# Trefftz Method for a Class of Time-Harmonic Two-Fields Friedrichs Systems

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

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

## To cite this version:

Sébastien Pernet, Matthias Rivet, Sébastien Tordeux. Trefftz Method for a Class of Time-Harmonic Two-Fields Friedrichs Systems. RR-9602, INRIA Bordeaux - Sud-Ouest. 2025. hal-05368909

HAL Id: hal-05368909

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

Submitted on 17 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inria

# Trefftz Method for a Class of Time-Harmonic Two-Fields Friedrichs Systems

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

## RESEARCH REPORT

### N° 9602

November 2025

Project-Team Makutu
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

# Trefftz Method for a Class of Time-Harmonic
Two-Fields Friedrichs Systems

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

Project-Team Makutu

Research Report n° 9602 — November 2025 — 47 pages

**Abstract:** This document presents a class of two-fields Friedrichs systems, allowing to encompass classic time-harmonic wave propagation problems into a unique formalism: notions of incoming and outgoing traces, in addition to a normal flux decomposition and a consistent numerical flux expression are defined in this general setting. Then, a Trefftz method is introduced for this class of problems: formulations based on the reciprocity formula or the Ultra Weak Variational Formulation are defined, allowing to prove weak-coercivity and contraction properties of the preconditioned system. Finally, we discuss the interpretation of the method from the point of view of Domain Decomposition methods, and recall classic error estimates results and discretisation by plane waves.

**Key-words:** Friedrichs systems, Trefftz method, Time-harmonic wave problems, Ultra Weak Variational Formulation, Numerical flux, Error estimate.

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

<footnote>* DTIS, ONERA, Université de Toulouse, 31000, Toulouse, France</footnote>

<footnote>† EPCI Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142, France</footnote>
# Méthode de Trefftz pour une classe de systèmes de Friedrichs harmoniques

**Résumé** : Ce document présente une classe de systèmes de Friedrichs à deux champs, permettant d'inclure les problèmes d'ondes harmoniques classiques dans un unique formalisme : des notions de traces entrantes et sortantes, ainsi qu'une décomposition du flux normal et une expression d'un flux numérique consistant sont ainsi définies de manière générale. Une méthode de Trefftz est alors introduite pour cette classe de problèmes : des formulations basées sur la formule de réciprocity ou la formulation ultra-faible sont ainsi définies, permettant de prouver des propriétés de coercivité faible et de contraction du système préconditionné. Enfin, nous discutons l'interprétation de cette méthode à l'aune des méthodes de Décomposition de Domaines, et rappelons des résultats d'estimation d'erreur classiques et de discrétisation par ondes planes.

**Mots-clés** : Systèmes de Friedrichs, Méthode de Trefftz, Problèmes d'onde harmoniques, Formulation Ultra Faible, Flux numérique, Estimation d'erreur.
# Contents

1 Introduction 4

2 Block Antidiagonal Two-Fields Friedrichs Systems 4

2.1 Introduction of the block antidiagonal two-fields Friedrichs systems 4

2.2 Explicit decomposition of the normal flux 6

2.2.1 Rewriting in terms of outgoing/incoming traces 9

2.3 Application in the framework of a mesh 9

2.3.1 Rewriting in terms of outgoing/incoming traces 11

2.4 Complete BATFFS Problem 12

2.5 A generic methodology to define consistent numerical fluxes for BATFFS 12

2.5.1 Continuity of the normal flux 13

2.5.2 Numerical flux computation 14

2.5.3 Numerical fluxes for BATFFS in a homogeneous configuration 15

2.5.4 Numerical fluxes for BATFFS in a heterogeneous configuration 16

3 Introduction to Trefftz methods for Block Antidiagonal Two-Fields Friedrichs Systems 18

3.1 Approximation spaces 20

3.2 Derivation of a Trefftz formulation from the reciprocity formula 21

3.3 Weak-coercivity property and well-posedness of the Trefftz formulation 23

3.4 Equivalence to an UWVF 27

3.4.1 Introduction of the traces and function spaces 27

3.4.2 An equivalent Ultra Weak Variational Formulation 30

3.5 The natural contracting preconditioner of Cessenat-Després 33

3.5.1 Contraction properties of the operators 33

3.5.2 Introduction of a Jacobi algorithm 35

3.5.3 Strict contraction of the discrete preconditioned system 36

3.6 Interpretation of the Trefftz method from a DDM point of view 38

3.7 Trefftz error estimate 40

3.7.1 A Céa's Lemma for the Trefftz formulation 40

3.7.2 Classic arguments of volumic bounding for wave propagation problems 42

3.8 Classic plane waves discretisation and numerical limits 43
# 1 Introduction

In this Report, we are interested in time-harmonic wave propagation problems (with a time convention in the form $t \mapsto e^{i\omega t}$), as the 3D Maxwell and 1D/2D Helmholtz equations, for which we aim at defining and studying a Trefftz formulation and its properties. Nevertheless, so as to avoid redundancy in the treatments of these different problems, we first propose a common and generic formalism, inspired by the Friedrichs systems, that encompasses all of them at the same time. This first-order point of view generally reads

$$ \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 one has $\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}$.

In this study, we suppose a piecewise constant distribution of the heterogeneity, in the sense that 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 $\mathcal{M}$ is constant in each $\Omega_j$. Then, we consider a conformal polygonal/polyhedral mesh $\mathcal{T}_h$, in the sense that each mesh cell is included in a single subdomain: $\mathcal{M}$ is then also constant in each cell (and the solution restriction is denoted as $\mathbf{y}^T := \mathbf{y}|_T$), while the flux matrices $(\mathcal{F}_j)_{j \in [1, d]}$ are fixed.

Then, being given this generic formalism (allowing to deal with the different time-harmonic wave problems at the same time), we introduce an adapted Trefftz method from the reciprocity formula, and following the Ultra Weak Variational Formulation. Its properties, interpretation in terms of Domain Decomposition Methods, classic error estimates and discretisation with plane waves are then proposed and detailed, allowing to recover previous results for the Helmholtz and Maxwell equations in this unified framework.

# 2 Block Antidiagonal Two-Fields Friedrichs Systems

The wave propagation problems we are interested in can be interpreted as time-harmonic Friedrichs systems (following the first-order formalism of Section 1). Yet, they specifically share a common crossed structure, related to their original formulation as second-order problems: the flux operators are all block 'antidiagonal', as they have a zero block diagonal. This case is part of the so-called two-fields Friedrichs systems, which were introduced and studied in [11, 10, 1]. This justifies the following denomination in the rest of the Report: Block Antidiagonal Two-Fields Friedrichs Systems (BATFFS).

This specific configuration will allow to develop explicit expressions of the normal flux decomposition, which justifies its dedicated introduction.

## 2.1 Introduction of the block antidiagonal two-fields Friedrichs systems

In the idea that the problems we consider are posed on a couple of two variables ($\mathbf{e}$ and $\mathbf{h}$ for the Maxwell's equations; $p$ and $\mathbf{v}$ for the Helmholtz equation), we suppose the unknown $\mathbf{y}$ can be split in 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}, \quad (1) $$
where the parameters $\varepsilon$ and $\mu$ satisfy a piecewise constant heterogeneity distribution. Thus, considering a mesh cell $T \subset \Omega$ conformal with the heterogeneity partition of the domain, these fields are constant in $T$ with values $\varepsilon_T$ and $\mu_T$. This allows to define the local matrix $M_T$ as

$$M_T := \begin{bmatrix} \varepsilon_T \mathcal{I}_{m_1} & \mathcal{O}_{m_1, m_2} \\ \mathcal{O}_{m_2, m_1} & \mu_T \mathcal{I}_{m_2} \end{bmatrix}.$$

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 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 has:

$$\mathcal{F}_{\mathbf{n}} = \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 (2)$$

Thus, time-harmonic Friedrichs systems showing this structure will be referred as *Block Antidiagonal Two-Fields Friedrichs Systems (BATFFS)* in the rest of the Report.

**Remark 2.1.** The three wave propagation problems we consider in this Report all follow this form, as:

* for the Maxwell system, one has: $\mathbf{y}_1 = \mathbf{e}$ ($m_1 = 3$), $\mathbf{y}_2 = \mathbf{h}$ ($m_2 = 3$) and $\mathcal{B}_{\mathbf{n}} = -\Gamma_{\times}[\mathbf{n}] \in \mathbb{R}^{3 \times 3}$, defined as

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

* for the 2D Helmholtz equation, one has: $\mathbf{y}_1 = \mathbf{p}$ ($m_1 = 1$), $\mathbf{y}_2 = \mathbf{v}$ ($m_2 = 2$) and $\mathcal{B}_{\mathbf{n}} = [-\mathbf{n}^*] \in \mathbb{R}^{1 \times 2}$.

* for the 1D Helmholtz equation, one has: $\mathbf{y}_1 = \mathbf{p}$ ($m_1 = 1$), $\mathbf{y}_2 = \mathbf{v}$ ($m_2 = 1$) and $\mathcal{B}_{\mathbf{n}} = [-\mathbf{n}] \in \mathbb{R}^{1 \times 1}$.

Finally, in the objective of defining numerical methods adapted to this formalism, we suppose a local unique continuation principle (so as to ensure a local well-posedness property).

**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 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 (3)$$

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 general case of BATFFS, when considering the complete cell boundary (i.e. $\Gamma = \partial T$), one can refer to [22]. Moreover, for a more general choice of $\Gamma$ and the wave problems of interest, Hypothesis 2.2 is verified for Maxwell [18] and Helmholtz [7] systems.
## 2.2 Explicit decomposition of the normal flux

In the specific case of BATFFS, we introduce a decomposition procedure, based on a generalised eigenvalue problem, which can be realised in an explicit way. This relies on the choice of a pivotal matrix $\tilde{\mathcal{M}}$ with a similar form to the matrix $\mathcal{M}$ (see (1))

$$ \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}_*^+$.

To do so, we consider a compact singular value decomposition of the matrix $\mathcal{B}_n$, 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 (4) $$

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** (Uniform bounding of $\Sigma$ and $\hat{U}^*$). Even if the previous decomposition (4) is relative to a given unit vector $\mathbf{n}$, its components can be bounded in an independent way of $\mathbf{n}$. Let $\mathbf{n} \in \mathbb{R}^d$ be a unit vector. Then, according to the compact SVD decomposition $\mathcal{B}_n = \mathcal{U}\Sigma\mathcal{V}^*$, one has

$$ \mathcal{B}_n \mathcal{B}_n^* = \hat{U} \Sigma^2 \hat{U}^* \quad \text{with} \quad \mathcal{B}_n = \sum_{j=1}^d n_j \mathcal{B}_j, $$

which implies that its norm can be bounded in an independent way of $\mathbf{n}$ as

$$ \|\mathcal{B}_n \mathcal{B}_n^*\| = \left\| \sum_{j,j'=1}^d n_j n_{j'} \mathcal{B}_j (\mathcal{B}_{j'})^* \right\| \le \sum_{j,j'=1}^d |n_j n_{j'}| \|\mathcal{B}_j\| \|\mathcal{B}_{j'}\| \le \left[ \sum_{j=1}^d \|\mathcal{B}_j\| \right]^2, $$

as $\mathbf{n}$ is unitary.

Moreover, relying on the definition of the norm $\|\mathcal{B}_n \mathcal{B}_n^*\|$ in terms of Rayleigh quotient, we get

$$ \|\mathcal{B}_n \mathcal{B}_n^*\| = \max_{\mathbf{x} \in \mathbb{R}^{m_1 \setminus \{0\}}} \frac{\mathbf{x}^* \mathcal{B}_n \mathcal{B}_n^* \mathbf{x}}{\mathbf{x}^* \mathbf{x}} = \max_{\tilde{\mathbf{x}} \in \mathbb{R}^r \setminus \{0\}} \frac{\mathbf{x}^* \hat{U} \Sigma^2 \hat{U}^* \mathbf{x}}{\mathbf{x}^* \mathbf{x}}. $$

Yet, $\hat{U}^* \mathbb{R}^{m_1} = \mathbb{R}^r$ as $\hat{U}^* \mathbb{R}^{m_1} \subset \mathbb{R}^r$, and for any $\tilde{\mathbf{x}} \in \mathbb{R}^r$, we can define $\mathbf{x} = \hat{U} \tilde{\mathbf{x}} \in \mathbb{R}^{m_1}$ such that $\hat{U}^* \mathbf{x} = \hat{U}^* \hat{U} \tilde{\mathbf{x}} = \tilde{\mathbf{x}}$, because of (5).

This allows to replace $\hat{U}^* \mathbf{x}$ by $\tilde{\mathbf{x}}$ in the previous expression, leading to

$$ \|\mathcal{B}_n \mathcal{B}_n^*\| = \max_{\mathbf{x} \in \mathbb{R}^{m_1 \setminus \{0\}}} \frac{\mathbf{x}^* \hat{U} \Sigma^2 \hat{U}^* \mathbf{x}}{\mathbf{x}^* \mathbf{x}} = \max_{\tilde{\mathbf{x}} \in \mathbb{R}^r \setminus \{0\}} \frac{\tilde{\mathbf{x}}^* \Sigma^2 \tilde{\mathbf{x}}}{\tilde{\mathbf{x}}^* \hat{U}^* \hat{U} \tilde{\mathbf{x}}} = \max_{\tilde{\mathbf{x}} \in \mathbb{R}^r \setminus \{0\}} \frac{\tilde{\mathbf{x}}^* \Sigma^2 \tilde{\mathbf{x}}}{\tilde{\mathbf{x}}^* \tilde{\mathbf{x}}} = \|\Sigma\|^2, $$

where we used, once again, the fact that $\hat{U}^* \hat{U} = \mathcal{I}_r$. Finally, we obtain

$$ \|\Sigma\| \le \sum_{j=1}^{d} \|\mathcal{B}_j\|, $$

which shows that, no matter the unitary vector $\mathbf{n} \in \mathbb{R}^d$, $\Sigma$ is bounded in an independent way.

Concerning the operator $\hat{U}^*$, we first remark that $\hat{U}^* \hat{U}$ is a projection and Hermitian, as (5) implies that

$$ (\hat{U}\hat{U}^*)^2 = \hat{U} \underbrace{\hat{U}^*\hat{U}}_{\mathcal{I}_r} \hat{U}^* = \hat{U}\hat{U}^*, $$
which induces the fact that $||\hat{u}\hat{u}^*|| = 1$. Then, one gets

$$ \forall \mathbf{x} \in \mathbb{C}^{m_1}, \quad \|\hat{u}^*\mathbf{x}\|^2 = \langle \hat{u}^*\mathbf{x}, \hat{u}^*\mathbf{x} \rangle = \langle \hat{u}\hat{u}^*\mathbf{x}, \mathbf{x} \rangle = \langle \hat{u}\hat{u}^*\mathbf{x}, \hat{u}\hat{u}^*\mathbf{x} \rangle = \|\hat{u}\hat{u}^*\mathbf{x}\|^2 = \|\mathbf{x}\|^2, $$

which shows that $||\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 (4), we define the following decomposition of $B_n$:

$$ B_n = uV^* \quad \text{with} \quad \left\{ \begin{array}{ll} u = \hat{u}\Sigma^{1/2} & \text{and } v = \hat{v}\Sigma^{1/2}, \\ u^*u = \Sigma & \text{and } v^*v = \Sigma, \end{array} \right. \qquad (5) $$

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 Report (even if (4) may have been preferred, as it would keep orthonormality for the singular vectors).

It is then possible to realise an explicit decomposition of the normal flux, leaning on the generalised eigenvalue problem

$$ \mathcal{F}_n \mathbf{w} = \lambda \tilde{\mathcal{M}} \mathbf{w}. $$

(i) We first introduce the symmetric positive definite square root $\tilde{\mathcal{M}}^{1/2}$ of $\tilde{\mathcal{M}}$, which leads to

$$ \begin{align*} \mathcal{F}_n \mathbf{w} = \lambda \tilde{\mathcal{M}} \mathbf{w} &\iff \mathcal{F}_n \tilde{\mathcal{M}}^{-1/2} \tilde{\mathcal{M}}^{1/2} \mathbf{w} = \lambda \tilde{\mathcal{M}}^{1/2} \tilde{\mathcal{M}}^{1/2} \mathbf{w} \\ &\iff \tilde{\mathcal{F}}_n \tilde{\mathbf{w}} = \lambda \tilde{\mathbf{w}}, \end{align*} $$

where we consider the intermediary variable $\tilde{\mathbf{w}} := \tilde{\mathcal{M}}^{1/2} \mathbf{w}$ and the real symmetric matrix

$$ \tilde{\mathcal{F}}_n := \tilde{\mathcal{M}}^{-1/2} \mathcal{F}_n \tilde{\mathcal{M}}^{-1/2} = \begin{bmatrix} O_{m_1, m_1} & \tilde{c} U V^* \\ \tilde{c} U V^* & O_{m_2, m_2} \end{bmatrix} \quad \text{with the velocity } \tilde{c} = (\tilde{\varepsilon} \tilde{\mu})^{-1/2}. $$

(ii) We are then led back to a classic eigenvalue problem, for which an explicit decomposition can be derived:

$$ \tilde{\mathcal{F}}_n = Q \Lambda Q^*, $$

where $Q \in \mathbb{R}^{m \times 2r}$ and $\Lambda \in \mathbb{R}^{2r \times 2r}$ respectively denote the matrix of the coordinates of the non-normalised eigenvectors and the diagonal matrix of the associated eigenvalues of $\tilde{\mathcal{F}}_n$, and are defined as

$$ Q := 2^{-1/2} \begin{bmatrix} u & u \\ v & -v \end{bmatrix} \in \mathbb{R}^{m \times 2r} \quad \text{and} \quad \Lambda := \begin{bmatrix} \tilde{c} \mathcal{I}_r & O_{r,r} \\ O_{r,r} & -\tilde{c} \mathcal{I}_r \end{bmatrix} \in \mathbb{R}^{2r \times 2r}. $$

Moreover, the matrix $Q$ is orthogonal, in the sense that

$$ Q^* Q = \frac{1}{2} \begin{bmatrix} U^* & V^* \\ U^* & -V^* \end{bmatrix} \begin{bmatrix} U & U \\ V & -V \end{bmatrix} = \frac{1}{2} \begin{bmatrix} U^* U + V^* V & U^* U - V^* V \\ U^* U - V^* V & U^* U + V^* V \end{bmatrix} = \begin{bmatrix} \Sigma & O_{r,r} \\ O_{r,r} & \Sigma \end{bmatrix}, $$
as we recall that $U^*U = V^*V = \Sigma$ (see (5)).

(iii) Finally, one can come back to the original matrix, leading to

$$ 
\begin{aligned}
\mathcal{F}_n &= \tilde{\mathcal{M}}^{1/2} \tilde{\mathcal{F}}_n \tilde{\mathcal{M}}^{1/2} \\
&= \tilde{\mathcal{M}}^{1/2} \mathcal{Q} \Lambda \mathcal{Q}^* \tilde{\mathcal{M}}^{1/2} \\
&= \tilde{\mathcal{M}} \underbrace{\tilde{\mathcal{M}}^{-1/2} \mathcal{Q}}_{\mathcal{P}} \Lambda \underbrace{\mathcal{Q}^* \tilde{\mathcal{M}}^{-1/2}}_{\mathcal{P}^*} \tilde{\mathcal{M}},
\end{aligned}
$$

where the matrix $\mathcal{P}$ is defined as

$$ \mathcal{P} := \tilde{\mathcal{M}}^{-1/2} \mathcal{Q} = 2^{-1/2} \begin{bmatrix} \tilde{\varepsilon}^{-1/2} \mathcal{U} & \tilde{\varepsilon}^{-1/2} \mathcal{U} \\ \tilde{\mu}^{-1/2} \mathcal{V} & -\tilde{\mu}^{-1/2} \mathcal{V} \end{bmatrix} \in \mathbb{R}^{m \times 2r}, $$

and is $\tilde{\mathcal{M}}$-orthogonal (but not orthonormal), in the sense that:

$$ \mathcal{P}^* \tilde{\mathcal{M}} \mathcal{P} = \mathcal{Q}^* \tilde{\mathcal{M}}^{-1/2} \tilde{\mathcal{M}} \tilde{\mathcal{M}}^{-1/2} \mathcal{Q} = \mathcal{Q}^* \mathcal{Q} = \begin{bmatrix} \mathbf{\Sigma} & \mathcal{O}_{r,r} \\ \mathcal{O}_{r,r} & \mathbf{\Sigma} \end{bmatrix}. \quad (6) $$

This leads to the following decomposition of the normal flux matrix:

$$ \mathcal{F}_n = \tilde{\mathcal{M}} \mathcal{P} \Lambda \mathcal{P}^* \tilde{\mathcal{M}} = \mathcal{F}_n^+ + \mathcal{F}_n^-, \quad (7) $$

where

$$ \mathcal{F}_n^{\pm} = \tilde{\mathcal{M}} \mathcal{P} \Lambda^{\pm} \mathcal{P}^* \tilde{\mathcal{M}} \quad \text{with} \quad \Lambda^{\pm} = \begin{bmatrix} \tilde{c} \mathcal{I}_r & \mathcal{O}_{r,r} \\ \mathcal{O}_{r,r} & \mathcal{O}_{r,r} \end{bmatrix} \quad \text{and} \quad \Lambda^{-} = \begin{bmatrix} \mathcal{O}_{r,r} & \mathcal{O}_{r,r} \\ \mathcal{O}_{r,r} & -\tilde{c} \mathcal{I}_r \end{bmatrix}, \quad (8) $$

or else, by denoting the admittance and impedance as $\tilde{Y} = (\tilde{\varepsilon}/\tilde{\mu})^{1/2}$ and $\tilde{Z} = (\tilde{\mu}/\tilde{\varepsilon})^{1/2}$ respectively

$$ \mathcal{F}_n^{\pm} = \frac{1}{2} \begin{bmatrix} \pm \tilde{Y} \mathcal{U} \mathcal{U}^* & \mathcal{U} \mathcal{V}^* \\ \mathcal{V} \mathcal{U}^* & \pm \tilde{Z} \mathcal{V} \mathcal{V}^* \end{bmatrix}. \quad (9) $$

The matrices $\mathcal{F}_n^{\pm}$ are then, by construction (8), positive and negative respectively.

**Remark 2.5** (Decomposition for the wave propagation problems). In the specific case of the wave propagation problems we study in this Report, an additional property allows to simplify this expression: the only non-zero singular value is 1, implying that $\Sigma = \mathcal{I}_r$. Then, this leads to

$$ \mathcal{U} \mathcal{V}^* = \hat{\mathcal{U}} \Sigma \hat{\mathcal{V}}^* = \mathcal{B}_n \quad \text{and} \quad \mathcal{V} \mathcal{U}^* = \hat{\mathcal{V}} \Sigma \hat{\mathcal{U}}^* = \mathcal{B}_n^*, $$

which can be inserted in (9) to obtain

$$ \mathcal{F}_n^{\pm} = \frac{1}{2} \begin{bmatrix} \pm \tilde{Y} \mathcal{B}_n \mathcal{B}_n^* & \mathcal{B}_n \\ \mathcal{B}_n^* & \pm \tilde{Z} \mathcal{B}_n^* \mathcal{B}_n \end{bmatrix}. \quad (10) $$

In particular, this leads to:
• for the Maxwell system, one has $r = 2$ and

$$ \mathcal{F}_{\mathbf{n}}^{\pm} = \frac{1}{2} \begin{bmatrix} \pm \tilde{Y} \Gamma_t[\mathbf{n}] & -\Gamma_{\times}[\mathbf{n}] \\ \Gamma_{\times}[\mathbf{n}] & \pm \tilde{Z} \Gamma_t[\mathbf{n}] \end{bmatrix}, $$

with

$$ \Gamma_{\times}[\mathbf{n}] = \begin{bmatrix} 0 & -n_3 & n_2 \\ n_3 & 0 & -n_1 \\ -n_2 & n_1 & 0 \end{bmatrix} \quad \text{and} \quad \Gamma_t[\mathbf{n}] = -(\Gamma_{\times}[\mathbf{n}])^2. $$

• for the 2D Helmholtz equation, one has $r = 1$ and

$$ \mathcal{F}_{\mathbf{n}}^{\pm} = \frac{1}{2} \begin{bmatrix} \pm \tilde{Y} & -\mathbf{n}^{*} \\ -\mathbf{n} & \pm \tilde{Z} \mathbf{n} \mathbf{n}^{*} \end{bmatrix}. $$

• for the 1D Helmholtz equation, one has $r = 1$ and

$$ \mathcal{F}_{\mathbf{n}}^{\pm} = \frac{1}{2} \begin{bmatrix} \pm \tilde{Y} & -n \\ -n & \pm \tilde{Z} \end{bmatrix}. $$

### 2.2.1 Rewriting in terms of outgoing/incoming traces

More precisely, the explicit decomposition (9) of the normal flux in the case of the BATFFS can be factorised as

$$ \mathcal{F}_{\mathbf{n}}^{\pm} = \frac{1}{2} \begin{bmatrix} \pm \tilde{Y} \mathcal{U} & \mathcal{O}_{m_1,r} \\ \mathcal{O}_{m_2,r} & \mathcal{V} \end{bmatrix} \begin{bmatrix} \mathcal{U}^{*} & \pm \tilde{Z} \mathcal{V}^{*} \\ \mathcal{U}^{*} & \pm \tilde{Z} \mathcal{V}^{*} \end{bmatrix}. $$

In particular, this expression highlights two quantities, that we refer as outgoing and incoming traces, as they allow to characterise the information going out and coming in respectively. They are respectively defined as

$$ \gamma_{\text{out}}(\mathbf{y}) := \mathcal{U}^* \mathbf{y}_1 + \tilde{Z} \mathcal{V}^* \mathbf{y}_2 \in \mathbb{C}^r \quad \text{and} \quad \gamma_{\text{in}}(\mathbf{y}) := \mathcal{U}^* \mathbf{y}_1 - \tilde{Z} \mathcal{V}^* \mathbf{y}_2 \in \mathbb{C}^r, \quad (11) $$

and will be thoroughly detailed in Section 3.4. Then, in this Section, we only consider them as intermediary quantities, providing a compact way to express the normal flux decomposition. This leads to the following decomposition of the normal flux in terms of outgoing and incoming traces:

$$ \mathcal{F}_{\mathbf{n}}^{+} \mathbf{y} = \frac{1}{2} \begin{bmatrix} \tilde{Y} \mathcal{U} & \mathcal{O}_{m_1,r} \\ \mathcal{O}_{m_2,r} & \mathcal{V} \end{bmatrix} \begin{pmatrix} \gamma_{\text{out}}(\mathbf{y}) \\ \gamma_{\text{out}}(\mathbf{y}) \end{pmatrix} \quad \text{and} \quad \mathcal{F}_{\mathbf{n}}^{-} \mathbf{y} = \frac{1}{2} \begin{bmatrix} -\tilde{Y} \mathcal{U} & \mathcal{O}_{m_1,r} \\ \mathcal{O}_{m_2,r} & \mathcal{V} \end{bmatrix} \begin{pmatrix} \gamma_{\text{in}}(\mathbf{y}) \\ \gamma_{\text{in}}(\mathbf{y}) \end{pmatrix}. \quad (12) $$

## 2.3 Application in the framework of a mesh

Being given this decomposition methodology, it can be adapted to a mesh, conformal with the parameters heterogeneity, by considering a cell $T \in \mathcal{T}_h$ and the outward unit normal vector $\mathbf{n}_T$ on its boundary.

**Remark 2.6.** We highlight the fact the compact singular value decomposition (5) of $\mathcal{B}_{\mathbf{n}_T}$ must be realised on each face $F \in \mathcal{F}(T)$ of the boundary cell $\partial T$ (where the unit outward normal vector $\mathbf{n}_T$ is constant), leading to

$$ \mathcal{B}_{\mathbf{n}_T} = \mathcal{U}_T^F (\mathcal{V}_T^F)^* \quad \text{on } F. $$
Nevertheless, so as to simplify the notations, the subscript ·F T will be presupposed without further mention.

Furthermore, the rank of $\mathcal{B}_{n_T}$, denoted by $r$, 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.

Thus, for a face $F$ of the boundary $\partial T$, the matrix $\tilde{\mathcal{M}}$ will naturally be chosen as $\tilde{\mathcal{M}} = \mathcal{M}_T$ (meaning that $\tilde{\varepsilon} = \varepsilon_T$, $\tilde{\mu} = \mu_T$, and then $\tilde{Z} = Z_T$ and $\tilde{Y} = Y_T$ in (8)), which leads to the following decomposition on $F$:

$$ \mathcal{F}_{\mathbf{n}_T}^{\pm} = \frac{1}{2} \begin{bmatrix} \pm Y_T \mathcal{U}\mathcal{U}^* & \mathcal{U}\mathcal{V}^* \\ \mathcal{V}\mathcal{U}^* & \pm Z_T \mathcal{V}\mathcal{V}^* \end{bmatrix}. $$

**Remark 2.7.** However, in case another cell $K$ shares this face $F$, such that $F = \partial T \cap \partial K$ and $n_T = -n_K$, the compact singular value decomposition operation of

$$ \mathcal{U}_T^F (\mathcal{V}_T^F)^* = \mathcal{B}_{n_T} = -\mathcal{B}_{n_K} = -\mathcal{U}_K^F (\mathcal{V}_K^F)^* $$

must be realised in each cell. This calls for distinct operators $\mathcal{U}_{T/K}^F$ and $\mathcal{V}_{T/K}^F$ in each cell, with (as an example where $\mathcal{U}$ is preserved):

$$ \mathcal{U}_T^F = \mathcal{U}_K^F, \quad \text{and} \quad \mathcal{V}_T^F = -\mathcal{V}_K^F. $$

In a simplification idea, we will then only refer to operators relative to $T$, allowing to rewrite the decomposition of $\mathcal{F}_{n_K}$ as:

$$ \mathcal{F}_{\mathbf{n}_K}^{\pm} = \frac{1}{2} \begin{bmatrix} \pm Y_K \mathcal{U}_T^F (\mathcal{U}_T^F)^* & -\mathcal{U}_T^F (\mathcal{V}_T^F)^* \\ -\mathcal{V}_T^F (\mathcal{U}_T^F)^* & \pm Z_K \mathcal{V}_T^F (\mathcal{V}_T^F)^* \end{bmatrix}. \quad (13) $$

Thus, in the rest of the Report, we propose to favour the cell $T$ (without loss of generality), as the operators $\mathcal{U}$ and $\mathcal{V}$ will systematically refer to the cell $T$ and the subscript ·F T will be presupposed without further mention (unless it is explicitly mentioned).

Nevertheless, in the case the face $F$ belongs to the boundary, we must consider an additional decomposition of the normal flux $\mathcal{F}_{n_{\partial\Omega}}$ adapted to the boundary condition $\mathcal{F}_{-n_{\partial\Omega}}\mathbf{y} = \mathbf{G}$ on $\partial\Omega$ (with $n_T = n_{\partial\Omega}$ in this case).

To do so, we introduce a specific pivotal boundary matrix

$$ \mathcal{M}_{\partial\Omega} = \begin{bmatrix} \varepsilon_{\partial\Omega} \mathcal{I}_{m_1} & \mathcal{O}_{m_1, m_2} \\ \mathcal{O}_{m_2, m_1} & \mu_{\partial\Omega} \mathcal{I}_{m_2} \end{bmatrix}, $$

with constant strictly positive values $\varepsilon_{\partial\Omega}$ and $\mu_{\partial\Omega}$ on $F$. Such boundary parameters may be associated to a physical domain lying outside of $\Omega$, or even correspond to a given impedance value which condensates the effects of the boundary.

Taking back the decomposition (9) with respect to these boundary parameters, we then obtain

$$ \mathcal{F}_{n_{\partial\Omega}}^{-} = \frac{1}{2} \begin{bmatrix} -Y_{\partial\Omega} \mathcal{U}\mathcal{U}^* & \mathcal{U}\mathcal{V}^* \\ \mathcal{V}\mathcal{U}^* & -Z_{\partial\Omega} \mathcal{V}\mathcal{V}^* \end{bmatrix}, $$

where we denote $Y_{\partial\Omega} = (\varepsilon_{\partial\Omega}/\mu_{\partial\Omega})^{1/2}$ and $Z_{\partial\Omega} = (\mu_{\partial\Omega}/\varepsilon_{\partial\Omega})^{1/2}$, which allows to rewrite the boundary condition as

$$ \mathcal{F}_{n_{\partial\Omega}}^{-} \mathbf{y} = \frac{1}{2} \begin{pmatrix} -Y_{\partial\Omega} \mathcal{U}\mathcal{U}^* \mathbf{y}_1 + \mathcal{U}\mathcal{V}^* \mathbf{y}_2 \\ \mathcal{V}\mathcal{U}^* \mathbf{y}_1 - Z_{\partial\Omega} \mathcal{V}\mathcal{V}^* \mathbf{y}_2 \end{pmatrix} = \mathbf{G} \quad \text{on } \partial\Omega. $$
### 2.3.1 Rewriting in terms of outgoing/incoming traces

In the framework of a mesh, the outgoing and incoming traces (11) are then defined as

$$ \gamma_{out}^{T}(\mathbf{y}^{T}) = \mathcal{U}^{*}\mathbf{y}_{1}^{T} + Z_{T}\mathcal{V}^{*}\mathbf{y}_{2}^{T} \quad \text{and} \quad \gamma_{in}^{T}(\mathbf{y}^{T}) = \mathcal{U}^{*}\mathbf{y}_{1}^{T} - Z_{T}\mathcal{V}^{*}\mathbf{y}_{2}^{T}, $$

where $\mathcal{U}$ and $\mathcal{V}$ are related to the compact SVD decomposition (5) of $\mathcal{B}_{n_T}$.

**Remark 2.8.** Considering two neighbouring cells T and K (separated by a face $F = \partial T \cap \partial K$) in a homogeneous domain (meaning that $Z_T = Z_K$), the traces on both sides are directly linked to each other.

Indeed, as $\mathbf{n}_T = -\mathbf{n}_K$ on F, by leaning on Remark 2.7 and specifically detailing the SVD operators for each cell, one gets on F

$$ \gamma_{out}^{T}(\mathbf{y}) = (\mathcal{U}_{T}^{F})^{*}\mathbf{y}_{1} + Z_{T}(\mathcal{V}_{T}^{F})^{*}\mathbf{y}_{2} = (\mathcal{U}_{K}^{F})^{*}\mathbf{y}_{1} + Z_{K}(-\mathcal{V}_{K}^{F})^{*}\mathbf{y}_{2} = \gamma_{in}^{K}(\mathbf{y}). $$

Thus, as one may expect, in a homogeneous configuration (where waves propagate without any reflection between the two cells), the outgoing trace from one cell is the one entering its neighbour.

Moreover, in case the face F belongs to the boundary, by relying on (12), the negative part of the normal flux can be given in terms of the incoming trace as

$$ \mathcal{F}_{\mathbf{n}_{\partial \Omega}}^{-} \mathbf{y} = \frac{1}{2} \begin{bmatrix} -Y_{\partial \Omega} & \mathcal{U} & \mathcal{O}_{m_1,r} \\ \mathcal{O}_{m_2,r} & \mathcal{V} & \end{bmatrix} \begin{pmatrix} \gamma_{in}^{\partial \Omega}(\mathbf{y}) \\ \gamma_{in}^{\partial \Omega}(\mathbf{y}) \end{pmatrix} \quad \text{where} \quad \gamma_{in}^{\partial \Omega}(\mathbf{y}) = \mathcal{U}^{*}\mathbf{y}_{1} - Z_{\partial \Omega} \mathcal{V}^{*}\mathbf{y}_{2}. \quad (14) $$

Thus, the usual boundary condition that we consider for general Friedrichs system can be condensed by only considering the condition

$$ \gamma_{in}^{\partial \Omega}(\mathbf{y}) = \mathbf{g}, \quad (15) $$

as the original one can be recovered by considering

$$ \mathbf{G} = \frac{1}{2} \begin{bmatrix} -Y_{\partial \Omega} & \mathcal{U} & \mathcal{O}_{m_1,r} \\ \mathcal{O}_{m_2,r} & \mathcal{V} & \end{bmatrix} \begin{pmatrix} \mathbf{g} \\ \mathbf{g} \end{pmatrix}. $$

**Remark 2.9.** More precisely, as its name may presuppose, imposing the value of the incoming trace $\gamma_{in}^{\partial \Omega}(\mathbf{y})$ allows to characterise the information coming in the domain or, in other words, the data going out of a fictitious cell K which would be out of $\Omega$ and share the face F with T.

Indeed, by considering $\mathcal{M}_K = \mathcal{M}_{\partial \Omega}$ and leaning on Remark 2.7 to only use operators relative to T, its outgoing trace would be defined as (considering that $\mathbf{n}_K = -\mathbf{n}_{\partial \Omega}$):

$$ \gamma_{out}^{K}(\mathbf{y}) = (\mathcal{U}_{K}^{F})^{*}\mathbf{y}_{1} + Z_{\partial \Omega} (\mathcal{V}_{K}^{F})^{*}\mathbf{y}_{2} = (\mathcal{U}_{T}^{F})^{*}\mathbf{y}_{1} - Z_{\partial \Omega} (\mathcal{V}_{T}^{F})^{*}\mathbf{y}_{2} = \gamma_{in}^{\partial \Omega}(\mathbf{y}). $$

**Remark 2.10 (Coherence with the usual BCs).** The condensed boundary condition (15) allows to generalise the usual impedance ones to the framework of BATFFS.

Indeed, by considering a boundary cell (such that $\mathbf{n} = \mathbf{n}_{\partial \Omega}$):

* for the Maxwell system, one gets:

$$ \gamma_{in}^{\partial \Omega}(\mathbf{y}) = \mathbf{g} \iff \gamma_t[\mathbf{n}_{\partial \Omega}](\mathbf{e}) + Z_{\partial \Omega} \gamma_{\times}[\mathbf{n}_{\partial \Omega}](\mathbf{h}) = \mathbf{g}, \quad (16) $$

where

$$ \gamma_t[\mathbf{n}_{\partial \Omega}](\mathbf{w}) = \mathbf{w} - (\mathbf{w} \cdot \mathbf{n}_{\partial \Omega}) \mathbf{n}_{\partial \Omega} = (\mathbf{n}_{\partial \Omega} \times \mathbf{w}) \times \mathbf{n}_{\partial \Omega} \quad \text{and} \quad \gamma_{\times}[\mathbf{n}_{\partial \Omega}](\mathbf{w}) = \mathbf{n}_{\partial \Omega} \times \mathbf{w}. \quad (17) $$
* for the 2D Helmholtz equation, one gets:

$$ \gamma_{in}^{\partial\Omega}(\mathbf{y}) = g \iff p + Z_{\partial\Omega} \mathbf{v} \cdot \mathbf{n}_{\partial\Omega} = g. $$

* for the 1D Helmholtz equation, one gets:

$$ \gamma_{in}^{\partial\Omega}(\mathbf{y}) = g \iff p + Z_{\partial\Omega} n_{\partial\Omega} v = g. $$

In particular, one will note that the term g appearing in (16) is of dimension r = 2 and d = 3 on both sides of the equivalence. Indeed, on the left side, g should be seen in terms of its projection $\tilde{g}$ over the space of eigenvectors of $B_n$ associated to non-zero eigenvalues (of dimension r). As an example, considering a face F orthogonal to $d_1$, one would have

$$ \mathbf{g} = \begin{pmatrix} 0 \\ g_2 \\ g_3 \end{pmatrix} \in L_t^2(\partial\Omega) \quad \text{and then} \quad \tilde{\mathbf{g}} = \begin{pmatrix} g_2 \\ g_3 \end{pmatrix} \in [L^2(\partial\Omega)]^r $$

as g is tangent to F.

Nevertheless, this distinction is likely to bring avoidable complexity in the notation, which leads us to stick with the same one for both g and $\tilde{g}$.

In a nutshell, the particular structure of the BATFFS flux operators allows to explicitly define a decomposition of the normal flux. In particular, it can be expressed in terms of outgoing and incoming traces, allowing to condense the boundary condition in terms of the data coming in the domain.

## 2.4 Complete BATFFS Problem

The previous developments then allow to write the generic form of a block antidiagonal two-fields Friedrichs system, that reads as:

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

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

that we equip with the boundary condition

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

for a source term $\mathbf{f} \in [L^2(\Omega; \mathbb{C})]^m$ and a boundary field $\mathbf{g} \in [L^2(\partial\Omega; \mathbb{C})]^r$.

## 2.5 A generic methodology to define consistent numerical fluxes for BATFFS

Some numerical methods usually rely on the definition of so-called numerical fluxes, allowing to generalise the definition of the trace to piecewise discontinuous functions. As the Trefftz method will also need such tools, we propose a general methodology to define consistent numerical fluxes in the framework of block antidiagonal two-fields Friedrichs systems, before specifying them for the three wave propagation problems of interest (for which they are consistent with classic expressions).
### 2.5.1 Continuity of the normal flux

First of all, we highlight an intrinsic continuity property of the normal flux for time-harmonic Friedrichs systems, which allows to recover (and prove) the expected ones for the wave propagation problems of interest.

**Proposition 2.12.** Let's consider $y \in H(\text{div}(\mathcal{F}), \Omega)$ smooth enough, such that one has

$$i\kappa \mathcal{M}y + \text{div}(\mathcal{F}y) = 0 \quad \text{a.e. in } \Omega.$$

Then, for any interface $\Gamma$ which would split $\Omega$ into two subdomains denoted as $\Omega_{\pm}$ (with their related variables associated to the superscripts $\pm$ respectively), one has:

$$\mathcal{F}_{n_+} y^+ + \mathcal{F}_{n_-} y^- = 0 \quad \text{a.e. on } \Gamma,$$

where $n_{\pm}$ stands for the unit outgoing normal to $\Omega_{\pm}$ on $\Gamma$.

*Proof.* This proof is adapted from [9].

Let's consider a test function $\psi \in [\mathcal{P}(\Omega; \mathbb{C})]^m$. Then, according to integration-by-parts Theorem (see [19, Section 3.1.2] and [22]) and the fact that $\psi = 0$ on $\partial\Omega$, one has:

$$\begin{aligned}
\int_{\Omega} \operatorname{div}(\mathcal{F}y) \cdot \bar{\psi} \,dx &= -\int_{\Omega} y \cdot \operatorname{div}(\overline{\mathcal{F}\psi}) \,dx \\
&= -\int_{\Omega_{+}} y^{+} \cdot \operatorname{div}(\overline{\mathcal{F}\psi}) \,dx - \int_{\Omega_{-}} y^{-} \cdot \operatorname{div}(\overline{\mathcal{F}\psi}) \,dx,
\end{aligned}
\qquad (20)$$

as $y \in [L^2(\Omega; \mathbb{C})]^m$. Moreover, by focusing on each side of $\Gamma$, one gets:

$$\begin{aligned}
-\int_{\Omega_{\pm}} y^{\pm} \cdot \operatorname{div}(\overline{\mathcal{F}\psi}) \,dx &= \int_{\Omega_{\pm}} \operatorname{div}(\mathcal{F}y) \cdot \bar{\psi} \,dx - \int_{\partial\Omega_{\pm}} \mathcal{F}_{n_{\pm}} y^{\pm} \cdot \bar{\psi} \,d\sigma \\
&= \int_{\Omega_{\pm}} \operatorname{div}(\mathcal{F}y) \cdot \bar{\psi} \,dx - \int_{\Gamma} \mathcal{F}_{n_{\pm}} y^{\pm} \cdot \bar{\psi} \,d\sigma,
\end{aligned}
\qquad (21)$$

as $\psi = 0$ on $\partial\Omega_{\pm}\setminus\Gamma \subset \partial\Omega$. Then, by adding the quantity

$$\int_{\Omega} i\kappa \mathcal{M}y \cdot \bar{\psi} \,dx = \int_{\Omega_{+}} i\kappa \mathcal{M}y^{+} \cdot \bar{\psi} \,dx + \int_{\Omega_{-}} i\kappa \mathcal{M}y^{-} \cdot \bar{\psi} \,dx,$$

to (20) and using (21), this leads to

$$\begin{aligned}
\int_{\Omega} (i\kappa \mathcal{M}y + \operatorname{div}(\mathcal{F}y)) \cdot \bar{\psi} \,dx &= \int_{\Omega_{+}} (i\kappa \mathcal{M}y^{+} + \operatorname{div}(\mathcal{F}y^{+})) \cdot \bar{\psi} \,dx \\
&\quad + \int_{\Omega_{-}} (i\kappa \mathcal{M}y^{-} + \operatorname{div}(\mathcal{F}y^{-})) \cdot \bar{\psi} \,dx \\
&\quad - \int_{\Gamma} (\mathcal{F}_{n_{+}} y^{+} + \mathcal{F}_{n_{-}} y^{-}) \cdot \bar{\psi} \,d\sigma.
\end{aligned}$$

Finally, as $y^{\pm}$ is solution of the Friedrichs system in $\Omega^{\pm}$, we obtain the continuities of $y$ in terms of normal flux. $\square$

Yet, in the specific case of the BATFFS, an additional simplification can be realised.
**Proposition 2.13.** Let's consider a sufficiently smooth $\mathbf{y} \in H(\mathrm{div}(\mathcal{F}), \Omega)$ such that one has

$$i\kappa \mathcal{M} \mathbf{y} + \mathrm{div}(\mathcal{F} \mathbf{y}) = \mathbf{0} \quad \text{in } \Omega,$$

where $\mathbf{y}$, $\mathcal{M}$ and $\mathcal{F}$ satisfy the conditions of a Block Antidiagonal Two-Fields Friedrichs System (see Subsection 2). Then, for any interface $\Gamma$ which would split $\Omega$ into two subdomains denoted as $\Omega_{\pm}$ (with their related variables associated to the superscripts $\pm$ respectively), one has:

$$\mathcal{U}^* \mathbf{y}_1^+ = \mathcal{U}^* \mathbf{y}_1^- \quad \text{and} \quad \mathcal{V}^* \mathbf{y}_2^+ = \mathcal{V}^* \mathbf{y}_2^- \quad \text{on } \Gamma,$$

where $\mathcal{U}$ and $\mathcal{V}$ are associated to the compact SVD (5) of one of the normal flux operators (associated to one of the outgoing normal vector $\mathbf{n}_{\pm}$ of $\Omega_{\pm}$ on $\Gamma$).

*Proof.* We take back Proposition 2.12 and consider the SVD of the operator relative to $\Omega_{+}$ for example. Then, by splitting the normal flux in the case of a BATFFS, one obtains

$$\mathcal{U}\mathcal{V}^*(\mathbf{y}_2^+ - \mathbf{y}_2^-) = \mathbf{0} \quad \text{and} \quad \mathcal{V}\mathcal{U}^*(\mathbf{y}_1^+ - \mathbf{y}_1^-) = \mathbf{0} \quad \text{on } \Gamma,$$

and by left multiplying respectively by $\Sigma^{-1}\mathcal{U}^*$ and $\Sigma^{-1}\mathcal{V}^*$ these two equalities, one gets the expected continuity relation (where we recall that $\Sigma$ is a strictly positive diagonal matrix such that $\mathcal{U}^*\mathcal{U} = \mathcal{V}^*\mathcal{V} = \Sigma$, see (5)). □

**Remark 2.14.** Taking back the classic wave propagation problems we considered previously and a unit vector $\mathbf{n} \in \mathbb{R}^d$, the continuity Proposition 2.12 reads:

* for the Maxwell problem, we have

$$\mathcal{F}_{\mathbf{n}\mathbf{y}} = \begin{bmatrix} \mathcal{O}_{3,3} & -\Gamma_{\times}[\mathbf{n}] \\ \Gamma_{\times}[\mathbf{n}] & \mathcal{O}_{3,3} \end{bmatrix} \begin{pmatrix} \mathbf{e} \\ \mathbf{h} \end{pmatrix} = \begin{pmatrix} -\gamma_{\times}[\mathbf{n}](\mathbf{h}) \\ \gamma_{\times}[\mathbf{n}](\mathbf{e}) \end{pmatrix},$$

which allows to recover the tangential continuity of both electric and magnetic fields.

* for the 2D Helmholtz problem, we have

$$\mathcal{F}_{\mathbf{n}\mathbf{y}} = \begin{bmatrix} 0 & -\mathbf{n}^* \\ -\mathbf{n} & \mathcal{O}_{2,2} \end{bmatrix} \begin{pmatrix} \mathbf{p} \\ \mathbf{v} \end{pmatrix} = - \begin{pmatrix} \mathbf{v} \cdot \mathbf{n} \\ p \mathbf{n} \end{pmatrix},$$

which allows to recover the continuity of both pressure $p$ and normal velocity $v_N$.

* for the 1D Helmholtz problem (where $n = \pm 1$), we have

$$\mathcal{F}_{\mathbf{n}\mathbf{y}} = \begin{bmatrix} 0 & -n \\ -n & 0 \end{bmatrix} \begin{pmatrix} p \\ v \end{pmatrix} = -n \begin{pmatrix} v \\ p \end{pmatrix},$$

which allows to recover the continuity of both pressure $p$ and velocity $v$.

### 2.5.2 Numerical flux computation

This continuity property should be recovered by the numerical methods chosen to approximate the solution of the Friedrichs system. Even if this can be realised by directly penalising the jumps of this quantity (as it is realised in IPDG or least-square formulations for example), it can be imposed in a weak way for methods relying on piecewise discontinuous basis functions (as DG or Trefftz) for which the values of the normal flux can mismatch on each side of an internal face.
To do so, considering the interface F between two cells T and K, the normal flux trace $\mathcal{F}_{n_T y}$ (which is then continuous) can be replaced in the variational formulations by the so-called numerical flux denoted as $\overline{\mathcal{F}_{n_T y}}$: this consistent quantity, computed from the values on each side of the mesh interface, allows to recover continuity at convergence.

The consistency property is absolutely needed, as it ensures the numerical scheme is (at least almost) satisfied by the actual solution and that its continuity is to be recovered at convergence. More precisely, the consistency property means that replacing the function $y$ by the actual solution $y_{exa}$ (satisfying then the continuity property 2.12), the numerical flux must match the actual value:

$$ \overline{\mathcal{F}_{n_T y_{exa}}} = \mathcal{F}_{n_T y_{exa}}. $$

To obtain such an expression, we decompose the normal flux operator $\mathcal{F}_{n_T}$ through a generalised eigenvalue problem defined with respect to the symmetric positive definite matrix $M_T$. This leads to the following decomposition

$$ \mathcal{F}_{n_T} = \mathcal{F}_{n_T}^{+} + \mathcal{F}_{n_T}^{-}, $$

where $\mathcal{F}_{n_T}^{\pm}$ are respectively positive and negative matrices.

In particular, $\mathcal{F}_{n_T}^{+}$ allows to select the information that goes out the cell T, while $\mathcal{F}_{n_T}^{-}$ corresponds to the quantities that come in. Intuitively, at the interface $F$, only the data going out of the two cells T and K is available. Hence, the numerical flux $\overline{\mathcal{F}_{n_T y}}$ will be searched as a combination of the outgoing fluxes from both sides, $\mathcal{F}_{n_T}^{+} y^T$ and $\mathcal{F}_{n_T}^{+} y^K$.

Such a procedure in the general case of time-harmonic Friedrichs systems in heterogeneous domains can be found in [9, Section 3.6] but leads to a quite complex expression which may be heavy to deal with. Thus, we choose to detail it in the specific case of the block antidiagonal two-fields Friedrichs systems, which encompasses the wave propagation problems of interest, in both homogeneous and heterogeneous configurations. In particular, this allows to lean on the explicit decompositions (9) and (12) of the normal flux, previously detailed.

Furthermore, we recall that the choice of a piecewise constant heterogeneity distribution, and a conformal associated mesh, lead to matrices $M$ constant in each cell T, and then referred as $M_T$.

### 2.5.3 Numerical fluxes for BATFFS in a homogeneous configuration

We first restrict to a homogeneous configuration, as a simpler case, meaning that:

$$ M_T = M_K \iff Y_T = Y_K \quad \text{and} \quad Z_T = Z_K. $$

In this case, the decompositions of the normal flux operators are identical in T and K, such that, lying on the decomposition (7) of the normal flux and the continuity Proposition 2.12, one has for an actual solution $y$:

$$ \mathcal{F}_{n_T}(y^T - y^K) = 0 \iff M_T P \Lambda P^* M_T(y^T - y^K) = 0, $$

which leads, by multiplying by $\Lambda^{-1}\text{diag}(\Sigma, \Sigma)^{-1}P^*$ and using the $M_T$-orthogonality of $P$ (see (6)), to

$$ P^* M_T(y^T - y^K) = 0 \quad \text{and then} \quad \mathcal{F}_{n_T}^{\pm}(y^T - y^K) = M_T P \Lambda^{\pm} P^* M_T(y^T - y^K) = 0. $$
Consequently, using the fact that $\mathcal{F}_{n_T}^{\pm} = \mathcal{F}_{n_K}^{\mp}$, the initial splitting $\mathcal{F}_{n_T} = \mathcal{F}_{n_T}^{+} + \mathcal{F}_{n_T}^{-}$ can then be rewritten as

$$ \mathcal{F}_{n_T} y^T = \mathcal{F}_{n_T}^{+} y^T + \mathcal{F}_{n_T}^{-} y^T = \mathcal{F}_{n_T}^{+} y^T + \mathcal{F}_{n_T}^{-} y^K = \mathcal{F}_{n_T}^{+} y^T + \mathcal{F}_{n_K}^{+} y^K. $$

Thus, leaning on this decomposition of the normal flux terms of the outgoing data from both cells, we introduce the numerical flux in a homogeneous configuration as

$$ \widehat{\mathcal{F}_{n_T} y} = \mathcal{F}_{n_T}^{+} y^T + \mathcal{F}_{n_K}^{+} y^K. $$

**Remark 2.15.** In particular, this result can also be obtained and interpreted by relying on a Riemann solver [9, Remark 3.10].

### 2.5.4 Numerical fluxes for BATFFS in a heterogeneous configuration

In the case of a heterogeneous configuration, the previous recipe cannot be followed anymore, as the normal flux is not ensured to share the same decomposition in the two neighbouring cells.

The numerical flux $\widehat{\mathcal{F}}_{n_T} y$ is then searched as a linear combination between the outgoing fluxes from both sides, $\mathcal{F}_{n_T}^{+} y^T$ and $\mathcal{F}_{n_K}^{+} y^K$, in the form:

$$ \widehat{\mathcal{F}_{n_T} y} = \begin{bmatrix} \alpha_T \mathcal{I}_{m_1} & \mathcal{O}_{m_1, m_2} \\ \mathcal{O}_{m_2, m_1} & \beta_T \mathcal{I}_{m_2} \end{bmatrix} \mathcal{F}_{n_T}^{+} y^T + \begin{bmatrix} \alpha_K \mathcal{I}_{m_1} & \mathcal{O}_{m_1, m_2} \\ \mathcal{O}_{m_2, m_1} & \beta_K \mathcal{I}_{m_2} \end{bmatrix} \mathcal{F}_{n_K}^{+} y^K, $$

or, taking back the expressions (9) and (13) of the normal flux decomposition for block antidiagonal two-fields Friedrichs systems, as

$$ \widehat{\mathcal{F}_{n_T} y} = \frac{1}{2} \begin{pmatrix} \mathcal{U} [\alpha_T Y_T \mathcal{U}^* \mathbf{y}_1^T + \alpha_T \mathcal{V}^* \mathbf{y}_2^T + \alpha_K Y_K \mathcal{U}^* \mathbf{y}_1^K - \alpha_K \mathcal{V}^* \mathbf{y}_2^K] \\ \mathcal{V} [\beta_T \mathcal{U}^* \mathbf{y}_1^T + \beta_T Z_T \mathcal{V}^* \mathbf{y}_2^T - \beta_K \mathcal{U}^* \mathbf{y}_1^K + \beta_K Z_K \mathcal{V}^* \mathbf{y}_2^K] \end{pmatrix}. \quad (22) $$

The coefficients $\alpha_T$, $\alpha_K$, $\beta_T$ and $\beta_K$ are then uniquely determined thanks to the only property that we have to ensure, consistency, which reads for an actual solution $\mathbf{y}$ (split in terms of $\mathbf{y}_1$ and $\mathbf{y}_2$) as

$$ \widehat{\mathcal{F}_{n_T} y} = \mathcal{F}_{n_T} y^T = \begin{pmatrix} \mathcal{B}_{n_T} \mathbf{y}_2^T \\ \mathcal{B}_{n_T}^* \mathbf{y}_1^T \end{pmatrix} = \begin{pmatrix} \mathcal{U} \mathcal{V}^* \mathbf{y}_2^T \\ \mathcal{V} \mathcal{U}^* \mathbf{y}_1^T \end{pmatrix}. \quad (23) $$

In particular, such a solution satisfies the continuity property 2.13: this implies that

$$ \mathcal{F}_{n_T} y^T + \mathcal{F}_{n_K} y^K = \mathbf{0} \iff \left\{ \begin{array}{l} \phi_1 := \mathcal{V}^* \mathbf{y}_2^T = \mathcal{V}^* \mathbf{y}_2^K \\ \phi_2 := \mathcal{U}^* \mathbf{y}_1^T = \mathcal{U}^* \mathbf{y}_1^K \end{array} \right. , $$

where we recall that $\mathbf{n}_K = -\mathbf{n}_T$ and we introduce the intermediary flux variables $\phi_1$ and $\phi_2$.

Then, Equation (22) can be rewritten as

$$ \widehat{\mathcal{F}_{n_T} y} = \frac{1}{2} \begin{pmatrix} \mathcal{U} [ (\alpha_T Y_T + \alpha_K Y_K) \phi_2 + (\alpha_T - \alpha_K) \phi_1 ] \\ \mathcal{V} [ (\beta_T - \beta_K) \phi_2 + (\beta_T Z_T + \beta_K Z_K) \phi_1 ] \end{pmatrix}, $$

and (23) then leads to

$$ \left\{ \begin{array}{l} 2 \phi_1 = (\alpha_T - \alpha_K) \phi_1 + (\alpha_T Y_T + \alpha_K Y_K) \phi_2 \\ 2 \phi_2 = (\beta_T Z_T + \beta_K Z_K) \phi_1 + (\beta_T - \beta_K) \phi_2 \end{array} \right. . $$
By identification, we obtain the systems

$$ \begin{cases} \alpha_T - \alpha_K = 2 \\ \alpha_T Y_T + \alpha_K Y_K = 0 \end{cases} \quad \text{and} \quad \begin{cases} \beta_T - \beta_K = 2 \\ \beta_T Z_T + \beta_K Z_K = 0 \end{cases}, $$

whose solutions can be written as

$$ \left\{ \begin{array}{l} \alpha_T = \frac{2Y_K}{Y_T + Y_K} = \frac{2Z_T}{Z_T + Z_K} \\ \\ \alpha_K = -\frac{2Y_T}{Y_T + Y_K} = -\frac{2Z_K}{Z_T + Z_K} \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} \beta_T = \frac{2Z_K}{Z_T + Z_K} \\ \\ \beta_K = -\frac{2Z_T}{Z_T + Z_K} \end{array} \right. . $$

Finally, the numerical flux can be expressed as

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \frac{1}{Z_T + Z_K} \left( \begin{array}{l} \mathcal{U} \mathcal{U}^*(\mathbf{y}_1^T - \mathbf{y}_1^K) + \mathcal{U} \mathcal{V}^*(Z_T \mathbf{y}_2^T + Z_K \mathbf{y}_2^K) \\ \mathcal{V} \mathcal{U}^*(Z_K \mathbf{y}_1^T + Z_T \mathbf{y}_1^K) + Z_T Z_K \mathcal{V} \mathcal{V}^*(\mathbf{y}_2^T - \mathbf{y}_2^K) \end{array} \right), \quad (24) $$

which can be simplified as

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \frac{1}{Z_T + Z_K} \left( \begin{array}{l} \mathcal{U} [\mathcal{U}^*(\mathbf{y}_1^T - \mathbf{y}_1^K) + \mathcal{V}^*(Z_T \mathbf{y}_2^T + Z_K \mathbf{y}_2^K)] \\ \mathcal{V} [\mathcal{U}^*(Z_K \mathbf{y}_1^T + Z_T \mathbf{y}_1^K) + Z_T Z_K \mathcal{V}^*(\mathbf{y}_2^T - \mathbf{y}_2^K)] \end{array} \right). \quad (25) $$

**Remark 2.16.** In the specific case of the wave propagation problems we study in this Report, Remark 2.5 allows to rewrite (24) as

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \frac{1}{Z_T + Z_K} \left( \begin{array}{l} \mathcal{B}_{n_T} \mathcal{B}_{n_T}^* (\mathbf{y}_1^T - \mathbf{y}_1^K) + \mathcal{B}_{n_T} (Z_T \mathbf{y}_2^T + Z_K \mathbf{y}_2^K) \\ \mathcal{B}_{n_T}^* (Z_K \mathbf{y}_1^T + Z_T \mathbf{y}_1^K) + Z_T Z_K \mathcal{B}_{n_T}^* \mathcal{B}_{n_T} (\mathbf{y}_2^T - \mathbf{y}_2^K) \end{array} \right). \quad (26) $$

Replacing the actual normal flux in this expression, the numerical flux for the wave propagation problems is defined as:

* for the Maxwell system, one gets

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \left( \frac{-\gamma_{\times}^T \mathbf{h}^T}{\gamma_{\times}^T \mathbf{e}^T} \right) = \frac{1}{Z_T + Z_K} \left( \begin{array}{l} \gamma_t (\mathbf{e}^T - \mathbf{e}^K) - \gamma_{\times}^T (Z_T \mathbf{h}^T + Z_K \mathbf{h}^K) \\ \gamma_{\times}^T (Z_K \mathbf{e}^T + Z_T \mathbf{e}^K) + Z_T Z_K \gamma_t (\mathbf{h}^T - \mathbf{h}^K) \end{array} \right), $$

where we compactly denote $\gamma_t = \gamma_t[\mathbf{n}_T]$ (without ambiguity as $\gamma_t[\mathbf{n}_T] = \gamma_t[-\mathbf{n}_T]$) and $\gamma_{\times}^T = \gamma_{\times}[\mathbf{n}_T]$, and then

$$ \left( \frac{\widehat{\gamma_t \mathbf{e}^T}}{\widehat{\gamma_{\times}^T \mathbf{h}^T}} \right) = \frac{1}{Z_T + Z_K} \left( \begin{array}{l} \gamma_t (Z_K \mathbf{e}^T + Z_T \mathbf{e}^K) - Z_T Z_K \gamma_{\times}^T (\mathbf{h}^T - \mathbf{h}^K) \\ -\gamma_t (\mathbf{e}^T - \mathbf{e}^K) + \gamma_{\times}^T (Z_T \mathbf{h}^T + Z_K \mathbf{h}^K) \end{array} \right), $$

which is coherent with the classic expressions [25].

* for the 2D Helmholtz equation, one gets

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \left( \frac{-v_N^T}{-\hat{p} \mathbf{n}_T} \right) = \frac{1}{Z_T + Z_K} \left( \begin{array}{l} (p^T - p^K) - (Z_T v_N^T - Z_K v_N^K) \\ -(Z_K p^T + Z_T p^K) \mathbf{n}_T + Z_T Z_K (v_N^T + v_N^K) \mathbf{n}_T \end{array} \right), $$

and then

$$ \left( \frac{\hat{p}}{v_N^T} \right) = \frac{1}{Z_T + Z_K} \left( \begin{array}{l} (Z_K p^T + Z_T p^K) - Z_T Z_K (v_N^T + v_N^K) \\ -(p^T - p^K) + (Z_T v_N^T - Z_K v_N^K) \end{array} \right). $$
• for the 1D Helmholtz equation, one gets

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \begin{pmatrix} -\hat{v} n_T \\ -\hat{p} n_T \end{pmatrix} = \frac{1}{Z_T + Z_K} \begin{pmatrix} (p^T - p^K) - (Z_T v^T + Z_K v^K) n_T \\ -(Z_K p^T + Z_T p^K) n_T + Z_T Z_K (v^T - v^K) \end{pmatrix}, $$

and then

$$ \begin{pmatrix} \hat{p} \\ \hat{v} \end{pmatrix} = \frac{1}{Z_T + Z_K} \begin{pmatrix} (Z_K p^T + Z_T p^K) - Z_T Z_K (v^T - v^K) n_T \\ -(p^T - p^K) n_T + (Z_T v^T + Z_K v^K) \end{pmatrix}. $$

Moreover, leaning on the expression (12) of the normal flux decomposition in terms of the outgoing and incoming traces, and (25), it can be expressed as

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \frac{1}{Z_T + Z_K} \begin{pmatrix} \mathcal{U}(\gamma_{out}^T(\mathbf{y}^T) - \gamma_{out}^K(\mathbf{y}^K)) \\ \mathcal{V}(Z_K \gamma_{out}^T(\mathbf{y}^T) + Z_T \gamma_{out}^K(\mathbf{y}^K)) \end{pmatrix}, \quad (27) $$

and by recalling the definitions (2) and (5) of $\mathcal{F}_n$ and its decomposition, one even obtains the following numerical traces

$$ \begin{pmatrix} \widehat{\mathcal{U}^* \mathbf{y}_1} \\ \widehat{\mathcal{V}^* \mathbf{y}_2} \end{pmatrix} = \frac{1}{Z_T + Z_K} \begin{pmatrix} Z_K \gamma_{out}^T(\mathbf{y}^T) + Z_T \gamma_{out}^K(\mathbf{y}^K) \\ \gamma_{out}^T(\mathbf{y}^T) - \gamma_{out}^K(\mathbf{y}^K) \end{pmatrix}. \quad (28) $$

Then, in the case $F$ is a boundary face, the cell $K$ does not have any physical meaning. Yet, interpreting it as a ghost cell on the other side of $F$, and referring to Remark 2.9, the boundary condition (19) stands as

$$ \gamma_{out}^K(\mathbf{y}) = \gamma_{in}^{\partial \Omega}(\mathbf{y}) = \mathbf{g}, $$

which leads to the following expression of the numerical flux on the boundary:

$$ \widehat{\mathcal{F}_{n_T} \mathbf{y}} = \frac{1}{Z_T + Z_{\partial \Omega}} \begin{pmatrix} \mathcal{U}(\gamma_{out}^T(\mathbf{y}^T) - \mathbf{g}) \\ \mathcal{V}(Z_{\partial \Omega} \gamma_{out}^T(\mathbf{y}^T) + Z_T \mathbf{g}) \end{pmatrix}, \quad (29) $$

and then

$$ \begin{pmatrix} \widehat{\mathcal{U}^* \mathbf{y}_1} \\ \widehat{\mathcal{V}^* \mathbf{y}_2} \end{pmatrix} = \frac{1}{Z_T + Z_{\partial \Omega}} \begin{pmatrix} Z_{\partial \Omega} \gamma_{out}^T(\mathbf{y}^T) + Z_T \mathbf{g} \\ \gamma_{out}^T(\mathbf{y}^T) - \mathbf{g} \end{pmatrix}. \quad (30) $$

## 3 Introduction to Trefftz methods for Block Antidiagonal Two-Fields Friedrichs Systems

The Trefftz approach can be interpreted as a Discontinuous Galerkin method for which the basis functions are local solutions of the actual problem. Expecting a reduced cost in dofs per wavelength while being adapted to iterative algorithms, it showed promising abilities to deal with wave propagation in large heterogeneous domains [25].

Thus, we dedicate this Section to a complete introduction of such a method in the framework of homogeneous Block Antidiagonal Two-Fields Friedrichs Systems (BATFFS), while recalling some classic results and highlighting the limits of its usual application environment.

It will mainly adapt the ideas, properties and proofs of [4], [25] and [9] to the generic framework of BATFFS.
Hence, in this Section, we are interested in BATFFS, in which framework we recall the functions $\mathbf{y} \in [L^2(\Omega; \mathbb{C})]^m$ are split in two contributions $\mathbf{y}_1$ and $\mathbf{y}_2$ of dimensions $m_1$ and $m_2 = m - m_1$ respectively. Concerning the real matrices, they follow the specific structure

$$ \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} \quad \text{and} \quad \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}, $$

while the parameters $\varepsilon$ and $\mu$ are strictly positive and piecewise constant (in each cell). Finally, considering a unit vector $\mathbf{n} \in \mathbb{R}^d$, the matrix $\mathcal{B}_n = \sum_{j=1}^d n_j \mathcal{B}_j$ admits the following compact singular value decomposition: $\mathcal{B}_n = \mathcal{U}\mathcal{V}^*$ where $\mathcal{U} \in \mathbb{R}^{m_1 \times r}$ and $\mathcal{V} \in \mathbb{R}^{m_2 \times r}$, while $r$ stands for the rank of $\mathcal{B}_n$. We refer the reader to Subsection 2 for a complete introduction.

In particular, the Trefftz method will be designed and studied in the absence of source term, leading to the following Problem.

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

$$ \mathcal{L}\mathbf{y} := i\kappa\mathcal{M}\mathbf{y} + \operatorname{div}(\mathcal{F}\mathbf{y}) = \mathbf{0} \quad \text{a.e. in } \Omega, \qquad (31) $$

that we equip with the boundary condition (see (14) and (15)):

$$ \gamma_{in}^{\partial\Omega}(\mathbf{y}) = \mathcal{U}^*\mathbf{y}_1 - Z_{\partial\Omega} \mathcal{V}^*\mathbf{y}_2 = \mathbf{g} \quad \text{a.e. on } \partial\Omega, \qquad (32) $$

for a boundary field $\mathbf{g} \in [L^2(\partial\Omega; \mathbb{C})]^r$ and a boundary impedance $Z_{\partial\Omega} : \partial\Omega \to \mathbb{R}^+$ that we suppose constant on each face of $\partial\Omega$ and non-uniformly zero.

**Remark 3.2.** Leaning on Remark 2.6, we recall the value of $r$ may change between mesh faces. Thus, so as to simplify the presentation, we suppose it is constant and strictly positive on every mesh face $F \in \mathcal{F}_h$.

Thus, this Section aims at defining a Trefftz method for these homogeneous block antidiagonal two-fields Friedrichs systems. Section 3.1 first introduces approximation spaces at stake, before three interpretations are provided, each offering a point of view highlighting a specific property of this method:

* a formulation issued from the reciprocity formula is derived in Section 3.2, while its well-posedness is analysed through a weak-coercivity argument in Section 3.3.

* then, this approach is proved to be equivalent to an Ultra-Weak formulation, whose function spaces and trace variables are detailed in Section 3.4. In particular, this formalism allows to highlight a natural contracting preconditioner in Section 3.5, making this method strongly adapted to iterative algorithms.

* finally, it is discussed within the scope of Domain Decomposition Methods in Section 3.6.

Last but not least, classic a priori error estimates are recalled in Section 3.7, with a special focus on the 2D Helmholtz equation, before classic plane waves discretisation is discussed in Section 3.8.
## 3.1 Approximation spaces

We introduce the function space on which the Trefftz method relies, which is made up of local solutions in each cell with regular traces on its boundaries.

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

$$ \mathcal{Y}^T := \{ \mathbf{y}^T \in [L^2(T; \mathbb{C})]^m : \mathcal{L} \mathbf{y}^T = \mathbf{0} \text{ in } T \text{ such that } \mathcal{F}_{n_T} \mathbf{y}^T \in [L^2(\partial T; \mathbb{C})]^m \}, $$

where we recall that $\mathcal{L}$ stands for the Friedrichs operator (see (31)). Considering the whole mesh, we then define the global space of piecewise solutions:

$$ \mathcal{Y} := \prod_{T \in \mathcal{T}_h} \mathcal{Y}^T, $$

which is equipped with the inner product (see Remark 3.5)

$$ \forall \mathbf{y}, \mathbf{y}' \in \mathcal{Y}, \langle \mathbf{y}, \mathbf{y}' \rangle_{\mathcal{Y}} := \sum_{T \in \mathcal{T}_h} \langle \mathbf{y}^T, \mathbf{y}'^T \rangle_{\mathcal{Y}^T}, \quad (33) $$

with

$$ \langle \mathbf{y}^T, \mathbf{y}'^T \rangle_{\mathcal{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}, $$

where we recall that $\mathcal{U}$ and $\mathcal{V}$ refer to the compact singular value decomposition (5) issued from $\mathcal{B}_{n_T}$.

**Remark 3.4.** We highlight the Friedrichs system naturally ensures its solution verifies

$$ \mathbf{y}^T \in H(\operatorname{div}(\mathcal{F}), T) := \{ \mathbf{y}^T \in [L^2(T, \mathbb{C})]^m : \operatorname{div}(\mathcal{F} \mathbf{y}^T) \in [L^2(T, \mathbb{C})]^m \}, \quad (34) $$

and then $\mathcal{F}_{n_T} \mathbf{y}^T \in [H^{-1/2}(\partial T; \mathbb{C})]^m$ on $\partial T$. Thus, the trace condition $\mathcal{F}_{n_T} \mathbf{y}^T \in [L^2(\partial T; \mathbb{C})]^m$ will force to consider more regular solutions.

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

$$ \|\mathbf{y}^T\|_{\mathcal{Y}^T} = 0 \implies \mathcal{F}_{n_T} \mathbf{y}^T = \mathbf{0} \text{ on } \partial T \implies \mathbf{y}^T = \mathbf{0} \text{ in } T. $$

Classically, as infinite-dimensional spaces cannot be handled by numerical solution, a finite-dimensional space $\mathcal{Y}_h$ will be needed. In particular, one usually prefers a conformal approximation of the original space $\mathcal{Y}$ in the sense that $\mathcal{Y}_h \subset \mathcal{Y}$: this may allow to conserve, at the discrete level, properties that are proved at the continuous one.

This conformal discretisation is realised by considering, in each cell $T \in \mathcal{T}_h$, a set of $N_{dof}^T$ boundary-regular local solutions $(\phi_j^T)_{j \in [1, N_{dof}^T]} \subset \mathcal{Y}^T$ and the space

$$ \mathcal{Y}_h^T := \operatorname{span}_{j \in [1, N_{dof}^T]} \{\phi_j^T\} \subset \mathcal{Y}^T. $$

From the point of view of the full mesh, this naturally allows to introduce the space

$$ \mathcal{Y}_h := \prod_{T \in \mathcal{T}_h} \mathcal{Y}_h^T, $$
<figure><img src="image_5.png" /><figcaption>Figure 1: Examples of global basis functions in 2D.</figcaption></figure>

whose $N_{dof} := \sum_{T \in \mathcal{T}_h} N_{dof}^T$ basis functions are referred as $(\phi_J)_{J \in [1, N_{dof}]} \subset \mathbb{Y}$ (see Figure 1).

In particular, any function $\mathbf{y}_h \in \mathbb{Y}_h$ can be equivalently represented by a vector $[\mathbf{y}_h] \in \mathbb{C}^{N_{dof}}$ which contains the contribution of each basis function in its definition. This allows to express $\mathbf{y}_h$ as

$$ \mathbf{y}_h = \sum_{J=1}^{N_{dof}} [\mathbf{y}_h]_J \phi_J. $$

**Remark 3.6.** We highlight the fact that each global basis function, identified by $J \in [1, N_{dof}]$, is associated to a unique cell $T_{dof2cell(J)} \in \mathcal{T}_h$ and a local dof number $j \in [1, N_{dof}^T]$ such that

$$ \phi_J = \begin{cases} \phi_j^{T_{dof2cell(J)}} & \text{in } T_{dof2cell(J)}, \\ \mathbf{0} & \text{in } \Omega \setminus T_{dof2cell(J)}. \end{cases} $$

Thus, each basis function is non-zero in a single mesh cell, leading to a reduced stencil and good sparsity properties of the associated matrices.

**Remark 3.7.** One will note that, in general, a function $\mathbf{w} \in \mathbb{Y}$ (and potentially the normal flux $\mathcal{F}_n \mathbf{w}$) is discontinuous across the faces of the mesh, in opposition to the results stated in Proposition 2.12: such a property will be recovered thanks to the introduction of numerical fluxes at interfaces.

## 3.2 Derivation of a Trefftz formulation from the reciprocity formula

The principle of the Trefftz approach is to consider basis functions which are local solutions in each cell and can be associated to a wide range of variational formulations (least-square, UWVF, Trefftz-DG...). Then, we present the derivation of a Trefftz formulation, and prove some of its properties in the case of BATFFS (even if the formulation may be introduced for time-harmonic Friedrichs systems, as long as numerical fluxes can be defined).

In this Section, we present the derivation of a formulation from the reciprocity formula, also referred as the virtual work formula, induced by the time-harmonic Friedrichs systems. It will allow to define well-posed variational formulations on this specific space.
**Proposition 3.8** (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, \quad (35)$$

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. \quad (36)$$

*Proof.* This proposition relies on the integration-by-parts formula (see [19, Section 3.1.2] and [22])

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

By using the fact that both $\mathbf{y}^T$ and $\mathbf{y}'^T$ satisfy

$$i\kappa \mathcal{M} \mathbf{y}^T + \mathrm{div}(\mathcal{F} \mathbf{y}^T) = i\kappa \mathcal{M} \mathbf{y}'^T + \mathrm{div}(\mathcal{F} \mathbf{y}'^T) = \mathbf{0} \quad \text{a.e. in } T,$$

one gets

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

Finally, the real symmetric character of $\mathcal{M}$ leads to

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

□

This naturally calls for the introduction of the following 'variational formulation'. Find $\mathbf{y} \in \mathbb{Y}$ such that for all $\mathbf{y}' \in \mathbb{Y}$, one has:

$$\sum_{T \in \mathcal{T}_h} r_T(\mathbf{y}^T, \mathbf{y}'^T) = \sum_{T \in \mathcal{T}_h} \langle \mathcal{F}_{n_T} \mathbf{y}^T, \mathbf{y}'^T \rangle_{0, \partial T} = 0.$$

Yet, such an approach is naturally ill-posed, as recovering the actual solution calls for the respect of boundary conditions and flux continuity properties between cells (see Remark 3.7), while no communication is allowed between neighbouring cells here.

Thus, we introduce the numerical flux at each face of the mesh to ensure these properties will be respected at convergence.

**Problem 3.9** (Trefftz formulation for time-harmonic Friedrichs system). Find $\mathbf{y} \in \mathbb{Y}$ such that for all $\mathbf{y}' \in \mathbb{Y}$, one has:

$$\tilde{a}(\mathbf{y}, \mathbf{y}', \mathbf{g}) := \sum_{T \in \mathcal{T}_h} \langle \overline{\mathcal{F}_{n_T} \mathbf{y}}, \mathbf{y}'^T \rangle_{0, \partial T} = 0. \quad (37)$$

Naturally, in the case of boundary faces, the components issued from the boundary conditions can be gathered as a linear form posed on $\mathbf{y}'$, which leads to the equivalent formulation

$$a(\mathbf{y}, \mathbf{y}') = \ell(\mathbf{y}'), \quad (38)$$

where $a$ is a sesquilinear form and $\ell$ is a conjugate-linear form which encompasses the boundary conditions.
In the specific case of BATFFS, the expression of the numerical flux on the boundary is explicitly known (see (29)), which allows to detail the forms a and $\ell$:

$$a(\mathbf{y}, \mathbf{y}') = \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_{\text{int}}(T)} \langle \overline{\mathcal{F}_{n_T} \mathbf{y}, \mathbf{y}'}^T \rangle_{0,F} \\ + \sum_{F \in \mathcal{F}_{\text{ext}}(T)} \frac{1}{Z_T + Z_{\partial\Omega}} \langle \mathcal{U}^* \mathbf{y}_1^T + Z_T \mathcal{V}^* \mathbf{y}_2^T, \mathcal{U}^* \mathbf{y}_1'^T \rangle_{0,F} \\ + \sum_{F \in \mathcal{F}_{\text{ext}}(T)} \frac{1}{Z_T + Z_{\partial\Omega}} \langle Z_{\partial\Omega} \mathcal{U}^* \mathbf{y}_1^T + Z_T Z_{\partial\Omega} \mathcal{V}^* \mathbf{y}_2^T, \mathcal{V}^* \mathbf{y}_2'^T \rangle_{0,F},$$

while, by recalling the boundary condition $\gamma_{in}^{\partial\Omega}(\mathbf{y}) = \mathbf{g}$, one has

$$\ell(\mathbf{y}') = \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_{\text{ext}}(T)} \frac{1}{Z_T + Z_{\partial\Omega}} \left( \langle \mathbf{g}, \mathcal{U}^* \mathbf{y}_1'^T \rangle_{0,F} - Z_T \langle \mathbf{g}, \mathcal{V}^* \mathbf{y}_2'^T \rangle_{0,F} \right).$$

The Trefftz formulation 3.9 relies on a specific characterisation of the exact solution: it is the only function that satisfies the boundary conditions and the normal flux continuity (see Proposition 2.12), while being a local solution in each cell. As the local solution characteristic is already verified thanks to the choice of basis functions, the variational formulation simply aims at imposing the continuity and BCs at interfaces.

Moreover, this Trefftz formulation can naturally be discretised by considering the finite-dimensional space $\mathcal{Y}_h$ instead of $\mathcal{Y}$.

**Problem 3.10** (Discrete Trefftz formulation for time-harmonic Friedrichs systems). Find $\mathbf{y}_h \in \mathcal{Y}_h$ such that for all $\mathbf{y}' \in \mathcal{Y}_h$, one has:

$$\tilde{a}(\mathbf{y}_h, \mathbf{y}', \mathbf{g}) = 0 \iff a(\mathbf{y}_h, \mathbf{y}') = \ell(\mathbf{y}').$$

By introducing the discretisation of the forms a and $\ell$ on $\mathcal{Y}_h$

$$\mathcal{A} = (a(\phi_J, \phi_I))_{I,J \in [1, N_{\text{dof}}]} \in \mathbb{C}^{N_{\text{dof}} \times N_{\text{dof}}} \quad \text{and} \quad [\mathbf{b}] = (\ell(\phi_J))_{J \in [1, N_{\text{dof}}]} \in \mathbb{C}^{N_{\text{dof}}},$$

the problem can equivalently be stated as

$$\mathcal{A}[\mathbf{y}_h] = [\mathbf{b}],$$

where we recall that $[\mathbf{y}_h] \in \mathbb{C}^{N_{\text{dof}}}$ is the vector that represents $\mathbf{y}_h$ in the basis of $\mathcal{Y}_h$.

## 3.3 Weak-coercivity property and well-posedness of the Trefftz formulation

A fundamental result for this formulation is that a weak-coercivity property holds, which ensures Problem 3.9 is well posed in the case of BATFFS.

**Proposition 3.11** (Weak-coercivity). We have the weak-coercivity property:

$$\forall \mathbf{y} \in \mathcal{Y}, \quad \Re[a(\mathbf{y}, \mathbf{y})] = \|\mathbf{y}\|_{DG}^2 \quad \text{with} \quad \|\mathbf{y}\|_{DG}^2 := |\mathbf{y}|_{\text{int}}^2 + |\mathbf{y}|_{\text{ext}}^2, \quad (39)$$
where $\| \cdot \|_{DG}$ is actually a norm on $\mathbb{Y}$, that we split in terms of contributions from internal and boundary interfaces as

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

and

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

*Proof.* (i) We start by proving (39), by relying on the reciprocity formula (35) and an expression of the sesquilinear form $a$ in a face-centred point of view: by gathering the terms coming from both sides of the internal faces, one gets

$$a(\mathbf{y}, \mathbf{y}) = \sum_{F \in \mathcal{F}_{\text{int}}} a_{\text{int}}^{F}(\mathbf{y}) + \sum_{F \in \mathcal{F}_{\text{ext}}} a_{\text{ext}}^{F}(\mathbf{y}),$$

with, by recalling that $\widehat{\mathcal{F}_{n_T} \mathbf{y}} = -\widehat{\mathcal{F}_{n_K} \mathbf{y}}$

$$a_{\text{int}}^{F}(\mathbf{y}) = \langle \widehat{\mathcal{F}_{n_T} \mathbf{y}}, \mathbf{y}^{T} \rangle_{0,F} + \langle \widehat{\mathcal{F}_{n_K} \mathbf{y}}, \mathbf{y}^{K} \rangle_{0,F} = \langle \widehat{\mathcal{F}_{n_T} \mathbf{y}}, \mathbf{y}^{T} - \mathbf{y}^{K} \rangle_{0,F}, \quad (40)$$

and

$$a_{\text{ext}}^{F}(\mathbf{y}) = \frac{Z_T}{Z_T + Z_{\partial\Omega}} \langle Y_T \mathcal{U}\mathcal{U}^* \mathbf{y}_1^T + \mathcal{B}_{n_{\partial\Omega}} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} \\ + \frac{Z_{\partial\Omega}}{Z_T + Z_{\partial\Omega}} \langle \mathcal{B}_{n_{\partial\Omega}}^* \mathbf{y}_1^T + Z_T \mathcal{V}\mathcal{V}^* \mathbf{y}_2^T, \mathbf{y}_2^T \rangle_{0,F}.$$

(a) Let $F \in \mathcal{F}_{\text{int}}$ be an internal face, intersection of two cells $T, K \in \mathcal{F}_h$ as $F = \partial T \cap \partial K$. Its contribution $a_{\text{int}}^F$ in the expression of the sesquilinear form $a$ (see (40)), can be split into three contributions by replacing the numerical flux by (25):

$$a_{\text{int}}^{F}(\mathbf{y}) = \frac{1}{Z_T + Z_K} (I_1 + I_2 + I_3),$$

with

$$I_1 = \langle \mathcal{U}\mathcal{U}^*(\mathbf{y}_1^T - \mathbf{y}_1^K), \mathbf{y}_1^T - \mathbf{y}_1^K \rangle_{0,F} + Z_T Z_K \langle \mathcal{V}\mathcal{V}^*(\mathbf{y}_2^T - \mathbf{y}_2^K), \mathbf{y}_2^T - \mathbf{y}_2^K \rangle_{0,F},$$

which corresponds to the term appearing in the definition of $|\mathbf{y}|_{\text{int}}^2$,

$$I_2 = Z_K \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + Z_T \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} + Z_T \langle \mathcal{B}_{n_K}^* \mathbf{y}_1^K, \mathbf{y}_2^K \rangle_{0,F} + Z_K \langle \mathcal{B}_{n_K} \mathbf{y}_2^K, \mathbf{y}_1^K \rangle_{0,F},$$

that is made up of crossed terms from the same cell, and

$$I_3 = Z_K \left( \langle \mathcal{B}_{n_T} \mathbf{y}_2^K, \mathbf{y}_1^T \rangle_{0,F} - \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^K \rangle_{0,F} \right) + Z_T \left( \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^K, \mathbf{y}_2^T \rangle_{0,F} - \langle \mathcal{B}_{n_T} \mathbf{y}_2^K, \mathbf{y}_1^K \rangle_{0,F} \right),$$

which comprises crossed terms between the neighbouring cells. In particular:

* $I_1$ is the sum of two positive terms, as $\mathcal{U}$ and $\mathcal{V}$ are real matrices.
• for $I_2$, as the inner product is a sesquilinear form, one has

$$ \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} = \overline{\langle \mathbf{y}_1^T, \mathcal{B}_{n_T} \mathbf{y}_2^T \rangle_{0,F}} = \overline{\langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F}}, \quad (41) $$

which implies that

$$ \Re \left[ Z_K \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + Z_T \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} \right] = (Z_T + Z_K) \Re \left[ \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} \right], $$

with

$$ \Re \left[ \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} \right] = \frac{1}{2} \left( \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} \right). $$

Thus, we obtain

$$ \begin{aligned} \Re[I_2] &= (Z_T + Z_K) \Re \left[ \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + \langle \mathcal{B}_{n_K}^* \mathbf{y}_1^K, \mathbf{y}_2^K \rangle_{0,F} \right] \\ &= \frac{Z_T + Z_K}{2} \left( \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} \right) \quad (42) \\ &\quad + \frac{Z_T + Z_K}{2} \left( \langle \mathcal{B}_{n_K}^* \mathbf{y}_1^K, \mathbf{y}_2^K \rangle_{0,F} + \langle \mathcal{B}_{n_K} \mathbf{y}_2^K, \mathbf{y}_1^K \rangle_{0,F} \right), \end{aligned} $$

where we recognise a part of the reciprocity formula (36).

• for $I_3$, a similar manipulation leads to

$$ \Re \left[ \langle \mathcal{B}_{n_T} \mathbf{y}_2^K, \mathbf{y}_1^T \rangle_{0,F} \right] = \Re \left[ \overline{\langle \mathbf{y}_1^T, \mathcal{B}_{n_T} \mathbf{y}_2^K \rangle_{0,F}} \right] = \Re \left[ \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^K \rangle_{0,F} \right], $$

showing that $I_3$ is purely imaginary, and then

$$ \Re[I_3] = 0. $$

(b) On the other side, by considering an external face $F \in \mathcal{F}_{ext}$, such that $F = \partial T \cap \partial \Omega$, the component $a_{ext}^F(\mathbf{y})$ can similarly be split as:

$$ a_{ext}^F(\mathbf{y}) = \frac{1}{Z_T + Z_{\partial\Omega}} (I_1 + I_2), $$

with

$$ I_1 = \langle \mathcal{U}\mathcal{U}^* \mathbf{y}_1^T, \mathbf{y}_1^T \rangle_{0,F} + Z_T Z_{\partial\Omega} \langle \mathcal{V}\mathcal{V}^* \mathbf{y}_2^T, \mathbf{y}_2^T \rangle_{0,F}, $$

which corresponds to the term appearing in the definition of $|\mathbf{y}|_{ext}^2$, and

$$ I_2 = Z_T \langle \mathcal{B}_{n_{\partial\Omega}} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} + Z_{\partial\Omega} \langle \mathcal{B}_{n_{\partial\Omega}}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F}. $$

In particular:

• $I_1$ is the sum of two positive terms, as $\mathcal{U}$ and $\mathcal{V}$ are real matrices.

• for $I_2$, the same manipulations as (41) and (42) lead to

$$ \Re[I_2] = \frac{Z_T + Z_{\partial\Omega}}{2} \left( \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} \right), $$

where we recognise the remaining part of the reciprocity formula (36).
(c) The last step consists in summing these contributions. In particular, one has:

$$ \Re[a(\mathbf{y}, \mathbf{y})] = \|\mathbf{y}\|_{DG}^2 + \Delta, $$

with the remainder

$$ \Delta = \frac{1}{2} \sum_{F \in \mathcal{F}_{\text{int}}} \left( \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F} \right. \\ \left. + \langle \mathcal{B}_{n_K}^* \mathbf{y}_1^K, \mathbf{y}_2^K \rangle_{0,F} + \langle \mathcal{B}_{n_K} \mathbf{y}_2^K, \mathbf{y}_1^K \rangle_{0,F} \right) \\ + \frac{1}{2} \sum_{F \in \mathcal{F}_{\text{ext}}} \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0,F} + \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0,F}. $$

Yet, as each face of each cell appears once and only once in this sum, it can be rearranged in the cell point of view as

$$ \Delta = \frac{1}{2} \sum_{T \in \mathcal{T}_h} \langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2^T \rangle_{0, \partial T} + \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1^T \rangle_{0, \partial T}. $$

The reciprocity formula (36) finally allows to conclude that

$$ \Delta = \frac{1}{2} \sum_{T \in \mathcal{T}_h} r_T(\mathbf{y}^T, \mathbf{y}^T) = 0. $$

(ii) It remains to prove that $\| \cdot \|_{DG}$ is a norm, by ensuring its definiteness (its positive character stems from the one of the impedances).

To do so, let's suppose we have $\mathbf{y} \in \mathcal{Y}$ such that $\|\mathbf{y}\|_{DG} = 0$. Then one has

$$ \left\{ \begin{array}{l} \forall F = \partial T \cap \partial K \in \mathcal{F}_{\text{int}}, \quad \mathcal{U}^*(\mathbf{y}_1^T - \mathbf{y}_1^K) = \mathcal{V}^*(\mathbf{y}_2^T - \mathbf{y}_2^K) = \mathbf{0}, \\ \forall F = \partial T \cap \partial \Omega \in \mathcal{F}_{\text{ext}}, \quad \mathcal{U}^* \mathbf{y}_1^T = \mathcal{V}^* \mathbf{y}_2^T = \mathbf{0}. \end{array} \right. $$

In particular, by relying on the form (2) of the normal flux for BATFFS and the decomposition (5) in terms of $\mathcal{U}$ and $\mathcal{V}$, we obtain

$$ \left\{ \begin{array}{ll} \mathcal{F}_{n_T} \mathbf{y}^T + \mathcal{F}_{n_K} \mathbf{y}^K = \mathbf{0} & \text{on } \partial T \cap \partial K, \\ \mathcal{F}_{n_{\partial \Omega}} \mathbf{y} = \mathbf{0} & \text{on } \partial \Omega. \end{array} \right. $$

Consequently, the function $\mathbf{y}$ verifies the time-harmonic Friedrichs system with homogeneous RHS and normal flux: the unique continuation Hypothesis 2.2 then ensures that $\mathbf{y} = \mathbf{0}$. $\square$

**Remark 3.12.** In the general case, without taking into account the fact that $\mathbf{y}^T$ is a local solution (and then satisfies the reciprocity formula), one has:

$$ \Re[a(\mathbf{y}, \mathbf{y})] = \|\mathbf{y}\|_{DG}^2 + \Delta \quad \text{with} \quad \Delta = \frac{1}{2} \sum_{T \in \mathcal{T}_h} r_T(\mathbf{y}^T, \mathbf{y}^T). $$

This weak-coercivity property then allows to prove the injectivity of the formulation, in both continuous and discrete frameworks.

**Corollary 3.13** (Injectivity of the formulation). The sesquilinear form $a$ is injective on $\mathcal{Y}$, as one has:

$$ \forall \mathbf{y} \in \mathcal{Y}, \quad a(\mathbf{y}, \mathbf{y}) = 0 \quad \implies \quad \mathbf{y} = \mathbf{0}. $$
Proof. According to Proposition 3.11, having $a(\mathbf{y}, \mathbf{y}) = 0$ implies that $\|\mathbf{y}\|_{DG} = 0$. Yet, as $\|\cdot\|_{DG}$ is a norm on $\mathbb{Y}$, this directly implies that $\mathbf{y} = \mathbf{0}$. $\square$

This injectivity property ensures the uniqueness of a solution to both continuous and discrete Trefftz Problems 3.9 and 3.10, due to the conformal character of the discrete approximation. In terms of existence, as the discrete space is finite-dimensional, the injectivity property directly implies the existence, and then the well-posedness of the discrete Problem 3.10. Considering the continuous one 3.9, this question is discussed in a general framework in [11, 10, 1] and supposed here with the well-posedness of the initial system (as we are looking for an approximation of such solution, it should at least exist). Moreover, it is proved in the cases of the wave propagation problems without referring to the Friedrichs formalism (see [18] for the Maxwell system and [2] for the Helmholtz equation).

## 3.4 Equivalence to an UWVF

This previous Trefftz approach can also be equivalently stated in the Ultra-Weak Variational Formulation (UWVF) framework, which was introduced by Cessenat and Després [4, 5]. The idea of this formulation is to introduce another set of variables on the local boundaries, the so-called outgoing traces, so as to condense the initial problem. In particular, as the functions we consider can be discontinuous between cells, these traces have to be defined independently on each side of an interface, in other words on the boundary of each mesh cell.

### 3.4.1 Introduction of the traces and function spaces

We begin by recalling the definition of the outgoing and incoming traces in the context of BATFFS, as they were previously introduced to define the normal flux decomposition and the numerical flux (see (12)). Indeed, they play a fundamental role in this ultra-weak formulation, as the latter is expressed in terms of these variables.

**Definition 3.14.** 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_{\mathrm{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_{\mathrm{in}}^T(\mathbf{y}^T) = \mathcal{U}^* \mathbf{y}_1^T - Z_T \mathcal{V}^* \mathbf{y}_2^T \end{array} \right. \qquad (43) $$

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

In particular, these traces correspond to the information going out and coming in $T$ respectively (see Remark 3.15).

Moreover, the expected regularity of $\mathbf{y}^T \in \mathbb{Y}^T$ ensures that

$$ \mathbf{x}^T, \mathcal{U}^T \mathbf{x}^T \in \mathcal{X}^T := [L^2(\partial T; \mathbb{C})]^r, $$

where we recall that $r$ is the rank of $\mathcal{B}_{n_T}$.

Then, we naturally introduce the space of square integrable traces on the faces of the mesh as

$$ \mathcal{X} := \prod_{T \in \mathcal{T}_h} \mathcal{X}^T = [L^2(\partial \mathcal{T}_h; \mathbb{C})]^r, $$

on which we consider the inner product

$$ \forall \mathbf{x}, \mathbf{x}' \in \mathcal{X}, \quad \langle \mathbf{x}, \mathbf{x}' \rangle_{\mathcal{X}} := \sum_{T \in \mathcal{T}_h} \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathcal{X}^T}, $$
with

$$ \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathcal{H}^T} := Y_T \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{0, \partial T}. \quad (44) $$

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

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

**Remark 3.15** (Interpretation of the outgoing/incoming traces with plane waves). The denomination of these variables can be illustrated by relying on plane waves, that will be thoroughly detailed in Section 3.8 in the context of time-harmonic Friedrichs systems. In a nutshell, such a function $\mathbf{y}_{\text{PW}}$ is a solution of the system associated to a specific unitary direction of propagation, $\mathbf{d} \in \mathbb{R}^d$, whose expression is of the form

$$ \forall \mathbf{x}, \quad \mathbf{y}_{\text{PW}}(\mathbf{x}) = \mathbf{Y}_{\text{PW}} e^{-i\kappa\lambda^{-1}\mathbf{d} \cdot \mathbf{x}} \quad \text{where} \quad \mathcal{F}_{\mathbf{d}}\mathbf{Y}_{\text{PW}} = \lambda \mathcal{M}\mathbf{Y}_{\text{PW}}, \quad (45) $$

with the couple of eigenvalue and eigenvector $(\lambda, \mathbf{Y}_{\text{PW}})$.

Then, taking back the temporal convention $e^{i\omega t}$, the associated time-dependent solution is

$$ \mathbf{y}_{\text{PW}}(\mathbf{x}, t) = \mathbf{Y}_{\text{PW}} e^{-i\kappa(\lambda^{-1}\mathbf{d} \cdot \mathbf{x} - ct)}. $$

Then, considering for example a propagation along $\mathbf{d} = \mathbf{d}_1$, an eigenvalue $\lambda > 0$ corresponds to a solution evolving in this direction (as a constant value is kept, when $t$ increases, by also increasing $x_1$).

Yet, in the case of a BATFFS where we split $\mathbf{Y}_{\text{PW}} \in \mathbb{C}^m$ in terms of $\mathbf{Y}_{\text{PW},1} \in \mathbb{C}^{m_1}$ and $\mathbf{Y}_{\text{PW},2} \in \mathbb{C}^{m_2}$, the eigenvalue problem (45) can be rewritten, thanks to the decompositions (2) and (5), as

$$ \begin{align*} \mathcal{F}_{\mathbf{d}}\mathbf{Y}_{\text{PW}} = \lambda \mathcal{M}\mathbf{Y}_{\text{PW}} &\iff \begin{cases} \mathcal{U}\mathcal{V}^*\mathbf{Y}_{\text{PW},2} = \lambda\varepsilon \mathbf{Y}_{\text{PW},1}, \\ \mathcal{V}\mathcal{U}^*\mathbf{Y}_{\text{PW},1} = \lambda\mu \mathbf{Y}_{\text{PW},2}, \end{cases} \\ &\iff \begin{cases} \mathbf{\Sigma}\mathcal{V}^*\mathbf{Y}_{\text{PW},2} = \lambda\varepsilon \mathcal{U}^*\mathbf{Y}_{\text{PW},1}, \\ \mathbf{\Sigma}\mathcal{U}^*\mathbf{Y}_{\text{PW},1} = \lambda\mu \mathcal{V}^*\mathbf{Y}_{\text{PW},2}, \end{cases} \end{align*} $$

by recalling that $\mathcal{U}^*\mathcal{U} = \mathcal{V}^*\mathcal{V} = \mathbf{\Sigma}$.

In particular, this leads to

$$ \mathcal{U}^*\mathbf{Y}_{\text{PW},1} = (\lambda\varepsilon)^{-1} \mathbf{\Sigma}\mathcal{V}^*\mathbf{Y}_{\text{PW},2} \quad \text{and} \quad \mathbf{\Sigma}^2\mathcal{V}^*\mathbf{Y}_{\text{PW},2} = \left(\frac{\lambda}{c}\right)^2 \mathcal{V}^*\mathbf{Y}_{\text{PW},2}, $$

and by using the diagonal positive definite character of $\mathbf{\Sigma}$ and the fact that $\lambda > 0$, one even gets $\mathbf{\Sigma}\mathcal{V}^*\mathbf{Y}_{\text{PW},2} = \frac{\lambda}{c}\mathcal{V}^*\mathbf{Y}_{\text{PW},2}$. Finally, one has

$$ \begin{align*} \gamma_{\text{in}}(\mathbf{y}_{\text{PW}}) &= (\mathcal{U}^*\mathbf{Y}_{\text{PW},1} - Z \mathcal{V}^*\mathbf{Y}_{\text{PW},2}) e^{-i\kappa\lambda^{-1}\mathbf{d} \cdot \mathbf{x}} \\ &= ((\lambda\varepsilon)^{-1} \mathbf{\Sigma}\mathcal{V}^*\mathbf{Y}_{\text{PW},2} - (\mu/\varepsilon)^{1/2} \mathcal{V}^*\mathbf{Y}_{\text{PW},2}) e^{-i\kappa\lambda^{-1}\mathbf{d} \cdot \mathbf{x}} \\ &= (\lambda\varepsilon)^{-1} e^{-i\kappa\lambda^{-1}\mathbf{d} \cdot \mathbf{x}} \left(\mathbf{\Sigma} - \frac{\lambda}{c} \mathcal{I}_r\right) \mathcal{V}^*\mathbf{Y}_{\text{PW},2} = \mathbf{0}. \end{align*} $$

This illustrates the fact, for a wave going to the right direction through the right boundary of a domain, the incoming trace of this function is zero, as all of the information goes out.
In particular, deriving this point of view through the notion of outgoing and incoming traces
leans on the well-posedness of a local problem, that we suppose in the rest of the Report (which
is verified for the wave propagation problems we are interested in).

**Hypothesis 3.16.** 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{T}^T$ (respectively $\mathcal{L}^T \mathbf{x}^T \in \mathcal{L}^T$), there exists a unique $S_{\text{out}}^T(\mathbf{x}^T) \in \mathcal{Y}^T$ (respectively $S_{\text{in}}^T(\mathcal{L}^T \mathbf{x}^T) \in \mathcal{Y}^T$) such that $\gamma_{\text{out}}^T(S_{\text{out}}^T(\mathbf{x}^T)) = \mathbf{x}^T$ (respectively $\gamma_{\text{in}}^T(S_{\text{in}}^T(\mathcal{L}^T \mathbf{x}^T)) = \mathcal{L}^T \mathbf{x}^T$) on $\partial T$.

Consequently, both traces are then directly linked as

$$ \mathcal{L}^T \mathbf{x}^T = \gamma_{\text{in}}^T(S_{\text{out}}^T(\mathbf{x}^T)) \quad \text{and} \quad \mathbf{x}^T = \gamma_{\text{out}}^T(S_{\text{in}}^T(\mathcal{L}^T \mathbf{x}^T)). $$

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

* for the Maxwell system, one gets

$$ \mathbf{x}^T = \gamma_t(\mathbf{e}) - Z_T \gamma_{\times}^T(\mathbf{h}) \quad \text{and} \quad \mathcal{L}^T \mathbf{x}^T = \gamma_t(\mathbf{e}) + Z_T \gamma_{\times}^T(\mathbf{h}), $$

which corresponds to classic expressions [4, 25].

* for the 2D Helmholtz equation, one gets

$$ \mathbf{x}^T = p^T - Z_T v_N^T \quad \text{and} \quad \mathcal{L}^T \mathbf{x}^T = p^T + Z_T v_N^T, $$

which corresponds to classic expressions [4, 5].

* for the 1D Helmholtz equation, one gets

$$ \mathbf{x}^T = p^T - n_T Z_T v^T \quad \text{and} \quad \mathcal{L}^T \mathbf{x}^T = p^T + n_T Z_T v^T. $$

**Remark 3.18.** In particular, as the outgoing and incoming traces are square integrable on $\partial T$, the inner product $\langle \mathbf{x}^T, \mathbf{x}^{T'} \rangle_{\mathcal{T}^T}$ can be decomposed as

$$ \langle \mathbf{x}^T, \mathbf{x}^{T'} \rangle_{\mathcal{T}^T} = \sum_{F \in \mathcal{F}(T)} \langle \mathbf{x}^T, \mathbf{x}^{T'} \rangle_{\mathcal{T}^T, F}, $$

with the restriction of the sesquilinear form to each face $F \in \mathcal{F}(T)$ that we naturally define as

$$ \langle \mathbf{x}^T, \mathbf{x}^{T'} \rangle_{\mathcal{T}^T, F} := Y_T \langle \mathbf{x}^T, \mathbf{x}^{T'} \rangle_{0, F}. $$

Moreover, in the case $F$ is an internal face separating the cells $T$ and $K$, the associated restricted sesquilinear forms can be related to each other.

Indeed, let's consider $\mathbf{x}^F, \mathbf{x}^{F'} \in [L^2(F; C)]^T$. Then, we obtain

$$ \begin{align*} \langle \mathbf{x}^F, \mathbf{x}^{F'} \rangle_{\mathcal{T}^T, F} &= Y_T \langle \mathbf{x}^F, \mathbf{x}^{F'} \rangle_{0, F} \\ &= Y_T Z_K \langle Y_K \langle \mathbf{x}^F, \mathbf{x}^{F'} \rangle_{0, F} \rangle_{\mathcal{T}^K, F} \\ &= Y_T Z_K \langle \mathbf{x}^F, \mathbf{x}^{F'} \rangle_{\mathcal{T}^K, F}. \end{align*} $$
Furthermore, a conformal discretisation of this trace space also has to be considered, with a view to handling numerical solution.

It is realised by directly taking the outgoing traces associated to the basis functions of the discrete space $\mathbb{Y}_h$ of boundary-regular local solutions. More precisely, in each cell $T \in \mathcal{T}_h$, we define the outgoing traces associated to the local basis functions of $\mathbb{Y}_h$

$$ \forall j \in [1, N_{dof}^T], \psi_j^T := \gamma_{out}^T(\phi_j^T), $$

and the space

$$ \mathcal{T}_h^T := \gamma_{out}^T (\mathbb{Y}_h^T) = \operatorname{span}_{j \in [1, N_{dof}^T]} \{\psi_j^T\} \subset \mathcal{T}^T, $$

which has the same dimension as $\mathbb{Y}_h^T$, due to the bijective character of the operator $\gamma_{out}^T$ (see Hypothesis 3.16).

From the point of view of the full mesh, this naturally allows to introduce the space

$$ \mathcal{T}_h := \prod_{T \in \mathcal{T}_h} \mathcal{T}_h^T, $$

whose $N_{dof}$ basis functions are referred as $(\psi_j)_{j \in [1, N_{dof}]} \subset \mathcal{T}$.

In particular, any outgoing trace $\mathbf{x}_h \in \mathbb{Y}_h$ can be equivalently represented by a vector $[\mathbf{x}_h] \in \mathbb{C}^{N_{dof}}$ which contains the contribution of each basis function in its definition. This allows to express $\mathbf{x}_h$ as

$$ \mathbf{x}_h = \sum_{J=1}^{N_{dof}} [\mathbf{x}_h]_J \psi_J. $$

**Remark 3.19.** In a similar way as we detailed in Remark 3.6, each basis function of $\mathcal{T}_h$ is non-zero on the boundary of a unique mesh cell. Indeed, considering a global basis function, identified by $J \in [1, N_{dof}]$, it is associated to a unique cell $T_{dof2cell(J)} \in \mathcal{T}_h$ and a local dof number $j \in [1, N_{dof}^{T_{dof2cell(J)}}]$ such that

$$ \psi_J = \begin{cases} \psi_j^{T_{dof2cell(J)}} & \text{on } \partial T_{dof2cell(J)}, \\ \mathbf{0} & \text{on } \partial \mathcal{T}_h \setminus \partial T_{dof2cell(J)}. \end{cases} $$

### 3.4.2 An equivalent Ultra Weak Variational Formulation

By relying on the definition of the incoming trace and the expressions (28) and (30) of the numerical flux in terms of traces, a numerical one $\widehat{\mathcal{U}^T \mathbf{x}}$ can be introduced.

**Definition 3.20.** In accordance with the definition (43) of the incoming trace, we introduce the numerical trace of the incoming trace on each face $F \in \mathcal{F}(T)$ relative to a cell $T \in \mathcal{T}_h$:

$$ \widehat{\mathcal{U}^T \mathbf{x}} = \widehat{\mathcal{U}^* \mathbf{y}_1} - Z_T \widehat{\mathcal{V}^* \mathbf{y}_2} $$

by relying on the numerical traces (28) and (30), which can be synthesised as

$$ \widehat{\mathcal{U}^T \mathbf{x}} = \begin{cases} \frac{Z_K - Z_T}{Z_T + Z_K} \mathbf{x}^T + \frac{2Z_T}{Z_T + Z_K} \mathbf{x}^K & \text{if } F = \partial T \cap \partial K \in \mathcal{F}_{\text{int}}, \\ \frac{Z_{\partial\Omega} - Z_T}{Z_T + Z_{\partial\Omega}} \mathbf{x}^T + \frac{2Z_T}{Z_T + Z_{\partial\Omega}} \mathbf{g} & \text{if } F = \partial T \cap \partial \Omega \in \mathcal{F}_{\text{ext}}. \end{cases} \quad (46) $$
**Remark 3.21.** The definition (46) of the numerical incoming trace can be interpreted in terms of reflection and transmission. Indeed, the expression presents similarities to the classic reflection coefficients at normal incidence, which allows to assimilate the incoming information as a weighted mean between the reflected outgoing one (which was leaving the cell) and the transmitted one (coming from the neighbouring cell).

This finally allows to introduce the following UWVF [4, 5]:

**Problem 3.22** (UWVF associated to 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 (47) $$

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 (48) $$

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 (49) $$

with

$$ 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 (50) $$

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} $$

Thanks to the Riesz representation theorem, both operators $k_{UWVF}$ and $\ell_{UWVF}$ can be defined as

$$ k_{UWVF}(\mathbf{x}, \mathbf{x}') := \langle \Pi_{\mathfrak{u}}\mathbf{x}, \mathfrak{u}\mathbf{x}' \rangle_{\mathcal{X}} = \langle \mathcal{K}\mathbf{x}, \mathbf{x}' \rangle_{\mathcal{X}} \quad \text{and} \quad \ell_{UWVF}(\mathbf{x}') = \langle \mathbf{g}_{UWVF}, \mathbf{x}' \rangle_{\mathcal{X}}. \quad (51) $$

Moreover, this UWVF formulation can be discretised by considering the finite-dimensional space $\mathcal{X}_h$ instead of $\mathcal{X}$.

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

$$ \tilde{a}_{UWVF}(\mathbf{x}_h, \mathbf{x}', \mathbf{g}) = 0 \iff a_{UWVF}(\mathbf{x}_h, \mathbf{x}') = \ell_{UWVF}(\mathbf{x}'). \quad (52) $$

By introducing the discretisation of the forms $a_{UWVF}$ and $\ell_{UWVF}$ on $\mathcal{X}_h$

$$ \mathcal{A}_{UWVF} = (a_{UWVF}(\psi_J, \psi_I))_{I,J \in [1, N_{\text{dof}}]} \in \mathbb{C}^{N_{\text{dof}} \times N_{\text{dof}}}, $$

and

$$ [\mathbf{b}_{UWVF}] = (\ell_{UWVF}(\psi_J))_{J \in [1, N_{\text{dof}}]} \in \mathbb{C}^{N_{\text{dof}}}, $$

the problem can equivalently be stated as

$$ \mathcal{A}_{UWVF}[\mathbf{x}_h] = [\mathbf{b}_{UWVF}], \quad (53) $$

where we recall that $[\mathbf{x}_h] \in \mathbb{C}^{N_{\text{dof}}}$ is the vector that represents $\mathbf{x}_h$ in the basis of $\mathcal{X}_h$.
In particular, this Ultra-Weak formulation is equivalent to the previous one, even though it offers another point of view, adapted to define a good preconditioner in Section 3.5.

**Proposition 3.24.** Both formulations (37) and (47) are equivalent, in the sense that for all $\mathbf{y}, \mathbf{y}' \in \mathbb{Y}$, by referring to their outgoing traces as $\mathbf{x} = \gamma_{out}(\mathbf{y})$ and $\mathbf{x}' = \gamma_{out}(\mathbf{y}')$, one has:

$$ \tilde{a}_{UWVF}(\mathbf{x}, \mathbf{x}', \mathbf{g}) = 2 \tilde{a}(\mathbf{y}, \mathbf{y}', \mathbf{g}). $$

*Proof.* Let $\mathbf{y}, \mathbf{y}' \in \mathbb{Y}$, and their associated outgoing traces $\mathbf{x} = \gamma_{out}(\mathbf{y}), \mathbf{x}' = \gamma_{out}(\mathbf{y}') \in \mathcal{X}$. Then, by referring to the definition (44) of the inner product on $\mathcal{X}^T$, one has:

$$ \begin{aligned} \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{\mathcal{U}^T \mathbf{x}^T}, \mathcal{U}^T \mathbf{x}'^T \rangle_{\mathcal{X}^T} \\ &= \sum_{T \in \mathcal{T}_h} Y_T \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{0, \partial T} - Y_T \langle \widehat{\mathcal{U}^T \mathbf{x}^T}, \mathcal{U}^T \mathbf{x}'^T \rangle_{0, \partial T}. \end{aligned} $$

Let's consider a cell $T \in \mathcal{T}_h$ and the two intermediate values

$$ I_1^T := \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{0, \partial T} \quad \text{and} \quad I_2^T := \langle \widehat{\mathcal{U}^T \mathbf{x}^T}, \mathcal{U}^T \mathbf{x}'^T \rangle_{0, \partial T}. $$

(i) First, the definition (43) of the outgoing traces leads to

$$ \begin{aligned} I_1^T &= \langle \mathbf{x}^T, \mathcal{U}^* \mathbf{y}_1^{T'} + Z_T \mathcal{V}^* \mathbf{y}_2^{T'} \rangle_{0, \partial T} \\ &= \frac{1}{Z_T + Z_K} \left( \langle (Z_T + Z_K) \mathcal{U} \mathbf{x}^T, \mathbf{y}_1^{T'} \rangle_{0, \partial T} + Z_T \langle (Z_T + Z_K) \mathcal{V} \mathbf{x}^T, \mathbf{y}_2^{T'} \rangle_{0, \partial T} \right). \end{aligned} $$

(ii) Considering the second term and the definition (46) of $\widehat{\mathcal{U}^T \mathbf{x}}$, it simplifies as

$$ \begin{aligned} I_2^T &= \langle \widehat{\mathcal{U}^T \mathbf{x}}, \mathcal{U}^T \mathbf{x}'^T \rangle_{0, \partial T} \\ &= \left\langle \frac{Z_K - Z_T}{Z_T + Z_K} \mathbf{x}^T + \frac{2Z_T}{Z_T + Z_K} \mathbf{x}^K, \mathcal{U}^* \mathbf{y}_1^{T'} - Z_T \mathcal{V}^* \mathbf{y}_2^{T'} \right\rangle_{0, \partial T} \\ &= \frac{1}{Z_T + Z_K} \left\langle (Z_K - Z_T) \mathcal{U} \mathbf{x}^T + 2Z_T \mathcal{U} \mathbf{x}^K, \mathbf{y}_1^{T'} \right\rangle_{0, \partial T} \\ &\quad - \frac{Z_T}{Z_T + Z_K} \left\langle (Z_K - Z_T) \mathcal{V} \mathbf{x}^T + 2Z_T \mathcal{V} \mathbf{x}^K, \mathbf{y}_2^{T'} \right\rangle_{0, \partial T}, \end{aligned} $$

where the cell $K$ (and its associated quantities) either corresponds to an actual cell, or a fictitious one out of the domain (such that $Z_K = Z_{\partial \Omega}$ and $\mathbf{x}^K = \mathbf{g}$): this allows to keep a unique and generic expression, as internal and boundary faces are treated in the same way.

(iii) Then by subtracting these two previous expressions, one obtains

$$ I_1^T - I_2^T = \frac{2Z_T}{Z_T + Z_K} \left( \langle \mathcal{U}(\mathbf{x}^T - \mathbf{x}^K), \mathbf{y}_1^{T'} \rangle_{0, \partial T} + \langle \mathcal{V}(Z_K \mathbf{x}^T + Z_T \mathbf{x}^K), \mathbf{y}_2^{T'} \rangle_{0, \partial T} \right). $$

Moreover, the definitions (27) and (29) of the numerical flux $\widehat{\mathcal{F}}_{n_T} \mathbf{y}$ lead to

$$ Y_T (I_1^T - I_2^T) = 2 \langle \widehat{\mathcal{F}}_{n_T} \mathbf{y}, \mathbf{y}'^T \rangle_{0, \partial T}. $$
We conclude by recalling that

$$ \tilde{a}_{UWVF}(\mathbf{y}, \mathbf{y}', \mathbf{g}) = \sum_{T \in \mathcal{T}_h} Y_T(I_1^T - I_2^T) \quad \text{and} \quad \tilde{a}(\mathbf{y}, \mathbf{y}', \mathbf{g}) = \sum_{T \in \mathcal{T}_h} \langle \widehat{\mathcal{F}_{n_T} \mathbf{y}}, \mathbf{y}'^T \rangle_{0, \partial T}. $$

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

In particular, this equivalence property does not use the fact that the functions are actual local solutions of the Friedrichs system, but only leans on the chosen variational formulation and the definition of the traces.

## 3.5 The natural contracting preconditioner of Cessenat-Després

The equivalent Ultra-Weak formulation we derived in Section 3.4 allows to highlight a natural preconditioner, first introduced by Cessenat and Després [4], which leads to a contracting problem. This justifies this Trefftz approach is particularly adapted to iterative solvers.

### 3.5.1 Contraction properties of the operators

We begin with a characterisation of the operators $\mathcal{U}$ and $\Pi_{\mathcal{U}}$.

**Proposition 3.25.** The operator $\mathcal{U}$ is an isometry and $\Pi_{\mathcal{U}}$ is a contraction, as one has

$$ \|\mathcal{U}\|_{\mathcal{X}} = 1 \quad \text{and} \quad \|\Pi_{\mathcal{U}}\|_{\mathcal{X}} \le 1. \qquad (54) $$

*Proof.* This proof is adapted from [5, Proposition 1.10] and [25, Proposition 3.1], and relies on two main steps.

(i) We first prove the operator $\mathcal{U}^T$ is an isometry. Indeed, for any cell $T \in \mathcal{T}_h$, one has:

$$ \begin{aligned} \langle \mathcal{U}^T \mathbf{x}^T, \mathcal{U}^T \mathbf{x}'^T \rangle_{\mathcal{X}^T} &= Y_T \langle \mathcal{U}^T \mathbf{x}^T, \mathcal{U}^T \mathbf{x}'^T \rangle_{0, \partial T} \\ &= Y_T \langle \mathcal{U}^* \mathbf{y}_1^T - Z_T \mathcal{V}^* \mathbf{y}_2^T, \mathcal{U}^* \mathbf{y}_1'^T - Z_T \mathcal{V}^* \mathbf{y}_2'^T \rangle_{0, \partial T}. \end{aligned} $$

Then, by recalling that

$$ \mathcal{B}_{n_T} = \mathcal{U}\mathcal{V}^* \quad \text{and} \quad \mathcal{B}_{n_T}^* = \mathcal{V}\mathcal{U}^*, $$

we can develop this expression as

$$ \begin{aligned} \langle \mathcal{U}^T \mathbf{x}^T, \mathcal{U}^T \mathbf{x}'^T \rangle_{\mathcal{X}^T} &= Y_T \langle \mathcal{U}\mathcal{U}^* \mathbf{y}_1^T, \mathbf{y}_1'^T \rangle_{0, \partial T} + Z_T \langle \mathcal{V}\mathcal{V}^* \mathbf{y}_2^T, \mathbf{y}_2'^T \rangle_{0, \partial T} \\ &\quad - \underbrace{\langle \mathcal{B}_{n_T}^* \mathbf{y}_1^T, \mathbf{y}_2'^T \rangle_{0, \partial T} + \langle \mathcal{B}_{n_T} \mathbf{y}_2^T, \mathbf{y}_1'^T \rangle_{0, \partial T}}_{r_T(\mathbf{y}^T, \mathbf{y}'^T)}. \end{aligned} \qquad (55) $$

On the other hand, similar manipulations lead to

$$ \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathcal{X}^T} = Y_T \langle \mathcal{U}\mathcal{U}^* \mathbf{y}_1^T, \mathbf{y}_1'^T \rangle_{0, \partial T} + Z_T \langle \mathcal{V}\mathcal{V}^* \mathbf{y}_2^T, \mathbf{y}_2'^T \rangle_{0, \partial T} + r_T(\mathbf{y}^T, \mathbf{y}'^T), \qquad (56) $$

and then

$$ \langle \mathcal{U}^T \mathbf{x}^T, \mathcal{U}^T \mathbf{x}'^T \rangle_{\mathcal{X}^T} = \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathcal{X}^T} - 2 r_T(\mathbf{y}^T, \mathbf{y}'^T). $$

The reciprocity formula (36) then ensures that $r_T(\mathbf{y}^T, \mathbf{y}'^T) = 0$. By summing on all the cells, we obtain

$$ \|\mathcal{U}\mathbf{x}\|_{\mathcal{X}}^2 = \sum_{T \in \mathcal{T}_h} \langle \mathcal{U}\mathbf{x}^T, \mathcal{U}\mathbf{x}^T \rangle_{\mathcal{X}^T} = \sum_{T \in \mathcal{T}_h} \langle \mathbf{x}^T, \mathbf{x}^T \rangle_{\mathcal{X}^T} = \|\mathbf{x}\|_{\mathcal{X}}^2, $$
and then $||\mathfrak{U}||_{\mathfrak{X}} = 1$.

(ii) Then, we prove the operator $\Pi_{\mathfrak{U}}$ is of norm smaller than 1. Indeed, let's split the contributions between internal and boundary faces:

$$ \begin{aligned} \|\Pi_{\mathfrak{U}}^T \mathbf{x}\|_{\mathfrak{X}}^2 &= \sum_{T \in \mathcal{F}_{\mathfrak{h}}} \sum_{F \in \mathcal{F}(T)} \langle \Pi_{\mathfrak{U}}^T \mathbf{x}^T, \Pi_{\mathfrak{U}}^T \mathbf{x}^T \rangle_{\mathfrak{X}^T, F} \\ &= \sum_{T \in \mathcal{F}_{\mathfrak{h}}} \sum_{F \in \mathcal{F}_{\text{int}}(T)} \underbrace{\|\mathfrak{U}^T \mathbf{x}\|_{\mathfrak{X}^T, F}^2}_{I_{\text{int}}^{T,F}} + \sum_{F \in \mathcal{F}_{\text{ext}}(T)} \underbrace{\|\Pi_{\mathfrak{U}}^T \mathbf{x}^T\|_{\mathfrak{X}^T, F}^2}_{I_{\text{ext}}^{T,F}}. \end{aligned} $$

Moreover, as each internal face appears twice, their contributions can be gathered such that

$$ \|\Pi_{\mathfrak{U}}^T \mathbf{x}\|_{\mathfrak{X}}^2 = \sum_{F \in \mathcal{F}_{\text{int}}} (I_{\text{int}}^{T,F} + I_{\text{int}}^{K,F}) + \sum_{F \in \mathcal{F}_{\text{ext}}} I_{\text{ext}}^{T,F}. \quad (57) $$

(a) We first consider an internal face $F = \partial T \cap \partial K \in \mathcal{F}_{\text{int}}$, for which we have

$$ \begin{aligned} I_{\text{int}}^{T,F} &= \frac{1}{(Z_T + Z_K)^2} \|(Z_K - Z_T) \mathbf{x}^T + 2Z_T \mathbf{x}^K\|_{\mathfrak{X}^T, F}^2 \\ &= \frac{(Z_T - Z_K)^2 \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 + 4Z_T^2 \|\mathbf{x}^K\|_{\mathfrak{X}^T, F}^2 + 4(Z_K - Z_T) \Re[Z_T \langle \mathbf{x}^T, \mathbf{x}^K \rangle_{\mathfrak{X}^T, F}]}{(Z_T + Z_K)^2}, \end{aligned} $$

and similarly, as $Z_K$ is real and $\Re[Z_K \langle \mathbf{x}^K, \mathbf{x}^T \rangle_{\mathfrak{X}^K, F}] = \Re[Z_K \langle \mathbf{x}^T, \mathbf{x}^K \rangle_{\mathfrak{X}^K, F}]$, one gets

$$ I_{\text{int}}^{K,F} = \frac{(Z_K - Z_T)^2 \|\mathbf{x}^K\|_{\mathfrak{X}^K, F}^2 + 4Z_K^2 \|\mathbf{x}^T\|_{\mathfrak{X}^K, F}^2 - 4(Z_K - Z_T) \Re[Z_K \langle \mathbf{x}^T, \mathbf{x}^K \rangle_{\mathfrak{X}^K, F}]}{(Z_T + Z_K)^2}. $$

Yet, relying on Remark 3.18, the two previous real parts are equal

$$ Z_T \langle \mathbf{x}^T, \mathbf{x}^K \rangle_{\mathfrak{X}^T, F} = Z_K \langle \mathbf{x}^T, \mathbf{x}^K \rangle_{\mathfrak{X}^K, F}, $$

and the inner products can be exchanged as

$$ \|\mathbf{x}^T\|_{\mathfrak{X}^K, F}^2 = Y_K Z_T \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 \quad \text{and} \quad \|\mathbf{x}^K\|_{\mathfrak{X}^T, F}^2 = Y_T Z_K \|\mathbf{x}^K\|_{\mathfrak{X}^K, F}^2, $$

which ensures that

$$ \begin{aligned} I_{\text{int}}^{T,F} + I_{\text{int}}^{K,F} &= \frac{(Z_T - Z_K)^2 + 4Z_T Z_K}{(Z_T + Z_K)^2} \left( \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 + \|\mathbf{x}^K\|_{\mathfrak{X}^K, F}^2 \right) \quad (58) \\ &= \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 + \|\mathbf{x}^K\|_{\mathfrak{X}^K, F}^2. \end{aligned} $$

(b) Then, we consider a boundary face $F = \partial T \cap \partial \Omega$. One has, thanks to the real positive character of the impedances:

$$ I_{\text{ext}}^{T,F} = \underbrace{\left( \frac{Z_{\partial \Omega} - Z_T}{Z_T + Z_{\partial \Omega}} \right)^2}_{\le 1} \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 \le \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2. \quad (59) $$

(c) Finally, by gathering (57), (58) and (59), we obtain

$$ \|\Pi_{\mathfrak{U}} \mathbf{x}\|_{\mathfrak{X}}^2 = \sum_{F \in \mathcal{F}_{\text{int}}} \left( \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 + \|\mathbf{x}^K\|_{\mathfrak{X}^K, F}^2 \right) + \sum_{F \in \mathcal{F}_{\text{ext}}} \left( \frac{Z_{\partial \Omega} - Z_T}{Z_T + Z_{\partial \Omega}} \right)^2 \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2, \quad (60) $$
and reorganising the sum with a cell point of view, we get

$$ \|\Pi_{\mathfrak{U}} \mathbf{x}\|_{\mathfrak{X}}^2 \le \sum_{F \in \mathcal{F}_{\text{int}}} \left( \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 + \|\mathbf{x}^K\|_{\mathfrak{X}^K, F}^2 \right) + \sum_{F \in \mathcal{F}_{\text{ext}}} \|\mathbf{x}^T\|_{\mathfrak{X}^T, F}^2 = \sum_{T \in \mathcal{T}_{\mathfrak{h}}} \|\mathbf{x}^T\|_{\mathfrak{X}^T}^2 = \|\mathbf{x}\|_{\mathfrak{X}}^2, \quad (61) $$

and then $\|\Pi_{\mathfrak{U}}\|_{\mathfrak{X}} \le 1$. $\square$

**Remark 3.26.** We highlight the fact that, without taking into account the fact that $\mathbf{y}^T$ satisfies the reciprocity formula, and lying on (56), one has:

$$ \left\{ \begin{array}{l} \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathfrak{X}^T} = \langle \mathbf{y}^T, \mathbf{y}'^T \rangle_{\mathfrak{Y}^T} + r_T(\mathbf{y}^T, \mathbf{y}'^T), \\ \langle \mathfrak{U}^T \mathbf{x}^T, \mathfrak{U}^T \mathbf{x}'^T \rangle_{\mathfrak{X}^T} = \langle \mathbf{y}^T, \mathbf{y}'^T \rangle_{\mathfrak{Y}^T} - r_T(\mathbf{y}^T, \mathbf{y}'^T), \end{array} \right. \quad (62) $$

and then

$$ \langle \mathfrak{U}^T \mathbf{x}^T, \mathfrak{U}^T \mathbf{x}'^T \rangle_{\mathfrak{X}^T} = \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathfrak{X}^T} - 2 r_T(\mathbf{y}^T, \mathbf{y}'^T). $$

In particular, this result allows to prove a kind of duality between the original variable and the outgoing traces, as one has

$$ \forall \mathbf{y} \in \mathcal{Y}, \quad \|\mathbf{y}\|_{\mathcal{Y}}^2 = \|\mathbf{x}\|_{\mathcal{X}}^2 \quad \text{with} \quad \mathbf{x} = \gamma_{\text{out}}(\mathbf{y}). $$

We then obtain a contraction property for $k_{UWVF}$ with respect to the $\mathfrak{X}$ scalar product.

**Corollary 3.27 (UWVF contraction).** We have the contraction property:

$$ \forall \mathbf{x}, \mathbf{x}' \in \mathcal{X}, \quad |k_{UWVF}(\mathbf{x}, \mathbf{x}')| = |\langle \Pi_{\mathfrak{U}} \mathbf{x}, \mathfrak{U} \mathbf{x}' \rangle_{\mathfrak{X}}| \le \|\mathbf{x}\|_{\mathfrak{X}} \|\mathbf{x}'\|_{\mathfrak{X}}. $$

*Proof.* Thanks to a Cauchy-Schwarz inequality and Proposition 3.25, one gets:

$$ \left| \sum_{T \in \mathcal{T}_{\mathfrak{h}}} \langle \Pi_{\mathfrak{U}}^T \mathbf{x}^T, \mathfrak{U}^T \mathbf{x}'^T \rangle_{\mathfrak{X}^T} \right| \le \|\Pi_{\mathfrak{U}} \mathbf{x}\|_{\mathfrak{X}} \|\mathfrak{U} \mathbf{x}'\|_{\mathfrak{X}} \le \|\mathbf{x}\|_{\mathfrak{X}} \|\mathbf{x}'\|_{\mathfrak{X}}. $$

$\square$

### 3.5.2 Introduction of a Jacobi algorithm

Corollary 3.27 induces the agreement of this Trefftz method with iterative solvers, as a contracting preconditioner can be highlighted. We propose to illustrate this by introducing the following Jacobi algorithm (leaning on the definition (49) of the sesquilinear form $a_{UWVF}$).

**Problem 3.28 (Jacobi UWVF).** Being given an initial guess $\mathbf{x}_{\mathfrak{h}}^0 \in \mathcal{X}_{\mathfrak{h}}$, determine the sequence $(\mathbf{x}_{\mathfrak{h}}^n)_{n \in \mathbb{N}^*} \subset \mathcal{X}_{\mathfrak{h}}$ such that

$$ \forall n \in \mathbb{N}^*, \quad \langle \mathbf{x}_{\mathfrak{h}}^n, \mathbf{x}' \rangle_{\mathfrak{X}} - k_{UWVF}(\mathbf{x}_{\mathfrak{h}}^{n-1}, \mathbf{x}') = \ell_{UWVF}(\mathbf{x}') \quad \text{for all } \mathbf{x}' \in \mathcal{X}_{\mathfrak{h}}, $$

which can be equivalently stated, considering the vector representations $([\mathbf{x}_{\mathfrak{h}}^n])_{n \in \mathbb{N}^*}$ of $(\mathbf{x}_{\mathfrak{h}}^n)_{n \in \mathbb{N}^*}$ in $\mathbb{C}^{N_{dof}}$, as

$$ \forall n \in \mathbb{N}^*, \quad \mathcal{D}[\mathbf{x}_{\mathfrak{h}}^n] = \mathcal{N}[\mathbf{x}_{\mathfrak{h}}^{n-1}] + [\mathbf{b}_{UWVF}], \quad (63) $$

where $\mathcal{D}, \mathcal{N} \in \mathbb{C}^{N_{dof} \times N_{dof}}$ represent the forms in the basis of $\mathcal{X}_{\mathfrak{h}}$:

$$ \mathcal{D} = (\langle \psi_J, \psi_I \rangle_{\mathfrak{X}})_{I,J \in [1, N_{dof}]} \quad \text{and} \quad \mathcal{N} = (k_{UWVF}(\psi_J, \psi_I))_{I,J \in [1, N_{dof}]} \quad (64) $$
The matrix $\mathcal{D}$ is naturally Hermitian positive definite, and then invertible: it seems particularly adapted as a left preconditioner for the linear system (63).

Indeed, as the Riesz representation $\mathcal{K}$ of $k_{UWVF}$ is represented by the matrix $\mathcal{D}^{-1}\mathcal{N}$, Corollary 3.27 and the conformal character of the approximation ensure that its spectral radius $\rho(\mathcal{D}^{-1}\mathcal{N})$ is bounded from above by 1.

**Corollary 3.29.** The spectral radius of the matrix $\mathcal{D}^{-1}\mathcal{N}$ is bounded from above by 1.

*Proof.* Let's consider generalised eigenvalue $\lambda \in \mathbb{C}$ and eigenvector $\mathbf{x}_h \in \mathcal{X}_h \setminus \{0\}$ (whose vector representation is denoted as $[\mathbf{x}_h] \in \mathbb{C}^{N_{dof}}$) associated to $\mathcal{D}$ and $\mathcal{N}$, such that $\mathcal{N}\mathbf{x}_h = \lambda \mathcal{D}\mathbf{x}_h$. In particular, one has

$$k_{UWVF}(\mathbf{x}_h, \mathbf{x}_h) = \langle \mathcal{N}[\mathbf{x}_h], [\mathbf{x}_h] \rangle = \lambda \langle \mathcal{D}[\mathbf{x}_h], [\mathbf{x}_h] \rangle = \lambda \|\mathbf{x}_h\|_{\mathcal{X}}^2,$$

and then, as $\mathcal{X}_h \subset \mathcal{X}$, one gets

$$|\lambda| = \frac{|k_{UWVF}(\mathbf{x}_h, \mathbf{x}_h)|}{\|\mathbf{x}_h\|_{\mathcal{X}}^2} \le \sup_{\mathbf{x}, \mathbf{x}' \in \mathcal{X} \setminus \{0\}} \frac{|k_{UWVF}(\mathbf{x}, \mathbf{x}')|}{\|\mathbf{x}\|_{\mathcal{X}} \|\mathbf{x}'\|_{\mathcal{X}}}.$$

Yet, the eigenvalues of $\mathcal{D}^{-1}\mathcal{N}$ coincide with the ones of this generalised eigenvalue problem. Thus, by referring to Corollary 3.27, the spectral radius can be bounded as

$$\rho(\mathcal{D}^{-1}\mathcal{N}) = \max_{\lambda \in \sigma(\mathcal{D}^{-1}\mathcal{N})} |\lambda| \le \sup_{\mathbf{x}, \mathbf{x}' \in \mathcal{X} \setminus \{0\}} \frac{|k_{UWVF}(\mathbf{x}, \mathbf{x}')|}{\|\mathbf{x}\|_{\mathcal{X}} \|\mathbf{x}'\|_{\mathcal{X}}} \le 1.$$

□

### 3.5.3 Strict contraction of the discrete preconditioned system

Yet, so as to ensure the convergence of the Jacobi algorithm 3.28, the spectral radius of the matrix $\mathcal{D}^{-1}\mathcal{N}$ should be strictly bounded from above by 1: this is the result of Proposition 3.30, which proves that the contraction property highlighted by Corollary 3.27 is strict in the finite-dimensional case. In particular, it is verified for the wave propagation problems we are interested in.

**Proposition 3.30.** We have the strict contraction property:

$$\exists \lambda_h \in [0, 1[, \quad \forall \mathbf{x}_h, \mathbf{x}_h' \in \mathcal{X}_h, \quad |k_{UWVF}(\mathbf{x}_h, \mathbf{x}_h')| \le \lambda_h \|\mathbf{x}_h\|_{\mathcal{X}} \|\mathbf{x}_h'\|_{\mathcal{X}}.$$

*Proof.* This proof is adapted from the one of [23, Proposition 2.3] in the case of the 3D time-harmonic Maxwell's equations, and relies on contradiction.

Leaning on Corollary 3.27, we suppose there exists an eigenvector $\mathbf{x}_h \in \mathcal{X}_h$ and $\lambda \in \mathbb{C}$ such that $\mathbf{x}_h \ne \mathbf{0}$, $|\lambda| = 1$ and

$$\forall \mathbf{x}_h' \in \mathcal{Y}_h, \quad k_{UWVF}(\mathbf{x}_h, \mathbf{x}_h') = \lambda \langle \mathbf{x}_h, \mathbf{x}_h' \rangle_{\mathcal{X}}. \qquad (65)$$

We will then show that $\mathbf{x}_h = \mathbf{0}$, which would be in contradiction with the hypothesis. First, relying on the definition (51) of $k_{UWVF}$, one gets

$$\begin{align*}
\|\Pi_{\mathcal{Y}}(\mathbf{x}_h - \lambda \mathfrak{U}\mathbf{x}_h)\|_{\mathcal{X}}^2 &= \|\Pi_{\mathcal{Y}}(\mathbf{x}_h)\|_{\mathcal{X}}^2 - 2\Re[\bar{\lambda} \langle \Pi_{\mathcal{Y}}(\mathbf{x}_h, \mathfrak{U}\mathbf{x}_h) \rangle_{\mathcal{X}}] + |\lambda|^2 \|\mathfrak{U}\mathbf{x}_h\|_{\mathcal{X}}^2 \\
&= \|\Pi_{\mathcal{Y}}(\mathbf{x}_h)\|_{\mathcal{X}}^2 - 2\Re[\bar{\lambda} k_{UWVF}(\mathbf{x}_h, \mathbf{x}_h)] + |\lambda|^2 \|\mathfrak{U}\mathbf{x}_h\|_{\mathcal{X}}^2 \\
&= \|\Pi_{\mathcal{Y}}(\mathbf{x}_h)\|_{\mathcal{X}}^2 - 2|\lambda|^2 \|\mathbf{x}_h\|_{\mathcal{X}}^2 + |\lambda|^2 \|\mathfrak{U}\mathbf{x}_h\|_{\mathcal{X}}^2,
\end{align*}$$

Inria
as $\mathbf{x}_h$ is an eigenvector associated to $\lambda$.

Moreover, leaning on (65) as well as the isometric and contracting characters of $\mathcal{U}$ and $\Pi_{\mathcal{U}}$ (see (54)), one gets

$$ \|\Pi_{\mathcal{U}}\mathbf{x}_h - \lambda \mathfrak{U}\mathbf{x}_h\|_{\mathcal{X}}^2 \le \|\mathbf{x}_h\|_{\mathcal{X}}^2 - 2|\lambda|^2 \|\mathbf{x}_h\|_{\mathcal{X}}^2 + |\lambda|^2 \|\mathbf{x}_h\|_{\mathcal{X}}^2 = 0, $$

as we recall that $|\lambda| = 1$, and then

$$ \Pi_{\mathcal{U}}\mathbf{x}_h = \lambda \mathfrak{U}\mathbf{x}_h. \qquad (66) $$

Moreover, let's introduce the space of cells of zero outgoing traces

$$ \tilde{\mathcal{F}}_h := \{T \in \mathcal{T}_h : \mathbf{x}_h^T = \mathbf{0} \text{ on } \partial T\}, $$

and then show that $\tilde{\mathcal{F}}_h = \mathcal{T}_h$.

(i) First, we prove that the boundary cells belong to $\tilde{\mathcal{F}}_h$. By relying on (54) and (66), we have

$$ \|\Pi_{\mathcal{U}}\mathbf{x}_h\|_{\mathcal{X}} = \|\mathfrak{U}\mathbf{x}_h\|_{\mathcal{X}} = \|\mathbf{x}_h\|_{\mathcal{X}}, $$

and the expression (60) of $\|\Pi_{\mathcal{U}}\mathbf{x}_h\|_{\mathcal{X}}^2$ leads to

$$ 0 = \|\mathbf{x}_h\|_{\mathcal{X}}^2 - \|\Pi_{\mathcal{U}}\mathbf{x}_h\|_{\mathcal{X}}^2 = \sum_{T \in \mathcal{T}_h} \sum_{F \in \mathcal{F}_{\text{ext}}(T)} \underbrace{\left[1 - \left(\frac{Z_{\partial\Omega} - Z_T}{Z_T + Z_{\partial\Omega}}\right)^2\right]}_{\ge 0} \|\mathbf{x}_h^T\|_{\mathcal{X}^T, F}^2. $$

Yet, as the field $Z_{\partial\Omega}$ is non-uniformly zero, there exists a boundary cell $T \in \mathcal{T}_h$ and $F \in \mathcal{F}_{\text{ext}}(T)$ such that $Z_{\partial\Omega} \ne 0$ and then $\mathbf{x}_h^T = \mathbf{0}$ on $F$. Thus, (66) and (50) imply that

$$ \mathbf{x}_h^T = \mathbf{0} \quad \text{and then} \quad \mathfrak{U}^T \mathbf{x}_h^T = \frac{1}{\lambda} \Pi_{\mathcal{U}}^T \mathbf{x}_h^T = \frac{1}{\lambda} \frac{Z_{\partial\Omega} - Z_T}{Z_T + Z_{\partial\Omega}} \mathbf{x}_h^T = \mathbf{0} \quad \text{on } F, $$

implying that $\mathcal{F}_{n_{\partial\Omega}}\mathbf{y}^T = \mathbf{0}$ on $F$, as $\lambda \ne 0$. Then, considering the unique continuation principle Hypothesis 2.2, it ensures that $\mathbf{y}^T = \mathbf{0}$ in $T$. In particular, we get $\mathbf{x}_h^T = \mathfrak{U}^T \mathbf{x}_h^T = \mathbf{0}$ on $\partial T$.

(ii) In a second time, let's suppose there exist two cells $T \in \mathcal{T}_h$ and $K \in \tilde{\mathcal{F}}_h$ such that they share a common face $F = \partial T \cap \partial K$. Yet, (66) implies that

$$ \Pi_{\mathcal{U}}^T \mathbf{x}_h^T = \lambda \mathfrak{U}^T \mathbf{x}_h^T \quad \text{and} \quad \Pi_{\mathcal{U}}^K \mathbf{x}_h^K = \lambda \mathfrak{U}^K \mathbf{x}_h^K \quad \text{on } F, $$

or, in other words

$$ \frac{Z_K - Z_T}{Z_T + Z_K} \mathbf{x}_h^T + \frac{2Z_T}{Z_T + Z_K} \mathbf{x}_h^K = \lambda \mathfrak{U}^T \mathbf{x}_h^T \quad \text{and} \quad \frac{Z_T - Z_K}{Z_T + Z_K} \mathbf{x}_h^K + \frac{2Z_K}{Z_T + Z_K} \mathbf{x}_h^T = \lambda \mathfrak{U}^K \mathbf{x}_h^K. $$

Having $\mathbf{x}_h^K = \mathfrak{U}^K \mathbf{x}_h^K = \mathbf{0}$ on $F$, we obtain $\mathbf{x}_h^T = \mathfrak{U}^T \mathbf{x}_h^T = \mathbf{0}$ on $F$. Similarly as (i), the unique continuation theorem ensures that $\mathbf{x}_h^T = \mathfrak{U}^T \mathbf{x}_h^T = \mathbf{0}$ on $\partial T$.

(iii) Finally, the completeness and connectivity of the neighbourhood graph allow to conclude that $\tilde{\mathcal{F}}_h = \mathcal{T}_h$ and $\mathbf{x}_h = \mathbf{0}$ in $\Omega$.

We conclude by contradiction that all the eigenvalues have their modules strictly smaller than 1. □
More generally, we showed that left-preconditioning the original discrete system (see (64))

$$ \mathcal{D}[\mathbf{x}_h] - \mathcal{N}[\mathbf{x}_h] = [\mathbf{b}_{UWVF}], $$

by $\mathcal{D}$ leads to a well-conditioned linear system

$$ (\mathcal{I} - \mathcal{D}^{-1}\mathcal{N})[\mathbf{x}_h] = \mathcal{D}^{-1}[\mathbf{b}_{UWVF}], $$

as we proved with Proposition 3.30 that the eigenvalues of $\mathcal{I} - \mathcal{D}^{-1}\mathcal{N}$ are strictly comprised in the circle of radius 1 and centre (1,0) (see Figure 2): the associated Jacobi algorithm is then (theoretically) ensured to converge, while other options may be considered.

**Remark 3.31.** The preconditioning process can be realised in a symmetric way by leaning on the square root $\mathcal{D}^{1/2}$ of the Hermitian positive definite matrix $\mathcal{D}$. Indeed, the Jacobi iteration (63) can be rewritten as

$$ \begin{align*} \mathcal{D}[\mathbf{x}_h^n] &= \mathcal{N}[\mathbf{x}_h^{n-1}] + [\mathbf{b}_{UWVF}] &\iff & \mathcal{D}^{1/2}(\mathcal{D}^{1/2}[\mathbf{x}_h^n]) &= \mathcal{N}\mathcal{D}^{-1/2}(\mathcal{D}^{1/2}[\mathbf{x}_h^{n-1}]) + [\mathbf{b}_{UWVF}] \\ &\iff & [\tilde{\mathbf{x}}_h^n] &= \mathcal{D}^{-1/2}\mathcal{N}\mathcal{D}^{-1/2}[\tilde{\mathbf{x}}_h^{n-1}] + \mathcal{D}^{-1/2}[\mathbf{b}_{UWVF}], \end{align*} $$

where we introduced the change of variable $[\tilde{\mathbf{x}}_h^n] = \mathcal{D}^{1/2}[\mathbf{x}_h^n]$.

Thus, we obtain a two-sided preconditioned linear system. In particular, this preconditioning process is systematically chosen in the numerical experiments presented in this Report.

As a summary, Proposition 3.30 has highlighted the adaptability of the Trefftz to iterative solution, as a natural preconditioner $\mathcal{D}$ can be designed to ensure a non-relaxed Jacobi algorithm will converge.

Moreover, more refined preconditioners can also be developed by separating the directional contributions in successive steps. This allows an additional reduction of the number of iterations needed to convergence [25, Section 4.3], at the price of an increased computation cost: the potential gain in terms of computation time must then be specifically studied according to the available material setting.

Nevertheless, even if this strict contraction property is theoretically proved, some of the actual eigenvalues may be really close to this bound, which strongly slows the convergence of the Jacobi algorithm [23].

Thus, more advanced iterative algorithms may be considered, leaning for example on Krylov subspaces as the Generalised Minimal Residual Method [24] (GMRES), with potential restart, deflation [6] or recycling [20]; or the biconjugate gradient stabilised [28] (BiCGSTAB).

Such approaches usually offer more robust convergence behaviours, even if they are sensitive to other eigenvalue distributions. Indeed, in the case of the GMRES, it is usually considered that eigenvalues close to the origin have to be avoided not to deteriorate the convergence [27, Theorem 35.2].

## 3.6 Interpretation of the Trefftz method from a DDM point of view

Finally, we offer a last point of view on this Trefftz approach, which may allow to allow interpretations and analysis from the Domain Decomposition Method community.

Indeed, relying on the unitary character (54) of $\mathfrak{U}$, one has for any cell $T \in \mathcal{T}_h$:

$$ \forall \mathbf{x}^T, \mathbf{x}'^T \in \mathcal{T}^T, \quad \langle \mathbf{x}^T, \mathbf{x}'^T \rangle_{\mathcal{T}^T} = \langle \mathfrak{U}^T \mathbf{x}^T, \mathfrak{U}^T \mathbf{x}'^T \rangle_{\mathcal{T}^T}, $$
<figure><img src="image_7.png" /><figcaption>Figure 2: Example of distribution of the eigenvalues of $\mathcal{I} - \mathcal{D}^{-1}\mathcal{N}$ in the complex plane, strictly comprised in the red circle of radius one and centre (1, 0).</figcaption></figure>

which allows to rewrite $\tilde{a}_{UWVF}$ as:

$$ \tilde{a}_{UWVF}(\mathbf{x}, \mathbf{x}', \mathbf{g}) = \sum_{T \in \mathcal{T}_h} \left\langle \mathfrak{U}^T \mathbf{x}^T - \widehat{\mathfrak{U}^T \mathbf{x}}, \mathfrak{U}^T \mathbf{x}'^T \right\rangle_{\mathcal{H}^T} \quad (67) $$

Yet, Hypothesis 3.16 ensures that $\mathfrak{U}^T \mathbf{x}'^T$ maps the whole space $\mathcal{H}^T = [L^2(\partial T; \mathbb{C})]^r$ when $\mathbf{x}'^T \in \mathcal{H}^T$. A solution to the Trefftz Problem 3.22 then satisfies

$$ \forall T \in \mathcal{T}_h, \quad \forall F \in \mathcal{F}(T), \quad \mathfrak{U}^T \mathbf{x}_h = \widehat{\mathfrak{U}^T \mathbf{x}_h} \quad \text{on } F. \quad (68) $$

**Remark 3.32.** In the specific case of the wave propagation problems we study in this Report, considering an internal face $F = \partial T \cap \partial K$ for $T, K \in \mathcal{T}_h$, we recover the classic Després transmission conditions [8]:

* for the Maxwell system, (68) can be rewritten as

$$ \left\{ \begin{array}{l} \gamma_t(\mathbf{e}^K) - Z_K \gamma_x^K(\mathbf{h}^K) = \gamma_t(\mathbf{e}^T) + Z_K \gamma_x^T(\mathbf{h}^T) \\ \gamma_t(\mathbf{e}^T) - Z_T \gamma_x^K(\mathbf{h}^T) = \gamma_t(\mathbf{e}^K) + Z_T \gamma_x^K(\mathbf{h}^K) \end{array} \right. , $$

which allows to recover the continuity of the tangential component and trace of the electric and magnetic fields.

* for the 2D Helmholtz equation, (68) can be rewritten as

$$ \left\{ \begin{array}{l} p^K - Z_K v_N^K = p^T + Z_K v_N^T \\ p^T - Z_T v_N^T = p^K + Z_T v_N^K \end{array} \right. , $$

which allows to recover the continuity of the pressure and normal velocity fields.

* for the 1D Helmholtz equation, (68) can be rewritten as

$$ \left\{ \begin{array}{l} p^K - Z_K n_K v^K = p^T + Z_K n_T v^T \\ p^T - Z_T n_T v^T = p^K + Z_T n_K v^K \end{array} \right. , $$

which allows to recover the continuity of the pressure and velocity fields.
Taking back the iterative problem 3.28 based on the Jacobi algorithm, it can be stated as

$$ \forall T \in \mathcal{T}_h, \forall F \in \mathcal{F}(T), \mathcal{U}^T \mathbf{x}_h^{n,T} = \mathcal{U}^T \mathbf{x}_h^{n-1} \quad \text{on } F. $$

Considering an internal face $F = \partial T \cap \partial K$, this can be rewritten as

$$ \mathcal{U}^T \mathbf{x}_h^{n,T} = \frac{Z_K - Z_T}{Z_T + Z_K} \mathbf{x}_h^{n-1,T} + \frac{2Z_T}{Z_T + Z_K} \mathbf{x}_h^{n-1,K}, $$

or else

$$ \mathcal{U}^* \mathbf{y}_1^{n,T} - Z_T \mathcal{V}^* \mathbf{y}_2^{n,T} = \mathcal{U}^* \left( \frac{(Z_K - Z_T) \mathbf{y}_1^{n-1,T} + 2Z_T \mathbf{y}_1^{n-1,K}}{Z_K + Z_T} \right) \\ + Z_T \mathcal{V}^* \left( \frac{(Z_K - Z_T) \mathbf{y}_2^{n-1,T} - 2Z_K \mathbf{y}_2^{n-1,K}}{Z_K + Z_T} \right), $$

which can be seen as a relaxed version of the Després algorithm usually considered in Domain Decomposition Methods.

## 3.7 Trefftz error estimate

Finishing with the theoretical introduction of this Trefftz method, we Report some classical tools and ideas used for proving a priori error estimates. This Section is separated in two parts:

* The first one presents a kind of Céa's Lemma, allowing to bound the Trefftz error in DG norm with respect to the error of the best-approximation $\mathbf{x}_{ba}$ in the stronger norm $\|\mathbf{x} - \mathbf{x}_{ba}\|_{X}$. Such a bound is inherent to the Trefftz formulation and can be derived for all BATFFS. Moreover, it must be followed by further estimates which are specific to the chosen discretisation, as the best approximation has to be characterised.

* Nevertheless, one may usually want to consider the error in a volumic norm, in opposition to the DG-norm which only consider jumps living on the mesh skeleton. Yet, such a relation does not seem to be derivable in a generic way for all BATFFS (as different arguments have been used in the literature for Maxwell and Helmholtz cases), due to the difference of regularity that such systems may ensure [12]. Then, we propose a presentation of the classic arguments and results that have been obtained in the case of the wave propagation problems we are interested in, with a focus on the 2D Helmholtz equation.

### 3.7.1 A Céa's Lemma for the Trefftz formulation

We first derive a kind of Céa's Lemma for the Trefftz formulation, which is derivable for any case of BATFFS.

**Proposition 3.33.** Under existence condition, let $\mathbf{y} \in \mathbb{Y}$ and $\mathbf{y}_h \in \mathbb{Y}_h$ denote the solutions of the continuous and discrete Trefftz Problems 3.9 and 3.10, in addition to an arbitrary $\mathbf{y}'_h \in \mathbb{Y}_h$. Then one has:

$$ \|\mathbf{y} - \mathbf{y}_h\|_{DG} \lesssim \|\gamma_{out}(\mathbf{y}) - \gamma_{out}(\mathbf{y}'_h)\|_{X}, \quad (69) $$

where we recall the DG norm (see Proposition 3.11) is defined for $\mathbf{w} \in \mathbb{Y}$ as

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

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

and

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

*Proof.* This proof is adapted from [4, 12, 25] and relies on the equivalent UWVF formulation of the Trefftz problem.

By first relying on the formulation equivalence Proposition 3.24 and the DG-norm definition (39), one has

$$\|\mathbf{y} - \mathbf{y}_h\|_{DG} = \Re[a(\mathbf{y} - \mathbf{y}_h, \mathbf{y} - \mathbf{y}_h)] \lesssim \Re[a_{UWVF}(\mathbf{x} - \mathbf{x}_h, \mathbf{x} - \mathbf{x}_h)] := \|\mathbf{x} - \mathbf{x}_h\|_{DG}, \quad (70)$$

where we naturally introduced the DG-norm on $\mathcal{X}$ (in agreement with the one on $\mathcal{Y}$, see (39)), $\mathbf{x} = \gamma_{\text{out}}(\mathbf{y}) \in \mathcal{X}$ and $\mathbf{x}_h = \gamma_{\text{out}}(\mathbf{y}_h) \in \mathcal{X}_h$ stand for the outgoing traces associated to $\mathbf{y}$ and $\mathbf{y}_h$ respectively.

Then, we obtain a Galerkin orthogonality property. As $\mathcal{X}_h \subset \mathcal{X}$, we have that for any $\mathbf{x}'_h \in \mathcal{X}_h$:

$$a_{UWVF}(\mathbf{x} - \mathbf{x}_h, \mathbf{x}'_h) = a_{UWVF}(\mathbf{x}, \mathbf{x}'_h) - a_{UWVF}(\mathbf{x}_h, \mathbf{x}'_h) = \ell_{UWVF}(\mathbf{x}'_h) - \ell_{UWVF}(\mathbf{x}_h) = 0.$$

Then, we obtain

$$\begin{aligned} \|\mathbf{x} - \mathbf{x}_h\|_{DG}^2 &= \Re[a_{UWVF}(\mathbf{x} - \mathbf{x}_h, \mathbf{x} - \mathbf{x}_h)] \\ &= \Re[a_{UWVF}(\mathbf{x} - \mathbf{x}_h, \mathbf{x} - \mathbf{x}'_h)], \end{aligned}$$

which leads, thanks to the Riesz representation Theorem (51), to

$$\|\mathbf{x} - \mathbf{x}_h\|_{DG}^2 = \Re[\langle (\mathcal{I} - \mathcal{K})(\mathbf{x} - \mathbf{x}_h), \mathbf{x} - \mathbf{x}'_h \rangle_{\mathcal{X}}],$$

and the Cauchy-Schwarz inequality then implies that

$$\|\mathbf{x} - \mathbf{x}_h\|_{DG}^2 \le \|(\mathcal{I} - \mathcal{K})(\mathbf{x} - \mathbf{x}_h)\|_{\mathcal{X}} \|\mathbf{x} - \mathbf{x}'_h\|_{\mathcal{X}}. \quad (71)$$

Yet, for any $\tilde{\mathbf{x}} \in \mathcal{X}$, one has

$$\|(\mathcal{I} - \mathcal{K})\tilde{\mathbf{x}}\|_{\mathcal{X}}^2 = \langle \tilde{\mathbf{x}}, \tilde{\mathbf{x}} \rangle_{\mathcal{X}} + \langle \mathcal{K}\tilde{\mathbf{x}}, \mathcal{K}\tilde{\mathbf{x}} \rangle_{\mathcal{X}} - 2 \Re[\langle \mathcal{K}\tilde{\mathbf{x}}, \tilde{\mathbf{x}} \rangle_{\mathcal{X}}].$$

By using the contraction Corollary 3.27 of $\mathcal{K}$, we get

$$\begin{aligned} \|(\mathcal{I} - \mathcal{K})\tilde{\mathbf{x}}\|_{\mathcal{X}}^2 &\le 2 \langle \tilde{\mathbf{x}}, \tilde{\mathbf{x}} \rangle_{\mathcal{X}} - 2 \Re[\langle \mathcal{K}\tilde{\mathbf{x}}, \tilde{\mathbf{x}} \rangle_{\mathcal{X}}] \\ &= 2 \Re[\langle (\mathcal{I} - \mathcal{K})\tilde{\mathbf{x}}, \tilde{\mathbf{x}} \rangle_{\mathcal{X}}] \\ &= 2 \Re[a_{UWVF}(\tilde{\mathbf{x}}, \tilde{\mathbf{x}})] \\ &= 2 \|\tilde{\mathbf{x}}\|_{DG}^2. \end{aligned}$$

Reporting this bounding in (71) for $\tilde{\mathbf{x}} = \mathbf{x} - \mathbf{x}_h \in \mathcal{X}$, we finally obtain

$$\|\mathbf{x} - \mathbf{x}_h\|_{DG} \lesssim \|\mathbf{x} - \mathbf{x}'_h\|_{\mathcal{X}}. \quad (72)$$

□
An *a priori* estimate could be derived from (69) in terms of the best-approximation in the context of BATFFS. To complete it, the quantity

$$ \inf_{\mathbf{x}'_h \in \mathcal{X}_h} \|\mathbf{x} - \mathbf{x}'_h\|_{\mathcal{X}}, $$

must be estimated in terms of a volumic (and mesh independent) norm, by making an explicit choice of the discretisation (some classic possibilities are discussed in Section 3.8).

In the classic case of Plane Waves, one may refer to [5, 17, 13, 14] for the Helmholtz equation and [12] for the Maxwell system, for complete *a priori* error estimates.

### 3.7.2 Classic arguments of volumic bounding for wave propagation problems

The previous part allowed to bound the Trefftz error in terms of the DG-norm, which seems unsatisfactory as it may not be natural to handle: the idea is then to try to consider this error with a volumic norm that should be bounded from above by the DG one.

Considering the classic results which have previously been derived for the Helmholtz [17, 13] and Maxwell systems [12], one may observe that a common idea to rely on an adjoint problem (first introduced by Monk and Wang [17]) cements both these approaches. Nevertheless, the noticeable difference of solution regularity between these problems does not allow to keep exactly the same arguments, and while an estimate in $L^2(\Omega)$-norm can be obtained for the Helmholtz equation, it is currently restricted to the weaker $\| \cdot \|_{H(div,\Omega)'}$-norm for the Maxwell system.

Thus, as these two specific cases already highlight differences in their arguments and final results, it does not seem obvious that a generic treatment may be considered in the more general framework of BATFFS. Then, we simply offer a description of the main arguments and results obtained for these two wave propagation problems.

In the case of the homogeneous Helmholtz equation (meaning that both $\varepsilon$ and $\mu$ are constant in $\Omega$), we rely on the arguments and proofs derived in [5, 17, 3, 13, 12]. So as to bound the volumic $L^2(\Omega)$-error in terms of the cells jumps (corresponding to the DG norm), we lean on the ideas of [17, 13] to introduce a global dual problem associated to a given function $\phi \in L^2(\Omega; \mathbb{C})$:

**Problem 3.34 (Dual Helmholtz problem).** Being given $\phi \in L^2(\Omega)$, find $\tilde{p} \in H^1(\Omega)$ such that one has in the sense of generalised functions

$$ \nabla \cdot (\mu^{-1} \nabla \tilde{p}) + \kappa^2 \varepsilon \tilde{p} = \phi \quad \text{in } \Omega, \qquad (73) $$

that we equip it with Fourier-Robin boundary conditions

$$ \tilde{p} - Z_{\partial\Omega} \tilde{v}_N = 0 \quad \text{with} \quad \tilde{v}_N = (i\kappa\mu)^{-1} \partial_{n_{\partial\Omega}} \tilde{p} \quad \text{on } \partial\Omega, $$

where $n_{\partial\Omega}$ stands for the outward normal from $\Omega$.

It is then possible to estimate the quantity $\langle \mathbf{y} - \mathbf{y}_h, \phi \rangle_{0,\Omega}$ in terms of $\|\phi\|_{0,\Omega}$ and $\|\mathbf{y} - \mathbf{y}_h\|_{DG}$ (thanks to the continuity of the dual solution, trace inequalities and stability estimates), which leads to:

**Corollary 3.35.** Let $\mathbf{y} = (p, \mathbf{v}) \in \mathbb{Y}$ and $\mathbf{y}_h = (p_h, \mathbf{v}_h) \in \mathbb{Y}_h$ denote the solutions of the continuous and discrete Trefftz Problems 3.9 and 3.10, in addition to the outgoing trace $\mathbf{x} = \gamma_{\text{out}}(\mathbf{y})$ and an arbitrary $\mathbf{x}'_h \in \mathcal{X}_h$. One has:

$$ \|p - p_h\|_{0,\Omega} \lesssim \text{diam}(\Omega) [h^{-1/2} + \kappa h^{1/2}] \|\mathbf{y} - \mathbf{y}_h\|_{DG} \lesssim \text{diam}(\Omega) [h^{-1/2} + \kappa h^{1/2}] \|\mathbf{x} - \mathbf{x}'_h\|_{\mathcal{X}}. $$
Nevertheless, a similar estimate as Corollary 3.35 cannot be obtained for the Maxwell system in such a strong norm as the L²(Ω) one. Following [12], this wave propagation problem does not offer comparable stability estimates for the dual problem with a non-globally divergence-free function $\phi$ on the right-hand side (choosing $\phi = y - y_h$ is only piecewise divergence-free). This usually calls for the introduction of a L²-orthogonal Helmholtz decomposition of $\phi$ as

$$ \phi = \phi_0 + \nabla \phi_1, $$

where $\phi_0 \in H(div^0; \Omega)$ and $\phi_1 \in H_0^1(\Omega)$ [18, Theorem 3.45]. Then, estimates for each term can be realised: following similar steps as in the case of the Helmholtz equation allows to derive a L² estimate for $\phi_0$, while this cannot be done for $\nabla \phi_1$.

Thus, the weaker control norm $\| \cdot \|_{H(div, \Omega)}$ has to be introduced, leading to an estimate in terms of mesh-skeleton norm [12, Theorem 4.9].

Facing this problematic, we tried to consider another angle of attack by introducing dual problems in each cell, allowing to recover a control of the L²(Ω)-norm of the error. Yet, this choice comes at the cost of the loss of a global continuity of the dual solution, as it is independently defined in each cell. Thus, the DG-norm of the error, naturally appearing in Céa's Lemma, is replaced by the stronger L²($\partial \mathcal{F}_h$) and may demand to faithfully characterise the h and $\kappa$-dependence between these two skeleton norms.

## 3.8 Classic plane waves discretisation and numerical limits

The main specificity of the Trefftz approach is to rely on basis functions which are actual solutions of the considered problem in each cell of the mesh. Then, its discretisation demands to select a finite number of such functions, which are classically chosen as analytical expressions that can be precisely estimated in the integral computations.

In the case of the wave propagation problems, a range of classic functions may be considered: Generalised Harmonic Polynomials [17, 15]; Bessel functions, which allow to take singularities into account, at the cost of specific integration techniques; or, classically, the propagative plane waves [4, 12, 25] (in opposition to the evanescent ones [26, 15, 21]) whose discretisation and numerical limits are detailed in this Section.

First, coming back to the generic point of view of time-harmonic Friedrichs systems, the plane waves can be seen as solutions associated to a specific propagation direction $\mathbf{d} \in \mathbb{R}^d$ which are 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}), $$

and satisfies

$$ \mathrm{div}(\mathcal{F}\mathbf{y}_{PW})(\mathbf{x}) = -i\kappa\lambda^{-1} \mathcal{F}_{\mathbf{d}}\mathbf{Y}_{PW} \exp(-i\kappa\lambda^{-1} \mathbf{d} \cdot \mathbf{x}). $$

Thus, ensuring that such a function is solution of the time-harmonic Friedrichs system (31) implies the following dispersion relation

$$ \mathcal{F}_{\mathbf{d}}\mathbf{Y}_{PW} = \lambda \mathcal{M}\mathbf{Y}_{PW}, $$

that we recognise as a generalised eigenvalue problem for the couple $(\mathbf{Y}_{PW}, \lambda)$.

Moreover, recalling the temporal convention $e^{i\omega t}$ in which the solutions are searched in the time-harmonic framework, their time-dependent counterpart are

$$ \mathbf{y}_{PW}(\mathbf{x}, t) = \mathbf{Y}_{PW} \exp(i\kappa(ct - \lambda^{-1} \mathbf{d} \cdot \mathbf{x})). $$

In particular, positive values of $\lambda$ corresponds to waves propagating in the direction $\mathbf{d}$, while plane waves associated to negative values go in the opposite direction.
**Remark 3.36.** In the specific case of the wave propagation problems we study in this manuscript, these plane waves are expressed as:

* for the Maxwell system, their expression can be synthesised as:

$$ \mathbf{y}(\mathbf{x}) = \begin{pmatrix} \mathbf{E} \\ \mathbf{H} \end{pmatrix} \exp (i\kappa\sqrt{\varepsilon\mu} \mathbf{d} \cdot \mathbf{x}), $$

where the amplitudes satisfy the relations $\mathbf{E} = Z \mathbf{d} \times \mathbf{H}$ and $\mathbf{H} = -Y \mathbf{d} \times \mathbf{E}$, and we recall that $Z = \sqrt{\mu/\varepsilon}$ and $Y = \sqrt{\varepsilon/\mu}$.

* for the 2D Helmholtz equation, their expression can be synthesised as:

$$ \mathbf{y}_{\text{PW}}(\mathbf{x}) = \begin{pmatrix} P \\ V \end{pmatrix} \exp (i\kappa\sqrt{\varepsilon\mu} \mathbf{d} \cdot \mathbf{x}), $$

where the amplitudes satisfy the relations $\mathbf{V} = Y P \mathbf{d}$.

Thus, plane waves allow to design actual and smooth solutions of the time-harmonic Friedrichs systems, which can easily be used as basis functions in a Treftz approach without refined quadra- ture tools. The last discretisation step then consists in describing a selection process of a finite number of propagation directions.

In the bidimensional configuration, the choice of propagation directions is naturally made by considering equispaced angles such that, being given a cell $T \in \mathcal{T}_h$ and a number of plane waves $N_{dof}^T$, the plane waves are chosen as

$$ \forall j \in [1, N_{dof}^T], \quad \phi_j^{T,PW}(\mathbf{x}) = \begin{pmatrix} 1 \\ Y_T \mathbf{d}_j^{T,pw} \end{pmatrix} \exp (i\kappa\sqrt{\varepsilon_T\mu_T} \mathbf{d}_j^{T,pw} \cdot \mathbf{x}), $$

with the propagation directions defined as (see Figure 3)

$$ \mathbf{d}_j^{T,pw} = \left( \cos \left( 2\pi \frac{j-1}{N_{dof}^T} \right), \sin \left( 2\pi \frac{j-1}{N_{dof}^T} \right) \right). $$

Yet, one will note these bases are not hierarchical, as considering two numbers $n_1 > n_2 > 0$ of dofs does not ensure the associated bases are accordingly embedded.

<figure><img src="image_8.png" /><figcaption>Figure 3: Illustration of the discretised propagation directions in 2D for $N_T = 5$.</figcaption></figure>
In the three-dimensional configuration, each propagation direction is associated to an am-
plitude space of dimension 2 (as **E** must be orthogonal to **d**). Concerning the choice of these
directions, it is not as obvious as in 2D, and different strategies may be considered:

* a discretisation strategy usually considered by neutron physicists consists in parameterising the direction by two angles belonging to $[0, \pi] \times [0, 2\pi]$ and splitting the phase space in polygons of equal area [4].

* another possibility is to consider the nodes of a surfacic mesh of the unit sphere [16] or the unit cube boundary [25] (see Figure 4). In particular, in case plane waves are chosen as basis functions in this manuscript, the latter is used to define their discretisation.

However, these methods (in particular the one we consider in the rest of this manuscript) usually
do not allow to define basis functions for any arbitrary number of dofs.

<figure><img src="image_9.png" /><figcaption>Figure 4: Illustration of the 26 discretised propagation directions in 3D for $N_T = 52$: the nodes of the cube mesh represent the 26 directions, while arrows are only drawn for the front face, for readability.</figcaption></figure>

References

[1] N. Antonić, K. Burazin, I. Crnjac, and M. Erceg. "Complex Friedrichs systems and applications". In: *Journal of Mathematical Physics* 58.10 (Oct. 2017), p. 101508. DOI: 10.1063/1.5005608.

[2] H. Barucq, T. Chaumont-Frelet, and C. Gout. "Stability analysis of heterogeneous Helmholtz problems and finite element solution based on propagation media approximation". In: *Mathematics of Computation* 86.307 (Sept. 2017), pp. 2129–2157. DOI: 10.1090/mcom/3165.

[3] A. Buffa and P. Monk. "Error estimates for the Ultra Weak Variational Formulation of the Helmholtz equation". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 42.6 (Nov. 2008), pp. 925–940. DOI: 10.1051/m2an:2008033.

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

[5] O. Cessenat and B. Despres. "Application of an Ultra Weak Variational Formulation of El-
liptic PDEs to the Two-Dimensional Helmholtz Problem". In: *SIAM Journal on Numerical
Analysis* 35.1 (Jan. 1998), pp. 255–299. DOI: 10.1137/S0036142995285873.
[6] A. Chapman and Y. Saad. "Deflated and Augmented Krylov Subspace Techniques". In: *Numerical Linear Algebra with Applications* 4.1 (Jan. 1997), pp. 43–66. DOI: 10.1002/(SICI)1099-1506(199701/02)4:1<43::AID-NLA99>3.0.CO;2-Z.

[7] D. Colton and R. Kress. *Inverse Acoustic and Electromagnetic Scattering Theory*. Vol. 93. Applied Mathematical Sciences. Cham, 2019. DOI: 10.1007/978-3-030-30351-8.

[8] B. Després. "Méthodes de décomposition de domaine pour la propagation d'ondes en régime harmonique. Le théorème de Borg pour l'équation de Hill vectorielle". Thèse de doctorat. Paris 9, Jan. 1991. URL: https://theses.fr/1991PA090032.

[9] I. Djiba. "Méthode de décomposition de domaine par méthode de Trefftz pour la propagation des ondes en géophysique." Thèse de doctorat. Université de Pau et des Pays de l'Adour, 2025. URL: https://theses.fr/s342267.

[10] A. Ern and J.-L. Guermond. "Discontinuous Galerkin Methods for Friedrichs' Systems. Part III. Multifield Theories with Partial Coercivity". In: *SIAM Journal on Numerical Analysis* 46.2 (Jan. 2008), pp. 776–804. DOI: 10.1137/060664045.

[11] A. Ern and J.-L. Guermond. "Discontinuous Galerkin Methods for Friedrichs' Systems. Part II. Second-order Elliptic PDEs". In: *SIAM Journal on Numerical Analysis* 44.6 (Jan. 2006), pp. 2363–2388. DOI: 10.1137/05063831X.

[12] R. Hiptmair, A. Moiola, and I. Perugia. "Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations". In: *Mathematics of Computation* 82.281 (July 2012), pp. 247–268. DOI: 10.1090/S0025-5718-2012-02627-5.

[13] R. Hiptmair, A. Moiola, and I. Perugia. "Plane Wave Discontinuous Galerkin Methods for the 2D Helmholtz Equation: Analysis of the p-Version". In: *SIAM Journal on Numerical Analysis* (Feb. 2011). DOI: 10.1137/090761057.

[14] R. Hiptmair, A. Moiola, and I. Perugia. "Plane Wave Discontinuous Galerkin Methods: Exponential Convergence of the h p-Version". In: *Foundations of Computational Mathematics* 16.3 (June 2016), pp. 637–675. DOI: 10.1007/s10208-015-9260-1.

[15] T. Luostari, T. Huttunen, and P. Monk. "Improvements for the ultra weak variational formulation". In: *International Journal for Numerical Methods in Engineering* 94.6 (2013), pp. 598–624. DOI: 10.1002/nme.4469.

[16] A. Moiola, R. Hiptmair, and I. Perugia. "Plane wave approximation of homogeneous Helmholtz solutions". In: *Zeitschrift für angewandte Mathematik und Physik* 62.5 (Oct. 2011), pp. 809–837. DOI: 10.1007/s00033-011-0147-y.

[17] P. Monk and D.-Q. Wang. "A least-squares method for the Helmholtz equation". In: *Computer Methods in Applied Mechanics and Engineering* 175.1 (June 1999), pp. 121–136. DOI: 10.1016/S0045-7825(98)00326-0.

[18] P. Monk. *Finite Element Methods for Maxwell's Equations*. Oxford University Press, Apr. 2003. DOI: 10.1093/acprof:oso/9780198508885.001.0001.

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

[20] M. L. Parks, E. de Sturler, G. Mackey, D. D. Johnson, and S. Maiti. "Recycling Krylov Subspaces for Sequences of Linear Systems". In: *SIAM Journal on Scientific Computing* 28.5 (Jan. 2006), pp. 1651–1674. DOI: 10.1137/040607277.

[21] E. Parolin, D. Huybrechs, and A. Moiola. "Stable approximation of Helmholtz solutions in the disk by evanescent plane waves". In: *ESAIM: Mathematical Modelling and Numerical Analysis* 57.6 (Nov. 2023), pp. 3499–3536. DOI: 10.1051/m2an/2023081.
[22] S. Pernet, M. Rivet, and 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.

[23] S. Pernet, M. Sirdey, and S. Tordeux. "Ultra-weak variational formulation for heterogeneous maxwell problem in the context of high performance computing". In: *ESAIM: Proceedings and Surveys* (2023). DOI: 10.1051/proc/202375096.

[24] Y. Saad and M. H. Schultz. "GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems". In: *SIAM Journal on Scientific and Statistical Computing* 7.3 (July 1986), pp. 856–869. DOI: 10.1137/0907058.

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

[26] R. Tezaur, L. Zhang, and C. Farhat. "A discontinuous enrichment method for capturing evanescent waves in multiscale fluid and fluid/solid problems". In: *Computer Methods in Applied Mechanics and Engineering*. Computational Methods in Fluid-Structure Interaction 197.19 (Mar. 2008), pp. 1680–1698. DOI: 10.1016/j.cma.2007.08.023.

[27] L. N. Trefethen and I. Bau David. *Numerical Linear Algebra*. Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics, Jan. 1997. DOI: 10.1137/1.9780898719574.

[28] H. A. van der Vorst. "Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems". In: *SIAM Journal on Scientific and Statistical Computing* 13.2 (Mar. 1992), pp. 631–644. DOI: 10.1137/0913035.
<figure><img src="image_10.png" /></figure>

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399