Complex-Shifted Helmholtz Operator

• Complex-Shifted Helmholtz operators are defined by shifting spectral parameters into the complex plane to introduce damping and enhance solver convergence.
• They encompass variants such as the shifted Laplacian, contour-based shifted systems, complex-scaled grids, and complex-frequency formulations used in preconditioning and spectral filtering.
• This approach improves numerical stability in iterative solvers like GMRES and multigrid through effective spectral relocation, optimal shift selection, and innovative preconditioning architectures.

The complex-shifted Helmholtz operator denotes a family of Helmholtz-type operators in which the spectral parameter, the wavenumber term, the grid, or the frequency is moved into the complex plane. In the formulations collected here, this includes shifted systems of the form $(A-zI)y=f$, complex-shifted Laplacians such as $-\Delta-(1+\beta i)k^2$, complex-scaled grids with $h\to he^{i\theta}$, and complex-frequency operators $\mathcal L u=-\Delta u+\zeta^2u$ with $\operatorname{Re}\zeta>0$. Across these variants, the common purpose is to introduce damping, move the spectrum away from the origin or the real axis, and thereby make indefinite high-frequency Helmholtz problems more amenable to multigrid, Krylov, contour-integral, and compression-based algorithms (1811.12378, Cools et al., 2011, Börm et al., 2019).

1. Operator forms and mathematical variants

The basic Helmholtz model in the cited work appears in several equivalent or closely related forms. For wave propagation, one standard form is

$(-\Delta-k(x)^2)u(x)=f(x),$

or, with variable coefficients and damping,

$$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$

In discretized notation this becomes $Au=f$, where $A$ is large, highly indefinite, and non-Hermitian or non-normal in the high-frequency regime (Kazemi, 2017).

A first major variant is the complex shifted Laplacian. In the notation of S. Cools and W. Vanroose, the modified coefficient is

$\tilde{\sigma}=\sigma(1+\beta i)=-k^2(1+\beta i),$

leading to an operator of the form

$-\Delta-(1+\beta i)k^2$0

or, equivalently,

$-\Delta-(1+\beta i)k^2$1

Here $-\Delta-(1+\beta i)k^2$2 is the complex shift parameter. The discretized one-dimensional example is

$-\Delta-(1+\beta i)k^2$3

with analogous higher-dimensional structures (Cools et al., 2011).

A second variant is the shifted linear system used in contour-based formulations:

$-\Delta-(1+\beta i)k^2$4

with complex shift $-\Delta-(1+\beta i)k^2$5, typically chosen with $-\Delta-(1+\beta i)k^2$6. In this setting, the shift is not introduced by modifying $-\Delta-(1+\beta i)k^2$7 directly, but by solving a family of resolvent systems at quadrature points on a contour in the complex plane (1811.12378).

A third variant is complex scaling or complex stretched grids. In the spacetree-based multilevel framework, the domain is rotated by

$-\Delta-(1+\beta i)k^2$8

so that the mesh width is effectively scaled as $-\Delta-(1+\beta i)k^2$9. This yields a discretized operator assembled on complex-scaled mesh regions and is explicitly described as inspired by complex shifted Laplacian preconditioning (Reps et al., 2015). A related complex stretched grid construction discretizes the Laplacian on a grid with complex-valued mesh widths, producing a preconditioning operator whose spectrum lies in a triangle with vertices determined by the interior and exterior complex mesh widths (Reps et al., 2010).

A fourth variant appears in boundary integral and compression contexts through complex frequency:

$h\to he^{i\theta}$0

In this case the fundamental solution decays exponentially for large arguments,

$h\to he^{i\theta}$1

so the complex shift is encoded in the real and imaginary parts of $h\to he^{i\theta}$2 rather than in a separate Laplacian preconditioner (Börm et al., 2019).

2. Spectral effect and analytical rationale

The principal analytical role of the complex shift is spectral relocation. In the contour-integral formulation for the 3D high-frequency Helmholtz equation, the shifted operators $h\to he^{i\theta}$3 are chosen so that the spectrum is moved away from the real axis and the origin; the paper states that the imaginary part of $h\to he^{i\theta}$4 is typically negative and that this improves definiteness and ensures invertibility or stability of $h\to he^{i\theta}$5 (1811.12378).

This spectral effect is central to the efficacy of Krylov methods. The heterogeneous Helmholtz study using shifted-Laplace preconditioners states that the Helmholtz operator is non-unitary and non-diagonalizable and that a complex shift alters the spectrum by shifting and clustering eigenvalues away from the origin and the negative real axis toward the positive real axis, ideally around the unit circle. In the reported eigenvalue plots, the preconditioned spectrum moves away from the negative real axis and the origin, and this is directly associated with faster GMRES convergence (Kazemi, 2017).

For multigrid, the same idea is expressed in smoothing language. Ira Livshits’ shifted Laplacian analysis characterizes the original Helmholtz operator as indefinite, with troublesome near-kernel components and Fourier modes for which Gauss–Seidel relaxation may diverge. Adding $h\to he^{i\theta}$6 moves eigenvalues away from the real axis, makes the shifted operator less indefinite, and regularizes highly oscillatory near-nullspace components so that relaxation and coarse-grid correction become more effective (Livshits, 2013).

In local Fourier analysis, the shift is quantified by an amplification factor $h\to he^{i\theta}$7. The minimal shift is defined by

$h\to he^{i\theta}$8

so that all Fourier modes contract. The dependence of $h\to he^{i\theta}$9 on wavenumber, mesh size, smoothing steps, and multigrid cycle structure is explicitly described as irregular, with the instability mechanism switching between coarse-grid correction and the smoother depending on the regime (Cools et al., 2011).

The complex-frequency formulation yields a related but distinct spectral mechanism. When $\mathcal L u=-\Delta u+\zeta^2u$0, the exponential factor $\mathcal L u=-\Delta u+\zeta^2u$1 produces decay, and the compression analysis is explicit in both $\mathcal L u=-\Delta u+\zeta^2u$2 and $\mathcal L u=-\Delta u+\zeta^2u$3. Higher values of $\mathcal L u=-\Delta u+\zeta^2u$4 reduce complexity and may even allow the discrete matrix to be replaced by its nearfield part (Börm et al., 2019).

3. Preconditioning architectures and iterative solvers

The complex-shifted Helmholtz operator is most often used as a preconditioner rather than as a replacement for the physical problem. In right-preconditioned form, one solves

$\mathcal L u=-\Delta u+\zeta^2u$5

with $\mathcal L u=-\Delta u+\zeta^2u$6 and complex $\mathcal L u=-\Delta u+\zeta^2u$7 in one heterogeneous Helmholtz implementation. The numerical examples summarized there use $\mathcal L u=-\Delta u+\zeta^2u$8, $\mathcal L u=-\Delta u+\zeta^2u$9 and report that preconditioned GMRES rapidly converges in 1D and yields results comparable to direct LU on a 2D Marmousi example at $\operatorname{Re}\zeta>0$0 (Kazemi, 2017).

Multigrid inversion of the shifted operator is a standard architecture, but the details vary. One matrix-free parallel two-level deflation method combines a Complex Shifted Laplacian preconditioner with matrix-free Galerkin coarsening or high-order re-discretization on the coarse grid and approximates the CSLP by one parallel geometric multigrid V-cycle. The matrix-vector multiplications and interpolation or restriction operators are implemented on finite-difference grids without constructing any coefficient matrix, and the reported experiments show wavenumber independence for medium wavenumbers together with satisfactory weak and strong parallel scalability (Chen et al., 2023).

Several papers replace classical point smoothers by more robust inner iterations on the shifted system. In the complex stretched-grid framework, GMRES is used as a smoother on each multigrid level because stationary smoothers can become unstable for indefinite spectra. Very few GMRES iterations per level are reported as sufficient to build a robust multigrid solver for the complex shifted or stretched problem (Reps et al., 2010). A different line of development uses Hermitian Skew-Hermitian Splitting. The shifted HSS approach applies HSS iteration to a shifted operator

$\operatorname{Re}\zeta>0$1

and proves $\operatorname{Re}\zeta>0$2- and mesh-robustness when $\operatorname{Re}\zeta>0$3 HSS iterations are performed. The shifted systems arising inside HSS are described as suitable for approximation by multigrid using standard smoothers and transfer operators, leading to a fully scalable algorithm (Cotter et al., 23 Jun 2025).

Contour integration leads to another solver architecture. The spectral projector identity

$\operatorname{Re}\zeta>0$4

is discretized by quadrature, so that the outer stage approximates $\operatorname{Re}\zeta>0$5 through a sum of shifted solves and the inner stage computes a correction in the complementary spectral range. In the cited 3D study, the shifted systems are solved by polynomial fixed-point iteration, and the reduced spectrum of the inner problem yields faster GMRES iterations (1811.12378).

4. Shift selection, optimality, and scalability trade-offs

Choosing the shift is a central problem because the same parameter controls multigrid stability, Krylov efficiency, and fidelity to the original Helmholtz problem. The local Fourier analysis of Cools and Vanroose makes the trade-off explicit: a large $\operatorname{Re}\zeta>0$6 improves multigrid stability, but preconditioning becomes less effective for the Krylov method. Their numerical experiments support the claim that $\operatorname{Re}\zeta>0$7 is the minimal requirement for a multigrid V-cycle to converge and is near-optimal in terms of Krylov iteration count (Cools et al., 2011).

This trade-off reappears in data-driven shift selection for a semi matrix-free two-grid shifted Laplacian method. There, the shifted problem is

$\operatorname{Re}\zeta>0$8

and the paper seeks a near-optimal exponent $\operatorname{Re}\zeta>0$9 in $(-\Delta-k(x)^2)u(x)=f(x),$0, with $(-\Delta-k(x)^2)u(x)=f(x),$1. The reported workflow generates samples over $(-\Delta-k(x)^2)u(x)=f(x),$2, optimizes the average residual convergence rate with respect to $(-\Delta-k(x)^2)u(x)=f(x),$3, and then fits a nonlinear regression map $(-\Delta-k(x)^2)u(x)=f(x),$4 so that the near-optimal shift can be obtained by simple evaluation. The benchmarks in two and three dimensions report that the map-based shift yields the smallest FGMRES iteration counts in nearly all scenarios tested and provides substantial speedups over fixed shift choices (Drzisga et al., 2021).

The literature also stresses that increasing the shift is not a monotone route to scalability. A recent multigrid study states that the standard complex shift prevents scalability, that convergence of the classical complex-shifted approach requires $(-\Delta-k(x)^2)u(x)=f(x),$5 as the wavenumber increases, and that shifting by more than $(-\Delta-k(x)^2)u(x)=f(x),$6 destroys scalability. Its proposed real-shifted coarse-grid correction therefore avoids a complex shift on the main operator and uses a real shift only on the coarsest grid; for 10 points per wavelength, the method combined with a modest complex shift is reported to outperform the standard complex shifted Laplacian method by an order of magnitude (Yovel et al., 21 Apr 2026).

A related 2025 development replaces point-wise smoothing by additive Vanka smoothing in shifted Laplacian multigrid. The stated objective is to require a much lower shift than point-wise smoothers and thereby enhance scalability. Local Fourier analysis and numerical experiments are reported to show that deep V-cycles with a small and bounded shift are feasible, even on challenging geophysical media in 2D and 3D (Yovel et al., 20 Nov 2025).

A common misconception is that the “best” shift is simply the largest one for which multigrid converges. The cited results consistently reject that interpretation: over-damping degrades the preconditioner, may increase Krylov iterations, and can destroy scalability (Cools et al., 2011, Drzisga et al., 2021, Yovel et al., 21 Apr 2026).

5. Contour, multi-shift, and boundary-integral formulations

The complex-shifted Helmholtz operator is also a vehicle for reformulating Helmholtz computations into families of shifted solves. In contour-integration methods, the quadrature approximation

$(-\Delta-k(x)^2)u(x)=f(x),$7

acts as a rational filter that restricts the solution to a favorable spectral subspace. Because each quadrature point produces an independent shifted system, the method is described as embarrassingly parallel. The same paper reports small storage requirements, applicability to both dense and sparse linear systems, and $(-\Delta-k(x)^2)u(x)=f(x),$8 matrix-vector products for solving a high-frequency 3D problem of size $(-\Delta-k(x)^2)u(x)=f(x),$9 with high accuracy (1811.12378).

The multi-shift structure is not confined to contour integration. A general non-Hermitian family

$$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$0

arises in PDEs, FDEs, control theory, and rational approximations to matrix functions. The shifted CMRH literature exploits Krylov subspace shift-invariance,

$$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$1

so that all shifted systems can share a single Krylov basis. Restarted and flexible variants enforce residual collinearity across cycles and allow variable preconditioning while preserving the shifted structure (Gu et al., 2016).

In boundary element formulations, the complex shift may appear as a complex wavenumber rather than an additive Laplacian perturbation. For the EFIE, one cited preconditioner uses the shifted Helmholtz operator

$$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$2

with

$$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$3

The negative complex part is stated to ensure damping, avoid spurious resonances, and ensure invertibility. The preconditioned discrete operator

$$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$4

is motivated by pseudo-differential order, and the paper reports stabilization of the number of GMRES iterations in low-frequency, dense-refinement, and high-frequency regimes together with quasi-linear complexity when proper acceleration strategies are used (Ciciriello et al., 5 Jun 2026).

Complex frequency is equally important for matrix compression. Directional $$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$5-matrix methods for $$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$6 with $$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$7 use a new admissibility condition explicit in $$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$8 and $$-\nabla\cdot \rho \nabla \mathbf u-k^2(1-\hat{\jmath}\alpha)\mathbf u=\mathbf f.$$9, a variable-order expansion

$Au=f$0

and error estimates whose decay is explicit in $Au=f$1. The reported complexity is $Au=f$2 in the sectorial regime, while high damping may eliminate all farfield expansions so that only the nearfield remains (Börm et al., 2019).

6. Limitations, alternatives, and ongoing questions

Although complex shifting is one of the standard cures for indefiniteness, the cited work presents several limitations. One recurring limitation is that the optimal parameter choice remains problem-dependent. The heterogeneous Helmholtz data-modelling study explicitly describes the analysis of suitability and the optimal choice of shifted-Laplace preconditioners as an ongoing field of study, with incomplete understanding of spectral and pseudospectral behavior for highly heterogeneous $Au=f$3 and high-frequency regimes (Kazemi, 2017).

A second limitation is scalability at very high frequency. The real-shifted coarse-grid correction work argues that standard multigrid does not converge for high-frequency Helmholtz problems, that the common cure is a complex shift, and that this cure itself prevents scalability. Its alternative is to correct numerical dispersion by real-shifting only the coarsest Galerkin operator, not by globally complex-shifting the original problem. The reported results show wavenumber independent convergence for heterogeneous geophysical media in 2D and 3D when there are 12 grid points per wavelength on the fine grid, and a convergent cycle with very few iterations for 11 grid points per wavelength (Yovel et al., 21 Apr 2026).

A third limitation is that complex shifting alone does not remove coarse-grid dispersion mismatch. A two-grid finite-element solver with dispersion matching uses a domain-decomposition smoother applied to a complex-shifted Helmholtz operator

$Au=f$4

but its local Fourier analysis concludes that good convergence of the two-grid method still requires the fine- and coarse-level dispersion relations to closely match. This suggests that the complex shift is effective as a smoother ingredient, yet not sufficient by itself to guarantee optimal coarse-grid correction (Stolk, 22 Sep 2025).

These limitations motivate several research directions already present in the literature. One direction reduces the amount of required shift through improved smoothers such as additive Vanka patches (Yovel et al., 20 Nov 2025). Another direction avoids global complex shifting through real-shifted coarse-grid correction (Yovel et al., 21 Apr 2026). A third direction embeds the shift in nested or split iterations, such as HSS, so that multigrid sees a strongly damped inner operator while the outer preconditioner remains robust with respect to $Au=f$5 (Cotter et al., 23 Jun 2025). A plausible implication is that the complex-shifted Helmholtz operator is no longer best regarded as a single preconditioner, but rather as a design principle that reappears in spectral filtering, multilevel relaxation, contour quadrature, compression, and boundary-integral regularization.