Lava Tidal Waves on Exoplanets

• Lava tidal waves are oscillatory flows in molten-surface magma oceans on highly irradiated rocky exoplanets, driven by time-dependent stellar tidal forces.
• The dynamics are modeled using linearized Laplace-tidal equations and Galerkin projection, revealing complex multi-modal interference within a hemispherical basin.
• Tidal heating from these waves creates variable surface heat patterns detectable in thermal phase curves and secondary eclipse measurements.

Lava tidal waves are oscillatory flows in the molten-surface magma oceans of highly irradiated rocky exoplanets—so-called “lava worlds”—driven by time-dependent tidal forcing from their host stars. Sustained by moderate eccentricities ($e_p \sim 0.01$–0.05), these waves constitute multi-modal dynamical responses within a shallow, viscous, hemispherical basin, producing marked variability in surface heat patterns and planetary thermal emission. Their presence distinguishes the tidal energetics and observable signatures of short-period rocky exoplanets from those of planets with solid or only weakly molten surfaces. Heat dissipated by lava tidal waves is typically regulated in the planet's interior by a combination of fluid, mushy, and solid-state convection, yielding complex equilibrium structures with potentially global-scale partial melting (Farhat et al., 11 Jan 2026).

1. Governing Equations and Dynamical Framework

The hydrodynamics of lava tidal waves on a planet of radius $R_p$ and magma ocean thickness $H$ are described by linearized Laplace-tidal (shallow-water) equations in a frame rotating at $\Omega_p$: $$\frac{\partial \mathbf{u}}{\partial t} + 2\,\boldsymbol{\Omega}_p \times \mathbf{u} + \sigma_R \mathbf{u} = -g\nabla\eta - \frac{1}{r}\nabla\Phi_{\rm tide}$$

$$\frac{\partial \eta}{\partial t} + \nabla \cdot (H\mathbf{u}) = 0$$

where $\mathbf{u}$ is the horizontal velocity field, $\eta$ the free-surface displacement, $g$ the surface gravity, $\Phi_{\rm tide}$ the tidal potential, and $\sigma_R$ the Rayleigh drag coefficient (incorporating viscous and Darcy drag, and mush dynamics). The tidal potential results from host star (mass $M_{\star}$), with harmonic decomposition capturing orbital eccentricity $e_p$ via Hansen coefficients $X_k^{n,m}(e_p)$, with forcing frequencies $\sigma_q^k = q\Omega_p - k n_{\rm orb}$.

Boundary conditions are imposed at the day–night terminator (no-flow/impermeability: $\mathbf{u} \cdot \hat{n} = 0$), effectively confining wave activity to a dynamic hemispherical basin. Dissipation localizes through Rayleigh drag, with the local dissipation per area given by

$$Q_{\rm tide}(\theta,\lambda) = \rho_m H \sigma_R \Big\langle |\mathbf{u}(\theta,\lambda,t)|^2 \Big\rangle_t,$$

where $\rho_m$ is the magma density.

2. Modal Structure, Boundary Conditions, and Eigenvalue Problem

The solution structure relies on Helmholtz decomposition of the horizontal displacement: $$\mathbf{x} = \nabla \Phi + \nabla \times (\Psi \hat{r}),$$ splitting into irrotational ($\Phi$) and solenoidal ($\Psi$) components. Both scalar fields are expanded in surface spherical harmonics within the “basin” $\mathcal{O}$, subject to:

• Neumann boundary condition on $\Phi$ $(\hat{n} \cdot \nabla \Phi = 0)$,
• Dirichlet boundary on $\Psi$ $(\Psi = 0)$. Eigenmodes satisfy

$$\nabla^2 \phi_r + \mu_r \phi_r = 0, \quad \nabla^2 \psi_r + \nu_r \psi_r = 0,$$

with indices truncating at $n_{\rm max}$ for numerical treatment. In the non-dissipative, non-rotating limit, modal frequencies are

$\omega_n^2 = gH \frac{n(n+1)}{R_p^2}$

with $k_n \approx \sqrt{n(n+1)/R_p^2}$.

3. Galerkin Solution and Wave Interference Phenomena

The linear system is advanced by Galerkin projection, expanding $\Phi$, $\Psi$ in a finite basis and projecting governing equations onto this set. In the frequency domain, each mode index $r$ yields coupled algebraic equations involving $\sigma_R$, Coriolis couplings ($\beta_{r,s}$), eigenvalues ($\mu_r, \nu_r$), and the spectral content of the tidal potential. The system is solved at each tidal forcing frequency $\sigma$, and physical velocity fields reconstructed via inverse transform and mode summation: $$\mathbf{u}(\theta,\lambda,t) = \operatorname{Re} \left\{ \sum_{r} i\sigma \left( \nabla \Phi + \nabla \times (\Psi \hat{r}) \right) \right\}.$$ The surface displacement $\eta(\theta, \lambda, t)$ emerges as a sum of oscillatory eigenmodes across $(n, m)$, and constructive or destructive interference among these “lava tidal waves” yields highly irregular, temporally variable heat patterns spatially distributed across the dayside magma ocean.

4. Tidal Heating, Wave Regimes, and Scaling Laws

Tidal dissipation rates are determined by the spatially and temporally resolved velocity fields. The dissipated power per area through Rayleigh drag,

$$\mathcal{P}_T(\theta,\lambda,t) = \rho_m H \sigma_R |\mathbf{u}(\theta,\lambda,t)|^2,$$

integrates over the basin and averages over the orbit to yield total tidal heating $\overline{\mathcal{P}_T}$. Parameter regimes are delineated via dimensionless groups:

• Eccentricity: $e_p \sim 0.01$–0.05 suffices for significant tidal forcing.
• Drag: $\widetilde{\sigma}_R = \sigma_R/\sqrt{gH/R_p}$; heavily damped ($\sigma_R \gg \sqrt{gH}k$) regimes suppress wave amplitude.
• Lamb parameter: $\epsilon_L = \Omega_p^2 R_p/g$.
• Depth: $\delta = \Omega_p^2 H / g$.
• Quality factor proxy: $Q \sim gHk^2/\sigma_R$ (very small for creep-flow magma oceans).

Regime transitions are summarized as follows:

Ocean Depth	Tidal Dissipation Scaling	Wave Phenomenology
Thin ($H \to 0$)	$\overline{\mathcal{P}_T} \propto H \to 0$	Weak flow, little dissipation
Intermediate $H$	Maximal when $\sigma \approx c_g k$ or $(c_g k)^2/\sigma_R$	Complex, multi-modal waves
Deep ($H \to R_p$)	$\overline{\mathcal{P}_T} \propto H^{-1} \to 0$	Flow confined to near-bottom

Maintained equilibrium between tidal dissipation and interior convective cooling—governed by Nusselt-Rayleigh scaling—supports self-consistent solutions for magma ocean depths of hundreds of kilometers for $P_{\rm orb} \lesssim 3$ d and $e_p \sim 0.03$.

5. Thermal Mapping, Surface Variability, and Light Curves

Tidal heat is removed through a combination of fluid convective, mushy, and solid-state mechanisms in the mantle. Local surface temperatures are established by balancing radiative losses (assuming a gray Lambertian blackbody, emissivity $\epsilon$) with stellar irradiation and internal tidal heating: $$\epsilon \sigma_{SB} T_s^4(\theta, \lambda, t) = (1-A_B) F_{\star}(t) \max[0, \cos Z] + Q_T(\theta, \lambda, t).$$ Resultant dayside heat maps exhibit marked spatial irregularity—hotspots can track east or west of the substellar point, and their locations and intensities vary aperiodically on both sub-orbital and orbit-to-orbit timescales. The observable global light curve, computed by integrating $$I(\theta,\lambda,\phi) \approx \epsilon \sigma_{SB} T_s^4/\pi$$ as a function of orbital phase, reflects amplitude and phase modulations induced by lava tidal wave interference. Measured metrics include dayside-average temperature variance, longitude drift of the thermal maximum, and flux spike irregularity.

6. Implications for Exoplanet Observations and Planetary Structure

For Earth-sized, ultra-short period planets, even modest eccentricities can induce global-scale partial or complete mantle liquefaction via sustained tidal heating. Observable photometric variability in secondary eclipse depth and phase curves may offer diagnostics of tidal heating regimes, inferable magma ocean properties, and the dynamical state of the planet. The lava tidal wave framework predicts that, contrary to models emphasizing insolation-dominated surface thermodynamics, enhanced and spatially variable tidal heating can dominate heat budgets and control surface temperature distributions on these worlds (Farhat et al., 11 Jan 2026). A plausible implication is that observed aperiodic thermal signatures could serve as indirect indicators of ongoing deep tidal liquefaction and magma ocean dynamics.

7. Theoretical and Computational Developments

The current theoretical architecture for lava tidal waves employs multi-modal Laplace-tidal equation solvers with terminator-confined basin eigenmodes, Rayleigh drag dissipation, and frequency-domain Galerkin projection. This enables systematic exploration of parameter dependencies for various tidal and rotational states, ocean depths, and drag regimes. Further development may integrate non-linear effects, feedbacks with atmospheric flows (if present), and three-dimensional mantle convection-coupled models to bridge the interface between observed surface signatures and dynamical planetary interiors (Farhat et al., 11 Jan 2026).