High-Order Accurate Scattering Matrix

• High-Order Accurate Scattering Matrix is a computational construct that maps incoming to outgoing wave coefficients with uniform high-order accuracy even in complex, singular geometries.
• It employs techniques like operator regularization, graded change-of-variable, and adaptive quadrature to resolve singularities and achieve spectral convergence.
• Hierarchical matrix compression and multipole expansions enable efficient simulation in applications ranging from computational electromagnetics to condensed matter physics.

A high-order accurate scattering matrix is a computational or analytical construct that encodes the exact mapping between incoming and outgoing wave coefficients in scattering phenomena, with accuracy and convergence properties that persist uniformly even in the presence of complex microstructure, geometric singularities, and multi-particle interactions. Such matrices serve a central role in fields ranging from computational electromagnetics to condensed matter physics, as they provide a compact summary of the full wave response of a system while supporting highly efficient simulation and analysis.

1. Mathematical Definition and Construction

A scattering matrix (S-matrix) relates the amplitudes of incoming modes (such as spherical harmonics, partial waves, or waveguide eigenmodes) to those of outgoing modes after an incident field interacts with a scatterer. For electromagnetic scattering, this takes the general form

$$\mathbf{b}_{\text{out}} = S \cdot \mathbf{a}_{\text{in}}$$

where $\mathbf{a}_{\text{in}}$ are incoming expansion coefficients, $\mathbf{b}_{\text{out}}$ are scattered (outgoing) coefficients, and $S$ is the scattering matrix.

High-order accuracy is pursued through:

• High-order quadrature for numerical integration of the surface/volume integral equations (e.g., using piecewise linear or higher order polynomials on curved patches [(Gimbutas et al., 2011); (Bremer et al., 2013); (Chorsi et al., 2016)]).
• Representation of solutions in rapidly convergent multipole expansions (e.g., truncating spherical harmonics at sufficiently high degree $p$, yielding expansions of size $(p+1)^2$ for 3D problems (Gimbutas et al., 2011)).
• Regularization of hypersingular integral operators by operator-based techniques, ensuring that numerical discretization preserves spectral convergence even when the underlying density has singularities at corners or edges (Sideris et al., 30 Jan 2025).
• Systematic precomputation of the mapping from any incoming expansion to the corresponding outgoing (scattered) expansion, not only for simple geometries (using Debye–Mie theory), but for arbitrarily shaped inclusions via integral equation solution followed by projection onto spherical harmonic bases [(Gimbutas et al., 2011); (Demésy et al., 2018)].

2. Numerical Strategies for High-Order Accuracy

High-order accurate scattering matrices are underpinned by:

• Operator regularization and change of variables: When singularities (e.g., corner blowup of the density) are present, a graded change of variables combined with a change of unknowns (absorption of the known singularity into the representation) renders the transformed density smooth, allowing high-order quadrature to remain effective (Sideris et al., 30 Jan 2025).
• Hierarchical matrix compression: Hierarchically Block Separable (HBS) and related formats, as well as butterfly and HODLR structures, expose and exploit low-rank structure in the off-diagonal coupling blocks of discretized integral operators, reducing the cost of S-matrix construction to O(N^{3/2}) or better while preserving the formal order of accuracy [(Gopal et al., 2020); (Bremer et al., 2013); (Liu et al., 2016)].
• Adaptive quadrature and discretization: Adaptive quad-trees, high-order Nystrom discretizations, and Chebyshev expansions permit local adaptation of discretization order to the regularity of the solution, ensuring uniform convergence even for heterogeneous inclusions or domains with multi-scale features (Ambikasaran et al., 2015, Sideris et al., 30 Jan 2025).

A representative formula for the regularized combined field integral equation (CFIE–R) for the scattered field is

$$u = -i\eta\int_\Gamma G_k(r,r')\,\phi(r')\,d\ell' + \int_\Gamma \frac{\partial G_k(r,r')}{\partial n(r')}\,R[\phi(r')]\,d\ell'$$

with $$R[\phi](r') = \int_\Gamma G_K(r',r'')\,\phi(r'')\,d\ell(r'')$$, allowing high-order numerical approximation even when $\phi$ diverges near a corner (Sideris et al., 30 Jan 2025).

3. Handling Geometric Singularities and Field Blowup

Field enhancement and the divergence of solution densities at corners and edges of conductors present severe obstacles to high-order schemes.

• Change of variable (CoV): For a patch containing a singular corner, a mapping $w(\theta)$ is employed so that the discretization nodes become clustered with the Jacobian and its first $p-1$ derivatives vanishing at the corner, transforming a singular density $\phi(s)$ into a smooth function $$\psi_q(\theta) = \phi(s_q(\theta))\,\widetilde{L}_q(\theta)$$ which is highly amenable to Chebyshev expansion and quadrature (Sideris et al., 30 Jan 2025).
• Patchwise parameterization: The boundary is partitioned such that each "corner-carrying" patch is parameterized to bring the corner to a fixed location, standardizing mesh grading and enabling integration without local reanalysis.
• Chebyshev and adaptive quadrature: Integration weights for near-singular interactions are accurately computed by projecting densities onto Chebyshev polynomials and precomputing kernel integrals on reference domains, yielding stable and accurate evaluation even in the vicinity of blowup (Sideris et al., 30 Jan 2025).

This dual strategy of regularization and high-order adaptive quadrature ensures that even for divergence rates as severe as $d^{-\nu}$ near corners, machine-precision solutions are attainable at minimal computational expense.

4. Generalization to Three Dimensions and Maxwell Problems

The described methodologies are directly extensible to full 3D Maxwell boundary integral problems, where tangential current densities typically diverge at edges:

• The same change-of-unknown and graded CoV regularizes the density associated with diverging edge-parallel current components, so that high-order quadrature remains valid even on Lipschitz or polyhedral non-smooth geometries.
• Operator regularization in the Maxwell context is expected to maintain favorable spectral properties and enable S-matrix construction with O(N^{3/2}) or better scaling, provided the kernel admits a suitable representation (Sideris et al., 30 Jan 2025).
• Unlike 2D cases where certain singularities can sometimes be avoided by choice of formulation, in 3D the tangential blowup is inherent and must be addressed via these regularization-based strategies.

5. Numerical Performance and Robustness

Numerical experiments validate the sharpness and robustness of these high-order accurate scattering matrix schemes:

• Relative errors as low as $10^{-12}$ are reported even for observation points at $10^{-8}$ distance from non-smooth corners, with convergence exceeding 9th- or 10th-order.
• In a comparison between regularized (CFIE–R) and unregularized formulations (such as MFIE), only the former achieves rapid convergence universally, as MFIE stagnates or exhibits resonance-related instability (Sideris et al., 30 Jan 2025).
• The eigenvalue spectra for regularized operators are tightly clustered and bounded, resulting in rapid convergence for Krylov methods (e.g., typically < 15 GMRES iterations even near resonance).
• Omission of any core element (change-of-unknown, precise quadrature, or mesh grading) causes a catastrophic loss of convergence order and attainable accuracy, underscoring the necessity of each component.

6. Practical Implications and Applications

The practical impact of high-order accurate scattering matrix construction with singularity-resolving capabilities is profound:

• It leads to simulation tools capable of tackling microstructured materials with thousands to millions of inclusions, arbitrary geometrical complexity, and material heterogeneity while maintaining efficiency and accuracy (Gimbutas et al., 2011).
• The techniques obviate the need for a priori basis enrichment or ad hoc local mesh refinement around corners—previously standard strategies for integrating over non-smooth geometries.
• Because the S-matrix is compact and universally applicable, it can be precomputed for inclusions of arbitrary shape and reused in reduced-order or large-scale multiple scattering computations.
• Extension to Maxwell's equations in 3D is immediate in form, promising broad applicability in computational photonics, antenna analysis, metamaterial design, and electromagnetic compatibility studies.

7. Summary Table: Key Elements for High-Order Accurate Scattering Matrices

Technique	Purpose	Effect on S-matrix Construction
Operator regularization	Smooths hypersingular terms	Well-conditioned, resonance-free system
Graded change of variable	Mesh clustering, singularity absorption	Bounded density, supports high-order rule
Change of unknown	Transforms unbounded densities to smooth	Enables spectral Chebyshev expansion
Chebyshev/adaptive quadrature	Precision-preserving integration	Machine-precision near singularities
Hierarchical matrix compression	Efficient storage and inversion	O(N)–O(N^{3/2}) scaling, reusability
Patchwise parameterization	Standardizes singularity location	Uniform application of graded coV

The collective adoption of these methods, as established in (Sideris et al., 30 Jan 2025), establishes a foundation for robust, high-order accurate computation of scattering matrices for a broad class of singular and non-smooth geometries, with direct generalization to full 3D Maxwell problems and wide relevance across computational science and engineering disciplines.