Elastic Wave Equation Fundamentals

• Elastic wave equation is a fundamental partial differential equation that models wave propagation in elastic materials by linking displacement, stress, and strain.
• It decouples into distinct P (compressional) and S (shear) wave modes, each with unique propagation characteristics essential for seismic and material analysis.
• Advanced numerical and inverse techniques apply the elastic wave equation to efficiently simulate elastic systems and accurately recover material properties.

The elastic wave equation is the fundamental partial differential equation describing the propagation of deformations—specifically, mechanical waves—in elastic continua. Its precise mathematical structure encodes the interplay between translation, rotation, stress, and strain in response to local and remote forces, and underlies the emergence of both longitudinal (P) and transverse (S) waves. The equation is central to seismology, nondestructive testing, medical elastography, and a wide range of solid mechanics and material science applications.

1. Kinematics, Constitutive Laws, and Dynamic Equations

The elastic wave equation arises by combining kinematic relations, constitutive laws, and fundamental balance principles. In the classical (Cauchy) framework for isotropic, linear elasticity, the displacement field $$u: \mathbb{R}^n \times [0, T] \rightarrow \mathbb{R}^n$$ satisfies

$\rho\,\partial_t^2 u = \nabla \cdot \sigma(u) + f,$

where $\rho(x) > 0$ is the mass density and $f$ is the body force. The Cauchy stress tensor $\sigma(u)$ is related to strain via

$$\sigma_{ij}(u) = \lambda\, (\nabla \cdot u)\,\delta_{ij} + 2\mu\,\varepsilon_{ij}(u),$$

with the symmetric strain tensor

$$\varepsilon_{ij}(u) = \tfrac{1}{2} (\partial_i u_j + \partial_j u_i).$$

Here $\lambda, \mu$ are Lamé parameters satisfying $\mu > 0$, $\lambda + 2\mu > 0$ for ellipticity. In fully anisotropic media with smooth domain $M \subset \mathbb{R}^n$, the law generalizes to

$$\rho(x)\,\partial_t^2 u_i(x,t) - \partial_j\big(C_{ijkl}(x)\,\partial_l u_k(x,t)\big) = 0,$$

where $C_{ijkl}$ is the stiffness tensor with the usual major/minor symmetries (Ilmavirta et al., 20 Nov 2025).

Recent work revisits this standard structure by modeling each material element not as a simple particle with only translational degrees of freedom, but as a local rigid body with both translational and rotational components (Shi, 2020). In this approach, the full displacement gradient $\xi_{ij} = \partial_j u_i$ is decomposed into symmetric (classical strain) and antisymmetric (local rotation) parts,

$$\xi_{ij} = \varepsilon_{ij} + \omega_{ij}, \quad \varepsilon_{ij} = \tfrac{1}{2}(\partial_j u_i + \partial_i u_j), \quad \omega_{ij} = \tfrac{1}{2}(\partial_j u_i - \partial_i u_j).$$

The constitutive relation becomes

$$\sigma_{ij} = \lambda \delta_{ij} \xi_{kk} + 2\mu \xi_{ij},$$

making $\sigma$ generally asymmetric and restoring conservation of total angular momentum by explicitly incorporating local spin dynamics.

2. Wave Mode Decomposition and Physical Interpretation

Diagonalization of the elastodynamic operator in the Fourier domain (Lamé symbol) reveals two distinct wave families:

• Longitudinal (P) waves: Associated with the eigenvalue $(\lambda+2\mu)|\xi|^2$ and eigenvector parallel to $\xi$. These correspond to compressional, dilatational motion.
• Transverse (S) waves: With eigenvalue $\mu|\xi|^2$ and eigenvectors orthogonal to $\xi$. These make up shear, rotational motions (Kim et al., 2021, Brytik et al., 2010).

In isotropic media, the principal symbol decomposes as

$$\sigma_2(A)(x,\xi) = (2\mu+\lambda)|\xi|^2 \Pi_P(\xi) + \mu|\xi|^2 \Pi_S(\xi),\qquad \Pi_P = \frac{\xi \otimes \xi}{|\xi|^2}.$$

It is possible to construct a (microlocally) diagonalizing operator such that, up to smoothing remainders, the elastic evolution decouples into two systemically independent wave equations—one for each mode (Brytik et al., 2010).

Physically, the symmetric (translational) strain produces P-waves, while the antisymmetric (rotational) part is the source of S-waves (Shi, 2020). Recent advances clarify that conservation of angular momentum mandates explicit inclusion of these rotational degrees of freedom, removing the necessity for a priori stress symmetry.

3. Extensions: Anisotropy, Heterogeneity, and Nonlinearity

Anisotropy modifies the principal symbol and dramatically alters dispersion and wavefront geometry. The stiffness tensor $C_{ijkl}(x)$ introduces direction-dependent phase velocities via the generalized Christoffel matrix

$\Gamma_{i\ell}(x,p) = C_{ijkl}(x) p_j p_k.$

Propagation reduces to the paper of Finsler geometry via the qP-branch of $\Gamma_{i\ell}$: if the leading eigenvalue $\lambda_{qP}$ is globally simple and smooth, associated Finsler-geometric inverse problems (such as boundary rigidity and X-ray transforms) can be rigorously posed (Ilmavirta et al., 20 Nov 2025).

Heterogeneity in density and elasticity coefficients mandates careful treatment of regularity. With coefficients in $C^{1,1}$, full decoupling (up to smoothing operators) and energy estimates remain valid (Brytik et al., 2010). In 1D, coordinate transformations require precise adjustment of pre-stress gradients to preserve form-invariance and time-synchronization, with only high-frequency limits allowing further reduction to classical forms (Yao et al., 2014).

Nonlinearity further extends the equation via higher-order strain terms, as in the Green–Saint-Venant (finite strain) expansion: $$\epsilon_{ij}(u) = \tfrac{1}{2}\big(\partial_j u_i + \partial_i u_j + \partial_k u_i \partial_k u_j\big).$$ The stress tensor admits third-order moduli $𝒜, 𝒝, 𝒞$ in addition to the Lamé parameters, leading to nonlinear wave interactions. Inverse boundary value problems in this broader context demonstrate that full recovery of all coefficients (including density) from Dirichlet-to-Neumann data is possible under appropriate geometric control (Uhlmann et al., 2023).

4. Dispersion, Attenuation, and Fractional Models

In real media, elastic waves suffer frequency-dependent attenuation and phase shifts. The fractional Zener elastic wave equation replaces the integer-order stress-strain relation by a fractional Caputo derivative, modeling viscoelasticity: $$\sigma(t) + \tau_s^\alpha D^\alpha \sigma(t) = E_0 \left[\epsilon(t) + \tau_\epsilon^\alpha D^\alpha \epsilon(t)\right],\quad 0<\alpha\le1.$$ Proceeding from mass and momentum conservation yields

$$\nabla^2 u - \frac{1}{c_0^2} \frac{\partial^2 u}{\partial t^2} + \tau_s^\alpha \frac{\partial^\alpha}{\partial t^\alpha} \nabla^2 u - \frac{\tau_\epsilon^\alpha}{c_0^2}\frac{\partial^{\alpha+2} u}{\partial t^{\alpha+2}} = 0.$$

This model reproduces frequency power-law attenuation with three distinct regimes and offers a concrete link to multiple-relaxation mechanical analogs (Maxwell-Wiechert models) (Nasholm et al., 2012).

5. Analytical Properties, Estimates, and Regularity

Sharp dispersive, Strichartz, and Morawetz estimates have been established for the elastic wave equation, enhancing control over space-time regularity, local energy decay, and global-in-time existence:

• Dispersive and Strichartz bounds: The elastic propagator inherits the full classical wave admissibility, modulo the spectral decomposition into P- and S-waves (Kim et al., 2021). Even with singular (inverse-square) matrix-valued potentials, global Strichartz and local energy estimates persist under optimal (Fefferman–Phong) regularity conditions (Kim et al., 2020).
• Morawetz-type inequalities: For power-type singular weights in space or space-time, refined Morawetz estimates quantify solution integrability by sharp Sobolev exponents, enabling scattering/decay analysis in elastic systems (Kim et al., 8 Oct 2025, Lai et al., 20 Oct 2025). Weighted multipliers have also yielded exponential stabilization for anisotropic and/or nonlinear systems (Ning et al., 2020).

These results are complemented by uniform Sobolev inequalities for the system, with applications to unique continuation and Carleman-type inequalities (Kim et al., 2021).

6. Numerical Discretizations and Computational Methods

Multiple advanced discretization frameworks have been developed for the elastic wave equation:

• Staggered-grid finite differences: Optimal for first-order forms, the latest schemes reduce computational load by combining short (second-order) and long (high-order) stencils for different spatial derivatives, preserving accuracy and reducing time per simulation step by ≈40% (Liang et al., 2017).
• High-order cut finite elements (cut-FEM): Robustly handle geometries where the domain boundary or interfaces cut arbitrarily through the background mesh. Symmetric Nitsche's enforcement combined with ghost-penalty stabilization guarantees unconditional energy stability and optimal a priori convergence, even with small cut elements (Sticko et al., 2018).
• Virtual element methods (VEM) via scalar potentials: Decompose the displacement into decoupled P- and S-wave scalar equations, allowing separate mesh refinement and polynomial order for each, matching the wave-number-dependent numerical challenges (Falletta et al., 2023).
• Fully discrete boundary integral methods: Exploit the Calderón calculus with precise Petrov–Galerkin discretization, staggered grids, and "look-around" quadrature for high order, uniform convergence in both time and space, well suited for exterior and crack problems in 2D (Dominguez et al., 2014).
• Frozen Gaussian approximation (FGA): High-frequency regime is treated by representing the wavefield as a superposition of Gaussian wavepackets propagated along rays, efficiently capturing P and S waves, polarization coupling (Berry-phase effects), and transmission/reflection at interfaces, drastically reducing computational cost in seismic tomography and inversion (Hateley et al., 2018).

7. Inverse Problems and Geometric Methods

The interplay of geometric analysis and the elastic wave equation enables powerful inverse problem statements:

• Boundary rigidity and X-ray problems: Under certain convexity and simplicity conditions on domains and wavespeeds, boundary measurements (displacement-to-traction map) determine all linear and nonlinear elastic parameters pointwise (Uhlmann et al., 2023).
• Finsler geometry and slowness surfaces: For anisotropic elasticity, the leading eigenbranch ('qP' slowness) defines a Finsler metric whose geodesics model P-wave rays. Minimal smoothness and algebraic separation correspond to explicit, checkable conditions on $C_{ijkl}(x)$, facilitating uniqueness in travel time tomography and related inverse problems (Ilmavirta et al., 20 Nov 2025).

These analyses clarify the minimal regularity and algebraic assumptions necessary for the geometric framework to apply and directly impact global uniqueness results and the stability of recovery algorithms.

---

References:

• New motion modeling and derivation of rotationally-aware elastic wave equation: (Shi, 2020)
• Fractional Zener equation, attenuation, and viscoelastic modeling: (Nasholm et al., 2012)
• High order CutFEM for elastic wave equation: (Sticko et al., 2018)
• Inverse boundary determination in nonlinear elasticity: (Uhlmann et al., 2023)
• Stabilization and multiplicity of decay mechanisms using Morawetz multipliers: (Ning et al., 2020)
• Morawetz-type estimates and singular weights: (Kim et al., 8 Oct 2025, Lai et al., 20 Oct 2025)
• Strichartz and uniform Sobolev estimates, spectral decomposition: (Kim et al., 2021)
• Local energy decay and singular potential estimates: (Kim et al., 2020)
• Fully discrete Calderón calculus for elastic waves: (Dominguez et al., 2014)
• Frozen Gaussian Approximation and seismic applications: (Hateley et al., 2018)
• Efficient staggered grid finite difference methods: (Liang et al., 2017)
• P/S mode decoupling for limited smoothness coefficients: (Brytik et al., 2010)
• Horizontal/vertical regularity in Finsler-geometry-based inverse theory: (Ilmavirta et al., 20 Nov 2025)
• Covariant 1D elastic wave equations under coordinate change: (Yao et al., 2014)
• Virtual element method for time-harmonic elastic waves: (Falletta et al., 2023)