# Quasi-Trefftz Method for Aeroacoustics: Part I - Model Problem and Application to the Helmholtz Equation

Lise-Marie Imbert-Gérard, Andréa Lagardère, Guillaume Sylvand, Sébastien Tordeux

► To cite this version:

Lise-Marie Imbert-Gérard, Andréa Lagardère, Guillaume Sylvand, Sébastien Tordeux. Quasi-Trefftz Method for Aeroacoustics: Part I - Model Problem and Application to the Helmholtz Equation. RR-9587, Inria Bordeaux - Sud Ouest. 2025. hal-05053699v2

HAL Id: hal-05053699

https://inria.hal.science/hal-05053699v2

Submitted on 21 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
informatics mathematics
Inria

# Quasi-Trefftz Method for Aeroacoustics: Part I – Model Problem and Application to the Helmholtz Equation

L.-M. Imbert-Gérard, A. Lagardère, G. Sylvand, S. Tordeux

## RESEARCH REPORT

Nº 9587

April 2025

Project-Team Makutu

ISSN 0249-6399 ISRN INRIA/RR-9587-FR+ENG
<figure><img src="image_2.png" /></figure>
informatics mathematics
Inria

Quasi-Trefftz Method for
Aeroacoustics: Part I – Model
Problem and Application to the
Helmholtz Equation

L.-M. Imbert-Gérard*, A. Lagardère†‡, G. Sylvand†,
S. Tordeux‡

Project-Team Makutu

Research Report n° 9587 — version 2 — initial version April
2025 — revised version July 2025 — 45 pages

**Abstract:** This work presents a quasi-Trefftz Discontinuous Galerkin method designed for the numerical resolution of problems governed by the Helmholtz equation. The method relies on local approximate solutions of the equation as test and trial functions. The formulation is based on a first-order system expressed within the Friedrichs framework. Three families of basis functions are considered: plane waves modulated in phase or amplitude, and polynomials. These are extended to the vector-valued setting—a key novelty—and used in a variational formulation compatible with the Friedrichs structure. Numerical experiments confirm the expected convergence rate and show that the method achieves comparable accuracy with fewer degrees of freedom than classical DG approaches. This work lays the foundation for future developments, including the extension to the convected Helmholtz equation for efficient aeroacoustic simulations.

**Key-words:** Quasi-Trefftz method – Discontinuous Galerkin – Helmholtz equation – Friedrichs’ systems – Aeroacoustic – Wave propagation

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

<footnote>* University of Arizona, Tucson, USA.</footnote>

<footnote>† Airbus Central R&T, Issy-les-Moulineaux, France.</footnote>

<footnote>‡ Équipe Makutu, Inria, UPPA, UMR CNRS 5142, TotalEnergies, France.</footnote>
# Méthode quasi-Trefftz pour l'aéroacoustique :
## Partie I – Problème modèle et application à l'équation d'Helmholtz

**Résumé** : Ce travail présente une méthode de Galerkin discontinue de type quasi-Trefftz conçue pour la résolution numérique de problèmes modélisés par l'équation de Helmholtz. Cette méthode utilise comme fonctions test des solutions approchées locales de l'équation. La formulation repose sur un système de premier ordre exprimé dans le cadre des systèmes de Friedrichs. Trois familles de fonctions de base sont considérées : des ondes planes généralisées, dont la phase ou l'amplitude est perturbée, ainsi que des polynômes. Ces fonctions sont étendues au cas vectoriel-une nouveauté majeure dans ce travail-puis utilisées dans une formulation variationnelle adaptée à la structure des systèmes de Friedrichs. Les simulations numériques confirment l'ordre de convergence attendu et montrent que cette approche atteint une précision comparable à celle des méthodes classiques de Galerkin discontinue, tout en nécessitant moins de degrés de liberté. Ce travail constitue une base pour de futurs développements, notamment l'adaptation à l'équation de Helmholtz convectée, et contribue à la conception d'outils numériques efficaces pour la simulation aéroacoustique.

**Mots-clés** : Méthode de quasi-Trefftz – Galerkin Discontinu – Équation de Helmholtz – Système de Friedrichs – Aéroacoustique – Propagation d'ondes
# Contents

<table><tr><td>1 Introduction</td><td>4</td></tr><tr><td>2 Model problem</td><td>5</td></tr><tr><td>2.1 Notation and physical parameters</td><td>6</td></tr><tr><td>2.2 Derivation from Fundamental Principles</td><td>7</td></tr><tr><td>2.2.1 Mass conservation equation</td><td>7</td></tr><tr><td>2.2.2 Momentum equation</td><td>8</td></tr><tr><td>2.2.3 The Helmholtz equation</td><td>9</td></tr><tr><td>2.2.4 The convected Helmholtz equation</td><td>10</td></tr><tr><td>2.3 Boundary conditions</td><td>11</td></tr><tr><td>2.3.1 Prandlt-Glauert-Lorentz transformation</td><td>11</td></tr><tr><td>2.3.2 Transformation of boundary conditions</td><td>12</td></tr><tr><td>2.4 A well-posed formulation</td><td>13</td></tr><tr><td>3 Mesh</td><td>13</td></tr><tr><td>4 Quasi-Trefftz method for Friedrichs systems</td><td>14</td></tr><tr><td>4.1 Friedrichs' formulation</td><td>14</td></tr><tr><td>4.2 Numerical hyperbolic flux</td><td>16</td></tr><tr><td>4.3 Discontinuous Galerkin variational formulation</td><td>17</td></tr><tr><td>4.4 Trefftz Discontinuous Galerkin variational formulation</td><td>19</td></tr><tr><td>4.5 Quasi-Trefftz Discontinuous Galerkin variational formulation</td><td>20</td></tr><tr><td>5 Application to the Helmholtz equation</td><td>21</td></tr><tr><td>5.1 Friedrichs' system formulation</td><td>21</td></tr><tr><td>5.2 Quasi-Trefftz functions</td><td>23</td></tr><tr><td>5.2.1 Scalar quasi-Trefftz functions</td><td>24</td></tr><tr><td>5.2.2 Vector-valued quasi-Trefftz spaces</td><td>34</td></tr><tr><td>6 Numerical results</td><td>35</td></tr><tr><td>6.1 Definition of the test case</td><td>35</td></tr><tr><td>6.2 Verification of the quasi-Trefftz Property</td><td>36</td></tr><tr><td>6.3 Quasi-Trefftz approximate solutions</td><td>37</td></tr><tr><td>7 Conclusion</td><td>39</td></tr><tr><td>A Horner's method</td><td>43</td></tr><tr><td>A.1 Horner's method in 1D</td><td>43</td></tr><tr><td>A.2 Horner's method in 2D</td><td>43</td></tr></table>
# 1 Introduction

Noise generated by reactors propagating inside the cabin and around an aircraft during the different flight phases (take-off, cruise, landing, taxying) is a significant factor in evaluating acoustic performance of an aircraft. At Airbus, a primary objective is to accurately model this sound propagation, with simulations taking into account air flow and variations in temperature, density and pressure. Indeed, there is an important heterogeneity in the jet-flow of the engines and in the boundary layer near the surface of the aircraft. Currently, several programs perform such simulations, primarily using a Finite Element Method coupled with a Boundary Element Method (FEM-BEM) [AJRT11, CES14]. The classical FEM requires a very refined mesh to account for all fluid variations. This requires a large amount of memory to solve the problem, along with significant computational time. We propose replacing the conventional FEM with a quasi-Trefftz Discontinuous Galerkin (DG) method, which should approximate the solution as accurately, but with fewer degrees of freedom.

First, we briefly recall the concept of Trefftz methods, from which quasi-Trefftz methods originate. Trefftz methods have been extensively studied for the approximation of acoustic problems, particularly in the homogeneous setting, which reduces to the standard Helmholtz equation. A survey of these methods is presented in [HMP16]. Their particularity is to take trial and test functions belonging to discrete spaces of local solutions of the governing PDE. This approach gives different variational formulations, like the so-called Ultra Weak Variational Formulation (UWVF) introduced by O. Cessenat and B. Després [Ces96, CD98]. One of the main benefits of the UWVF is that it only involves integrals on the skeleton of the mesh, which greatly reduces memory and computational costs compared to the standard Finite Element Method (FEM). Additionally, we can construct various discrete spaces for basis functions. For example, A. Buffa and P. Monk compare two of them [BM08]; one using Plane Waves (PW) and one using Bessel functions, which are fundamental solutions of the Helmholtz equation in 2D. In the homogeneous case, i.e. with constant coefficients, numerous results have been demonstrated. Multiple convergence proofs were derived [BM08, CD98, MW99] with some convergence constants determined explicitly [Moi11]. However, the problem we are interested in is characterized by heterogeneous coefficients in which case there are in general no exact solutions available to define explicitly and construct discrete spaces. This means that we cannot apply Trefftz methods here, and instead quasi-Trefftz methods can be introduced. Instead of using exact local solutions, these methods rely on approximate local solutions, referred to as quasi-solutions [IGS24].

The long term goal is to tackle an aeroacoustic problem, modeled by the convected Helmholtz equation, in an infinite three-dimensional domain. In this report, we
first introduce the physical origins of this equation. Then, for a specific and simpler case—the heterogeneous Helmholtz equation in a two-dimensional bounded domain—we describe the DG variational formulation we work with. The formulation we choose fits within the Friedrichs’ framework with some specific numerical fluxes. The construction of quasi-Trefftz basis functions was implemented and tested for their known high-order approximation properties. The corresponding quasi-Trefftz method was also implemented to illustrate high-order convergence. Some numerical results are presented in Section 6, including comparison with a classical Discontinuous Galerkin (DG) method. This is promising for reducing memory costs, and will form the foundation of future work on the convected Helmholtz equation in three dimensions.

## 2 Model problem

The purpose of this work is to introduce the quasi-Trefftz method and to apply it to an aeroacoustic problem that represents the noise propagation in a heterogeneous domain. In this setting, the fluid flow is assumed to be given and smooth, and the convected Helmholtz equation is used as a model for acoustic wave propagation. To model this problem, we make several assumptions about the flow characteristics. These assumptions are detailed below, followed by the mathematical formulation of the aeroacoustic problem.

**Assumption 1.** The convected flow is heterogeneous, isentropic everywhere and irrotational, with isentropic perturbations.

**Problem 1 (Convected Helmholtz equation).** We consider the propagation of sound through a heterogeneous flow, as described in Assumption 1. The spatial dimension is $d \in \{2, 3\}$, and the unknown $\varphi$ denotes the velocity potential. The wave number $\kappa_0$, the Mach vector $M_0$, and the fluid density $\rho_0$ are assumed to be smooth functions. Under these assumptions, the velocity potential satisfies the convected Helmholtz equation:

$$ \rho_0(\kappa_0^2\varphi + i\kappa_0 M_0 \cdot \nabla\varphi) + \nabla \cdot (\rho_0(\nabla\varphi - (M_0 \cdot \nabla\varphi)M_0 + i\kappa_0\varphi M_0)) = 0 \quad \text{in } \mathbb{R}^d. \quad (1) $$

Throughout this report, we consider only the subsonic case, as stated in the following assumption.

**Assumption 2 (Subsonic flow).** The absolute ratio of the carrier flow velocity field to the speed of sound is less than 1, i.e., $|M_0| < 1$ everywhere.

In the next subsection, we provide a detailed explanation of all the variables involved in the convected Helmholtz equation (1).
## 2.1 Notation and physical parameters

The convected Helmholtz equation is generally defined in $\mathbb{R}^d$ with $d = 3$, but in this work, we often restrict our attention to the case of $d = 2$ for the sake of simplicity. Yet many results can be easily extended to 3D.

Rather than in the time-domain, the problem of interest is expressed in the frequency-domain, and $\omega$ denotes the time-frequency. Depending on the formulations, the various frequency-domain unknowns describing the noise signal are functions on the space variable:

* $p : \boldsymbol{x} \in \mathbb{R}^d \mapsto p(\boldsymbol{x}) \in \mathbf{C}$, the pressure perturbation relative to the atmospheric pressure,

* $v : \boldsymbol{x} \in \mathbb{R}^d \mapsto v(\boldsymbol{x}) \in \mathbf{C}^d$, the velocity perturbation relative to the carrier flow velocity field,

* $\varphi : \boldsymbol{x} \in \mathbb{R}^d \mapsto \varphi(\boldsymbol{x}) \in \mathbf{C}$, the potential velocity relative to the reference state, such as $v(\boldsymbol{x}) = \nabla\varphi(\boldsymbol{x})$.

The noise signal of interest propagates in a given air flow, modeled by a given fluid flow surrounding an aircraft. The characteristics of the flow are considered to be problem data.

*Notation* 1. The density of the fluid is $\rho_0$, $v_0$ is the velocity of the fluid, $c_0$ is the speed of sound in the fluid, $\kappa_0$ is the wave number and $M_0$ is the Mach vector. The open bounded Lipschitz domain $\Omega \subset \mathbb{R}^d$ is the domain surrounding the aircraft within which the fluid flow is not necessarily constant.

This fluid flow is assumed to be constant far enough from the aircraft, and inhomogeneous only in a neighborhood of the aircraft.

*Assumption* 3. For $\Omega^\infty := \mathbb{R}^d \setminus \Omega$, there exists constants ($\rho_\infty$, $v_\infty$, $c_\infty$) $\in \mathbb{R} \times \mathbb{R}^d \times \mathbb{R}$, denoting for convenience $M_\infty = \frac{v_\infty}{c_\infty}$, such that

* $\rho_0(\boldsymbol{x}) = \rho_\infty \in \mathbb{R}, \forall \boldsymbol{x} \in \Omega^\infty,$

* $v_0(\boldsymbol{x}) = v_\infty \in \mathbb{R}^d, \forall \boldsymbol{x} \in \Omega^\infty,$

* $c_0(\boldsymbol{x}) = c_\infty \in \mathbb{R}, \forall \boldsymbol{x} \in \Omega^\infty,$

* $\kappa_0(\boldsymbol{x}) = \kappa_\infty = \frac{\omega}{c_\infty} \in \mathbb{R}, \forall \boldsymbol{x} \in \Omega^\infty,$

* $M_0(\boldsymbol{x}) = M_\infty \in \mathbb{R}^d, \forall \boldsymbol{x} \in \Omega^\infty.$
## 2.2 Derivation from Fundamental Principles

The Helmholtz convected equation refers to many different expressions [Cam07]. The derivation of (1) in homogeneous case presented here follows [RH04]. The model considers coupled time-dependent density $ρ$, pressure $p$ and velocity $v$, and later the potential $φ$. The derivation then relies on considering a small perturbation of a stationary mean flow, respectively defined by $ρ'$, $p'$, $v'$, $φ'$ and $ρ_0, p_0, v_0, φ_0$. To linearize the various equations of interest, consider

$$ρ = ρ_0 + ρ', \quad p = p_0 + p', \quad v = v_0 + v' \quad \text{and} \quad φ = φ_0 + φ'. \quad (2)$$

Moreover, the time-harmonic model is obtained by assuming that each perturbation quantity varies in the time-domain quantities as

$$\left\{
\begin{aligned}
p'(\boldsymbol{x}, t) &= \Re(p'(\boldsymbol{x})e^{-i\omega t}), \\
v'(\boldsymbol{x}, t) &= \Re(v'(\boldsymbol{x})e^{-i\omega t}), \\
\varphi'(\boldsymbol{x}, t) &= \Re(\varphi'(\boldsymbol{x})e^{-i\omega t}), \\
\rho'(\boldsymbol{x}, t) &= \Re(\rho'(\boldsymbol{x})e^{-i\omega t}).
\end{aligned}
\right.
\qquad (3)$$

### 2.2.1 Mass conservation equation

We start by considering a very important and common equation, the time-dependent mass conservation equation:

$$\frac{\partial \rho}{\partial t} + \text{div}(\rho v) = 0, \qquad (4)$$

and the mean flow is assumed to satisfy

$$\frac{\partial \rho_0}{\partial t} + \text{div}(\rho_0 v_0) = 0.$$

Plugging in (2) into (4) and neglecting quadratic terms for the perturbation quantities then gives

$$\frac{\partial \rho'}{\partial t} + \text{div}(\rho_0 v') + \text{div}(\rho' v_0) = 0. \qquad (5)$$

The first hypothesis, stating that *the perturbations are isentropic*, can be modeled [RH04, F.19] by

$$p' = c_0^2 \rho'$$

Substituting $ρ'$ then gives a time-dependent equation coupling $p'$ and $v'$:

$$\frac{1}{c_0^2} \frac{\partial p'}{\partial t} + \text{div}(\rho_0 v') + \text{div}\left(\frac{v_0 p'}{c_0^2}\right) = 0.$$
Finally, (3) leads to the linearized time-harmonic mass conservation law

$$ \frac{-i\omega p'}{c_0^2} + \operatorname{div}(\rho_0 v') + \operatorname{div}\left(\frac{v_0 p'}{c_0^2}\right) = 0. \qquad (6) $$

### 2.2.2 Momentum equation

The momentum equation,

$$ \frac{\partial v}{\partial t} + (v \cdot \nabla)v = -\frac{\nabla p}{\rho}, $$

governs the evolution of fluid velocity in terms of pressure and density.

The first hypothesis, stating that the flow is isentropic everywhere can be modeled [RH04] by

$$ \frac{p}{\rho^{\gamma}} = C, \quad \text{with } C \in \mathbb{R}. $$

where $\gamma$ represents the adiabatic index, which is a real constant for the specific gas under consideration. For air, $\gamma$ is typically taken as 1.4.

In particular, this implies that on the one hand

$$ \forall j = \{1, \dots, d\} \quad \partial_j \left( \frac{p}{\rho} \right) = \partial_j \left( C^{\frac{1}{\gamma}} p^{1-\frac{1}{\gamma}} \right) = C^{\frac{1}{\gamma}} \left( 1 - \frac{1}{\gamma} \right) p^{-\frac{1}{\gamma}} \partial_j p = \left( \frac{\gamma - 1}{\gamma} \right) \frac{\partial_j p}{\rho}, $$

and on the other hand from (2)

$$ C^{\frac{1}{\gamma}} = \frac{p^{\frac{1}{\gamma}}}{\rho} = \frac{p_0^{\frac{1}{\gamma}}}{\rho_0} \Rightarrow \frac{p}{\rho} = \frac{p_0^{\frac{1}{\gamma}} p^{1-\frac{1}{\gamma}}}{\rho_0} = \frac{p_0^{\frac{1}{\gamma}} (p_0 + p')^{1-\frac{1}{\gamma}}}{\rho_0}, $$

and with the Taylor expansion of the function $p \mapsto p^{1-\frac{1}{\gamma}}$, we obtain

$$ \frac{p_0^{\frac{1}{\gamma}} p^{1-\frac{1}{\gamma}}}{\rho_0} = \frac{p_0^{\frac{1}{\gamma}}}{\rho_0} \left( p_0^{1-\frac{1}{\gamma}} + p' \left(1 - \frac{1}{\gamma}\right) p_0^{-\frac{1}{\gamma}} + \underset{p' \to 0}{O}(p'^2) \right) = \frac{p_0}{\rho_0} + \frac{p' \gamma - 1}{\rho_0 \gamma} + \underset{p' \to 0}{O}(p'^2). \quad (7) $$

So the momentum equation becomes

$$ \frac{\partial v}{\partial t} + (v \cdot \nabla)v + \frac{\gamma}{\gamma - 1} \nabla \left(\frac{p}{\rho}\right) = 0, \qquad (8) $$

and with (7) and (2), this implies

$$ \frac{\partial(v_0 + v')}{\partial t} + \left( (v_0 + v') \cdot \nabla \right) (v_0 + v') + \frac{\gamma}{\gamma - 1} \nabla \left( \frac{p_0}{\rho_0} + \frac{p' \gamma - 1}{\rho_0 \gamma} + \underset{p' \to 0}{O}(p'^2) \right) = 0. $$
After linearizing and subtracting the steady-state equation, we have

$$ \frac{\partial \boldsymbol{v}'}{\partial t} + (\boldsymbol{v}_0 \cdot \nabla)\boldsymbol{v}' + (\boldsymbol{v}' \cdot \nabla)\boldsymbol{v}_0 + \frac{\nabla p'}{\rho_0} - \frac{\nabla \rho_0}{\rho_0^2} p' = 0. $$

Finally, (3) leads to the linearized time harmonic momentum equation

$$ -i\omega\boldsymbol{v}' + (\boldsymbol{v}_0 \cdot \nabla)\boldsymbol{v}' + (\boldsymbol{v}' \cdot \nabla)\boldsymbol{v}_0 + \frac{\nabla p'}{\rho_0} - \frac{\nabla \rho_0}{\rho_0^2} p' = 0. \quad (9) $$

### 2.2.3 The Helmholtz equation

From here the Helmholtz equation can be derived under the hypothesis on a constant equilibrium flow.

Combining the mass conservation equation (6) and the momentum equation (9) leads to the following time-harmonic system coupling $p'$ and $\boldsymbol{v}'$

$$ \left\{ \begin{array}{l} \displaystyle \frac{-i\omega p'}{\rho_0 c_0^2} + \frac{\nabla \rho_0}{\rho_0} \cdot \boldsymbol{v}' + \operatorname{div}(\boldsymbol{v}') + \frac{1}{\rho_0} \operatorname{div}\left(\frac{\boldsymbol{v}_0}{c_0^2}\right) p' + \frac{\boldsymbol{v}_0}{\rho_0 c_0^2} \cdot \nabla p' = 0, \\[1.5ex] \displaystyle -i\omega \rho_0 \boldsymbol{v}' + \rho_0(\boldsymbol{v}_0 \cdot \nabla)\boldsymbol{v}' + \rho_0(\boldsymbol{v}' \cdot \nabla)\boldsymbol{v}_0 + \nabla p' - \frac{\nabla \rho_0}{\rho_0} p' = 0. \end{array} \right. \quad (10) $$

Under the assumption that the equilibrium characteristics $\rho_0$, $\boldsymbol{v}_0$ and $c_0$ are constant, the system (10) becomes

$$ \left\{ \begin{array}{l} \displaystyle \frac{-i\omega p'}{\rho_0 c_0^2} + \operatorname{div}(\boldsymbol{v}') + \frac{\boldsymbol{v}_0}{\rho_0 c_0^2} \cdot \nabla p' = 0, \\ \displaystyle -i\omega \rho_0 \boldsymbol{v}' + \rho_0(\boldsymbol{v}_0 \cdot \nabla)\boldsymbol{v}' + \nabla p' = 0. \end{array} \right. $$

Further assuming that $\boldsymbol{v}_0 = 0$, we obtain the time-harmonic system:

$$ \left\{ \begin{array}{l} \displaystyle \frac{-i\omega p'}{\rho_0 c_0^2} + \operatorname{div}(\boldsymbol{v}') = 0, \\ \displaystyle -i\omega \rho_0 \boldsymbol{v}' + \nabla p' = 0. \end{array} \right. \quad (11) $$

This system is equivalent to the classical heterogeneous Helmholtz equation

$$ \Delta p'(\boldsymbol{x}) + \kappa_0^2(\boldsymbol{x}) p'(\boldsymbol{x}) = 0, \quad \forall \boldsymbol{x} \in \mathbb{R}^d, \quad (12) $$

with $\kappa_0(\boldsymbol{x}) = \frac{\omega}{c_0(\boldsymbol{x})}$ for all $\boldsymbol{x} \in \mathbb{R}^d$.

**Remark 2.2.1.** We observe that when $\rho_0$, $\boldsymbol{v}_0$ and $c_0$ are constant in (10), then $p'$ solves the convected Helmholtz equation (1) under the restrictive assumption of a homogeneous equilibrium flow. However, this equation in an inhomogeneous equilibrium flow requires additional work, and is the topic of the next section.
### 2.2.4 The convected Helmholtz equation

We now aim to recover the convected Helmholtz equation (1) in inhomogeneous medium from the pressure-velocity system (10). Additional assumptions on the fluid flow lead to a scalar formulation involving only the velocity potential. Under the hypothesis that the convected flow is irrotational, the velocity field can be expressed as $v = \nabla\varphi$, with $\varphi$ the velocity potential and $v_0 = \nabla\varphi_0$. Then the linearized time-harmonic mass conservation law simply becomes

$$ \frac{-i\omega p'}{c_0^2} + \operatorname{div}(\rho_0 \nabla\varphi') + \operatorname{div}\left(\frac{v_0 p'}{c_0^2}\right) = 0. \quad (13) $$

To analyze the influence of this potential velocity on the momentum equation, we return to equation (8) component-wise before linearizing. For $j \in [1, d]$, we have

$$ \frac{\partial v_j}{\partial t} + \sum_{i=1}^{d} v_i \partial_i v_j + \frac{\gamma}{\gamma-1} \partial_j \left( \frac{p}{\rho} \right) = 0 \Rightarrow \frac{\partial(\partial_j \varphi)}{\partial t} + \sum_{i=1}^{d} v_i \partial_i \varphi + \frac{\gamma}{\gamma-1} \partial_j \left( \frac{p}{\rho} \right) = 0. $$

As a consequence we obtain

$$ \partial_j \frac{\partial \varphi}{\partial t} + \sum_{i=1}^{d} v_i \partial_j v_i + \frac{\gamma}{\gamma-1} \partial_j \left( \frac{p}{\rho} \right) = 0. $$

Finally, the fact that $\partial_j(v_i^2) = 2v_i \partial_j v_i$ leads to

$$ \partial_j \left( \frac{\partial \varphi}{\partial t} + \frac{1}{2} \sum_{i=1}^{d} v_i^2 + \frac{\gamma}{\gamma-1} \frac{p}{\rho} \right) = 0. $$

We simplify this equation, which gives

$$ \frac{\partial \varphi}{\partial t} + \frac{1}{2} |v|^2 + \frac{\gamma}{\gamma-1} \frac{p}{\rho} = C, \quad \text{with } C \in \mathbb{R}. $$

For the perturbed unknown and using the Taylor expansion computed in (7), we obtain

$$ \frac{\partial \varphi_0}{\partial t} + \frac{\partial \varphi'}{\partial t} + \frac{1}{2} |v_0|^2 + v_0 \cdot v' + \frac{1}{2} |v'|^2 + \frac{\gamma}{\gamma-1} \frac{p_0}{\rho_0} + \frac{p'}{\rho_0} + O(p'^2) = C. $$

Linearizing and subtracting the steady-state equation for the stationary quantities gives

$$ \frac{\partial \varphi'}{\partial t} + v_0 \cdot v' + \frac{p'}{\rho_0} = 0. $$
Together with (13), in the harmonic regime this gives the following $p-\varphi$ system

$$ \left\{ \begin{array}{l} \displaystyle \frac{-i\omega p'}{c_0^2} + \operatorname{div}(\rho_0 \nabla\varphi') + \operatorname{div}\left(\frac{\mathbf{v}_0 p'}{c_0^2}\right) = 0, \\ \displaystyle -i\omega\rho_0\varphi' + \rho_0\mathbf{v}_0 \cdot \nabla\varphi' + p' = 0. \end{array} \right. \qquad (14) $$

This system is equivalent to the convected Helmholtz equation (1), as can be seen by simple elimination of $p'$ to obtain a second order equation for $\varphi'$ with the wave number $\kappa_0 = \frac{\omega}{c_0}$ and the Mach vector $\mathbf{M}_0 = \frac{\mathbf{v}_0}{c_0}$.

## 2.3 Boundary conditions

As explained in subsection 2.1, generally, for aeroacoustic simulations, we assume that there exists a domain $\Omega$ in which all parameters vary, and outside this domain they are assumed to remain constant. In this paper we focus on solving the problem inside the bounded domain $\Omega$. To ensure that the problem is well-posed on $\Omega$, appropriate boundary conditions will be defined. To do so, we use the Prandtl-Glauert-Lorentz (PGL) transformation [BRT22]. Indeed, we know that outside of $\Omega$, the parameters are constant. Thus, the PGL transformation allows us to switch from the classical Helmholtz equation to the convected Helmholtz equation. This will facilitate the definition of boundary conditions for the convected equation. We now provide a brief introduction to this transformation.

### 2.3.1 Prandlt-Glauert-Lorentz transformation

Let $p$ be a function satisfying the convected Helmholtz equation (1) on $\Omega^\infty = \mathbb{R}^d \setminus \Omega$, such that all the parameters involved remain constant. The PGL transformation changes $p$ into the transformed function $\tilde{p}$ given by

$$ \tilde{p}(\tilde{\mathbf{x}}) = \exp\left(i\frac{\kappa_0}{\alpha^2}\mathbf{M}_0 \cdot \mathbf{x}\right)p(\mathbf{x}), \quad \text{with } \tilde{\mathbf{x}} = B^{-1}\mathbf{x}, \qquad (15) $$

where $\kappa_0$ and $\mathbf{M}_0$ are defined in equation (1), the Lorentz factor and the anisotropy matrix are given by

$$ \alpha := \sqrt{1 - |\mathbf{M}_0|^2} \quad \text{and} \quad B := \mathrm{Id} - \frac{1}{\alpha + 1} \mathbf{M}_0 \mathbf{M}_0^T. $$

Remark 2.3.1. Matrix $B$ can also be written as $B = \mathrm{Id} + \frac{\alpha - 1}{|\mathbf{M}_0|^2} \mathbf{M}_0 \mathbf{M}_0^T$, due to the expression of $\alpha$. Moreover, we have the following expressions for the inverse matrix: $B^{-1} = \mathrm{Id} + \frac{1}{\alpha(\alpha + 1)} \mathbf{M}_0 \mathbf{M}_0^T = \mathrm{Id} + \frac{1}{|\mathbf{M}_0|^2} \left(\frac{1}{\alpha} - 1\right) \mathbf{M}_0 \mathbf{M}_0^T$.
The goal of this transformation is that $\tilde{p}$ satisfies the classical Helmholtz equation

$$ \tilde{\Delta}\tilde{p} - \tilde{\kappa}^2\tilde{p} = 0, \quad \text{on } \tilde{\Omega}^{\infty}, $$

with $\tilde{\Omega}$ the PGL transform of $\Omega$, $\tilde{\Delta}$ the Laplacian operator in PGL coordinates, and $\tilde{\kappa} = \frac{\kappa_0}{\alpha}$.

### 2.3.2 Transformation of boundary conditions

Now, we know that the PGL transformation converts the convected Helmholtz equation into the standard Helmholtz, allowing us to define some boundary conditions for the convected Helmholtz equation more easily. Indeed, since we want to find the outgoing solution of this convected equation, we apply the classical absorbing boundary condition for the Helmholtz equation:

$$ \partial_{\tilde{n}}\tilde{p} - i\tilde{\kappa}\tilde{Z}\tilde{p} = g, \quad (16) $$

where $\tilde{n}$ is the outgoing normal vector of $\tilde{\Omega}$, $\tilde{Z}$ a local operator and $g$ a known function. Now, we rewrite this boundary condition according to the variables $p$ and $v$. We start by expressing $\tilde{\nabla}\tilde{p}$ in terms of $\nabla p$. Indeed we have

$$ \tilde{\nabla}\tilde{p}(\tilde{x}) \cdot \tilde{n} = \left( B\nabla p(\mathbf{x}) + i\frac{\kappa_0 M_0}{\alpha} p(\mathbf{x}) \right) \cdot \frac{B^{-1}\mathbf{n}}{|B^{-1}\mathbf{n}|}, $$

due to the PGL transformation (15), and by using $\tilde{n} = \frac{B^{-1}\mathbf{n}}{|B^{-1}\mathbf{n}|}$.

Finally, by using the symmetry of $B^{-1}$, the boundary condition (16) becomes

$$ \nabla p(\mathbf{x}) \cdot \mathbf{n} - i \frac{\kappa_0}{\alpha} Z(\mathbf{x}) p(\mathbf{x}) = g(\mathbf{x}), \quad (17) $$

with

$$ Z(\mathbf{x}) = |B^{-1}\mathbf{n}|\tilde{Z}(\tilde{x}) - \frac{M_0 \cdot \mathbf{n}}{\alpha}. $$

To conclude, we replace $\nabla p$ with the fact that $(p, v)$ satisfies system (10) with constant coefficients, which leads to the following boundary condition for the convected Helmholtz problem

$$ (i\omega\rho_0 v - \rho_0(v_0 \cdot \nabla)v) \cdot \mathbf{n} - i\frac{\kappa_0}{\alpha}Z(\mathbf{x})p(\mathbf{x}) = g(\mathbf{x}). \quad (18) $$

*Remark 2.3.2.* In the case where $v_0 = 0$, i.e., $M_0 = 0$, this boundary condition simplifies to:

$$ v \cdot n - \frac{Z}{\rho_0 c_0} p = g. $$
## 2.4 A well-posed formulation

In this subsection, we briefly describe the outline of the proof of well-posedness for the Helmholtz problem on a bounded domain, equipped with the boundary conditions introduced previously.

This section focuses on the following boundary-value problem governed by Helmholtz equation:

$$ egin{cases} 
\Delta p + \frac{\omega^2}{c_0^2} p = 0 & \text{in } \Omega, \\
\partial_n p - \frac{i\omega}{c_0} p = g & \text{on } \partial\Omega,
\end{cases} 
\qquad (19) $$

where $\Omega \in \mathbb{R}^d$ denotes a bounded domain with spatially varying sound speed $c_0$. Well-posedness of this problem can be proved following standard techniques: applying the Fredholm alternative [Bre05] thanks to the compactness and injectivity of certain operators, and concluding thanks to the unique continuation principle [Mü54].

## 3 Mesh

From this section onward, we restrict our attention to dimension $d = 2$. In order to develop a quasi-Trefftz Discontinuous Galerkin method, the domain is discretized into several elements. In this section, we provide all necessary details about the mesh.

Notation 2. We introduce the following notations, illustrated in Figure 1:

* the set of mesh elements is denoted by $\mathcal{T}$,

* a generic element is denoted by $T \in \mathcal{T}$,

* $\partial T$ is the boundary of the element $T$; it can be decomposed into faces $f$, the set of faces of the element $T$ is denoted by $\mathcal{F}_T$,

* $\mathcal{F}_{int}$ is the set of all interior faces, i.e., faces shared by two adjacent elements,

* the element $V \in \mathcal{T}$ is the neighbor of $T$ through the face $f \in \mathcal{F}_T$,

* $\mathcal{F}_{ext}$ is the set of exterior faces, i.e., all faces for which $f \subset \partial\Omega$,

* for any function $g$, $g_T$ denotes its restriction to the element $T$.

We assume that the domain $\Omega$ is a bounded polygonal subset of $\mathbb{R}^2$. Each element $T \in \mathcal{T}$ is a polygon, and the mesh is conforming, meaning that each interior face belongs to exactly two elements, while each boundary face belongs to a single element and lies on $\partial\Omega$.
<figure><img src="image_3.png" /></figure>

$f \in \mathcal{F}_{\text{ext}}$

$f \in \mathcal{F}_{\text{int}}$

<caption>Figure 1: Domain mesh visualization</caption>

# 4 Quasi-Trefftz method for Friedrichs systems

This section presents a general framework for a quasi-Trefftz method to tackle boundary-value problems governed by Friedrichs systems. The next section shows the application of this framework to the Helmholtz equation.

## 4.1 Friedrichs' formulation

Friedrichs' formulation for hyperbolic systems [Fri58] describes a family of first-order linear differential systems with properties of symmetry and positivity.

**Definition 4.1.1 (Hyperbolic system).** A hyperbolic system is a system of partial differential equations that can be written in the form

$$M(\boldsymbol{x}) \partial_t U(\boldsymbol{x}, t) + \sum_{j=1}^{d} F_j \partial_j U(\boldsymbol{x}, t) = 0, \quad \forall \boldsymbol{x} \in \mathbb{R}^d, \forall t \in \mathbb{R}^+, \quad (20)$$

with $U \in (C^\infty(\mathbb{R}^d))^m$ the real unknown for $m \in \mathbb{N}^*$, $M \in (C^\infty(\mathbb{R}^d))^{m \times m}$ a symmetric positive definite matrix field, and $F_j$ symmetric constant matrices of size $m$ by $m$ for $j = 1, 2, \dots, d$.
With the objective of solving complex problems, like the convected Helmholtz equation (1) in the future, this section will introduce a generalized time-harmonic version of such systems.

Many numerical methods have been developed for problems governed by Friedrichs systems, notably by Lesaint [Les73] in the context of symmetric hyperbolic systems using Galerkin methods. The Friedrichs framework also serves as the foundation of the unified approach proposed by Ern and Guermond [EG06a, EG06b] for the design and analysis of Galerkin methods. A wide range of applications of DG methods to problems governed by Friedrichs systems – covering elliptic, parabolic, hyperbolic, and mixed-type equations – can also be found in the work of Jensen [Jen04], further emphasizing the flexibility and unifying aspects of this formulation. Guided by our interest for time-harmonic regime, we turn to a time-harmonic version of Friedrichs systems.

**Definition 4.1.2.** (Time-harmonic Friedrichs' system) We call a time-harmonic Friedrichs' system a problem that can be written in the form

$$ -i\omega M(\boldsymbol{x})\boldsymbol{U}(\boldsymbol{x}) + \sum_{j=1}^{d} F_j \partial_j \boldsymbol{U}(\boldsymbol{x}) = 0, \quad \forall \boldsymbol{x} \in \Omega, \qquad (21) $$

where $\boldsymbol{U}$ is the complex unknown, $M \in (C^\infty(\bar{\Omega}))^{m \times m}$ is a symmetric positive definite matrix field, and $F_j$ are symmetric constant matrices of size $m$ by $m$ for $j = 1, 2, \dots, d$.

An important notion for the development of DG methods for Friedrichs system is the notion of flux. It relies on the so-called hyperbolic symbol, a matrix field that depends only on the space derivative terms of Friedrichs' systems.

**Definition 4.1.3.** For any vector $\boldsymbol{n} \in \mathbb{R}^d$, we denote by $F[\boldsymbol{n}]$ the matrix field of size $m$ by $m$, called the hyperbolic symbol. It is defined by

$$ F[\boldsymbol{n}] = \sum_{j=1}^{d} F_j n_j, \quad \forall \boldsymbol{n} \in \mathbb{R}^d. \qquad (22) $$

Interestingly, this notion is the same for both the time-domain and time-harmonic version of Friedrichs' systems.

While the first order system equivalent to Helmholtz equation naturally fits in this framework, it is not the case for the convected Helmholtz equation. As a consequence we introduce a new formulation inspired by Friedrichs' approach for future work on the convected Helmholtz equation.
**Definition 4.1.4** (Generalized Friedrichs' system). A Generalized Friedrichs' problem is a system that can be written in the form

$$-i\omega M(\boldsymbol{x})U(\boldsymbol{x}) + \sum_{j=1}^{d} F_j(\boldsymbol{x})\partial_j U(\boldsymbol{x}) + G(\boldsymbol{x})U(\boldsymbol{x}) = 0, \quad \forall \boldsymbol{x} \in \Omega, \quad (23)$$

with the complex unknown $U \in (C^\infty(\Omega))^m$, $m \in \mathbb{N}^*$, $M \in (C^\infty(\bar{\Omega}))^{m \times m}$ a symmetric positive definite matrix field, and $G, F_j \in (C^\infty(\Omega))^{m \times m}$ for $j \in [1, d]$.

## 4.2 Numerical hyperbolic flux

In quasi-Trefftz methods in particular and in DG methods in general, discontinuous functions are defined piecewise on each element of the mesh. The notion of their traces at the interfaces between elements is critical, and the weak formulation relies on the definition of a corresponding numerical flux relying on the hyperbolic symbol $F[n]$, required to evaluate functions at the boundaries of elements.

The proposed decomposition of the matrix field $F[n]$ into positive and negative parts represents the fluxes involved at the interfaces. This decomposition, which has also been explored by [RPT23] and [Dji25], can be interpreted in terms of energy flux when returning to the notion of work done by external forces. The decomposition will rely on a generalized eigenvalue problem with respect to the matrix $M$ defining the Friedrichs' system, see Definition 4.1.2. Indeed, $M$ being symmetric positive definite defines a metric related to total energy of the problem.

**Proposition 4.2.1.** Given $M \in (C^\infty(\bar{\Omega}))^{m \times m}$ a symmetric positive definite matrix field and $\{F_j \in \mathbb{R}^{m \times m}, j = 1, 2, ..., d\}$ symmetric constant matrices, consider the corresponding Friedrichs' system. Then for all $n \in \mathbb{R}^d$, the matrix field $F[n]$ defined in (22) can be decomposed in the form

$$F[n] = F_+[n] + F_-[n],$$

where $F_+[n]$ represents the outgoing flux operator, and $F_-[n]$ the incoming flux operator.

*Proof.* The matrices $F[n]$ and $M$ are symmetric by definition. Consider the generalized eigenvalue problem: find $w \in \mathbb{R}^m$ and $\lambda \in \mathbb{R}$ such that

$$F[n]w = \lambda Mw.$$

Then denote by $W$ the matrix of generalized eigenvectors and $\Lambda$ the diagonal matrix containing the corresponding generalized eigenvalues, which are real. In particular, they satisfy $F[n] = MW\Lambda W^{-1}$ and $W^T MW = Id$, and $F[n] = MW\Lambda W^T M$.
The decomposition then relies on splitting the positive and negative eigenvalues, and accordingly defining $F^{\pm}[n] = MW\Lambda^{\pm}W^T M$, where $\Lambda^{\pm}$ contains only the positive eigenvalues, and $\Lambda^{-}$ only the negative ones. $\square$

These flux operators can then be leveraged, first to express the boundary condition, second to define numerical fluxes.

Given $U_{inc}$ an incident wave, we can reformulate the boundary conditions (16) as follows

$$F_{-}[n]U_{T} = F_{-}[n]U_{inc}, \quad \text{on } \partial T \cap \partial \Omega, \qquad (24)$$

for all elements $T \in \mathcal{T}$. Therefore, problem (19) can then be written under the standard form: with $M \in (C^{\infty}(\bar{\Omega}))^{m \times m}$ a symmetric positive definite matrix field and $\{F_j \in \mathbb{R}^{m \times m}, j = 1, 2, \dots, d\}$ symmetric constant matrices;

$$\left\{ \begin{array}{l} -i\omega M U + \sum_{j=1}^{d} F_j \partial_j U = 0 \quad \text{in } \Omega, \\ F_{-}[n]U - F_{-}[n]U_{inc} = 0, \quad \text{on } \partial\Omega, \end{array} \right. \qquad (25)$$

To define the numerical flux on each face, we add the flux incoming from the neighboring element with the flux outgoing from the current element. This is detailed in the following definition.

**Definition 4.2.2** (Numerical hyperbolic flux). Consider a function $U: x \in \mathbb{R}^d \to U(x) \in \mathbb{C}^m$ and a mesh $\mathcal{T}$. The hyperbolic numerical flux that represents $F[n]U$ on any face $f$ is given by

$$\widehat{F[n]U} = \begin{cases} F_{+}[n]U_{T} + F_{-}[n]U_{V}, & \text{on } f = \partial T \cap \partial V, \\ F_{+}[n]U_{T} + F_{-}[n]U_{inc}, & \text{on } f = \partial T \cap \partial \Omega, \end{cases}$$

with $T$ and $V$ two neighboring elements, $n$ the outgoing normal vector of $T$, and $U_{inc}$ a given incident wave.

Based on this numerical flux, the next step is to define the DG variational formulation.

## 4.3 Discontinuous Galerkin variational formulation

Discontinuous Galerkin methods are widely used for solving partial differential equations. There are numerous explanations of these methods; including the derivation of the variational formulation and the analysis of convergence, notably for hyperbolic problems [Coc03], elliptic problems [ABCM02], and a wide range of problems such as advection-reaction and diffusion and Friedrichs' systems [DPE11].
We detail in this subsection how to derive the DG variational formulation for problems governed by a standard Friedrichs' system, before turning to the Trefftz and quasi-Trefftz formulations

In a standard way, the DG variational formulation is derived by multiplying the PDE (21) by a complex-valued test function $U' \in \mathcal{H}$ where $\mathcal{H} := \prod_{T \in \mathcal{T}} \mathcal{H}_T$ with

$$ \mathcal{H}_T := \{ \mathbf{U}'_T \in (L^2(T))^m \mid \operatorname{div}(F\mathbf{U}'_T) \in (L^2(T))^m \text{ and } F[\mathbf{n}][\mathbf{U}'_T] \in (L^2(\partial T))^m \}, $$

and where $\operatorname{div}(F\mathbf{U}')$ denotes the compact form of $\sum_{j=1}^d \partial_j(F_j \mathbf{U}'_T)$. We then integrate over a generic mesh element $T$ as follow:

$$ \int_T -i\omega M \mathbf{U}_T \cdot \overline{\mathbf{U}'_T} + \sum_{j=1}^d F_j \partial_j \mathbf{U}_T \cdot \overline{\mathbf{U}'_T} = 0. $$

By integrating by parts, we obtain

$$ \int_T -i\omega M \mathbf{U}_T \cdot \overline{\mathbf{U}'_T} - \mathbf{U}_T \cdot \overline{\sum_{j=1}^d \partial_j(F_j \mathbf{U}'_T)} + \int_{\partial T} \mathbf{U}_T \cdot \overline{F[\mathbf{n}]\mathbf{U}'_T} = 0, $$

where $\mathbf{n}$ is the outgoing normal vector of $T$, or equivalently

$$ \int_T \mathbf{U}_T \cdot \left( \overline{i\omega M \mathbf{U}'_T - \sum_{j=1}^d F_j \partial_j \mathbf{U}'_T} \right) + \int_{\partial T} F[\mathbf{n}] \mathbf{U}_T \cdot \overline{\mathbf{U}'_T} = 0. $$

Then $F[\mathbf{n}]\mathbf{U}_T$ at the element interfaces is replaced with the hyperbolic numerical flux denoted by $\overline{F[\mathbf{n}]\mathbf{U}}$, giving

$$ \int_{\partial T} \overline{F[\mathbf{n}]\mathbf{U} \cdot \overline{\mathbf{U}'_T}} = \int_{\partial T} F_+[\mathbf{n}]\mathbf{U}_T \cdot \overline{\mathbf{U}'_T} + \sum_{f \in \mathcal{F}_T \cap \mathcal{F}_{\text{int}}} \int_f F_-[\mathbf{n}]\mathbf{U}_V \cdot \overline{\mathbf{U}'_T} \\ + \sum_{f \in \mathcal{F}_T \cap \mathcal{F}_{\text{ext}}} \int_f F_-[\mathbf{n}]\mathbf{U}_{\text{inc}} \cdot \overline{\mathbf{U}'_T}. $$

The classical variational formulation is then obtained by adding over all elements $T \in \mathcal{T}$ with numerical fluxes. The resulting formulation can then be written compactly as: find $\mathbf{U} \in \mathcal{H}$ such that

$$ a(\mathbf{U}, \mathbf{U}') = \ell(\mathbf{U}'), \quad \forall \mathbf{U}' \in \mathcal{H}, \qquad (26) $$
in terms of the following sesquilinear and antilinear forms:

$$a(\mathbf{U}, \mathbf{U}') := \sum_{T \in \mathcal{T}} a_T(\mathbf{U}_T, \mathbf{U}'_T) + \sum_{T \in \mathcal{T}} \sum_{f \in \mathcal{F}_T} r_{T,f}^+(\mathbf{U}_T, \mathbf{U}'_T) + r_{T,f}^-(\mathbf{U}_V, \mathbf{U}'_T),$$

$$a_T(\mathbf{U}_T, \mathbf{U}'_T) := \int_T \mathbf{U}_T \cdot \left( i\omega M \mathbf{U}'_T - \sum_{j=1}^d F_j \partial_j \mathbf{U}'_T \right),$$

$$r_{T,f}^+(\mathbf{U}_T, \mathbf{U}'_T) := \int_f F_+[\mathbf{n}]\mathbf{U}_T \cdot \overline{\mathbf{U}'_T},$$

$$r_{T,f}^-(\mathbf{U}_V, \mathbf{U}'_T) := \begin{cases} \int_f F_-[\mathbf{n}]\mathbf{U}_V \cdot \overline{\mathbf{U}'_T}, & \text{for } f \in \mathcal{F}_{\text{int}}, f = \partial T \cap \partial V, \\ 0, & \text{for } f \in \mathcal{F}_{\text{ext}}, \end{cases}$$

$$\ell(\mathbf{U}') := \sum_{T \in \mathcal{T}} \sum_{f \in \mathcal{F}_T} \ell_{T,f}(\mathbf{U}'_T)$$

$$\ell_{T,f}(\mathbf{U}'_T) := \begin{cases} 0, & \text{for } f \in \mathcal{F}_{\text{int}}, \\ -\int_f F_-[\mathbf{n}]\mathbf{U}_{\text{inc}} \cdot \overline{\mathbf{U}'_T}, & \text{for } f \in \mathcal{F}_{\text{ext}}. \end{cases} \qquad (27)$$

Trefftz and quasi-Trefftz methods will then be defined in terms of some specific function spaces.

## 4.4 Trefftz Discontinuous Galerkin variational formulation

We first extend this formulation to Trefftz methods before moving on to quasi-Trefftz methods. In the case of self-adjoint problems, which we consider throughout this work, Trefftz methods take both trial and test functions as local solutions of the main equation. This approach eliminates all the volume terms in the variational formulation, retaining only the terms on the mesh interfaces.

**Definition 4.4.1 (Local Trefftz space).** On a given mesh element $T$, we define $\mathcal{X}_T$ as the local Trefftz space, i.e., the space of all functions $\mathbf{U}'_T \in \mathcal{H}_T$ that satisfy

$$-i\omega M \mathbf{U}'_T + \sum_{j=1}^d F_j \partial_j \mathbf{U}'_T = 0, \quad \text{on } T. \qquad (28)$$

Moreover, we define $\mathcal{X}$ as the Trefftz space, given by

$$\mathcal{X} = \prod_{T \in \mathcal{T}} \mathcal{X}_T$$
Remark 4.4.2. We observe that in the standard Friedrichs’ framework, the problem is self-adjoint, then trial and test functions both solve the same PDE.

The Trefftz Discontinuous Galerkin variational formulation is then defined from (26) by replacing the standard DG spaces by the Trefftz spaces: find $U \in \mathcal{X}$ such that

$$a^{\mathrm{TDG}}(U, U') = \ell^{\mathrm{TDG}}(U'), \quad \forall U' \in \mathcal{X}, \qquad (29)$$

with the following sesquilinear and antilinear forms

$$a_T^{\mathrm{TDG}}(U, U') = \sum_{T \in \mathcal{T}} \sum_{f \in \mathcal{F}_T} r_{T,f}^+(U_T, U'_T) + r_{T,f}^-(U_V, U'_T),$$

$$\ell^{\mathrm{TDG}}(U') = \ell(U'),$$

where $r_{T,f}^+$, $r_{T,f}^-$ and $\ell$ are defined in (27).

## 4.5 Quasi-Trefftz Discontinuous Galerkin variational formulation

The quasi-Trefftz methods are based on the same principle as the Trefftz methods, but instead of local solutions of the adjoint problem, they use local approximate solutions. This allows to handle differential equations with variable coefficients, for which exact analytical solutions are most often not available.

The quasi-Trefftz space $\mathcal{Q}$, as any DG space, is defined as a broken space:

$$\mathcal{Q} = \prod_{T \in \mathcal{T}} \mathcal{Q}_T,$$

but each local space $\mathcal{Q}_T$ consists of functions $U' \in \mathcal{H}_T$ satisfying a so-called quasi-Trefftz property:

$$-i\omega M U'_T + \sum_{j=1}^{d} F_j \partial_j U'_T \approx 0, \quad \text{on } T.$$

In previous works on quasi-Trefftz methods [IG15, IG21, IGS21, IGS24], the precise meaning of this approximation has been interpreted in the sense of a truncated Taylor expansion.

Since the quasi-Trefftz functions do not solve exactly the governing equation, the volume terms in the DG variational formulation are no longer cancelled. Consequently, the quasi-Trefftz variational formulation reads: fund $U \in \mathcal{Q}$ such that

$$a(U, U') = \ell(U'), \quad \forall U' \in \mathcal{Q}, \qquad (30)$$

with $a$ and $\ell$ defined in (27).
# 5 Application to the Helmholtz equation

In this section we consider $\Omega$ as a 2-dimensional domain with

* a constant density $\rho_0$ and a zero velocity field of the carrier flow, $v_0 = 0$, leading to the classical heterogeneous Helmholtz equation,

* equipped with Fourier-Robin boundary conditions.

This section defines the quasi-Trefftz Discontinuous variational for this problem. First it focuses on the Friedrichs' system framework, including the explicit flux decomposition. Second it turns to the quasi-Trefftz space, including the construction of a basis based on previous work [IG15, IG21, IGS21, IGS24].

## 5.1 Friedrichs' system formulation

This subsection is dedicated to determining the matrices of the standard Friedrichs' system (21) for the heterogeneous Helmholtz equation (12), and computing the associated numerical flux as introduced in subsection 4.2. As a reminder, the Helmholtz equation is equivalent to the following first order system

$$ \left\{ \begin{array}{l} -\frac{i\omega}{\rho_0 c_0^2} p + \operatorname{div}(v) = 0, \\ -i\omega\rho_0 v + \nabla p = 0. \end{array} \right. $$

This system can be written in the formalism of Friedrichs' problem presented in Definition 4.1.2 for the complex unknown $U$ defined as $U(x) = \begin{pmatrix} p(x) \\ v(x) \end{pmatrix} \in \mathbb{C}^3$ for all $x \in \Omega$, thanks to the particular choice of matrices

$$ M = \begin{pmatrix} \chi & 0 & 0 \\ 0 & \rho_0 & 0 \\ 0 & 0 & \rho_0 \end{pmatrix}, \quad F_1 = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad F_2 = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix}, $$

with the notation $\chi(x) = \frac{1}{\rho_0 c_0^2(x)}$ for all $x \in \Omega$.

The corresponding hyperbolic symbol defined in (22) is then given by

$$ F[n] = \begin{pmatrix} 0 & n^T \\ n & 0_{2 \times 2} \end{pmatrix}, \quad \forall n \in \mathbb{R}^2. \qquad (31) $$

In this case, the flux decomposition is provided explicitly in the following proposition.
**Proposition 5.1.1.** The hyperbolic symbol of the heterogeneous Helmholtz equation, given in (31), can be expressed as follows

$$F[n] = F_{+}[n] + F_{-}[n], \quad \text{with } F_{\pm}[n] = \frac{1}{2} \begin{pmatrix} \pm Y & n^T \\ n & \pm \frac{1}{Y} n n^T \end{pmatrix},$$

where $nn^T \in \mathbb{R}^{2 \times 2}$, and $Y = \sqrt{\frac{\chi}{\rho_0}}$.

*Proof.* The matrix $F[n]$, for all $n \in \mathbb{R}^2$ is decomposed following the procedure outlined in subsection 4.2. We solve the generalized eigenvalue problem:

seek $w = \begin{pmatrix} w_p \\ w_v \end{pmatrix} \in \mathbb{R}^3$ and $\lambda \in \mathbb{R}$ satisfying $F[n]w = \lambda Mw$.

or equivalently

$$\begin{cases} w_v \cdot n = \lambda \chi w_p, \\ w_p n = \lambda \rho_0 w_v. \end{cases}$$

Considering $\|n\| = 1$, we multiply the second line by $n$: combining it with the first line we have

$$w_p = w_p n \cdot n = \lambda \rho_0 w_v \cdot n = \lambda^2 \rho_0 \chi w_p.$$

Now we have two possibilities:

* either $w_p \neq 0$, which implies $\lambda^2 \rho_0 \chi = 1$. It leads to two eigenvalues:

$$\lambda_+ = \frac{1}{\sqrt{\chi \rho_0}} \quad \text{and} \quad \lambda_- = \frac{-1}{\sqrt{\chi \rho_0}}.$$

* or $w_p = 0$, which leads to $\lambda_0 = 0$.

We compute an eigenvector associated with $\lambda_+$: if $w_p = 1$, then $w_v = \frac{1}{\lambda_+ \rho_0} n = Y n$. To normalize this vector, we first compute its M-weighted norm:

$$\left\| \begin{pmatrix} 1 \\ Y n \end{pmatrix} \right\|_M = \sqrt{\left( \begin{pmatrix} 1 \\ Y n \end{pmatrix} \right)^T M \begin{pmatrix} 1 \\ Y n \end{pmatrix}} = \sqrt{2\chi}.$$

Therefore, we obtain the eigenvector: $w_+ = \frac{1}{\sqrt{2\chi}} \begin{pmatrix} 1 \\ Y n \end{pmatrix} = \begin{pmatrix} (2\chi)^{1/2} \\ (2\rho_0)^{1/2} n \end{pmatrix}$. Repeating the above steps with $\lambda_-$ gives $w_- = \frac{1}{\sqrt{2\chi}} \begin{pmatrix} 1 \\ -Y n \end{pmatrix} = \begin{pmatrix} (2\chi)^{1/2} \\ -(2\rho_0)^{1/2} n \end{pmatrix}$.
Then, for $\lambda_0$, $w_p = 0$ implies $w_v \cdot n = 0$. In two dimensions, we define $n^\perp = (-n_2, n_1)^T$ as a unit vector orthogonal to $n$. Therefore, $w_v = n^\perp$, and we compute the norm

$$ \left\| \begin{pmatrix} 0 \\ n^\perp \end{pmatrix} \right\|_M = \sqrt{\begin{pmatrix} 0 \\ n^\perp \end{pmatrix}^T M \begin{pmatrix} 0 \\ n^\perp \end{pmatrix}} = \sqrt{\rho_0}. $$

to obtain the unit eigenvector $w_0 = \frac{1}{\sqrt{\rho_0}} \begin{pmatrix} 0 \\ n^\perp \end{pmatrix}$.

Hence the matrix field $F[n]$ can be expressed as

$$ F[n] = MW\Lambda W^T M \quad \text{where} \quad W = (w_+ w_- w_0) \quad \text{and} \quad \Lambda = \begin{pmatrix} \lambda_+ & 0 & 0 \\ 0 & \lambda_- & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$

To conclude the proof, we compute the matrices $F[n]_+$ and $F[n]_-$:

$$ F_{\pm}[n] = MW\Lambda^{\pm}W^T M, \quad \text{with} \quad \Lambda_{+} = \begin{pmatrix} \lambda_{+} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \text{and} \quad \Lambda_{-} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \lambda_{-} & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$

□

Now that we have established the expressions for all the matrices, the only remaining element to obtain the well-posed variational formulation (30) is the definition of the quasi-Trefftz space $\mathcal{Q}$.

## 5.2 Quasi-Trefftz functions

We seek to define quasi-Trefftz spaces to solve the Helmholtz problem (19). As stated in subsection 4.5, quasi-Trefftz functions are local quasi-solutions of the Friedrichs' system (21). Explicit definitions of quasi-Trefftz functions have been provided in several articles [IG15, IG21, IGS21, IGS24] for linear scalar differential operators. However, in the Friedrichs' framework we look for a vector-valued unknown. Therefore, in order to obtain vector-valued quasi-Trefftz functions for the problem of interest, we use the definitions and procedures from these articles to compute functions expressed in terms of the scalar unknown $p$, that quasi-solve the Helmholtz equation (12). Then, we deduce functions involving the unknown $(p, v)$ from the system (11), which is equivalent to the Helmholtz equation. These functions will form the quasi-Trefftz space used in the Discontinuous Galerkin variational formulation.
### 5.2.1 Scalar quasi-Trefftz functions

We now give the definition of a quasi-Trefftz function that we choose to work with in the context of the heterogeneous Helmholtz equation. Indeed, we seek to solve the following equation

$$ \mathcal{L}p = 0, \quad \text{in } \Omega \subset \mathbb{R}^2, \qquad (32) $$

with $\mathcal{L}$ the linear differential operator defined as

$$ \mathcal{L}: H^1(\Omega) \to H^{-1}(\Omega), \\ p \mapsto \Delta p + \kappa^2 p, $$

where $p$ is the scalar complex-valued unknown and $\kappa \in C^\infty(\Omega)$ is a smooth variable coefficient.

As briefly mentioned in subsection 4.5, we need quasi-solutions of equation (32), i.e., smooth functions $\phi \in C^\infty(\Omega)$ satisfying the local quasi-Trefftz property:

$$ \mathcal{L}\phi \approx 0. \qquad (33) $$

There is more than one way to define this quasi-Trefftz property (33). We adopt here the definition presented in [IGS21], using Taylor expansion. Therefore, we introduce our notation for the Taylor expansion: Let $x_0 \in \Omega$ be such that $\kappa^2(x_0) \neq 0$ and let $q \in \mathbb{N}$ and then, for all $x = (x, y)$ in a neighborhood of $x_0$, we have:

$$ \kappa^2(\mathbf{x}) = \sum_{k=0}^{q-1} \kappa^{2,k}(\mathbf{x}) + O(\|\mathbf{x} - \mathbf{x}_0\|^q) \qquad (34) $$

with

$$ \forall k < q, \quad \kappa^{2,k}(\mathbf{x}) = \sum_{\ell=0}^{k} \frac{\partial_x^\ell \partial_y^{k-\ell} \kappa^2(\mathbf{x}_0)}{\ell!(k-\ell)!} (x-x_0)^\ell (y-y_0)^{k-\ell}. $$

In a more general way, we denote by $\mathcal{T}^k(\phi)$ the order-$k$ Taylor polynomial of $\phi$ centered at $x_0 \in \Omega$.

Then, we reformulate the approximation (33) into the equation:

$$ \mathcal{L}\phi = O(\|\mathbf{x} - \mathbf{x}_0\|^q), \quad \text{with } q \in \mathbb{N}, \qquad (35) $$

or equivalently

$$ \forall k < q, \quad \mathcal{T}^k(\mathcal{L}\phi) = 0, \qquad (36) $$

in a neighborhood of $x_0$. Due to the variability of the coefficient $\kappa$, the quasi-Trefftz basis functions must be local. Several types of quasi-Trefftz functions are available. We present three of them: phase-based Generalized Plane Waves (GPWs) [IG15], amplitude-based GPWs [IG21] and polynomial functions [IGS24].
Notation 3. We introduce some useful notation concerning polynomials:

• $\boldsymbol{x}_0 = (x_0, y_0)$ a point in $\Omega$,

• $\boldsymbol{X} = (X, Y) = \boldsymbol{x} - \boldsymbol{x}_0$, for all $\boldsymbol{x} = (x, y) \in \Omega$,

• $P$ is a complex polynomial of degree $dP \in \mathbb{N}$ described by its Taylor series expansion centered at point $\boldsymbol{x}_0$

$$P(\boldsymbol{x}) = \sum_{k=0}^{dP} P^k(\boldsymbol{x})$$

with $\lambda_{\ell}^k \in \mathbb{C}$ complex coefficients for all $\ell \in [0, k]$, such as

$$P^k(\boldsymbol{x}) = \sum_{\ell=0}^{k} \lambda_{\ell}^k X^{\ell} Y^{k-\ell},$$

• $\boldsymbol{\Lambda}^k \in \mathbb{C}^{k+1}$ complex vectors defined by $\boldsymbol{\Lambda}^k = (\lambda_0^k, \dots, \lambda_k^k)^T$.

Phase-based Generalized Plane Waves

Plane waves (PWs) of the form

$$pw(\boldsymbol{x}) = e^{i\boldsymbol{k}_0 \cdot \boldsymbol{X}}, \quad \forall \boldsymbol{x} \in \mathbb{R}^2 \quad \text{where} \quad \boldsymbol{k}_0 \in \mathbb{R}^2 \text{ and } |\boldsymbol{k}_0| = \kappa_0,$$

are exact solutions to (32) when the Helmholtz coefficient is constant, i.e., $\kappa^2 = \kappa_0^2$, with $\kappa_0 \in \mathbb{R}$.

The initial idea to adapt this method to variable coefficients is to build functions that mimic PWs. We begin by perturbing their phase, introducing a complex phase function: Phase-based GPWs, which can be defined for all $\boldsymbol{x} \in \mathbb{R}^2$, as

$$\phi(\boldsymbol{x}) = e^{P(\boldsymbol{x})},$$

with the notations introduced in Notation 3.

As explained in [IG15], we must choose the polynomial degree; indeed, for each polynomial $P$ associated to a quasi-Trefftz function of approximation order $q$ we have

• $\frac{(dP + 1)(dP + 2)}{2}$ coefficients of $P$ to determine,

• $\frac{q(q + 1)}{2}$ equations, to cancel the terms of degree lower than $q$ in the Taylor expansion of $\mathcal{L}\phi$.
First, we do not want this system to be overdetermined; this implies that $dP > q - 1$. Second, we seek a solvable system, i.e., $dP$ must be large enough to ensure at least one non-zero linear term in each equation. However, the larger $dP$ is, the more unknowns we have, therefore, we choose the smallest $dP$ that satisfies these conditions. By analyzing the terms involved in this system–details of which will be shown in the next proof–we arrive at $dP = q + 1$, as explained in [IGS21]. This quantity is fixed for the remainder of this paper. Now, for each quasi-Trefftz function $\phi = e^P$, we must determine $\lambda_\ell^k$ for all $k \in [0, q+1]$ and $\ell \in [0, k]$.

The procedure to compute these coefficients is presented in the following proposition.

**Proposition 5.2.1.** Functions $\phi = e^P$, characterized by the following properties, are quasi-solutions of the PDE (32).

For an angle $\theta \in [0, 2\pi]$ and $N \neq 0 \in \mathbf{C}$,

* $\lambda_0^0 = 0$,

* $(\lambda_0^1, \lambda_1^1) = N(\sin \theta, \cos \theta)$,

* for all $k \in [2, q+1]$, $\lambda_0^k = \lambda_1^k = 0$.

* for all $k = 2, \dots, q+1$, and for all $\ell = 0, \dots, k-2$,

$$ \lambda_{\ell+2}^k = \frac{\beta_\ell^k - (k-\ell)(k-\ell-1)\lambda_\ell^k}{(\ell+2)(\ell+1)}, \quad (37) $$

where

$$ \beta_{\ell}^{k} = -\kappa_{\ell}^{2,k-2} - \sum_{m=0}^{k-1} \sum_{i=\max\{0,m+l+2-k\}}^{\min\{m,\ell+2\}} \lambda_{i}^{m} \lambda_{\ell-i+2}^{k-m} i(\ell-i+2) \\ - \sum_{m=0}^{k-1} \sum_{i=\min\{0,m+\ell-k\}}^{\min\{m,\ell\}} \lambda_{i}^{m} \lambda_{\ell-i}^{k-m} (m-i)(k-m-\ell+i). \quad (38) $$

*Proof.* We seek the coefficients of $P$ to obtain a quasi-solution to the Helmholtz equation (32), that is to say we want, for a fixed $q \in \mathbb{N}$,

$$ \Delta e^P + \kappa^2 e^P = (\Delta P + \nabla P \cdot \nabla P + \kappa^2)e^P = O(\|\mathbf{x} - \mathbf{x}_0\|^q). $$

Moreover, $e^P$ is bounded and non-zero near $\mathbf{x}_0$, so we simplify this equation and obtain the truncated Taylor expansion as follows:

$$ \mathcal{T}^{q-1} \left( \Delta \left( \sum_{k=0}^{q+1} P^k \right) + \nabla \left( \sum_{k=0}^{q+1} P^k \right) \cdot \nabla \left( \sum_{k=0}^{q+1} P^k \right) + \sum_{k=0}^{q-1} \kappa^{2,k} \right) = 0. $$

Inria
To determine the polynomial coefficients, we consider each homogeneous term of degree less than or equal to $q-1$; for all $k \in [2, q+1]$ we have

$$ \Delta P^k = -\kappa^{2,k-2} - \sum_{m=1}^{k-1} \nabla(P^m) \cdot \nabla(P^{k-m}). \quad (39) $$

First, we focus on the left-hand side term: we compute the matrix associated with the Laplace operator. To do this, we need to introduce the operator $W^k$, that builds a complex polynomial $P^k$ from its coefficients contained in $\Lambda^k$. It is defined by

$$ W^k : \begin{array}{l} \mathbf{C}^{k+1} \rightarrow \mathbb{P}^k[\mathbf{C}] \\ \mathbf{\Lambda}^k \mapsto W^k \mathbf{\Lambda}^k = P^k \end{array} \quad \text{where} \quad P^k : (X,Y) \mapsto \sum_{\ell=0}^{k} \lambda_{\ell}^{k} X^{\ell} Y^{k-\ell}, \quad (40) $$

with $\mathbb{H}^k[\mathbf{C}]$ the space of complex polynomials of degree $k$. Then, we seek a matrix $T^k$ such that

$$ \Delta P^k = W^{k-2} T^k \Lambda^k. \quad (41) $$

By expending $\Delta P^k$, we obtain

$$ \begin{aligned} \Delta P^k &= \sum_{\ell=0}^{k} \lambda_{\ell}^{k} \Delta(X^{\ell} Y^{k-\ell}), \\ &= \sum_{\ell=2}^{k} \ell(\ell-1)\lambda_{\ell}^{k} X^{\ell-2} Y^{k-\ell} + \sum_{\ell=0}^{k-2} (k-\ell)(k-1-\ell)\lambda_{\ell}^{k} X^{\ell} Y^{k-2-\ell}, \\ &= \sum_{\ell=0}^{k-2} (\ell+2)(\ell+1)\lambda_{\ell+2}^{k} X^{\ell} Y^{k-\ell-2} + \sum_{\ell=0}^{k-2} (k-\ell)(k-1-\ell)\lambda_{\ell}^{k} X^{\ell} Y^{k-2-\ell}. \end{aligned} $$

Finally, we express the real matrix $T^k$ of size $(k-1) \times (k+1)$ by setting all its coefficients to 0 except:

$$ T_{\ell,\ell}^{k} = (k-\ell)(k-\ell-1) \quad \text{and} \quad T_{\ell,\ell+2}^{k} = (\ell+1)(\ell+2), \quad \forall 0 \le \ell \le k-2, \quad (42) $$

and thus we have (41).

Now that we know the matrix $T^k$ for all $k=2, \dots, d$, we consider the right-hand side of equation (39). Thus, we seek for $\beta^k$ such as

$$ W^{k-2} T^k \Lambda^k = -\kappa^{k-2} - \sum_{m=1}^{k-1} \nabla(P^m) \cdot \nabla(P^{k-m}) = W^{k-2} \beta^k, \quad \text{i.e.} \quad T^k \Lambda^k = \beta^k. $$
Moreover, we introduce the notation $\beta_{\ell}^k$: the $\ell^{th}$ component of $\beta^k$, i.e., the coefficient associated to $X^{\ell} Y^{k-2-\ell}$, and we denote $B_1^k = -\kappa^{2,k-2}$ and $B_2^k = -\sum_{m=1}^{k-1} \nabla(P^m) \cdot \nabla(P^{k-m})$. Then, we identify all the components of $\beta^k$.

• In $B_1^k$ we have the term $-\kappa_{\ell}^{2,k-2}$ defined by $\kappa_{\ell}^{2,k-2}(\boldsymbol{x}) = \sum_{\ell=0}^{k} \kappa_{\ell}^{2,k-2} X^{\ell}Y^{k-2-\ell}$.

• We develop the term $B_2^k$ which gives with $m' = k - m$

$$ 
\begin{aligned}
B_2^k &= -\sum_{m=1}^{k-1} \left( \sum_{i=0}^{m} \lambda_i^m i X^{i-1} Y^{m-i} \right) \left( \sum_{j=0}^{m'} \lambda_j^{m'} j X^{j-1} Y^{m'-j} \right) \\
&\quad -\sum_{m=1}^{k-1} \left( \sum_{i=0}^{m} \lambda_i^m (m-i) X^i Y^{m-i-1} \right) \left( \sum_{j=0}^{m'} \lambda_j^{m'} (m'-j) X^j Y^{m'-1-j} \right), \\
&= -\sum_{m=1}^{k-1} \left( \sum_{i=0}^{m} \sum_{j=0}^{m'} \lambda_i^m \lambda_j^{m'} ij X^{i+j-2} Y^{k-i-j} \right) \\
&\quad -\sum_{m=1}^{k-1} \left( \sum_{i=0}^{m} \sum_{j=0}^{m'} \lambda_i^m \lambda_j^{m'} (m-i)(m'-j) X^{i+j} Y^{k-2-i-j} \right).
\end{aligned}
$$

In order to swap the sums, we set $\lambda_{\ell}^k = 0$ if $\ell > k$ or $\ell, k < 0$ or $k > dP$, we obtain

$$ 
\begin{aligned}
B_2^k &= - \sum_{m=-\infty}^{+\infty} \left( \sum_{i=-\infty}^{+\infty} \sum_{j=-\infty}^{+\infty} \lambda_i^m \lambda_j^{m'} ij X^{i+j-2} Y^{k-i-j} \right) \\
&\quad - \sum_{m=-\infty}^{+\infty} \left( \sum_{i=-\infty}^{+\infty} \sum_{j=-\infty}^{+\infty} \lambda_i^m \lambda_j^{m'} (m-i)(m'-j) X^{i+j} Y^{k-2-i-j} \right).
\end{aligned}
$$

Then, we denote $\ell = i + j - 2$ for the first sum and $\ell' = i + j$ for the second one

$$ 
\begin{aligned}
B_2^k &= - \sum_{m=-\infty}^{+\infty} \left( \sum_{i=-\infty}^{+\infty} \sum_{\ell=-\infty}^{+\infty} \lambda_i^m \lambda_{\ell-i+2}^{m'} i(\ell-i+2) X^{\ell} Y^{k-2-\ell} \right) \\
&\quad - \sum_{m=-\infty}^{+\infty} \left( \sum_{i=-\infty}^{+\infty} \sum_{\ell'=-\infty}^{+\infty} \lambda_i^m \lambda_{\ell'-i}^{m'} (m-i)(m'-\ell'+i) X^{\ell'} Y^{k-2-\ell'} \right), \\
&= - \sum_{\ell=-\infty}^{+\infty} X^{\ell} Y^{k-2-\ell} \left( \sum_{m=-\infty}^{+\infty} \sum_{i=-\infty}^{+\infty} \lambda_i^m \lambda_{\ell-i+2}^{m'} i(\ell-i+2) \right) \\
&\quad - \sum_{\ell'=-\infty}^{+\infty} X^{\ell'} Y^{k-2-\ell'} \left( \sum_{m=-\infty}^{+\infty} \sum_{i=-\infty}^{+\infty} \lambda_i^m \lambda_{\ell'-i}^{m'} (m-i)(m'-\ell'+i) \right).
\end{aligned}
$$

Inria
Therefore, in terms of $k$ and $\ell$ the coefficient associated to $X^\ell Y^{k-\ell-2}$ is the sum of:

$$ -\sum_{m=1}^{k-1} \sum_{i=\max\{0, m+\ell+2-k\}}^{\min\{m, \ell+2\}} \lambda_i^m \lambda_{\ell-i+2}^{k-m} i(\ell - i + 2) $$

and

$$ -\sum_{m=1}^{k-1} \sum_{i=\min\{0, m+\ell-k\}}^{\min\{m, \ell\}} \lambda_i^m \lambda_{\ell-i}^{k-m} (m-i)(k-m-\ell+i). $$

Then by summing terms found in $B_1^k$ and $B_2^k$ we obtain $\beta_\ell^k$ introduced in Proposition 5.2.1.

The matrix $T^k$ is defined only for $k \in [2, dP]$. Therefore, we first fix the two missing levels:

* $\lambda_0^0 = 0$, to facilitate the computations

* $(\lambda_0^1, \lambda_1^1) = N(\sin\theta, \cos\theta)$, with $N \in \mathbb{C}$ and $\theta \in [0, 2\pi]$, in order to mimic plane waves.

Then, for each $k$ the value of $\Lambda^k$ is determined by solving the linear system $T^k \Lambda^k = \beta^k$. We follow the procedure in [IG15], and compute the coefficients of $P$ layer by layer; we compute $\Lambda^k$ by incrementing $k$ from 2 to $dP$. Therefore, to determine $\beta^k$ we need to know $\Lambda^j$ for $j < k$. Moreover, the matrix $T^k$ given by (42), solve $T^k \Lambda^k = \beta^k$ is equivalent to use the recurrence relation (37), in which we must know $\lambda_i^k$ and $\lambda_{i+1}^k$ to compute $\lambda_{i+2}^k$. That is why, for each level $k \in [2, dP]$, we simply fix the two first components of $\Lambda^k$:

* $\lambda_0^k = \lambda_1^k = 0$.

This completes the proof of Proposition 5.2.1. □

*Remark 5.2.2.* These quasi-Trefftz functions are defined locally, i.e., in a neighborhood of $x_0 \in \mathbb{R}^2$. In general we choose the value of $N$ as follows: $N = i\kappa(x_0)$ if $\kappa^2(x_0) > 0$ and $N = \kappa(x_0)$ if $\kappa^2(x_0) < 0$.

## Amplitude-based Generalized Plane Waves

Another idea to build quasi-Trefftz functions is to perturb the amplitude of a PW, instead of its phase as we have just seen. Amplitude-based GPWs can be defined for all $x \in \mathbb{R}^2$, as

$$ \phi(x) = P(x)e^{ik_0 \cdot x}, $$
with $k_0 \in \mathbb{R}^2$ such that $|k_0|^2 = \kappa^2(\mathbf{x}_0)$ and $P$ defined in Notation 3.

As for the Phase-based GPWs, we fix the degree of the polynomial to be $dP = q+1$, where $q \in \mathbb{N}$ is the order in the quasi-Trefftz property (35). Next, we determine all the coefficients of this polynomial $P$ for a given $k_0 \in \mathbb{R}^2$. The procedure to compute these coefficients is very similar to the one for phase-based GPWs.

**Proposition 5.2.3.** Functions $\phi = Pe^{ik_0 \cdot \mathbf{X}}$, characterized by the following properties, are quasi-solutions of the PDE (32).

For an angle $\theta \in [0, 2\pi]$ and $N \neq 0 \in \mathbf{C}$,

* $k_0 = (k_1, k_2) = N(\sin \theta, \cos \theta)$,

* $\lambda_0^0 = 1$,

* for all $k \in [1, q+1]$, $\lambda_0^k = \lambda_1^k = 0$.

* for all $k = 2, \dots, q+1$, and for all $\ell = 0, \dots, k-2$,

$$ \lambda_{\ell+2}^k = \frac{\beta_\ell^k - (k-\ell)(k-\ell-1)\lambda_\ell^k}{(\ell+2)(\ell+1)}, \quad (43) $$

where

$$ \begin{aligned} \beta_{\ell}^{k} = & |\mathbf{k}_{0}|^{2} \lambda_{\ell}^{k-2} - 2i\mathbf{k}_{1} \lambda_{\ell+1}^{k-1} (\ell+1) - 2i\mathbf{k}_{2} \lambda_{\ell}^{k-1} (k-1-\ell) \\ & - c - \sum_{m=0}^{k-2} \sum_{i=\max\{0, m+\ell+2-k\}}^{\min\{m, \ell\}} \kappa_{i}^{2,m} \lambda_{\ell-i}^{k-2-m}. \end{aligned} \quad (44) $$

*Proof.* We proceed in the same way as in the proof of Proposition 5.2.1.

We aim to determine coefficients of $P$ such that $\phi = Pe^{ik_0 \cdot \mathbf{X}}$ satisfies the quasi-Trefftz property (35):

$$ \Delta(\mathrm{Pe}^{ik_0 \cdot \mathbf{X}}) + \kappa^2 \mathrm{Pe}^{ik_0 \cdot \mathbf{X}} = (\Delta P + 2i\mathbf{k}_0 \cdot \nabla P - |\mathbf{k}_0|^2 P + \kappa^2 P) e^{ik_0 \cdot \mathbf{X}} = O(\|\mathbf{x} - \mathbf{x}_0\|^q). $$

Since $e^{ik_0 \cdot \mathbf{X}}$ is bounded and non-zero near $\mathbf{x}_0$, we simplify the equation to the truncated Taylor expansion as follows:

$$ \mathcal{T}^{q-1} \left( \Delta \left( \sum_{k=0}^{q+1} P^k \right) + \left( \sum_{k=0}^{q-1} \kappa^{2,k} \right) \left( \sum_{k=0}^{q+1} P^k \right) + 2i \mathbf{k}_0 \cdot \nabla \sum_{k=0}^{q+1} P^k - |\mathbf{k}_0|^2 \sum_{k=0}^{q+1} P^k \right) = 0. $$

To determine the polynomial coefficients, we consider each homogeneous term of degree less than or equal to $q-1$; for all $k \in [2, q+1]$ we have

$$ \Delta P^k = - \sum_{m=0}^{k-2} \kappa^{2,m} P^{k-2-m} - 2i \mathbf{k}_0 \cdot \nabla P^{k-1} + |\mathbf{k}_0|^2 P^{k-2}. $$
We use the matrix $T^k$ associated with the Laplace operator that we have already computed (42) in the previous proof, as well as the operator $W^k$ introduced in (40) for all $k \in \mathbb{N}$. Again, we know that we want to solve $T^k \Lambda^k = \beta^k$ with this time $\beta^k$ such that

$$ -\sum_{m=0}^{k-2} \kappa^{2,m} P^{k-2-m} - 2i\mathbf{k}_0 \cdot \nabla P^{k-1} + |\mathbf{k}_0|^2 P^{k-2} = W^{k-2} \beta^k. $$

We identify for $\ell = 0, \dots, k$, the $\ell$th component of $\beta^k$: the coefficient associated with $X^\ell Y^{k-2-\ell}$. We denote $B_1^k = |\mathbf{k}_0|^2 P^{k-2}$, $B_2^k = -2i\mathbf{k}_0 \cdot \nabla P^{k-1}$ and $B_3^k = -\sum_{m=0}^{k-2} \kappa^{2,m} P^{k-2-m}$.

* In $B_1^k$ we simply have the term $|\mathbf{k}_0|^2 \lambda_\ell^{k-2}$.

* We rewrite $B_2^k$ in order to find the coefficients multiplied by $X^\ell Y^{k-2-\ell}$:

$$ B_2^k = -2ik_1 \sum_{m=1}^{k-1} \lambda_m^{k-1} m X^{m-1} Y^{k-1-m} - 2ik_2 \sum_{m=0}^{k-2} \lambda_m^{k-1} (k-1-m) X^m Y^{k-2-m}. $$

So we have the term: $-2ik_1 \lambda_{\ell+1}^{k-1}(\ell+1) - 2ik_2 \lambda_\ell^{k-1}(k-1-\ell)$.

* Finally, we develop $B_3^k$ where we recall that $\kappa^{2,m}(\mathbf{x}) = \sum_{i=0}^m \kappa_i^{2,m} X^i Y^{m-i}$.

$$ \begin{aligned} B_3^k &= -\sum_{m=0}^{k-2} \left( \sum_{i=0}^{m} \kappa_i^{2,m} X^i Y^{m-i} \right) \left( \sum_{j=0}^{k-2-m} \lambda_j^{k-2-m} X^j Y^{k-2-m-j} \right), \\ &= -\sum_{m=0}^{k-2} \left( \sum_{i=0}^{m} \sum_{j=0}^{k-2-m} \kappa_i^{2,m} \lambda_j^{k-2-m} X^{i+j} Y^{k-2-i-j} \right). \end{aligned} $$

To swap the sums, we set $\lambda_\ell^k = 0$ if $\ell > k$ or $\ell, k < 0$ or $k > dP$, we obtain

$$ B_3^k = - \sum_{m=-\infty}^{+\infty} \left( \sum_{i=-\infty}^{+\infty} \sum_{j=-\infty}^{+\infty} \kappa_i^{2,m} \lambda_j^{k-2-m} X^{i+j} Y^{k-2-i-j} \right). $$

We then set $\ell = i + j$, which gives

$$ \left\{ \begin{aligned} B_3^k &= - \sum_{m=-\infty}^{+\infty} \left( \sum_{i=-\infty}^{+\infty} \sum_{\ell=-\infty}^{+\infty} \kappa_i^{2,m} \lambda_{\ell-i}^{k-2-m} X^\ell Y^{k-2-\ell} \right) \\ &= - \sum_{\ell=-\infty}^{+\infty} X^\ell Y^{k-2-\ell} \left( \sum_{m=-\infty}^{+\infty} \sum_{i=-\infty}^{+\infty} \kappa_i^{2,m} \lambda_{\ell-i}^{k-2-m} \right). \end{aligned} \right. $$
Finally, the term associated with $X^\ell Y^{k-2-\ell}$ in $B_3^k$ is:

$$ -\sum_{m=0}^{k-2} \sum_{i=\max\{0, m+\ell+2-k\}}^{\min\{m, \ell\}} \kappa_i^{2,m} \lambda_{\ell-i}^{k-2-m} $$

By summing all these terms we recover the expression for $\beta_\ell^k$ as stated in (44). Next, following the same approach as for the Phase-based GPWs, we first fix the values of $\Lambda^0$ and $\Lambda^1$, by taking

* $\lambda_0^0 = 1$,

* $\lambda_0^1 = \lambda_1^1 = 0$.

Then for each $k$ from 2 to $q+1$, we compute $\Lambda^k$ with the recurrence relation (43), deduced from the definition of $T^k$ (42), which is equivalent to solving $T^k\Lambda^k = \beta^k$. For this, we use coefficients of inferior degrees already computed and fix $\lambda_0^k = \lambda_1^k = 0$. $\square$

**Remark 5.2.4.** As for the Phase-based GPWs, these quasi-Trefftz functions are defined locally, i.e., in a neighborhood of $x_0 \in \mathbb{R}^2$, and the value of $N$ is generally chosen as follows: $N = i\kappa(x_0)$ if $\kappa^2(x_0) > 0$ and $N = \kappa(x_0)$ if $\kappa^2(x_0) < 0$.

## Polynomial functions

The last family of quasi-Trefftz functions is the family of polynomial functions. Indeed, quasi-Trefftz functions are not exact solutions of the problem of interest, therefore we can take a simpler form than GPWs. We seek to build functions as follows: for all $x \in \mathbb{R}^2$:

$$ \phi(x) = P(x), $$

that satisfy the quasi-Trefftz property (35) of order $q$. As for the GPWs, we choose $dP$, the degree of the polynomial $P$ equal to $q+1$. The following proposition explains how to determine all the coefficients of these functions.

**Proposition 5.2.5.** Functions $\phi = P$, characterized by the following properties, are quasi-solutions of the PDE (32).

$$ \text{Let } (i_0, j_0) \in I_0 = \{(0,0)\} \cup \{(0,1) \times [1, q+1]\} $$

* $\lambda_{i_0}^{j_0} = 1$,

* for all $(\ell, k) \neq (i_0, j_0) \in I_0$, $\lambda_\ell^k = 0$,
• for all $k = 2, \dots, q + 1$, and for all $\ell = 0, \dots, k - 2$,

$$ \lambda_{\ell+2}^k = \frac{\beta_\ell^k - (k-\ell)(k-\ell-1)\lambda_\ell^k}{(\ell+2)(\ell+1)}, \quad (45) $$

where

$$ \beta_{\ell}^{k} = - \sum_{m=0}^{k-2} \sum_{i=\max\{0, m+\ell+2-k\}}^{\min\{m, \ell\}} \kappa_{i}^{2,m} \lambda_{\ell-i}^{k-2-m}. \quad (46) $$

*Proof.* In order to find coefficients of $P$, we proceed as in the proofs of Proposition 5.2.1 and 5.2.3.

We know that $\phi = P$ satisfies the quasi-Trefftz property (35), so we have

$$ \Delta P + \kappa^2 P = O(\|\boldsymbol{x} - \boldsymbol{x}_0\|^q), $$

and equivalently

$$ \mathcal{T}^{q-1} \left( \Delta \left( \sum_{k=0}^{q+1} P^k \right) + \left( \sum_{k=0}^{q-1} \kappa^{2,k} \right) \left( \sum_{k=0}^{q+1} P^k \right) \right) = 0. $$

To determine the polynomial coefficients, we consider each homogeneous term of degree less than or equal to $q-1$; for all $k \in [2, q+1]$ we have

$$ \Delta P^k = - \sum_{m=0}^{k-2} \kappa^{2,m} P^{k-2-m}. $$

Then, the proof is immediate. Indeed, in the proof of Proposition 5.2.1 we dealt with the left-hand side by determining the matrix $T^k$ associated with the Laplace operator (42). Moreover, the right-hand side has also been handled in the term $B_3^k$ of the proof of Proposition 5.2.3. The remainder of the proof follows similarly. To obtain the recurrence formula on the polynomial coefficient, we compute them layer by layer, handling each homogeneous polynomial in order of increasing degree. In the case of polynomial functions we no longer have directional constraints. This is why we fix the coefficients in a simpler way: one of the $\lambda_{i_0}^{j_0}$, with $(i_0, j_0) \in I_0$ equals to 1, and all the others with indices in $I_0$ equal to 0. The reminding polynomial coefficients are deduced by solving the linear system obtained, i.e., equation (45). $\square$

Now that we have described how to compute three families of quasi-Trefftz functions for the scalar Helmholtz equation (12), we adapt this construction to fit within Friedrichs' framework. The next subsection presents the vector-valued quasi-Trefftz space required for our variational formulation. We will then move on to numerical results.
### 5.2.2 Vector-valued quasi-Trefftz spaces

In the previous subsections, we have described three families of quasi-Trefftz functions for the scalar Helmholtz equation (12). However, in subsection 4.5 where we presented the quasi-Trefftz Discontinuous Galerkin variational formulation, we introduced quasi-Trefftz spaces composed of vector-valued functions that approximate the solution of a first-order system. To adapt this approach to vector-valued equations, we define the quasi-Trefftz space as follows. This definition is specific to the Helmholtz equation in two dimensions, but it can easily be extended to three dimensions.

**Definition 5.2.6 (Local quasi-Trefftz polynomial space).** For $q \in \mathbb{N}$, $q \ge 1$, we denote by $Q_T^q$ the local quasi-Trefftz space defined by

$$Q_T^q := \left\{ \Phi = \begin{pmatrix} \phi_p \\ \phi_v \end{pmatrix} \in \mathbb{P}^{q+1} \times (\mathbb{P}^q)^2 \mid \mathcal{T}^{q-1}(\Delta\phi_p + \kappa^2\phi_p) = 0, \phi_v = \frac{\nabla\phi_p}{i\omega\rho_0} \right\}.$$

This definition is stated for polynomial functions, but it can also be applied to another ansatz previously presented in this section. Indeed, the local vector-valued quasi-Trefftz space can be defined using functions where $\phi_p$ is either an amplitude-based or phase-based GPW, both characterized by polynomials of degree $q+1$. To clarify, the scalar component $\phi_p$ in $\Phi$ satisfies the quasi-Trefftz property for the scalar Helmholtz equation, and the second unknown $\phi_v$ is simply deduced from the system (11) which is equivalent to the Helmholtz equation.

One of the main benefits of the quasi-Trefftz method is that it requires fewer degrees of freedom than the classical Discontinuous Galerkin method. This is because, in the scalar case, we showed that only $2dP + 1$ coefficients are fixed for polynomials of degree $dP$. Thus, this method has $2dP + 1$ degrees of freedom per element, since the vector-valued component $\phi_v$ is fully determined from $\phi_p$. By comparison, classical methods using polynomials have $(dP+1)(dP+2)/2$ degrees of freedom for the scalar unknown. We must also add $dP(dP+1)$ (in two dimensions) for the additional component $\phi_v$, which typically consists of polynomials of degree $dP-1$.

Finally, to construct basis functions in these quasi-Trefftz spaces, in the case of GPWs, we just have to take $2dP + 1$, equivalently $2q + 3$, linearly independent directions $\theta$ during the construction of quasi-Trefftz functions as described in Proposition 5.2.1 and 5.2.3, to obtain linearly independent functions. Whereas for the polynomial functions, there are no directions involved. Thus, in the construction procedure we successively select distinct indices $(i_0, j_0) \in I_0$ as described in Proposition 5.2.5.

In future work, we aim to propose a way to compute vector-valued quasi-Trefftz
functions using a similar procedure. We retain the layer-by-layer construction, but this time, we work directly with the first-order system and determine all the unknowns simultaneously. Furthermore, for the Helmholtz equation, it leads to the same quasi-Trefftz spaces.

Now that the quasi-Trefftz variational formulation has been fully established and appropriate approximation spaces have been defined, we proceed to numerical experiments. We first verify that the proposed functions satisfy the quasi-Trefftz property (35) for which they are designed for. Then, we examine the errors associated with the corresponding approximate solutions.

# 6 Numerical results

We now turn to the illustration of the implementation of the quasi-Trefftz Galerkin method for the Helmholtz problem presented in Section 5. To that end, we consider a test case with explicit exact solution throughout the domain. On the one hand, for various types of basis functions, the high-order convergence of the remainder of the quasi-Trefftz property is illustrated. On the other hand, for these bases, the quasi-Trefftz approximation to the solution of the BVP is computed and the local error is displayed.

## 6.1 Definition of the test case

Consider the heterogeneous Helmholtz equation (12): $\Delta p + \kappa^2 p = 0$, associated with the boundary condition: $\partial_n p - i\kappa p = \partial_n p_{inc} - i\kappa p_{inc}$, with $p_{inc}$ an incident wave.

The test case of interest in terms of the scalar problem is defined as follows.

* the domain of interest is $\Omega = [0, 1] \times [0, 1] \in \mathbb{R}^2$,

* the constant density is $\rho_0 = \frac{4}{\pi^2}$,

* the harmonic frequency is $\omega = \frac{\pi}{2}$,

* the coefficient is $\kappa^2(x, y) = \frac{\omega^2}{c^2(x, y)} = 24 - y^4 + 6y^2 + 6y$, for all $(x, y) \in \Omega$,

* the exact solution is $p_{ex}(x, y) = e^{5ix+y-y^3}$ for all $(x, y) \in \Omega$,

* the input for the boundary condition is $p_{inc} = p_{ex}$ for all $(x, y) \in \Omega$.
Correspondingly the solution to the first-order system formulation is

$$ \forall (x, y) \in \Omega, \text{ and hence } \mathbf{U}_{ex}(\mathbf{x}) = \begin{pmatrix} p_{ex}(\mathbf{x}) \\ \frac{\nabla p_{ex}(\mathbf{x})}{i\omega\rho_0} \end{pmatrix}. $$

## 6.2 Verification of the quasi-Trefftz Property

As a reminder, all quasi-Trefftz basis functions satisfy the local quasi-Trefftz property (35): for any basis function $\phi$ defined around a point $x_0 \in \Omega$, and with an approximation order $q \in \mathbb{N}$, then $\Delta\phi + \kappa^2\phi = O(\|x - x_0\|^q)$.

This local property is here investigated at the particular point $x_0 = (0.75, 0.25) \in \Omega$, for the three types of quasi-Trefftz functions introduced earlier. For each type of basis, we estimate the maximal residuals for all elements of a given basis. Figure 2 displays the maximal residual, defined as:

$$ res_{\infty}(r) = \max_{j \in J} \max_{\|x-x_0\|=r} \frac{|\Delta\phi_j(\mathbf{x}) + \kappa^2(\mathbf{x})\phi_j(\mathbf{x})|}{\|\phi_j\|_{\infty}}, \quad \text{with } \mathbf{x} \in \Omega, $$

where $J = [0, 2dP - 1]$ is the set of all the indices of basis functions and $r$ is the distance from the center $x_0$. We choose to study the case $dP = 5$, so we expect a residual of order $q = 4$.

<figure><img src="image_4.png" /><figcaption>Figure 2: Maximal residual of basis functions in terms of the distance from $x_0$.</figcaption></figure>
For all three families of functions, we observe the expected approximation order.

## 6.3 Quasi-Trefftz approximate solutions

In this subsection, we compare several quasi-Trefftz approximate solutions by numerically solving the Helmholtz problem using the quasi-Trefftz variational formulation (30) from Section 5. We compare three numerical solutions, each based on a different quasi-Trefftz space defined on a unique mesh:

* $U_{ph}^h = \begin{pmatrix} p_{ph}^h \\ v_{ph}^h \end{pmatrix}$ built with phase-based GPWs from $Q_{ph}$,

* $U_{am}^h = \begin{pmatrix} p_{am}^h \\ v_{am}^h \end{pmatrix}$ built with amplitude-based GPWs from $Q_{am}$,

* $U_p^h = \begin{pmatrix} p_p^h \\ v_p^h \end{pmatrix}$ built with quasi-Trefftz polynomials from $Q_p$,

where $h$ denotes the mesh size. Although the numerical solution is vector-valued, we first consider the scalar component $p^h$ to assess the accuracy of the pressure field alone, which is often of particular interest in acoustics. We then proceed to evaluate the error of the full vector field.

Figure 3 shows the relative error between the numerical and the exact solutions over the domain. This pointwise relative error is defined by:

$$err(\boldsymbol{x}) = \frac{|p^h(\boldsymbol{x}) - p_{ex}(\boldsymbol{x})|}{||p_{ex}||_{\infty}}, \quad \forall \boldsymbol{x} \in \Omega.$$

For this simulation, we set the polynomial degree of the basis functions to $dP = 5$, and the mesh size to $h = 0.1$.

<figure><img src="image_5.png" /><figcaption>Figure 3: Relative error of the approximate solution for each type of basis functions.</figcaption></figure>
We observe that the relative error remains below $10^{-5}$ for all types of basis functions. However, during these simulations, we observed that the total computational time, including the computation of basis functions, the creation of the matrices associated with the variational formulation, and the solution of the linear system, was approximately twice as fast for the phase-based GPWs as for amplitude-based GPWs, and also twice as fast for the polynomial bases compared to the phase-based GPWs. To speed up the computation, we used the Horner's method, explained in appendix A.

We now turn our attention to the full vector-valued quasi-Trefftz solution $U^h$. We focus on the $L^2$ norm–with the metric induced by the matrix $M$ from the Friedrichs' formulation (21)–to express the relative errors defined by

$$err_{L^2}^h = \frac{\|U^h - U_{ex}\|_{L^2}}{\|U_{ex}\|_{L^2}}, \quad \text{with} \quad \|U\|_{L^2} = \int_{\Omega} M U \cdot \bar{U} \, dx. \quad (47)$$

Figure 4 illustrates the relative $L^2$-error between approximate and exact solutions as a function of the ratio of $h$, as the mesh is refined, and $\lambda_{min}$, the smallest wavelength in $\Omega$. This is more representative than only considering $h$; a mesh can be either under- or over-refined depending on the problem characteristics. For this simulation we recall that the variable wave number $\kappa$ is given by $\kappa^2(x, y) = 24-9y^4+6y^2+6y$. Thus, we compute the smallest wavelength over $\Omega = [0, 1] \times [0, 1]$:

$$\lambda_{min} = \frac{2\pi}{\kappa_{max}} \approx 1.17, \quad \text{where} \quad \kappa_{max} = \max_{x \in \Omega} \kappa(x).$$

<figure><img src="image_6.png" /><figcaption>Figure 4: Relative $L^2$-error of approximate solutions as a function of $h/\lambda_{min}$.</figcaption></figure>
We observe that for each family of functions constructed from polynomials of degree $dP = 5$, we obtain a convergence rate of order 5. We can compare this result for the polynomial family with theoretical convergence of classical Discontinuous Galerkin method. For the classical method, as for the quasi-Trefftz method, when the scalar components $\phi_p$ of basis functions are polynomials of degree $dP$, then components $\phi_v$ are polynomials of degree $dP - 1$. Therefore, Discontinuous Galerkin methods give a convergence rate $dP$, see [ABCM02]. While the three quasi-Trefftz methods display the same order on this example, they require considerably fewer degrees of freedom on each element. Classical methods typically involve 21 degrees of freedom per element for polynomial of degree 5, whereas the quasi-Trefftz method uses only 11 degrees of freedom on each element of the mesh. In summary, on this example, the quasi-Trefftz method successfully captures the behavior of the Helmholtz solution using significantly fewer degrees of freedom than classical methods, while maintaining convergence rates matching the theoretical order of classical methods.

# 7 Conclusion

To summarize, after explaining the physical context of the convected Helmholtz equation (1) and forming a well-posed problem with appropriate boundary conditions, we introduced a convenient Friedrichs' formulation (21). We then presented a Discontinuous Galerkin variational formulation adapted to this framework, based on a particular choice of numerical fluxes. Next, we proposed to leverage this formulation with quasi-Trefftz spaces. These spaces can be constructed following the procedure from [IGS24], and we extended this notion to vector-valued unknowns through a first-order formulation of the problem of interest. In order to tackle problems governed by the convected Helmholtz equation rather than the standard Helmholtz equation, we plan to investigate possible Friedrichs' system formulations of the PDE. This will also require the definition of a corresponding notion of quasi-Trefftz spaces. As for the classical Helmholtz equation, there are papers about the scalar formulation of the convected equation but none, to our knowledge, addressing a vector-valued one. Finally, we will implement this method to demonstrate the practical performance and validate the theoretical predictions of high-order convergence, comparing the results with those obtained using the classical Discontinuous Galerkin method.

This study lays the foundation for a general quasi-Trefftz approach adapted to complex aeroacoustic problems. Further goals include extension to 3D and coupling with BEM to integrate the method into Airbus' computational framework.
References

[ABCM02] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems. *SIAM Journal on Numerical Analysis*, 39(5):1749–1779, January 2002. Publisher: Society for Industrial and Applied Mathematics.

[AJRT11] T. Abboud, P. Joly, J. Rodriguez, and I. Terrasse. Coupling discontinuous Galerkin methods and retarded potentials for transient wave propagation on unbounded domains. *Journal of Computational Physics*, 230(15):5877–5907, July 2011.

[BM08] A. Buffa and P. Monk. Error estimates for the ultra weak variational formulation of the Helmholtz equation. *ESAIM: Mathematical Modelling and Numerical Analysis*, 42(6):925–940, 2008.

[Bre05] H. Brezis. *Analyse fonctionnelle: théorie et applications*. Mathématiques appliquées pour le master. Dunod, Paris, 2005.

[BRT22] H. Barucq, N. Rouxelin, and S. Tordeux. Low-order Prandtl-Glauert-Lorentz based Absorbing Boundary Conditions for solving the convected Helmholtz equation with Discontinuous Galerkin methods. *Journal of Computational Physics*, 468:111450, November 2022.

[Cam07] L. M. B. C. Campos. On 36 Forms of the Acoustic Wave Equation in Potential Flows and Inhomogeneous Media. *Applied Mechanics Reviews*, 60(4):149–171, July 2007.

[CD98] O. Cessenat and B. Despres. Application of an Ultra Weak Variational Formulation of Elliptic PDEs to the Two-Dimensional Helmholtz Problem. *SIAM Journal on Numerical Analysis*, 35(1):255–299, January 1998.

[Ces96] O. Cessenat. Application d'une nouvelle formulation variationnelle aux équations d'ondes harmoniques : problèmes de Helmholtz 2D et de Maxwell 3D. These de doctorat, Paris 9, January 1996.

[CES14] F. Casenave, A. Ern, and G. Sylvand. Coupled BEM–FEM for the convected Helmholtz equation with non-uniform flow in a bounded domain. *Journal of Computational Physics*, 257:627–644, January 2014.

[Coc03] B. Cockburn. Discontinuous Galerkin methods. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte
Mathematik und Mechanik, 83(11):731–754, 2003. eprint:
https://onlinelibrary.wiley.com/doi/pdf/10.1002/zamm.200310088.

[Dji25] I. Djiba. Méthode de décomposition de domaine par méthode de Trefftz pour la propagation des ondes en géophysique. These en préparation, Pau, 2025.

[DPE11] D. A. Di Pietro and A. Ern. Mathematical Aspects of Discontinuous Galerkin Methods. Springer Science & Business Media, November 2011.

[EG06a] A. Ern and J. L. Guermond. Discontinuous Galerkin Methods for Friedrichs' Systems. I. General theory. SIAM Journal on Numerical Analysis, 44(2):753–778, January 2006. Publisher: Society for Industrial and Applied Mathematics.

[EG06b] A. Ern and J.-L. Guermond. Discontinuous Galerkin Methods for Friedrichs' Systems. Part II. Second-order Elliptic PDEs. SIAM Journal on Numerical Analysis, 44(6):2363–2388, January 2006. Publisher: Society for Industrial and Applied Mathematics.

[Fri58] K. O. Friedrichs. Symmetric positive linear differential equations. Communications on Pure and Applied Mathematics, 11(3):333–418, 1958.

[HMP16] R. Hiptmair, A. Moiola, and I. Perugia. A Survey of Trefftz Methods for the Helmholtz Equation. In Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, pages 237–279. Springer International Publishing, Cham, 2016.

[IG15] L.-M. Imbert-Gérard. Interpolation properties of generalized plane waves. Numerische Mathematik, 131(4):683–711, December 2015.

[IG21] L.-M. Imbert-Gérard. Amplitude-based Generalized Plane Waves: New Quasi-Trefftz Functions for Scalar Equations in two dimensions. SIAM Journal on Numerical Analysis, 59(3):1663–1686, January 2021.

[IGS21] L.-M. Imbert-Gérard and G. Sylvand. A roadmap for Generalized Plane Waves and their interpolation properties. Numerische Mathematik, 149(1):87–137, September 2021.

[IGS24] L.-M. Imbert-Gérard and G. Sylvand. Three types of quasi-Trefftz functions for the 3D convected Helmholtz equation: construction and
approximation properties. *IMA Journal of Numerical Analysis*, page drae060, October 2024.

[Jen04] M. Jensen. *Discontinuous Galerkin methods for Friedrichs systems with irregular solutions*. Doctoral Thesis, University of Oxford, 2004.

[Les73] P. Lesaint. Finite element methods for symmetric hyperbolic equations. *Numerische Mathematik*, 21(3):244–255, June 1973.

[Moi11] A. Moiola. *Trefftz-Discontinuous Galerkin Methods for Time-Harmonic Wave Problems*. Doctoral Thesis, ETH Zurich, 2011.

[MW99] P. Monk and D-Q. Wang. A least-squares method for the Helmholtz equation. *Computer Methods in Applied Mechanics and Engineering*, 175(1):121–136, June 1999.

[Mü54] C. Müller. On the behavior of the solutions of the differential equation $\Delta U = F(x,u)$ in the neighborhood of a point. *Communications on Pure and Applied Mathematics*, 7(3):505–515, August 1954.

[RH04] S. W. Rienstra and A. Hirschberg. *An Introduction to Acoustics*. Technische Universiteit Eindhoven, January 2004.

[RPT23] M. Rivet, S. Pernet, and S. Tordeux. Flux reconstruction method for time-harmonic linear propagation problems: 1d a priori error analysis, December 2023. HAL preprint, https://hal.science/hal-04384657.

[Smo08] C. Smorynski. *History of Mathematics: A supplement*. Springer Science & Business Media, 2008.
# A Horner's method

Horner's method is a way to efficiently evaluate polynomials by reducing the number of computational operations. It rewrites the polynomial, making calculations faster and more accurate. It's commonly used in numerical methods to improve both speed and precision in polynomial evaluations.

## A.1 Horner's method in 1D

We first introduce the technique in 1 dimension, see [Smo08].

Let $P$ be a complex polynomial in the variable $x \in \mathbb{R}$ of degree $N$, expressed as

$$P(x) = \sum_{i=0}^{N} a_i x^i, \quad \text{where } \forall i \in [0, N], a_i \in \mathbb{C}.$$

The principle of Horner's method is to reformulate the polynomial $P$ as:

$$P(x) = p^H(x) = a_0 + x(a_1 + x(a_2 + \cdots x(a_{N-1} + a_N x) \cdots)).$$

Here is a description of how the algorithm works:

**Algorithm 1:** Horner's method in 1D

**Data:** $N \ge 0$: degree of $P$, $(a_n)_n$: coefficients of $P$.

1. $P \leftarrow 0$

2. **for** $i = 0$ **to** $N$ **do**

3.   $\bar{P} \leftarrow P * x + a_{N-i}$

4.   **return** $P$

## A.2 Horner's method in 2D

Now that we have introduced the principle of this method in one dimension, we extend it to the two-dimensional case.

Let $P$ be a complex polynomial of $(x, y) \in \mathbb{R}^2$ of total degree $N$, expressed as

$$P(x, y) = \sum_{i=0}^{N} \sum_{j=0}^{N-i} a_{i,j} x^i y^j, \quad \text{with } \forall (i, j) \in \mathbb{N}^2, i + j \le N, a_{i,j} \in \mathbb{C}.$$
In order to reformulate $P$, we start by grouping the terms with respect to the powers of $x$, and then apply the 1D Horner's method to $y$:

$$P(x, y) = \sum_{i=0}^{N} x^{i} \left( \sum_{j=0}^{N-i} a_{i,j} y^{j} \right) = \sum_{i=0}^{N} x^{i} p_{i}^{H}(y),$$

with

$$p_{i}^{H}(y) = a_{i,0} + y(a_{i,1} + \cdots y(a_{i,N-i-1} + a_{i,N-i}y) \cdots).$$

Thus, we apply Horner's method to $x$ and get

$$P(x, y) = p_{0}^{H}(y) + x(p_{1}^{H}(y) + \cdots x(p_{N-1}^{H}(y) + p_{N}^{H}(y)x) \cdots).$$

As in the one-dimensional case, we present the associated algorithm :

<caption>**Algorithm 2:** Horner's method in 2D</caption>

**Data:** $N \ge 0$: degree of $P$, $(a_{n,m})_{n,m}$: coefficients of $P$.
1 &nbsp; $P \leftarrow 0$
2 &nbsp; **for** $i = 0$ **to** $N$ **do**
3 &nbsp;&nbsp;&nbsp;&nbsp; $S \leftarrow 0$
4 &nbsp;&nbsp;&nbsp;&nbsp; **for** $j = 0$ **to** $i$ **do**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor S \leftarrow S * y + a_{N-i,i-j} \rfloor$
6 &nbsp;&nbsp;&nbsp;&nbsp; $P \leftarrow P * x + S$
7 &nbsp; **return** $P$

We also adapt this method to compute the derivatives of $P$. Indeed, $P$ is a polynomial, so we know the expression for all its derivatives: for all $(k, \ell) \in \mathbb{N}^2$ such that $k + \ell \le N$:

$$\frac{\partial^{k+\ell} P}{\partial x^k \partial y^\ell}(x, y) = \sum_{i=k}^{N} \sum_{j=\ell}^{N-i} [i]_k [j]_\ell a_{i,j} x^{i-k} y^{j-\ell}, \quad \text{where } [i]_k = \frac{i!}{(i-k)!}.$$

In the same way as we did for $P$, we start by grouping the terms according to the power of $x$ and apply Horner's method to $y$

$$\frac{\partial^{k+\ell} P}{\partial x^k \partial y^\ell}(x, y) = \sum_{i=k}^{N} [i]_k x^{i-k} \left( \sum_{j=\ell}^{N-i} [j]_\ell a_{i,j} y^{j-\ell} \right) = \sum_{i=k}^{N} [i]_k x^{i-k} dp_i^H(y),$$

with

$$dp_i^H(y) = [\ell]_\ell a_{i,\ell} + y([\ell+1]_\ell a_{i,\ell+1} + \cdots y([N-i-1]_\ell a_{i,N-i-1} + [N-i]_\ell a_{i,N-i} y) \cdots).$$

Inria
Thus we apply Horner's method to $x$ and get

$$ \frac{\partial^{k+\ell} P}{\partial x^k \partial y^\ell}(x, y) = [k]_k dp_k^H(y) + x([k+1]_k dp_{k+1}^H(y) + \cdots + x([N-1]_k dp_{N-1}^H(y) + [N]_k dp_N^H(y)) \cdots). $$

To conclude we provide the associated algorithm:

**Algorithm 3:** Horner's method in 2D for the derivatives of a polynomial

**Data:** $N \ge 0$: degree of $P$, $(a_{n,m})_{n,m}$: coefficients of $P$, $k$ and $\ell$ order of the derivative according to $x$ and $y$.

1. $dP \leftarrow 0$

2. **for** $i = \ell$ **to** $N-k$ **do**

3. &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 0$ **to** $i-\ell$ **do**

4. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $k = 0$ **to** $j-i$ **do**

5. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $l = 0$ **to** $j-i$ **do**

6. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$dP \leftarrow dP * x + [N-i]_k * S$

7. **return** $dP$
informatics mathematics
Inria

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