Jean-Luc Akian
email: jean-luc.akian@onera.fr

Radjesvarane Alexandre
email: alexandreradja@gmail.com

Salma Bougacha
email: salma.bougacha@ecp.fr

A Gaussian beam approach for computing Wigner measures in convex domains

Keywords: Mathematics Subject Classification: 35L05, 35L20, 81S30 wave equation, Gaussian beam summation, Wigner measures, FBI transform, reflection

A Gaussian beam method is presented for the analysis of the energy of the high frequency solution to the mixed problem of the scalar wave equation in an open and convex subset Ω of R n , with initial conditions compactly supported in Ω, and Dirichlet or Neumann type boundary condition. The transport of the microlocal energy density along the broken bicharacteristic flow at the high frequency limit is proved through the use of Wigner measures. Our approach consists first in computing explicitly the Wigner measures under an additional control of the initial data allowing to approach the solution by a superposition of first order Gaussian beams. The results are then generalized to standard initial conditions.

Introduction

We are interested in the high frequency limit of the initial-boundary value problem (IBVP) for the wave equation

P u ε = ∂ 2 t u ε - n j=1 ∂ xj c 2 (x)∂ xj u ε = 0 in [0, T ] × Ω, (1.1a
)

Bu ε = 0 in [0, T ] × ∂Ω, (1.1b)
u ε | t=0 = u I ε , ∂ t u ε | t=0 = v I ε in Ω, (1.1c)
where B stands for a Dirichlet or Neumann type boundary operator. Above, T > 0 is fixed, Ω is a bounded domain of R n with a C ∞ boundary and the wave propagation velocity c is in C ∞ ( Ω), though this assumption may be relaxed.

The initial data depend on a small wavelength parameter ε > 0 and we assume that u I ε and v I ε are uniformly bounded w.r.t. ε respectively in H 1 (Ω) and L 2 (Ω).

(H1)

We are interested in the description of the behavior of the local energy density

1 2 |∂ t u ε | 2 + 1 2 n j=1 c 2 |∂ xj u ε | 2
, at the high frequency limit ε → 0, in which case, it is well known that this quantity can be computed through the use of Wigner measures. The Wigner transform is a phase space distribution introduced by E. Wigner [START_REF] Wigner | On the quantum correction for thermodynamic equilibrium[END_REF] in 1932 to study quantum corrections to classical statistical mechanics. In the 90's, mathematicians became increasingly interested by these transforms and related measures, see for example [START_REF] Lions | Sur les mesures de Wigner[END_REF][START_REF] Markowich | The classical limit of a self-consistent Quantum-Vlasov equation in 3D[END_REF][START_REF] Markowich | A Wigner-function approach to (semi)classical limits: electrons in a periodic potential[END_REF][START_REF] Markowich | Weak limits of finite difference schemes of Schrödinger-type equations[END_REF] for the semiclassical limit of Schrödinger equations. A general theory for their use in the homogenization of energy densities of dispersive equations was laid out by Gérard et al. in [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF], see also [START_REF] Gérard | Mesures semi-classiques et ondes de Bloch[END_REF][START_REF] Gasser | Quantum hydrodynamics, Wigner transform and the classical limit[END_REF]. Wigner measures are also related to the H-measures and microlocal defect measures introduced in [START_REF] Tartar | H-measures, a new approach for studying homogenization, oscillations and concentration effects in partial differential equations[END_REF] and [START_REF] Gérard | Microlocal defect measures[END_REF], see also [START_REF] Burq | Mesures semi-classiques et mesures de défaut, (French) [Semiclassical measures and defect measures[END_REF][START_REF] Bougacha | Oscillations in PDE with singularities of codimension one. Part I : review of the symbolic calculus and basic definitions[END_REF]. Whereas there is no notion of scale for the latter measures, Wigner transforms are associated to a small parameter tending to zero. In quantum mechanics, this parameter is the rescaled Planck constant, while it will be typically the distance between two points of the medium's periodic structure for homogenization problems.

The Wigner transform, at the scale ε, is defined for a given sequence (a ε , b ε ) in S ′ (R n ) p × S ′ (R n ) p as the tempered distribution

w ε (a ε , b ε )(x, ξ) = (2π) -n R n e -iv•ξ a ε (x + ε 2 v)b * ε (x - ε 2 v)dv.
If a ε is uniformly bounded w.r.t. ε in L 2 (R n ) p , then w ε [a ε ] := w ε (a ε , a ε ) converges as ε goes to 0 in M p S ′ (R n x × R n ξ ) to a positive hermitian matrix measure (modulo the extraction of a subsequence), which is called a Wigner measure associated to (a ε ) and denoted w[a ε ]. The Wigner measures associated to the solution of the wave equation (and hyperbolic problems in general, see e.g. [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF][START_REF] Papanicolaou | Waves and Transport, in "Hyperbolic Equations and Frequency Interactions[END_REF]) are related to the energy density in the high frequency limit. More precisely, under suitable hypotheses (see Proposition 1.7 in [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF]), the density of energy associated to the solution u C ε of the Cauchy problem for the scalar wave equation converges as ε → 0 in the sense of measures to

R n E u C ε (t, .) (x, dξ),
where

E u C ε (t, .) = 1 2 w[∂ t u C ε (t, .)] + 1 2 n j=1 w[c∂ xj u C ε (t, .)].
Moreover, E u C ε (t, .) is the sum of two measures satisfying transport equations of Liouville type (see e.g. [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF]).

For the Dirichlet or Neumann initial boundary value problem connected with the wave equation, we shall study the same quantity after extending ∂ t u ε (t, .) and c∂ xj u ε (t, .), j = 1, . . . , n, to functions of L 2 (R n ) by setting

∂ t u ε = 1 Ω ∂ t u ε , ∂ xj u ε =
1 Ω ∂ xj u ε and extending c outside Ω in a smooth way. Hence, we call microlocal energy density of u ε the distribution

1 2 w ε ∂ t u ε (t, .) + 1 2 n j=1
w ε c∂ xj u ε (t, .)

and its high frequency limit the measure

E (u ε (t, .)) = 1 2 w ∂ t u ε (t, .) + 1 2 n j=1
w c∂ xj u ε (t, .) . (H3)

v I ε = 1 Ω v I ε and ∂ xj u I ε = 1 Ω ∂ xj u I ε (j =
Our present study will be restricted to the case where the rays starting from the support of the initial data do not face diffraction on the boundary, nor do they glide along ∂Ω. Therefore, we also assume that u I ε and v I ε have supports contained in a fixed compact set of Ω independent of ε, (H4)

Ω is convex with respect to the bicharacteristics of the wave operator, that is every ray originating from Ω hits the boundary twice and transversally, and the boundary has no dead-end trajectories, that is infinite number of successive reflections cannot occur in a finite time.

These geometric hypotheses insure that the only phenomena occurring at the boundary is the reflection according to the geometrical optics laws.

Wigner measures for the wave equation in presence of a boundary or an interface have been studied by Miller [START_REF] Miller | Refraction of high-frequency waves density by sharp interfaces and semiclassical measures at the boundary[END_REF] who proved refraction results for sharp interfaces and Burq [START_REF] Burq | Contrôlabilité exacte des ondes dans des ouverts peu réguliers, (French) [Exact controllability of waves in nonsmooth domains[END_REF] who described their support for a Dirichlet boundary condition. Similar results have been established for other problems [START_REF] Burq | Mesures de défaut de compacité, application au système de Lamé, (French) [Microlocal defect measures and application to the Lame system[END_REF][START_REF] Duyckaerts | Stabilization of the linear system of magnetoelasticity[END_REF][START_REF] Fouassier | High frequency limit of Helmholtz equations: refraction by sharp interfaces[END_REF], in particular the eigenfunctions for the Dirichlet problem [START_REF] Zelditch | Uniform distribution of eigenfunctions on compact hyperbolic surfaces[END_REF][START_REF] Gérard | Ergodic properties of eigenfunctions for the Dirichlet problem[END_REF] and for the Neumann and Robin problems [START_REF] Burq | Quantum ergodicity of boundary values of eigenfunctions: a control theory approach[END_REF]. All these works are based on pseudo-differential calculus, and in particular the use of a tangential pseudo differential calculus.

In this paper, we present an approach to compute Wigner measures based on the Gaussian beam formalism. Therefore, we avoid any use of adapted pseudodifferential calculus. Though a Gaussian beam technic requires much more work, compared to the above mentioned papers, one advantage is that we are able to give asymptotic estimates for remainders terms, which could be useful for numerical purposes for instance.

Let us recall that Gaussian beams (or the related coherent states) are waves with a Gaussian shape at any instant, localized near a single ray [START_REF] Babič | Eigenfunctions concentrated in a neighborhood of a closed geodesic[END_REF][START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF]. They play the role of a basis of fundamental solutions of wave motion and furthermore can be used to study general solutions of partial differential equations (PDEs). For example, they can help for the understanging of propagation of singularities [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF], to prove lack of observability [START_REF] Macià | On the lack of observability for wave equations: a Gaussian beam approach[END_REF] and to study semiclassical measures [START_REF] Paul | On the pointwise behavior of semi-classical measures[END_REF] and trace formulas [START_REF] Wilkinson | A semiclassical sum rule for matrix elements of classically chaotic systems[END_REF][START_REF] Combescure | A proof of the Gutzwiller semiclassical trace formula using coherent states decomposition[END_REF].

To describe non localized solutions of PDEs, one can use the Gaussian beam summation method [START_REF] Katchalov | The application of the Gaussian beam summation method to the computation of high-frequency wave fields[END_REF][START_REF] Červený | Computation of wave fields in inhomogeneous media -Gaussian beam approach[END_REF][START_REF] Klimeš | Expansion of a high-frequency time-harmonic wavefield given on an initial surface into Gaussian beams[END_REF]. The initial field is expanded as a sum of Gaussian beams. Each individual beam is computed and the solution is then obtained at an observation point by superposing a selection of Gaussian beams. The summation strategies are numerous. The sum can be discrete [START_REF] Motamed | Taylor expansion and discretization errors in Gaussian beam superposition[END_REF][START_REF] Tanushev | Gaussian beam decomposition of high frequency wave fields[END_REF][START_REF] Ariel | Gaussian beam decomposition of high frequency wave fields using expectation-maximization[END_REF] or continuous [START_REF] Liu | Recovery of high frequency wave fields for the acoustic wave equation[END_REF][START_REF] Liu | Error Estimates for Gaussian Beam Superpositions[END_REF], the selection of the beams to be superposed can be done according to several criteria. In [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF], a weighted integral of Gaussian beams was designed to build an approximate solution of the IBVP (1.1) under an additional assumption (H5) on the initial data (see p.9). See also [START_REF] Leung | Eulerian Gaussian beams for Schrödinger equations in the semi-classical regime[END_REF][START_REF] Leung | The backward phase flow and FBI-transform-based Eulerian Gaussian beams for the Schrödinger equation[END_REF][START_REF] Qian | Fast multiscale Gaussian wavepacket transforms and multiscale Gaussian beams for the wave equation[END_REF] for recent numerical implementations related to this method.

Gaussian beams seem to be very well suited for the study of Wigner measures. Indeed, the Wigner transform of two different beams vanishes when ε goes to zero. Even better, the Wigner measure of one individual beam is a Dirac mass localized on the corresponding bicharacteristic. Thus Gaussian beams act as an orthogonal family for the Wigner measure. Using these elementary solutions for studying Wigner measures is not new, see for example in the whole space domain the work by Robinson [START_REF] Robinson | Semiclassical mechanics for time-dependent Wigner functions[END_REF] for the Schrödinger equation, and more recently the paper by Castella [START_REF] Castella | The radiation condition at infinity for the high frequency Helmholtz equation with source term: a wave packet approach[END_REF] who used a coherent states approach for the Helmholtz equation.

As the microlocal energy density of one individual beam is concentrated near its associated bicharacteristic, one would expect that the Wigner measure of a summation of weighted Gaussian beams will yield easily that the associated weights are transported along the broken bicharacteristic flow (see p.8 for the construction of reflected flows and p.29 for the definition of the broken flow). Unfortunately this result is not immediate as even different beams become infinitely close to each other. However, we shall show by elementary computations that this intuition is indeed true and that the microlocal energy density of the considered approximate solution is transported at the high frequency limit along the broken bicharacterisitic flow. Since the asymptotic solution is close to the exact solution u ε , we may deduce the same consequence for E (u ε (t, .)).

The additional hypothesis (H5) consists in assuming that in the frequency space, the initial data are supported in a compact that does not contain 0 (modulo infinitely small residues). When studying Wigner measures by the pseudo-differential calculus techniques, the frequency behavior of the initial conditions is only controlled by the less restrictive hypotheses (H2) and (H3). Hence, the assumption (H5) is artificial (though not for numerical purposes) and is required only by the Gaussian beam summation method we have chosen. However, for ε-oscillatory initial data with Wigner measures not charging the set R n × {ξ = 0}, a truncation of their frequency support at infinity and at zero does not affect the energy density of the solution as ε → 0. By achieving such a truncation, we succeed to derive the transport property of the energy density under the traditional hypotheses (H2) and (H3): Theorem 1.1 Assume the hypotheses (H1)-(H4) on the initial conditions hold true. Let E ± = 1 2 w v I ε ± ic|D|u I ε and denote by ϕ t b the broken bicharacteristic flow associated to -i∂ t -c|D| obtained after successive reflections on the boundary ∂Ω.

Then E (u ε (t, .)) = 1 2 E + o (ϕ -t b ) -1 + E -o (ϕ t b ) -1 in Ω × (R n \{0}).
As mentioned already in our Introduction, this result is well known. But our method of proof is able to give more precise estimations than those stated above for the Wigner measures. In particular, we have estimations on the Wigner transforms of the solutions.

The rest of the paper is organized as follows. In Section 2, we recall the construction of first order Gaussian beams and the structure of the asymptotic solutions obtained as an infinite sum of such beams. The derivatives of the asymptotic solutions are then expressed using Gaussian type integrals. We simplify the expression of the Wigner transform of such integrals in Section 3, following initial computations of [START_REF] Robinson | Semiclassical mechanics for time-dependent Wigner functions[END_REF] in the Schrödinger case. We then compute the microlocal energy density of the asymptotic solution by exploiting the expressions of the beams phases and amplitudes and using the dominated convergence theorem. We prove the propagation along the broken flow of E (u ε (t, .) at the high frequency limit, with the help of assumptions (H2) and (H3) on the initial data. Some complementary results are collected in an Appendix, Section 4.

Let us end this Introduction with a few notations which will be used hereafter.

A vector x ∈ R d will be denoted by (x 1 , . . . , x d ), the inner product of two vectors a, b ∈ R d by a • b, and the transpose of a matrix A by A T . If E is a subset of R d , we denote E c its complementary and 1 E its characteristic function. For a function f ∈ L 2 (Ω), we let f = 1 Ω f . For r > 0, χ r denotes a cut-off function in

C ∞ 0 (R n , [0, 1]) such that χ r (x) = 1 if |x| ≤ r/2 and χ r (x) = 0 if |x| ≥ r.
We use the following definition of the Fourier transform

F x u(ξ) = R d u(x)e -ix•ξ dx for u ∈ L 2 (R d ).
If no confusion is possible, we shall omit the reference to the lower index x.

We keep the standard multi-index notations. For a scalar function

f ∈ C ∞ (R d x , C), ∂ x f will denote its gradient vector (∂ xj f ) 1≤j≤d and ∂ 2 x f will denote its Hessian matrix (∂ xj ∂ x k f ) 1≤j,k≤d . For a vector function g ∈ C ∞ (R d , C p ), the notation Dg is used for its Jacobian matrix (Dg) j,k = ∂ x k g j . If g is a function in C ∞ (R n y × R n η , C p ), we denote (D y g) j,k = ∂ y k g j and (D η g) j,k = ∂ η k g j .
We use the letter C to denote a (possible different at each occurence) positive constant. For (y ε ) and (z ε ) sequences of R + with ε ∈]0, ε 0 ], we use the notation y ε z ε if there exists C > 0 independent of ε such that y ε ≤ Cz ε for ε small enough. We write

y ε ε ∞ or y ε = O(ε ∞ ) if for any s ≥ 0 there exists C s > 0 s.t. for ε small enough y ε ≤ C s ε s . Finally, if E is in an open subset of R 2n and ν ε , ν ′ ε are two distributions s.t. lim ε→0 (ν ε -ν ′ ε ) = 0 in E, we shall write ν ε ≈ ν ′ ε in E.
2 Tool-box and construction of the asymptotic solution

We recall the construction made in [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF] of an asymptotic solution as a superposition of Gaussian beams and give the expression of its time and spatial derivatives with the help of so called Gaussian integrals.

First order Gaussian beams

Beams in the whole space

Let h + (x, ξ) = c(x)|ξ| and (x t , ξ t ) be a Hamiltonian flow for h + , that is a solution of the system

d dt x t = ∂ ξ h + (x t , ξ t ) = c(x t ) ξ t |ξ t | , d dt ξ t = -∂ x h + (x t , ξ t ) = -∂ x c(x t )|ξ t |.
The curves (t, x ±t ) of R n+1 are called the rays of P . An individual first order (Gaussian) beam for the wave equation associated to a ray (t, x t ) has the form ω ε (t, x) = a 0 (t, x)e iψ(t,x)/ε , with a complex phase function ψ real-valued on (t, x t ), an amplitude function a 0 null outside a neighborhood of (t, x t ), and such that sup

t∈[0,T ] P ω ε (t, .) L 2 (Ω) = O(ε m ),
for some m > 0.

The construction of such a beam is achieved by making the amplitudes of P ω ε vanish on the ray up to fixed suitable orders [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF][START_REF] Katchalov | Inverse Boundary Spectral Problems[END_REF][START_REF] Macià | On the lack of observability for wave equations: a Gaussian beam approach[END_REF]]

P ω ε = ε -2 p(x, ∂ t ψ, ∂ x ψ)a 0 +ε -1 i 2∂ t ψ∂ t a 0 -2c 2 ∂ x ψ∂ x a 0 + P ψa 0 +h.o.t. e iψ/ε , (2.1) where p(x, τ, ξ) = c 2 (x)|ξ| 2 -τ 2
is the principal symbol of P and h.o.t. denotes higher order terms. The first equation is then the eikonal equation

p (x, ∂ t ψ(t, x), ∂ x ψ(t, x)) = 0 (2.2)
on x = x t up to order 2 (see Remark 2.1 in [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF] for an explanation of the choice of this specific order), which means

∂ α x [p (x, ∂ t ψ(t, x), ∂ x ψ(t, x))] | x=x t = 0 for |α| ≤ 2.
Orders 0 and 1 of the previous equation are fulfilled on the ray by setting

∂ t ψ(t, x t ) = -h + (x t , ξ t ) and ∂ x ψ(t, x t ) = ξ t . (2.3)
Choosing ψ(0, x 0 ) as a real quantity, it follows that

ψ(t, x t ) is real. (2.4)
Order 2 of eikonal (2.2) on the ray may be written as a Riccati equation given an initial symmetric matrix ∂ 2 x ψ 0, x 0 with a positive definite imaginary part (see the proof of Lemma 2.56 p.101 in [START_REF] Katchalov | Inverse Boundary Spectral Problems[END_REF]).

d dt ∂ 2 x ψ(t, x t ) + H 21 (x t , ξ t )∂ 2 x ψ(t, x t ) + ∂ 2 x ψ(t, x t )H 12 (x t , ξ t ) + ∂ 2 x ψ(t, x t )H 22 (x t , ξ t )∂ 2 x ψ(t, x t ) + H 11 (x t , ξ t ) = 0, (2.5
The phase is defined beyond the ray as a polynomial of order 2 w.r.t. (xx t ) [START_REF] Tanushev | Mountain waves and Gaussian beams[END_REF]

ψ(t, x) = ψ(t, x t ) + ξ t • (x -x t ) + 1 2 (x -x t ) • ∂ 2 x ψ(t, x t )(x -x t ). (2.7)
Next, we make the term associated to the power ε -1 in the expansion (2.1) vanish on (t,

x t ) 2∂ t ψ∂ t a 0 -2c 2 ∂ x ψ∂ x a 0 + P ψa 0 = 0 on (t, x t ), (2.8)
which leads to a linear ordinary differential equation (ODE) on a 0 (t, x t ). The amplitude is then chosen under the form

a 0 (t, x) = χ d (x -x t )a 0 (t, x t ),
where d is a positive parameter. The constructed beams are thus defined for all (t, x) ∈ R n+1 and they satisfy the estimate

ε -n 4 +1 P ω ε (t, .) L 2 (Ω) = O( √ ε) uniformly w.r.t. t ∈ [0, T ].
Note that Gaussian beams for P associated to the ray (t, x -t ) are ω ε (-t, x).

Incident and reflected beams in a convex domain

Assume that c(x) is constant for dist(x, Ω) larger than some constant C > 0. Given a point (y, η) in the phase space

o T * R n , where o T * U denotes U × (R n \{0}) if U is an open set of R n , the Hamiltonian flow ϕ t 0 (y, η) = (x t 0 (y, η), ξ t 0 (y, η)) satisfying: d dt x t 0 = c(x t 0 ) ξ t 0 |ξ t 0 | , d dt ξ t 0 = -∂ x c(x t 0 )|ξ t 0 |, x t 0 | t=0 = y, ξ t 0 | t=0 = η, η = 0,
is called incident flow. A beam associated to the incident ray (t, x t 0 ) is denoted ω 0 ε and called an incident beam. Since we have dependence w.r.t. the initial conditions (y, η), we write the incident beam as ω 0 ε (t, x, y, η) = a 0 (t, x, y, η)e iψ0(t,x,y,η)/ε .

Let R be the reflection involution

R : o T * R n | ∂Ω → o T * R n | ∂Ω (X, Ξ) → X, (Id -2ν(X)ν(X) T )Ξ , (2.9)
where ν denotes the exterior normal field to ∂Ω. We shall only consider initial points (y, η) ∈ B = ∪ t∈R ϕ t 0 ( o T * Ω) giving rise to rays that enter the domain Ω at some instant. Each associated flow ϕ t 0 (y, η) hits the boundary twice. Reflection of ϕ t 0 (y, η) at the exit time t = T 1 (y, η) s.t.

x T1(y,η) 0 (y, η) ∈ ∂Ω and ẋT1(y,η)

0 (y, η) • ν x T1(y,η) 0 (y, η) > 0
gives birth to the reflected flow ϕ t 1 (y, η) = (x t 1 (y, η), ξ t 1 (y, η)) defined by the condition ϕ

T1(y,η) 1 (y, η) = R oϕ T1(y,η) 0 (y, η).
Similarly, we also define the reflection time T -1 (y, η) and the flow ϕ t -1 (y, η) by reflecting ϕ t 0 (y, η) as follows

x

T-1(y,η) 0 (y, η) ∈ ∂Ω and ẋT-1(y,η)

0 (y, η) • ν x T-1(y,η) 0 (y, η) < 0, ϕ T-1(y,η) -1 (y, η) = R o ϕ T-1(y,η) 0 (y, η).
We denote, for k = ±1, the reflected beams by ω k ε (t, x, y, η) = a k 0 (t, x, y, η)e iψ k (t,x,y,η)/ε .

These beams are associated to the reflected bicharacteristics ϕ t k . Let us introduce, for k = 0, ±1, the boundary amplitudes d k -mB +j s.t.

Bω k ε = mB j=0 ε -mB +j d k -mB +j e iψ k /ε .
Above, m B denotes the order of B (m B = 0 for Dirichlet and m B = 1 for Neumann).

The construction of the reflected phases and amplitudes is achieved by imposing that 1. the time and tangential derivatives of

ψ k equal at (T k , x T k 0 ) those of ψ 0 up to order 2, 2. d 0 -mB + d k -mB (T k , x T k 0 ) = 0, for k = ±1.
These constraints uniquely determine the reflected phases and amplitudes, once the incident ones are fixed [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF]. If T is sufficiently small, at most one reflection occurs in the interval [0, T ] and in the interval [-T, 0] for a fixed starting position and vector speed (y, η) ∈ o T * Ω, and the following boundary estimates are satisfied [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF]

B ε -n 4 +1 ω 0 ε (., y, η) + ε -n 4 +1 ω 1 ε (., y, η) H s ([0,T ]×∂Ω) = O(ε -mB -s+ 3 2 ),
and

B ε -n 4 +1 ω 0 ε (., y, η) + ε -n 4 +1 ω -1 ε (., y, η) H s ([-T,0]×∂Ω) = O(ε -mB -s+ 3 2 ),
for s ≥ 0.

Gaussian beam summation

The construction of asymptotic solutions to the IBVP (1.1a)-(1.1b) with initial conditions (1.1c') having a suitable frequency support (see below) is recalled, through the Gaussian beam summation introduced in [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF]. We focus on a superposition of first order beams, for which exact expressions of the phases and amplitudes are displayed in Subsection 2.2.2. These beams lead to a first order approximate solution, close to the exact one up to √ ε. Then, the derivatives of the first order solution will be approximated by some Gaussian type integrals.

Construction of the approximate solution

In [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF], we have constructed a family of asymptotic solutions to the IBVP for the wave equation for initial data satisfying (H1), (H4) and an additional hypothesis (H5) concerning their FBI transforms.

Let us recall here that the FBI transform (see [START_REF] Martinez | An Introduction to Semiclassical and Microlocal Analysis[END_REF]) is, for a given scale ε, the operator

T ε : L 2 (R n ) → L 2 (R 2n ) defined by T ε (a)(y, η) = c n ε -3n 4 R n a(x)e iη•(y-x)/ε-(y-x) 2 /(2ε) dx, c n = 2 -n 2 π -3n 4 , a ∈ L 2 (R n ), (2.10
) with adjoint operator given by

T * ε (f )(x) = c n ε -3n 4 R 2n f (y, η)e iη•(x-y)/ε-(x-y) 2 /(2ε) dydη, f ∈ L 2 (R 2n ).
As the Fourier transform, the FBI transform is an isometry, satisfying T * ε T ε = Id. The extra assumption on the initial data needed in [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF] is

T ε u I ε L 2 (R n ×R c η ) = O(ε ∞ ) and T ε v I ε L 2 (R n ×R c η ) = O(ε ∞ ), ( H5
)
where R c η denotes the complementary in R n of some ring

R η = {η ∈ R n , r 0 ≤ |η| ≤ r ∞ }, 0 < r 0 ≪ r ∞ .
In general, this assumption may be not satisfied. Therefore, we construct a family of initial data (u I ε,r0,r∞ , v I ε,r0,r∞ ) close to (u I ε , v I ε ), satisfying the same assumptions as (H1), (H4) and having FBI transforms small in L 2 (R n × R c η ). Letting r 0 go to 0 and r ∞ go to +∞ makes these data approach (u I ε , v I ε ) in a sense that will be specified in Section 3.3. In any case, the needed convergence is weaker than a L 2 convergence since we are interested in the study of Wigner measures.

Let us first truncate T ε u I ε and T ε v I ε outside R η by multiplying them by a cut-off

γ r0,r∞ ∈ C ∞ 0 (R n , [0, 1]) supported in the interior of R η γ r0,r∞ = χ r∞/2 (1 -χ 4r0 ).
(2.11) Lemma 4.5 from the Appendix (Section 4) yields

T ε T * ε γ r0,r∞ (η)T ε u I ε L 2 (R n ×R c η ) = O(ε ∞ ), and T ε T * ε γ r0,r∞ (η)T ε v I ε L 2 (R n ×R c η ) = O(ε ∞ ).
In order to have data supported in fixed compact sets of Ω independent of ε, we multiply

T * ε γ r0,r∞ (η)T ε u I ε , T * ε γ r0,r∞ (η)T ε v I ε by a cut-off ρ ∈ C ∞ 0 (R n , [0, 1]
) supported in Ω, and consider

u I ε,r0,r∞ = ρT * ε γ r0,r∞ (η)T ε u I ε and v I ε,r0,r∞ = ρT * ε γ r0,r∞ (η)T ε v I ε . (1.1c') It is assumed that ρ(x) = 1 if dist(x, suppu I ε ∪ suppv I ε ) ≤ C for some C > 0. The required estimates T ε u I ε,r0,r∞ L 2 (R n ×R c η ) = O(ε ∞ ) and T ε v I ε,r0,r∞ L 2 (R n ×R c η ) = O(ε ∞ ) (H5')
are fulfilled since Lemma 4.4 from the Appendix implies that

(1 -ρ)T * ε γ r0,r∞ (η)T ε u I ε L 2 x e -C/ε and (1 -ρ)T * ε γ r0,r∞ (η)T ε v I ε L 2 x e -C/ε .
(2.12) Using the boundedness of the operator

T * ε γ r0,r∞ T ε from L 2 (R n ) to L 2 (R n ) and the relations ∂ yj T ε = T ε ∂ xj , ∂ xj T * ε = T * ε ∂ yj , (2.13)
obtained by integrations by parts in the expressions of T ε and T * ε , one can show that the new initial data (u I ε,r0,r∞ , v I ε,r0,r∞ ) is also uniformly bounded w.r.t. ε in

H 1 (Ω) × L 2 (Ω). Let ρ ′ be a cut-off of C ∞ 0 (R n , [0, 1]) supported in a compact K y ⊂ Ω and satis- fying ρ ′ (y) = 1 if dist(y, suppρ) < C for some C > 0,
and

γ ′ a cut-off of C ∞ 0 (R n , [0, 1]) supported in K η ⊂ R n \{0} s.t. γ ′ ≡ 1 on R η .
Without loss of generality, we assume that either the incident ray or the reflected one propagating in the positive sense is in the interior of the domain at the instant T (x T 0 (y, η) ∈ Ω or x T 1 (y, η) ∈ Ω) when y varies in K y and η in R n \{0}. This is always possible upon reducing T because the number of reflections for initial position and vector speed varying in K y × (R n \{0}) is uniformly bounded (see Section 2.3 of [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF] for similar arguments). And similarly for the instant -T for rays propagating in the negative sense. Then, the IBVP (1.1a)-(1.1b) with initial conditions (1.1c') has a family of approximate solutions u appr ε,r0,r∞ in C 0 ([0, T ], H 1 (Ω))∩C 1 ([0, T ], L 2 (Ω)) obtained as a summation of first order beams. A general result using a superposition of beams of any order was proven in [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF], and it reads for first order beams as follows:

Proposition 1 ([4], Theorem 1.1). Denote for t ∈ [0, T ] and x ∈ R n the following superposition of Gaussian beams u appr ε,r0,r∞ (t, x) = 1 2 ε -3n 4 +1 c n R 2n ρ ′ (y)γ ′ (η)T ε v I ε,r0,r∞ (y, η) k=0,1 ω k ε ′ (t, x, y, η) - k=0,-1 ω k ε ′ (-t, x, y, η) dydη + 1 2 ε -3n 4 +1 c n R 2n ρ ′ (y)γ ′ (η)ε -1 T ε u I ε,r0,r∞ (y, η) k=0,1 ω k ε (t, x, y, η) + k=0,-1 ω k ε (-t, x, y, η) dydη.
Above, ω 0 ε , ω 0 ε ′ are incident Gaussian beams with the same phase ψ 0 satisfying at

t = 0 ψ 0 (0, x, y, η) = η • (x -y) + i 2 (x -y) 2 (2.14)
and different amplitudes a 0 0 , a 0 0

′ satisfying a 0 0 (0, x, y, η) = χ d (x -y), i∂ t ψ 0 a 0 0 ′ (0, x, y, η) = χ d (x -y) + O(|x -y|). (2.15) ω ±1 ε and ω ±1 ε ′ denote the associated reflected beams. Then u appr ε,r0,r∞ is asymptotic to u ε,r0,r∞ the exact solution of the problem (1.1a)-(1.1b) with initial conditions (1.1c') in the sense that sup t∈[0,T ] u ε,r0,r∞ -u appr ε,r0,r∞ H 1 (Ω) ≤ C(r 0 , r ∞ , Ω, T ) √ ε,
and

sup t∈[0,T ] ∂ t u ε,r0,r∞ -∂ t u appr ε,r0,r∞ L 2 (Ω) ≤ C(r 0 , r ∞ , Ω, T ) √ ε.
We refer to [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF] for further details, and just mention that the proof relies on the use of a family of approximate operators acting from

L 2 (R 2n ) to L 2 (R n ).
A simple version of the estimate of these operators norms is recalled in Section 4.2 of the Appendix.

Expression of the phases and amplitudes

In order to compute the first order beams, we begin by analyzing the relationship between the incident phase and amplitudes, and the Jacobian matrix of the incident flow. A similar relationship involving the reflected phases and amplitudes and the reflected flows will be also given.

The requirement (2.3) for the incident phase implies that

d dt ψ 0 (t, x t 0 ) = ∂ t ψ 0 (t, x t 0 ) + ∂ x ψ 0 (t, x t 0 ) • ẋ0 t = 0.
Taking into account the initial null value ψ 0 (0, y) = 0 chosen in (2.14), one gets a null phase on the ray ψ 0 (t, x t 0 ) = 0. With the aim of computing ∂ 2

x ψ 0 (t, x t 0 ), we note that the Jacobian matrix of the bicharacteristic F t 0 = Dϕ t 0 satisfies the linear ordinary differential system

d dt F t 0 = JH(x t 0 , ξ t 0 )F t 0 , F 0 0 = Id,
where J = 0 Id -Id 0 is the standard symplectic matrix. Writing F t 0 as

F t 0 = D y x t 0 D η x t 0 D y ξ t 0 D η ξ t 0 leads to the following ordinary differential system on (U t 0 , V t 0 ) = (D y x t 0 + iD η x t 0 , D y ξ t 0 + iD η ξ t 0 ) d dt U t 0 = H 21 (x t 0 , ξ t 0 )U t 0 + H 22 (x t 0 , ξ t 0 )V t 0 , (2.16
)

d dt V t 0 = -H 11 (x t 0 , ξ t 0 )U t 0 -H 12 (x t 0 , ξ t 0 )V t 0 .
(2.17)

Note that F t 0 is a symplectic matrix, i.e.

(F t 0 ) T JF t 0 = J.
Using the symmetry of the following matrices

(D y x t 0 ) T D y ξ t 0 , (D η x t 0 ) T D η ξ t 0 , D y x t 0 (D η x t 0 ) T , and D y ξ t 0 (D η ξ t 0 ) T
and the relations

(D y x t 0 ) T D η ξ t 0 -(D y ξ t 0 ) T D η x t 0 = Id and D y x t 0 (D η ξ t 0 ) T -D η x t 0 (D y ξ t 0 ) T = Id, one has (U t 0 ) T V t 0 = (V t 0 ) T U t 0 , (V t 0 ) T Ū t 0 -(U t 0 ) T V t 0 = 2iId and U t 0 is invertible. (2.18)
Putting together (2.16), (2.17) and (2.18) shows that V t 0 (U t 0 ) -1 is a symmetric matrix with a positive definite imaginary part and fulfills the Riccati equation (2.5) with initial value iId. Since this is the initial condition for ∂ 2

x ψ 0 (t,

x t 0 ) given in (2.14), it follows that ∂ 2 x ψ 0 (t, x t 0 ) = V t 0 (U t 0 ) -1 . (2.19)
The incident beams amplitudes are computed as follows. Using (2.3) and the Hamiltonian system satisfied by (x t 0 , ξ t 0 ), the equation (2.8) at order zero yields the following transport equation for the value of the amplitude on the ray [START_REF] Katchalov | Inverse Boundary Spectral Problems[END_REF]

d dt a 0 0 ( ′ ) (t, x t 0 ) + 1 2 Tr H 21 (x t 0 , ξ t 0 ) + H 22 (x t 0 , ξ t 0 )∂ 2 x ψ 0 (t, x t 0 ) a 0 0 ( ′ ) (t, x t 0 ) = 0, ( 2
.20) which may be written using the matrices U t 0 and V t 0 as

d dt a 0 0 ( ′ ) (t, x t 0 ) + 1 2 Tr H 21 (x t 0 , ξ t 0 )U t 0 + H 22 (x t 0 , ξ t 0 )V t 0 (U t 0 ) -1 a 0 0 ( ′ ) (t, x t 0 ) = 0.
The time evolution for U t 0 , see (2.16), combined with the choice of the initial values a 0 0 (0, y) = 1 and a 0 0

′ (0, y) = (-ic(y)|η|) -1 from (2.15), yields a 0 0 (t, x t 0 ) = det U t 0 -1 2 and a 0 0 ′ (t, x t 0 ) = i(c(y)|η|) -1 det U t 0 -1 2 .
Above the square root is defined by continuity in t from 1 at t = 0. The expression of the reflected phases

ψ k , k = ±1, is similar to the incident phase. In fact, since d dt (ψ k (t, x t k )) = 0 and ψ k (T k , x T k 0 ) = ψ 0 (T k , x T k 0 ) because of the requirement 1 p.8, we get ψ k (t, x t k ) = 0.
We then apply the general relation between incident and reflected beams phases given in Lemma 4.1 from the Appendix, to compute the Hessian matrices of ψ ±1 on the rays. As far as we know, the result stated in this Lemma is particularly simple enough so that we stated it in the Appendix 4.1. The matrices ∂ 2

x ψ ±1 (t, x ±t ±1 ) can also be computed by solving the Riccati equations with the proper values at the instants of reflections t = T ±1 (see eg. [START_REF] Norris | Elastic Gaussian wave packets in isotropic media[END_REF][START_REF] Tanushev | Gaussian beam decomposition of high frequency wave fields[END_REF]). One gets (see Appendix 4.1)

∂ 2 x ψ k (t, x t k ) = V t k (U t k ) -1 where U t k = D y x t k + iD η x t k and V t k = D y ξ t k + iD η ξ t k . As ϕ t k is symplectic, (U t k , V t k ) share the same properties (2.18) as (U t 0 , V t 0 ) (U t k ) T V t k = (V t k ) T U t k , (V t k ) T Ū t k -(U t k ) T V t k = 2iId and U t k is invertible. (2.21)
The reflected amplitudes evaluated on the rays have an expression similar to the incident amplitudes (see Appendix 4.1)

a k 0 (t, x t k ) = -si det U t k -1 2 and a k 0 ′ (t, x t k ) = s(c(y)|η|) -1 det U t k -1 2 for k = ±1,
where the square root is defined by continuity from i det

U T k 0 -1 2 at t = T k , s = -1
for the Dirichlet boundary condition and s = 1 for the Neumann condition. We summarize the previous form of the beams in the following result:

Lemma 2.1 For k = 0, ±1, the incident and reflected beams ω k ε have the form

ω k ε ( ′ ) (t, x) = β k χ d (x -x t k )a ( ′ ) k (t)e iψ k /ε , with β 0 = 1, β 1 = β -1 = -si, a k (t) = [det U t k ] -1 2 , a ′ k (t) = i(c(y)|η|) -1 [det U t k ] -1 2 , ψ k = ξ t k • (x -x t k ) + i 2 (x -x t k ) • Λ k (t)(x -x t k ), and Λ k (t) = -iV t k (U t k ) -1 .

Gaussian integrals

It follows that the approximate solution u appr ε,r0,r∞ has the form (recall the dependence of Gaussian beams w.r.t. variables (y, η))

u appr ε,r0,r∞ (t, x) = 1 2 ε -3n 4 +1 c n R 2n ρ ′ (y)γ ′ (η) k=0,1 χ d (x -x t k )β k p ε,k (t, y, η)e iψ k (t,x,y,η)/ε dydη + 1 2 ε -3n 4 +1 c n R 2n ρ ′ (y)γ ′ (η) k=0,-1 χ d (x -x -t k )β k q ε,k (-t, y, η) e iψ k (-t,x,y,η)/ε dydη, with p ε,k (t, y, η) = a k (t, y, η)ε -1 T ε u I ε,r0,r∞ (y, η) + a ′ k (t, y, η)T ε v I ε,r0,r∞ (y, η), and q ε,k (t, y, η) = a k (t, y, η)ε -1 T ε u I ε,r0,r∞ (y, η) -a ′ k (t, y, η)T ε v I ε,r0,r∞ (y, η).
Because of the phases expression given in (2.7), time and spatial derivatives of u appr ε,r0,r∞ may be written as a sum of integrals of the form

ε -3n 4 R 2n ρ ′ (y)γ ′ (η)f ε (y, η)ε j (x -x t k ) α r k j,α (t, x, y, η) e iψ k (t,x,y,η)/ε dydη, j, k = 0, 1, |α| ≤ 2,
arising from differentiation of ω 0 ε (t, .) and ω 1 ε (t, .). Other terms of the same form originate from derivatives of ω 0 ε (-t, .) and ω -1 ε (-t, .). f ε stands for ε -1 T ε u I ε,r0,r∞ or T ε v I ε,r0,r∞ and r k j,α are smooth functions vanishing for |x -

x t k | ≥ d. For a function f depending on (t, x, z, θ) ∈ R n+1 × B and k = 0, ±1, let f k (t, x, z, θ) = f (t, x, ϕ t k -1 (z, θ)). Set K k z,θ (t) = ϕ t k (K y × K η ). Let Π k (t) be a cut-off of C ∞ 0 (R 2n , [0, 1]) supported in B and satisfying Π k (t) ≡ 1 on K k z,θ (t). The volume preserving change of variables (z, θ) = ϕ t k (y, η)
transforms the previous integrals as

ε -3n 4 R 2n Π k (t) ρ ′ ⊗ γ ′ k f ε k ε j (x -z) α r k j,α k e i ψ k k (t,x,z,θ)/ε dzdθ, j, k = 0, 1, |α| ≤ 2.
(2.22) We can write the leading terms obtained for j = 0 and α = 0 using Gaussian type integrals I ε (h, Φ) defined as

I ε (h, Φ)(t, x) = ε -3n 4 c n R 2n h(t, z, θ)e iΦ(t,x,z,θ)/ε dzdθ,
for a given phase function Φ ∈ C ∞ (R n+1 t,x × B, C) polynomial of order 2 in xz and satisfying, for t ∈ [0, T ] and (z, θ) ∈ B

Φ(t, z, z, θ) is real, ∂ x Φ(t, z, z, θ) = θ, Im ∂ 2 x Φ(t, z, z, θ) is positive definite, (2.23)
and a given amplitude function h ∈ C 0 ([0, T ], L 2 (R 2n z,θ )) supported for every fixed t ∈ [0, T ] in a compact of B. By Proposition 3 in the Appendix, one has

R 2n h(t, z, θ)χ(x -z)e iΦ(t,x,z,θ)/ε dzdθ L 2 x h(t, .) L 2 z,θ .
Noticing that e iΦ/ε is exponentially decreasing for |x -z| ≥ 1, one can use the following crude estimate

|x-z|≥a h(t, z, θ)e iΦ(t,x,z,θ)/ε dzdθ L 2 x e -C/ε h(t, .) L 2 z,θ for a > 0 (2.24)
to deduce that I ε (h, Φ)(t, .) is uniformly bounded w.r.t. ε in L 2

x . The same notation I ε (h, Φ) will be also used for a vector valued function h.

The contribution of the terms (2.22) with j = 1 or |α| ≥ 1 to the derivatives of u appr ε,r0,r∞ is of order √ ε as stated in the following Lemma, whose proof is given Appendix 4.2 and relies on the approximation operators defined therein.

Lemma 2.2 ∂ t u appr ε,r0,r∞ (t, .) is uniformly bounded w.r.t. ε in L 2 (R n ) and satisfies ∂ t u appr ε,r0,r∞ (t, x) = 1 2 v + t,ε (t, x) -v - t,ε (-t, x) + O( √ ε) in L 2 (R n ) uniformly w.r.t. t ∈ [0, T ],
where

(v + t,ε ) and (v - t,ε ) are sequences of L 2 (R n ) uniformly bounded w.r.t. ε given by v + t,ε = k=0,1 β k I ε (-ic(z)|θ|Π k ρ ′ ⊗ γ ′ k p ε,k k , ψ k k ), v - t,ε = k=0,-1 β k I ε (-ic(z)|θ|Π k ρ ′ ⊗ γ ′ k q ε,k k , ψ k k ).
Likewise, ∂ x u appr ε,r0,r∞ (t, .) is uniformly bounded w.r.t. ε in L 2 (R n ) n and satisfies

∂ x u appr ε,r0,r∞ (t, x) = 1 2 v + x,ε (t, x) + v - x,ε (-t, x) + O( √ ε) in L 2 (R n ) n uniformly w.r.t. t ∈ [0, T ],
where

(v + x,ε ) and (v - x,ε ) are sequences of L 2 (R n ) n uniformly bounded w.r.t. ε given by v + x,ε = k=0,1 β k I ε (iθΠ k ρ ′ ⊗ γ ′ k p ε,k k , ψ k k ), v - x,ε = k=0,-1 β k I ε (iθΠ k ρ ′ ⊗ γ ′ k q ε,k k , ψ k k ).

Wigner transforms and measures

We now compute the scalar measures associated to the sequences ∂ t u appr ε,r0,r∞ (t, .) and c∂ x u appr ε,r0,r∞ (t, .) . As |β k | = 1, the Wigner transform associated to v + t,ε (t, .) is a finite sum of terms of the form

w ε I ε (f k t,ε , Φ k )(t, .), I ε (f l t,ε , Φ l )(t, .) ,
where k, l = 0, 1,

f k t,ε = c|θ|Π k ρ ′ ⊗ γ ′ k p ε,k k and Φ k = ψ k k .
As regards the Wigner transforms associated to cv +

x,ε (t, .) , since c is uniformly continuous on R n , one has by a classical result ( [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF], p.8)

w ε cv + x,ε (t, .), cv + x,ε (t, .) ≈ c 2 w ε v + x,ε (t, .), v + x,ε (t, .) in R 2n , (3.1)
and therefore the involved quantities have the form

c 2 w ε I ε (f k x,ε , Φ k )(t, .), I ε (f l x,ε , Φ l )(t, .) , with f k x,ε = θΠ k ρ ′ ⊗ γ ′ k p ε,k k .
Similarly, we define for k = 0, -1 the sequences

g k t,ε = c|θ|Π k ρ ′ ⊗ γ ′ k q ε,k k , which
are needed when considering the Wigner transform associated with cv - t,ε (-t, .) and the cross Wigner transform between cv + t,ε (t, .) and cv - t,ε (-t, .) , as well as

g k x,ε = θΠ k ρ ′ ⊗ γ ′ k q ε,k k .
Then, forgetting the powers of ε factors, all the previous Wigner transforms tested on cut-off functions have the form

R 6n T ε κ ε k (z, θ) T ε τ ε l (z ′ , θ ′ )b k,l 1 (z, θ, z ′ , θ ′ , x, v)e iΨ k,l 1 (z,θ,z ′ ,θ ′ ,x,v)/ε dzdθdz ′ dθ ′ dxdv, (3.2)
with κ ε , τ ε = ε -1 u I ε,r0,r∞ , v I ε,r0,r∞ and k, l = 0, ±1, or after expanding the FBI transforms

R 8n κ ε (w)τ ε (w ′ )b k,l 2 (z, θ, z ′ , θ ′ , x, v)e iΨ k,l 2 (w,w ′ ,z,θ,z ′ ,θ ′ ,x,v)/ε dwdw ′ dzdθdz ′ dθ ′ dxdv. (3.
3) This type of oscillating integrals is traditionally estimated by the stationary phase theorem. For example, this method was successfully used in [START_REF] Castella | The radiation condition at infinity for the high frequency Helmholtz equation with source term: a wave packet approach[END_REF] for the computation of a Wigner measure for smooth data. There the phase was complex and its Hessian matrix restricted to the stationary set was assumed to be non-degenerate in the normal direction to this set. However, in our case, the amplitude is not smooth as no such assumption was made on u I ε and v I ε , and we cannot estimate immediately the global integral (3.3) by the same techniques. One possibility of solving this issue would be to resort to the stationary phase theorem with a complex phase depending on parameters for estimating

R 6n b 2 (z, θ, z ′ , θ ′ , x, v)e iΨ2(w,w ′ ,z,θ,z ′ ,θ ′ ,x,v)/ε dzdθdz ′ dθ ′ dxdv,
and then study the whole integral involving κ ε (w)τ ε (w ′ ).

An alternative method was used in [START_REF] Robinson | Semiclassical mechanics for time-dependent Wigner functions[END_REF], where an integral of the form (3.2) associated to the Wigner transform for the Schrödinger equation with a WKB initial condition was simplified by elementary computations into an integral over R 4n .

Though the method therein faced difficulties in deducing the exact relation between the Wigner measure of the solution and of the initial data, we adapt the result of [START_REF] Robinson | Semiclassical mechanics for time-dependent Wigner functions[END_REF] to our problem in Section 3.1 and complete the analysis to prove the propagation along the flow of the microlocal energy density of u appr ε,r0,r∞ as ε → 0 in Section 3.2. The proof is simple and elementary and the computations are made in an explicit way. Section 3.3 is devoted to the Wigner measures associated to the derivatives of u ε the exact solution of (1.1).

Wigner transform for Gaussian integrals

The sequences (f k t,ε ), (f k x,ε ), (g l t,ε ) and (g l x,ε ) are uniformly bounded w.r.t. ε in L 2 (R 2n ) and their supports are contained in a fixed compact independent of ε. Slight modifications of the computations of [START_REF] Robinson | Semiclassical mechanics for time-dependent Wigner functions[END_REF] lead to the following more general result: Lemma 3.1 Let (f ε ) and (g ε ) be sequences uniformly bounded in L 2 (R 2n ) and having their supports contained in a fixed compact independent of ε. Let F be an open set containing suppf ε ∪ suppg ε and Φ, Ψ be phase functions in

C ∞ (R n x × F, C) satisfying Φ(x, z, θ) = r Φ (z, θ) + θ • (x -z) + i 2 (x -z) • H Φ (z, θ)(x -z), Ψ(x, z ′ , θ ′ ) = r Ψ (z ′ , θ ′ ) + θ ′ • (x -z ′ ) + i 2 (x -z ′ ) • H Ψ (z ′ , θ ′ )(x -z ′ ), for x ∈ R n and (z, θ), (z ′ , θ ′ ) ∈ F , with r Φ , r Ψ ∈ C ∞ (F, R
) and the matrices

H Φ , H Ψ ∈ C ∞ (F, M n (C)) having positive definite real parts. Then for φ ∈ C ∞ 0 (F, R) < w ε (I ε (f ε , Φ), I ε (g ε , Ψ)) , φ > = R 4n φ(s, σ)f ε (s + √ εr, σ + √ εδ)g * ε (s - √ εr, σ - √ εδ) A(Φ, Ψ)(s, σ)e iΘε(Φ,Ψ)(s,σ,r,δ) drdδdsdσ + o(1), where A(Φ, Ψ)(s, σ) = c 2 n 2 5n 2 π n 2 det H Φ (s, σ) + H Ψ (s, σ) -1 2 ,
and

Θ ε (Φ, Ψ)(s, σ, r, δ) =r Φ (s + √ εr, σ + √ εδ)/ε -r Ψ (s - √ εr, σ - √ εδ)/ε -2σ • r/ √ ε + i(r, δ) • Q H Φ (s, σ), H Ψ (s, σ) (r, δ).
The matrix Q H Φ (s, σ), H Ψ (s, σ) and the square root are defined in Lemma 3.2.

Proof: It consists in two steps. Firstly, the Fourier transform of a Gaussian type function is computed explicitly. Then, a Gaussian approximation is used for several smooth functions appearing in the Wigner transform integral.

For simplicity we denote u(x, z, θ) by u and u(x, z ′ , θ ′ ) by u ′ when integrating w.r.t. z, θ, z ′ , θ ′ . We also omit the index ε in the notation of f ε and g ε .

Step 1. Fourier transform. We note that the Wigner transform at point (x, ξ) ∈ R 2n may be written as

w ε (I ε (f, Φ), I ε (g, Ψ)) (x, ξ) =π -n c 2 n ε -5n 2 R 5n f g * ′ e irΦ/ε-ir ′ Ψ /ε+ix•(θ-θ ′ )/ε+i(θ ′ •z ′ -θ•z)/ε F v e -(v+x-z)•HΦ(v+x-z)/(2ε) × e -(v-x+z ′ )•HΨ ′ (v-x+z ′ )/(2ε) ((2ξ -θ -θ ′ )/ε) dvdzdz ′ dθdθ ′ .
The Fourier transform of a Gaussian functions product is given by the following Lemma, whose proof is postponed to the end of this Section:

Lemma 3.2 Let a, b ∈ R d and M, N ∈ M d (C
) symmetric matrices with positive definite real parts, then

F x e -(x-a)•M(x-a)/2 e -(x-b)•N (x-b)/2 (ξ) =(2π) d 2 (det(M + N )) -1 2 e -iξ•(b+a)/2-(b-a,ξ)•Q(M,N )(b-a,ξ)/4 ,
where Q(M, N ) is the symmetric symplectic matrix given by

Q(M, N ) = 2M (M + N ) -1 N i(N -M )(M + N ) -1 i(M + N ) -1 (N -M ) 2(M + N ) -1 ,
and the square root is defined as explained in Section 3.4 of [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I. Distribution Theory and Fourier Analysis[END_REF].

Moreover, Q(M, N )A(M, N ) = B(M, N ) with A(M, N ) = Id Id -iN iM and B(M, N ) = N M -iId iId
, and Q(M, N ) has a positive definite real part

Re Q(M, N ) = 2A(M, N ) * -1 Re N 0 0 Re M A(M, N ) -1 .
Hence

w ε (I ε (f, Φ), I ε (g, Ψ)) (x, ξ) =c 2 n 2 n 2 π -n 2 ε -2n R 4n f g * ′ det(H Φ + H Ψ ′ ) -1 2 e irΦ/ε-ir ′ Ψ /ε e i(θ+θ ′ -2ξ)•(z-z ′ )/(2ε)+i(θ-θ ′ )•x/ε+i(θ ′ •z ′ -θ•z)/ε e -(2x-z-z ′ ,2ξ-θ-θ ′ )•Q(HΦ,HΨ ′ )(2x-z-z ′ ,2ξ-θ-θ ′ )/(4ε) dzdz ′ dθdθ ′ .
Making the changes of variables

(z, z ′ ) = (s + √ ǫr, s - √ ǫr), (θ, θ ′ ) = (σ + √ ǫδ, σ - √ ǫδ),
and writing f + for f (s + √ εr, σ + √ εδ) and g -for g(s -√ εr, σ -√ εδ) leads to

w ε (I ε (f, Φ), I ε (g, Ψ)) (x, ξ) =c 2 n 2 5n 2 π -n 2 ε -n R 4n f + g * -det(H Φ+ + HΨ -) -1 2 e irΦ + /ε-irΨ -/ε+2iδ•(x-s)/ √ ε e -2iξ•r/ √ ε-(x-s,ξ-σ)•Q(HΦ + ,HΨ -)(x-s,ξ-σ)/ε drdδdsdσ.
Step 2. Gaussian approximations. Taking the duality product of the Wigner transform with a test function φ ∈ C ∞ 0 (F, R), and after setting (x ′ , ξ ′ ) = (xs, ξσ)/ √ ε, one has

< w ε (I ε (f, Φ), I ε (g, Ψ)) , φ > =c 2 n 2 5n 2 π -n 2 R 6n φ(s + √ εx ′ , σ + √ εξ ′ )f + g * -det(H Φ+ + H Ψ-) -1 2 e irΦ + /ε-irΨ -/ε e -2iσ•r/ √ ε+2i(x ′ ,ξ ′ )•(δ,-r)-(x ′ ,ξ ′ )•Q(HΦ + ,HΨ -)(x ′ ,ξ ′ ) dx ′ dξ ′ drdδdsdσ. (3.4)
Let ρ ′ f and ρ ′ g be cut-off functions supported in F s.t. ρ ′ f ≡ 1 on a fixed compact containing suppf and ρ ′ g ≡ 1 on a fixed compact containing suppg, and consider

b ε : (x ′ , ξ ′ , s, σ, r, δ) → φ(s + √ εx ′ , σ + √ εξ ′ ) -φ(s, σ) ρ ′ f + ρ ′ g -e -(x ′ ,ξ ′ )•Q(HΦ + ,HΨ -)(x ′ ,ξ ′ ) .
The r.h.s. of (3.4) may be written as

< w ε (I ε (f, Φ), I ε (g, Ψ)) , φ > =c 2 n 2 5n 2 π -n 2 R 6n φ(s, σ)f + g * -det(H Φ+ + H Ψ-) -1 2 e irΦ + /ε-irΨ -/ε-2iσ•r/ √ ε e 2i(x ′ ,ξ ′ )•(δ,-r)-(x ′ ,ξ ′ )•Q(HΦ + ,HΨ -)(x ′ ,ξ ′ ) dx ′ dξ ′ drdδdsdσ + c 2 n 2 5n 2 π -n 2 R 4n det(H Φ+ + H Ψ-) -1 2 f + g * -e irΦ + /ε-irΨ -/ε-2iσ•r/ √ ε F (x ′ ,ξ ′ ) b ε (-2δ, 2r, s, σ, r, δ)drdδdsdσ. (3.5)
Leibnitz formula yields for a multiindex α

∂ α x ′ ,ξ ′ b ε (x ′ , ξ ′ , s, σ, r, δ) =ρ ′ f + ρ ′ g -φ(s + √ εx ′ , σ + √ εξ ′ ) -φ(s, σ) ∂ α x ′ ,ξ ′ e -(x ′ ,ξ ′ )•Q(HΦ + ,HΨ -)(x ′ ,ξ ′ ) + ρ ′ f + ρ ′ g - β+γ=α,β =0 C(β, γ)ε |β| 2 ∂ β x ′ ,ξ ′ φ(s + √ εx ′ , σ + √ εξ ′ ) × ∂ γ x ′ ,ξ ′ e -(x ′ ,ξ ′ )•Q(HΦ + ,HΨ -)(x ′ ,ξ ′ ) . As (s + √ εr, σ + √ εδ) varies in suppρ ′ f and (s - √ εr, σ - √ εδ) varies in suppρ ′ g , one can find by continuity a constant C > 0 s.t. Re Q(H Φ+ , H Ψ-) ≥ CId on supp(ρ ′ f + ρ ′ g -).
Since (s, σ) and √ ε(r, δ) are bounded on supp(ρ

′ f + ρ ′ g -), (3.6)
it follows that there exists a constant C ′ > 0 s.t.

|∂ α x ′ ,ξ ′ b ε (x ′ , ξ ′ , s, σ, r, δ)| √ εe -C ′ (x ′ ,ξ ′ ) 2 for all (x ′ , ξ ′ , s, σ, r, δ),
which leads to

|F (x ′ ,ξ ′ ) b ε (-2δ, 2r, s, σ, r, δ)| √ ε(1 + (r, δ) 2 ) -n-1 for all (s, σ, r, δ).
The second integral in the r.h.s. of (3.5) is then dominated by

√ ε R 4n |f + ||g -|(1 + (r, δ) 2 ) -n-1 drdδdsdσ.
We deduce by Cauchy-Schwartz inequality w.r.t. s, σ that

< w ε (I ε (f, Φ), I ε (g, Ψ)) , φ > -c 2 n 2 5n 2 π n 2 R 4n φ(s, σ) det(H Φ+ + H Ψ-) -1 2 f + g * -e irΦ + /ε-irΨ -/ε-2iσ.r/ √ ε e -(δ,-r)•Q(HΦ + ,HΨ -) -1 (δ,-r) drdδdsdσ √ ε f L 2 g L 2 ,
where we used det

Q(H Φ+ , H Ψ-) = 1 since Q(H Φ+ , H Ψ-) is symplectic. Next, we extend H Φ and H Ψ outside F as λH Φ + (1 -λ)Id and λH Ψ + (1 -λ)Id by using a cut-off λ ∈ C ∞ 0 (R 2n , [0, 1]) supported in F s.t. λ ≡ 1 on the compact set suppρ ′ f ∪ suppρ ′ g ∪ suppφ,
the extended matrices having positive definite real parts. The smoothness of these matrices implies by the mean value theorem and (3.6) that

det(H Φ+ + H Ψ-) -1 2 -det H Φ (s, σ) + H Ψ (s, σ) -1 2 √ ε|(r, δ)| on supp(φf + g * -)
.

By symplecticity and symmetry of Q(H Φ+ , H Ψ-), its inverse is -JQ(H Φ+ , H Ψ-)J.

Thus the quantity e -(δ,-r)•Q(HΦ + ,HΨ -) -1 (δ,-r)e -(r,δ)•Q(HΦ(s,σ),HΨ(s,σ))(r,δ)

is dominated by

(r, δ) • Q(H Φ+ , H Ψ-) -Q H Φ (s, σ), H Ψ (s, σ) (r, δ) × sup u∈[0,1]
e -u(r,δ)•Q(HΦ + ,HΨ -)(r,δ)-(1-u)(r,δ)•Q(HΦ(s,σ),HΨ(s,σ))(r,δ) .

The positivity of Re Q(H Φ+ , H Ψ-) and Re Q H Φ (s, σ), H Ψ (s, σ) and the mean value theorem for the matrix function

Q λH Φ + (1 -λ)Id, λH Ψ + (1 -λ)Id give by (3.6) e -(δ,-r)•Q(HΦ + ,HΨ -) -1 (δ,-r) -e -(r,δ)•Q(HΦ(s,σ),HΨ(s,σ))(r,δ) √ ε|(r, δ)| 3 e -C(r,δ) 2 for (s, σ) ∈ suppφ, (s + √ εr, σ + √ εδ) ∈ suppρ ′ f and (s - √ εr, σ - √ εδ) ∈ suppρ ′ g . It follows that < w ε (I ε (f, Φ), I ε (g, Ψ)) , φ > -c 2 n 2 5n 2 π n 2 R 4n φ(s, σ)(det[H Φ + H Ψ ]) -1 2 (s, σ)f + g * -e irΦ + /ε-irΨ -/ε-2iσ•r/ √ ǫ e -(r,δ)•Q(HΦ,HΨ)(s,σ)(r,δ) drdδdsdσ √ ε f L 2 g L 2 .
2 Proof: [Proof of Lemma 3.2] The matrix M + N has a positive definite real part and is thus non-singular. By elementary calculus we have

(x -a) • M (x -a) + (x -b) • N (x -b) =(b -a) • M (M + N ) -1 N (b -a) + x -(M + N ) -1 (M a + N b) • (M + N ) x -(M + N ) -1 (M a + N b) .
Using the value of the Fourier transform of a Gaussian function (see Theorem 7.6.1 of [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I. Distribution Theory and Fourier Analysis[END_REF]), it follows that

F x e -(x-a)•M(x-a)/2 e -(x-b)•N (x-b)/2 (ξ) =(2π) d 2 (det[M + N ]) -1 2 e -(b-a)•M(M+N ) -1 N (b-a)/2 e -iξ•(M+N ) -1 (Ma+N b)-ξ•(M+N ) -1 ξ/2 . Writing M = 1/2(M + N ) + 1/2(M -N ) and N = 1/2(M + N ) -1/2(M -N ),
we get the expression with the matrix Q(M, N ) and the relation

Q(M, N )A(M, N ) = B(M, N ).
One can easily show that

B(M, N ) T JB(M, N ) = 0 i(M + N ) -i(M + N ) 0 = A(M, N ) T JA(M, N ),
from which follows the symplecticity of Q(M, N ). Then write

Q(M, N ) + Q(M, N ) =A(M, N ) * -1 (A(M, N ) * B(M, N ) + B(M, N ) * A(M, N )) A(M, N ) -1
to obtain the value of Re Q(M, N ).

2 From now on, we drop the index ε in the notation of v ± t,ε , v ± x,ε , f k t,ε etc. for simplicity. We fix t ∈ [0, T ] and apply Lemma 3.1 with F = B on the sequences (f k t ), (f l t ) (respectively (f k x ), (f l x )) and the phase functions Φ k , Φ l for the Wigner transforms associated to v + t (t, .) (respectively (v + x (t, .))). To evaluate the cross Wigner transforms between v + t (t, .) and v - t (-t, .) (respectively (v + x (t, .)) and (v -

x (-t, .))) , we use this Lemma on the sequences (f k t ), (g l t ) (respectively (f k x ), (g l x )).

Wigner measures for superposed Gaussian beams

We shall prove that the cross Wigner transforms

w ε v + t (t, .), v - t (-t, .) , w ε v + x (t, .), v - x (-t, .) and w ε I ε (f k t,x , Φ k ), I ε (f l t,x , Φ l ) , w ε I ε (g k t,x , Φ k ), I ε (g l t,
x , Φ l ) with k = l do not contribute to the microlocal energy density limit E u appr ε,r0,r∞ (t, .) in o T * Ω. We compute Θ ε (Φ k , Φ k ) and A(Φ k , Φ k ) and analyze the transported FBI transforms at points (s± √ εr, σ ± √ εδ), which will complete the study of the Wigner measures for superposed Gaussian beams.

Firstly, we note that (1

-ρ ′ ⊗ γ ′ )Π k p ε,k L 2 y,η = O(ε ∞ ) for k = 0, ±1.
Indeed, γ ′ ≡ 1 on R η so one gets from (H5') that T ε u I ε,r0,r∞ , T ε v I ε,r0,r∞ have infinitely small contributions in L 2 (R n × supp(1γ ′ )).

On the other hand, dist supp(1ρ ′ ), suppu I ε,r0,r∞ ∪ suppv I ε,r0,r∞ > C. Then, Lemma 4.3 implies that T ε u I ε,r0,r∞ , T ε v I ε,r0,r∞ have infinitely small contributions in L 2 (supp(1ρ ′ ) × R n ). Therefore

w ε I ε (f k t , Φ k ), I ε (f l t , Φ l ) ≈A(Φ k , Φ l ) R 2n c|σ|Π k p ε,k k + c|σ|Π l pε,l l - e iΘε(Φ k ,Φ l ) drdδ in o T * Ω,
and a similar relation holds true for w ε I ε (f k x , Φ k ), I ε (f l x , Φ l ) . We start by approaching (c(s)|σ|) + (c(s)|σ|) -by c(s) 2 |σ| 2 in the previous integral

w ε I ε (f k t , Φ k ), I ε (f l t , Φ l ) ≈A(Φ k , Φ l )c(s) 2 |σ| 2 R 2n Π k p ε,k k + Π l pε,l l - e iΘε(Φ k ,Φ l ) drdδ in o T * Ω, (3.7)
and σ + σ * -by σσ * in the integral giving

w ε I ε (f k x , Φ k ), I ε (f l x , Φ l ) w ε I ε (f k x , Φ k ), I ε (f l x , Φ l ) ≈A(Φ k , Φ l )σσ * R 2n Π k p ε,k k + Π l pε,l l - e iΘε(Φ k ,Φ l ) drdδ in o T * Ω. (3.8)
Indeed, these approximations are proved with the help of the following Lemma

w ε (I ε (αf ε , Φ), I ε (βg ε , Ψ)) ≈ α βw ε (I ε (f ε , Φ), I ε (g ε , Ψ)) in F.
Proof: The proof relies on the use of Taylor's formula on ρ ′ f α and ρ ′ g β, where ρ ′ f and ρ ′ g are the cut-offs used in the proof of Lemma 3.1 (supported in F and equal to 1 on suppf ε and suppg ε respectively).

2 It follows by using (3.7) and (3.8) that

c 2 Tr w ε I ε (f k x , Φ k ), I ε (f l x , Φ l ) ≈ w ε I ε (f k t , Φ k ), I ε (f l t , Φ l ) in o T * Ω,
which leads to

w ε v + t (t, .), v + t (t, .) ≈ c 2 Tr w ε v + x (t, .), v + x (t, .) in o T * Ω.
Similarly

w ε v - t (-t, .), v - t (-t, .) ≈ c 2 Tr w ε v - x (-t, .), v - x (-t, .) in o T * Ω,
and

w ε v + t (t, .), v - t (-t, .) ≈ c 2 Tr w ε v + x (t, .), v - x (-t, .) in o T * Ω. (3.9)
The approximations linking the derivatives of u appr ε,r0,r∞ to v ± t,x given in Lemma 2.2 and equation (3.1) lead to 4E u appr ε,r0,r∞ (t, .)

≈w ε [v + t (t, .)] + c 2 Tr w ε [v + x (t, .)] + w ε [v - t (-t, .)] + c 2 Tr w ε [v - x (-t, .)] -w ε v + t (t, .), v - t (-t, .) + c 2 Tr w ε v + x (t, .), v - x (-t, .) -w ε v - t (-t, .), v + t (t, .) + c 2 Tr w ε v - x (-t, .), v + x (t, .) in R 2n
by using the standard estimate (see Proposition 1.1 in [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF])

| < w ε (a ε , b ε ), φ > | a ε L 2 (R n ) b ε L 2 (R n ) , (3.10)
for sequences (a ε ), (b ε ) in L 2 (R n ) and φ ∈ C ∞ 0 (R 2n , R). The cross terms between v + t,x and v - t,x cancel in o T * Ω by using (3.9), leading to

E u appr ε,r0,r∞ (t, .) ≈ 1 2 w ε [v + t (t, .)] + 1 2 w ε [v - t (-t, .)] in o T * Ω. (3.11)
Thus, we are left with the computation of the Wigner measure associated to (v + t ), computations being similar for (v - t ). One has

w ε [v + t ] ≈ k,l=0,1 c(s) 2 |σ| 2 A(Φ k , Φ l ) R 2n Π k p ε,k k + Π l pε,l l - e iΘε(Φ k ,Φ l ) drdδ in o T * Ω.
(3.12)

Moreover the inverse of the reflected/incident flow in

o T * Ω is a reflected/incident flow {ϕ t k } -1 = ϕ -t -k , k = 0, 1.
Thus, for (s, σ) ∈ o T * Ω, at most one of the points x -t -k (s, σ) and x -t -l (s, σ) is in Ω. Consequently, the contribution of cross terms between different Gaussian beams in (3.12) vanishes in o T * Ω, and we need to compute only the limits when ε goes to zero of the following two distributions:

µ t ε,k = c(s) 2 |σ| 2 w ε [I ε (Π k p ε,k k , Φ k )], k = 0, 1. (3.13) Remember that p ε,k k = a k k ε -1 T ε u I ε,r0,r∞ k + a ′ k k T ε v I ε,r0,r∞ k
, so µ t ε,k may be written as

µ t ε,k =c 2 (s)|σ| 2 w ε I ε (Π k a k k ε -1 T ε u I ε,r0,r∞ k , Φ k ) + w ε I ε (Π k a k k T ε v I ε,r0,r∞ k , Φ k ) -ic(s)|σ|w ε I ε (Π k a k k ε -1 T ε u I ε,r0,r∞ k , Φ k ), I ε (Π k a k k T ε v I ε,r0,r∞ k , Φ k ) + ic(s)|σ|w ε I ε (Π k a k k T ε v I ε,r0,r∞ k , Φ k ), I ε (Π k a k k ε -1 T ε u I ε,r0,r∞ k , Φ k ) . (3.14)
In the remainder of this Section we prove the following Proposition, compute µ t ε,k and the limit when ε → 0 of the microlocal energy density of u appr ε,r0,r∞ . Proposition 2 Let (κ ε ), (τ ε ) be uniformly bounded sequences in L 2 (R n ). Then

w ε I ε (Π k a k k T ε κ ε k , Φ k ), I ε (Π k a k k T ε τ ε k , Φ k ) ≈ Π 2 k w ε (κ ε , τ ε ) o ϕ t k -1 in o T * Ω.
Above ϕ t k is extended outside B as the identity. Proof: We simplify the integral

A(Φ k , Φ k ) R 2n Π k T ε κ ε k + Π k T ε τ ε k - e iΘε(Φ k ,Φ k ) drdδ
obtained when applying Lemma 3.1 in o T * Ω by firstly computing the phase Θ ε and the amplitude A and then analyzing the transported FBI transforms. Computa-

tion of Θ ε (Φ k , Φ k ) and A(Φ k , Φ k ). We consider (s, σ) ∈ o T * Ω and start from Θ ε (Φ k , Φ k )(s, σ, r, δ) = -2σ • r/ √ ε + i(r, δ) • Q Λ k k (t, s, σ), Λk k (t, s, σ) (r, δ).
The particular form of Λ k (t) = -iV t k (U t k ) -1 , see Lemma 2.1, induces a similar form for the matrix

Q Λ k k (t), Λk k (t) Q Λ k k (t), Λk k (t) Y t k = -iZ t k ,
where Y t k and Z t k are the 2n × 2n matrices

Y t k =   Ū t k k U t k k V t k k V t k k   and Z t k =   -V t k k V t k k Ū t k k -U t k k   .
Replacing U t k and V t k by their definitions links Y t k and Z t k to the Jacobian matrix

F t k Y t k = -i F t k k J -Id Id iId iId and Z t k = J F t k k -Id Id iId iId , so that Q Λ k k (t), Λk k (t) = -J F t k k J F t k k -1
.

As ϕ t k o ϕ -t -k = Id, one has F t k k F -t -k = Id.
Combining this relation with the symplecticity of F t k , one gets the following relation for the matrix

Q Λ k k (t), Λk k (t) Q Λ k k (t), Λk k (t) = (F -t -k ) T F -t -k . Therefore Θ ε (Φ k , Φ k )(s, σ, r, δ) = -2σ • r/ √ ε + i F -t -k (s, σ)(r, δ) 2 .

Moving to the amplitude

A(Φ k , Φ k ) = c 2 n 2 5n 2 π n 2 det( Λ k k + Λk k ) - 1 2
, one gets by using (2.18) and (2.21)

Λ k (t) + Λk (t) = 2 ( Ū t k ) -1 T (U t k ) -1 .
Hence

A(Φ k , Φ k ) = c 2 n 2 2n π n 2 det U t k k .
Plugging the form of the incident and reflected amplitudes in Lemma 2.1 and using the C 1 smoothness of a

( ′ )
k on B yields by Lemmas 3.1 and 3.3

w ε I ε (Π k a k k T ε κ ε k , Φ k ), I ε (Π k a k k T ε τ ε k , Φ k ) ≈c 2 n 2 2n π n 2 R 2n Π k T ε κ ε k + Π k T ε τ ε k - e -i2σ•r/ √ ε-(F -t -k (r,δ)) 2 drdδ =:J t ε,k (κ ε , τ ε ).
Analysis of the transported FBI transforms. It remains to analyze the most difficult terms in the amplitude, which involve transported FBI transforms

Π k T ε κ ε k + = Π k T ε κ ε o ϕ -t -k (s + √ εr, σ + √ εδ),
and

Π k T ε τ ε k - = Π k T ε τ ε o ϕ -t -k (s - √ εr, σ - √ εδ).
Let φ be a test function in

C ∞ 0 ( o T * Ω, R) and ϑ -t -k a map of C ∞ 0 (R 2n , R 2n ) that coincides with ϕ -t -k on K k z,θ (t) ∪ suppφ (see Theorem 1.4.1 of [21]
). We use Taylor's formula for this map to get for (s ± √ εr, σ ± √ εδ) ∈ K k z,θ (t) and (s, σ) ∈ suppφ

x -t -k ± = x -t -k ± √ ε D y x -t -k r ± √ ε D η x -t -k δ + εr x± ε , ξ -t -k ± = ξ -t -k ± √ ε D y ξ -t -k r ± √ ε D η ξ -t -k δ + εr ξ± ε , with r x± ε (s, σ, r, δ) = |α|=2 2 α! (r, δ) α 1 0 (1 -u)∂ α y ϑ -t -k (s, σ) ± u √ ε (r, δ) du, r ξ± ε (s, σ, r, δ) = |α|=2 2 α! (r, δ) α 1 0 (1 -u)∂ α η ϑ -t -k (s, σ) ± u √ ε (r, δ) du.
The change of variables (r

′ , δ ′ ) = F -t -k (s, σ)(r, δ) in J t ε,k (κ ε , τ ε )(s, σ) is thus appro- priate. Notice that for (s, σ) ∈ o T * Ω one has the following relations [26] D y x -u -k (s, σ) T ξ -u -k (s, σ) -σ = 0 and D η x -u -k (s, σ) T ξ -u -k (s, σ) = 0 for u ∈ R.
In fact, one can show that the derivatives of the previous equations w.r.t. u are zero. Besides, the equalities clearly hold true at u = 0 for k = 0, and at u = T k (s, σ) for k = ±1, as a consequence of (4.9). Hence, it follows that r,δ). In order to use the change of variables (s, σ) = ϕ t k (y, η) for < J t ε,k (κ ε , τ ε ), φ >, we extend ϕ t k outside B by the identity and still denote it ϕ t k , making ϕ t k a one to one map from R 2n to ϕ t k (R 2n ). Then Π k o ϕ t k and φ o ϕ t k belong to C ∞ 0 (R 2n , R) and are supported in B. Expanding the FBI transforms gives

σ • r = ξ -t -k (s, σ) • D y x -t -k (s, σ)r + D η x -t -k (s, σ)δ = ξ -t -k (s, σ) • r ′ , which leads in o T * Ω to J t ε,k (κ ε , τ ε ) =c 2 n 2 2n π n 2 R 2n (Π k ) + T ε κ ε (x -t -k + √ εr ′ + εr x+ ε ′ , ξ -t -k + √ εδ ′ + εr ξ+ ε ′ ) (Π k ) -T ε τ ε (x -t -k - √ εr ′ + εr x- ε ′ , ξ -t -k - √ εδ ′ + εr ξ- ε ′ ) e -2iξ -t -k •r ′ / √ ε-r ′2 -δ ′2 dr ′ dδ ′ , where (r x ε ′ , r ξ ε ′ )(s, σ, r ′ , δ ′ ) = (r x ε , r ξ ε )(s, σ,
< J t ε,k (κ ε , τ ε ), φ > =c 4 n 2 2n π n 2 ε -3n 2 R 6n φ o ϕ t k (y, η)κ ε (z)τ ε (z ′ ) (Π k o ϕ t k )(y + √ εr ′ + εR x+ ε , η + √ εδ ′ + εR ξ+ ε ) (Π k o ϕ t k )(y - √ εr ′ + εR x- ε , η - √ εδ ′ + εR ξ- ε ) e iη•(2 √ εr ′ +εR x+ ε -εR x- ε -z+z ′ )/ε+iδ ′ •(2y-z-z ′ +εR x+ ε +εR x- ε )/ √ ε e iR ξ+ ε •(y+ √ εr ′ +εR x+ ε -z)-iR ξ- ε •(y- √ εr ′ +εR x- ε -z ′ ) e -(y+ √ εr ′ +εR x+ ε -z) 2 /(2ε)-(y- √ εr ′ +εR x- ε -z ′ ) 2 /(2ε) e -2iη•r ′ / √ ε-r ′2 -δ ′2 dr ′ dδ ′ dzdz ′ dydη, where (R x ε , R ξ ε )(y, η, r ′ , δ ′ ) = (r x ε ′ , r ξ ε ′ )(s, σ, r ′ , δ ′ ).
We perform the changes of variables

(x, u) = ( z + z ′ 2 , z -z ′ ε ) and y ′ = (y - z + z ′ 2 )/ √ ε to obtain < J t ε,k (κ ε , τ ε ), φ > =c 4 n 2 2n π n 2 R 6n κ ε (x + ε 2 u)τ ε (x - ε 2 u)d ε e iγε-iη•u dr ′ dδ ′ dxdudy ′ dη,
where

d ε (x, y ′ , η, r ′ , δ ′ ) =φ o ϕ t k (x + √ εy ′ , η)(Π k o ϕ t k )(x + √ εy ′ + √ εr ′ + εR x+ ε ′ , η + √ εδ ′ + εR ξ+ ε ′ ) (Π k o ϕ t k )(x + √ εy ′ - √ εr ′ + εR x- ε ′ , η - √ εδ ′ + εR ξ- ε ′ ), γ ε (x, y ′ , η, r ′ , δ ′ , u) =η • (R x+ ε ′ -R x- ε ′ ) + δ ′ • (2y ′ + √ εR x+ ε ′ + √ εR x- ε ′ ) + √ εR ξ+ ε ′ • (y ′ + r ′ + √ εR x+ ε ′ - √ ε u 2 ) - √ εR ξ- ε ′ • (y ′ -r ′ + √ εR x- ε ′ + √ ε u 2 ) + ir ′ 2 + iδ ′ 2 + i(y ′ + r ′ + √ εR x+ ε ′ - √ εu/2) 2 /2 + i(y ′ -r ′ + √ εR x- ε ′ + √ εu/2) 2 /2, and (R x ε ′ , R ξ ε ′ )(x, y ′ , η, r ′ , δ ′ ) = (R x ε , R ξ ε )(x + √ εy ′ , η, r ′ , δ ′ ).
Notice that d ε (x, y ′ , η, r ′ , δ ′ ) converges when ε → 0 to ,η). On the other hand, since εr ±

d 0 (x, η) = φ o ϕ t k (x, η)(Π k o ϕ t k ) 2 (x
x are the remainder terms in the Taylor expansions of

x -t -k (s ± √ εr, σ ± √ εδ) at order 2, r x+ ε -r x- ε is of order √ ε and so is R x+ ε ′ -R x- ε ′ , leading to γ ε (x, y ′ , η, r ′ , δ ′ , u) → ε→0 γ 0 (y ′ , r ′ , δ ′ ) = 2δ ′ • y ′ + iy ′ 2 + 2ir ′ 2 + iδ ′ 2 .
One has

< J t ε,k (κ ε , τ ε ), φ > -c 4 n 2 2n π n 2 R 6n κ ε (x + ε 2 u)τ ε (x - ε 2 u)d 0 e iγ0 e -iη•u dr ′ dδ ′ dudy ′ dxdη R 4n R n |κ ε |(x + ε 2 u)|τ ε |(x - ε 2 u)dx sup x F η d ε e iγε -d 0 e iγ0 (x, y ′ , u, r ′ , δ ′ , u) dr ′ dδ ′ dudy ′ . (3.15)
Cauchy-Schwartz inequality w.r.t. dx insures that the bracket integral is less than

κ ε L 2 τ ε L 2 . Let us examine the term R 4n sup x F η d ε e iγε -d 0 e iγ0 (x, y ′ , u, r ′ , δ ′ , u) dr ′ dδ ′ dudy ′ .
For fixed y ′ , r ′ , δ ′ the functions d ε and d 0 are compactly supported w.r.t. (x, η) so On the other hand, successive integrations by parts give

sup x F η d ε e iγε -d 0 e iγ0 (x, y ′ , u, r ′ , δ ′ , u) sup (x,η) d ε e iγε -d 0 e iγ0 (x, y ′ , η, r ′ , δ ′ , u) .
R n d ε e iγε e -iη•u dη = (1 + u 2 ) -n R n L d ε e iγε e -iη•u dη,
with L a differential operator w.r.t. η, of order 2n. Thus,

sup x F η d ε e iγε (x, y ′ , u, r ′ , δ ′ , u) (1 + u 2 ) -n sup (x,η) max |α|≤2n ∂ α η d ε e iγε (x, y ′ , η, r ′ , δ ′ , u) , (3.16)
for every y ′ , r ′ , δ ′ , u. The quantities (x + √ εy ′ , η) and

√ ε(r ′ , δ ′ ) are bounded on the support of d ε , so R x± ε ′ , R ξ± ε ′
and their derivatives w.r.t. η are dominated by (r ′ , δ ′ ) 2 . Hence for a given multiindex α, there exists C > 0 s.t.

∂ α η d ε ≤C, ∂ α η γ ε ≤C|(r ′ , δ ′ )| |(r ′ , δ ′ )| + |y ′ + r ′ + √ εR x+ ε ′ - √ εu/2| + |y ′ -r ′ + √ εR x- ε ′ - √ εu/2| if |α| ≥ 1,
for all (x, y ′ , η, r ′ , δ ′ ) ∈ suppd ε and u ∈ R n . Thus, there exists C, C ′ > 0 s.t.

∂ α η d ε e iγε ≤ Ce -C ′ (y ′ +r ′ + √ εR x+ ε ′ - √ εu/2) 2 -C ′ (y ′ -r ′ + √ εR x- ε ′ + √ εu/2) 2 -C ′ r ′2 -C ′ δ ′ 2 ≤ Ce -C ′ (2y ′ + √ εR x+ ε ′ + √ εR x- ε ′ ) 2 -C ′ r ′2 -C ′ δ ′ 2 ,
for all (x, y ′ , η, r ′ , δ ′ ) ∈ suppd ε and u ∈ R n . On the support of d ε ,

√ εR x± ε ′ are dominated by |(r ′ , δ ′ )|, which implies for some C 0 > 0 that (2y ′ + √ εR x+ ε ′ + √ εR x- ε ′ ) 2 ≥ 4y ′ 2 -C 0 |(r ′ , δ ′ )||y ′ |. Hence, if |y ′ | ≥ C 0 |(r ′ , δ ′ )|, e -C ′ (2y ′ + √ εR x+ ε ′ + √ εR x- ε ′ ) 2 ≤ e -C ′′ y ′2 . Otherwise, e -C ′ r ′2 -C ′ δ ′ 2 ≤ e -C ′′ y ′2 -C ′′ r ′2 -C ′′ δ ′2
. In all cases, there exists C ′ , C ′′ > 0 s.t.

∂ α η d ε e iγε ≤ C ′ e -C ′′ y ′2 -C ′′ r ′2 -C ′′ δ ′2 ,
for every x, y ′ , η, r ′ , δ ′ , u and ε ∈]0, ε 0 ] with some ε 0 > 0. Using this in (3.16) leads to

sup x F η d ε e iγε (x, y ′ , u, r ′ , δ ′ , u) (1 + u 2 ) -n e -Cy ′2 -Cr ′2 -Cδ ′ 2 ,
and repeating the same arguments for sup

x

F η d 0 e iγ0 gives sup x F η d ε e iγε -d 0 e iγ0 (x, y ′ , u, r ′ , ′ , u) (1 + u 2 ) -n e -Cy ′2 -Cr ′2 -Cδ ′ 2 ,
for every y ′ , u, r ′ , δ ′ and ε ∈]0, ε 0 ]. By the dominated convergence theorem, one obtains

R 4n sup x |F η d ε e iγε -d 0 e iγ0 (x, y ′ , u, r ′ , δ ′ , u)|dy ′ dudr ′ dδ ′ → ε→0 0.
From the inequality (3.15) concerning the distribution J t ε,k (κ ε , τ ε ), one finally has by plugging the expressions of d 0 and γ 0

< J t ε,k (κ ε , τ ε ), φ >= c 4 n 2 2n π n 2 R 6n κ ε (x + ε 2 u)τ ε (x - ε 2 u) e 2iδ ′ •y ′ -y ′2 -2r ′ 2 -δ ′ 2 e -iη•u dr ′ dδ ′ dxdudy ′ dη + o(1).
Integration w.r.t. r ′ , δ ′ , y ′ , η yields

< J t ε,k (κ ε , τ ε ), φ > =(2π) -n R 2n F η Π 2 k o ϕ t k φ o ϕ t k (x, u)κ ε (x + ε 2 u)τ ε (x - ε 2 u)dxdu + o(1).
The integral in the r.h.s. is exactly the Wigner transform of (κ ε , τ ε ) tested on

Π 2 k o ϕ t k φ o ϕ t k . 2
We are now able to compute the measure µ t ε,k given in (3.14) by using the previous Proposition and the Lemma 4.7

µ t ε,k ≈ Π 2 k w ε v I ε,r0,r∞ -ic|D|u I ε,r0,r∞ o ϕ t k -1 in o T * Ω.
Recalling the relation between the Wigner measure and the FBI transform (see Proposition 1.4 of [START_REF] Gérard | Ergodic properties of eigenfunctions for the Dirichlet problem[END_REF])

R 2n |T ε a ε | 2 θdydη → ε→0 < w[a ε ], θ >, (3.17)
for θ ∈ C ∞ 0 (R 2n , R) and (a ε ) a uniformly bounded sequence in L 2 (R n ), it follows that w ε v I ε,r0,r∞ -ic|D|u I ε,r0,r∞ ≈ 0 in (K y × K η ) c or equivalently

w ε v I ε,r0,r∞ -ic|D|u I ε,r0,r∞ o ϕ t k -1 ≈ 0 in (K k z,θ (t)) c . Since Π k ≡ 1 on K k z,θ (t), one deduces µ t ε,k ≈ w ε v I ε,r0,r∞ -ic|D|u I ε,r0,r∞ o ϕ t k -1 in o T * Ω.
By summing over k = 0, 1 and letting ε → 0, we get

w[v + t (t, .)] = k=0,1 w v I ε,r0,r∞ -ic|D|u I ε,r0,r∞ o ϕ t k -1 in o T * Ω.
For u ∈ [-T, T ] and (y, η) ∈ K y × (R n \{0}), the incident and reflected flows are related to the broken bicharacteristic flow associated to -i∂ t -c|D| as follows:

ϕ u b (y, η) =    ϕ u -1 (y, η) if u < T -1 (y, η), ϕ u 0 (y, η) if T -1 (y, η) < u < T 1 (y, η), ϕ u 1 (y, η) if u > T 1 (y, η).
We extend ϕ u b at times of reflections arbitrary. We define ϕ u b in (Ω\K y ) × (R n \{0}) by successively reflecting the rays at the boundary.As only one incident/reflected ray can be in the interior of the domain at a fixed time t ∈ [-T, T ]

φ o ϕ t b = k=0,1 φ o ϕ t k in K y × R n \{0}.
It follows that

w[v + t (t, .)] = w v I ε,r0,r∞ -ic|D|u I ε,r0,r∞ o ϕ t b -1 in o T * Ω.
The computations for v - t are similar. One has just to replace the index k = 1 by k = -1 and p ε,k k by q ε,k k in (3.13) and to repeat the same techniques. If we denote Υ ± ε,r0,r∞ = v I ε,r0,r∞ ± ic|D|u I ε,r0,r∞ , then one gets

w[v - t (-t, .)] = w Υ + ε,r0,r∞ o ϕ -t b -1 in o T * Ω.
Using these results in (3.11) as ε → 0 leads to

E u appr ε,r0,r∞ (t, .) = 1 2 w Υ + ε,r0,r∞ o ϕ -t b -1 + 1 2 w Υ - ε,r0,r∞ o ϕ t b -1 in o T * Ω.
(3.18)

Proof of the main Theorem

A consequence of the estimate (3.10) is

| < w(a ε , b ε ), θ > | lim sup ε→0 a ε L 2 (Ω) lim sup ε→0 b ε L 2 (Ω) , (3.19)
for (a ε ), (b ε ) uniformly bounded sequences in

L 2 (R n ) and θ ∈ C ∞ 0 ( o T * Ω, R).
Applying this estimate to the difference between the derivatives of the exact and approximate solutions of the IBVP (1.1a)-(1.1b) with initial conditions (1.1c'), one deduces the measures associated to ∂ t u ε,r0,r∞ and ∂ x u ε,r0,r∞ and gets by (3.18)

E (u ε,r0,r∞ (t, .)) = 1 2 w Υ + ε,r0,r∞ o ϕ -t b -1 + 1 2 Υ - ε,r0,r∞ o ϕ t b -1 in o T * Ω.
Remark 1 Gaussian beam summation of first order beams allows to compute the microlocal energy density of the solution of the IBVP (1.1) as ε → 0, under the hypotheses (H1),(H4) and (H5) on initial conditions. Summation of higher order beams may imply asymptotic formulas for the Wigner transforms and thus for the energy density. Higher order terms in the expansion of the Wigner transform were studied for instance in [START_REF] Filippas | Semiclassical Wigner function and geometrical optics[END_REF] and [START_REF] Pulvirenti | Semiclassical expansion of Wigner functions[END_REF] for WKB initial data.

Let us now study the microlocal energy density for the problem (1.1) when ε → 0, by making the data (u

I ε,r0,r∞ , v I ε,r0,r∞ ) approach (u I ε , v I ε ). The contribution of the sets {η ∈ R n , |η| ≥ r ∞ /4} and {η ∈ R n , |η| ≤ 4r 0 } where γ r0,r∞ ≡ 1 (remember the definition of γ r0,r∞ in (2.11)) to T ε u I ε , T ε v I ε is controlled asymptotically by the assumptions (H2) and (H3). Set Υ ± ε = v I ε ± ic|D|u I ε and denote φ t = φoϕ t b . Then φ t ∈ C ∞ 0 (R 2n , R
) and one has

< E (u ε (t, .)) , φ > - 1 2 < w Υ + ε , φ -t > - 1 2 < w Υ - ε , φ t > ≤ 1 2 < w Υ + ε,r0,r∞ -w Υ + ε , φ -t > + 1 2 < w Υ - ε,r0,r∞ -w Υ - ε , φ t > + < w ∂ t u ε (t, .) -w ∂ t u ε,r0,r∞ (t, .) , φ > + n j=1 < w c∂ xj u ε (t, .) -w c∂ xj u ε,r0,r∞ (t, .) , φ > + < E (u ε,r0,r∞ (t, .)) , φ > - 1 2 < w Υ + ε,r0,r∞ , φ -t > - 1 2 < w Υ - ε,r0,r∞ , φ t > . (3.20)
We use (3.10) to get

| < w Υ + ε,r0,r∞ -w Υ + ε , φ -t > | lim sup ε→0 Υ + ε,r0,r∞ -Υ + ε L 2 (R n ) lim sup ε→0 Υ + ε,r0,r∞ L 2 (R n ) + Υ + ε L 2 (R n ) lim sup ε→0 v I ε -v I ε,r0,r∞ L 2 (Ω) + lim sup ε→0 u I ε -u I ε,r0,r∞ H 1 (Ω) .
Similarly, by (3.19)

< w ∂ t u ε (t, .) -w ∂ t u ε,r0,r∞ (t, .) , φ > lim sup ε→0 ∂ t u ε (t, .) -∂ t u ε,r0,r∞ (t, .) L 2 (Ω) lim sup ε→0 ∂ t u ε (t, .) L 2 (Ω) + lim sup ε→0 ∂ t u ε,r0,r∞ (t, .) L 2 (Ω) ,
and for j = 1, . . . , n

< w ∂ xj u ε (t, .) -w ∂ xj u ε,r0,r∞ (t, .) , φ > lim sup ε→0 ∂ xj u ε (t, .) -∂ xj u ε,r0,r∞ (t, .) L 2 (Ω) lim sup ε→0 ∂ xj u ε (t, .) L 2 (Ω) + lim sup ε→0 ∂ xj u ε,r0,r∞ (t, .) L 2 (Ω) .
The solution of the IBVP for the wave equation is given by a continuous unitary evolution group on the space H 1 (Ω, dx) × L 2 (Ω, dx). Hence

∂ t u ε (t, .) -∂ t u ε,r0,r∞ (t, .) L 2 (Ω) v I ε -v I ε,r0,r∞ L 2 (Ω) + u I ε -u I ε,r0,r∞ H 1 (Ω) , ∂ xj u ε (t, .) -∂ xj u ε,r0,r∞ (t, .) L 2 (Ω) v I ε -v I ε,r0,r∞ L 2 (Ω) + u I ε -u I ε,r0,r∞ H 1 (Ω) ,
for j = 1, . . . , n. Finally, by using (3.18), the estimate (3.20) is simplified into

| < E (u ε (t, .)) , φ > - 1 2 < w Υ + ε , φ -t > - 1 2 < w Υ - ε , φ t > | lim sup ε→0 v I ε -v I ε,r0,r∞ L 2 (Ω) + lim sup ε→0 u I ε -u I ε,r0,r∞ H 1 (Ω) . (3.21)
We therefore need to estimate the difference between initial data (1.1c) and (1.1c'). We start by the initial speed. By the exponential decrease of T * ε γ r0,r∞ T ε v I ε on the support of 1ρ (see (2.12)), one has

v I ε -v I ε,r0,r∞ L 2 (Ω) ε ∞ + v I ε -T * ε γ r0,r∞ T ε v I ε L 2 (Ω) . Because T * ε is bounded on L 2 (R 2n ) → L 2 (R n ) and T * ε T ε = Id v I ε -T * ε γT ε v I ε L 2 (R n ) ≤ (1 -χ r∞/2 )T ε v I ε L 2 (R 2n ) x + χ r∞/2 χ 4r0 T ε v I ε L 2 (R 2n ) y .
Firstly, Lemma 4.2 yields

(x) 2 = c n (2π) -n 2 ε -n 4 1 -χ r∞/2 (η) R 2n F v I ε (ξ)e iξ•y-(η-εξ) 2 /(2ε) dξ 2 L 2 (R 2n y,η ) .
It follows by Parseval equality that

(x) 2 = c 2 n ε -n 2 |εξ|≤r∞/8 1 -χ r∞/2 (η) 2 |F v I ε (ξ)| 2 e -(η-εξ) 2 /ε dξdη +c 2 n ε -n 2 |εξ|≥r∞/8 1 -χ r∞/2 (η) 2 |F v I ε (ξ)| 2 e -(η-εξ) 2 /ε dξdη.
The first integral in the r.h.s. is exponentially decreasing, which leads to lim sup ε→0

x lim sup

ε→0 |εξ|≥r∞/8 |F v I ε (ξ)| 2 dξ 1 2
.

Secondly, as dist suppv

I ε , supp(1 -ρ) > 0, one gets (1 -ρ)T ε v I ε L 2 (R 2n
) ≤ e -C/ε by Lemma 4.3 and thus lim sup

ε→0 ( y) 2 = lim sup ε→0 ρ(y)χ r∞/2 (η)χ 4r0 (η)T ε v I ε 2 L 2 (R 2n ) .
It results from the relation (3.17) applied with

a ε = v I ε that (y) 2 → ε→0 < w v I ε , ρ 2 ⊗ χ 2 r∞/2 χ 2 4r0 > .
Because w v I ε is a regular measure, assumption (H3) yields

∀α > 0, ∃l 0 (α) > 0 s.t. w v I ε ({|ξ| ≤ l 0 (α)}) ≤ α.
One deduces, for 4r 0 ≤ l 0 (α), that lim sup

ε→0 y √ α,
which leads to lim sup

ε→0 v I ε -v I ε,r0,r∞ L 2 (Ω) lim sup ε→0 |εξ|≥r∞/8 |F v I ε (ξ)| 2 dξ 1 2 + √ α.
For the analysis of u I εu I ε,r0,r∞ in H 1 (Ω), we begin by estimating the spatial derivatives of the difference. It follows by using the relation (2.13) when differentiating the inverse FBI transform that

∂ xj u I ε -∂ xj u I ε,r0,r∞ = ∂ xj u I ε -(∂ xj ρ)T * ε γ r0,r∞ T ε u I ε -ρT * ε γ r0,r∞ ∂ yj T ε u I ε .
< E (u ε (t, .)) , φ > - 1 2 < w Υ + ε , φ -t > - 1 2 < w Υ - ε , φ t > √ α + lim sup ε→0 ε|ξ|≥r∞/8 |F v I ε (ξ)| 2 dξ 1 2 + n j=1 lim sup ε→0 ε|ξ|≥r∞/8 |F ∂ xj u I ε (ξ)| 2 dξ 1 2 . (3.22)
The assumption (H2) of ε-oscillation means by definition that lim sup Since the l.h.s. of the estimate (3.22) does not depend on α nor r ∞ , one deduces by taking the limits α → 0 and r ∞ → ∞ that

ε→0 ε|ξ|≥R |F v I ε (ξ)| 2 dξ → R→+∞ 0, ( 3
E (u ε (t, .)) = 1 2 w Υ + ε o ϕ -t b -1 + 1 2 w Υ - ε o ϕ t b -1 in o T * Ω.

Appendix

4.1 Reflected first order and higher order beams

Higher order beams

Higher order beams, possibly with more than one amplitude, can be constructed to satisfy better interior and boundary estimates. In this case, the eikonal equation (2.2) must be satisfied up to order R ≥ 2 on the rays. If r ≥ 3, the equations

∂ α x (p(x, ∂ t ψ, ∂ x ψ)) (t, x t ) = 0, |α| = r, (4.1)
give systems of linear ODEs of order 1 on (∂ α x ψ(t, x t )) |α|=r with second members involving lower order spatial derivatives of the phase. In fact, the key observation is the equality

∂ τ p(ϕ t )∂ t ∂ α x ψ(t, x t ) + ∂ ξ p(ϕ t ) • ∂ x ∂ α x ψ(t, x t ) =2c(x t )|ξ t |∂ t ∂ α x ψ(t, x t ) + 2c 2 (x t )ξ t • ∂ x ∂ α x ψ(t, x t ) =2c(x t )|ξ t | d dt ∂ α x ψ(t, x t ) ,
used for |α| = r to eliminate the r + 1-th order derivatives of ψ in equation (4.1).

To summarize, the requirements

∂ t ψ(t, x t ) = -c(x t )|ξ t |, ∂ x ψ(t, x t ) = ξ t , p (x, ∂ t ψ(t, x), ∂ x ψ(t, x)) = 0 on x = x t up to order R,
uniquely determine the spatial derivatives of ψ on the ray up to the order R under the knowledge of their initial values on (0, x 0 ). We refer to [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF] for further details.

A general relation between incident and reflected beams phases

By (2.19), the Hessian matrix of the incident beam's phase is related to the Jacobian matrix of the incident flow. One can prove that its higher order derivatives are also related to the higher order derivatives of the incident flow. Computations exhibiting such relations can be found for instance in the Appendix of [START_REF] Norris | Elastic Gaussian wave packets in isotropic media[END_REF]. We shall give a nice relation between an incident phase ψ inc and the associated reflected phase ψ ref for beams of any order. This relation is intuitive true on geometrical grounds and it provides with the derivatives of the reflected phase up to order R, which might be useful in applications of Gaussian beams. Consider the following auxiliary function linking ϕ t 1 to ϕ t 0 for any fixed time t

s 1 : B → B (x, ξ) → ϕ -T1(x,ξ) 0 o R o ϕ T1(x,ξ) 0 (x, ξ).
For a given point (x, ξ) ∈ B, s 1 (x, ξ) is its "image by the mirror" ∂Ω. For instance, Chazarain used this type of auxiliary functions in [START_REF] Chazarain | Paramétrix du problème mixte pour l'équation des ondes à l'intérieur d'un domaine convexe pour les bicaractéristiques[END_REF] to show propagation of regularity for wave type equations in a convex domain.

By the Implicit functions theorem, T 1 is C ∞ on the open set B and so is s 1 . Since ϕ t 0 o s 1 satisfies the same Hamiltonian equations as ϕ t 1 and ϕ

T1(x,ξ) 1 (x, ξ) = ϕ T1(x,ξ) 0 o s 1 (x, ξ) for (x, ξ) ∈ B, one has ϕ t 1 = ϕ t 0 o s 1 .
Besides, noticing that T 1 (ϕ t 0 ) = T 1t, one has also

ϕ t 1 = s 1 o ϕ t 0 . ( 4
s 1 = (r, λ). For m ∈ N, f , g functions in C ∞ R n u × (R n ξ \{0}), C p , u 0 ∈ R n a fixed point and V ∈ C ∞ (R n u , C n ξ ) a phase function s.t. V (u 0 ) ∈ R n ξ \{0}, we introduce the notation f (u, V (u)) m ≍ u=u0 g (u, V (u)) ,
to denote that the formal partial derivatives of f (u, V (u)) and g (u, V (u)) up to the order m coincide on u 0 . The differentiation here is viewed formally, since V may be complex valued out of u 0 , which makes f (u, V (u)) and g(u, V (u)) not defined for u = u 0 . However, on the exact point u 0 , one can always use the formula of composite functions derivatives to get a formal expression of the derivatives. We will use the same notation

f (t, x, V (t, x)) m ≍ x=x t g (t, x, V (t, x)) , for functions f, g ∈ C ∞ R t × R n x × (R n ξ \{0}), C p and phase function V ∈ C ∞ (R t × R n x , C n ξ ) s.t. for t ∈ R, V (t, x t ) ∈ R n ξ \{0}
to denote that the formal partial derivatives of f (t, x, V (t, x)) and g (t, x, V (t, x)) w.r.t. x up to order m coincide on (t, x t ) for all t ∈ R. We will be sloppy with respect to the notation of the dependence of the phase V on its variables.

Consider an integer R ≥ 2 and an incident phase ψ inc satisfying

∂ t ψ inc (t, x t 0 ) = -c(x t 0 )|ξ t 0 |, ∂ x ψ inc (t, x t 0 ) = ξ t 0 and p(x, ∂ t ψ inc , ∂ x ψ inc ) R ≍ x=x t 0 0.
As a particular case, the phase ψ 0 is obtained by setting R = 2 and choosing its initial value on the ray as zero and its initial Hessian matrix on the ray as iId.

Let

ψ ref ∈ C ∞ (R t × R n
x , C) be the reflected phase associated to ψ inc , that is the phase satisfying

∂ t ψ ref (t, x t 1 ) = -c(x t 1 )|ξ t 1 |, ∂ x ψ ref (t, x t 1 ) = ξ t 1 and p(x, ∂ t ψ ref , ∂ x ψ ref ) R ≍ x=x t 1 0,
and having the same time and tangential derivatives as ψ inc at the instant and the point of reflection (T 1 , x T1 0 ) up to the order R. Since ϕ t 0 and the reflection R conserve c(x)|ξ| (see (2.9)), one has for every (x, ξ) ∈ B and τ ∈ R * p (r(x, ξ), τ, λ(x, ξ)) = p(x, τ, ξ).

Thus p (r(x, ∂ x ψ inc ), ∂ t ψ inc , λ(x, ∂ x ψ inc )) ∞ ≍ x=x t 0 p(x, ∂ t ψ inc , ∂ x ψ inc ),
which implies, by construction of

ψ inc p (r(x, ∂ x ψ inc ), ∂ t ψ inc , λ(x, ∂ x ψ inc )) R ≍ x=x t 0 0. (4.3)
Compare this with the equation

p (r(x, ∂ x ψ inc ), ∂ t ψ ref (t, r(x, ∂ x ψ inc )) , ∂ x ψ ref (t, r(x, ∂ x ψ inc ))) R ≍ x=x t 0 0
resulting from the construction of ψ ref and (4.2). This suggests the following Lemma Lemma 4.1

∂ t ψ ref (t, r(x, ∂ x ψ inc )) R-1 ≍ x=x t 0 ∂ t ψ inc and ∂ x ψ ref (t, r(x, ∂ x ψ inc )) R-1 ≍ x=x t 0 λ(x, ∂ x ψ inc ).
A similar result linking the reflected phase associated to the ray (t, x -t -1 ) to ψ inc can be established. Proof: The strategy of the proof is the following: we consider a phase function θ satisfying the relations announced in Lemma 4.1 and we prove that θ fulfills the eikonal equation on the reflected ray up to order R and has the correct derivatives at the instant and point of reflection. This proves that θ coincides with the reflected phase on the reflected ray up to the order R.

Denote r (x, ∂ x ψ inc (t, x)) by ̺(t, x) or simply by ̺ if no confusion arises and let us first verify that for a fixed k ≥ 1 there exists a phase function

θ ∈ C ∞ (R t ×R n x , C) s.t. ∂ x θ(t, ̺) k ≍ x=x t 0 λ(x, ∂ x ψ inc ). (4.4) Let A(t, x, ξ) = D x r(x, ξ) + D ξ r(x, ξ)∂ 2 x ψ inc (t, x) and B(t, x, ξ) = D x λ(x, ξ) + D ξ λ(x, ξ)∂ 2 x ψ inc (t, x). Then D x ̺(t, x) = A(t, x, ∂ x ψ inc ), D x [λ(x, ∂ x ψ inc (t, x))] = B(t, x, ∂ x ψ inc ) and for v ∈ C ∞ (R t × R n x , C p ) one has D x (v(t, ̺)) ∞ ≍ x=x t 0 D x v(t, ̺)A(t, x, ∂ x ψ inc ).
Hence, θ exists if A(t, x t 0 , ξ t 0 ) is non singular and

B(t, x, ∂ x ψ inc )A(t, x, ∂ x ψ inc ) -1 k-1 ≍ x=x t 0 A(t, x, ∂ x ψ inc ) T -1 B(t, x, ∂ x ψ inc ) T . (4.5)
From (4.2) one gets

A(t, x t 0 , ξ t 0 )(D y x t 0 + iD η x t 0 ) = D y x t 1 + iD η x t 1 .
Since ϕ t 1 is symplectic, the matrix This proves that D y x t 1 + iD η x t 1 is invertible and so is A(t, x t 0 , ξ t 0 ). On the other hand,

A B = D x r D ξ r D x λ D ξ λ Id ∂ 2 x ψ inc . Let M (x, ξ) = D x r(x, ξ) D ξ r(x, ξ) D x λ(x, ξ) D ξ λ(x, ξ) . Then [A T B -B T A] = Id ∂ 2 x ψ inc T M T JM Id ∂ 2 x ψ inc . Since M T JM = Ds T 1 JDs 1 , the symplecticity of s 1 leads to M T JM = J. Hence [A T B -B T A] = Id ∂ 2 x ψ inc T J Id ∂ 2
x ψ inc = 0 and the requirement (4.5) is fulfilled. The relation (4.4) fixes the derivatives of ∂ t ∂ x θ on (t, x t 1 ) up to order k -1. Indeed, using the compatibility condition

d dt f (t, x, ∂ x ψ inc (t, x)) | x=x t 0 =∂ t [f (t, x, ∂ x ψ inc (t, x))] | x=x t 0 + ∂ x [f (t, x, ∂ x ψ inc (t, x))] | x=x t 0 • ẋt 0 on the maps (t, x, ξ) → ∂ x θ (t, r(x, ξ)), (x, ξ) → λ(x, ξ)
and their derivatives yields recursively by (4.4)

∂ t [∂ x θ(t, ̺)] k-1 ≍ x=x t 0 D ξ λ(x, ∂ x ψ inc )∂ t ∂ x ψ inc . Thus ∂ t ∂ x θ(t, ̺) + ∂ 2 x θ(t, ̺)D ξ r(x, ∂ x ψ inc )∂ t ∂ x ψ inc k-1 ≍ x=x t 0 D ξ λ(x, ∂ x ψ inc )∂ t ∂ x ψ inc . Using the relations ∂ 2 x θ(t, ̺) k-1 ≍ x=x t 0 BA -1 (t, x, ∂ x ψ inc
) and (4.5) in the previous equation yields

∂ t ∂ x θ(t, ̺) k-1 ≍ x=x t 0 D ξ λ(x, ∂ x ψ inc ) -A T -1 B T (t, x, ∂ x ψ inc )D ξ r(x, ∂ x ψ inc ) ∂ t ∂ x ψ inc . Since A T D ξ λ -B T D ξ r = Id ∂ 2 x ψ inc T M T JM 0 Id = Id, it follows that A(t, x, ∂ x ψ inc ) T ∂ t ∂ x θ(t, ̺) k-1 ≍ x=x t 0 ∂ t ∂ x ψ inc . Note that ∂ x (u(t, ̺)) ∞ ≍ x=x t 0 A(t, x, ∂ x ψ inc ) T ∂ x u(t, ̺) for u ∈ C ∞ (R t × R n x , C), (4.6)
so one gets

∂ x (∂ t θ(t, ̺)) k-1 ≍ x=x t 0 ∂ t ∂ x ψ inc . Setting ∂ t θ(t, x t 1 ) = ∂ t ψ inc (t, x t 0 ) implies then that ∂ t θ(t, ̺) k ≍ x=x t 0 ∂ t ψ inc . ( 4
p(x, ∂ t θ, ∂ x θ) R ≍ x=x t 1 0.
To compare the time and tangential derivatives of θ and

ψ inc at (T 1 , x T1 0 ), let us introduce a C ∞ parametrization of a neighborhood U of x T1 0 in ∂Ω σ : N → R n , where N is an open subset of R n-1 , σ(N ) = U and σ a diffeomorphism from N to U. For x ∈ R n close to x T1 0 , we may write x = σ(v) + v n ν (σ(v)) , with v ∈ N and v n ∈ R. Denote σ(v 1 ) = x T1 0 and set θ b (t, v) = θ (t, σ(v)) and (ψ inc ) b (t, v) = ψ inc (t, σ(v)) the phases at the boundary near x T1 0 . Since r(X, Ξ) = X for (X, Ξ) ∈ o T * R n | ∂Ω , it follows that ̺ (t, σ(v)) ∞ ≍ (t,v)=(T1,v1) σ(v),
which implies by (4.7) that have the same time and tangential derivatives at (T 1 , v1 ) from the order 1 to the order k + 1.

∂ t θ b k ≍ (t,v)=(T1,v1) ∂ t (ψ inc ) b . Similarly λ(X, Ξ) = Ξ -2 (Ξ • ν(X)) ν(X) for (X, Ξ) ∈ o T * R n | ∂Ω , leading to Dσ(v) T λ (σ(v), ∂ x ψ inc (t, σ(v))) ∞ ≍ (t,v)=(T1,v1) Dσ(v) T ∂ x ψ inc (t, σ(v)) . (4.8) Since ∂ v θ b (t, v) = Dσ(v) T ∂ x θ (t, σ (
If we assume that θ(T

1 , x T1 0 ) = ψ inc (T 1 , x T1 0 ), then θ b k+1 ≍ (t,v)=(T1,v1) (ψ inc ) b ,
and θ satisfies all the requirements that determine the reflected phase associated to ψ inc and concentrated on (t, x t 1 ). The phases θ and ψ ref are thus equal on (t, x t 1 ) up to the order R.

∂ 2 x ψ 1 (t, x t 1 ) D x r(x t 0 , ξ t 0 ) + D ξ r(x t 0 , ξ t 0 )∂ 2 x ψ 0 (t, x t 0 ) =D x λ(x t 0 , ξ t 0 ) + D ξ λ(x t 0 , ξ t 0 )∂ 2 x ψ 0 (t, x t 0 ).

One obtains by plugging the expression

(2.19) of ∂ 2 x ψ 0 (t, x t 0 ) ∂ 2 x ψ 1 (t, x t 1 ) D x r(x t 0 , ξ t 0 )U t 0 + D ξ r(x t 0 , ξ t 0 )V t 0 = D x λ(x t 0 , ξ t 0 )U t 0 + D ξ λ(x t 0 , ξ t 0 )V t 0 .
From (4.2), it follows that

∂ 2 x ψ k (t, x t k ) = V t k (U t k ) -1 where U t k = D y x t k + iD η x t k and V t k = D y ξ t k + iD η ξ t k ,
and a similar relation holds true for ∂ 2 x ψ -1 (t, x t -1 ). The reflected amplitudes evaluated on the associated rays satisfy transport equations which are similar to (2.20) and may be written as

d dt a k 0 ( ′ ) (t, x t k ) + 1 2 Tr H 21 (x t k , ξ t k )U t k + H 22 (x t k , ξ t k )V t k (U t k ) -1 a k 0 ( ′ ) (t, x t k ) = 0.
One can obtain a similar equation to (

a k 0 ( ′ ) (t, x t k ) = a k 0 ( ′ ) (T k , x T k 0 ) det U t k det U T k k -1 2 , k = ±1,
where the square root is obtained by continuity from 1 at t = T k . On the other hand, for k = ±1

d 0 -mB + d k -mB = b(x, ∂ x ψ 0 )a 0 0 + b(x, ∂ x ψ k )a k 0
, where b denotes the principal symbol of B. Thus, the condition 2 p.8 required for the construction of the reflected amplitudes implies that a k 0

( ′ ) (T k , x T k 0 ) = sa 0 0 ( ′ ) (T k , x T k 0 )
, with s = -1 for Dirichlet condition and s = 1 for Neumann condition.

In order to find the relationship between U T k k and U T k 0 for k = ±1, we differentiate the equality

x T k k = x T k 0 D y,η x T k k + ẋT k k (∂ y,η T k ) T = D y,η x T k 0 + ẋT k 0 (∂ y,η T k ) T , and compute the derivatives of T k from the condition x T k 0 ∈ ∂Ω ∂ y,η T k = - 1 ẋ0 T k • ν(x T k 0 ) (D y,η x T k 0 ) T ν(x T k 0 )
to get after elementary computations

U T k k = Id -2ν(x T k 0 )ν(x T k 0 ) T U T k 0 . (4.9)
Hence

a k 0 (t, x t k ) = -si det U t k -1 2 and a k 0 ′ (t, x t k ) = s(c(y)|η|) -1 det U t k -1 2 for k = ±1,
where the square root is defined by continuity from i[det U T k 0 ] -1 2 at t = T k .

Approximation operators

We briefly recall a simple version of the integral operators with complex phases used in [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF] and the estimates established therein. We then use these results to prove Lemma 2.2. For t ∈ [0, T ], let K z,θ (t) be a compact of R 2n and consider the set h(z, θ)l ε (t, x, z, θ)(xz) α e iΦ(t,x,z,θ)/ε dzdθ, h ∈ L 2 (R 2n ), then, under the previous hypotheses on Φ and l ε , we have the following estimate: Hence, only (r k 0,0 ) k contributes to ∂ t,x u appr ε,r0,r∞ , the residue being of order √ ε. One has

E 1 = {(t, x, z, θ) ∈ [0, T ] × R 3n , (z, θ) ∈ K z,θ (
r
If u ∈ L 2 (R n ) is supported in K then T * ε θ(η)T ε u L 2 (E) e -C/ε u L 2 (R n ) .
Proof: The kernel of

1 E T * ε θ(η)T ε 1 K : L 2 (R n x ) → L 2 (R n w ) is k ε (w, x) = ε -3n 2 c 2 n 1 E (w)1 K (x) R 2n
θ(η)e iη•(w-x)/ε-(y-x) 2 /(2ε)-(w-y) 2 /(2ε) dydη = 1 E (w)1 K (x)ε -n (2π) -n F θ( xw ε )e -(x-w) 2 /(4ε) .

For w ∈ R n , one has by Cauchy-Schwartz inequality

R n |k ε (w, x)|dx ≤ F θ L 2 (R n ) (2π) -n ε -n 2 R n
1 E (w)1 K (x)e -(x-w) 2 /(2ε) dx

T ε T * ε θ(η)T ε L 2 (R n )→L 2 (R n ×E) e -C/ε .
Proof: Consider the operator H ε : L 2 (R 2n y,η ) → L 2 (R 2n x,ξ ) defined by

H ε f (x, ξ) = 1 E (ξ)T ε T * ε (1 K (η)f (y, η)) (x, ξ).
It is easy to compute its kernel h ε h ε (x, ξ, y, η) = c 2 n π n 2 ε -n 1 E (ξ)1 K (η)e i(ξ+η)•(x-y)/(2ε)-(x-y) 2 /(4ε)-(ξ-η) 2 /(4ε) .

Hence, R 2n |h ε (x, ξ, y, η)|dxdξ e -C/ε and R 2n |h ε (x, ξ, y, η)|dydη e -C/ε . For u ∈ L 2 (R n ) , it follows by Schur's Lemma that

H ε T ε u L 2 (R 2n x,ξ ) = T ε T * ε θ(η)T ε u L 2 (R n x ×E ξ ) e -C/ε T ε u L 2 (R 2n y,η )
e -C/ε u L 2 (R n ) . Proof: Differentiating (2.10) w.r.t. y j , 0 ≤ j ≤ n, yields

ε 1 2 ∂ yj (T ε u I ε,r0,r∞ ) = iη j ε -1 2 T ε u I ε,r0,r∞ -c n ε -3n 4 R n u I ε,r0,r∞ (x)ε -1 2 (y j -x j )
e iη•(y-x)/ε-(y-x) 2 /(2ε) dx.

The l.h.s. is bounded in L 2 y,η because ∂ yj (T ε u I ε,r0,r∞ ) = T ε (∂ xj u I ε,r0,r∞ ). The second term of the r.h.s. is the Fourier transform of a bounded function in L 2

x , thus it can be estimated using Parseval equality. One gets

ε -3n 4 R n u I
ε,r0,r∞ (x)ε -1 2 (y jx j ) e iη•(y-x)/ε-(y-x) 2 /(2ε) dx L 2 y,η u I ε,r0,r∞ L 2

x .

Thus ε - Proof: Let φ be a test function in C ∞ 0 (R n × (R n \{0}), R) and denote c ε = |D|b ε . We use another expression of < w ε (a ε , c ε ), φ > exhibiting the Fourier transform of c ε :

< w ε (a ε , c ε ), φ >= (2π) -n R 2n F ξ φ(x - ε 2 v, v)a ε (x)c ε (x -εv)dvdx.
Since F ξ φ is rapidly decreasing

sup x F ξ φ(x - ε 2 v, v) -F ξ φ(x, v) ε(1 + v 2 ) -n-1 .
By Cauchy-Schwartz inequality w.r.t. dx

R 2n | F ξ φ(x - ε 2 v, v) -F ξ φ(x, v) a ε (x)c ε (x -εv)|dvdx ε a ε L 2 c ε L 2 .
It follows that

< w ε (a ε , c ε ), φ >= (2π) -n R 3n
φ(x, ξ)e -iv•ξ a ε (x)c ε (xεv)dvdxdξ + o(1).

  1, . . . , n) will satisfy the usual assumptions needed in the general context of the study of Wigner measures: their Wigner measures are supposed unique and v I ε and ∂ xj u I ε , j = 1, . . . , n, are ε-oscillatory (see (3.23) and (3.24)), (H2) the Wigner measures of v I ε and ∂ xj u I ε , j = 1, . . . , n, do not charge the set R n × {ξ = 0}.

Lemma 3 . 3

 33 Let (f ε ), (g ε ) and Φ, Ψ satisfy the hypotheses of Lemma 3.1. If α and β are in C 1 (F, C) then

F

  Note that d ε e iγεd 0 e iγ0 is dominated by |d εd 0 | + |d 0 | e iγε-iγ0 -1 . The convergence of d ε when ε → 0 to its limit d 0 is uniform w.r.t. (x, η) and so is the convergence of γ ε to γ 0 on the support of d 0 . Thus d ε e iγε converges to d 0 e iγ0 uniformly w.r.t. (x, η). It follows that sup x η d ε e iγεd 0 e iγ0 (x, y ′ , u, r ′ , δ ′ , u) → ε→0 0 for every y ′ , u, r ′ , δ ′ .

  The term involving the derivative of ρ is exponentially decreasing by Lemma 4.4. Since the FBI transform of a derivative is the derivative of the FBI transform by (2.13), one has to estimate ∂ xj u I ε -ρT * ε γ r0,r∞ T ε ∂ xj u I ε L 2 (Ω) . Employing the same previous techniques yields for j = 1, . . . , n lim sup ε→0 ∂ xj u I ε -∂ xj u I ε,r0,r∞ L 2 (Ω) 4r 0 ≤ l j (α) and w ∂ xj u I ε ({|ξ| ≤ l j (α)}) ≤ α. Set r 0 = 1 4 min 0≤j≤n l j (α), then the Poincaré inequality yields the same bound for lim sup ε→0 u I εu I ε,r0,r∞ L 2 (Ω) . Coming back to (3.21) we deduce that

  .23) lim sup ε→0 ε|ξ|≥R |F ∂ xj u I ε (ξ)| 2 dξ → R→+∞ 0 for j = 1, . . . , n. (3.24)

  v)) and a similar relation holds true for ∂ v (ψ inc ) b , one gets from (4.4) and (4.8) that ∂ v θ b k ≍ (t,v)=(T1,v1) ∂ v (ψ inc ) b . Hence θ b and (ψ inc ) b

2 4. 1 . 3

 213 First order reflected beams' phases and amplitudes Lemma 4.1 gives at order one

  t), |x -z| ≤ 1}, which we assume compact. Let Φ be a phase function smooth on an open set containing E 1 and satisfying (2.23) for t ∈ [0, T ] and (z, θ) ∈ K z,θ (t). Then there exists r[Φ] ∈]0, 1] s.t. Im Φ(t, x, z, θ) ≥ C(xz) 2 for t ∈ [0, T ], (z, θ) ∈ K z,θ (t) and |x -z| ≤ r[Φ]. Let l ε ∈ C ∞ ([0, T ] × R 3n , C) satisfying for t ∈ [0, T ], l ε (t, x, z, θ) = 0 if (z, θ) / ∈ K z,θ (t) or |x -z| > r[Φ], ε k 2 ∂ k xj l ε is uniformly bounded in L ∞ ([0, T ] × R 3n )for every 1 ≤ j ≤ n and k ∈ N. (4.10) If O α (l ε (t, .), Φ(t, .)/ε) denotes, for a given multiindex α and t ∈ [0, T ], the operator [O α (l ε (t, .), Φ(t, .)/ε) h] (x) = R 2n

= 1 B

 1 6) the properties (2.23) for t ∈ [0, T ] and (z, θ) ∈ K k z,θ (t). We fix some r[ ψ k k ] ∈]0, 1] so that Im ψ k k (t, x, z, θ) ≥ C(xz) 2 for t ∈ [0, T ], (z, θ) ∈ K k z,θ (t) and |x -z| ≤ r[ ψ k k ]. For t ∈ [0, T ], Π k ρ ′ ⊗ γ ′ k (t, z, θ) r k j,α k (t,x, z, θ) depends smoothly on its variables and vanishes for |x -z| > d or (z, θ) / ∈ K k z,θ (t). Hence, upon choosing d ≤ r[ ψ k k ], the amplitude Π k ρ ′ ⊗ γ ′ k r k j,α ksatisfies the properties formulated in (4.10). Let uscheck if 1 B f ε k T ε v I ε,r0,r∞ k , 1 B ε -1 T ε u I ε,r0,r∞ k is uniformly bounded in L 2 (R 2n ). Clearly T ε v I ε,r0,r∞ is, and the property holds true for ε -1 T ε u I ε,r0,r∞ by Lemma 4.6. One can then use the approximation operators O α to write the integral (2.22) asε -z) α e i ψ k k /ε dzdθ =ε -3n 4 +j O α Π k ρ ′ ⊗ γ ′ k z) α e i ψ k k /ε dzdθ L 2

Lemma 4 . 2 -n 4 R n F 2 Lemma 4 . 3 ([ 4 ], Lemma 2 . 4 ) 2 Lemma 4 . 4

 424n243424244 For u in L 2 (R n ) T ε u(y, η) = c n (2π) -n 2 ε u(ξ)e -(η-εξ) 2 /(2ε) e iξ•y dξ.Proof:The equality is proven by Parseval formula. Let a be a positive real, E a measurable subset ofR n and K ⊂ R n a compact set s.t. dist(K, E) ≥ a. If u ∈ L 2 (R n x ) is supported in K then T ε u L 2 (E×R n η ) = c n ε -n 4 1 E (y)u(x)e -(x-y) 2 /(2ε) L 2 y,xe -a 2 /(4ε) u L 2x . Proof: The proof consists of writing the FBI transform as the Fourier transform w.r.t. x of some auxiliary function and using Parseval equality. Let θ be a cut-off of C ∞ 0 (R n η , R), E a measurable subset of R n and K ⊂ R n a compact set s.t. dist(K, E) > 0.

Lemma 4 . 6 (

 46 [START_REF] Bougacha | Gaussian beams summation for the wave equation in a convex domain[END_REF], Lemma 3.4)ε -1 T ε u I ε,r0,r∞ L 2 (R 2n ) 1.

2 Lemma 4 . 7

 247 Reproducing the same arguments on the equality∂ yj (T ε u I ε,r0,r∞ ) = iη j ε -1 T ε u I ε,r0,r∞c n ε -3n 4 R n ε -1 2 u I ε,r0,r∞ (x)ε -1 2 (y jx j ) e iη•(y-x)/ε-(y-x) 2 /(2ε) dx leads to u I ε,r0,r∞ L 2 (R n ) ε. Let (a ε ) and (b ε ) be two sequences uniformly bounded in L 2 (R n ) and H 1 (R n ) respectively. If ε -1 b ε is uniformly bounded in L 2 (R n ), then w ε (a ε , |D|b ε ) ≈ |ξ|w ε (a ε , ε -1 b ε ) on R n × (R n \{0}).

  Integrating w.r.t. v leads to< w ε (a ε , c ε ), φ >= (2π) -n ε -n R 2n φ(x, ξ)e -ix•ξ/ε a ε (x)F c ε (ξ/ε)dxdξ + o(1),and replacing F c ε (ξ/ε) by ε -1 |ξ|F b ε (ξ/ε) ends the proof.2

  .2) ϕ t 0 and ϕ t 1 are symplectic C ∞ diffeomorphisms from B to B [22], and so is s 1 . One can define a similar auxiliary function s -1 : B → B s.t. ϕ t -1 = ϕ t 0 o s -1 and ϕ t -1 = s -1 o ϕ t 0 for t ∈ R. Let us introduce the components of s 1 as

  For m ≤ R, the l.h.s. is zero so one can show recursively on |β| ≤ R that ∂ β x π(t, x t 1 ) = 0. One thus has the following eikonal equation on θ

	t 0 ) = 0 and 1 ) is zero. More generally, 0 ) is non singular, it follows by (4.6) that ∂ x π(t, x t 0 , ξ t A(t, x t
	for m ≥ 1, the formula of composite functions' high derivatives yields
	n	n
	∂ xi 1 . . . ∂ xi m [π (t, ̺(t, x))] (t, x t 0 ) =	∂ xj 1 . . . ∂ xj m π(t, x t 1 )	A j k i k (t, x t 0 , ξ t 0 )
	j1,...,jm=1	k=1
	+ z i1...im (t),

.7)

Putting together (4.3), (4.4) and (4.7) shows that the phase θ satisfies

p (̺, ∂ t θ(t, ̺), ∂ x θ(t, ̺)) R ≍ x=x t 0 0 under the further assumption k ≥ R. Let π(t, x) = p (x, ∂ t θ(t, x), ∂ x θ(t, x)). Since ∂ x (π(t, ̺)) (t, x

where z i1...im depends on derivatives of π on (t, x t 1 ) of order lower than m.

  Proposition 3 ([4], Lemma 3.3) O α (l ε (t, .), Φ(t, .)/ε) L 2 (R 2n )→L 2 (R n ) εThis estimate allows to prove Lemma 2.2.Proof:[Proof of Lemma 2.2] Consider the integrals (2.22) giving the derivatives of u appr ε,r0,r∞ and fix j, k and α. The transported phase ψ k k is smooth and satisfies by (2.3), (2.4) and (2.

	3n 4 +	|α| 2	uniformly w.r.t. t ∈ [0, T ].

  ∂ t ψ k (t, x t k ) = -c(x t k )|ξ t k |, ∂ x ψ k (t, x t k ) = ξ t k . (-i)β k c(z)|θ|χ d (xz)Π k (t, z, θ) ρ ′ ⊗ γ ′ iβ k θ j χ d (xz)Π k (t, z, θ) ρ ′ ⊗ γ ′ iβ k θ j χ d (xz)Π k (-t, z, θ) ρ ′ ⊗ γ ′ (R n ),uniformly w.r.t. t ∈ [0, T ]. One can get rid of the cut-off χ d (xz) appearing in ∂ t u appr ε,r0,r∞ (t, x) and ∂ xj u appr ε,r0,r∞ (t, x) by using the estimate (2.24). 2 4.3 Results related to the FBI and the Wigner transforms

	in L 2 (R n ), uniformly for t ∈ [0, T ] and
			∂ xj u appr ε,r0,r∞ (t, x)
		=	1 2	ε -3n 4 c n	R 2n k=0,1	k	(t, z, θ)
							p ε,k	k (t, z, θ)e i ψ k	k (t,x,z,θ)/ε dzdθ
				+	1 2	ε -3n 4 c n	R 2n k=0,-1	k	(-t, z, θ)
				+ O( √	ε)	q ε,k	k (-t, z, θ)e i ψ k	k (-t,x,z,θ)/ε dzdθ
	in L 2
							( ′ ) k (t, y, η),
	and by (2.3)
	It follows that
		∂ t u appr ε,r0,r∞ (t, x)
	=	1 2	ε -3n 4 c n	R 2n k=0,1	k	(t, z, θ)
							p ε,k	k (t, z, θ)e i ψ k	k (t,x,z,θ)/ε dzdθ
		+	1 2	ε -3n 4 c n	R 2n k=0,-1
		+ O( √	ε)

k 0,0 (t, x, y, η) = i 2 c n β k ∂ t,x ψ k (t, x t k )χ d (xx t k )a iβ k c(z)|θ|χ d (xz)Π k (-t, z, θ) ρ ′ ⊗ γ ′ k (-t, z, θ) q ε,k k (-t, z, θ)e i ψ k k (-t,x,z,θ)/ε )dzdθ

  Similarly, R n |k ε (w, x)|dw is dominated by e -C/ε , so one gets by Schur's LemmaT * ε θ(η)T ε u L 2 (Ew) e -C/ε u L 2 (R n x ) . Lemma 4.5 Let E be a measurable subset of R n and K ⊂ R n a compact set s.t. dist(K, E) > 0. If θ is a cut-off of C ∞ 0 (R n η , R) supported in K then

	2
	1
	2

e -C/ε .