pith. sign in

arxiv: 2602.16591 · v2 · submitted 2026-02-18 · 🧮 math.NA · cs.NA

Fast Ewald Summation using Prolate Spheroidal Wave Functions

Erik Bostr\"om , Anna-Karin Tornberg , Ludvig af Klinteberg This is my paper

Pith reviewed 2026-05-15 21:04 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords Ewald summationprolate spheroidal wave functionsfast Fourier transformCoulomb interactionsmolecular dynamicserror estimateswindow functionsmollifiers
0
0 comments X

The pith

Prolate spheroidal wave functions enable Ewald summation with fewer Fourier modes and smaller windows than Gaussian or B-spline methods.

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

The paper shows that the first prolate spheroidal wave function, chosen for its optimal concentration in both real and Fourier space, can define both the mollifier that splits the Coulomb kernel into short- and long-range parts and the window used for spreading charges and interpolating forces onto an FFT grid. This choice produces rigorous bounds on truncation and aliasing errors together with closed-form approximations that let parameters be set so the realized error matches a user-specified tolerance. A reader cares because molecular dynamics and other particle simulations must recompute long-range forces many times; lowering the number of Fourier modes and the support size of the window directly reduces the dominant FFT cost while preserving accuracy. Experiments confirm the theoretical advantage holds in practice.

Core claim

The first prolate spheroidal wavefunction is used to define both the mollifier in the kernel split and the window function in the spreading and interpolation steps of fast Ewald summation. This yields rigorous error estimates and closed-form approximations of the Fourier truncation and aliasing errors, permitting parameter choices where the achieved error closely matches the prescribed tolerance. As a result, the PSWF-based approach achieves a given accuracy with significantly fewer Fourier modes and smaller window supports than Gaussian- and B-spline-based methods.

What carries the argument

The first prolate spheroidal wave function (PSWF), which possesses optimal joint concentration in physical and Fourier space and is applied both as the mollifier for the interaction kernel and as the window for FFT spreading and interpolation.

If this is right

  • Parameters for any target tolerance can be obtained directly from the closed-form error approximations without iteration.
  • The FFT step operates on a coarser grid, cutting its arithmetic cost.
  • Spreading and interpolation steps operate with smaller supports, reducing their work per particle.
  • The resulting scheme supplies a practical alternative for repeated long-range force evaluations in particle simulations.

Where Pith is reading between the lines

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

  • The same replacement of mollifier and window could be tested on other long-range kernels that admit an Ewald-type split.
  • In large-scale molecular dynamics the reduction in FFT size may permit either larger system sizes or longer trajectories at fixed wall-clock time.
  • Existing Ewald codes could adopt the approach by swapping only the two functions while retaining the rest of the infrastructure.

Load-bearing premise

The theoretical concentration optimality of the PSWF translates directly into lower total cost in the complete Ewald pipeline for arbitrary particle configurations without hidden stability or implementation penalties.

What would settle it

A numerical test on a collection of charged particles in which the PSWF-based method requires as many Fourier modes as a Gaussian-based method to reach the same accuracy.

Figures

Figures reproduced from arXiv: 2602.16591 by Anna-Karin Tornberg, Erik Bostr\"om, Ludvig af Klinteberg.

Figure 1
Figure 1. Figure 1: Comparison of the Gaussian and the first PSWF of order zero in real space (left) and Fourier space (right). For equal effective support, the PSWF mollifier requires roughly half as many Fourier modes to achieve the same accuracy. Until recently, PSWF analogues, such as the KB and ES functions, have primarily been considered as window functions in fast Ewald summation. However, employing them also as mollif… view at source ↗
Figure 2
Figure 2. Figure 2: The tail of |ψ c 0 (s)| for c = √ 3 and c = 10, shown together with the bound (32), the approximation obtained from (34) by neglecting Ac (s), and its corresponding upper bound. By Lemma 4, the bound |ψ c 0 (s)| ≤ ψ c 0 (1)/s holds for c ≥ √ 3 (left). For large s, the approximation 2ψ c 0 (1) λ0 | sin(cs)| cs accurately captures the oscillatory tail (right). Lemma 3. Let c > 0. Then the smallest eigenvalue… view at source ↗
Figure 3
Figure 3. Figure 3: Split function Φ cs rc (x) (left) and residual kernel R(x) (right) for rc = 0.5 and cs ∈ {1, 6, 11, . . . , 49}. By construction, the split function satisfies Φ cs rc (rc) = 1, and the residual kernel is compactly supported on [−rc, rc] [PITH_FULL_IMAGE:figures/full_fig_p015_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Resolution requirements for direct Ewald summation using PSWF and Gaussian splits. Relative ℓ2 errors are shown for cutoff radii rc ∈ {0.1, 0.3, 0.5} with L = 1 and n = 100. Dashed lines indicate reference slopes log(1/ε) (PSWF) and 2 log(1/ε) (Gaussian). Expressing the resolution as ωmaxrc causes the results to collapse and confirms the expected scaling with respect to rc. For practically relevant toleran… view at source ↗
Figure 5
Figure 5. Figure 5: Absolute RMS errors for fast Ewald summation (Algorithm 2) using the PSWF split and PSWF window (L = 1, rc = 0.1, n = 100). Curves with symbols show numer￾ical results, while black curves indicate theoretical error models. Dotted lines indicate the simple asymptotic exponential model e −{cs,cw} ; solid lines correspond to the split-error model ε = √ 1 V √ rc ∥ρ∥2 Asc −1/2 s e −cs (89) with As ≈ 5; and dash… view at source ↗
Figure 6
Figure 6. Figure 6: Dependence of the absolute RMS error on the system size n and the cutoff radius rc for fast Ewald summation with PSWF split and PSWF mollifier. (a) Window support vs. absolute RMS error for varying n. (b) Same data as in (a), rescaled by ∥ρ∥2. (c) Window support vs. absolute RMS error for varying rc. (d) Same data as in (c), rescaled by √ rc. Curves with symbols show numerical results, while black curves i… view at source ↗
Figure 7
Figure 7. Figure 7: Contour plots of the absolute RMS error (L = 1, rc = 0.1, n = 100). In (a), error in (m, P) space. In (b), error in (cs, cw) space. The two parameter spaces are related linearly by cs = (πrc/L) m and cw = (π/2) P. The green curve in each panel shows the theoretical optimal relation (126); the computed minimizers (m, P) (or equivalently (cs, cw)) for each error level align closely with this curve. For refer… view at source ↗
Figure 8
Figure 8. Figure 8: Measured absolute RMS error as a function of the requested tolerance ε for different values of n and rc [PITH_FULL_IMAGE:figures/full_fig_p034_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Total number of Fourier modes m3 versus relative ℓ2 error for PSWF/B-spline and Gaussian/B-spline (SPME) configurations at L = 1, rc = 0.1, n = 100, and B-spline order P ∈ {4, 6, 8}. The split parameters cs (PSWF) and σ (Gaussian) were chosen to yield the same direct-sum split error of 10−5 . For small and moderate m, the PSWF split yields faster error decay and reaches a given tolerance with fewer modes; … view at source ↗
Figure 10
Figure 10. Figure 10: Curve-fitted approximations of key PSWF quantities [PITH_FULL_IMAGE:figures/full_fig_p041_10.png] view at source ↗
read the original abstract

Fast Ewald summation efficiently evaluates Coulomb interactions and is widely used in molecular dynamics simulations. It is based on a split into a short-range and a long-range part, where evaluation of the latter is accelerated using the fast Fourier transform (FFT). The accuracy and computational cost depend critically on the mollifier in the kernel split and the window function used in the spreading and interpolation steps that enable the use of the FFT. The first prolate spheroidal wavefunction (PSWF) has optimal concentration in real and Fourier space simultaneously, and is used when defining both a mollifier and a window function. We provide a complete description of the method and derive rigorous error estimates. In addition, we obtain closed-form approximations of the Fourier truncation and aliasing errors, yielding explicit parameter choices for the achieved error to closely match the prescribed tolerance. Numerical experiments confirm the analysis: PSWF-based Ewald summation achieves a given accuracy with significantly fewer Fourier modes and smaller window supports than Gaussian- and B-spline-based approaches, providing a superior alternative to existing Ewald methods for particle simulations.

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

2 major / 2 minor

Summary. The manuscript introduces a fast Ewald summation method for Coulomb interactions in particle simulations that uses the first prolate spheroidal wave function (PSWF) to define both the mollifier in the short-range/long-range kernel split and the window function for spreading and interpolation steps prior to the FFT. It derives rigorous error estimates together with closed-form approximations to the Fourier truncation and aliasing errors that permit explicit parameter choices matching a prescribed tolerance, and reports numerical experiments showing that the PSWF approach attains a target accuracy with substantially fewer Fourier modes and smaller window supports than Gaussian- or B-spline-based Ewald schemes.

Significance. If the error bounds prove tight and the reported reductions in modes and support translate into lower overall computational cost without offsetting precomputation or stability penalties, the method would constitute a meaningful improvement over existing Ewald techniques for molecular-dynamics workloads. The exploitation of the known optimal concentration properties of the first PSWF is a natural and potentially powerful choice; the provision of closed-form error approximations and machine-checkable parameter selection would be a concrete practical strength.

major comments (2)
  1. [Abstract and §4] Abstract and §4 (error analysis): the rigorous error estimates and closed-form truncation/aliasing approximations must be shown to remain uniform for the nonuniform particle distributions typical of molecular dynamics; the current derivation appears to rely on periodic uniform-grid assumptions that may not carry over directly to the full spreading-interpolation pipeline.
  2. [Numerical experiments] Numerical experiments section: the reported reductions in Fourier modes and window support size are quantified, yet no breakdown of total flop count—including the cost of evaluating or tabulating high-order PSWFs and computing prolate eigenvalues for tolerances below 10^{-8}—is provided, leaving open whether the headline savings produce net wall-clock or memory gains.
minor comments (2)
  1. Notation for the PSWF mollifier and window should be introduced with a single consistent symbol set early in the manuscript to avoid later confusion between the two uses.
  2. Figure captions comparing error versus number of modes should explicitly state the tolerance and particle configuration used for each curve.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the thorough review and valuable suggestions. We address the major comments below and will revise the manuscript accordingly where appropriate.

read point-by-point responses

  1. Referee: [Abstract and §4] Abstract and §4 (error analysis): the rigorous error estimates and closed-form truncation/aliasing approximations must be shown to remain uniform for the nonuniform particle distributions typical of molecular dynamics; the current derivation appears to rely on periodic uniform-grid assumptions that may not carry over directly to the full spreading-interpolation pipeline.

    Authors: The error analysis in §4 is based on the properties of the PSWF mollifier and window, which provide uniform bounds on the truncation and aliasing errors independent of the particle positions. The spreading and interpolation steps use the same window function, and the error bounds derived from the optimal concentration of the PSWF hold for arbitrary charge distributions because they rely on the Fourier decay properties rather than grid uniformity. The periodic uniform-grid assumption is only for the FFT computation itself, but the error estimates carry over to the nonuniform case via the standard analysis of the particle-mesh method. We will add a clarifying paragraph in §4 to explicitly state this uniformity for molecular dynamics applications. revision: partial

  2. Referee: [Numerical experiments] Numerical experiments section: the reported reductions in Fourier modes and window support size are quantified, yet no breakdown of total flop count—including the cost of evaluating or tabulating high-order PSWFs and computing prolate eigenvalues for tolerances below 10^{-8}—is provided, leaving open whether the headline savings produce net wall-clock or memory gains.

    Authors: We agree that including a detailed computational cost analysis would be beneficial. In the revised manuscript, we will add a subsection in the numerical experiments providing flop counts for the main operations, precomputation costs for PSWFs and eigenvalues (which are computed once per tolerance and can be tabulated), and wall-clock timings on representative MD systems. This will demonstrate that the reductions in modes and support size lead to overall efficiency gains even after accounting for precomputation. revision: yes

Circularity Check

0 steps flagged

No significant circularity: claims rest on classical PSWF optimality and independent error derivations

full rationale

The paper selects the first PSWF for mollifier and window based on its established optimal concentration properties (a classical result from Slepian and others, external to this work). It then derives rigorous error estimates, closed-form truncation/aliasing approximations, and explicit parameter choices directly from Fourier analysis and the kernel split. These steps are self-contained and do not reduce by construction to fitted inputs, self-definitions, or load-bearing self-citations. Numerical experiments serve only as confirmation, not as definitional inputs. No ansatz smuggling, renaming of known results, or uniqueness theorems imported from the authors' prior work appear in the derivation chain.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The method relies on the known optimal concentration properties of the first prolate spheroidal wave function and on standard Fourier-analysis error bounds; no new free parameters, invented entities, or ad-hoc axioms are introduced in the abstract.

axioms (2)
  • domain assumption The first prolate spheroidal wave function has optimal simultaneous concentration in real and Fourier space
    Invoked to justify its use as both mollifier and window function
  • standard math Standard Fourier truncation and aliasing error analysis applies to the smoothed kernel
    Used to derive the rigorous bounds and closed-form approximations

pith-pipeline@v0.9.0 · 5493 in / 1453 out tokens · 21114 ms · 2026-05-15T21:04:28.859645+00:00 · methodology

discussion (0)

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

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

What do these tags mean?
matches
The paper's claim is directly supported by a theorem in the formal canon.
supports
The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends
The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses
The paper appears to rely on the theorem as machinery.
contradicts
The paper's claim conflicts with a theorem or certificate in the canon.
unclear
Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

36 extracted references · 36 canonical work pages · 1 internal anchor

  1. [1]

    af Klinteberg, L

    L. af Klinteberg, L. Greengard, S. Jiang, and A.-K. Tornberg. Fast summa- tion of Stokes potentials using a new kernel-splitting in the DMK framework. arXiv:2509.21471 [math.NA], 2025. 37

    work page arXiv 2025

  2. [2]

    af Klinteberg, D

    L. af Klinteberg, D. S. Shamshirgar, and A.-K. Tornberg. Fast Ewald summation for free-space Stokes potentials.Research in the Mathematical Sciences, 4(1):1, 2017

    work page 2017

  3. [3]

    Bagge and A.-K

    J. Bagge and A.-K. Tornberg. Fast Ewald summation for Stokes flow with arbitrary periodicity.Journal of Computational Physics, 493:112473, 2023

    work page 2023

  4. [4]

    A. H. Barnett. Aliasing error of the exp(β √ 1−z 2) kernel in the nonuniform fast fourier transform.Applied and Computational Harmonic Analysis, 51:1–16, 2021

    work page 2021

  5. [5]

    Exponential of Semicircle

    A. H. Barnett, J. Magland, and L. Af Klinteberg. A Parallel Nonuniform Fast Fourier Transform Library Based on an “Exponential of Semicircle” Kernel.SIAM Journal on Scientific Computing, 41(5):C479–C504, 2019

    work page 2019

  6. [6]

    G. Beylkin. On the Fast Fourier Transform of Functions with Singularities.Applied and Computational Harmonic Analysis, 2(4):363–381, 1995

    work page 1995

  7. [7]

    J. W. Cooley and J. W. Tukey. An Algorithm for the Machine Calculation of Complex Fourier Series.Mathematics of Computation, 19(90):297–301, 1965

    work page 1965

  8. [8]

    Darden, D

    T. Darden, D. York, and L. Pedersen. Particle mesh Ewald: AnN·log(N) method for Ewald sums in large systems.The Journal of Chemical Physics, 98(12):10089– 10092, 1993

    work page 1993

  9. [9]

    T. M. Dunster. Uniform Asymptotic Expansions for Prolate Spheroidal Functions with Large Parameters.SIAM Journal on Mathematical Analysis, 17(6):1495–1524, 1986

    work page 1986

  10. [10]

    Dutt and V

    A. Dutt and V. Rokhlin. Fast Fourier Transforms for Nonequispaced Data.SIAM Journal on Scientific Computing, 14(6):1368–1393, 1993

    work page 1993

  11. [11]

    Essmann, L

    U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Peder- sen. A smooth particle mesh Ewald method.The Journal of Chemical Physics, 103(19):8577–8593, 1995

    work page 1995

  12. [12]

    P. P. Ewald. Die Berechnung optischer und elektrostatischer Gitterpotentiale.An- nalen der Physik, 369(3):253–287, 1921

    work page 1921

  13. [13]

    W. H. J. Fuchs. On the eigenvalues of an integral equation arising in the theory of band-limited signals.Journal of Mathematical Analysis and Applications, 9(3):317– 330, 1964

    work page 1964

  14. [14]

    Greengard and J.-Y

    L. Greengard and J.-Y. Lee. Accelerating the Nonuniform Fast Fourier Transform. SIAM Review, 46(3):443–454, 2004

    work page 2004

  15. [15]

    R. W. Hockney and J. W. Eastwood.Computer Simulation Using Particles. CRC Press, Boca Raton, 2021

    work page 2021

  16. [16]

    Hofmann, F

    M. Hofmann, F. Nestler, and M. Pippig. NFFT based Ewald summation for electro- static systems with charges and dipoles.Applied Numerical Mathematics, 122:39– 65, 2017

    work page 2017

  17. [17]

    Jiang and L

    S. Jiang and L. Greengard. A dual-space multilevel kernel-splitting framework for discrete and continuous convolution.Communications on Pure and Applied Mathematics, 78(5):1086–1143, 2025. 38

    work page 2025

  18. [18]

    Keiner, S

    J. Keiner, S. Kunis, and D. Potts. Using NFFT 3—A software library for vari- ous nonequispaced fast Fourier transforms.ACM Transactions on Mathematical Software (TOMS), 36(4):1–30, 2009

    work page 2009

  19. [19]

    Accelerating Molecular Dynamics Simulations using Fast Ewald Summation with Prolates

    J. Liang, L. Lu, A. Barnett, L. Greengard, and S. Jiang. Accelerating Fast Ewald Summation with Prolates for Molecular Dynamics Simulations.arXiv:2505.09727 [math.NA], 2025

    work page internal anchor Pith review Pith/arXiv arXiv 2025

  20. [20]

    Lindbo and A.-K

    D. Lindbo and A.-K. Tornberg. Spectral accuracy in fast Ewald-based methods for particle simulations.Journal of Computational Physics, 230(24):8744–8761, 2011

    work page 2011

  21. [21]

    F. Nestler. Automated parameter tuning based on RMS errors for nonequispaced FFTs.Advances in Computational Mathematics, 42(4):889–919, 2016

    work page 2016

  22. [22]

    Nestler, M

    F. Nestler, M. Pippig, and D. Potts. Fast Ewald summation based on NFFT with mixed periodicity.Journal of Computational Physics, 285:280–315, 2015

    work page 2015

  23. [23]

    Osipov and V

    A. Osipov and V. Rokhlin. On the evaluation of prolate spheroidal wave functions and associated quadrature rules.Applied and Computational Harmonic Analysis, 36(1):108–142, 2014

    work page 2014

  24. [24]

    Osipov, V

    A. Osipov, V. Rokhlin, and H. Xiao.Prolate Spheroidal Wave Functions of Order Zero: Mathematical Tools for Bandlimited Approximation, volume 187 ofApplied Mathematical Sciences. Springer US, Boston, MA, 2013

    work page 2013

  25. [25]

    Plonka, D

    G. Plonka, D. Potts, G. Steidl, and M. Tasche.Numerical Fourier Analysis. Applied and Numerical Harmonic Analysis. Springer International Publishing, Cham, 2018

    work page 2018

  26. [26]

    Potts, G

    D. Potts, G. Steidl, and M. Tasche. Fast Fourier Transforms for Nonequispaced Data: A Tutorial. In J. J. Benedetto and P. J. S. G. Ferreira, editors,Modern Sam- pling Theory: Mathematics and Applications, pages 247–270. Birkh¨ auser, Boston, MA, 2001

    work page 2001

  27. [27]

    I. J. Schoenberg.Cardinal Spline Interpolation. Society for Industrial and Applied Mathematics, 1973

    work page 1973

  28. [28]

    D. S. Shamshirgar, J. Bagge, and A.-K. Tornberg. Fast Ewald summation for electrostatic potentials with arbitrary periodicity.The Journal of Chemical Physics, 154(16):164109, 2021

    work page 2021

  29. [29]

    D. Slepian. Some comments on Fourier analysis, uncertainty and modeling.SIAM Review, 25(3):379–393, 1983

    work page 1983

  30. [30]

    Slepian and H

    D. Slepian and H. O. Pollak. Prolate Spheroidal Wave Functions, Fourier Analysis and Uncertainty - I.Bell System Technical Journal, 40(1):43–63, 1961

    work page 1961

  31. [31]

    G. Steidl. A note on fast Fourier transforms for nonequispaced grids.Advances in Computational Mathematics, 9(3):337–352, 1998

    work page 1998

  32. [32]

    Szeg¨ o.Orthogonal polynomials, volume 23

    G. Szeg¨ o.Orthogonal polynomials, volume 23. American Mathematical Soc., 1939

    work page 1939

  33. [33]

    T. A. Driscoll, N. Hale, and L. N. Trefethen, editors.Chebfun Guide. Pafnuty Publications, Oxford, 2014. 39

    work page 2014

  34. [34]

    A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.Computer Physics Communi...

    work page 2022

  35. [35]

    Van Der Spoel, E

    D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen. GROMACS: Fast, flexible, and free.Journal of Computational Chem- istry, 26(16):1701–1718, 2005. A PSWF asymptotics and approximations In this section, we summarize the asymptotic behavior ofψ c 0, along with corresponding approximations. These results are used to deriv...

    work page 2005

  36. [36]

    Thus,Fis nonincreasing on (1,∞)

    Hence, (pq) ′(s)≥0 for alls >1, and thereforeF ′(s)≤0 on (1,∞). Thus,Fis nonincreasing on (1,∞). Therefore, u(s)2 ≤F(s)≤lim t→1+ F(t) =u(1) 2 = ψc 0(1) 2.(159) It follows that |ψc 0(s)| ≤ ψc 0(1) s .(160) D Auxiliary standard results D.1 Radial extension (Fourier-side definition) Letf:R→Rbe even and belong toL 1(R), with Fourier transform bf. Define the r...

    work page