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

This paper presents an $O(M)$ algorithm that efficiently and accurately evaluates the azimuthal Fourier modes of the 3D Helmholtz Green's function and their derivatives for any real wavenumber by combining contour deformation with stable recurrence relations, enabling high-performance axisymmetric acoustic scattering simulations.

The Gist

Imagine you are trying to listen to a complex symphony played by a giant, rotating drum. The sound isn't just a single note; it's a mix of thousands of different "modes" or layers of vibration, each spinning at a different speed. In the world of physics and engineering, calculating how sound (or light, or radio waves) bounces off a round object is like trying to figure out exactly what each of those thousands of layers sounds like.

This paper introduces a new, super-fast way to calculate those layers, along with how they change (their "derivatives"), without getting bogged down by the math usually required.

Here is the breakdown of what the authors did, using everyday analogies:

1. The Problem: The "Oscillating Wave" Nightmare

Usually, to figure out how a wave behaves around a round object, you have to do a massive amount of math involving integrals (summing up tiny pieces).

• The Catch: If you want to calculate many layers (modes), the old methods get slower and slower. It's like trying to count every grain of sand on a beach one by one.
• The Difficulty: Sometimes the waves are huge, and sometimes they are so tiny they are practically invisible (exponentially small). Standard math tools often lose accuracy when the numbers get that small, like trying to weigh a feather on a scale meant for elephants.
• The Geometry: The math gets even messier when the source of the sound and the target are very close together, creating a "near-singular" situation where the numbers blow up.

2. The Solution: A Two-Step "Magic Trick"

The authors created an algorithm that solves this in linear time ($O(M)$). This means if you double the number of layers you want to calculate, the time it takes only doubles, rather than exploding into a massive calculation.

They achieved this by combining two clever strategies:

Strategy A: The "Steep Slide" (Contour Deformation)

Imagine you are trying to walk across a bumpy, oscillating field to get from point A to point B. Walking straight across is exhausting because you have to step up and down thousands of times.

• The Trick: Instead of walking on the surface, the authors found a secret "slide" (a path in the complex plane) that goes underneath the bumps. On this slide, the wavy, bumpy terrain turns into a smooth, straight slope that goes downhill.
• The Benefit: You can slide down this path very quickly and accurately, regardless of how wavy the original terrain was. They use this only for a few "boundary" layers (the very first and very last ones you need).

Strategy B: The "Domino Chain" (Recurrence Relations)

Once they have the first and last layers calculated using the "slide," they don't calculate the middle ones one by one.

• The Trick: They realized the layers are connected like a chain of dominoes. If you know the first and last domino, you can figure out all the ones in the middle by solving a giant, structured puzzle (a linear system).
• The Benefit: This avoids the instability of trying to push the dominoes from just one end (which often causes the chain to fall over or become inaccurate). By pinning both ends, the whole chain stands up perfectly.

3. Handling the "Tiny" and the "Messy"

• The Tiny Layers: In the "decay regime," the layers get so small they vanish into the noise. The authors use a special technique (similar to Miller's algorithm) where they pretend the very far-away layers are zero and work backward. This ensures that even the tiniest, almost invisible layers are calculated with high precision, not lost to rounding errors.
• The Messy Neighbors: When the source and target are right next to each other, the math gets "singular" (blows up). The authors use a special type of calculator (Generalized Gaussian Quadrature) designed specifically to handle these sharp spikes without losing accuracy.

4. The "Bonus" Feature: Derivatives

In physics, you often need not just the sound level, but how fast it's changing (first derivative) or how the rate of change is changing (second derivative).

• The Paper's Claim: Usually, calculating these extra details takes a lot of extra work. The authors show that once you have the main layers, you can get all these extra details using stable "recurrence" formulas.
• The Cost: It only adds a tiny, constant amount of time (about 30% more) to get all these extra details. It's like getting a full report card (grades, attendance, and behavior) for the same price as just getting the grades.

5. The Result: Speed and Independence

The most impressive claim is that this method is independent of the wavenumber (how fast the wave vibrates) and the distance between the source and target.

• Analogy: Imagine a delivery service. Usually, if the package is heavy (high frequency) or the distance is tricky (close proximity), the delivery takes longer. This new algorithm delivers the package in the exact same amount of time, whether it's a feather or a boulder, and whether it's next door or across town.

Summary

The paper presents a mathematical "shortcut" that allows computers to calculate how waves interact with round objects. By using a "slide" to get the start and end points and a "domino chain" to fill in the middle, they can calculate thousands of wave layers and their changes in a flash. This makes it possible to simulate complex acoustic and electromagnetic scattering (like radar or sound bouncing off a submarine) much faster and more accurately than before, without the computer getting confused by tiny numbers or close distances.

Technical Summary

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

Problem Statement
The paper addresses the numerical evaluation of the azimuthal Fourier modes, $G_{k,m}$, of the three-dimensional Helmholtz Green's function, along with their first- and second-order derivatives with respect to cylindrical source and target coordinates. These modes arise naturally in boundary integral equation (BIE) methods for problems with rotational symmetry, such as acoustic scattering and electromagnetic metasurface design.

Evaluating these quantities poses significant challenges:

1. Oscillation: The integrands oscillate at frequencies determined by both the mode number $m$ and the wavenumber $k$, requiring a number of quadrature nodes that grows with both parameters in naive approaches.
2. Near-Singularity: When source and target points are close, the integrand becomes near-singular, causing standard smooth quadrature to fail.
3. Exponential Decay: For large $m$, the magnitude of $|G_{k,m}|$ decays exponentially. Direct quadrature loses relative accuracy in this regime due to catastrophic cancellation between the large integrand and the small result.
4. Derivatives: BIE applications often require derivatives, which introduce even stronger near-singular factors.
5. Computational Cost: Existing methods either scale as $O(M^2)$ to compute all modes $m=0, \dots, M$ or have costs dependent on the wavenumber $k$. There is a lack of an $O(M)$ algorithm that is independent of $k$ and maintains high relative accuracy across all regimes, including the decay regime.

Methodology
The proposed algorithm achieves $O(M)$ complexity with cost independent of the wavenumber $k$ and source-target separation by combining contour deformation with a boundary-value formulation of recurrence relations. The method operates in three distinct regimes based on the geometric parameter $\alpha$ and the mode number $m$:

1. Near-Axis Regime ($\alpha \ll 1$): Handled via precomputed Taylor-expansion tables.
2. Non-Decay and Decay Regimes:
   • Contour Deformation: For a small number of boundary modes (specifically $m=0, 1$ and $m=M-1, M$), the algorithm employs contour deformation in the complex plane. The integration path is deformed into steepest descent contours ($\gamma_1, \gamma_2$) and a Bernstein ellipse arc ($E_c$). This renders the spherical wave term non-oscillatory and exponentially decaying, allowing for accurate evaluation with a bounded number of quadrature nodes independent of $k$.
   • Handling Derivatives: The same deformed contours are used for derivative integrals. To handle the stronger near-singular factors in derivative integrands (particularly near $z=1$), the method utilizes precomputed Generalized Gaussian Quadratures rather than monomial recurrences, avoiding separate recurrences for each derivative type.
   • Recurrence Relations: The interior modes are not computed via direct integration. Instead, the algorithm utilizes a five-term recurrence relation satisfied by $G_{k,m}$ (and associated recurrences for derivatives). Rather than using this recurrence as a forward or backward iteration (which is unstable in certain regimes), it is cast as a banded boundary-value problem. The boundary values are supplied by the contour deformation method, and the interior modes are recovered by solving a pentadiagonal linear system.
   • Decay Regime Strategy: When $m$ exceeds the transition point $m^*$ where modes decay exponentially, the boundary values $G_{M-1}$ and $G_M$ become smaller than machine epsilon. In this case, the algorithm employs a Miller-type strategy: it sets the boundary values at a sufficiently high index $M'$ (where the values are negligible) to zero and solves the linear system for the interior modes. This preserves relative accuracy in the exponentially decaying tail.
   • Derivative Recovery: Once $G_{k,m}$ is known, the first- and second-order derivatives are obtained via stable upward or downward recurrences (depending on the regime) and pointwise identities. Cancellation-free combinations of derivatives are computed directly to ensure stability in the near-singular regime.

Key Contributions

• $O(M)$ Complexity: The algorithm computes all modes $m=0, \dots, M$ and their derivatives in linear time with respect to $M$.
• Wavenumber Independence: The computational cost is independent of the wavenumber $k$ and the source-target separation distance.
• High Relative Accuracy: The method maintains high relative accuracy (close to double precision) even for modes whose magnitudes are exponentially small, a regime where direct quadrature fails.
• Derivative Efficiency: First- and second-order derivatives are computed with only a small constant factor increase in cost relative to evaluating $G_{k,m}$ alone, utilizing stable recurrences and cancellation-free formulations.
• Robustness: The combination of contour deformation for boundary values and the pentadiagonal boundary-value formulation for interior modes avoids the instabilities associated with pure forward or backward recursion.

Numerical Results
Numerical experiments demonstrate:

• Accuracy: Relative errors remain uniform across modes $m=0, \dots, M$ (up to $M=3000$) at the level of $10^{-11}$ to $10^{-12}$, even in high-frequency and near-singular geometries. In the decay regime, relative accuracy is preserved over more than 15 orders of magnitude of decay.
• Scaling: Runtime scales linearly with $M$. The inclusion of first-order derivatives adds negligible overhead, while second-order derivatives increase runtime by approximately 30% due to additional contour evaluations and solves.
• Independence: Runtimes are essentially independent of the wavenumber $k$ (tested from $k=10$ to $5000$) and the source-target separation (tested from well-separated to near-singular).
• Application: The evaluator was integrated into modal BIE solvers for acoustic scattering from bodies of revolution (using CFIE and Müller's formulations). The kernel evaluation cost became negligible compared to the dense linear algebra costs, confirming that the evaluator successfully removes the kernel evaluation bottleneck.

Significance
The paper claims that this algorithm provides the first $O(M)$ evaluation of all Helmholtz modal Green's functions and their derivatives with cost independent of the wavenumber. By decoupling the kernel evaluation cost from the wavenumber and source-target geometry, the method shifts the dominant computational cost in axisymmetric BIE solvers to the dense per-mode linear algebra. This enables the application of fast direct solvers based on hierarchical compression to axisymmetric problems, potentially allowing for substantially larger frequency-domain scattering calculations while retaining the accuracy and geometric flexibility of boundary integral formulations. The authors note that while the componentwise accuracy of the pentadiagonal solve is observed experimentally, a rigorous relative-error analysis for this specific formulation remains an open direction for future work.