Iterative Born solver for the acoustic Helmholtz equation with heterogeneous sound speed and density

Efficient numerical solution of the acoustic Helmholtz equation in heterogeneous media remains challenging, particularly for large-scale problems with spatially varying density-a limitation that restricts applications in biomedical acoustics and seismic imaging. A fast iterative solver that extends the convergent Born series [Osnabrugge, Leedumrongwatthanakun, and Vellekoop, J. Comput. Phys. 322, 113-124 (2016)] method to handle arbitrary variations in sound speed, density, and absorption simultaneously is presented. This approach reformulates the Helmholtz equation as a first-order system and applies the universal split-preconditioner from Vettenburg and Vellekoop [arXiv:2207.14222v2 (2022)], yielding a matrix-free algorithm that leverages Fast Fourier Transforms for computational efficiency. Unlike existing Born series methods, this solver accommodates heterogeneous density without requiring expensive matrix decompositions or preprocessing steps, making it suitable for large-scale three-dimensional problems with minimal memory overhead. The method provides forward and adjoint solutions, enabling its application for inverse problems. Accuracy is validated through comparison against an analytical solution and the solver's practical utility is demonstrated through transcranial ultrasound simulations. The solver achieves convergence for strong scattering scenarios, offering a computationally efficient alternative to time-domain methods and matrix-based Helmholtz solvers for applications ranging from medical ultrasound treatment planning to seismic exploration.