Dynamical tidal Love numbers of rapidly rotating planets and stars

Tidal interactions play an important role in many astrophysical systems, but uncertainties regarding the tides of rapidly rotating, centrifugally distorted stars and gaseous planets remain. We have developed a precise method for computing the dynamical, non-dissipative tidal response of rotating planets and stars, based on summation over contributions from normal modes driven by the tidal potential. We calculate the normal modes of isentropic polytropes rotating at up to $\simeq90\%$ of their critical breakup rotation rates, and tabulate fits to mode frequencies and tidal overlap coefficients that can be used to compute the frequency-dependent, non-dissipative tidal response (via potential Love numbers $k_{\ell m}$). Although fundamental modes (f-modes) possess dominant tidal overlap coefficients at (nearly) all rotation rates, we find that the strong coupling of retrograde inertial modes (i-modes) to tesseral ($\ell>|m|$) components of the tidal potential produces resonances that may be relevant to gas giants like Jupiter and Saturn. The coupling of f-modes in rapid rotators to multiple components of both the driving tidal potential and the induced gravitational field also affect the tesseral response, leading to significant deviations from treatments of rotation that neglect centrifugal distortion and high-order corrections. For very rapid rotation rates ($\gtrsim 70\%$ of breakup), mixing between prograde f-modes and i-modes significantly enhances the sectoral ($\ell=|m|$) tidal overlap of the latter. The tidal response of very rapidly rotating, centrifugally distorted planets or stars can also be modified by resonant sectoral f-modes that are secularly unstable via the Chandrasekhar-Friedman-Schutz (CFS) mechanism.


INTRODUCTION
The question of how a self-gravitating fluid body (such as a star or a planet) responds to the gravitational potential of an orbiting satellite becomes complicated, especially at a quantitative level, when the body in question rotates at a significant fraction of the "breakup" rotation rate, Ω d = (GM/R 3 eq ) 1/2 (where M and R eq are the body's mass and equatorial radius). Recent measurements of the non-spherical gravity fields of Jupiter and Saturn provided by Juno and Cassini promise insight into this problem, since the two planets possess relatively rapid bulk rotation rates of Ω 0.3Ω d and 0.4Ω d (respectively). The satellite data have facilitated the placement of observational constraints on the (dissipative as well as non-dissipative) tidal response for both planets (Lainey et al. 2009(Lainey et al. , 2012(Lainey et al. , 2017(Lainey et al. , 2020Durante et al. 2020). In particular, gravity field measurements allow for the inference of the so-called 'potential Love numbers' k m -response functions that measure the ratios between the harmonic components of the gravitational potentials associated with tidal deformation of the planet, and components of the potentials imposed by planetary satellites (e.g., Ogilvie 2014). These unprecedented measurements have the potential to constrain planetary internal structures, and invite further theoretical investigations into the dynamical tides of rapidly rotating planets and stars.
Previous efforts in the planetary sciences community have primarily adopted 3D 'Concentric MacLaurin Spheroid' (CMS) calculations of distorted planetary structure (Wahl et al. 2017(Wahl et al. , 2020Nettelmann 2019). These calculations rely on a hydrostatic approximation, neglecting the finite tidal forcing frequencies associated with Jupiter and Saturn's moons. Significant discrepancies between the observed Love numbers and those calculated with the CMS method for Jupiter indicate that dynamical, non-hydrostatic effects modify the planet's tidal response. Idini & Stevenson (2021) showed through direct computations incorporating the Coriolis force that dynamical tides can explain the Love number discrepancies, in particular for the quadrupolar k 22 . Lai (2021) demonstrated that the same results can be achieved more efficiently with a phase space expansion in the normal mode oscillations of the planet. However, Lai (2021) employed a perturbative treatment of rotation (to linear order in Ω), and both Idini & Stevenson (2021) and Lai (2021) ignored the ∼ O(Ω 2 ) effects of the centrifugal distortion of the equilibrium state on the fundamental 'f-modes' that typically dominate the gravitational response. Jupiter's rapid rotation makes the accuracy of perturbative treatments unclear.
In this paper, we use a non-perturbative method to compute the normal modes of rapidly rotating fluid bodies. Focusing on n = 1 polytropes with rotation rates Ω 0.9Ω d , we tabulate fits to the mode properties required to calculate the (frequency-dependent) potential Love numbers k m with azimuthal wave numbers m = 1, 2, 3, 4 and spherical harmonic degrees = 2, 3, 4, 5, 6 (for even − m). In addition to giant planets, whose structures can be well approximated by n = 1 polytropes, our results can also be applied approximately to rapidly rotating neutron stars.
As in the non-rotating limit, f-modes dominate the non-dissipative tidal response at all rotation rates, except near resonances where the tidal forcing frequency matches the frequency of a mode. In particular, inertial modes (i-modes) restored by the Coriolis force have frequencies that can resonate with the tidal forcing. The i-modes also possess a strong tesseral ( > |m|) overlap that is often not considered in studies of tidal interactions focused on, e.g., the quadrupolar component of the perturbing potential (but see, e.g., Ogilvie 2013;Braviner & Ogilvie 2015). This tesseral coupling may be relevant to discrepant observed and predicted values of tesseral Love numbers (such as k 42 ; Durante et al. 2020). In addition, our non-perturbative inclusion of f-modes' overlap with multiple spherical harmonics due to centrifugal distortion also leads to significant (non-resonant) deviations of the tesseral Love numbers k m+2,m from the calculations of Lai (2021). This paper is structured as follows. In Section 2 we describe our methods for calculating centrifugally distorted stellar models, normal mode oscillations and po-tential Love numbers. In Section 3 we review the general properties of the oscillation modes considered. We present results in Section 4, and conclude in Section 5.

METHODS
In this section we outline our methods for calculating oblate models of rotating polytropes (Section 2.1), and their normal mode oscillations (Section 2.2). We then introduce the relevant quantities required to characterize the response of a fluid body to an imposed tidal potential (Section 2.3). Throughout the paper, we refer both to spherical and cylindrical polar coordinates denoted by (r, θ, φ) and (R, φ, z), respectively.

Stellar and planetary models
We consider rapidly rotating, neutrally stratified (i.e., fully and efficiently convective) polytropes characterized by the barotropic equation of state P ∝ ρ 1+1/n , focusing on the polytropic index n = 1. The convective envelopes in current models of Jupiter and Saturn match relatively well with n = 1 polytropes, although they may differ both in details and in the deep interior (e.g., . Newtonian polytropes with n = 1 are also frequently used to approximate the interiors of neutron stars (e.g., Passamonti et al. 2009b;Xu & Lai 2017), since this index reproduces predictions by equations of state for nuclear matter of an approximately constant neutron star radius over a range of masses around 1.4M (see, e.g., Lattimer & Prakash 2007). The inclusion of super-fluidity (Passamonti et al. 2009a) and general relativity (De Pietri et al. 2018) are, however, required for more realistic representations.
We use a two-dimensional pseudospectral method to calculate axisymmetric but oblate polytropic models. It is reasonable to start from an axisymmetric basic state when computing the linear tidal response; the tidal deformations of Jupiter and Saturn, for instance, are very weak in comparison with their rotational deformations. Appendix A provides the details of these calculations. In brief, we start by calculating one-dimensional, nonrotating solutions to the Lane-Emden equation, and then use an iterative scheme to compute centrifugally distorted models. We verify the accuracy of our model calculations by considering virial errors.
In addition to the polytropic index n and rotation rate Ω, different models are characterized by the equatorial and polar radii R eq and R pol (resp.), the central value for the pseudo-enthalpy H = ρ −1 dP = (1+n)P/ρ, the total mass, and the ratio between kinetic and potential energy. Unless otherwise stated, we use units with G = M = R eq = 1 so that angular velocities take on units of Ω d = (GM/R 3 eq ) 1/2 (Ω d is frequently referred to as the "dynamical frequency").

Mode calculations
In the limit of zero rotation, the partial differential equations (PDEs) governing linear perturbations to a self-gravitating fluid in hydrostatic equilibrium are separable in terms of spherical harmonics. Rotation disrupts this separability, since both the Coriolis force and centrifugal distortion of the equilibrium state break spherical symmetry. Perturbative treatments, in which eigenfunctions and frequencies are expanded in powers of Ω, prove adequate when rotation is slow (see, e.g., Unno et al. 1989). However, the normal mode oscillations of fluid bodies rotating at an appreciable fraction of the equatorial Keplerian frequency Ω d require more complete treatments, as do modes with frequencies comparable to the rotation rate (even when Ω Ω d ). Non-perturbative methods (e.g., Lignières et al. 2006;Reese et al. 2006Reese et al. , 2009Reese et al. , 2013Reese et al. , 2021Ouazzani et al. 2012) treat the linearized PDEs as fundamentally nonseparable, eschewing assumptions of spherical symmetry and instead solving an inherently two-dimensional problem. Having calculated equilibrium pressure and density profiles P 0 (r, θ) and ρ 0 (r, θ) for a given polytropic index n and rotation rate Ω, we use the spectral method described in Dewberry et al. (2021) to solve for adiabatic normal modes with the time dependence exp[−iωt] in the rotating frame. We refer interested readers to that work for the details of these calculations, but make note of a few points salient to this investigation: In a rapidly rotating planet or star, normal mode eigenfunctions involve Eulerian perturbations v, ρ , P , Φ to the fluid velocity, density, pressure and gravitational field (resp.) that depend non-trivially on radius and polar angle; given the non-separability of the governing PDEs, mode eigenfunctions cannot be associated with a single spherical harmonic Y m . Instead, the non-perturbative approach involves a more agnostic expansion in series of spherical harmonics. For example, we represent the gravitational potential perturbation Φ (r, θ, φ, t) with an expansion of the form Note that because the equilibrium states are axisymmetric, the equations are still separable in φ, and each oscillation mode can be identified with a unique azimuthal wavenumber m. We adopt the convention of strictly positive m, with prograde (retrograde) propagation in the rotating frame differentiated by positive (negative) ω. The coefficients Φ (ζ) in this expansion depend on a quasi-radial coordinate ζ associated with a nonorthogonal coordinate system originally introduced by Bonazzola et al. (1998), which is constructed to match the oblate surface r s (θ) of the rotating fluid body (at ζ/R eq = 1), and to relax to spherical coordinates both at the origin (ζ = 0) and at the edge of the exterior vacuum region included in the computational domain (defined by ζ/R eq = 2).
We solve Laplace's equation ∇ 2 Φ = 0 in the external vacuum, which we include for the purposes of boundary condition application: at the oblate surface r s we impose the continuity of Φ and its gradient, in addition to the free surface condition of a vanishing Lagrangian pressure perturbation. Following Reese et al. (2006), we apply the actual boundary condition for the gravitational potential on the spherical surface ζ = r = 2R eq , matching harmonic coefficients Φ to the analytical solutions that vanish at infinity. We enforce the standard regularity conditions at the origin (see Chapter III in Unno et al. 1989) implicitly.

Tides and Love numbers
We work in a frame centred on, and rotating with the tidally disturbed planet or star (hereafter "planet") under consideration. In the vicinity of this planet, a satellite with mass M and orbital separation a on a circular and co-planar orbit produces a gravitational potential that can be expanded as where the tidal forcing frequency is and Ω o = [G(M + M )/a 3 ] 1/2 is the orbital frequency. The harmonic coefficients U m are then given by (Jackson 1962;Press & Teukolsky 1977) where W m = 0 for non-integer ( +m)/2, and otherwise The tidal potential U drives density perturbations in the planet, which then generate an external potential perturbation that can be expanded as The potential Love number k m is then defined as the ratio (e.g., Ogilvie 2014) In the absence of any dissipative processes, k m is realvalued; we ignore the Love numbers' imaginary parts, but note that they can be significant near resonances, or in rotating planets with solid cores.
In the linear theory, the tidal response of the planet can be written as an expansion in normal modes (e.g., Schenk et al. 2002;Lai & Wu 2006). Specifically, we write where k α m describes the contribution of a mode labelled by α. Following the procedure of (Lai 2021), with the inclusion of higher-order centrifugal effects we find Defining the inner product ξ, ξ = V ρ 0 (r, θ)ξ * · ξdV for the Lagrangian displacement ξ = (−iω) −1 v, Equation 9 involves the quantity along with the tidal overlap coefficient Note that all quantities in Equation 9 are dimensionless, under our adoption of units with G = M = R eq = 1. Each coefficient Q α m characterizes the spatial overlap of the mode α with the m component of the tidal potential. In the last equality of Equation 11, the quantity in square brackets gives the coefficient of degree in a harmonic expansion in spherical coordinates of the mode's gravitational potential perturbation (at radius r = R eq ). This can differ significantly from the coefficient Φ α (ζ)| ζ=1 in the expansion of Equation 1, since the rotating body is non-spherical.
Equation 9 is more general than Eq. 11 in Lai (2021), in that it includes an additional sum over spherical harmonic degrees indexed by . This sum is necessary in order to fully account for ∼ O(Ω 2 ) centrifugal effects on the planet's structure and modes. In a spherically symmetric planet, all terms with = would vanish, and the ratio U m /U m would collapse to one. More generally, where q = M /M is the mass ratio [recall that a − = (a/R eq ) − in our units].
, which increases as the separation a increases (Ω o decreases) for − > 0. This leads to a divergence in some k m (specifically tesseral Love numbers with > m) as a → ∞. For example, consider k 42 ( = 4, m = 2): in a centrifugally flattened planet, modes excited by the quadrupolar U 22 potential contribute to δΦ 42 . As a → ∞, the quadrupolar contribution to δΦ 42 can in fact dominate the contribution from U 42 , simply because |U 22 | |U 42 |. The ratio k 42 = δΦ 42 /U 42 then "diverges," even though the gravitational response δΦ 42 to the satellite remains finite. The unphysical nature of this divergence, which we discuss further in Section 4.3, suggests that Love numbers as defined by Equation 7 to relate like-coefficients (i.e., those with the same ) in harmonic expansions may not provide the best description of the tidal response of rotationally deformed bodies.
Computing k α m therefore requires knowledge of i) the mode frequency ω α , ii) α , and iii) Q α m (for multiple ). The need for α can be eliminated by normalizing ξ α such that α /ω α = 1 (see, e.g., Dewberry et al. 2021), but we adopt the normalization ξ α , ξ α = 1 (as in Lai 2021) because in the perturbative regime, α is simply the mode frequency in the nonrotating limit. This normalization additionally ensures that the rotational correction to the f-modes' Q α m is of order O(Ω 2 ).

MODE COMPENDIUM
In this paper we focus on the properties of fundamental modes (f-modes), inertial modes (i-modes), and to a lesser extent acoustic modes (p-modes), deferring consideration of stable stratification and internal gravity modes (g-modes) to later works. We additionally focus on modes with even equatorial parity (i.e., even − m), which are the most relevant to planets with spin-aligned satellites. The following subsections briefly review the relevant properties of these different types of modes.

Fundamental modes
Fundamental modes (f-modes), also known as surface gravity modes, manifest even in homogeneous and incompressible stellar models. They usually possess the simplest eigenfunctions of all non-radial oscillations, with no nodes in the radial direction, and have the largest magnitude tidal overlap coefficients Q α m . We label f-modes with the notation f ± m , where superscript +/− denote prograde/retrograde propagation, and subscript is the spherical harmonic degree uniquely associated with each mode in the non-rotating limit. With increasing rotation the f-mode eigenfunctions come to involve a large number of spherical harmonic degrees, but generally maintain a clear correspondence to the non-rotating .

Acoustic modes
Compressible (n > 0) polytropic models also support acoustic oscillations (p-modes). Their contributions to the potential Love numbers are usually unimportant, because of their high frequencies and small |Q α m | values in comparison with f-modes. We nonetheless include a few for completeness; employing the notation p ± mnp , where n p is the number of radial nodes in the nonrotating limit, we track even-parity p-modes with < 7, ( −m)/2+n p ≤ 3 for Ω = 0. Despite minimal frequency shifts due to the Coriolis force, strong centrifugal distortion at high rotation rates can complicate p-mode eigenfunctions and frequency spectra (Lignières & Georgeot 2009;Reese et al. 2009Reese et al. , 2013. However, the low-degree, low-order p-modes that we consider generally maintain a regular spectrum with increasing rotation.

Inertial modes
Rotating isentropic stellar models also support global inertial modes (i-modes) that are restored by the Coriolis force, and occupy the low-frequency regime |ω| < 2Ω (Bryan 1889;Greenspan 1968). In the absence of a solid core and/or stable stratification, i-modes form a regular, dense spectrum that can be characterized analytically in the simplest cases: ignoring centrifugal distortion, the i-modes of a rotating fluid body with a homogeneous or power-law density profile have ρ ∝ P m (x 1 )P m (x 2 ), where P m are Legendre polynomials, and (x 1 , x 2 ) are ellipsoidal coordinates that depend on the parameter ω/(2Ω) (Wu 2005). Lindblom & Ipser (1999) and Braviner & Ogilvie (2014) provided similar analytical solutions for the i-modes (and f-modes) of non-spherical Maclaurin spheroids, while  and  investigated inertial modes via asymptotic and numerical approaches (respectively).
For each and m in eigensolutions ∝ P m (x 1 )P m (x 2 ), Wu (2005) found − m eigenfrequencies. The different frequencies correspond to meridional eigenfunctions with differing numbers of nodes with respect to ellipsoidal generalizations of both cylindrical radius R (n 1 ), and cylindrical z (n 2 ). These mode numbers satisfy − m = 2(n 1 + n 2 ) for i-modes with even equatorial parity, and − m = 2(n 1 + n 2 ) − 1 for odd equatorial parity. We adopt the labelling scheme i ± mn1n2 , 1 where + and − again refer to prograde and retrograde propagation. Our focus on modes with even equatorial parity excludes so-called 'r-modes,' a subset of odd-parity, retrograde inertial modes (e.g., Schenk et al. 2002). These r-modes cannot be excited in co-planar systems, but may be relevant when there is a significant spin-orbit misalignment (Ho & Lai 1999;Braviner & Ogilvie 2015;Xu & Lai 2017).
The meridional cross sections in Figure 1 show ρ (r, θ) for even-parity inertial modes with n 1 , n 2 ∈ [0, 2], calculated for an n = 1, Ω/Ω d 0.3 polytrope. We see that n 1 (n 2 ) can also be counted as the number of surface zeros between the pole (equator) and a 'critical latitude' at which the mode wavelengths shortens, and the amplitude approaches a maximum. For each mode, the black dashed line indicates the WKB approximation θ = arccos[|ω|/(2Ω)] to this critical latitude.
The i-modes with the largest Q α m are those with the longest wavelengths (i.e., smallest n 1 +n 2 ), in particular the oscillations i − m01 and i + m10 with n 1 + n 2 = 1 (the first and third cross sections in the top row of Figure 1). For convenience, we sometimes refer to these oscillations as the retrograde and prograde inertial modes, respectively (for comparison, i − 201 and i + 210 are the same as the j = 3, m = 2 modes described in table IV of Xu & Lai 2017). For completeness, for each m we additionally consider the four shorter-wavelength modes with n 1 + n 2 = 2.
We note that the presence of a rigid core complicates inertial wave propagation (Rieutord & Valdettaro 1997;Ogilvie 2005;Goodman & Lackner 2009), by introducing singularities on the inner boundary that disrupt the regular spectrum illustrated by Figure 1. In the (apparent) absence of the global inertial modes of coreless models, inertial wave attractors do appear as spatially periodic features. However, these attractors do not generally correspond to the largest peaks in (frequency-dependent) tidal dissipation (Ogilvie 2009;Rieutord & Valdettaro 2010;Ogilvie 2013).

RESULTS
In this section we present the results of our mode and Love number calculations for rapidly rotating, n = 1 polytropes. Section 4.1 demonstrates the relative importance of the different modes described in Section 3 to the tidal response, Section 4.2 characterizes the variation of mode frequencies and tidal coefficients with increasing rotation rate, and Section 4.3 presents representative Love number calculations relevant to tides in rotating gas giants. Lastly, Section 4.4 describes additional results for more rapid rotation.

Mode relevance
The first question to ask is which modes contribute most significantly to which Love numbers. As might be gleaned from Equation 9, mode contributions k α m to a given k m are largest for i) large Q α m , ii) small α , and iii) frequencies ω α close to resonance with ω m (i.e., small ω α − ω m ). Since α is approximately equal to the zerorotation mode frequency (for f and p-modes), i) and iii) prove to be the differentiating factors. Figure 2 demonstrates the relative importance of different oscillations to an n = 1 polytrope rotating wwith Ω/Ω d 0.3. The left and right-hand panels plot |Q α 22 | and |Q α 42 | (resp.) vs. |ω α | for a selection of f-modes (blue), p-modes (green) and i-modes (orange). The filled circles (plus signs) indicate prograde (retrograde) modes with positive (negative) ω α .
The plots show that for a rotation rate typical of gas giants (Ω/Ω d 0.3, 0.4 respectively for Jupiter and Saturn) the f-modes f ± m (associated with spherical harmonic Y m in the absence of rotation) produce Q α m orders of magnitude larger than any of the other modes. P-modes possess reasonably large Q α m , but have large frequencies that are unlikely to come anywhere near resonance with satellites. Inertial modes, on the other hand, have smaller Q α m that can nevertheless be compensated by frequencies falling much closer to resonance with realistic tidal forcing frequencies.
Figure 2 also demonstrates coupling across spherical harmonics by rotation, which produces non-zero Q α m for degrees other than that associated with each mode's non-rotating . For instance, when Ω = 0 the only 3. Filled circles and plus signs denote prograde and retrograde modes, respectively. We label each point according to the conventions described in Section 3: f ± m for f-modes, p ± mnp for p-modes, and i ± mn 1 n 2 for i-modes.
non-zero tidal coupling factor for f ± 22 is Q 22 , but with rotation the modes have non-zero Q 42 , Q 62 , ... that decrease in magnitude with increasing (see Appendix C for a more detailed breakdown of Q α m calculations for different modes). As the f-mode eigenfunctions come to involve more spherical harmonics, they overlap more strongly with multiple components of the tidal potential. For moderate to rapid rotation rates, a given k m can no longer be reasonably approximated through sole consideration of the corresponding f ± m , since every mode overlaps non-negligibly with multiple components of both the tidal potential and the induced response.
Comparison of the left and right panels in Figure 2 reveals an important aspect of inertial mode overlap with the tidal potential: i-modes of a given m couple most strongly to tesseral harmonics with degrees > m. For instance, the longest wavelength m = 2 i-modes (i − 201 and i + 210 ) produce larger |Q 42 | than |Q 22 |. This preferential coupling with tesseral harmonics is exact for homogeneous, slowly rotating fluid bodies (Ogilvie 2013). The i-modes of the distorted n = 1 polytropes considered here have nonzero sectoral overlap, but the preferential coupling with tesseral harmonics still holds approximately, leading to tesseral overlap coefficients that are larger by at least an order of magnitude.

Frequency and overlap variation with rotation
The top panels in Figure 3 plot frequencies ω (left), dominant tidal overlap coefficients Q m (middle), and coefficients (right) as a function of rotation rate for several f-modes and the longest wavelength (n 1 + n 2 = 1) imodes of the n = 1 polytrope. In this section we present analytic formulae (polynomial fits) to these quantities as a function of rotation rate.

Fits for f-modes
The frequencies and parameters of f-modes (f ± m ) can be expanded to order Ω 2 as where ω 0 is the non-rotating frequency. Note that has no linear dependence on Ω. The linear term AΩ arises due to the Coriolis force, and can be obtained by treating rotation as a perturbation. Defining the Lagrangian displacement for a mode α in the absence of rotation as ξ 0 α = (ξ rr + rξ ⊥ ∇ ⊥ )Y m (here ∇ ⊥ is the angular part of the gradient operator), with the normalization ξ 0 α , ξ 0 α = 1, the perturbative approach yields (e.g., Unno et al. 1989) The coefficient describing the overlap between f ± m and the m-tidal potential can be expanded as where Q 0 is the tidal overlap in the absence of rotation. Table 1 gives values obtained Table 1 and Table 2). Bottom panels: residuals associated with the polynomial fits. Due to the dominance of inertial mode coupling to tesseral (rather than sectoral) harmonics, in the center panels we plot Qm+2,m (rather than Qmm) for the i-modes. Table 1. Fits to f-mode frequencies ω ω0 + AΩ + BΩ 2 , tidal overlap coefficients Q α m Q0 + CΩ 2 , and coefficients ω0 + DΩ 2 . The mode notation is f ± m (see Section 3.1). All quantities are in units with G = M = Req = 1.  Barker et al. (2016), who found that for homogeneous spheroids the relative importance of centrifugal distortion to changes in f-mode frequencies (i.e., the ratio B/A) increases with (see also Section 3.1 of Ho & Lai 1999).
As noted in Section 4.1, for Ω > 0 a given f-mode f ± m also couples to tidal potentials associated with different spherical harmonic degrees ( ± 2, ± 4, ..). In the non-rotating limit the sign of Q 0 is unimportant, since only its square appears in Equation 9. For a rapidly rotating planet, however, the (relative) signs of the different Q α m associated with a mode α must be retained, since the sum over degrees in Equation 9 includes cross terms such as Q α m Q α ±2,m . Adopting the convention that Q 0 > 0, we find that the most important of the subdominant tidal overlap coefficients are the tesseral Q m+2,m associated with the sectoral f-modes f ± mm . Fitting our numerical results, we find The negative values of these overlap coefficients (relative to Q 0 ) are important to include, along with the signs of the W m coefficients. The additional overlap coefficients for non-sectoral f-modes also scale with Ω 2 [e.g., , but we find that their contributions to the relevant Love numbers are negligible. Table 2 provides coefficients for polynomial fits to the frequencies, overlap coefficients Q mm and Q m+2,m , and the -coefficients of the longest-wavelength inertial modes with even equatorial parity (i − m01 and i + m10 ). We Table 2. Fits to i-mode frequencies ω AΩ + BΩ 3 , tidal overlap coefficients Qmm CΩ 2 + DΩ 4 and Qm+2,m EΩ 2 + F Ω 4 , and coefficients GΩ. The mode notation is i ± mn 1 n 2 (see Section 3.3). We do not list Qmm for the two m = 1 i-modes because the tidal potential does not include an = m = 1 component.

Fits for i-modes
GΩ.
Note that the coefficients in these expansions are distinct from those introduced in Equations (13), (14), and (16). We choose the (larger amplitude) Q m+2,m to be positive. As indicated by the table, this leads to negative Q mm for the prograde inertial modes, and positive Q mm for the retrograde inertial modes.
Our results for the m = 2 i-modes agree with those given in Xu & Lai (2017) (compare with the j = 3, m = 2 modes in their Table IV), although they did not consider the (somewhat larger) Q 42 coefficients. To further compare with Xu & Lai (2017), we have computed a few m = 1 inertial modes with odd equatorial parity, namely the m = 1 r-mode (labelled i − 1r ), and the three m = 1, odd-parity i-modes corresponding to their j = 3, m = 1 modes. Note that because these modes have odd equatorial parity, they cannot couple to the tidal potential for systems with spin-orbit alignment (the relevant case for Jupiter and Saturn).
Our results, given in Table 3, show agreement with those of Xu & Lai (2017). Note that the m = 1 r-mode has a frequency exactly equal to −Ω (to the precision of our calculations) for all rotation rates considered in this paper. As explained by Xu & Lai (2017), this implies that the mode has zero inertial-frame frequency, constituting a spin-over perturbation. Table 3. Same as Table 3, but for odd-parity, m = 1 i-mode frequencies ω AΩ + BΩ 3 , tidal coupling factors Q21 CΩ 2 + DΩ 4 , and coefficients EΩ. The so-called 'hydrostatic' Love numbers can be computed by taking the limit of zero tidal forcing frequency ω m = m(Ω o − Ω) = 0, i.e., k hs m = k m (Ω 0 /Ω = 1).
This quantity is satellite independent in the limit of negligible mass ratio q = M /M , characterizing the intrinsic (linear) response of the planet to a hypothetical perturber with infinitesimal mass and an orbit synchronized to the spin rate.  This definition differs from that used in studies employing the CMS method that report hydrostatic Love numbers for multiple satellites; we reserve the term hydrostatic for the response that is stationary in the corotating frame, which for circular and co-planar systems with vanishing q can only be produced by a single orbital frequency Ω o and associated satellite separation a. A more direct comparison of our results to those of, e.g., Wahl et al. (2017Wahl et al. ( , 2020 could be made by assuming ω m = 0 in Equation 9, but then introducing multiple values of a in Equation 12 (this is inconsistent within our framework because fixing ω m = 0, and hence Ω o = Ω, in turn fixes a).
The analytically known k hs 22 = 15/π 2 − 1 of a nonrotating, n = 1 polytrope provides a zeroth order test of our computations; including only contributions from f ± 22 and p ± 221 , we find an error of 0.0001%. Table 4 lists hydrostatic Love numbers taken from the calculations shown in Figure 4 for a variety of and m, in addition to Table 4. Love numbers k m computed for an n = 1 polytrope with Ω/Ω d 0.30 (∼ Jupiter's rotation rate), both in the hydrostatic limit of zero forcing frequency (i.e., Ωo = Ω), and at the orbital frequencies associated with the Galilean satellites.  Table 4, the tesseral Love numbers with > m and m > 1 (e.g., k hs 42 ) exhibit the largest deviations from the hydrostatic values found by Lai (2021)  (25) The top panel shows dynamical corrections for k 31 and the sectoral Love numbers k mm , while the bottom panel shows corrections for the tesseral k m+2,m . In the following subsections we break down the contributions of different modes to different Love numbers.

The effect of f-modes on Love numbers
As indicated by the dashed colored lines in Figure 4 and Figure 5, the f-modes contribute relatively gradual frequency variation to k 31 , k 22 , k 33 and k 44 . The top panel of Figure 5 shows dynamical corrections that follow the same general trends found by Idini & Stevenson (2021) and Lai (2021); the inclusion of f-mode modification by centrifugal distortion in our computations affects the overall values of these Love numbers (i.e., k hs m ) more significantly than their dynamical corrections.
In contrast, the f-modes contribute dramatic frequency variation to the tesseral Love numbers with > m (and m > 1) that cannot be captured by O(Ω) treatments of rotation that neglect centrifugal distortion. The curves for k 42 , k 53 and k 64 in the log-log plot shown in Figure 4 indicate a power-law dependence ∝ Ω −4/3 o (black dotted line) as Ω o → 0. The origin of this power law can be traced back to Equation 12. Indeed, the expression for an arbitrary tesseral Love number k m+2,m is given by The second term in curly brackets indicates a crucial dependence on the amplitude and sign of the "crosscorrelation" Q α mm Q α m+2,m . Numerically, we find that this cross-correlation is both negative and non-negligible for the sectoral f-modes f ± mm , but negligible for most of the other modes (see Appendix C). The driving of sectoral f-modes by the sectoral component of the tidal potential then produces the dominant ∝ Ω −4/3 o dependence (which derives from the ratio U mm /U m+2,m ) illustrated for k 42 , k 53 and k 64 in Figure 4. Even steeper power laws prevail at low Ω o for larger > m + 2 (omitted for clarity), while k 31 does not show the same power law because the m = 1 sectoral f-mode is a trivial mode with zero frequency.
This inverse dependence on the orbital frequency of the satellite (similarly found in CMS calculations, and Juno observations: Wahl et al. 2017Wahl et al. , 2020Nettelmann 2019;Durante et al. 2020) implies that an infinitely distant satellite would produce an infinite tesseral Love number. However, this divergence does not have any physical consequence; it merely indicates that in a centrifugally distorted planet the response δΦ m+2,m to a distant satellite with small Ω o will be primarily driven by U mm , since |U mm | |U m+2,m | for large separation. δΦ m+2,m becomes large relative to U m+2,m , but not in an absolute sense: δΦ m+2,m /U mm remains finite, and in fact δΦ m+2,m decreases with increasing a. For very rapid rotation, it may be more useful to consider the full complement of potential Love numbers (see, e.g., Eq. 37 in Ogilvie 2013), rather than just the subset relating the components of the response and the tidal potential with the same and m.

The effect of i-modes on Love numbers
As illustrated in Figure 4 and Figure 5, the i-modes encounter resonances where their frequencies coincide with the tidal frequency. In the absence of dissipation, these resonances produce infinite peaks in the Love numbers; the peaks in Figure 4 and Figure 5 have a finite height because of the finite grid-spacing in Ω o . Dissipation would introduce physically meaningful peak amplitudes, and could additionally broaden the resonances. The resonances are narrow for the sectoral Love numbers k mm , since the inertial modes couple weakly to the sectoral tides, and wider for the tesseral k m+2,m . For instance, for m = 2 the longest wavelength retrograde mode i − 201 produces a relatively narrow resonance in k 22 at Ω o /Ω 0.45, and a more significant dynamical correction to k 42 at the same frequency. The retrograde inertial modes i − 301 and i − 401 similarly affect the tesseral Love numbers k 53 and k 64 (respectively).
This stronger tesseral response of inertial modes was absent from the calculations of Lai (2021), who considered the effects of the longest wavelength m = 2 i-modes on k 22 but not k 42 . The effect is interesting in light of the significant deviation between the k 42 inferred from Juno observations, and the theoretical k 42 computed via the CMS method. We find k 42 1.54 at Io's orbital frequency; discrepancies between this and the observed value of k 42 (measured as 1.29 by Durante et al. 2020) could be explained by resonances involving retrograde i-modes with frequencies shifted by, e.g., a stably stratified interior.
The sub-inertial mode spectrum of course becomes more complicated for realistic planetary interior models including discrete regions of stable stratification, as the regular spectrum of pure inertial modes present in isentropic polytropes is replaced by mixed gravito-inertial modes (e.g., Ouazzani et al. 2020). However, even inertial waves confined to a convective envelope maintain a strong tesseral response (Ogilvie 2013). We therefore expect the strong tesseral response of retrograde gravitoinertial modes close to resonance with the tidal forcing frequency to remain relevant in more realistic planetary models with stably stratified interiors.

More rapid rotation
In this subsection we describe our results for more rapid rotation (Ω/Ω d 0.5 − 0.9). While these results are primarily of academic interest when it comes to planetary science, they may bear relevance to very rapidly rotating neutron stars.

Mixing between f-modes and i-modes
As the rotation rate increases, the prograde f-modes f + 22 , f + 33 and f + 44 cross into the sub-inertial frequency range ω < 2Ω at Ω/Ω d 0.4 (as found in MacLauren spheroids by Barker et al. 2016). This leads to mode mixing or "avoided crossings"-mode interactions in which oscillations come close to one another in frequency and trade physical character-between f-modes and i-modes. Figure 6 illustrates two such avoided crossings, between f + 22 and both i + 211 (at Ω/Ω d 0.51), and i + 210 (at Ω/Ω d 0.81). The left panels show the frequencies ω (top) and overlap coefficient amplitudes |Q 22 | (bottom) as a function of Ω for all three modes. Meanwhile, the cross-sections on the right compare the gravitational perturbations for f + 22 (top) and i + 210 (bottom), for Ω/Ω d 0.30 (left) and for Ω/Ω d 0.81 (right).
As shown in the top left panel of Figure 6, the avoided crossing between f + 22 and i + 211 takes place over a narrow range of frequency. Consequently, the eigenfunction overlap associated with even the small frequency separation near Ω/Ω d 0.51 is insufficient to push the value of |Q 22 | for i + 211 beyond 0.1. In contrast, the mixing between f + 22 and i + 210 involves strong frequency repulsion, such that the two oscillations are physically very similar for rotation rates Ω/Ω d ∼ 0.7 − 0.9. The interaction imbues this particular prograde inertial mode with a strong coupling to the = m = 2 tidal potential, so that for these rapid rotation rates there are effectively two sectoral, prograde f-modes. As a result, it is important to include i + 210 in Love number computations for Ω/Ω d 0.7.
We note that  identified an avoided crossing between i + 210 and another i-mode in their calculations, which excluded both centrifugal distortion and gravitational perturbations. We similarly find that i-modes with different wavelengths can undergo avoided crossings, but these are generally less important to the tidal response than the mixing between i-modes and f-modes demonstrated in Figure 6. Figure 7 shows our Love number calculations at multiple rotation rates for an n = 1 polytrope, as indicated by color schemes that extend from Ω/Ω d 0.01 (dark) to Ω/Ω d 0.90 (light). The left panel set shows k m vs. Ω o /Ω (like Figure 4), while the right panel set shows dynamical corrections to k hs m (like Figure 5). In both the left and right sets of panels, the separate columns show Love numbers computed with (right) and without (left) the contributions from the inertial modes included in this study. Although we have checked that adding additional modes does not significantly alter the results summarized in Figure 7, for the most rapid rotation rates (Ω 0.8Ω d ) centrifugal distortion is strong enough (see Figure 6, right) that our truncated mode expansion should be treated with caution, since additional, shorter-wavelength modes might play a role in the aggregate. Figure 7 demonstrates some qualitative differences between the tidal response with/without the inclusion of i-modes, in particular the longest wavelength i − m01 and i + m10 . First, the strong overlap of i − 201 , i − 301 and i − 401 with the tesseral ( m) = (42), (53), and (64) components of the tide persists at all rotation rates, as seen in the panels showing k 42 , k 53 , and k 64 (resp.). The tesseral coupling of i − 101 with the ( m) = (31) tidal potential also affects k 31 , but the mode frequency is too negative to resonate with the tidal forcing frequencies considered. As in the case of Ω/Ω d 0.3, the shorter wavelength, retrograde i-modes with n 1 + n 2 = 2 contribute much narrower resonances to all the other Love numbers. and Ω/Ω d 0.81 (right). At moderate rotation, the eigenfunction for i + 210 is qualitatively distinct from that of f + 22 , and its gravitational perturbation is small under the normalization ξ, ξ = 1. During the mixing around Ω/Ω d 0.7, however, the two oscillations lose distinction.

Love numbers for very rapid rotators
Secondly, the avoided crossing between the prograde modes f + 22 and i + 210 (see Figure 6) results in significant differences between the k 22 and k 42 values calculated with and without the inclusion of the latter mode, even far from resonance (see the ∆k 22 and ∆k 42 panels). This considerably alters the zero-frequency response (k hs m ) with/without i-modes. For very rapid rotation (Ω/Ω d 0.87), mixing with the m = 2 f-mode causes even the shorter-wavelength, prograde i-mode i + 220 to have a larger impact on the hydrostatic Love numbers; we find that excluding it leads to a decrease in k hs 22 and k hs 42 . At all rotation rates, forcing of the sectoral f-modes by the = m component of the tide, along with nonzero cross-correlations Q mm Q m+2,m due to centrifugal distortion, leads to a divergence in the tesseral Love numbers (k 42 , k 53 , and k 64 ) as Ω o → 0. However, as discussed in Section 4.3.1, this divergence is not particularly meaningful in that δΦ 42 , δΦ 53 , and δΦ 64 remain finite (for moderate rotation rates), only becoming large relative to U 42 , U 53 , and U 64 (respectively).
On the other hand, at very rapid rotation rates k 33 , k 53 , k 44 and k 64 exhibit more meaningful divergence at orbital frequencies Ω o /Ω > 0, due to resonances: as indicated by Figure 3  These resonances arise because f-modes that are retrograde in the rotating frame (ω < 0) can become prograde in the inertial frame for rotation rates large enough that their inertial-frame frequency is positive (σ = ω + mΩ > 0). We note that the property of propagating backward in the rotating frame, but forward in an inertial frame renders these modes unstable to secular excitation via the Chandrasekhar-Friedman-Schutz (CFS) instability (Chandrasekhar 1970;Friedman & Schutz 1978). The CFS instability involves the enhancement of retrograde modes' negative angular momentum, due to the removal of positive angular momentum by gravitational wave radiation.
The required rotation rates for the onset of the CFS instability for sectoral f-modes are quite large for small m, and decrease for larger m. Ipser & Lindblom (1990) showed that for an n = 1 polytrope, all the retrograde sectoral f-modes save for f − 22 can experience CFS instabiity before the mass-shedding (break-up) rotation limit. In any case, the required rotaton rates are unlikely to be . Love numbers (left panels) and dynamical Love number corrections (right panels), both with and without the inclusion of inertial modes (with n1 + n2 ≤ 2). The prograde inertial mode i + 210 contributes more strongly to the hydrostatic Love number k hs 22 = k22(Ωo/Ω = 1) and k hs 42 for rotation rates Ω 0.7, as it mixes with f + 22 (see Figure 6). In all but k31, retrograde i-modes produce resonant features that are wider for the tesseral Love numbers with > m. The tesseral Love numbers additionally diverge as Ωo → 0, due to contributions from sectoral f-modes driven by the sectoral tidal potential (see Section 4.3.1). On the other hand, for rapid rotation rates k33, k53, k44 and k64 exhibit more meaningful resonances at orbital frequencies Ωo > 0, as the retrograde sectoral f-modes f − 33 and f − 44 pass into the CFS-unstable regime (this occurs when σ = ω + mΩ = mΩo). realized in giant planets. Such rapid rotation is not excluded for neutron stars, however; as found by Ho & Lai (1999), resonant f-mode excitation in inspiralling neutron star binaries may affect the observed gravitational waveform.

SUMMARY AND CONCLUSIONS
We have developed a new, precise method for computing the dynamical, non-dissipative tidal response of rapidly rotating planets and stars, modeled in this paper by isentropic polytropes. The method involves summing over contributions to the potential Love numbers k m from normal mode oscillations, which we have computed using a non-perturbative treatment of rotation that fully accounts for the effects of both the Coriolis force and centrifugal distortion.

Moderate rotation rates
Focusing on the polytropic index n = 1 most relevant to giant planets and neutron stars, we have evaluated the relative importance of fundamental modes (fmodes), acoustic modes (p-modes), and inertial modes (i-modes), as dictated by their frequencies ω and tidal overlap coefficients Q m (see Figure 2). Identifying the most significant modes for rotation rates Ω 0.4Ω d , where Ω d = (GM/R 3 eq ) 1/2 is the breakup rotation rate, we have provided analytical expressions (polynomial fits; see Figure 3; Table 1; Table 2) to the quantities required to compute a variety of Love numbers.
For most rotation rates in this range, a perturbative treatment of rotation (to linear order in Ω) focused on f-modes is adequate for capturing the sectoral ( = |m|) dynamical tidal Love numbers, to which p-modes and i-modes contribute negligibly (Lai 2021). Such a perturbative calculation can, for instance, satisfactorily explain the discrepancy between the value of Jupiter's k 22 inferred from Juno observations, and predicted by calculations ignoring dynamical effects. However, the tesseral ( > |m|) coupling of retrograde i-modes to the tidal potential is more important to consider than their sectoral coupling; we show that the "leading" m = 2 inertial modes can induce broad resonant features in (e.g.) k 42 , significantly modifying the dynamical corrections to the hydrostatic value (see Figures 4-5).
Additionally, we show that sectoral f-modes driven by the sectoral components of the tidal potential (U m with = m) can produce a finite tesseral gravitational response (δΦ m , with > m)-an order O(Ω 2 ) effect that is not captured by perturbative treatments of rotation that ignore centrifugal distortion of the planet or star. Driving of the tesseral tidal response by sectoral components of the tidal potential results in tesseral Love numbers that diverge with increasing companion separation (see Section 4.3.1). This divergence suggests that alternative response functions (e.g., potential Love numbers relating unlike components of the tidal potential and the response) may be more useful for characterizing the response of significantly centrifugally distorted bodies.

Rapid rotation rates
We have also considered the tidal response of much more rapid rotators. For Ω/Ω d 0.7, mixing between prograde, sectoral f-modes and prograde i-modes can enhance the overlap of the latter with the sectoral components of the tidal potential ( Figure 6; panels showing k 22 in Figure 7). These "avoided crossings" make the contributions of prograde i-modes important to consider when characterizing the tidal response of rapidly rotating planets or stars.
Lastly, the criterion of negative (positive) frequencies in the rotating frame (inertial frame) for secular CFS instability also demarcates a regime in which retrograde, sectoral f-modes with m > 2 may come into resonance with tidal forcing frequencies associated with distant, spin-orbit aligned companions (bottom four panels of Figure 7). Such resonant f-modes may significantly alter the dynamical corrections to both the sectoral and tesseral Love numbers k mm and k m+2,m in very rapidly rotating, centrifugally distorted bodies.
We thank the anonymous referee who reviewed this manuscript, and provided thorough and constructive comments that significantly improved the quality of the paper. J. W. D. gratefully acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) [funding reference #CITA 490888-16], and from the Sloan Foundation through grant FG-2018-10515. R pol /R eq M/( c R 3 eq ) /(4 ) 10 × T/|W| Figure 8. Left: relative changes in Φ(r, µ) (blue) and virial errors V (orange) with successive iterations for calculations of a near-maximally rotating, n = 1 polytrope with Ω/Ω d 0.99 (note that less rapidly rotating models converge with fewer iterations). Right: variation of model parameters for the n = 1 polytropic models considered in this paper.
We solve this equation in the radial domain r/R eq ∈ [0, 2], using Chebyshev collocation with two radial domains split between r/R eq ∈ (0, 1] and [1,2]. We enforce the regularity of each Φ at the origin, and the continuity of each Φ and its radial derivative at r/R eq = 1. At the outer boundary r/R eq = 2, we impose the boundary condition [∂ r +( +1)/r]Φ = 0. We then reconstruct the 2D field Φ(r, µ) from the coefficient Φ (r), before calculating subsequent values of H c , ρ, and λ from the algebraic Equations A6-A8 [setting ρ = 0 exterior to the surface r s (µ) defined by H = 0]. Figure 8 (left) illustrates the convergence of the iteration scheme for a model with polytropic index n = 1 and nearly maximal Ω/Ω d 0.99. The plot shows relative changes in gravitational field with each iteration (blue), and also a virial error defined by V = 1 + (2T + U )/W (orange), where The virial error assesses the degree to which the stellar model satisfies the virial theorem 2T + U + W = 0; we achieve v 10 −6 for all models considered in this paper (with much smaller errors for more centrally condensed polytropes). Figure 8 (right) plots the ratio between polar and equatorial radius (blue), total mass (orange), eigenvalue λ/(4π) = Gρ c R 2 eq /H c (green), and ten times the ratio of kinetic to gravitaional potential energy, 10T /|W |. Accounting for differences in dimensionalization, we find values in agreement with Reese et al. (2006) and Passamonti et al. (2009b).

B. NUMERICAL DETAILS
The mode (model) calculations presented in this work were performed with N ζ = 200 (N r = 400) Chebyshev collocation points ζ = −x ∈ (0, 1], where x = cos[πi/(2N ζ − 1)], i = N ζ , ..., 2N ζ − 1 constitute half of a Gauss-Lobatto (endpoint-extrema) grid covering the stellar/planetary interior. We have halved the grid and constructed spectral derivative matrices from exclusively even or odd Chebyshev polynomials in order to apply regularity boundary conditions at ζ = 0 implicitly (see Chapter 8 in Boyd 2001), but we note that we recover the same results using a full Gauss-Lobatto grid and explicit boundary condition enforcement via boundary bordering. In the exterior vacuum, our mode calculations used a (full) Gauss-Lobatto grid with N ζ = 31, and spherical harmonic expansions varying from N = 8 to N = 50 (depending on rotation rate). Both frequencies and eigenfunctions were validated against calculations at lower resolution (in both N ζ and N ).  Figure 2, providing a more complete breakdown of tidal coupling coefficients for different f-modes (top), and i-modes (bottom) calculated for an n = 1, Ω/Ω d 0.3 polytrope. For each mode labelled on the x-axis, the colormaps show Q m for different spherical harmonic degrees (and the single m associated with each mode). The color-plots for f-modes (top), and i-modes (bottom) are all saturated identically, and transition from log to linear scale near zero.
Fundamental modes produce the largest Q m by far. Inertial modes' Q m values are smaller, but the frequencies of these oscillations generally lie much closer to resonance with the tidal forcing. Inertial modes also have a much stronger tesseral than sectoral overlap (e.g., the m = 2 i-modes i − 201 and i + 210 have much larger Q 42 than Q 22 ). Additionally, with non-zero rotation the f-modes f ± m invariably gain overlap coefficients Q +2,m with the opposite sign of Q m . As described in Section 4.3.1, this sign difference produces the cross-correlations Q +2,m Q m < 0 that are responsible for the divergence of tesseral Love numbers with increasing satellite separation (decreasing orbital frequency).