🔢 mathematics

A Hierarchical direct solver for Lippmann-Schwinger equations based on boundary traces

This paper presents a high-order hierarchical direct solver for two-dimensional Lippmann-Schwinger equations that constructs a reusable scattering operator by operating entirely on boundary traces above the leaf level, achieving $O(N^{3/2})$ setup cost and $O(N)$ application cost for variable-media scattering problems that challenge iterative methods.

Original author: Yunhui Cai, Simon Dirckx, Joar Bagge, Per-Gunnar Martinsson

Published 2026-11-03

📖 6 min read🧠 Deep dive

CC BY 4.0

Original author: Yunhui Cai, Simon Dirckx, Joar Bagge, Per-Gunnar Martinsson

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the author. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine sound or light traveling through a medium that is not uniform, perhaps passing through a cloud of varying density or a lens with changing thickness. In physics, predicting how these waves bend, scatter, and emerge is a fundamental challenge. The equations that describe this behavior are well known, but solving them for complex, real-world shapes is notoriously difficult. When the waves encounter an object with a complicated internal structure, the standard mathematical approach requires breaking the entire space into a massive grid of tiny points. To find the answer, a computer must solve a giant system of equations where every single point is connected to every other point. This creates a computational bottleneck: if the object changes, or if the incoming wave comes from a different angle, the entire massive calculation must be started over from scratch. For scientists studying everything from medical imaging to the design of advanced optical materials, this repeated, heavy lifting has long been a barrier to efficiency.

A team of researchers at the University of Texas at Austin and the KTH Royal Institute of Technology has developed a new way to handle this problem. They have created a high-precision method that acts like a reusable map for how waves interact with a specific object. Instead of recalculating the entire interior of an object every time a new wave hits it, their approach builds a single, comprehensive description of the object's scattering behavior. Once this description is constructed, it can be used instantly to predict the outcome for any number of incoming waves, regardless of their direction or shape. The method is designed for two-dimensional scenarios, such as waves moving across a flat surface, and it achieves a level of accuracy that allows it to handle difficult cases, including objects with sharp edges, high-frequency waves, and complex internal cavities, which often confuse older methods.

The core idea behind this new solver is to stop thinking about the problem as a single, monolithic block of data and instead view it as a hierarchy of smaller, manageable pieces. The researchers divide the area containing the object into a grid of small squares, which are then grouped into larger squares, and those into even larger ones, forming a tree-like structure. At the very bottom of this tree, the smallest squares are treated individually. For each tiny square, the researchers solve the wave equation locally to determine how that specific piece responds to waves entering its boundaries. They do not keep track of the wave's behavior at every single point inside the square; instead, they record only how the wave enters the square and how it leaves it. This creates a compact "scattering operator" for that small piece, which acts as a rulebook for its interaction with the outside world.

Once these rules are established for the smallest pieces, the researchers begin to merge them. They take two neighboring small squares and combine their individual rulebooks to create a single rulebook for the larger square they form together. This merging process is clever: it calculates how the two pieces talk to each other across their shared boundary without ever needing to look inside them again. The new, larger rulebook describes how the combined area scatters waves based solely on what happens at its outer edges. This process repeats up the tree, merging pairs of boxes into larger and larger ones, until the entire object is represented by a single, master rulebook at the top of the hierarchy. The entire construction relies on the fact that the complex interactions inside the object can be fully summarized by what happens at the boundaries between these regions.

The result is a powerful tool that separates the work into two distinct phases. The first phase is the "build" stage, where the researchers construct the master rulebook for a specific object and frequency. This stage is computationally intensive and takes time, scaling with the size of the problem in a way that is manageable for modern computers. However, once this rulebook is built, the second phase, the "solve" stage, becomes incredibly fast. If a scientist wants to know what happens when a wave comes from a new angle, they do not need to rebuild the model or re-solve the interior equations. They simply feed the new wave into the existing rulebook, and the answer emerges almost instantly. The researchers demonstrated that for a problem with millions of data points, the initial construction might take a few minutes, but processing hundreds of different incoming waves afterward takes only seconds.

In their tests, the team applied this method to a variety of challenging scenarios. They simulated waves hitting a smooth, rounded bump, a sharp-edged pyramid, a hollow cavity, and even a complex grid of tiny scatterers resembling a photonic crystal. In every case, the method produced highly accurate results, matching the precision of much slower, traditional techniques. The accuracy remained high even when the waves had very short wavelengths or when the object had irregular shapes that caused waves to bounce back and forth in complex ways. The researchers noted that while the method works best for smooth objects, it still performs well for objects with sharp corners, though the convergence to the exact answer is slightly slower in those cases.

The significance of this work lies in its ability to turn a problem that usually requires a fresh, heavy calculation for every new scenario into a one-time investment with unlimited reuse. By reducing the problem to boundary interactions and organizing the solution hierarchically, the researchers have found a way to bypass the computational wall that has limited the study of variable-media scattering. This approach allows for the rapid exploration of how waves behave in complex environments, opening the door to more efficient designs in optics, acoustics, and imaging. The method does not rely on approximations that sacrifice accuracy for speed; rather, it achieves speed by organizing the mathematics in a way that respects the physical structure of the problem, ensuring that the heavy lifting is done only once, while the rest of the work flows effortlessly through the established hierarchy.

Technical Summary: A Hierarchical Direct Solver for Lippmann–Schwinger Equations Based on Boundary Traces

Problem Statement
The paper addresses the numerical solution of the variable-coefficient Helmholtz equation in two dimensions, modeling acoustic scattering from inhomogeneous media. The problem is formulated as a Lippmann–Schwinger integral equation:
$q(\mathbf{x}) + \kappa^2 b(\mathbf{x}) \int_{\Omega} G_\kappa(\mathbf{x}, \mathbf{y}) q(\mathbf{y}) d\mathbf{y} = -\kappa^2 b(\mathbf{x}) u_{in}(\mathbf{x}), \quad \mathbf{x} \in \Omega$
where $q$ is the source density, $\kappa$ is the wavenumber, $b$ is a compactly supported scattering potential, and $u_{in}$ is an incoming field.

Standard high-order discretizations of this equation result in dense linear systems. While iterative solvers accelerated by Fast Multipole Methods (FMM) or FFTs are common, they require a new global solve for every distinct incident field. This becomes computationally prohibitive for applications involving many incident fields (e.g., varying angles of incidence) or regimes where global iteration converges slowly (e.g., high frequencies, strong scattering, or limited regularity). The authors aim to construct a reusable direct solver that pre-computes a scattering operator for the medium, allowing subsequent incident fields to be processed rapidly without revisiting the volume discretization.

Methodology
The proposed method is a high-order hierarchical direct solver that operates by constructing a field-to-field scattering operator. The core innovation is a dimensional reduction where volume unknowns appear only at the finest level of a hierarchical tree; all higher-level operations are performed entirely on boundary traces.

Hierarchical Domain Decomposition: The computational domain $\Omega$ is partitioned into a binary tree of square or rectangular subdomains (leaf boxes).
Continuum Formulation (Field-to-Field):
- Local Operators: For each leaf box, a local scattering operator is defined as a map from an incoming trace (boundary values of the incoming field) to an outgoing trace (boundary values of the scattered field). This is constructed by solving a local Lippmann–Schwinger problem within the leaf.
- Merging: The scattering operator for a parent box is constructed by merging the operators of its two children. This involves:
  - Continuation: Extending the incoming field from the parent boundary to the children's boundaries.
  - Interaction Resolution: Solving a system that accounts for multiple scattering between the children using a solution-transfer operator $T_\tau = (I - E_\tau F)^{-1} U_\tau$ , where $E_\tau$ propagates fields between siblings and $F$ represents the block-diagonal child operators.
  - Assembly: Combining the outgoing fields of the children to form the parent's outgoing field.
- Crucially, the merging process relies only on boundary-to-boundary operators ( $U_\tau, E_\tau, R_\tau$ ) and does not require re-solving the interior volume equations of the children.
Discretization:
- Leaves: High-order local volume solves are performed on tensor-product Chebyshev grids using corrected quadrature rules for the Lippmann–Schwinger equation. Boundary traces are represented via collocation on Gauss–Legendre nodes.
- Merging: The universal operators ( $U, E, R$ ) are constructed using high-order boundary integral equation (BIE) techniques (specifically Brakhage–Werner formulations) to handle the propagation of Helmholtz fields between boundaries.
Algorithm Flow:
- Upward Pass: Construct local leaf scattering operators and merge them recursively up the tree to form the global scattering operator for the entire domain.
- Downward Pass: For a specific incident field, propagate the incoming trace down the tree to determine the incoming traces on all leaves.
- Recovery: Solve the local volume problems on the leaves to recover the full interior field if required.

Key Contributions

Continuum-First Derivation: Unlike previous hierarchical methods that discretize the volume integral first and then exploit matrix rank structure, this work derives the hierarchical reduction at the continuum level. The discretization is introduced only after the field-to-field scattering operators are defined. This distinction is highlighted as critical for maintaining high-order accuracy, where discrete rank structures can become large.
Boundary Trace Reduction: The method ensures that volume unknowns are confined to the leaf level. All interactions between subdomains at higher hierarchy levels are represented solely by boundary traces, significantly reducing the dimensionality of the merge operations.
Reusable Direct Solver: The construction of the scattering operator is independent of the incident field. Once the operator is built for a specific medium and frequency, it can be applied to any number of incoming fields with linear or near-linear cost relative to the number of degrees of freedom.
High-Order Accuracy: The solver employs high-order polynomial representations ( $p=10$ or $p=20$ ) and spectral methods, achieving high accuracy even for challenging problems with cavities, high wavenumbers, and limited regularity.

Results
Numerical experiments were conducted on six scattering potentials (Gaussian bump, smooth mollifier, cavity, lens, photonic crystal, and pyramid) with wavenumbers ranging from $\kappa=25$ to $\kappa=200$ .

Accuracy: The method demonstrates high-order convergence. For smooth potentials, errors reach $O(10^{-10})$ at fine resolutions. Potentials with limited regularity (e.g., the pyramid case) or high-frequency features (lens, photonic crystal) require finer discretizations but still converge to high accuracy.
Complexity:
- Build Time: Constructing the hierarchical scattering representation requires $O(N^{3/2})$ work and $O(N \log N)$ storage (where $N$ is the total degrees of freedom).
- Evaluation Time: Once built, processing a new incident field requires $O(N)$ work if only the outgoing field is needed, and $O(N \log N)$ if the full interior field is recovered.
Multiple Right-Hand Sides: The solver efficiently handles multiple incident fields. In a test with 500 incidence angles, the time to solve for the interior field was comparable to a single solve, and the exterior field could be computed orders of magnitude faster using a precomputed exterior evaluation matrix.
Memory: The full hierarchy requires significant memory ( $O(N \log N)$ ), but a variant storing only the root operator reduces storage to $O(N)$ .

Significance and Claims
The authors claim that this work provides a robust, high-order accurate direct solver for 2D variable-media scattering that overcomes the limitations of iterative methods in specific regimes (high frequency, strong scattering, multiple sources). By separating the computation into a reusable precomputation stage and a rapid application stage, the method is particularly suited for inverse scattering and parametric studies where the medium is fixed but the incident field varies.

The paper modestly notes that the current implementation does not yet exploit the internal low-rank or Hierarchically Semiseparable (HSS) structures present in the dense merge matrices. The authors identify this as a primary avenue for future acceleration, which could further reduce the $O(N^{3/2})$ build complexity and storage requirements. The current work establishes the theoretical framework and a high-order numerical realization of the boundary-trace hierarchical approach.

Technical Summary: A Hierarchical Direct Solver for Lippmann–Schwinger Equations Based on Boundary Traces

More like this