Vibration eigenfrequencies of an elastic sphere with large radius

An estimation is given for the free vibration eigenfrequencies (normal modes) of a homogeneous solid sphere with a large radius, with application to Earth's free vibrations. The free vibration eigenfrequencies of a fluid sphere are also derived as a particular case. Various corrections arising from static and dynamic gravitation, rotation, and inhomegeneities are estimated, and a tentative notion of an earthquake temperature is introduced.

Introduction

The long-standing interest in the vibrations of a solid sphere is related to the seismic vibrations of the Earth [1,2]. After a relatively short burst of energy in an earthquake the Earth continues to vibrate freely for a long time. Though with a liquid outer core and a viscous mantle, the Earth is still approximated by a solid sphere. Great progress in studying the vibrations of a homogeneous and isotropic elastic sphere has been made since the early days when Lamb introduced the vector spherical harmonics (Hansen vectors) [3-5]. The relevant eigenfrequencies were computed numerically as early as 1898 [6]. We discuss in this paper a natural simplification of this problem, which arises from the fact that a large radius of the sphere is a natural cutoff. Apart from giving formally the general solution of vibrations generated by the seismic tensorial force, we show that a large radius simplifies appreciably the boundary conditions, leading readily to the estimation of the eigenfrequencies (normal modes). The particular case of a fluid sphere is treated to a larger extent.

As it is well known, the vibrations of the Earth following an earthquake are of great importance in revealing the inner structure of the crust, mantle, and even the inner cores of the Earth. These vibrations imply a large number of modes, usually classified as spheroidal and toroidal, with periods in a wide range from 10⁻³ − 10⁻⁴s to hours. They attenuate slowly in time, leading to the thermalization of the residual energy of the earthquake. Usually, they are studied numerically, from the recorded data. We give in this paper a thorough description of the vibrations of a sphere in the limit of a large radius, which simplifies greatly the problem. This simplification allows us to perform analytical calculations to a great extent. First, we present the eigenvibrations for a homogeneous and isotropic elastic solid sphere, with general boundary conditions. Second, we introduce the assumption of a large radius, as appropriate for Earth's vibrations. We show that the toroidal vibrations are easily amenable to analytical calculations, while useful quantitative estimations can be made for the spheroidal vibrations. As a useful example, we include the analysis of the vibrations of a fluid sphere. Further on, we investigate the effects of static and dynamic gravitation, a problem with a higher degree of difficulty. Also, the effects of the Coriolis and centrifugal forces are analyzed, with emphasis on their well-known frequency splitting. Finally, we discuss the possibility of estimating the temperature gained by the Earth, as a consequence of the thermalization energy released in an earthquake.

Solid sphere

The elastic vibrations of a homogeneous and isotropic solid are described by the equation

$\mu curlcurlu-\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}graddivu-\rho {\omega }^{2}u\mathrm{<mo>=</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\text{ (1)}$

where u is the local displacement, ρ is the density, µ and λ are the Lame elastic moduli, ω is the frequency, and F(ω) is the force [7]. The components of the seismic tensorial force are

$$F_i\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{M}_{ij}\mathrm{<mo\; stretchy="false">(</mo>}\omega \mathrm{<mo\; stretchy="false">)</mo>}{\partial }_j\delta \mathrm{<mo\; stretchy="false">(</mo>}r-r_0\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\text{ (2)}$$

where M_ij(ω) is the Fourier transform of the seismic moment, r₀ is the position of the point where the force is placed and i, j, … = 1,2,3 are cartesian labels [8,9]. An equivalent form of equation (1) is

$c_2^2\Delta u+\mathrm{<mo\; stretchy="false">(</mo>}{c}_1^2-c_2^2\mathrm{<mo\; stretchy="false">)</mo>}graddivu+{\omega }^{2}u\mathrm{<mo>=</mo>}-f\mathrm{,}\text{ (3)}$

where $c_1\mathrm{<mo>=</mo>}\sqrt{\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}\rho }$ is the velocity of the longitudinal elastic waves, $c_2\mathrm{<mo>=</mo>}\sqrt{\mu \mathrm{<mo>/</mo>}\rho }$ is the velocity of the transverse elastic waves, and f(ω) = F(ω) / ρ (also, m_ij(ω) = M_ij (ω) / ρ). As it is well known, equation (3) is separated into two inhomogeneous Helmholtz equations

$\begin{array}{c}{c}_1^2\Delta \Phi +{\omega }^2\Phi \mathrm{<mo>=</mo>}-\phi \mathrm{,}{c}_2^2\Delta A+{\omega }^{2}A\mathrm{<mo>=</mo>}-h\mathrm{,}\end{array}\text{ (4)}$

by $u\mathrm{<mo>=</mo>}grad\Phi +curlA\text{, }divA\mathrm{<mo>=</mo>0}$ , $f\mathrm{<mo>=</mo>}grad\phi +curlh\text{, }divh\mathrm{<mo>=</mo>0}$ , where ϕ and h are given by $\Delta \phi \mathrm{<mo>=</mo>}divf\text{, }\Delta h\mathrm{<mo>=</mo>}-curlf$ (Helmholtz potentials). We get

$$\begin{array}{c}\phi \mathrm{<mo>=</mo>}-\frac{1}{4\pi }{m}_{ij}{\partial }_i{\partial }_j\frac{1}{|r-r_0|}\mathrm{,}{h}_i\mathrm{<mo>=</mo>}\frac{1}{4\pi }{\epsilon }_{ijk}{m}_{kl}{\partial }_j{\partial }_l\frac{1}{|r-r_0|}\end{array}\text{ (5)}$$

(where εijk is the antisymmetric tensor of rank three), such that we are led to consider the equation

$c^2\Delta F+{\omega }^{2}F\mathrm{<mo>=</mo>}\frac{1}{r}\text{ (6)}$

with solution

$F\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{1-\cos kr}{{\omega }^{2}r}\mathrm{,}{k}^2\mathrm{<mo>=</mo>}{\omega }^2\mathrm{<mo>/</mo>}{c}^2\mathrm{<mo>;</mo>}\text{ (7)}$

This solution results immediately from the vibration Green function $G\mathrm{<mo>=</mo>}-\frac{\cos kr}{4\pi c^{2}r}$ of the Helmholtz equation $c^2\Delta G+{\omega }^{2}G\mathrm{<mo>=</mo>}\delta \mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>.}$ We get a particular solution of equation (3)

$$\begin{array}{c}{u}_i^p\mathrm{<mo>=</mo>}\frac{1}{4\pi }{m}_{ij}{\partial }_j\Delta F_2\mathrm{<mo\; stretchy="false">(</mo>}|r-r_0|\mathrm{<mo\; stretchy="false">)</mo>}+\\ \\ +\frac{1}{4\pi }{m}_{jk}{\partial }_i{\partial }_j{\partial }_k\left[F_1\mathrm{<mo\; stretchy="false">(</mo>}|r-r_0|\mathrm{<mo\; stretchy="false">)</mo>}-F_2\mathrm{<mo\; stretchy="false">(</mo>}|r-r_0|\mathrm{<mo\; stretchy="false">)</mo>}\right]\mathrm{.}\end{array}\text{ (8)}$$

For a fluid, where c₂ = 0 (µ = 0) and m_ij = −m_δij, this solution becomes

$$u^p\mathrm{<mo>=</mo>}-\frac{m}{4\pi c_1^2}grad\frac{\cos k_1|r-r_0|}{|r-r_0|}\mathrm{.}\text{ (9)}$$

To apply these results to a sphere we need to use expansions in a series of (orthogonal) vector spherical harmonics, defined by[10]

$\begin{array}{c}{R}_{lm}\mathrm{<mo>=</mo>}{Y}_{lm}{e}_r\mathrm{,}\\ \\ S_{lm}\mathrm{<mo>=</mo>}\frac{\partial Y_{lm}}{\partial \theta }{e}_{\theta }+\frac{1}{\sin \theta }\frac{\partial Y_{lm}}{\partial \phi }{e}_{\phi }\mathrm{,}\\ \\ T_{lm}\mathrm{<mo>=</mo>}\frac{1}{\sin \theta }\frac{\partial Y_{lm}}{\partial \phi }{e}_{\theta }-\frac{\partial Y_{lm}}{\partial \theta }{e}_{\phi }\mathrm{,}\end{array}\text{ (10)}$

l ≠ 0, where Y_lm are spherical harmonics and e_r,θ,ϕ are the spherical unit vectors. The functions R_lm S_lm are called spheroidal functions, while the functions T_lm are called toroidal functions. The series expansion reads

$u^p\mathrm{<mo>=</mo>}{\displaystyle \sum _{lm}}\mathrm{<mo\; stretchy="false">(</mo>}{f}_{lm}^p{R}_{lm}+g_{lm}^p{S}_{lm}+h_{lm}^p{T}_{lm}\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\text{ (11)}$

where $f_{lm}^p$ , $g_{lm}^p$ and $h_{lm}^p$ are functions only of the radius r. A similar series holds also for the free solution uf of equation (1).

The explicit form of the coefficients $f_{lm}^p$ , $g_{lm}^p$ and $h_{lm}^p$ is extremely cumbersome. We prefer to work formally with equation (1) and series expansions of the full solution u = u^p + u^f and the force F(ω), with coefficients f_lm, g_lm, and h_lm and F^r,s,t, respectively. Making use of such series expansions and the properties of the vector spherical harmonics,[10] we get the equations

$\begin{array}{c}{f}^{{\text{'}}^{\prime }}+\frac{2}{r}{f}^{\mathrm{<mo>\text{'}</mo>}}+\frac{\rho {\omega }^2}{\lambda +2\mu }f-\left[2+\frac{\mu l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{\lambda +2\mu }\right]\frac{1}{r^2}f+\\ \\ +\frac{\mathrm{<mo\; stretchy="false">(</mo>}\lambda +3\mu \mathrm{<mo\; stretchy="false">)</mo>}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}{r}^2}g-\frac{\mathrm{<mo\; stretchy="false">(</mo>}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}r}{g}^{\mathrm{<mo>\text{'}</mo>}}\mathrm{<mo>=</mo>}-\frac{F^r}{\lambda +2\mu }\mathrm{,}\\ \\ g^{{\text{'}}^{\prime }}+\frac{2}{r}{g}^{\mathrm{<mo>\text{'}</mo>}}+\frac{\rho {\omega }^2}{\mu }g-\frac{\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{\mu r^2}g+\\ \\ +\frac{\mathrm{2<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}}{\mu r^2}f+\frac{\lambda +\mu }{\mu r}{f}^{\mathrm{<mo>\text{'}</mo>}}\mathrm{<mo>=</mo>}-\frac{F^s}{\mu }\mathrm{,}\\ \\ h^{{\text{'}}^{\prime }}+\frac{2}{r}{h}^{\mathrm{<mo>\text{'}</mo>}}+\frac{\rho {\omega }^2}{\mu }h-\frac{l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{r^2}h\mathrm{<mo>=</mo>}-\frac{F^t}{\mu }\mathrm{,}\end{array}\text{ (12)}$

where, for the sake of simplicity, we dropped out the suffixes lm.

We turn now to the boundary conditions. The force P acting (inwards) on the surface r = R of the sphere, where R is the radius of the sphere, with the spherical components P_α = (α = r, θ ,ϕ) is P_α = nβσαβ = σ_αr, where the stress tensor is given by ${\sigma }_{\alpha \beta }\mathrm{<mo>=</mo>2}\mu u_{\alpha \beta }+\lambda u_{\gamma \gamma }{\delta }_{\alpha \beta }$ ; we get

$\begin{array}{c}2\mu u_{\theta r}\mathrm{<mo>=</mo>}{P}_{\theta }\mathrm{,}2\mu u_{\phi r}\mathrm{<mo>=</mo>}{P}_{\phi }\mathrm{,}\\ \\ 2\mu u_{rr}+\lambda divu\mathrm{<mo>=</mo>}{P}_r\mathrm{,}\end{array}\text{ (13)}$

where divu is written in spherical coordinates,

$divu\mathrm{<mo>=</mo>}{\displaystyle \sum _{lm}}\left[\frac{1}{r^2}\frac{d}{dr}\mathrm{<mo\; stretchy="false">(</mo>}{r}^2{f}_{lm}\mathrm{<mo\; stretchy="false">)</mo>}-\frac{g_{lm}}{r}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\right]Y_{lm}\text{ (14)}$

(by using the properties of the vector spherical harmonics equations[10]). We compute the strain tensor u_αβ in spherical coordinates[7]

$\begin{array}{c}{u}_{rr}\mathrm{<mo>=</mo>}\frac{\partial u_r}{\partial r}\mathrm{,}{u}_{\theta \theta }\mathrm{<mo>=</mo>}\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_r}{r}\mathrm{,}{u}_{\phi \phi }\mathrm{<mo>=</mo>}\frac{1}{r\sin \theta }\frac{\partial u_{\phi }}{\partial \phi }+\frac{u_{\theta }}{r}\cot \theta +\frac{u_r}{r}\mathrm{,}\\ \\ 2u_{\theta \phi }\mathrm{<mo>=</mo>}\frac{1}{r}\left(\frac{\partial u_{\phi }}{\partial \theta }-u_{\phi }\cot \theta \right)+\frac{1}{r\sin \theta }\frac{\partial u_{\theta }}{\partial \phi }\mathrm{,}2u_{r\theta }\mathrm{<mo>=</mo>}\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}+\frac{1}{r}\frac{\partial u_r}{\partial \theta }\mathrm{,}\\ \\ 2u_{\phi r}\mathrm{<mo>=</mo>}\frac{1}{r\sin \theta }\frac{\partial u_r}{\partial \phi }+\frac{\partial u_{\phi }}{\partial r}-\frac{u_{\phi }}{r}\end{array}\text{ (15)}$

by using the spherical components

$\begin{array}{c}{u}_r\mathrm{<mo>=</mo>}{\displaystyle \sum _{lm}}{f}_{lm}{Y}_{lm}\mathrm{,}\\ \\ u_{\theta }\mathrm{<mo>=</mo>}{\displaystyle \sum _{lm}}{g}_{lm}\frac{\partial Y_{lm}}{\partial \theta }+{\displaystyle \sum _{lm}}\frac{h_{lm}}{\sin \theta }\frac{\partial Y_{lm}}{\partial \phi }\mathrm{,}\\ \\ u_{\phi }\mathrm{<mo>=</mo>}{\displaystyle \sum _{lm}}\frac{g_{lm}}{\sin \theta }\frac{\partial Y_{lm}}{\partial \phi }-{\displaystyle \sum _{lm}}{h}_{lm}\frac{\partial Y_{lm}}{\partial \theta }\end{array}\text{ (16)}$

of the expansion of the displacement vector and the definition of the vector spherical functions (equations (10)). Similarly, we decompose the force P in vector spherical harmonics (with coefficients P^r,s,t) and identify its spherical components. The boundary conditions given by equations (13) lead to

$\begin{array}{c}2\mu f^{\mathrm{<mo>\text{'}</mo>}}+\lambda \left[\frac{2}{r}f+f^{\mathrm{<mo>\text{'}</mo>}}-\frac{g}{r}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\right]{|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>}{P}^r\mathrm{,}\\ \\ \mu \left(\frac{g}{r}-g^{\mathrm{<mo>\text{'}</mo>}}-\frac{f}{r}\right){|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>}-P^s\mathrm{,}\\ \\ \mu \left(\frac{h}{r}-h^{\mathrm{<mo>\text{'}</mo>}}\right){|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>}-P^t\mathrm{,}\end{array}\text{ (17)}$

where we dropped the subscripts lm.

Vibration eigenfrequencies for large radius

The solutions f, g and h of equations (12) consist of free solutions (solutions of the homogeneous equations (12)) plus particular solutions. The homogeneous third equation (12), which describes toroidal vibrations, is the equation of the spherical Bessel functions $j_l\mathrm{<mo\; stretchy="false">(</mo>}kr\mathrm{<mo\; stretchy="false">)</mo>},k\mathrm{<mo>=</mo>}\sqrt{\rho {\omega }^2\mathrm{<mo>/</mo>}\mu }\mathrm{<mo>=</mo>}\omega \mathrm{<mo>/</mo>}{c}_2$ . For F^t = 0 and P^t = 0 (a free surface) the third equation in the boundary conditions (17) gives

$j_l\mathrm{<mo\; stretchy="false">(</mo>}kR\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}kR{j}_l^{\mathrm{<mo>\text{'}</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}kR\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>;</mo>}\text{ (18)}$

this equation has an infinity of solutions β_ln, labeled by an integer n, such that we get the eigenfrequencies

${\omega }_{ln}\mathrm{<mo>=</mo>}\frac{c_2}{R}{\beta }_{ln}\mathrm{.}\text{ (19)}$

We can get an approximate estimate of the numbers β_ln by using the asymptotic expression of the spherical Bessel functions[11]

$j_l\mathrm{<mo\; stretchy="false">(</mo>}kr\mathrm{<mo\; stretchy="false">)</mo>}\simeq \frac{1}{kr}\cos \left[kr-\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\frac{\pi }{2}\right]\mathrm{,}kr\gg 1\mathrm{<mo>;</mo>}\text{ (20)}$

For kR ≫ 1 equation (18) becomes

$\tan \left[kR-\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\frac{\pi }{2}\right]\mathrm{<mo>=</mo>}-\frac{2}{kR}\mathrm{,}\text{ (21)}$

which have the approximate zeroes

${\beta }_{ln}\simeq n\pi +\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\frac{\pi }{2}\mathrm{,}\text{ (22)}$

where n is any (large) integer. We can see that the frequencies are dense for large $R\left(\Delta {\omega }_{ln}\mathrm{<mo>=</mo>}\pi c_2\mathrm{<mo>/</mo>}R\right)$ . The free toroidal solution is a superposition of $j_l\mathrm{<mo\; stretchy="false">(</mo>}{k}_{ln}r\mathrm{<mo\; stretchy="false">)</mo>}$ , where $k_{ln}\mathrm{<mo>=</mo>}{\omega }_{ln}\mathrm{<mo>/</mo>}{c}_2\mathrm{<mo>=</mo>}{\beta }_{ln}\mathrm{<mo>/</mo>}R$ , with undetermined coefficients.

In general, for F^t ≠ 0 and P^t ≠ 0 the free toroidal solution is C_ljl(kr), where the constants C_l are determined from the boundary condition. It is easy to see that these coefficients include singular factors proportional to $\sim \frac{1}{\omega -{\omega }_{ln}}$ , such that, the integration over frequencies leads to toroidal vibrations governed by the eigenfrequencies ω_ln. In general, the solution h_lm depends on two integration constants, which are determined by the boundary condition and the condition of a finite solution at the origin.

We pass now to the spheroidal components which involve the functions f and g in equations (12) and (17). We note that the two coupled equations (12) of for the functions f and g include Bessel operators for spherical Bessel functions. We can get a simplified picture of these equations for large values of r. Indeed, it is easy to see that in the limit wr / c_1,2 ≫ l² the free solutions are

$$\begin{array}{c}f\simeq \frac{A}{r}\cos \mathrm{<mo\; stretchy="false">(</mo>}\omega r\mathrm{<mo>/</mo>}{c}_1+{\phi }_{1l}\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}g\simeq \frac{B}{r}\cos \mathrm{<mo\; stretchy="false">(</mo>}\omega r\mathrm{<mo>/</mo>}{c}_2+{\phi }_{2l}\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\end{array}\text{ (23)}$$

where the coefficients A and B and the phases ϕ_1l,2l remain undetermined. In this limit the boundary conditions are $f^{\prime }{|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>}{g}^{\prime }{|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>0}$ and the eigenfrequencies are given by $${\omega }_{nl}R\mathrm{<mo>/</mo>}{c}_{\mathrm{1,2}}+{\phi }_{1l\mathrm{,2}l}\mathrm{<mo>=</mo>}n\pi$$ , where n is any integer (the roots of the equation $$\sin \mathrm{<mo\; stretchy="false">(</mo>}\omega R\mathrm{<mo>/</mo>}{c}_{\mathrm{1,2}}+{\phi }_{1l\mathrm{,2}l}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0}$$ ).

The condition wr / c_1,2 ≫ l² is satisfied for a large r and a reasonably large range of frequencies and parameter l. For instance, for r close to the Earth's surface, which is the spatial region of interest, we get wR / c of the order ≃10³w for the mean radius of the Earth R = 6370 km and the mean velocity of the elastic waves c = 5 km / s. Indeed, frequencies as low as w = 10^-3 s^-1 are known for Earth's seismic vibrations [12-14].

We can see that there are two branches of spheroidal eigenfrequencies (corresponding to the velocities c_1,2), which are dense (continuous) for large R, very similar to the infinite space (as expected for large R); the ω⁽²⁾ -branch, although close to the toroidal branch, is distinct (there is a total of three branches of eigenfrequencies, corresponding to the three degrees of freedom; in the limit of the rotations of the sphere as a whole their frequencies go to zero (acoustic modes)). For non-vanishing forces, we have spheroidal vibrations driven by these forces, as discussed in the previous cases. The set of all eigenfrequencies is called the (seismic) spectrum. Earth's eigenmodes with eigenfrequencies of the order 10⁻³ − 10⁻⁴ s⁻¹, excited by earthquakes, have been discussed in Refs [12-14].

The numerical solution of equations (12) indicates that the lowest mode (the fundamental mode) is S_lm with l = 2 and n = 0 (therefore, we may denote it as) $S_{l\mathrm{<mo>=</mo>2,}m}^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo>=</mo>0<mo\; stretchy="false">)</mo>}}$ ;[6] it is denoted by ₀ S₂, and its eigenfrequency is denoted ω₂₀; the corresponding period is approximately 1 an hour. Much later, the Earth's crust was modeled as a series of superposed layers, with welded interfaces; the vibrations of such a stack of layers can be computed and long periods of the fundamental modes have been obtained; the dispersion relation of these modes (i.e., the dependence of the frequency on their label n) can give information about the inner crustal structure [15,16]. The first observation of " free oscillations of the Earth as a whole" was made for the Kamchatka earthquake of November 4, 1952;[17] they were followed by many observations of the Earth's vibrations caused by the great Chile earthquake of May 22, 1960 [18-20] (with magnitude greater than 8, which saturated the scales[21]). Today, eigenoscillations of the Earth can be recorded even for small earthquakes [22].

From studies of propagation of the seismic waves, it was inferred the Earth's solid inner core [23,24] of ≃1000km and the outer liquid core of radius radius ≃2000km. The inner-outer core discontinuity is called the Bullen, or Lehmann, discontinuity. The temperature of the inner core is ≃3000K (iron and nickel) and the pressure is ≃1012 dyn /cm2. The buoyancy at this boundary could be the source of convection currents that generate the Earth's magnetic field (geodynamo effect). The next layers are a viscous mantle of thickness ≃3000km and the solid crust of thickness≃70km. The boundary between mantle and crust is known as the Mohorovic discontinuity.

Fluid sphere

For a fluid sphere, the shear modulus µ is zero (µ= 0); equations (12) become

$\begin{array}{c}{f}^{{\text{'}}^{\prime }}+\frac{2}{r}{f}^{\mathrm{<mo>\text{'}</mo>}}+k^{2}f-\frac{2}{r^2}f-\frac{d}{dr}\left[\frac{l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}g}{r}\right]\mathrm{<mo>=</mo>}-\frac{F^r}{\lambda }\\ \\ \frac{1}{r}{f}^{\mathrm{<mo>\text{'}</mo>}}+\frac{2}{r^2}f-\frac{l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{r^2}g+k^{2}g\mathrm{<mo>=</mo>}-\frac{F^s}{\lambda }\mathrm{,}\end{array}\text{ (24)}$

where k² = ρω² / λ = ω² / c²; the boundary condition reads

$\left[\frac{2}{r}f+f^{\mathrm{<mo>\text{'}</mo>}}-\frac{g}{r}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\right]{|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>}\frac{P^r}{\lambda }\mathrm{.}\text{ (25)}$

Let us introduce divu, given by equation (14), which includes

$d\mathrm{<mo>=</mo>}{f}^{\mathrm{<mo>\text{'}</mo>}}+\frac{2}{r}f-\frac{g}{r}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\mathrm{.}\text{ (26)}$

Then the boundary condition becomes

$d{|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>}\frac{P^r}{\lambda }\mathrm{,}\text{ (27)}$

the second equation (24) reads

$\frac{d}{r}+k^{2}g\mathrm{<mo>=</mo>}-\frac{F^s}{\lambda }\text{ (28)}$

and the first equation (24) is

$d^{\mathrm{<mo>\text{'}</mo>}}+k^{2}f\mathrm{<mo>=</mo>}-\frac{F^r}{\lambda }\mathrm{.}\text{ (29)}$

Hence, we have

$g\mathrm{<mo>=</mo>}-\frac{d}{k^{2}r}-\frac{F^s}{\lambda k^2}\mathrm{,}f\mathrm{<mo>=</mo>}-\frac{d^{\mathrm{<mo>\text{'}</mo>}}}{k^2}-\frac{F^r}{\lambda k^2}\mathrm{.}\text{ (30)}$

Now we introduce these functions in equation (26) and get

$d^{{\text{'}}^{\prime }}+\frac{2d^{\mathrm{<mo>\text{'}</mo>}}}{r}+k^{2}d-\frac{l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{r^2}d\mathrm{<mo>=</mo>}-\frac{{\mathrm{<mo\; stretchy="false">(</mo>}{F}^r\mathrm{<mo\; stretchy="false">)</mo>}}^{\mathrm{<mo>\text{'}</mo>}}}{\lambda }-\frac{2F^r}{\lambda r}+\frac{F^s}{\lambda r}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\mathrm{.}\text{ (31)}$

For free vibrations this is the Bessel equation for spherical Bessel functions $d\mathrm{<mo>=</mo>}{j}_l\mathrm{<mo\; stretchy="false">(</mo>}kr\mathrm{<mo\; stretchy="false">)</mo>}$ ; the boundary condition (27) leads to the eigenfrequencies ${\omega }_{ln}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}c\mathrm{<mo>/</mo>}R\mathrm{<mo\; stretchy="false">)</mo>}{\beta }_{ln}$ , $j_l\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_{ln}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0}$ . In a fluid, we have only pressure p, and the stress tensor is ${\sigma }_{ij}\mathrm{<mo>=</mo>}-p{\delta }_{ij}$ ( ${\sigma }_{ij}\mathrm{<mo>=</mo>2}\mu u_{ij}+\lambda u_{kk}{\delta }_{ij}$ with $\mu \mathrm{<mo>=</mo>0}$ ); therefore, for a fluid $p\mathrm{<mo>=</mo>}-\lambda u_{ii}\mathrm{<mo>=</mo>}-\lambda divu$ ; the equations written above for d are in fact equations for the pressure p. It is convenient to introduce the decomposition in Helmholtz potentials $u\mathrm{<mo>=</mo>}grad\Phi +curlA\mathrm{,}$ $divA\mathrm{<mo>=</mo>0}$ and $F\mathrm{<mo>=</mo>}grad\phi +curlh$ $divh\mathrm{<mo>=</mo>0}$ , ; then, $p\mathrm{<mo>=</mo>}-\lambda \Delta \Phi$ and the equation of motion $\rho \ddot{u}\mathrm{<mo>=</mo>}\lambda grad\cdot divu+F\mathrm{<mo>=</mo>}-gradp+F$ becomes $\rho \ddot{\Phi }\mathrm{<mo>=</mo>}\lambda \Delta \Phi +\phi$ , where the potential $\phi$ is given by $\Delta \phi \mathrm{<mo>=</mo>}divF$ and $h\mathrm{<mo>=</mo>0}$ , $A\mathrm{<mo>=</mo>0}$ . For vibrations this equation reads $c^2\Delta \Phi +{\omega }^2\Phi \mathrm{<mo>=</mo>}-\frac{1}{\rho }\phi$ and for $F\mathrm{<mo>=</mo>}-Mgrad\delta \mathrm{<mo\; stretchy="false">(</mo>}r-r_0\mathrm{<mo\; stretchy="false">)</mo>}$ we get $\phi \mathrm{<mo>=</mo>}-M\delta \mathrm{<mo\; stretchy="false">(</mo>}r-r_0\mathrm{<mo\; stretchy="false">)</mo>}$ , $\Phi \mathrm{<mo>=</mo>}-m\frac{\cos k|r-r_0|}{4\pi c^2|r-r_0|}$ ( $m\mathrm{<mo>=</mo>}M\mathrm{<mo>/</mo>}\rho$ ) and the solution $$u^p$$ given by equation (9).

Static self-gravitation

A gravitational force

$FdV\mathrm{<mo>=</mo>}G\rho dV\frac{m}{r^2}\mathrm{<mo>=</mo>}\frac{4\pi }{3}G{\rho }^{2}rdV\text{ (32)}$

acts upon a volume element dV placed at a distance r from the center of a sphere, where G = 6.67 ´ 10-8 cm3/ g × s2 is the universal constant of gravitation, r is the density of the sphere (assumed incompressible), and m = (4p / 3)rr3 is the mass of the sphere with radius r. If the sphere is compressible, the gravitational potential j is given by the Poisson equation Dj = 4pGr and the gravitational force per unit mass is F = -gradj; the condition of (hydrostatic) equilibrium (for a non-rotating sphere) reads gradr = rF = -rgradj, such that div[(gradr)/r] = -4pGr; the dependence of the pressure on the density is given by the equation of state; for a constant density the pressure for a self-gravitating sphere of radius R at rest with free surface is p = (2p/3) Gr2(R2-r2) (it seems that the pressure in the inner Earth's (solid) core is ≃300GPa = 3 ´ 1012dyn/cm2). Making use of equation (32), the equation of the elastic motion reads

$\rho \ddot{u}-\mu \Delta u-\mathrm{<mo\; stretchy="false">(</mo>}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}graddivu\mathrm{<mo>=</mo>}F\mathrm{<mo>=</mo>}-\gamma r\mathrm{,}\text{ (33)}$

where $\gamma \mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>4}\pi \mathrm{<mo>/</mo>3<mo\; stretchy="false">)</mo>}G{\rho }^2$ . Since $Y_{00}\mathrm{<mo>=</mo>1<mo>/</mo>}\sqrt{4\pi }$ we may write

$F\mathrm{<mo>=</mo>}-\gamma r\mathrm{<mo>=</mo>}-\sqrt{4\pi }\gamma r{Y}_{00}{e}_r\mathrm{,}\text{ (34)}$

whence we can see that F has an expansion is series of spheroidal and toroidal functions with all the coefficients zero, except the coefficient $F_{00}^r\mathrm{<mo>=</mo>}-\sqrt{4\pi }\gamma r$ of the function ₀₀; it follows that the motion may include all the eigenmodes S_lm and T_lm, as well as all the eigenmodes Rlm, the latter with l ≠ 0; for l = 0 , m = 0 the motion, described by f = f₀₀, is driven by the gravitational force. We note also that the force in equation (33) is static, which means that its Fourier transform is proportional to δ(ω). For l = 0 the first equation (12) includes only the function f, i.e. f δ(ω); this equation reads

$\begin{array}{c}{f}^{{\text{'}}^{\prime }}+\frac{2}{r}{f}^{\mathrm{<mo>\text{'}</mo>}}-\frac{2}{r^2}f\mathrm{<mo>=</mo>}\frac{\sqrt{4\pi }\gamma }{\lambda +2\mu }r\mathrm{.}\end{array}\text{ (35)}$

It is easy to see that a particular solution of this equation is $\mathrm{<mo\; stretchy="false">[</mo>}\sqrt{4\pi }\gamma \mathrm{<mo>/</mo>10<mo\; stretchy="false">(</mo>2}\mu +\lambda \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}{r}^3$ , while the homogeneous part of this equation has the solution C₁r +C₂/ r², where C_1,2 there are constants of integration; we must take C₂ = 0 , because the solution is finite at the origin. We are left with the solution

$u_r\mathrm{<mo>=</mo>}A{r}^3+C_{1}r\mathrm{,}A\mathrm{<mo>=</mo>}\frac{\gamma }{\mathrm{10<mo\; stretchy="false">(</mo>2}\mu +\lambda \mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{.}\text{ (36)}$

This solution must satisfy the boundary conditions at the surface of the sphere; making use of equations (16), we have the strain tensor $u_{rr}\mathrm{<mo>=</mo>}{u}_r^{\mathrm{<mo>\text{'}</mo>}}$ and $u_{\theta \theta }\mathrm{<mo>=</mo>}{u}_{\phi \phi }\mathrm{<mo>=</mo>}{u}_r\mathrm{<mo>/</mo>}r$ ; the force on the surface is $-{\sigma }_{\alpha r}{|}_R$ , where the stress tensor is given by ${\sigma }_{\alpha \beta }\mathrm{<mo>=</mo>2}\mu u_{\alpha \beta }+\lambda u_{\gamma \gamma }{\delta }_{\alpha \beta }$ ; for a free surface we get the boundary condition

$\mathrm{<mo\; stretchy="false">(</mo>2}\mu +\lambda \mathrm{<mo\; stretchy="false">)</mo>}{u}_r^{\mathrm{<mo>\text{'}</mo>}}+2\lambda \frac{u_r}{r}{|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>0}\text{ (37)}$

( ${\sigma }_{\alpha r}{|}_R\mathrm{<mo>=</mo>0}$ ), whence we determine the constant $C_1\mathrm{<mo>=</mo>}-\mathrm{<mo\; stretchy="false">[</mo><mo\; stretchy="false">(</mo>6}\mu +5\lambda \mathrm{<mo\; stretchy="false">)</mo><mo>/</mo><mo\; stretchy="false">(</mo>2}\mu +3\lambda \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}A{R}^2$ and, finally, the radial displacement

$u_r\mathrm{<mo>=</mo>}Ar\left(r^2-\frac{6\mu +5\lambda }{2\mu +3\lambda }{R}^2\right)\mathrm{<mo>=</mo>}\frac{\gamma }{\mathrm{10<mo\; stretchy="false">(</mo>2}\mu +\lambda \mathrm{<mo\; stretchy="false">)</mo>}}r\left(r^2-\frac{6\mu +5\lambda }{2\mu +3\lambda }{R}^2\right)\mathrm{<mo>;</mo>}\text{ (38)}$

we note that the radial displacement u, is negative, as expected. It is worth estimating the radial displacement at the surface due to gravitation

$u_r{|}_{r\mathrm{<mo>=</mo>}R}\mathrm{<mo>=</mo>}-\frac{\gamma }{\mathrm{5<mo\; stretchy="false">(</mo>2}\mu +3\lambda \mathrm{<mo\; stretchy="false">)</mo>}}{R}^3\mathrm{<mo>;</mo>}\text{ (39)}$

making use of ρ = 5g /cm³, λ,µ ≃ 10¹¹dyn/cm² (parameters for Earth), we get γ ≃ 10⁻⁶ g/cm³s² and $u_r{|}_R\simeq {10}^{-18}{R}^{3}cm\simeq {10}^{8}cm{\mathrm{<mo>=</mo>10}}^{3}km$ , for the Earth's radius R ≃ 6 ×10⁸ cm; this is a distance of the order of the Earth's radius. Moreover, the strain is of the order 1/6, which may cast doubts on the validity of the linear elasticity used in this estimation. In addition, we note that the density suffers an important change due to the static gravitational field. Indeed, the change in density is dr = -div(ru) = -r0divu, where r0 is the uniform initial density; with u given by equation (38) we get

$$\frac{\delta \rho }{{\rho }_0}\mathrm{<mo>=</mo>}A\mathrm{<mo\; stretchy="false">(</mo>3}\alpha R^2-5r^2\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\alpha \mathrm{<mo>=</mo>}\frac{6\mu +5\lambda }{2\mu +3\lambda }\mathrm{,}\text{ (40)}$$

which is of the order unity. The proper estimation of the static effect of the self-gravitational field on the elastic sphere is to solve simultaneously the equation of elastic equilibrium (33) with F = -rgradj and the Poisson equation for the gravitational field j, Dj =4pGr. With spherical symmetry we have

$F\mathrm{<mo>=</mo>}-\frac{4\pi }{3r^2}G\rho {\displaystyle {\int }_{\mathrm{0<mo><</mo>}{r}^{\mathrm{<mo>\text{'}</mo>}}\mathrm{<mo><</mo>}r}}d{r}^{\mathrm{<mo>\text{'}</mo>}}\rho \frac{r}{r}\mathrm{<mo>;</mo>}\text{ (41)}$

the Poisson equation for the gravitational potential may be written as D(F/r) = -4pGgradr, such that the problem involves two equations and unknowns, u and r. Since this is a more difficult problem it is preferable to consider the density r as an empirically known function of r (a parametrization in powers of r can be used for r and a variational approach can be applied to the problem). Even so, the equations governing the influence of the gravitational field upon the elasticity of a self-gravitating sphere are difficult.

Dynamic self-gravitation

Let us assume a spheric, non-rotating, homogeneous, elastic Earth at equilibrium under the action of its gravitational field; we consider small elastic deformations of this equilibrium state; in the first approximation, we have a small change denoted by K in the gravitational potential as a consequence of the small changes in density div(ρu), i.e., we have

$\Delta K\mathrm{<mo>=</mo>}-4\pi Gdiv\mathrm{<mo\; stretchy="false">(</mo>}\rho u\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\text{ (42)}$

where ρ is a known function of r. The equation of elastic motion reads

$\rho \ddot{u}-\mu \Delta u-\mathrm{<mo\; stretchy="false">(</mo>}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}graddivu\mathrm{<mo>=</mo>}-\rho gradK\mathrm{.}\text{ (43)}$

These two coupled (vectorial) equations are difficult to be treated by an analytical method, due to the non-uniformity of the density. For a uniform density, taking the div in equation (43) and using equation (42) we get for D = divu

$\rho \ddot{D}-\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}\Delta D\mathrm{<mo>=</mo>4}\pi G{\rho }^{2}D\mathrm{,}\text{ (44)}$

an equation which indicates that the frequency ω changes by

$\Delta \mathrm{<mo\; stretchy="false">(</mo>}{\omega }^2\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}-4\pi G\rho \mathrm{<mo>;</mo>}\text{ (45)}$

for frequencies as low as ω = 10⁻⁴ S⁻¹ the variation given by equation (45) is large. Let us use the Helmholtz decomposition $u\mathrm{<mo>=</mo>}grad\Phi +curlA$ divA = 0; then, from equation (42) we have K = −4πGρΦ and from equation (43) we get $\Delta \Phi +k_1^2\Phi \mathrm{<mo>=</mo>0},\Delta A+k_2^{2}A\mathrm{<mo>=</mo>0}$ . These are the same equations as those which hold in the absence of the gravitational field, except that $k_1^2$ is changed into $k_1^2\to k_1^2+4\pi \rho G\mathrm{<mo>/</mo>}{c}_1^2$ . Moreover, we can see that only the spheroidal modes are affected by gravitation (since divT_lm = 0). It follows that the spheroidal frequencies (i.e., the branches ω^(1,2)) are given by the same relations of the type ω = (c/ R)β, where β denote the zeros the spherical Bessel functions in the limit of largeR; for c =c1, this relation reads ${\omega }^2+4\pi \rho G\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{c}_1^2\mathrm{<mo>/</mo>}{R}^2\mathrm{<mo\; stretchy="false">)</mo>}{\beta }^2$ . Hence, we may see that we should have the inequality $\mathrm{<mo\; stretchy="false">(</mo>}{c}_1^2\mathrm{<mo>/</mo>}{R}^2\mathrm{<mo\; stretchy="false">)</mo>}{\beta }^2\mathrm{<mo>></mo>4}\pi \rho G$ , or $\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}{\beta }^2\mathrm{<mo>></mo>4}\pi {\rho }^{2}G{R}^2$ . The term on the right side of this inequality is, up to an immaterial numerical factor, the pressure due to the gravitation at the origin; it is much larger than the elastic pressure λ + 2µ. The inequality is not satisfied for small values of β (as required by experimental observations). It follows that the model of an elastic solid Earth is not valid for the interior of the Earth. In those central regions, the elasticity is not able to sustain the gravitational pressure. Likely, an additional pressure exists there, which compensates for the gravitational pressure. The large dimensions of the mantle and liquid outer core complicate the matter, and such an Earth's model may exhibit very low frequencies (undertones);[25] If so, we may leave aside the effects of the gravitation in estimating the elastic vibrations of the Earth. In this case, c= 5km/s we get a period T ≃ (2.2/β) of hours; the smallest zero of j2 (corresponding approximately to the mode 0 S2) is β= 3.6;[11] we get T ≃ 37 minutes (for a velocity c = 3km/s the period is T = 61 minutes, which agrees with the experimental observations).

Rotation effect

If a vector a rotates, its change is $\delta a+\alpha \times a$ , where $\delta \alpha$ is the infinitesimal rotation angle; therefore, its velocity is $$\dot{a}+\Omega \times a$$ , where Ω is the angular velocity; its acceleration is $\ddot{a}+\dot{\Omega }\times a+2\Omega \times \dot{a}+\Omega \times \mathrm{<mo\; stretchy="false">(</mo>}\Omega \times a\mathrm{<mo\; stretchy="false">)</mo>}$ . Let us apply this relation to the displaced position a = r + u; we get the acceleration $\ddot{u}+\dot{\Omega }\times (r+u\mathrm{<mo\; stretchy="false">)</mo>}+2\Omega \times \dot{u}+\Omega \times \mathrm{<mo\; stretchy="false">[</mo>}\Omega \times (r+u\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}$ ; we can see that additional forces appear in rotation: is the Coriolis acceleration and $$-\Omega \times \mathrm{<mo\; stretchy="false">[</mo>}\Omega \times (r+u\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}$$ is the centrifugal acceleration. The Earth rotates with a constant angular velocity Ω=2π/T, T = 24 hours, oriented along the z -axis. We write the equation of elastic motion as

$\rho \ddot{u}+2\rho \Omega \times \dot{u}\mathrm{<mo>=</mo>}F\mathrm{,}\text{ (46)}$

where F includes the elastic force (i.e., F_i = ∂jσij) and other external forces, and the centrifugal force is omitted since Ω is much smaller than the eigenfrequencies of the Earth (an estimation of the longest periods of the Earth's eigenmodes gives an order of magnitude 2πR/cβ≃37 minutes, for the wave velocity c = 5km/s and β = 3.6 where R is the Earth's radius).

In the absence of the Coriolis force in equation (46) we decompose the force F and the displacement u in normal modes by using the spheroidal and toroidal functions. Let us focus on one normal mode, for instance, a toroidal mode $u_{lm}^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{h}_l^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}{T}_{lm}$ , corresponding to the eigenfrequency ω_ln = (c₂ / R)β_ln, where β_ln is, approximately, a zero of the function $j_l\mathrm{<mo\; stretchy="false">(</mo>}{k}^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}R\mathrm{<mo\; stretchy="false">)</mo>}$ ; the eigenfunctions $h_l^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}$ are given by the spherical Bessel functions $j_l\mathrm{<mo\; stretchy="false">(</mo>}{k}^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}r\mathrm{<mo\; stretchy="false">)</mo>}$ ; it is preferable to multiply these functions by constants and fix these constants such as

${\displaystyle \int }dr\cdot r^2{h}_l^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{h}_l^{\mathrm{<mo\; stretchy="false">(</mo>}{n}^{\mathrm{<mo>\text{'}</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\delta }_{n{n}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{<mo>;</mo>}\text{ (47)}$

we recall that the toroidal functions are orthogonal, i.e.

${\displaystyle \int }do{T}_{lm}{T}_{l^{\mathrm{<mo>\text{'}</mo>}}{m}^{\mathrm{<mo>\text{'}</mo>}}}^{\mathrm{<mo>*</mo>}}\mathrm{<mo>=</mo>}{\delta }_{l{l}^{\mathrm{<mo>\text{'}</mo>}}}{\delta }_{m{m}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{.}\text{ (48)}$

Since Ω/ω_{ln << 1} we solve equation (46) by a perturbation-theory method. First, we drop the labels l,m and n and use the notations $u_{lm}^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{u}_0$ , ω_lm = ω₀; we seek the solution as a series in powers of Ω/ω0

$u\mathrm{<mo>=</mo>}{u}_0+\frac{\Omega }{{\omega }_0}{u}_1+\mathrm{...}\mathrm{,}\text{ (49)}$

where u₁ to be determined, is assumed orthogonal on u₀,[26] with respect to the scalar product defined as the integration over the whole space, i.e.

${\displaystyle \int }dr{u}_1{u}_0\mathrm{<mo>=</mo>0}\mathrm{.}\text{ (50)}$

A similar series is valid for the frequency

$\omega \mathrm{<mo>=</mo>}{\omega }_0+\frac{\Omega }{{\omega }_0}{\omega }_1+\mathrm{...}\mathrm{.}\text{ (51)}$

Introducing these series in equation (46), with time Fourier transforms, we get

$\begin{array}{c}-\rho {\omega }_0^2{u}_0\mathrm{<mo>=</mo>}F\mathrm{,}\\ \\ -\rho {\omega }_0\Omega u_1-2\rho \Omega {\omega }_1{u}_0-2i\rho {\omega }_0\Omega \times u_0\mathrm{<mo>=</mo>0}\mathrm{<mo>;</mo>}\end{array}\text{ (52)}$

the first equation (52) defines the function u₀; in the second equation (52) we take the scalar product with u₀ and use the orthogonality of u₀ with u₁ we get

${\omega }_1\mathrm{<mo>=</mo>}-\frac{i{\omega }_0}{l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}{\displaystyle \int }dr{e}_z\mathrm{<mo\; stretchy="false">(</mo>}{u}_0\times u_0^{\mathrm{<mo>*</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\text{ (53)}$

where we put Ω = Ωe_z, e_z being the unit vector along the z -axis. Here we use e_z = cosθe_r − sinθe_θ, $u_0\mathrm{<mo>=</mo>}{h}_l^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}{T}_{lm}$ and T_lm from equations (10); we get immediately

${\omega }_1\mathrm{<mo>=</mo>}{\omega }_0\frac{m}{l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}\text{ (54)}$

where m denotes all integers from -l to l. It follows that the frequencies wln, which are degenerate with respect to m, are split into 2 l + l branches

${\omega }_{ln}\to {\omega }_{ln}+\Omega \frac{m}{l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>;</mo>}\text{ (55)}$

using ω₁ thus determined, we can get u₁ from the second equation (52). Higher-order contributions can be obtained similarly. An m -band occurs for each ω_ln, of the widths 2Ω/(l+1), with the separation frequency Ω/(l+1). For a typical eigenperiod 60 minutes the ratio Ω/ω₀ is approximately 1/20 ≪ 1.

Centrifugal force

The equation of the elastic motion for a body in rotation with a (constant) angular velocity Ω reads

$$\rho \ddot{u}+2\rho \Omega \times \dot{u}+\rho \Omega \times \mathrm{<mo\; stretchy="false">[</mo>}\Omega \times (r+u\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo><mo>=</mo>}\mu \Delta u+\mathrm{<mo\; stretchy="false">(</mo>}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}graddivu+F\mathrm{,}\text{ (56)}$$

where F is an external force. We note that the centrifugal term ρΩ×(Ω×r) is static, so we can write it as

$F_c\mathrm{<mo>=</mo>}\rho \Omega \mathrm{<mo\; stretchy="false">(</mo>}\Omega \text{r}\mathrm{<mo\; stretchy="false">)</mo>}-\rho {\Omega }^{2}r\mathrm{,}\text{ (57)}$

where we denoted by F_c the centrifugal force and removed any other external force (F=0); we may neglect u in the centrifugal force, since it is very small in comparison to r. The angular velocity is oriented along the z -axis, Ω = Ωez. Making use of ez = cos er − sin θe_θ and the spherical harmonics

$Y_{00}\mathrm{<mo>=</mo>}\frac{1}{\sqrt{4\pi }}\mathrm{,}{Y}_{20}\mathrm{<mo>=</mo>}\sqrt{\frac{5}{16\pi }}\mathrm{<mo\; stretchy="false">(</mo>1}-3{{\displaystyle \cos }}^2\theta \mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}\text{ (58)}$

it is easy to see that we can write F_c as a series expansion

$F_c\mathrm{<mo>=</mo>}-\rho {\Omega }^{2}r\left(\alpha R_{00}+2\beta R_{20}-\beta S_{20}\right)\text{ (59)}$

in spheroidal functions, where $\alpha \mathrm{<mo>=</mo>2}\sqrt{4\pi }\mathrm{<mo>/</mo>3}$ and $\beta \mathrm{<mo>=</mo>}\sqrt{16\pi \mathrm{<mo>/</mo>5}}$ . We seek a similar expansion for the displacement u,

$u\mathrm{<mo>=</mo>}{f}_1{R}_{00}+f_2{R}_{20}+g{S}_{20}\mathrm{<mo>;</mo>}\text{ (60)}$

equations (12) lead to

$\begin{array}{c}{f}_1^{{\text{'}}^{\prime }}+\frac{2}{r}{f}^{\mathrm{<mo>\text{'}</mo>}}-\frac{2}{r^2}{f}_1\mathrm{<mo>=</mo>}-\frac{\rho {\Omega }^2\alpha }{\lambda +2\mu }r\mathrm{,}\\ \\ f_2^{{\text{'}}^{\prime }}+\frac{2}{r}{f}_2^{\mathrm{<mo>\text{'}</mo>}}-\frac{\mathrm{2<mo\; stretchy="false">(</mo>}\lambda +5\mu \mathrm{<mo\; stretchy="false">)</mo>}}{\lambda +2\mu }\frac{1}{r^2}{f}_2+\frac{\mathrm{6<mo\; stretchy="false">(</mo>}\lambda +3\mu \mathrm{<mo\; stretchy="false">)</mo>}}{\lambda +2\mu }\frac{1}{r^2}g-\frac{\mathrm{6<mo\; stretchy="false">(</mo>}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}}{\lambda +2\mu }\frac{1}{r}{g}^{\mathrm{<mo>\text{'}</mo>}}\mathrm{<mo>=</mo>}-\frac{2\rho {\Omega }^2\beta }{\lambda +2\mu }r\mathrm{,}\\ \\ g^{{\text{'}}^{\prime }}+\frac{2}{r}{g}^{\mathrm{<mo>\text{'}</mo>}}-\frac{\mathrm{6<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}}{\mu }\frac{1}{r^2}g+\frac{\mathrm{2<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}}{\mu }\frac{1}{r^2}{f}_2+\frac{\lambda +\mu }{\mu }\frac{1}{r}{f}_2^{\mathrm{<mo>\text{'}</mo>}}\mathrm{<mo>=</mo>}-\frac{\rho {\Omega }^2\beta }{\mu }r\mathrm{.}\end{array}\text{ (61)}$

We seek solutions of these equations of the form f_1,2,g = Arⁿ; the solution of the homogeneous equations (regular in the origin) corresponds to n =1; we get

$f_1\mathrm{<mo>=</mo>}-\frac{\rho {\Omega }^2\alpha }{\mathrm{10<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}}{r}^3+C_{1}r\text{ (62)}$

and

$f_2\mathrm{<mo>=</mo>}{C}_{2}r\mathrm{,}g\mathrm{<mo>=</mo>}\frac{\rho {\Omega }^2\beta }{6\lambda }{r}^3+C_{3}r\mathrm{,}\text{ (63)}$

where C_1,2,3 are constants of integration. These constants are determined by the boundary conditions given by equations (17) for a free surface. Finally, we get the displacement

$\begin{array}{c}u\mathrm{<mo>=</mo>}-\frac{\rho {\Omega }^2}{3\lambda }\left[\frac{\lambda }{\mathrm{5<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}}r\left(r^2-\frac{5\lambda +2\mu }{3\lambda +2\mu }{R}^2\right)-R^{2}r\mathrm{<mo\; stretchy="false">(</mo>1}-3{{\displaystyle \cos }}^2\theta \mathrm{<mo\; stretchy="false">)</mo>}\right]e_r+\\ \\ +\frac{\rho {\Omega }^2}{3\lambda }r\left[r^2-\frac{\mathrm{2<mo\; stretchy="false">(</mo>3}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}}{3\lambda }{R}^2\right]\sin \theta \cos \theta e_{\theta }\mathrm{.}\end{array}\text{ (64)}$

It is worth estimating the equatorial displacement (q = p/2) for the Earth radius R = 6370km; with r= 5g/cm³ and l,m = 10¹¹dyn/cm² we get u = u_r, ≃ 10km.

Earthquake " temperature"

Let us multiply by u the equation of the elastic motion,

$\rho \ddot{u}+\mu curlcurlu-\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}graddivu\mathrm{<mo>=</mo>}F\mathrm{<mo>;</mo>}\text{ (65)}$

integrating by parts, we get the law of energy conservation

$\frac{\partial ℰ}{\partial t}\mathrm{<mo>=</mo>}-divS+w\mathrm{,}\text{ (60)}$

where

$ℰ\mathrm{<mo>=</mo>}\frac{1}{2}\rho {\dot{u}}^2+\frac{1}{2}\mu {\mathrm{<mo\; stretchy="false">(</mo>}curlu\mathrm{<mo\; stretchy="false">)</mo>}}^2+\frac{1}{2}\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}divu{\mathrm{<mo\; stretchy="false">)</mo>}}^2\text{ (67)}$

is the energy density,

$S_i\mathrm{<mo>=</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{\dot{u}}_j{\partial }_j{u}_i-{\dot{u}}_j{\partial }_i{u}_j\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}\lambda +2\mu \mathrm{<mo\; stretchy="false">)</mo>}{\dot{u}}_i{\partial }_j{u}_j\text{ (68)}$

are the components of the energy flux density and $w\mathrm{<mo>=</mo>}\dot{u}F$ is the density of mechanical work done by the external force per unit time. It is worth noting that the energy density given by equation (67) differs from the energy density derived from the other form of the equation of motion, e.g.,

$\rho \ddot{u}-\mu \Delta u-\mathrm{<mo\; stretchy="false">(</mo>}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}graddivu\mathrm{<mo>=</mo>}F\mathrm{,}\text{ (69)}$

by the divergence of a vector; it follows that the energy density and the energy flux density are not unique (well defined).

Making use of equations (6), (7) and (10) we can write symbolically

$\begin{array}{c}curlu\mathrm{<mo>=</mo>}\frac{h}{r}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}R+\frac{1}{r}\frac{d}{dr}\mathrm{<mo\; stretchy="false">(</mo>}rh\mathrm{<mo\; stretchy="false">)</mo>}S+\left[\frac{f}{r}-\frac{1}{r}\frac{d}{dr}\mathrm{<mo\; stretchy="false">(</mo>}rg\mathrm{<mo\; stretchy="false">)</mo>}\right]T\mathrm{,}\\ \\ divu\mathrm{<mo>=</mo>}\frac{1}{r^2}\frac{d}{dr}\mathrm{<mo\; stretchy="false">(</mo>}{r}^{2}f\mathrm{<mo\; stretchy="false">)</mo>}-\frac{g}{r}l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}\mathrm{.}\end{array}\text{ (70)}$

We compute the total energy E by introducing these expressions for $curlu$ and divu in equation (67), integrating over the solid angle and integrating by parts over the radius r; for large values of R the boundary conditions given by equations (17) for free vibrations ensure the vanishing of the " surface terms" in the r-integration by parts; in addition, for large values of R we may neglect the f -term in $curlu$ and the g -term in divu; making use of the equations of motion (12), we get finally

$E\simeq \frac{2l+1}{8\pi }{\displaystyle \int }dr\rho {\omega }^2\left[f^2+l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}{g}^2+l\mathrm{<mo\; stretchy="false">(</mo>}l+\mathrm{1<mo\; stretchy="false">)</mo>}{h}^2\right]\mathrm{,}\text{ (71)}$

where the summation over l is omitted (the factor 2l + 1 arises from the summation over m). The functions ωf, ωg and ωh in equation (71) are superpositions of their own normal modes (labeled by n); for large values of R all these eigenmodes may be taken as the spherical Bessel functions, and the eigenfrequencies are given by the zeros of the derivatives of the spherical Bessel functions; we note that these eigenmodes are orthogonal concerning the r -integration; the f -part in equation (71) is related to the velocity c₁ (the combination of $\lambda +2\mu$ of the elastic moduli), while the g - and h -parts are related to the velocity c₂ (modulus µ).

Let us write the energy given by equation (71) for the normal modes as

$E\simeq \frac{1}{8\pi }{\displaystyle \sum _{lmn}}{\displaystyle \int }dr\left[\rho {\omega }_{ln}^{\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>2}}{f}_{ln}^2+\rho {\omega }_{ln}^{\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>2}}{g}_{ln}^2+\rho {\omega }_{ln}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>2}}{h}_{ln}^2\right]\mathrm{,}\text{ (72)}$

where the summation over m is restored and the coefficients l(l+1) are included in g_ln and h_ln. We may use approximately the asymptotic expressions for the functions f_ln, g_ln h_ln of the form f_ln = a_ln cos[kr−(l+1)π / 2]/kr (spherical Bessel functions), with amplitudes a_ln; and, similarly, for g_ln and h_ln with amplitudes b_ln and, respectively, c_ln Effecting the integral, we get

$E\simeq \frac{1}{4}R{\displaystyle \sum _{lmn}}\left[\rho c_1^2{a}_{ln}^2+\rho c_2^2{b}_{ln}^2+\rho c_2^2{c}_{ln}^2\right]\mathrm{,}\text{ (73)}$

where R is the radius of the sphere and c_1,2 are the wave velocities. This is a simple expression, of the form

$E\mathrm{<mo>=</mo>}{\displaystyle \sum _s}\rho R{c}^2{a}_s^2\mathrm{,}\text{ (74)}$

where is a generic notation for the normal modes.

Let us assume that energy E is given to the vibrating sphere; we ask how it is distributed among the normal modes. It is reasonable to assume that, after many reflections from the surface, the distribution of energy reaches an equilibrium state, in the sense that it does not depend anymore on time. This state is characterized by a probability density w, which is multiplicative for different spheres; In w is additive, and the function

$S\mathrm{<mo>=</mo>}-w\ln w\text{ (75)}$

should have a maximum value in the equilibrium state, corresponding to a maximal " disorder"; this represents our idea of equilibrium. Obviously, the function S given by equation (75) is the entropy. Its maximum value for constant energy is reached for the extremum of the function S −βwE, where β is a Lagrange multiplier; we get the Boltzmann (canonical) distribution

$w\mathrm{<mo>=</mo>}const\cdot e^{-\beta E}\mathrm{,}\text{ (76)}$

or, for one mode,

$w\mathrm{<mo>=</mo>}\sqrt{\beta \rho R{c}^2\mathrm{<mo>/</mo>}\pi }{e}^{-\beta \rho R{c}^2{a}^2}\mathrm{.}\text{ (77)}$

The mean energy per mode is

$\overline{e}\mathrm{<mo>=</mo>}\frac{1}{2}T\mathrm{,}\text{ (78)}$

as expected, and the mean value of the square amplitude is

$\overline{a^2}\mathrm{<mo>=</mo>}\frac{T}{2\rho R{c}^2}\mathrm{,}\text{ (79)}$

where we introduced the temperature T = 1/β. The total mean energy is $\overline{E}\mathrm{<mo>=</mo>}N\overline{e}\mathrm{<mo>=</mo>}NT\mathrm{<mo>/</mo>2}$ , where N is the total number of modes; this equality gives the temperature parameter.

Making use of the asymptotic expressions of the spherical Bessel functions (for the radial functions) we get the normal modes given by $k_{ln}R\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>2}n+l+\mathrm{1<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo>/</mo>2}$ ; hence, we see that the normal modes are equidistant; the corresponding wavelengths are ${\lambda }_{ln}\mathrm{<mo>=</mo>4}R\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>2}n+l+\mathrm{1<mo\; stretchy="false">)</mo>}$ . We may take, tentatively, a cutoff of short wavelengths of the order 10⁻⁴ cm (1µm, corresponding to a frequency ≃5GHz, for velocity 5km/s); it is reasonable to admit that below this distance the homogeneous elastic qualities of the Earth do not hold anymore. For this cutoff, we get a maximum number 2n+l+1 of the order N_c = 10¹³ and several modes of the order $N\mathrm{<mo>=</mo>}{N}_c^3{\mathrm{<mo>=</mo>10}}^{39}$ . Also, it is reasonable to assume that the earthquake energy is spent mainly on mechanical work (like fracture of the rocks, structural damage, etc). The largest part of the energy released in an earthquake is spent in mechanical work associated with the motion of the rocks, soil, and the damage produced at the Earth's surface; the remaining is dissipated as heat, after a long while. We may assume tentatively that only the energy $\overline{E}\simeq {10}^{15.65}erg$ is spent in thermalization (corresponding to magnitude M_wM = 0) in the Gutenberg-Richter law $\lg E\mathrm{<mo>=</mo>}\frac{3}{2}{M}_w+15.65$ (in erg). This way, we get a temperature T = 10⁻²³ erg (i.e., ≃10⁻⁷ K, since $1.38\times {10}^{-16}erg\mathrm{<mo>=</mo>1}K,1erg=1dyn\cdot cm$ ); the inner Earth's temperature is ≃6000K. The quantity ρRc² in equation (78) is ρRc² ≃ 10²⁰ g/s² (for ρ = 5g/cm³, R≃6×10⁸ cm and c=5km/s). The estimation of the temperature is very sensitive to the number of eigenmodes N, which may be much lower; also, we may allow for the thermalization of higher energies, corresponding to higher magnitudes. In both cases the temperature increases. We note that the cutoff wavelength, which affects essentially the numerical estimation of the temperature, corresponds to the mean inter-atomic distance in the Debye estimation of the statistical equilibrium of the elastic vibrations (phonons) in crystals.

Concluding remarks

A systematic analysis of vibrations of an elastic sphere is presented in this paper, to get analytical results. Such results are useful, on one hand, in analyzing the recorded data, and, on the other, in comparing them with the current numerical investigations. It is shown that the hypothesis of a large radius, appropriate for Earth's vibrations, simplifies the analysis to a great extent. We discuss the general formulation of the problem, the use of spherical harmonics, the approximation of a large radius, and the example of a fluid sphere. Specific results are given for toroidal and spheroidal vibrations. Also, the self-gravitation and rotation effects are analyzed in detail.

Apart from self-gravitation and rotation, the inhomogeneities may have an important effect on the vibrations of the solid sphere. For instance, from equation (1), a (uniform) change δρ in density causes a change δω /ω = −δρ/2ρ in frequency. The effect of similar changes in the elastic moduli λ and µ can be estimated by using the changes in the wave velocities c in the relation ω_ln ≃ (c/R)β_ln.

An approximate procedure is given in this paper for estimating the spectrum of eigenfrequencies (and eigenfunctions) of the vibrations of a solid sphere, with application to Earth's vibrations, as those produced by an earthquake. The procedure is sufficiently convenient to apply to other, more complex situations involving the vibrations of a solid sphere, as, for instance, the corrections brought about by self-gravitation, rotation, and inhomogeneities. The distribution of the energy among the vibrations eigenmodes is also estimated here and the concept of the earthquake " temperature" is tentatively introduced, as another means of characterizing earthquakes and estimating the earthquake's effects.

The author is indebted to his colleagues in the Department of Engineering Seismology, Institute of Earth’s Physics, Magurele-Bucharest, for many enlightening discussions, and to the members of the Laboratory of Theoretical Physics at Magurele-Bucharest for many useful discussions and a throughout checking of this work. This work was carried out within the Program Nucleu SOL4RISC, contract #24N/03.01.2023, funded by the Romanian Ministry of Research, Innovation and Digitalization, project # PN23360202.

1. Aki K, Richards PG. Quantitative Seismology. University Science Books, Sausalito, CA. 2009.
2. Ben-Menahem A, Singh JD. Seismic Waves and Sources. Springer, New York. 1981.
3. Lamb H. On the vibrations of an elastic sphere. Proc. London Math. Soc. 1882; 13:189-212.
4. Lamb H. On the oscillations of a viscous spheroid. Proc. London Math. Soc. 1881; 13:51-66.
5. Hansen WW. A new type of expansion in radiation problems. Phys. Rev. 1935; 47:139-143.
6. Bromwich TJI'A. On the influence of gravity on elastic waves, and, in particular, on the vibrations of an elastic globe. Proc. London Math. Soc. 1898; 30:98-120.
7. Landau L, Lifshitz E. Course of Theoretical Physics, Theory of Elasticity, Elsevier, Oxford. 1986; 7.
8. Apostol BF. Elastic waves inside and on the surface of an elastic half-space. Q. J. Mech. Appl. Math. 2017; 70:281-308.
9. Apostol BF. Elastic displacement in a half-space under the action of a tensor force. A general solution for the half-space with point forces. J. Elast. 2017; 126:231-244.
10. Barrera RG, Estevez GA, Giraldo J. Vector spherical harmonics and their applications to magnetostatics. Eur. J. Phys. 1985; 6:287-294.
11. Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. National Bureau of Standards, USA Government Printing Office, Washington. 1964.
12. Mendiguren J. Identification of free oscillation spectral peaks for 1970 July 31 Colombian deep shock using the excitation criterion. Geophys. J. Roy. Astron. Soc. 1973; 33:281-321.
13. Mendiguren J. High-resolution spectroscopy of the Earth's free oscillations, knowing the earthquake source mechanism. Science. 1973; 179:179-180.
14. Gilbert F, Dziewonski AM. An application of the normal mode theory to the retrieval of the structural parameters and source mechanisms from seismic spectra. Phil. Trans. Roy. Soc. (London). 1975; A278:187-269.
15. Haskell NA. The dispersion of surface waves in multilayered media, Bul. Seism. Soc. Am. 1953; 43:17-34.
16. Ewing M, Jardetzky W, Press F. Elastic Waves in Layered Media. McGraw-Hill, NY. 1957.
17. Benioff H, Gutenberg B, Richter CF. Progr. Report. Trans. Am. Geophys. Union. 1954; 35:979-987.
18. Benioff H, Press F, Smith SW. Excitation of the free oscillations of the Earth by earthquakes. J. Geophys. Res. 1961; 66:605-620.
19. Ness NF, Harrison JC, Slichter LB. Observations of the free oscillations of the Earth. J. Geophys. Res. 1961; 66:621-629.
20. Alsop LE, Sutton GH, Ewing M. Free oscillations of the Earth observed on strain and pendulum seismographs. J. Geophys. Res. 1961; 66:631-641.
21. Kanamori H. The energy release in great earthquakes. J. Geophys. Res. 1977; 82:2981-2987.
22. Block B, Dratler J, Moore RD. Earth and normal modes from a  magnitude earthquake. Nature. 1970; 226:343-344.
23. Lehmann I. Publications du Bureau Central Seismologique International, Serie A, Travaux Scientifiques. 1936; 14:87-115.
24. Dziewonski AM, Gilbert F. Solidity of the inner core of the Earth inferred from normal mode observations. Nature. 1971; 234:465-466.
25. Pekeris CL, Accad Y. Dynamics of the liquid core of the Earth. Phil. Trans. Roy. Soc. London.1972; A273:237-260.
26. Dahlen FA, Smith ML. The influence of rotation on the free oscillations of the Earth. Phil. Tras. Roy. Soc. London. 1975; A279:583-629.