We have derived a convergent scattering series solution for the frequency-domain wave equation in acoustic media with variable density and velocity. The convergent scattering series solution is based on the homotopy analysis of a vectorial integral equation of the Lippmann-Schwinger type. By using the Green's function and partial integration, we have derived the vectorial integral equation of the Lippmann-Schwinger type that involves the pressure gradient field as well as the pressure field from the wave equation. The vectorial Lippmann-Schwinger equation can in principle be solved via matrix inversion, but the computational cost of matrix inversion scales like N 3 , where N is the number of grid blocks. The computational cost can be significantly reduced if one solves the vectorial Lippmann-Schwinger equation iteratively. A simple iterative solution is the Born series, but it is only convergent when the scattering potential is sufficiently small. In this study, we have used the so-called homotopy analysis method to derive an iterative solution for the vectorial Lippmann-Schwinger equation which can be made convergent even in strongly scattering media. The computational cost of our convergent scattering series scales as N 2 . Our algorithm, which is based on the homotopy analysis method, involves a convergence control operator that we select using hierarchical matrices. We use a three-layer model and a resampled version of the SEG/EAGE salt model to show the performance of the developed convergent scattering series.