A finite ice shell, from boundary tractions to Love numbers

Choose a world with one tractable interface

Use a spherical world of radius five hundred kilometres and constant density one thousand kilograms per cubic metre. Its solid shell and inviscid liquid interior have the same density. Both are incompressible. The shell has uniform Maxwell properties, the liquid has no shear rigidity, and the tide changes slowly enough that inertia is omitted from the internal force balance.

Equal density does real work in this simplification. Moving the internal boundary creates no density contrast and hence no additional interfacial gravitational mass. The outer displaced surface produces the entire induced gravity potential. For a degree-two tide its gravity Love number is three-fifths of the radial Love number. Density stratification would remove this fixed proportionality.

Define the liquid-radius fraction as one minus thickness divided by total radius. At zero liquid radius the body is fully solid. The solid branch uses a regular displacement field all the way to the centre. No rock core, compressibility, ocean resonance, rotation-induced internal dynamics or depth-dependent ice temperature is present.

g=\frac{4\pi G\rho R}{3},\quad x=1-\frac dR,\quad k_2=\frac35h_2

The surface gravity is 0.139786 metres per square second in the saved synthetic example.

Reading: Beuthe, Tides on Europa: the membrane paradigm, arXiv:1410.4735v1

Solve the shape before its overall amplitude

Separate the angular degree-two pattern from its radial dependence. Normalize the surface radial shape to one and write the internal radial function as a combination of four powers. Incompressibility fixes the tangential function. The elastic equilibrium equation supplies a pressure term, leaving four coefficients that can be determined from the boundary conditions. Here the radial coordinate is divided by the total radius.

The outer surface has zero tangential traction. At the liquid interface, tangential traction also vanishes and the elastic perturbation of normal traction vanishes after the common fluid hydrostatic load is removed. Together with surface normalization these are four linear equations. The Lab solves that system directly. At the fully solid limit regularity excludes the negative powers, and the two outer conditions give the separate central-regular solution.

The dimensionless outer normal traction defines the geometric stiffness. Its value rises from approximately 24 times thickness divided by eleven radii for a thin shell to 19 divided by five for a solid sphere. For the default fifty-kilometre shell it is 0.232691. The tangential surface displacement is obtained from the same coefficients, so the visible shell boundary and the energy integral share one mechanical solution.

U(s)=a s+b s^3+c s^{-2}+d_0 s^{-4},\quad V(s)=\frac{sU^\prime+2U}{6},\quad p(s)=7b s^2+2c s^{-3}

The symbol for the last coefficient is separate from physical shell thickness. Pressure and traction share the omitted displacement/modulus normalization.

T_t=V^\prime+\frac{U-V}{s},\quad T_r=-p+2U^\prime,\quad U(1)=1,\quad T_t(1)=T_t(x)=T_r(x)=0

The four conditions apply for a nonzero liquid radius. For a solid sphere the regular radial function has linear and cubic coefficients eight-fifths and minus three-fifths.

z_h=\frac{24}{5}\frac{19-75x^3+112x^5-75x^7+19x^{10}}{24+40x^3-45x^7-19x^{10}}

The published closed form checks a displacement field obtained independently from the numerical boundary solve.

Reading: Beuthe, Tides on Europa: the membrane paradigm, arXiv:1410.4735v1

Let self-gravity oppose the elastic response

A radial displacement changes gravitational potential and hydrostatic pressure. At the surface, inward hydrostatic restoration is opposed by the applied and induced tidal potentials, with an additional elastic traction. Divide the normal balance by density times the applied potential. The result contains the radial Love number, minus one, minus the gravity Love number, plus the geometric elastic resistance.

Insert the equal-density relation between the Love numbers and solve. The fluid gravitational coefficient is two-fifths. Replacing elastic shear modulus by its complex Maxwell value is the linear correspondence step. The real restoring coefficient stays positive, while the passive material’s positive imaginary stiffness gives a negative imaginary displacement response.

A zero-frequency Maxwell body has enough time to relax its shear stress. Its radial Love number tends to two-and-a-half and its gravity Love number to one-and-a-half. An ideal elastic body instead retains finite stiffness and has no imaginary response. These are different limiting operations. The liquid limit and elastic zero-loss limit are tested separately.

h_2-1-k_2+z_h\frac{\mu^*}{\rho gR}h_2=0\quad\Longrightarrow\quad h_2=\left[\frac25+z_h\frac{\mu^*}{\rho gR}\right]^{-1},\qquad k_2=\frac35h_2

Density times gravity times radius is a pressure; the entire denominator is dimensionless.

Reading: Beuthe, Tides on Europa: the membrane paradigm, arXiv:1410.4735v1

Calculate the default shell from start to finish

Take shell thickness fifty kilometres, eccentricity 0.005, elastic modulus 3.5 gigapascals and viscosity ten to the power 14.5 pascal-seconds. The liquid-radius fraction is 0.9. Solving the four traction equations gives geometric stiffness 0.232691; the world’s density, gravity and radius multiply to a pressure scale about 69.8931 million pascals.

At the orbital frequency the material chapter gives complex modulus approximately 3.18439 plus 1.00251 imaginary gigapascals. Divide it by the pressure scale, multiply by geometric stiffness, and add the fluid coefficient two-fifths. The complex denominator is about 11.0016 plus 3.33760 imaginary. Invert it by multiplying numerator and denominator by its conjugate, then take three-fifths for the gravity response.

The surface potential amplitude varies with direction. Maximizing its magnitude over the sphere gives the expression below. Multiply it by the radial Love magnitude and divide by gravity to obtain 3.08319 metres. Its phase lag is 16.8764 degrees. Finally the full work formula gives 100.009 gigawatts; the independent volume dissipation check uses the internal strain field, not this final scalar power.

h_2=\frac1{11.0016+3.33760i}\simeq0.0832352-0.0252513i,\qquad k_2\simeq0.0499411-0.0151508i

The pressure, geometric stiffness and modulus all appear in this explicit numerical substitution.

|\widehat W|_{\max}=\frac32\sqrt{\frac{32}{7}}\frac{GM_pR^2}{a^3}e,\qquad |\widehat u_r|_{\max}=\frac{|h_2|}{g}|\widehat W|_{\max}

The maximum is the largest harmonic radial amplitude among all surface directions. The scene samples its real displacement at a selected phase.

Reading: Beuthe, Tides on Europa: the membrane paradigm, arXiv:1410.4735v1

Keep both parts of the changing eccentricity tide

Set the body-fixed first axis toward the mean planet direction and the third axis along the spin axis. Synchronous mean rotation removes the average orbital angle. To first order in eccentricity, distance variation changes the tidal amplitude while optical libration rotates its direction. The complex symmetric trace-free tensor below includes both effects. The steady tide is removed.

For a unit radial direction, sandwiching the tensor between that direction gives the complex surface potential. Its radial response is the Love number times that potential divided by gravity. Inside the shell the radial and tangential functions from the traction solve supply the full vector displacement. In the liquid, a regular incompressible potential-flow continuation matches radial motion and allows tangential slip.

The orbital drawing uses a Kepler ellipse while the forcing keeps only the leading eccentricity term. Their discrepancy starts at second order in eccentricity. The allowed eccentricities are small; the experiment does not add nonlinear orbital harmonics merely by drawing a more eccentric ellipse. Larger eccentricity would require upgrading the forcing model and internal response together.

\widehat W(\boldsymbol n)=\boldsymbol n^{\mathsf T}\widehat{\boldsymbol S}\boldsymbol n,\qquad \widehat{\boldsymbol S}=\frac{GM_pR^2}{2a^3}e\begin{pmatrix}6&-6i&0\\-6i&-3&0\\0&0&-3\end{pmatrix}

The imaginary off-diagonal terms encode optical libration under the positive-time-exponential convention.

\widehat{\boldsymbol u}(r,\boldsymbol n)=\frac{h_2}{g}\left[U(r/R)\widehat W\boldsymbol n+V(r/R)\boldsymbol\nabla_\Omega\widehat W\right]

At the outer surface radial normalization is one. The liquid continuation is proportional to radius and is regular at the centre.

Reading: Beuthe, Tides on Europa: the membrane paradigm, arXiv:1410.4735v1

Integrate work rather than painting a heat map

Incompressibility converts tidal work in the volume to a surface integral of potential times radial velocity. Averaging the product of two real harmonics yields half the real part of one amplitude times the conjugate of the other. Substituting the displacement response leaves minus the imaginary Love number. A passive lag therefore gives positive power.

The angular integral must include the full tensor. For a symmetric trace-free tensor the integral of the squared potential amplitude is eight pi over fifteen times its squared tensor norm. The diagonal and optical-libration terms together give the familiar factor of twenty-one over two in the final power. Dropping the direction oscillation changes that coefficient.

At the saved default, external-work integration gives about 100.009 gigawatts. A separate volume calculation takes numerical Cartesian derivatives of the displacement, builds the symmetric strain tensor and integrates viscous shear loss through the shell. It agrees with the external work. Neither calculation interprets the scene’s blue and white materials as temperature.

\mathcal P=\rho R^2\int\left\langle W\dot u_r\right\rangle d\Omega=-\frac{\rho R^2\omega}{2g}\operatorname{Im}h_2\int|\widehat W|^2d\Omega

The sign follows directly from the declared harmonic convention and radial displacement response.

\int|\widehat W|^2d\Omega=\frac{84\pi}{5}\left(\frac{GM_pR^2}{a^3}\right)^2e^2,\quad \mathcal P=-\frac{21}{2}\frac{GM_p^2R^5}{a^6}\omega e^2\operatorname{Im}k_2

For the physical synchronous eccentricity tide the angular frequency is the orbital mean motion. Harmonic diagnostic probes hold spatial forcing amplitude fixed while changing frequency.

Reading: Beuthe, Tides on Europa: the membrane paradigm, arXiv:1410.4735v1 · Petricca et al., Titan’s strong tidal dissipation precludes a subsurface ocean, Nature 648 (2025)

Saved numerical example, units and method

Synthetic incompressible equal-density world. The parent mass is five hundred trillion trillion kilograms and the orbital semimajor axis is three hundred thousand kilometres. A four-condition elastic traction solve fixes the shell displacement; complex Maxwell correspondence fixes lag and heating. These values are predictions, with no observational noise except in the explicitly labelled fitted rows.

Radius, shell thickness and eccentricity	$R=500\,\mathrm{km},\quad d=50\,\mathrm{km},\quad e=0.00500$
Density in both layers and surface gravity	$\rho=1000\,\mathrm{kg\,m^{-3}},\quad g=0.139786\,\mathrm{m\,s^{-2}}$
Finite-shell geometric stiffness	$z_h=0.232691$
Radial and gravity response	$h_2=0.0832352-0.0252513i,\quad k_2=0.0499411-0.0151508i$
Largest radial amplitude and total power	$u_{r,\max}=3.08319\,\mathrm{m},\quad\mathcal P=100.009\,\mathrm{GW}$

The displayed scene removes the steady tide and magnifies the cyclic radial shape. Its gain is normalized separately at the reference eccentricity for each material and thickness; physical metre-valued amplitudes remain available for absolute comparisons. Independent external-work and volume-dissipation integrals test the power.

Test a boundary condition by moving the boundary

Drag from a thin shell toward a thick shell. Explain the changing liquid radius and compare the physical metre-valued displacement with the displayed gain.
Double eccentricity and predict the factor-of-four power change before checking the readout. Then switch to the circular preset.

Exercises with explanations

Why does the solid limit need a separate regular solution rather than evaluating negative radial powers at the centre?
Read the explanation
Negative powers diverge at zero radius. Their coefficients must vanish before taking the limit. The remaining regular solution gives geometric stiffness 19/5, and reproduces the homogeneous incompressible solid Love response.
Double eccentricity while retaining orbit scale and material. How should displacement and power change?
Read the explanation
The leading tidal potential and displacement double, while power quadruples because it contains the product of forcing and velocity. The display uses the same reference-eccentricity gain for that material, so its visible cyclic deformation doubles too.
A real moon has a denser rocky core. Can the fixed three-fifths relation between gravity and radial Love numbers still be imposed?
Read the explanation
No. Density contrasts add displaced internal mass and alter the gravitational boundary problem. One must solve the stratified structure before comparing the two observations; the proportionality is an assumption of this equal-density example.

Sources and calculation scope

Beuthe, Tides on Europa: the membrane paradigm, arXiv:1410.4735v1
Fixed 17 October 2014 preprint; Icarus publication DOI 10.1016/j.icarus.2014.10.027. Appendix A equations 105–116 give the finite-thickness homogeneous incompressible shell; Appendix C gives Maxwell correspondence; Appendix E derives work and power; Appendix G fixes the eccentricity/libration tide. Figures 8–9 and 12–13 expose approximation and thermal-model dependence. Full relevant sections and appendices were read. Our uniform-density Maxwell example is much simpler than a real moon.
Petricca et al., Titan’s strong tidal dissipation precludes a subsurface ocean, Nature 648 (2025)
Version of record, 17 December 2025, DOI 10.1038/s41586-025-09818-x, CC BY 4.0; full publisher-deposited text and Methods, PMC12711566, read. Figure 1 reports gravity-response estimates, Figure 2 compares conditional interior posteriors and Figure 3 is an inferred-interior schematic. Figure 1 error bars are three-sigma; Figure 2 ellipses are one-sigma. Positive reported imaginary response denotes their loss convention; our positive-time-exponential convention uses negative imaginary Love numbers for passive loss.