pith. sign in

arxiv: 2604.07170 · v1 · submitted 2026-04-08 · 🧮 math.NA · cs.NA

A spectral method for the rapid evaluation of hyperbolic potentials in two dimensions using windowed Fourier projection

Nour G. Al Hassanieh , Leslie Greengard , Alex H. Barnett This is my paper

Pith reviewed 2026-05-10 16:56 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords wave equationhyperbolic potentialsfast algorithmsFourier projectionsum-of-exponentialsnumerical methodsscattering solvers
0
0 comments X

The pith

A splitting algorithm evaluates 2D wave potentials from point sources in quasi-linear time by separating local, near-history, and far-history contributions.

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

The paper presents a fast method to compute the non-smooth solution of the free-space 2D scalar wave equation when many point sources each carry a high-frequency band-limited time signature. Direct evaluation costs O(M squared N_t squared) for M sources and N_t time steps because of the full history dependence. The method divides the solution at each time into a non-smooth local part handled by direct quadrature, a smooth near-history part treated with truncated windowed Fourier projection, and a very smooth far-history part captured by a large-time sum-of-exponentials approximation of the wave kernel. This decomposition yields quasi-linear scaling overall. Tests on problems with a million sources and targets inside a 300-by-300 wavelength domain reach six-digit accuracy at five orders of magnitude less cost than direct summation.

Core claim

The central claim is that the free-space 2D wave kernel solution with many band-limited point sources can be evaluated in quasi-linear time by splitting the history at each observation time into a local non-smooth piece, a controlled-oscillation near-history piece amenable to truncated windowed Fourier projection, and a far-history piece that admits an accurate large-time sum-of-exponentials representation thanks to the weak Huygens principle.

What carries the argument

The decomposition of the solution at each time into local, near-history, and far-history components, with TK-WFP applied to the near part and a large-time sum-of-exponentials approximation applied to the far part of the wave kernel.

If this is right

  • Enables practical time-domain scattering solvers that rely on layer potentials for large numbers of sources and long simulation intervals.
  • Achieves five orders of magnitude speedup while retaining six-digit accuracy on domains hundreds of wavelengths across.
  • Allows the near-history window to remain of fixed length O(1) domain traversals, keeping the Fourier-projection cost independent of total simulation length.
  • The local part remains a small, fixed-cost direct quadrature independent of total history length.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same history-splitting idea could be tested on other hyperbolic equations whose kernels satisfy a weak Huygens principle.
  • Refining the sum-of-exponentials fit for even longer times might support simulations over thousands of domain traversals without restart.
  • Coupling the method to an adaptive choice of near-history length could balance accuracy and cost when source frequencies vary.

Load-bearing premise

The large-time sum-of-exponentials approximation must represent the far history with small enough error that inaccuracies do not accumulate noticeably over many time steps.

What would settle it

Run the algorithm on a sequence of increasing total simulation times while comparing the output error against a direct quadrature reference on the same sources; error that remains bounded rather than growing linearly or faster with time steps would support the claim.

Figures

Figures reproduced from arXiv: 2604.07170 by Alex H. Barnett, Leslie Greengard, Nour G. Al Hassanieh.

Figure 2
Figure 2. Figure 2 [PITH_FULL_IMAGE:figures/full_fig_p007_2.png] view at source ↗
Figure 4
Figure 4. Figure 4 [PITH_FULL_IMAGE:figures/full_fig_p011_4.png] view at source ↗
Figure 7
Figure 7. Figure 7 [PITH_FULL_IMAGE:figures/full_fig_p018_7.png] view at source ↗
Figure 7
Figure 7. Figure 7 [PITH_FULL_IMAGE:figures/full_fig_p019_7.png] view at source ↗
Figure 7
Figure 7. Figure 7 [PITH_FULL_IMAGE:figures/full_fig_p021_7.png] view at source ↗
Figure 7
Figure 7. Figure 7 [PITH_FULL_IMAGE:figures/full_fig_p022_7.png] view at source ↗
read the original abstract

We present a fast algorithm for evaluating the (non-smooth) solution of the free-space two-dimensional (2D) scalar wave equation with many point sources, each with a high-frequency band-limited time signature. Such an algorithm is key to an efficient time-domain scattering solver using spatially-discretized hyperbolic layer potentials. Given $M$ sources/targets and $N_t$ time steps, direct evaluation costs $O(M^2N_t^2)$, due to the history dependence. We develop a quasi-linear scaling algorithm that splits the solution at a given time into (a) a non-smooth time-local part, (b) a (smooth) near history involving sources up to ${\mathcal O}(1)$ domain traversal times into the past, plus (c) a (very smooth) far history comprising all waves emitted before the near history. The local part is computed directly via high-order quadrature. A naive spatial Fourier transform for (b) plus (c) would be both slowly converging and arbitrarily oscillatory as time progresses. Yet in (b) the oscillations are controlled, so we use the recent truncated windowed Fourier projection (TK-WFP) method to give rapid convergence. For (c) -- present due to the weak Huygens' principle -- we exploit a new large-time sum-of-exponentials approximation of the free-space wave kernel. Numerical examples with up to a million sources and targets, a domain of $300\times 300$ wavelengths, and 6-digit accuracy, show an acceleration of five orders of magnitude relative to direct evaluation.

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.

Referee Report

1 major / 2 minor

Summary. The manuscript presents a fast algorithm for evaluating the non-smooth solution of the free-space 2D scalar wave equation with many high-frequency point sources. It achieves quasi-linear scaling in M sources/targets and Nt time steps by splitting the solution at each time into a non-smooth local part (direct high-order quadrature), a smooth near-history part (truncated windowed Fourier projection, TK-WFP), and a very smooth far-history part (new large-time sum-of-exponentials approximation to the wave kernel, exploiting the weak Huygens principle). Numerical examples with up to 10^6 sources/targets, domains of 300x300 wavelengths, and 6-digit accuracy report five orders of magnitude speedup over direct O(M^2 Nt^2) evaluation.

Significance. If the far-history approximation error remains controlled without accumulation, the method would provide a practical acceleration for time-domain scattering solvers based on hyperbolic layer potentials, enabling large-scale 2D simulations that are currently prohibitive. The work builds on existing TK-WFP and sum-of-exponentials techniques with a new application and large-scale numerical validation.

major comments (1)
  1. [Abstract and numerical examples section] The accuracy and quasi-linear scaling claims rest on the far-history sum-of-exponentials approximation not accumulating significant error over many time steps. The abstract and numerical examples report 6-digit accuracy for the tested regimes, but no a priori bound is provided on the approximation error per step or demonstration that the total far-history error remains bounded (rather than growing with Nt) when the same exponential fit is used across increasing numbers of time steps, given the slowly decaying tail from the weak Huygens principle in 2D. This is load-bearing for the central performance claim.
minor comments (2)
  1. [Abstract] The local-part quadrature is described only as 'high-order' without specifying the rule, order, or error analysis for the non-smooth time-local contribution.
  2. [Abstract] Notation for the splitting parameters (e.g., the O(1) domain traversal time separating near and far history) could be defined more explicitly with a figure or equation.

Simulated Author's Rebuttal

1 responses · 1 unresolved

We thank the referee for the positive summary and significance assessment, as well as for identifying a key point regarding error control in the far-history approximation. We respond to the major comment below.

read point-by-point responses

  1. Referee: [Abstract and numerical examples section] The accuracy and quasi-linear scaling claims rest on the far-history sum-of-exponentials approximation not accumulating significant error over many time steps. The abstract and numerical examples report 6-digit accuracy for the tested regimes, but no a priori bound is provided on the approximation error per step or demonstration that the total far-history error remains bounded (rather than growing with Nt) when the same exponential fit is used across increasing numbers of time steps, given the slowly decaying tail from the weak Huygens principle in 2D. This is load-bearing for the central performance claim.

    Authors: We appreciate the referee highlighting this important consideration for the reliability of the performance claims. The manuscript's numerical examples, spanning up to 10^6 sources/targets and domains of 300x300 wavelengths, report stable 6-digit accuracy across the tested time-step counts, with no observed growth in error as Nt increases. This provides empirical support that the large-time sum-of-exponentials approximation to the wave kernel controls the far-history contribution without significant accumulation in the regimes considered. We agree, however, that the manuscript does not include an explicit study isolating the far-history error versus Nt or a rigorous a priori bound on the per-step approximation error. Deriving such a bound would require a detailed analysis of the approximation error for the 2D wave kernel under the weak Huygens principle, which lies outside the scope of the present algorithmic and validation-focused work. In the revised manuscript we will add a dedicated numerical study in the examples section that varies Nt while holding the exponential fit fixed, to explicitly demonstrate boundedness of the far-history error. We will also insert a short remark in the discussion noting that error control for this component is established numerically rather than via a priori estimates. revision: partial

standing simulated objections not resolved
  • Derivation of a rigorous a priori bound on the per-step approximation error and its accumulation for the far-history sum-of-exponentials approximation

Circularity Check

0 steps flagged

No significant circularity in derivation chain

full rationale

The paper introduces a new splitting of the wave solution into local (direct quadrature), near-history (TK-WFP), and far-history (new sum-of-exponentials approximation) components, with quasi-linear scaling demonstrated via numerical experiments on up to 10^6 sources. The TK-WFP reference is to a recent external method and is not used to justify any uniqueness theorem or to force the central performance claim by construction. The sum-of-exponentials fit is presented as novel here and validated numerically rather than renamed or fitted to the target result. No step reduces the claimed acceleration or accuracy to a self-referential definition, fitted input, or author-overlapping citation chain.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Based solely on the abstract, no free parameters, axioms, or invented entities are explicitly detailed or required for the central claim. The method relies on standard numerical analysis techniques and approximations.

pith-pipeline@v0.9.0 · 5596 in / 1260 out tokens · 78013 ms · 2026-05-10T16:56:31.456071+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

8 extracted references · 8 canonical work pages

  1. [1]

    fast-hybrid

    [1]N. G. Al Hassanieh, A. H. Barnett, and L. Greengard, Truncated kernel windowed Fourier projection: a fast algorithm for the 3D free-space wave equation, 2025, https://arxiv.org/ abs/2511.20824. [2]N. G. Al Hassanieh, A. H. Barnett, and L. Greengard, A fast algorithm for the wave equation using time-windowed Fourier projection, J. Comput. Phys., 548 (20...

    work page doi:10.1016/j.jcp.2025.114570 2025

  2. [2]

    Betcke, N

    [11]T. Betcke, N. Salles, and W. Smigaj, Overresolving in the Laplace domain for convolution quadrature methods, SIAM J. Sci. Comput., 39 (2017), https://doi.org/10. 1137/16M106474X. [12]O. P. Bruno and M. A. Santana, Efficient time-domain scattering synthesis via frequency-domain singularity subtraction, 2025, https://arxiv.org/abs/2505.06189. [13]W. C. ...

    work page arXiv 2017

  3. [3]

    Engquist and A

    [14]B. Engquist and A. Majda, Absorbing boundary conditions for numerical simulation of waves, Proceedings of the National Academy of Sciences, 74 (1977), pp. 1765–1766, https://doi. org/10.1073/pnas.74.5.1765. [15]I. Gradshteyn and I. Ryzhik, Table of Integrals, Series and Products, New York: Academic, 8th ed.,

    work page doi:10.1073/pnas.74.5.1765 1977

  4. [4]

    Greengard and P

    [16]L. Greengard and P. Lin, Spectral approximation of the free-space heat kernel, Applied and Computational Harmonic Analysis, 9 (2000), pp. 83–97, https://doi.org/10.1006/acha. 2000.0310. [17]L. Greengard and J. Strain, A fast algorithm for the evaluation of heat potentials, Comm. Pure Appl. Math., 43 (1990), pp. 949–963. [18]T. Hagstrom, D. Givoli, D. ...

    work page doi:10.1006/acha 2000

  5. [5]

    [24]J. Kaye, A. Barnett, and L. Greengard, A high-order integral equation-based solver for the time-dependent Schr¨ odingerequation, Communications on Pure and Applied Mathematics, 75 (2022), pp. 1657–1712, https://doi.org/10.1002/cpa.21959. [25]Y. Liu, A. C. Y ¨ucel, H. Ba ˘gcı, A. C. Gilbert, and E. Michielssen, A wavelet-enhanced PWTD-accelerated time-...

    work page doi:10.1002/cpa.21959 2022

  6. [6]

    Soffer and C

    [34]A. Soffer and C. Stucchio, Open boundaries for the nonlinear Schr¨ odingerequation, J. Comput. Phys., 225 (2007), pp. 1218–1232, https://doi.org/10.1016/j.jcp.2007.01.020. [35]A. Soffer, C. Stucchio, and M.-B. Tran, Time Dependent Phase Space Filters: A Stable Absorbing Boundary Condition, Springer Briefs on PDEs and Data Science, Springer Na- ture,

    work page doi:10.1016/j.jcp.2007.01.020 2007

  7. [7]

    Stucchio and A

    [37]C. Stucchio and A. Soffer, Multiscale resolution of shortwave-longwave interaction, Com- munications on Pure and Applied Mathematics, 62 (2009), pp. 82–124. [38]T. Takahashi, An interpolation-based fast-multipole accelerated boundary integral equation method for the three-dimensional wave equation, Journal of Computational Physics, 258 (2014), pp. 809...

    work page doi:10.1016/j.jcp.2013.11.008 2009

  8. [8]

    Wilber, W

    [41]H. Wilber, W. Vaes, A. Gopal, and G. Martinsson, A time-frequency method for acoustic scattering with trapping, 2025, https://arxiv.org/abs/2506.15165. [42]A. Yilmaz, J.-M. Jin, and E. Michielssen, Time domain adaptive integral method for surface integral equations, IEEE Transactions on Antennas and Propagation, 52 (2004), pp. 2692– 2708, https://doi....

    work page doi:10.1109/tap.2004.834399 2025