Convergence et stabilité d'une méthode Galerkin discontinue basée sur un schéma saute-mouton d'ordre élevé pour la résolution des équations de Maxwell en maillage non-conformes

Résumé : Nous étudions dans ce rapport la stabilité et la convergence d'une méthode Galerkin discontinue d'ordre élevé pour la résolution des équations de Maxwell en domaine temporel en maillages non-conformes. La méthode repose sur des fonctions de base nodales pour approcher le champ électromagnétique dans un simplexe, un schéma centré pour calculer les flux sur les interfaces entre éléments voisins du maillage et un schéma saute-mouton d'ordre élevé pour l'intégration temporelle. De plus, l'ordre d'interpolation est défini au niveau d'un élément et le maillage est raffiné d'une manière non-conforme autorisant un nombre arbitraire de noeuds flottants. On prouve que la méthode est stable sous une condition de type CFL et qu'un équivalent discret de l'énergie électromagnétique est exactement conservée. On réalise aussi une étude de convergence hp a priori ainsi qu'une étude de convergence de l'erreur sur la divergence. Des résultats numériques pour des ordres élevés montrent le potentiel de la méthode.

Mots-clés : équations de Maxwell, domaine temporel, méthode Galerkin discontinue, schéma sautemouton, maillages non-conformes,

Introduction

Time domain solutions of Maxwell's equations find applications in the applied sciences and engineering problems such as the design and optimization of antennas and radars, the design of emerging technologies (high speed electronics, integrated optics, etc.), the study of human exposure to electromagnetic waves [START_REF] Clatz | Modélisation numérique réaliste de l'exposition des tissus de la tête à un champ électromagnétique issu d'un téléphone mobile[END_REF], to name a few. These problems require high fidelity approximate solutions with a rigorous control of the numerical errors. Even for linear problems such conditions force one to look beyond standard computational techniques and seek new numerical frameworks enabling the accurate, efficient, and robust modeling of wave phenomena over long simulation times in settings of realistic geometrical complexity. Recently, discontinuous Galerkin methods have attracted much research to solve electromagnetic wave propagation problems. Being higher order versions of traditional finite volume methods [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Timedomain solution of Maxwell's equations[END_REF], discontinuous Galerkin time-domain (DGTD) methods based on discontinuous finite element spaces, easily handle elements of various types and shapes, irregular non-conforming meshes [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF]- [START_REF] Fahs | Preliminary investigation of a non-conforming discontinuous Galerkin method for solving the time domain Maxwell equations[END_REF], and even locally varying polynomial degree. They hence offer great flexibility in the mesh design, but also lead to (block-) diagonal mass matrices and therefore yield fully explicit, inherently parallel methods when coupled with explicit time stepping [START_REF] Bernacki | Parallel unstructured mesh solvers for heterogeneous wave propagation problems[END_REF]. Moreover, continuity is weakly enforced across mesh interfaces by adding suitable bilinear forms (so-called numerical fluxes) to the standard variational formulations. Whereas high-order discontinuous Galerkin time-domain methods have been developed on conforming hexahedral [START_REF] Cohen | A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell's equations in time domain[END_REF] and tetrahedral [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the heterogeneous Maxwell equations on unstructured meshes[END_REF] meshes, the design of non-conforming discontinuous Galerkin time-domain methods is still in its infancy. In practice, the non-conformity can result from a local refinement of the mesh (i.e. h-refinement), of the interpolation degree (i.e. p-enrichment) or of both of them (i.e. hp-refinement). This work is concerned with the study of high-order leap-frog schemes that are extensions of the second-order leap-frog scheme adopted in the DGTD methods that are studied in [START_REF] Fahs | Numerical evaluation of a non-conforming discontinuous Galerkin method on triangular meshes for solving the time-domain Maxwell equations[END_REF]- [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF]- [START_REF] Fahs | Preliminary investigation of a non-conforming discontinuous Galerkin method for solving the time domain Maxwell equations[END_REF]- [START_REF] Fahs | Etude de stabilité d'une méthode Galerkin discontinu pour la résolution numérique des équations de Maxwell 2D en domaine temporel sur des maillages triangulaires non-conformes[END_REF]- [START_REF] Fahs | A hp-like discontinuous Galerkin method for solving the 2D time-domain Maxwell's equations on non-conforming locally refined triangular meshes[END_REF]- [START_REF] Fahs | Development of a non-conforming discontinuous Galerkin method on simplex meshes for electromagnetic wave propagation[END_REF]. The motivation behind this study is to improve the overall accuracy for the same mesh resolution and/or to improve convergence when the mesh resolution is increased. Not surprisingly, the arbitrary high-order DGTD methods discussed in this work are consistently more accurate than the DGTD methods based on the second-order leap-frog scheme. The high-order leap-frog schemes require more computational operations to update a cell. Fortunately, this can be alleviated by the ability to use discretization meshes with fewer points per wavelength for the same level of accuracy.

The high-order leap-frog schemes considered in this work were initially proposed by Young [START_REF] Young | High-order, leapfrog methodology for the temporally dependent Maxwell's equations[END_REF], and further studied by Spachmann et al. [START_REF] Spachmann | High order explicit time integration schemes for Maxwell's equations[END_REF]. The chief attributes of these integrators are that the memory requirements are small and the algorithmic complexity is not significantly increased, with respect to the second-order leap-frog scheme. In the case of Fang's approach [START_REF] Fang | Time domain finite difference computation for Maxwell's equations[END_REF]- [START_REF] Tuomela | Fourth-order schemes for the wave equation, Maxwell's equations, and linearized elastodynamic equations[END_REF], high-order temporal derivatives are replaced with high-order mixed spatial derivatives. Although the memory requirements of Fang's approach are the same as the one of the second-order leap-frog scheme, the resulting scheme is far too complex to implement for most practical problems. One can note here that to our knowledge, these high-order leap-frog schemes have not been used yet in the context of the Maxwell's equations discretized by a discontinuous Galerkin method.

The existing high-order discontinuous Galerkin methods are mostly based on high-order Runge-Kutta (RK) schemes. The low storage RK schemes introduced in [START_REF] Carpenter | Fourth-order 2N-storage Runge-Kutta schemes[END_REF] are among the most popular choices for time integration of the DG space-discretized Maxwell equations. High-order RKDG schemes have been used by Monk and Richter [START_REF] Monk | A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media[END_REF] for solving linear symmetric hyperbolic problems, Hesthaven and Warburton [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Timedomain solution of Maxwell's equations[END_REF], Chen et al. [START_REF] Chen | High-order RKDG methods for computational electromagnetics[END_REF] and Lu et al. [START_REF] Lu | Discontinuous Galerkin methods for dispersive and lossy Maxwell's equations and PML boundary conditions[END_REF] for solving time-domain electromagnetics. A dispersion and dissipation study for a high-order DG method for solving Maxwell's equations have been conducted in [START_REF] Sármány | Dispersion and dissipation error in highorder Runge-Kutta discontinuous Galerkin discretisations of the Maxwell equations[END_REF] using several high-order RK schemes. Recently, Kanevsky et al. [START_REF] Kanevsky | Application of implicitexplicit high-order Runge-Kutta methods to discontinuous Galerkin schemes[END_REF] have applied a hybrid implicit-explicit high-order RK scheme to DG methods for solving conservation laws.

This report is structured as follows. In Sec. 2, we introduce the high-order non-conforming DGTD method for solving the system of Maxwell equations. Our two main results, the stability and the hp-convergence of the proposed method, are stated and proved in Sec. 3. In this section we also establish bounds on the behavior of the divergence error. In Sec. 4 we verify our theoretical results through numerical experiments. Finally, some concluding remarks are presented in Sec. 5.

2 An arbitrary high-order non-conforming DGTD method

We consider the Maxwell equations in three space dimensions for heterogeneous anisotropic linear media with no source. The electric permittivity tensor ¯ (x) and the magnetic permeability tensor μ(x) are varying in space, time-invariant and both symmetric positive definite. The electric field E = t (E x , E y , E z ) and the magnetic field H = t (H x , H y , H z ) verify:

¯ ∂ t E = curl H, μ∂ t H = -curl E, (1)
where the symbol ∂t denotes a time derivative. These equations are set and solved on a bounded polyhedral domain Ω of R 3 . For the sake of simplicity, a metallic boundary condition is set everywhere on the domain boundary ∂Ω, i.e. n × E = 0 (where n denotes the unitary outwards normal).

Space discretization

We consider a partition Ω h of Ω into a set of tetrahedra τ i of size h i with boundaries ∂τ i such that h = max τi∈Ω h h i . To each τ i ∈ Ω h we assign an integer p i ≥ 0 (the local interpolation degree) and we collect the p i in the vector p = {p i : τ i ∈ Ω h }. Of course, if p i is uniform in all element τ i of the mesh, we have p = p i . Within this construction we admit meshes with possibly hanging nodes i.e. by allowing non-conforming (or irregular) meshes where element vertices can lie in the interior of faces of other elements. However, we assume that the local mesh sizes and approximation degrees are of bounded variation, that is, there exist a constant κ 1 > 0, depending only on the shape-regularity of the mesh, and a constant κ 2 > 0, such that:

κ -1 1 h i ≤ h k ≤ κ 1 h i , κ -1 2 p i ≤ p k ≤ κ 2 p i , (2)
for all neighboring elements τ i and τ k in Ω h . Nevertheless, the above hypothesis is not restrictive in practice and allows, in particular for geometric refinement and linearly increasing approximation degrees. We also assume that Ω h is shape regular in the sense that there is a constant η > 0 such that:

∀ τ i ∈ Ω h , h i ρ i ≤ η, (3)
where ρ i is the diameter of the insphere of τ i . Each tetrahedron τ i is assumed to be the image, under a smooth bijective (diffeomorphic) mapping, of a fixed reference tetrahedron τ = {x, ŷ, ẑ| x, ŷ, ẑ ≥ 0; x + ŷ + ẑ ≤ 1}. Assuming that τ i is a straight faced tetrahedron defined through the coordinates of the four vertices v i 1 , v i 2 , v i 3 and v i 4 (see Fig. 1), the correspondence between the two tetrahedra τ and τ i is established through the use of the barycentric coordinates (λ 1 , λ 2 , λ 3 , λ 4 ). We recall that any point x i ∈ τ i can be expressed as a convex combination of the vertices of τ i and the mapping is simply given by χ i : (x, ŷ, ẑ) ∈ τ → x i , such that:

x i (x, ŷ, ẑ) = λ 1 v i 1 + λ 2 v i 2 + λ 3 v i 3 + λ 4 v i 4 , where λ 1 + λ 2 + λ 3 + λ 4 = 1 and 0 ≤ (λ 1 , λ 2 , λ 3 , λ 4 ) ≤ 1 with λ 1 = 1 -x -ŷ -ẑ, λ 2 = x, λ 3 = ŷ and λ 4 = ẑ.
For each τ i , V i denotes its volume, and ¯ i and μi are respectively the local electric permittivity and magnetic permeability tensors of the medium, which could be varying inside the element τ i . For two distinct tetrahedra τ i and τ k in Ω h , the intersection τ i ∩ τ k is an (oriented) triangle a ik which we will call interface. For each internal interface a ik , we denote by S ik the measure of a ik and by n ik the unitary normal vector, oriented from τ i towards τ k . For the boundary interfaces, the index k corresponds to a fictitious element outside the domain. We denote by F I h the union of all interior INRIA faces of Ω h , by F B h the union of all boundary faces of Ω h , and by F h = F I h ∩ F B h . Furthermore, we identify F B h to ∂Ω since Ω is a polyhedron. Finally, we denote by V i the set of indices of the elements which are neighbors of τ i (having an interface in common). We also define the perimeter P i of τ i by P i = k∈Vi S ik . We have the following geometrical property for all elements:

k∈Vi S ik n ik = 0.

v v v v i i i i 1 2 3 4 x y z τ i (0,0,0) (0,1,0) (1,0,0) (0,0,1) τ z y x (x,y,z) = χ i (x,y,z) (x,y,z) = χ i -1 (x,y,z)
Figure 1: Mapping between the physical tetrahedron τ i and the reference tetrahedron τ .

In the following, for a given partition Ω h and for a vector p, we seek approximate solutions to [START_REF] Babuska | The hp-version of the finite element method with quasiuniform meshes[END_REF] in the finite dimensional subspace

V p (Ω h ) = { v ∈ L 2 (Ω) 3 : v |τi ∈ P pi (τ i ), ∀τ i ∈ Ω h }, where P pi (τ i )
denotes the space of nodal polynomials of degree at most p i inside the element τ i . Note that the polynomial degree p i may vary from element to element in the mesh. By non-conforming interface we mean an interface a ik which is such that at least one of its vertices is a hanging node, or/and such that p i|a ik = p k|a ik .

Following the discontinuous Galerkin approach, the electric and magnetic fields inside each finite element are seeked for as linear combinations ( E i , H i ) of linearly independent basis vector fields ϕ ij , 1 ≤ j ≤ d i , where d i denotes the local number of degrees of freedom (DOF) inside τ i . We denote by

P i = Span( ϕ ij , 1 ≤ j ≤ d i )
. The approximate fields ( E h , H h ), defined by (∀i, E h|τi = E i , H h|τi = H i ) are allowed to be completely discontinuous across element boundaries. For such a discontinuous field U h , we define its average { U h } ik through any internal interface a ik , as { U h } ik = ( U i|a ik + U k|a ik )/2 . Note that for any internal interface a ik , { U h } ki = { U h } ik . Because of this discontinuity, a global variational formulation cannot be obtained. However, dot-multiplying (1) by any given vector function ϕ ∈ P i , integrating over each single element τ i and integrating by parts, yields:

         τi ϕ • ¯ i ∂ t E = τi curl ϕ • H - ∂τi ϕ • ( H × n), τi ϕ • μi ∂ t H = - τi curl ϕ • E + ∂τi ϕ • ( E × n). (4)
In Eq. ( 4), we now replace the exact fields E and H by the approximate fields E h and H h in order to evaluate volume integrals. For integrals over ∂τ i , a specific treatment must be introduced since the approximate fields are discontinuous through element faces. We choose to use a fully centered numerical flux, i.e. ∀i, ∀k ∈

V i , E |a ik { E h } ik , H |a ik { H h } ik .
The metallic boundary condition on a boundary interface a ik (where k is the element index of a fictitious neighboring element) is dealt with weakly, in the sense that traces of fictitious fields E k and H k are used for the computation of numerical fluxes for the boundary element τ i . In the present case, where all boundaries are metallic, we simply take E k|a ik = -E i|a ik and H k|a ik = H i|a ik . Replacing surface integrals using the centered numerical flux in (4) and re-integrating by parts yields:

           τi ϕ • ¯ i ∂ t E i = 1 2 τi (curl ϕ • H i + curl H i • ϕ) - 1 2 k∈Vi a ik ϕ • ( H k × n ik ), τi ϕ • μi ∂ t H i = - 1 2 τi (curl ϕ • E i + curl E i • ϕ) + 1 2 k∈Vi a ik ϕ • ( E k × n ik ). (5)
We can rewrite this formulation in terms of scalar unknowns. Inside each element, the fields being recomposed according to

E i = 1≤j≤di E ij ϕ ij and H i = 1≤j≤di H ij ϕ ij .
Let us denote by E i and H i respectively the column vectors (E il ) 1≤l≤di and (H il ) 1≤l≤di . Eq. ( 5) can be rewritten as:

       M i ∂ t E i = K i H i - k∈Vi S ik H k , M µ i ∂ t H i = -K i E i + k∈Vi S ik E k , (6)
where the symmetric positive definite mass matrices M σ i (σ stands for or µ) and the symmetric stiffness matrix K i (all of size d i × d i ) are given by :

(M σ i ) jl = τi t ϕ ij • σi ϕ il , (K i ) jl = 1 2 τi t ϕ ij • curl ϕ il + t ϕ il • curl ϕ ij .
For any interface a ik , the d i × d k rectangular matrix S ik is given by:

(S ik ) jl = 1 2 a ik t ϕ ij • ( ϕ kl × n ik ), 1 ≤ j ≤ d i , 1 ≤ l ≤ d k . (7)

Time discretization

The choice of the time discretization is a crucial step for the global efficiency of the numerical method.

The temporal integration methods are divided into two major families: implicit and explicit schemes. Implicit schemes require the solution of large matrix systems resulting in a high computational effort per time step. The quality of the scheme depends strongly on the efficiency of the used linear system solver. The advantage of implicit schemes is their flexibility regarding the choice of the time step since usually, these time schemes are unconditionally stable. Thus, an analysis requires only a small number of solver runs but every time step is burdened by a high numerical effort. Explicit schemes in contrast are easy to implement, produce greater accuracy with less computational effort than implicit methods, but are restricted by a stability criterion enforcing a close linkage of the time step to the spatial discretization parameter. This restriction may result in a large number of iterations per analysis, each iteration with a low computational effort. In this study, we focus on explicit time integration schemes and our objective is to design an arbitrary high-order DGTD method which combines the spatial discretization features of the methods discussed in [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF]- [START_REF] Fahs | Preliminary investigation of a non-conforming discontinuous Galerkin method for solving the time domain Maxwell equations[END_REF]- [START_REF] Fahs | A hp-like discontinuous Galerkin method for solving the 2D time-domain Maxwell's equations on non-conforming locally refined triangular meshes[END_REF], with a familty of high-oder explicit leap-frog schemes. The ordinary differential system (5) can be formally seen as a system of the form (see Sec. 3.1 for more details):

M∂ t E = UH, N∂ t H = VE. (8)
The one-step explicit time integration methods like Runge-Kutta or leap-frog schemes imply the convenience of storing just one old solution followed by a single update step. Therefore, they are INRIA computationally efficient per update cycle and easy to implement. Moreover, the leap-frog scheme has the advantage to be free of time dissipation. We can introduce the N th-order explicit leap-frog (LF N ) integrator as an approximation of the solution of the first-order ODE:

ẏ(t) = Ay(t) ⇒ y(t) = e At y(t 0 ), (9)
with y(t 0 ) as initial value. The time discrete equivalent of ( 9) is given by:

y(n∆t) = e A∆t y((n -1)∆t). ( 10
)
System ( 8) can be rewritten as:

∂ t H E = 0 N -1 V M -1 U 0 A H E Y(t) . ( 11
)
Note that the system matrix A depends only on the spatial configuration. Seeking a time discrete solution of [START_REF] Fahs | A hp-like discontinuous Galerkin method for solving the 2D time-domain Maxwell's equations on non-conforming locally refined triangular meshes[END_REF], a discretization in time with a global time step is introduced. The time discrete solution of the first-order system of ODEs ( 11) is a discretized version of the exponential solution according to its scalar equivalent [START_REF] Fahs | Etude de stabilité d'une méthode Galerkin discontinu pour la résolution numérique des équations de Maxwell 2D en domaine temporel sur des maillages triangulaires non-conformes[END_REF]:

Y(n∆t) = Φ(∆t)Y((n -1)∆t), (12) with:
Φ(∆t) = ∞ i=0 ∆t i i! A i := e A∆t . (13)
Finally, the solution of ( 11) is written as:

H(n∆t) E(n∆t) = Φ 11 Φ 12 Φ 21 Φ 22 Φ H((n -1)∆t) E((n -1)∆t) Y((n-1)∆t) . ( 14
)
Introducing a staggered temporal grid, as in the case of LF schemes, we obtain the general LF update equation:

H n+1 = [Φ 2 11 -Φ 11 Φ 12 Φ -1 22 Φ 21 ]H n-1 + [Φ 12 + Φ 11 Φ 12 Φ -1 22 ]E n , (15a)
E n+2 = [Φ 21 + Φ 22 Φ 21 Φ -1 11 ]H n+1 + [Φ 2 22 -Φ 22 Φ 21 Φ -1 11 Φ 12 ]E n , (15b)
where the electric field quantities are located at even time steps and magnetic quantities at odd time steps as illustrated on Fig. 2. In the case of a non-conducting material, the relation:

[Φ 2 11 -Φ 11 Φ 12 Φ -1 22 Φ 21 ] = [Φ 2 22 -Φ 22 Φ 21 Φ -1 11 Φ 12 ] = I + O(∆t N +1 ), (16)
holds, which is a characteristic of a LF scheme.

In the sequel, superscripts refer to time stations and ∆t is the global time step. The unknowns related to the electric field are approximated at integer time-stations t n = n∆t and are denoted by E n i . The unknowns related to the magnetic field are approximated at half-integer time-stations t n+1/2 = (n + 1/2)∆t and are denoted by Note that here, the used time step ∆t is in comparison to the time step defined in [START_REF] Fahs | Etude de stabilité d'une méthode Galerkin discontinu pour la résolution numérique des équations de Maxwell 2D en domaine temporel sur des maillages triangulaires non-conformes[END_REF] twice as large.

H n+ 1 2 i . RR n 6699 (n-2) (n-1) (n+1) n ∆ ∆ ∆ ∆t t t t H E t leap-frog time axis
The LF N for N = 2 and 4 integrators are constructed as follows [START_REF] Young | High-order, leapfrog methodology for the temporally dependent Maxwell's equations[END_REF]- [START_REF] Spachmann | High order explicit time integration schemes for Maxwell's equations[END_REF]:

                                       T 1 = ∆t(M i ) -1 curl H n+ 1 2 i , T 1 = -∆t(M µ i ) -1 curl E n+1 i , T 2 = -∆t(M µ i ) -1 curlT 1 , T 2 = ∆t(M i ) -1 curlT 1 , T 3 = ∆t(M i ) -1 curlT 2 , T 3 = -∆t(M µ i ) -1 curlT 2 . LF 2 : E n+1 i = E n i + T 1 , H n+ 3 2 i = H n+ 1 2 i + T 1 . LF 4 : E n+1 i = E n i + T 1 + T 3 /24, H n+ 3 2 i = H n+ 1 2 i + T 1 + T 3 /24. ( 17
)
Here the T i and T i (i = 1, 2, 3) can be combined. Moreover, T 3 (resp. T 3 ) is a temporary vector. However, in a practical implementation, this vector is not used since the calculations for T 3 and E n+1 i (resp.T 3 and H

n+ 1 2 i
) can be combined. Thus, the LF 4 scheme requires 2 times more memory storage and 3 times more arithmetic operations than the LF 2 scheme. In general, the LF N scheme requires N/2 times more memory storage and (N -1) times more arithmetic operations than the LF 2 scheme. A pictorial representation of the extended LF 4 integrator is shown on Fig. 3.

CURL CURL CURL CURL CURL CURL Σ Σ E E H H n n n n+1 -1/2 +1/2
Figure 3: A pictorial representation of the extended LF 4 scheme (multiplicative constants omitted). Figure taken from [START_REF] Young | High-order, leapfrog methodology for the temporally dependent Maxwell's equations[END_REF].

For the treatment of the boundary condition on an interface a ik ∈ F B h , we use:

E n k|a ik = -E n i|a ik and H n+ 1 2 k|a ik = -H n+ 1 2 i|a ik (18)

Stability and convergence analysis

In this section we study the stability and convergence properties of the high-order discontinuous Galerkin method introduced previously.

INRIA

Stability

Here, we aim at giving and proving a sufficient condition for the L 2 -stability of the proposed high-order DGTD method. We use the same kind of energy approach as in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the heterogeneous Maxwell equations on unstructured meshes[END_REF] where a quadratic form plays the role of a Lyapunov function of the whole set of numerical unknowns. To this end, we suppose that all electric (resp. magnetic) unknowns are gathered in a column vector

E (resp. H) of size d = i d i .
Then the space discretized system (6) can be rewritten as:

M ∂ t E = KH -AH -BH, M µ ∂ t H = -KE + AE -BE, (19)
where we have the following definitions and properties:

M , M µ and K are d × d block diagonal matrices with diagonal blocks equal to M i , M µ i and K i respectively. Therefore M and M µ are symmetric positive definite matrices, and K is a symmetric matrix.

A is also a d × d block sparse matrix, whose non-zero blocks are equal to S ik when a ik is an internal interface of the mesh. Since n ki = -n ik , it can be checked from ( 7) that (S ik ) jl = (S ki ) lj and then S ki = t S ik ; thus A is a symmetric matrix.

B is a d × d block diagonal matrix, whose non-zero blocks are equal to S ik when a ik is a metallic boundary interface of the mesh. In that case, (S ik ) jl = -(S ik ) lj ; thus B is a skew-symmetric matrix.

The discontinuous Galerkin DGTD-P pi method using centered fluxes combined with a N th order leap-frog (LF N ) time scheme can be written, using the matrix S = K -A -B, in the general form:

       M E n+1 -E n ∆t = S N H n+ 1 2 , M µ H n+ 3 2 -H n+ 1 2 ∆t = -t S N E n+1 , (20)
where the matrix S N (N being the order of the leap-frog scheme) verifies:

S N =                  S if N = 2, S(I - ∆t 2 24 M -µ t SM -S) if N = 4, S I + N/2-1 i=1 (-1) i (2i + 1)!2 2i (∆t 2 M -µ t SM -S) i ∀ N > 2, even. (21)
We now define the following discrete version of the electromagnetic energy.

Definition 3.1 We consider the following electromagnetic energies inside each tetrahedron and in the whole domain Ω:

• the local energy : ∀i, E n i = 1 2 ( t E n i M i E n i + t H n-1 2 i M µ i H n+ 1 2 i ), (22)
• the global energy :

E n = 1 2 ( t E n M E n + t H n-1 2 M µ H n+ 1 2 ). ( 23
)
In the following, we shall prove that the global energy ( 23) is conserved through a time step and that it is a positive definite quadratic form of all unknowns under a CFL-like condition on the global time step ∆t.

Lemma 3.1 Using the DGTD-P pi method ( 20)-( 21) for solving [START_REF] Babuska | The hp-version of the finite element method with quasiuniform meshes[END_REF] with metallic boundaries only, the global discrete energy ( 23) is exactly conserved, i.e. E n+1 -E n = 0, ∀ n.

Proof. We denote by

E n+ 1 2 = E n+1 + E n 2 .
We have :

E n+1 -E n = t E n+ 1 2 M (E n+1 -E n ) + 1 2 t H n+ 1 2 M µ (H n+ 3 2 -H n-1 2 ) = ∆t t E n+ 1 2 S N H n+ 1 2 - 1 2 ∆t t H n+ 1 2 ( t S N E n+1 + t S N E n ) = ∆t t H n+ 1 2 ( t S N -t S N )E n+ 1 2 = 0.
This concludes the proof.

Lemma 3.2 Using the DGTD-P pi method ( 20)-( 21), the global discrete electromagnetic energy E n (23) is a positive definite quadratic form of all unknowns if:

∆t ≤ 2 d N , with d N = M -µ 2 t S N M - 2 , ( 24
)
where . denote the canonical norm of a matrix (∀X, AX ≤ A X ), and the matrix M

-σ 2
is the inverse square root of M σ . Also, for a given mesh, the stability limit of the LF 4 scheme is roughly 2.85 times larger than that of the LF 2 scheme.

Proof. The mass matrices M and M µ are symmetric positive definite and we can construct in a simple way their square root (also symmetric positive definite) denoted by M 2 and M µ 2 respectively. Moreover:

E n = 1 2 t E n M E n + 1 2 t H n-1 2 M µ H n+ 1 2 = 1 2 t E n M E n + 1 2 t H n-1 2 M µ H n-1 2 - ∆t 2 t H n-1 2 t S N E n ≥ 1 2 M 2 E n 2 + 1 2 M µ 2 H n-1 2 2 - ∆t 2 | t H n-1 2 M µ 2 M -µ 2 t S N M - 2 M 2 E n | ≥ 1 2 M 2 E n 2 + 1 2 M µ 2 H n-1 2 2 - d N ∆t 2 M µ 2 H n-1 2 M 2 E n .
At this point, we choose to use an upper bound for the term

M µ 2 H n-1 2
M 2 E n which might lead to suboptimal lower bounds for the energy (and then to a slightly too severe stability limit for the DGTD method). Anyway, this stability limit is only sufficient, and not really close to necessary. We use the inequality:

M µ 2 H n-1 2 M 2 E n ≤ 1 2 ( M µ 2 H n-1 2 2 + M 2 E n 2 ).
We then sum up the lower bounds for E n to obtain:

E n ≥ 1 2 (1 - d N ∆t 2 ) M 2 E n 2 + 1 2 (1 - d N ∆t 2 ) M µ 2 H n-1 2 2 .
Then, under the condition proposed in Lemma 3.2, the electromagnetic energy E n is a positive definite quadratic form of all unknowns. For a given mesh, using the definition (21) of S N , the LF 4 scheme is stable if:

INRIA ∆t M -µ 2 t S 4 M - 2 ≤ 2, ⇒ ∆t M -µ 2 (S 2 - ∆t 2 24 S 2 M -µ t S 2 M -S 2 )M - 2 ≤ 2, ⇒ |∆td 2 - ∆t 3 24 d 3 2 | ≤ 2.
This inequality is verified if and only if

d 2 ∆t ≤ 2( 3 √ 2 + 3 √ 4) 2(2.847
). This concludes the proof. Now, we denote by ν N = CFL N /CFL 2 the ratio between the stability limit of the LF N scheme and the LF 2 scheme, and by r N = ν N /(N/2) the ratio between ν N and the additional memory storage between the LF N and LF 2 schemes. Tab. 1 lists the values of ν N and r N for several values of N . As it can be seen from Tab. 1, the choice of the LF 4 scheme is advantageous with respect to the r N ratio. Now, our objective is to give an explicit CFL condition on ∆t under which the local energy ( 22) is a positive definite quadratic form of the numerical unknowns E n i and H n-

1 2 i
. We first need some classical definitions. Definition 3.2 We assume that the media is isotropic and such that the tensors ¯ i and μi are piecewise constant, i.e. ¯ i = i and μi = µ i . We denote by c i = 1/ √ i µ i the propagation speed in the element τ i . We also assume that there exist dimensionless constants α i and β ik (k ∈ V i ) such that:

∀ X ∈ P i ,        curl X τi ≤ α i P i V i X τi , X 2 a ik ≤ β ik S ik V i X 2 τi , (25)
where X τi and X a ik denote the L 2 -norms of the vector field X over τ i and the interface a ik respectively.

Lemma 3.3 Using the scheme ( 6)-( 17)- [START_REF] Monk | A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media[END_REF], under assumptions of Definition 3.2, the local discrete energy

E n i (22) is a positive definite quadratic form of all unknowns (E n i , H n-1 2 i
) and the scheme is stable if the time step ∆t is such that:

∀i, ∀k ∈ V i , c i ∆t[2α i + β ik ] < 4V i P i , (26)
(with the convention that, in the above formula, k should be replaced by i for a metallic boundary interface a ik ).

Proof. Using the scheme (5) to replace the occurrences of H

n+ 1 2 i
in the definition of E i , and using the boundary fluxes given in [START_REF] Monk | A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media[END_REF], we get:

E n i = i 2 E n i 2 τi + i 2 H n-1 2 i 2 τi - ∆t 4 X n i , with X n i = Ti curl H n-1 2 i • E n i + curl E n i • H n-1 2 i - k∈Vi a ik ( H n-1 2 i × E n k ) • n ik .

RR n 6699

For any metallic or internal interface a ik , we have:

a ik ( H n-1 2 i × E k ) • n ik ≤ 1 √ µ i i a ik √ µ i H n-1 2 i √ i E n k ≤ 1 2
µ i i H n-1 2 i 2 a ik + 1 2 i µ i E n k 2 a ik ≤ 1 2
µ i i β ki S ik V i H n-1 2 i 2 τi + 1 2 i µ i β ki S ik V k E n k 2 τi .
In the remainder of the proof, we omit the superscripts n and n-1/2 in the electric and magnetic variables respectively. We have that:

|X n i | ≤ curl H i τi E i τi + curl E i τi H i τi + 1 2 k∈Vi µ i i H i 2 a ik + i µ i E k 2 a ik ≤ 2α i P i V i H i τi E i τi + 1 2 k∈Vi µ i i β ik S ik V i H i 2 τi + i µ i β ki S ik V k E k 2 τi .
Noticing that

H i τi E i τi ≤ c i 2 (µ i H i 2 τi + i E i 2 τi
), gathering all lower bounds for terms in the expression of E n i and using P i = k∈Vi S ik leads to:

E n i ≥ k∈Vi S ik 1 2P i - α i c i ∆t 4V i ( i E i 2 τi + µ i H i 2 τi ) - ∆t 8 k∈Vi S ik µ i i β ik V i H i 2 τi + i µ i β ki V k E k 2 τi .
Then, summing up these inequalities in order to obtain a lower bound for i E i leads to an expression that we reorganize using sums over interfaces a ik . We find that

i E i ≥ a ik
S ik W ik with:

W ik = i E i 2 τi 1 2P i - α i c i ∆t 4V i - β ik c i ∆t 8V i + µ i H i 2 τi 1 2P i - α i c i ∆t 4V i - β ik c i ∆t 8V i + k E k 2 τ k 1 2P k - α k c k ∆t 4V k - β ki c k ∆t 8V k + µ k H k 2 τ k 1 2P k - α k c k ∆t 4V k - β ki c k ∆t 8V k .
Then, under the conditions of Lemma 3.3, W ik is a positive definite quadratic form and the local energy is a positive definite quadratic form of all unknowns. This concludes the proof. Note that, the existence of the constants α i and β ik (k ∈ V i ) is always ensured. The values of α i only depend on the local polynomial order p i while the values of β ik depend on p i and on the number of hanging nodes on the interface a ik . For instance, for orthogonal polynomials on a dsimplex β ik = (p i + 1)(p i + d)/d (see [START_REF] Warburton | On the constants in hp-finite element trace inequalities[END_REF]), and for arbitrary basis functions these values are given by:

INRIA ( α 2 i P 2 i V 2 i ; β ik S ik V i ) = ( M -1/2 S 1 M -1/2 ; M -1/2 S 2 M -1/2 ),
where M is the mass matrix without material parameter, S 2 = 2S ik , and

S 1 = τi curl ϕ ij •curl ϕ il , 1 ≤ j, l ≤ d i .
Moreover, the value of β ik verifies the properties of Lemma 2 of [START_REF] Fahs | Etude de stabilité d'une méthode Galerkin discontinu pour la résolution numérique des équations de Maxwell 2D en domaine temporel sur des maillages triangulaires non-conformes[END_REF].

Convergence

A convergence analysis of the LF 2 based DGTD-P pi method is conducted in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the heterogeneous Maxwell equations on unstructured meshes[END_REF] in the case of conforming simplicial meshes and p i = p everywhere. In this section, our objective is to obtain an a priori error estimates depending on h and p, which establishes the rate of convergence of the proposed hp-like DGTD method. We consider again the Maxwell problem: find ( E, H)

: Ω×]0, T [→ R 3 × R 3 such that,    ¯ ∂ t E = curl H in Ω, μ∂ t H = -curl E in Ω, (27a)
         n × E = 0 on ∂Ω, E( x, 0) = E 0 ( x) in Ω, H( x, 0) = H 0 ( x) in Ω, (27b)
   ∇ • E = 0 in Ω, ∇ • H = 0 in Ω. ( 27c
)
We assume that ¯ , μ ∈ [L ∞ (Ω)] 3×3 and ∃ C 1 , C 2 > 0 such that:

∀ ξ ∈ R 3 :    C 1 | ξ| 2 ≤ ¯ ξ • ξ ≤ C 2 | ξ| 2 , C 1 | ξ| 2 ≤ μ ξ • ξ ≤ C 2 | ξ| 2 . ( 28
)
The problem (27) admits a unique solution ( E, H) ∈ [C 1 (0, T ; [L 2 (Ω)] 3 ) ∩ C 0 (0, T ; H 0 (curl, Ω))] 2 (see [START_REF] Remaki | A new finite volume scheme for solving Maxwell system[END_REF] for more details), where

H 0 (curl, Ω) = { u ∈ H(curl, Ω) such that n × u = 0}.
For a real s ≥ 0, we define the classical broken space :

H s (Ω h ) = {v ∈ L 2 (Ω) : ∀τ i ∈ Ω h , v |τi ∈ H s (τ i )}. ( 29
)
The space H s (Ω h ) is equipped with the natural norm, for v ∈ H s (Ω h ) :

v s,h = τi∈Ω h v 2 s,τi 1 2 , ( 30
)
where . s,τi is the usual Sobolev norm of H s on τ i . For s > 1 2 , the elementwise traces of functions in H s (Ω h ) belongs to tr(F h ) = Π τi∈Ω h L 2 (∂τ i ). We denote by H s (Ω h ) the vectorial broken space [H s (Ω h )] 3 and the associated norm defined by :

v s,h = 3 j=1 v j 2 s,h 1 2 , ( 31
)
where v = (v 1 , v 2 , v 3 ) ∈ H s (Ω h ).
We define the jump of a function v ∈ H s (Ω h ):

∀a ik ∈ F I h , [[ v]] i ik = [[ v]] τi a ik = ( v k|a ik -v i|a ik ) × n ik , ∀a ik ∈ F B h , [[ v]] i ik = -v i|a ik × n ik . (32)
We associate to the continuous problem (27a) the following space discretized problem: find ( E(., t), H(., t))

∈ H 1 (Ω h ) × H 1 (Ω h ) such that, ∀ τ i ∈ Ω h and ∀ φ, ψ ∈ H 1 (Ω h ),                          τi φ i • ¯ i ∂ t E i - τi H i • curl φ i + k ∈ V i a ik ∈ F I h a ik φ i • ( H |a ik × n ik ) + k ∈ V i a ik ∈ F B h a ik φ i • ( H |a ik × n ik ) = 0, τi ψ i • μi ∂ t H i + τi E i • curl ψ i - k ∈ V i a ik ∈ F I h a ik ψ i • ( E |a ik × n ik ) = 0, ( 33
)
where φ i = φ |τi and ψ i = ψ |τi . Summing up the identities in (33) with respect to i, we consider the following semi-discrete discontinuous Galerkin problem : find ( E h (., t),

H h (., t)) ∈ V p (Ω h ) × V p (Ω h ) such that, ∀ τ i ∈ Ω h and ∀ φ h , ψ h ∈ V p (Ω h ),                    i τi φ hi • ¯ i ∂ t E i - i τi H i • curl φ hi + a ik ∈F h a ik [[ φ h ]] i ik • { H h } ik = 0, i τi ψ hi • μi ∂ t H i + i τi E i • curl ψ hi - a ik ∈F h a ik [[ ψ h ]] i ik • { E h } ik = 0, E h (0) = Π p h E 0 and H h (0) = Π p h H 0 . (34)
Here Π p h :

L 2 (Ω) → V p (Ω h ) is the L 2 -projection onto V p (Ω h ).
The problem (34) can be rewritten in the following form: find

U h = ( E h , H h ) ∈ V p (Ω h ) × V p (Ω h ) such that: J(∂ t U h , U h ) + a( U h , U h ) + b( U h , U h ) = 0, ∀ U h = ( φ h , ψ h ) ∈ V p (Ω h ) × V p (Ω h ). ( 35
)
For W = ( u, v) and W = ( u , v ), the bilinear forms J, a and b defined on V p (Ω h ) × V p (Ω h ) are given by:

                     J( W, W ) = i τi ¯ u • u + μ v • v , a( W, W ) = i τi u • curl v -v • curl u , b( W, W ) = a ik ∈F h a ik { v} • [[ u ]] -{ u} • [[ v ]] , (36)
taking into account that, for boundary faces a ik ∈ F B h we have { v} = v. The semi-discrete discontinuous Galerkin formulation (35) is consistent with the original continuous problem (27) in the following sense: if U = ( E, H) is the exact solution of (27), such that ∀ h, ∀ t ∈ [0, T ], ( E(., t), H(., t)) ∈ H s (Ω) × H s (Ω), then we have:

J(∂ t U, U ) + a( U, U ) + b( U, U ) = 0, ∀ U ∈ V p (Ω h ) × V p (Ω h ). ( 37
)
The following approximation results will be used to bound the error [START_REF] Babuska | The hp-version of the finite element method with quasiuniform meshes[END_REF]- [START_REF] Schwab | p-and hp-finite element methods. Theory and applications to solid and fluid mechanics[END_REF] INRIA Lemma 3.4 (Babuska and Suri [START_REF] Babuska | The hp-version of the finite element method with quasiuniform meshes[END_REF]) Let τ i ∈ Ω h and suppose that u ∈ H s+1 (τ i ) for s ≥ 0. Let Π be a linear continuous operator from H s+1 (τ i ) onto P pi (τ i ), p i ≥ 1, such that Π( u) = u, ∀ u ∈ P pi (τ i ).

Then we have:

u -Π( u) 0,τi ≤ C h νi+1 i p s+1 i u s+1,τi , (38)
u -Π( u) 0,∂τi ≤ C h νi+ 1 2 i p s+ 1 2 i u s+1,τi , (39)
where ν i = min{s, p i } and C is a positive constant independent of u, h i and p i , but dependent on s and on the shape regularity of the mesh parameter η.

Lemma 3.5 (Schwab [21]) For all q ∈ P pi (τ i ), p i ≥ 1, we have:

q 2 0,∂τi ≤ C inv p 2 i h i q 2 0,τi ,
where C inv is a positive constant depending only on the shape regularity of the mesh parameter η.

Let U = ( E, H) and U h = ( E h , H h ). We denote by ε ε ε τi (t) the local error and by ε ε ε(t) = τi∈Ωj ε ε ε τi (t) the global error. Then we have:

ε ε ε τi (t) = E -Π p h E + Π p h E -E h 2 0,τi + H -Π p h H + Π p h H -H h 2 0,τi ≤ 2( E -Π p h E 2 0,τi + H -Π p h H 2 0,τi ) + 2( Π p h E -E h 2 0,τi + Π p h H -H h 2 0,τi ) = 2 U -Π p h U 2 0,τi + 2 Π p h U -U h 2 0,τi = 2ε ε ε a τi + 2ε ε ε b τi ,
where ε ε ε a τi is due to the error introduced by the polynomial approximation of the exact solution while ε ε ε b τi measures the errors associated with the semi-discrete approximation of Maxwell's equations.

To bound ε ε ε a τi we need only recall Lemma 3.4 to state

Lemma 3.6 Assume that U ∈ H s+1 (τ i ) × H s+1 (τ i ).
Then there exists a constant C, dependent on s and on the shape regularity of the mesh η, but independent of U, h i and p i , such that:

U -Π p h U 0,τi ≤ C h νi+1 i p s+1 i U s+1,τi , (40)
where ν i = min{s, p i } and s ≥ 0.

Theorem 3.1 Assume that a solution ( E(t), H(t)) ∈ H s+1 (τ i ) × H s+1 (τ i ) with s ≥ 1/2 to Maxwell's equations in Ω h = i τ i exists. Then the numerical solution, ( E h (t), H h (t)) ∈ V p (Ω h ) × V p (Ω h
), to the semi-discrete approximation (34) converges to the exact solution and the global error is bounded as:

E -E h 2 0,Ω + H -H h 2 0,Ω 1 2 ≤ C h ν+1 p s+1 min + T h ν p s-1 2 min max t∈[0,T ] E(t), H(t) s+1,Ω , (41)
where ν = min{s, p min } and p min = min{p i , τ i ∈ Ω h }, p i ≥ 1. The constant C > 0 depends on the material properties and on the shape regularity of the mesh parameter η, but not on p min and h.

Proof. Let q = U -U h . Since Π p h U h = U h , we have i ε ε ε b τi = Π p h q 2 0,Ω .
To obtain a bound for Π p h q 0,Ω , we introduce σ(t) = 1 2 J(Π p h q(t), Π p h q(t)) with Π p h q(., t) belongs to V p (Ω h ) × V p (Ω h ). Using the discrete initial conditions of (34), we have σ(0) = 0 and then, for 0 < t ≤ T ,

σ(t) = 1 2 t 0 d ds J(Π p h q(s), Π p h q(s))ds = t 0 J(∂ s Π p h q(s), Π p h q(s))ds.
For any

U h ∈ V p (Ω h ) × V p (Ω h ), we have a( U h , U h ) + b( U h , U h ) = 0
, and we get:

σ(t) = t 0 ( J(∂ s Π p h q(s), Π p h q(s)) + a(Π p h q(s), Π p h q(s)) + b(Π p h q(s), Π p h q(s)) ) ds. (42)
Subtracting (35) from the consistency result (37) with U = U h = Π p h q(s) yields:

J(∂ s q(s), Π p h q(s)) + a( q(s), Π p h q(s)) + b( q(s), Π p h q(s)) = 0. ( 43
)
Now, subtracting the above equality (43) from (42) leads to:

σ(t) = t 0 J([Π p h ∂ s U -∂ s U](s), Π p h q(s)) + a([Π p h U -U](s), Π p h q(s)) + b([Π p h U -U](s), Π p h q(s)) ds. Since Π p h is a projector onto V p (Ω h ) and Π p h q(., t) belongs to V p (Ω h ) × V p (Ω h ), then J(Π p h ∂ s U - ∂ s U, Π p h q) = 0 and a(Π p h U -U, Π p h q) = 0.
Using the lower bound C 1 > 0 of ¯ and μ (28), we thus get:

C 1 2 Π p h q(t) 2 0,Ω ≤ t 0 b([Π p h U -U](s), Π p h q(s))ds. (44)
Now, we bound the surface integrals deriving from the definition of b(., .). We assume that q = ( q E , q H ), where q E and q H denote the error in E and H respectively. Let a ik ∈ F I h be an internal interface shared by the tetrahedra τ i and τ k . We denote by

I E = a ik {Π p h H -H} ik • [[Π p h q E ]
] ik , we have, using the Cauchy-Schwarz-Buniakovsky (CSB) inequality:

I E ≤ 1 2 a ik (Π p h H i -H i ) + (Π p h H k -H k ) 2 1 2 a ik n ik (Π p h q E ) i|a ik + n ki (Π p h q E ) k|a ik 1 2 ≤ Π p h H i -H i 2 0,a ik + Π p h H k -H k 2 0,a ik 1 2 (Π p h q E ) i 2 0,a ik + (Π p h q E ) k 2 0,a ik 1 2 .
Using Lemma 3.4 and Lemma 3.5, yields:

I E ≤ C h νi+ 1 2 i p s+ 1 2 i H s+1,τi + C h ν k + 1 2 k p s+ 1 2 k H s+1,τ k 1 2 C p 2 i h i Π p h q E 2 0,τi + C p 2 k h k Π p h q E 2 0,τ k 1 2 .
According to the assumptions (2), we finally get:

I E ≤ K(κ 1 , κ 2 ) h ν i p s-1 2 i H 2 s+1,τi + H 2 s+1,τ k 1 2 Π p h q E 2 0,τi + Π p h q E 2 0,τ k 1 2 , (45)

INRIA

where K > 0 does not depend on h i and p i , but depends on κ 1 and κ 2 , and on the local material properties ( ¯ i/k , μi/k ) associated to τ i and τ k .

The term

I H = a ik {Π p h E -E} ik • [[Π p h q H ]
] ik is treated in the same way, yielding the result:

I H ≤ K(κ 1 , κ 2 ) h ν i p s-1 2 i E 2 s+1,τi + E 2 s+1,τ k 1 2 Π p h q H 2 0,τi + Π p h q H 2 0,τ k 1 2 . ( 46
)
For boundary interfaces a ik ∈ F B h , we obtain the same upper bounds as (45) and ( 46) but without the norms on τ k . Summing up with respect to all τ i ∈ Ω h , and using the CSB inequality, yields:

b([Π p h U -U](s), Π p h q(s)) ≤ K(κ 1 , κ 2 ) h ν p s-1 2 min Π p h q(s) 0,Ω ( E(s), H(s)) s+1,Ω . (47)
Integrating in t ∈ [0, T ] and combining this with Lemma 3.6 establishes the result and proves convergence on weak assumptions of local, elementwise smoothness of the solution.

Corollary 3.1 Under the assumptions of Theorem 3.1 and assuming that (

E(t), H(t)) ∈ H s (τ i ) × H s (τ i ), s ≥ 3/2
, the global error is bounded as:

E -E h 2 0,Ω + H -H h 2 0,Ω 1 2 ≤ C h ν p s min + T h ν-1 p s-3 2 min max t∈[0,T ] E(t), H(t) s,Ω , (48)
where ν = min{s, p min + 1} and p min = min{p i , τ i ∈ Ω h }, p i ≥ 1. The constant C > 0 depends on the material properties and on the shape regularity of the mesh η, but not on p min and h.

We have hence established the semi-discrete result that the error cannot grow faster than linearly in time and that we can control the growth rate by adapting the resolution parameters h and p accordingly. As we shall verify in Sec. 4 this linear growth is a sharp result. However, the numerical experiments will also show that we can expect that the growth rate approaches zero spectrally fast when increasing the approximation order p provided that the solution is sufficiently smooth. Note that the convergence result of Corollary 3.1 is different from the one obtained by Fezoui et al [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the heterogeneous Maxwell equations on unstructured meshes[END_REF]. The convergence result in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the heterogeneous Maxwell equations on unstructured meshes[END_REF] considers only the case of a conforming discontinuous Galerkin formulation where the interpolation degree is constant. The result presented here remains valid on any kind of mesh and discontinuous elements, including hp-type or non-conformal refinement. Now, we give the consistency order of the time-discretized problem. The discretized scheme (6) can be formally seen as the discretization in time of a system of ODEs [START_REF] Sármány | Dispersion and dissipation error in highorder Runge-Kutta discontinuous Galerkin discretisations of the Maxwell equations[END_REF]. The estimation of the consistency error comes directly from Taylor expansions.

If ( E n+1 h , H n+ 3 2 h ) is computed from E n h = E h (t n ) and H n+ 1 2 h = H h (t n+ 1
2 ) by [START_REF] Cohen | A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell's equations in time domain[END_REF] where ( E h (.), H h (.)) denotes the semi-discrete solution of (34), then there exists a constant C independent of ∆t and h, but dependent on the order of the leap-frog scheme N , such that:

E n+1 h -E n h 0,τi + H n+ 3 2 h -H n+ 1 2 h 0,τi ≤ C∆t N +1 N/2+1 m=1 (odd) ( ∂ m t E h 2 0,τi ) 1 2 + N/2+1 m=1 (odd) ( ∂ m t H h 2 0,τi ) 1 2
.

(49)

Provided that the exact solution of the Maxwell system (27) is regular enough, (∂ m t E h , ∂ m t H h ) are some discrete approximation of (∂ m t E, ∂ m t H) which is also a solution of Maxwell equations and, ∂ m t E h 0,τi and ∂ m t H h 0,τi can be bounded independently from h, which proves that the consistency error is of order O(∆t N ).

Convergence of the divergence error

In the absence of sources, it is well known that the electric and the magnetic fields must remain solenoidal1 throughout the computation. Indeed, taking the divergence of Eqs. (27a) and applying Eqs. (27c) in combination with Gauss' law for charge conservation immediately confirms that if the initial conditions satisfy Eqs. (27c), and the fields are evolved according to Maxwell's equations Eqs. (27a), the solution will satisfy Eqs. (27c) at all times. Hence, one can view Eqs. (27c) as a consistency condition on the initial conditions and limit the solution to the time-dependent part of Maxwell's equations, Eqs. (27a). The scheme in Eqs. [START_REF] Cohen | A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell's equations in time domain[END_REF] does not solve Eqs. (1), however, but rather an approximation to it. Hence one needs to consider the question of how well Eqs. ( 6) conserve the divergence. Using the results of Sec. 3.2 we can state the following result. Theorem 3.2 Assume that a solution U = ( E(t), H(t)) ∈ H s (τ i ) × H s (τ i ) with s ≥ 7/2 to Maxwell's equations in Ω h = i τ i exists. Then there exist a constant C dependent on s and the shape regularity of the mesh parameter η, but independent of U, h, and p, such that the divergence of the numerical solution U h to the semi-discrete approximation (34) is bounded as:

∇ • ( E -E h ) 2 0,Ω + ∇ • ( H -H h ) 2 0,Ω 1 2 ≤ C h ν-1 p s-1 min + T h ν-2 p s-7 2 min max t∈[0,T ] E(t), H(t) s,Ω , (50)
where ν = min{s, p min + 1} and p min = min{p i , τ i ∈ Ω h }, p i ≥ 1.

Proof. Consider the local divergence of H on any τ i ∈ Ω h we have:

∇ • ( H -H h ) 2 0,τi ≤ 2 ∇ • ( H -Π p h H) 2 0,τi + 2 ∇ • (Π p h H -H h ) 2 0,τi . (51)
The first term can be bounded using Lemma 3.4 as:

∇ • ( H -Π p h H) 0,τi ≤ C h νi-1 i p s-1 i H s,τi , (52)
where ν i = min{s, p i + 1} and s ≥ 1.

Using the inverse inequality [START_REF] Schwab | p-and hp-finite element methods. Theory and applications to solid and fluid mechanics[END_REF]:

∇ • u h 0,τi ≤ C p 2 i h i u h s,τi , (53)
for all u h ∈ P pi (τ i ), we can bound the second term as:

∇ • (Π p h H -H h ) 0,τi ≤ C p 2 i h i Π p h H -H h 0,τi ≤ CT p 2 i h i h ν-1 i p s-3 2 i ( E, H) s,τi ≤ CT h ν-2 i p s-7 2 i ( E, H) s,τi , (54)
by combining (44) with ( 47). An equivalent bound can be obtained for the divergence of E h in the case of a source free medium which, combined with the above, yields the result.

As could be expected, the result inherits the temporal linear growth from the convergence result and confirms the possibility of recovering spectral convergence of the divergence under the assumption of sufficient smoothness of the solutions. It should be noted that while the result confirms high-order accuracy and convergence, the estimate for the actual convergence rate is certainly suboptimal and leaves room for improvement.

Numerical validation

In the following, we shall discuss the validity of the main theoretical results of the previous sections through the numerical solution of the two-dimensional Maxwell's equations in the TM polarization, i.e. we solve for (H x , H y , E z ). We consider the propagation of an eigenmode which is a standing wave of frequency f = 212 MHz and wavelength λ = 1.4 m in a unitary metallic cavity with = µ = 1 in normalized units. Owing to the existence of an exact analytical solution, this problem allows us to appreciate the numerical results at any point and time in the cavity. Numerical simulations make use of a non-conforming locally refined triangular meshes of the square [0, 1] × [0, 1] as shown on Fig. As a first verification of the theoretical estimates, we consider a non-conforming mesh consists of 782 triangles and 442 nodes (36 of them are hanging nodes). All simulations are carried out for t = 150 which corresponds to 106 periods. We plot on Fig. 5 the time evolution of the L 2 error of the DGTD-P p and DGTD-P pc :P p f methods using the LF 2 and LF 4 schemes. Tab. 3 gives the final L 2 error, the number of degrees of freedom and the CPU time to reach t = 150. It can be observed from Fig. 5 that the gain in the L 2 error is notable when the accuracy in space and time are increased. Moreover, it is clear from [START_REF] Lu | Discontinuous Galerkin methods for dispersive and lossy Maxwell's equations and PML boundary conditions[END_REF] and Lemma 3.2 that, for a given mesh, each time step of LF 4 scheme requires 2

INRIA

On Fig. 9 we show the numerical h-and p-convergence of the divergence of H using the LF 2 and LF 4 schemes. Consistent with the theoretical result in Theorem 3.2, the divergence error vanishes spectrally as we increase the approximation order p. Corresponding asymptotic convergence orders of the divergence of H are given in Tab. 7. One can observe that the convergence order is bounded by N + 2 contrary to what we have observed for the h-convergence of the DGTD methods which confirms that the estimate (50) is suboptimal and leaves room for improvement.

Figure 2 :

 2 Figure 2: Temporal allocation of electric (×) and magnetic ( ) fields in the LF scheme.

Figure 4 :

 4 Figure 4: Example of a non-conforming locally refined triangular mesh.

Figure 6 :

 6 Figure 6: h-convergence of the DGTD-P p (top) and DGTD-P pc :P p f (bottom) methods.L 2 error as a function of the square root of #DOF.

  scheme h=1/2 h=1/3 h=1/4

Figure 7 :

 7 Figure 7: p-convergence of the DGTD-P p (top) and DGTD-P pc :P p f (bottom) methods.L 2 error as a function of the approximation order p.

Figure 8 :

 8 Figure 8: Global L 2 error of the divergence of H as a function of time and p. DGTD-P p (top) and DGTD-P pc :P p f (bottom) methods using LF 2 (left) and LF 4 (right) schemes.

10

 10

Figure 9 :

 9 Figure 9: h-and p-convergence of the divergence of H. DGTD-P p (top) and DGTD-P pc :P p f (bottom) methods.

Table 1 :

 1 The values of ν N and r N for several LF N schemes.

	N 2	4	6	8	10	12	14	16	18	20
	ν N 1 2.847 3.681 3.793 5.272 4.437 6.422 7.534 7.265 8.909
	r N 1 1.424 1.227 0.948 1.05 0.739 0.917 0.942 0.807 0.891

  4. In the sequel, we compare the LF 2 and LF 4 schemes using the DGTD-P p and DGTD-P pc :P p f methods previously studied in[START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplicial meshes for electromagnetic wave propagation[END_REF]-[START_REF] Fahs | Preliminary investigation of a non-conforming discontinuous Galerkin method for solving the time domain Maxwell equations[END_REF]-[START_REF] Fahs | A hp-like discontinuous Galerkin method for solving the 2D time-domain Maxwell's equations on non-conforming locally refined triangular meshes[END_REF]. In Tab. 2, we summarize the CFL values of the LF 2 based DGTD-P p and DGTD-P pc :P p f methods. The CFL values of the LF 4 schemes are given by CFL(LF 4 ) = 2.847 CFL(LF 2 ). If p c = p f , the DGTD-P pc :P p f method has the same stability limit as the DGTD-P min{pc,p f } method, as long as the mesh is actually refined.

	1 1
	0.8 0.8
	0.6 0.6
	0.4 0.4
	0.2 0.2
	0 0
	0 0	0.2 0.2	0.4 0.4	0.6 0.6	0.8 0.8	1 1

Table 3 :

 3 # DOF, L 2 -error and CPU time using the LF 2 and LF 4 based DGTD methods.

	DGTD-P p method		LF 2		LF 4
	p	# DOF	Error	CPU (min)	Error	CPU (min)
	2	4692	1.8E-03	11	5.5E-04	8
	3	7820	3.1E-04	39	2.4E-05	28
	4	11730	1.9E-04	98	1.5E-05	70
	5	16422	1.5E-04	220	1.3E-05	155
	DGTD-P pc :P p f method		LF 2		LF 4
	p c :p f	# DOF	Error	CPU (min)	Error	CPU (min)
	3:2	6668	1.3E-03	17	2.3E-05	12
	4:2	9138	1.3E-03	27	1.5E-05	19
	4:3	10290	3.2E-04	61	1.5E-05	44
	5:4	14694	2.0E-04	134	1.4E-05	95

Table 4 :

 4 Asymptotic convergence orders of the LF 2 and LF 4 based DGTD methods.

	DGTD-P p method, p =	0	1	2	3	4	5	6
	LF 2 scheme	1.06 1.19 2.18 2.37 2.29 2.25 2.26
	LF 4 scheme	1.06 1.14 2.23 3.03 4.30 4.50 4.50
	DGTD-P pc :P p f method, p c :p f = 1:0	2:1 3:2 4:3 5:4 6:5
	LF 2 scheme	1.30 2.23 2.08 2.27 2.13 2.17
	LF 4 scheme	1.05 2.20 3.01 4.21 4.50 4.48

Table 5 :

 5 p-convergence of the DGTD-P p method. LF 4 1 3.0E-01 2.6E-01 1.8E-01 1.7E-01 1.1E-01 9.8E-02 2 4.9E-02 4.5E-02 1.8E-02 2.2E-02 9.2E-03 9.0E-03 3 9.1E-03 8.9E-03 2.8E-03 2.4E-03 1.4E-03 8.9E-04 4 3.2E-02 1.0E-03 1.3E-03 2.0E-04 7.4E-04 6.9E-05 5 1.6E-02 1.7E-04 7.4E-04 2.2E-05 4.1E-04 5.1E-06 6 9.3E-04 2.0E-05 4.1E-04 2.1E-06 2.3E-04 7.3E-07 7 1.7E-04 2.5E-06 4.4E-05 3.3E-07 2.7E-05 1.1E-07

				DGTD-P p method
		h = 1/2	h = 1/3	h = 1/4
	p	LF 2	LF 4	LF 2	LF 4	LF 2

Table 6 :

 6 p-convergence of the DGTD-P pc :P p f method.

				DGTD-P p method
		h = 1/2	h = 1/3	h = 1/4
	p c :p f	LF 2	LF 4	LF 2	LF 4	LF 2	LF 4
	1:0	8.3E-01 8.5E-01 1.3E-01 1.9E-01 7.4E-02 7.6E-02
	2:1	1.9E-01 1.8E-01 3.9E-02 3.9E-02 2.1E-02 1.3E-02
	3:2	1.8E-02 1.8E-02 3.2E-03 2.7E-03 1.4E-03 8.5E-04
	4:3	2.3E-03 2.3E-03 5.6E-04 2.5E-04 2.9E-04 7.6E-05
	5:4	8.2E-04 2.0E-04 3.2E-04 1.7E-05 1.8E-04 4.4E-06
	6:5	4.1E-04 2.3E-05 1.8E-04 1.9E-06 1.0E-04 4.7E-07
	7:6	2.3E-04 3.4E-06 1.0E-04 2.8E-07 5.7E-05	-
	8:7	1.0E-04 6.1E-07 4.6E-05	-	2.1E-05	-
	9:8	4.9E-05	-	1.9E-05	-	7.9E-06	-

Table 7 :

 7 Asymptotic convergence orders of the divergence of H. LF 2 scheme 2.33 2.81 3.84 3.24 3.46 LF 4 scheme 2.26 2.73 3.94 4.40 5.50

	DGTD-P p method, p =	1	2	3	4	5	6
	LF 2 scheme	0.89 2.10 2.94 4.07 3.49 3.45
	LF 4 scheme	0.97 2.05 3.00 4.09 4.58 5.66
	DGTD-P

pc :P p f method, p c :p f = 2:1 3:2 4:3 5:4 6:5

RR n 6699

A solenoidal vector is a vector field v with zero divergence, ∇ • v = 0.INRIA

Fahs & Lanteri

times more memory than the LF 2 , but its stability limit is almost 2.85 times less restrictive. Then, the LF 4 schemes requires almost 1.5 times less CPU time and is roughly 15 times more accurate than the LF 2 scheme based on the observed L 2 errors. Furthermore, for a given accuracy, the LF 4 based DGTD-P pc :P p f method requires less CPU time than the LF 4 based DGTD-P p method. Fig. 6 illustrates the numerical h-convergence of the DGTD-P p and DGTD-P pc :P p f methods. Corresponding asymptotic convergence orders are summarized in Tab. 4. As it could be expected from the use of a N th accurate time integration scheme, the asymptotic convergence order is bounded by N independently of the approximation order p. On Fig. 7 we show the numerical p-convergence of the DGTD-P p and DGTD-P pc :P p f methods for different approximation orders p and different mesh resolutions h. Corresponding L 2 errors are given in Tab. 5 and 6. Following the main result, Theorem 3.1, we expect that the error grows at most linearly in time and that the growth rate should vanish spectrally for smooth solution. The results on Fig. 7 and in Tab. 5 and 6 not only confirm the validity of both statements but also illustrate that Theorem 3.1 is sharp, i.e. we cannot in general guarantee slower than linear growth, although we can control the growth rate by the approximation order p. We conclude this experimental study by considering the numerical behavior of the divergence error. For this purpose, we still consider the eigenmode problem. The computational domain is discretized by a non-conforming locally refined mesh with 48 triangles (32 of them in the refined region) and 37 nodes (16 of them are hanging nodes), which corresponds to a grid resolution of 5 points per wavelength. Simulations are carried out for t = 30 which corresponds to 20 periods. Fig. 8 shows the global L 2 error of the divergence of H as a function of time and the approximation order p using respectively the DGTD-P p and DGTD-P pc :P p f methods. One can observe that for p ≤ N + 2, the error vanishes faster than for p > N + 2, N being the order of the leap-frog (LF N ) scheme.

INRIA 5 Concluding remarks

The main purpose of this report has been to study both theoretically and numerically an arbitrarily high-order DGTD method for the discretization of the time-domain Maxwell equations on nonconforming simplicial meshes. The central element which distinguishes the current work from previous attempts to develop such DGTD methods is that a high-order leap-frog time integration scheme is adopted here instead of a high-order Runge-Kutta method. We have proved that the resulting DGTD method conserves a discrete equivalent of the electromagnetic energy and is stable under some CFLtype condition. Also, we have developed a complete, if not optimal, convergence theory. We have confirmed the results of the analysis thorough numerical experiments using an eigenmode problem in two space dimensions, illustrating the flexibility, versatility, and efficiency of the proposed arbitrarily high-order DGTD method.


Shreya Nandi

Jiangxiang Xiang Wu

Pardis Simon

Nicolas Nuns

Martine Trentesaux

Asma Tougerti

Emiliano Fonda

Jean-Sébastien Girardon

Jean-François Paul

Anne-Sophie Mamede
email: anne-sophie.mamede@ensc-lille.fr

Elise Berrier
email: elise.berrier@univ-lille.fr

Mn-or Cu-substituted LaFeO 3

published or not. The documents may come

Introduction

In response to rising concerns on air quality, gasoline-based vehicles are nowadays equipped with efficient catalytic converters, which ensure the concurrent abatement of chief pollutants such as carbon monoxide (CO), nitrogen oxides (NO x ) and unburnt hydrocarbons (HC). Within the catalytic converter, a three-way catalyst (TWC) plays the simultaneous role of oxidizing the CO and HC and reducing the NO x into molecular nitrogen via a redox process. [1] The interest towards LaFeO 3 (LFO) perovskite-based materials for TWC application rose immensely with the discovery of their so-called "self-regenerating properties" by Nishihata et al. [2]. By means of an in situ investigation of a Pd-doped perovskite, LaFe 0.57 Co 0.38 Pd 0.05 O 3 , using X-ray diffraction and X-ray absorption spectroscopy, the authors have postulated that the palladium oscillates between lattice position under oxidative conditions and out of the lattice in the form of Pd 0 nanoparticles under reducing conditions. This particular behaviour was proposed to be responsible for the standing catalytic activity of the perovskite-based catalysts in comparison with their alumina-supported counterparts, by preventing the sintering of noble metal particles and their subsequent deactivation [3]. Since then, this so-called self-regeneration phenomenon has been extensively investigated [4] and exploited [5,6], revealing overall benefits of the use of perovskite-based oxide formulations for TWC applications while drastically lowering the content of platinum-group metals [7,8].

Previous studies on LaFeO 3 -based TWC by Schoen et al. have evidenced that purposely prepared Ladeficient formulations, La 1-x FeO 3 (0.1< x <0.33), show better catalytic activity than the stoichiometric LaFeO 3 one, not only in CO oxidation but also in NO reduction [9]. The authors have correlated this catalytic enhancement with the restoration of a Fe/La ratio (XPS) close to unity at the surface of the material while LFO typically features a Fe/La atomic ratio of 0.4 to 0.6 at the surface [9,10]. This trend is intrinsic to the perovskites prepared using citrate complexation route in combination with mild annealing temperatures (≤ 600 °C), typically preferred over solid-state reaction, for preparing highest possible specific surface area catalysts while keeping the perovskite structure. As a result of this compromise between structural and textural properties, the as-obtained materials show the expected bulk composition and perovskite structure, but feature a clear lanthanum enrichment at the surface. The latter has been assigned to the high reactivity of lanthanum oxides with air to form surface lanthanum hydroxides and carbonates [9][10][11][12].

Further, substitution at A-and/or B-site of the perovskite can effectively tune the catalytic properties by varying the valency of the embedded cations and introducing lattice defects and ionic vacancies [13].

In parallel, the co-existence of two different cations at B-site can bring synergetic effects and hence, can enhance the overall reducibility [14]. It was shown recently by Wu et al. that a 20% substitution of Fe by Cu at B-site of a La-deficient La 0.67 FeO 3 , further denoted as dLFO, significantly lowers the temperature required for reaching 50% of CO conversion, θ 50 (CO), under TWC-relevant conditions1 [15]. Indeed, while the θ 50 (CO) measured on dLFO was 385 °C, the θ 50 (CO) for Cu-doped dLFO with La 0.67 Fe 0.8 Cu 0.2 O 3 formulation was 260 °C. In general, the same trend is reported for Cu-doped stoichiometric LaFeO 3 , which all show enhanced TWC activity in comparison with the parent LaFeO 3 material [7,16,17]. Similarly, although less spectacular, 20% substitution of Fe 3+ by Mn 3+ yielding La 0.67 Fe 0.8 Mn 0.2 O 3 stoichiometry also led to a noticeable lowering of the θ 50 (CO) by 45 °C under similar conditions. In the case of Cu-doped LaFeO 3 , such beneficial effect has been proposed to originate from a strong interaction between copper and iron-related species [16]. From a structural point of view, it is well established that at ambient pressure, La-Cu mixed oxides rather crystallize in the Ruddlesden-Popper phase (La 2 CuO 4 ) instead of the simple ABO 3 perovskite phase [18][19][20][21][22]. This makes the ocurrence of LaCuO 3 phase unlikely in the present case. Moreover, Cu-substitution has been widely carried out for yielding a series of AB 1-x Cu x O 3 perovskite-based oxides, which can lead to the formation of the segregated CuO phase with a higher Cu content, typically for x>0.4 [23][24][25][26][27][28][29]. On the other hand, Mn-substitution in ferrite perovskites, LaFe 1-x Mn x O 3 with x ranging between 0 and 1, mainly results into a solid solution owing to close crystal parameters of LaFeO 3 and LaMnO 3 perovskites [30,31].

In this paper, we present a systematic comparative study of Cu or Mn substituted La 0.67 FeO 3 catalysts using in situ/operando spectroscopies to understand the beneficial role of each substituted B-site element under TWC-relevant model reaction. Particularly, we have used operando Raman spectroscopy to follow any structure-activity variations in the catalysts in real-time during CO-TPR. Moreover, quasiin situ XPS was drawn upon to detect the chemical changes after a sequential exposure to CO and NO, while the structure and surface composition was refined using Fe K-edge XANES and ToF-SIMS analysis, respectively.

Materials and Methods

Sample preparation. The sample preparation was detailed elsewhere [15]. In brief, powders with La 0.67 Fe 0.8 Mn 0.2 O 3 and La 0.67 Fe 0.8 Cu 0.2 O 3 formulation -further denoted as Mn-dLFO and Cu-dLFO for conciseness -were prepared using the citrate complexation method. Therefore, weighted amounts of nitrate salts of La, Fe, Cu and Mn were dissolved in de-ionised water and allowed to react with citric acid at 60 °C under reflux before being dried at 80 °C, decomposed at 150 °C and further calcined in air at 600 °C overnight.

Operando Raman analyses. The Raman spectra were collected in confocal mode using the 488 nm line of a Ar + laser (Melles Griot), focused on the sample placed in an atmospheric chamber (Harrick Scientific) equipped with a dedicated flat silica dome by means of a 50X long working distance objective (NA: 0.6, Olympus). The scattered light was guided through a 150 µm pinhole, dispersed and collected using a Peltier-cooled CCD implemented within a HORIBA HR equipment. The laser power was reduced using neutral density filters in order to avoid sample alteration. The gas feed composition was tuned using mass flow controllers (Brooks), while the gaseous products were analysed by means of a three-columns micro-GC (Agilent). Prior to the CO-TPR, both Mn-dLFO and Cu-dLFO powders were stabilised and dehydrated by means of a thermal treatment under dry air up to 450 °C with a ramp of 5 °C per minute and a dwell time of 1 hour.

Quasi-in situ XPS analyses were performed using a Kratos Analytical Axis UltraDLD spectrometer.

The excitation was ensured by a monochromatic aluminum Kα source at 1486.6 eV , while a pass energy of 20 eV for an analysed surface area of 700 µm x 300 µm was used. Charge compensation was applied to neutralize any charging effects. The measurements were done at a pressure of 10 -9 mbar in the analysis chamber. All spectra were calibrated against the adventitious C 1s spectrum at 285 eV. Quasi in situ measurements were made using the catalysis cell, where the samples were thermally pre-treated under gaseous reactant flow at atmospheric pressure and then transferred to the analysis chamber under ultra-high vacuum (UHV). Both powders were stabilised by a thermal treatment up to 500 °C for 1 hour under O 2 before the CO (5% in He) or NO (2% in He) treatment.

XRD. The X-ray diffraction structural characterization of the powders were carried out using a Bruker AXS D8 Advance diffractometer using CuK α (1.5418 Å, 40kV, 30mA) radiation operated in the Bragg-Brentano reflection geometry.

ToF-SIMS. The samples were analysed on a ToF.SIMS 5 spectrometer (ION TOF GmbH, Germany) equipped with a pulsed bismuth liquid metal ion gun with the primary ion cluster of Bi + 3 (25 keV, 1 pA current) as heavier ions can promote formation of heavier fragments. The analysed area was 500 µm x 500 µm for spectral analyses while maintaining static conditions.

XANES. measurements were made at SAMBA beamline in SOLEIL synchrotron facility (France) [START_REF] Briois | UVX 2010 -10e Colloque sur les Sources Cohérentes et Incohérentes UV[END_REF]. The beam size at the sample was 200 µm x 300 µm with a focusing monochromator between two bendable cylindrical mirrors with 5-43 keV energy range. Spectra were collected at the Fe K-edge in the transmission mode. Reference powder compounds were recorded accordingly after proper dilution in boron nitride.

Results and Discussion

Ex situ physico-chemical properties

X-Ray diffraction (XRD).

The XRD pattern of Cu-dLFO powder is presented in Fig. 1a). The diffractogram shows the characteristic reflections of the orthorhombic LaFeO 3 perovskite with pnma space group (JCPDS PDF 37-1493).

The lattice constants are a=5.545 Å, b=7.857 Å, and c=5.550 Å. In addition, patterns of weak intensity are visible at 35.6 • and 38.7 • , which are assigned to the (-111) and (111) reflections of a segregated CuO phase, respectively.

The diffraction pattern of Mn-dLFO shown in Fig. 1b) displays diffraction peaks at 32.2 • , 40

X-Ray Photoelectron spectroscopy (XPS).

Both Fe 2p and Cu 2p photopeaks of Cu-dLFO are shown in Fig. 2. The Fe 2p signal at the binding en- ergy of 710.2 eV displays a shake-up satellite feature at 718.8 eV, typical of Fe 3+ ions [33]. On the other hand, the shape of the Cu 2p photopeak envelope mostly refers to Cu 2+ with a binding energy of 934.0 eV together with an intense and broad satellite around 943 eV [34]. From a quantitative point of view, a La/(Fe+Cu) atomic ratio of 0.68 is estimated, which is in line with the theoretical bulk composition.

The atomic ratio of the B-site cations, Cu/Fe, is of 0.33 instead of 0.25 as per the La 0.67 Fe 0.8 Cu 0.2 O 3 formula, indicating a slight surface over-representation of the copper content. A Fe/La ratio of 1.1 is estimated, which is identical to the one previously reported for La 0.67 FeO 3 (dLFO) [9].

The Fe 2p and Mn 2p photopeaks of the as-prepared Mn-dLFO are presented in Fig. 3. The Fe 2p

Figure 3: XPS Fe 2p and Mn 2p spectra of Mn-substituted perovskite powder, Mn-dLFO.

signal is present at a binding energy of 710.4 eV also showing the characteristic Fe 3+ satellite feature at 719 eV. The Mn 2p photopeak is detected at a binding energy of 641.8 eV, which is characteristic of oxidized manganese. Since the quality of the Mn 3s photopeak was too weak to allow in-depth analysis, the Mn 2p orbital has been used for refining the oxidation state. The Mn 2p spectrum has been fitted by relying on the multiplet splitting parameters given by Biesinger et al. as result of their refining of best parameters from various standard references of Mn oxides with fixed peak positions, relative areas and peak widths [33]. This systematic study makes it possible to quantify the average oxidation state(s) of manganese with reasonable reliability when Mn 3s is not available, since it provides robust multiplet splitting parameters to assist the peak fitting of Mn 2p orbital. Hence, the peak decomposition for the Mn 2p photopeak shown in Fig. 3 indicates that 100% of manganese stands as Mn 3+ . Thereby, the as-received sample of Mn-dLFO perovskite contains both B-site cations in +III oxidation state. A La/(Fe+Mn) atomic ratio of 0.7 has been obtained, which stands in line with the nominal La 0.67 Fe 0.8 Mn 0.2 O 3 formula, as is the case for the as obtained B-site ionic ratio, Mn/Fe, of 0.26.

Raman analysis.

The ex-situ Raman spectrum of the Cu-dLFO catalyst is shown in Fig. 4, together with that of the parent LaFeO 3 (LFO) perovksite for structural comparison. All the Raman lines characteristic of the modes of LaFeO 3 previously discussed in literature are detected at the expected wavenumbers in the Raman spectrum of Cu-dLFO [35,36]. In brief, the La-centered vibrations, denoted as A modes are observed below 170 cm -1 , the FeO 6 octahedra tilt modes (T) are visible between 260 and 290 cm -1 while the bending modes (B) are responsible for the Raman lines around 430 cm -1 . The Fe-O stretching modes (S) are observed around 500 (anti-symmetric) and 625 cm -1 (symmetric). Noticeably, the scattering in the region of T and B modes is not as well resolved as in the parent LFO material, while the S modes region feature a clearly sharper line of larger intensity. In general, the Raman signature of CuO is composed of two sharp lines observed at 290-300 cm -1 (strong, Ag symmetry) and 340-350 cm -1 (weak, Bg) together with a broader component at ∼ 625-636 cm -1 (medium, Bg) [37,38], while Cu(OH) 2

features Raman lines at ∼292 and ∼488 cm -1 [38][39][40]. In addition, the mixed-valence phase Cu 4 O 3 is generally identified from its broad Raman lines at 320, 540 and 650 cm -1 [41]. In the present case, although CuO was detected by XRD, the Raman spectrum of Cu-dLFO showed no clear evidence of any copper oxide species and the poor signal quality has discouraged us to attempt a more detailed analysis.

Note, however, that if any, the Raman lines of the copper oxide and related compounds would largely overlap with those of the LFO in the present case, making the technique possibly blind to such a minor phase.

The Raman spectrum of Mn-dLFO is shown in Fig. 5. On the very first look, this Raman spectrum appears completely different from the parent LFO perovskite structure. Indeed, the overall spectral intensity and the signal-to-noise ratio are significantly higher as compared to the spectra shown in Fig. 4.

Moreover, the A, T, B and S modes of LFO discussed earlier are not clearly visible, rather the Raman spectrum of Mn-dLFO is dominated by two intense peaks at ∼500 and ∼650 cm -1 together with a weak band at around 300 cm -1 . Actually, this spectrum is more alike the signature of the orthorhombic LaMnO 3 perovskite [42,43] showing modes at 493 and 611 cm -1 related to oxygen motions. The intensity of those modes was proposed to be enhanced by the Jahn-Teller (JT) distortion [43] induced by Mn 3+ cations having a d 4 electronic configuration in high spin state (E 1 g ) [44]. In comparison with LaMnO 3 , all three modes shown in Fig. 5 feature a noticeable frequency shift towards higher frequencies, as was already reported for LaFe 1-x Mn x O 3 solid solutions with perovskite structure [45,46].

Hence, the Raman signature of Mn-dLFO suggests that a large share of Mn 3+ cations are inserted at B site of a La 0.67 Fe 0.8 Mn 0.2 O 3 solid solution.

Figure 5: Raman spectrum of Mn-substituted perovskite powder, Mn-dLFO (λ exc. = 488 nm).

X-Ray absorption spectroscopy.

The Fe K-edge XANES spectra of Cu-dLFO and Mn-dLFO are shown as solid lines in Fig. 6a) and b), respectively. The XANES fingerprints of Cu-dLFO and Mn-dLFO powders show absorption edges of 7128.1 eV and 7127.6 eV, respectively, signifying the presence of Fe 3+ ions, as expected for the LFO perovskite [47]. This result is in full agreement with the XPS results, indicating that iron stands as Fe 3+

cations at both the surface and the bulk of the catalysts. A structural refinement using LaFeO 3 and dLFO as reference materials was made, and the XANES spectrum of the most relevant reference was plotted as dotted lines together with the one of the catalysts in Fig. 6. Hence, the Fe-K edge XANES spectrum of Cu-dLFO is found to be similar to the one of LFO, while the Mn-dLFO is close to that of dLFO.

Looking more in detail, the spectrum of Cu-dLFO actually shows a minor decrease of the white line intensity as compared to LaFeO 3 . A linear least-square fitting (LCF) analysis of Cu-dLFO, given in Fig. 6c) performed using reference spectra of iron oxides (Fig. S2 in supporting information) suggested that Cu-dLFO may contain some Fe 3+ embedded in a minority α-Fe 2 O 3 phase besides the majority LaFeO 3 one. On the other hand, the similarity of the XANES spectrum of Mn-dLFO to the one of the La-deficient dLFO material (Fig. 6b)) indicates that a more significant share of Fe 3+ ions are involved in the α-Fe 2 O 3 phase, as was previously evidenced by our group for the La 0.67 FeO 3 (dLFO) catalyst using Raman spectroscopy [48]. This is further confirmed by the LCF results shown in Fig. 6d), where no less than 33% of the Fe 3+ ions stands as α-Fe 2 O 3 and rest 67% as LaFeO 3 . The higher content of iron oxide in Mn-dLFO compared to the one of Cu-dLFO suggests that fewer Fe 3+ are embedded within the perovskite lattice after Mn substitution.

To sum up the results discussed above, Cu-dLFO powder roughly consists of a LaFeO 3 supported CuO catalyst. Our XANES results suggests a possible minority hematite phase segregated on the surface, while we could not completely rule out the insertion of copper at B site of the perovskite as trace amount. On the other hand, Mn-dLFO is more alike the La-deficient material, dLFO, which embeds a mixture of orthorhombic perovskite together with a segregated iron oxide phase. Moreover, our results also suggests that most Mn 3+ actually substitute Fe 3+ at B site of the perovskite. This results in a large exsolution of Fe 3+ as Fe 2 O 3 . In order to support this assumption, we have tracked the presence of mixed oxide fragments in the mass spectra of both catalysts using ToF-SIMS analysis.

ToF-SIMS detection of surface mixed oxide fragments. The detection of the mixed oxide ions produced by ToF-SIMS is a direct indication of the presence of a mixed oxide phase and unlikely from ionic recombination [49], which was further confirmed by us in the present case. The most representative mixed oxide fragments, namely La 2 FeMO + 5 and La 3 FeMO + 6 , where M = Mn or Cu, are shown in Fig. 7 for both catalysts. In the case of Mn-dLFO, such ionic fragments are clearly detected, supporting the hypothesis of a single LaFe 1-x Mn x O 3±δ phase at the surface. In contrast, the intensity of similar mixed oxide ionic fragments is negligible in the case of Cu-dLFO. Similarly, their isotopes with 65 Cu and/or 54 Fe are completely absent in the mass spectrum.

This highlights the fact that, if any, only a tiny share of Cu 2+ ions stands at the B-site of the perovskite.

Starting from the assumptions made above, a simple reasoning made by linking both the qualitative (from XRD and Raman) and quantitative (LCF of XANES) assessments with the nominal bulk compositions gives the structural picture of both catalysts shown in Table 1. Concerning the surface compo- Operando Raman monitoring of CO -temperature programmed reaction (TPR)

Cu-dLFO

The operando Raman spectra recorded under a flow of 5% CO in He during the TPR are shown in Fig.

8.

• Between 25 and 100 °C, minor changes were detected in the Raman spectrum of Cu-dLFO, which include the transient detection of a sharp signal around 2135 cm -1 assigned to physisorbed CO together with an increase of the relative intensity of the S modes region.

• Between 100 and 250 °C, a broad line assigned to poly-aromatic hydrocarbons (PAHs) gradually emerges at ∼1595 cm -1 . Meanwhile, the T and S modes of the perovskite decrease in intensity whereas the intensity in the region of A and B modes remains stable.

• Between 350 and 400 °C, the PAHs were detected yet again -after having vanished at 300 °Cwhile the intensity of T and S modes of the perovskite further decreased. Interestingly, the intensity of scattering between 100 and 200 cm -1 (A modes) was kept quite high. A similar deposition of PAHs was observed in the copper-free LaFeO 3 and La 0.67 FeO 3 counterparts under similar conditions [48].

The corresponding CO consumption measured on-line during this experiment is shown in Fig. 9. During TPR, the CO consumption was accompanied by the formation of CO 2 with an associated carbon balance kept between 96.5 and 100%. As shown in Fig. 9 (left), the consumption of CO to CO 2 , marginal below 100 °C, sharply increased up to 0.105 mmol.min -1 g -1 cat at 170 °C. The CO consumption was found to rapidly drop during the plateau at 200 °C. This trend is not surprising in absence of a co-reactant, such as NO, to participate in the re-oxidation and regeneration of the active phase after reduction. Similarly, still less spectacular spurt-drop sequences were further detected between 250 and 350 °C, with another noticeable increase of CO consumption rate up to 0.091 mmol.min -1 g -1 cat at 400 °C. Fig. 9 (right) displays the CO consumption under isothermal conditions at 100, 200, 300 and 400 °C plotted against the same time range for allowing a direct visualisation of the deactivation kinetics. At 100, 200 and 300 °C, the CO consumption rate stabilised rapidly to 24.4, 31.3 and 35.3 µmol.min -1 .g -1 cat. , respectively. On the other hand, the CO consumption rate was still not stabilised after 1 hour reaction at 400 °C. More precisely, the CO consumption rate follows an exponential decay over time with a constant λ=0.28±0.04

min -1 at 200 °C, while it is almost linearly decreasing at 400 °C (λ=0.028±0.001 min -1 ). Keeping in mind that the regeneration of the active sites cannot be mediated by the gas phase and solely relies on the solid phase, this trend suggests that at mild temperatures (100-200 °C), a limited number of highly active surface oxygen species sustains the decent catalytic activity observed for only a short period of time. On the other hand, a more sizeable pool of less active oxygen moieties are engaged at 400 °C.

These observations are in agreement with the commonly agreed description of adsorbed surface oxygen moieties, which drive oxidation reactions at milder temperatures, while bulk oxygens are involved at higher temperatures on migration to the surface [13,50].

Besides the formation of PAHs and the subsequent increase of the broad lines at 1360 and 1595 cm -1 , the return spectrum recorded at 100 °C under He/CO after the CO-TPR is quite similar to the one recorded before under the same conditions with the exception of a slight decrease of the scattered intensity around 500 cm -1 (Fig. S1 in supporting information). This indicates that the changes observed in 100-800 cm -1

range by the operando Raman spectroscopy are mostly related to temperature-driven reversible processes.

Our operando Raman investigation has allowed to get only fragmentary data on the structural evolution of the active phase of Cu-dLFO during the model reaction. Hence, while Cu-dLFO was reaching a maximal CO consumption rate between 100 and 200 °C, the Raman signature of the catalyst barely evolved (see b and c spectra in Fig. 8), confirming that the technique was not capable of monitoring the structural evolution of the active phase. However, interestingly the first detection of PAHs at 200 °C coincides with the sharp increment of the CO consumption rate. Two regimes were detected for CO consumption, one fast process involving surface oxygen species present at mild temperatures (150-200 °C), and another one which requires higher temperatures. This higher temperature regime may indicate that oxygen diffusion from the sub-surface to the surface has become the rate-limiting factor.

Mn-dLFO. structure. This trend is further confirmed by the comparison of Raman spectra recorded at 100 °C under CO at beginning and end of CO-TPR shown in Fig. 11. Moreover, some PAHs are also detected around 1595 cm -1 at 350 °C, while already detected at 200 °C in the case of Cu-dLFO. Although evidenced by XANES, the fingerprints of hematite (212, 280 and ∼425 cm -1 ) were never detected along CO-TPR, as was, in fact, already the case for La 0.67 FeO 3 [48]. The apparent recovering of the LaFeO 3 modes suggests that the J-T enhancement typical of the LaFe 1-x Mn x O 3 phase is not at play any more. This could point towards a change of Mn 3+ cations from d 4 electronic configuration featuring strong JT effect towards a JT-inactive one under reducing conditions.

The CO consumption rate recorded on-line during CO-TPR is shown in Fig. 12. In general, the CO consumption rate of Mn-dLFO is lower by at least one half as compared to Cu-dLFO. The maximal CO consumption rate per unit mass of catalyst is observed at 350 °C with a carbon balance maintained above 97% along the experiment. As for Cu-dLFO, the CO consumption rate rapidly decreased over time at each temperature plateau, but the pool of both surface and diffusible bulk oxygen species is significantly smaller as compared to Cu-dLFO. Quasi-in situ XPS investigation of CO-NO sequence.

By using the quasi-in situ mode, it was possible to detect the chemical changes occurring under CO, especially for the Cu-related species, which was not accessible by operando Raman spectroscopy. The perovskites were initially thermally treated under 20% O 2 /He for stabilization prior to the CO exposure.

After the CO treatment, the catalysts were also subjected to NO to observe the redox behaviour in each catalyst. Here, the La 3d signal around 834 eV representing the La 3+ state remains unchanged during the redox treatment, as it is not expected to show any redox activity and therefore is not discussed here.

Cu-dLFO. Catalyst reduction under 5% CO/He (500 °C). The XPS spectra of the B-site elements after CO exposure are shown in Fig. 13 (black line). The binding energy of Fe 2p, the intensity of the Fe 3+ satellite as well as the FWHM values remains unchanged after the CO treatment, indicating that iron do not majorly participate in the redox process at play in Cu-dLFO. In contrast, the Cu 2p spectrum signal [34]. From our experimental data, α is estimated to be 1851.2 eV corresponding to metallic copper (Cu 0 ). From a quantitative point of view, the A/B ratio increased severely after CO exposure from 0.68 to 1.14 (Table 2). This trend is mostly due to the drastic decrease (-78%) of the Cu/La atomic ratio at the surface, while the Fe/La ratio was reduced to a lesser extent by 32%. As a result, the Cu/Fe atomic ratio changed remarkably from 0.34 to 0.11. This evolution accompanies the reduction of Cu(II)

to metallic copper, which can denote the reduction of formerly dispersed CuO moieties into larger Cu 0 nanoparticles.

Catalyst re-oxidation under 2% NO/He (500 °C). Upon reaction with NO, the Cu 2+ satellite appears in the Cu 2p signal again, along with a clear shoulder of the Cu 2p 3/2 at 933.9 eV. At this point, the copper is present in two oxidation states, made evident from the presence of the photopeak at 932.1

eV. The calculated value of the α parameter was 1849.4 eV, corresponding to the one expected from

Cu + [34]. Thereby, the Cu 0 formed under the CO step have partly re-oxidized to a mixture of CuO and Cu 2 O after the NO exposure. The A/B site and Cu/Fe atomic ratios are only partially shifted back from 1.14 to 0.82 and 0.11 to 0.17, respectively, without recovering their initial levels before the CO or NO treatment. More precisely, while the Fe/La atomic ratio is almost recovered, the Cu/La one is only half restored.

Mn-dLFO. The quasi-in situ Fe 2p and Mn 2p photopeaks under each treatment step are shown in Fig. 14.

After activation under oxygen, the surface chemistry did not change for Fe or Mn as compared to the as-received state. Both the Mn and Fe ions are present in the +III valence state. from its satellite structure is not so obvious due to its overlapping with the predominant Fe 3+ fingerprint. The spectral decomposition reveals less than 5-10% of Fe 2+ component falling in the same range of uncertainty as the spectral decomposition, while the remaining exists as +III, leading to an average oxidation state of 2.7±0.1 of iron. The Mn 2p orbital has been decomposed according to the Biesinger parameters [33] as mentioned previously, and the corresponding fitting results (Fig. S3 in supporting information) showed the presence of 67% Mn 2+ besides 33% Mn 3+ (average oxidation state=2.3). From a quantitative point of view, Mn 3+ cations show more extensive reduction at the catalyst's surface than Fe 3+ , when reduced under CO. On the basis of our XANES results stating that Mn-dLFO embeds substantial amounts of α-Fe 2 O 3 , and in comparison with the absence of Fe reduction under same conditions in Cu-dLFO, containing only limited hematite content, we assume that the reduction of Fe observed here mostly occurs at the α-Fe 2 O 3 phase. The reduction of hematite to metal iron involves the mixed II-III valence oxide Fe 3 O 4 (magnetite) as intermediate [51].

Re-oxidation with 2% NO in He (500 °C). After the NO treatment, the Fe 2p photopeak returned to the very same shape observed after stabilisation under O 2 /He, indicating that iron was fully re-oxidised into 100% of Fe 3+ . On the other hand, Mn 2p signal, denoted that Mn remained partially reduced with 47% Mn 2+ and 53% Mn 3+ (average oxidation state=2.7±0.1). This means that only 30% of the Mn 2+ re-oxidised back to Mn 3+ when heated up to 500 °C under NO.

In brief, our quasi-in situ XPS investigation has made it possible to highlight different redox paths in Cu-dLFO and Mn-dLFO catalysts. In Cu-dLFO, the CuO phase participates in the redox process under CO/NO reaction completely, while the iron majorly remains as Fe 3+ and embedded in the LaFeO 3 phase. The copper is reduced to metal Cu 0 nanoparticles under CO and re-oxidised to Cu 2+ as CuO and Cu + as Cu 2 O under NO. The CO/NO cycling is accompanied by a sintering of the copper phase, as revealed by the lowering of both A/B and Cu/Fe ratios. In Mn-dLFO, both Fe and Mn are reduced under CO. Two third of Mn 3+ involved in the Mn-dLFO phase is reduced to Mn 2+ on surface while a minor part of Fe 3+ is reduced to metallic iron and likely some Fe 2+ . The latter are fully re-oxidised to Fe 3+ under NO, whereas the re-oxidation of Mn to Mn 3+ is not complete and thus, led to a balanced mixture of Mn 2+ and Mn 3+ . The different reduction paths observed for Mn 3+ and Fe 3+ cations can be further made clear using the iron-oxygen thermal diagram. Indeed, it is widely reported that below 570 °C, the nominal Fe(+II) oxide FeO (Wüstite) decomposes into metal iron and magnetite according to the formal equation [52][53][54]:

4F eO -→ F e 0 + F e 3 O 4 (1)

Conclusion

Our present investigations have highlighted some major differences in the catalytic operating modes of two formally B-site, 20% substituted perovskites, namely Cu-dLFO and Mn-dLFO. Under dehydrated conditions, Cu-dLFO consists of a LaFeO 3 -supported CuO phase together with a minority iron oxide phase on the catalyst's surface. Under decoupled CO/NO model conditions, CuO is evidenced as being the only redox-active phase. We have also highlighted that all copper moieties are efficiently reduced from Cu 2+ to metal Cu 0 under CO at 500 °C, while only partly re-oxidised to a mixture of +II and +I cations. Moreover, the Cu-dLFO catalyst retains a noticeable activity in CO oxidation at mild temperatures (around 200 °C), which can indicate the presence of a large pool of highly active surface oxygen moieties.

In contrast, the Mn-dLFO perovskite mostly consists of the perovskite phase together with an XRDand Raman-blind surface Fe 2 O 3 nano-phase, further confirmed by the XANES analysis. The Mn 3+ ion causes strong JT effect that results in a complete change in the Raman fingerprint as compared to the parent structure of LaFeO 3 , while the reduction of Mn 3+ to Mn 2+ under CO allows the partial recovery of LaFeO 3 signature. The reduction to Mn 2+ in presence of CO was confirmed by the evolution of the Mn 2p photopeak during the quasi-in situ XPS experiment. Interestingly, the Fe 2 O 3 phase was also found to undergo redox behaviour, which could also play a role in the as-observed activity of the Mn-dLFO catalyst.

Figure 1 :

 1 Figure 1: X-ray diffraction pattern of a) Cu-dLFO and b) Mn-dLFO.

Figure 2 :

 2 Figure 2: XPS Fe 2p and Cu 2p spectra of Cu-substituted perovskite powder, Cu-dLFO.

Figure 4 :

 4 Figure 4: Ex-situ Raman spectra of a) Cu-dLFO and b) LFO; λ exc. =488 nm (colours online).

Figure 6 :

 6 Figure 6: Compared XANES spectra of (a) Cu-dLFO (blue) and LFO (red dotted), (b) Mn-dLFO (green) and LFO (yellow dotted) recorded at the Fe K-edge in transmission mode. Full XANES regions are shown in the insets. LCF analysis of the XANES spectra of (c) Cu-dLFO and (d) Mn-dLFO using Febased references.

Figure 7 :

 7 Figure 7: ToF-SIMS mass spectra of the mixed ionic fragments La 2 FeMO + 5 and La 3 FeMO + 6 (M = Cu or Mn) from Mn-dLFO (a and c) and Cu-dLFO (b and d), respectively, analyzed in positive mode.

Figure 8 :

 8 Figure 8: In situ Raman analysis of Cu-dLFO during TPR under 5% CO/He at a) 25 °C, b) 100 °C, c) 200 °C, d) 250 °C, e) 300 °C, f) 350 °C and g) 400 °C. The Raman spectra were collected after 30 min stabilisation at a given temperature.

Figure 9 :

 9 Figure 9: CO consumption rate as a function of temperature during operando Raman analysis of La 0.67 Fe 0.8 Cu 0.2 O 3 under 5% CO/He; CO-TPR (left) and time representation (right).

Figure 10 :

 10 Figure 10: Raman spectra of Mn-dLFO during CO-TPR under 5% CO/He at a) 100 °C, b) 200 °C, c) 300 °C, d) 350 °C, e) 400 °C and f) 450 °C; Intensities of spectra c-f were magnified as indicated.

Figure 11 :

 11 Figure 11: Raman spectra of Mn-dLFO recorded at 100 °C at beginning (black), after CO-TPR (red) under 5% CO/He and their 10X magnified difference spectrum (green).

Figure 12 :

 12 Figure 12: CO consumption rate as a function of temperature during operando Raman analysis of Mn-dLFO under 5% CO/He; CO-TPR (left) and time representation (right).

Figure 13 :

 13 Figure 13: Quasi-in situ XPS spectra of Cu-dLFO powder recorded under UHV as-received (blue) and after reaction at 500 °C under 20% O 2 /He (green), 5% CO/He (black) and 2% NO/He (pink).

Figure 14 :

 14 Figure 14: Quasi-in situ XPS spectra of Mn-dLFO powder recorded under UHV as-received (blue) and after reaction at 500 °C under 20% O 2 /He (green), 5% CO/He (black) and 2% NO/He (pink).

Table 1 :

 1 Experimental bulk composition in Cu-dLFO and Mn-dLFO catalysts.

	Sample	Composition
	Cu-dLFO 0.7 LaFeO 3 0.2 CuO 0.05 Fe 2 O 3
	Mn-dLFO	0.7 LaFe 0.7 Mn 0.3 O 3	0.15 Fe 2 O 3

sition, our XPS results suggest that the surface of Cu-dLFO is largely covered by the segregated CuO phase, whereas there is no hint to support the surface localisation of the exsolved α-Fe 2 O 3 phase in Mn-dLFO.

Table 2 .

 2 After the O 2 step, the binding energy of Fe 2p is 710.4 eV, featuring a +III valence state, similar to the as-received state. In case of the Cu 2p, the binding energy is slightly shifted to 934.2 eV still retaining the intense satellite peak typical of Cu 2+ species. Besides, the La/(Fe+Cu) and Cu/Fe atomic ratios remains unchanged at 0.64 and 0.34, respectively.

Table 2 :

 2 Summary of the data extracted from the quasi-in situ XPS spectra of Cu-dLFO; : oxidation state. CuO, Cu 2 O undergoes drastic changes under CO: the Cu 2+ satellite disappeared completely, the binding energy is shifted from 934.2 eV to 932.1 eV and the peak width is reduced. The latter is indicative of reduced copper species such as Cu 0 or Cu + . Since binding energy for both are experimentally expected at the same value, the modified Auger parameter (α) can be used for accurate assignment of the oxidation state for copper, given by, α = E k(LM M ) + E b(Cu) , where E k(LM M ) is the kinetic energy of the LMM Auger peak of copper (see supporting information) and E b(Cu) is the binding energy of the Cu 2p 3/2

		Fresh	After O 2 After CO	After NO
	A/B	0.68	0.64	1.14	0.82
	Cu/Fe	0.33	0.34	0.11	0.17
	Cu os +II, (+I/0)	+II	0	+I, +II
	Fe os	+III	+III	+III	+III
	Fe as	LaFeO 3	LaFeO 3	LaFeO 3	LaFeO 3
	Cu as	CuO	CuO	Cu 0

Table 3 :

 3 Summary of the data extracted from the quasi-in situ XPS and XANES spectra of Mn-dLFO; : oxidation state. After CO exposure, the Fe 2p signal exhibited the presence of metallic iron at 706.8 eV, corresponding to 6% of the total iron probed by XPS. The presence of Fe 2+

		Fresh After O 2	After CO	After NO
	A/B	0.72	0.69	1.37	1.08
	Mn/Fe	0.26	0.23	0.31	0.34
	Mn o.s.	+III	+III	+II, +III	+II, +III
	Fe o.s.	+III	+III	0, (+II) +III	+III
	Reduction under 5% CO/He (500 °C).

This complex feed involves CO and NO together with hydrocarbons and H

O.

Acknowledgements

This work was funded by the European Commission as part of the PARTIAL-PGMs project under Grant Agreement #686086.

Chevreul Institute (FR 2638), Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation, Hauts-de-France Region and FEDER are acknowledged for supporting and funding partially this work.

The authors are thankful to the SAMBA beamline staff members at Synchrotron SOLEIL for their precious help in running the XAS experiments.