WaveHoltz: Time-Domain Helmholtz Solver

• WaveHoltz is a time-domain iterative method that solves the Helmholtz equation by filtering periodic wave solutions, yielding a positive definite system ideal for Krylov methods.
• It leverages one-period wave propagation combined with harmonic filtering to isolate the desired frequency component, enhancing convergence and numerical stability.
• The method extends to multiple applications, including multi-frequency solves, time-harmonic Maxwell's and elastic waves, and GPU-optimized domain decomposition for scalable computations.

WaveHoltz is a time-domain iterative method for solving the Helmholtz equation by filtering solutions of a related wave equation over one period. In its original formulation, the method was introduced as an iterative solution of the Helmholtz equation via the wave equation, with the key structural property that the iteration corresponds to a coercive operator, or to a positive definite matrix in the discretized case (Appelo et al., 2019). Subsequent work extended the framework to constant-coefficient problems on all of $\mathbb{R}^d$, multi-frequency Helmholtz solves, Maxwell and elastic waves, eigenvalue computation, low-rank compression, heterogeneous multiscale methods, and GPU-oriented domain decomposition, while also developing convergence theory, discrete error analysis, and acceleration strategies (Runborg et al., 17 Oct 2025).

1. Origins and conceptual basis

The motivating observation behind WaveHoltz is that the solution of a Helmholtz problem can be identified with a time-periodic solution of a wave equation driven at angular frequency $\omega$. Rather than solving the indefinite frequency-domain problem directly, WaveHoltz repeatedly propagates a wave equation for one period and applies a filter that isolates the $\omega$-frequency component. In the original formulation, the Helmholtz equation

$\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$

is related to the wave equation

$$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$

and the time-periodic solution $w(t,x)=u(x)e^{i\omega t}$ is the target object extracted by filtering (Appelo et al., 2019).

This construction addresses a standard obstacle in Helmholtz computation: direct discretizations are highly indefinite, whereas the WaveHoltz reformulation produces an operator $I-\mathcal S$ that is self-adjoint and positive definite away from resonance. That feature makes the method amenable to CG, GMRES, and related Krylov accelerations without ever assembling the indefinite Helmholtz operator explicitly (Appelo et al., 2019).

A common misunderstanding is to identify WaveHoltz with simple long-time marching to steady state. The method is more specific: it is a fixed-point iteration built from one-period wave propagation and harmonic filtering. This distinction is central to the convergence theory and to later extensions in which the filtered time-domain operator, rather than asymptotic time evolution alone, is the primary numerical object (Garcia et al., 2022).

2. Core mathematical formulation

A particularly explicit setting is the constant-coefficient Helmholtz equation on all of $\mathbb{R}^d$,

$\Delta u+\omega^2u=f(x),\qquad x\in\mathbb{R}^d,$

with the physical solution selected by the Sommerfeld radiation condition

$$\lim_{|x|\to\infty}|x|^{(d-1)/2}\left(\partial_r-i\omega\right)u(x)=0.$$

Because such solutions generally are not in standard Sobolev spaces, the analysis is formulated in weighted Sobolev spaces

$\omega$0

with corresponding spaces $\omega$1 and $\omega$2 (Runborg et al., 17 Oct 2025).

In the real-valued constant-coefficient formulation, the WaveHoltz operator is

$\omega$3

where $\omega$4 solves

$\omega$5

and the filter is

$\omega$6

The iteration is

$\omega$7

In this setting, the fixed point of $\omega$8 is the real part of the outgoing solution to the Helmholtz equation (Runborg et al., 17 Oct 2025).

For the more general energy-conserving formulation, the filter acts on both displacement and velocity,

$\omega$9

and the fixed-point relation is

$\omega$0

For real data with energy-conserving boundary conditions, the iteration simplifies to a displacement-only update (Appelo et al., 2019).

3. Spectral structure and convergence theory

The spectral description of WaveHoltz is organized around a filter transfer function. In the original eigenfunction-based analysis,

$\omega$1

with $\omega$2 and $\omega$3 for $\omega$4. If $\omega$5 is the operator with eigenvalues $\omega$6, then the WaveHoltz operator is $\omega$7, and the error satisfies

$\omega$8

Away from resonance, this yields geometric contraction in the underlying modal basis (Appelo et al., 2019).

The semi-discrete theory replaces eigenfunction expansions by assumptions on the spectrum of a stable spatial discretization. For stable semi-discretizations of the wave equation, the WaveHoltz iteration is guaranteed to converge to an approximate solution of the corresponding frequency domain problem, if it exists. For certain classes of frequency domain problems, the unaccelerated iteration converges in $\omega$9 iterations, with the constant factor depending logarithmically on the desired tolerance. It is conjectured that Helmholtz problems in open domains with no trapping waves belong to this class, and finite difference and discontinuous Galerkin experiments in one and two dimensions support that scaling (Rotem et al., 2024).

Separate extensions cover non-energy-conserving settings. Convergence has been proved for impedance boundary conditions in a single spatial dimension, and for interior Dirichlet or Neumann problems with damping in any spatial dimension. For a sufficient level of damping, the WaveHoltz iteration converges in a number of iteration independent of the frequency. The same line of work also proves that the fixed-point of the discrete WaveHoltz iteration converges to the discrete Helmholtz solution with the order of the time-stepper chosen (Garcia et al., 2022).

The most explicit whole-space convergence result currently available is the constant-coefficient analysis on $\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$0. There, the error $\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$1 satisfies

$\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$2

with Fourier representation

$\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$3

If $\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$4, $\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$5, and $\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$6, then

$\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$7

and if $\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$8,

$\nabla\cdot(c^2(x)\nabla u)+\omega^2u=f(x)$9

These estimates guarantee convergence of the real parts of the iterates to the real part of the outgoing solution in weighted Sobolev norms (Runborg et al., 17 Oct 2025).

4. Discretization, accuracy, and optimal-complexity implementations

Discrete analysis is a major part of the WaveHoltz literature because the practical fixed point depends on both time stepping and quadrature. A central result is that, for a family of higher order time-stepping schemes, the discrete fixed-point converges to the discrete Helmholtz solution with the order of the time-stepper chosen. Moreover, time discretization error can be completely removed through careful analysis of the discrete iteration together with updated quadrature formulas, and the paper reports this for centered modified equation schemes of order $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$0, $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$1, and $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$2 (Garcia et al., 2022).

This principle is developed algorithmically in the Multi-Frequency WaveHoltz (MFWH) method, which solves multiple Helmholtz problems simultaneously by solving a single wave equation combined with multiple time filters. MFWH defines a fixed-point iteration that can be accelerated with Krylov methods such as GMRES. The associated wave equation can be solved with either explicit time-stepping or implicit time-stepping using as few as five time-steps per period, and when combined with an $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$3 solver for the implicit equations, such a multigrid, the scheme has an $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$4 solution cost when the frequencies are fixed and the number of grid points $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$5 increases. The method uses high-order accurate approximations in space together with second-order accurate approximations in time, and includes an error-removal mechanism so that the MFWH solutions converge to the corresponding solutions to the discretized Helmholtz problems (Appelö et al., 31 Jul 2025).

For complex geometry, WaveHoltz has been combined with overset grids and implicit or implicit-explicit time stepping. The overset-grid Helmholtz solver uses Cartesian grids throughout most of the domain together with curvilinear grids near boundaries, accelerates the basic fixed-point iteration with GMRES and eigenmode deflation, and solves the wave equation with implicit time-stepping using as few as five time-steps per period, independent of the mesh size. With multigrid for the implicit equations, the resulting scheme scales linearly with the total number of grid points $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$6 at fixed frequency in both CPU-time and memory usage (Appelo et al., 3 Apr 2025). Closely related implicit and partitioned implicit-explicit modified-equation schemes were developed specifically for wave equations on overset grids, with the fully implicit schemes identified as useful for WaveHoltz applications where very large time-steps are desired (Carson et al., 2024).

5. Variants and generalizations

WaveHoltz has generated a family of related methods that preserve the time-filtering viewpoint while changing the target PDE, the spectral task, or the data representation.

Variant	Target problem	Distinctive feature
EM-WaveHoltz	time-harmonic Maxwell's equations	positive definite linear system; CG or GMRES
El-WaveHoltz	time-harmonic elastic wave equations	explicit and implicit schemes with time-discretization error removal
MFWH	multiple frequencies and right-hand sides	single wave solve with multiple time filters
EigenWave	eigenvalues and eigenvectors of elliptic boundary value problems	target frequency selection and matrix-free Arnoldi
LR-WaveHoltz	low-rank Helmholtz solution	SVD in two dimensions; tensor trains in three dimensions
WaveHoltz HMM	Helmholtz problems with rapidly varying coefficients	HMM wave solver combined with WaveHoltz filtering

EM-WaveHoltz reformulates time-harmonic Maxwell's equations through time-domain simulations and produces a positive definite system of equations amenable to conjugate gradient or GMRES. Theoretical results guarantee convergence away from resonances, and the framework can be wrapped around standard time-domain solvers such as FDTD or DGTD. It also extends to multiple frequencies in a single time-domain simulation (Peng et al., 2021).

El-WaveHoltz applies the same filtered fixed-point structure to time-harmonic elastic waves with energy conserving boundary conditions. As in the acoustic case, the fixed-point iteration is recast as a positive definite linear system solved by a Krylov method. The method includes one explicit and one novel implicit time-stepping scheme that completely remove time discretization error from the WaveHoltz solution after a simple modification of the initial data and time-stepping scheme. Numerical experiments indicate iteration scaling similar to that of the original WaveHoltz method, with the convergence rate dictated by the shortest, shear wave speed of the problem (Appelö et al., 2022).

EigenWave shifts the target from forced frequency response to spectral extraction. It computes eigenvalues and eigenvectors of elliptic boundary value problems by time-filtering the wave equation, allows the choice of an arbitrary target frequency, and can be embedded within a matrix-free Arnoldi algorithm to compute multiple eigenpairs near the target frequency. With implicit time-stepping and about $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$7 time-steps per period, and with multigrid for the implicit equations, the cost scales linearly with the number of grid points $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$8 as the mesh is refined (Appelo et al., 24 Jul 2025).

LR-WaveHoltz adapts the Helmholtz solver itself to low-rank structure. In two dimensions it uses the singular value decomposition, in three dimensions tensor trains, and to control rank growth it uses step-truncation during time stepping together with a low-rank Anderson acceleration for the WaveHoltz fixed-point iteration. Numerical experiments were reported for free- and half-space problems in two and three dimensions with constant and piecewise constant wave speeds (Granath et al., 10 Oct 2025).

The WaveHoltz Heterogeneous Multiscale Method combines a finite difference Heterogeneous Multiscale Method for the wave equation with the WaveHoltz iteration for the time-periodic Helmholtz solution. A notable property of this construction is that the time-domain solver does not artificially impose boundary conditions on the micro-scale problems, so no boundary errors from the micro-scale problems are present in the homogenized frequency domain solution (Rotem et al., 7 Jul 2026).

6. Acceleration, deflation, and high-performance computing

Although WaveHoltz is often used together with Krylov methods, WaveHoltz itself is not a Krylov method. This distinction becomes especially important in high-performance implementations. In GPU-oriented domain decomposition for Helmholtz problems, WaveHoltz is used as the subdomain solver because it is a fixed-point iteration that is uniquely well-suited to the GPU execution model due to its minimal memory footprint and no reduction operations. In the reported CUDA implementation on an NVIDIA A100, WaveHoltz achieves $$w_{tt}=\nabla\cdot(c^2(x)\nabla w)-f(x)e^{i\omega t},$$9-$w(t,x)=u(x)e^{i\omega t}$0 speedup over MINRES, with the advantage growing with subdomain size, and single-precision subdomain solves yield an additional $w(t,x)=u(x)e^{i\omega t}$1-$w(t,x)=u(x)e^{i\omega t}$2 speedup (Rotem, 19 Jun 2026).

Another acceleration direction is spectral deflation. For energy-conserving Dirichlet or Neumann boundary conditions, the WaveHoltz fixed-point iteration converges slowly at high frequency, requiring approximately $w(t,x)=u(x)e^{i\omega t}$3 iterations in $w(t,x)=u(x)e^{i\omega t}$4 dimensions. A numerical study of eigenvector deflation shows that deflating the eigenvectors whose eigenvalues lie nearest the driving frequency substantially reduces iteration counts. The study considers direct eigenvector deflation and augmented-Krylov eigenvector deflation using deflated conjugate gradient, augmented GMRES, and augmented BiCGSTAB. In two dimensions, when the number of deflation vectors grows quadratically with $w(t,x)=u(x)e^{i\omega t}$5, the asymptotic convergence rate remains essentially constant; the same study reports that the deflated solver breaks even against the undeflated solver after as few as two right-hand sides when the cost of precomputing the eigenvectors is included (Appelo et al., 30 Jun 2026).

This combination of operator reformulation, explicit spectral analysis, and architecture-aware implementation suggests a distinctive computational profile for WaveHoltz: the method is not defined by a single solver kernel, but by a filtered wave-propagation operator that can be accelerated by GMRES, CG, Anderson acceleration, deflation, multigrid-backed implicit stepping, low-rank compression, and block-level GPU concurrency, depending on the target problem and discretization regime (Appelo et al., 3 Apr 2025).