A fullwave model of the nonlinear wave equation with multiple relaxations and relaxing perfectly matched layers for high-order numerical finite difference solutions

Large-scale acoustic simulation is essential for ultrasound imaging and therapy but faces computational bottlenecks when modeling nonlinearity, frequency-dependent attenuation, and absorbing boundaries in heterogeneous media. To address this, we present Fullwave 2, a unified time-domain formulation that integrates power-law tissue attenuation and Perfectly Matched Layers (PMLs) via complex coordinate stretching. This framework preserves the structure of the d'Alembertian operator, enabling the use of high-order staggered-grid finite difference stencils optimized for long distance propagation. Since the attenuation in the interior of the domain uses the same formulation as the boundaries, sophisticated implementations of the PMLs are easily implementable and do not add any computational burden. Here a two-stage Convolutional Perfectly Matched Layer (C-PML) with a transition region is demonstrated to be highly stable. The PML attains high absorption efficiency, achieving reflection coefficients below -50dB with a compact 4λ footprint. The domain-wide multiple relaxation mechanisms achieve <5% attenuation error and <0.5% phase-velocity error over a 1--20 MHz bandwidth. Validation against a 1D Burgers solution confirms nonlinear accuracy up to the 7th harmonic. The capabilities of the simulation tool are finally demonstrated in ultrasound imaging of 2D abdominal phantoms, and in 3D transcranial applications demonstrating how the method accurately captures complex scattering and aberration artifacts, and provides an efficient and acoustical wave propagation tool for medical ultrasound research.