We investigate the linear response to longitudinal and latitudinal libration of a rapidly rotating fluid-filled sphere. Asymptotic methods are used to explore the structure of resonant modes in both cases, provided that the nondimensional libration frequency is in the range $ \omega \u2208 [ 0 , 2 ]$. High-resolution numerics are then used to map out this entire frequency range, picking out both the resonant peaks as well as the non-resonant troughs in between. The kinetic energy is independent of the Ekman number *E* at the peaks and scales as $ E 1 / 2$ at the troughs. As the Ekman number is reduced, down to $ E = 10 \u2212 10$ for longitudinal libration and $ E = 10 \u2212 9$ for latitudinal libration, the frequency response also exhibits an increasingly fractal structure, with more and more peaks and troughs emerging. The spacing between peaks is seen to follow an $ E 1 / 2$ self-similarity factor. However, detailed examinations of some of the more prominent troughs shows that their widths follow an $ E \u223c 0.23$ self-similarity factor.

## I. INTRODUCTION

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é,^{4} 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^{5} and Aldridge and Toomre,^{6} 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.*^{21} and Chap. 4 of Le Bars and Lecoanet.^{22}

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.

## II. GOVERNING EQUATIONS

*R*that is undergoing libration. The mean rotation axis of the cavity points in the

*z*-direction, and the mean rotation rate is $ \Omega 0$. We will adopt the convention that for latitudinal libration, the libration axis is parallel to the

*x*-direction. In a frame of reference that is attached to the librating sphere, the so-called mantle frame, the Navier–Stokes equation becomes

*t*,

*p*,

*ρ*, and

*ν*denote the velocity field, time, pressure, mass density, and kinematic viscosity, respectively. The last term on the right-hand side of Eq. (1) is termed the Poincaré acceleration and is a consequence of the non-inertial character of the librating mantle frame. The total rotation vector $\Omega $ is given by

*θ*, the instantaneous inclination with respect to the mean rotation axis. More precisely, we have

_{L}*ω*is the libration frequency.

*R*as a length scale, we find

*E*,

## III. ASYMPTOTIC LINEAR THEORY

*E*justifies the use of the asymptotic linear theory as developed by Greenspan

^{26}and Zhang and Liao.

^{27}In the rest of the paper, we follow a similar approach to these authors recalling the important steps. We perform asymptotic analysis in the reference frame of mean rotation following Greenspan,

^{26}while Zhang and Liao

^{27}used the mantle frame for the asymptotic theory. We, thus, focus on the linearized equations of motion, expressed with respect to the frame of mean rotation,

*p*

_{0}are the leading-order solutions within the interior of the sphere, whereas $ u \u0303$ and $ p \u0303$ denote boundary layer quantities that are asymptotically small everywhere except in the Ekman layer of thickness $ O ( E 1 / 2 )$ near the boundary. As it will turn out, this Ekman layer drives a small mass flux into the interior region, the so-called Ekman pumping, with related quantities $ u \u0302$ and $ p \u0302$.

*E*and can be neglected at leading order. The leading order interior solution is, thus, governed by the system

^{28,29}Based on this result, we decompose the velocity $ u 0$ and the pressure

*p*

_{0}in terms of inertial modes as follows:

*p*are solutions of

_{lmk}*l*and

*m*refer to spherical harmonic degree and order;

*k*corresponds to the radial eigenfunctions. Further properties of these analytically obtainable inertial modes in spherical enclosures are discussed in Appendix A.

*A*vanishes at the leading order of approximation, except if the libration frequency

_{lmk}*ω*equals the eigenfrequency

*ω*. Inertial mode eigenfrequencies being restricted to the interval $ \omega lmk \u2208 [ 0 , 2 ]$, and being dense within this interval, we can conclude the following:

_{lmk}-
If the libration frequency $ \omega \u2265 2$,

*A*is necessarily 0, which implies that it is not possible to drive an_{lmk}*O*(1) flow. -
If, on the other hand, $ \omega < 2$, it may be possible to excite an inertial mode whose eigenfrequency matches the libration frequency. The modal amplitudes

*A*, however, cannot be found directly from Eq. (19), but, as we shall see in the next paragraphs, require taking into account viscous effects, reflecting the purely viscous nature of the coupling between the fluid and the shell._{lmk}

^{26}that we discuss in detail in Appendix B. The final result is

*A*.

_{lmk}*θ*only and are given by (24) and (25). Given a triplet of indices {

*l*,

*m*,

*k*}, Eq. (30) permits the computation of the modal amplitude

*A*given a driving frequency $ \omega = \omega lmk$. Since longitudinal and latitudinal libration are characterized, respectively, by an

_{lmk}*m*= 0 and

*m*= 1 azimuthal dependency, it follows from (30) that these mechanisms can only drive modes with azimuthal wavenumbers

*m*= 0 and

*m*= 1. Taking into account the expression (A4) for the pressure eigenmode at the spherical surface, the

*θ*component of $ \u2207 p lmk$ is of parity $ ( l \u2212 m \u2212 1 )$ with respect to the equator, and the $\varphi $ component is of parity $ ( l \u2212 m )$. Combining this with the specific symmetry of the boundary forcing, we come to the conclusion that the numerator of (30) vanishes if

*l*is odd. In other words, both longitudinal and latitudinal libration can only drive modes for which

*l*is even.

*l*,

*m*,

*k*):

*l*,

*m*,

*k*), and by means of high-resolution numerical simulations. For comparison, we also present the values of the kinetic energy from Zhang and Liao

^{27}(p. 243); noticing that their kinetic energy is volume-averaged, their results must also be multiplied by $ 4 \pi / 3$, and one should only account for the mode contribution and not the apparent fluid motion due to their choice to look at the dynamics from the mantle frame of reference. For instance, the kinetic energy for the mode at $ \omega = 1.3093$, labeled (0, 1, 2) in Zhang and Liao

^{27}and (0, 4, 1) with our conventions, the kinetic energy in their manuscript is

(m, l, k)
. | ω
. _{lmk} | S
. _{lmk} | $ E lmk / \epsilon \varphi 2$ . | $ E num / \epsilon \varphi 2$ . | $ E Z L / \epsilon \varphi 2$ . |
---|---|---|---|---|---|

(0,4,1) | 1.309 307 341 | 0.426 775 579 0 | $ 4.437 \u2009 333 \u2009 386 \xd7 10 \u2212 2$ | $ 4.421 \xd7 10 \u2212 2$ | $ 4.440 \xd7 10 \u2212 2$ |

(0,6,1) | 0.937 697 587 | 0.318 769 489 4 | $ 2.852 \u2009 976 \u2009 237 \xd7 10 \u2212 2$ | $ 2.834 \xd7 10 \u2212 2$ | $ 2.848 \xd7 10 \u2212 2$ |

(0,6,2) | 1.660 447 793 | 0.555 430 038 4 | $ 8.244 \u2009 384 \u2009 738 \xd7 10 \u2212 3$ | $ 8.009 \xd7 10 \u2212 3$ | $ 7.959 \xd7 10 \u2212 3$ |

(0,8,1) | 1.799 515 990 | 0.610 750 963 0 | $ 2.289 \u2009 522 \u2009 846 \xd7 10 \u2212 3$ | $ 2.252 \xd7 10 \u2212 3$ | ⋯ |

(0,8,2) | 0.726 234 926 | 0.253 129 413 2 | $ 1.827 \u2009 816 \u2009 203 \xd7 10 \u2212 2$ | $ 1.834 \xd7 10 \u2212 2$ | ⋯ |

(0,8,3) | 1.354 372 557 | 0.471 011 105 6 | $ 9.463 \u2009 982 \u2009 745 \xd7 10 \u2212 3$ | $ 9.649 \xd7 10 \u2212 3$ | ⋯ |

(m, l, k)
. | ω
. _{lmk} | S
. _{lmk} | $ E lmk / \epsilon \varphi 2$ . | $ E num / \epsilon \varphi 2$ . | $ E Z L / \epsilon \varphi 2$ . |
---|---|---|---|---|---|

(0,4,1) | 1.309 307 341 | 0.426 775 579 0 | $ 4.437 \u2009 333 \u2009 386 \xd7 10 \u2212 2$ | $ 4.421 \xd7 10 \u2212 2$ | $ 4.440 \xd7 10 \u2212 2$ |

(0,6,1) | 0.937 697 587 | 0.318 769 489 4 | $ 2.852 \u2009 976 \u2009 237 \xd7 10 \u2212 2$ | $ 2.834 \xd7 10 \u2212 2$ | $ 2.848 \xd7 10 \u2212 2$ |

(0,6,2) | 1.660 447 793 | 0.555 430 038 4 | $ 8.244 \u2009 384 \u2009 738 \xd7 10 \u2212 3$ | $ 8.009 \xd7 10 \u2212 3$ | $ 7.959 \xd7 10 \u2212 3$ |

(0,8,1) | 1.799 515 990 | 0.610 750 963 0 | $ 2.289 \u2009 522 \u2009 846 \xd7 10 \u2212 3$ | $ 2.252 \xd7 10 \u2212 3$ | ⋯ |

(0,8,2) | 0.726 234 926 | 0.253 129 413 2 | $ 1.827 \u2009 816 \u2009 203 \xd7 10 \u2212 2$ | $ 1.834 \xd7 10 \u2212 2$ | ⋯ |

(0,8,3) | 1.354 372 557 | 0.471 011 105 6 | $ 9.463 \u2009 982 \u2009 745 \xd7 10 \u2212 3$ | $ 9.649 \xd7 10 \u2212 3$ | ⋯ |

(m, l, k)
. | ω
. _{lmk} | S
. _{lmk} | $ E lmk / \epsilon \theta 2$ . | $ E num / \epsilon \theta 2$ . |
---|---|---|---|---|

(1,2,1) | 1.000 000 000 | 0.258 463 392 8 | $ 2.094 \u2009 395 \u2009 103 \xd7 10 \u2212 1$ | $ 2.097 \xd7 10 \u2212 1$ |

(1,4,1) | 1.708 023 904 | 0.504 571 758 0 | $ 4.869 \u2009 232 \u2009 497 \xd7 10 \u2212 2$ | $ 4.827 \xd7 10 \u2212 2$ |

(1,4,2) | 0.820 008 840 | 0.316 519 274 5 | $ 7.972 \u2009 175 \u2009 759 \xd7 10 \u2212 3$ | $ 1.015 \xd7 10 \u2212 2$ |

(1,4,3) | 0.611 984 936 | 0.180 786 895 0 | $ 3.807 \u2009 897 \u2009 735 \xd7 10 \u2212 2$ | $ 3.867 \xd7 10 \u2212 2$ |

(1,6,1) | 0.440 454 452 | 0.138 816 148 5 | $ 1.465 \u2009 674 \u2009 965 \xd7 10 \u2212 2$ | $ 1.408 \xd7 10 \u2212 2$ |

(1,6,2) | 1.306078717 | 0.429 554 895 2 | $ 2.188 \u2009 323 \u2009 255 \xd7 10 \u2212 2$ | $ 2.198 \xd7 10 \u2212 2$ |

(1,6,3) | 1.861 684 240 | 0.560 603 944 5 | $ 1.724 \u2009 853 \u2009 435 \xd7 10 \u2212 2$ | $ 1.736 \xd7 10 \u2212 2$ |

(1,6.4) | 0.537 333 891 | 0.202 556 729 8 | $ 5.284 \u2009 148 \u2009 108 \xd7 10 \u2212 3$ | $ 5.706 \xd7 10 \u2212 3$ |

(1,6,5) | 1.404 216 852 | 0.526 006 866 5 | $ 7.314 \u2009 363 \u2009 353 \xd7 10 \u2212 4$ | $ 1.005 \xd7 10 \u2212 3$ |

(m, l, k)
. | ω
. _{lmk} | S
. _{lmk} | $ E lmk / \epsilon \theta 2$ . | $ E num / \epsilon \theta 2$ . |
---|---|---|---|---|

(1,2,1) | 1.000 000 000 | 0.258 463 392 8 | $ 2.094 \u2009 395 \u2009 103 \xd7 10 \u2212 1$ | $ 2.097 \xd7 10 \u2212 1$ |

(1,4,1) | 1.708 023 904 | 0.504 571 758 0 | $ 4.869 \u2009 232 \u2009 497 \xd7 10 \u2212 2$ | $ 4.827 \xd7 10 \u2212 2$ |

(1,4,2) | 0.820 008 840 | 0.316 519 274 5 | $ 7.972 \u2009 175 \u2009 759 \xd7 10 \u2212 3$ | $ 1.015 \xd7 10 \u2212 2$ |

(1,4,3) | 0.611 984 936 | 0.180 786 895 0 | $ 3.807 \u2009 897 \u2009 735 \xd7 10 \u2212 2$ | $ 3.867 \xd7 10 \u2212 2$ |

(1,6,1) | 0.440 454 452 | 0.138 816 148 5 | $ 1.465 \u2009 674 \u2009 965 \xd7 10 \u2212 2$ | $ 1.408 \xd7 10 \u2212 2$ |

(1,6,2) | 1.306078717 | 0.429 554 895 2 | $ 2.188 \u2009 323 \u2009 255 \xd7 10 \u2212 2$ | $ 2.198 \xd7 10 \u2212 2$ |

(1,6,3) | 1.861 684 240 | 0.560 603 944 5 | $ 1.724 \u2009 853 \u2009 435 \xd7 10 \u2212 2$ | $ 1.736 \xd7 10 \u2212 2$ |

(1,6.4) | 0.537 333 891 | 0.202 556 729 8 | $ 5.284 \u2009 148 \u2009 108 \xd7 10 \u2212 3$ | $ 5.706 \xd7 10 \u2212 3$ |

(1,6,5) | 1.404 216 852 | 0.526 006 866 5 | $ 7.314 \u2009 363 \u2009 353 \xd7 10 \u2212 4$ | $ 1.005 \xd7 10 \u2212 3$ |

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^{27} 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

^{26}such that $ \omega lmk ( E ) = \omega lmk + S lmk E 1 / 2$, also reported in Tables I and II. Note that Zhang and Liao

^{27}did not consider the viscous correction to the eigen-frequency, whereas Greenspan

^{26}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.

## IV. NUMERICAL METHOD

*e*and

*f*are expanded in terms of associated Legendre functions as

*m*= 0, 1 for longitudinal and latitudinal libration, respectively. Inserting Eq. (33) into (32), taking the

*r*-components of the first and second curls of Eq. (32), and using properties of the associated Legendre functions then yields a block tri-diagonal matrix structure,

^{30,31}where each block contains the details of the radial expansion.

*r*and

*r*

^{2}ensures that the leading-order regularity conditions at the origin are satisfied. The difference between

*r*for

*e*vs

*r*

^{2}for

*f*is due to the odd vs even Legendre functions in Eqs. (34) and (35), which induces different symmetries in

*r*, as shown, for example, in the appendix of Ref. 33. The Chebyshev variable

*x*and the actual radial coordinate

*r*are related by either

*r*=

*x*, or a slightly more complicated formula discussed below.

The boundary conditions associated with Eq. (32) are Eqs. (10) or (11) as before, except that the previous amplitude factors $ \Delta \epsilon \varphi $ and $ \Delta \epsilon \theta $ are now absent since the problem is linear. Working through the details, and recalling also that the $ cos \u2009 ( \omega t )$ factors there are already accounted for by the $ e i \omega 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 \u2032 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*, 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.

_{nk}First, the relationship *r* = *x* was altered to $ r = x + \delta ( x \u2212 x 3 )$. For $ 0 < \delta < 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 $ \delta = 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*^{2}, so reducing *K* by 30% allows a substantial increase in the angular resolution *N*, for a given core memory size.

*ω*enters as a parameter in the matrix $A$, since it appears in Eq. (32). Thus, having solved the problem for a given

*ω*, if we then change

*ω*by some tiny increment, do we really have to start all over, or can we somehow use the fact that we already have the LU decomposition of $A$ at the original frequency? To see how we can use this information, let $ A 0$ be the matrix at frequency

*ω*

_{0}, and then consider the problem at frequency $ \omega 0 + \Delta \omega $. Noting how

*ω*enters Eq. (32), the new matrix will clearly be of the form $ A 0 + \Delta \omega \u2009 A \u0303$, where $ A \u0303$ is simply the discretized version of the term $ i u$. We, thus, wish to solve the problem

The final procedure is thus as follows: first compute the LU decomposition of $ A 0$, corresponding to some frequency *ω*_{0}. This is accomplished using NAG routine F07BRF, and becomes computationally very expensive at high resolutions, scaling as *NK*^{3}. Next iterate according to Eq. (38), and for some $ \xb1 \Delta \omega $ range centered around *ω*_{0}. Each iteration is accomplished using NAG routine F07BSF and is computationally much cheaper, scaling only as *NK*^{2}. How large $ | \Delta \omega |$ can be taken before the iteration fails to converge, and a new LU decomposition corresponding to a new *ω*_{0} 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 *ω*_{0} 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 $ \u223c 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*^{2}) vs a new LU decomposition (*NK*^{3}) means that the iterative procedure is far faster, even if on average a dozen or so iterations might be required for each value of $ \Delta \omega $.

With these improvements in place, the code was then run at resolutions $ N \xd7 K = 80 \xd7 170$ for Ekman number $ E \u2265 10 \u2212 8$, 100 × 270 for $ E = 10 \u2212 9$, and 130 × 430 for $ E = 10 \u2212 10$, where according to Eqs. (34) and (35), the maximum spherical harmonic degree is then 2*N* in each case. Recall also that the memory requirements scale as *NK*^{2} and the CPU time for each LU decomposition scales as *NK*^{3}. As *E* is reduced, this iterative procedure (38) also converges for increasingly small $ | \Delta \omega |$ 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*^{3}. The final result is that $ E \u2265 10 \u2212 8$ is quite straightforward on a standard PC, $ E = 10 \u2212 9$ still fits on a PC in terms of memory, but requires considerably more CPU time, whereas $ E = 10 \u2212 10$ requires a machine with larger memory, and becomes so CPU-intensive that only a small fraction of the full range $ \omega \u2208 [ 0 , 2 ]$ could be scanned.

## V. RESULTS

### A. Frequency response and flow structure

Figures 3 and 4 show the time-averaged kinetic energy $ E num / \epsilon \varphi , \theta 2$ as a function of the libration frequency *ω*. The upper panels in each figure show scans over the entire frequency range $ \omega \u2208 [ 0 , 2 ]$, computed with a spacing $ \Delta \omega = 2 \xd7 10 \u2212 3$, and down to Ekman numbers $ E = 10 \u2212 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 $ \omega = 2 \u2009 sin \u2009 ( \pi / 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 , \u2026 , 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 \u2212 7$ in the top rows with $ E = 10 \u2212 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, $ ( \omega n + 1 \u2212 \omega 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 \u2192 0$, the entire interval $ \omega \u2208 [ 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 $ \omega 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 $ \theta = arcsin ( \omega / 2 )$ with the vertical axis. As an example, we show results for longitudinal libration at $ \omega = 2 \u2009 sin \u2009 ( \pi / 8 ) = 0.765$ in Fig. 7. As explained in the study by Rieutord *et al.*,^{34} for frequencies of the form $ \omega = 2 \u2009 sin \u2009 ( \pi \u2009 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,^{35} 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.^{6}

### B. Convergence toward the asymptotic regime

Although numerical simulations such as the ones performed here allow us to compute Ekman numbers as small as $ 10 \u2212 10$, the values in planetary cores can be several orders of magnitude smaller still, as small as $ O ( 10 \u2212 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 $ \omega lmk + S lmk E 1 / 2$ proposed by Greenspan.^{26} 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 $ \Delta \omega lmk = \omega num \u2212 \omega lmk \u2212 S lmk E 1 / 2$. As shown in Fig. 9, these errors are bounded by $ \Delta \omega lmk < O ( E )$, as expected from the theory proposed by Greenspan.

^{26}

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.

### C. The trough regions

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 $ \omega = 2 \u2009 sin \u2009 ( \pi \u2009 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 $ \omega = 2$. The nearest peaks at $ \omega \u2212$ and $ \omega +$ define the width of the trough as $ \delta = | \omega + \u2212 \omega \u2212 |$. The characteristic energy of a trough is defined at the central frequency $ f = 1 2 ( \omega + + \omega \u2212 )$, 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, $ \delta \u221d E \alpha $ and $ E num / \epsilon \varphi , \theta 2 \u221d E \alpha \u2032$, 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 / \epsilon \varphi , \theta 2 \u221d 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 $ \delta \u221d E 0.23 \xb1 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 $ \omega t = 2$. That is, we rescale the energy by $ E 1 / 2$, and the frequency as $ ( \omega \u2212 \omega t ) / E 0.23$. The troughs for the three different values $ E = 10 \u2212 8 , \u2009 10 \u2212 9 , \u2009 10 \u2212 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.

## VI. CONCLUSIONS AND PERSPECTIVES

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,^{27} which form a dense set on the frequency interval $ \omega \u2208 [ 0 , 2 ]$. Including the first-order viscous correction to the eigenfrequencies, $ \omega = \omega 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 \u2212 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;^{27} 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 \u223c 0.23 \u223c 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 [ $ \xi e q = ( a 2 \u2212 b 2 ) / a 2$] and polar [ $ \xi pol = ( a 2 \u2212 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;^{27} other conventions exist in the literature leading to different inequalities, for instance adopting the conventions of Grannan *et al.*^{19} [ $ \beta = | a 2 \u2212 b 2 | / ( a 2 + b 2 )$] or the polar flattening $ \eta = ( a \u2212 b ) / a$ (Refs. 37 and 38) the condition becomes $ \beta , \eta \u226a 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 $ \Delta \theta / E 1 / 2$ (Refs. 25 and 27) near its resonant frequency, i.e., $ | \omega \u2212 2 / ( 2 \u2212 \xi pol 2 ) | \u226a 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,^{23} similarly to what is found for precession.^{39}

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 $ \epsilon \u223c 10 \u2212 4$. A typical velocity is obtained from our Tables I and II as $ U rms = \epsilon \Omega 0 R ( E lmk / \epsilon 2 ) / ( 4 \pi / 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.^{34} 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.^{40}

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: INERTIAL MODES IN A SPHERE

^{26,41,42}There is an infinite set of eigenvalues

*ω*, which are the real roots of the transcendental equation,

_{lmk}*l*and order

*m*. For each combination of $ l > 0 , \u2212 l \u2264 m \u2264 l$, the characteristic equation (A2) has a finite number of zeros. Hence, we can label each mode by an integer triplet {

*l*,

*m*,

*k*}. The corresponding pressure eigenfunction is given by

*r*= 1, the pressure eigenmode takes the form

^{26}Eq. 2.7.2)

### APPENDIX B: SOLUTION OF THE BOUNDARY LAYER EQUATIONS

In this appendix, we derive the boundary layer flow $ u \u0303$ that matches a given interior flow $ u 0$ to a boundary condition $ u \Sigma $ of the form (10) and (11). Following the discussion in Sec. III, the libration frequency *ω* lies within the inertial mode spectrum $ \omega \u2208 [ 0 , 2 ]$ and the leading-order interior flow is taken to be of the form $ u 0 = A lmk u lmk$ with $ u lmk$ an inertial mode whose eigenfrequency is *ω*.

*O*(1), and this implies that $ u \u0303 r$ is $ O ( E )$. Thus, the leading-order term of the corrective velocity $ u \u0303$ is purely tangential.

*r*= 1, i.e.,

*ξ*= 0 eventually gives

^{26}Kudlick,

^{43}these singularities can be ignored when computing the Ekman pumping.

## REFERENCES

*Solar System Dynamics*

*Fluid Mechanics of Planets and Stars*

*Theory and Modeling of Rotating Fluids: Convection, Inertial Waves and Precession*