Matthias Rivet

Sébastien Pernet

Sébastien Tordeux

Mathematical Modelling and Numerical Analysis Will be set by the publisher Modélisation Mathématique et Analyse Numérique FLUX RECONSTRUCTION METHOD FOR TIME-HARMONIC LINEAR PROPAGATION PROBLEMS: 1D A PRIORI ERROR ANALYSIS

Keywords: 2020 Mathematics Subject Classification. 65N12, 65N35, 35L02, 35Q61. 22/12/2023 Hyperbolic problem, Time-harmonic equations, Flux Reconstruction method, High-order approximation, Maxwell's equations, A priori error estimate

The Flux Reconstruction (FR) method is well established for hyperbolic equations in the Computational Fluid Dynamics field, but has barely been studied for electromagnetism. In this article, we propose to describe the FR formulation for time-harmonic linear hyperbolic problems. In particular, this formalism includes the Maxwell's equations and the unidimensional wave equations, for which the method is detailed. We then focus on the wave equations for incoming boundary conditions, and prove the well-posedness of the associated FR method and quasi-optimal a priori error estimates, which are explicit in terms of the flux correction polynomials and discretisation parameters. Numerical experiments finally validate the main behaviours of the estimates, and confirm the good properties of the method for the Maxwell problem.

Résumé. La méthode Flux Reconstruction (FR) est largement établie pour les équations hyperboliques de la mécanique des fluides, mais n'a été encore que peu étudiée dans le cadre de l'électromagnétisme. Dans cet article, nous nous proposons donc de décrire la formulation FR pour des problèmes hyperboliques linéaires harmoniques. Ce formalisme général contient en particulier le problème de Maxwell et l'équation des ondes unidimensionnelle, pour lesquelles la méthode est détaillée. Nous nous concentrons ensuite sur l'équation des ondes pour des conditions de bord entrantes, et prouvons le caractère bien posé de la méthode FR associée, ainsi que des estimations d'erreur a priori quasi-optimales et explicites en les polynômes de correction de flux et les paramètres de discrétisation. Finalement, des expériences numériques permettent de valider les principaux comportements de ces estimations, et confirment les bonnes propriétés de la méthode pour les équations de Maxwell.

Introduction

Electromagnetic waves propagation is a fundamental physical phenomenon at the heart of many science and engineering fields: their numerical simulation therefore represents a considerable stake and has led to the development and adaptation of numerous methods. In particular, the Finite Difference Method (FDM) [START_REF] Gary | Higher-Order Numerical Methods for Transient Wave Equations[END_REF][START_REF] Taflove | Computational electrodynamics: the finite-difference time-domain method[END_REF][START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF] is the first approach which has been introduced to deal with Maxwell's equations. Then, the search of highorder methods gave rise to the Finite Element Method (FEM) [START_REF] Ainsworth | Dispersive and Dissipative Properties of Discontinuous Galerkin Finite Element Methods for the Second-Order Wave Equation[END_REF][START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF][START_REF] Cohen | Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations[END_REF][START_REF] Gary | Higher-Order Numerical Methods for Transient Wave Equations[END_REF][START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF][START_REF] Jin | The finite element method in electromagnetics[END_REF][START_REF] Monk | Finite element methods for Maxwell's equations[END_REF], the Discontinuous Galerkin (DG) [START_REF] Cohen | Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations[END_REF][START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF][START_REF] Houston | Mixed discontinuous Galerkin approximation of the Maxwell operator: The indefinite case[END_REF][START_REF] Perugia | The hp-local discontinuous Galerkin method for low-frequency time-harmonic Maxwell equations[END_REF] or the Boundary Element Method (BEM) [START_REF] Buffa | Boundary element methods for Maxwell's equations on non-smooth domains[END_REF][START_REF] Sauter | Boundary Element Methods[END_REF] for example: these classic approaches are thus widely used in both academic and industrial contexts. This inexhaustive list of methods and articles embodies the broad spectrum of strategies which have been considered to face time-harmonic problems, with two main point of views: weak-formulations for DG and FEM in opposition to strong ones for FDM.

At the same time, the Computational Fluid Dynamics (CFD) community particularly focused on Finite Volume Method (FVM), which are commonly used in industrial configurations at relatively low order. The recent increase in computational abilities has also elicited interest for high-order methods. This tends to accelerate the spreading of the Flux Reconstruction (FR) method, which is the topic of this article, and of the DG ones, which have intensively been studied for wave propagation problems. The FR method was first introduced by Huynh [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF] and relies on the strong formulation with a hyperbolic point of view: the solution and flux approximations are searched as piecewise polynomial functions of degree k and k `1 respectively. To take Boundary Conditions (BCs) and continuity properties into account, the flux variables result from the use of correction polynomial functions. Such high-order methods relying on the strong formulation had particularly gained interest with the introduction of the staggered-grid method [START_REF] Kopriva | A Conservative Staggered-Grid Chebyshev Multidomain Method for Compressible Flows[END_REF], which was then extended to triangular cells by Liu et al. [START_REF] Liu | Spectral difference method for unstructured grids I: Basic formulation[END_REF] and renamed Spectral Difference (SD). Such a generalisation for FR has been introduced in [START_REF] Wang | A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids[END_REF] through a method named Lifting Collocation Penalty: both methods were lately gathered under the denomination of Correction Procedure via Reconstruction (CPR). As this article will only consider the classic framework of FR, we will stick to this designation in what follows. Further developments were then led for Euler and Navier-Stokes equations in 2 and 3 dimensions on various mesh types. More generally, a class of energy-stable schemes, referred as Energy Stable Flux Reconstruction (ESFR) or Vincent-Castonguay-Jameson-Huynh (VCJH), tends to bring a pertinent approach to the construction of high-order stable methods for these time-dependent problems [START_REF] Jameson | A Proof of the Stability of the Spectral Difference Method for All Orders of Accuracy[END_REF][START_REF] Vincent | A New Class of High-Order Energy Stable Flux Reconstruction Schemes[END_REF]. Finally, entropy stability considerations start to be taken into account in recent studies [START_REF] Marchal | Entropy-Stable Spectral Difference and Flux Reconstruction Methods for Discontinuous Flows[END_REF]. The reader may consult the reviews [START_REF] Huynh | High-Order Methods for Computational Fluid Dynamics: A Brief Review of Compact Differential Formulation on Unstructured Grids[END_REF][START_REF] Wang | A review of flux reconstruction or correction procedure via reconstruction method for the Navier-Stokes equations[END_REF] for an overview of the main developments of the method up to 2015. In terms of noteworthy properties, the FR method relies on a strong formulation which does not depend on any quadrature rule and offers a generic way to create high-order accurate methods (by choosing different definitions for the flux correction polynomial functions). Moreover, classic methods can be retrieved with this construction scheme in usually less sophisticated and expensive ways, such as particular SD and nodal DG methods [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF][START_REF] Wang | A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids[END_REF]: their good sparsity and conditioning properties can then be obtained in this framework. In particular, the close relationships between a class of FR and filtered DG have been investigated by Allaneau and Jameson [START_REF] Allaneau | Connections between the filtered discontinuous Galerkin method and the flux reconstruction approach to high order discretizations[END_REF], offering an alternative view of the stability criterion of Vincent et al [START_REF] Vincent | A New Class of High-Order Energy Stable Flux Reconstruction Schemes[END_REF]. Last but not least, this formalism allows an easy understanding of the method effects (by being able to compute the flux correction polynomial functions and their properties) while intrinsically respecting the hyperbolic structure of the equation.

The objective is then to adapt the FR approach to time-harmonic wave propagation problems. Indeed, due to the non-linear properties of CFD equations, most of the FR developments have been realised in a timedependent context for transport, Euler and Navier-Stokes equations: to our knowledge, only a few adaptations to the transient Maxwell equations have already been realised for transverse electric and magnetic waves [START_REF] Hazra | Globally constraint-preserving FR/DG scheme for Maxwell's equations at all orders[END_REF][START_REF] Liu | Spectral difference method for unstructured grids I: Basic formulation[END_REF], and in the specific case of SD [START_REF] Ritzenthaler | Spectral Difference method On Structured-Grids for Maxwell's Equations in Time Domain[END_REF]. We will then detail the method in a general framework of time-harmonic linear hyperbolic problems, which includes acoustic, elastic and electromagnetic wave equations. In particular, this paper deals with both Maxwell's equations and the 1D wave equations, for a Cartesian mesh. The question of the well-posedness of the associated FR methods then naturally arises, which leads us to derive such conditions for a specific case of the 1D wave equations, by leaning on the uncoupling property of one-way variables. Moreover, while most of the theoretical study of the FR methods relies on Fourier (Von Neumann) analysis, we derive a priori error estimates in Sobolev norms in this framework to explicit the impact of flux correction polynomial functions on the method.

Introduction of the problems of interest

First, we introduce the two problems on which we will focus and detail the FR approach.

The homogeneous time-harmonic 1D wave equations

The first model problem we consider in this article is the 1D time-harmonic wave problem written in its hyperbolic form. Its interest is purely theoretical and does not lie in its discretisation.

Problem 1.1 (Time-harmonic 1D wave problem). Let Ω :" r0, Ls Ă R be an interval and κ ą 0 be the wavenumber. Find y " pu, vq T P " H 1 pΩq ‰ 2 such that iκy `dϕ dx " 0 with ϕ " Fpyq in Ω, [START_REF] Ainsworth | Dispersive and Dissipative Properties of Discontinuous Galerkin Finite Element Methods for the Second-Order Wave Equation[END_REF] equipped with Fourier-Robin Boundary Conditions (BCs) written as up0q ´Z1 vp0q " g 1 and upLq `Z2 vpLq " g 2 ,

with the linear flux operator F " ˆ0 ´1 ´1 0 ˙, Z 1 , Z 2 P tZ P C : ℜpZq ą 0u, and pg 1 , g 2 q P C 2 .

Remark 1.2. The simplest way to obtain this system is to consider the time-harmonic solution Y " pU, V q, with U px, tq " ℜ `upxqe iωt ˘and V px, tq " ℜ `upxqe iωt ˘, of the d'Alembert equation written in a hyperbolic way

BV
equipped with the Fourier-Robin boundary conditions U p0, tq ´Z1 V p0, tq " ℜ pg 1 exppiωtqq and U pL, tq `Z2 V pL, tq " ℜ pg 2 exppiωtqq ,

with ω ą 0 the angular frequency and c ą 0 the velocity, such that κ " ω c .

with the determinant δ " 2 cospκLqpZ 1 `Z2 q `2i sinpκLqp1 `Z1 Z 2 q ‰ 0. This decomposition will find its usefulness in the estimates proof of Section 4.

1.1.2. The time-harmonic 3D Maxwell problem Then, we introduce the time-harmonic 3D Maxwell problem, as a more complex problem (not being in one dimension anymore will call on a few subtleties): its importance for industrial problematics justifies the design of adapted FR methods. Let's consider an orthogonal parallelepiped domain Ω Ă R 3 , whose boundary is denoted as BΩ. The vacuum permittivity and permeability are denoted as ϵ 0 and µ 0 , while the relative ones are respectively ϵ r and µ r , and we suppose they verify the conditions outlined in Section 4.2 of [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF]. Then, in the absence of charges and currents, the Maxwell problem for the electric and magnetic fields E and H takes the form @ px, tq P Ω ˆR`,

$ ' & ' % ∇ ˆEpx, tq " ´BB Bt px, tq, ∇ ˆHpx, tq " BD Bt px, tq, (6)
with, in the case of a linear and isotropic material,

D " ϵ 0 ϵ r E and B " µ 0 µ r H. (7)
We obtain the time-harmonic version by looking for E and H as

Epx, tq " E 0 ℜ `epxq e iωt ˘and Hpx, tq " H 0 ℜ `hpxq e iωt ˘, (8)
denoting their amplitudes as E 0 and H 0 " c ϵ 0 µ 0 E 0 , the speed of light in a vacuum as c 0 " pϵ 0 µ 0 q ´1{2 and the wavenumber as κ " ω c 0 . This finally leads to the following normalised time-harmonic Maxwell equations that we equip with impedance BCs.

Problem 1.4 (Time-harmonic Maxwell problem). Find the electromagnetic field E :" pe, hq P "

L 2 pΩq ‰ 6 s.t. $ ' ' ' ' & ' ' ' ' % iκME `Bϕ x Bx `Bϕ y By `Bϕ z Bz " 0 in Ω, ϕ x " F x E, ϕ y " F y E and ϕ z " F z E in Ω, γ t rn BΩ se `ZBΩ γ ˆrn BΩ sh " g on BΩ, (9)
by denoting

E " pe x e y , e z , h x , h y , h z q T , M " " ε r I 3 0 3 0 3 µ r I 3 ȷ , F x " " 0 3 ´Γˆr e x s Γ ˆre x s 0 3 ȷ , F y " " 0 3 ´Γˆr e y s Γ ˆre y s 0 3 ȷ and F z " " 0 3 ´Γˆr e z s Γ ˆre z s 0 3 ȷ , (10)
where Z BΩ P L 8 pBΩq denotes the positive boundary impedance that we suppose constant on each face, g P " L 2 pBΩq ‰ 3 and I 3 and 0 3 stand for the identity and null matrices of C 3,3 . The tangential component and trace of w are denoted on BΩ as γ t rn BΩ sw " w ´pw ¨nBΩ q n BΩ " pn BΩ ˆwq ˆnBΩ and γ ˆrn BΩ sw " n BΩ ˆw, [START_REF] Ern | An Intrinsic Criterion for the Bijectivity of Hilbert Operators Related to Friedrich' Systems[END_REF] where n BΩ stands for the outgoing normal from the domain Ω.

The matrices associated to these operators for a vector n " pn 1 , n 2 , n 3 q are denoted as Γ t rns and Γ ˆrns with

Γ ˆrns " » - 0 ´n3 n 2 n 3 0 ´n1 ´n2 n 1 0 fi fl and Γ t rns " ´pΓ ˆrnsq 2 " » - pn 2 q 2 `pn 3 q 2 ´n1 n 2 ´n1 n 3 ´n1 n 2 pn 1 q 2 `pn 3 q 2 ´n2 n 3 ´n1 n 3 ´n2 n 3 pn 1 q 2 `pn 2 q 2 fi fl , (12)
such that one has γ t rnsw " Γ t rnsw and γ ˆrnsw " Γ ˆrnsw.

Theorem 1.5 (Monk). Problem 1.4 is well-posed for g P L 2

t pBΩq, with

L 2 t pBΩq " ! w P " L 2 pBΩq ‰ 3 : w ¨nBΩ " 0 ) . ( 13
)
Proof. By referring to Theorem 4.17 of [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF], Ω, ϵ r , µ r and Z BΩ all satisfy its conditions of application, which ensures the well-posedness of the problem. □

The abstract hyperbolic setting

After developing the two specific systems on which we will specify the FR method, we propose a framework for general linear hyperbolic problems. Let p, q P N ˚and Ω Ă R p be a p-hyperrectangular domain, whose 2p faces are denoted as pBΩ s q sP 1,2p .

For a given Right Hand Side (RHS) ψ P " L 2 pΩq ‰ q and κ ą 0, we consider a time-harmonic hyperbolic system iκMy `p ÿ j"1

Bϕ j Bx j " ψ in Ω, (14)
governing the primal variable y : Ω Ñ C q related to dual variables ϕ j : Ω Ñ C q , for j P 1, p , by the flux operator F j P R q,q as ϕ j " F j y.

(15) In this paper, we suppose that M P R q,q denotes a symmetric positive definite metric operator and that the F j are all symmetric, which implies the classical hyperbolic hypothesis that Frξs " ξ 1 F 1 `¨¨¨`ξ p F p is diagonizable in R for any ξ P R p . We assume this equation is equipped with boundary conditions

By " By BC , (16)
for boundary operator B : BΩ Ñ C q,q and field y BC P " L 2 pBΩq ‰ q , which finally leads to the frequency problem:

Problem 1.6 (p-dimensional time-harmonic hyperbolic problem). Find y P " L 2 pΩq ‰ q such that we have

$ ' & ' % iκMy `p ÿ j"1
Bϕ j Bx j " ψ in Ω with ϕ j " F j y, @ j P 1, p , By " By BC on BΩ.

(

) 17
Remark 1.7. We highlight the fact this part does not aim at giving general conditions on M, pF j q jP 1,p and B to ensure the well-posedness of Problem 1.6 (even if it has specifically been done for both Problems 1.1 and 1.4). This general presentation presupposes they are all given (as well as flux correction polynomial functions and a numerical trace computation definition, as detailed in Subsection 2.3), and it will introduce the associated FR formulation.

Admissible boundary conditions through flux operator decomposition

In this subsection, we precise some generic boundary conditions leading to a uniquely solvable Problem 1.6.

Definition 1.8. For a boundary face of normal n BΩ , let's introduce a symmetric positive definite surface metric Ă M P R q,q and consider the generalised eigenvalue problem

Frn BΩ sw " λ Ă Mw, (18)
where we recall that Frξs " ξ 1 F 1 `¨¨¨`ξ p F p for any ξ P R p . This leads to the decomposition

Frn BΩ s " Ă MPΛP ´1, (19)
for a real diagonal matrix Λ P R q,q of eigenvalues, and a matrix P P R q,q of eigenvector coordinates. Then, we can introduce the decomposition of Frn BΩ s into positive and negative parts as

Frn BΩ s " F Ñ Ă M rn BΩ s `FÐ Ă M rn BΩ s with F Ñ Ă M rn BΩ s " Ă MPΛ `P´1 and F Ð Ă M rn BΩ s " Ă MPΛ ´P´1 , ( 20
)
where Λ `(respectively Λ ´) denotes the positive (respectively strictly negative) part of Λ.

Thus, most of the BCs can be written in terms of the incoming flux as

By " F Ð Ă M rn BΩ sy " By BC , (21)
which will be illustrated in Section 3 for our two model problems.

Remark 1.9. In particular, the definition B " F Ð Ă M rn BΩ s falls within the framework of the admissible BCs for Friedrichs systems [START_REF] Ern | An Intrinsic Criterion for the Bijectivity of Hilbert Operators Related to Friedrich' Systems[END_REF][START_REF] Friedrichs | Symmetric positive linear differential equations[END_REF], and it is possible to show the solution uniqueness of the associated Problem 1.6 by leaning on the unique continuation theorem.

FR Method principle: a strong polynomial PDE

After defining the general hyperbolic framework, one naturally would like to introduce the discrete counterpart of [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF]. The approximations y h and ϕ j h of the primal and dual variables are then defined, for a given mesh, by: Find y h P V h,k such that we have in every cell of the mesh

iκMy h `p ÿ j"1 Bϕ j h Bx j " ψ h with ϕ j h " F j y h , @ j P 1, p , (22)
where V h,k stands for a piecewise discontinuous polynomial function space defined on the mesh, and ψ h is a piecewise polynomial approximation of the RHS ψ. This formulation is unfortunately ill-posed since it contains neither information on the solution continuity nor the boundary conditions. The idea of the FR method (as initially introduced by Huynh [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF]) is to introduce a non-trivial extension r F j of the flux operator F j to piecewise polynomial functions. This leads to the FR formulation: Problem 2.1 (Discrete problem). Find y h P V h,k such that we have in every cell of the mesh

iκMy h `p ÿ j"1 B r ϕ j h Bx j " ψ h with r ϕ j h " r F j y h @ j P 1, p . (23)
Thus, this section is devoted to the introduction of the mesh, the approximation space V h,k and the flux operator extension r F j .

X 1 0 X 1 n1-2 X 1 n1-1 X 1 n1 X 1 n1+1 X 1 N1 X 2 0 X 2 n2-1 X 2 n2 X 2 N2 K n 1,-1 K n K n 1,1 I 2 n2 I 1 n1-1 I 1 n1 I 1 n1+1 (i) Main mesh characteristics X 1 0 X 1 n1-2 X 1 n1-1 X 1 n1 X 1 n1+1 X 1 N1 X 2 0 X 2 n2-1 X 2 n2 X 2 N2 S 1 n2 K n 1,-1 K n K n 1,1 S 2 n1-1 F 1 (n1+1,n2) F 2 (n1+1,N2)
(ii) Faces and strips Figure 1. Mesh introduction in the case p " 2.

Mesh characterisation

We introduce a p-hyperrectangular domain Ω Ă R p and a direction-wise uniform mesh (see Figure 1i) whose cells K n are numbered by a multiindex n " pn j q jP 1,p P N with

Ω " p ź j"1 r0, L j s, K n " p ź j"1 I j nj , N " p ź j"1 1, N j , I j nj " rX j nj ´1, X j nj s, X j n " nh j , h j " L j N j . ( 24
)
Definition 2.2. We introduce the set F j " ´Fj n F ¯nF PN j F of faces which are orthogonal to e j (see Figure 1ii)

F j n F " I 1 n F 1 ˆ¨¨¨ˆI j´1 n F j´1 ˆ!X j n F j ) ˆIj`1 n F j`1 ˆ¨¨¨ˆI p n F p , (25)
being numbered by n F P N j F with

N j F " 1, N 1 ˆ¨¨¨ˆ 1, N j´1 ˆ 0, N j ˆ 1, N j`1 ˆ¨¨¨ˆ 1, N p . (26)
Then, in order to shorten the notations, we also introduce dedicated notations to specific behaviours relative to a chosen direction j P 1, p . Definition 2.3 (Direction-wise notations). For all a, b P N ˚, a, b P rN ˚sp and x " px 1 , . . . , x p q P R p , we denote

x zj " px 1 , . . . , x j´1 , x j`1 , . . . , x p q P R p´1 , a, b zj " a, b ztju, a, b "

ź j 1 P 1,p a j 1 , b j 1 Ă N p , and a, b zj " ź j 1 P 1,p zj a j 1 , b j 1 Ă N p´1 . (27)
Finally, the following definition makes it easy to refer to the cells adjacent to the one of interest.

Definition 2.4. For any c " pc 1 , . . . , c p q P N p , let c j,˘1 denote the shifted vector along direction j c j,˘1 " pc 1 , . . . , c j´1 , c j ˘1, c j`1 , . . . , c p q. (28)

Primal and dual variable approximations

In this subsection, we define the function spaces for the approximation of the primal variable y P " L 2 pΩq ‰ q and the dual variables ϕ j P " L 2 pΩq ‰ q satisfying (17), i.e.

p ÿ j"1 Bϕ j Bx j " ψ ´iκMy P " L 2 pΩq ‰ q . ( 29
)
The primal variable y is approximated by a piecewise multivariate polynomial function y h P V h,k .

Definition 2.5 (Primal approximation space). We introduce the discrete primal approximation space (in 1 and q dimensions) of complex valued piecewise multivariate polynomial functions of degree at most k " pk j q jP 1,p P N p in each variable

V h,k " ␣ w P L 2 pΩq : @ n P N , w| Kn P Q k pK n q ( and V h,k " rV h,k s q , (30)
where

Q k pK n q " $ & % w P L 2 pK n q : wpxq " ÿ mP 0,k α m x m , α m P C , .
-

with x m " ź jP 1,p px j q mj . (31)
Definition 2.6. We also introduce the space H 1 pT h q of piecewise H 1 functions on the mesh, which will give us a general framework for further estimates H 1 pT h q " ␣ w P L 2 pΩq : @ n P N , w| Kn P H

1 pK n q ( . (32)
This allows, for any n P N and w P " H 1 pT h q ‰ q , to denote its restriction to K n as w n :" w| Kn P " H 1 pK n q ‰ q .

Since p ÿ j"1

Bϕ j Bx j P " L 2 pΩq ‰ q , it is natural to look for the dual variable approximation r ϕ j h as a continuous function in the direction e j . This introduces the space C 0 j pΩq "

" w P L 2 pΩq : @ n zj P 1, N zj , w| S j n zj P C 0 ´Sj n zj ¯* (33)
of functions which are continuous on every strip S j n zj of direction j defined by (see Figure 1ii)

S j n zj " I 1 n1 ˆ¨¨¨ˆI j´1 nj´1 ˆr0, L j s ˆIj`1 nj`1 ˆ¨¨¨ˆI p np . (34)
Thus, with the aim of having all the terms of (23) with the same polynomial degree, each dual variable ϕ j will be approximated by r ϕ j h P W j h,k , with:

Definition 2.7 (Dual approximation spaces). For a direction j P 1, p , we introduce the discrete dual approximation space of complex valued piecewise multivariate polynomial functions of degree at most k j,1 P N p in each variable, and continuous in direction e j W j h,k " V h,k j,1 X " C 0 j pΩq ‰ q with k j,1 " pk 0 , . . . , k j´1 , k j `1, k j`1 , . . . , k p q P N p .

Finally, the idea of the FR method is to approximate the solution (primal variable), the fluxes (dual variable) and the RHS by piecewise polynomials which solve the initial PDE on each cell K n , for n P N

iκMy h,n `p ÿ j"1 B r ϕ j h,n Bx j " pΠ h,k ψq n , (36)
where Π h,k : " L 2 pΩq ‰ q Ñ V h,k denotes a projection from " L 2 pΩq ‰ q to the discrete space V h,k , whose choice belongs to the user (such as the L 2 -projection or an interpolation on a set of points for example).

Local approximation of the fluxes

The flux operator F j is a local linear operator. It can naturally be extended to V h,k and associates to the approximation y h of the primal variable an approximation ϕ j h of the dual variable given by

ϕ j h " F j y h . (37)
Yet, the FR formulation relies on a nontrivial and nonlocal extension r

F j : V h,k Ñ W j h,k , y h Þ Ñ r ϕ j h of the flux operator F j which is obtained by perturbing it as r F j " F j `δF j . ( 38
)
through a polynomial correction of the flux ϕ j h at the interfaces which will be detailed in this Subsection. Moreover, so as to ensure the consistency of the method, we will only consider perturbations which vanish for the exact solution, meaning that δF j pyq " 0 for the exact solution y.

First, as y h and ϕ j h may be discontinuous across interfaces, we need to introduce their 'direct' traces on any face from the neighbouring cells (two cells in the case of an internal interface, and one otherwise), as illustrated in Figure 2. Definition 2.8 (Direct traces). Let F P F j be an internal face orthogonal to e j such that F " K n j,´1 X K n . The direct traces of y h and ϕ j h on F (respectively 'just after' and 'just before' it) are defined by

@ x P F, $ ' & ' %
y h,F px zj q " y h,n pxq and ϕ j,h ,F " F j " y h,F ı , y h,F px zj q " y h,n j,´1 pxq and ϕ j,h ,F "

F j " y h,F ı . (40)
In particular, y h,F , ϕ j,h ,F P " Q k zj pFq ‰ q , and this definition naturally extends to boundary interfaces.

Then, as we want to impose the BCs and the continuity of the flux approximation along the interfaces, we introduce flux numerical traces on them: the flux on both sides of the faces will be corrected to fit with these values and thus be single-valued on these interfaces (as illustrated in Figure 2). Therefore, Definitions 2.9 and 2.11 will introduce generic ways to define such flux numerical traces on internal and boundary faces, even if absolutely general methods could be used.

First, for an internal face F P F j , we presuppose the reader is able to give a decomposition of F j that will divide the influence of both sides in terms of the flux numerical trace. Definition 2.9 (Internal flux numerical traces). Let F P F j be an internal face orthogonal to e j . Assume the following decomposition of the directional flux operator: F j " F j,Ñ `Fj,Ð . The numerical trace γ F r ϕ j h of the

0 1 2 0 1 2 x 1 h 1 x 2 h 2 ϕ 1 h ϕ 1,h ,F 1 p1,2q ϕ 1,h ,F 1 p1,2q γ F 1 p1,2q r ϕ 1 h Figure 2.
Illustration of the piecewise discontinuous flux in direction 1, ϕ 1 h " F 1 y h , with its direct and numerical traces, for p " 2 and k " p2, 1q. reconstructed flux r ϕ j h on F is defined by:

γ F r ϕ j h " F j,Ñ " y h,F ı `Fj,Ð " y h,F ı , (41)
where y h,F were introduced in Definition 2.8. In particular, γ F r

ϕ j h P " Q k zj pFq ‰ q .
Remark 2.10. As detailed in Subsection 1.3 and illustrated in Section 3, such a decomposition can be obtained by solving a generalised eigenvalue problem for a given symmetric positive definite matrix Ă M P R q,q , taking F j,Ñ and F j,Ð as the associated positive and negative parts of F j (see [START_REF] Houston | hp-Discontinuous Galerkin finite element methods for hyperbolic problems: error analysis and adaptivity[END_REF] for example). Moreover, one classically makes the choice Ă M " M for the internal interfaces.

In the same idea, such a numerical trace has to be defined on the boundary faces to take the BCs into account. To do so, we presuppose the reader is able to express the incoming flux with respect to a given surface metric.

Definition 2.11 (Boundary flux numerical traces). Let us consider a direction j P 1, p , a boundary face F P F j orthogonal to e j , with n P N the index of its neighbouring cell and s P 1, 2p the index of the corresponding boundary interface, such that F " K n X BΩ s . Moreover, we denote the incoming and outgoing directions as id and od, which stand for Ñ and Ð (respectively Ð and Ñ) if n j " 1 (respectively n j " N j ). Furthermore, is stands for the internal side of the boundary, corresponding to `(respectively ´) if n j " 1 (respectively n j " N j ). Finally, we suppose that, for this face F, we are given a decomposition of the directional flux operator F j " F j,Ñ `Fj,Ð , such that the BC can be written as

F j,id y " g F on F, (42)
where g F : F Ñ im pF j,id q. This allows to define a flux numerical trace polynomial γ F r ϕ j h P " Q k zj pFq ‰ q on F in the same way as [START_REF] Van Den Abeele | On the Stability and Accuracy of the Spectral Difference Method[END_REF], by replacing the incoming (and then undefined) component thanks to the previous expression:

γ F r ϕ j h " F j,od " y is h,F ‰ `ΠF h,k zj " g F ‰ , (43)
where Π F h,k zj :

" L 2 pFq ‰ q Ñ " Q k zj pFq
‰ q stands for a given projection on multivariate polynomial functions defined on F. Remark 2.12. We highlight the fact that Definitions 2.9 and 2.11 offer specific ways to define the flux numerical traces in the consistency framework of (39): for both the equations of Subsection 1.1, their applications will be detailed in Section 3. Yet, they are not the only consistent possibilities, let alone the only choice to define a FR method. Indeed, one may choose completely general ways to define γ F r ϕ j h on internal and boundary faces (not even ensuring consistency), and Definition 2.13 will still allow to define the associated FR formulation. Finally, in any case, neither well-posedness nor convergence properties of the FR method are naturally ensured: this paper will prove them only for the specific Problem 4.4, and it is not destined to characterise a class of flux numerical traces that should be used for a general problem. This finally allows to define the flux approximation r ϕ j h . To do so, we impose the flux numerical traces at each interface thanks to flux correction polynomial functions of degree k j `1: this will satisfy the expected continuity and polynomial degree properties. Definition 2.13 (Flux approximation). For j P 1, p , let's introduce two polynomials P j,Ñ , P j,Ð P Q kj `1pr0, 1sq such that " P j,Ñ p0q " 1, P j,Ñ p1q " 0, P j,Ð p0q " 0, P j,Ð p1q " 1.

We define r ϕ j h as the perturbation of ϕ j h " F j y h P V h,k whose values at any interface F P F j of normal e j match the numerical trace polynomial functions γ F r ϕ j h , thanks to corrections realised by P j,Ñ and P j,Ð . More precisely, let's consider a cell n P N and its two faces along e j , respectively F `" F j n j,´1 and F ´" F j n on the 'left' and 'right' (in the j-th direction). Thus, we define r ϕ j h,n as

@ x P K n , r ϕ j h,n pxq " ϕ j h,n pxq `δÑ h,F `pxq `δÐ h,F ´pxq, (45)
where the perturbations from the 'left' and 'right' interfaces (in the j-th direction) are respectively

δ Ñ h,F `pxq " ´γF `r ϕ j h ´ϕj,h ,F
`¯px zj q P j,Ñ pp x j q and δ Ð h,F ´pxq " ´γF ´r ϕ j h ´ϕj,h ,F

´¯px zj q P j,Ð pp x j q,

where x j " X j nj ´1 `hj p x j (see Figure 3).

Remark 2.14. One will observe that Definition 2.13 naturally ensures the increase to k j `1 of the flux polynomial degree with respect to x j , in addition to the continuity with respect to the j-th component (as the values at the interfaces of normal e j are then single-valued with the numerical trace polynomials):

@ n F P N j F , @ x P F j n F , r ϕ j h,n F pxq " r ϕ j h,pn F q j,1 pxq " γ F j n F r ϕ j h px zj q. ( 47
)
This allows to complete the introduction of the Flux Reconstruction method with the strong polynomial equation [START_REF] Di | The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications[END_REF], the flux numerical traces γ F r ϕ j h definition (with the examples of Definitions 2.9 and 2.11, or a general choice) and the flux approximation Definition 2.13: in any cell K n , with n P N , we get

iκMy h,n `p ÿ j"1 B r ϕ j h,n Bx j " pΠ h,k ψq n with r ϕ j h,n " ϕ j h,n `δÑ h,F j n j,´1 `δÐ h,F j n . ( 48
)
Remark 2.15. We highlight the fact we did not introduce any set of solution points, as it is used to being done in classic definitions of the FR method (see [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF][START_REF] Wang | A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids[END_REF]). Such nodes are usually used as collocation points, where the polynomial equation ( 48) is evaluated: the approximate solution is then reconstructed thanks to Lagrange polynomials. Yet, in this linear case, it has been shown thewe solution nodes have no influence on the method [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF][START_REF] Van Den Abeele | On the Stability and Accuracy of the Spectral Difference Method[END_REF] (without paying attention to the RHS projection which may depend on it, especially in case of interpolation), which avoids their introduction and additional notations.

(i) Flux correction in direction 1 (ii) Flux correction in direction 2
Remark 2.16. The Spectral Difference (SD) method can be seen as a specific case of the FR approach in this linear framework [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF]. Indeed, for each direction j, let's consider a set of k j `2 distinct flux points tp x F,j 0 , . . . , p x F,j kj `1u Ă r0, 1s such that p x F,j 0 " 0 and p x F,j kj `1 " 1, with the associated Lagrange polynomials ´LF,j l ¯lP 0,kj `1 . Thus, by considering the specific flux correction polynomials

P j,Ñ " L F,j 0 and P j,Ð " L F,j kj `1, (49)
the associated FR formulation allows to retrieve the SD one. Similarly to Remark 2.15, SD definitions (see [START_REF] Kopriva | A Conservative Staggered-Grid Chebyshev Multidomain Method for Compressible Flows[END_REF][START_REF] Liu | Spectral difference method for unstructured grids I: Basic formulation[END_REF]) rely on sets of solution points which have no influence in this linear case: they are only introduced as collocation points for the numerical implementation. Equivalently, in this case, r ϕ j h can be defined as the unique element of V h,k j,1 such that its restriction r ϕ j h,n to any K n is equal to:

' γ F j n j,´1 r ϕ j
h on F j n j,´1 (the 'left' interface in the j-th direction). ' ϕ j h on the internal faces of the cell associated to the flux points, defined by

´I1 n1 ˆ¨¨¨ˆI j´1 nj´1 ˆ!h j ´nj ´1 `p x j,F l ¯) ˆIj`1 nj`1 ˆ¨¨¨ˆI p np ¯lP 1,kj . (50)
' γ F j n r ϕ j h on F j n (the 'right' interface in the j-th direction).

Application of the FR approach to the problems of interest

We will then specify the applications of Definitions 2.9 and 2.11 to the Problems 1.1 and 1.4: this will allow to characterise the FR methods we introduce for them.

3.1. The FR formulation for the 1D wave equations Note 3.1. In the case of a single dimension, we will abandon the superscript of the direction j and the faces, reduced to a single point, will directly be referred as their number n F P 0, N (see (51)). Moreover, as the cells K n and intervals I n are identical, we will only refer to I n in this Subsection.

We take back Problem 1.1. Then, giving the set of numerical traces at the interfaces

´γn F r ϕ h ¯nF P 0,N , (51)
the corrected flux (see [START_REF] Wang | A review of flux reconstruction or correction procedure via reconstruction method for the Navier-Stokes equations[END_REF] and ( 46)) can be expressed on any I n , for n P 1, N , as

r ϕ h,n " ϕ h,n `´γ n´1 r ϕ h ´ϕh,n pX n´1 q ¯P Ñ n `´γ n r ϕ h ´ϕh,n pX n q ¯P Ð n , (52)
where P Ñ n and P Ð n stand for the transported versions of the flux correction polynomial functions on I n : @ x P I n , P Ñ n pxq " P Ñ pp xq and P Ð n pxq " P Ð pp xq where x " X n´1 `hp x.

Then, the last thing we need to completely define the FR method is the computation of these numerical traces.

To do so, as mentioned in Remark 2.10, let's diagonalize the flux operator F with respect to M " I 2 :

F " PΛP ´1 with P "

ˆ1 ´1 1 1 ˙and Λ " ˆ´1 0 0 1 ˙. (54)
Thus, we introduce the flux operator decomposition as the positive and negative parts of

F $ ' ' ' & ' ' ' % F Ñ " PΛ `P´1 " 1 2 ˆ1 ´1 ´1 1 ˙with Λ `" ˆ0 0 0 1 ˙, F Ð " PΛ ´P´1 " 1 2 ˆ´1 ´1 ´1 ´1˙w ith Λ ´" ˆ´1 0 0 0 ˙, (55)
which allows to define the flux numerical traces for internal faces according to Definition 2.9.

Concerning the boundaries, we refer to Subsection 1.3: we will highlight a surface metric Ă M whose associated decomposition of the flux operator F encompasses the BCs (2). Indeed, one has, for given complex values Z 0 ‰ 0 and Y 0 " pZ 0 q ´1 ‰ 0:

F " Ă MQΛQ ´1 with Ă M " ˆY0 0 0 Z 0 ˙and Q " ˆpZ 0 q 1{2 ´pZ 0 q 1{2 pY 0 q 1{2 pY 0 q 1{2 ˙, (56)
which allows to introduce

F Ñ Ă M " Ă MQΛ `Q´1 " 1 2 ˆY0 ´1 ´1 Z 0 ˙and F Ð Ă M " Ă MQΛ ´Q´1 " 1 2 ˆ´Y 0 ´1 ´1 ´Z0 ˙. (57)
Then, replacing pZ 0 , Y 0 q by pZ 1 , Y 1 q and pZ 2 , Y 2 q for the left and right boundaries respectively (and denoting the associated surface metrics as Ă M 1 and Ă M 2 ) that we suppose to be non-zero, one has

$ ' ' ' & ' ' ' % F Ñ Ă M1 yp0q " g 1 2 ˆY1 ´1˙ð ñ ˆY1 rup0q ´Z1 vp0q ´g1 s ´rup0q ´Z1 vp0q ´g1 s ˙" 0 ðñ up0q ´Z1 vp0q " g 1 , F Ð Ă M2 ypLq " g 2 2 ˆ´Y 2 ´1 ˙ðñ ˆ´Y 2 rupLq `Z2 vpLq ´g2 s ´rupLq `Z2 vpLq ´g2 s ˙" 0 ðñ upLq `Z2 vpLq " g 2 . ( 58
)
Thus, the framework of Definition 2.11, with

F Ñ Ă M1 , F Ð Ă M2
and (58), allows to take the BCs into account, and then defines an associated numerical trace.

Remark 3.2 (BCs for Dirichlet and Neumann conditions). Finally, one will note that the previous condition does not allow to take into account Dirichlet (Z 0 " 0) and Neumann (Y 0 " 0) conditions. To solve this limitation, let's first highlight that (55) leads, on the left and right boundaries respectively, to

$ ' ' ' & ' ' ' % F Ñ Ă M1 " 1 2 "ˆ1 `Y1 0 0 1 `Z1 ˙FÑ `ˆ1 ´Y1 0 0 1 ´Z1 ˙FÐ ȷ , F Ð Ă M2 " 1 2
"ˆ1 ´Y2 0 0 1 ´Z2 ˙FÑ `ˆ1 `Y2 0 0 1 `Z2 ˙FÐ ȷ , (59)
which allows to rewrite the BCs (58) as the expression of the incoming natural flux in terms of the reflection of the outgoing one:

$ ' ' & ' ' % F Ñ yp0q " R ãÑ F Ð yp0q `G1 with G 1 " g 1 Z 1 `1 ˆ1 ´1˙a nd R ãÑ " ˆ´R 1 0 0 R 1 ˙, F Ð ypLq " R Ðâ F Ñ ypLq `G2 with G 2 " g 2 Z 2 `1 ˆ´1 ´1˙a nd R ãÑ " ˆ´R 2 0 0 R 2 ˙, (60)
with

R 1 " Z 1 ´1 Z 1 `1 and R 2 " Z 2 ´1 Z 2
`1 the classic expressions of the reflection coefficients at normal incidence. Moreover, replacing the undefined term of Definition 2.9 with these rewritings allows to define a numerical trace in more general conditions.

The FR formulation for the 3D Maxwell problem

In what follows, we suppose that ε r and µ r are constant in the whole domain (which is a subcase of the more general hypotheses of Problem 1.4). In the framework of Subsection 2.3, we finally need to define the numerical traces. To do so, we consider the same mechanism as the one introduced in Subsection 3.1: the idea is to diagonalize each flux operator with respect to the scalar product induced by M (defined in [START_REF] Antonio | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]), before using their positive and negative parts as a decomposition. For example, we have F x " MP x Λ x pP x q ´1, with

P x " ¨0 0 0 0 0 pε r q ´1{2 ´pε r q ´1{2 0 pε r q ´1{2 0 0 0 0 pε r q ´1{2 0 ´pε r q ´1{2 0 0 0 0 0 0 pµ r q ´1{2 0 0 pµ r q ´1{2 0 pµ r q ´1{2 0 0 pµ r q ´1{2 0 pµ r q ´1{2 0 0 0 ‹ ‹ ‹ ‹ ‹ ‹ ' and Λ x " ¨´cI 2 0 2 0 2 0 2 cI 2 0 2 0 2 0 2 0 2 ' (61)
where Λ x is the diagonal matrix composed of the eigenvalues, while P x is made of the eigenvector coordinates.

The speed is also denoted as c " 1 ? ε r µ r .

Thus, we introduce the flux operator decomposition F x " F x,Ñ `Fx,Ð as its positive and negative parts.

This leads to F x,Ñ " MP x Λ x,`p P x q ´1 and F x,Ð " MP x Λ x,´p P x q ´1 with

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % F x,Ñ " MP x Λ x,`p P x q ´1 " 1 2 ˆY Γ t re x s ´Γˆr e x s Γ ˆre x s ZΓ t re x s ˙with Λ x,`" ¨02 0 2 0 2 0 2 cI 2 0 2 0 2 0 2 0 2 ', F x,Ð " MP x Λ x,´p P x q ´1 " 1 2 ˆ´Y Γ t re x s ´Γˆr e x s Γ ˆre x s ´ZΓ t re x s ˙with Λ x,´" ¨´cI 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 ', (62)
where Z " c µ r ε r and Y " Z ´1 stand for the internal impedance and admittance, and Γ t re x s and Γ ˆre x s denote the matrices respectively associated to the operators γ t re x s and γ ˆre x s (see [START_REF] Ernst | Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods[END_REF]). This allows to define the flux numerical traces for internal faces according to Definition 2.9, in a similar way to classic ones [START_REF] Hesthaven | High-order nodal discontinuous Galerkin methods for the Maxwell eigenvalue problem[END_REF][START_REF] Alireza | Computation of electromagnetic scattering and radiation using a time-domain finite-volume discretization procedure[END_REF].

Then, for the boundaries, we introduce a surface metric Ă M whose associated decomposition will allow to take into account the BCs, as detailed in Subsection 1.3. Thus, for given values Z 0 ą 0, Y 0 " pZ 0 q ´1 ą 0 and the surface metric

Ă M " ˆZY 0 ε r I 3 0 3 0 3 Z 0 Y µ r I 3 ˙(63)
the generalised eigenvalue problem leads to the decomposition

F x " F x,Ñ Ă M `Fx,Ð Ă M with F x,Ñ Ă M " 1 2
ˆY0 Γ t re x s ´Γˆr e x s Γ ˆre x s Z 0 Γ t re x s ˙and F x,Ð Ă M " 1 2
ˆ´Y 0 Γ t re x s ´Γˆr e x s Γ ˆre x s ´Z0 Γ t re x s ˙. (64)
Then, considering left and right boundary faces of respective outgoing normals n BΩ " ´ex and n BΩ " e x , and their impedance Z BΩ and admittance Y BΩ " pZ BΩ q ´1 that we suppose strictly positive (with the associated surface metric Ă M BΩ ), one has

$ ' ' ' & ' ' ' % F Ñ Ă M BΩ E " 1 2
ˆYBΩ Γ t r´e x s ´Γˆr ´ex s ˙g ðñ γ t r´e x se `ZBΩ γ ˆr´e x sh " g,

F Ð Ă M BΩ E " 1 2
ˆ´Y BΩ Γ t re x s Γ ˆre x s ˙g ðñ γ t re x se `ZBΩ γ ˆre x sh " g, (65)
where we use the fact that g P L 2 t pBΩq implies that g " γ t re x sg " γ t r´e x sg on these boundary faces. Thus, the framework of Definition 2.11 allows to take the BCs into account, and then defines an associated numerical trace. Similar considerations can be applied for the y and z axes. Remark 3.3 (BCs for perfect conductor conditions). Finally, one will note that the previous condition does not allow to take into account perfect electric (Z BΩ " 0) and magnetic (Y BΩ " 0) conductor conditions. To solve this limitation, let's first highlight that (62) leads, on left and right boundaries respectively, to

$ ' ' ' & ' ' ' % F x,Ñ Ă M BΩ " 1 2 "ˆ1 `YBΩ Z 0 0 1 `Y Z BΩ ˙Fx,Ñ `ˆ1 ´YBΩ Z 0 0 1 ´Y Z BΩ q ˙Fx,Ð ȷ , F x,Ð Ă M BΩ " 1 2
"ˆ1 ´YBΩ Z 0 0 1 ´Y Z BΩ ˙Fx,Ñ `ˆ1 `YBΩ Z 0 0 1 `Y Z BΩ ˙Fx,Ð ȷ , (66)
which allows to rewrite the BCs (65) as the expression of the incoming natural flux (associated to M) in terms of the reflection of the outgoing one:

$ ' ' & ' ' % F x,Ñ E " R ãÑ F x,Ð E `Gã Ñ with G ãÑ " Z Z BΩ `Z ˆY Γ t r´e x s ´Γˆr ´ex s ˙g and R ãÑ " ˆ´RI 3 0 3 0 3 RI 3 ˙, F x,Ð E " R Ðâ F x,Ñ E `GÐâ with G Ðâ " Z Z BΩ `Z ˆ´Y Γ t re x s Γ ˆre x s ˙g and R Ðâ " ˆ´RI 3 0 3 0 3 RI 3 ˙, ( 67
)
with the classic expression of the reflection coefficient at normal incidence R " Z BΩ ´Z Z BΩ `Z . Moreover, replacing the undefined term of Definition 2.9 with these rewritings allows to define a numerical trace in more general conditions.

Remark 3.4. These numerical traces can also be obtained by considering a 1D Riemann solver or an upwinding mechanism (see respectively 2.3.2 and 2.3.3 from [START_REF] Sirdey | Méthode itérative de Trefftz pour la simulation d'ondes électromagnétiques en trois dimensions[END_REF]). For the latter, the idea to look for the numerical trace as a linear combination of output traces from the neighbouring cells of the face (associating a virtual cell out of Ω in the boundary case).

A priori error estimates and well-posedness character

After detailing the construction of the Flux Reconstruction method, we focus on the homogeneous 1D timeharmonic wave Problem 4.4, with the numerical traces introduced in Subsection 3.1 and incoming BCs. Thus, this Section is devoted to show the well-posedness of the associated FR formulation, in addition to explicit a priori error estimates.

Main results

First, we detail the problem on which we focus, before summarising the main results.

Conjecture 4.1. For the rest of the paper, we assume the flux correction polynomials P Ñ and P Ð are of degrees exactly k `1 and that the parameters are such that k`1 ÿ l"0 pP Ñ q plq p0q p´iκhq ´l ‰ 0 and

k`1 ÿ l"0 pP Ð q plq p1q piκhq ´l ‰ 0, (68)
where pP Ñ q plq and pP Ð q plq denote the l th derivatives of P Ñ and P Ð respectively.

Remark 4.2. We highlight the fact that: 1. For any given flux correction polynomial functions P Ñ and P Ð , by introducing the real polynomials of degrees k

`1 T Ñ pXq " k`1 ÿ l"0 pP Ñ q plq p0q X l and T Ð pXq " k`1 ÿ l"0 pP Ð q plq p1q X l , (69)
Conjecture 4.1 may be refuted for at most Z k `1 2 ^values of κh (which is the maximum number of their couples of non-zero conjugated purely imaginary roots). 2. Condition (68) is verified for h small enough, as pP Ñ q pk`1q p0q ‰ 0 and pP Ð q pk`1q p1q ‰ 0. 3. All the numerical tests that we have implemented (for various flux correction polynomial functions, polynomial orders, domain lengths and wavenumbers) have verified this condition (see Subsection 5.1).

Claim 4.3.

In what follows, we restrict the study to the homogeneous case ψ " p0, 0q with the Boundary Conditions associated to Z 1 " 1 and Z 2 " 1, the numerical traces associated to (55) and (58), and flux correction polynomial functions of degrees exactly k `1 verifying Conjecture 4.1:

Problem 4.4. Find y h " pu h , v h q P V h,k and r ϕ h P W h,k such that for all n P 1,

N iκy h,n `d r ϕ h,n dx " 0 in I n , (70)
associated to the Boundary Conditions

# upX 0 q ´vpX 0 q " g 1 , upX N q `vpX N q " g 2 , (71)
and the flux numerical traces from Subsection 3.1.

Just before presenting the error estimates, we introduce two semi-norms for which they are developed.

Definition 4.5. For w P H 1 pT h q and w P " H 1 pT h q ‰ 2 (see ( 32)), we define the semi-norms |w| ˝and |w| ˝by

|w| 2 ˝" |w 1 pX 0 q| 2 `N´1 ÿ n"1 | w n | 2 `|w N pX N q| 2 and |w| 2 ˝" |w 1 pX 0 q| 2 `N´1 ÿ n"1 | w n | 2 `|w N pX N q| 2 , ( 72
)
where we have denoted the modulus and the usual euclidean norm on C 2 as | ¨| and the jumps of w and w as @ n P 1, N ´1 , w n " w n`1 pX n q ´wn pX n q and w n " w n`1 pX n q ´wn pX n q.

(73)

Moreover, the semi-norms |w| 1,h and |w| 1,h denote the broken H 1 semi-norms as

|w| 2 1,h " N ÿ n"1 ż In |pw n q p1q pxq| 2 dx and |w| 2 1,h " N ÿ n"1 ż In |pw n q p1q pxq| 2 dx. (74)
Finally, we denote the L 2 -norm on Ω as }w} 0 " }w} L 2 pΩq and }w} 0 " }w} L 2 pΩq .

Thus, we can introduce the main result, which ensures the well-posedness of the FR approach, in addition to explicit asymptotic estimates of the error. Moreover, for κh small enough, we have the following estimates of the error ε h " y ´yh (for the norms introduced in Definition 4.5):

piq |ε h | 2 ˝À pκLq pκhq 2k`1 `CÑ ˝|g 1 | 2 `CÐ ˝|g 2 | 2 ˘, piiq }ε h } 2 0 À κ ´1 pκLq pκhq 2k`2 `CÑ L 2 |g 1 | 2 `CÐ L 2 |g 2 | 2 ˘, piiiq |ε h | 2 1,h À κ pκLq pκhq 2k `CÑ H 1 |g 1 | 2 `CÐ H 1 |g 2 | 2 ˘, (75)
where the notation A À B means that there exists a constant C ą 0, independent of κ, h, k and L, such that we have A ď B C. Moreover, we have the explicit expressions of the constants

$ ' ' ' & ' ' ' % C Ñ ˝" 1 |T Ñ k`1 | 2 C Ð ˝" 1 |T Ð k`1 | 2 , $ ' ' ' ' & ' ' ' ' % C Ñ L 2 " |B Ñ k | 2 `1 3 pκLq 2 |A Ñ k | 2 |T Ñ k`1 | 2 C Ð L 2 " |B Ð k | 2 `1 3 pκLq 2 |A Ð k | 2 |T Ð k`1 | 2 and $ ' ' ' ' & ' ' ' ' % C Ñ H 1 " |C Ñ k | 2 |T Ñ k`1 | 2 C Ð H 1 " |C Ð k | 2 |T Ð k`1 | 2 , ( 76
)
which depend on $ ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' % A Ñ k " ż 1 0 P Ñ pp xq dp x and A Ð k " ż 1 0 P Ð pp xq dp x, B Ñ k " }P Ñ } L 2 p0,1q and
B Ð k " }P Ð } L 2 p0,1q , C Ñ k " › › ›pP Ñ q p1q › › › L 2 p0,1q
and

C Ð k " › › ›pP Ð q p1q › › › L 2 p0,1q , T Ñ k`1 " pP Ñ q pk`1q p0q and T Ð k`1 " pP Ð q pk`1q p1q . ( 77
)
Remark 4.7. One will note that Theorem 4.6 ensures quasi-optimal convergence orders, such as DG methods [START_REF] Antonio | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF][START_REF] Montseny | Dissipative terms and local time-stepping improvements in a spatial high order discontinuous galerkin scheme for the time-domain maxwell's equations[END_REF] (as it will be highlighted in the numerical Section 5). In particular, this behaviour is allowed by the regularity of the exact solution y (see Remark 4.10). Moreover, the dependencies are explicit in terms of the flux correction polynomial functions, the number of wavelengths in the domain κL, and the size of a cell relatively to a wavelength κh.

Finally, the dependence of the estimate piiq on κL highlights two asymptotic behaviours in ? κL and pκLq 3{2 , according to the value of A k , and then the flux correction polynomial functions (see Subsection 5.5).

Proof of the well-posedness result and error estimates

In what follows, we will detail the proof of Theorem 4.6. It will mainly rely on the introduction of oneway variables, for which the associated approximations and errors verify local strong equations. Then, error estimates will be derived going from one cell to its neighbour, and propagating the error along it in an iterative process: summing up on all the cells will thus allow to conclude.

Strong equations for the one-way variables

The interest of the Flux Reconstruction approach is that strong equations can be derived on each interval for the approximation of the solution, and then for the local error. First, in every cell I n , we define local flux correction polynomial functions P Ñ n and P Ð n thanks to the reference flux correction polynomial functions P Ñ and P Ð (see Definition 2.13) as @ x P I n , P Ñ n pxq " P Ñ pp xq and P Ð n pxq " P Ð pp xq with x " X n´1 `hp x.

Then, to drastically simplify the continuous and discrete systems, we decompose the solution in the one-way variables y Ñ " u ´v and y Ð " u `v,

with their discrete counterparts y Ñ h " u h ´vh and y Ð h " u h `vh .

Remark 4.8. In agreement with the BCs (71), we introduce the virtual values out of the boundaries y Ñ h,0 pX 0 q " g 1 and y Ð h,N `1pX N q " g 2 .

(81)

Then, these changes of variables allow to decouple the systems of Problems 1.1 and 4.4.

Proposition 4.9 (One-way reformulation). y Ñ and y Ð are solutions of @ x P Ω, " py Ñ q p1q pxq " ´iκy Ñ pxq y Ñ pX 0 q " g 1 and " py Ð q p1q pxq " iκy Ð pxq,

y Ð pX N q " g 2 , (82)
while y Ñ h and y Ð h are solutions of

@ n P 1, N , @ x P I n , $ ' & ' % ´yÑ h,n ¯p1q pxq " ´iκy Ñ h,n pxq ` y Ñ h n´1 pP Ñ n q p1q pxq, ´yÐ h,n ¯p1q pxq " iκy Ð h,n pxq ´ y Ð h n pP Ð n q p1q pxq, (83)
where we recall the jump of w P H 1 pT h q at any interface is denoted as @ n P 0, N , w n " w n`1 pX n q ´wn pX n q.

(84)

Proof. Equation ( 82) is derived from (1) with Z 1 " Z 2 " 1 and the definition (79) of y Ñ and y Ð . Equation ( 83) is due to the strong equation (70) verified by y h , the definition (80) of y Ñ h and y Ð h , and the following expressions of the interface corrections (which rely on the flux numerical trace operators defined thanks to (55) and ( 58)):

' For an internal interface n F P 1, N ´1 , one obtains

γ n F r ϕ h ´ϕh,n F `1pX n F q " 1 2 ˜´ u h ´vh n F u h ´vh n F ¸and γ n F r ϕ h ´ϕh,n F pX n F q " ´1 2 ˜ u h `vh n F u h `vh n F ¸. (85)
' For the boundary interfaces, one obtains by referring to Remark 4.8

γ 0 r ϕ h ´ϕh,1 pX 0 q " 1 2 ˜´ u h ´vh 0 u h ´vh 0 ¸and γ N r ϕ h ´ϕh,N pX N q " ´1 2 ˜ u h `vh N u h `vh N ¸, (86)
which can be seen as the direct application of (85) with fictive cells associated to the BCs.

□

Remark 4.10. We highlight the fact that the exact solutions of (82) can be easily obtained as @ x P Ω, y Ñ pxq " g 1 e ´iκpx´X0q and y Ð pxq " g 2 e iκpx´X N q , (87) which ensures that y P rC 8 pΩqs 2 .

Proposition 4.11 (Strong error equations in one-way variables). The one-way errors

ε Ñ h " y Ñ ´yÑ h and ε Ð h " y Ð ´yÐ h , (88)
satisfy the strong local equations

@ n P 1, N , @ x P I n , $ ' & ' % ´εÑ h,n ¯p1q pxq " ´iκε Ñ h,n pxq ` ε Ñ h n´1 pP Ñ n q p1q pxq, ´εÐ h,n ¯p1q pxq " iκε Ð h,n pxq ´ ε Ð h n pP Ð n q p1q pxq, (89)
and the boundary conditions ε Ñ h,0 pX 0 q " 0 and ε Ð h,N `1pX N q " 0. (90)

Proof. Equation ( 89) is due to (82), (83) and the continuity of y Ñ and y Ð , while (90) relies on (81) and (87). □

Well-posedness of the one-way discrete problem and local error estimates

In this Subsection, we will focus on the propagating errors, to derive preliminary local estimates. The proofs will only be detailed for ε Ñ h , as similar steps can be followed for ε Ð h . We begin with the proof that y Ñ h and y Ð h are well defined on each interval if Conjecture 4.1 is verified, and we detail a bound on the error jump realised on ε Ñ h and ε Ð h at the interface between two cells. Proposition 4.12. In the framework of Conjecture 4.1, the one-way variables y Ñ h and y Ð h are well-defined in V h,k and the error jumps are bounded as

@ n P 0, N ´1 , | ε Ñ h n | ď α Ñ κh,k pκhq k`1 |T Ñ k`1 | `ˇε Ñ h,n pX n q ˇˇ`|g 1 | ˘, @ n P 1, N , | ε Ð h n | ď α Ð κh,k pκhq k`1 |T Ð k`1 | `ˇε Ð h,n`1 pX n q ˇˇ`|g 2 | ˘, (91)
where we denote the derivatives at 0 and 1 of the reference flux correction polynomial functions as

@ l P 0, k `1 , T Ñ l " pP Ñ q plq p0q and T Ð l " pP Ð q plq p1q , (92) and α
Ñ κh,k " pκhq ´pk`1q |T Ñ k`1 | ˇˇˇˇk `1 ÿ l"0 p´iκhq ´l T Ñ l ˇˇˇˇa nd α Ð κh,k " pκhq ´pk`1q |T Ð k`1 | ˇˇˇˇk `1 ÿ l"0 piκhq ´l T Ð l ˇˇˇˇ, (93)
with the asymptotic behaviour α Ñ κh,k ÝÑ (94)

Proof. First, let's consider an interval n P 1, N for which we suppose y Ñ h,n´1 pX n´1 q is well-defined (which is the case for n " 1, as y Ñ h,0 pX 0 q " g 1 ). By differentiating (83) l ´1 times, we obtain for l P 1, k `1 `yÑ h,n ˘plq pX n´1 q " ´iκ `yÑ h,n ˘pl´1q pX n´1 q ``y Ñ h,n pX n´1 q ´yÑ h,n´1 pX n´1 q ˘pP Ñ n q plq pX n´1 q,

and by induction `yÑ h,n ˘plq pX n´1 q " p´iκq l y Ñ h,n pX n´1 q ``y Ñ h,n pX n´1 q ´yÑ h,n´1 pX n´1 q ˘l ÿ m"1 p´iκq l´m pP Ñ n q pmq pX n´1 q. (96)

Thus, the last thing we need to characterise the polynomial function y Ñ h,n (of degree at most k) is the expression of y Ñ h,n pX n´1 q, as the derivatives ˆ´y Ñ h,n ¯plq pX n´1 q ˙lP 1,k can be deduced from (96). To do so, by taking l " k `1 in this equation and as ´yÑ h,n ¯pk`1q pX n´1 q " 0, we get

p´iκq k`1 y Ñ h,n pX n´1 q ``y Ñ h,n pX n´1 q ´yÑ h,n´1 pX n´1 q ˘k`1 ÿ m"1 p´iκq k`1´m pP Ñ n q pmq pX n´1 q " 0. ( 97
)
The definition (78) of transported flux correction polynomial functions leads to @ m P 0, k `1 , @ x P I n , pP Ñ n q pmq pxq " h ´m pP Ñ q pmq pp xq with x " X n´1 `hp x,

and by using the fact that T Ñ 0 " P Ñ p0q " 1, (97) can be rewritten as

« k`1 ÿ l"0 p´iκhq ´l T Ñ l ff y Ñ h,n pX n´1 q " « k`1 ÿ l"1 p´iκhq ´l T Ñ l ff y Ñ h,n´1 pX n´1 q, (99)
where we recall that T Ñ l " pP Ñ q plq p0q, for l P 0, k `1 . Conjecture 4.1 ensures that

k`1 ÿ l"0 p´iκhq ´l T Ñ l ‰ 0:
y Ñ h,n pX n´1 q, and then all the derivatives ˆ´y Ñ h,n ¯plq pX n´1 q ˙lP 0,k , are well-defined with respect to y Ñ h,n´1 pX n´1 q. The local solution approximation y Ñ h,n being a polynomial function of degree at most k, it is then perfectly determined. Global well-posedness is deduced by induction. Concerning the jump estimates, we can write

´ y Ñ h n´1 " y Ñ h,n´1 pX n´1 q k`1 ÿ l"0 p´iκhq ´l T Ñ l . ( 100
)
Artificially introducing y Ñ pX n´1 q " g 1 e ´iκXn´1 (see (87)) to make ε Ñ h,n´1 pX n´1 q " py Ñ ´yÑ h,n´1 qpX n´1 q appear, we get

ε Ñ h n´1 " ´ y Ñ h n´1 " y Ñ pX n´1 q ´εÑ h,n´1 pX n´1 q k`1 ÿ l"0 p´iκhq ´l T Ñ l , (101)
where we have used the continuity of y Ñ at x " X n´1 . Then, the triangular inequality and the definition (93) of α Ñ κh,k lead to

| ε Ñ h n´1 | ď α Ñ κh,k pκhq k`1 |T Ñ k`1 | `|g 1 | `|ε Ñ h,n´1 pX n´1 q| ˘. (102)
Finally, the asymptotic behaviour of α Ñ κh,k results from T Ñ k`1 ‰ 0, as P Ñ is of degree k `1 (see Claim 4.3). □ Proposition 4.12 is the key argument to propagate the error estimate between two consecutive cells. The next Proposition estimates |ε Ñ h,n pX n q| and |ε Ð h,n`1 pX n q|, respectively at the right and left sides of each cell (see Figure 4). Proposition 4.13. In the framework of Conjecture 4.1, the propagation of the error on each interval leads to

@ n P 0, N , |ε Ñ h,n pX n q| ď "`1 `γÑ κh,k ˘n ´1‰ |g 1 | and |ε Ð h,n`1 pX n q| ď " `1 `γÐ κh,k ˘N´n ´1ı |g 2 |, ( 103
)
where we denote, thanks to the definition (93

) of α Ñ κh,k and α Ð κh,k γ Ñ κh,k " α Ñ κh,k |J Ñ κh | |T Ñ k`1 | pκhq k`1 and γ Ð κh,k " α Ð κh,k |J Ð κh | |T Ð k`1 | pκhq k`1 , (104)
J Ñ κh " 1 `ż 1 0 pP Ñ q p1q pp xq e iκhp x dp x and J Ð κh " 1 ´ż 1 0 pP Ð q p1q p1 ´p xq e iκhp x dp x. ( 105
)
with the asymptotic behaviour

γ Ñ κh,k " κhÑ0 `O `pκhq k`2 ˘and γ Ð κh,k " κhÑ0 `O `pκhq k`2 ˘. ( 106
)
Proof. By considering a variation of parameters on the differential equation (89) solved by ε Ñ h , we get

@ n P 1, N , ε Ñ h,n pX n q " " ε Ñ h n´1 J Ñ κh `εÑ h,n´1 pX n´1 q ‰ e ´iκh , (107)
where we have denoted, and simplified thanks to the change of variable x " hp x `Xn´1 and (98)

J Ñ κh " 1 `żIn pP Ñ n q p1q pxq e iκpx´Xn´1q dx " 1 `ż 1 0 pP Ñ q p1q pp xq e iκhp x dp x. ( 108
)
Thus, triangular inequality and Proposition 4.12 lead to

|ε Ñ h,n pX n q| ď p1 `γÑ κh,k q|ε Ñ h,n´1 pX n´1 q| `γÑ κh,k |g 1 |, ( 109
)
where we have used the definition (104) of γ Ñ κh,k . Thus, by considerations on an arithmetico-geometric sequence, we obtain

@ n P 0, N ´1 , `|ε Ñ h,n pX n q| `|g 1 | ˘ď p1 `γÑ κh,k q n `|ε Ñ h,0 pX 0 q| `|g 1 | ˘, (110)
and ( 90) allows to conclude. Finally, the asymptotic behaviour of γ Ñ κh,k comes out of Remark 4.14 and Lemma 4.15. □ Remark 4.14. We highlight two specific behaviours of γ Ñ κh,k and γ Ð κh,k with respect to

A Ñ k " ż 1 0 P Ñ pp xq dp x and A Ð k " ż 1 0 P Ð pp xq dp x. (111)
Indeed, by using Lemma 4.15, one has:

' If A Ñ k ‰ 0 and A Ð k ‰ 0, then γ Ñ κh,k and γ Ð κh,k can be rewritten as γ Ñ κh,k " α Ñ κh,k β Ñ κh,k |A Ñ k | |T Ñ k`1 | pκhq k`2 and γ Ð κh,k " α Ð κh,k β Ð κh,k |A Ð k | |T Ð k`1 | pκhq k`2 , ( 112
)
where we denote

β Ñ κh,k " |J Ñ κh | |A Ñ k |κh and β Ð κh,k " |J Ð κh | |A Ð k |κh , (113)
with the asymptotic behaviour

β Ñ κh,k ÝÑ κhÑ0 `1 and β Ð κh,k ÝÑ κhÑ0 `1. (114)
' Otherwise, by using the asymptotic behaviour (94) of α Ñ κh,k and α Ð κh,k , one has

γ Ñ κh,k " κhÑ0 `O `pκhq k`3 ˘and γ Ð κh,k " κhÑ0 `O `pκhq k`3 ˘. ( 115
)
Lemma 4.15. We have the asymptotic estimates

J Ñ κh " κhÑ0 `´iA Ñ k κh `O `pκhq 2 ˘and J Ð κh " κhÑ0 `´iA Ð k κh `O `pκhq 2 ˘. ( 116
)
Proof. Let's first note that ż 1 0 pP Ñ q p1q pp xq dp x " P Ñ p1q ´P Ñ p0q " ´1,

which allows to rewrite J Ñ κh as

J Ñ κh " ż 1 0 pP Ñ q p1q pp xq ´eiκhp x ´1¯d p x. ( 118
)
By writing e iκhp x ´1 " κhÑ0 `iκhp x `O `pκhq 2 ˘for p x P r0, 1s, an integration by parts leads to

J Ñ κh " κhÑ0 `iκh " ´ż 1 0 P Ñ pp xq dp x `rp x P Ñ pp xqs 1 0 ȷ `ż 1 0 pP Ñ q p1q pp xq O `pκhq 2 ˘dp x. ( 119
)
According to the fact that P Ñ p1q " 0 and the definition (111) of A Ñ k , we get iκh

" ´ż 1 0 P Ñ pp xq dp x `rp x P Ñ pp xqs 1 0 ȷ " ´iκhA Ñ k . ( 120
)
Finally, a Cauchy-Schwarz inequality ensures that, for κh small enough

ż 1 0 pP Ñ q p1q pp xq O `pκhq 2 ˘dp x À pκhq 2 › › ›pP Ñ q p1q › › › L 2 p0,1q . ( 121
) □ 0 0.2 0.4 0.6 0.8 1 ´2 ´1 0 1 2 ¨10 ´2 x ℜpε Ñ h q ℑpε Ñ h q |ε Ñ h | (i) ε Ñ h 0 0.2 0.4 0.6 0.8 1 ´2 ´1 0 1 2 ¨10 ´2 x ℜpε Ð h q ℑpε Ð h q |ε Ð h | (ii) ε Ð h Figure 4. Visualisation of ε Ñ h and ε Ð h along Ω for L " 1, N " 10
, the parameters of Table 1 and configuration SD CLo (see Section 5): error 'absorption' in each cell and global accumulation with propagation towards right and left respectively.

Global error estimates in the one-way framework

Note 4.16. For the rest of this Section, we recall that the notation A À B means that there exists a constant C ą 0, independent of κ, h, k and L, such that we have A ď B C. Moreover, the following estimates will be derived in an asymptotic regime, for κh small enough, which ensures that Conjecture 4.1 is verified (see Remark 4.2): all the local estimates of Subsection 4.2.2 are then valid. This finally leads to a bound on the jump at each interface, and then on the weak semi-norm | ¨|˝.

Proposition 4.17. The jump error at any internal interface is bounded, for κh small enough, by

@ n P 0, N ´1 , | ε Ñ h n | À pκhq k`1 |T Ñ k`1 | |g 1 | and | ε Ð h n`1 | À pκhq k`1 |T Ð k`1 | |g 2 |, (122)
which ensures the following bounds of the semi-norm of propagating errors

|ε Ñ h | 2 ˝À pκLq pκhq 2k`1 |T Ñ k`1 | 2 |g 1 | 2 and |ε Ð h | 2 ˝À pκLq pκhq 2k`1 |T Ð k`1 | 2 |g 2 | 2 . ( 123
)
Proof. First, by merging Propositions 4.12 and 4.13, we directly get

@ n P 1, N ´1 , | ε Ñ h n | ď α Ñ κh,k pκhq k`1 |T Ñ k`1 | p1 `γÑ κh,k q n |g 1 |. ( 124
)
Yet, as N " Lh ´1 and by using the asymptotic behaviour (106) of γ Ñ κh,k , one has for any n P 1, N ´1

`1 `γÑ κh,k ˘n " κhÑ0 `1 `O `pκLqpκhq k`1 ˘, (125)
which justifies the bound of (122) for κh small enough, with the asymptotic behaviour (94) of α Ñ κh,k . Moreover, recalling that N " Lh ´1, we have

N ´1 ÿ n"0 `1 `γÑ κh,k ˘2n " κhÑ0 `κL κh `O `pκLq 2 pκhq k ˘, (126)
which implies, for κh small enough, that

N ´1 ÿ n"0 | ε Ñ h n | 2 " κhÑ0 `" κL κh `O `pκLq 2 pκhq k ˘ȷ " α Ñ κh,k pκhq k`1 |T Ñ k`1 | |g 1 | ȷ 2 . ( 127
)
For the boundary terms, we have thanks to (90)

|ε Ñ h,1 pX 0 q| 2 " | ε Ñ h 0 | 2 , ( 128
)
and according to Proposition 4.13

|ε Ñ h,N pX N q| 2 " κhÑ0 `"O `pκLqpκhq k`1 ˘|g 1 | ‰ 2 . ( 129
)
Finally, by recalling the definition (72) of |ε Ñ h | ˝, the right boundary term can be neglected for κh small enough: this gives the expected semi-norm bounding as α Ñ κh,k has limit 1 as κh tends to 0 `. □

Then, let's consider the L 2 norm of the error.

Proposition 4.18. We have the following bound of the L 2 norm of the error for κh small enough

}ε Ñ h } 2 0 À L |B Ñ k | 2 `1 3 pκLq 2 |A Ñ k | 2 |T Ñ k`1 | 2 pκhq 2k`2 |g 1 | 2 and }ε Ð h } 2 0 À L |B Ð k | 2 `1 3 pκLq 2 |A Ð k | 2 |T Ð k`1 | 2 pκhq 2k`2 |g 2 | 2 , ( 130
)
where we denote

B Ñ k " }P Ñ } L 2 p0,1q and B Ð k " }P Ð } L 2 p0,1q . (131)
Proof. Let's take back the expression of ε Ñ h obtained by variation of parameters: for any n P 1, N , we have

@ x P I n , ε Ñ h,n pxq " « ε Ñ h n´1 ˜1 `ż x Xn´1 pP Ñ n q p1q ptq e iκpt´Xn´1q dt ¸`ε Ñ h,n´1 pX n´1 q ff e ´iκpx´Xn´1q . ( 132
)
Thus, by triangular and Young's inequalities, one has

@ n P 1, N , }ε Ñ h,n } 2 L 2 pInq À K Ñ n phq | ε Ñ h n´1 | 2 `h |ε Ñ h,n´1 pX n´1 q| 2 , ( 133
)
where we have denoted

K Ñ phq " ż In ˇˇˇˇ1 `ż x Xn´1 pP Ñ n q p1q ptq e iκpt´Xn´1q dt ˇˇˇˇ2 dx. ( 134
)
Then, by using Propositions 4.13 and 4.17 and Lemma 4.19, we obtain for κh small enough

@ n P 1, N , }ε Ñ h } 2 L 2 pInq À h " |B Ñ k | 2 |T Ñ k`1 | 2 pκhq 2k`2 ``p1 `γÑ κh,k q n´1 ´1˘2 ȷ |g 1 | 2 . (135)
Thus, by summing over n from 1 to N " Lh ´1, we obtain

}ε Ñ h,n } 2 0 À " L |B Ñ k | 2 |T Ñ k`1 | 2 pκhq 2k`2 `hS Ñ κh,k ȷ |g 1 | 2 , (136)
where we have denoted

S Ñ κh,k " N ÿ n"1 " p1 `γÑ κh,k q n´1 ´1‰ 2 . (137)
Finally, reporting Lemma 4.20 in this bounding gives the expected estimate for κh small enough:

' If A Ñ k ‰ 0
, we obtain the expected estimate. ' Otherwise, hS Ñ κh,k can be absorbed by the other term, and the estimate is still valid (with A Ñ k " 0). □ Lemma 4.19. For any n P 1, N , we have the bound

K Ñ n phq " ż In ˇˇˇˇ1 `ż x Xn´1 pP Ñ n q p1q ptq e iκpt´Xn´1q dt ˇˇˇˇ2 dx À δ Ñ κh,k h|B Ñ k | 2 , (138)
where we denote

δ Ñ κh,k " 1 `|B Ñ k | ´2 ż 1 0 ˇˇˇˇż p z 0 pP Ñ q p1q pp rq ´eiκhp r ´1¯d p r ˇˇˇˇ2 dp z, (139)
with the asymptotic behaviour

δ Ñ κh,k ÝÑ κhÑ0 `1. (140)
Proof. First, let's observe that K Ñ n phq is independent of n, as the the changes of variable x " hp z `Xn´1 and t " hp r `Xn´1 and (98) ensure that

K Ñ n phq " h ż 1 0 ˇˇˇˇ1 `ż p z 0 pP Ñ q p1q pp rq e iκhp r dp r ˇˇˇˇ2 dp z (141)
that we denote K Ñ phq for the rest of the proof. Then, let's note that @ p z P r0, 1s,

ż p z 0 pP Ñ q p1q pp rq dp r " P Ñ pp zq ´1, (142)
which allows to rewrite K Ñ phq as

K Ñ phq " h ż 1 0 ˇˇˇˇP Ñ pp zq `ż p z 0 pP Ñ q p1q pp rq ´eiκhp r ´1¯d p r ˇˇˇˇ2 dp z. (143)
Thus, Young's inequality allows to separate the bounding of K Ñ phq in two terms which will be treated separately

K Ñ phq À h |B Ñ k | 2 `h ż 1 0 ˇˇˇˇż p z 0 pP Ñ q p1q pp rq ´eiκhp r ´1¯d p r ˇˇˇˇ2 dp z. (144)
Focusing on the second term denoted as M Ñ phq, the Cauchy-Schwarz inequality leads, for κh small enough, to

M Ñ phq " h ż 1 0 ˇˇˇˇż p z 0 pP Ñ q p1q pp rq ´eiκhp r ´1¯d r ˇˇˇˇ2 dp z À κ 2 h 3 › › ›pP Ñ q p1q › › › 2 L 2 p0,1q (145)
Finally, by introducing

δ Ñ κh,k " 1 `M Ñ phq h |B Ñ k |
2 , this allows to conclude. □ Lemma 4.20. Since N " Lh ´1, we have the following asymptotic estimate

S Ñ κh,k " N ÿ n"1 " p1 `γÑ κh,k q n´1 ´1‰ 2 " κhÑ0 `pκLq 3 3 |A Ñ k | 2 |T Ñ k`1 | 2 pκhq 2k`1 `o `pκhq 2k`1 ˘. (146)
Proof. We develop the expression of S Ñ κh,k

S Ñ κh,k " N ´1 ÿ n"0 " p1 `γÑ κh,k q 2n ´2p1 `γÑ κh,k q n `1‰ " p1 `γÑ h,k q 2L{h ´1 2γ Ñ h,k `pγ Ñ h,k q 2 ´2 p1 `γÑ h,k q L{h ´1 γ Ñ h,k `L h , (147)
and realise series expansions of the different terms using the fact that, for κh small enough, γ Ñ κh,k À pκhq k`2 . □ Finally, a similar estimate can be derived for the H 1 -broken semi-norm.

Proposition 4.21. We have the following bound of the H 1 -broken semi-norm of the error for κh small enough

|ε Ñ h | 2 1,h À κpκLqpκhq 2k |C Ñ k | 2 |T Ñ k`1 | 2 |g 1 | 2 and |ε Ð h | 2 1,h À κpκLqpκhq 2k |C Ð k | 2 |T Ð k`1 | 2 |g 2 | 2 , (148)
where we denote

C Ñ k " › › ›pP Ñ q p1q › › › L 2 p0,1q
and

C Ð k " › › ›pP Ð q p1q › › › L 2 p0,1q . (149)
Proof. Let's take back the differential equation verified by ε Ñ h on each cell

@ n P 1, N , @ x P I n , `εÑ h,n ˘p1q pxq " ´iκε Ñ h,n pxq ` ε Ñ h n´1 pP Ñ n q p1q pxq. (150)
Thus, we obtain by triangular and Young's inequalities

@ n P 1, N , › › › `εÑ h,n ˘p1q › › › 2 L 2 pInq À κ 2 }ε Ñ h,n } 2 L 2 pInq `| ε Ñ h n´1 | 2 › › ›pP Ñ n q p1q › › › 2 L 2 pInq . (151)
Yet, (98) ensures that

› › ›pP Ñ n q p1q › › › L 2 pInq " 1 h 1{2 › › ›pP Ñ q p1q › › › L 2 p0,1q , (152)
and summing over n from 1 to N " h ´1L with the jump estimate of Proposition 4.17 leads to

|ε Ñ h | 2 1,h À κ 2 }ε Ñ h } 2 0 `κ pκLqpκhq 2k |T Ñ k`1 | 2 › › ›pP Ñ q p1q › › › 2 L 2 p0,1q |g 1 | 2 . (153)
Finally, the L 2 -broken norm estimate of Proposition 4.18 allows to neglect the first term when κh is small enough, and then to conclude. □

FR formulation well-posedness and global error estimates

Finally, we gather the precedent results to prove Theorem 4.6.

Proof. Proposition 4.12 ensures that y Ñ h and y Ð h are well-defined on each interval, in the framework of Conjecture 4.1. Then, by inverting their definitions (80), we get

@ n P 1, N , u h,n " y Ñ h,n `yÐ h,n 2 and v h,n " y Ð h,n ´yÑ h,n 2 , (154)
which ensures that y h is also uniquely defined in V h,k . Concerning the error estimates, the proof is exactly the same in the three cases, and we will then only present it for (i), the semi-norm | ¨|˝.

By definition, we have

|ε h | 2 ˝" |u ´uh | 2 ˝`|v ´vh | 2 ˝, (155)
and by using the definition (88) of ε Ñ h and ε Ð h , we obtain thanks to the triangular and Young's inequalities

|ε h | 2 ˝À |ε Ñ h `εÐ h | 2 ˝`|ε Ð h ´εÑ h | 2 À |ε Ñ h | 2 ˝`|ε Ð h | 2 ˝. (156)
Thus, Proposition 4.17 allows to conclude. Similarly, Proposition 4.18 is used for (ii), while Proposition 4.21 corresponds to (iii). □

Numerical validations for the time-harmonic 1D wave equations

After the introduction of the error estimates of Theorem 4.6, this section will aim at illustrating them through numerical experiments.

To do so, we take back Problem 4.4 in a domain Ω " r0, Ls, with a uniform mesh of N P N ˚cells of size h " L N . In all the following experiments, we consider the parameters summarised in Table 1 with uniform polynomial degree k and flux correction polynomial functions in all the cells.

Concerning the flux correction polynomial functions, we will focus on 'symmetrical couples' of degree exactly k `1 (as proposed in [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF][START_REF] Wang | A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids[END_REF]), in that @ x P r0, 1s, P Ð pxq " P Ñ p1 ´xq,

κ g 1 g 2 2π 2.3 `0.4i ´1.2i (157)
Table 1. Common characteristics to the different simulation configurations. which avoids asymmetrical behaviours for waves coming from the left and the right. Indeed, exchanging the BCs will have no effect on the solution precision and no direction will be favored. Moreover, referring to the constants of (77), this implies

A k " |A Ñ k | " |A Ð k |, B k " |B Ñ k | " |B Ð k |, C k " |C Ñ k | " |C Ð k | and T k`1 " |T Ñ k`1 | " |T Ð k`1 |, (158)
where A k , B k , C k and T k`1 are referred as 'flux correction polynomial constants'. Similarly, the 'estimate constants' of (76) will be denoted as C ˝" C Ñ ˝" C Ð ˝and so on. In particular, the numerical experiments will be realised for four different families, for which P Ñ is defined as: ' The first Lagrange polynomial associated to the Chebyshev-Lobatto nodes

@ l P 0, k `1 , x CLo l " 1 2 ˆ1 ´cos ˆl k `1 π ˙˙and @ x P r0, 1s, P Ñ pxq " ź lP 1,k`1 x ´xCLo l x CLo 0 ´xCLo l . (159)
According to Remark 2.16, the associated FR formulation can then be seen as a SD one in this configuration: it will then be referred as SD CLo. ' The first Lagrange polynomial associated to the Internal Gauss nodes px IG l q lP 0,k`1 . px IG l q lP 1,k are defined as the roots of the Legendre polynomial of degree k (translated to r0, 1s), referred as P Legendre k , while x IG 0 " 0 and x IG k`1 " 1. This then leads to @ x P r0, 1s, P Ñ pxq "

ź lP 1,k`1 x ´xIG l x IG 0 ´xIG l . ( 160
)
Similarly to the previous case, the associated FR formulation will be referred as SD IG. ' The right Radau polynomial of degree k `1 defined as @ x P r0, 1s, P Ñ pxq "

P Radau k`1 pxq " p´1q k`1 2 ´P Legendre k`1 pxq ´P Legendre k pxq ¯. (161)
In particular, it is orthogonal to all the polynomial functions of degree at most k ´1, denoted as P k´1 pr0, 1sq. The associated FR formulation will then be referred as F R Radau.

In this linear configuration, this method can also be seen as a nodal DG formulation using the same flux numerical traces [START_REF] Allaneau | Connections between the filtered discontinuous Galerkin method and the flux reconstruction approach to high order discretizations[END_REF][START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF]. ' The G2 polynomial (according to Huynh's original denomination in [START_REF] Huynh | A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[END_REF]) of degree k `1 defined as @ x P r0, 1s, P Ñ pxq " P G2 k`1 pxq "

k 2k `1 P Radau k`1 pxq `k `1 2k `1 P Radau k pxq. (162)
In particular, it is orthogonal to all the polynomial functions of P k´2 pr0, 1sq and verifies `P G2 k`1 ˘p1q p1q " 0. The associated FR formulation will then be referred as F R G2.

Remark 5.1. One will note that in the case k " 0, the flux correction polynomial functions must be defined as @ x P r0, 1s, P Ñ pxq " 1 ´x and P Ð pxq " x,

because of [START_REF] Wang | A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids[END_REF]. Then, we will only focus on the cases k ě 1 for which the four families differ.

First, this Section will numerically verify Conjecture 4.1 for these flux correction polynomial functions of different degrees. Then, it will present the evolution of the associated estimate constants before illustrating error behaviours relatively to various parameters: a study of the h-convergence, an overview of the k-convergence and a focus on the κL-dependence. 2. Numerical approximation of the roots of polynomial T Ñ for the different flux correction polynomial families and polynomial degrees k P 1, 3 .

To do so, Table 2 presents the roots of the polynomial T Ñ (defined in (69)) for the different flux correction families and various polynomial degrees k: the key point is that none of them is purely imaginary. Indeed, according to Remark 4.2, the Conjecture is verified if and only if p´iκhq ´1 and piκhq ´1 are not roots of T Ñ and T Ð respectively. Moreover, as correction symmetry (157) implies that T Ð pXq " T Ñ p´Xq and none of them admits a purely imaginary root, the FR method is well-posed for any mesh in these configurations. We also insist on the fact that all the numerical experiments we intended have fallen within the framework of the Conjecture 4.1.

Flux correction polynomial constants

We present in Table 3 the evolution of the constants A k , B k , C k and T k`1 according to the polynomial degree k, for the four flux correction polynomial families we introduced. In particular, one will note the decreasing behaviours of A k and B k , while C k and T k`1 show an opposite behaviour. Thus, the constants C ˝and C L 2 will naturally decrease (and then improve their respective convergence constants) with k, while the behaviour of C H 1 will depend on the considered correction. Moreover, an important remark is that the value of A k for SD IG, F R Radau and F R G2 is 0 for k ě 2 (due to their orthogonality properties to polynomial spaces). Coming back to the definition (76) of C L 2 and Remark 4.7, one will then expect them to show a ? κL-dependence for their }ε h } 0 estimate (see Subsection 5.5).

SD CLo

SD IG FR Radau FR G2 k A k B k C k T k`1 A k B k C k T k`1 A k B k C k T k`1 A k B k C k T k`1 1

h-convergence: mesh refinement influence

In this Subsection, we focus on the h-convergence of the method. Then, we fix the value L " 1 and we make the number of cells N " Lh ´1 vary: the objective is to highlight the behaviour of the method for different polynomial orders k when the mesh is refined. respectively (with respect to the norm of the exact solution), when the mesh is refined.

In accordance with the estimates (75) of Theorem 4.6, quasi-optimal convergence orders are shown for the three error estimates when κh tends to 0. Moreover, we highlight the fact that the SD CLo and FR G2 method have shown inferior abilities to FR Radau in this harmonic framework, as it had already been observed for CFD problems [START_REF] Castonguay | Application of High-Order Energy Stable Flux Reconstruction Schemes to the Euler Equations[END_REF][START_REF] Vincent | Insights from von Neumann analysis of high-order flux reconstruction schemes[END_REF]. 4. Evolution of the relative error |ε h | rel ˝according to the number of mesh cells N and the polynomial degree k.

SD

k-convergence: polynomial degree influence

Then, we focus on the k-convergence of the method. Then, we fix the values L " 1 and N " 4, and we make the polynomial degree k vary: the objective is to highlight the behaviour of the method, depending on the flux correction polynomial functions, when k increases. Then, Figure 5 presents the relative values of the errors |ε h | rel ˝, }ε h } rel 0 and |ε h | rel 1,h respectively, when the polynomial degree evolves. One will note comparable behaviours, as the error decreases in similar ways for all the flux correction polynomial families, as it could be expected.

κL-dependence of the errors

Finally, we focus on the error dependence on κL: the objective is to highlight the behaviour of the method when the number of wavelengths in the domain becomes larger and larger. So as to get comparable experiments in this subsection, we fix the number of degrees of freedom (dof) per wavelength, N dof {λ " 2πpk`1q κh , to a constant value N dof {λ " 600. Then, κ " 2π is let constant, while various 5.89ˆ10 ´8 -5.0 4.62ˆ10 ´8 -5.0 2.82ˆ10 ´8 -5.0 5.88ˆ10 ´8 -5.0 100 3.03ˆ10 ´11 -5.01 2.38ˆ10 ´11 -5.0 1.46ˆ10 ´11 -5.0 3.03ˆ10 ´11 -5.0 Table 5. Evolution of the relative error }ε h } rel 0 according to the number of mesh cells N and the polynomial degree k. In accordance with the estimates (75) of Theorem 4.6, Remark 4.7 and the flux correction polynomial constants of Table 3, different behaviours can be highlighted:

SD CLo SD IG FR Radau FR G2 k N }ε h } rel 0 r N }ε h } rel 0 r N }ε h } rel 0 r N }ε h } rel 0 r N 1

SD CLo

SD IG FR Radau FR G2 k N |ε h | rel 1,h r N |ε h | rel 1,h r N |ε h | rel 1,h r N |ε h | rel 1,h r N 1

SD CLo

SD IG FR Radau FR G2 k L |ε h | ˝rL |ε h | ˝rL |ε h | ˝rL |ε h | ˝rL 1 0.1 1.12ˆ10 ´3 0.5 1.12ˆ10 ´3 0.5 7.43ˆ10 ´4 0.
' Tables 7 and9 show a

? κL-dependence of |ε h | ˝and |ε h | 1,h
on the number of wavelengths in the domain for all the flux correction polynomial families. ' For }ε h } 0 , two cases can be detailed.

First, if A k " 0 (for F R Radau, and SD IG and F R G2 when k ě 2), the same ? κL-dependence is observed. On the contrary, if A k ą 0 (for SD CLo), the ? κL-dependence only happens when κL tends to 0, while having κL diverging to `8 implies a pκLq 3{2 -dependence. Moreover, this transition seems to happen for larger values of L when the polynomial degree k increases. Finally, the estimates (75) highlight the potential interest of flux correction polynomial functions for which A k " 0, as their dependence on the number of wavelengths in the domain is optimised. 1.9ˆ10 ´6 0.5 1.29ˆ10 ´6 0.5 2.41ˆ10 ´6 0.5 1.0 6.42ˆ10 ´6 0.5 6.0ˆ10 ´6 0.5 4.08ˆ10 ´6 0.5 7.64ˆ10 ´6 0.5 10.0 2.04ˆ10 ´5 0.51 1.9ˆ10 ´5 0.5 1.29ˆ10 ´5 0.5 2.41ˆ10 ´5 0.5 4 0.1 1.61ˆ10 ´8 0.5 1.43ˆ10 ´8 0.5 9.17ˆ10 ´9 0.5 1.74ˆ10 ´8 0.5 1.0 5.1ˆ10 ´8 0.5 4.52ˆ10 ´8 0.5 2.9ˆ10 ´8 0.5 5.5ˆ10 ´8 0.5 10.0 1.62ˆ10 ´7 0.5 1.43ˆ10 ´7 0.5 9.17ˆ10 ´8 0.5 1.74ˆ10 ´7 0.5 Table 9. Evolution of the non-relative error |ε h | 1,h according to the domain size L, for a constant number of dof per wavelength.

SD

CLo SD IG FR Radau FR G2 k L }ε h } 0 r L }ε h } 0 r L }ε h } 0 r L }ε h } 0 r L 1
G2 k L |ε h | 1,h r L |ε h | 1,h r L |ε h | 1,h r L |ε h | 1,h r L 1 0.1 2.95ˆ10 ´2 0.5 2.95ˆ10 ´2 0.5 2.57ˆ10 ´2 0.

Time-harmonic 3D Maxwell problem: error analysis through numerical experiments

Finally, we come back to the Maxwell Problem 1.4. Contrary to the 1D wave problem, we did not derive a priori error estimates for this case. Yet, we present numerical results for the implementation of the FR formulation associated to the numerical traces introduced in Remark 3.3 of Subsection 3.2. We consider a domain Ω " r0, Ls 3 , with a uniform mesh of N cells of size h " L N in each direction. In all the following experiments, we consider the parameters summarised in Table 10 with uniform polynomial degree k and flux correction polynomial functions in all the cells and every direction (even if direction-wise choices could have been made). Finally, we take back the four flux correction polynomial families introduced in Section 5 to compare their influence in this 3D configuration. Moreover, we consider impedance boundary conditions associated to an exact solution E exa " pe exa , h exa q of κ L ε r µ r Z BΩ k pw e pw h pw P d X d 2π 1 1 1 1 2 1{2 πp1, ´1, 0q p0, 0, 1q ´2´1{2 p1, 1, 0q 2 ´1{2 p1, ´1, 0.q ´0.5p1, 1, 1q

where k pw ˆepw " ´κµ r h pw and k pw ˆhpw " κϵ r e pw . ' An electric dipole of dipole moment P d and position X d , defined in the framework of the dipole approximation [START_REF] Labat | Multi-scale modeling of the electromagnetic wave scattering by small obstacles[END_REF] by

@
where we denote r " }x ´Xd } and r x " x ´Xd r .

Thus, we focus on the h-convergence of the formulation. We fix all the parameters and we make the number of cells in each direction, N , vary: the goal is to highlight the effect of mesh refinement on the error ε h . Moreover, the convergence rate is evaluated in a similar way as in Subsection 5.3 (see (164)). Tables 11 and12 present the relative L 2 -norm error }ε h } rel 0 , for the two exact solutions, when the number of cells per direction increases. For all the flux correction polynomial functions and polynomial degrees k, it shows an asymptotic quasi-optimal behaviour, as it was previously proved in the 1D case, and highlights properties of a high-order method. Moreover, the flux correction polynomial family hierarchy seems to be the same as in the 1D wave equations case: FR Radau clearly presents the best abilities in these contexts. Remark 6.1. One will note that the finest mesh refinements are not presented for the highest polynomial degrees in Tables 11 and12: this is due to the fact these cases imply a prohibitive memory cost (due to the direct resolution of the linear system) which prevented us from computing the associated results.

Conclusion

We introduced a Flux Reconstruction method in the framework of generic time-harmonic linear hyperbolic problems, for Cartesian meshes of hyperrectangular domains. It relies on the piecewise polynomial approximation of the solution and fluxes, thanks to interface corrections to impose boundary conditions and flux continuity. Such perturbations are realised thanks to user-dependent flux correction polynomial functions, whose choice allows to recover some classical methods (as nodal Discontinuous Galerkin or Spectral Difference). A general flux numerical trace definition through flux operator decomposition is also detailed, and the associated FR methods have been specified for the homogeneous 1D time-harmonic wave equations and the 3D Maxwell problem. Then, we focused on the 1D wave equations with ingoing boundary conditions, for which well-posedness conditions of the associated FR method have been derived. This also allowed to determine quasi-optimal a priori error estimates in the asymptotic limit for various error norms, with explicit expressions of the constants in terms of the polynomial degree, flux correction polynomials, wavenumber, domain length and mesh properties. Numerical experiments were realised for some classical flux correction polynomial families, and confirmed the expected behaviours of the estimates in terms of the mesh refinement, polynomial degree and dependence on the number of wavelengths in the domain. Finally, the FR method has been implemented for the 3D time-harmonic Maxwell problem. Even if a priori estimates were not derived in this case, numerical experiments highlighted quasi-optimal orders for two different actual solutions and allow to expect this method to show interesting properties for wave equations.

The presented approach relied on the possibility to uncouple the 1D wave equations thanks to the one-way variables, and further developments may be needed to extend this work to general boundary conditions or heterogeneous materials: in case this strong framework would not be applicable anymore, one could contemplate classic methodology for the associated weak formulation. The same remark holds in the 3D framework of the Maxwell problem, for which it seems difficult to introduce uncoupling variables (which would rely on the simultaneously diagonalisable character of the flux operators). Furthermore, optimisations may be applied to such a method: one can think about Domain Decomposition Methods [START_REF] Collino | Exponentially convergent non overlapping domain decomposition methods for the Helmholtz equation[END_REF][START_REF] Ernst | Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods[END_REF][START_REF] Modave | An efficient domain decomposition method with cross-point treatment for Helmholtz problems[END_REF] or static condensation [START_REF] Carlos | Hybridized formulations of flux reconstruction schemes for advection-diffusion problems[END_REF] as for example Hybrid Discontinuous Galerkin (HDG) [START_REF] Cockburn | Unified Hybridization of Discontinuous Galerkin, Mixed, and Continuous Galerkin Methods for Second Order Elliptic Problems[END_REF][START_REF] Moro | A hybridized discontinuous petrov-galerkin scheme for scalar conservation laws[END_REF][START_REF] Nguyen | Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell's equations[END_REF] or Hybrid High Order (HHO) [START_REF] Montseny | Dissipative terms and local time-stepping improvements in a spatial high order discontinuous galerkin scheme for the time-domain maxwell's equations[END_REF][START_REF] Di | The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications[END_REF], allowing to reduce the amount of degrees of freedom. This work allowed to explicit the influence of the error norms with respect to the flux correction polynomial functions and discretisation parameters. Numerically, the Radau polynomials (corresponding to a nodal DG method) led to the best results in both these time-harmonic frameworks, in comparison with the SD and G2 corrections. Then, the explicit character of the 1D estimates allows to get a glimpse of the possibility to optimise their constants according to the flux correction polynomial, which will be investigated in further work in order to define an 'optimised' FR scheme.

Figure 3 .

 3 Figure 3. Illustration, for p " 2 and k " p2, 1q, of the flux correction with the numerical traces, with respect to direct flux approximation: r ϕ 1 h P V h,p3,1q and r ϕ 2 h P V h,p2,2q .

Theorem 4 . 6 .

 46 For κL fixed and κh verifying Conjecture 4.1, the Flux Reconstruction formulation associated to Problem 4.4 is well-posed.

Table 3 .

 3 18 3.57 6144 0.0 0.2 4.33 8400 0.0 0.22 5.0 15120 0.0 0.21 4.22 6720 Evolution of the flux correction polynomial constants for the different flux correction families and polynomial degrees k P 1, 4 .

  For a given error norm |ε h | (along the three ones we introduced in Definition 4.5) and a cell number N ě 2, the estimate of the convergence rate r N is given thanks to the error norm values for N ˘1 cells, respectively denoted |ε h | N ˘1, as r N " log 10 p|ε h | N `1q ´log 10 p|ε h | N ´1q log 10 pN `1q ´log 10 pN ´1q . (164) Then, Tables 4, 5 and 6 present the values of r N and the relative error norms |ε h | rel ˝, }ε h } rel 0 and |ε h | rel 1,h

  ´2 -3.44 7.95ˆ10 ´3 -3.43 4.54ˆ10 ´3 -3.43 1.06ˆ10 ´2 -3´4 -4.45 5.75ˆ10 ´4 -4.45 3.2ˆ10 ´4 -4.45 7.18ˆ10 ´4 -4.44 22 1.02ˆ10 ´6 -4.5 7.49ˆ10 ´7 -4.5 4.16ˆ10 ´7 -4.5 9.37ˆ10 ´7 -4.5 100 1.13ˆ10 ´9 -4.5 8.24ˆ10 ´10 -4.5 4.58ˆ10 ´10 -4.5 1.03ˆ10 ´9 -4.5 Table

4 5 1 .Table 6 .Figure 5 .

 4165 Figure 5. Evolution of the relative errors according to the polynomial degree k for the four flux correction families: SD CLo ( ), SD IG ( ), FR Radau ( ) and FR G2 ( ).

  ´3 0.51 10.0 1.14ˆ10 ´2 0.52 1.14ˆ10 ´2 0.52 7.43ˆ10 ´3 0.5 2.39ˆ10 ´2 0´10 0.5 3.01ˆ10 ´10 0.5 1.67ˆ10 ´10 0.5 3.77ˆ10 ´10 0.5 1.0 1.3ˆ10 ´9 0.5 9.53ˆ10 ´10 0.5 5.29ˆ10 ´10 0.5 1.19ˆ10 ´9 0.5 10.0 4.12ˆ10 ´9 0.5 3.01ˆ10 ´9 0.5 1.67ˆ10 ´9 0.5 3.77ˆ10 ´9 0.5

5 4 . 3 0.1 2 .

 432 47ˆ10 ´2 0.5 1.0 9.37ˆ10 ´2 0.5 9.37ˆ10 ´2 0.5 8.14ˆ10 ´2 0´4 0.5 2.43ˆ10 ´4 0.5 1.82ˆ10 ´4 0.5 3.32ˆ10 ´4 0.5 1.0 7.88ˆ10 ´4 0.5 7.7ˆ10 ´4 0.5 5.75ˆ10 ´4 0.5 1.05ˆ10 ´3 0.5 10.0 2.52ˆ10 ´3 0.51 2.43ˆ10 ´3 0.5 1.82ˆ10 ´3 0.5 3.32ˆ10 ´3 0.5 03ˆ10 ´6 0.5

  5.1. Numerical validation of Conjecture 4.1As highlighted in the proof Subsection 4.2, the well-posedness of the FR method relies on Conjecture 4.1. We will then numerically verify it for the four families of flux correction polynomial functions. ˘0.18i 0.24 ˘0.1i 0.09 ˘0.17i 0.22 ˘0.08i 0.1 ˘0.14i 0.19 ˘0.06i 0.1 ˘0.19i 0.23 ˘0.09i

	k	SD CLo		SD IG		FR Radau	FR G2
	1 0.37 ˘0.33i		0.38 ˘0.33i		0.33 ˘0.24i		0.5 ˘0.5i
	2 0.15 ˘0.25i	0.36	0.17 ˘0.24i	0.33	0.16 ˘0.18i	0.27	0.18 ˘0.27i	0.38
	3 0.08 Table

Table 7 .

 7 Evolution of the non-relative error |ε h | ˝according to the domain size L, for a constant number of dof per wavelength.

Table 8 .

 8 Evolution of the non-relative error }ε h } 0 according to the domain size L, for a constant number of dof per wavelength.

		0.1	2.38ˆ10 ´5 0.53 2.38ˆ10 ´5 0.53 1.57ˆ10 ´5 0.5	5.97ˆ10 ´5 0.57
		1.0	1.44ˆ10 ´4 1.23 1.44ˆ10 ´4 1.23 4.96ˆ10 ´5 0.5	5.25ˆ10 ´4 1.38
		10.0	3.9ˆ10 ´3	1.5	3.9ˆ10 ´3	1.5	1.59ˆ10 ´4 0.53 1.56ˆ10 ´2 1.5
		0.1	1.52ˆ10 ´7 0.51	1.4ˆ10 ´7	0.5	8.88ˆ10 ´8 0.5	2.22ˆ10 ´7 0.5
	2	1.0	6.02ˆ10 ´7 0.86 4.41ˆ10 ´7 0.5	2.81ˆ10 ´7 0.5	7.02ˆ10 ´7 0.5
		10.0 1.16ˆ10 ´5 1.48 1.41ˆ10 ´6 0.52 8.88ˆ10 ´7 0.5	2.26ˆ10 ´6 0.54
		0.1 9.94ˆ10 ´10 0.5 8.75ˆ10 ´10 0.5 5.42ˆ10 ´10 0.5	1.2ˆ10 ´9	0.5
	3	1.0	3.62ˆ10 ´9 0.75 2.77ˆ10 ´9 0.5	1.71ˆ10 ´9 0.5	3.79ˆ10 ´9 0.5
		10.0 5.78ˆ10 ´8 1.47 8.75ˆ10 ´9 0.5	5.42ˆ10 ´9 0.5	1.2ˆ10 ´8	0.5
		0.1 6.59ˆ10 ´12 0.5 5.63ˆ10 ´12 0.5 3.43ˆ10 ´12 0.5 7.15ˆ10 ´12 0.5
	4	1.0 2.26ˆ10 ´11 0.65 1.78ˆ10 ´11 0.5 1.09ˆ10 ´11 0.5 2.26ˆ10 ´11 0.5
		10.0 2.82ˆ10 ´10 1.5 5.63ˆ10 ´11 0.5 3.44ˆ10 ´11 0.51 7.15ˆ10 ´11 0.5
			SD CLo	SD IG		FR Radau	FR

Table 10 .

 10 Common characteristics and parameters of the different simulation configurations.the Maxwell problem, which leads to:g " γ t e exa `ZBΩ n BΩ ˆγt h exa on BΩ,(166)and the numerical approximation E h is thus compared to E exa thanks to the error ε h " E exa ´Eh . Finally, we consider two different exact solutions, whose parameters are summarised in Table10:' A Plane Wave defined by @ x P Ω, e exa pxq " e pw e ikpw¨x and h exa pxq " h pw e ikpw¨x ,

  x P Ω,

	$ ' ' ' &	e exa pxq "	e ´iκr 4πr	"ˆ´1 r 2	´iκ r	`κ2 ˙pr x ˆpP ˆr xqq	`2 ˆ1 r 2	`iκ r	˙pP ¨r xq r x	ȷ	,
	' ' ' %	h exa pxq "	e ´iκr 4πr	ˆ´iκ r	`κ2 ˙pr x ˆPq ,