A Spectral Method to Compute the Tides of Laterally Heterogeneous Bodies

Body tides reveal information about planetary interiors and affect their evolution. Most models to compute body tides rely on the assumption of a spherically symmetric interior. However, several processes can lead to lateral variations of interior properties. We present a new spectral method to compute the tidal response of laterally heterogeneous bodies. Compared to previous spectral methods, our approach is not limited to small-amplitude lateral variations; compared to finite element codes, this approach is more computationally efficient. While the tidal response of a spherically symmetric body has the same wavelength as the tidal force; lateral heterogeneities produce an additional tidal response with a spectra that depends on the spatial pattern of such variations. For Mercury, the Moon, and Io, the amplitude of this signal is as high as 1%–10% of the main tidal response for long-wavelength shear modulus variations higher than ∼10% of the mean shear modulus. For Europa, Ganymede, and Enceladus, shell-thickness variations of 50% of the mean shell thickness can cause an additional signal of ∼1% and ∼10% for the Jovian moons and Encelaudus, respectively. Future missions, such as BepiColombo and JUICE, might measure these signals. Lateral variations of viscosity affect the distribution of tidal heating. This can drive the thermal evolution of tidally active bodies and affect the distribution of active regions.


INTRODUCTION
The deformation of planets and moons under a tidal force depends on their internal properties -e.g., density stratification, elastic and anelastic properties.As a consequence, observations of body tides can be used to infer the internal properties of planetary objects (e.g., Konopliv & Yoder 1996;Wahr et al. 2006;Steinbrüugge et al. 2018;Williams & Boggs 2015;Sohl et al. 2014;Hoppa et al. 1999;Segatz et al. 1988;Hamilton et al. 2013).
Traditionally, body tides have been studied under the assumption of spherical symmetry, i.e., the internal properties of the body do only vary radially and not laterally.However, planets and moons are not spherically symmetric.Seismic tomography evidences that Earth's lithosphere and mantle are not spherically symmetric but present lateral structures associated with tectonics and deep mantle dynamics (Woodhouse & Dziewonski 1984;Ritsema et al. 2011); the Moon shows a farside-nearside the spectral-finite-element method of Martinec (2000); Tanaka et al. (2021), angular integrals are performed analytically using the properties of tensor-spherical harmonics and not numerically; and the equations are solved in the Fourier domain, rather than in the time domain (as done for lateral variations of viscosity in Martinec (2000)) or iteratively (as done for lateral variations of elastic properties in Tanaka et al. (2021)).
The paper is structured as follows.Section 2 introduces the viscoelastic-gravitational problem, in Sections 3 and 4 the governing equations are transformed to the spectral domain and solved.In Section 5 we compare the new spectral method with a spectral-perturbation method and a FEM code; and in Sections 6 and 7 we obtain the tidal response of elastic and viscoelatic bodies to illustrate with various types of lateral heterogeneities.The paper finishes with a summary of the main results in Section 8.
Alongside with this paper, we release the LOV3D software package, which uses the method described here to compute the tidal response of laterally-heterogeneous bodies (Rovira Navarro & Matsuyama 2024).

Governing Equations
The tidal response of a self-gravitating body can be obtained by solving the mass, momentum and Poisson's equation.The equations are linearised around a state of hydrostatic equilibrium with gravitational potential ϕ 0 and pressure p 0 given by ∇p 0 = −ρ 0 ∇ϕ 0 and ∇ 2 ϕ 0 = 4πGρ 0 , where G is the universal gravitational constant and ρ 0 the density of the unperturbed body.The resulting equations are then (e.g., Sabadini et al. 2016): ∇ • σ ′ − ρ 0 ∇(gu • e r ) + gρ 0 χe r − ρ 0 ∇ϕ ′ = 0, (1b) u is the displacement vector.σ ′ is the incremental material stress tensor, and ρ ′ and ϕ ′ are the incremental local density and gravitational potential, respectively -with the latest including both the tidal potential and the potential arising from the deformation of the body.χ = ∇ • u is the divergence of the displacement vector, e r is the radial vector and g the gravity of the unperturbed body.To close the system, a constitutive equation relating stress and displacements is required.For an elastic body, we have where † indicates transpose and I is the identity matrix.λ is the second Lame parameter (λ = κ− 2 3 µ), and µ and κ are the shear and bulk modulus.If the body is viscoelastic, the correspondence principle can be used to relate stress and displacement (e.g., Peltier 1974;Sabadini et al. 2016).Eq. (2) still holds but now µ and λ are the Fourier-transformed shear modulus and Lame parameter μ, λ.For Maxwellian rheology, we have: ω is the forcing frequency and τ M is the Maxwell time τ M = η/µ, with η being the viscosity.Different expressions for the Fourier-transformed shear modulus and Lame parameter can be obtained for other rheological models (e.g., Andrade, Voigt-Kelvin, Burgers, Sundberg-Cooper) (Renaud & Henning 2018).
The tidal response is obtained by solving Eqs.(1,2) for a given tidal load, interior structure and boundary conditions.
(4) Y m n are spherical harmonics of degree n and order m (see Appendix A.2), and θ and φ are the co-latitude and longitude, respectively.
As in previous work (e.g., Qin et al. 2016), we account for the effect of thickness variations ∆H of a layer using using effective shear modulus variations: ∆µ/µ 0 = ∆H/H 0 (e.g., Qin et al. 2016).Yet this is an approximation that deserves further comment.The tidal response of a body whose outer layer is a thin solid shell floating above a liquid layer (i.e., thin-shell body) is controlled by the extension and the bending rigidity: 2(1 + ν)µH and µH 3 /6(1 − ν) (Beuthe 2018).The effect of shear modulus and shell thickness variations in the extension rigidity is equivalent.However, this is not true for the bending rigidity.Bending effects are small provided (H/R) 2 ≪ 1 and the response has a characteristic long wavelength; when this is the case, they can be neglected (membrane approximation).As bending effects become more relevant (e.g., thicker outer shells and shorter wavelength), the use of an effective shear modulus to mimic thickness variations becomes less accurate.Beuthe (2018) applied thin shell theory to study the deformation of Enceladus with lateral shell thickness variations and showed that the membrane approximation renders good results provided shell variations are of long wavelength.As bending effects are expected to be smaller in larger icy moons (e.g., Europa and Ganymede), we expect this approximation to hold better for these moons.
While the method described below can be applied to any layered interior structure, we will consider a simplified 3-layer model consisting of a non-deformable solid core of density ρ 1 surrounded by an incomprehensible liquid layer in hydrostatic equilibrium with outer radius R 2 and density ρ 2 ; and a solid shell with outer radius R, density ρ 3 and mean shear modulus, bulk modulus and viscosity µ 0 , κ 0 and η 0 .This simple model can be used to represent rocky and icy bodies.For the former, the innermost two layers correspond to the inner and the outer core and the outermost layer represents both the rocky mantle and crust.For the latter, the core corresponds to the rocky core, the liquid layer to a subsurface ocean and the outer layer to the ice shell.This simplified model neglects the gravitational coupling between the inner solid core and the surrounding shell, which typically have a small effect on the tidal response, as well as dynamic liquid tides in the liquid layer, which under certain circumstances can play a relevant role in the tidal response (e.g., Rovira-Navarro et al. 2023, 2019).
Eq. ( 28 For the moon, we consider a core radius of 390 km and adjust the core and mantle densities to meet the density density and M oI constraint.For Mercury the interior structure follows form the average density, and the total and outer shell moment of inertia (Genova et al. 2019).Io's interior model is from Steinke et al. (2020a).Ice shell thickness for Europa, Ganymede and Enceladus are from Beuthe et al. (2016), Hussmann et al. (2002) and Vance et al. (2018); we assume ice and ocean densities of 1000 kg m −3 and the core density is set to meet the mean density constraint.Io's viscosity is adjusted to match the ℑ(k 2 ) from astrometric observations (Lainey et al. 2009).Lateral variations are given in terms of maximum peak-topeak variations with respect to the mean value.For Io, representative viscosity and shear modulus variations follow from expected variations in tidal heating (Steinke 2021); for Enceladus and Europa expected shell thickness variations are from (Beuthe et al. 2016) and Nimmo et al. (2007), respectively.The gravitational and radial displacement Love numbers (k u 2 and h u 2 ) for the spherically symmetric elastic case is also given.
The tidal response is controlled by a set of non-dimensional parameters.The tidal response of an incompressible homogeneous elastic body only depends on the effective rigidity µ eff .A body with high effective rigidity exhibits small tidal deformations.The role of compressibility depends on the Poisson's ratio ν.For ν ≈ 0.5, the body behaves as an incompressible body.Viscous effects depend on the non-dimensional Maxwell time τ .Bodies with τ → ∞ behave elastically.The introduction of a solid core and a liquid layer adds three non-dimensional parameters: the ratio between the outer radius of the liquid layer and surface radius, r R ; the ratio between the outer layer density and the mean density, r ρ ; and the density contrast between the liquid and solid shells ,r ∆ρ .For icy moons with subsurface oceans r ∆ρ ≈ 0, while for rocky worlds we will assume that the inner and outer core have the same density.Lateral rheology variations introduce additional non-dimensional parameters.The effect of lateral variations of rheology properties depend on the parameters κ m n , µ m n and η m n .Characteristic non-dimensional parameters for tidally-active Solar System bodies are listed in Table 1.

SPECTRAL METHOD
Different approaches can be used to solve Equation (1).We use a spectral method that employs tensor spherical harmonics as basis functions (James & Cook 1976).The properties of spherical harmonics allows to transform the three-dimensional governing equations to a set of partial differential equations that depends only on radial distance.In this section, we describe the new spectral method.First, we explain how the different fields are expanded using tensor spherical harmonics (Section 3.1); we then provide the governing equations and boundary conditions in the spectral domain (Sections 3.2 and 3.3); finally, we extend the definition of Love numbers traditionally used for sphericallysymmetric bodies to laterally-heterogeneous bodies and obtain expressions to compute tidal energy dissipation for the anelastic case (Sections 3.4 and 3.5).The appendices contain further details about the method.Appendix A provides definitions of tensor spherical harmonics; Appendix B gives explicit expressions to compute various integrals of products of tensor spherical harmonics required for this approach; and Appendix C provides explicit expressions for the governing equations.

Spectral Expansion
As we are considering a periodic forcing, it is convenient to work in the Fourier domain.The tidal potential, ϕ T , can be written as n k , m k and ω k are respectively the degree, order and frequency of the tidal force.ϕ T should be real, implying that for each , where indicates complex-conjugate.The amplitude of ϕ (T ) m k n k ,ω k depends on the tidal component considered -Appendix D provides the expression corresponding to a synchronous satellite in an eccentric orbit.As Eq. ( 1) is linear, the complete tidal response can be obtained by linearly adding the solution for each tidal component.In what follows, we will consider the solution to a tidal force of frequency ω of the form The + and − components correspond to the +ω and −ω components of the tidal force.We expand rank 0 tensors, such as the perturbing potential ϕ ′ and the divergence of the displacement vector ∇ • u, using spherical harmonics of rank 0 where summation is over n The rheology parameters are also expanded in spherical harmonics (Eq.4).For an anelastic body, the complex shear modulus is obtained using the correspondence principle (Eq.( 3)) with τ M,0 = η 0 /µ 0 the average Maxwell time.The Fourier-transformed shear modulus can be decomposed into spherical harmonics ℜ and ℑ denote the real and imaginary components, respectively.For an elastic body we have μ0 = µ 0 and μm n = (−1) m μ−m n .The displacement vector is expanded using tensor spherical harmonics of rank 1 (Appendix A.3) where summation is over n The decomposition in tensor spherical harmonics automatically separates between spheroidal (n 1 = n±1) and toroidal components n 1 = n.This expansion is different to that traditionally used to solve the viscoelastic-gravitational problem, which uses the radial-poloidal-toroidal decomposition where c.c means complex conjugate.U m n , V m n , W m n and u m n,n−1 , u m n,n , u m n,n+1 are related via a linear transformation (Eq.A14).
(14) This projection has the advantage of separating the trace l = 0, anti-symmetric l = 1 and symmetric trace-free l = 2 components.For a symmetric rank-2 tensor there are 6 non-zero components of the stress and strain tensor for each degree and order: the trace (l = 0, n 2 = 0) and the symmetric spheroidal (l = 2, n 2 = n and n 2 = n ± 2) and toroidal (l = 2, n 2 = n ± 1) components.The Zerilli and rank-2 spherical harmonics tensor spherical harmonics are linearly related (see Appendix A.4).
To apply the boundary conditions, the radial components of the stress tensor are required.The radial component of the stress tensor can be written in the rank-1 spherical harmonics basis or alternatively in the radial-poloidal-toroidal basis often used in literature, The spheroidal, poloidal and toroidal components can be written in terms of the σ (l) n,n 2 ;m components of the stress tensor (Eq.(C34)).
The previous tensors must be real, which implies that the + and − components are related.Once the solution to the + component is obtained the solution to the − component follows immediately from the complex-conjugate properties of tensor spherical harmonics (see Appendix A.1)

Spectral Form of the Governing Equations
Using the spectral expansions defined in Section 3.1 and the properties of tensor spherical harmonics (James & Cook 1976), we can transform the mass and momentum conservation equations and Poisson equation (Eq. 1) as well as the constitutive law (Eq.2) to the spectral domain.
The Poisson equation (Eq.( 1c)) is where ∂ r stands for the derivative with respect to the radial distance r, and D n is an operator defined in Eq. (A20).The right-hand side is 0 if the density is uniform within the layer and the solid is incompressible.The explicit forms of the momentum and the Poisson equation are given in Appendix C.
To obtain the spectral form of the constitutive law (Eq.( 2)), it is useful to rewrite it as T and S stand for the trace and trace-free symmetric parts of the tensor and we have used the definition of the strain tensor An explicit expression for the strain tensor, which can be obtained using the properties listed in §3 of James & Cook (1976), is given in Appendix C. Plugging Eq. (4) into Eq.( 18) and using the spectral expansion of the strain tensor, we can obtain the stress tensor Â0 and A 0 Âm n stand for the laterally-uniform and laterally-varying components of A, respectively.Finally, the components of the stress tensor σ (l) n,n 2 ;m are obtained by projecting the stress tensor into a basis of Zerilli's tensors where : is the inner-product.Using Eq. ( 20), we find are the coupling coefficients, which indicate the coupling of mode (n α , m α ) and (n, m) due to rheology variations of degree and order (n β , m β ).The coupling coefficients with non-zero value determine which modes are coupled.Their explicit form and properties are given in Appendix B.

Boundary Conditions
To obtain the tidal response, the previous set of equations must be solved under appropriate boundary conditions.We impose boundary conditions at the surface R and at the core-mantle boundary R 2 .The derivation of the boundary conditions can be found in Sabadini et al. (2016).At the surface, R, the stress vanishes and the potential stress, Q,is given by with n T and m T the degree and order of the considered tidal forcing.At the boundary between the liquid layer and the outermost solid layer (R 2 ), there is zero tangential and toroidal stresses and the radial stress is given by the difference between the radial displacement and the geoid The gradient of the gravitational potential is Degree 0 and 1 require a special treatment (Farrell 1972;Qin et al. 2014).For degree 0 the displacement and radial stress vector only have radial components.Because of compressibility, the body can experience non-zero radial displacements.However, as the mass of the body does not change, the perturbing potential should be 0, ϕ 0 = 0.The degree 1 solution automatically satisfies gQ m 1 4πG = 0, which implies that two of the boundary conditions are a linear combination of each other.This is because the degree 1 solution includes a translation of the center of mass that does not introduce stress.As we are working in the center of mass reference center we constrain this translation to be 0, . This is because the toroidal mode contains a net-rotation that does not introduce stress.Instead of using T m n (R 2 ) = 0, we impose a 0 toroidal displacement at R 2 without introducing extra stress W m 1 (R 2 ) = 0, as done in Qin et al. (2014).

Love Numbers of an Aspherical Body
The tidal response of a body can be expressed in terms of Love numbers, a set of dimensionless proportionality constants that relate the tidal force and the body's response.For a sphericallysymmetric body, a tidal forcing of a given degree and order results in a response of the same degree and order.Moreover, due to the spherical-symmetry, the Love numbers are independent of the order of the forcing.Because of this, there exist a set of frequency-dependent Love numbers ( radial h n , poloidal l n , and gravitational k n ) per degree that can be evaluated at any radial point.If there are lateral-heterogenities, a forcing of a given wavelength might excite a mode with different wavelength.Because of this the definition of Love numbers must be extended.More generally, we can write the tidal response as n β e iωt + c.c.
(28) α and β indicate the tidal response at degree and order n β , m β due to a forcing at degree and order n α , m α .t are the toroidal Love number, which are 0 for a spherically-symmetric body.The only non-zero numbers Love numbers for a spherically-symmetric body are h nα,mα nα,mα , l nα,mα nα,mα , k nα,mα nα,mα , and as they are independent of order, can be simply written as h n , l n , and k n .

Energy Dissipation
If the body is anelastic, its deformation is not adiabatic.The volumetric energy dissipation is given by ė The previous integral can be computed either in the temporal-spatial domain using ϵ(t, r, θ, φ) and σ(t, r, θ, φ) or alternatively in the spectral domain using their spectral expansions.Plugging Eq. ( 14) into Eq.( 29) and performing the time integral, we find For a given radius, we can project the previous expression into spherical harmonics to obtain the spectra of energy dissipation, i.e., is an integral defined in Appendix B. The average energy dissipation for radius r is given by ė0 0 (r) and can be obtained using the orthogonality and complex conjugate properties of tensor spherical harmonics (Appendix A.1) The total energy dissipation can be found by radially integrating ė0 0 , Alternatively, the total tidal dissipation can be obtained from the work done by the tidal force (Love 1927) which can be transformed into a surface integral (e.g., Zschau 1978;Peale & Cassen 1978;Platzman 1984) where ϕ ′ is the perturbing potential (i.e, ϕ = ϕ ′ + ϕ T ), and the subscript n indicates that only the component of degree n of the potential is considered.Using the spectral expansion of the tidal potential (Eq. ( 5)), the definition of the Love numbers (Eq.28) and the orthogonality of spherical harmonics, we find: If we consider a spherically-symmetric body in an eccentric orbit, the previous expression reduces to the classic expression (e.g., Peale et al. 1979): 4. NUMERICAL APPROACH The system of differential equations of Section 3.2 is written in terms of the 3 components of the displacement vector (u m n,n−1 , u m n,n , u m n,n+1 ), the 6 components of the strain (ϵ n,n−2;m , ϵ (2) n,n+2;m ) and stress (σ n,n+2;m ) tensors, and the gravitational potential (ϕ m n ) for each degree and order.However, it is more convenient to rewrite the equations in terms of the variables traditionally employed to solve the viscoelastic deformation of a self-gravitating body (e.g., Love 1911;Farrell 1972;Sabadini et al. 2016): the radial, poloidal and toroidal components of the displacement vector and of the radial component of the stress tensor (U m n , V m n , W m n , R m n , S m n , T m n ) and the gravitational potential and its gradient (ϕ, ∂ r ϕ): where we have also introduced the y radial functions traditionally used in the viscoealstic-gravitational problem (e.g., Sabadini et al. 2016).We note that the problem is formulated in terms of the gradient of the gravitational potential y 6 = ∂ r ϕ rather than the potential stress Q.The equations can be then be cast in the form (Appendix C) where D is a linear operator that depends on the interior properties and the radial distance r.Given a forcing of degree n T and order m T not all the modes (n T , m T ) are excited.If there are no lateral rheology variations, equations of different degree and order are decoupled.. Additionally, if we only consider a zonal forcing m T = 0 and no longitudinal rheology variations (i.e., ) and toroidal (W m n , T m n ) modes are decoupled.More generally, the modes (n, m) involved in the tidal response depend on the forcing spectra (n T , m T ) and that of the lateral variations (n LV , m LV ) via the coupling terms.Given a mode (n α , m α ) and lateral variations of the form (n LV , m LV ) the coupled modes (n α , m α ) (n LV , m LV ) ⇒ (n, m) are given by the non-zero coupling coefficients.The coupled modes can be obtained recursively by using the selection rules listed in Appendix B. Starting from the tidal force (n T , m T ) the modes involved in the solution can be obtained by recursively applying the selection rules: first order modes, (n As opposed to a spherically symmetric body, for which there is only one mode involved in the tidal response per forcing harmonic (n T , m T ), an infinite set of modes are excited by the tidal force when lateral variations are considered.This makes it impossible to obtain an exact solution to the problem using the spectral method.An approximate solution is obtained by setting a maximum cut-off order N p .The reduced system of differential equations is given by with y a containing only the considered modes and D a their corresponding dynamics.y a is a 8N modes vector or 8(N modes − 1) + 4 vector if the degree 0 response is excited.To obtain the tidal response, Eq. ( 41) is integrated numerically from the boundary between the liquid layer and the outermost solid layer (R 2 ) to the surface using a Runge-Kutta scheme.The solution is given by C m n,k are a set of 8N modes , 8(N modes − 1) + 4 if the degree 0 is excited, integration constants; and y m n,k are the corresponding set of solution vectors with The integration constants are obtained by applying the boundary conditions listed in Section 3.3 and the solution follows from Eq. ( 42).Once the radial functions y are obtained, the components of the displacement vector, and strain and stress tensors in tensor spherical harmonic basis are obtained using Eqs.(A14), (C26) and (C29).The solution can be transformed from the spectral domain to the spatial and time domain using Eqs.( 7), ( 11) and ( 14) and the definition of the rank 0 (Eq.( A9)), rank 1 (Eq.( A11)) and rank 2 (Eq.( A16)) tensor spherical harmonics.
Figure 1 provides an overview of the method described above, which is implemented in the LOV3D software repository.

MODEL BENCHMARK AND COMPARISON WITH PREVIOUS METHODS
We compare the results obtained with the spectral method presented here with those obtained using the spectral-perturbation method of Qin et al. (2014Qin et al. ( , 2016) ) and the FEM code of Berne et al. (2023a,b).
The spectral method of Qin et al. (2014Qin et al. ( , 2016) ) relies on perturbation theory.As opposed to the method presented above, the solution is computed recursively.The tidal response of the sphericallysymmetric body is used to obtain the first-order modes response, which is then employed to compute second-order modes, etc.For each mode, only the radial functions corresponding to that mode are considered unknown.The coupling terms (second term in the RHS of Eq. ( 23)) are considered known and follow from a lower order mode, effectively acting as a forcing.This way, the equations corresponding to each of the modes participating in the tidal response are decoupled.This means that for a mode of a given perturbation order, the effect of higher order modes is ignored.
The spectral-perturbation method of Qin et al. (2014Qin et al. ( , 2016) ) provides an approximated solution whose accuracy is expected to decrease as the amplitude of lateral variations increases.To test the accuracy of the perturbation method, we consider an elastic Io-like body (Table 1 (2, 0)-and compute the tidal response.We note that similar results are obtained for different interior parameters and forcing.For each of the considered lateral heterogeneities, we obtain the gravitational Love numbers using the two methods and compute the difference between them.The tidal response can be written as the tidal response of a spherically symmetric body, given by k u 2 , plus an additional response arising from lateral variations ∆k n,m 2,0 .Figure 2 shows the tidal Love number spectra for the three sets of lateral heterogeneities; first, second and third order modes are indicated.The additional tidal response increases from being 0.01% of the tidal response of the uniform body for peak-to-peak shear modulus variations of ∼ 0.1% the mean shear modulus to about 10% for variations of the same order of magnitude as the mean shear modulus.The additional tidal response is clearly dominated by first-order modes, which for small Eq. ( 23) Eq. ( 28) Eq. ( 41) Section 3.3 Figure 1.Flowchart of the methodology employed to obtain the tidal response lateral variations are orders of magnitude higher than second-order modes.However, the difference between first and second order modes decreases as the amplitude of lateral variations increases.
The results obtained with the spectral and the spectral-perturbation method show good agreement.Differences between the results obtained with the two models are as small as ∼ 10 −5 % for lateral variations of less than ∼ 10 −3 %.The difference is likely caused by differences in how the equations are radially-integrated; for instance, the difference decreases as the number of radial points used in the Runge-Kutta integration increases.This difference increases as the amplitude of lateral heterogeneities increases.However, it remains small even when substantial lateral heterogenites are considered.For 10% peak-to-peak variations of the shear modulus the difference remains below 1%, discrepancies of 10% or more are only attained when peak-to-peak variations of the shear modulus have the same order of magnitude as the mean shear modulus.The remarkable agreement is due to the fast decrease of mode amplitude with increasing mode order, as shown in Figure 2a,c and d.This means that the effect of modes of higher order modes on modes of lower order remains small even with high-amplitude lateral heterogeneities.
We also consider the case of an icy moon.We obtain the tidal response for the reference Enceladus model (Table 1) for various types of lateral variations using both spectral methods and the FEM code.For the FEM approach, we adapt the methodology outlined in Berne et al. (2023a) and Berne et al. (2023b) to develop a FEM capable of simulating tidal deformation on laterally heterogeneous ocean worlds.FEMs solve the equation of motion for quasi-static problems by formulating 3D displacements (i.e., in response to applied forces and boundary conditions) as a series of linear shape functions across a mesh domain.For this work, we discretize mesh domains using tetrahedra with a maximum edge length of 1km and consider between 18•10 6 and 33•10 6 elements for each simulation.To simulate tidal loading, we upload mesh geometries to a modified version of the geodynamic software package Pylith (Aagaard et al. 2007) which can self-consistently consider forces arising from external tidal potentials, self-gravity, and radial displacements at boundaries between internal density layers (for additional details, see Supplementary S1.1 of Berne et al. (2023a)).To incorporate lateral variations in elastic properties, we sample analytic basis functions (i.e., spherical harmonics) for a given heterogeneity at mesh node locations.Following simulations, we expand displacements into spherical harmonics and compute Love numbers for comparison to semi-analytic solutions presented in this work.
Figure 3 compares the results obtained with the three methods (i.e., spectral method, spectralperturbation method and FEM).With few exceptions, the results obtained with the FEM and the spectral method show the best agreement.As explained above, the difference between the spectral and spectral perturbation method grows with the amplitude of lateral variations reaching values as high as 5% for peak-to-peak shear modulus variations of 50%.In contrast, the difference between results obtained with the spectral and the FEM remains below 1% and is often one order of magnitude smaller than for the former.We ascribe small discrepancies between results with the spectral method and the FEM code to errors associated with resolution.We find that the agreement between both methods improves when the FEM mesh size decreases from 5 to 1 km.
To benchmark the viscoelastic component, we compute the Love numbers of a spherically-symmetric uniform body obtained using analytical expressions (e.g., Matsuyama et al. 2018, Appendix B) and LOV3D and find excellent agreement (relative differences of < 10 −6 %).For a viscoelastic body with lateral variations we check that energy dissipation is computed self-consistently by comparing the total energy dissipation obtained using the two approaches outline in Section 3.5 (i.e., Eqs. ( 34) and ( 37)), finding good agreement (see Appendix E).  2014) (lower panels).An elastic Io (Table 1) with shear modulus variations of different wavelengths (n LV , m LV ) is assumed.The amplitude of lateral variations is given in terms of the peak-to-peak variation of the shear modulus with respect to its mean.The solution is cut-off at perturbation order 3.The line-width indicates the order of the mode, with thicker lines indicating lower order modes.In the upper panels the dashed lines corresponds to the solution obtained with the perturbation method of Qin et al. (2014).
, with R and Q standing for the method presented here and the method of Qin et al. (2014), respectively.

THE TIDES OF A LATERALLY-HETEROGENEOUS ELASTIC BODY
Figure 2 evidences that different spatial patterns of lateral variations results in a distinct tidal response.This is further illustrated in Figure 4, which shows the degree 2 gravitational Love spectra for monochromatic (single total and zonal wave-number) shear modulus variations.Each type of lateral variations leads to a unique Love number spectra.Hence, if the full tidal response was measured the inverse problem could be solved and the spatial pattern and amplitude of lateral variations inferred, this technique is known as tidal tomography.
The use of tidal tomography presents several challenges.The amplitude of the additional tidal response arising from lateral variations is small, making it challenging to measure.Because of this, we can expect that only the terms with highest amplitudes will be measured.Moreover, lateral variations in different regions and of different properties (e.g., elastic properties, layer thickness) can  2 but for the Enceladus model (Table 1).Differences are obtained as 100 ∆k n,m , where X can either be the results obtained using the result in a similar tidal spectra.Finally, lateral variations are likely not monochromatic but feature various wavelengths.This makes the solution of the inverse problem non-unique.As an example, both degree 1 and 3 lateral variations produce a degree 3 tidal signal (see Figure 4).If no other terms of the spectra are measured, a degree 3 response would be indicative of a hemispherical dichotomy but could not be used to distinguish between degree 3 and 1 lateral variations.Distinguishing between the two would require measuring the degree 5 tidal response, which is more prominent for degree 3 lateral variations.Lateral variations also alter the tidal response at the degree and order of the forcing, most prominently for degree 2 variations.This makes it challenging to distinguish the effects of lateral heterogeneity from the tidal response and might even result in errors when estimating the mean properties of the body.In such cases, a way to detect lateral variations is by comparing the k 2 Love numbers at different orders, which are only equal if the body is spherically symmetric.The tidal potential of a moon in an eccentric orbit has both order 0 and 2 components, making this approach attractive.So far, tidal tomography has just been used for Earth (Lau et al. 2017), for which there is high quality geodetic data.However, future space missions might make it possible to use tidal tomography for other Solar System bodies.Figure 5 shows the gravitational Love number spectra for Io, the Moon and Mercury (Table 1).In all cases, the additional tidal response (∆k) is normalized by the tidal response of the spherically-symmetric body (k u 2 ), listed in Table 1.The Gravity Recovery and Interior Laboratory (GRAIL) mission provided very accurate lunar gravity data, which can be used to constrain the amplitude of lateral heterogenities in the Lunar interior (Qin et al. 2014).We caution that a comprehensive inversion can only be done if an expression for the Love numbers that account for the effect of lateral variations (i.e., Eq (28)) rather than the classic expression often used (e.g., Konopliv et al. 2013, Eq. ( 17)) is used when obtaining the gravity solution from raw GRAIL data, and the model parameter space is thoroughly sampled.However, the results shown in Figure 5 combined with the available GRAIL gravity field give a qualitative impression of the amplitude of lateral variations.For instance, Williams et al. (2014) showed that the measured degree 3 response of the Moon is consistent with a spherically-symmetric model, this suggests that there are not high amplitude odd-degree lateral variations, as otherwise one would measure a high degree 3 response.Additionally, the degree 2 orders 0 and 2 Love numbers of the Moon differ by ∼ 1%, but are within the uncertainty of the measurements (Konopliv et al. 2013;Lemoine et al. 2013).From Figure 5, this suggests that degree two variations are less than ∼ 10%.
Of the three bodies, Mercury is the one for which the additional tidal response is the smallest relative to the tidal response of a spherically-symmetric.Nevertheless, as Mercury is the body with the highest k u 2 , lateral variations lead to the highest gravity signal in absolute terms.The MESSENGER mission measured Mercury's k 2 with an accuracy of approximately 5%, insufficient to observe lateral variations.In contrast, BepiColombo, scheduled to start its science operations in 2026, is expected to improve the accuracy to ∼ 0.1% (Genova et al. 2021).The accuracy at which the non-diagonal Love numbers can be obtained (e.g., k n β ,m β nα,mα , with n α ̸ = n β , m α ̸ = m β ) can only be determined via a simulation of BepiColombo's gravity science experiment.Assuming the accuracy inferred for k 2 extends to the other Love numbers, variations of the shear modulus or of the thickness of Mercury's rocky envelope as small as 2% might be detected.As demonstrated by Figure 5, the detection threshold depends on the spatial pattern of lateral variations.
We also consider the role of lateral heterogeneities in icy moons.As noted in Section 2.2, shell thickness variations can be approximately mapped to shear modulus variations provided the shell is thin compared to the moon's radius.This is the case of icy moons -for Europa and Ganymede the shell to radius ratio is ∼ 0.02 and 0.05, respectively; for Enceladus it is higher ∼ 0.1, making the approximation less accurate.Figure 6 shows the additional tidal response for zonal variations of the shear modulus for Europa, Ganymede and Enceladus.The relative additional tidal response resulting from lateral variations (∆k n,m n T ,m T /k u n T ) is approximately one order of magnitude higher for Enceladus than for the two Jovian moons.For Enceladus, lateral variations of 50% cause an additional tidal response ∼ 10% the amplitude of the main tidal response (see also Berne et al. (2023a) and Bêhounková et al. ( 2017)), while for the latter the additional tidal response is ∼ 1%.Nevertheless, as k 2 of Enceladus is more than one order of magnitude smaller (Table 1), ∆k n,m n T ,m T is similar for the three icy moons.
The possibility to observe lateral variations depends on their amplitude and the accuracy to which the Love number spectra can be obtained.For Enceladus, gravity and shape data indicates that shell thickness varies from 29 km at the equator to 7 km at the south pole (Beuthe et al. 2016;Čadek et al. 2016;Hemingway & Mittal 2019).In contrast, shell thickness variations for Europa are not expected to exceed 7 km (Nimmo et al. 2007).Ganymede's shell thickness variations are not constrained, yet ocean circulation models predict shell thickness variations to be smaller for big icy moons (Kang & Jansen 2022).This makes Enceladus a prime candidate for the use of tidal tomography in the future.With the JUpiter Icy Moons Explorer (JUICE ) on its way to the Jovian system, it is interesting to consider if it can pick up the signal arising from lateral variations of shell properties.JUICE will ).An elastic Io is assumed (Table 1), and the solution is cut-off at perturbation order 4. Color intensity indicates the strength of the mode, and the terms with highest amplitude are indicated with a cross.For all cases, we assume peak-to-peak shear modulus variations of 10% of the mean shear modulus.measure Ganymede's k 2 with an expected accuracy of 10 −4 (Cappuccio et al. 2020), which is less than 0.1% of Ganymede's expected k 2 .This accuracy is sufficient to detect differences in k 2,0 2,0 and k 2,2 2,2 caused by lateral variations of ice shell properties, making tidal tomography a promising tool for constraining lateral variations in Ganymede's shell.

THE TIDES OF A LATERALLY-HETEROGENEOUS MAXWELL BODY
In the previous section, we considered the tidal response of an elastic body with lateral variations.We will now consider how lateral variations affect the tidal response of a viscoelastic body.If a body is not perfectly elastic, part of the tidal energy is converted into heat.Even for a sphericallysymmetric body, tidal heating is not homogeneously distributed within the body's interior.Here we explore the effect that lateral variations of the viscosity has on tidal heating patterns.
We employ our reference Io model (see Table 1) and consider various patterns of viscosity variations.The complex shear modulus spectra is obtained using Eqs.( 8)-( 10).Even in the case of monochromatic viscosity variations the complex shear modulus spectra is not monochromatic but contains infinite terms of varying amplitude.To use the spectral method, the complex shear modulus spectra needs to be cut-off.We consider terms up to 2 orders of magnitude smaller than the leading non-degree-0 term of the complex shear modulus.The heating pattern depends on the tidal potential spectra.We focus on the tides experienced by a synchronous body in an eccentric orbit with zero obliquity -for which we give the tidal potential in Appendix D. We obtain the tidal response for each component of the tidal potential and compute the energy spectra as described in Section 3.5.
Figures 7 and 8 show the spatial distribution of tidal heating for various types of monochromatic lateral viscosity variations and its corresponding spectra (Eq.31), respectively.The tidal heating pattern characteristic of a spherically-symmetric body consists of terms of degree and order (0, 0), (2, 0), (2, ±2), (4, 0), (4, ±2) and (4, ±4) (e.g., Beuthe 2013).For the three layer model considered here (solid core, liquid core and rocky envelope), tidal heating is maximum at the poles; along the equator, tidal heating exhibits minima at the sub-planet (0 • ) and anti-planet (180 • ) points and maxima at the center of the trailing and leading hemispheres (90 • and −90 • , respectively).Viscosity lateral variations introduce additional terms.The difference between the tidal heating pattern of a uniform body and a body with lateral variations is dominated by the same pattern as the considered viscosity variations.For example, a polar dichotomy of degree 1 and order 0 in viscosity translates into the same dichotomy in tidal heating (see second panel in Figure 7).This also follows from the tidal heating spectra, where the most prominent terms are those with same wavelength as the considered viscosity variations (Figure 8).Apart from this dominant term, other terms with smaller amplitude also arise.
Some terms of the tidal heating spectra are not symmetric with respect to the 0 • meridian.This is most evident by looking at the tidal heating map corresponding to the degree 2 order 0 lateral variations, where we observe a westward shift in the longitude at which the minimum heat flux is attained compared to the spherically-symmetric case.The appearance of a trailing-leading hemispheres asymmetry is remarkable since the considered patterns of lateral variations do not contain such an asymmetry.A similar phenomenon was observed by Steinke et al. (2021), who used FEM to obtain tidal heating patterns for a multi-layered, laterally-heterogeneous Io.The asymmetry arises Figure 6.Effect of lateral heterogeneities of different wavelengths for Europa (first row), Ganymede (second row) and Enceladus (third row).The highest amplitude component of the tidal response arising from zonal lateral variations of degrees 1 (first column), 2 (second column) and 3 (third column) is depicted.For reference, the degree two Love number of a spherically symmetric body, k u 2 , is indicated (dashed black lines).
because the tidal potential is not symmetric with respect to the east-west direction.As shown in Appendix D, the tidal potential can be broken into a degree 2 order 0 standing wave and degree 2 order 2 westward and eastward propagating waves of different amplitudes.The amplitude of the eastward component is greater than that of the westward component, breaking the symmetry of the problem.No asymmetry is observed if a single standing-wave is considered.Tidal heating patterns can help constraining the interior properties of rocky and icy worlds (e.g., Breuer et al. 2022).The distribution of volcanoes in Io has been used as a proxy of tidal heating.In Io, the concentration of volcanoes is greater towards mid-to low-latitudes and bimodal, exhibiting two maxima 30 • −60 • eastward of the sub-Jovian and anti-Jovian points (see Figure 7).Moreover, the distribution of volcanoes also contains a statistically-significant degree-6 component (Hamilton et al. 2013;Kirchoff et al. 2011;Steinke et al. 2020b).The link between tidal heating and volcanic patterns is not simple: convection and melt transport affects how the two patterns relate (Steinke et al. 2020a), and the inferred volcanic patterns are affected by biases in observations.Recent observations by JUNO have discovered new volcanoes in previously poorly-covered polar regions and even hinted that the north pole could have a higher concentration of hot-spots (Zambon et al. 2023;Davies et al. 2024).
The tidal heating pattern characteristic of a 2-layer (core+rocky envelope) spherically-symmetric Io does not match the observed distribution of volcanoes.As illustrated in Figure 7, the surface heat flux peaks at the poles, does not exhibit a 30 • − 60 • eastward shift with respect to the tidal axis, does not have degree-6 terms nor a polar dichotomy.The heating pattern is affected by both radial and lateral variations of internal properties.The presence of a low-viscosity asthenosphere can account for the concentration of volcanoes at mid-to low-latitudes.However, for spherically-symmetric cannot reproduce the observed eastward shift of volcanic activity with respect to the sub-Jovian point.Tyler et al. (2015) showed that tidal heating in a magma ocean can cause this shift.The results discussed here, as well as the FEM model of Steinke et al. (2020a), indicate that lateral viscosity variations can introduce a longitudinal shift.The magnitude and direction of the longitudinal shift depends on the particular form of lateral and radial variations of internal properties.A comprehensive exploration of the joint effect of radial and lateral variations is left for future work.As evidenced in Figure 8, lateral variations of internal properties can also introduce a degree 6 heat flux pattern and cause a polar-dichotomy (terms with odd m).
Interior properties and heating patterns are tightly related.Viscosity and shear modulus depend on temperature and melt fraction.Lateral variations of tidal heating affect both quantities (Steinke et al. 2020a) and hence the heating pattern.This feedback can have important implications for the interior evolution of tidally-active bodies.As shown in Figure 8, even a body with a uniform interior exhibits non-uniform tidal heating.Such a pattern will alter material properties and feedback into the heating pattern, driving interior evolution.Alternatively, primordial lateral heterogeneities might be amplified by the feedback in a similar way.The extent to which this happens depends on how heat is transported within the body.Tackling this complex problem requires coupling a tidal and a thermal model.The lower computational cost of spectral method compared to FEM methods makes them an attractive tool to approach the problem.

SUMMARY AND OUTLOOK
In this work, we presented a spectral method to compute the tidal response of bodies with lateral heterogeneities and applied it to elastic and viscoelastic bodies.Below, we summarize the key points: • The new spectral method shows good agreement with results obtained with the spectralperturbation approach of Qin et al. (2014) and the FEM model of Berne et al. (2023a,b).Differences between the results obtained with the new method and the perturbation method grow as the amplitude of lateral variations increase.Predictions obtained using the perturbation method differ by less than 1% for shear modulus variations of ≲ 10% and only reach ∼ 10% for peak-to-peak variations of the same order of magnitude as the mean shear modulus.This makes the perturbation method a powerful tool to compute the tidal response of bodies unless For the spherically-symmetric body the total heating pattern is shown while for the rest the difference between their corresponding heating pattern and that of a spherically-symmetric body is shown.All plots are normalized by the average tidal heating of the spherically-symmetric body.The reference Io model is employed (see Table 1) with peak to peak viscosity variations of 50%.The sub-planet (SP) and anti-planet (AP) points as well as the centers of the trailing and leading hemispheres (TH and LH) are indicated.
The longitudinally-and longitudinally-averaged tidal heating are also shown.In the first panel, Io's volcanoes are indicated (blue dots), with contours corresponding to the density of Io's volcanoes (red and blue contours indicating higher and lower than the global average) (Steinke et al. 2020b).such variations have high amplitude.Differences between the FEM approach and the presented spectral method remain small even for high amplitude lateral variations.Spectral methods are more computationally-efficient: a model that takes ∼ 10 days to run using the FEM code in a two-core computer is solved in less than a minute using the spectral method presented here; in terms of memory, the grids employed in the FEM runs presented here required ∼ 1 GB of memory, while the coupling matrices employed in the spectral model required ∼ 1 MB.This makes spectral methods more suitable to tackle the inverse problem.However, unlike FEMs, the spectral code presented here cannot include faults and cracks, non-linear rheology, and layer thickness variations can only be treated by mapping them to shear modulus variations (Berne et al. 2023a).
• Each set of lateral variations results in distinct Love number spectra.Thus, measurements of the complete tidal response of a body can be used to constrain lateral variations.Nevertheless, the solution of the inverse problem is non-unique as it is unlikely that the complete Love number spectra can be measured and lateral variations at different depths might result in similar spectra (Section 6).Other measurements, such as static gravity and topography, might help to solve this degeneracy.
• For the Moon, Io and Mercury, shear modulus variations of the same order of magnitude as the mean shear modulus can cause an additional tidal response as high as ∼ 1 − 10% the main tidal response.BepiColombo should observe the fingerprint of lateral variations provided they are higher than approximately 2% the mean shear modulus (Section 6).
• In icy worlds, lateral variations of ice shell thickness modify the tidal response.Due to the expected amplitude of shell-thickness variations, Enceladus is a prime candidate to use tidal tomography.For Europa and Ganymede, shell-thickness variations can lead to an additional tidal response in the order of 0.1 − 1% the main tidal response depending on their amplitude.The accuracy of JUICE makes it possible to detect this signal (Section 6).
• Lateral variations modify the distribution of tidal heating for viscoelastic bodies.The additional tidal heating pattern due to lateral variations of the viscosity is dominated by the pattern of such variations.Lateral viscosity variations can cause a trailing-leading hemisphere asymmetry in tidal heating.This could explain the eastward shift of volcanic activity with respect to the sub-Jovian and anti-Jovian points observed in Io (Section 7).
• The dependence of interior properties on temperature and melt fraction, which in turn depend on the distribution of tidal heating, gives rise to a complex feedback that can drive interior evolution.The computational efficiency of the spectral method makes it a good candidate to study this feedback (Section 7).
• For this work, we considered a simplified interior structure consisting of a solid non-deformable core, overlaid by an hydrostatic liquid and a solid envelope with radially-uniform properties (Section 2.2).These assumptions can be relaxed -i.e., the deformation of the inner core can be considered, dynamic liquid tides included and radial variations of interior properties introduced.For problems with spherical geometry, expanding the tensors in tensor spherical harmonics allows to employ the properties of tensor spherical harmonics to eliminate longitude and latitude from the equations of motion.A rank k tensor a can be written as (James & Cook 1976) Below, we define tensor spherical harmonics and introduce some of their main properties.We also give explicit expressions of tensor spherical harmonics of ranks 0 , 1 and 2; explain how to transform between them and other basis often used in literature and in this manuscript; and list operators used in the text.An exhaustive list of properties can be found in James & Cook (1976).

TENSOR SPHERICAL HARMONICS INTEGRALS
To obtain both the coupling coefficients and the energy dissipation spectra the following integral should be evaluated T (lα) nα,n 2α ;mα : T Eq. ( B21a) is required to compute the coupling coefficients and (B21b) to compute the energy dissipation integrals.The two expressions can be obtained from ) using the complex conjugate properties of tensor spherical harmonics (Eq.( A7)).To evaluate T , we write it in terms of rank-2 spherical harmonics Y m n(2) (Eq.A17): which can be obtained using the expressions provided in James & Cook (1976) §4 We have used the abbreviation Λ(n(k)) = (2n + 1)(2n 1 + 1) . . .(2n k + 1) The amount of integrals that need to be evaluated is greatly reduced by noting that many are 0. T 3. n 2α , n 2β , and n ν satisfy the triangular inequality The previous set of properties can be used to deduce a set of selection rules for the excited modes.A mode of degree and orders (n 0 , m 0 ) together with lateral variations of degree and order (n 1 , m 1 ) results in modes of degree (n 2 , m 2 ) if: where x or y indicate if the mode is spheroidal or toroidal and the cases x = s, y = t, and x = t, y = s can be considered.

C. EXPLICIT FORM OF THE EQUATIONS
The components of the strain tensor follow from the gradient of the displacement vector (Eq.( 19)) χ m n is the divergence of the displacement vector Eq. ( C26) can be written more concisely as n,n 2 ,m being a vector containing the elements of the strain tensor, ϵ n,n+2,m . ..) † .The components of the stress tensor can be written in terms of the strain tensor and the material properties , or in matrix form σ with σ (l) n,n 2 ,m being a vector containing the components of the stress tensor, σ The explicit form of the momentum (Eq.16) and Poisson equations (Eq.17) are To cast the governing equations in terms U m n , V m n , W m n , R m n , S m n , T m n , ϕ m n , ∂ r ϕ m n and their derivatives, we first need to obtain the radial components of the stress tensor (James & Cook 1976, Eq. 3.22) The terms with R v arise from lateral variations and couple equations involving different modes (see Eq. C32).Without lateral variations, Eq. (C36) reduces to the equations typical of a sphericallysymmetric body.Using linear algebra, the previous set of equations can be cast as ∂ r y = Dy.The ϕ (T ) 0,+ 2 , ϕ (T ) 0,− 2 ; ϕ (T ) 2,+ 2 , ϕ (T ) 2,− 2 ; and ϕ (T ) −2,+ 2 , ϕ (T ) −2,− 2 correspond to standing, westward propagating and eastward propagating waves, respectively.As the tidal potential is real, we have that ϕ (T ) −m,− n = (−1) m ϕ (T ) m,+ n .The amplitude of the different components of the tide are given in Table 2.

E. TEST ON THE TOTAL ENERGY DISSIPATION
We test that the total energy dissipation obtained using the stress and strain tensors (Eqs.( 29)-( 34)) and considering the work done by the tidal potential (Eqs.( 35)-( 37)) is the same.We consider a viscoelastic Io that is spherically-symmetric and one that has zonal degree 2 lateral viscosity variations with peak to peak amplitude of 50% the mean viscosity and compute tidal heating using Eq. ( 34) ( Ė1 ) and (37) ( Ė2 ).We find excellent agreement between the total energy dissipation obtained with the two methods, with the difference decreasing as the number of radial points used for the radial integration increases (Figure 9).Equivalent results are obtained for lateral variations of different amplitude, and degree and order.35), Ė2 ).We assume the viscoelastic Io model of Table 1 and zonal degree 2 lateral viscosity variations with peak to peak amplitude of 50% the mean viscosity.

Figure 2 .
Figure 2. Gravitational Love numbers for different types of lateral shear modulus variations (upper panels); and difference with the Love numbers obtained with the perturbation approach of Qin et al. (2014) (lower panels).An elastic Io (Table1) with shear modulus variations of different wavelengths (n LV , m LV ) is assumed.The amplitude of lateral variations is given in terms of the peak-to-peak variation of the shear modulus with respect to its mean.The solution is cut-off at perturbation order 3.The line-width indicates the order of the mode, with thicker lines indicating lower order modes.In the upper panels the dashed lines corresponds to the solution obtained with the perturbation method ofQin et al. (2014).For (n LV , m LV ) = (1, 1) only the +|m| modes are shown, −|m| modes of amplitude k Figure 2. Gravitational Love numbers for different types of lateral shear modulus variations (upper panels); and difference with the Love numbers obtained with the perturbation approach of Qin et al. (2014) (lower panels).An elastic Io (Table1) with shear modulus variations of different wavelengths (n LV , m LV ) is assumed.The amplitude of lateral variations is given in terms of the peak-to-peak variation of the shear modulus with respect to its mean.The solution is cut-off at perturbation order 3.The line-width indicates the order of the mode, with thicker lines indicating lower order modes.In the upper panels the dashed lines corresponds to the solution obtained with the perturbation method ofQin et al. (2014).For (n LV , m LV ) = (1, 1) only the +|m| modes are shown, −|m| modes of amplitude k n,−|m| 2,0 = (−1) |m| k n,|m| 2,0 are also excited.The difference between

FEM
model or the spectral perturbation method.Solutions obtained using the FEM of Berne et al. (2023a) and the perturbation method of Qin et al. (2014) are indicated with filled circles and dashed lines, respectively.

Figure 4 .
Figure 4. Love number spectra for a degree two order zero (upper panel), one (center panel) and two (lower panel) tidal forcing for lateral variations of different degrees and orders (n LV , m LV).An elastic Io is assumed (Table1), and the solution is cut-off at perturbation order 4. Color intensity indicates the strength of the mode, and the terms with highest amplitude are indicated with a cross.For all cases, we assume peak-to-peak shear modulus variations of 10% of the mean shear modulus.

Figure 5 .
Figure 5. Gravitational Love numbers, k n,m 2,0 , for lateral shear modulus variations of different wavelengths (n LV , m LV ) for Io (solid lines), the Moon (dashed lines) and Mercury (dotted lines).

Figure 7 .
Figure7.Pattern of tidal heating for a spherically-symmetric body (upper row) and for bodies with different types of lateral viscosity variations.For the spherically-symmetric body the total heating pattern is shown while for the rest the difference between their corresponding heating pattern and that of a spherically-symmetric body is shown.All plots are normalized by the average tidal heating of

Figure 8 .
Figure 8. Tidal heating spectra (Eq.(31)) corresponding to the cases of Figure 7.The energy spectra is given in terms of amplitude (upper panel) and longitudinal eastward shift in fraction of longitudinal wavelength with respect to ℜ(Y m n ) (lower panel).

M
.R. and I.M. were supported by the National Aeronautics and Space Administration (NASA) under grant No. 80NSSC20K0570 issued through the NASA Solar System Workings program.A.B was supported by the Future Investigators in NASA Earth and Space Science and Technology (FI-NESST) Program (80NSSC22K1318).The authors thank Allard Veenstra for his contributions to the LOV3D software and for providing feedback on the original manuscript.The LOV3D software can be downloaded and run at https://github.com/mroviranavarro/LOV3Dopen.Software: LOV3D (Rovira Navarro & Matsuyama 2024), Wigner 3j-6j-9j (Sovkov 2024), cmocean (Thyng et al. 2016), M Map (Pawlowicz 2020), export fig Altman (2024), Tidal-Response Qin (2016) APPENDIX A. TENSOR SPHERICAL HARMONICS Tensor spherical harmonics are a generalization of spherical harmonics to tensors of rank > 0.
with the coupling coefficients.If there are no lateral variations, R v is 0.

Figure 9 .
Figure9.Difference in tidal heating computed using the stress and strain tensors (Eq.(29), Ė1 ) and the work done by the tidal force (Eq.(35), Ė2 ).We assume the viscoelastic Io model of Table1and zonal degree 2 lateral viscosity variations with peak to peak amplitude of 50% the mean viscosity.