Hélène Barucq

Julien Diaz

Véronique Duprat

Long-term stable acoustic absorbing boundary conditions for regular-shaped surfaces

Keywords: Absorbing boundary conditions, Acoustic wave equation, Well-posedness Conditions aux limites absorbantes, Équation des

When dealing with regular shaped boundaries, the most widely used Absorbing Boundary Condition is a first-order condition including the curvature of the artificial surface. It is easy to include in variational numerical methods and it performs well. In this paper, we prove that this condition is equivalent to a second-order condition. Actually, we establish that it belongs to a one-parameter family of ABC that lead to well-posed and long-term stable problems. We then propose a 2D performance analysis of the different ABCS, including BGT-like conditions also and we conclude that the first-order condition is the one which performs better.

Introduction

The numerical simulation of acoustic waves is often carried out after truncating the propagation domain. The size of the computational domain is then reduced and the spatial approximation can be done with finite element spaces. The corresponding boundary value problem is defined by the coupling of the initial system of equations with a boundary condition set on the external boundary limiting the computational domain. The boundary condition ends the problem and represents the behavior of any wave impinging the external boundary. A "good" condition should be the one for which the numerical waves are not reflected by the external boundary. The computed field would thus be exactly the restriction of the exact wave that propagates inside the larger domain for which the external boundary has been pushed away. The construction of "good" conditions, that are generally called Absorbing Boundary Conditions (ABCs), has been considered in many works. Regarding time-dependent problems, Engquist and Majda [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] have shown that the micro-diagonalization of the wave equation, and more generally of hyperbolic systems, let construct ABCs on arbitrarily-shaped surfaces. The reference work [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF] has been followed by a series of papers in which the issue was to derive high-order ABCs to improve the accuracy of the numerical waves by diminishing the amplitude of the reflected waves. For instance, Higdon [START_REF] Higdon | Numerical absorbing boundary conditions for the wave equation[END_REF] has derived high-order ABCs which are obtained as the combination of low-order ABCs. Grote and Keller [START_REF] Grote | Exact nonreflecting boundary conditions for the time dependent wave equation[END_REF] have proposed non reflecting ABCs which are obtained from the approximation of the Dirichlet-to-Neumann (DtN) operator of the sphere. There are many other very interesting papers on the subject and we refer to the review article [START_REF] Tsynkov | Numerical solution of problems on unbounded domains[END_REF] for a more exhaustive list of references than herein. Most of the papers addressed the question of constructing ABCs for flat boundaries. It is quite surprising because there exist plenty of applications where the boundary should be curved. For instance, as far as sonar is concerned, it seems more convenient to embed a submarine in an ellipsoid than in a cube. The size of the computational region is then reduced and the computational costs are thus optimized. In the case of harmonic problems, Antoine et al. [START_REF] Antoine | Bayliss-turkel like radiation conditions on surfaces of arbitrary shape[END_REF] have proposed a second-order ABC that outperforms the BGT [START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF] one and obviously any first-order condition. In the case of time-dependent problems, to the best of our knowledge, the same issue has not be addressed and it is not clear that a second-order ABC can perform better than a first-order one. The most widely used first-order condition that is employed for the acoustic wave equation reads as:

B n u `Bt u `κ 2 u " 0 on Σ ˆr0, T s. (1)
It involves the curvature κ of the boundary Σ. It thus takes the geometry of Σ into account and it is well-known that it performs better than the simplest condition:

B n u `Bt u " 0 on Σ ˆr0, T s.

(

) 2
while both are easy to implement in a variational formulation. Now it is not obvious that condition (1) is outperformed by higher-order conditions in the case of the acoustic wave equation.

In [START_REF] Trefethen | Well-posedness of one-way wave equations and absorbing boundary conditions[END_REF], it has even been shown that ABCs of order greater than 3 can lead to ill-posed mixed problems. Thus the question of which ABC performs better for the acoustic wave equation remains posed.

Another point to consider is the computational burden that are required by the model. Conditions of high-order are quite easy to implement if we introduce additional unknowns [START_REF] Givoli | High-order non-reflecting boundary scheme for time-dependent waves[END_REF]. These latter are only defined on the external boundary. Thus, they do not induce significant computational costs if the number of auxiliary unknowns is low. However, the introduction of a high number of auxiliary unknowns can become problematic by impeding the implementation of parallel computations because it cause hamper the load balancing. Finally, it was recently shown that we can significantly improve the performance of the ABCs themselves by combining the ABCs with a condition capable of reproducing the evanescent and/or creeping waves [START_REF]Complete Radiation Boundary Conditions: Minimizing the Long Time Error Growth of Local Methods. T. hagstrom and t. warburton[END_REF][START_REF] Barucq | Micro-differential boundary conditions modelling the absorption of acoustic waves by 2d arbitrarily-shaped convex surfaces[END_REF]. Following these works, it would not be wise to use high-order ABCs. That is why we investigate in this paper the performance of condition 1 compared to higher conditions. More precisely, our goal is to decide if condition 1 is really outperformed by a second order ABC. Regarding the second order ABC, we have decided to use a BGT-like condition. Our choice is guided by the fact that this condition written for monochromatic waves is more efficient compared to other second order conditions [START_REF] Antoine | Bayliss-turkel like radiation conditions on surfaces of arbitrary shape[END_REF].

The paper is organized as follows. In section 2, we show how to derive a family of ABCs for the wave equation by applying a micro-diagonalization process that was formerly introduced by M.E. Taylor in [START_REF] Taylor | Pseudodifferential Operators[END_REF]. We then get ABCs depending on a parameter γ and at Section 3, we establish that the corresponding mixed problem is well-posed. Moreover, we prove that if γ " κ{4, the condition is equivalent to condition [START_REF] Agut | Stability analysis of the interior penalty discontinuous galerkin method for the wave equation[END_REF], which indicates that condition (1) is more than of order 1. In Section 4, we show that if γ ą" κ{4, the corresponding numerical solution behaves like the exact solution, that is the solution converges to zero when the time t goes to infinity. Finally, we address the question of evaluating the numerical performances of each condition. We conduct a performance analysis by comparing the BGT-like condition with the family of ABCs knowing that for γ " γ{4, it corresponds to condition [START_REF] Agut | Stability analysis of the interior penalty discontinuous galerkin method for the wave equation[END_REF]. The numerical experiments are carried out with a Discontinuous Galerkin code using a penalized formulation of the wave equation. Our conclusion (1) provides equivalent or better results when compared with all the other conditions, including the BGT-like one.

A new family of second-order ABCs for the acoustic wave equation

In this section, we construct a new absorbing boundary condition using the micro-diagonalization method developed by M.E. Taylor [START_REF] Taylor | Pseudodifferential Operators[END_REF]. This condition is written for regular domains and takes propagating waves into account. In the next Subsection, we present the main steps of the micro-diagonalization method and refer to [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF][START_REF] Antoine | Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism[END_REF][START_REF] Barucq | A new family of second-order absorbing boundary conditions for the acoustic wave equation -Part I: Construction and mathematical analysis[END_REF] for more details. For the sake of simplicity, the velocity c is supposed to be equal to 1, but our work can be easily reproduced for any velocity.

The micro-diagonalization method applied to the wave equation

We tackle the construction of ABCs by applying a factorization theorem proposed by M.E. Taylor [START_REF] Taylor | Pseudodifferential Operators[END_REF] to study the propagation of singularities of strictly hyperbolic systems. Because we want to build low-order conditions, we limit our work to the application of the first step of factorization. We are thus only dealing with the diagonalization of the principal symbol of the wave equation, according to the mathematical analysis in [START_REF] Antoine | Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism[END_REF], which shows that the next steps necessarily involve differential operators with order higher than two.

The ABCs that we consider are derived from the micro-local approximation of the DtN operator related to the artificial surface Σ. We thus begin with rewriting the acoustic wave equation in a local coordinate system pr, sq. The couple pr, sq describes a point in the neighborhood of Σ in such a way that Σ " tr " 0u. We use the same coordinate system as in [START_REF] Antoine | Bayliss-turkel like radiation conditions on surfaces of arbitrary shape[END_REF] and the acoustic wave equation reads then as

B 2 t u ´B2 r u ´κr B r u ´h´1 B s `h´1 B s u ˘" 0, (3)

Inria

where κ is the curvature of Σ, h " 1 `rκpsq and κ r " h ´1κ.

Next, to apply Taylor's method, we rewrite (3) as a first-order system. We thus introduce an auxiliary unknown v which satisfies B t v `Br u " 0 in a neighborhood of Σ and such that the field U " pv, uq is solution to the first-order system B r U " LU. The entries of L are first-order pseudo-differential operators and σpLq " L " L 1 `L0 is given by

L 1 " ˆ0 ´h´2 ξ 2 ´ω2 iω ´iω 0 ˙P S 1 and L 0 " ˜´κ r ´h´3 Bsphqξ ω 0 0 ¸P S 0 ,
where ξ and ω are the dual variables respectively associated to s and t. We use standard notations such as σpLq for the symbol of L. Thereafter, as in [START_REF] Taylor | Pseudodifferential Operators[END_REF], we note more concisely by S m the symbol class S m 1,0 introduced by [START_REF] Hörmander | Pseudodifferential operators and hypoelliptic equations[END_REF]. We also denote by τ ´mpσpLqq the truncation at the order m of the asymptotic expansion of the symbol σpLq. Let λ 1 denote the symbol λ 1 " `h´2 ξ 2 ´ω2 ˘1{2 . Then, when λ 1 ‰ 0, the principal symbol of L admits two single eigenvalues λ 1 and ´λ1 . When h ´2ξ 2 ´ω2 ą 0, pω, ξq belongs to the elliptic region and λ 1 is real. When h ´2ξ 2 ´ω2 ă 0, λ 1 is imaginary and corresponds to propagating modes for which pω, ξq covers the hyperbolic region. Following [START_REF] Antoine | Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism[END_REF][START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF], we introduce

V 0 " 1 ? 2 ˆiω λ1 1 1 ´λ1 iω ˙, so that V 0 L 1 V ´1 0 " ˆλ1 0 0 ´λ1 ˙" M 1 .
Now, a first-order ABC can be derived from the approximation and the localization of the global boundary condition rpI `K´1 q V 0 Us 2 " 0 on Σ,

where r¨s 2 represents the second component of the vector, V 0 is the Fourier Integral Operator (FIO) with symbol V 0 and K ´1 is a regularizing operator of order -1.

According to [START_REF] Taylor | Pseudodifferential Operators[END_REF], the principal symbol of K ´1 is fixed by imposing that W 1 " pI `K´1 q V 0 U is solution to B r W 1 " pΛ 1 `Λ0 q W 1 `R´1 W 1 , where Λ 1 is the FIO with the symbol M 1 and Λ 0 is a diagonal FIO also (see [START_REF] Antoine | Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism[END_REF][START_REF] Sesquès | Conditions aux limites artificielles pour le système de Maxwell[END_REF]). Indeed, by setting (4) on Σ, we impose that the in-going component of W 1 vanishes on Σ. By this way, we select only the out-going component of W 1 and an ABC is obtained after localizing (4). Therefore, the principal symbol of K ´1, denoted by K ´1, is determined by the condition : Λ 0 is a diagonal operator. Observe that this condition fixes only the off-diagonal entries of K ´1 (see [START_REF] Antoine | Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism[END_REF]) and most of the works using this approach consider

K ´1 " ˜0 ´iκω 3 4λ 4 1 ´iκω 4λ 2 1 0 ¸. (5)
Then, after using a truncated Taylor expansion of λ 1 in the hyperbolic region with ω 2 " h ´2ξ 2 , (4) gives rise to B t u `Bn u `κ 2 u " 0 on Σ. The resulting condition involves differential operators but it should be called micro-differential since it is justified in the propagating cone ω 2 " h ´2ξ 2 . This condition is widely used in case of curved surfaces. We will refer to this condition as the C-ABC for "Curvature ABC". Now, it might be interesting to see how the condition changes by considering non-zero diagonal entries for K ´1. We address this issue in the following Subsection.

A new family of ABCs

The idea of modifying K ´1 by introducing non-zero diagonal entries has been formerly applied in [START_REF] Barucq | A new family of first-order boundary conditions for the maxwell system : derivation, well-posedness and long-time behavior[END_REF] for the 2D Maxwell system, but only from a theoretical point of view. To the best of our knowledge, the numerical impact of using these conditions has never been investigated. Herein, we propose to modify K ´1 as follows

K ´1 " ¨0 ´iκω 3 4λ 4 1 ´iκω 4λ 2 1 γpsq λ 1 ‹ ‹ ' , (6)
where γ is a parameter depending only on the curvilinear abscissa s. Let us note that we do not change the first diagonal entry since it is not involved in condition (4). We then get Theorem 2.1 A family of first-order condition is given by

B t pB n u `Bt uq " ´κ 4 ´γ¯B n u ´´κ 4 `γ¯B t u on Σ, ( 7
)
where γ is a regular function defined on Σ.

Proof : We recall that the first-order boundary condition is given by rpI `K´1q V0Us 2 " 0 on Σ. The symbol of the corresponding operator reads as σ ppI `K´1q V0q " pI2 `K´1q V0 `R´2, where R´2 P S ´2. Therefore, the truncation of σ ppI `K´1q V0q in S ´1 is given by

τ´1 ppI `K´1q V0q " pI2 `K´1q V0 " 1 ? 2 ¨iω λ1 ´iκω 3 4λ 4 1 1 `κω 2 4λ 3 1 κω 2 4λ 3 1 `1 `γpsq λ1 ´iκω 4λ 2 1 ´λ1 iω ´γpsq iω ‹ ‹ ' . (8)
Using a first-order Taylor expansion for ω " h ´1ξ, we then obtain

τ´1 pI2 `K´1q V0 " 1 ? 2 ¨1 ´iκ 4ω 1 `iκ 4ω 1 `iκ 4ω `γ iω iκ 4ω ´1 ´γ iω '. (9)
Since U " pv, uq, with Btv `Bnu " 0 on Σ, is supposed to satisfy (4), we get Bt pBnu `Btuq " ´κ 4 ´γ¯B nu ´´κ 4 `γ¯B tu on Σ, When γ " κ 4 , we have

τ ´1 pI 2 `K´1 q V 0 " 1 ? 2 ¨1 ´iκ 4ω 1 `iκ 4ω 1 iκ 2ω ´1 ‹ '
and the boundary condition reduces to B n u `Bt u `κ 2 u " 0 on Σ. We thus recover the C-ABC for the particular choice γ " κ 4 .

Inria

Before we focus on the property (3) defined in the introduction to this paper, we wish to go back over the construction method that we used. Indeed it might seem more natural to apply the Nirenberg factorization technique [START_REF] Nirenberg | Pseudodifferential operators and some applications[END_REF] in particular because it is well-suited for scalar equations like the acoustic wave equation. But, if we set the same degree of approximations on the symbols, this process will only lead to the construction of the C-ABC, whereas the one used herein leads to a family of ABCs. Our approach is thus more productive for the same fixed level of approximation of the symbol and the same amount of computations.

Mathematical analysis

We consider the solution of the scattering by a sound soft obstacle when using the family of ABCs that have been constructed in the previous section. The corresponding BVP reads as : find u solution to

B 2 t u ´∆u " 0 in Ω ˆp0, `8q ; up0, xq " u 0 pxq, B t up0, xq " u 1 pxq in Ω;
u " 0 on Γ ˆp0, `8q ;

B t pB n u `Bt uq " ´κ 4 ´γ¯B n u ´´κ 4 `γ¯B t u on Σ ˆp0, `8q . (10)
Now, we aim at numerically solving [START_REF] Barucq | Micro-differential boundary conditions modelling the absorption of acoustic waves by 2d arbitrarily-shaped convex surfaces[END_REF] and, as a matter of fact, the condition on Σ is not easy to include in a variational formulation of the acoustic wave equation. Indeed, the boundary condition on Σ steps in the numerical scheme through the formal integrand

¯´1

. The sparsity of the finite element matrix is thus violated. This is why we propose an equivalent formulation involving an auxiliary unknown ψ defined on Σ. It is

B 2 t u ´∆u " 0 in Ω ˆp0, `8q ; ´Bt ´κ 4 `γ¯ψ " B t u, on Σ ˆp0, `8q ; B n u `Bt u `κ 2 ψ " 0 on Σ ˆp0, `8q ; u " 0 on Γ ˆp0, `8q ; up0, xq " u 0 pxq, B t up0, xq " u 1 pxq in Ω; ψp0, xq " ψ 0 pxq on Σ. (11)
The initial datum ψ 0 is given as a function of u 0 and u 1 on Σ satisfying the compatibility relation

B n u 0 `u1 `κ 2 ψ 0 " 0 on Σ.
Remark 3.1 Observe that in the context of an ABC method, Σ must be chosen in such a way that the supports of u 0 and u 1 are fully embedded inside Ω. In that case, we thus have ψpx, 0q " ψ 0 " 0 on Σ.

In the following, the function κ always denotes the curvature of Σ and γ is a regular parameter defined on Σ. The domain Ω is a bounded domain and its boundary BΩ " Γ Y Σ is assumed to be regular, with Γ X Σ " ∅ (see Fig. 1).

We do not consider the case γ " κ 4 because it corresponds to a particular case where the ABC is the C-ABC. Then, the BVP (10) can be studied directly, without considering the auxiliary problem [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF]. We thus make do with setting existing results that are proved in [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF]. We first transform [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF] in a first order system in time. We introduce an auxiliary unknown v defined by v " B t u. The vector U " pu, v, ψq is thus solution to

dU dt " AU, A " ¨0 Id 0 ∆ 0 0 0 1 κ 4 ´γ ' (13)
with the boundary conditions u " 0 on Γ ˆp0, `8q and B n u `Bt u `κ 2 ψ " 0 on Σ ˆp0, `8q .

In the following, we will concentrate on formulation (13) and we will interest ourselves on its solution in suitable Hilbert spaces. As usual, H s denotes the Sobolev space of order s P R. Let us first introduce H as the product space defined by H " H 1 Γ pΩq ˆL2 pΩq ˆL2 pΣq, where H 1 Γ pΩq " th 1 P H 1 pΩq , h 1 " 0 on Γu. We equip H with the Hilbertian graph norm

} ph 1 , h 2 , h 3 q } H " ´}h 1 } 2 L 2 pΩq `}∇h 1 } 2 L 2 pΩq `}h 2 } 2 L 2 pΩq `}h 3 } 2 L 2 pΣq ¯1{2 .
Let V be the product space defined by

V " tpv 1 , v 2 , ϕq P H, Apv 1 , v 2 , ϕq P H, B n v 1 `v2 `κ 2 ϕ " 0 on Σu.
The space V corresponds to the domain of A. By enforcing Apv 1 , v 2 , ϕq P H, the regularity of the unknown is improved. Indeed, we then have v 2 P H 1 Γ pΩq and ∆v 1 P L 2 pΩq. Then,

v 2 P H 1 Γ pΩq implies that v 2| Σ is defined in H 1{2 pΣq and ∆v 1 P L 2 pΩq implies that B n v 1| Σ P H ´1{2 pΣq, knowing that v 1 P H 1 pΩq. Moreover, the relation B n v 1 `v2 `κ 2 ϕ " 0 on Σ improves the regularity of B n v 1| Σ since v 2 `κ 2 ϕ P L 2 pΣq. Hence, to sum up, we have V " tpv 1 , v 2 , ϕq P H, ∆v 1 P L 2 pΩq, v 2 P H 1 Γ pΩq, B n v 1| Σ P L 2 pΣq, B n v 1 `v2 `κ 2 ϕ " 0 on Σu.

Inria

We recall the Green formula we will use: for all pu, vq P H 1 pΩq ˆH1 pΩq such that ∆u P L 2 pΩq, we have ż

Ω ∆uv dx " ´żΩ ∇u ¨∇v dx `xB n u, vy H ´1{2 pBΩq,H 1{2 pBΩq . ( 14
)
Lemma 3.2 Let κ be given in L 8 pΣq and such that min xPΣ κpxq " κ 0 ą 0. Then, for all h P H, the quantity ~h~"

ˆżΩ |∇h 1 | 2 `|h 2 | 2 dx `żΣ κ 2 |h 3 | 2 dσ

˙1{2

is a norm on H equivalent to the norm }h} H .

Proof : It is well-known that }∇ ¨}L 2 pΩq defines a norm on H 1 Γ pΩq which is equivalent to the standard norm in H 1 pΩq, as a consequence of the Poincaré inequality. Hence, since the curvature κ is supposed to be in L 8 pΣq with min xPΣ κpxq ą 0, it is straightforward that ~¨d efines a norm on H equivalent to the conventional norm } ¨}H .

In the following, we denote by p , q the scalar product derived from the norm ~~. Lemma 3.3 Let κ and γ be given such that ( 12) is checked. Then, for all v P V , we have pAv, vq ď 0.

Proof : Let v " pv1, v2, ϕq in V . By definition of A, Av " `v2,

Lemma 3.4

The operator A, with domain V , is maximal.

RR n°8203

Proof : Given f " pf1, f2, f3q in H, we consider the following mixed problem: find v P V such that pA ´Iq v " f . We thus seek v " pv1, v2, ϕq P V such that

v2 ´v1 " f1 in Ω; ∆v1 ´v2 " f2 in Ω; v2 `´κ 4 ´γ ´1¯ϕ " f3, Bnv1 `v2 `κ 2 ϕ " 0 on Σ; v1 " 0 on Γ; (15)
Assume now that the problem (15) has a solution in V . Then, by removing v2 thanks to the first equation, v2 " f1 `v1, in Ω, and ϕ, thanks to the third equation, ϕ " f 3 ´f1 ´v1 κ 4 ´γ´1 , we obtain that v1 is solution to the boundary-value problem

´∆v1 `v1 " f in Ω;
v1 " 0 on Γ;

Bnv1 `αv1 " g on Σ.

(

) 16
with f :" ´pf2 `f1q in L 2 pΩq, α " 1 `κ 2p1`γ´κ 4 q ą 0 and

g :" ´˜1 `κ 2 `1 `γ ´κ 4 ˘¸f1 `κ 2 `1 `γ ´κ 4 ˘f3 in L 2 pΣq.
It is obvious that, since κ P L 8 pΣq and α ą 0 by hypothesis, the functional |v|1,α " ˆ}v}

ż Ω f φ dx. Moreover, φ | Γ " 0. Therefore, @φ P T `Ω˘, ż Ω ∇v1∇φ dx ´żΣ Bnv1φ dσ `żΩ v1φ dx " ż Ω f φ dx.
Then, by using that Bnv1 " g ´αv1 on Σ, we obtain

@φ P T `Ω˘, ż Ω ∇v1∇φ dx `żΣ αv1φ dσ `żΩ v1φ dx " ż Ω f φ dx `żΣ gφ dσ. ( 17
)
Let ap¨, ¨q be the bilinear form defined by apv1, φq "

ż Ω ∇v1∇φ dx `żΣ αv1φ dσ `żΩ v1φ dx.
It is obvious that ap¨, ¨q is continuous on H

@φ P DpΩq Ă H 1 Γ pΩq, ż Ω ∇v1∇φ dx `żΩ v1φ dx " ż Ω f φ dx.
We then deduce that @φ P DpΩq, xv1 ´∆v1 ´f , φy " 0, where x , y denotes the duality pairing between D 1 pΩq and DpΩq, and we thus obtain v1 ´∆v1 " f in D 1 pΩq.

This identity allows us to give a sense to ∆v1 in L 2 pΩq. Therefore, Bnv1 | BΩ P H ´1{2 pBΩq and we also have

@φ P T `Ω˘, ż Ω ∇v1∇φ dx `żΣ αv1φ dσ `żΩ v1φ dx " ż Ω f φ dx `żΣ gφ dσ.
Using the Green formula ( 14), we get, @φ P T `Ω˘: The existence of v1 as a solution to ( 16) is thus proved. Next, since v1 and f1 are in H 1 pΩq, we deduce the existence of v2 " f1 `v1 in H 1 pΩq. Moreover since v2 | Σ and f3 are in L 2 pΣq, ϕ " f3 ´f1 ´v1 κ 4 ´γ ´1 exists in L 2 pΣq. To complete the proof, we have to check that pv1, v2, ϕq P V , which is obvious, using the relations v2 " f1 `v1 in Ω and αϕ " f3 ´v2 on Σ and that f1 P H 1 Γ pΩq and f3 P L 2 pΣq.

´żΩ ∆v1φ dx `xBnv1,
We are now willing to set the main result of this section:

Theorem 3.5 Let pu 0 , u 1 , ϕ 0 q P V . The problem (11) admits a unique solution u such that pu, B t u, ϕq P C 1 pr0, `8r; V q X C 0 pr0, `8r; Hq.

Proof : The two previous lemmas show that the operator A is a maximal dissipative operator in its domain V . According to the Hille-Yosida theorem [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF], the problem [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF] has one and only one solution U " pu, v, ψq such that pu, v, ψq P C 1 pr0, `8r; V q X C 0 pr0, `8r; Hq. Remark 3.6 A is thus the infinitesimal generator of a semi-group of contraction Zptq. We can thus define a finite energy solution of [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF] with initial data pu 0 , u 1 , ψ 0 q in H only. In that case, pu, v, ψq " Zptqpu 0 , u 1 , ψ 0 q P C 1 pr0, `8r; V q X C 0 pr0, `8r; Hq.

When γ " κ 4 , we have [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF] :

Theorem 3.7 Let pu 0 , u 1 q P W 1 with W 1 " ! pu, vq P H 1 pΩq ˆH1 pΩq, ∆u P L 2 pΩq, B n u `v `κ 2 u " 0 on Σ ) .
Then the problem (10) admits a unique solution u such that pu, B t u, ϕq P C 1 pr0, `8r; W 1 q X C 0 pr0, `8r;

W 0 q, ( 18
)
where W 0 " H 1 Γ pΩq ˆL2 pΩq.

We have established that the problem is well-posed if κ ą 0 and γ ě κ 4 . The problem admits thus a solution when the external boundary Σ is convex. In the framework of ABCs, this is not restrictive since Σ is chosen by the engineer. Nevertheless, the well-posedness can also be obtained without any assumption on κ and γ, except suitable regularity hypotheses. Indeed, by multiplying u by e ´δt with δ ą maxp0, max xPΣ κpxq 4 ´γpxq, ´min xPΣ κpxq 4 `γpxqq, we get an auxiliary and equivalent problem for which a proof of the existence and uniqueness is achieved for any κ and γ. We refer to [START_REF] Barucq | A new family of second-order absorbing boundary conditions for the acoustic wave equation -Part I: Construction and mathematical analysis[END_REF] for a detailed proof for the general case. Let us also mention that the well-posedness of the initial formulation [START_REF] Barucq | Micro-differential boundary conditions modelling the absorption of acoustic waves by 2d arbitrarily-shaped convex surfaces[END_REF] has been established in [START_REF] Barucq | A new family of second-order absorbing boundary conditions for the acoustic wave equation -Part I: Construction and mathematical analysis[END_REF] and the case of a Neumann boundary condition is also considered. Herein, we decided to select a set of suitable values for κ and γ that allows to consider a large rank of applications and that provides simplified but informative proofs.

Regarding applications, it is interesting to consider the case where the acoustic wave is generated by a source field S. Then the initial data vanish and, in practice, S is compactly supported in a time interval r0, T s. We can apply the Hille-Yosida theory providing S P C 1 pr0, `8r, Hq. The wave field can then be written thanks to the Duhamel formula in r0, T s and for t ě T , the previous study is still valid by considering the solution pu, B t u, ψq at t " T as an initial wave field.

Long-term behavior

The results of Section 3 can be enhanced by introducing the functional defined on H by

Eph 1 , h 2 , h 3 q " 1 2 ż Ω `|∇h 1 | 2 `|h 2 | 2 ˘dx `1 2 ż Σ κ 2 |h 3 | 2 dσ. (19)
If κpxq ě 0 for all x P Σ, E is positive on H and if κ ą 0 on Σ, E 1{2 is obviously a norm in H equivalent to the norm ~~, according to Lemma 3.2. Then, Epu, B t u, ψq defines an energy on H. Moreover, Lemma 4.1 For all pu 0 , u 1 , ψ 0 q P V, t Þ ÝÑ Epu, B t u, ψq is differentiable and is decreasing under the condition (12): γpxq ě κpxq 4 , @x P Σ.

Proof : In the previous section, we have seen that if the initial conditions pu0, u1, ψ0q are in V , Epu, Btu, ψq P C 1 pr0, `8rq. Moreover, Remark 4.2 Lemma 4.1 has been proved in [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF][START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF] when γ " κ 4 .

d
In the following, we will assume that ψ 0 " 0 on Σ, without being restrictive on the application domain of our work. Indeed, as it was previously discussed in Remark 3.1, ψ 0 must vanish on Σ in the framework of ABCs methods. Theorem 4.3 Under ( 12), for all pu 0 , u 1 , 0q P V , lim tÑ`8

Epu, B t u, ψq " 0.

Proof : We have already seen that A is the generator of a continuous contraction semigroup Zptq. As it is sufficient to prove the theorem on a dense subspace of V " DpAq, we consider the initial data pu0, u1, 0q in DpA 2 q, where DpA 2 q " tpv1, v2, ψq P V, Apv1, v2, ψq P V u is equipped with the graph norm }pv1, v2, ψq} DpA 2 q " }pv1, v2, ψq}V `}Apv1, v2, ψq}V `}A 2 pv1, v2, ψq}V .

For any solution to [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF], we have }pu, Btu, ψq} DpA 2 q " }Zptqpu0, u1, 0q} DpA 2 q " }Zptqpu0, u1, 0q}V `}ApZptqpu0, u1, 0qq}V `}A 2 pZptqpu0, u1, 0qq}V .

As A, A 2 and Zptq are commuting on DpA 2 q, }pu, Btu, ψq} DpA 2 q " }Zptqpu0, u1, 0q}V `}ZptqApu0, u1, 0q}V `}ZptqA 2 pu0, u1, 0q}V .

Since Zptq is continuous in V , we easily deduce that there exists a positive constant C such that }pu, Btu, ψq} DpA 2 q ď C}pu0, u1, 0q} DpA 2 q . We thus have a bounded sequence of solutions in DpA 2 q which implies that we can extract a subsequence denoted by Zpt k qpu0, u1, 0q which weakly converges to pu8, v8, ψ8q in DpA 2 q. Now, let us denote by pupt k q, vpt k qq the sequence that is converging to pu8, v8q. By definition of DpA 2 q, pupt k q, vpt k qq is bounded in H 3{2 pΩq Ĥ3{2 pΩq and ∆upt k q is bounded in H 1 pΩq. Indeed, any pu, vq in DpA 2 q satisfies u P H 1 pΩq, ∆u P H 1 pΩq, u " 0 on Γ, Bnu | Σ P L 2 pΣq, v P H 1 pΩq, ∆v P H 1 pΩq, v " 0 on Γ, Bnv | Σ P L 2 pΣq, and, since Ω is regular,

w P H 1 pΩq, ∆w P H 1 pΩq, w " 0 on Γ, Bnw | Σ P L 2 pΣq ( Ă H 3{2 pΩq. (21)
We can thus deduce that pupt k q, vpt k qq strongly converges to pu8, v8q in H 1 pΩq ˆH1 pΩq according to the compact injection of H 3{2 pΩq into H 1 pΩq. Moreover, ∆upt k q strongly converges in L 2 pΩq since it is bounded in H 1 pΩq. Since ∆upt k q converges to ∆u8 in D 1 pΩq, we have that ∆upt k q strongly converges to ∆u8 in L 2 pΩq, according to the limit uniqueness.

We then deduce that Bnupt k q | Σ strongly converges to Bnu8 | Σ in H ´1{2 pΣq and that vpt k q | Σ strongly converges to v8 | Σ in H 1{2 pΣq. This implies that ψpt k q strongly converges to ψ8 in H ´1{2 pΣq.

Nevertheless, this convergence result is not sufficient to prove that pupt k q, vpt k q, ψpt k qq strongly converges to pu8, v8, ψ8q in V , because ψpt k q should strongly converge in L 2 pΣq. That is why we go back to the equation satisfied by ψpt k q. By construction, ψ is solution to: which implies }ψ} H 1{2 pΣq ď 1 α min p1 ´e´α mint q}v} H 1{2 pΣq . Then, ψptq is uniformly bounded in H 1{2 pΣq and we can then deduce that ψpt k q strongly converges to ψ8 in L 2 pΣq, as a consequence of the compact injection from H 1{2 pΣq into L 2 pΣq. As a conclusion, pupt k q, vpt k q, ψpt k qq strongly converges to pu8, v8, ψ8q in V . Since t Þ ÝÑ Epu, Btu, ψq is continuous, we thus have lim tÑ`8 Epu, Btu, ψq " lim tÑ`8 EpZptqpu0, u1, 0qq

Btψ
" lim t k Ñ`8
EpZpt k qpu0, u1, 0qq " Epu8, v8, ψ8q.

Moreover, for all positive s, lim tÑ`8

EpZpt `sqpu0, u1, 0qq " lim

t k Ñ`8
EpZpsqZpt k qpu0, u1, 0qq " EpZpsqpu8, v8, ψ8qq.

Inria

Then, if pw, Btw, ϕq " Zptqpu8, v8, ψ8q denotes the solution to problem [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF], with initial data pu8, v8, ψ8q in DpAq, we have Epw, Btw, ϕq " Epu8, v8, ψ8q for all positive t. Hence, according to the proof of Lemma 4.1, since

d dt Epw, Btw, ϕq " ´"ż Σ κ 2 ´γ ´κ 4 ¯|ϕ| 2 dσ `żΣ |Btw| 2 dσ  ,
we necessarily have Btw " 0 on Σ and ϕ " 0 on Σ. We then deduce that w is solution to the following problem

B 2 t w ´∆w " 0 in Ω ˆr0, ` 8r
wpx, 0q " u8, Btwpx, 0q " v8 in Ω

w " 0 on Γ ˆr0, `8r

Bnw " 0 on Σ ˆp0, `8q and z :" Btw is solution to

B 2 t z ´∆z " 0 in Ω ˆr0, `8r zpx, 0q " v8, Btzpx, 0q " ∆u8 in Ω z " 0 on Γ ˆr0, ` 8r
Bnz " z " 0 on Σ ˆp0, `8q .

Since Bnz " z " 0 on Σ, we deduce that z " 0 in Ω ˆr0, `8r, as a consequence of the Holmgren theorem (see e.g. Lions [START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués, T.1: Contrôlabilité exacte[END_REF]). Therefore, w is solution to ∆w " 0 in Ω ˆr0, `8r wpx, 0q " u8 in Ω

w " 0 on Γ ˆr0, `8r

Bnw " 0 on Σ ˆp0, `8q which implies that w " 0 in Ω ˆr0, `8r since Ω is connected. The pair pu8, v8q is thus equal to zero. Consequently, we also get that ψ8 is also equal to zero and a fortiori, we have Epu8, v8, ψ8q " 0.

We are now willing to set Theorem 4.4 Let pu 0 , u 1 , 0q in V . If [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier-stokes equations[END_REF] is satisfied, the solution u to [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF] satisfies

lim tÑ`8
pu, B t u, ψq " p0, 0, 0q in H.

Proof : If the pair pu0, u1, 0q is in V , we know that lim tÑ`8

Epu, Btu, ψq " 0 (cf. Theorem 4.3) and that E is a norm on H, equivalent to } }. Therefore, lim tÑ`8 pu, Btu, ψq " p0, 0, 0q in H.

Theorem 4.4 is still valid when γ " κ 4 . The result has been proved in [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF].

Numerical Experiments

In this section, we analyze numerically the impact of the parameter γ on the accuracy of the ABC (Subsection 5.2), we compare the performance of the ABC to the performance of the BGT2 condition (Subsection 5.3), and we illustrate the exponential decay of the energy E (Subsection 5.4).

To perform the experiments, we have discretized [START_REF] Barucq | Etude asymptotique du système de maxwell en dimension deux d'espace avec la condition aux limites absorbante de silver-müller[END_REF] by the Interior Penalty Discontinuous Galerkin (IPDG) method, which is described in the Subsection 5.1.

General setting of the IPDG method

The IPDG method, also known as Symmetric Interior Penalty (SIP) method, has been introduced in [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF][START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier-stokes equations[END_REF] for general elliptic problems. It is a discontinuous Galerkin approximation method in which a penalization term is introduced to impose the weak continuity of the solution through each element of the mesh. As a Discontinuous Galerkin method, it allows to use block-diagonal (and thus easily invertible) matrices. Moreover, it has been shown in [5] that IPDG is one of the only Discontinuous Galerkin method that is both stable and consistant, which guarantees an optimal order of convergence. We consider a triangulation T h of Ω composed of elements K, triangles in 2D or tetrahedra in 3D. We denote by Ω h the set of elements, by Σ abs the set of the edges on the absorbing boundary Σ, by Σ D the set of the edges on Γ and by Σ i the set of the edges inside the domain. We thus have Σ i X pΣ D Y Σ abs q " H. For each σ P Σ i , two elements, denoted arbitrarily by K `and K

´, share σ. We define the jump and the average over each edge by:

rrvss :" v `ν``v´ν´a nd ttvuu :" v ``v 2 ,
where v `and v ´respectively refer to the restriction of v in K `and K ´and ν ˘stands for the unit outward normal vector to K ˘.

Regarding the external boundary condition, it reads as B n u " ´Bt u ´κ 2 ψ on Σ, with `Bt ´κ 4 `γ˘ψ " B t u on Σ. We seek an approximation of the solution u in the finite element space V k h defined for k P N by V k h " v P L 2 pΩq; v |K P P k pKq, @K P T h ( , k P N and an approximation of ψ in the finite element space W k h defined by W k h " w P L 2 pΣq; w |σ P P k pσq, @σ P Σ abs ( where P k pKq (respectively P k pσq) is the set of polynomials of degree at most k on K (respectively on σ).

The discrete problem is then given by Find u h P V k h ˆp0, `8q and ψ P W k h ˆp0, `8q such that, @pv h , w

h q P V k h ˆW k h , ÿ KPT h ż K B 2 t u h v h dx `apu h , v h q `ÿ σPΣ abs ż σ ´Bt u h `κ 2 ψ h ¯vh dσ " ÿ KPT h ż K Sv h dx, ÿ σPΣ abs ż σ ´Bt ´κ 4 `γ¯ψ h w h dσ " ÿ σPΣ abs ż σ B t u h w h dσ, with apu, vq " ÿ K ż K ∇u∇vdx ´ÿ σPΣi ż σ
´tt∇uuu rrvss `tt∇vuu rruss ´α h σ rruss rrvss ¯dσ.

The function h σ is the smallest radius of the inscribed circles of the two elements sharing σ and α P R `˚is a penalization coefficient [2]. We refer to [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF]5] for more details on IPDG methods, to [2,[START_REF] Agut | Stability analysis of the interior penalty discontinuous galerkin method for the wave equation[END_REF] for an analysis of the penalization coefficient and to [START_REF] Grote | Discontinuous galerkin finite element method for the wave equation[END_REF] for the application of the IPDG method to the discretization of the wave equation.

Inria

We consider in the following two dimensional problems. Then, for each element K (resp. for each edge σ) the dimension of P k pKq (resp. P k pσq) is n " C k k`2 (resp. m " k `1). If the mesh is composed of N e elements and N f edges, we then have N :" N e ˆn degrees of freedom to interpolate u h and M :" N f ˆm degrees of freedom to interpolate ψ h .

If tv i , 1 ď i ď N u denotes a basis of V k h and tw i ,

1 ď i ď M u is a basis of W k h , we set u h px, tq " N ÿ i"1 U i ptq v i pxq and ψ h px, tq " M ÿ i"1 Ψ i ptq w i pxq .
The algebraic form of the problem is thus given by

M d 2 U dt 2 `B dU dt `1 2 B κ Ψ `KU " S, C dΨ dt `Cκ, γ Ψ " D dU dt (22)
where U and Ψ are the vectors of unknowns and

M " ˜ÿ KPT h ż K v i v j dx ¸1ďi,jďN , K " papv i , v j qq 1ďi,jďN B " ˜ÿ σPΣ abs ż σ v i v j dσ ¸1ďi,jďN , B κ " ˜ÿ σPΣ abs ż σ κw j v i dσ ¸1 ď i ď N 1 ď j ď M C " ˜ÿ σPΣ abs ż σ w i w j dσ ¸1ďi,jďM , C κ, γ " ˜ÿ σPΣ abs ż σ ´γ ´κ 4 ¯wi w j dσ ¸1ďi,jďM , D " ˜ÿ σPΣ abs ż σ v j w i dσ ¸1 ď i ď M 1 ď j ď N , S " ˜ÿ KPT h ż K S v i dx ¸1ďiďN .
For the time discretization, we use a second-order finite difference scheme with a time step ∆t and we obtain

MX n`1 " ˆ∆t 2 S n 0 ˙`ˆ´∆ t 2 K `2M 0 ˙Xn `MX n´1 , (23)
where

M " ¨M `∆t 2 B ∆t 2 2 B κ ´D C `∆tC κ, γ ', M " ¨´M `∆t 2 B ´∆t 2 2 B κ ´D C ´∆tC κ, γ ' and X i " ˆUi Ψ i ˙.
In the following numerical experiments, we use P 3 elements and we choose α " 20, ∆t " 0.05h, where h denotes the radius of the circumscribed circle of the smallest cell of the mesh.

Remark 5.1 In practice, we do not invert the matrix M and we solve [START_REF] Sesquès | Conditions aux limites artificielles pour le système de Maxwell[END_REF] element by element.

To simplify the presentation, we assume, without loss of generality, that the degrees of freedom are numbered such that 1. the element K i is the support of functions `vnpi´1q`j ˘j"1,n ; 2. the edge σ i is the support of functions `wmpi´1q`j ˘j"1,m .

We then denote by U Ki (resp. Ψ σi ) the restrictions of U (resp. Ψ) to the element K i (resp. the edge σ i ): U Ki " `Unpi´1q`j ˘j"1,n , Ψ σi " `Ψmpi´1q`j ˘j"1,m , by M Ki and B Ki the restrictions of M and B to the element K i :

M Ki " `Mnpi´1q`j,npi´1q`k ˘j,k"1,n , B Ki " `Bnpi´1q`j,npi´1q`k ˘j,k"1,n ,
and by B κσ i , C σi , C κ,γ σi and D σi the restrictions of B κ , C, C κ,γ and D to the edge σ i .

B κσ i " `Bκnpj i´1q`k,mpi´1q`l ˘k"1,n l"1,m , C σi " `Cmpi´1q`k,mpi´1q`l ˘k"1,m l"1,m C κ,γ σi " ´Cκ,γ mpi´1q`k,mpi´1q`l ¯k"1,m l"1,m , D σi " `Dmpi´1q`k,npji´1q`l ˘k"1,m l"1,n
(if we suppose that the edge i belongs to the element j i ). Then, 1. if K i is an interior element (i.e. an element with no edge on Σ), we have to solve the system M Ki U n`1 Ki " G where G is a vector that only depends on U n , U n´1 and S n ;

2. if K i is a boundary element, with for instance two external edges σ j1 and σ j2 (in practice a boundary element has one or two external edge), we have to solve the system

¨MKi `∆t 2 B Ki ∆t 2 2 B κσ j 1 ∆t 2 2 B κσ j 2 ´Dσj 1 C σj 1 `∆tC κ, γ σj 1 0 ´Dσj 2 0 C σj 2 `∆tC κ, γ σj 2 ‹ ‹ ‹ ' ¨U n`1 Ki Ψ n`1 σj 1 Ψ n`1 σj 2 ‹ ‹ ' " G
where G is a vector that only depends on U n , U n´1 , Ψ n , Ψ n´1 and S n .

If K i has two external edges (resp. one), we thus have to invert a pn `2mq ˆpn `2mq matrix (resp. pn `mq ˆpn `mq.

In all cases, we just have to solve small systems easily invertible.

Note that, in the particular case of the curvature condition (i.e γ " κ{4), we have u " ψ and the system can be rewritten as

M d 2 U dt 2 `B dU dt `Bκ U `KU " S, (24)
so that there is no need to compute the auxiliary variable.

Accuracy of the ABC

Now, we want to analyze the performances of the ABC (7) compared to the ones of the C-ABC which corresponds to γ " κ 4 . To this aim, we will consider two configurations denoted respectively by Configuration 1 and Configuration 2. In Configuration 1, the two-dimensional domain Ω 1 (see Fig. 2) is delimited by a circular exterior boundary Σ 1 and the boundary Γ 1 represents the boundary of the obstacle. Both circles are centered at the origin. The radius of Computational domain -Configuration 2 Σ 1 is R ext " 3m and the radius of Γ 1 is R int " 1m. The point source S is set at (0m, 1.3m) and is defined as a second-derivative of a Gaussian with a dominant frequency f 0 of 1Hz:

Inria

S " δ x0 2λ ´λ pt ´t0 q 2 ´1¯e ´λpt´t0q 2 ,
with x 0 " p0m, 1.3mq, λ " π 2 f 2 0 and t 0 " 1{f 0 . Since the source decays exponentially in time for t ě t 0 , we consider that Sptq " 0 for t ě 5t 0

In Configuration 2, Ω 2 (see Fig. 3) is a two-dimensional domain delimited by an exterior boundary Σ 2 and by an interior boundary Γ 2 . Σ 2 is an ellipse of semi-major axis a ext " 6m and semi-minor axis b ext " 3m centered in the origin. Γ 2 is the boundary of an elliptical obstacle of semi-major axis a "2m and semi-minor axis b "1m centered at the origin. To compare the efficiency of the ABC [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF] for different values of γ, we first compute the relative L 2 -error in time at a given receiver set near the exterior boundary. To evaluate this relative error, we have to compute the exact solution at each receiver. We did not actually compute the exact solution, but we computed an approximate solution in a larger domain including Ω 1 and Ω 2 . Indeed, we multiplied the dimensions of Σ 1 and Σ 2 by 3. The relative L 2 px,yq pr0, T sq error at point px, yq is defined by ´şT 0 pu app pt, px, yqq ´uex pt, px, yqqq 2 dt

¯1{2 ´şT 0 pu ex pt, px, yqqq 2 dt ¯1{2 , (25)
where u app is the numerical solution and u ex is the exact solution computed inside the large domain. The error is computed after 6000 time steps (with a time step equal to 3.6 ˚10 ´3s). We consider three receivers set near the boundary with coordinates p0, 2.85mq, p0.7m, 2.75mq and p1.4m, 2.45mq for Configuration 1 and four receivers whose coordinates are p6m, 2.05mq, p5m, 2.3mq, p6m, ´2.05mq and p5m, ´2.3mq for Configuration 2. The results are given in Tab. 1 for Configuration 1 and Tab. 2 for Configuration 2. We can see that the relative L 2 -error does not vary significantly with the values of γ and is similar to the errors obtained with the C-ABC pγ " κ 4 q. The error seems to increase when γ is greater than κ but does not grow fast. For the next numerical tests, we fix γ " κ 4 . Now, we interest ourselves on the evolution of the errors when moving the absorbing boundary. That is to find out where the boundary should be to obtain suitable results without high computational burdens. We first consider Configuration 1. We have tested six different radius R ext : 1.5m, 2m, 3m, 4m, 5m and 6m in order to find the critical value of the radius of the exterior circle we should take to obtain accurate results. First, we set three receivers of coordinates p1.025m, 1.025mq, p1.256m, 0.725mq and p1.4m, 0.375mq near the obstacle and we evaluate the relative L 2 -error in time defined in [START_REF]Complete Radiation Boundary Conditions: Minimizing the Long Time Error Growth of Local Methods. T. hagstrom and t. warburton[END_REF]. The results we have obtained are presented in Tab. [START_REF] Antoine | Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism[END_REF] . By this way, the impact of measures near the external boundary is averaged and we thus get a better illustration of the accuracy of the solution. In order to compare accurately this error, we used the same time step for computing both the approximate and the reference solution. The results we have obtained are presented in Fig. 4. For external radius smaller than 3m, the computed relative errors are between 5 and 25%. When the external radius is greater than 4m, we obtain small relative errors around 2.5% and 1%. Hence, choosing an external radius of 4m when the interior radius is equal to 1m seems to be again the best choice to obtain good approximations of the solution of the wave equation without high computation burdens. Now, We have set three receivers near the obstacle of coordinates p0.75m, 1.4mq, p2.51m, 0.725mq and p2.8m, 0.375mq and we evaluate the relative L 2 -error in time defined in [START_REF]Complete Radiation Boundary Conditions: Minimizing the Long Time Error Growth of Local Methods. T. hagstrom and t. warburton[END_REF]. The results we have obtained are displayed in Tab. 4. In this case, we obtain really good results when a ext ě 4a and b ext ě 4b. If we take a ext " 6m and b ext " 3m we obtain suitable results but the difference with the case when a ext " 8m and b ext " 4m is important. Next, we compute the relative L 2 pr0, ts ˆΩq-error. The results are depicted in Fig. 5. We deduce that in such a configuration, we should take a ext " 3a and b ext " 3b to obtain small relative errors as in the case of a circular scatterer. In these numerical experiments, we have not observed any instabilities even in long-term simulations. Global relative error for pa ext , b ext q " p4, 2q using BGT2 condition on a larger time interval -Configuration 2 Inria able to compute the solution with BGT2 condition. One solution could be to use an implicit time-scheme instead of the Leap-Frog scheme in the whole domain. This would require the inversion of a huge matrix at each time step, which, added to the fact that we have to compute two auxiliary variables on the interface, would increase too strongly the computational costs.

As a conclusion, we recommend to use the curvature condition instead of the BGT2 condition, since this latter increases the computational costs and is unstable when coupled with an IPDG scheme. However, the experiments we have presented have been performed on relatively short time, and they do not prove the long-term stability of the curvature condition. This is the object of the next Subsection. 3 ensures that the energy E is decreasing. In practice, its discrete equivalent is actually exponentially decreasing. To illustrate this property, we have computed the evolution of the logarithm of the discrete energy (see Fig. 16). Let us remark that the logarithm of the energy decreases linearly until t " 250s and is constant for t ą 250s. This is due to the fact that the energy becomes smaller than 10 ´16 and thus smaller than the round-off error. Using a linear regression method, we found that the logarithm of the energy could be approximated, for t ă 250s by the function hptq " ´3 ´0.14t, so that the energy could be approximated by the function gptq " expp´0.14 ˚t ´3qq. These coefficients depends probably on the curvature Figure 16: Log(energy) vs time for γ " κ of the artificial boundary and of its distance to the obstacle. A more detailed numerical study will be necessary to analyze this dependence, but it is outside the scope of this paper. In Fig. 17, we compare the evolution of the discrete energy obtained with γ " κ (black curve) with the function g (blue dashed curve), the y scale is magnified by 100 in the second picture and by 10 8 in the third one. It is clear that g is always greater than the discrete energy which seems to indicate that the energy could be exponentially decreasing.

Conclusion

We have constructed a new family of ABCs depending on a parameter γ. When γ ě κ 4 , the ABCs satisfy each criterion defined in the introduction. In particular, we have proved that there exists a functional of time which is an energy when γ ě κ 4 and which is decreasing to zero. The numerical experiments that we have performed did not allow us to define that there exists a best ABC among the family of conditions. Nevertheless, our study let conclude that the C-ABC, which involves one-order differential operators, performs as well as the family of ABCs with γ ą κ 4 . Moreover, the C-ABC performs as well as the BGT2 condition and is long-term stable when coupled with IPDG method and Leap-Frog scheme, contrary to BGT2 condition. It would be now interesting to prove that the energy is actually exponentially decreasing. This result is already known for the C-ABC [START_REF] Bardos | Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary[END_REF][START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation de systèmes distribués, T.1: Contrôlabilité exacte[END_REF] and it might be satisfied for the full family of ABC with γ ě κ 4 .

A Implementation of the BGT2 condition in the IPDG formulation

We describe here how we obtain condition (28) with auxiliary variables from condition 26 and how we implement it in the IPDG formulation. The main issue is the term B s ˆp1`B t κpsq q ´1 2κpsq B s ˙u.

In the harmonic case, when B t is replaced by ´ik, this term can be easily taken into account into a variational formulation. This is unfortunately not possible for the transient wave equation, since the inverse of the operator 1 `Bt κpsq has to be computed via an auxiliary unknown. One solution could be to rewrite the above term as However, the first term would break the symmetry of the IPDG formulation and we prefer to approximate the curvature by a function which is constant on each edge. To do so, for each system reads as

M U n`1 ´2U n `Un´1 ∆t 2 `B U n`1 ´Un´1 ∆t ´DT ˆΨ1 n`1 `Ψ1 n´1 2 `Ψ2 n`2 `Ψ2 n´1 2 ˙`KU n " S n , C κ Ψ 1 n`1 ´Ψ1 n´1 ∆t `C Ψ 1 n`1 `Ψ1 n´1 2 `Kκ U n " 0, C Ψ n`1 2 ´2Ψ n 2 `Ψn´1 2 ∆t 2 `Un`1 `Un´1 2 " 0. ( 33
)
Note that this system is implicit on the interface. However, since matrices M , B, D, C κ are block-diagonal, we do not have to invert a global linear system but a linear system on each element of the boundary. However, the second penalty parameter hampers the CFL condition of the global system and we observed instabilities, even in the case of a circular boundary. To circumvent this problem, we decided to approximate the term K κ U by its mean value and we obtained the system

M U n`1 ´2U n `Un´1 ∆t 2 `B U n`1 ´Un´1 ∆t ´DT ˆΨ1 n`1 `Ψ1 n´1 2 `Ψ2 n`2 `Ψ2 n´1 2 ˙`KU n " S n , C κ Ψ 1 n`1 ´Ψ1 n´1 ∆t `C Ψ 1 n`1 `Ψ1 n´1 2 `Kκ U n`1 `Un´1 2 " 0, C Ψ n`1 2 ´2Ψ n 2 `Ψn´1 2 ∆t 2 `Un`1 `Un´1 2 " 0. (34)
Since the matrix K κ is not block-diagonal, we can not solve this system element by element. However, the linear system involves only the unknowns on the absorbing boundary.

żΣB

  n uϕdσ, where ϕ is a test-function. Hence, to include the boundary condition on Σ requires to multiply by the pseudo-differential operator ´Bt `γ ´κ 4

Figure 1 :

 1 Figure 1: Studied domain

Figure 2 : 1 Figure 3 :

 213 Figure 2: Computational domain -Configuration 1 Figure 3: Computational domain -Configuration 2

Figure 4 :

 4 Figure 4: Global relative error for different distances to the obstacle -Configuration 1

Figure 5 : 2 InriaFigure 6 : 1 Figure 7 : 2 Figure 8 :Figure 9 :

 52617289 Figure 5: Global relative error for different distances to the obstacle -Configuration 2

Figure 10 :

 10 Figure 10: Energy vs time for γ " κ

Figure 11 :

 11 Figure 11: Energy vs time for γ " 0.1κ

Figure 12 :Figure 14 :

 1214 Figure 12: Energy for γ " 0.249κ Figure 13: Energy for γ " 0.25κ

Figure 15 : 4 5. 4 . 3

 15443 Figure 15: Energy for γ " 0.25κ -y-scale magnified by 10 4

Figure 17 :

 17 Figure 17: Exponential decay of the energy

  ∆v1, v2 ``κ

													4	´γ˘ϕ˘.
	Then, we have
	pAv, vq "	ż Ω	∇v2 ¨∇v1 dx	`żΩ	∆v1v2 dx	`żΣ	κ 2	´v2	`´κ 4	´γ¯ϕ¯ϕ dσ.
	Using the Green formula (14), we get
				ż
	pAv, vq "	∇v2 ¨∇v1 dx	´żΩ	∇v2 ¨∇v1 dx	`żΣ	Bnv1v2 dσ
					Ω
				`xBnv1, v2y H ´1{2 pΓq,H 1{2 pΓq	`żΣ	κ 2	´v2	`´κ 4	´γ¯ϕ¯ϕ dσ.
	pAv, vq "	´żΣ	|v2| 2 dσ	´żΣ	κ 2	ϕv2 dσ	`żΣ	κ 2	ϕv2 dσ	`żΣ	κ 2 ´κ 4	´γ¯| ϕ| 2 dσ.
	Now, since κ 2 `κ 4	´γ˘ď 0 on Σ, we get that for all v in V
				pAv, vq "	´żΣ	|v2| 2 dσ	`żΣ	κ 2 ´κ 4	´γ¯| ϕ| 2 dσ ď 0,
	which completes the proof of Lemma 3.3.

Moreover, in V , we have v1 | Γ " v2 | Γ " 0 and on Σ, Bnv1 " ´v2 ´κ 2 ϕ. Hence,

  φy H ´1{2 pBΩq,H 1{2 pBΩq

						`żΣ	αv1φ dσ	`żΩ	v1φ dx "
		ż
		f φ dx	`żΣ	gφ dσ,
		Ω
						ż
	i.e., @φ P T	`Ω˘,	xBnv1, φy H ´1{2 pBΩq,H 1{2 pBΩq	`żΣ	αv1φ dσ "

Σ gφ dσ. Now, we have φ | Γ " 0, which implies that xBnv1 `˜1 `κ 2 `1 `γ ´κ 4 ˘¸v1 ´g, φy H ´1{2 pΣq,H 1{2 pΣq " 0, and we then have Bnv1 `˜1 `κ 2 `1 `γ ´κ 4 ˘¸v1 " g on Σ.

  Btu, ψq ď 0, if γ ě κ 4 and κ ě 0, and completes the proof of Lemma 4.1.

	Using the Green formula (14) and the relation B 2 t u " ∆u in Ω, we obtain
	d dt	Epu, Btu, ψq "	´żΩ	Btu∆u dx	`żBΩ	BnuBtu dσ	`żΩ	Btu∆u dx	`żΣ	κ 2	ψBtψ dσ.
	Since Btu | Γ " 0 and Bnu " ´Btu ´κ 2 ψ on Σ,
						d dt	Epu, Btu, ψq "	´żΣ	|Btu| 2 dσ	´żΣ	κ 2	ψBtudσ.
	At last, since Btψ " Btu ´`γ ´κ 4 ˘ψ, we get
		d dt	Epu, Btu, ψq "	´"ż	Σ	κ 2	´γ	´κ 4	¯|ψ| 2 dσ	`żΣ	|Btu| 2 dσ		.
	which implies that	d dt	Epu,
											ż Ω	∇pBtuq ¨∇u dx	`żΩ	BtuB 2 t u dx	`żΣ	κ 2	ψBtψ dσ.	(20)
													Inria

dt Epu, Btu, ψq "

Table 2 :

 2 Relative L 2 error ( in % ) -Source in (1.5,1) -Configuration 2

		p0, 2.85mq p0.7m, 2.75mq p1.4m, 2.45mq
	γ " κ 4 γ " κ 3		4.92 4.92	13.11 13.10	13.30 13.31
	γ " κ		4.89	13.19	13.37
	γ " 3κ		4.92	13.42	13.83
	γ " 10κ	5.33	14.04	14.22
	Table 1: Relative L 2 error ( in % ) -Source in (1.5,1) -Configuration 1
	p6m, 2.05mq p5m, 2.3mq p6m, ´2.05mq p5m, ´2.3mq
	γ " κ 4 γ " κ 3	2.49 2.49	3.65 3.65	6.77 6.79	11.53 11.56
	γ " κ	2.50	3.66	6.93	11.75
	γ " 3κ	2.52	3.67	7.35	12.31
	γ " 10κ	2.60	3.74	7.89	13.10

Table 3 :

 3 . We see that, when the artificial boundary is set Relative L 2 error ( in % ) -Source in (1.3,0) -Configuration 1 near the boundary of the obstacle, the relative errors are very important because the amplitude of reflected waves is high. When the radius of the artificial boundary is greater than 3m, we obviously obtain small relative errors because the reflected waves are generated by propagating waves impinging the external boundary with a low amplitude. We can observe that the error level stabilizes around 1% as soon as R ext ě 3m. Therefore, when the obstacle is a circle we recommend to take for the artificial boundary a circle whose radius is equal to three times the radius of the obstacle. Let us precise that the same conclusion holds for any value of γ greater than κ 4 . Now, we perform the same kind of analysis but looking at the global error in space and in Inria time, i.e. we evaluate the global relative L 2 pr0, ts ˆΩq-error defined by

		p1.025m, 1.025mq p1.256m, 0.725mq p1.4m, 0.375mq
	R ext " 1.5m	29.87	35.08	41.75
	R ext " 2m	6.49	10.84	18.17
	R ext " 3m	1.46	2.46	4.52
	R ext " 4m	0.96	1.42	2.03
	R ext " 5m	0.94	1.33	1.11
	R ext " 6m	0.93	0.87	1.07

Table 4 :

 4 Relative L 2 error ( in % ) -Source in (1.3,0) -Configuration 2

		p0.75m, 1.4mq p2.51m, 0.725mq p2.8m, 0.375mq
	a ext " 3m, b ext " 1.5m	35.03	61.61	69.34
	a ext " 4m, b ext " 2m	4.27	22.69	30.77
	a ext " 6m, b ext " 3m	4.1	8.75	10.52
	a ext " 8m, b ext " 4m	3.72	2.48	3.27
	a ext " 10m, b ext " 5m	0.59	1.15	1.77
	a ext " 12m, b ext " 6m	0.44	0.55	1.38

RR n°8203

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

Comparison with BGT2 condition

In this Subsection, we compare the performance of the curvature ABC with the performance of a BGT2 condition for arbitrarily convex surfaces. This condition reads as

B s ‹ 'u " 0, [START_REF] Trefethen | Well-posedness of one-way wave equations and absorbing boundary conditions[END_REF] where s denotes the curvilinear abscissa of the artificial boundary. We call it BGT2 condition since it is obtained from a classical BGT2 condition [START_REF] Antoine | Bayliss-turkel like radiation conditions on surfaces of arbitrary shape[END_REF] for the Helmholtz equation :

where k denotes the frequency of the monochromatic wave.

The time-domain formulation of this condition involves the inversions of time derivative operators. To avoid these inversions, we reformulate (26) by introducing two auxiliary unknowns defined only on the artificial boundary Σ :

We explain in the appendix how we derived this condition and how we implemented it in the IPDG formulation. Because of the second order derivative with respect to the curvature, the implementation of this condition requires the introduction of an additional penalty parameter on the absorbing boundary.

Note that, in the particular case where the boundary is a circle of radius R, the term B 2 s κ vanishes and we only need to compute the auxiliary variable ψ 1 .

In Fig. 6, we represent the global relative error obtained for the first configuration (the circle) for R ext " 1.5, 2, 3, 4, 5 and 6.

We remark that, compared to Fig. 4, the BGT2 condition improves the accuracy of the solution for R ext " 1.5 and 2. However, the errors are still greater than 10%. For R ext ě 3, there is no significant improvement. This means that, in order to obtain an accuracy around 1%, it is still necessary to choose R e xt ě 4, and the curvature condition should be preferred since it does not involve an auxiliary variable.

In Fig. 7, we represent the global relative error obtained for the second configuration (the ellipse) for pa ext , b ext q " p3, 1.5q, p4, 2q, p6, 3q, p8, 4q, p10, 5q and p12, 6q. In this configuration, the BGT2 condition improves the accuracy of the solution whatever the distance between the boundary and the obstacle is, even if this improvement is less important when the distance increases. However, the solution is long-term unstable in the case pa ext , b ext q " p3, 1.5q. To emphasize this fact, we have represented this solution on a larger time interval in Fig. 8. We have also represented the solution obtained with pa ext , b ext q " p4, 2q in Fig. 9 It does not necessarily mean that the continuous problem corresponding to the BGT2 condition is unstable. It may also mean that the numerical scheme that we employ may be not

Behavior of the discrete energy

In this Subsection, we illustrate the time behavior of the energy E defined at [START_REF] Hörmander | Pseudodifferential operators and hypoelliptic equations[END_REF]. To this aim, we first introduce the discrete energy corresponding to E.

For n P N, we set

We consider the Configuration 1. We first want to illustrate that the discrete energy is decreasing when the source is switched off. In a second part, we illustrate the fact that the condition γ ě κ 4 is an optimal condition to ensure the stability of the problem. In the last part, we show that the decay rate of the discrete energy seems to be exponential.

Energy decay

In Fig. 10, we represent the evolution of the discrete energy along the time when γ " κ. We see that the energy first increases, as long as the source is injecting some energy inside the system. Once Sptq " 0 (i.e. for t ě 5t 0 ), the energy is constant as long as the wave has not reached the boundary of the domain. Finally, the energy starts decreasing as predicted by Theorem 4.3.

Stability of the ABC

Theorem 4.3 states that the problem is long-term stable if γ is greater than κ 4 . We have performed some numerical experiments when the condition γ ě κ 4 . For example, if we take γ " 0.1κ, the system is unstable as we can see in Fig. 11. To illustrate the property that κ 4 is a critical value for γ, we compare the evolution of the discrete energy for γ " 0.249κ and γ " 0.25κ. In Fig. 12 (resp. in Fig. 13), we depict the evolution of the discrete energy during 1 000 000 iterations when γ " 0.249κ (resp. when γ " 0.25κ). If we give these figures a cursory glance, both systems seems to be stable for long-term, even when γ ă κ 4 . But if we look at Fig. 14 (resp. Fig. 15) where we have magnified the y-scale by a factor 10 4 , for γ " 0.249κ (resp. when γ " 0.25κ), we see that the scheme is not stable when γ ă κ 4 . The numerical tests thus indicate that κ 4 could be critical value for γ.

Barucq & Diaz & Duprat

edge σ on the boundary, we denote by s 1 σ and s 2 σ the curvilinear coordinates of the two nodes of σ and we define κ σ " pκps 1 σ q `κps 2 σ qq{2. This simplification will also be useful for the numerical computation of the integrals on the boundary. In the same way, we define B 2 s κ σ " pB 2 s κps 1 σ q `B2 s κps 2 σ qq{2 and the approximated BGT2 condition on an external edge σ reads as

Now, following the methodology we used above, we introduce two auxiliary variables ψ 1 and ψ 2 and we rewrite the condition as

Note that when the boundary is a circle, there is no need to introduce the second auxiliary variable ψ 2 . To discretize the operator B 2 s u by the IPDG method, we introduce the bilinear form

where rruss s 1 σ and ttuuu s 1 σ denote respectively the jump and the mean value of u at point xps 1 σ q. We obtain the algebraic system .

and the other matrices have been defined above.

We first chose to use a classical Leap-Frog scheme in the whole domain and the fully discretized Inria


Marc-Antoine Bouzigues
email: marcantoine.bouzigues@orange.com

Isabelle Siaud
email: isabelle.siaud@orange.com

Maryline Hélard
email: maryline.helard@insa-rennes.fr

Anne-Marie Ulmer-Moll
email: annemarie.ulmermoll@orange.com

Anne-Marie Ulmer-Moll Time

Member IEEE A.-M Ulmer-Moll

Time Reversal for Green Radio Communications

Keywords: Time Reversal, Green Communications, Fast Session Transfer, Cooperative Multi-Point

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

OVER the past few years, the carbon footprint of Information and Communication Technologies (ICTs) has grown. From 2% of global carbon emissions in 2007 it may reach 15% in 2025 as traffic demand grows exponentially while energy efficiency of network infrastructures is linearly increasing due to hardware natural improvement [START_REF] Gati | Green mobile access networks[END_REF]. Therefore, interest for energy-aware (green) techniques is growing and Time Reversal (TR) has been identified as a good candidate for green communications [START_REF] Wang | Green wireless communications: a time-reversal paradigm[END_REF].

An evidence of the burgeoning interest for green communications is the growth of energy-aware consortium and international projects such as GreenTouch [START_REF] Kilper | GreenTouch Consortium: Building the Roadmap[END_REF], and EARTH (Energy Aware Radio and neTwork tecHnologies) [START_REF] Bohn | Most Promising Tracks of Green Radio Technologies[END_REF]. GreenTouch is a consortium of international ICTs actors. Its mission is to deliver the architecture, specifications and roadmap needed to increase network energy efficiency by a factor 1000 from current levels. Similarly, international project EARTH goal is to significantly reduce the power consumption of cellular networks delivering different approaches applicable on component, node, and link levels.

In order to determine if TR is an efficient green wireless communication technique, theoretical gains are faced to link level simulations performed in relevant green scenarios. Dedicated use cases have been then defined and mapped on the novel realistic green scenarios of multi-Radio Access Technology (RAT) [START_REF] Siaud | EE multi-RAT mobile backhauling[END_REF] that are promoted by French operator Orange in the GreenTouch consortium. The first scenario is the Long Term Evolution-Advanced Coordinated Multi-Point transmission/reception (LTE-Advanced CoMP) scheme [START_REF] Taoka | MIMO and CoMP in LTE-Advanced[END_REF], [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA); Further advancements for E-UTRA Physical layer aspects[END_REF] which aims to reduce border-cell interference. The second scenario is an extension [START_REF] Cariou | Fast Session Transfer[END_REF], [9] of the IEEE 802.1X Fast Session Transfer (FST) mechanism [START_REF]Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications: Enhancements for Very High Throughput in the 60GHz Band[END_REF] which allows stations to transfer a session from a channel to another channel in the same or different frequency bands (2.4/5 GHz and 60 GHz) by considering different WLAN/WPAN air interfaces benefiting from Wi-Fi network backward compatibility. The third scenario is a multi-RAT heterogeneous network aiming to provide a seamless connectivity under power and Quality of Service (QoS) constraints considering a switching between independent networks using a Multiple Interface Management (MIM) utilizing inter-Medium Access Control (MAC) layer and dedicated multi-RAT Channel Quality Indicator (CQI) metrics. Such studies and results have been performed in the ICT-OMEGA (hOME Gigabit Access) project for the Home Networking 1 [START_REF] Javaudin | OMEGA FP7 project[END_REF]1 [START_REF] Caldera | Seventh Framework programme theme 3: Information & Communication Technologies (ICT)[END_REF].

The aim of this paper is to describe and evaluate TR benefits in the green communications context. The paper is organized as follows. Section II introduces TR principle and properties. Section III exposes novel uses cases for green communications turned towards convergent and heterogeneous networks. TR generic use cases are defined in section IV. Section V gives a propagation analysis carried out to predict advantages and limits of RT signal processing. A conclusion is drawn in section VI.

II. TIME REVERSAL PRINCIPLE

TR is a signal processing technique adapted from the acoustic field [START_REF] Fink | Time reversal of ultrasonic fields -part i basic principles[END_REF] toward electromagnetic field [START_REF] Lerosey | Time Reversal of electromagnetic waves[END_REF]. It has been applied to Ultra-Wideband (UWB) systems [START_REF] Kaleghi | Demonstration of time-reversal in indoor ultra-wideband communication: time domain measurement[END_REF] and, more recently, to narrow band communications systems based on non-impulsive waveforms [START_REF] Dubois | On the use of timereversal for digital communications with non-impulsive waveforms[END_REF]. A proof of the growing interest for TR applied to green communications is the project TRIMARAN: Time Reversal MIMO (Multiple-Input Multiple-Output) OFDM (Orthogonal Frequency Division Multiplexing) Green communications based on Micro Structured Antenna Arrays [START_REF] Phan Huy | TRIMARAN 2011 -Time Reversal MIMO OFDM Green communications base don Micro Structured Antenna Arrays[END_REF].

TR can be considered as a special case of beamforming pre-filtering technique and enables to focus waves in the temporal and spatial domain. TR spatial focusing property can be used to reduce interferences [START_REF] Hieu | On the MSI mitigation for MIMO UWB time reversal systems[END_REF] and its time compression property simplifies the receiver complexity [START_REF] Oestges | Time Reversal Techniques for Broadband Wireless Communication Systems[END_REF]- [START_REF] Strohmer | Application of Time-Reversal with MMSE Equalizer to UWB Communications[END_REF], thus, using TR allows reducing both power emission and power consumption of systems.

The TR principle lies in two steps: the sounding (fig. 1) and the focusing (fig. 2). The sounding step consists in an estimation of the impulse response of the multipath channel from the receiver to the transmitter. Then, the transmitter time reverses the Channel Impulse Response (CIR) and uses the result to pre-filter the signal. Electromagnetic waves reexperience the paths and add coherently at the target location.

Let us denote 1234 the transmitted signal, 5234 the received signal, 6237 84 the CIR where 8 is the delay variable, 9234 the Additive White Gaussian Noise (AWGN),

The multipath propagation channel is expressed as a sum of echoes:

6237 84 D F 637 8 A 234 8 8 A 234 A (2)
The complex coefficient 6237 8 A 2344 is the gain the B echo of the baseband CIR of the propagation channel at the time3. When TR is used, the signal is pre-filter by

The equivalent channel has then a peak which corresponds to the center of the autocorrelation function. Width of this central peak can be used to estimates the quality of the TR focalization properties.

III. MULTI-RAT AND CONVERGENT BASED USE CASES FOR GREEN COMMUNICATIONS

The aim of green communications is to increase the energy-efficiency of ICTs. One way to achieve this goal is to reduce the complexity and transmission power of systems under QoS constraint, a convergent system should be able to switch toward the greenest technology available according to quality metrics. Switching mechanisms between different WLAN/WPAN interfaces and associated quality metrics to perform Multiple Interface Management are designed and evaluated in the project OMEGA (hOME Gigabit Access) 1 [START_REF] Caldera | Seventh Framework programme theme 3: Information & Communication Technologies (ICT)[END_REF] .Green communications uses cases extended to mobile radio communications and convergent indoor/mobile networks have been designed by Orange Labs [START_REF] Siaud | EE multi-RAT mobile backhauling[END_REF] in the Green Touch consortium. These use cases are described below.

A. LTE-Advanced CoMP

Inter-cell interference is a recurring problem for wireless communications. LTE-Advanced CoMP is a radio-access technology intended to reduce border-cell interference: The station (eNode B) and Remote Radio Equipments (RREs) cooperate together on the downlink transmission in order to mitigate interference. It can be divided in two categories: Coordinated Scheduling / Coordinated Beamforming (CS/CB) and joint processing [START_REF] Karakayali | Network coordination for spectrally efficient communications in cellular systems[END_REF], [START_REF] Andrews | Overcoming interference in spatial multiplexing MIMO Cellular Networks[END_REF] (fig. 3). Joint processing means that several eNode-B/RREs address the user (either with joint transmission or instantaneous cell selection) while with CS/CB a given station focus on the user while other stations try to null their interferences.

When beamforming is used the main lobe focus on the user, but secondary lobes can create interferences toward other users of the same cell. CS/CB intends to deal with this interference but it may be hard to perform perfect CS/CB with a high number of users. With TR, the focalization is not made toward the user location but at the user location and without secondary lobes. Therefore, the main benefit of Time Reversal applied to CoMP is to efficiently achieve spatial focalization without interferences toward other users.

Joint processing may benefit from the ability of TR to address a user even in a highly scattering medium and in Non-Line of Sight (NLOS) situation.

B. Fast Session Transfer

Fast Session Transfer (FST) is a mechanism that allows a station to transfer a session from a channel to another channel in the same or different frequency band [START_REF]Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications: Enhancements for Very High Throughput in the 60GHz Band[END_REF] by considering several Wi-Fi certified air interfaces. It can be applied to IEEE 802.1X systems in OBand (e.g., 2.4 and 5 GHz) and DBand (e.g., 60 GHz) if stations are multi-band capable. Stations can become noticeable, or not, for other stations by changing carrier frequency, FST can then be use to perform multi-hop relay and backhauling. This mechanism was first described within the IEEE 802.11ad standard [START_REF]Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications: Enhancements for Very High Throughput in the 60GHz Band[END_REF] but compatibility with other IEEE802.1X systems (e.g., 802.15 or 802.16) is intended.

Transition from one band/channel toward another band/channel is done according to the protocol defined in fig. 4. From the initial state, the stations exchange FST setup request and response. Once the Link Loss Timeout (LLT) reaches the value 0, stations operate a switch toward a new bandwidth/channel. If the switching if confirmed by acknowledgments (Ack), the transfer is completed. Fig. 5 shows an example of its utilization to target and lock a station at 60 GHz thanks to the 5 GHz range.

As TR performance depends on the propagation channel, transmission bandwidth size and Radio-Frequency (RF) carrier, selection related to FST mechanisms may be efficiently combined with TR processing.

C. Simple and multi-RAT heterogeneous network

The last use case includes the previous ones. It consists in a simple and multi-RAT heterogeneous network. Given Channel Quality Indicators (CQIs), the network must be able to select the greenest access combination (access points, technology and modulation and coding scheme) under QoS constraints. Therefore, it must include CQIs evaluation mechanisms and switching protocols enabling seamless connectivity as investigated in the ICT OMEGA project.

LTE-Advanced CoMP can be viewed as a simple-RAT heterogeneous network where the best scheme of stations combination is selected to increase the Signal to Interference plus Noise Ratio (SINR) at the user location. FST is a switching protocol combining RF multi-band processing to MIM between Wi-Fi certified air interfaces.

IV. GENERIC TR USE CASES

This section is dedicated to the identification of limits and advantages of TR by defining TR use cases to be applied to multi-RAT and convergent use cases (cf. III). Impact of the frequency selectivity and bandwidth size, Doppler effect impact, antenna mismatch alignment, synchronization problems and number of transmission antennas are studied.

A. Frequency selectivity and Bandwidth size influence

The ideal TR equivalent CIR is the Dirac delta function. To strive toward this ideal case, the digital channel needs to have sufficiently uncorrelated taps. Additional taps can be seen as additional information about the user location and the more information we have about the user location, the more we can focus on it. An easy way of getting this extra information is to increase the transmission bandwidth as it increases the resolution of the channel, i.e., the number of taps we can detect. Another way is to transmit over a highly scattering medium. Root Mean Square (RMS) delay spread8 is a common metric that quantifies the multipath richness of the medium. For a digital channel it can be expressed as:

8 D ! " 28 A 8 4 # $ % A A " % A A & '(# ( 5
)
where % A is the average power of echo B and 8 is the average delay of the channel given by:

8 D " 8 A $ % A A " % A A (6)
On fig. 6 curves show that, for a given channel, focalization increases with transmission bandwidth. Moreover, focalization (i.e. first drop of the delay profile) is better with channel E which has a RMS delay spread of 100ns than with channel C which has a RMS delay spread of 30ns

Wide bandwidth and/or high delay spread (hence, high frequency selectivity) are then prerequisites to obtain good performances with TR as shown on fig. 6.

B. Doppler effect impact

One of the main drawbacks of TR is the need for slow channel variations. In order to get good performances, the channel must remain constant between the sounding and focusing step. If the estimated channel used for pre-filtering is outdated, the instantaneous correlation function becomes a cross-correlation instead of an autocorrelation and focalization properties are lost. A high Doppler effect will then induce strong uncorrelation between the estimated channel and the effective channel, causing a drop of performances.

As shown by fig. 7, if the system can estimates properly the CIR or if the channel variations are too fast, it is better not to use TR. Uncorrelated Time-Reversal means that the estimated CIR used for pre-filtering is outdated. Synchronization is made at the central peak location but as the CIR is outdated the main tap is no more at the peak location, the equivalent channel seems in NLOS propagation.

C. Antenna mismatch alignment

Systems using directional antennas can suffer from mismatching, i.e. an unwanted transition from LOS to NLOS. With TR, this critical drop of performances can be avoided if FST is used to switch to a wider channel or a different carrier (with different propagation properties) when SNR drops.

D. Synchronisation problems

Ideally, systems using TR should synchronize on the central peak of the equivalent channel, thus, systems should be able to deal with the previous taps of the equivalent channel (e.g., to avoid interferences between consecutive OFDM symbols). Moreover, if the synchronization is biased, it can cause a drop of performance when the central peak of the equivalent channel is very sharp.

E. Number of antennas

In order to increase TR performances, it is also possible to increase the number of transmit antennas and to apply TR before each. The resulting equivalent channel is then the sum of the autocorrelation functions of the different CIRs between each transmit antennas and the receive antenna. On fig. 8, curves shows that increasing the number of transmit antennas increases the focalization. The difference between the two first taps is increase but the remaining of the power delay profile is unchanged. Then, TR benefits from spatial diversity if channels are enough uncorrelated.

V. PROPAGATION ANALYSIS

In order to evaluate benefits of TR, a comparison of TR performances is done through a propagation analysis of equivalent channels. We used IEEE 802.11 High Throughput (HT) and Very High Throughput (VHT) transmission modes of IEEE 802.11 standards [START_REF]Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications: Enhancements for Very High Throughput in the 60GHz Band[END_REF], [START_REF]Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications[END_REF] over TGn propagation channel models C and E [START_REF] Erceg | TGn Channel Models[END_REF]. These models are valid for simulations with a communication bandwidth up to 40MHz. As we intend to simulate bandwidths up to 80MHz, we linearly interpolated TGn channels C and E. RMS delay spreads of channels stay unchanged: 30ns for channel C and 100ns for channel E.

As we can see on fig. 9, thanks to its focalization property which is in enhanced by bandwidth size, TR helps to increase the difference between the 2 main taps when the receiver synchronizes on the central peak. The main tap of channel E is not the first tap (as it is the case in NLOS situations), but when TR is used, if the receiver synchronizes on the central peak the multipath effect of the channel E is softened. Improvements are better with this channel which corresponds to a more scattering medium. This analysis corroborates the conclusions of part IV A: TR improvements are greater in a highly scattering medium and with a large transmission bandwidth.

VI. CONCLUSION

In this article we carry out a first analysis of benefits and limits of TR for green communications applied to Wi-Fi certified technologies. Practical Green use cases derived from Orange Labs contribution to GreenTouch were first exposed. Then we defined generic TR uses cases. Such cases can be mapped on the practical Green use cases in order to determine the benefits and limits of TR for green communication. A propagation analysis was also performed upon channels used upon WLAN/WPAN systems operating in 2.4/5 GHz bands with different channel bandwidth sizes, typically ranged from 20 MHz to 160 MHz.

The focalization property of TR allows interference reduction for LTE-Advanced CoMP and can simplify the addressing of selected user because the cancellation of interferences toward other users is inherent to TR. Moreover, its compatibility with highly scattering media can be use with FST to ensure a seamless connectivity even in case of antennas mismatching. TR improves systems robustness while reducing complexity of receiver thanks to its time compression property. Green communications can therefore highly benefit from TR. But we should however be careful with its utilization when transmission bandwidth is limited, or when the medium is not enough scattering. Moreover, if applying TR with a perfect knowledge of the channel improves performances of systems, a bad or outdated estimation of the channel can severely drop performances.

8 A

 8 the delay of path B and C the convolution operation.Without TR, the received signal is expressed as follows:

6 C

 6 237 8 34, where * denominates the conjugate. Then, the received signal becomes: 5234 D 1234 C 237 84 E 9234 (3) where 237 84 is the instantaneous correlation function of the channel associated with the equivalent TR channel 6 237 84 expressed as follow: 6 237 84 D 6 C 237 8 84 C 6237 84 D 237 84

Fig. 1 .

 1 Fig. 1. Sounding step: the user sends a probe signal used by the base station to estimates the CIR.

Fig. 2 .

 2 Fig. 2. Focusing step: the base station pre-filters the signal with a timereversed version of the CIR. From the user side, the equivalent channel is then the autocorrelation of the CIR.

Fig. 4 .

 4 Fig. 4. States of the FST setup protocol.

Fig. 5 .

 5 Fig. 5. An example of the use of FST: the station (STA) is located by the access point (AP) in the 5 GHz range, its location is then used to perform beamforming (BF) and reach it with 60 GHz technology.

Fig. 3 .

 3 Fig. 3. Downlink CoMP transmission categories adapted from [6]. With CS/CB, stations null their interference toward other users while addressing desired user. With joint processing user can either use joint coherent/noncoherent transmission or make instantaneous cell selection.

Fig. 8 .

 8 Fig. 8. Average Power Delay profiles (in dB) of channel C with 20 MHz and 80 MHz bandwidths for 1,2 and 4 transmit antennas and 1 receive antenna.

Fig. 7 .

 7 Fig. 7. Average Power Delay profile (in dB) of channel C with 20 MHz bandwidth.

Fig. 6 .

 6 Fig. 6. First 100ns of the average delay profiles (in dB) of TGn channel C and E [25] with receiver synchronization on the central peak of the equivalent CIR.