Resonant and non-resonant flows in longitudinally and latitudinally librating spheres

Many planetary satellites undergo a phenomenon called libration. This term refers to a harmonic oscillation of the boundary of a rotating body with respect to its mean spin axis. One can distinguish between two types of librations: east–west and north–south oscillations are referred to as librations in longitude and latitude, respectively (see Fig. 1). Another distinction can be made based on the physical mechanism giving rise to the librations. First, forced librations are caused by time-dependent gravitational torques that are exerted on an orbiting planet or moon. The forced longitudinal librations arise from the tidal bulge resulting from a low spin–orbit resonance and an eccentric orbit leading to a time-dependent orbital velocity, i.e., Kepler's third law.^1,2 The periodic misalignment of the bulge generates an axial torque, leading to a harmonic oscillation of the rotation rate of the planets and moons around their mean value, as illustrated in Fig. 2(a). Forced latitudinal librations are due to the obliquity of the orbit combined with the polar flattening of the mantle, leading to an equatorial torque that tends to realign the equator of the planets or moons with their orbital planes. As a result, the celestial objects exhibit a small oscillation along an equatorial axis that is fixed in the frame rotating with the body, see Fig. 2(b). Second, free librations both in longitude and latitude may follow from a large impact; in such a case the oscillations are expected to decay over time due to dissipation both in the solid and liquid shells. Note also that nutations and polar motion can be thought of as types of latitudinal librations, and variations in length of day can be thought of as a type of longitudinal libration, although in all these cases the frequencies involved are much smaller than for the librations we are interested in.

Starting with the pioneering work of Poincaré,⁴ many scientists have used observations of changes in the spin rate of a planet or moon to infer its interior structure. In most of these studies, the couplings between fluid and solid parts of a planet are often simplified under the ad hoc assumption that the fluid is essentially behaving as a rigidly rotating solid body. Testing the validity of this approach requires detailed fluid mechanical investigation of the response of a fluid body to a change in the rotation of its overlaying mantle. Roughly speaking, one can identify three physical mechanisms that can transmit the motion of the mantle to the fluid: (i) topographic coupling due to a deviation of sphericity, (ii) viscous coupling emanating from the core–mantle boundary, and (iii) electromagnetic coupling in the presence of a magnetic field and if the mantle is electrically conducting (at least weakly).

The first fluid mechanical investigations of libration-driven flows in spherical enclosures were carried out by Greenspan⁵ and Aldridge and Toomre,⁶ who used a theoretical and experimental approach, respectively, to study the linear response of the fluid to longitudinal libration forcing (see also Ref. 7). Later studies have addressed the genesis of boundary-layer instabilities, the presence of an inner shell, and the zonal flows produced by nonlinear effects.^8–16 As far as topographic coupling is concerned, most studies have focused on the instabilities and turbulence that ensue from an elliptical or higher-order multipolar deformation of the boundary.^17–20 For a general discussion, we refer also to the two reviews on precession, libration and tidally driven flows, Le Bars et al.²¹ and Chap. 4 of Le Bars and Lecoanet.²² 

Much less attention has been devoted to flows driven by libration in latitude. References 23 and 24 investigated the spheroidal geometry and found that a topographic resonance can occur when the libration frequency matches the eigenfrequency of the spin-over inertial mode. This was extended toward triaxial ellipsoids by Ref. 25, who also investigated parametric instabilities driven by the shear and strain induced by the ellipsoidal geometry.

In this work, we focus on viscous coupling mechanisms. Using theoretical and numerical methods, we study both longitudinal and latitudinal libration-driven flows in the limit of vanishing viscosity and small forcing amplitude. First, we introduce the governing equations in Sec. II. In Sec. III, we use the classical boundary layer theory of Greenspan to compute the linear response of the fluid to the libration forcing; in Secs. IV and V, we compare our theoretical results against high-resolution numerical simulations.

Our results and Zhang and Liao are in quantitative agreement, as expected from the similarity of the derivation but for the choice of frame of reference. Zhang and Liao²⁷ also treated latitudinal librations only in spheroids rather than spheres, so no direct comparison with our results is available in that case. We find that the amplitude of the kinetic energy at the peak is effectively independent of the Ekman number, and in particular, there is no divergence of the amplitude of the modes as the Ekman number vanishes, in agreement with Refs. 7 and 27.

To go beyond the first order approximation, we account for the viscous correction to the inviscid eigen-frequency S_lmk derived by Greenspan²⁶ such that $ω lmk ( E ) = ω lmk + S lmk E 1 / 2$⁠, also reported in Tables I and II. Note that Zhang and Liao²⁷ did not consider the viscous correction to the eigen-frequency, whereas Greenspan²⁶ considered both the real and imaginary parts of the viscous correction to the eigen-frequency. Nevertheless, the viscous correction of $O ( E 1 / 2 )$ to the inviscid eigen-frequency only produces small discrepancies in the total energy as we can see from Table I.

The boundary conditions associated with Eq. (32) are Eqs. (10) or (11) as before, except that the previous amplitude factors $Δ ε ϕ$ and $Δ ε θ$ are now absent since the problem is linear. Working through the details, and recalling also that the $cos ( ω t )$ factors there are already accounted for by the $e i ω t$ factors here, the final result is remarkably simple, namely, just $e 1 ( 1 ) = 1$⁠, for both longitudinal and latitudinal libration. For all other modes, the original no-slip boundary conditions on $u$ convert to simply $e n = f n = f ′ n = 0$⁠.

The original equation and its associated boundary conditions therefore convert to a matrix problem $A v = b$⁠, where $A$ is this block tri-diagonal matrix discretising the left-hand side of Eq. (32) and the boundary conditions, $v$ is a column vector containing the spectral coefficients e_nk and f_nk, and $b$ is a column vector containing all zeros, except for a single non-zero entry corresponding to the librational forcing term $e 1 ( 1 ) = 1$⁠. The NAG routines F07BRF and F07BSF (https://www.nag.com/numeric/nl/nagdoc_24/nagdoc_fl24/html/f07/f07brf.html and https://www.nag.com/numeric/nl/nagdoc_24/nagdoc_fl24/html/f07/f07bsf.html) are then used to first compute the LU decomposition of $A$⁠, and then do the back-substitution to obtain the desired solution $v$⁠. In principle one could then do all numerical calculations with this code as just outlined. There are though two further improvements that were implemented.

First, the relationship r = x was altered to $r = x + δ ( x − x 3 )$⁠. For $0 < δ < 0.5$⁠, this has the property of concentrating the Chebyshev collocation points even more toward the boundary at r = 1 than they naturally already are. It was found empirically that taking $δ = 0.45$ allows the radial resolution K to be up to 30% less than the original δ = 0 version, while still yielding identical results. The memory requirements to store the matrix $A$ scale as NK², so reducing K by 30% allows a substantial increase in the angular resolution N, for a given core memory size.

The final procedure is thus as follows: first compute the LU decomposition of $A 0$⁠, corresponding to some frequency ω₀. This is accomplished using NAG routine F07BRF, and becomes computationally very expensive at high resolutions, scaling as NK³. Next iterate according to Eq. (38), and for some $± Δ ω$ range centered around ω₀. Each iteration is accomplished using NAG routine F07BSF and is computationally much cheaper, scaling only as NK². How large $| Δ ω |$ can be taken before the iteration fails to converge, and a new LU decomposition corresponding to a new ω₀ has to be computed, is determined by simply trying it out. It was found that each of the peaks and troughs seen in the figures to follow typically required a few ω₀ values, but once the corresponding LU decompositions had been computed, this iterative solution method then allowed the precise details of the curves to be mapped out in far greater detail at relatively little additional cost. The final result is then identical to what the original direct solver would have accomplished, but allows a very fine scan of ω for $∼ 1 / 30$ the CPU time that a separate LU decomposition at each new frequency would have required. Specifically, at the highest resolutions, the factor of K saved by an iteration (NK²) vs a new LU decomposition (NK³) means that the iterative procedure is far faster, even if on average a dozen or so iterations might be required for each value of $Δ ω$⁠.

With these improvements in place, the code was then run at resolutions $N × K = 80 × 170$ for Ekman number $E ≥ 10 − 8$⁠, 100 × 270 for $E = 10 − 9$⁠, and 130 × 430 for $E = 10 − 10$⁠, where according to Eqs. (34) and (35), the maximum spherical harmonic degree is then 2N in each case. Recall also that the memory requirements scale as NK² and the CPU time for each LU decomposition scales as NK³. As E is reduced, this iterative procedure (38) also converges for increasingly small $| Δ ω |$ ranges as the peak/trough structures become increasingly narrow, thus requiring more and more LU decompositions to map out a given ω interval. The overall computational effort, therefore, increases even more rapidly than NK³. The final result is that $E ≥ 10 − 8$ is quite straightforward on a standard PC, $E = 10 − 9$ still fits on a PC in terms of memory, but requires considerably more CPU time, whereas $E = 10 − 10$ requires a machine with larger memory, and becomes so CPU-intensive that only a small fraction of the full range $ω ∈ [ 0 , 2 ]$ could be scanned.

Figures 3 and 4 show the time-averaged kinetic energy $E num / ε ϕ , θ 2$ as a function of the libration frequency ω. The upper panels in each figure show scans over the entire frequency range $ω ∈ [ 0 , 2 ]$⁠, computed with a spacing $Δ ω = 2 × 10 − 3$⁠, and down to Ekman numbers $E = 10 − 8$⁠. The lower panels zoom in on subsets covering only 1/10 of the full frequency range and include even smaller Ekman numbers on those subsets. In both figures, we observe a complicated series of peaks and troughs and exhibiting finer and finer structures for smaller E. As E is reduced the values at the peaks become independent of E, in agreement with the previous theoretical analysis; the main peaks can also be associated with some of the inertial modes given in Tables I and II.

In contrast, at the troughs between the peaks, the energies scale as $E 1 / 2$⁠. Some of the most prominent troughs are also seen to correspond to the frequencies $ω = 2 sin ( π / n )$ for integer n, for which the internal shear layers form particularly simple reflection patterns (periodic orbits). The upper panels in Figs. 3 and 4 indicate these frequencies for $n = 3 , … , 10$⁠.

For longitudinal libration, all of these values correspond to troughs, but for latitudinal libration, some of them are actually peaks. For example, ω = 1 is the n = 6 reflection pattern, but according to Table II is also one of the resonant modes (the spin-over mode in this case). In cases such as this the resonant mode behavior dominates, and the result is a peak rather than a trough.

The other striking feature about Figs. 3 and 4 is how more and more peaks emerge as the Ekman numbers are reduced. Indeed, the pattern appears to be fractal: if we compare $E = 10 − 7$ in the top rows with $E = 10 − 9$ in the bottom rows, and temporarily ignoring different ω ranges, the curves are qualitatively very similar. If we then recall that the frequency ranges in the bottom rows are in each case 1/10 of the full ranges in the top rows, this suggests $E 1 / 2$ as the self-similarity factor. Figure 5 further quantifies this by showing the distributions of the distance between two consecutive peaks, $( ω n + 1 − ω n )$⁠, re-scaled by $E 1 / 2$⁠. The fact that distributions for different Ekman numbers collapse on top of one another confirms that $E 1 / 2$ is the correct self-similarity factor. In the limit $E → 0$⁠, the entire interval $ω ∈ [ 0 , 2 ]$ would then become densely filled with peaks but always also still containing troughs between the peaks.

Returning briefly to the peaks, we saw previously that the numerically computed energies are in excellent agreement with some of the theoretical values for the inertial modes in Tables I and II. Figure 6 shows an example illustrating that the corresponding spatial structures are also in near-perfect agreement with the theoretical mode profiles given by expressions (A5)–(A7)). The top row shows the $( m , l , k ) = ( 0 , 4 , 1 )$ free eigenmode (which by itself has an arbitrary amplitude) matched up with the forced librational response at its frequency $ω mlk = 12 / 7$ in the bottom row. If the amplitude of the free mode is suitably adjusted, the two are seen to be virtually identical, confirming that the libration is indeed picking out this particular free mode.

Regarding the troughs, these are (relatively large) frequency bands with no inertial mode resonances at a given Ekman number. In these regions, the velocity field is dominated by conical shear layers emanating from the critical latitudes making an angle $θ = arcsin ( ω / 2 )$ with the vertical axis. As an example, we show results for longitudinal libration at $ω = 2 sin ( π / 8 ) = 0.765$ in Fig. 7. As explained in the study by Rieutord et al.,³⁴ for frequencies of the form $ω = 2 sin ( π q )$⁠, with q a rational number, the ray patterns of the inertial waves form periodic orbits, which has the consequence that regular inviscid inertial modes, in general, do not exist at these frequencies. Therefore, no inertial mode would be excited at these particular frequencies in the limit of vanishing viscosity, hence explaining the troughs. The detailed scalings within these conical shear layers were also predicted by Stewartson and Roberts,³⁵ with the thickness scaling as $E 1 / 5$ and their characteristic velocities as $E 1 / 5$⁠. By examining cuts along the blue line in Fig. 7, we obtained results in excellent agreement with these scalings.

Note finally that all the features obtained in this section confirm the theoretical inviscid spectrum predicted by Fig. 2 of Aldridge and Toomre.⁶ 

Although numerical simulations such as the ones performed here allow us to compute Ekman numbers as small as $10 − 10$⁠, the values in planetary cores can be several orders of magnitude smaller still, as small as $O ( 10 − 15 )$⁠. In this section, we therefore quantify the convergence of the asymptotic model as we decrease the Ekman number.

In Fig. 8, we show the frequency-shift of the peak corresponding to the resonance with the modes (0,8,3) and (1,2,1) excited by librations in longitude and latitude, respectively. The numerical results (full-circles) converge toward the asymptotic values $ω lmk + S lmk E 1 / 2$ proposed by Greenspan.²⁶ For all peaks in Tables I and II, we next compute the residual errors between the numerical values ω_num and the asymptotic values, that is, we evaluate the quantity $Δ ω lmk = ω num − ω lmk − S lmk E 1 / 2$⁠. As shown in Fig. 9, these errors are bounded by $Δ ω lmk < O ( E )$⁠, as expected from the theory proposed by Greenspan.²⁶ 

The second quantity that an asymptotic model will allow us to estimate in the case of a resonance at planetary settings is the kinetic energy, from which one can derive a typical rms velocity, $U = E kin$⁠. From a theoretical standpoint we expect the amplitude of the peaks in the energy spectrum to be Ekman independent in the asymptotic regime, except for small contributions from the boundary layer at $O ( E 1 / 2 )$ and the interior conical shear layers at $O ( E 3 / 5 )$⁠. In Fig. 10, we compare the energy at the peak frequencies from our numerical simulations with the asymptotic values in Tables I and II for libration in longitude and latitude, respectively. As before for the frequency, we observe a good quantitative agreement between the numerical simulations and the asymptotic model, even at moderate values of the Ekman number for the large-scale modes, having small l, k values.

We now focus on the trough regions introduced in Sec. V A, where no resonant modes are excited, and the low-energy flow instead resembles the periodic orbit corresponding to the particular frequencies $ω = 2 sin ( π q )$⁠, with q a rational number. The trough regions cannot be captured by the asymptotic model as they correspond to flows driven by the boundary layer critical latitude, thus requiring a small but finite Ekman number. As we showed in Fig. 5, the distribution of the distance between peaks is at first order self-similar when re-scaled by $E 1 / 2$⁠, from which one might conjecture that the width of the troughs would also follow the same scaling. We shall see this is not the case though.

We define a trough in our numerical spectrum as the region around one of the frequencies associated with a periodic orbit bounded by the two nearest peaks, as shown in Fig. 11 for $ω = 2$⁠. The nearest peaks at $ω −$ and $ω +$ define the width of the trough as $δ = | ω + − ω − |$⁠. The characteristic energy of a trough is defined at the central frequency $f = 1 2 ( ω + + ω − )$⁠, indicated by the red circles in the left panel of Fig. 11. We perform the same analysis for all the troughs near the frequencies of the periodic orbits represented by vertical lines in Figs. 3 and 4.

For each trough, we fit a power law using linear regression in log–log space for the width and energy, $δ ∝ E α$ and $E num / ε ϕ , θ 2 ∝ E α ′$⁠, respectively. Figure 12 shows the distribution of the exponents for the width and energy scaling laws for all periodic orbit troughs driven by longitudinal and latitudinal librations. Quite clearly, the energy scales as $E 1 / 2$⁠, consistent with an energy dominated by the flow in the Ekman boundary layer $E num / ε ϕ , θ 2 ∝ E 1 / 2$ with an additional smaller contribution from the conical shear layers in the bulk contributing at the order $E 3 / 5$⁠.

What is more surprising is the scaling of the width $δ ∝ E 0.23 ± 0.05$⁠, close to $E 1 / 4$⁠. Indeed, from the self-similar distribution of the distance between peaks re-scaled by $E 1 / 2$ (see Sec. V A), one would have guessed that the width of the troughs should also follow the same scaling. However, one should note that the widest troughs correspond to the tails of the distributions in Fig. 5, which do not exhibit any self-similarity when re-scaled by $E 1 / 2$⁠. To the best of our knowledge there is no theory that explains these observations regarding the widths of the trough regions.

In Fig. 13, we collapse the longitudinal libration energy spectra in the vicinity of the trough at $ω t = 2$⁠. That is, we rescale the energy by $E 1 / 2$⁠, and the frequency as $( ω − ω t ) / E 0.23$⁠. The troughs for the three different values $E = 10 − 8 , 10 − 9 , 10 − 10$ are seen to almost completely coincide under these rescalings. Similar collapsing is observed for all the troughs, whether driven by libration in longitude or in latitude.

We have used asymptotics and numerics to systematically study libration in both longitude and latitude in a rotating, fluid-filled sphere. The asymptotic theory allows a clear prediction of the resonant eigenmodes, in agreement with previous studies for longitudinal libration by Zhang and Liao,²⁷ which form a dense set on the frequency interval $ω ∈ [ 0 , 2 ]$⁠. Including the first-order viscous correction to the eigenfrequencies, $ω = ω lmk + S lmk E 1 / 2$ yields results that are then always within O(E) of the corresponding numerical results. This indicates that extrapolation to Ekman numbers found in planetary settings, with Ekman numbers as small as $E = O ( 10 − 15 )$⁠, is possible with very high accuracy.

For increasingly small Ekman numbers, the numerically computed frequency response tends toward a fractal structure, with increasingly finely spaced peaks and troughs. The energies at the peaks become independent of E, a well-known results in spheres;²⁷ they correspond to resonant conditions of eigenmodes. The energies at the troughs scale as $E 1 / 2$⁠, and correspond to frequencies where the internal reflection patterns form periodic orbits. The fractal self-similarity factor for the spacing between the peaks is $E 1 / 2$ in general, but the self-similarity factor for the widths of the widest troughs is $E ∼ 0.23 ∼ E 1 / 4$⁠. An asymptotic theory that explains these aspects of the fractal pattern, and especially why peaks and troughs have different self-similarity factors, does not yet exist, but would certainly be of considerable interest.

Librations arise mainly in tidally locked satellites such as the Earth's moon, Callisto, Ganymede, Io and Titan, for which the main librations are at ω = 1. For those satellites with equatorial [ $ξ e q = ( a 2 − b 2 ) / a 2$] and polar [ $ξ pol = ( a 2 − c 2 ) / a 2$] ellipticity small with respect to $E 1 / 4$ (Refs. 23 and 36), our spherical theory provides a realistic estimate of the velocities. Note that here we adopt the convention of Zhang and Liao;²⁷ other conventions exist in the literature leading to different inequalities, for instance adopting the conventions of Grannan et al.¹⁹ [ $β = | a 2 − b 2 | / ( a 2 + b 2 )$] or the polar flattening $η = ( a − b ) / a$ (Refs. 37 and 38) the condition becomes $β , η ≪ E 1 / 2$⁠. For finite ellipticities, i.e., in true ellipsoids, longitudinal librations still do not lead to direct resonances of any inertial modes, resulting in amplitudes scaling comparable to what is obtained in the full sphere, see Ref. 27 for more details. In contrast, for latitudinal librations, a finite polar flattening will result in a pressure coupling leading to a scaling of the spin-over mode amplitude proportional to $Δ θ / E 1 / 2$ (Refs. 25 and 27) near its resonant frequency, i.e., $| ω − 2 / ( 2 − ξ pol 2 ) | ≪ O ( E 1 / 2 )$⁠, potentially very large at planetary settings. Meanwhile, for all other frequencies, it has been shown that latitudinal librations cannot resonate with any inertial mode,²³ similarly to what is found for precession.³⁹ 

A detailed application of our results to any particular celestial object and specific libration frequency would require an extensive discussion of the validity of the assumptions underlying our model, which is beyond the scope of this study. However, one can estimate an order of magnitude of the velocity at the peaks in typical planetary settings. To do so, we shall use a core radius R = 500 km, a rotation period of 1 day and a libration amplitude $ε ∼ 10 − 4$⁠. A typical velocity is obtained from our Tables I and II as $U rms = ε Ω 0 R ( E lmk / ε 2 ) / ( 4 π / 3 )$ leading to motions in the sub-mm/s range.

Including an inner core remains analytically challenging as the inertial modes in a spherical shell form an ill-posed mathematical problem.³⁴ Our approach could be extended to include an inner shell and provide an efficient way to scrutinize the linear response induced by the librations in a spherical shell at low Ekman numbers.⁴⁰ 

Y.L. was supported by NSFC grants (Nos. 42174215 and 12250012) and the National Key R&D Program of China (No. 2022YFF0503200). R.H. was supported by STFC Grant Nos. ST/S000275/1 and ST/W000873/1. R.H. also thanks the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme “Frontiers in dynamo theory: from the Earth to the stars” where work on this paper was undertaken. This work was supported by EPSRC Grant No. EP/R014604/1. R.H.'s visit to the Newton Institute was supported by a grant from the Heilbronn Institute. J.N. was supported by SNSF Grant Nos. 200021_140708 and 200021_185088. J.N. and S.V. are grateful for the provision of computing resources by the Swiss National Supercomputer Centre (CSCS) under allocation s1111. J.N. would like to dedicate this paper to Dr. Hugo Perfettini for 30 years of loyal friendship.

The authors have no conflicts to disclose.

Yufeng Lin: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Rainer Hollerbach: Conceptualization (equal); Formal analysis (lead); Writing – original draft (equal); Writing – review & editing (equal). Jerome Noir: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Stijn Vantieghem: Conceptualization (equal); Formal analysis (lead); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.