A reduced viscoelastic FDTD formulation for ultrasound-driven shear wave propagation in soft tissue

Ultrasound-driven shear wave propagation in soft tissue underlies shear wave elastography (SWE) and emerging elastomechanical hypotheses of ultrasonic neuromodulation, both of which require accurate, efficient modeling of radiation-force--induced tissue motion. General-purpose finite-element elastodynamic solvers are often computationally expensive and unnecessarily broad for shear-dominant applications. We derive a reduced viscoelastic formulation by applying near-incompressibility, small-strain linearization, Helmholtz decomposition, and solenoidal force projection to the full Navier equations, yielding a Kelvin--Voigt shear wave equation that retains only the transverse dynamics relevant to radiation-force--induced motion. An explicit finite-difference time-domain (FDTD) implementation with second-order spatial and temporal accuracy enforces the solenoidal body-force constraint via a matrix-free conjugate-gradient Poisson solve. For separable radiation-force sources, a pre-computed projection reduces this cost by one to two orders of magnitude. In homogeneous media the solver recovers the theoretical shear wavespeed to within $<$0.1\%, exhibits clear second-order grid convergence (trace-$L_2$ self-convergence, $p\gtrsim2$), and matches analytical Kelvin--Voigt attenuation and phase speed to within ${\sim}3\%$ and $<$1\% over a 16-point viscosity sweep ($\eta\in[0,1.5]$~Pa$\cdot$s). An end-to-end RSNA QIBA phantom benchmark recovers shear wave speeds to within ${\sim}1\%$ across a tenfold shear-modulus range ($G=1$--$10$~kPa). The framework accommodates spatial heterogeneity in shear modulus, density, and viscosity, and integrates with acoustic simulators. Transcranial demonstrations through micro-CT skull geometries produce shear displacements of 1.7--5~$\mu$m consistent with clinical ARFI.