Plane-Wave Discontinuous Galerkin Method

• The plane-wave discontinuous Galerkin method is a Trefftz-type discretization that uses exact local plane-wave solutions to satisfy the underlying homogeneous PDE, enhancing computational efficiency.
• It employs weakly enforced interelement and boundary conditions via DG fluxes, which concentrates the coupling on faces and supports flexible mesh designs.
• Recent extensions include DtN coupling, generalized and anisotropic bases, hybridizations with FEM and spline spaces, and neural network approaches to improve accuracy and condition numbers in complex settings.

The plane-wave discontinuous Galerkin method (PWDG) is a Trefftz-type discontinuous Galerkin discretization for time-harmonic wave problems in which the local approximation space on each element is spanned by plane waves that satisfy the underlying homogeneous partial differential equation exactly inside the element. Global conformity is not imposed strongly; instead, interelement transmission conditions and boundary conditions are enforced weakly through DG fluxes on the mesh skeleton. In the Helmholtz setting, this replaces the global $H^1(\Omega)$ trial space by a broken plane-wave Trefftz space and leads to formulations whose essential coupling is concentrated on faces, an architecture that has been developed for acoustic scattering, Maxwell systems, periodic diffraction problems, and several hybrid wave-based discretizations (Congreve et al., 2017, Kapita et al., 2016)

1. Core mathematical structure

In a standard two-dimensional Helmholtz model, one seeks $u$ on a bounded Lipschitz domain $\Omega\subset\mathbb{R}^2$ such that

$$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$

with wave number $k>0$. The corresponding weak formulation is well posed, but PWDG replaces the conforming space by a piecewise Trefftz space and shifts the coupling to jumps and averages on the mesh skeleton (Congreve et al., 2017).

Let $\mathcal{T}_h$ be a partition of $\Omega$ into convex polygonal elements $K$. On each element, the local plane-wave space of degree $p$ is

$$\mathrm{PW}_p(K)=\left\{v(x)=\sum_{j=1}^p \alpha_j \exp\big(ik\mathbf{d}_j\cdot(x-x_K)\big): \alpha_j\in\mathbb{C}\right\},$$

where $u$0 is the mass center of $u$1, the directions $u$2 are unit vectors, and in the numerical investigations they are chosen uniformly on the unit circle. Every basis function is an exact local solution of the homogeneous Helmholtz equation, so every $u$3 satisfies the Trefftz property

$u$4

The global discrete space

$u$5

is in general discontinuous across interfaces (Congreve et al., 2017).

Across an interior face, standard DG averages and jumps are used, for example

$u$6

The resulting PWDG sesquilinear form involves traces of $u$7 and $u$8 on faces, together with penalty-type terms on jumps of the field and its gradient. A defining structural feature is that volume integrals of $u$9 do not appear, because the local basis is Trefftz. This face-based construction supports general polygonal meshes and is closely connected to the ultra-weak variational formulation tradition in Trefftz discretization (Congreve et al., 2017, Kretzschmar et al., 2013).

2. Fluxes, artificial boundaries, and periodicity

The DG character of PWDG is expressed through numerical fluxes. For interior edges in Helmholtz discretizations, typical fluxes use averages of the field and gradient plus penalties on jumps. On Dirichlet boundaries, the boundary data are enforced weakly. In exterior scattering, this flux machinery is extended to artificial truncation boundaries, where the main difficulty is the faithful representation of outgoing waves (Kapita et al., 2016).

A prominent construction is the DtN-PWDG method for acoustic scattering outside a bounded obstacle. The unbounded exterior is truncated by a circular artificial boundary $\Omega\subset\mathbb{R}^2$0, and a Dirichlet-to-Neumann map $\Omega\subset\mathbb{R}^2$1 is introduced through a Fourier–Hankel expansion. In computation, the truncated map

$\Omega\subset\mathbb{R}^2$2

is inserted directly into new consistent numerical fluxes on $\Omega\subset\mathbb{R}^2$3. The corresponding formulation yields error estimates with respect to both the truncation order $\Omega\subset\mathbb{R}^2$4 and the mesh width, and the reported numerical results show that DtN boundary conditions can significantly improve accuracy relative to first-order absorbing conditions, especially in shadow regions and in multi-obstacle or resonant configurations (Kapita et al., 2016).

Periodic scattering requires a different boundary architecture. In periodic diffraction-grating problems, the computational domain is reduced to a single quasi-periodic cell, with quasi-periodic conditions on the lateral sides and DtN operators on top and bottom truncation boundaries. The Trefftz DG formulation again uses plane-wave discrete spaces, but now the boundary operators are mode-wise in the quasi-periodic Fourier basis. For polygonal meshes, all linear-system entries can be computed analytically. Using a Rellich identity, an explicit stability estimate is proved that is robust in the small material jump limit, which is a notable feature in periodic penetrable–impenetrable configurations (Moiola et al., 29 May 2025).

3. Approximation spaces beyond standard elementwise plane waves

Classical PWDG relies on exact local plane waves, which is natural when the local operator has constant coefficients. When coefficients vary smoothly, exact plane waves are no longer local Trefftz functions. One response is the generalized plane wave construction, in which local basis functions have the form

$\Omega\subset\mathbb{R}^2$5

with $\Omega\subset\mathbb{R}^2$6 a complex polynomial chosen so that $\Omega\subset\mathbb{R}^2$7 vanishes to high order at the element centroid. This produces a modified Trefftz DG scheme for piecewise smooth coefficients that retains the high-order convergence behavior of the original TDG framework and keeps the same number of degrees of freedom per element as standard plane-wave approximations. The trade-off is explicit: the skeleton-only advantage is lost, because additional volume residual stabilization terms are required (Imbert-Gerard et al., 2015).

Anisotropic media require a different generalization. For the Helmholtz equation and time-harmonic Maxwell equations in three-dimensional anisotropic media, novel plane-wave bases are obtained by scaling and coordinate transformations built from the spectral decomposition of the coefficient matrices. The transformed problem becomes isotropic, standard plane waves are defined there, and the basis is pulled back to the physical domain. The resulting PWDG analysis yields error estimates with explicit dependence on the condition number of the coefficient matrices, under a new shape-regularity assumption imposed on the transformed mesh rather than on the physical mesh. The reported numerical results verify the theory and indicate high accuracy (Yuan et al., 2019).

A more algebraic route is the embedded Trefftz DG method. Here a standard DG method is projected onto a Trefftz-type subspace through an embedding operator constructed from the kernel of a local operator matrix, without explicitly computing Trefftz basis functions. In the simplest cases this recovers established Trefftz DG methods, while also permitting extensions to inhomogeneous sources and non-constant coefficient operators. For Helmholtz, the reported numerical behavior shows improved accuracy similar to Trefftz DG methods based on plane waves (Lehrenfeld et al., 2022).

A distinct but related extension is the space-time Trefftz DG formulation of wave propagation. In one dimension it uses transport-polynomial bases that satisfy the PDE exactly on space-time elements, and for higher-dimensional Maxwell systems the authors propose plane-wave-like space-time bases. The resulting method exhibits spectral convergence in the $\Omega\subset\mathbb{R}^2$8-norm over the entire space-time domain, not just in space at fixed time, which places time-domain Trefftz discretization in direct conceptual continuity with frequency-domain PWDG (Kretzschmar et al., 2013).

4. Conditioning, basis orthogonalization, and iterative solution

Conditioning is a central practical issue in PWDG. Congreve, Gedicke, and Perugia carried out a numerical investigation of both local basis conditioning and global system conditioning for the Helmholtz problem. On a single element, the spectral condition number of the local plane-wave mass matrix depends algebraically on mesh size and wave number and exponentially on the number of plane-wave directions; it also depends on element shape. For regular polygons, conditioning improves as the number of sides increases, whereas for anisotropic rectangles the condition number increases exponentially with the aspect ratio. On a square element, extensive numerical experiments support an empirical model of the form

$\Omega\subset\mathbb{R}^2$9

capturing algebraic dependence on $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$0 and exponential dependence on $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$1 (Congreve et al., 2017).

The same study shows that the conditioning of the assembled global system can be improved by orthogonalization of the local basis functions with the modified Gram–Schmidt algorithm. The orthogonalization is performed elementwise with respect to the Hermitian part of the local system matrix,

$$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$2

leading to a block-diagonal basis transformation $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$3 and a transformed global matrix $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$4. The paper stresses that this is a basis change, not a preconditioner. Numerically, it reduces global condition numbers and produces significantly fewer GMRES iterations: on the quadrilateral mesh with $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$5, the GMRES count for the original basis rises from 60 at $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$6 to 217 at $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$7, whereas for the orthogonalized basis it rises from 47 to 73 over the same range. The effective range of this procedure is limited by finite precision: in double precision it remains effective up to about $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$8, and in single precision up to about $$-\Delta u-k^2u=0 \quad \text{in }\Omega,\qquad \nabla u\cdot \mathbf{n}+iku=g \quad \text{on }\partial\Omega,$$9 (Congreve et al., 2017).

These observations undercut two common simplifications. First, increasing the number of plane-wave directions is not purely an approximation issue; it is also a severe conditioning issue. Second, element regularity matters even in Trefftz schemes: nearly round polygonal elements are beneficial, while strong anisotropy is harmful at the local basis level (Congreve et al., 2017).

5. Hybridizations and application domains

PWDG ideas have been coupled systematically with other discretizations. One example is the nonconforming coupling of a wave-based discontinuous Galerkin method with a standard finite element model for the Helmholtz equation. The coupling is formulated without Lagrange multipliers, allows incompatible meshes across the interface, and is reported to be optimal in the sense that the convergence rates of the individual methods are maintained. The intended use is geometric decomposition: standard FEM is employed near complex surfaces and geometric details, while the wave-based DG discretization covers large, simply shaped regions where plane-wave approximations are more efficient (Gaborit et al., 2018).

A second major application area is periodic full-potential electronic structure. In one formulation, the periodic cell is decomposed into atomic spheres and an interstitial region. Plane waves are used in the interstitial region, radial basis functions times spherical harmonics are used inside the atomic spheres, and the two are coupled by a symmetric interior penalty DG method. The scheme follows the philosophy of APW/LAPW methods while possessing systematically spectral convergence rates, and the paper provides rigorous a priori error analysis for the corresponding linear eigenvalue problems (Li et al., 2019).

A more recent variant replaces radial–spherical expansions by tensor-product B-splines in localized atomic patches and couples them to interstitial plane waves through a symmetric interior penalty DG formulation. This is motivated by multiscale regularity in full-potential calculations: orbitals exhibit cusp singularities near nuclei and are smooth in the interstitial region. The method develops a combined fast Fourier transform and Chebyshev correction strategy for restricted plane-wave integrals, constructs a trace-block DG preconditioner to alleviate conditioning deterioration caused by SIPG penalty terms, and proves a priori error estimates showing algebraic convergence near the nuclei and superalgebraic convergence in the interstitial region (Lai et al., 24 May 2026).

Taken together, these hybrid formulations show that plane-wave DG is not restricted to a homogeneous wave field over a single uniform discretization. It also functions as a coupling technology for multiscale decompositions in which oscillatory regions are described by plane waves and singular or geometrically intricate regions by alternative local bases (Gaborit et al., 2018, Lai et al., 24 May 2026).

6. Relation to adjacent DG methods and recent directions

The label “plane-wave DG” is sometimes used loosely. A useful counterexample is the weight-adjusted discontinuous Galerkin method for the poroelastic wave equation. That method is high-order DG and is tested against analytical plane-wave solutions, but its trial and test spaces are polynomial $k>0$0 spaces on simplices, not plane-wave Trefftz spaces. The paper itself explicitly notes that it is not a PWDG in the classical sense: plane waves are used only as analytical benchmarks, not as basis functions. This distinction matters because flux design, heterogeneity treatment, and conditioning mechanisms are different in polynomial DG and in genuine plane-wave Trefftz DG (Shukla et al., 2019).

Recent work pushes the PWDG idea toward adaptive nonlinear approximants. A discontinuous Galerkin plane wave neural network method has been proposed for Helmholtz and Maxwell equations. It defines an elliptic-type variational problem as in the plane wave least square method with $k>0$1-refinement and introduces adaptive recursively augmented discontinuous Galerkin subspaces whose basis functions are element-wise neural network realizations with $k>0$2-refinement. The activation function is chosen as a complex-valued exponential function like the plane wave function, and convergence is established without assuming boundedness of the neural-network parameters (Yuan et al., 11 Jun 2025).

This suggests a broad but still technically coherent picture. In its classical form, PWDG is a Trefftz DG method with exact elementwise plane waves. In contemporary variants, that core idea is being extended in several directions: DtN coupling for exterior and periodic problems, hybridization with finite elements and spline spaces, transformed bases for anisotropy, algebraic embeddings for non-constant operators, and neural constructions that preserve the local oscillatory ansatz while relaxing the rigid linear basis structure (Moiola et al., 29 May 2025, Yuan et al., 11 Jun 2025)