A fullwave model of the nonlinear wave equation with multiple relaxations and relaxing perfectly matched layers for...

Large-scale acoustic simulations underpin the development of ultrasound imaging and therapy, but modeling nonlinearity, frequency-dependent attenuation, and absorbing boundaries in heterogeneous tissue remains computationally demanding. We present Fullwave 2, a unified time-domain formulation with arbitrary power-law tissue attenuation and perfectly matched layers (PMLs) within a single framework suited to high-order finite difference solution. Attenuation and dispersion are encoded directly into complex coordinate-stretched spatial derivatives through multiple relaxation mechanisms. Because the same mechanism describes both interior tissue attenuation and the absorbing boundary, the convolutional PML (C-PML) becomes a special case of the domain-wide model and adds no extra computational burden. The formulation preserves the structure of the d'Alembertian operator, which allows high-order staggered-grid finite difference stencils optimized for long-distance propagation, and a two-stage C-PML with a transition region is introduced to ensure numerical stability in the presence of multiple relaxations. The domain-wide multiple relaxation model reproduces power-law attenuation with less than 5% attenuation error and 0.5% phase-velocity error over 1-20 MHz. The two-stage C-PML reaches reflection coefficients below -49 dB with a compact 4 lambda footprint. Nonlinear propagation is validated against a 1D Burgers solution, with agreement up to the 7^(th) harmonic. The framework is demonstrated on 2D abdominal wall imaging and 3D transcranial rat skull simulations, where it accurately captures complex scattering and aberration artifacts. Fullwave 2 unifies nonlinear propagation, arbitrary power-law attenuation, and absorbing boundaries in a single, computationally efficient time-domain formulation, providing an accurate and scalable wave propagation tool for medical ultrasound research.