# A Quasi-Trefftz Method Based on Local Polynomial Impedance Boundary Conditions for Time-Harmonic Wave Propagation Problems

Matthias Rivet, Sébastien Pernet, Sébastien Tordeux

► To cite this version:

Matthias Rivet, Sébastien Pernet, Sébastien Tordeux. A Quasi-Trefftz Method Based on Local Polynomial Impedance Boundary Conditions for Time-Harmonic Wave Propagation Problems. 2026. hal-05452348

HAL Id: hal-05452348

https://inria.hal.science/hal-05452348v1

Preprint submitted on 10 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Quasi-Trefftz Method Based on Local Polynomial Impedance Boundary Conditions for Time-Harmonic Wave Propagation Problems

Matthias Rivet¹,², Sébastien Pernet¹, and Sébastien Tordeux²

¹DTIS, ONERA, Université de Toulouse, 31000, Toulouse, France

²EPI Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142

January 10, 2026

## Abstract

Facing the need to solve time-harmonic wave problems in wide propagation domains, one needs to resort to iterative algorithms, to which classic numerical methods are usually poorly suited. The Trefftz method, based on the use of actual local solutions as basis functions, can then appear as a viable option for this kind of configurations.

Thus, we first propose a general formalism allowing to gather different wave problems, so that a general Ultra Weak Variational Formulation can be defined for them at the same time. Yet, as its classic discretisation by plane waves is used to causing independence and accuracy limitations, we suggest to characterise local solutions by an impedance boundary conditions, which is discretised thanks to piecewise polynomial fields. Unfortunately, analytic expressions of such functions are not practically derivable, and we resort to a local numerical method so as to provide approximations: these quasi-solutions are finally chosen as basis functions in the Trefftz formulation, which is then referred as a quasi-Trefftz method.

In the end, numerical experiments allow to highlight robust independence, convergence and iterative behaviour of this method, under calibration conditions of the local solver with respect to the polynomial BC discretisation.

## 1 Introduction

The propagation of time-harmonic electromagnetic waves is a physical phenomenon which knows a various range of applications, from telecommunication and imaging to stealth. Even though dedicated numerical methods have long been specifically designed for the stakes of this class of problems (as the Finite Difference in Frequency Domain [1, 2, 3], the Finite Element Method [4, 5, 6, 7] or the ones lying on integral equations [8, 9]), the growing complexity of the practical configurations highlights the actual limits of these classic approaches.

Indeed, considering the use of high-frequency devices in large and complex propagation environments (with the positioning of a 4G/5G antenna in a city for example), current applications call for the globally accurate (in the whole domain at the same time) solution of such problems in wide domains (whose size can reach hundreds of wavelengths). However, this problematic is used to giving rise to the so-called 'pollution error' when dealing with wave problems [10], as a dispersion error is likely to appear when considering multiple wavelengths, unless a sufficiently large number of degrees of freedom is dedicated to represent the oscillation: this constraint then becomes more and more demanding when the domain size increases. Thus, the linear systems induced by numerical methods for these applications can quickly become unsuitable for direct solvers, while these classic approaches are usually poorly suited to iterative algorithms [11].

In this context, Trefftz approaches [12, 13] may offer an opportunity to deal with these stakes. In the continuity of the recent surge of interest offered by the introduction of the Ultra Weak Variational Formulation by Cessenat and Després [14, 15], this method relies on the use of local actual solutions as basis functions, instead of the classic polynomial approximations of Discontinuous Galerkin (DG) [16]. Its usual discretisation then lies on plane waves, as they are natural, analytic and smooth solutions of such wave problems, leading to thorough theoretical studies of their approximation and convergence properties [17, 18, 19, 20, 21, 22]. Nevertheless, these favourable characteristics come with numerical dependence phenomena, and then limited representation abilities
of such bases [23, 24]. Reduction techniques have then been developed to get rid of the spurious modes [24, 25], while other kinds of analytic solutions (as evanescent plane waves [26, 27, 28] or Generalised Harmonic Polynomials [17, 27]) were also assessed as basis functions.

In addition to these considerations, another point of view has also gained interest: considering that the constraint of using **actual** local solutions may be too demanding, it may be possible to loosen it by only resorting to **quasi-solutions**. Such a method is then referred as quasi-Trefftz.

One of the first kind of quasi-Trefftz method, referred as generalised plane waves, has been introduced by Imbert-Gérard and Després [29]. Its principle is, as the name suggests, to consider a generalisation of the classic plane waves $x \mapsto e^{i\kappa x}$ in the case of smooth non-constant coefficients. For example, in 1D, considering the equation

$$-p^{(2)} + \alpha p = 0 \quad \text{in } \Omega, \qquad (1)$$

this approach leans on two polynomial approximations in each cell $T$:

* the varying coefficient $\alpha$ is approximated thanks to a Taylor expansion $\tilde{\alpha} \in \mathbb{Q}_{q-1}(T)$ with respect to $x_G \in T$.

* as the classic phase of the plane wave can be seen as a polynomial function of degree 1 in $x$, it is more generally searched as a polynomial $P \in \mathbb{Q}_I(T)$, such that the associated basis function is $x \mapsto e^{P(x-x_G)}$. In this unidimensional case, the two basis functions are naturally searched as independent.

Replacing these forms in (1) then leads to a linear system which can be solved for judicious choices of the polynomial degrees $q$ and $I$ (and using the independence hypothesis of the basis functions) following an inductive procedure [30], while this principle can be extended to higher-dimensional cases.

More generally, this approach relies on a local approximation property with respect to a position $x_G \in T$: the idea is that, for any function $y$ satisfying the PDE $\mathcal{L}y = 0$, there exists a discrete approximation $y_h$ such that, for a given order $q \in \mathbb{N}^*$:

$$\forall \mathbf{x} \in T, \quad \mathcal{L}y_h(\mathbf{x}) = O(|\mathbf{x} - x_G|^q) \quad \text{or} \quad |y - y_h|(\mathbf{x}) = O(|\mathbf{x} - x_G|^q). \qquad (2)$$

Such a philosophy was then extended by considering a high-order polynomial development of the plane wave amplitude [31] of the form $x \mapsto P(x - x_G)e^{i\kappa(x-x_G)}$, or even purely polynomial functions as $x \mapsto P(x - x_G)$ [32]. In particular, these approaches share a common structure, highlighted in [32], leading to a generic methodology of computation of these polynomial functions thanks to an iterative process.

Nevertheless, the quasi-solution character of these functions can also be ensured by other means than a local property with respect to a given position, as detailed in (2). For example, Lehrenfeld and Stocker propose an embedded Trefftz DG method where this property is characterised in terms of orthogonality of the residual to a polynomial space [33], as $\mathbb{Q}_q$ for example. This idea then allows to keep the usual polynomial spaces used for discretisation in DG, which are then projected to only keep a subspace of quasi-solutions such that

$$\Pi \mathcal{L}y_h = 0 \quad \text{or} \quad \forall y'_h \in \mathbb{Q}_q(T), \quad \langle \mathcal{L}y_h, y'_h \rangle = 0, \qquad (3)$$

where $\Pi$ stands for the projection operator to the inner product $\langle \cdot, \cdot \rangle$ (which is usually chosen as the $L^2$ one).

In the vein of the stance of quasi-Trefftz methods, we propose to extend the principle proposed by Fure, Pernet, Sirdey and Tordeux in [34] for the 2D Maxwell's equations, to a general class of time-harmonic wave propagation problems. Characterising the local solutions by an impedance boundary conditions, their discretisation is then realised by resorting to a piecewise polynomial basis of this boundary field. Moreover, due to the practical unavailability of such basis functions (as no analytic expression can be derived), they are locally approximated thanks to a numerical method (under the form of **quasi-solutions**).

This article is organised as follows. We first introduce a common formalism to a class of time-harmonic wave problems, based on the two-fields Friedrichs systems, so as to tackle all of them at the same time. An associated Trefftz method, lying on the Ultra Weak Variational Formulation, is then defined while its classic discretisation with plane waves is discussed, and its independence and accuracy limitations are highlighted.

This leads us to formally introduce a discretisation of local solutions through local polynomial impedance boundary conditions, while an *a priori* error estimate is derived for the associated Trefftz method. Yet, the need to rely on an approximate local solver to get usable (quasi-solution) basis functions finally leads to a quasi-Trefftz method. This approach is then numerically discussed in terms of basis independence, in addition to $h$ and $p$-convergences,
allowing to calibrate the local solver with respect to the BC discretisation (so as to keep the best performances). In the end, the influence of the quasi-solution character of the basis functions on the iterative solution of the method is assessed, in order to verify its adequation to this mandatory resolution process when dealing with wide propagation domains. A final illustration, with the scattering of a plane wave by a metallic structure, is proposed, showing such an approach can overtake the accuracy saturation induced by the plane waves discretisation.

## 2 A Trefftz method for time-harmonic wave propagation problems

In this article, we are interested in time-harmonic source-free wave propagation problems, as the Maxwell and Helmholtz equations, for which we consider the temporal convention $t \mapsto e^{i\omega t}$, being given an angular frequency $\omega \in \mathbb{R}_*^+$. So as to deal with this full range of problems at the same time, we first introduce a common first-order formalism

$$ \mathcal{L}\mathbf{y} := i\kappa\mathcal{M}\mathbf{y} + \operatorname{div}(\mathcal{F}\mathbf{y}) = \mathbf{0} \quad \text{with} \quad \operatorname{div}(\mathcal{F}\mathbf{y}) := \sum_{j=1}^{d} \partial_{x_j}(\mathcal{F}_j\mathbf{y}) \quad \text{in } \Omega \subset \mathbb{R}^d, $$

where $\mathbf{y} : \Omega \to \mathbb{C}^m$, $\mathcal{M} : \Omega \to \mathbb{R}^{m \times m}$ and $\mathcal{F}_j \in \mathbb{R}^{m \times m}$. A Trefftz method is then adapted and described to this framework, lying on a polygonal/polyhedral mesh $\mathcal{T}_h$ of the domain $\Omega$, conformal with a piecewise heterogeneity distribution: we suppose there exists a partition of $\Omega$ in a finite set of $N_\Omega$ disjointed connected non-degenerate polygonal/polyhedral open subdomains $(\Omega_j)_{j \in [1, N_\Omega]}$ of Lipschitz boundary, such that each mesh cell is completely included in a unique subdomain. Thus, $\mathcal{M}$ is supposed to be constant in each mesh cell $T \in \mathcal{T}_h$ (while the solution field restriction is denoted as $\mathbf{y}^T := \mathbf{y}|_T$), while the flux matrices $(\mathcal{F}_j)_{j \in [1, d]}$ are fixed.

### 2.1 A unique formalism for time-harmonic wave propagation problems

In the idea that the problems we consider are posed on a couple of two variables (e and h for the Maxwell's equations; p and v for the Helmholtz equation), we propose to resort to a unique formalism, allowing to tackle all these problems at the same time (we refer to [35] for additional details).

#### 2.1.1 The Block Antidiagonal Two-Fields Friedrichs Systems

Following an adaptation of the Friedrichs systems [36, 37, 38] to the time-harmonic framework [39, 40], and relying on a subcase of the two-fields Friedrichs systems [41, 42], we suppose the unknown $\mathbf{y}$ can be split into two elements $\mathbf{y}_1$ and $\mathbf{y}_2$ of respective dimensions $m_1$ and $m_2$ such that $m_1 + m_2 = m$:

$$ \mathbf{y} = \begin{pmatrix} \mathbf{y}_1 \\ \mathbf{y}_2 \end{pmatrix}, \quad \mathbf{y}_1 : \Omega \to \mathbb{C}^{m_1} \quad \text{and} \quad \mathbf{y}_2 : \Omega \to \mathbb{C}^{m_2}. $$

Moreover, we suppose the matrix $\mathcal{M}$ can be written in a diagonal way as

$$ \forall \mathbf{x} \in \Omega, \quad \mathcal{M}(\mathbf{x}) = \begin{bmatrix} \varepsilon(\mathbf{x}) \mathcal{I}_{m_1} & \mathcal{O}_{m_1, m_2} \\ \mathcal{O}_{m_2, m_1} & \mu(\mathbf{x}) \mathcal{I}_{m_2} \end{bmatrix}, \qquad (4) $$

where $\mathcal{I}_{m_1}$ and $\mathcal{O}_{m_1, m_2}$ respectively stand for the identity and null matrices of $\mathbb{R}^{m_1 \times m_1}$ and $\mathbb{R}^{m_1 \times m_2}$, while the parameters $\varepsilon$ and $\mu$ satisfy a piecewise constant heterogeneity distribution (i.e. they are constant in each mesh cell) and are supposed to be strictly positive: being given a cell $T \in \mathcal{T}_h$, their associated values are referred as $\varepsilon(\mathbf{x}) = \varepsilon_T$ and $\mu(\mathbf{x}) = \mu_T$ respectively, such that $\mathcal{M}(\mathbf{x}) = \mathcal{M}_T$ for $\mathbf{x} \in T$. On the other side, the flux matrices take a symmetric block antidiagonal form

$$ \forall j \in [1, d], \quad \mathcal{F}_j = \begin{bmatrix} \mathcal{O}_{m_1, m_1} & \mathcal{B}_j \\ \mathcal{B}_j^* & \mathcal{O}_{m_2, m_2} \end{bmatrix}, $$

for fixed space-independent real matrices $\mathcal{B}_j \in \mathbb{R}^{m_1 \times m_2}$, and then, for any unit vector $\mathbf{n} \in \mathbb{R}^d$, one obtains the normal flux:

$$ \mathcal{F}_{\mathbf{n}} := \sum_{j=1}^{d} n_j \mathcal{F}_j = \begin{bmatrix} \mathcal{O}_{m_1, m_1} & \mathcal{B}_{\mathbf{n}} \\ \mathcal{B}_{\mathbf{n}}^* & \mathcal{O}_{m_2, m_2} \end{bmatrix} \quad \text{with} \quad \mathcal{B}_{\mathbf{n}} := \sum_{j=1}^{d} n_j \mathcal{B}_j \in \mathbb{R}^{m_1 \times m_2}. \qquad (5) $$

Thus, systems showing this structure are hereafter referred as Block Antidiagonal Two-Fields Friedrichs Systems (BATFFS).
**Remark 2.1.** In the specific case of the wave propagation problems we study in this article, they can be recovered as:

• for the Maxwell system, one has the electric $\mathbf{y}_1 = \mathbf{e} : \Omega \to \mathbb{C}^3$ and magnetic fields $\mathbf{y}_2 = \mathbf{h} : \Omega \to \mathbb{C}^3$, with $m_1 = m_2 = 3$. Indeed, one has $\mathcal{B}_j = -\Gamma_\times[\mathbf{d}_j] \in \mathbb{R}^{3 \times 3}$ where $(\mathbf{d}_1, \mathbf{d}_2, \mathbf{d}_3)$ denotes the canonical basis of $\mathbb{R}^3$ and, being given a vector $\mathbf{n} := (n_1, n_2, n_3) \in \mathbb{R}^3$, the matrix $\Gamma_\times[\mathbf{n}]$ denotes

$$ \Gamma_\times[\mathbf{n}] = \begin{bmatrix} 0 & -n_3 & n_2 \\ n_3 & 0 & -n_1 \\ -n_2 & n_1 & 0 \end{bmatrix}. \qquad (6) $$

• for the 2D Helmholtz equation, one has the pressure $y_1 = p : \Omega \to \mathbb{C}$ and velocity $\mathbf{y}_2 = \mathbf{v} = (i\kappa\mu)^{-1} \nabla p : \Omega \to \mathbb{C}^2$ fields, with $m_1 = 1$ and $m_2 = 2$. Indeed, one has $\mathcal{B}_j = -\mathbf{d}_j^* \in \mathbb{R}^{1 \times 2}$.

Moreover, so as to ensure the well-posedness of the PDE and the following numerical methods, we suppose that a unique continuation principle stands.

**Hypothesis 2.2.** We suppose the following unique continuation principle in any cell $T \in \mathcal{T}_h$:

let $\mathbf{y}^T \in [L^2(T)]^m$ such that

$$ \left\{ \begin{array}{ll} i\kappa \mathcal{M}_T \mathbf{y}^T + \operatorname{div}(\mathcal{F} \mathbf{y}^T) = \mathbf{0} & \text{in } T, \\ \mathcal{F}_{n_T} \mathbf{y}^T = \mathbf{0} & \text{on } \Gamma, \end{array} \right. \qquad (7) $$

where $\Gamma \subset \partial T$ is a subset of $\partial T$ of non-zero measure. Then, one has $\mathbf{y}^T = \mathbf{0}$ a.e. in $T$.

**Remark 2.3.** In the case one has $\Gamma = \partial T$ for general BATFFS, one can refer to [40]. For a general $\Gamma$, considering the wave problems of interest, Hypothesis 2.2 is verified for Maxwell [7] and Helmholtz [43] systems.

Then, it remains to detail boundary conditions adapted to the BATFFS formalism we have just introduced.

To do so, being given a unit vector $\mathbf{n} \in \mathbb{R}^d$, we consider a compact singular value decomposition of the matrix $\mathcal{B}_n$ (defined in (5)), such that

$$ \mathcal{B}_n = \hat{U}\Sigma\hat{V}^* \quad \text{with} \quad \hat{U}^*\hat{U} = \mathcal{I}_r \quad \text{and} \quad \hat{V}^*\hat{V} = \mathcal{I}_r, \qquad (8) $$

where $r \le \min\{m_1, m_2\}$ is the rank of $\mathcal{B}_n$, the orthonormal left and right-singular vectors are denoted as $\hat{U} \in \mathbb{R}^{m_1 \times r}$, $\hat{V} \in \mathbb{R}^{m_2 \times r}$, and $\Sigma \in \mathbb{R}^{r \times r}$ is a diagonal matrix which contains the strictly positive singular values of $\mathcal{B}_n$.

**Remark 2.4.** In particular, even if such operators are defined relatively to the unit vector $\mathbf{n}$, they can be bounded in an independent way to $\mathbf{n}$. Indeed, by relying on the definition of the norm $\|\mathcal{B}_n\mathcal{B}_n^*\|$ in terms of Rayleigh quotient, and by using the fact $\hat{U}\hat{U}^*$ is Hermitian and a projection, one gets

$$ \|\Sigma\|^2 = \|\mathcal{B}_n\mathcal{B}_n^*\| \le \left[ \sum_{j=1}^d \|\mathcal{B}_j\| \right]^2 \quad \text{and} \quad \|\hat{U}^*\| = 1. $$

Moreover, so as to simplify the notations, we rather introduce the symmetric positive definite square root $\Sigma^{1/2}$ of $\Sigma$, allowing to report the singular values into the singular vectors. Thus, leaning on (8), we define the following compact singular value decomposition of $\mathcal{B}_n$:

$$ \mathcal{B}_n = U\Sigma^* \quad \text{with} \quad \left\{ \begin{array}{ll} U = \hat{U}\Sigma^{1/2} & \text{and} \quad V = \hat{V}\Sigma^{1/2}, \\ U^*U = \Sigma & \text{and} \quad V^*V = \Sigma, \end{array} \right. \qquad (9) $$

where $U \in \mathbb{R}^{m_1 \times r}$ and $V \in \mathbb{R}^{m_2 \times r}$ are then simply orthogonal (as we recall that $\Sigma$ is a diagonal matrix). As it avoids to call for $\Sigma$, we will stick to this decomposition for the rest of the article (even if (8) may have been preferred, as it would keep orthonormality for the singular vectors).

**Remark 2.5.** In the specific case of the wave propagation problems we study in this article, one gets $\Sigma = \mathcal{I}_r$ and:
* for the Maxwell system, one has $U^*\mathbf{y}_1 = \gamma_t[\mathbf{n}](\mathbf{e})$ and $V^*\mathbf{y}_2 = -\gamma_\times[\mathbf{n}](\mathbf{h})$ (with $r=2$), where

$$ \gamma_t[\mathbf{n}](\mathbf{w}) = \mathbf{w} - (\mathbf{w} \cdot \mathbf{n}) \mathbf{n} = (\mathbf{n} \times \mathbf{w}) \times \mathbf{n} \quad \text{and} \quad \gamma_\times[\mathbf{n}](\mathbf{w}) = \mathbf{n} \times \mathbf{w}, \qquad (10) $$

or more specifically, their projection over the two-dimensional space of tangent vectors.

* for the 2D Helmholtz equation, one has $U^*\mathbf{y}_1 = p$ and $V^*\mathbf{y}_2 = -v_N := -\mathbf{v} \cdot \mathbf{n} = -(i\kappa\mu)^{-1}\partial_\mathbf{n}p$, with $r=1$.

**Remark 2.6.** We highlight the fact the compact singular value decomposition (9) of $B_{n_T}$ must be realised on each face F of the boundary cell $\partial T$ (where the unit outward normal vector $\mathbf{n}_T$ is constant), leading to $B_{n_T} = U_T^F(V_T^F)^*$ on F. Nevertheless, so as to simplify the notations, the subscript $\cdot_T^F$ will be presupposed without further mention. Similarly, the rank $r$ of $B_{n_T}$ must also be computed for each face F, and may evolve along the cell boundary $\partial T$. Yet, so as to simplify the notations, we stick to the simple denotation $r$, presupposing it is relative to each face of each cell.

Leaning on these variables, it is then possible to define notions of outgoing and incoming traces (see Subsection 2.2.1) and to recover classic impedance boundary conditions: all these notions allow to completely define a general class of time-harmonic problems, in the formalism of BATFFS.

**Problem 2.7 (BATFFS).** Being given a wavenumber $\kappa > 0$, find $\mathbf{y} \in [L^2(\Omega)]^m$ such that

$$ \mathcal{L}\mathbf{y}(\mathbf{x}) := i\kappa M\mathbf{y}(\mathbf{x}) + \operatorname{div}(F\mathbf{y})(\mathbf{x}) = \mathbf{0} \quad \text{for } \mathbf{x} \in \Omega, \qquad (11) $$

that we equip with an impedance boundary condition, for a non uniformly zero impedance field $Z_{\partial\Omega} : \partial\Omega \to \mathbb{C}$

$$ \gamma_{in}^{\partial\Omega}(\mathbf{y})(\mathbf{x}) = U^*\mathbf{y}_1(\mathbf{x}) - Z_{\partial\Omega}(\mathbf{x}) V^*\mathbf{y}_2(\mathbf{x}) = \mathbf{g}(\mathbf{x}) \quad \text{for } \mathbf{x} \in \partial\Omega, \qquad (12) $$

for a boundary field $\mathbf{g} \in [L^2(\partial\Omega)]^r$.

**Remark 2.8.** In the specific case of the wave propagation problems we study in this article, we recover the classic impedance BCs by referring to Remark 2.5 and (12):

* for the Maxwell system, one gets $\gamma_{in}^{\partial\Omega}(\mathbf{y}) = \gamma_t[\mathbf{n}_{\partial\Omega}](\mathbf{e}) + Z_{\partial\Omega} \gamma_\times[\mathbf{n}_{\partial\Omega}](\mathbf{h})$.

* for the Helmholtz equation, one gets $\gamma_{in}^{\partial\Omega}(\mathbf{y}) = p + Z_{\partial\Omega} v_N$.

### 2.1.2 Explicit decomposition of the normal flux and numerical traces

The main interest of the BATFFS formalism is to allow an explicit derivation of a decomposition of the normal flux into incoming and outgoing ones (and then a consistent numerical flux), similarly to what is classically realised through Riemann solvers [44, Section 2.3.2]. The classic expressions for Maxwell and Helmholtz equations can then be recovered.

First, the decomposition of the normal flux $\mathcal{F}_n$ can be realised thanks to the generalised eigenvalue problem

$$ \mathcal{F}_n \mathbf{x} = \lambda \tilde{\mathcal{M}} \mathbf{x} \quad \text{with} \quad \tilde{\mathcal{M}} = \begin{bmatrix} \tilde{\varepsilon} \mathcal{I}_{m_1} & \mathcal{O}_{m_1, m_2} \\ \mathcal{O}_{m_2, m_1} & \tilde{\mu} \mathcal{I}_{m_2} \end{bmatrix}, $$

where we consider generic values $\tilde{\varepsilon}, \tilde{\mu} \in \mathbb{R}_*^+$. This leads, by denoting the admittance and impedance as $\tilde{Y} = \sqrt{\tilde{\varepsilon}/\tilde{\mu}}$ and $\tilde{Z} = \sqrt{\tilde{\mu}/\tilde{\varepsilon}}$ respectively, to

$$ \mathcal{F}_n = \mathcal{F}_n^+ + \mathcal{F}_n^- \quad \text{with} \quad \mathcal{F}_n^{\pm} = \frac{1}{2} \begin{bmatrix} \pm \tilde{Y} \mathbf{u}\mathbf{u}^* & \mathbf{u}\mathbf{v}^* \\ \mathbf{v}\mathbf{u}^* & \pm \tilde{Z} \mathbf{v}\mathbf{v}^* \end{bmatrix}. \qquad (13) $$

where the matrices $\mathcal{F}_n^\pm$ are, by construction, positive and negative respectively.
Lying on this explicit decomposition, one can look for a consistent numerical flux $\overline{\mathcal{F}_{n_T}y}$ (at the interface between two cells T and K) as a linear combination of outgoing data $\mathcal{F}_{n_T}^+y^T$ and $\mathcal{F}_{n_K}^+y^K$. Its coefficients are uniquely determined by imposing consistency and using continuity properties of the exact solution [45]:

$$ \mathcal{F}_{n_T}y^T + \mathcal{F}_{n_K}y^K = 0, \quad \text{such that} \quad U^*y_1^T = U^*y_1^K \quad \text{and} \quad V^*y_2^T = V^*y_2^K \quad \text{a.e. on } \partial T \cap \partial K, $$

where $U$ and $V$ are associated to the compact SVD (9) of the normal flux operator $\mathcal{B}_{n_T}$ (which is arbitrarily favoured hereafter). This finally leads to

$$ \overline{\mathcal{F}_{n_T}y} = \begin{pmatrix} U\overline{V^*y_2} \\ V\overline{U^*y_1} \end{pmatrix} = \frac{1}{Z_T + Z_K} \begin{pmatrix} U[\mathcal{U}^*(y_1^T - y_1^K) + \mathcal{V}^*(Z_T y_2^T + Z_K y_2^K)] \\ V[\mathcal{U}^*(Z_K y_1^T + Z_T y_1^K) + Z_T Z_K \mathcal{V}^*(y_2^T - y_2^K)] \end{pmatrix}, \quad (14) $$

where $Z_T := \sqrt{\mu_T/\epsilon_T}$ and $Y_T := \sqrt{\epsilon_T/\mu_T}$ respectively denote the relative impedance and admittance of T.

**Remark 2.9.** In the case of a boundary face F, K can be interpreted as a ghost cell lying outside of $\Omega$ that shares F with T. By choosing $Z_K = Z_{\partial\Omega}$, the quantity $U^*y_1^K - Z_K V^*y_2^K = \gamma_{in}^{\partial\Omega}(y^K)$ appearing in the numerical flux (14), can then be replaced by the boundary condition (12): this allows to naturally define a numerical flux for boundary faces.

In particular, relying on Remark 2.5, classic expressions of the numerical fluxes are recovered (see [25] for the Maxwell's equations for example).

## 2.2 Introduction of a Trefftz method for BATFFS

Once the formalism of BATFFS has been introduced, it is possible to define an adapted Trefftz method. Such an approach can be obtained by relying on an Ultra Weak Variational Formulation (UWVF): we detail the approximation spaces it lies on, the formulation in the case of BATFFS, in addition to its noteworthy properties.

In particular, the proofs of the latter are not derived here for the sake of conciseness but are natural extensions of the ones for Maxwell and Helmholtz equations, to which we refer. In particular, a thorougher introduction and proofs are all gathered in [35].

### 2.2.1 Approximation spaces

The Trefftz method can be interpreted as a Discontinuous Galerkin method, whose specificity is to call for basis functions which are actual solutions of the PDE in each cell. Then, we begin with the introduction of the associated spaces, on which relies the study of this method.

**Definition 2.10.** For a cell $T \in \mathcal{T}_h$, we define the local and global spaces of boundary-regular local solutions

$$ Y^T := \{ y^T \in [L^2(T)]^m : \mathcal{L}y^T = 0 \text{ in } T \text{ such that } \mathcal{F}_{n_T}y^T \in [L^2(\partial T)]^m \} \quad \text{and} \quad Y := \prod_{T \in \mathcal{T}_h} Y^T, $$

where we recall that $\mathcal{L}$ stands for the Friedrichs operator (see (11)). They are equipped with the inner products (see Remark 2.11)

$$ \forall \mathbf{y}, \mathbf{y}' \in \mathbb{Y}, \quad \langle \mathbf{y}, \mathbf{y}' \rangle_{\mathbb{Y}} := \sum_{T \in \mathcal{T}_h} \langle \mathbf{y}^T, \mathbf{y}'^T \rangle_{\mathbb{Y}^T} \quad \text{with} \quad \langle \mathbf{y}^T, \mathbf{y}'^T \rangle_{\mathbb{Y}^T} := Y_T \langle \mathcal{U}^* \mathbf{y}_1^T, \mathcal{U}^* \mathbf{y}_1'^T \rangle_{0, \partial T} + Z_T \langle \mathcal{V}^* \mathbf{y}_2^T, \mathcal{V}^* \mathbf{y}_2'^T \rangle_{0, \partial T}, \quad (15) $$

where $\langle \cdot, \cdot \rangle_{0, \partial T}$ stands for the $[L^2(\partial T)]^T$ scalar product, and we recall that $U$ and $V$ refer to the compact singular value decomposition (9) issued from $\mathcal{B}_{n_T}$.

**Remark 2.11.** In each cell $T \in \mathcal{T}_h$, the symmetric positive sesquilinear form $\langle \cdot, \cdot \rangle_{\mathbb{Y}^T}$ is also definite on $\mathbb{Y}^T$. Indeed, having $\|\mathbf{y}^T\|_{\mathbb{Y}^T} = 0$ implies that $\mathcal{F}_{n_T}y^T = 0$ on $\partial T$. Then, relying on the unique continuation principle (see Hypothesis 2.2) finally ensures that $\mathbf{y}^T = 0$ in T, which can be summarised as:

$$ \|\mathbf{y}^T\|_{\mathbb{Y}^T} = 0 \quad \Rightarrow \quad \mathcal{F}_{n_T}y^T = 0 \text{ on } \partial T \quad \Rightarrow \quad \mathbf{y}^T = 0 \text{ in } T. $$
In particular, functions belonging to such a space satisfy the following reciprocity formula.

**Proposition 2.12** (Reciprocity formula). Let's consider a cell $T \in \mathcal{T}_h$, in addition to $\mathbf{y}^T, \mathbf{y}'^T \in \mathbb{Y}^T$. Then, one has:

$$r_T(\mathbf{y}^T, \mathbf{y}'^T) := \langle \mathcal{F}_{n_T} \mathbf{y}^T, \mathbf{y}'^T \rangle_{0, \partial T} = 0, \qquad (16)$$

which can be rewritten for BATFFS as

$$r_T(\mathbf{y}^T, \mathbf{y}'^T) = \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1'^T \rangle_{0, \partial T} + \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2'^T \rangle_{0, \partial T} = 0. \qquad (17)$$

*Proof.* This proposition relies on the integration-by-parts formula (see [46, Section 3.1.2] and [40]), stating that

$$r_T(\mathbf{y}^T, \mathbf{y}'^T) = \int_T \mathrm{div}(\mathcal{F}\mathbf{y}^T) \cdot \overline{\mathbf{y}^T} \, d\mathbf{x} + \int_T \mathbf{y}^T \cdot \mathrm{div}(\overline{\mathcal{F}\mathbf{y}^T}) \, d\mathbf{x}.$$

By using the fact that both $\mathbf{y}^T$ and $\mathbf{y}'^T$ are local solutions, one gets

$$r_T(\mathbf{y}^T, \mathbf{y}'^T) = \int_T \left(-i\kappa \mathcal{M} \mathbf{y}^T\right) \cdot \overline{\mathbf{y}^T} \, d\mathbf{x} + \int_T \mathbf{y}^T \cdot \left(\overline{-i\kappa \mathcal{M} \mathbf{y}^T}\right) \, d\mathbf{x},$$

and the real symmetric character of $\mathcal{M}$ allows to conclude. $\square$

Instead of relying on such functions that live in the full cell, the UWVF point of view leans on their outgoing and incoming traces, which are restricted to the cell boundary and used as solving variables, allowing to condense the problem on the mesh skeleton: we introduce their definition in this generalised framework.

**Definition 2.13.** For a cell $T \in \mathcal{T}_h$, we introduce the outgoing and incoming traces $\mathbf{x}^T$ and $\mathcal{U}^T \mathbf{x}^T$, associated to $\mathbf{y}^T \in \mathbb{Y}^T$, as

$$\left\{ \begin{array}{l} \mathbf{x}^T := \gamma_{out}^T(\mathbf{y}^T) = \mathcal{U}^* \mathbf{y}_1^T + Z_T \mathcal{V}^* \mathbf{y}_2^T \quad \text{on } \partial T, \\ \mathcal{U}^T \mathbf{x}^T := \gamma_{in}^T(\mathbf{y}^T) = \mathcal{U}^* \mathbf{y}_1^T - Z_T \mathcal{V}^* \mathbf{y}_2^T \end{array} \right. \qquad (18)$$

where we refer to (9) for a definition of the decomposition of $\mathcal{B}_{n_T}$ in terms of $\mathcal{U}$ and $\mathcal{V}$.

Moreover, the expected regularities of $\mathbf{y}^T \in \mathbb{Y}^T$ and $\mathbf{y} \in \mathbb{Y}$ ensure that

$$\mathbf{x}^T, \mathcal{U}^T \mathbf{x}^T \in \mathbb{X}^T := \gamma_{out}^T(\mathbb{Y}^T) = \left[L^2(\partial T)\right]^T \quad \text{and} \quad \mathbf{x}, \mathcal{U}\mathbf{x} \in \mathbb{X} := \prod_{T \in \mathcal{T}_h} \mathbb{X}^T = \gamma_{out}(\mathbb{Y}) = \left[L^2(\partial \mathcal{T}_h)\right]^T, \qquad (19)$$

where we recall that $r$ is the rank of $\mathcal{B}_{n_T}$, and the global outgoing $\mathbf{x} = \gamma_{out}(\mathbf{y}) \in \mathbb{X}$ and incoming traces $\mathcal{U}\mathbf{x} = \gamma_{in}(\mathbf{y}) \in \mathbb{X}$ are defined through their restrictions

$$\forall T \in \mathcal{T}_h, \quad \mathbf{x}^T = \gamma_{out}^T(\mathbf{y}^T) \quad \text{and} \quad \mathcal{U}^T \mathbf{x}^T = \gamma_{in}^T(\mathbf{y}^T).$$

Concerning the whole domain $\Omega$, this allows to define, for $\mathbf{y} \in \mathbb{Y}$, the global outgoing $\mathbf{x} = \gamma_{out}(\mathbf{y}) \in \mathbb{X}$ and incoming $\mathcal{U}\mathbf{x} = \gamma_{in}(\mathbf{y}) \in \mathbb{X}$ traces such that

$$\forall T \in \mathcal{T}_h, \quad \mathbf{x}^T = \gamma_{out}^T(\mathbf{y}^T) \quad \text{and} \quad \mathcal{U}^T \mathbf{x}^T = \gamma_{in}^T(\mathbf{y}^T).$$

These spaces are equipped with the inner products

$$\forall \mathbf{x}, \mathbf{x}' \in \mathbb{X}, \quad \langle \mathbf{x}, \mathbf{x}' \rangle_{\mathbb{X}} := \sum_{T \in \mathcal{T}_h} \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathbb{X}^T} \quad \text{with} \quad \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathbb{X}^T} := Y_T \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{0, \partial T}. \qquad (20)$$

**Remark 2.14.** These notions can be illustrated in the framework of plane waves (whose general definition is detailed in Subsection 2.3). Considering a propagation direction going out of a cell $T \in \mathcal{T}_h$, the incoming trace of the associated plane wave is then zero, as no information enters the domain but everything goes out: $\gamma_{in}^T(\mathbf{y}_{PW}) = \mathbf{0}$.

In particular, lying on the definition (18) of the incoming trace and the numerical traces (14), one can define the following numerical incoming trace:

$$\widehat{\mathcal{U}^T \mathbf{x}} = \widehat{\mathcal{U}^* \mathbf{y}_1} - Z_T \widehat{\mathcal{V}^* \mathbf{y}_2} = \begin{cases} \frac{Z_K - Z_T}{Z_T + Z_K} \mathbf{x}^T & + \frac{2Z_T}{Z_T + Z_K} \mathbf{x}^K \quad \text{on } \partial T \cap \partial K, \\ \frac{Z_{\partial \Omega} - Z_T}{Z_T + Z_{\partial \Omega}} \mathbf{x}^T & + \frac{2Z_T}{Z_T + Z_{\partial \Omega}} \mathbf{g} \quad \text{on } \partial T \cap \partial \Omega. \end{cases} \qquad (21)$$

Naturally, characterising a function by its outgoing/incoming trace calls for the well-posedness of the associated problem, that we suppose in the rest of the article.
**Hypothesis 2.15.** We suppose that, in every cell $T \in \mathcal{T}_h$, the local problem associated to the time-harmonic Friedrichs system and the outgoing/incoming traces is well posed.

More precisely, we suppose that for any $\mathbf{x}^T \in \mathcal{X}^T$ (respectively $\mathfrak{U}^T \mathbf{x}^T \in \mathcal{X}^T$), there exists a unique $S_{out}^T(\mathbf{x}^T) \in \mathcal{Y}^T$ (respectively $S_{in}^T(\mathfrak{U}^T \mathbf{x}^T) \in \mathcal{Y}^T$) such that $\gamma_{out}^T(S_{out}^T(\mathbf{x}^T)) = \mathbf{x}^T$ (respectively $\gamma_{in}^T(S_{in}^T(\mathfrak{U}^T \mathbf{x}^T)) = \mathfrak{U}^T \mathbf{x}^T$) on $\partial T$.

**Remark 2.16.** In the specific case of the wave propagation problems we study in this article, such local solvers and traces are well posed as:

* for the Maxwell system, one recovers the classic expressions [15, 44]:

$$ \mathbf{x}^T = \gamma_t(\mathbf{e}) - Z_T \gamma_{\times}^T(\mathbf{h}) \quad \text{and} \quad \mathfrak{U}^T \mathbf{x}^T = \gamma_t(\mathbf{e}) + Z_T \gamma_{\times}^T(\mathbf{h}). $$

* for the 2D Helmholtz equation, one recovers the classic expressions [15, 47]:

$$ \mathbf{x}^T = p^T - Z_T v_N^T \quad \text{and} \quad \mathfrak{U}^T \mathbf{x}^T = p^T + Z_T v_N^T. $$

### 2.2.2 A Trefftz formulation for BATFFS

Lying on the notions of outgoing/incoming traces, an Ultra Weak Variational Formulation [15, 47, 25] can then be adapted to BATFFS.

**Problem 2.17 (UWVF for BATFFS).** Find $\mathbf{x} \in \mathcal{X}$ such that for all $\mathbf{x}' \in \mathcal{X}$, one has:

$$ \tilde{a}_{UWVF}(\mathbf{x}, \mathbf{x}', \mathbf{g}) := \sum_{T \in \mathcal{T}_h} \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathcal{X}^T} - \langle \widehat{\mathfrak{U}^T \mathbf{x}^T}, \mathfrak{U}^T \mathbf{x}'^T \rangle_{\mathcal{X}^T} = 0. \quad (22) $$

Naturally, by separating the terms calling for the boundary condition, this problem can equivalently be stated as

$$ a_{UWVF}(\mathbf{x}, \mathbf{x}') = \ell_{UWVF}(\mathbf{x}'), \quad (23) $$

where $a_{UWVF}$ is a sesquilinear form defined for $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$ as

$$ a_{UWVF}(\mathbf{x}, \mathbf{x}') := \langle \mathbf{x}, \mathbf{x}' \rangle_{\mathcal{X}} - k_{UWVF}(\mathbf{x}, \mathbf{x}'), \quad (24) $$

with (see (21) for the definition of the numerical incoming trace)

$$ k_{UWVF}(\mathbf{x}, \mathbf{x}') := \langle \Pi_{\mathfrak{U}} \mathbf{x}, \mathfrak{U} \mathbf{x}' \rangle_{\mathcal{X}} \quad \text{and} \quad \Pi_{\mathfrak{U}}^T \mathbf{x}^T = \begin{cases} \widehat{\mathfrak{U}^T \mathbf{x}}|_F & \text{if } F \in \mathcal{F}_{\text{int}} \\ \frac{Z_{\partial\Omega} - Z_T}{Z_T + Z_{\partial\Omega}} \mathbf{x}^T & \text{if } F \in \mathcal{F}_{\text{ext}} \end{cases}, \quad (25) $$

where $\mathcal{F}_{\text{int}}$ and $\mathcal{F}_{\text{ext}}$ respectively denote the internal and boundary interfaces of the mesh, and $\ell_{UWVF}$ is a conjugate-linear form which encompasses the boundary conditions

$$ \ell_{UWVF}(\mathbf{x}') := \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_{\text{ext}}(T)} \left\langle \frac{2Z_T}{Z_T + Z_{\partial\Omega}} \mathbf{g}, \mathfrak{U}^T \mathbf{x}'^T \right\rangle_{\mathcal{X}^T, F}, $$

where $\mathcal{F}_{\text{ext}}(T)$ stands for the faces of $\partial T$ which also belong to the boundary $\partial\Omega$ (and then $\mathcal{F}_{\text{ext}}$).

In terms of noteworthy properties, this UWVF formulation can be proved to be equivalent to another one (see [25, Theorem 2.1] for the Maxwell system), derived from the reciprocity formula, and which allows one to prove a weak-coercivity result (see [25, Proposition 1.3] for the Maxwell system), with respect to a control norm based on the jumps on the mesh skeleton. In particular, this ensures the well-posedness of the method.

Moreover, this UWVF point of view allows one to naturally highlight a preconditioner (the mass matrix associated to the X-scalar product) which leads to a (strictly in the discrete case) contracting linear system to solve (see [47, Proposition 1.10 and Lemma 2.5] for the Helmholtz equation and [25, Proposition 3.1] for the Maxwell system). This property justifies the strong adaptability of the Trefftz approach to iterative algorithms.
## 2.3 Classic discretisation and numerical limits

In terms of discretisation, we classically resort to conformal finite-dimensional subspaces $\mathbb{Y}_h \subset \mathbb{Y}$ and $\mathbb{X}_h \subset \mathbb{X}$ of the ones introduced in Subsection 2.2.1, such that one has: $\mathbb{X}_h = \gamma_{out}(\mathbb{Y}_h)$, just as $\mathbb{X} = \gamma_{out}(\mathbb{Y})$ (see (19)). It remains to replace $\mathbb{X}$ by $\mathbb{X}_h$ in Problem 2.17.

Nevertheless, one must recall the specificity of the approximation spaces of the Trefftz method: they call for functions which are **actual local solutions**. Then, so as to discretise the formulation, one needs to define, in each cell $T \in \mathcal{T}_h$, a discrete set of such local solutions $(\phi_j^T)_{j \in [1, N_{dof}^T]} \subset \mathbb{Y}^T$ and their associated outgoing traces $(\psi_j^T = \gamma_{out}^T(\phi_j^T))_{j \in [1, N_{dof}^T]} \subset \mathbb{X}^T$.

In the case of the wave propagation problems, a range of classic functions may be considered: Generalised Harmonic Polynomials [17, 27]; Bessel functions, which allow to take singularities into account, at the cost of specific integration techniques; or, classically, the propagative plane waves [15, 21, 44] (in opposition to the evanescent ones [26, 27, 28]). Indeed, such functions are natural smooth solutions of wave propagation problems, with an easily computable analytical expression which does not require sophisticated quadrature rules for integral approximations. In the general case of time-harmonic Friedrichs systems, a plane wave $\mathbf{y}_{PW}$ can be seen as a solution associated to a specific propagation direction $\mathbf{d} \in \mathbb{R}^d$ which is searched, for $\lambda \in \mathbb{R}^*$ and an amplitude $\mathbf{Y}_{PW} \in \mathbb{C}^m$, in the form

$$ \forall \mathbf{x}, \quad \mathbf{y}_{PW}(\mathbf{x}) := \mathbf{Y}_{PW} \exp(-i\kappa\lambda^{-1}\mathbf{d} \cdot \mathbf{x}), \quad \text{such that} \quad \mathcal{F}_{\mathbf{d}}\mathbf{Y}_{PW} = \lambda \mathcal{M}\mathbf{Y}_{PW}. $$

Unfortunately, their use in numerical methods is used to implying stability issues because of numerical dependence phenomena. Indeed, considering distinct propagation directions theoretically ensures the associated plane waves should be linearly independent [48, 49]. Thus, the local mass matrix $\mathcal{D}^T$ (which is supposed to serve as preconditioner) defined as

$$ \forall i, j \in [1, N_{dof}^T], \quad \mathcal{D}_{ij}^T = \langle \psi_j^T, \psi_i^T \rangle_{\mathbb{X}^T} = \langle \gamma_{out}^T(\phi_j^T), \gamma_{out}^T(\phi_i^T) \rangle_{\mathbb{X}^T}, \quad (26) $$

should be invertible for any value of $N_{dof}^T$.

Yet, the numerical discretisation of the plane waves is usually responsible for the loss of this theoretical independence [27, 23, 24]: being given a large enough number of plane waves, adding another one does not allow to get new information but only redundant one. This phenomenon is illustrated in Figure 1, where the condition number of the mass matrix is presented with respect to the number $N_{dof}^T$ of plane waves in a cell of size 1 wavelength: it skyrockets when the number of basis functions increases, until preventing the numerical invertibility of $\mathcal{D}^T$.

Facing this phenomenon, some reduction techniques have been developed to withdraw the spurious modes, by only retaining the largest eigenvalues (and their associated eigenvectors) [24, 25]: even if this principle allows to reduce the size of the basis with minimum accuracy loss, it does not prevent the precision saturation, as no new information is brought when $N_{dof}^T$ increases. This assessment calls for the introduction of naturally well-conditioned Trefftz basis functions, that we discuss in Section 3.

# 3 A new set of (quasi-)Trefftz functions characterised by local polynomial impedance boundary conditions

Facing the dependence limitation of plane waves as Trefftz basis functions, we propose and analyse a new discretisation process allowing to define a naturally well-conditioned basis, lying on the extension of the approach introduced by Fure, Pernet, Sirdey and Tordeux [34] for the 2D Maxwell's equations. The idea is then to characterise the local solutions by a local impedance boundary condition, that is discretised thanks to polynomial functions. Furthermore, due to the fact they lean on a local solver that cannot be used in practice, we suggest to resort to an approximate one, leading to a quasi-Trefftz approach.

In particular, for the sake of presentation, we introduce this method in the framework of a Cartesian mesh, whose cell size is bounded by a quantity $h$, even if the following principles and proofs could be straightforwardly adapted to more general configurations.

## 3.1 Local problem and polynomial BCs discretisation

We begin with the introduction, for each mesh cell, of a local problem equipped with specific impedance boundary conditions, in addition to their polynomial discretisation: this point of view then allows to define a finite-dimensional
<figure><img src="image_2.png" /><figcaption>Figure 1: Distribution of the eigenvalues of the local mass matrix $D^T$ for 36 plane waves (2D Helmholtz) and a cell of size 1 wavelength; evolution of the condition number of the local mass matrix $D^T$ with respect to the number of plane waves, for a cell of size 1 wavelength.</figcaption></figure>

space of Trefftz functions.

### 3.1.1 Local impedance field and problem

Let's consider a cell $T \in \mathcal{T}_h$, whose set of faces is denoted $\mathcal{F}(T)$. We introduce a local impedance boundary field $Z_{\partial T}: \partial T \to \mathbb{C}$, that we suppose to be non uniformly zero and constant on each face of $\mathcal{F}(T)$.

Then, we define the local solution operator $S_T$ which solves the local Friedrichs system for parameters $\varepsilon_T$ and $\mu_T$, with an impedance boundary condition leaning on the boundary field $Z_{\partial T}$:

**Problem 3.1.** We consider the local solution operator $S_T$ which associates, to any $\mathbf{g}^T \in [L^2(\partial T)]^r$, the solution $\mathbf{y}^T = S_T(\mathbf{g}^T) \in \mathbb{Y}^T$ such that one has

$$ \left\{ \begin{array}{ll} \mathcal{L}\mathbf{y}^T = i\kappa M_T \mathbf{y}^T + \operatorname{div}(\mathcal{F}\mathbf{y}^T) = \mathbf{0} & \text{a.e. in } T, \\ \gamma_{in}^{\partial T}(\mathbf{y}^T) = \mathcal{U}^*\mathbf{y}_1^T - Z_{\partial T} \mathcal{V}^*\mathbf{y}_2^T = \mathbf{g}^T & \text{a.e. on } \partial T. \end{array} \right. \qquad (27) $$

Naturally, being given $\mathbf{g} = (\mathbf{g}^T)_{T \in \mathcal{T}_h} \in [L^2(\partial \mathcal{T}_h)]^r$, we denote the global solution as $\mathbf{y} = S(\mathbf{g}) \in \mathbb{Y}$ such that

$$ \forall T \in \mathcal{T}_h, \quad \mathbf{y}^T = S_T(\mathbf{g}^T). $$

**Remark 3.2.** As we recall the impedance $Z: \Omega \to \mathbb{R}_*^+$ is constant in each cell (and is denoted as $Z_T$), we propose three exemplary choices of this boundary impedance field:

(a) first, it may be chosen as the cell one for example, meaning that $Z_{\partial T} = Z_T = \sqrt{\mu_T/\varepsilon_T}$. The impedance boundary condition will then match the incoming trace (18).

(b) then, another possibility is to consider, on each face $F \in \mathcal{F}(T)$, the impedance $Z_K$ of the neighbouring cell $K$ (and $Z_{\partial \Omega}$ in the case of a boundary cell): $Z_{\partial T} = Z_K$ on $F = \partial T \cap \partial K$.

(c) finally, one may rely on a constant and uniform value for the cell, and then independent of the local parameters, as $Z_{\partial T} = C_T \in \mathbb{C}$.

In the framework of a Cartesian mesh where all cell dimensions are identical, options (a) and (c) then lead the following local problem to be identical in all homogeneous subdomains.

**Remark 3.3.** Generally, the impedance field $Z_{\partial T}$ is asked to satisfy $\Re[Z_{\partial T}] \ge Z_{min}$ for $Z_{min} > 0$ (on a part of $\partial T$ of non-zero measure), so as to ensure the well-posedness of the classic wave propagation problems [7].

Naturally, we ask for the well-posedness of the local problem in every cell:
**Hypothesis 3.4.** In the rest of the article, we suppose the boundary impedance field $Z_{\partial T}$, the time-harmonic Friedrichs system and the mesh are chosen in such a way that Problem 3.1 is well-posed in every cell.

**Remark 3.5.** In the specific case of the wave propagation problems we study in this manuscript, such a result holds in the usual framework of Remark 3.3:

* for the Maxwell system, we refer to [7, Theorem 4.17] that ensures the well-posedness of the local problem with impedance boundary conditions.

* for the 2D Helmholtz equation, we refer to [50, Section 8.8] and [51] that ensures the well-posedness of the local problem with impedance boundary conditions.

### 3.1.2 Polynomial BC discretisation

After characterising the Trefftz function of $\mathbb{Y}^T$ thanks to their boundary impedance conditions, we propose to discretise them, thanks to piecewise polynomial functions. To do so, we begin with the introduction of a **micro-mesh** of the boundary $\partial T$.

**Definition 3.6 (Micro-mesh).** We consider a **micro-mesh** $\mathcal{T}^\mu$ (whose relative parameters will be referred with the superscript '$\mu$') of the boundary $\partial T$, made up of $N^\mu$ micro-faces in each direction of each face $F \in \mathcal{F}(T)$ (see Figure 2).

The micro-mesh of each face is then referred as $\mathcal{T}^\mu(F)$, such that:

$$ \mathcal{T}^\mu = \bigcup_{F \in \mathcal{F}(T)} \mathcal{T}^\mu(F). $$

In particular, for any micro-face $f \in \mathcal{T}^\mu$, $w : \partial T \to \mathbb{C}$ and $\mathbf{w} : \partial T \to \mathbb{C}^r$, we denote their restrictions to $f$ as

$$ w^f := w|_f \quad \text{and} \quad \mathbf{w}^f := \mathbf{w}|_f. $$

<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the micro-mesh for $N^\mu = 3$ in 2D (3 segments per face, identified with red dots) and $N^\mu = 2$ in 3D (2 micro-faces per direction of each face, identified with red dotted lines).</figcaption></figure>

Then, we introduce the function spaces on which the discretisation of the boundary condition relies, which are made up of discontinuous and continuous piecewise polynomial functions on each face.

**Definition 3.7.** For a cell $T \in \mathcal{T}_h$, one of its faces $F \in \mathcal{F}(T)$ and the polynomial degree $k^\mu \in \mathbb{N}$, we define the space of piecewise polynomial functions on its micro-faces of degree at most $k^\mu$ in each tangential direction as:

$$ \mathbb{Q}_{k^\mu} (\mathcal{T}^\mu(F)) := \{ w \in L^2(F) : \forall f \in \mathcal{T}^\mu(F), w^f \in \mathbb{Q}_{k^\mu}(f) \}. $$

The continuity of the functions between the micro-faces can also be imposed by considering

$$ \mathbb{Q}_{k^\mu}^0 (\mathcal{T}^\mu(F)) := \mathbb{Q}_{\mathbf{k}^\mu} (\mathcal{T}^\mu(F)) \cap \mathcal{C}^0(F). $$

Finally, considering all the faces $F \in \mathcal{F}(T)$, we then define:

$$ \mathbb{Q}_{k^\mu} (\mathcal{T}^\mu) := \prod_{F \in \mathcal{F}(T)} \mathbb{Q}_{k^\mu} (\mathcal{T}^\mu(F)) \quad \text{and} \quad \mathbb{Q}_{k^\mu}^0 (\mathcal{T}^\mu) := \prod_{F \in \mathcal{F}(T)} \mathbb{Q}_{k^\mu}^0 (\mathcal{T}^\mu(F)), \quad (28) $$

whose examples are presented on Figure 3 in 2D and 3D configurations.
<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the piecewise polynomial boundary conditions for $N^\mu = 3$ and $k^\mu = 2$ in 2D (quadratic in each micro-face), and $N^\mu = 2$ and $k^\mu = 1$ in 3D ($\mathbb{Q}_1$ in each micro-face) in the continuous case.</figcaption></figure>

**Remark 3.8.** Technically, we chose to discretise the polynomial space $\mathbb{Q}_{k^\mu}(f)$, on each micro-face $f \in \mathcal{T}_T^\mu(F)$, by relying on $(k^\mu + 1)^{d-1}$ tensorised Gauss-Lobatto nodes $(k^\mu + 1$ in each direction) and the associated Lagrange polynomials.

**Remark 3.9.** One will observe that even if it may sound reasonable to consider continuity of the boundary conditions on a single face, Definition 3.7 does not enforce continuity between the faces. Indeed, as the unit outward normal vector $\mathbf{n}_T$ is discontinuous between the different faces $F$ of $\partial T$, it is the same for $\mathcal{F}_{\mathbf{n}_T}$ and its SVD decomposition. Thus, only continuity along each face can be expected from the local boundary conditions.

**Remark 3.10.** To simplify the notations, we only considered that all the cells and faces were using the same parameters $N^\mu$ and $k^\mu$. Nevertheless, one may consider different numbers $N_{T,F}^\mu$ of micro-faces for each face $F \in \mathcal{F}(T)$ of $T \in \mathcal{T}_h$, in addition to polynomial orders $k_f^\mu$ relative to each micro-face $f \in \mathcal{T}_T^\mu$. This may allow specific and local control of the basis functions (in terms of their number and regularity for example), and then p-refinement.

### 3.1.3 Definition of the new Trefftz functions

Finally, we introduce our methodology of Trefftz functions discretisation through a polynomial approximation of local impedance boundary conditions.

For each cell, we consider a mesh of its boundary and piecewise polynomial fields on it. The $r$ components of the local impedance boundary condition are then searched in this polynomial boundary space and the associated solutions of the local Friedrichs system are taken as local basis functions for this cell: we then denote these Trefftz functions as **Local Polynomial Impedance Trefftz** ones (LPIT).

**Definition 3.11 (LPIT functions).** Let's consider a cell $T \in \mathcal{T}_h$, a micro-mesh $\mathcal{T}_T^\mu$ of its boundary $\partial T$ and a polynomial degree $k^\mu \in \mathbb{N}$. Then, let $(\mathbf{g}_j^T)_{j \in [1, N_{\mathrm{dof}}^T]}$ be a basis of $[\mathbb{Q}_{k^\mu}(\mathcal{T}_T^\mu)]^r$: each one of them has piecewise polynomial components defined on $\partial T$, which are interpreted as a local boundary condition.

Finally, referring to Problem 3.1 for the definition of the local solver, we denote the associated solutions in the cell $T$ as:

$$ \forall j \in [1, N_{\mathrm{dof}}^T], \quad \phi_j^T = S_T(\mathbf{g}_j^T). $$

We can finally define a discrete space of local solutions $\mathbb{Y}_h^T \subset \mathbb{Y}^T$ and associated outgoing traces $\mathbb{X}_h^T \subset \mathbb{X}^T$ as

$$ \mathbb{Y}_h^T := \operatorname{span}_{j \in [1, N_{\mathrm{dof}}^T]} \{\phi_j^T\} = \{\mathbf{y}^T \in [\mathbb{L}^2(T; \mathbb{C})]^m : \exists \mathbf{g}^T \in [\mathbb{Q}_{k^\mu}(\mathcal{T}_T^\mu)]^r, \mathbf{y}^T = S_T(\mathbf{g}^T)\}, $$

and

$$ \mathbb{X}_h^T := \gamma_{\mathrm{out}}^T (\mathbb{Y}_h^T) = \{\mathbf{x}^T \in [\mathbb{L}^2(\partial T; \mathbb{C})]^r : \exists \mathbf{y}^T \in \mathbb{Y}_h^T, \mathbf{x}^T = \gamma_{\mathrm{out}}^T(\mathbf{y}^T)\}. $$
Of course, the same definition stands by replacing $\mathbb{Q}_{k^\mu} (\mathcal{T}_T^\mu)$ by $\mathbb{Q}_{k^\mu}^0 (\mathcal{T}_T^\mu)$, so as to consider local boundary conditions which are continuous on each face.

We can naturally extend it on the whole mesh by defining

$$ \mathbb{Y}_h := \prod_{T \in \mathcal{T}_h} \mathbb{Y}_h^T \quad \text{and} \quad \mathbb{X}_h := \prod_{T \in \mathcal{T}_h} \mathbb{X}_h^T. $$

**Remark 3.12.** This methodology does not allow to choose explicitly the number of basis functions which are taken in each cell. Indeed, according to the fact the local polynomial BCs are chosen in $\mathbb{Q}_{k^\mu}$ or $\mathbb{Q}_{k^\mu}^0$, the number of local basis functions for a Cartesian cell is given by:

$$ N_{dof}^T = \begin{cases} 2rd(N^\mu(k^\mu + 1))^{d-1} & \text{if the BCs are \textbf{discontinuous} between micro-faces,} \\ 2rd(N^\mu k^\mu + 1)^{d-1} & \text{if the BCs are \textbf{continuous} between micro-faces,} \end{cases} $$

as each Cartesian cell has 2d faces of dimension $d-1$ and the local BCs have $r$ components.

Choosing continuous BCs then allows to reduce the number of basis functions, while it seems natural, for regular enough exact solutions, that the quantity $\mathcal{U}^*\mathbf{y}_1^T - Z_{\partial T} \mathcal{V}^*\mathbf{y}_2^T$ is continuous on each face. This configuration will then be favoured and systematically used in the rest of the article.

**Remark 3.13.** No analytic expressions of such basis functions exist, and one would have to rely on a local solver to compute approximations. This approximation process and its analysis will be detailed in Subsection 3.3 and Section 4. In practice, these LPIT functions are evaluated thanks to a strongly refined Discontinuous Galerkin solver.

In a nutshell, a regular-enough solution of the Friedrichs system in a cell $T$ can be characterised by its local impedance boundary conditions (due to the well-posedness of the local Problem 3.1, thanks to Hypothesis 3.4). Thus, we propose to introduce a piecewise polynomial approximation of this field on any face of $\partial T$ (see Figure 3). A basis of such polynomial boundary fields can then naturally be defined, and their associated solutions in $T$ are taken as basis functions of the discrete Trefftz space $\mathbb{Y}_h$.

## 3.2 Trefftz error estimate

Another justification of this definition of Trefftz basis functions leans on the following classic a priori error estimate (which is adapted from [15, 21, 44]): letting $\mathbf{y} \in \mathbb{Y}$ and $\mathbf{y}_h \in \mathbb{Y}_h$ denote the solutions of the continuous and discrete Trefftz Problem 2.17, one has

$$ \|\mathbf{y} - \mathbf{y}_h\|_{DG} \lesssim \inf_{\mathbf{x}'_h \in \mathbb{X}_h} \|\mathcal{U}\mathbf{x} - \mathcal{U}\mathbf{x}'_h\|_{\mathbb{X}} \quad \text{with} \quad \mathbf{x} = \gamma_{\text{out}}(\mathbf{y}), \qquad (29) $$

where the DG-norm is defined for $\mathbf{w} \in \mathbb{Y}$ as

$$ \|\mathbf{w}\|_{DG}^2 = \Re[a(\mathbf{w}, \mathbf{w})] = |\mathbf{w}|_{\text{int}}^2 + |\mathbf{w}|_{\text{ext}}^2, $$

with

$$ |\mathbf{w}|_{\text{int}}^2 = \sum_{F \in \mathcal{F}_{\text{int}}} \frac{1}{Z_T + Z_K} \|\mathcal{U}^*(\mathbf{w}_1^T - \mathbf{w}_1^K)\|_{0,F}^2 + \frac{Z_T Z_K}{Z_T + Z_K} \|\mathcal{V}^*(\mathbf{w}_2^T - \mathbf{w}_2^K)\|_{0,F}^2, $$

and

$$ |\mathbf{w}|_{\text{ext}}^2 = \sum_{F \in \mathcal{F}_{\text{ext}}} \frac{1}{Z_T + Z_{\partial \Omega}} \|\mathcal{U}^* \mathbf{w}_1^T\|_{0,F}^2 + \frac{Z_T Z_{\partial \Omega}}{Z_T + Z_{\partial \Omega}} \|\mathcal{V}^* \mathbf{w}_2^T\|_{0,F}^2, $$

while $\mathcal{F}_{\text{int}}$ and $\mathcal{F}_{\text{ext}}$ respectively stand for the sets of internal and boundary faces of the mesh.

Then, as the error of the Trefftz scheme is controlled by the best approximation of the incoming trace on the boundary of each mesh cell, it seems reasonable to approximate it (or a similar quantity, when $Z_{\partial T} \neq Z_T$, see Remark 3.2) with polynomials.

This Subsection then aims at closing this a priori error estimate in the case of these specific Trefftz functions, for which the best-approximation has to be characterised.

To complete the estimate (29), we have to bound the error of the best-approximation of the local impedance BC in terms of a mesh-independent norm of the exact solution $\mathbf{y}$. In particular, all the following proofs are realised with constant real local impedance fields $Z_{\partial T}$ on $\partial T$ (which can be distinct for different cells):
**Hypothesis 3.14.** In this Subsection, we suppose that, in each cell $T \in \mathcal{T}_h$ of the mesh, the impedance field is constant:

$$ \forall T \in \mathcal{T}_h, \exists \tilde{Z}_{\partial T} \in \mathbb{R}_+^*, \forall \mathbf{x} \in \partial T, Z_{\partial T}(\mathbf{x}) = \tilde{Z}_{\partial T}, $$

and uniformly bounded over the mesh:

$$ \exists \tilde{Z} \in \mathbb{R}_+^*, \forall T \in \mathcal{T}_h, |Z_{\partial T}| \le \tilde{Z}. $$

To do so, we first rely on the Bramble-Hilbert Lemma, as these Trefftz functions lean on a polynomial approximation of the boundary fields.

**Proposition 3.15.** Let's consider, in a cell $T \in \mathcal{T}_h$ of diameter $h$, a constant local boundary impedance $Z_{\partial T} \in \mathbb{R}_+^*$ on $\partial T$. Then, for a regular solution $\mathbf{y}^T \in [\mathbb{H}^{k^\mu+2}(T)]^m$ and its outgoing trace $\mathbf{x}^T = \gamma_{\text{out}}^T(\mathbf{y}) \in \mathbb{X}^T$, there exists $\tilde{\mathbf{x}}_h^T \in \mathbb{X}_h^T$ such that

$$ \|\mathcal{U}^T \mathbf{x}^T - \mathcal{U}^T \tilde{\mathbf{x}}_h^T\|_{\mathbb{X}^T} \lesssim \left(\frac{h}{N^\mu}\right)^{k^\mu+1} \sum_{F \in \mathcal{F}(T)} |\mathcal{U}^* \mathbf{y}_1^T|_{k^\mu+1,F} + |\mathcal{V}^* \mathbf{y}_2^T|_{k^\mu+1,F}. $$

*Proof.* We consider a cell $T \in \mathcal{T}_h$ and the local BC $\mathbf{g}^T := \mathcal{U}^* \mathbf{y}_1^T - Z_{\partial T} \mathcal{V}^* \mathbf{y}_2^T$. Then, we choose a face $F \in \mathcal{F}(T)$ and one of its micro-faces $f \in \mathcal{T}_h^\mu(F)$ (whose dimension is bounded by $h/N^\mu$). In particular, due to the fact that $\mathbf{y}^T \in [\mathbb{H}^{k^\mu+2}(T)]^m$, one has $\mathbf{y}^f \in [\mathbb{H}^{k^\mu+3/2}(f)]^m$: the constant character of $\mathcal{U}^*$ and $\mathcal{V}^*$ on $f$ then ensures that $\mathbf{g}^f := \mathbf{g}^T|_f \in [\mathbb{H}^{k^\mu+3/2}(f)]^r$.

(i) Then, the Bramble-Hilbert Lemma [52] ensures there exists a polynomial approximation $\tilde{\mathbf{g}}_h^f \in [\mathbb{Q}_{k^\mu}(f)]^r$ of the local BC $\mathbf{g}^f$ on $f$, such that

$$ \|\mathbf{g}^f - \tilde{\mathbf{g}}_h^f\|_{0,f} \lesssim \left(\frac{h}{N^\mu}\right)^{k^\mu+1} |\mathbf{g}^f|_{k^\mu+1,f}. \quad (30) $$

It is then possible to define the boundary field $\tilde{\mathbf{g}}_h^T \in [\mathbb{Q}_{k^\mu}(\mathcal{T}_T^\mu)]^r$ on $\partial T$ by joining the restrictions on each one of its micro-faces:

$$ \forall f \in \mathcal{T}_T^\mu, \quad \tilde{\mathbf{g}}_h^T|_f := \tilde{\mathbf{g}}_h^f, $$

which satisfies

$$ \|\mathbf{g}^T - \tilde{\mathbf{g}}_h^T\|_{0,\partial T}^2 = \sum_{f \in \mathcal{T}_T^\mu} \|\mathbf{g}^f - \tilde{\mathbf{g}}_h^f\|_{0,f}^2 \lesssim \left(\frac{h}{N^\mu}\right)^{2(k^\mu+1)} \sum_{f \in \mathcal{T}_T^\mu} |\mathbf{g}^f|_{k^\mu+1,f}^2 = \left(\frac{h}{N^\mu}\right)^{2(k^\mu+1)} \sum_{F \in \mathcal{F}(T)} |\mathbf{g}^T|_{k^\mu+1,F}^2, $$

where we used the fact that $\mathbf{g}^T|_F \in [\mathbb{H}^{k^\mu+1}(F; \mathbb{C})]^r$, in addition to $\partial T = \bigcup_{f \in \mathcal{T}_T^\mu} f$ and (30).

Naturally, recalling the fact that $\mathbf{g}^T = \mathcal{U}^* \mathbf{y}_1^T - Z_{\partial T} \mathcal{V}^* \mathbf{y}_2^T$ and the uniformly bounded character of $Z_{\partial T}$ (see Hypothesis 3.14), one gets

$$ |\mathbf{g}^T|_{k^\mu+1,F}^2 \lesssim |\mathcal{U}^* \mathbf{y}_1^T|_{k^\mu+1,F}^2 + |\mathcal{V}^* \mathbf{y}_2^T|_{k^\mu+1,F}^2, $$

and then

$$ \|\mathbf{g}^T - \tilde{\mathbf{g}}_h^T\|_{0,F}^2 \lesssim \left(\frac{h}{N^\mu}\right)^{2(k^\mu+1)} \left( |\mathcal{U}^* \mathbf{y}_1^T|_{k^\mu+1,F}^2 + |\mathcal{V}^* \mathbf{y}_2^T|_{k^\mu+1,F}^2 \right). $$

(ii) Moreover, by Definition 3.11 of the discrete Trefftz spaces, there exist unique $\tilde{\mathbf{y}}_h^T \in \mathbb{Y}_h^T$ and $\tilde{\mathbf{x}}_h^T = \gamma_{\text{out}}^T(\tilde{\mathbf{y}}_h^T) \in \mathbb{X}_h^T$ such that

$$ \mathcal{U}^* \tilde{\mathbf{y}}_1^T - Z_{\partial T} \mathcal{V}^* \tilde{\mathbf{y}}_2^T = \tilde{\mathbf{g}}_h^T \quad \text{on } \partial T, $$

and by introducing the difference to the exact solution as $\tilde{\mathbf{\epsilon}}^T = \mathbf{y}^T - \tilde{\mathbf{y}}^T \in \mathbb{Y}^T$, one has

$$ \mathbf{g}^T - \tilde{\mathbf{g}}_h^T = \mathcal{U}^* \tilde{\mathbf{\epsilon}}_1^T - Z_{\partial T} \mathcal{V}^* \tilde{\mathbf{\epsilon}}_2^T \quad \text{and} \quad \mathcal{U}\mathbf{x}^T - \mathcal{U}\tilde{\mathbf{x}}_h^T = \mathcal{U}^* \tilde{\mathbf{\epsilon}}_1^T - Z_T \mathcal{V}^* \tilde{\mathbf{\epsilon}}_2^T. \quad (31) $$
Nevertheless, leaning on the reciprocity formula (17), the error in terms of incoming traces can be bounded thanks to the local BC. Indeed, one gets

$$
\begin{align*}
\|\mathfrak{U}\mathbf{x}^T - \mathfrak{U}\tilde{\mathbf{x}}_h^T\|_{\mathbf{x}^T}^2 &= Y_T \langle \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T - Z_T \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T, \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T - Z_T \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T \rangle_{0,\partial T} \\
&= Y_T \langle \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T, \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T \rangle_{0,\partial T} + Z_T \langle \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T, \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T \rangle_{0,\partial T} - \underbrace{\left\langle \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T, \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T \right\rangle_{0,\partial T} + \left\langle \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T, \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T \right\rangle_{0,\partial T}}_{r_T(\tilde{\mathbf{\epsilon}}^T, \tilde{\mathbf{\epsilon}}^T)=0},
\end{align*}
$$

as $\tilde{\mathbf{\epsilon}}^T \in \mathbb{Y}^T$. Yet, due to the fixed heterogeneity distribution (which is independent of the mesh), the local parameters are uniformly bounded, and the two first positive terms can be bounded by replacing $Z_T$ by $Z_{\partial T}$ (which is also constant over $\partial T$) so as to make $\mathbf{g}^T - \tilde{\mathbf{g}}_h^T$ appear. This leads to

$$
\| \mathfrak{U} \mathbf{x}^T - \mathfrak{U} \tilde{\mathbf{x}}_h^T \|_{\mathbf{x}^T}^2 \lesssim \langle \mathcal{U}^* \tilde{\mathbf{\epsilon}}_1^T, \mathcal{U}^* \tilde{\mathbf{\epsilon}}_1^T \rangle_{0, \partial T} + Z_{\partial T}^2 \langle \mathcal{V}^* \tilde{\mathbf{\epsilon}}_2^T, \mathcal{V}^* \tilde{\mathbf{\epsilon}}_2^T \rangle_{0, \partial T},
$$

and subtracting $Z_{\partial T} r_T(\tilde{\mathbf{\epsilon}}^T, \tilde{\mathbf{\epsilon}}^T) = 0$ to this quantity leads to

$$
\begin{align*}
\|\mathfrak{U}\mathbf{x}^T - \mathfrak{U}\tilde{\mathbf{x}}_h^T\|_{\mathbf{x}^T}^2 &\lesssim \langle \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T, \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T \rangle_{0,\partial T} + Z_{\partial T}^2 \langle \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T, \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T \rangle_{0,\partial T} - Z_{\partial T} r_T(\tilde{\mathbf{\epsilon}}^T, \tilde{\mathbf{\epsilon}}^T) \\
&= \langle \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T - Z_{\partial T} \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T, \mathcal{U}^*\tilde{\mathbf{\epsilon}}_1^T - Z_{\partial T} \mathcal{V}^*\tilde{\mathbf{\epsilon}}_2^T \rangle_{0,\partial T}.
\end{align*}
$$

Finally, by using (31), one gets

$$
\| \mathfrak{U} \mathbf{x}^T - \mathfrak{U} \tilde{\mathbf{x}}_h^T \|_{\mathbf{x}^T}^2 \lesssim \| \mathcal{U}^* \tilde{\mathbf{\epsilon}}_1^T - Z_{\partial T} \mathcal{V}^* \tilde{\mathbf{\epsilon}}_2^T \|_{0, \partial T}^2 = \| \mathbf{g}^T - \tilde{\mathbf{g}}_h^T \|_{0, \partial T}^2.
$$

(iii) By gathering the results of (i) and (ii), we can conclude.

Yet, the quantity $|\mathcal{U}^*\mathbf{y}_1^T|_{k^\mu+1,F} + |\mathcal{V}^*\mathbf{y}_2^T|_{k^\mu+1,F}$ is still mesh-dependent: we come back to a volumic estimate to get rid of this dependency.

**Proposition 3.16.** Let's consider a cell $T \in \mathcal{T}_h$, $\mathbf{y}^T \in [H^{k^\mu+2}(T)]^m$ and its associated outgoing trace $\mathbf{x}^T = \gamma_{\text{out}}^T(\mathbf{y}^T)$. Then one has

$$
\sum_{F \in \mathcal{F}(T)} |\mathcal{U}^* \mathbf{y}_1^T|_{k^\mu+1,F} + |\mathcal{V}^* \mathbf{y}_2^T|_{k^\mu+1,F} \lesssim h^{-1/2} |\mathbf{y}^T|_{k^\mu+1,T} + h^{1/2} |\mathbf{y}^T|_{k^\mu+2,T}.
$$

*Proof.* We first get rid of the operators $\mathcal{U}^*$ and $\mathcal{V}^*$. To do so, we consider a cell $T \in \mathcal{T}_h$ and one of its faces $F \in \mathcal{F}(T)$ (on which both $\mathcal{U}$ and $\mathcal{V}$ are constant). We recall the first compact SVD (8) of $\mathcal{B}_{n_T}$ leads to $\mathcal{U} = \hat{\mathcal{U}}\mathcal{\Sigma}^{1/2}$, where Remark 2.4 ensures these operators are bounded in an independent way of $\mathbf{n}_T$: $\|\hat{\mathcal{U}}\| = 1$ and $\|\mathcal{\Sigma}\| \le \sum_{j=1}^d \|\mathcal{B}_j\|$. Thus, there exists a constant $C > 0$ (independent of $\mathbf{n}_T$) such that

$$
\| \mathcal{U}^* \mathbf{y}_1^T \|_{0,F} \le C \| \mathbf{y}_1^T \|_{0,F}.
$$

As this manipulation can be realised on each face (on which $\mathcal{U}^*$ is constant) and also for $\mathcal{V}^*\mathbf{y}_2^T$, one obtains

$$
\sum_{F \in \mathcal{F}(T)} |\mathcal{U}^* \mathbf{y}_1^T|_{k^\mu+1,F}^2 + |\mathcal{V}^* \mathbf{y}_2^T|_{k^\mu+1,F}^2 \lesssim \sum_{F \in \mathcal{F}(T)} |\mathbf{y}_1^T|_{k^\mu+1,F}^2 + |\mathbf{y}_2^T|_{k^\mu+1,F}^2 = \sum_{F \in \mathcal{F}(T)} |\mathbf{y}^T|_{k^\mu+1,F}^2.
$$

Yet, coming back to the definition of the semi-norm $|\cdot|_{k^\mu+1,F}$, it only calls for tangential derivatives (that we denote $\mathcal{D}_F^{\tilde{\alpha}}$ for $\tilde{\alpha} \in \mathbb{N}^{d-1}$) which are included in the 'volumic' ones (that we denote $\mathcal{D}^{\alpha}$ for $\alpha \in \mathbb{N}^d$, where the derivative in the normal direction to $F$ is added).

In particular, considering a Cartesian cell and one of its faces $F \in \mathcal{F}(T)$ which is orthogonal to $\mathbf{d}_1$, any tangential derivative $\mathcal{D}_F^{\tilde{\alpha}}$ (with $\tilde{\alpha} \in \mathbb{N}^{d-1}$) can be seen as a 'volumic' one for $\alpha = (0, \tilde{\alpha}_1, \dots, \tilde{\alpha}_{d-1}) \in \mathbb{N}^d$: this simply corresponds to the subcase where there is no differentiation along the normal direction to $F$.
Thus, one has

$$
\sum_{F \in \mathcal{F}(T)} |\mathbf{y}^T|_{k^\mu+1,F}^2 = \sum_{F \in \mathcal{F}(T)} \sum_{\substack{\tilde{\alpha} \in \mathbb{N}^{d-1} \\ |\tilde{\alpha}| = k^\mu+1}} \|\mathcal{D}_F^{\tilde{\alpha}} \mathbf{y}^T\|_{0,F}^2 \lesssim \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha| = k^\mu+1}} \sum_{F \in \mathcal{F}(T)} \|\mathcal{D}^{\alpha} \mathbf{y}^T\|_{0,F}^2 = \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha| = k^\mu+1}} \|\mathcal{D}^{\alpha} \mathbf{y}^T\|_{0,\partial T}^2,
$$
where we used the fact that $\mathcal{D}^{\alpha} \mathbf{y}^T \in [\mathrm{H}^1(T)]^m$, and then its trace satisfies $\mathcal{D}^{\alpha} \mathbf{y}^T|_{\partial T} \in [\mathrm{L}^2(\partial T)]^m$.

Then, by using a trace inequality [53, Theorem 1.6.6], we obtain

$$ \|\mathcal{D}^{\alpha} \mathbf{y}^T\|_{0,\partial T}^2 \lesssim h^{-1} \|\mathcal{D}^{\alpha} \mathbf{y}^T\|_{0,T}^2 + h \|\nabla \mathcal{D}^{\alpha} \mathbf{y}^T\|_{0,T}^2, $$

which implies that

$$ \sum_{F \in \mathcal{F}(T)} |\mathbf{y}^T|_{k^{\mu}+1,F}^2 \lesssim h^{-1} |\mathbf{y}^T|_{k^{\mu}+1,T}^2 + h |\mathbf{y}^T|_{k^{\mu}+2,T}^2. $$

□

Finally, gathering these different estimates, we obtain an a priori error estimate for this Trefftz formulation in DG-norm:

**Theorem 3.17.** Let $\mathbf{y} \in \mathbb{Y}$ and $\mathbf{y}_h \in \mathbb{Y}_h$ denote the solutions of the continuous and discrete Trefftz Problem 2.17 obtained for LPIT functions associated to a constant local boundary field of $Z_{\partial T} \in \mathbb{R}_*^+$ in each cell. Moreover, we suppose an additional regularity on the solution, as $\mathbf{y} \in [\mathrm{H}^{k^{\mu}+2}(\Omega; \mathbb{C})]^m$. Then, one has the a priori error estimate

$$ \|\mathbf{y} - \mathbf{y}_h\|_{DG} \lesssim \left(\frac{h}{N^{\mu}}\right)^{k^{\mu}+1} \left(h^{-1/2} |\mathbf{y}|_{k^{\mu}+1,\Omega} + h^{1/2} |\mathbf{y}|_{k^{\mu}+2,\Omega}\right). \quad (32) $$

*Proof.* This estimate naturally results from (29), in addition to Propositions 3.15 and 3.16, as the regular character of $\mathbf{y}$ in $\Omega$ ensures that

$$ \sum_{T \in \mathcal{T}_h} |\mathbf{y}^T|_{k^{\mu}+1,T}^2 = |\mathbf{y}|_{k^{\mu}+1,\Omega} \quad \text{and} \quad \sum_{T \in \mathcal{T}_h} |\mathbf{y}^T|_{k^{\mu}+2,T}^2 = |\mathbf{y}|_{k^{\mu}+2,\Omega}. $$

□

**Remark 3.18.** In this finite dimensional setting for which the DG sesquilinear form is a norm (see [35, Section 3.3]), the error estimate (32) also allows to ensure the well-posedness of the numerical method. Indeed, as existence and uniqueness are equivalent in this finite-dimensional case, we can prove solution uniqueness by considering the boundary condition $\mathbf{g} = \mathbf{0}$ on $\partial \Omega$ and its associated solution $\mathbf{y} = \mathbf{0}$. Then, by denoting the Trefftz solution as $\mathbf{y}_h \in \mathbb{Y}_h$, one gets

$$ \|\mathbf{y}_h\|_{DG} = \|\mathbf{y} - \mathbf{y}_h\|_{DG} \lesssim \left(\frac{h}{N^{\mu}}\right)^{k^{\mu}+1} \left(h^{-1/2} |\mathbf{y}|_{k^{\mu}+1,\Omega} + h^{1/2} |\mathbf{y}|_{k^{\mu}+2,\Omega}\right) = 0, $$

which implies that $\mathbf{y}_h = \mathbf{0}$, as DG is a norm on $\mathbb{Y}_h \subset \mathbb{Y}$.

## 3.3 Introduction of the quasi-Trefftz basis functions

Nevertheless, considering such polynomial boundary conditions, one cannot derive an analytic expression of the associated solutions: the exact local solver 3.1 is then not explicitly available, which compels us to introduce an approximate one. The associated functions will then no longer be actual solutions (as the LPIT functions), but only quasi-solutions: this leads to a numerical method referred as a quasi-Trefftz one.

### 3.3.1 Approximate local solver

Facing the practical unavailability of the exact local solver, we resort to an approximate one (which can be chosen as a classic numerical method for example), so as to supply approximations of the LPI Trefftz functions introduced in Section 3.1.

**Problem 3.19.** We consider the approximate local solution operator $\tilde{S}_T$ which associates, to any $\mathbf{g}^T \in [\mathrm{L}^2(\partial T)]^r$, the approximate solution $\tilde{\mathbf{y}}^T = \tilde{S}_T(\mathbf{g}^T) \in [\mathrm{L}^2(T; \mathbb{C})]^m$ to the original time-harmonic Friedrichs Problem 3.1. Naturally, being given a set of local boundary conditions $\mathbf{g} = (\mathbf{g}^T)_{T \in \mathcal{T}_h} \in [\mathrm{L}^2(\partial \mathcal{T}_h)]^r$, we denote the set of associated approximate local solutions as $\tilde{\mathbf{y}} = \tilde{S}(\mathbf{g}) \in [\mathrm{L}^2(\Omega; \mathbb{C})]^m$ such that

$$ \forall T \in \mathcal{T}_h, \quad \tilde{\mathbf{y}}^T = \tilde{S}_T(\mathbf{g}^T). $$
In particular, this means that we expect the approximate solution almost verifies the PDE and boundary conditions:

$$ \forall T \in \mathcal{T}_h, \quad \mathcal{L}_T \tilde{\mathbf{y}}^T \approx \mathbf{0} \quad \text{in } T \quad \text{and} \quad \gamma_{in}^{\partial T}(\tilde{\mathbf{y}}^T) \approx \mathbf{g}^T \quad \text{on } \partial T. $$

This highlights the 'quasi' character of this approach, while the expected meaning of the symbol '$\approx$' will be discussed in Section 4.

**Remark 3.20.** Due to the fact the whole solution $\mathbf{y}^T$ is required by the Trefftz formulation (see Problem 2.17), solving the first-order problem may be pertinent in this case. Indeed, considering the Maxwell's equations for example, solving the second-order problem on the electric field $\mathbf{e}$ usually leads to precision losses in the post-processing estimation of $\mathbf{h} = (i\kappa\mu)^{-1} \nabla \times \mathbf{e}$.

For the following numerical experiments, we then rely on an auxiliary solver based on the Discontinuous Galerkin method. Nevertheless, we highlight the fact that any other kind of numerical method may have been considered as the FEM [34] or integral equations [54] for example.

**Remark 3.21.** Naturally, if two cells share the same parameters (shape, local impedance field $Z_{\partial T}$ and dielectric parameters $(\epsilon_T, \mu_T)$), then it has to be inverted only once. In particular, this computation needs to be done a single time in any homogeneous part of a Cartesian mesh (no matter how large it is).

### 3.3.2 Discretisation parameters of the approximate local solver

The approximate solver $\tilde{S}_T$, being obtained thanks to a numerical method, is usually characterised by discretisation parameters in each cell $T \in \mathcal{T}_h$.

In particular, considering the framework of the DG method, it relies on two characteristics:

* a submesh $\mathcal{T}_T^n$ of the cell $T$, that we denote as **nano-mesh** (illustrated in Figure 4). All the characteristics relative to it will then be referred with the superscript '$n$'.

In the case of a Cartesian cell $T$, we also choose a Cartesian nano-mesh which is made up of $N^n$ nano-cells in each direction. Moreover, it is naturally defined in a conformal way with the micro-mesh $\mathcal{T}_T^\mu$ of $\partial T$: the separation between two micro-faces coincides with the one of two nano-cells. In other words, it means there exists a whole number of nano-cells per micro-face, $N^{n/\mu} \in \mathbb{N}$, such that

$$ N^n = N^{n/\mu} \times N^\mu. \qquad (33) $$

* a polynomial degree $k^n \in \mathbb{N}$ which characterises the approximation in each nano-cell.

**Remark 3.22.** We highlight the fact that the micro-mesh $\mathcal{T}_T^\mu$ and the nano-mesh $\mathcal{T}_T^n$ are absolutely independent (as long as the nano-mesh is conformal with the micro-mesh). Indeed, once the micro-mesh is fixed, the associated polynomial boundary conditions do not evolve: it is then possible to refine the nano-mesh so as to improve the approximation of the associated solution, which is independent of this nano-mesh.

On the contrary, the original approach proposed in [34] relies on a FEM local solver, whose nano-mesh coincides with the micro-mesh on the boundaries. In our case, a nano-mesh refinement will be possible, without any effect on the number of basis functions.

**Remark 3.23.** This whole procedure can be realised for general polygonal and polyhedral cells without any additional difficulty, as long as an adapted local solver is available for the cells geometry: this Cartesian framework only aims at facilitating the presentation.

### 3.3.3 Definition of the quasi-Trefftz functions

In the end, as the exact solver $S_T$ is not available in practice, we rely on the approximate one $\tilde{S}_T$ of Subsubsection 3.3.1 to define Local Polynomial Impedance Quasi-Trefftz (LPIQT) functions as approximations of the LPIT ones.

**Definition 3.24 (LPIQT functions).** Let's consider a cell $T \in \mathcal{T}_h$, a micro-mesh $\mathcal{T}_T^\mu$ of its boundary $\partial T$ and a polynomial degree $k^\mu \in \mathbb{N}$. Then, let $(\mathbf{g}_j^T)_{j \in [1, N_{dof}^T]}$ be a basis of $[\mathbb{Q}_{k^\mu}(\mathcal{T}_T^\mu)]^T$: each one of them has piecewise polynomial components defined on $\partial T$, which are interpreted as a local boundary condition. Thus, we denote the associated quasi-solutions in the cell $T$ as:

$$ \forall j \in [1, N_{dof}^T], \quad \tilde{\phi}_j^T = \tilde{S}_T(\mathbf{g}_j^T). $$
<figure><img src="image_5.png" /><figcaption>(i) 2D configuration with $\mathbf{N} = (3, 2)$, $N^\mu = 3$ and $N^n = 6$ (i.e. $N^{n/\mu} = 2$)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(ii) 3D configuration with $\mathbf{N} = (3, 1, 2)$, $N^\mu = 2$ and $N^n = 2$ (i.e. $N^{n/\mu} = 1$)</figcaption></figure>

Figure 4: Illustration of the different meshes in 2D and 3D: mesh $\mathcal{T}_h$ of the domain $\Omega$ (in black on the left, made up of $\mathbf{N} = (N_1, \dots, N_d)$ cells in each direction) and zoom on one cell $T \in \mathcal{T}_h$, with both the micro-mesh $\mathcal{T}_T^\mu$ of its boundary $\partial T$ (in red on the right, made up of $N^\mu$ micro-faces in each direction) and the nano-mesh $\mathcal{T}_T^n$ of $T$ (in blue on the right, made up of $N^n$ nano-cells in each direction).

We can finally define a discrete space of local quasi-solutions $\tilde{\mathbb{Y}}_h^T \subset \mathbb{Y}^T$ and associated outgoing traces $\tilde{\mathbf{X}}_h^T \subset \mathbf{X}^T$ as

$$ \tilde{\mathbb{Y}}_h^T := \operatorname{span}_{j \in [1, N_{\mathrm{dof}}^T]} \{\tilde{\phi}_j^T\} = \{\tilde{\mathbf{y}}^T \in [L^2(T; \mathbb{C})]^m : \exists \mathbf{g}^T \in [\mathbb{Q}_{k^\mu} (\mathcal{T}_T^\mu)]^r, \tilde{\mathbf{y}}^T = \tilde{S}_T(\mathbf{g}^T)\}, $$

and

$$ \tilde{\mathbf{X}}_h^T := \gamma_{\mathrm{out}}^T (\tilde{\mathbb{Y}}_h^T) = \{\tilde{\mathbf{x}}^T \in [L^2(\partial T; \mathbb{C})]^r : \exists \tilde{\mathbf{y}}^T \in \tilde{\mathbb{Y}}_h^T, \tilde{\mathbf{x}}^T = \gamma_{\mathrm{out}}^T(\tilde{\mathbf{y}}^T)\}. $$

Of course, the same definition stands by replacing $\mathbb{Q}_{k^\mu} (\mathcal{T}_T^\mu)$ by $\mathbb{Q}_{k^\mu}^0 (\mathcal{T}_T^\mu)$, so as to consider local boundary conditions which are continuous on each face.

We can naturally extend it on the whole mesh by defining

$$ \tilde{\mathbb{Y}}_h := \prod_{T \in \mathcal{T}_h} \tilde{\mathbb{Y}}_h^T \quad \text{and} \quad \tilde{\mathbf{X}}_h := \prod_{T \in \mathcal{T}_h} \tilde{\mathbf{X}}_h^T. $$

Once these discretisation spaces are defined, it is then natural to introduce the quasi-Trefftz formulation associated to them, as a counterpart of Problem 2.17:

**Problem 3.25** (Discrete quasi-Trefftz UWVF for BATFFS). Find $\tilde{\mathbf{x}}_h \in \tilde{\mathbf{X}}_h$ such that for all $\tilde{\mathbf{x}}' \in \tilde{\mathbf{X}}_h$, one has:

$$ a_{\mathrm{UWVF}}(\tilde{\mathbf{x}}_h, \tilde{\mathbf{x}}') = \ell_{\mathrm{UWVF}}(\tilde{\mathbf{x}}'), $$

where we refer to Problem 2.17 for the definitions of $a_{\mathrm{UWVF}}$ and $\ell_{\mathrm{UWVF}}$.

In terms of noteworthy properties, the ones of the original Trefftz method (recalled in Subsection 2.2.2) cannot be expected to be recovered, as they strongly relied on the **actual solution character** of the associated basis functions: this is no longer the case, as quasi-Trefftz ones are only **quasi-solutions**.

Nevertheless, it can be expected that, if the approximate local solver is refined enough, the LPIQT functions will be
really close to the LPIT ones, allowing to recover the properties of the Trefftz method. This feeling can be proved
to be right: for a sufficiently refined approximate local solver, the quasi-Trefftz method satisfies a weak-coercivity
property and the natural mass preconditioner leads to a contracting linear system.

**Remark 3.26.** A parallel can be drawn between the approach we propose and the CHDG method (where the letter 'C' refers to 'characteristic variable' and 'HDG' corresponds to Hybridisable Discontinuous Galerkin), where the local problem in each subdomain is condensed on its boundary by relying on a characteristic variable, which can be interpreted as a transmission one [55, 56]: indeed, the Local Polynomial Impedance Trefftz functions are also associated to a polynomial representation of a Fourier-Robin kind of trace.

Nevertheless, in our case, the associated formulation is considered and studied in a Trefftz point of view, while HDG usually enforces a direct relation between the local scattering DG solver and the polynomial degree of the hybrid variable. In the case of our quasi-Trefftz approach, both the local traces and approximation of the associated basis functions are decoupled, allowing to adapt this latter to get the best performances.

# 4 Numerical study and calibration of the quasi-Trefftz approach

The previous Sections allowed to introduce a quasi-Trefftz formulation adapted to BATFFS, for which the basis
functions are approximate local solutions. Being obtained resorting to a numerical method in each cell, we numer-
ically discuss two fundamental questions in this last Section: is this quasi-Trefftz approach robust and accurate
to deal with the solution of such problems, and how does the local solver have to be calibrated (in terms of its
discretisation parameters, see Subsection 3.3.2, with respect to the ones of the BCs, see Subsection 3.1.2)?
In particular, we propose a numerical illustration of the abilities of the quasi-Trefftz approach in terms of basis
independence, *h* and *p*-convergences, in addition to its adaptability to iterative solution.

## 4.1 Validation of the basis independence

First of all, we recall the Local Polynomial Impedance (Quasi-)Trefftz functions have initially been introduced to
tackle the numerical dependence limitations induced by plane waves (see Subsection 2.3). Thus, considering a
cell $T \in \mathcal{T}_h$, we take back the condition number of the mass matrix $D^T$ (see (26)) which is supposed to serve as
preconditioner, allowing to assess the dependence of local basis functions.

More precisely, we consider a unique Cartesian cell $T$ of size 1 wavelength in each direction and compute the
local mass matrix $D^T$ for the LPIQT basis functions associated to different BC polynomial degrees $k^\mu \in \{1, 2, 3\}$. In
particular, the influence of the local solver is discussed in terms of its polynomial approximation degree $k^n \in [0, 5]$ and the refinement of the nano-mesh $\mathcal{T}_h^n$ (through the number of nano-cells per micro-face, $N^{n/\mu}$).

Then, Figure 5 presents the evolution of the condition number of the mass matrix $D^T$ with respect to the number $N_{\text{dof}}^T$ of basis functions per cell (which is increased thanks to the number of micro-faces $N^\mu$), for the 2D Helmholtz equation (referred as 2DH) and 3D Maxwell system (referred as 3DM) respectively. In every configuration, one can notice that the condition number is bounded, in opposition to the behaviour of plane waves observed in Subsection 2.3. Moreover, there seems to exist a limit curve towards which the LPIQT functions converge, when the polynomial degree $k^n$ increases or the nano-mesh is refined (meaning the local solver becomes more accurate): this obviously corresponds to the case of an 'exact' local solver, and then to the LPI Trefftz functions.

Nevertheless, considering the coarsest nano-mesh (for $N^{n/\mu} = 1$), the configurations associated to $k^n < k^\mu$ do not
appear, as the mass matrix is found to be not invertible. Thus, a minimum accuracy of the local solver (with
respect to the polynomial degree of the BCs) must be ensured so as to obtain decent enough approximations of
actual solutions. In particular, choosing a slightly more refined FR method, with $k^n = k^\mu + 1$, already delivers
conditioning evolutions close to the ones of the LPIT functions.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Figure 5: Evolution of the condition number of the mass matrix for LPIQT functions in 2D and 3D, with respect to the number of basis functions per cell, $N_{dof}^T$, for different polynomial degrees $k^\mu$ of the BCs and nano-mesh refinements (through the number of nano-cells per micro-face, $N^{n/\mu}$). We consider a FR local solver of degree $k^n = 0$ (■+), $k^n = 1$ (■+), $k^n = 2$ (■+), $k^n = 3$ (■+), $k^n = 4$ (■+), and $k^n = 5$ (■+).</figcaption></figure>
## 4.2 Calibration of the local solver through convergence studies

The previous Subsection allowed to confirm the LPIQT functions show robust independence properties, under the condition of a sufficiently refined local solver (with respect to the discretisation parameters of the boundary conditions). We propose, in this Section, to complete this first assessment with the development of a calibration rule of the approximate local solver so as to get optimal abilities, thanks to numerical studies of the $h$ and $p$-convergences of the quasi-Trefftz method.

In opposition to the previous Subsection which only considered the LPIQT basis functions, without coupling with the Trefftz methodology, this Subsection considers the full quasi-Trefftz method: both quasi-solution character of the LPIQT functions and their interaction with the variational formulation are going to affect the results.

More precisely, this Subsection aims at studying the quasi-Trefftz method for different configurations of the local solver, focusing on two distinct regimes for given domain $\Omega$ and polynomial degree $k^\mu$ of the approximation on each micro-face:

* $h$-convergence, where the polynomial discretisation of the local boundary conditions (through the number $N^\mu$ of micro-faces in each direction) and the approximation of the associated Trefftz functions (through the local solver parameters $k^n$ and $N^{n/\mu}$) are fixed, while the mesh $\mathcal{T}_h$ is refined. More specifically, the number of basis functions in each cell (and their approximation) does not evolve, but the number of cells is increased.

* $p$-convergence, where the mesh $\mathcal{T}_h$ (in addition to the local solver parameters $k^n$ and $N^{n/\mu}$) is fixed, while the local boundary mesh is refined (meaning that the number of micro-faces in each direction, $N^\mu$, and then the number of local basis function, $N_{dof}^T$ (see Remark 3.12), are increased). In particular, as only the number $N^{n/\mu}$ of nano-cells per micro-face is fixed, the nano-mesh is refined in accordance to the micro-mesh (see (33)). This corresponds to a convergence in terms of the number of basis functions.

In particular, the numerical results we present in this Subsection are realised in 2D and 3D, for the Helmholtz and Maxwell's equations respectively, for $\kappa = 2\pi$ in a domain $\Omega = [0, 1]^d$ (the domain is then of size 1 wavelength in each direction). Moreover, we consider a uniform Cartesian mesh for simplicity of implementation and description (with the same number of cells, $N$, in each direction), in addition to a constant impedance boundary field in each cell $T \in \mathcal{T}_h$.

Hence, the error is computed in the $L^2(\Omega)$ norm (for $p - p_h$ in 2D and $y - y_h$ in 3D) with respect to an actual solution made up of a plane wave in both cases: this smooth and regular configuration allows to highlight the best convergence properties that may be expected from this method. Moreover, the linear system is symmetrically preconditioned thanks to the mass matrix $D$, and solved thanks to a restarted GMRES [57, 58, 59] with Krylov space of size $N_{Kry} = 20$ with a residual target of $10^{-8}$.

### 4.2.1 $h$-convergence

We begin with the mesh-convergence of the method: we fix the BC discretisation (through the number of micro-faces per direction to $N^\mu = 2$, and the polynomial degree $k^\mu \in \{1, 2\}$) and the parameters of the local solver (through the polynomial degree $k^n$ of the approximation in each nano-cell, and the number $N^{n/\mu}$ of nano-cells per micro-face for the nano-mesh refinement), while the number $N$ of cells of $\mathcal{T}_h$ is increased.

Then, Figure 6 presents the evolution of the relative $L^2(\Omega)$-error with respect to $N$, for the 2D Helmholtz and 3D Maxwell's equations respectively. In particular, in case some curves are not drawn, this means that the GMRES solver failed to converge (in the sense the residual stagnates, preventing it from reaching convergence after a maximum number of iterations), or the local solver was too expensive (in the 3D configuration).

Consistently with what could be observed in Subsection 4.1, a common limit convergence is reached when the local solver is accurate enough (thanks to a refined nano-mesh or a high-order approximation): it corresponds to what can be obtained with LPIT functions, as the solver is then 'perceived' as 'perfect'. In particular, when the mesh $\mathcal{T}_h$ becomes fine enough (i.e. for a large enough value of $N$), one can observe that quasi-Trefftz curves separate from the 'ideal' one of the LPIT functions. Then, a different convergence rate (smaller than the 'ideal' one of LPIT) is obtained, depending on the polynomial degree $k^n$.

Nevertheless, this common convergence seems to be reached only in the case $k^n > k^\mu$, and the larger $k^n$ is, the larger is the domain of accordance with it. Moreover, refining the nano-mesh (by increasing $N^{n/\mu}$) also leads to the same effect. In a nutshell, when the local solver is accurate enough, then the quasi-solutions are sufficiently close to
the actual ones, allowing to recover the convergence of the Trefftz method with LPIT functions (which is coherent with the recovery of the Trefftz properties for an accurate enough local solver).

<figure><img src="image_19.png" /><figcaption>(i) $k^\mu = 1$ and $N^{n/\mu} \in \{1, 2, 3\}$ for 2DH</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(ii) $k^\mu = 1$ and $N^{n/\mu} \in \{3, 4, 5\}$ for 2DH</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(iii) $k^\mu = 1$ and $N^{n/\mu} \in \{5, 6, 7\}$ for 2DH</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(iv) $k^\mu = 2$ and $N^{n/\mu} \in \{1, 2, 3\}$ for 2DH</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(v) $k^\mu = 2$ and $N^{n/\mu} \in \{3, 4, 5\}$ for 2DH</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(vi) $k^\mu = 2$ and $N^{n/\mu} \in \{5, 6, 7\}$ for 2DH</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(vii) $k^\mu = 1$ and $N^{n/\mu} = 1$ for 3DM</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(viii) $k^\mu = 1$ and $N^{n/\mu} \in \{1, 2, 3\}$ for 3DM</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(ix) $k^\mu = 1$ and $N^{n/\mu} \in \{3, 4, 5\}$ for 3DM</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(x) $k^\mu = 2$ and $N^{n/\mu} = 1$ for 3DM</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(xi) $k^\mu = 2$ and $N^{n/\mu} \in \{1, 2, 3\}$ for 3DM</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(xii) $k^\mu = 2$ and $N^{n/\mu} \in \{3, 4, 5\}$ for 3DM</figcaption></figure>

Figure 6: Evolution of the $L^2(\Omega)$-error of the quasi-Trefftz method in 2D and 3D, with respect to the number $N$ of cells per direction, for different polynomial degrees $k^\mu$ of BCs and nano-mesh refinements (through the number of nano-cells per micro-face, $N^{n/\mu}$). We consider a FR local solver of degree $k^n = 0$ (+), $k^n = 1$ (+), $k^n = 2$ (+), $k^n = 3$ (x), $k^n = 4$ (x) and $k^n = 5$ (x).
**Remark 4.1.** We highlight the fact that the optimal asymptotic convergence rate seems to be $k^\mu + 2$ for the 2D Helmholtz equation, and $k^\mu + 1$ for the 3D Maxwell system. This improved convergence in the case of the Helmholtz equation is likely to be due to the usually stronger regularity of its solutions (in comparison to the ones of the Maxwell's equations), or the fact the presented error is $\|p - p_h\|_{0,\Omega}$: the first-order structure of the system is likely to ensure a robust behaviour of $v = (i\kappa\mu)^{-1}\nabla p$, and then offer an improved convergence of $p$.

In a more global point of view, considering higher-order polynomial boundary conditions is likely to allow a quicker convergence, with larger asymptotic rates (as the number of basis functions per cell is also increased). Yet, this comes with the need of a more accurate local solver, and then a costlier auxiliary method. Thus, a balance must be found between pre-computation cost and expected accuracy of the full quasi-Trefftz method, even if we recall the local solver is to be called only for different cell configurations (which can be really limited if the domain is constituted of large homogeneous regions).

### 4.2.2 $p$-convergence

Then, we discuss the convergence in basis functions of the quasi-Trefftz method: we fix the mesh $\mathcal{T}_h$ with $N = 4$ cells in each direction, the polynomial degree $k^\mu$ of the BCs and the parameters of the local solver (through the polynomial degree $k^n$ of the approximation in each nano-cell, and the number $N^{n/\mu}$ of nano-cells per micro-face for the nano-mesh refinement), while the number $N^\mu$ of micro-faces per direction (and then the number $N_{dof}^T$ of local basis functions, see Remark 3.12) is increased. Naturally, if $N^{n/\mu}$ is fixed, refining the micro-mesh $\mathcal{T}_h^\mu$ on the boundary implies that the nano-mesh $\mathcal{T}_h^n$ of the cell is adequately refined.

Then, Figure 7 presents the evolution of the relative $L^2(\Omega)$-error with respect to $N^\mu$, for the 2D Helmholtz and 3D Maxwell's equations respectively. Once again, in case some curves are not drawn, this means that the GMRES solver failed to converge (or the local solver was too expensive, in the 3D configuration).

The presented results are very close to the conclusions we drew in Subsection 4.2.1 for $h$-convergence. First, when the local solver is accurate enough, a common convergence is verified, no matter the parameters $k^n$ and $N^{n/\mu}$: as the quasi-solutions are close enough to actual ones, we recover the abilities of the LPIT functions. Nevertheless, reaching this 'ideal' convergence calls for $k^n > k^\mu$, while it is maintained up to a given micro-mesh refinement, after which a smaller convergence rate is obtained. Furthermore, this 'optimal' range can be widened by increasing either the order of the auxiliary solver (i.e. $k^n$) or by refining its associated nano-mesh (i.e. by increasing $N^{n/\mu}$).

**Remark 4.2.** In the case of $p$-convergence, the optimal asymptotic convergence rate seems to be $k^\mu + 2.5$ for the 2D Helmholtz equation, and $k^\mu + 1.5$ for the 3D Maxwell system. The fact it is improved of 0.5 with respect to the $h$-convergence is then coherent with the error estimate 3.17.

Moreover, it seems the common convergence regime is less demanding for the Maxwell's equations than the Helmholtz ones, as the separation from this optimal regime happens quicklier for the latter.

Globally, choosing high-order BCs allows to potentially get access to larger convergence rates, even if this comes at the cost of a more refined (and then costly) auxiliary solver.
<figure><img src="image_31.png" /><figcaption>(i) $k^\mu = 1$ and $N^{n/\mu} \in \{1, 2, 3\}$ for 2DH</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(ii) $k^\mu = 1$ and $N^{n/\mu} \in \{2, 3, 4\}$ for 2DH</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(iii) $k^\mu = 1$ and $N^{n/\mu} \in \{4, 5, 6\}$ for 2DH</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(iv) $k^\mu = 2$ and $N^{n/\mu} \in \{1, 2, 3\}$ for 2DH</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(v) $k^\mu = 2$ and $N^{n/\mu} \in \{2, 3, 4\}$ for 2DH</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(vi) $k^\mu = 2$ and $N^{n/\mu} \in \{4, 5, 6\}$ for 2DH</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(vii) $k^\mu = 1$ and $N^{n/\mu} = 1$ for 3DM</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(viii) $k^\mu = 1$ and $N^{n/\mu} \in \{1, 2\}$ for 3DM</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(ix) $k^\mu = 2$ and $N^{n/\mu} = 1$ for 3DM</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(x) $k^\mu = 2$ and $N^{n/\mu} \in \{1, 2\}$ for 3DM</figcaption></figure>

Figure 7: Evolution of the $L^2(\Omega)$-error of the quasi-Trefftz method in 2D and 3D, with respect to the number $N^\mu$ of micro-faces per direction, for different polynomial degrees $k^\mu$ of BCs and nano-mesh refinements (through the number of nano-cells per micro-face, $N^{n/\mu}$). We consider a FR local solver of degree $k^n = 0$ (+), $k^n = 1$ (-), $k^n = 2$ (+), $k^n = 3$ (±), $k^n = 4$ (±) and $k^n = 5$ (±).
### 4.2.3 Calibration of the local solver

Both numerical convergence studies realised in this Subsection allowed to highlight that, when the local solver is refined enough, the quasi-Trefftz method recovers the results obtained with the Trefftz method for LPIT basis functions. Nevertheless, this desirable behaviour calls for the consideration of polynomial approximations of order $k^n \ge k^\mu + 1$, and its validity range increases with $k^n$ and $N^{n/\mu}$.

Consequently, when using this quasi-Trefftz method, to our opinion, one should consider a local solver of order at least $k^\mu + 1$ (even if $k^n = k^\mu + 2$ seems appreciable to get the best from of the method) and with at least 2 nano-cells per micro-face, so as to ensure the local approximate solution is sufficiently close to an actual one. Moreover, if additional computational cost is available, it seems more favourable to increase the polynomial degree $k^n$ than refining the nano-mesh $\mathcal{T}^n$: a balance must be found between computational resources (both in terms of hardware and solving time) and the expected precision of the method.

## 4.3 Effects on the iterative convergence

Then, we propose an analysis of the influence of quasi-solutions on the iterative solving process, which is adapted to such a method. Indeed, as the quasi-Trefftz is likely to be less accurate than the LPIT one (out of the common optimal convergence regime), it should not additionally lead to an unpractical approach when it comes to iteratively solve the linear system.

<figure><img src="image_41.png" /><figcaption>Figure 8: Comparison of the iterative GMRES convergence of the preconditioned Trefftz method according to the FR solver refinement, in comparison to the LPIT (■): in terms of polynomial degree $k^n = 0$ (■), $k^n = 1$ (■), $k^n = 2$ (■), $k^n = 3$ (■), $k^n = 4$ (■) and $k^n = 5$ (■), and number of nano-cells per micro-face $N^{n/\mu} = 1$ (■), $N^{n/\mu} = 2$ (■■), $N^{n/\mu} = 3$ (■■■) and $N^{n/\mu} = 4$ (■■).</figcaption></figure>
In order to discuss this legitimate problematic, we consider the Helmholtz equations for a domain $\Omega = [0, L]^2$ (with $L \in \{10, 20, 40\}$) split in $N = 20$ cells in each direction, with $\kappa = 2\pi$ and an exact solution made up of a sum of 10 random plane waves. The local boundary conditions are discretised thanks to a boundary micro-mesh of $N^\mu = 3$ micro-faces per face, with piecewise polynomial BCs of degree $k^\mu \in \{1, 2\}$. The quasi-Trefftz basis functions are then obtained thanks to a DG solver of order $k^n \in [0, 5]$, associated to a nano-mesh of $N^{n/\mu}$ nano-cells per micro-face (leading to $N^\mu \times N^{n/\mu}$ nano-cells per direction).

Figure 8 then presents the evolution of the residual of a restarted GMRES($N_{Kry}$) method (with $N_{Kry} = 20$) applied to the preconditioned Trefftz matrix, for different local solver refinements. First of all, one can observe that considering larger cells (when $L$ increases, as the number $N$ cells is fixed) calls for a more accurate local solver so that the method converges: the approximations are then more demanding, as it could be expected.

Nevertheless, in case it converges (considering $k^n \ge 3$ here), it seems the quasi-Trefftz methods generally need less iterations in comparison to the LPIT (when an exact solver is considered). Consequently, even if the quasi-Trefftz method is likely to be less precise than the LPIT one (when it is out of the common optimal regime), the solution of its associated linear system is usually cheaper: this potential accuracy loss is 'counterbalanced' by an accelerated iterative solution.

## 4.4 An illustrative example: scattering of a plane wave by a metallic structure

As a more practical example, we propose a numerical illustration in the case of a realistic configuration: the scattering of a plane wave on a perfect-metal object in 2D (see Figure 9). The reference solution is computed thanks to a BEM method, and is used to define exact BCs.

<figure><img src="image_42.png" /><figcaption>Figure 9: Illustration of a scattering configuration (total pressure field with plane wave incident field) on a perfect-metal object (in black).</figcaption></figure>

The visualisation of the total solution field is presented in Figure 9 in a domain $\Omega = [-2, 2]^2$ with $N = 24$ cells per direction and $\kappa = 10\pi$, while Figure 10 shows the evolution of the volumic $L^2(\Omega)$-error of the quasi-Trefftz method for different bases, with respect to the number of basis functions per cell (by increasing $N^\mu$ with a fixed mesh). In particular, in this illustration, we consider a converged local solver.

As introduced in Subsection 2.3, one can observe that the plane waves basis induces a limited convergence, as the numerical dependency leads to an explosion of the error when the number of basis functions goes beyond 30. This phenomenon can be mitigated by considering a reduction procedure (which consists in only keeping the eigenvalues larger than a given threshold [24]) which avoids the increase of the error, but does not allow to recover convergence. On the contrary, the LPIT bases show robust convergence of the error, for the different polynomial orders $k^\mu$ of the local boundary conditions. Yet, in a configuration where the exact solution may present limited regularity because of corner/edge singularities, all values of the BCs polynomial order $k^\mu$ lead to similar convergence.
<figure><img src="image_43.png" /><figcaption>Figure 10: Evolution of the normalised volumic L²(Ω)-error of the Trefftz method, with respect to the number $N_{dof}^T$ of basis functions per cell, for the 2D scattering configuration: plane waves (⊗), reduced plane waves (⊕), and LPIT for $k^\mu = 1$ (⊕), $k^\mu = 2$ (⊕) and $k^\mu = 3$ (⊕).</figcaption></figure>

# 5 Conclusion and openings

In this article, we first introduced a formalism adapted from the two-fields Friedrichs systems to the time-harmonic framework, corresponding to systems where the solution variable can be split into two distinct parts which respectively depend on the other one. Thus, this allowed to encompass different wave propagation problems under a common form (referred as *Block Antidiagonal Two-Fields Friedrichs Systems*), for which an associated Trefftz method has been defined. This approach relies on the use of local solutions as basis functions and, leaning on the Ultra Weak Variational Formulation, very favourable properties were highlighted, in particular with a view to resorting to iterative algorithms for its solution. Nevertheless, its classic discretisation thanks to plane waves is used to leading to numerical dependence phenomena, and then to a saturated accuracy of the method.

Facing this limitation, we proposed a new discretisation principle, by relying on the characterisation of local solutions by their local impedance boundary conditions (through a well-posed local problem). These variables can then be discretised thanks to piecewise polynomial fields, leading to basis functions referred as *Local Polynomial Impedance Trefftz (LPIT) functions*: an *a priori* error estimate of the associated Trefftz method could thus be derived. Nevertheless, the exact local solver on which it relies is not practically available (as no analytic expression of such functions can be determined), which calls for the introduction of an approximate one in each cell (using a Discontinuous Galerkin method for example): the associated approximations of the LPIT functions (being no longer actual local solutions, but quasi-ones) are then referred as *Local Polynomial Impedance Quasi-Trefftz (LPIQT) ones*. The abilities of such a quasi-Trefftz method were then numerically assessed and characterised in terms of local solver parameters with respect to the discretisation of the boundary conditions. Globally, the basis independence, *h* and *p*-convergence experiments highlighted robust abilities of this quasi-Trefftz approach, when the local solver is refined enough. More precisely, a common optimal regime is then shared when the local solver is sufficiently accurate, allowing to recover the abilities of the LPIT functions (the local solver is then perceived as 'exact'). Last but not least, in case of precision loss out of this desirable regime, it seems to be counterbalanced by an improved iterative convergence.

In the end, this quasi-Trefftz approach offers a robust behaviour in terms of accuracy, basis independence and iterative behaviour for the solution of time-harmonic wave propagation problems.

However, we observed an influence of the choice of the local impedance field $Z_{\partial T}$ on the abilities of the method, paving the way to a potential optimisation of this parameter, which may lie on an *a priori* error estimate of the method. Moreover, other numerical methods may be considered as approximate local solvers and could offer improved performances for this kind of applications.
References

[1] T. Weiland, A discretization model for the solution of Maxwell's equations for six-component fields, Archiv Elektronik und Uebertragungstechnik 31 (1977) 116-120.
URL https://ui.adsabs.harvard.edu/abs/1977ArE1U..31..116W

[2] K. Beilenhoff, W. Heinrich, H. Hartnagel, Improved finite-difference formulation in frequency domain for three-dimensional scattering problems, IEEE Transactions on Microwave Theory and Techniques 40 (3) (1992) 540-546. doi:10.1109/22.121730.

[3] N. J. Champagne, J. G. Berryman, H. M. Buettner, FDFD: A 3D Finite-Difference Frequency-Domain Code for Electromagnetic Induction Tomography, Journal of Computational Physics 170 (2) (2001) 830-848. doi: 10.1006/jcph.2001.6765.

[4] J. C. Nedelec, Mixed Finite Elements in IR3, Numerische Mathematik 35 (1980) 315-342. doi:10.1007/BF01396415.

[5] J. C. Nédélec, A new family of mixed finite elements in $\mathbb{R}^3$, Numerische Mathematik 50 (1) (1986) 57-81.
doi:10.1007/BF01389668.

[6] P. Monk, A finite element method for approximating the time-harmonic Maxwell equations, Numerische Mathematik 63 (2) (1992) 243-262. doi:10.1007/BF01385860.

[7] P. Monk, Finite Element Methods for Maxwell's Equations, Clarendon Press, 2003. doi:10.1093/acprof:oso/9780198508885.001.0001.

[8] A. Buffa, R. Hiptmair, Galerkin Boundary Element Methods for Electromagnetic Scattering, in: M. Ainsworth, P. Davies, D. Duncan, B. Rynne, P. Martin (Eds.), Topics in Computational Wave Propagation: Direct and Inverse Problems, Springer, Berlin, Heidelberg, 2003, pp. 83-124. doi:10.1007/978-3-642-55483-4_3.

[9] A. Buffa, R. Hiptmair, T. v. Petersdorff, C. Schwab, Boundary Element Methods for Maxwell Transmission Problems in Lipschitz Domains, Numerische Mathematik 95 (3) (2003) 459-485. doi:10.1007/s00211-002-0407-z.

[10] I. M. Babuška, F. Ihlenburg, E. T. Paik, S. A. Sauter, A Generalized Finite Element Method for Solving the Helmholtz Equation in Two Dimensions with Minimal Pollution, Technical Note BN-1179, Institute for Physical Science and Technology, University of Maryland at College Park (1994). doi:10.1016/0045-7825(95) 00890-X.

[11] O. G. Ernst, M. J. Gander, Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods, in: I. G. Graham, T. Y. Hou, O. Lakkis, R. Scheichl (Eds.), Numerical Analysis of Multiscale Problems, Lecture Notes in Computational Science and Engineering, Springer, Berlin, Heidelberg, 2012, pp. 325-363.
doi:10.1007/978-3-642-22061-6_10.

[12] E. Trefftz, Ein gegenstuck zum ritzschen verfahren, Proceedings of the 2nd International Congress of Applied Mechanics, Orell Fussli Verlag (1926) 131—137.

[13] R. Hiptmair, A. Moiola, I. Perugia, A Survey of Trefftz Methods for the Helmholtz Equation, in: G. R. Barrenechea, F. Brezzi, A. Cangiani, E. H. Georgoulis (Eds.), Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Springer International Publishing, Cham, 2016, pp. 237-279. doi:10.1007/978-3-319-41640-3_8.

[14] B. Després, Sur une formulation variationnelle de type ultra-faible, C. R. Acad. Sci. Paris Sér. I Math. 318 (1994) 939-944.

[15] O. Cessenat, Application d'une nouvelle formulation variationnelle aux équations d'ondes harmoniques : problèmes de Helmholtz 2D et de Maxwell 3D, These de doctorat, Paris 9 (Jan. 1996).
URL https://www.theses.fr/1996PA090039

[16] P. Helluy, Résolution numérique des équations de Maxwell harmoniques par une méthode d'éléments finis discontinus, phdthesis, Ecole nationale superieure de l'aeronautique et de l'espace (Jan. 1994).
URL https://theses.hal.science/tel-00657828
[17] P. Monk, D.-Q. Wang, A least-squares method for the Helmholtz equation, Computer Methods in Applied Mechanics and Engineering 175 (1) (1999) 121-136. doi:10.1016/S0045-7825(98)00326-0.

[18] A. Buffa, P. Monk, Error estimates for the Ultra Weak Variational Formulation of the Helmholtz equation, ESAIM: Mathematical Modelling and Numerical Analysis 42 (6) (2008) 925-940, number: 6 Publisher: EDP Sciences. doi:10.1051/m2an:2008033.

[19] A. Moiola, Trefftz-Discontinuous Galerkin Methods for Time-Harmonic Wave Problems, Doctoral Thesis, ETH Zurich, accepted: 2017-08-10T11:43:11Z (2011). doi:10.3929/ethz-a-006698757.

[20] R. Hiptmair, A. Moiola, I. Perugia, Plane Wave Discontinuous Galerkin Methods for the 2D Helmholtz Equation: Analysis of the p-Version, SIAM Journal on Numerical Analysis Publisher: Society for Industrial and Applied Mathematics (Feb. 2011). doi:10.1137/090761057.

[21] R. Hiptmair, A. Moiola, I. Perugia, Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations, Mathematics of Computation 82 (281) (2012) 247-268. doi:10.1090/S0025-5718-2012-02627-5.

[22] R. Hiptmair, A. Moiola, I. Perugia, Plane Wave Discontinuous Galerkin Methods: Exponential Convergence of the h p -Version, Foundations of Computational Mathematics 16 (3) (2016) 637-675. doi: 10.1007/s10208-015-9260-1.

[23] S. Congreve, J. Gedicke, I. Perugia, Numerical Investigation of the Conditioning for Plane Wave Discontinuous Galerkin Methods, in: F. A. Radu, K. Kumar, I. Berre, J. M. Nordbotten, I. S. Pop (Eds.), Numerical Mathematics and Advanced Applications ENUMATH 2017, Springer International Publishing, Cham, 2019, pp. 493-500. doi:10.1007/978-3-319-96415-7_44.

[24] H. Barucq, A. Bendali, J. Diaz, S. Tordeux, Local strategies for improving the conditioning of the plane-wave Ultra-Weak Variational Formulation, Journal of Computational Physics 441 (2021) 110449. doi:10.1016/j.jcp.2021.110449.

[25] S. Pernet, M. Sirdey, S. Tordeux, Ultra-weak variational formulation for heterogeneous maxwell problem in the context of high performance computing, ESAIM: Proceedings and Surveys Publisher: EDP Sciences (2023). doi:10.1051/proc/202375096.

[26] R. Tezaur, L. Zhang, C. Farhat, A discontinuous enrichment method for capturing evanescent waves in multi-scale fluid and fluid/solid problems, Computer Methods in Applied Mechanics and Engineering 197 (19) (2008) 1680-1698. doi:10.1016/j.cma.2007.08.023.

[27] T. Luostari, T. Huttunen, P. Monk, Improvements for the ultra weak variational formulation, International Journal for Numerical Methods in Engineering 94 (6) (2013) 598-624. doi:10.1002/nme.4469.

[28] E. Parolin, D. Huybrechs, A. Moiola, Stable approximation of Helmholtz solutions in the disk by evanescent plane waves, ESAIM: Mathematical Modelling and Numerical Analysis 57 (6) (2023) 3499-3536, number: 6 Publisher: EDP Sciences. doi:10.1051/m2an/2023081.

[29] L.-M. Imbert-Gérard, B. Després, A generalized plane wave numerical method for smooth non constant coefficients, IMA Journal of Numerical Analysis 34 (Jul. 2014). doi:10.1093/imanum/drt030.

[30] L.-M. Imbert-Gérard, Interpolation properties of generalized plane waves, Numerische Mathematik 131 (4) (2015) 683-711. doi:10.1007/s00211-015-0704-y.

[31] L.-M. Imbert-Gérard, Amplitude-based Generalized Plane Waves: New Quasi-Trefftz Functions for Scalar Equations in two dimensions, SIAM Journal on Numerical Analysis 59 (3) (2021) 1663-1686, publisher: Society for Industrial and Applied Mathematics. doi:10.1137/20M136791X.

[32] L.-M. Imbert-Gérard, G. Sylvand, Three types of quasi-Trefftz functions for the 3D convected Helmholtz equation: construction and approximation properties, IMA Journal of Numerical Analysis 45 (4) (2025) 2274-2328. doi:10.1093/imanum/drae060.

[33] C. Lehrenfeld, P. Stocker, Embedded Trefftz discontinuous Galerkin methods, International Journal for Numerical Methods in Engineering 124 (17) (2023) 3637-3661. doi:10.1002/nme.7258.
[34] H. S. Fure, S. Pernet, M. Sirdey, S. Tordeux, A discontinuous Galerkin Trefftz type method for solving the two dimensional Maxwell equations, SN Partial Differential Equations and Applications 1 (23) (2020) 19. doi:10.1007/s42985-020-00024-0.

[35] S. Pernet, M. Rivet, S. Tordeux, Trefftz Method for a Class of Time-Harmonic Two-Fields Friedrichs Systems, Tech. Rep. RR-9602, INRIA Bordeaux - Sud-Ouest (Nov. 2025). URL https://inria.hal.science/hal-05368909

[36] K. O. Friedrichs, Symmetric positive linear differential equations, Communications on Pure and Applied Mathematics 11 (3) (1958) 333-418. doi:10.1002/cpa.3160110306.

[37] A. Ern, J.-L. Guermond, G. Caplain, An Intrinsic Criterion for the Biancectivity of Hilbert Operators Related to Friedrichs' Systems, Communications in Partial Differential Equations 32 (2007) 317-341. doi:10.1080/03605300600718545.

[38] D. A. Di Pietro, A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Vol. 69 of Mathématiques et Applications, Springer, Berlin, Heidelberg, 2012. doi:10.1007/978-3-642-22980-0.

[39] N. Antonić, K. Burazin, I. Crnjac, M. Erceg, Complex Friedrichs systems and applications, Journal of Mathematical Physics 58 (10) (2017) 101508. doi:10.1063/1.5005608.

[40] S. Pernet, M. Rivet, S. Tordeux, Unique continuation theorem for symmetric harmonic hyperbolic systems with constant coefficients, Tech. Rep. RR-9544, Inria Bordeaux - Sud-Ouest (Mar. 2024). URL https://inria.hal.science/hal-04504052

[41] A. Ern, J. Guermond, Discontinuous Galerkin Methods for Friedrichs' Systems. Part II. Second-order Elliptic PDEs, SIAM Journal on Numerical Analysis 44 (6) (2006) 2363-2388, publisher: Society for Industrial and Applied Mathematics. doi:10.1137/05063831X.

[42] A. Ern, J.-L. Guermond, Discontinuous Galerkin Methods for Friedrichs' Systems. Part III. Multifield Theories with Partial Coercivity, SIAM Journal on Numerical Analysis 46 (2) (2008) 776-804, publisher: Society for Industrial and Applied Mathematics. doi:10.1137/060664045.

[43] D. Colton, R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Vol. 93 of Applied Mathematical Sciences, Springer International Publishing, Cham, 2019. doi:10.1007/978-3-030-30351-8.

[44] M. Sirdey, Méthode itérative de Trefftz pour la simulation d'ondes électromagnétiques en trois dimensions., Theses, Université de Pau et des Pays de l'Adour (Dec. 2022). URL https://theses.hal.science/tel-04172930

[45] I. Djiba, Méthode de décomposition de domaine par méthode de Trefftz pour la propagation des ondes en géophysique., Theses, Université de Pau et des Pays de l'Adour (2025). URL https://theses.fr/s342267

[46] J. Nečas, Direct Methods in the Theory of Elliptic Equations, Springer Monographs in Mathematics, Springer, Berlin, Heidelberg, 2012. doi:10.1007/978-3-642-10455-8.

[47] O. Cessenat, B. Despres, Application of an Ultra Weak Variational Formulation of Elliptic PDEs to the Two-Dimensional Helmholtz Problem, SIAM Journal on Numerical Analysis 35 (1) (1998) 255-299, publisher: Society for Industrial and Applied Mathematics. doi:10.1137/S0036142995285873.

[48] P. R. S. Antunes, S. S. Valtchev, A meshfree numerical method for acoustic wave propagation problems in planar domains with corners and cracks, Journal of Computational and Applied Mathematics 234 (9) (2010) 2646-2662. doi:10.1016/j.cam.2010.01.031.

[49] A. Moiola, R. Hiptmair, I. Perugia, Plane wave approximation of homogeneous Helmholtz solutions, Zeitschrift für angewandte Mathematik und Physik 62 (5) (2011) 809-837. doi:10.1007/s00033-011-0147-y.

[50] F. J. Sayas, T. S. Brown, M. E. Hassell, Variational Techniques for Elliptic Partial Differential Equations: Theoretical Tools and Advanced Applications, CRC Press, Boca Raton, 2019. doi:10.1201/9780429507069.

[51] I. Graham, S. Sauter, Stability and finite element error analysis for the Helmholtz equation with variable coefficients, Mathematics of Computation 89 (321) (2020) 105-138. doi:10.1090/mcom/3457.
[52] G. Cohen, S. Pernet, Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations, Scientific Computation, Springer Netherlands, Dordrecht, 2017. doi:10.1007/978-94-017-7761-2.

[53] S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, Vol. 15 of Texts in Applied Mathematics, Springer New York, New York, NY, 2008. doi:10.1007/978-0-387-75934-0.

[54] H. Barucq, A. Bendali, M. Fares, V. Mattesi, S. Tordeux, A symmetric Trefftz-DG formulation based on a local boundary element method for the solution of the Helmholtz equation, Journal of Computational Physics 330 (2017) 1069–1092. doi:10.1016/j.jcp.2016.09.062.

[55] A. Modave, T. Chaumont-Frelet, A hybridizable discontinuous Galerkin method with characteristic variables for Helmholtz problems, Journal of Computational Physics 493 (2023) 112459. doi:10.1016/j.jcp.2023.112459.

[56] S. Pescuma, G. Gabard, T. Chaumont-Frelet, A. Modave, A hybridizable discontinuous Galerkin method with transmission variables for time-harmonic acoustic problems in heterogeneous media, Journal of Computational Physics 534 (2025) 114009. doi:10.1016/j.jcp.2025.114009.

[57] Y. Saad, M. H. Schultz, GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM Journal on Scientific and Statistical Computing 7 (3) (1986) 856–869, publisher: Society for Industrial and Applied Mathematics. doi:10.1137/0907058.

[58] R. Morgan, GMRES with deflation restarting, Siam Journal on Scientific Computing 24 (Aug. 2002). doi: 10.1137/S1064827599364659.

[59] Y. Saad, Iterative Methods for Sparse Linear Systems, Other Titles in Applied Mathematics, Society for Industrial and Applied Mathematics, 2003. doi:10.1137/1.9780898718003.