A Gaussian beam approach for computing Wigner measures in convex domains

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, with initial conditions compactly supported in this set, 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 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 [50] 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 [29,33,34,35] 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 [20], see also [17,16]. Wigner measures are also related to the H-measures and microlocal defect measures introduced in [49] and [18], see also [6,1]. 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 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. [20,40]) are related to the energy density in the high frequency limit. More precisely, under suitable hypotheses (see Proposition 1.7 in [20]), 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 Moreover, E u C ε (t, .) is the sum of two measures satisfying transport equations of Liouville type (see e.g. [20]).
v I ε = 1 Ω v I ε and ∂ xj u I ε = 1 Ω ∂ xj u I ε (j = 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)), the Wigner measures of v I ε and ∂ xj u I ε , j = 1, . . . , n, do not charge the set R n × {ξ = 0}. (H3) 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 ε, Ω 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 [37] who proved refraction results for sharp interfaces and Burq [5] who described their support for a Dirichlet boundary condition. Similar results have been established for other problems [8,13,15], in particular the eigenfunctions for the Dirichlet problem [52,19] and for the Neumann and Robin problems [7]. 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 [3,44]. 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 [44], to prove lack of observability [32] and to study semiclassical measures [41] and trace formulas [51,12].
To describe non localized solutions of PDEs, one can use the Gaussian beam summation method [24,10,25]. 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 [38,47,2] or continuous [30,31], the selection of the beams to be superposed can be done according to several criteria. In [4], 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 [27,28,43] 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 [45] for the Schrödinger equation, and more recently the paper by Castella [9] 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 ∂Ω.
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 [45] 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 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.
2 Tool-box and construction of the asymptotic solution We recall the construction made in [4] 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.

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 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 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 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 [44,23,32] is the principal symbol of P and h.o.t. denotes higher order terms. The first equation is then the eikonal equation on x = x t up to order 2 (see Remark 2.1 in [4] for an explanation of the choice of this specific order), which means Orders 0 and 1 of the previous equation are fulfilled on the ray by setting Choosing ψ(0, x 0 ) as a real quantity, it follows that 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 [23]).
The phase is defined beyond the ray as a polynomial of order 2 w.
Next, we make the term associated to the power ε −1 in the expansion (2.1) vanish which leads to a linear ordinary differential equation (ODE) on a 0 (t, x t ). The amplitude is then chosen under the form where d is a positive parameter. The constructed beams are thus defined for all (t, x) ∈ R n+1 and they satisfy the estimate 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 is an open set of R n , the Hamiltonian flow ϕ t 0 (y, η) = (x t 0 (y, η), ξ t 0 (y, η)) satisfying: 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 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.
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.
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, These constraints uniquely determine the reflected phases and amplitudes, once the incident ones are fixed [44]. 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 [44] B ε − n 4 +1 ω 0 ε (., y, η) 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 [4]. 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 [4], 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 [36]) is, for a given scale ε, the operator T ε : (2.10) with adjoint operator given by As the Fourier transform, the FBI transform is an isometry, satisfying T * ε T ε = Id. The extra assumption on the initial data needed in [4] is 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.
(2.11) Lemma 4.5 from the Appendix (Section 4) yields 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 are fulfilled since Lemma 4.4 from the Appendix implies that (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 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 Let ρ ′ be a cut-off of C ∞ 0 (R n , [0, 1]) supported in a compact K y ⊂ Ω and satisfying ρ ′ (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 [4] 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 [4], and it reads for first order beams as follows: x, y, η) dydη Above, ω 0 ε , ω 0 ε ′ are incident Gaussian beams with the same phase ψ 0 satisfying at and different amplitudes a 0 0 , a 0 0 ′ satisfying 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 and sup We refer to [4] 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 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 leads to the following ordinary differential system on Using the symmetry of the following matrices 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 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 [23] 20) which may be written using the matrices U t 0 and V t 0 as 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 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. [39,47]). One gets (see Appendix 4.1) The reflected amplitudes evaluated on the rays have an expression similar to the incident amplitudes (see Appendix 4.1) where the square root is defined by continuity from i det for the Dirichlet boundary condition and s = 1 for the Neumann condition. We summarize the previous form of the beams in the following result:

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) 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 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 The volume preserving change of variables transforms the previous integrals as (2.22) We can write the leading terms obtained for j = 0 and α = 0 using Gaussian type integrals I ε (h, Φ) defined as for a given phase function Φ ∈ C ∞ (R n+1 t,x × B, C) polynomial of order 2 in x − z and satisfying, for t ∈ [0, T ] and (z, θ) ∈ B 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 Noticing that e iΦ/ε is exponentially decreasing for |x − z| ≥ 1, one can use the following crude estimate 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. .

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 As regards the Wigner transforms associated to cv + x,ε (t, .) , since c is uniformly continuous on R n , one has by a classical result ( [20], p.8) and therefore the involved quantities have the form 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 with κ ε , τ ε = ε −1 u I ε,r0,r∞ , v I ε,r0,r∞ and k, l = 0, ±1, or after expanding the FBI transforms 3) This type of oscillating integrals is traditionally estimated by the stationary phase theorem. For example, this method was successfully used in [9] 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 and then study the whole integral involving κ ε (w)τ ε (w ′ ).
An alternative method was used in [45], 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 [45] 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 [45] lead to the following more general result: 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 for x ∈ R n and (z, θ), (z ′ , θ ′ ) ∈ F , with r Φ , r Ψ ∈ C ∞ (F, R) and the matrices 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 The Fourier transform of a Gaussian functions product is given by the following Lemma, whose proof is postponed to the end of this Section:  Hence Making the changes of variables , Step 2. Gaussian approximations. Taking the duality product of the Wigner transform with a test function φ ∈ C ∞ 0 (F, R), and after setting ( 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 The r.h.s. of (3.4) may be written as Leibnitz formula yields for a multiindex α 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.
The second integral in the r.h.s. of (3.5) is then dominated by We deduce by Cauchy-Schwartz inequality w.r.t. s, σ that 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 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 Using the value of the Fourier transform of a Gaussian function (see Theorem 7.6.1 of [21]), it follows that 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 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.

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 which leads to 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 by using the standard estimate (see Proposition 1.1 in [20]) 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, .) ≈ Thus, we are left with the computation of the Wigner measure associated to (v + t ), computations being similar for (v − t ). One has 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: , so µ t ε,k may be written as (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 Above ϕ t k is extended outside B as the identity. Proof: We simplify the integral 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- 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 where Y t k and Z t k are the 2n × 2n matrices 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 Combining this relation with the symplecticity of F t k , one gets the following relation Moving to the amplitude A(Φ k , Φ k ) = c 2 n 2 Hence 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 Analysis of the transported FBI transforms. It remains to analyze the most difficult terms in the amplitude, which involve transported FBI transforms [21]). We use Taylor's formula for this map to get for (s ± √ εr, σ ± √ εδ) ∈ K k z,θ (t) and (s, σ) ∈ suppφ is thus appropriate. Notice that for (s, σ) ∈ o T * Ω one has the following relations [26] 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 We perform the changes of variables Notice that d ε (x, y ′ , η, r ′ , δ ′ ) converges when ε → 0 to η). On the other hand, since εr ± x are the remainder terms in the Taylor expansions of One has Cauchy-Schwartz inequality w.r.t. dx insures that the bracket integral is less than κ ε L 2 τ ε L 2 . Let us examine the term For fixed y ′ , r ′ , δ ′ the functions d ε and d 0 are compactly supported w.
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 On the other hand, successive integrations by parts give with L a differential operator w.r.t. η, of order 2n. Thus, for all (x, y ′ , η, r ′ , δ ′ ) ∈ suppd ε and u ∈ R n . Thus, there exists C, C ′ > 0 s.t.
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 Recalling the relation between the Wigner measure and the FBI transform (see 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 By summing over k = 0, 1 and letting ε → 0, we get 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: 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 ] It follows that 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 Using these results in (3.11) as ε → 0 leads to E u appr ε,r0,r∞ (t, .) = (3.18)

Proof of the main Theorem
A consequence of the estimate (3.10) is 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) 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 [14] and [42] for WKB initial data.

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 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 used for |α| = r to eliminate the r + 1-th order derivatives of ψ in equation (4.1).
To summarize, the requirements 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 [44] 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 [39]. 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 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 [11] 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 ϕ Besides, noticing that T 1 (ϕ t 0 ) = T 1 − t, one has also ϕ 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 : Let us introduce the components of s 1 as 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 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 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 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, τ, ξ).
which implies, by construction of ψ inc Compare this with the equation resulting from the construction of ψ ref and (4.2). This suggests the following Lemma 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) Hence, θ exists if A(t, x t 0 , ξ t 0 ) is non singular and From (4.2) one gets Since ϕ t 1 is symplectic, the matrix is symplectic. This implies in particular the relation and at the same time, 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, Since M T JM = Ds T 1 JDs 1 , the symplecticity of s 1 leads to Hence 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 on the maps (t, x, ξ) → ∂ x θ (t, r(x, ξ)), (x, ξ) → λ(x, ξ) and their derivatives yields recursively by (4.4) Using the relations ∂ 2 x θ(t, ̺) x, ∂ x ψ inc ) and (4.5) in the previous equation yields so one gets Putting together (4.3), (4.4) and (4.7) shows that the phase θ satisfies under the further assumption k ≥ R. Let π(t, x) = p (x, ∂ t θ(t, x), ∂ x θ(t, x)). Since ∂ x (π(t, ̺)) (t, x t 0 ) = 0 and A(t, x t 0 , ξ t 0 ) is non singular, it follows by (4.6) that ∂ x π(t, x t 1 ) is zero. More generally, for m ≥ 1, the formula of composite functions' high derivatives yields where z i1...im depends on derivatives of π on (t, x t 1 ) of order lower than m. 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 θ To compare the time and tangential derivatives of θ and ψ inc at ( where N is an open subset of R n−1 , σ(N ) = U and σ is a diffeomorphism from N to U. For x ∈ R n close to x T1 0 , we may write x = σ(v) + v n ν (σ(v)) , withv ∈ 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 which implies by (4.7) that Similarly λ(X, Ξ) = Ξ − 2 (Ξ · ν(X)) ν(X) for (X, Ξ) ∈ o T * R n | ∂Ω , leading to ) and a similar relation holds true for ∂v (ψ inc ) b , one gets from (4.4) and (4. have the same time and tangential derivatives at (T 1 ,v 1 ) from the order 1 to the order k + 1.
If we assume that θ(T 1 , 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.
One obtains by plugging the expression (2.19) of ∂ 2 From (4.2), it follows that 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 One can obtain a similar equation to (2.16) on U t k involving H 21 (x t k , ξ t k ) and H 22 (x t k , ξ t k ), by using the relation ϕ t k = ϕ t 0 o s k . On the whole 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 to get after elementary computations 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 [4] 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 Im Φ(t, x, z, θ) ≥ C(x − z) 2 for t ∈ [0, T ], (z, θ) ∈ K z,θ (t) and |x − z| ≤ r[Φ].
2 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 T ε T * ε θ(η)T ε L 2 (R n )→L 2 (R n ×E) e −C/ε .