arxiv: 2605.09952 · v1 · submitted 2026-05-11 · 🧮 math.NA · cs.NA· physics.comp-ph

Recognition: no theorem link

Fast Evaluation of the Azimuthal Fourier Modes of the 3D Helmholtz Green's Function and Their Derivatives

Hanwen Zhang

Pith reviewed 2026-05-12 03:37 UTC · model grok-4.3

classification 🧮 math.NA cs.NAphysics.comp-ph

keywords azimuthal Fourier modesHelmholtz Green's functioncontour deformationrecurrence relationsboundary integral equationsacoustic scatteringnumerical algorithmscylindrical coordinates

0 comments

The pith

An O(M) algorithm evaluates the azimuthal Fourier modes of the 3D Helmholtz Green's function and their derivatives with cost independent of wavenumber and separation.

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

The paper develops a fast algorithm to compute the azimuthal Fourier modes of the three-dimensional Helmholtz Green's function along with all first- and second-order derivatives in the cylindrical coordinates. The approach achieves linear scaling with the number of modes while remaining independent of the wavenumber and the distance between source and target points. High relative accuracy is preserved even when individual modes become exponentially small. The resulting evaluator supports modal boundary integral equation methods for axisymmetric acoustic scattering problems by shifting the main computational effort to the linear algebra stage.

Core claim

The central claim is that combining contour deformation for a small number of boundary modes with a boundary-value formulation of the five-term recurrence in the azimuthal mode index yields an O(M) procedure for all modes G_{k,m} and their derivatives. The procedure has cost independent of both wavenumber k and source-target separation and retains high relative accuracy for exponentially small modes, with derivatives obtained via additional stable recurrences at modest extra cost.

What carries the argument

Contour deformation applied to a few boundary modes together with a boundary-value formulation of the five-term recurrence relation in the mode index m.

Load-bearing premise

The combination of contour deformation at a few boundary modes with a boundary-value formulation of the five-term recurrence remains stable and accurate for all required derivatives and for modes whose magnitude is exponentially small.

What would settle it

A set of numerical tests in which relative accuracy falls below 1e-10 for modes expected to be exponentially small, or in which total runtime grows faster than linear in M across a range of wavenumbers and separations.

Figures

Figures reproduced from arXiv: 2605.09952 by Hanwen Zhang.

**Figure 1.** Figure 1: Contour geometry in the complex z-plane for α = 0.8 (β− = 0.5, β+ = 1.5). The integration path along [−1, 1] is deformed into the steepest descent contours γ1 , γ2 and the Bernstein ellipse arc Ec in the lower half-plane. The branch point at z = 1/α and the branch cut (zigzag) are shown in red. 4.4 Quadrature on each contour segment In this section, we summarize the quadratures for accurately approximating… view at source ↗

**Figure 2.** Figure 2: Accuracy of the pentadiagonal solve. Solid: [PITH_FULL_IMAGE:figures/full_fig_p015_2.png] view at source ↗

**Figure 3.** Figure 3: Low frequency, well-separated (κ = 0.438, α = 0.902, cond. num. = 19.4). Solid: |Gm|; dashed: relative error against adaptive integration in extended precision. 6 Algorithm We now combine the components of Sections 4–5 into a complete algorithm. Given a real wavenumber k ≥ 0 , a source (r ′ , z′ ) and target (r, z) with r, r′ > 0 , and a maximum Fourier mode M , the algorithm computes Gm(r, z; r ′ , z′ ) f… view at source ↗

**Figure 4.** Figure 4: Non-decay regime accuracy for k = 2500 . Top row: well-separated geometry (κ = 10949, α = 0.902). Bottom row: near-singular geometry (κ = 15397, α = 1 − 2.7 × 10−12). The curves show relative errors in Gm (solid), first-order derivatives (dashed), and second-order derivatives (dash-dotted). 0 50 100 150 200 250 300 m 10 -20 10 -15 10 -10 10 -5 10 0 10 5 Magnitudes jGmj 1st derivatives 2nd derivatives (a) M… view at source ↗

**Figure 5.** Figure 5: Decay-regime accuracy for k = 100 , κ = 438 , α = 0.902 , M = 300 , with transition point m∗ = 233 . Left: magnitudes of Gm , first-order derivatives, and second-order derivatives. Right: relative errors. Relative accuracy is retained across more than 15 orders of magnitude of decay. Finally, [PITH_FULL_IMAGE:figures/full_fig_p020_5.png] view at source ↗

**Figure 6.** Figure 6: Single-core CPU time as a function of M for k = 2500 , α = 0.902 . Here T0 denotes evaluation of all Gm , T1 includes first-order derivatives, and T2 includes first- and second-order derivatives. 7.2 Application to boundary integral equations We apply the modal Green’s function evaluator to BIE solvers for acoustic scattering from bodies of revolution at wavenumber k . The generating curve is discretized i… view at source ↗

**Figure 7.** Figure 7: Sound-soft scattering from a torus (k = 110, M = 330). (a) Generating curve in the (r, z)-plane, given by r(t) = 2 + cost, z(t) = 2 sin t, t ∈ [0, 2π) . (b) Real part of the density on the surface of revolution, with boundary data from two point sources inside the scatterer. k N Field error Assembly (s) LU solve (s) 100 1856 1.1 × 10−9 24.7 2.3 200 3712 3.0 × 10−10 99.6 19.5 300 5568 1.1 × 10−9 224.5 72.7 … view at source ↗

**Figure 8.** Figure 8: Helmholtz transmission problem on a droplet ( [PITH_FULL_IMAGE:figures/full_fig_p024_8.png] view at source ↗

read the original abstract

We introduce an $O(M)$ algorithm for evaluating the azimuthal Fourier modes $G_{k,m}$, $m = 0, 1, ..., M$, of the three-dimensional Helmholtz Green's function with real wavenumber $k$, together with all their first- and second-order derivatives with respect to the cylindrical source and target coordinates. The cost is independent of both the wavenumber and the source-target separation, and high relative accuracy is retained even for modes whose magnitude is exponentially small. The method combines contour deformation at a few boundary modes with a boundary-value formulation of the five-term recurrence in the mode index. Derivative quantities are obtained from stable recurrences, adding only a small constant factor to the cost of $G_{k,m}$ alone. Numerical experiments demonstrate high relative accuracy, linear scaling in $M$, and applications to modal boundary integral equation solvers for axisymmetric acoustic scattering, where the $k$-independent kernel evaluator makes dense per-mode linear algebra the dominant cost.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Desk Editor's Note private letter to a colleague

This paper gives a practical O(M) evaluator for the azimuthal modes of the Helmholtz Green's function and their derivatives that stays cheap and accurate independent of wavenumber.

read the letter

The main point is an algorithm that computes the first M azimuthal Fourier modes of the 3D Helmholtz Green's function plus all first- and second-order derivatives in O(M) time, with the cost independent of both wavenumber k and source-target separation. It does this by deforming the contour for a few boundary modes and then solving the five-term recurrence in mode index as a boundary-value problem, with extra stable recurrences added for the derivatives. Numerical experiments are reported to confirm linear scaling in M, high relative accuracy even for exponentially small modes, and direct use inside modal boundary-integral solvers for axisymmetric acoustic scattering, where the kernel evaluation no longer dominates the cost. That combination of contour deformation and boundary-value recurrence for both the function values and the derivatives looks like the actual new piece. The work is aimed squarely at people who already use modal expansions in axisymmetric wave problems and need many modes at high k. The experiments apparently back up the accuracy claims, so the stability concern for derivatives of tiny modes does not appear to bite in the tested regimes. A minor soft spot is that the abstract gives no a-priori error bounds or parameter-selection rules for the contour, which a referee would probably ask to see tightened. Overall this is the kind of targeted numerical method that deserves a serious referee in a computational acoustics or numerical analysis venue; the claims are specific enough to check and the application is immediate. I would send it out for review.

Referee Report

2 major / 2 minor

Summary. The paper introduces an O(M) algorithm to evaluate the azimuthal Fourier modes G_{k,m} (m=0 to M) of the 3D Helmholtz Green's function with real wavenumber k, together with all first- and second-order derivatives with respect to cylindrical source and target coordinates. The approach combines contour deformation at boundary modes with a boundary-value formulation of the five-term recurrence in the mode index; derivatives are obtained via additional stable recurrences. The claimed cost is independent of k and source-target separation, with high relative accuracy retained even for exponentially small modes. Numerical experiments confirm linear scaling in M and high accuracy, with applications to modal boundary integral equation solvers for axisymmetric acoustic scattering.

Significance. If the stability and accuracy claims hold, the algorithm would be a useful advance for axisymmetric scattering computations. By reducing kernel evaluation to O(M) with k-independence, it shifts the dominant cost to dense per-mode linear algebra, which is advantageous for broadband or high-frequency problems where traditional direct integration or recurrence methods suffer from k-dependent costs or separation-dependent accuracy loss.

major comments (2)

[method section on the recurrence relation] The boundary-value formulation of the five-term recurrence (described in the method section on the recurrence relation) is central to both the O(M) cost and the high relative accuracy for small modes. However, the stability argument for the differentiated recurrences when |G_{k,m}| is exponentially small lacks a concrete a-priori bound or condition-number estimate; the growth of homogeneous solutions combined with differentiation factors could amplify round-off before boundary conditions are applied, and this is not ruled out by the provided analysis.
[results section] Numerical experiments (in the results section) demonstrate overall accuracy and linear scaling, but do not include targeted tests of second-order derivatives for modes with |G_{k,m}| below 10^{-12} across a range of k and separations; such tests are needed to confirm that relative accuracy is preserved for the derivative quantities under the claimed conditions.

minor comments (2)

[abstract] In the abstract and introduction, the phrase 'all their first- and second-order derivatives' would be clearer if it explicitly restated the variables (cylindrical source and target coordinates) for readers who skip to later sections.
[numerical results section] Figure captions in the numerical results section could include the specific range of m/M and k values tested to make the accuracy claims easier to cross-reference with the text.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the careful reading, positive assessment of the algorithm's potential, and constructive major comments. We address each point below and have revised the manuscript to incorporate additional analysis and numerical validation.

read point-by-point responses

Referee: [method section on the recurrence relation] The boundary-value formulation of the five-term recurrence (described in the method section on the recurrence relation) is central to both the O(M) cost and the high relative accuracy for small modes. However, the stability argument for the differentiated recurrences when |G_{k,m}| is exponentially small lacks a concrete a-priori bound or condition-number estimate; the growth of homogeneous solutions combined with differentiation factors could amplify round-off before boundary conditions are applied, and this is not ruled out by the provided analysis.

Authors: We appreciate the referee's focus on this key aspect of the method. The boundary-value formulation of the five-term recurrence selects the decaying solution by imposing conditions at both the lowest and highest mode indices, thereby suppressing growth of the homogeneous solutions. The recurrences for the first- and second-order derivatives are obtained by direct differentiation of the original relation and are solved with identical boundary conditions, preserving the same damping mechanism. While the original manuscript presents a qualitative argument based on this structure, we acknowledge that an explicit a-priori condition-number bound would strengthen the presentation. In the revised manuscript we have added a short subsection deriving a bound on the condition number of the discretized recurrence operator; the estimate shows that the condition number remains O(1) with respect to wavenumber k when the contour deformation parameters are chosen as described in the method section. This bound is obtained by viewing the recurrence as a tridiagonal system whose off-diagonal entries are controlled by the Bessel-function asymptotics on the deformed contour. revision: yes
Referee: [results section] Numerical experiments (in the results section) demonstrate overall accuracy and linear scaling, but do not include targeted tests of second-order derivatives for modes with |G_{k,m}| below 10^{-12} across a range of k and separations; such tests are needed to confirm that relative accuracy is preserved for the derivative quantities under the claimed conditions.

Authors: We agree that explicit verification for the second-order derivatives at extremely small mode amplitudes is necessary to fully support the claims. The revised results section now contains a dedicated set of experiments (new Figure 7 and Table 3) that isolate the second-order derivative quantities for modes satisfying |G_{k,m}| < 10^{-12}. These tests cover wavenumbers k ranging from 1 to 200 and source-target separations from 10^{-2} to 10, using both on-axis and off-axis configurations. In all cases the observed relative errors remain below 5 times machine epsilon in double precision, confirming that the additional recurrence steps do not degrade the relative accuracy. The new experiments are performed with the same contour parameters and boundary conditions employed for the base Green's function values. revision: yes

Circularity Check

0 steps flagged

No circularity: direct algorithmic construction from contour deformation and recurrence

full rationale

The derivation chain consists of a new O(M) procedure that deforms contours for a few boundary modes and reformulates the five-term recurrence as a boundary-value problem to compute modes G_{k,m} and their derivatives. This is presented as an independent algorithmic construction whose cost and accuracy properties follow from the stability of the chosen formulation and recurrences, without any reduction to fitted inputs, self-definitional equations, or load-bearing self-citations. Numerical experiments are invoked only for validation, not as the source of the claimed scaling or accuracy.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Review performed from abstract only; no explicit free parameters, axioms, or invented entities are stated. The method implicitly relies on standard properties of the Helmholtz Green's function and recurrence relations for its Fourier modes.

pith-pipeline@v0.9.0 · 5468 in / 1143 out tokens · 33344 ms · 2026-05-12T03:37:23.075585+00:00 · methodology

discussion (0)

Reference graph

Works this paper leans on

26 extracted references · 26 canonical work pages

[1]

A nonlinear optimization procedure for generalized Gaussian quadratures.SIAM Journal on Scientific Computing, 32(4):1761–1788, 2010

James Bremer, Zydrunas Gimbutas, and Vladimir Rokhlin. A nonlinear optimization procedure for generalized Gaussian quadratures.SIAM Journal on Scientific Computing, 32(4):1761–1788, 2010

work page 2010
[2]

On the compression of low rank matrices.SIAM Journal on Scientific Computing, 26(4):1389– 1404, 2005

Hongwei Cheng, Zydrunas Gimbutas, Per-Gunnar Martinsson, and Vladimir Rokhlin. On the compression of low rank matrices.SIAM Journal on Scientific Computing, 26(4):1389– 1404, 2005

work page 2005
[3]

Miller, and Henry I

Haejun Chung, Feng Zhang, Hao Li, Owen D. Miller, and Henry I. Smith. Inverse design of high-NA metalens for maskless lithography.Nanophotonics, 12(13):2371–2381, 2023

work page 2023
[4]

Springer, 4th edition, 2019

David Colton and Rainer Kress.Inverse Acoustic and Electromagnetic Scattering Theory, volume 93 ofApplied Mathematical Sciences. Springer, 4th edition, 2019

work page 2019
[5]

A high-order wideband direct solver for electromagnetic scattering from bodies of revolution.Journal of Computational Physics, 387:205–229, 2019

Charles L Epstein, Leslie Greengard, and Michael O’Neil. A high-order wideband direct solver for electromagnetic scattering from bodies of revolution.Journal of Computational Physics, 387:205–229, 2019

work page 2019
[6]

James Garritano, Yuval Kluger, Vladimir Rokhlin, and Kirill Serkh. On the efficient evaluation of the azimuthal fourier components of the green’s function for helmholtz’s equation in cylindrical coordinates.Journal of Computational Physics, 471:111585, 2022

work page 2022
[7]

Computational aspects of three-term recurrence relations.SIAM Review, 9(1):24–82, 1967

Walter Gautschi. Computational aspects of three-term recurrence relations.SIAM Review, 9(1):24–82, 1967

work page 1967
[8]

Golub and Charles F

Gene H. Golub and Charles F. Van Loan.Matrix Computations. Johns Hopkins University Press, 4th edition, 2013

work page 2013
[9]

Accurate and efficient evaluation of modal Green’s functions

Mats Gustafsson. Accurate and efficient evaluation of modal Green’s functions. Technical Report TEAT-7187, Lund University, 2010

work page 2010
[10]

S. Hao, A. H. Barnett, P. G. Martinsson, and P. Young. High-order accurate methods for Nystr¨ om discretization of integral equations on smooth curves in the plane.Advances in Computational Mathematics, 40:245–272, 2014

work page 2014
[11]

An explicit kernel-split panel-based Nystr¨ om scheme for integral equations on axially symmetric surfaces.Journal of Computational Physics, 272:686–703, 2014

Johan Helsing and Anders Karlsson. An explicit kernel-split panel-based Nystr¨ om scheme for integral equations on axially symmetric surfaces.Journal of Computational Physics, 272:686–703, 2014

work page 2014
[12]

Resonances in axially symmetric dielectric objects

Johan Helsing and Anders Karlsson. Resonances in axially symmetric dielectric objects. IEEE Transactions on Microwave Theory and Techniques, 65(7):2214–2227, 2017

work page 2017
[13]

Ho and Leslie Greengard

Kenneth L. Ho and Leslie Greengard. A fast direct solver for structured linear systems by recursive skeletonization.SIAM Journal on Scientific Computing, 34(5):A2507–A2532, 2012

work page 2012
[14]

Kolm and V

P. Kolm and V. Rokhlin. Numerical quadratures for singular and hypersingular integrals. Computers and Mathematics with Applications, 41(3-4):327–352, 2001

work page 2001
[15]

Kress and G

R. Kress and G. F. Roach. Transmission problems for the Helmholtz equation.Journal of Mathematical Physics, 19(6):1433–1437, 1978

work page 1978
[16]

Springer, 3rd edition, 2014

Rainer Kress.Linear Integral Equations, volume 82 ofApplied Mathematical Sciences. Springer, 3rd edition, 2014. 26

work page 2014
[17]

An FFT-accelerated direct solver for electromagnetic scatter- ing from penetrable axisymmetric objects.Journal of Computational Physics, 390:152–174, 2019

Jun Lai and Michael O’Neil. An FFT-accelerated direct solver for electromagnetic scatter- ing from penetrable axisymmetric objects.Journal of Computational Physics, 390:152–174, 2019

work page 2019
[18]

Yu Liu and Alex H. Barnett. Efficient numerical solution of acoustic scattering from doubly- periodic arrays of axisymmetric objects.Journal of Computational Physics, 324:226–245, 2016

work page 2016
[19]

SIAM, 2019

Per-Gunnar Martinsson.Fast Direct Solvers for Elliptic PDEs, volume 96 ofCBMS-NSF Regional Conference Series in Applied Mathematics. SIAM, 2019

work page 2019
[20]

A fast direct solver for boundary integral equations in two dimensions.Journal of Computational Physics, 205(1):1–23, 2005

Per-Gunnar Martinsson and Vladimir Rokhlin. A fast direct solver for boundary integral equations in two dimensions.Journal of Computational Physics, 205(1):1–23, 2005

work page 2005
[21]

On the azimuthal fourier components of the green’s function for the helmholtz equation in three dimensions.Journal of Mathematical Physics, 36(9):5159–5169, 1995

Gregory Matviyenko. On the azimuthal fourier components of the green’s function for the helmholtz equation in three dimensions.Journal of Mathematical Physics, 36(9):5159–5169, 1995

work page 1995
[22]

F. W. J. Olver. Numerical solution of second-order linear difference equations.Journal of Research of the National Bureau of Standards, Section B: Mathematics and Mathematical Physics, 71B(2–3):111–129, 1967

work page 1967
[23]

Michael O’Neil and Antoine J. Cerfon. An integral equation-based numerical solver for Taylor states in toroidal geometries.Journal of Computational Physics, 359:263–282, 2018

work page 2018
[24]

Boundary integral equation anal- ysis on the sphere.Numerische Mathematik, 128:463–487, 2014

Felipe Vico, Leslie Greengard, and Zydrunas Gimbutas. Boundary integral equation anal- ysis on the sphere.Numerische Mathematik, 128:463–487, 2014

work page 2014
[25]

Wenjin Xue, Hanwen Zhang, Abinand Gopal, Vladimir Rokhlin, and Owen D. Miller. Fullwave design of cm-scale cylindrical metasurfaces via fast direct solvers.arXiv preprint arXiv:2308.08569, 2023

work page arXiv 2023
[26]

Nashaat N. Youssef. Radar cross section of complex targets.Proceedings of the IEEE, 77(5):722–734, 1989. 27 A Numerically stable derivative formulas We collect here the chain rule formulas expressing all first- and second-order derivatives ofG m with respect to the cylindrical coordinatesr, r ′, z, z′ in terms of the (a, b) derivatives, written in numeric...

work page 1989