---
abstract: 'Approximations by Trefftz functions are rapidly gaining popularity in the numerical solution of boundary value problems of mathematical physics. By definition, these functions satisfy locally, in weak form, the underlying differential equations of the problem, which often results in high-order or even exponential accuracy with respect to the size of the basis set. We highlight two separate examples in applied electromagnetics and photonics: (i) homogenization of periodic structures, and (ii) numerical simulation of electromagnetic waves in slab geometries. Extensive numerical evidence and theoretical considerations show that Trefftz approximations can be applied much more broadly than is traditionally done: they are effective not only in physically homogeneous regions but also in complex inhomogeneous ones. Two mechanisms underlying the high accuracy of Trefftz approximations in such complex cases are pointed out. The first one is related to trigonometric interpolation and the second one – somewhat surprisingly – to well-posedness of random matrices.'
address:
- |
    Department of Electrical and Computer Engineering, The University of Akron\
    Akron, OH 44325-3904, USA\
    igor@uakron.edu
- |
     School of Physical & Mathematical Sciences, Nanyang Technological University\
    21 Nanyang Link, Singapore 637371\
    yidong@ntu.edu.sg
- |
     Department of Radiology, University of Pennsylvania, Philadelphia, PA 19104, USA\
    vmarkel@pennmedicine.upenn.edu 
author:
- Igor Tsukerman
- 'Shampy Mansha, Y. D. Chong'
- 'Vadim A. Markel'
bibliography:
- 'Igor\_reference\_dbase.bib'
- 'flame-slab-references.bib'
title: |
    Trefftz Approximations in Complex Media:\
    Accuracy and Applications
---

Trefftz approximations ,convergence ,Maxwell’s equations ,homogenization ,photonic devices ,wave scattering ,interpolation ,finite difference schemes ,random matrices 65M06 ,76M20 ,35B27 ,76M50 ,74Q15 ,35Q61 ,15B52

Introduction {#sec:Intro}
============

Many classical numerical methods for partial differential equations rely on polynomial or piecewise-polynomial approximations of the solution. Examples include traditional finite difference (FD) schemes, the finite element method (FEM), and the boundary element method (BEM). But a strong incentive to achieve qualitatively higher accuracy of the numerical solution has led, over several decades of research, to the development of Trefftz-based methods. By definition, Trefftz functions satisfy locally (in weak form) the underlying differential equations of the problem, which often results in high-order algebraic or even exponential convergence with respect to the dimension of the basis. This qualitative accuracy improvement has been demonstrated in a large variety of mathematical methods and engineering applications: Domain Decomposition [@Herrera00; @Farhat-DD-DG09], Generalized FEM [@Melenk96; @Babuska97; @Babuska-GFEM2004; @Plaks03; @Proekt02; @Strouboulis-GFEM-Helmholtz2006], Discontinuous Galerkin [@Farhat-DD-DG09; @Cockburn00; @Arnold02; @Buffa-Monk-ultraweak08; @Gittelson09; @Gabard-wave-based-DG-UW-LS11; @Hiptmair-PWDG-Helmholtz2011; @Kretzschmar-IMA16], and finite difference (“Flexible Local Approximation MEthods,” FLAME) [@Tsukerman-JCP10; @Tsukerman05; @Tsukerman06; @Tsukerman-PBG08].

It is not our intention to review all, or even some, of these Trefftz-oriented methods; several good reviews are already available: [@Deckers-wave-method-overview2014; @Qin-Trefftz-FEM2005] and especially [@Hiptmair2016]. Rather, our focus is on one question central in these methods: *why are Trefftz approximations so effective*?

A simplified intuitive picture is shown in [Fig. \[fig:test-waves\]]{}, left panel. Several incident waves, schematically indicated with solid arrows, are impinging on an object (in general, physically inhomogeneous) and give rise to the respective total fields inside and to scattered fields outside that object. For visual clarity, only the incident waves are sketched in the figure, and their number is limited to three.

The total fields inside the scatterer, by definition, form a Trefftz set. One may view it as a “database” or “training set,” which can be precomputed and then used to approximate the field induced by another wave, indicated with a dashed arrow and a question mark in [Fig. \[fig:test-waves\]]{}. This approximation of one physically meaningful solution by other physically meaningful solutions (as opposed to, say, generic polynomials) certainly makes intuitive sense but is not trivial from the mathematical perspective.

\[fig:test-waves\] ![Several incident waves (schematically indicated with solid arrows) give rise to the respective total fields inside an inhomogeneous scatterer. These total fields, by definition, form a Trefftz set. This “training set” can be precomputed and then used to approximate the field induced by another wave, indicated with a dashed arrow and a question mark. For visual clarity, only three Trefftz waves are sketched, and only the incident components. Left: one inhomogeneous object is present. Right: additional scatterers (such as S1, S2, S3) may be present in the case of an unknown field (dashed arrow) whose approximation is sought within a given small subdomain $\Omega_h$.](test_waves.png "fig:"){width="0.4\linewidth"} 0.2in ![Several incident waves (schematically indicated with solid arrows) give rise to the respective total fields inside an inhomogeneous scatterer. These total fields, by definition, form a Trefftz set. This “training set” can be precomputed and then used to approximate the field induced by another wave, indicated with a dashed arrow and a question mark. For visual clarity, only three Trefftz waves are sketched, and only the incident components. Left: one inhomogeneous object is present. Right: additional scatterers (such as S1, S2, S3) may be present in the case of an unknown field (dashed arrow) whose approximation is sought within a given small subdomain $\Omega_h$.](test_waves_complex.png "fig:"){width="0.4\linewidth"}

The right panel of [Fig. \[fig:test-waves\]]{} illustrates a more interesting, and more complicated, case. Suppose that the Trefftz training set has been generated for the original inhomogeneous scatterer – same as in the left panel. However, the unknown “dashed arrow” solution *may involve additional objects* – such as S1, S2, S3 – in the computational domain. Obviously, under this complication, little can be inferred about the unknown solution from the Trefftz set *in the whole domain*, especially in the regions around the additional scatterers. One may hope, however, that the field *within a given small subdomain* $\Omega_h$ *inside the original scatterer* can still be approximated accurately as a superposition of the known Trefftz waves. This setup is the central issue of Sections \[sec:Auxiliary-basis\], \[sec:Random-matrices\], and is inspired by our numerical experiments with pseudorandom structures of [Section \[sec:FLAME-slab\]]{}, as well as by our earlier work on multiparticle problems [@Dai-Webb11].

The overall motivation for the paper is to highlight applications of Trefftz functions to problems involving complex, inhomogeneous media. Much of mathematical analysis so far has revolved around the homogeneous case (that is, equations with constant coefficients), where cylindrical, spherical or plane waves serve as Trefftz functions for the Helmholtz equation, while harmonic polynomials are used for the Laplace equation. One can refer, for example, to papers by Melenk, Hiptmair, Moiola, Perugia *et al.* cited above, to the references in these papers, and to Perrey-Debain’s paper [@Perrey-Debain-PW-convergence2006]. Much less attention has been paid to the inhomogeneous case [@Melenk95 Chapter IV], [@Melenk99 Section 3], [@Laghrouche-jumps2005; @Imbert-Gerard-interp-PW2015], which is substantially more complicated but at the same time more rewarding in practice.

For illustration, in Sections \[sec:Trefftz-homogenization\] and \[sec:FLAME-slab\] we consider two application examples where Trefftz approximations prove to be effective for two different variations of the generic setup shown in [Fig. \[fig:test-waves\]]{}. The first example is non-asymptotic and nonlocal two-scale homogenization. Instead of a single scatterer, in this case one deals with a periodic structure; Trefftz functions on the fine scale are Bloch waves traveling in different directions, and on the coarse scale – the corresponding plane waves.

The second example involves a common setup in metasurface and nanophotonics research: a patterned finite-thickness slab. This problem is especially challenging computationally when the pattern is non-periodic and the slab is geometrically large relative to the vacuum wavelength. One possible simulation procedure relies on high-order Trefftz difference schemes (FLAME). The Trefftz bases are computed “locally,” i.e. over relatively small segments of the structure ([Section \[sec:FLAME-slab\]]{}).

Sections \[sec:Trig-interpolation\] and \[sec:FD-Trefftz\] provide background information needed in the application examples of Sections \[sec:Trefftz-homogenization\] and \[sec:FLAME-slab\]. The underlying mechanisms for the accuracy of Trefftz approximations are discussed in Section \[sec:Trefftz-approximation\].

Preliminaries: Trigonometric Projection and Interpolation {#sec:Trig-interpolation}
=========================================================

Trigonometric approximation of periodic functions is a well-established subject. Here we summarize the key mathematical results that will be needed in Section \[sec:Interpolation-argument\].

For any Lipschitz-continuous periodic function $g$ on $[-\pi, \pi]$, one may consider its best possible approximation by a trigonometric polynomial $T_n$ in the maximum norm: $$\label{eqn:En-best-approx-trig-poly}
E_n^T(g) \,=\, \min_{\alpha, \beta} \max_{\phi \in [-\pi, \pi]}
\left| g(\phi) - T_n(\alpha, \beta, \phi) \right|, ~~~$$ where $$\label{eqn:trig-poly-defined}
T_n(\alpha, \beta, \phi) \, \equiv \, \alpha_0 + \sum_{\nu=1}^n 
(\alpha_\nu \cos \nu \phi + \beta_\nu \sin \nu \phi),$$ $$\alpha \equiv \{\alpha_0, \alpha_1, \ldots, \alpha_n\},
~~~ \beta \equiv \{\beta_1, \ldots, \beta_n\}$$ A slightly modified notation of [@Meinardus-Approximation-theory1967] is used here. Note that the total number of coefficients $\alpha$, $\beta$ in the trigonometric series is $N = 2n+1$.

It follows from Jackson’s theorem [@Jackson-Theory-approximation1930], or [@Meinardus-Approximation-theory1967 Theorem 41], that if the derivative $g^{(l+1)} (\phi)$ exists and is bounded, i.e. $$\label{g-bounded-derivative}
\left| g^{(l+1)} (\phi) \right|  \, \leq \, M_{l+1},
~~~ l = 0,1, \ldots$$ then $$\label{Eg-leq-nk1}
E_n^T(g) \, \leq\, \frac{c^{l+1} M_{l+1}} {n^{l+1}}, ~~~
c = 1 + \frac{\pi^2}{2}$$ For reasons that will become apparent in Section \[sec:Interpolation-argument\], we are interested primarily in trigonometric *interpolation* rather than the best approximation, and thus need to relate the two. The interpolant $\tilde{T}_N(\phi)$ of a given function $g(\phi)$ over a set of $N = 2n + 1$ equidistant knots $\{\phi_m \}$ is defined in a standard way, by requiring that $$\label{eqn:interpolant-defined}
   \tilde{T}_N(g, \phi_m) = g(\phi_m), ~~ \phi_m = \frac{2\pi m}{N},
   ~~~m = 0,1, \ldots, N-1$$ It is known that this interpolant exists and is unique. Furthermore, there is an upper bound for the interpolation error: $$\label{eqn:interp-error-vs-approximation-error}
   \| g - \tilde{T}_N(g) \|_\infty ~ \leq ~ (1 + \Lambda_N) \,
   \| g - T_N \|_\infty ~\equiv~ (1 + \Lambda_N) \, E_n^T(g)$$ where $\Lambda_N$ is the Lebesgue constant, which itself has an upper bound [@Cheney1975] $$\label{eqn:Lebesgue-const-upper-bound}
\Lambda_N \, \leq \, 2 \pi^{-1} \log N + \frac{5}{3}$$ All the above information can be found in a variety of sources, including very recent ones [@Austin-Trefethen-trig-interp2017; @Austin-PhD-thesis-interpolation2016], [@Trefethen-Weideman-exp-convergent-trapezoidal2014 Section 7].

Combining , , and , one has $$\label{eqn:interp-error-algebraic-convergence-vs-N}
\| g - \tilde{T}_N(g) \| \, \leq \, \left(2 \pi^{-1} \log N + \frac{8}{3} \right) 
\left(1 + \frac{\pi^2}{2} \right)^{l+1} \frac{M_{l+1}}  {n^{l+1}}$$ This indicates fast uniform algebraic convergence of the interpolant with respect to the number of knots. Moreover, under additional assumptions of analyticity of $g(\theta)$ in a strip of the complex plane Re$\, \theta \in (0, 2\pi)$, $|\mathrm{Im} \, \theta | < \delta$, convergence becomes exponential [@Trefethen-Weideman-exp-convergent-trapezoidal2014 (7.19)]: $$\label{eqn:interp-error-exp-convergence-vs-N}
\| g - \tilde{T}_N(g) \| \, \leq \, \frac{4M \exp[-\delta(N+1)/2]}
{1 - \exp(-\delta)}$$ We are also interested in the approximation of the integral $$\label{eqn:integral-of-g}
I \,=\, \int_0^{2\pi} g(\theta) \, d \theta$$ using the values of $g$ at the equispaced knots: $$\label{eqn:num-quadrature-of-g}
I_N \,=\, \frac{2\pi}{N} \, \sum_{m=1}^{N-1} g(\theta_m),
~~~ \theta_m = \frac{2\pi m}{N},$$ which is the trapezoidal rule for the numerical quadrature. Under the same analyticity assumptions as above, the error of this quadrature can be bounded as [@Trefethen-Weideman-exp-convergent-trapezoidal2014 (7.20)] $$\label{eqn:quadrature-error-exp-convergence-vs-N}
\| I_N - I \| \, \leq \, \frac{8\pi M \exp[-\delta(N+1)/2]}
{1 - \exp(-\delta)}$$ A similar result can be found in [@Austin-Trefethen-trig-interp2017 Theorem 1]. Adapted to our needs and notation, it states:

0.1in If $f$ is $l$ times continuously differentiable and $f^{(l)}$ is Lipschitz continuous, then $$\label{eqn:Austin-Trefethen-trig-interp-accuracy-power}
| I - \tilde I_N^{} | ,  \| f - \tilde t_N^{} \| 
= \mathcal{O}(N^{-(l+1)}).$$ If $f$ can be analytically continued to a $2\pi$-periodic function for $-\delta < \mathrm{Im} \, x < \delta$ for some $\delta > 0$, then for any $\hat \delta < \delta$, $$\label{eqn:Austin-Trefethen-trig-interp-accuracy-exp}
| I - \tilde I_N^{} | , \| f - \tilde t_N^{} \| = 
\mathcal{O}(\exp(-\hat \delta N)).$$ 0.1in

The qualitative conclusion of this section is that *trigonometric interpolation of a smooth periodic function provides a very accurate approximation of the function and its integrals.*

Preliminaries: Finite Difference Trefftz Schemes {#sec:FD-Trefftz}
================================================

Another preliminary subject, which will be needed in [Section \[sec:FLAME-slab\]]{}, is FLAME [@Tsukerman-JCP10; @Tsukerman05; @Tsukerman06; @Tsukerman-PBG08; @Tsukerman-book07; @Pinheiro07]. Recall that classical FD schemes are typically derived from Taylor expansions; but this is problematic if the solution is not sufficiently smooth – e.g. at material interfaces. That is the root cause of the notorious “staircase” effect at slanted or curved interface boundaries that do not conform geometrically to the grid lines. FLAME replaces Taylor polynomials with Trefftz functions, which often produces high-order schemes.

The key ideas of FLAME are as follows. Let a boundary value problem be defined in a computational domain $\Omega$ and consider a small subdomain $\Omega_h$ within which a difference scheme is to be formed. In $\Omega_h$, introduce a set of $m$ degrees of freedom (DoF). These DoF are, by definition, linear functionals, $l_\beta(u)$ ($\beta = 1,2,\dots, m$), each mapping any admissible field $u$ to a number (real or complex, depending on the problem). The simplest example of DoF for a scalar field $u$ is as set of nodal values $l_\beta(u) \equiv u(\mathbf{r}_\beta)$, where $\mathbf{r}_1,\dots,\mathbf{r}_m$ are a set of grid nodes in $\Omega_h$. In the case of vector fields, one may also consider fluxes, circulations, etc. as other examples of DoF .

Locally, within $\Omega_h$, the solution $u$ is approximated by a linear combination of Trefftz functions $\psi_\alpha$ ($\alpha = 1,2,\dots, n$) : $$\label{eqn:uh-eq-c-psi}
   u(\mathbf{r}) \,\approx\, u_h(\mathbf{r}) \equiv 
   \sum\nolimits_{\alpha=1}^n c_\alpha \psi_\alpha (\mathbf{r}) \,=\,
   \underline{c}^T \underline{\psi}(\mathbf{r}),$$ where $\underline{c} \in \mathbb{C}^n$ is a coefficient vector and $\underline{\psi}$ is a vector of basis functions (both generally complex). In $\Omega_h$, we seek an FD equation of the form $$\label{eqn:s-beta-l-beta-eq-0}
    \sum\nolimits_{\beta=1}^m s_\beta l_\beta (u) = 0,$$ where $\underline{s} = (s_1, s_2,\dots, s_m)^T$ is a vector of complex coefficients (a “scheme”) to be determined. In the simplest version of FLAME, the scheme is required to be exact for any linear combination of basis functions. Then, after straightforward algebra, one obtains [@Tsukerman05; @Tsukerman06] $$\label{eqn:s-in-null-Nt}
\underline{s} \in \mathrm{Null}(N^T), ~~~
\mathrm{where}~~  N^T_{\alpha \beta} = l_\beta (\psi_\alpha).$$ There are also least-squares versions of this idea [@Boag94; @Tsukerman-JCP10].

Many illustrative examples are given in [@Tsukerman-book07; @Tsukerman05; @Tsukerman06]. Here we mention just one of them, closely related to the construction of FLAME schemes in [Section \[sec:FLAME-slab\]]{}. For the 2D Helmholtz equation, one may consider a Trefftz basis set of eight plane waves traveling at the angles $\phi_0 + m \pi/4$ ($m = 0,1, \ldots, 7$), where $\phi_0$ is a given angle; practical choices are $\phi_0 = 0$ or $\phi_0 = \pi / 8$. Evaluating these plane waves over a standard $3 \times 3$ grid “molecule,” one obtains an $8 \times 9$ matrix $N^T$ whose null vector is the FLAME scheme. The result for $\phi_0 = 0$ is a nine-point ($3 \times 3$) order-six scheme [@Tsukerman06]. For $\phi_0 = \pi/8$, one arrives at a scheme derived by Babuška *et al*. in 1995 [@Babuska-Ihlenburg95] from very different considerations.

Trefftz Homogenization of Electromagnetic Structures {#sec:Trefftz-homogenization}
====================================================

We consider Trefftz-based homogenization of electromagnetic periodic structures (photonic crystals and metamaterials). The general description of the problem in this section follows [@Tsukerman-PLA17; @Tsukerman-Markel14] closely; but our focus here is on *Trefftz approximation*, the importance of other aspects of the problem notwithstanding.

The physical essence of the problem is as follows. A sample of a periodic material is illuminated by incoming monochromatic electromagnetic waves at a given frequency $\omega$ and the corresponding free-space wavenumber $k_0 = \omega / c$. To sidestep the complicated problem of field behavior at corners, the sample is assumed to be a finite-thickness slab contained between the planes $z = 0$ and $z = L$, and infinite in the $x$ and $y$ directions. The periodic medium in the sample is to be replaced with a homogeneous material in such a way that the scattering wave pattern would be preserved as accurately as possible.

Following [@Tsukerman-PLA17; @Tsukerman-Markel14], let us define the problem more precisely. Assume that the intrinsic dielectric permittivity $\tilde{\epsilon}({\mathbf{r}})$ within the slab is lattice-periodic, and that all material constituents are nonmagnetic, $\tilde{\mu}({\bf r}) = 1$. Let all constitutive relationships be local and linear, and let the sample be illuminated by monochromatic waves with a given far-field pattern; these waves are reflected by the metamaterial.

The problem has two principal scales (levels). *Fine-level* fields are the exact solutions of Maxwell’s equations for given illumination conditions for a given sample. These fields are denoted with small letters ${\mathbf{e}}$, ${\mathbf{d}}$, ${\mathbf{h}}$ and ${\mathbf{b}}$. In general, their variation in space is rapid and consistent with the microstructure of metamaterial cells. *Coarse-level* fields ${\mathbf{E}}$, ${\mathbf{D}}$, ${\mathbf{H}}$, ${\mathbf{B}}$ vary on a characteristic scale greater that the cell size. They represent some smoothed (averaged) versions of the fine-level fields and are auxiliary mathematical constructions rather than measurable physical quantities. The coarse-level fields are sought to satisfy Maxwell’s equations *and all interface boundary conditions* as accurately as possible.

Importantly, effective magnetic properties of metamaterials cannot be determined from the bulk behavior alone *as a matter of principle*. This is due, in particular, to the fact that the Maxwell equation $\nabla \times {\mathbf{H}}= -i k_0 {\mathbf{D}}$ is invariant with respect to an arbitrary simultaneous rescaling of vectors **H** and **D**. Loosely speaking, bulk behavior defines the dispersion relation only, while magnetic characteristics depend on the boundary impedance as well.

The fine-level fields satisfy macroscopic Maxwell’s equations of the form $$\label{eq:Maxwell_exact}
   \nabla \times {\mathbf{h}}({\mathbf{r}}) = -i k_0 \tilde{\varepsilon}({\mathbf{r}}) \, {\mathbf{e}}({\mathbf{r}}) \ , \ \ ~~
   \nabla \times {\mathbf{e}}({\mathbf{r}}) = i k_0 \, {\mathbf{h}}({\mathbf{r}})$$ everywhere in space, supplemented by the usual radiation boundary conditions at infinity. Outside the slab, the most general solution of can be written as a superposition of incident, transmitted and reflected waves. For the electric field, we can write these in the form of angular-spectrum expansions [@Tsukerman-Markel14]:

\[eq:Efine\] $$\begin{aligned}
\label{eq:Ei}
&& {\mathbf{e}}_i({\mathbf{r}}) = \int {\mathbf{s}}_i(k_x, k_y) e^{i\left(k_x x + k_y y +
    k_z z \right) } dk_x dk_y \ ,    \\
\label{eq:Et}
&& {\mathbf{e}}_t({\mathbf{r}}) = \int {\mathbf{s}}_t(k_x, k_y)
e^{i\left(k_x x + k_y y + k_z z \right)} d k_x d k_y \ , \ \ z > L \ ,    \\
\label{eq:Er}
&& {\mathbf{e}}_r({\mathbf{r}}) = \int {\mathbf{s}}_r(k_x, k_y) e^{i \left(k_x x + k_y y - k_z
    z \right)} d k_x d k_y \ , \ \ z < 0 \ ,\end{aligned}$$

where $$\label{eq:disp_vac}
   k_z = \sqrt{k_0^2 - k_x^2 - k_y^2} \ ,$$ and the square root branch is defined by the condition $0 \leq \arg k_z < \pi$. Expressions for the magnetic field are obtained from by using the second Maxwell equation in . In , ${\mathbf{s}}_i(k_x,k_y)$, ${\mathbf{s}}_t(k_x,k_y)$ and ${\mathbf{s}}_r(k_x,k_y)$ are the angular spectra of the incident, transmitted and reflected fields. Waves included in these expansions can be evanescent or propagating. For propagating waves, $k_x^2 + k_y^2 < k_0^2$, otherwise the waves are evanescent.

Everywhere in space, the total electric field ${\mathbf{e}}({\mathbf{r}})$ can be written as a superposition of the incident and scattered fields, viz, $$\label{eqn:ei_es}
{\mathbf{e}}({\mathbf{r}}) = {\mathbf{e}}_i({\mathbf{r}}) + {\mathbf{e}}_s({\mathbf{r}}) \ .$$ Outside the material, the reflected and transmitted fields form the scattered field: $$\label{eqn:es-eq-er-et}
{\mathbf{e}}_s({\mathbf{r}}) = \left\{
\begin{array}{ll}
{\mathbf{e}}_r({\mathbf{r}}) \ , & z < 0 \ , \\
{\mathbf{e}}_t({\mathbf{r}}) \ , & z > L \ .
\end{array} \right.$$ The scattered field *inside* the material is also formally defined by .

It is natural to approximate fine-level fields via a basis set of Bloch waves traveling in different directions: $$\label{eqn:eh-Bloch}
{\mathbf{e}}_{m \alpha}({\mathbf{r}}) = \tilde{{\mathbf{e}}}_{m \alpha}({\mathbf{r}}) \exp(i {\mathbf{q}}_{m \alpha} \cdot {\mathbf{r}}) \ , \ \
{\mathbf{h}}_{m \alpha} = \tilde{{\mathbf{h}}}_{m \alpha}({\mathbf{r}}) \exp(i {\mathbf{q}}_{m \alpha}
\cdot {\mathbf{r}}) \ ,$$ where index $\alpha$ labels both the wave vector and the polarization state of the Bloch wave in a lattice cell $m$; $\tilde{{\mathbf{e}}}_{m \alpha}({\bf r})$, $\tilde{{\mathbf{h}}}_{m \alpha}({\bf r})$ are the respective lattice-periodic factors. As the notation indicates, the basis is defined cell-wise; different bases in different lattice cells could be used. This makes the homogenization problem tractable and reducible to a single cell, rather than global and encompassing the whole sample.

On the coarse scale, a natural counterpart of the fine-scale Bloch basis is a set of generalized plane waves $$\label{eqn:Psi-eq-E0-H0-exp}
\Psi_{m \alpha} = \{ \mathbf{E}_{m \alpha}, \mathbf{H}_{m \alpha} \} = 
\{ \mathbf{E}_{0 m \alpha}, \mathbf{H}_{0m \alpha} \} \exp(i {\mathbf{q}}_{m \alpha} \cdot {\mathbf{r}})$$ which satisfy Maxwell’s equations in a homogeneous but possibly anisotropic medium; subscript ‘0’ indicates the field amplitudes to be determined. Further technical details of the procedure can be found in [@Tsukerman-Markel14; @Tsukerman-PLA17]. The final result is as follows. First, the coarse-level wave vector for each plane wave is taken to be the same as its counterpart for the corresponding Bloch wave, which is already reflected in our notation above , . Secondly, the amplitudes $\{ \mathbf{E}_{0 m \alpha}, \mathbf{H}_{0m \alpha} \}$ of each plane wave are the *boundary* average of the *tangential* components of the respective fine-scale Bloch wave: $$\label{eqn:EH0-eq-face-avrg-eh}
{\mathbf{E}}_{0m \alpha} \,=\, {\mathcal{A}}_{m} ^{\, \tau} \tilde{{\mathbf{e}}}_{m \alpha}, ~~~
{\mathbf{H}}_{0m \alpha} \,=\, {\mathcal{A}}_{m}^{\, \tau} \tilde{{\mathbf{h}}}_{m \alpha}$$ The averaging operator ${\mathcal{A}}_{m}^{\, \tau}$ for tangential components of a generic vector field $\mathbf{f}$ is defined, in the case of an orthorhombic cell $\mathbb{C}_{m}$, as $$\label{eqn:bdry-avrging-operator-tau}
({\mathcal{A}}_m^{\, \tau})_{\gamma} \, \mathbf{f} \,\equiv\, 
\frac{
\int_{\partial \mathbb{C}_{m}} f_{\gamma} \, |{\hat{\bf n}}\times \hat{{\mathbf{r}}}_{\gamma}| \, dS
}
{
\int_{\partial \mathbb{C}_{m}} |{\hat{\bf n}}\times \hat{{\mathbf{r}}}_{\gamma}| \, dS
}
, ~~~ \gamma = 1,2,3;  ~~ \hat{{\mathbf{r}}}_{1,2,3} = \hat{x}, \hat{y}, \hat{z}$$ Here $|{\hat{\bf n}}\times \hat{{\mathbf{r}}}_{\gamma}|$ acts simply as the Kronecker delta for the faces of the cell parallel to a given coordinate direction $\hat{{\mathbf{r}}}_{\gamma}$, $\gamma = 1,2,3$. Note that the averages in involve the *periodic factor* of the Bloch wave. The amplitudes ${\mathbf{E}}_{0m \alpha}$, ${\mathbf{H}}_{0m \alpha}$, along with the Bloch wave vector, define the coarse-level basis function $\alpha$ in a lattice cell $m$.

The homogenization procedure of [@Tsukerman-PLA17; @Tsukerman-Markel14] leads to a system of algebraic equations of the form $$\label{eqn:Psi-DB-eq-M-Psi-EH}
\Psi_{DB} \overset{l.s.}{=} {\mathcal{M}}\Psi_{EH}$$ Here ‘l.s.’ stands for ‘least squares’. Each column of the rectangular matrix $\Psi_{EH}$ corresponds to a given coarse-level basis function $\alpha$, and the entries of that column are the $xyz$-components of the wave amplitudes ${\mathbf{E}}_{0m \alpha}$, ${\mathbf{H}}_{0m \alpha}$. The number of columns $n$ is equal to the chosen number of basis functions; the number of rows is, in general, six, unless some of the field components are known to be zero (e.g. for $s$- or $p$-polarized waves). The $\Psi_{DB}$ matrix is completely analogous and contains the $\mathbf{DB}$ amplitudes derived from Maxwell’s curl equations: $$\label{B0-D0-eq-q-times-EH}
{\mathbf{B}}_{0m \alpha} = k_0^{-1} {\mathbf{q}}_{m \alpha} \times {\mathbf{E}}_{0m \alpha}, ~~
{\mathbf{D}}_{0m \alpha} = -k_0^{-1} {\mathbf{q}}_{m \alpha} \times {\mathbf{H}}_{0m \alpha}$$ The (local) material tensor is represented, in general, by a $6 \times 6$ matrix. Since the number of columns in matrix $\Psi_{EH}$ is typically greater than the number of rows, the matrix equation for the material tensor is solved in the least squares sense: $$\label{eqn:M-eq-DB-over-EH}
{\mathcal{M}}\,=\, \Psi_{DB} \Psi_{EH}^+;    
~~~~~ \delta_{\mathrm{l.s.}} \,=\, \| \Psi_{DB} - {\mathcal{M}}\Psi_{EH} \|_2$$ where $\Psi_{EH}^+$ is the Moore-Penrose pseudoinverse of $\Psi_{EH}$, and $\delta_{\mathrm{l.s.}}$ is the associated least-squares error.

![Example A of a layered medium from [@Markel-Tsukerman-PRB2013; @Tsukerman-Markel14]. The real part of $R$ (left) and $T$ (right) vs. the sine of the angle of incidence; non-asymptotic and nonlocal homogenization. The lattice cell contains three layers of widths $a/4$, $a/2$ and $a/4$, with scalar permittivities $\epsilon_1$, $\epsilon_2$, and $\epsilon_1$, respectively. ($\epsilon_1 = 4 + 0.1i$ and $\epsilon_2 = 1$.) Fine-level basis: $2n_{\mathrm{dir}}$ Bloch modes traveling at $n_{\mathrm{dir}} = 7$ different angles in $(-\pi/2, \pi/2)$; $n_{\mathrm{dir}} = 7.$ The kernel width parameter $\tau_0 = a$. The reflection and transmission coefficients from nonlocal homogenization are visually indistinguishable from the exact ones. The nonlocal procedure includes two additional DoF: the convolution integrals of the tangential components of the electric and magnetic fields.[]{data-label="fig:Re-RT-vs-angle-ExA-PRB2013"}](Re_R_vs_angle_Ex_A_PRB2013 "fig:"){width="0.85\linewidth"} ![Example A of a layered medium from [@Markel-Tsukerman-PRB2013; @Tsukerman-Markel14]. The real part of $R$ (left) and $T$ (right) vs. the sine of the angle of incidence; non-asymptotic and nonlocal homogenization. The lattice cell contains three layers of widths $a/4$, $a/2$ and $a/4$, with scalar permittivities $\epsilon_1$, $\epsilon_2$, and $\epsilon_1$, respectively. ($\epsilon_1 = 4 + 0.1i$ and $\epsilon_2 = 1$.) Fine-level basis: $2n_{\mathrm{dir}}$ Bloch modes traveling at $n_{\mathrm{dir}} = 7$ different angles in $(-\pi/2, \pi/2)$; $n_{\mathrm{dir}} = 7.$ The kernel width parameter $\tau_0 = a$. The reflection and transmission coefficients from nonlocal homogenization are visually indistinguishable from the exact ones. The nonlocal procedure includes two additional DoF: the convolution integrals of the tangential components of the electric and magnetic fields.[]{data-label="fig:Re-RT-vs-angle-ExA-PRB2013"}](Re_T_vs_angle_Ex_A_PRB2013 "fig:"){width="0.85\linewidth"}

![Example A of a layered medium from [@Markel-Tsukerman-PRB2013; @Tsukerman-Markel14]. Absolute error in $R$ (left) and $T$ (right) vs. $a / \lambda$; non-asymptotic and nonlocal homogenization. The lattice cell contains three layers of widths $a/4$, $a/2$ and $a/4$, with scalar permittivities $\epsilon_1$, $\epsilon_2$, and $\epsilon_1$, respectively. ($\epsilon_1 = 4 + 0.1i$ and $\epsilon_2 = 1$.) Fine-level basis: $2n_{\mathrm{dir}}$ Bloch modes traveling at $n_{\mathrm{dir}} = 7$ different angles in $(-\pi/2, \pi/2)$; $n_{\mathrm{dir}} = 7.$ The kernel width parameter $\tau_0 = a$. The accuracy of the nonlocal procedure is, by far, the highest. The nonlocal procedure includes two additional DoF: the convolution integrals of the tangential components of the electric and magnetic fields. []{data-label="fig:error-R-vs-angle-ExA-PRB2013"}](error_R_vs_lambda_Ex_A_PRB2013 "fig:"){width="0.85\linewidth"} ![Example A of a layered medium from [@Markel-Tsukerman-PRB2013; @Tsukerman-Markel14]. Absolute error in $R$ (left) and $T$ (right) vs. $a / \lambda$; non-asymptotic and nonlocal homogenization. The lattice cell contains three layers of widths $a/4$, $a/2$ and $a/4$, with scalar permittivities $\epsilon_1$, $\epsilon_2$, and $\epsilon_1$, respectively. ($\epsilon_1 = 4 + 0.1i$ and $\epsilon_2 = 1$.) Fine-level basis: $2n_{\mathrm{dir}}$ Bloch modes traveling at $n_{\mathrm{dir}} = 7$ different angles in $(-\pi/2, \pi/2)$; $n_{\mathrm{dir}} = 7.$ The kernel width parameter $\tau_0 = a$. The accuracy of the nonlocal procedure is, by far, the highest. The nonlocal procedure includes two additional DoF: the convolution integrals of the tangential components of the electric and magnetic fields. []{data-label="fig:error-R-vs-angle-ExA-PRB2013"}](error_T_vs_lambda_Ex_A_PRB2013 "fig:"){width="0.85\linewidth"}

As demonstrated in [@Tsukerman-PLA17], the homogenization accuracy can be further improved by including, in addition to the $EH$ amplitudes, integral DoF of the form $$\label{eqn:D-eq-Eps-star-E}
{\mathbf{D}}({\mathbf{r}}) \,=\, \int_{{\mathrm{\Omega}}} \mathcal{E}({\mathbf{r}}, {\mathbf{r}}') \, {\mathbf{E}}({\mathbf{r}}') \, d {\mathrm{\Omega}}$$ where $\mathcal{E}$ is a convolution kernel depending only on the coordinates tangential to the boundary of the sample: $$\mathcal{E}({\mathbf{r}}, {\mathbf{r}}') = 
\mathcal{E}({\hat{\bf n}}\times {\mathbf{r}}, {\hat{\bf n}}\times {\mathbf{r}}')$$ A natural (but certainly not unique) choice for this kernel is a Gaussian $$\mathcal{E}({\mathbf{r}}, {\mathbf{r}}') = 
\mathcal{E}_0 \exp(-\tau_0^{-2} \, |{\hat{\bf n}}\times ({\mathbf{r}}- {\mathbf{r}}')|^2 )$$ where the amplitude $\mathcal{E}_0$ and width $\tau_0$ are adjustable parameters, and ${\hat{\bf n}}$ is the unit normal vector.

Since our focus is on the approximation properties of Trefftz functions and not on the homogenization procedure per se, we do not discuss the physics of the problem here, or the merits and demerits of nonlocal vs. local theory.[^1] We also omit further technical details and limit ourselves to just one illustration example.

Shown in Figs. \[fig:Re-RT-vs-angle-ExA-PRB2013\] and \[fig:error-R-vs-angle-ExA-PRB2013\] are the reflection $R$ and transmission $T$ coefficients for electromagnetic waves propagating through a layered slab. These coefficients are defined in a standard way, as the ratio of the complex amplitudes of the reflected/transmitted waves to that of the incident wave. The geometric and physical parameters correspond to Example A of [@Markel-Tsukerman-PRB2013]: the lattice cell of a width $a$ contains three layers of widths $a/4$, $a/2$ and $a/4$, with scalar permittivities $\epsilon_1$, $\epsilon_2$, and $\epsilon_1$, respectively; $\epsilon_1 = 4 + 0.1i$ and $\epsilon_2 = 1$. The fine-level Trefftz basis contains $2n_{\mathrm{dir}}$ Bloch modes traveling at $n_{\mathrm{dir}} = 7$ equispaced angles in $(-\pi/2, \pi/2)$; $n_{\mathrm{dir}} = 7$. In nonlocal homogenization, the additional DoF are the integrals of the form , with the Gaussian kernel of width $\tau_0 = a$.

Fig. \[fig:Re-RT-vs-angle-ExA-PRB2013\] shows the real part of $R$ and $T$ as a function of the angle of incidence, for $a / \lambda = 0.2$. (The imaginary parts are not plotted to save space but are qualitatively similar). Since analytical solutions for wave propagation in layered media are fairly simple and well known, one may easily calculate the errors in $R$ and $T$; those are plotted in Fig. \[fig:error-R-vs-angle-ExA-PRB2013\].

The figures show that our numerical results, especially for nonlocal homogenization, are highly accurate. In fact, we are not aware of any alternative methods that could produce a comparable level of accuracy at a comparable computational cost.[^2]

What explains this high accuracy? Plausible mechanisms are presented in [Section \[sec:Trefftz-approximation\]]{}.

Electromagnetic Waves in Slab Geometries {#sec:FLAME-slab}
========================================

Formulation of the problem {#sec:Formulation}
--------------------------

![Schematic for the structure used in the sample FLAME-slab calculation. The structure consists of 10 dielectric pillars positioned aperiodically on a dielectric slab. Light is incident from the top, with wavenumber $k$ and incidence angle $\theta_{\mathrm{inc}}$.[]{data-label="fig:schematic"}](schematic_CAMWA.pdf){width="11cm"}

The general description of the problem in this section closely follows the recently published paper [@Mansha-OpEx17], which explores a new computational method, “FLAME-slab,” for electromagnetic wave scattering problems in aperiodic photonic structures – specifically, structures possessing short-range regularity but lacking long-range order, such as amorphous or quasicrystalline lattices. Structures of this type can exhibit a variety of interesting properties, e.g. highly isotropic band gaps and fractal photonic spectra, but are difficult to study numerically . FLAME-slab exploits the short-range regularity of the structure by generating a Trefftz basis in a relatively small segment of the structure.

As an example, we consider a slab substrate patterned with aperiodically placed but geometrically identical pillars (Fig. \[fig:schematic\]). The slab has thickness $d$, and there are 10 pillars of height $d$ and width $w = 0.8 \, d$. Both the substrate and the pillars have dielectric constant $\varepsilon=12$. The surrounding medium is air. In our calculations, we adopt computational units where the vacuum constants and the speed of light are all set to unity: $\varepsilon_{0} = 1$, $\mu_{0} = 1$, $c = 1$. Then the frequency $f$ has the units of $1/\lambda$, where $\lambda$ is the free space wavelength.

Light is incident from the top, as shown in Fig. \[fig:schematic\], with a wavenumber $k$ and incidence angle $\theta_{\mathrm{inc}}$ relative to the $z$-axis. We take the entire structure to be a supercell of length $L_{x}$, with quasi-periodic boundary conditions (see below). The electric and magnetic fields in the structure are governed by Maxwell’s equations: $$\begin{aligned}
\begin{aligned}
\nabla \times \mathbf{E} &= ik\mathbf{H},\\
\nabla \times \mathbf{H} &= -ik\varepsilon\mathbf{E}.
\end{aligned}
\label{eq:Maxwell}\end{aligned}$$ We consider the case where the electric field is $s$-polarized, $\mathbf{E}=E\hat{y}$, so that the magnetic field has the form $\mathbf{H}=H_{x}\hat{x}+H_{z}\hat{z}$. The quasi-periodic boundary conditions are: $$\begin{aligned}
\begin{aligned}
\mathbf{E}(L_{x}/2,z) &= \mathbf{E}(-L_{x}/2,z)\exp(ik_{x}L_{x}), \\
\mathbf{H}(L_{x}/2,z) &= \mathbf{H}(-L_{x}/2,z)\exp(ik_{x}L_{x}),
\end{aligned}\end{aligned}$$ where $k_x = k\sin \theta_{\mathrm{inc}}$ is the $x$-component of the incident wave vector ${\mathbf{k}}$.

The scattered electric field is defined as $$E_s(\mathbf{r}) = E_{\mathrm{tot}}(\mathbf{r}) - E_{\mathrm{inc}}(\mathbf{r}), \hspace{1cm} [\,\mathbf{r}\equiv (x,z)\,],$$ where $E_{\mathrm{tot}}$ and $E_{\mathrm{inc}}$ are the total and incident electric fields, respectively. The magnetic field is split similarly. The scattered field is purely outgoing on both the upper side (towards the negative $z$-direction) and the lower side (towards the positive $z$-direction) of the structure.

Fig. \[fig:discretization\](a) shows the discretization scheme for FLAME. The structure is discretized into $N_x$ grid points in the horizontal direction. In the vertical direction, the number of layers is deliberately limited to three ($z_{-}$, $z_{0}$, $z_{+}$), to demonstrate that FLAME-slab can work well on very coarse grids. The electric fields in these three layers are denoted with $E^{\alpha}_m$, $\alpha=\{-,0,+\}$, $m = 1,2, \ldots, N_x$. Similarly, the magnetic fields in the upper and bottom layers are denoted with $H^{\beta}_m$, $\beta=\{-,+\}$.

We define three distinct types of patches, with their corresponding grid “molecules” and FD stencils. The first is a standard 9-point stencil containing just the electric field degrees of freedom (DoF), as shown in the left panel of Fig. \[fig:discretization\](b). The second is a 6-point stencil over the middle and top layers, containing both the electric and magnetic fields (middle panel). The third is a 6-point stencil over the middle and bottom layers, containing both the electric and magnetic fields (right panel of Fig. \[fig:discretization\](b)).

Each type of patch thus contains 9 degrees of freedom. FLAME uses 8 basis functions, to be determined by solving Maxwell’s equations for “Trefftz cells” matching the local dielectric environment in each patch. Each Trefftz cell contains a segment of length $L_{i}$ with a single pillar on the substrate; quasi-periodic boundary conditions are imposed. We choose $L_{i} \ll L_x$, so that Maxwell’s equations can be solved much more rapidly for the Trefftz cell than for the entire aperiodic structure. We generate 8 different Trefftz basis functions by picking two different segment lengths ($L_1$ and $L_2$), and four different angles of incidence for each $L_i$. To compute the fields in the Trefftz cell, we use the existing rigorous coupled wave analysis (RCWA) solver $S^{4}$ [@Liu2012].

![(a) Discretization of the structure into $N_x \times 3$ nodes ($N_x$ in the horizontal direction and 3 layers in the vertical direction). (b) Variations of the 9-point stencils. Left: 9 nodes with a single DoF (the values of the electric field). Center and right: three nodes with a single DoF (the electric field) and another three nodes (double circles) with double DoF (electric and magnetic fields.)[]{data-label="fig:discretization"}](discretization_CAMWA.pdf){width="10cm"}

The FLAME procedure now yields a matrix equation of the form $$\mathbf{A}_{\mathrm{FL}} \, \boldsymbol{\psi}_{tot} = 0,
\label{eq:Aflame_psi}$$ where $\mathbf{A}_{\mathrm{FL}}$ is a matrix of stencil coefficients and $\boldsymbol{\psi}_\mathrm{tot}$ is a column vector containing the nodal values of the total electric and magnetic fields. In our examples, $\mathbf{A}_{\mathrm{FL}}$ has the size $3N_{x}\times 5N_{x}$, and $\boldsymbol{\psi}_\mathrm{tot}$ has the size $5N_{x}\times 1$; we emphasize that this is just one possible choice of discretization, and other choices can be handled in a completely analogous way. Details about the calculation of $\mathbf{A}_{\mathrm{FL}}$ can be found in [@Mansha-OpEx17].

FLAME schemes need to be supplemented with radiation boundary conditions. One way of implementing such conditions is via the Dirichlet-to-Neumann (DtN) maps in the semi-infinite air strips above and below the slab. DtN maps can be efficiently calculated via Fast Fourier Transforms (FFTs). More specifically, from Maxwell’s equations in free space, $$H_s(x,z) = \frac{i}{\omega} \frac{\partial E_s}{\partial z}
\label{eq:EandHMaxwell}$$ The operating frequency $\omega=2\pi f = k$, under the assumed normalization $c=1$. We expand the scattered electric field into its Fourier series: $$\label{eq:EscatBloch}
   E_s(x,z) = \sum_{n} c_{n} \exp[i(k_{nz}z + k_{nx}x)] \exp(iqx),$$ where the factor of $\exp(iqx)$ comes from the quasiperiodic boundary conditions in the $x$ direction, with $q = k\sin \theta_{\mathrm{inc}}$. The summation $n$ runs over the integer values, $k_{nx}=2\pi n/ L_x$ is the horizontal wavenumber, and $$k_{nz} = \pm \sqrt{k^2 - (k_{nx} + q)^2}.
\label{eq:knzexpression}$$ In the above equation, the choice of $\pm$ depends upon the layer we are dealing with ($-$ for the upper layer and $+$ for the bottom layer), so that the scattered field is outgoing. Eqs.  and give $$\label{eq:HscatBloch}
   H_s(x,z) = -\frac{1}{\omega}\sum_{n} c_{n} \, k_{nz} \, \exp\left[i(k_{nx}x + k_{nz}z)\right] \exp(iqx).$$ The coefficients $c_n$ in can be efficiently computed via a Fast Fourier Transform, and then the scattered magnetic field can be obtained via the respective inverse transform (detailed expressions can be found in [@Mansha-PhDthesis18]). This leads to equations in the following matrix form: $$\label{eq:Adef}
    \begin{pmatrix} 
       \mathbf{A}_{\mathrm{FL}} \\ \mathbf{A}_{\mathrm{BC}} 
    \end{pmatrix} 
    \boldsymbol{\psi}_s
    = 
    \begin{pmatrix} 
       -\mathbf{A}_{\mathrm{FL}} \boldsymbol{\psi}_{\mathrm{inc}} \\ \boldsymbol{0}
     \end{pmatrix},$$ where $\mathbf{A}_{\mathrm{FL}}$ is a sparse sub-matrix obtained using FLAME, and $\mathbf{A}_\mathrm{BC}$ is sub-matrix obtained from the boundary relations [@Mansha-OpEx17]. In our 2D examples, standard direct solvers in Matlab were sufficient for finding $\boldsymbol{\psi}_s$. In 3D, iterative solvers will need to be used, but this issue is completely beyond the scope of the present paper.

Results
-------

Fig. \[fig:EHfields\] compares the fields calculated using FLAME-slab to a reference RCWA calculation. The structure is the one shown in Fig. \[fig:schematic\], with frequency $f=0.25$ and incidence angle $\theta_{\mathrm{inc}} = 30\degree$. For the FLAME-slab calculation, we take a horizontal discretization of $N_x=101$, and precompute the Trefftz basis functions with $N_G=150$ (the number of expansion terms used in the RCWA subroutine [@Liu2012]) and $N_T=800$ (the cell dicretization used for storing the Trefftz basis functions). The pure RCWA reference solution is computed using $N_{G}^{ref}=1000$ – an “overkill” setting meant to produce a highly accurate solution. The figure shows two representative field components: the real part of the scattered electric field ($E^{0}_s$) in the middle layer ($z_0$) in Fig. \[fig:EHfields\](a), and the scattered magnetic field ($H^{+}_s$) in the bottom layer ($z_{+}$) in Fig. \[fig:EHfields\](b). The FLAME-slab solution is seen to be in excellent agreement with the RCWA solution.

![(a) Real part of the scattered electric field $E_{scat}$ in the middle layer ($z_0$). (b) Real part of the scattered magnetic field $H_{scat}$ in the bottom layer ($z_{+}$). The calculations were done for the slab shown in Fig. \[fig:schematic\], with $f=0.25$ and $\theta_{\mathrm{inc}}=30\degree$. The FLAME-slab parameters are $N_x=101, $$N_G=150$, and $N_T=800$. Blue dots show the FLAME-slab results and the red curve shows the result from RCWA obtained by setting $N_{G}^{ref}=1000$.[]{data-label="fig:EHfields"}](EH_field_plots_lambda_040_theta_30.pdf){width="90.00000%"}

![Consistency error ($\xi$) vs. angle of incidence $\theta_{\mathrm{inc}}$ for the 10 pillar system as shown Fig. \[fig:schematic\]. The value of the parameters used are: $N_x=101$, $N_G=150$, $N_T=800$ and $f=0.25$.[]{data-label="fig:cons_error_vs_theta"}](cons_error_vs_theta_freq025_10rods_Nx101_NT800_NG150.pdf){width="90.00000%"}

The central issue of this paper is approximation, and the finite-difference measure most closely related to it is the (normalized) consistency error $$\xi = \frac{\| \mathbf{A}_{\mathrm{FL}} \, 
   \boldsymbol{\psi}_{\mathrm{tot}}^{\mathrm{ref}} \|}
   {\| \mathbf{A}_{\mathrm{FL}} \| \, 
    \|\boldsymbol{\psi}_{\mathrm{tot}}^{\mathrm{ref}} \|}.
\label{eq:consistency}$$ where Euclidean vector norms and the Frobenius matrix norm are implied.

In , $\boldsymbol{\psi}_{\mathrm{tot}}^{\mathrm{ref}}$ should ideally be the exact solution, which is not available; hence an overkill RCWA solution with $N_{G}^{ref}=1000$ is used in its stead.

Since FLAME-slab contains a few adjustable parameters, we study the dependence of the consistency error on these parameters separately.

Fig. \[fig:cons\_error\_vs\_theta\] displays the consistency error versus the incidence angle $\theta_{\mathrm{inc}}$. For this calculation, we set $N_x=101$, $N_G=150$, $N_T=800$ and $f=0.25$. The consistency error oscillates but remains bounded by $\lesssim 10^{-5}$ over the entire range of $\theta_{\mathrm{inc}}$.

![Consistency error ($\xi$) vs. $\lambda_{\mathrm{vac}}$ (vacuum wavelength) for the slab structure shown in Fig. \[fig:schematic\], with $\theta_{\mathrm{inc}} = 30\degree$. The FLAME-slab parameters are $N_x=101$, $N_G=150$, and $N_T=800$.[]{data-label="fig:cons_error_vs_lambda"}](cons_error_vs_lambda_theta_30_freq025.pdf){width="10cm"}

![Consistency error ($\xi$) vs. $N_x$ for the slab shown in Fig. \[fig:schematic\], with $f = 0.25$ and $\theta_{\mathrm{inc}} = 0$. The FLAME-slab parameters are $N_G=150$ and $N_T=800$.[]{data-label="fig:cons_error_vs_Nx"}](cons_error_vs_Nx_NG150_cachenx_800.pdf){width="10cm"}

Fig. \[fig:cons\_error\_vs\_lambda\] shows the consistency error versus the vacuum wavelength $\lambda_{\mathrm{vac}}$ for the 10 pillar system, with fixed incidence angle $\theta_{\mathrm{inc}} = 30\degree$. The FLAME-slab parameters are fixed at $N_x=101$, $N_G = 150$, and $N_T = 800$. As $\lambda_{\mathrm{vac}}$ is increased, $\xi$ decreases from $10^{-4}$ to around $10^{-6}$. Past this point, $\xi$ saturates.

Fig. \[fig:cons\_error\_vs\_Nx\] shows the consistency error versus spatial discretization $N_x$, for $f = 0.25$ and normal incidence $\theta_{\mathrm{inc}}=0$. The other FLAME-slab parameters are $N_G=150$ and $N_T=800$. The consistency error decreases with $N_x$, saturating at $\approx 10^{-7}$ for $N_x \gtrsim 500$.

For the purposes of the paper, the main qualitative conclusion of this section is that *Trefftz functions*, on which FLAME-slab is based, *provide an accurate approximation of the electromagnetic field in a geometrically and physically complex structure*.

The Accuracy of Trefftz Approximations {#sec:Trefftz-approximation}
======================================

An Interpolation Argument {#sec:Interpolation-argument}
-------------------------

The numerical results for the two application examples of the previous sections show that Trefftz approximations are surprisingly effective. What explains their high accuracy?

As noted in [Section \[sec:Intro\]]{}, in the mathematical literature this question has been studied primarily for homogeneous subdomains (harmonic polynomials, plane/cylindrical/spherical wave expansions) but needs to be posed much more broadly, because complex inhomogeneous media are of great theoretical and practical interest. This section is an attempt to understand the general mechanisms of high accuracy of Trefftz approximations. Due to the complexity of this subject, some of the material, especially that of [Section \[sec:Random-matrices\]]{}, is speculative and intended to stimulate further analysis and discussion.

In the case of Trefftz homogenization (Section \[sec:Trefftz-homogenization\]), one can apply an interpolation argument using the summary in Section \[sec:Trig-interpolation\]. Indeed, the key parameters in our homogenization methodology are the boundary averages of the Bloch fields . Each of these averages is, trivially, a periodic function of the angle (direction) of propagation of the respective Bloch wave and, as such, can be accurately approximated by the trigonometric interpolant over a set of equispaced knots. But these knots correspond precisely to the basis set of Bloch waves chosen in our procedure. Per [Section \[sec:Trefftz-homogenization\]]{}, the accuracy of this interpolation is $\mathcal{O}(N^{l+1})$ if the respective Bloch average is $l$ times continuously differentiable, or, under additional analyticity assumptions, even $\mathcal{O}(\exp(-\alpha N))$, where $N$ is the size of the Bloch basis set (which is the same as the number of interpolation knots).

![An inhomogeneous scatterer ${\mathrm{\Omega}}$ (solid red) is enclosed in a shell $\tilde{{\mathrm{\Omega}}}$ (textured). The material parameters are fixed within $\Omega$ in all cases. However, in $\tilde{{\mathrm{\Omega}}}$ these parameters for the unknown field (dashed arrow) may differ from the parameters used in the construction of the Trefftz training set (solid arrows). A local Trefftz approximation of the unknown field is sought in a small subdomain $\Omega_h \subset \Omega$.[]{data-label="fig:test-waves-random"}](test_waves_random.png){width="0.65\linewidth"}

In our second example of wave propagation and scattering in a slab geometry, the interpolation argument is not sufficient. This is because our Trefftz functions are defined over *a segment* of the structure, whereas the full electromagnetic problem is defined over the *whole structure*. Hence a more sophisticated explanation for the accuracy of Trefftz approximations in this case is needed.

We start with a slightly more abstract physical setup than that of Fig. \[fig:test-waves\]. Namely, let us assume, as before, that an inhomogeneous scatterer occupies a Lipschitz domain ${\mathrm{\Omega}}$ (solid red in Fig. \[fig:test-waves-random\]) which is enclosed in a shell $\tilde{{\mathrm{\Omega}}}$ (textured area). As previously, we consider a Trefftz “training set” corresponding to several incident waves, and are interested in approximating a different, generally unknown, solution in a small subdomain ${\mathrm{\Omega}}_h \subset {\mathrm{\Omega}}$. This approximation can be used, for example, to generate a high-order difference scheme in ${\mathrm{\Omega}}_h$, as was done in [Section \[sec:FLAME-slab\]]{}.

The Trefftz training set is also generated for fixed position-dependent parameters in ${\mathrm{\Omega}}\cup \tilde{{\mathrm{\Omega}}}$. Importantly, however, the unknown solution may correspond to material parameters which *differ* in $\tilde{{\mathrm{\Omega}}}$ from those assumed for the training set (but are the same in ${\mathrm{\Omega}}\supset {\mathrm{\Omega}}_h$). The presence of the variable layer $\tilde{\Omega}$ makes this case peculiar. The following section examines why accurate *local* Trefftz approximations can still be expected.

An Auxiliary “Reference” Basis {#sec:Auxiliary-basis}
------------------------------

Let us assume that in $\Omega_h$ there is *an auxiliary basis* $\zeta_\alpha$ ($\alpha = 1,2, ..., n_\zeta$) which can provide an accurate approximation of a (generic) solution of the wave equation: $$\label{eqn:E-vs-modes}
     u(\mathbf{r}) = \sum_{\alpha}
     \gamma_\alpha \zeta_\alpha (\mathbf{r}) + \delta({\mathbf{r}}),
     ~~ {\mathbf{r}}\in \Omega_h$$ $$\| \underline{\gamma} \|_2 \equiv \| \{\gamma_\alpha \} \|_2
     \leq C(\Omega_h, n_\zeta, k) \| u \|_{H^1(\Omega_h)},
     ~~ \| \delta \|_{H^1(\Omega_h)} \leq c(\Omega_h, n_\zeta, k) \| u \|_{H^1(\Omega_h)}$$ Here $\delta$ is an error term, $\underline{\gamma}$ is a coefficient vector, $C$ and $c$ are some generic constants, the latter being “small” in some sense (see Theorems below). In the specific example of $s$-wave scattering in Section \[sec:FLAME-slab\], the unknown is the $E$-field; but here we use the “generic” symbol $u$ as an indication that our analysis could be applied more broadly.

Assuming that holds, one applies it to the training set $\underline{\psi}_T$ of $n_T$ Trefftz waves, and arrives at the linear transformation $$\underline{\psi}_T({\mathbf{r}}) = P_{\zeta \rightarrow T} \underline{\zeta}({\mathbf{r}}) 
   + \underline{\delta}_T({\mathbf{r}})$$ where column vectors are underlined; $P_{\zeta \rightarrow T}$ is the $ n_T \times n_\zeta $ transformation matrix, and $\underline{\delta}_T$ is the approximation error for the Trefftz functions in terms of the local $\zeta$ basis. If $n_\zeta \leq n_T$, and if matrix $P_{\zeta \rightarrow T}^* P_{\zeta \rightarrow T}$ is invertible, then $$\underline{\zeta}({\mathbf{r}}) = P_{\zeta \rightarrow T}^+ \underline{\psi}_T({\mathbf{r}}) + 
   P_{\zeta \rightarrow T}^+ \underline{\delta}_T({\mathbf{r}})$$ where the ‘+’ subscript indicates the Moore-Penrose pseudoinverse.

\[fig:Trefftz-vs-reference-basis\] ![A schematic illustration of the role of the reference basis. If the Trefftz basis and the solution of a given boundary value problem can be approximated via the reference basis, and if the reference-to-Trefftz transformation has a bounded pseudoinverse, then one can approximate the solution via the Trefftz basis (by following, conceptually, the two solid arrows in the sketch). ](Trefftz_vs_reference_basis "fig:"){width="0.4\linewidth"}

The solution in ${\mathrm{\Omega}}_h$ can therefore be expressed as $$\label{eqn:Eexact-Trefftz-approximation}
    u({\mathbf{r}}) = \underline{\gamma}^T \underline{\zeta}({\mathbf{r}}) = 
    \underline{\gamma}^T P_{\zeta \rightarrow T}^+ \underline{\psi}_T({\mathbf{r}}) + 
    \underline{\gamma}^T P_{\zeta \rightarrow T}^+ \underline{\delta}_T({\mathbf{r}})
   +  \delta_u({\mathbf{r}})$$ where $\delta_u$ is the approximation error of this solution via the $\zeta$ basis. Thus the smallness of the Trefftz approximation error hinges on the smallness of the norm of the pseudoinverse $P_{\zeta \rightarrow T}^+$ – that is, on the inverse of its minimum singular value $\sigma_{\min}$; we discuss that below.

The transformations above are schematically illustrated in Fig. \[fig:Trefftz-vs-reference-basis\]. If the Trefftz basis and the solution $u$ can be approximated via the reference basis as in , and if $\sigma_{\min} (P)$ is bounded from below, then one can approximate $u_{\mathrm{exact}}$ via the Trefftz basis (by following, conceptually, the two solid arrows in the sketch).

An example of this auxiliary basis is, in the special case of a homogeneous domain $\Omega_h$, a set of cylindrical harmonics $\zeta_{\mathrm{cyl}}(r,k,\theta,n)
= J_n (kr) \exp(in \theta)$, $n = 0, \pm 1, \pm 2, \ldots$; $J_n$ is the Bessel function of the first kind. Detailed error analyses have been carried out by Melenk, Hiptmair, Moiola and Perugia [@Melenk95; @Hiptmair2016; @Moiola-PhD11]. For our purposes, the most convenient final results can be found in [@Babuska97; @Melenk99].

**[@Babuska97 Theorem 4].** Let ${\mathrm{\Omega}}\subset \mathbb{R}^2$ be a simply connected, bounded Lipschitz domain. Let $\tilde{{\mathrm{\Omega}}} \supset {\mathrm{\Omega}}$ and assume that $u \in L^2(\tilde{{\mathrm{\Omega}}})$ solves the homogeneous Helmholtz equation on $\tilde{{\mathrm{\Omega}}}$. ¹ Then $$\label{eqn:Babuska-Melenk-cyl-harm-estimate1}
    \inf_{u_p \in V_p} \| u - u_p \|_{H^1({\mathrm{\Omega}})}  \,\leq\, C \exp(-\gamma p) 
    \, \| u \|_{L^2{(\tilde{{\mathrm{\Omega}}})}}$$ where $V_p \equiv \mathrm{span} \{\zeta_{\mathrm{cyl}}(r,k,\theta,n) \}$, $n = 0,1, \ldots, p$; $C, \gamma$ depend only on ${\mathrm{\Omega}}$, $\tilde{{\mathrm{\Omega}}}$, and the wavenumber $k$.

Under the assumptions of this theorem, the presence of a “buffer region” $\tilde{{\mathrm{\Omega}}} - {\mathrm{\Omega}}$ ensures that high-order harmonics from the boundary of $\tilde{{\mathrm{\Omega}}}$ die out sufficiently. If this assumption is not made, an alternative error estimate, dependent on the level of smoothness of the solution, reads:

**[@Babuska97 Theorem 5].** Let $\Omega \subset \mathbb{R}^2$ be a simply connected, bounded Lipschitz domain, star-shaped with respect to a ball. Let the exterior angle of ${\mathrm{\Omega}}$ be bounded from below by $\lambda \pi$, $0 < \lambda < 2$. Assume that $u \in H^s({\mathrm{\Omega}})$, $s > 1$, satisfies the homogeneous Helmholtz equation. Then[^3] $$\label{eqn:Babuska-Melenk-cyl-harm-estimate2}
   \inf_{u_p \in V_p} \| u - u_p \|_{H^j({\mathrm{\Omega}})}  \,\leq\, 
    C_j \, \left( \frac{\ln^2 p}{p} \right)^{\lambda(s-j)}
    \, \| u \|_{H^s{(\tilde{{\mathrm{\Omega}}})}},
    ~~~ j = 0,1, \ldots, [s]$$ Obviously, in our case ${\mathrm{\Omega}}_h$ plays the role of the generic ${\mathrm{\Omega}}$ in the estimates above. These estimates of the error term $\delta$ in are valid for the 2D Helmholtz equation in a physically *homogeneous medium* within ${\mathrm{\Omega}}_h$.

Also in the special case of a homogeneous domain $\Omega_h$, and the Trefftz set consisting of plane waves traveling in $n_T$ equispaced angular directions, the norm of the pseudoinverse $P_{\mathrm{cyl} \rightarrow \mathrm{PW}}^+$ can be evaluated explicitly. From the Jacobi-Anger expansion, the entries of the matrix $\hat{P} \equiv P_{\mathrm{cyl} \rightarrow \mathrm{PW}}$ are $$\label{eqn:P-cyl-to-PW}
  \hat{P}_{ml} \,=\, i^l \exp \left(-iml \, \frac{2\pi}{n_T} \right),
  ~~ 0 \leq m \leq n_T-1, ~~ 0 \leq l \leq n_\zeta - 1$$ This matrix corresponds to a discrete Fourier transform, and its columns are easily shown to be orthogonal, so that $$\label{eqn:P-P-diag}
    \hat{P}^* \hat{P} \,=\, n_T I_{n_\zeta},
    ~~~ n_\zeta \leq n_T$$ where $I_{n_\zeta}$ is the identity matrix of dimension $n_\zeta$. It then immediately follows that $$\label{eqn:sigma-min-P-cyl-PW}
   \| P_{\mathrm{cyl} \rightarrow \mathrm{PW}}^+ \|_2 \,=\,
   \sigma_{\min}^{-1}( P_{\mathrm{cyl} \rightarrow \mathrm{PW}}) = n_T^{-\frac12}$$ so in this case stability of the transformation is guaranteed.

A Connection with Random Matrix Theory {#sec:Random-matrices}
--------------------------------------

A natural, and critical, question is whether the well-posedness of the transformation noted above is accidental and valid in special cases only, or whether it has broader applicability. Practical experience with multiparticle problems, random and quasi-random structures of different kind [@Tsukerman05; @Tsukerman06; @Tsukerman-PBG08; @Dai-Webb11; @Mansha-OpEx17], [@Tsukerman-book07 Chapters 4, 6] strongly suggests the latter. Rigorous mathematical analysis is so far available only for a narrow subset of cases [@Melenk95 Chapter IV], [@Melenk99 Section 3], [@Laghrouche-jumps2005; @Imbert-Gerard-interp-PW2015], and may constitute an interesting direction of future research.

In the remainder of this section, we outline – on physical grounds – a curious connection between the accuracy of Trefftz approximations and the theory of random matrices. This theory dates back to von Neumann and Wigner [@vonNeumann-inverting-matrices-1947; @Wigner-random-matrices-1955] and is now quite mature [@Akemann-handbook-random-matrices-2011; @Tao2010; @Rudelson-Vershynin-least-sing-value2008; @Rudelson-Vershynin-smallest-sing-value2009; @Rudelson-Vershynin-nonasymptotic-sing-value2010]. Particularly relevant to us is the following result.

0.1in

**Rudelson & Vershynin [@Rudelson-Vershynin-nonasymptotic-sing-value2010 Theorem 3.3]**.\
Let $A$ be an $N \times n$ random matrix whose entries are independent and identically distributed (i.i.d.) subgaussian random variables with zero mean and unit variance. Then $$P \left( \sigma_{\min}(A) \le \epsilon (\sqrt{N} - \sqrt{n-1}) \right) 
    \le (C\epsilon)^{N-n+1} + c^N,
    ~~~ \epsilon \ge 0$$ where $C > 0$ and $c \in (0,1)$ depend only on the subgaussian moment of the entries.

0.1in

The connection of this theorem with the previous subsection can be outlined as follows.

-   The Trefftz “training set” can be viewed as a particular realization of some random distribution (e.g. angles of incidence randomly chosen and/or random properties of the “shell” $\tilde{\Omega}$). A notable feature of random matrix theory is *universality*: only mild dependence of the spectral bounds on the distribution of the random variables.

-   One major restrictive condition, however, is that the matrix entries be i.i.d. variables. Strictly speaking, this condition can be immediately ascertained only under additional symmetry assumptions, e.g. the bases being invariant under rotation by a given angle. It is hoped that such strong assumptions can be relaxed.

-   The assumption of zero mean is less restrictive and valid if the probability distribution of each function in the Trefftz training set $\underline{\psi}_T({\mathbf{r}})$, for all ${\mathbf{r}}$, is invariant with respect to the sign change of that function.

-   Clearly, the theorem is applied with $N \equiv n_T$, $n \equiv n_\zeta$.

-   The assumption that the distribution is subgaussian is satisfied, in particular, by all bounded random variables and hence is not restrictive.[^4][^5]

-   Complex bases and matrices can be decomplexified by the substitutions of the form $\psi \rightarrow (\mathrm{Re} \, \psi, \mathrm{Im} \, \psi)^T$, $P \rightarrow 
        \begin{pmatrix}
        \mathrm{Re} \, P & -\mathrm{Im} \, P \\
        \mathrm{Im} \, P & \mathrm{Re} \, P 
        \end{pmatrix}
        $. This preserves the relevant norms and hence does not affect the spectral bounds.

-   The assumption of unit variance is obviously a matter of scaling only.

The theorem affirms that stability of the transformation is not accidental. In fact, with a probability close to one, $\sigma_{\min}(P)$ is not small, for any reasonable choice of the Trefftz basis.

Conclusion {#sec:Conclusion}
==========

The key argument of this paper is that Trefftz approximations – that is, approximations by functions satisfying (locally) a given differential equation – deserve to be studied and applied more broadly than is traditionally done. Conventionally, these approximations are used in homogeneous subdomains, where the underlying differential equation has constant coefficients; this is done in various contexts (GFEM, DG, FD).

As an illustration of a much broader use of Trefftz functions, the paper reviews two disparate but representative examples: (i) non-asymptotic and nonlocal two-scale homogenization of periodic electromagnetic media, and (ii) special Trefftz FD (FLAME) schemes for wave scattering from photonic structures with slab geometries. In both cases, Trefftz approximations are applied in complex inhomogeneous domains and prove to be quite effective.

We discuss possible mechanisms engendering the high accuracy of Trefftz approximations. One such mechanism is trigonometric interpolation, which itself is known to be surprisingly accurate for smooth periodic functions, in comparison with other typical forms of interpolation. We also outline, on physical grounds, a curious connection of Trefftz approximations with the theory of random matrices.

It is hoped that these considerations will stimulate further mathematical research and practical applications of Trefftz-based methods.

Acknowledgment {#acknowledgment .unnumbered}
==============

The work of IT was supported in part by the US National Science Foundation Grants DMS-1216927 and DMS-1620112. The research of SM and YC was supported by the Singapore MOE Academic Research Fund Tier 2 Grant MOE2016-T2-1-128, the Singapore MOE Academic Research Fund Tier 2 Grant MOE2015-T2-2-008, and the Singapore MOE Academic Research Fund Tier 3 Grant MOE2016-T3-1-006. The work of VM was supported in part by the US National Science Foundation Grants DMS-1216970.

IT thanks Ralf Hiptmair, Andrea Moiola, Lise-Marie Imbert-Gérard and Ben Schweizer for discussions.

References {#references .unnumbered}
==========

[^1]: It should, however, be noted that our nonlocal procedure operates in real space, in contrast with $k$-space techniques that we critiqued elsewhere [@Markel-Tsukerman-PRB2013].

[^2]: The latter provision is needed to exclude from consideration “brute force” numerical optimization of the material tensor.

[^3]: There is an apparent misprint in [@Babuska97]: $H^k$ instead of $H^s$ in the norm on the right hand side.

[^4]: A random variable $X$ is called subgaussian if there exists a positive constant $w$ such that $P (|X| > x) \le 2 \exp(-x^2/w^2)$ for $x > 0$.

[^5]: There is an apparent misprint in [@Rudelson-Vershynin-nonasymptotic-sing-value2010]: $n \times n$ instead of $N \times n$.
