Pith. sign in

REVIEW 3 major objections 4 minor 28 references

Iterative Born Solver for the Acoustic Helmholtz Equation with Heterogeneous Sound Speed and Density

T0 review · 3 major / 4 minor · reviewed 2026-08-06 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read This paper claims that the acoustic Helmholtz equation with heterogeneous sound speed, density, and absorption can be solved by a matrix-free, FFT-based iterative Born-series scheme, and validates it on transcranial ultrasound.

desk verdict Extends convergent Born series to variable density via a first-order staggered-grid formulation; accuracy is well validated, but the headline efficiency claim rests on unmeasured 3D iteration counts. read the letter →

arxiv 2507.16087 v1 pith:NVV6QG5J submitted 2025-07-21 physics.comp-ph

Antonio Stanziola , Simon R. Arridge , Bradley E. Treeby , Benjamin T. Cox This is my paper · ORCID
classification physics.comp-ph MSC 65F1065N2265T5076Q05
keywords HelmholtzequationheterogeneousmediumBornseriestranscranialultrasounditerativesolvertimeharmonicacousticsdensitycontrastmatrix-freemethod
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper tries to establish that the acoustic Helmholtz equation in heterogeneous media can be solved quickly and at large scale even when density varies in space. Existing Born-series solvers handle sound-speed and absorption contrasts but not density without expensive matrix preprocessing. The proposed solver rewrites the Helmholtz equation as a first-order system of pressure and particle velocity and applies the universal split-preconditioner, so each iteration only needs diagonal scalings and FFTs. If correct, this gives a fast, low-memory forward and adjoint solver for problems like transcranial ultrasound and seismic imaging.

What carries the argument

The central object is the universal split-preconditioner for the first-order acoustic system. The forward operator is split as $A = L + V$, where $L$ contains the gradient and divergence blocks and $V$ is a diagonal operator holding the material properties; the preconditioner $\Gamma^{-1} = \nu^{-1}(I - V)(L + I)^{-1}$ is then iterated as a Neumann series. The load-bearing part is that $(L + I)^{-1}$ is computed exactly and cheaply by FFTs on a staggered grid, so the method stays matrix-free and memory-light. The staggered grid also preserves the skew-Hermitian structure of $L$ at the discrete level, which is what keeps the iteration convergent.

What would settle it

A concrete check: run the algorithm on a fixed 3D problem with strong density contrast while sweeping the two shift parameters around the median values; if the iteration count varies strongly and no plateau appears near the median, the heuristic selection is not robust.

Watch

Extended reading notes

Core claim

On its own terms, the paper's central claim is that reformulating the second-order Helmholtz equation as a first-order system makes the operator split into a skew-Hermitian differential part and an accretive diagonal material part, and the universal split-preconditioner then guarantees Neumann-series convergence for arbitrary sound speed, density, and absorption contrasts. Because the preconditioner's inverse is evaluated in the Fourier domain on a staggered grid, the whole iteration is matrix-free and uses eight 3D FFTs per step, with no LU-type decomposition or hierarchical matrix preprocessing. The paper demonstrates sub-1% error at the focus of a transcranial ultrasound benchmark compared with an established time-domain solver and agreement with an analytical fluid-sphere solution.

Load-bearing premise

The solver's practical speed rests on a heuristic choice of two tuning constants inside the preconditioner; the paper does not prove that this choice keeps iteration counts low for all strong contrasts.

Editorial extensions

If this is right

  • Forward and adjoint time-harmonic solutions are computed by the same matrix-free iteration, so inverse problems such as full-waveform inversion can reuse the solver directly.
  • Large 3D domains become feasible because memory scales with the grid and per-iteration cost is eight FFTs, with no setup-phase matrix factorization.
  • Clinical transcranial ultrasound planning can use a frequency-domain solver with sub-1% focus error relative to a time-domain reference, avoiding long steady-state time stepping.
  • Strong contrasts, such as a sound speed of 3500 m/s against a 1500 m/s background and a density of 3000 kg/m^3 against 1000 kg/m^3, converge, covering bone-water interfaces.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Editorial inference: the same split-preconditioner recipe could apply to other first-order coupled systems, such as visco-elastic or poroelastic wave equations, whenever the differential block is skew-Hermitian and the material block is accretive.
  • Editorial inference: the two shift parameters $a_1$ and $a_2$ could be selected by a minimax or spectral-radius optimization instead of the median heuristic; this is a direct, testable improvement the paper leaves open.
  • Editorial inference: because convergence degrades with contrast, one could expect a practical contrast ceiling per problem size; the paper reports convergence but does not quantify a threshold.
  • Editorial inference: the FFT-based preconditioner ties the method to translation-invariant free-space boundaries; extending to arbitrary Dirichlet or Neumann surfaces would require replacing $(L + I)^{-1}$ with a local sparse representation, sacrificing the matrix-free structure.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper presents an iterative, matrix-free solver for the first-order acoustic Helmholtz system with heterogeneous sound speed, density, and absorption. The method reformulates the Helmholtz equation as a first-order system and applies the universal split-preconditioner of Vettenburg and Vellekoop, yielding an FFT-based iteration whose per-iteration cost is dominated by eight 3D FFTs. The authors prove accretiveness of the continuous and staggered-grid discrete operators, provide Fourier-domain expressions for (L+I)^{-1}, and validate accuracy against Anderson's analytical fluid-sphere solution and a transcranial ultrasound benchmark solved with k-Wave. The paper also claims forward and adjoint capability for inverse problems and emphasizes suitability for large-scale 3D problems with minimal memory overhead.

Significance. If the efficiency claim is substantiated, the paper extends convergent Born-series methods to variable density without matrix factorizations or preprocessing, which is genuinely useful for transcranial ultrasound, full-waveform inversion, and similar large-scale acoustic problems. The accuracy validation against independent benchmarks (Anderson's solution and k-Wave) is credible, and the authors provide a public implementation. The derivation of the discrete skew-Hermitian operator on a staggered grid is a concrete contribution. However, the headline 'fast' and 'large-scale 3D' claims currently rest on per-iteration FFT cost rather than measured total iteration counts or runtimes, and the convergence rate depends on heuristically chosen shifts.

major comments (3)
  1. [Sections 3.2 and 4] The central efficiency claim ('fast', 'suitable for large-scale 3D problems') is not supported by the reported experiments. The 3D benchmarks in Sections 3.1 and 3.2 report only accuracy (relative error vs. the Anderson solution and the k-Wave reference), and Figure 4 shows residual curves only for a 2D disc. Since Section 4 states that each iteration costs eight 3D FFTs, the total runtime is determined by the number of iterations, and no iteration counts, wall-clock times, or memory measurements are reported for the 120x70x70 mm transcranial case. Please report these quantities for the 3D benchmarks, including the iteration count to reach each tolerance, and, if possible, a comparison with the k-Wave reference runtime.
  2. [Section 2.3.2 and Discussion] The convergence rate depends on the heuristic choice of a1 and a2 (Algorithm 1, Eq. (22)), and the Discussion admits that no systematic optimization strategy has been developed. This is not merely a theoretical concern: Figure 4 shows that for 2D disc contrasts the solver needs up to about 2000 iterations to reach a residual of 1e-4, and the corresponding iteration counts for the 3D transcranial case are not given. Please provide a sensitivity study of the iteration count to the choice of a1 and a2 for a representative 3D case, or at least report the chosen values and the resulting iteration counts so the reader can assess the practical impact of this heuristic.
  3. [Section 2.3, Eq. (19)] The displayed adjoint operator appears to be missing a factor 1/rho0 in the second diagonal block: conjugating the second diagonal entry of Eq. (7), (i*omega + gamma)/(rho0*c^2), gives (-i*omega + gamma)/(rho0*c0^2)*(1 + 2*i*omega*alpha0*c0), whereas Eq. (19) writes 1/c0^2 without rho0. As written, Eq. (19) is not the adjoint of Eq. (7), which undermines the stated inverse-problem capability. Please correct the formula or explain any additional scaling used.
minor comments (4)
  1. [Section 2.3.1] The amplitude scaling factor 2c0/dx is stated without derivation; since it enters the source normalization, please provide the derivation or a reference for this factor.
  2. [Eq. (22)] The definition of C mixes a 2x2 matrix with scalar max expressions; please write lambda1 and lambda2 explicitly so that the diagonal scaling is unambiguous.
  3. [Appendix B] The inverse formulas use lambda0, lambda1, a, and b, while the main text uses lambda1, lambda2, a1, and a2; please align the notation across the paper.
  4. [Throughout] There are several typographical errors, including 'choosen', 'rearrangment', and inconsistent capitalization of 'Vettenburg & Vellekoop'; a careful proofread would improve readability.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: derivation is from an independently published split-preconditioner theory and validation is against external analytical and time-domain benchmarks.

full rationale

The paper's central contribution—a matrix-free iterative solver for the first-order Helmholtz system with heterogeneous density—is derived from an independently published universal split-preconditioner (Ref. [2], Vettenburg and Vellekoop), not from a self-citation. The derivation in Section 2.3 establishes accretiveness of the diagonal operator and skew-Hermiticity of the differential operator, then applies the cited preconditioner; the inverse (L+I)^{-1} is computed analytically in the Fourier domain (Appendix B). Accuracy is validated against Anderson's analytical solution (Section 3.1, Figure 2) and against the externally established k-Wave time-domain solver on a transcranial benchmark (Section 3.2, Figure 3); neither benchmark uses fitted parameters from the solver. Free parameters a1, a2, and beta are chosen from material-property medians and a spectral-norm bound, not calibrated to the validation data, and the paper explicitly discloses the absence of a systematic optimization theory (Discussion, page 9). The efficiency claim is a performance risk because iteration counts for 3D cases are not reported, but that is a completeness and measurement concern, not circular reasoning.

Assumptions & free parameters 4 free parameters · 5 assumptions · 0 invented entities

The method introduces no new physical entities; it relies on the universal split-preconditioner theory from prior work, the accretiveness of the acoustic operator, and the discrete skew-Hermitian structure of a staggered k-space gradient. The free parameters are numerical tuning choices, not fitted data parameters.

free parameters (4)
  • a1 = not specified (median of rho0(iw+gamma))
    Complex shift for the velocity block, chosen by heuristic to minimize the norm of V0. Convergence rate depends on it; no optimization is provided.
  • a2 = not specified (median of (iw+gamma)/(rho0*c^2))
    Complex shift for the pressure block, chosen by heuristic. Same impact on convergence rate.
  • beta = not specified, <1
    Scaling factor in C that enforces ||V|| < 1, a necessary condition for the convergence of the split-preconditioned Neumann series.
  • nu = not specified, 0<nu<1
    Relaxation parameter in the fixed-point update (eq. 13); affects iteration stability and speed.
assumptions (5)
  • standard math Universal split-preconditioner convergence theorem (Ref. [2])
    Guarantees Neumann-series convergence for accretive A with ||V||<1; the entire iteration rests on this theorem, cited in Section 2.2.
  • domain assumption Material parameters are positive and the diagonal operator is accretive (rho0>0, c0>0, gamma>0, alpha>=0)
    Invoked in Section 2.3 to show the diagonal part of A has positive real part, making A accretive.
  • standard math Staggered-grid finite-difference operator with shift operators preserves discrete skew-Hermitianness
    Shown in Appendix B.1 via the shift theorem; this is what makes the discrete L accretive and the universal preconditioner applicable.
  • domain assumption Free-space radiation conditions implemented through a diagonal absorbing layer gamma with translation-invariant boundary conditions
    Assumed in Section 2.1 and Discussion; limits the method to open domains and enables FFT-based inversion of (L+I).
  • standard math Equivalence between the first-order system (5) and the second-order Helmholtz equation (3)
    Derived in Appendix A; the solver targets the first-order system under the complex sound speed definition (4).

how reviews work

0 comments
Cite this review

Pith. "Pith review of Iterative Born Solver for the Acoustic Helmholtz Equation with Heterogeneous Sound Speed and Density." pith.science (2026). https://pith.science/paper/NVV6QG5J

@misc{pith2026250716087,
  author       = {Pith},
  title        = {Pith review of: Iterative Born Solver for the Acoustic Helmholtz Equation with Heterogeneous Sound Speed and Density},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/NVV6QG5J}},
  note         = {Machine review of arXiv:2507.16087}
}
read the original abstract

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. We present a fast iterative solver that extends the Convergent Born Series method to handle arbitrary variations in sound speed, density, and absorption simultaneously. Our approach reformulates the Helmholtz equation as a first-order system and applies Vettenburg and Vellekoop's universal split-preconditioner, yielding a matrix-free algorithm that leverages Fast Fourier Transforms for computational efficiency. Unlike existing Born series methods, our solver accommodates heterogeneous density without requiring expensive matrix decompositions or pre-processing steps, making it suitable for large-scale 3D problems with minimal memory overhead. The method provides both forward and adjoint solutions, enabling its application for inverse problems. We validate accuracy through comparison against an analytical solution and demonstrate the solver's practical utility 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.

Figures

Figures reproduced from arXiv: 2507.16087 by the authors.

Figure 1
Figure 1. Setup and results for the heterogeneous experiments. The source is a planar line source at the top of the [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. Comparison of the solutions obtained from the proposed numerical solver and an analytic solution obtained [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗
Figure 3
Figure 3. Comparison of the solutions obtained from a reference solve (k-Wave time domain simulation run to steady [PITH_FULL_IMAGE:figures/full_fig_p008_3.png] view at source ↗
Figures from the paper (1 more)
Figure 4
Figure 4. Figure 4: Normalized residual error as a function of iteration count, illustrating the convergence behavior of the iterative [PITH_FULL_IMAGE:figures/full_fig_p009_4.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

28 extracted references · 23 canonical work pages

  1. [1]

    A convergent Born series for solving the inhomogeneous Helmholtz equation in arbitrarily large media

    Gerwin Osnabrugge, Saroch Leedumrongwatthanakun, and Ivo M Vellekoop. A convergent Born series for solving the inhomogeneous Helmholtz equation in arbitrarily large media. Journal of computational physics, 322:113–124, 2016

    work page 2016

  2. [2]

    A universal matrix-free split preconditioner for the fixed-point iterative solution of non-symmetric linear systems

    Tom Vettenburg and Ivo M Vellekoop. A universal matrix-free split preconditioner for the fixed-point iterative solution of non-symmetric linear systems. arXiv preprint arXiv:2207.14222v2, 2022

    work page Pith review arXiv 2022

  3. [3]

    Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models

    Jean-Francois Aubry, Oscar Bates, Christian Boehm, Kim Butts Pauly, Douglas Christensen, Carlos Cueto, Pierre Gélat, Lluis Guasch, Jiri Jaros, Yun Jing, et al. Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models. The Journal of the Acoustical Society of America, 152(2):1003– 1019, 2022

    work page 2022

  4. [4]

    Waveholtz: Iterative solution of the helmholtz equation via the wave equation

    Daniel Appelo, Fortino Garcia, and Olof Runborg. Waveholtz: Iterative solution of the helmholtz equation via the wave equation. SIAM Journal on Scientific Computing, 42(4):A1950–A1983, 2020. 14

    work page 2020

  5. [5]

    Advances in iterative methods and preconditioners for the Helmholtz equation

    Yogi A Erlangga. Advances in iterative methods and preconditioners for the Helmholtz equation. Archives of Computational Methods in Engineering, 15:37–66, 2008

    work page 2008

  6. [6]

    An improved two-grid preconditioner for the solution of three-dimensional Helmholtz problems in heterogeneous media

    Henri Calandra, Serge Gratton, Xavier Pinel, and Xavier Vasseur. An improved two-grid preconditioner for the solution of three-dimensional Helmholtz problems in heterogeneous media. Numerical Linear Algebra with Applications, 20(4):663–688, 2013

    work page 2013

  7. [7]

    Martin J Gander and Hui Zhang. A class of iterative solvers for the helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized schwarz methods. Siam Review, 61(1):3–76, 2019

    work page 2019

  8. [8]

    A helmholtz equation solver using unsupervised learning: Application to transcranial ultrasound

    Antonio Stanziola, Simon R Arridge, Ben T Cox, and Bradley E Treeby. A helmholtz equation solver using unsupervised learning: Application to transcranial ultrasound. Journal of computational physics, 441:110430, 2021

    work page 2021

Show all 28 references
  1. [9]

    A learned born series for highly-scattering media

    Antonio Stanziola, Simon Arridge, Ben T Cox, and Bradley E Treeby. A learned born series for highly-scattering media. JASA Express Letters, 3(5), 2023

    2023

  2. [10]

    A neural multigrid solver for helmholtz equations with high wavenumber and heterogeneous media

    Chen Cui, Kai Jiang, and Shi Shu. A neural multigrid solver for helmholtz equations with high wavenumber and heterogeneous media. arXiv preprint arXiv:2404.02493, 2024

    2024 arXiv

  3. [11]

    Solving 2-d helmholtz equation in the rectangular, circular, and elliptical domains using neural networks

    D Veerababu and Prasanta K Ghosh. Solving 2-d helmholtz equation in the rectangular, circular, and elliptical domains using neural networks. Journal of Sound and Vibration, page 119022, 2025

    2025

  4. [12]

    Non-invasive transcranial ultrasound therapy based on a 3d ct scan: protocol validation and in vitro results

    Fabrice Marquet, M Pernot, Jean-Francois Aubry, G Montaldo, L Marsac, M Tanter, and M Fink. Non-invasive transcranial ultrasound therapy based on a 3d ct scan: protocol validation and in vitro results. Physics in Medicine & Biology, 54(9):2597, 2009

    2009

  5. [13]

    Homotopy scattering series for seismic forward modelling with variable density and velocity

    Kui Xiang, Kjersti Solberg Eikrem, Morten Jakobsen, and Geir Nævdal. Homotopy scattering series for seismic forward modelling with variable density and velocity. Geophysical Prospecting, 70(1):3–18, 2021

    2021

  6. [14]

    Efficient Helmholtz equation solver for frequency domain waveform inversion based on the decomposition into one-way wave equations

    Rehman Ali, Feiyu Wang, Trevor Mitcham, and Nebojsa Duric. Efficient Helmholtz equation solver for frequency domain waveform inversion based on the decomposition into one-way wave equations. In Medical Imaging 2024: Ultrasonic Imaging and Tomography, volume 12932, pages 262–26...

    2024

  7. [15]

    Iterative solution of the Lippmann–Schwinger equation in strongly scattering acoustic media by randomized construction of preconditioners

    Kjersti Solberg Eikrem, Geir Nævdal, and Morten Jakobsen. Iterative solution of the Lippmann–Schwinger equation in strongly scattering acoustic media by randomized construction of preconditioners. Geophysical Journal International, 224(3):2121–2130, 10 2020

    2020

  8. [16]

    A k-space method for coupled first-order acoustic propagation equations

    Makoto Tabei, T Douglas Mast, and Robert C Waag. A k-space method for coupled first-order acoustic propagation equations. The Journal of the Acoustical Society of America, 111(1):53–63, 2002

    2002

  9. [17]

    An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6):WCC1–WCC26, 2009

    Jean Virieux and Stéphane Operto. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6):WCC1–WCC26, 2009

    2009

  10. [18]

    Physics- informed machine learning

    George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics- informed machine learning. Nature Reviews Physics, 3(6):422–440, 2021

    2021

  11. [19]

    Sound scattering from a fluid sphere

    Victor C Anderson. Sound scattering from a fluid sphere. The Journal of the Acoustical Society of America , 22(4):426–431, 1950

    1950

  12. [20]

    Sound scattering from two concentric fluid spheres

    Jared McNew, Roberto Lavarello, and William D O’Brien. Sound scattering from two concentric fluid spheres. The Journal of the Acoustical Society of America, 125(1):1–4, 2009

    2009

  13. [21]

    Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method

    Bradley E Treeby, Jiri Jaros, Alistair P Rendell, and Benjamin T Cox. Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method. The Journal of the Acoustical Society of America, 131(6):4324–4336, 2012

    2012

  14. [22]

    Computational frameworks for the fast Fourier transform

    Charles Van Loan. Computational frameworks for the fast Fourier transform. SIAM, 1992

    1992

  15. [23]

    Vellekoop

    Swapnil Mache and Ivo M. Vellekoop. Domain decomposition of the modified Born series approach for large-scale wave propagation simulations. arXiv:2410.02395v2, 2024

    2024

  16. [24]

    Ultrasonic neuromodulation by brain stimulation with transcranial ultrasound

    Yusuf Tufail, Anna Yoshihiro, Sandipan Pati, Monica M Li, and William J Tyler. Ultrasonic neuromodulation by brain stimulation with transcranial ultrasound. nature protocols, 6(9):1453–1470, 2011

    2011

  17. [25]

    Trans-spinal focused ultrasound stimulation selectively modulates descending motor pathway

    Evgenii Kim, Jeungeun Kum, and Hyungmin Kim. Trans-spinal focused ultrasound stimulation selectively modulates descending motor pathway. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 30:314–320, 2022

    2022

  18. [26]

    Three-dimensional frequency-domain full-waveform inversion with an iterative solver

    René-Édouard Plessix. Three-dimensional frequency-domain full-waveform inversion with an iterative solver. Geophysics, 74(6):WCC149–WCC157, 2009. 15

    2009

  19. [27]

    Neural operator: Graph kernel network for partial differential equations

    Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485, 2020

    2003 arXiv

  20. [28]

    Beyond regular grids: Fourier-based neural operators on arbitrary domains

    Levi Lingsch, Mike Y Michelis, Emmanuel de Bezenac, Sirani M Perera, Robert K Katzschmann, and Sid- dhartha Mishra. Beyond regular grids: Fourier-based neural operators on arbitrary domains. arXiv preprint arXiv:2305.19663, 2023. 16

    2023

Pith tools

Reviewed August 6, 2026 · model on record in the stance chip above.