Lukas Rummelhard

Amaury Nègre

Mathias Perrollaz

Christian Laugier

Probabilistic Grid-based Collision Risk Prediction for Driving Application

Keywords: Dynamic Occupancy Grid, Risk assessment, Time To Collision, Intelligent Vehicle

In the recent years, more and more modern cars have been equipped with perception capabilities. One of the key applications of such perception systems is the estimation of a risk of collision. This is necessary for both Advanced Driver Assistance Systems and Autonomous Navigation. Most approach for risk estimation propose to detect and track the dynamic objects in the scene. Then the risk is estimated as a Time To Collision (TTC) by projecting the object's trajectory in the future. In this paper, we propose a new grid-based approach for collision risk prediction, based on the Hybrid-Sampling Bayesian Occupancy Filter framework. The idea is to compute an estimation of the TTC for each cell of the grid, instead of reasoning on objects. This strategy avoids to solve the difficult problem of multi-objects detection and tracking and provides a probabilistic estimation of the risk associated to each TTC value. After promising initial results, we propose in this paper to evaluate the relevance of the method for real on-road applications, by using a real-time implementation of our method in an experimental vehicle.

Introduction

In the development of ADAS (Advanced Driver Assistance Systems) that can prevent collision to append, a basic requirement is an assessment of the collision risk. This paper presents an original method for collision risk estimation. The idea is to compute an estimation of the TTC for each cell of a grid, instead of reasoning on objects. This strategy avoids to solve the difficult problem of multi-objects detection and tracking and provides a probabilistic estimation of the risk associated to each TTC value. The method has been implemented for real time operation on GPU in our experimental vehicle. Then the vehicle is used for real on-road testing in order to evaluate the method in real conditions. The paper is organized as follow: Section II presents related work. Section III describes the method. Section IV provides experimental results. Section V concludes.

Related work

In the literature, various approaches have been proposed in order to detect potential collisions in advance and trigger an assistance like a driver alarm or an automatic braking. Most of them rely on the detection and tracking of the moving objects in the scene. This tracking is used to estimate a risk metric, used for a decision layer. The simplest approach consists in computing a simple Time to Collision (TTC) by predicting the trajectories of both the ego-vehicle and the other objects using a simple linear motion model [START_REF] Labayrade | Experimental assessment of the rescue collision-mitigation system[END_REF]. More advanced approach propose to generate a set of possible trajectories for all the objects, including the ego-vehicle, and compute the number of possible intersections between those trajectories [START_REF] Kaempchen | Situation assessment of an autonomous emergency brake for arbitrary vehicle-to-vehicle collision scenarios[END_REF]. The trajectory generation is performed by applying a set of possible controls to the objects. This approach has the advantage of providing a probabilistic estimation of the risk.

Considering that a control can change over time and is dependant from the intention of the driver, the authors in [START_REF] Laugier | Probabilistic Analysis of Dynamic Scenes and Collision Risk Assessment to Improve Driving Safety[END_REF] propose to constrain the possible trajectories using the estimated maneuver intention of the driver. This allows to consider longer time prediction.

All this methods rely on the detection and tracking of the moving objects [START_REF] Petrovskaya | Awareness of road scene participants for autonomous driving[END_REF]. This stage is a difficult problem to solve, which can be computationally costly and can generate errors in the process. Alternatively to methods based on objects trajectory prediction, other methods propose to estimate time to collision based on visual features [START_REF] Alenya | A Comparison of Three Methods for Measure of Time to Contact[END_REF], but such methods are only designed to work with computer vision.

Another alternative to the object-based approaches are the grid-based approaches, like the Bayesian Occupancy Filter [START_REF] Tay | An efficient formulation of the bayesian occupation filter for target tracking in dynamic environments[END_REF], or its extension Hybrid -Sampling Bayesian Occupancy Filter (HSBOF) [START_REF] Negre | Hybrid sampling bayesian occupancy filter[END_REF], which represents the environment as a probabilistic occupancy grid. This framework allows to model both the static and the dynamic environment, by estimating velocity probability distributions for each cell in the grid. It is efficient for short term prediction, but since there are no notions of objects or behaviors, longer term prediction requires to integrate additional prior knowledge like map information [START_REF] Gindele | Bayesian occupancy grid filter for dynamic environments using prior map knowledge[END_REF].

Grid based methods are designed to be very efficiently implemented on parallel architecture and avoid to deal with object-level data association for both tracking and sensor fusion. Consequently, they are able to provide in real time a robust representation of the environment. However, in order to use such approaches for driving applications, the authors generally add a subsequent clustering stage to retrieve the notion of objects before any application [START_REF] Danescu | Modeling and tracking the driving environment with a particle-based occupancy grid. Intelligent Transportation Systems[END_REF][START_REF] Laugier | Probabilistic Analysis of Dynamic Scenes and Collision Risk Assessment to Improve Driving Safety[END_REF].

An hybrid approach is presented in [START_REF] Laugier | Dynamic obstacle avoidance in uncertain environment combining pvos and occupancy grids[END_REF], with velocity obstacles. The authors combine a notion of objects with an estimation of velocities within a dynamic occupancy grid.

Based on the occupancy grid framework, we propose in this paper a method that allows to estimate both a probabilistic collision risk and a TTC for each cell of a grid, without any detection and tracking of objects. This dense representation of the risk is intended to be used to take decisions about the short term evolution of the scene. In this work, we present the approach and propose to evaluate how it can be used for automotive applications.

3 Description of the method

HSBOF presentation

The Hybrid Sampling Bayesian Occupancy Filter (HSBOF) [START_REF] Negre | Hybrid sampling bayesian occupancy filter[END_REF] is a Bayesian filtering perception technique which models the environment at a sub-object level, in term of spatial occupancy and dynamics. The surrounding of the subject is divided into cells, to which are associated random variables, symbolizing their occupancies and velocities. The scene is interpreted through the estimated distributions of those variables, estimations that are recursively updated according to the observations. The motion field of the scene is inferred through the spatial occupancy evolution, and is described as the combination of a static part, depicted by likelihood values in a regular grid, and a dynamic part, sampled by moving particles which transfer occupancy between cells. All the distributions are jointly generated and updated. The distinction in dynamics allows to apply specific motion representations, and then to efficiently allocate memory and computation power. Indeed, while the static part of the scene does not require an important number of samples to symbolize its velocity distribution, the dynamic one necessitate a substantial precision and reactivity. The figure 1 summarizes the HSBOF scene representation. This data structure also allows to conveniently change the reference frame, which is needed in case of a moving subject. Using motion sensors embedded in the mobile robot (Inertial Measurement Unit, GPS, Wheel speed and steering sensor, visual odometry, etc.), the displacement of the grid between two updates is estimated, then applied on the models. The global framework of the HSBOF is summarized on the figure 2.

Grid-based collision assessment

As an output of the HSBOF, are available estimates of the robot surrounding occupancy and dynamics. Thanks to the already mentioned embedded motion sensors, the ego-motion is estimated. Using those and a robot model, which specifies its form, size and motion model, a time projection is achieved to assess collision risk in the future and localize it in the grid. Indeed, to each current occupancy component is associated a velocity and motion model, in such a way that its position over time can be predicted. Simultaneously, the space occupied by the robot over time can also be assessed. The risk evaluation consists in a succession of scene configuration and robot position predictions, and so potential overlaps, in other words collisions can be assessed in time and space (figure 3).

In practice, to each static cell and dynamic particle of the HSBOF representation are associated a "Time To Contact" (TTC) value and a risk weight. Until a time horizon, both vehicle and elementary components are iteratively moved according to their motion models. The time step used is chosen according to the grid resolution and the maximal reachable speed, to prevent any failed collision detection due to time resolution. The intersection between the vehicle and each basic component is evaluated. If a collision is confirmed, the TTC value is set and a weight is associated to the element, according to its associated quantity of occupancy. Otherwise, the next time step is processed. Beyond the chosen time range, if no collision is observed, a null weight is assigned (since there is no predicted collision, the TTC value is meaningless, by convention it is set to the time limit). Those computations performed, for each cell in the current occupancy grid a risk value over the time range is calculated, corresponding to the likelihood of the cell to be occupied, and for its content to collide with the robot during the selected period. The risk calculation simply consists in a sum over this time range of the related weights.

As a result at this stage, a risk grid presents for each cell the probability to hit an obstacle coming from that cell over a given period if the scene dynamics remain unchanged (conservative model, relevant for short-term prediction). Given that grid, a single global collision value is computed, and used to inform the system of any imminent danger. Two approaches are proposed to evaluate this general risk :

The first one proposes to evaluate the likelihood of colliding with a selected surface size. According to this area, and to the area of a cell, the number of n cells to be struck is thus given, the computed risk value being the probability to collide with n cells of the grid. This evaluation presents the advantage of rightfully integrate risk when the occupancy estimation is diffuse. However, in the case of high speed robot displacement, the integration of occluded spaces, yet with unknown occupancy, can lead to systematic high risk scores. Although crossing unknown areas is rightfully dangerous, the estimation becomes meaningless under such circumstances.

The second evaluation selects, after applying a median filter on the risk grid to remove noisy data, the maximum value of collision over the cells. This simpler approach turns out to be more effective in practice.

By comparing this value to a threshold, a risk detector is thus defined. The system is then used to generate alerts, for the driver if the timing makes it relevant, and for a driving assistance device.

t t+1 t+2 t t+1 t+2

Collision Fig. 3. Collision risk estimation over time for a specific cell. The cell position is predicted according to its velocity, along with the mobile robot. This risk profile is computed for every cell, and then used to integrate over time the global collision risk.

Experimental Results

Experimental platform

For the experiments, we equipped a Lexus LS600h car with two Ibeo Lux multilayer lidars under the two front lights (see Fig. 4). The horizontal field of view covered by the two lidars is almost 160 degrees. For the odometry calculation we use vehicle velocity and steering data collected from the CAN bus system. The input of our algorithm is a sequence of occupancy grids computed from the two lidars. To merge the eight laser scan layers acquired by the two lidars and to compute the input observation occupancy grid, we use a method similar to [START_REF] David Adarve | Computing Occupancy Grids from Multiple Sensors using Linear Opinion Pools[END_REF].

The vehicle is used to collect real on-road data for both downtown and highway scenarios, in order to assess the performances of the proposed approach. The collected data are also processed in real time using a computer embedded in the car trunk. As the method is designed to be efficiently implemented on parallel, it is implemented on GPU, using nvidia CUDA. The grid resolution used is of 750x300, and the number of particles is 262144 (0.86 particle/ cell on average). The complete chain for perception (including lidar to grid mapping, HSBOF and TTC estimation) runs in real time at 15 FPS on a nVidia Quadro 2000 graphic board.

Methodology

The experiments focus on validating the approach in different contexts. On the one side, real on-road data is used to asses the realism of the approach in complex scenarios. One the other side, it is not possible to generate real risky situations on open road. Therefore, collision are performed with soft objects in a closed environment, to simulate risky driving situations.

Beyond the various contexts, we propose to use our risk estimation strategy for various applications. Therefore, we consider 3 different time horizons corresponding to 3 different TTC values: A1 -pre-crash: a collision is likely to appear in less than t 1 = 500ms. In this context, the driver is no longer able to react and the collision detection is used to prepare the vehicle to the impact; A2 -collision mitigation: a collision is likely to appear between t 1 and t 2 = 1s

In this case, the vehicle can perform an automatic action (e.g. automatic braking) in order to mitigate the danger or avoid the impact; A3 -dangerous driving: a collision is likely to appear between t 2 and t 3 = 1.5s.

In this last case, the assistance system can warn the driver of a potential risk and the driver may react in a way or another

First experiment : dynamic collisions

The first set of experiments was performed in a mostly controlled environment. On a parking lot, the experimental vehicle was used to model its environment and predicts potential collisions. To simulate obstacles, 3 bouncy balls of 55, 65 and 75 centimeter of diameter were thrown around the vehicle and on its path. The main advantage of the use of balls is the immediate validation of the existence of a risk. Indeed, on real road data, as actual collisions are not conceivable, the ground truth for risk is a hard task as it implies a subjective assessment. In this experiment, the vehicle could hit the obstacles, the collisions being annotated, confirming a preceding risk without ambiguity. Many scenarios were staged, including:

immobile or moving vehicle immobile or moving obstacles various number of obstacles various trajectory options -impacts or close crossings limited visibility (late appearance of obstacles, occlusions, intersection during curves)

Estimated risk values Figure 5-a shows an example of the estimated risk values while the balls are thrown toward the vehicle. The impacts are annotated as the vertical purple lines. It appears that for all the impacts, a high risk value is properly predicted in advance. Figure 5-b shows in more detail an example of risk prediction before the impact. As expected, the risk corresponding to T T C ≤ 1.5s rises first (between 1.5s and 2s before the collision), followed by the risk values corresponding to T T C ≤ 1.0s, and T T C ≤ 0.5s. Figure 5-c shows an example where the ball crosses close to the car without hitting it. In this case, the risk corresponding to T T C ≤ 1.5s is high and then decreases as soon as the uncertainty about a possible collision becomes very low. Alerts generation The experiment with the balloons has been repeated about a hundred times, including cases where the car hits the balloons and other where the balls come very close to the car without hitting it. This provides the opportunity to evaluate possible thresholds for generating alerts for the different values It appears on the precision/recall curves that efficient detectors can be made by selectionning the adequate threshold for every time horizon.

Second experiment : the mannequin

In order to go further more realistic scenarios, we designed a fake pedestrian in the form of a tissue mannequin filled with bubble wrap. This allows to make collisions with a human-shaped objects. Figure 7 shows the mannequin and the system used for hanging it. This system is equipped with a runner, in order to permit lateral displacements of the pedestrian. Over the different scenarios of car velocity, pedestrian motion and occlusion experimented, the system successfully generated risk alarms every time. Figure 8 shows examples of those experiments with the mannequin. The first examples corresponds to the sudden appearance of a pedestrian which was hidden by the preceding car. The risk map then shows immediately a important risk at the position of the pedestrian. It means that the system is reactive enough to detect a risk appearing at the last moment before the collision. The second one corresponds to the road crossing of a moving pedestrian. As the system estimate its motion, a risk appears on the risk map even before the pedestrian is on the car trajectory. Figure 9 shows the output of the HSBOF and the extraction of cells with a high shot-term probability of collision. On the first example, the driver arrives too quickly at the traffic light, the car ahead then generates a collision risk corresponding to a A3 alert. On the second example, we can see that the system rightly detects the risk of collision with the pedestrian crossing the road, even through he is not already on the trajectory of the vehicle. Indeed, as the occupancy dynamics are inferred by the HSBOF, the method predicts the likely impact.

Figure 10 shows an example of the risk values estimated during a urban driving sequence. Since we are driving safely, the risk corresponding to T T C ≤ 1.0s remains very low. There are a few peaks for the risk corresponding to T T C ≤ 1.5s. These happen at the end of turning maneuvers where the predicted trajectory of the car is supposed to be circular, while the driver is about to go back to a straight trajectory. Thus, these cases do not generate high values of risk for shorter TTC values. Alarms generation on-road In order to evaluate the relevance of the approach for pre-crash applications, we focus on the A1 alerts. The car is driven in various road environment (downtown, highway, national roads) to collect data. A dataset containing 27000 frames, including complex situations is recorded and the algorithm is applied, with a variable alarm threshold. The curve on figure 11 shows the number of A1 alarms with respect to the risk threshold. With a risk threshold set to 0.8 (which is relevant with respect to the curves shown on figure 6) only 4 alarms were generated.

Conclusion

In this paper, we presented a method for short-term collisions-risk estimation. The method takes advantage of the ability of the HSBOF algorithm to provide both an occupancy grid and accurate estimation of velocities over an occupancy grid. The use of this information, combined with the prediction of the vehicle's trajectory, allows to estimate Time To Collision for each cell of the grid. The methods shows several advantages, like the ability to have probabilistic estimation of collision risk, the ability to be implemented efficiently in parallel, and also to estimate TTC in highly dynamic environment, without any object tracking algorithm. A comprehensive evaluation of the algorithm has been proposed, based on both real road data and controlled situation. This evaluation showed that the method was able to perform well. In controlled environments the system is very reactive and every collision is properly detected while in difficult scenarios the system do not generate abusive alarms. There is still some work to do to deal with the remaining errors and to validate more extensively the algorithm on hundreds of kilometers, but the results obtained here from the experiments are promising. Therefore, we project to couple this perception layer with prototypes of Human-Machine-Interfaces for testing there use in Advanced Driver Assistance Systems.

Fig. 1 .

 1 Fig. 1. HSBOF data representation. A dense grid is used to represent the occupancy while for each cell velocity is split among a static field and a set of particles.

Fig. 2 .

 2 Fig.2. HSBOF algorithm summary. From sensor data instantaneous occupancy grids are successively computed. Those observations are integrated in a Bayesian filter in which coexist and jointly adapt two models, a static grid and a dynamic set of moving particles. The result is obtain by their combination, which provides a filtered occupancy grid as well as inferred motion distributions for cells.

Fig. 4 .

 4 Fig. 4. Experimental platform : Lexus LS600h car equipped with two Ibeo Lux Lidars and cameras.

Fig. 5 .

 5 Fig. 5. (a) Evolution of the collision risk for TTC ≤ 0.5s, 1.0s, 1.5s, in the case of multiple balls thrown toward the vehicle. (b) detail when a collision happens. (c) detail in case of the ball passing very near without collision. For 1.5s, a significant risk is estimated, while for shorter TTC nothing is detected, as expected.

Fig. 7 .Fig. 8 .

 78 Fig. 7. (a) Fake pedestrian used for experiments. (b) The mannequin is attached to a system with a runner, in order to allow lateral displacements.

Fig. 9 .

 9 Fig. 9. First HSBOF and risk estimation results. Left image: output of the HSBOF process with the filtered occupancy grid (white for empty space, black for occupied space) and the velocity field (from each moving cell are drawn red rays representing the velocities). Right image: dangerous cells are reprojected in the camera image.

4. 5

 5 Third experiment: On-road evaluation HSBOF and risk estimation We have performed experiments on various recorded sequences of road journey, in highway and downtown environments.

Fig. 10 .

 10 Fig. 10. Example of a sequence of driving downtown.

Fig. 11 .

 11 Fig. 11. Alarms generated for risky situations at T T C ≤ 0.5s while driving in crowded urban situations, with respect to the risk threshold (tested on 27000 frames).

Acknowledgment

The authors would like to thank Toyota Motor Europe for their collaboration and their continuous support on the Lexus car. The authors acknowledge the contribution of the "Institut de Recherche Technologique NanoElec" which has been founded by the french program "Investissement d'Avenirs" ANR-10-AIRT-05. Our thanks are also given to Nicolas Turro, Laurence Boissieux and Jean-François Cuniberto for their assistance in setting up the experiments.

of TTC. Figure 6 shows how the precision and recall rates change with respect to the threshold. Fig. 6. Precision (red dots) and recall (green triangles) of risk detection according to the threshold selected. From left to right, those are given for a estimation over 0.5s, 1s and 1.5s.


Abderemane Morame
email: abderemane.morame@univ-nantes.fr

Francoise Truc
email: francoise.truc@ujf-grenoble.fr

Semiclassical analysis for a Schrödinger operator with a U(2) artificial gauge: the periodic case

Keywords: semiclassical asymptotic, spectrum, eigenvalues, Schrodinger, periodic potential, BKW method, width of the first band, magnetic field AMS MSC 2000 : 35J10, 35P15, 47A10, 81Q10, 81Q20

We consider a Schrödinger operator with a Hermitian 2x2 matrix-valued potential which is lattice periodic and can be diagonalized smoothly on the whole R n . In the case of potential taking its minimum only on the lattice, we prove that the well-known semiclassical asymptotic of first band spectrum for a scalar potential remains valid for our model.

Semiclassical analysis for a Schrödinger operator with a U(2) artificial gauge: the periodic case 1 Introduction

Schrödinger operators with periodic matrix-valued potentials appear in many models in physics. Such models have been used recently to describe the motion of an atom in optical fields ( [Co], [Co-Da], [Da-al]), see also [Ca-Yu]. The aim of this paper is to investigate their spectral properties using semiclassical analysis. We focus on the first spectral band and assume that the potential has a non degenerate minimum. The Schrödinger operators with a non-Abelian gauge potential are Hamiltonian operators on L 2 (R n ; C m ) of the following form :

H h = h 2 n k=1 (D x k I -A k ) 2 + V + hQ + h 2 R = P h (x, hD) .
(1.1)

The classical symbol of P h (x, hD), P h (x, ξ), for (x, ξ) ∈ R n × R n , is given by

P h (x, ξ) = n k=1 (ξ k I -hA k (x)) 2 + ih 2 ∂ x k A k (x) + V (x) + hQ(x) + h 2 R(x) ,
(1.2) I is the identity m × m matrix, V, Q, R and the A k are hermitian m × m matrix with smooth coefficients and Γ periodic:

A k = (a k,ij (x)) 1≤i,j≤m , V = (v ij (x)) 1≤i,j≤m , Q = (q ij (x)) 1≤i,j≤m , R = (r ij (x)) 1≤i,j≤m , a k,ij , v ij , q ij , r ij ∈ C ∞ (R n ; C), a k,ji = a k,ij , v ji = v ij , q ji = q ij , r ji = r ij a k,ij (x -γ) = a k,ij (x), v ij (x -γ) = v ij (x), q ij (x -γ) = q ij (x) and r ij (x -γ) = r ij (x) ∀γ ∈ Γ ;                (1.3) Γ is a lattice of R n , Γ = { n k=1 m k β k ; m k ∈ Z},
β 1 , β 2 , . . . , β n ∈ R n form a basis, det(β 1 , β 2 , . . . , β n ) = 0. We use the notation D = (D x 1 , . . . , D xn ) where

D x k = -i∂ x k , k = 1...n, so D 2 = -∆ is the Laplacian operator on L 2 (R n ).
The dual basis {β ⋆ 1 , . . . , β ⋆ n } of the reciprocal lattice Γ ⋆ , is the basis of R n defined by the relations

β ⋆ j .β k = 2πδ jk : Γ ⋆ = { n k=1 m k β ⋆ k ; m k ∈ Z}.
The fundamental cell, the Wigner-Seitz cell,

W n = { n k=1 x k β k ; x k ∈] - 1 2 , 1 2 [} ,
will be identified with the n-dimensional torus T n = R n /Γ and the dual cell, the Brillouin zone, is defined by

B n = { n k=1 θ k β ⋆ k ; θ k ∈] - 1 2 , - 1 2 [} .
We will identify L 2 (T n ; C m ) with Γ periodic functions of L 2 loc (R n ; C m ) provided with the norm of L 2 (W n ; C m ). In the same way the Sobolev space W k (T n ; C m ), with k ∈ N, may be identified with Γ periodic functions of

W k loc (R n ; C m ) provided with the norm of W k (W n ; C m ).
By Floquet theory, (see [Ea] or [Re-Si] ), we have

H h = ⊕ B n H h,θ dθ , with H h,θ the partial differential operator P h (x, h(D -θ)) on L 2 (T n ; C m ). The ellipticity of P h (x, h(D -θ)) implies that the spectrum of H h,θ is discrete sp(H h,θ ) = {λ h,θ j ; j ∈ N ⋆ }, λ h,θ 1 ≤ λ h,θ 2 ≤ . . . ≤ λ h,θ j ≤ λ h,θ j+1 ≤ . . . (1.4) each λ h,θ j
is an eigenvalue of finite multiplicity and each eigenvalue is repeated according to its multiplicity.

(When m = 1 and V = Q = R = A k = 0, ( 1 |T n | e iω.x ) ω∈Γ ⋆ is the Hilbert basis of L 2 (T n ) which is composed of eigenfunctions of h 2 (D -θ) 2 ).
The Floquet theory guarantees that

sp(H h ) = θ∈B n sp(H h,θ ) = ∞ j=1 b h j , (1.5)
where b h j denotes the j-th band b h j = {λ h,θ j , θ ∈ B n }. In the sequel h 0 will be a non negative small constant, h will be in ]0, h 0 [, and any non negative constant which doesn't depend on h will invariably be denoted by C.

Preliminary: the artificial gauge model

We will be interested in the model of artificial gauge considered in [Co], [Co-Da] and [Da-al]

m = 2, V = vI + W, A k = Q = R = 0, ∀k, W = w.σ, with w = (w 1 , w 2 , w 3 ), v and the w j are in C ∞ (R n ; R), (2.1)
we denote σ = (σ 1 , σ 2 , σ 3 ), where the σ j are the Pauli matrices.

Let us remark that

V = vI + W, W = w.σ, W 2 = |w| 2 I . (2.2)
In the sequel we will assume that

|w| > 0 v(x) -|w(x)| has a unique non degenerate minimum on T n . (2.3)
Due to the invariance of the Laplacian by translation and by the action of O(n), we can assume, up to a composition by a translation of the potentials, that

v(γ) -|w(γ)| < v(x) -|w(x)|, ∀x ∈ R n \ Γ and ∀γ ∈ Γ, v(x) -|w(x)| = E 0 + n k=1 τ 2 k x 2 k + O(|x| 3 ), as |x| → 0,      (2.4) (τ k > 0, ∀k).
There exists U ∈ U(2), ( a unitary 2 × 2 matrix), such that

U ⋆ V U = V = v -|w| 0 0 v + |w| . (2.5)
As |w| never vanishes, U = U(x) can be chosen smooth and Γ periodic:

U = u 11 u 12 u 21 u 22 ∈ C ∞ (T n ; U(2));
for example u 11 = 1 2 |w|(|w| -Re((w 1 + iw 2 )e -iθ )) (w 3 -|w| + e iθ (w 1 -iw 2 )),

u 21 = 1 2 |w|(|w| -Re(w 1 + iw 2 )e -iθ )) (w 1 + iw 2 -e iθ (w 3 + |w|)), u 12 = u 21 , u 22 = -u 11 and θ = χ( w 2 2 + w 2 3 |w| 2 ) π 2 ,
where χ(t) is a smooth function on the real line, 0 ≤ χ(t) ≤ 1,

χ(t) = 1 when |t| ≤ 1/4 and χ(t) = 0 when |t| ≥ 1/2. So U = (α, β, ρ).σ + iδσ 0 , with (α, β, ρ, δ) ∈ C ∞ (T n ; S 3 ); (2.6)
σ 0 is the 2 × 2 identity matrix and S 3 is the unit sphere of R 4 . When w 1 + iw 2 = 0 or when w 3 < 0, one can choose U such that δ = 0 by

taking (α, β, ρ) = 1 2|w| (- w 1 |w| -w 3 , - w 2 |w| -w 3 , |w| -w 3 ).
Firstly let us expand the formula of the operator

H h = U ⋆ H h U = h 2 D 2 I + U ⋆ V U -2ih 2 n k=1 (U ⋆ ∂ x k U)D x k -h 2 U ⋆ ∂ 2 x k U
which can be rewritten as

H h = U ⋆ H h U = h 2 n k=1 (D x k I -A k ) 2 + V + h 2 R , (2.7)
where

A k = iU ⋆ ∂ x k U : A k = [(∂ x k α, ∂ x k β, ∂ x k ρ)∧(α, β, ρ)+(δ∂ x k α-α∂ x k δ, δ∂ x k β-β∂ x k δ, δ∂ x k ρ-ρ∂ x k δ)].σ , (2.8) and R = n k=1 (U ⋆ ∂ x k U) 2 + (∂ x k U ⋆ ).(∂ x k U) .
(2.9)

So we can assume that H h is of the form (1.1) with m = 2, Q = 0, A k and R given by (2.8) and (2.9), with U defined by (2.6), and V = V a diagonal matrix given by (2.5).

Theorem 2.1 Under the above assumptions, the first bands b h j , j = 1, 2, . . . , of H h are concentrated around the value hµ j + E 0 j = 1, 2, . . . , in the sense that, there exist N 0 > 1 and h 0 > 0 such that

distance(hµ j + E 0 , b h j ) ≤ Ch 2 , ∀j < N 0 and ∀h, 0 < h < h 0 ,
where

µ j = n k=1 (2j k + 1)τ k , j k ∈ N, the (µ ℓ ) ℓ∈N ⋆ is the increasing sequence of the eigenvalues of the harmonic oscillator -∆ + n k=1 τ 2 k x 2 k .
3 Proof of Theorem 2.1

Proof. According to the above discussion, we can assume that

H h = P h (x, hD) , with P h (x, hD) = P h 11 (x, hD) P h 12 (x, hD) P h 21 (x, hD) P h 22 (x, hD) , (3.1) with P h 11 (x, hD) = h 2 (D -a .,11 (x)) 2 + v(x) -|w(x)| + h 2 r 11 (x) P h 22 (x, hD) = h 2 (D + a .,11 (x)) 2 + v(x) + |w(x)| + h 2 r 22 (x) P h 12 (x, hD) = -h 2 a .,12 (x).(D + a .,11 (x)) -h 2 a .,12 (x).(D -a .,11 (x)) + ih 2 div(a .,12 (x)) + h 2 r 12 (x) P h 21 (x, hD) = -h 2 a .,21 (x)(D -a .,11 (x)) -h 2 a .,21 (x).(D + a .,11 (x)) + ih 2 div(a .,21 (x)) + h 2 r 21 (x) .                (3.2) (D = (D x 1 , D x 2 , . . . , D xn ) and a .,ij (x) = (a 1,ij (x), a 2,ij (x), . . . , a n,ij (x)) .) (We used that a .,22 = -a .,11 ).
Let us denote by H h,θ 11 and H h,θ 22 the operators associated with

P h 11 (x, h(D -θ)) and P h 22 (x, h(D -θ)) on L 2 (T n ; C). Then, if c 0 = min |w(x)| and c 1 = max R(x) , sp(H h,θ 11 ) ⊂ [E 0 -h 2 c 1 , +∞[ and sp(H h,θ 22 ) ⊂ [E 0 -h 2 c 1 + 2c 0 , +∞[.
To prove the theorem it is then enough to prove the proposition below.

Proposition 3.1 Let us consider a constant c, 0 < c < c 0 . Then there exists C 0 > 0 such that, for any E h ∈] -∞, E 0 + 2c[, we have E h ∈ sp(H h,θ ) ⇒ distance(E h , sp(H h,θ 11 ) ≤ C 0 h 2 E h ∈ sp(H h,θ 11 ) ⇒ distance(E h , sp(H h,θ )) ≤ C 0 h 2 . (3.3) Proof. For such E h , (H h,θ
22 -E h ) -1 exists and, thanks to semiclassical pseudodifferential calculus of [Ro] (see also [Di-Sj] ), for

h 0 > 0 small, if 0 < h < h 0 then (H h,θ 22 -E h ) -1 L 2 (T n ) + h(D -θ)(H h,θ 22 -E h ) -1 L 2 (T n ) + (H h,θ 22 -E h ) -1 h(D -θ) L 2 (T n ) + h(D -θ)(H h,θ 22 -E h ) -1 h(D -θ) L 2 (T n ) ≤ C, and then P h 12 (x, h(D -θ))(H h,θ 22 -E h ) -1 P h 21 (x, h(D -θ)) L 2 (T n ) ≤ h 2 C. So if E h ∈ sp(H h,θ ), then u h = (u h 1 , u h 2 ) = (0, 0) is an eigenfunction of H h,θ associated with E h iff H h,θ 11 u h 1 + P h 12 (x, h(D -θ))u h 2 = E h u h 1 u h 2 = -(H h,θ 22 -E h I) -1 P h 21 (x, h(D -θ))u h 1 .
(3.4)

In fact E h ∈] -∞, E 0 + c[ will be an eigenvalue of H h,θ iff there exists

u h 1 in the Sobolev space W 2 (T n ; C), u h 1 L 2 (T n ) = 0, such that H h,θ 11 u h 1 -P h 12 (x, h(D -θ))(H h,θ 22 -E h I) -1 P h 21 (x, h(D -θ))u h 1 = E h u h 1 , then we get the first part of Proposition 3.1. If E h is an eigenvalue of H h,θ
11 satisfying the assumption of Proposition 3.1, and u h 1 an associated eigenfunction, then with

u h = (u h 1 , -(H h,θ 22 -E h ) -1 P h 21 (x, h(D -θ))u h 1 ), one has (H h,θ -E h I)u h L 2 (T n ;C 2 ) = P h 12 (x, h(D -θ))(H h,θ 22 -E h ) -1 P h 21 (x, h(D -θ))u h 1 L 2 (T n ;C) ≤ h 2 C u h L 2 (T n ;C 2 )
, we get the second part of Proposition 3.1.

Theorem 2.1 follows from Proposition 3.1 and , , [He-Sj-1] and [He-Sj-2] results, (see also [He]), which guarantee that the sequence of eigenvalues of

H h,θ 11 , (λ j (H h,θ 11 )) j∈N ⋆ satisfies ∀N 0 > 1, ∃h 0 > 0, C 0 > 0 s.t. ∀h, 0 < h < h 0 and ∀j ≤ N 0 , |λ j (H h,θ 11 ) -(hµ h j + E 0 )| ≤ C 0 h 2

Asymptotic of the first band

For any real Lipschitz Γ periodic function φ, and for any u ∈ W 2 (T n ; C 2 ), we have the identity

Re < P h (x, h(D -θ))u | e 2φ/h u > L 2 (T n ;C 2 ) = n k=1 h 2 ((D x k -θ k )I -A k )e φ/h u 2 L 2 (T n ;C 2 ) + < ( V -|∇φ| 2 I + h 2 R)u | e 2φ/h u > L 2 (T n ;C 2 ) .          (4.1)
This identity enables us to apply the method used in [He-Sj-1], (see also [He] and [Ou]). We define the Agmon [Ag] distance on

R n d(y, x) = inf γ 1 0 v(γ(t)) -|w(γ(t))| -E 0 | γ(t)|dt , (4.2)
the inf is taken among paths such that γ(0) = y and γ(1) = x.

For common properties of the Agmon distance, one can see for example [Hi-Si].

We will use that, for any fixed

y ∈ R n , the function d(y, x) is a Lipschitz function on R n and |∇ x d(y, x)| 2 ≤ v(x) -|w(x)| -E 0 almost everywhere on R n .
Using that the zeros of v(x) -w(x) -E 0 are the elements of Γ and are non degenerate, we get that the real function

d 0 (x) = d(0, x) satisfies, (see [He-Sj-1]), |∇d 0 (x)| 2 = v(x) -|w(x)| -E 0 in a neighbourhood of 0.
We summarize the properties of the Agmon distance we will need:

i) ∃ R 0 > 0 s.t. d 0 (x) ∈ C ∞ (B 0 (R 0 )) ii) |∇d 0 (x)| 2 = v(x) -|w(x)| -E 0 , ∀x ∈ B 0 (R 0 ) iii) |∇d 0 (x)| 2 ≤ v(x) -|w(x)| -E 0 iv) |∇d Γ (x)| 2 ≤ v(x) -|w(x)| -E 0        (4.3) where d 0 (x) = d(0, x), B 0 (r) = {x ∈ R n ; d 0 (x) < r} and d Γ (x) = d(Γ, x) = min ω∈Γ d(ω, x).
The least Agmon distance in Γ is

S 0 = inf 1≤k≤n d 0 (β k ) = inf ρ =ω, (ω,ρ)∈Γ 2 d(ω, ρ) .
(4.4)

The Agmon distance on T n , d T n (., .), is defined by its Γ-periodic extension on

(R n ) 2 d T n (y, x) = min ω∈Γ d(y, x + ω).
Then

S 0 2 = sup r {r > 0 s.t. {x ∈ T n ; d T n (x 0 , x) < r} is simply connected} , (4.5)
where x 0 is the single point in T n such that v(x 0 ) -|w(x 0 )| = E 0 . The Γ-periodic function on R n , d Γ (x) is the one corresponding to the extension of d T n (x 0 , x). If λ h,θ is an eigenvalue of H h,θ and if u h,θ is an associated eigenfunction, then by (4.1) one gets as in the scalar case considered in ) and ,

n k=1 h 2 ((D x k -θ k )I -A k )e φ/h u h,θ 2 L 2 (T n ;C 2 ) + < [ V -|∇φ| 2 I + h 2 R -λ h,θ I] + u h,θ | e 2φ/h u h,θ > L 2 (T n ;C 2 ) =< [ V -|∇φ| 2 I + h 2 R -λ h,θ I] -u h,θ | e 2φ/h u h,θ > L 2 (T n ;C 2 ) ,          (4.6) so, when φ(x) = d T n (x 0 , x), necessarily λ h,θ -E 0 + O(h 2 ) > 0, and if h/C < λ h,θ -E 0 < hC, then u h,θ is localized in energy near x 0 , for any η ∈]0, 1[, ∃C η > 0 such that n k=1 h 2 ((D x k -θ k )I -A k )e ηd T n (x 0 ,x)/h u h,θ 2 L 2 (T n ;C 2 ) +(1 -η 2 ) < (v -|w| -E 0 )u h,θ | e 2ηd T n (x 0 ,x)/h u h,θ > L 2 (T n ;C 2 ) ≤ hC η {x∈T n ; d T n (x 0 ,x)< √ hC} |u h,θ (x)| 2 dx.          (4.7)
Let Ω ⊂ T n an open and simply connected set with smooth boundary satisfying, for some η, 0 < η < S 0 /2,

{x ∈ T n ; d T n (x 0 , x) < S 0 -η 2 } ⊂ Ω ⊂ {x ∈ T n ; d T n (x 0 , x) < S 0 /2} (4.8)
Let H h Ω be the selfadjoint operator on L 2 (Ω; C 2 ) associated with P h (x, hD) with Dirichlet boundary condition. We denote by (λ j (H h Ω )) j∈N ⋆ the increasing sequence of eigenvalues of H h Ω . Using the method of [He-Sj-1], we get easily the following results Theorem 4.1 For any η, 0 < η < S 0 /2, there exist h 0 > 0 and N 0 > 1 such that, if 0 < h < h 0 and j ≤ N 0 , then 2h) . For the first band, we have the following improvement

∀θ ∈ B n , 0 < λ j (H h Ω ) -λ h,θ j ≤ Ce -(S 0 -η)/(2h) ; (4.9) so the length of the band b h j satisfies |b h j | ≤ Ce -(S 0 -η)/(
|b h 1 | ≤ Ce -(S 0 -η)/h . (4.10)

Sketch of the proof.

As Ω is simply connected and the one form θdx is closed, there exists a smooth real function ψ θ (x) on Ω such that e -iψ θ (x) P h (x, hD)e iψ θ (x) = P h (x, h(D-θ)) : the Dirichlet operators on Ω, H h Ω and H h,θ Ω associated to P h (x, hD) and P h (x, h(Dθ)) are gauge equivalent, so they have the same spectrum.

Therefore the min-max principle says that

0 < λ j (H h,θ Ω ) -λ h,θ j = λ j (H h Ω ) -λ h,θ j .
But the exponential decay of the eigenfunction ϕ h,θ j (x) associated with λ h,θ j , given by (4.7) implies that

(P h (x, h(D -θ)) -λ h,θ j )χϕ h,θ j (x) L 2 (Ω;C 2 ) ≤ Ce -(S 0 -η+ǫ)/(2h) ,
for some ǫ > 0, and for a smooth cut-off function χ supported in Ω and 2h) . This achieves the proof of (4.9).

χ(x) = 1 if d T n (x 0 , x) ≤ (S 0 -η + ǫ)/2. So distance(λ h,θ j , sp(H h Ω )) ≤ Ce -(S 0 -η+ǫ)/(
Let us denote E h,θ , (respectively E h,Ω ), the first eigenvalue λ h,θ 1 , (respectively λ 1 (H h Ω )), and ϕ h,θ (x), (respectively ϕ h,Ω (x)), the associated normalized eigenfunctions. Let χ be a cut-off function satisfying the same properties as before. Then

P h (x, h(D-θ))(e -iψ θ (x) χ(x)ϕ h,Ω (x)) = λ 1 (H h
Ω )e -iψ θ (x) χ(x)ϕ h,Ω (x) + e -iψ θ (x) r h 0 (x) with, thanks to the same identity as (4.7) for Dirichlet problem on Ω,

r h 0 L 2 (T n ;C 2 ) ≤ Ce -S 0 -η 2h .
The same argument used in [He-Sj-1], (see also [He]), gives this estimate

|E h,θ -E h,Ω -< r h 0 | χϕ h,Ω > L 2 (T n ;C 2 ) | ≤ Ce -(S 0 -η)/h . As τ h =< r h 0 | χϕ h,Ω > L 2 (T n ;C 2 ) does not depend on θ, so ∀θ ∈ B n , |E h,θ -E h,Ω -τ h | ≤ Ce -(S 0 -η)/h
this estimate ends the proof of (4.10) As for the tunnel effect in [He-Sj-1] and , we have an accuracy estimate for the first band, like the scalar case in and in [Ou] (see also [He]).

Theorem 4.2 There exists h 0 > 0 such that, if 0 < h < h 0 then

|b h 1 | ≤ Ce -S 0 /h .
Sketch of the proof. Instead of comparing H h,θ with an operator defined in a subset of T n , we have to work on the universal cover R n of T n . We take Ω ⊂ R n an open and simply connected set with smooth boundary satisfying, for some η 0 , 0 < η 0 < η 1 < S 0 /2, B 0 ((S 0 + η 0 )/2) ⊂ Ω ⊂ B 0 ((S 0 + η 1 )/2).

(4.11)

So Ω contains the Wigner set W n , more precisely

W n ⊂ Ω ⊂ 2W n and Ω ∩ Γ = {0}.
We let also denote H h Ω the Dirichlet operator on L 2 (Ω; C 2 ) associated with P h (x, hD), and E h Ω its first eigenvalue. The associated eigenfunction is also denoted by ϕ h,Ω (x).

In the same way as to get (4.7), we have

n k=1 h 2 (D x k I -A k )e d 0 (x)/h ϕ h,Ω 2 L 2 (Ω;C 2 ) ≤ hC B 0 ( √ hC)
|ϕ h,Ω (x)| 2 dx, (4.12)

then the Poincaré estimate gives

Ω e 2d 0 (x)/h |ϕ h,Ω (x)| 2 dx ≤ h -1 C B 0 ( √ hC) |ϕ h,Ω (x)| 2 dx. (4.13)
Let χ a smooth cut-off function satisfying

χ(x) = 1 if d 0 (x) ≤ (S 0 + η 0 )/2 and χ(x) = 0 if x / ∈ Ω.
Then the function

ϕ h,θ (x) = ω∈Γ e iθ(x-ω) χ(x -ω)ϕ h,Ω (x -ω)
is Γ-periodic and satisfies

(P h (x, h(D -θ)) -E h Ω )ϕ h,θ (x) = r h,θ and r h,θ L 2 (W n ;C 2 ) ≤ Ce -(S 0 +η 0 )/(2h) ϕ h,θ L 2 (W n ;C 2 ) (4.14) and < r h,θ | ϕ h,θ > L 2 (W n ;C 2 ) = ω,ρ∈Γ 0 e iθ(ρ-ω) W n ([P h (x, hD); χ]ϕ h,Ω )(x -ω).(χϕ h,Ω )(x -ρ)dx with Γ 0 = {0, ±β 1 , . . . , ±β n } and [P h (x, hD); χ] = -2h 2 i n k=1 ∂ x k χ(D x k I -A k ) -h 2 ∆χI. So 1 ϕ h,θ 2 L 2 (W n ;C 2 ) < r h,θ | ϕ h,θ > L 2 (W n ;C 2 ) ≤ Ce -S 0 /h . (4.15)
The proof comes easily from (4.12) and (4.13) as in [Ou] or in [He].

Using the same argument of as in the proof of (4.10), we get that

|E h,θ -E h Ω -τ h,θ | ≤ Ce -(S 0 +η 0 )/h ) , (4.16) with τ h,θ = 1 ϕ h,θ 2 L 2 (W n ;C 2 ) < r h,θ | ϕ h,θ > L 2 (W n ;C 2 ) .
Theorem 4.2 follows from (4.15) and (4.16)

B.K.W. method for the Dirichlet ground state

Let Ω be an open set satisfying (4.8). more precisely Ω ⊂ R n an open, bounded and simply connected set with smooth boundary satisfying, for some

η 1 and η 2 , 0 < η 1 < η 2 < S 0 /2, {x ∈ R n ; d 0 (x) < S 0 -η 2 2 } ⊂ Ω ⊂ {x ∈ R n ; d 0 (x) < S 0 -η 1 2 } (5.1) Theorem 5.1 The first eigenvalue E h,Ω = λ 1 (H h Ω ) of the Dirichlet operator H h Ω
admits an asymptotic expansion of the form

E h,Ω ≃ ∞ j=0 h j e j ,
and if S 0 -η 1 is small enough, the associated eigenfunction ϕ h,Ω has also an asymptotic expansion of the form

ϕ h,Ω = e -φ/h (f + h , f - h ) , f ± h ≃ ∞ j=0 h j f ± j , (f - 0 = 0) .
As usual Proof. When the gauge potential matrix is identified with the one form

e 0 = E 0 , e 1 = τ 1 , e 2 = r 11 (0) + n k=1 |a k,11 (0)| 2 , (5.2
Adx = n k=1 A k (x)dx k ,
its curvature form appears to be the related magnetic field B = dA + A ∧ A :

B = 1≤j<k≤n (∂ x j A k (x)-∂ x k A j (x))dx j ∧dx k + 1≤j<k≤n (A j (x)A k (x)-A k (x)A j (x))dx j ∧dx k .
For our purpose, only the vector magnetic potential a .,11 is significant. We will work with Coulomb vector gauge a .,11 : div(a We look for an eigenvalue E h ≃ ∞ j=0

       (5.6) Let us write e -φ/h P h (x, hD)(e -φ/h f h ) = W 0 (x)f h + hW 1 (x, D)f h + h 2 W 2 (x, D)f h (5.7) with W 0 (x) = V (x) -|∇φ(x)| 2 I W 1 (x, D) = ∆φI + 2i n k=1 ∂ x k φ(D x k I -A k ) W 2 (x, D) = n k=1 (D x k I -A k ) 2 + R(x). So W 1 (x, D) = ∆φ -2i∇φ.(i∇ + a .,
h j e j and an associated eigenfunction

f h ≃ ∞ j=0 h j f j , so e φ/h (P h (x, hD) -E h I)(e -φ/h f h ) ≃ ∞ j=0 h j κ j (We used that f + 0 (0) = 1). So κ 2 = 0 implies e 2 = (i∇ + a .,11 ) 2 f + 0 (0) + r 11 (0) , 2∇φ(x).∇f + 1 (x) = -(∆φ(x) -e 1 -2ia .,11 (x).∇φ(x))f + 1 (x) -( (i∇ + a .,11 ) 2 + e 2 -r 11 (x))f + 0 (x) + 2i(a .,12 (x).∇φ(x))f - 1 (x) f - 2 (x) = 1 2|w(x)| [2i(a .,21 (x).∇φ(x))f + 1 (x) -2∇φ(x).∇f - 1 (x) -2i(a .,11 (x).∇φ(x))f - 1 (x) -2ia .,21 (x).∇f + 0 (x) -i(div(a .,21 (x))f + 0 (x) -r 21 (x)f + 0 (x)] .
4) Term of order j> 2. We assume that e ℓ for ℓ = 0, 1, . . . , j-1, the functions f ± ℓ (x) for ℓ = 0, 1, . . . , j-2, and the one f - j-1 (x) are well-known, f + ℓ (0) = 0 when 0 < ℓ < j -1. The equation κ

+ j = 0 becomes 2∇φ(x).∇f + j-1 (x) + (∆φ(x) -e 1 -2ia .,11 (x).∇φ(x))f + j-1 (x) = 2i(a .,21 (x).∇φ(x))f - j-1 (x) + j-3 ℓ=0 e j-ℓ f + ℓ (x) -( (i∇ + a .,11 ) 2 -e 2 + r 11 (x) )f + j-2 (x) -( r 12 (x) + idiv(a 21 ) )f - j-2 (x) -2ia 21 .∇f - j-2 (x)              (5.11)
This equation has a unique solution f

+ j-1 (x) such that f + j-1 (0) = 0 iff e j = ( (i∇ + a .,11 ) 2 + r 11 (0) )f + j-2 (0) + r 12 (0))f - j-2 (0) -2ia .,21 (0).∇)f - j-2 (0) -j-3 ℓ=0 e j-ℓ f + ℓ (0) .
(5.12)

The equation κ - j = 0 gives

f - j (x) = 1 2|w(x)| × [-2∇φ(x).∇f - j-1 (x) + (e 1 -2ia .,11 (x).∇φ(x))f - j-1 (x) + 2ia .,21 (x).∇φ(x)f + j-1 (x) +(e 2 -(i∇-a .,11 ) 2 )f - j-2 (x)-idiv(a .,21 (x))f + j-2 (x)-2ia .,21 (x).∇f + j-2 (x)+ j-3 ℓ=0 e j-ℓ f - ℓ (x)]
5) End of the proof. Let χ(x) be a cut-off function equal to 1 in a neighbourhood of 0 and supported in Ω. Then taking χ(x)f j (x) instead of f j (x), we get a function ϕ h,Ω satisfying Dirichlet boundary condition and Theorem 5.1. The self-adjointness of the related operator ensures that the computed sequence (e j ) is real 6 Sharp asymptotic for the width of the first band

Returning to the proof of Theorem 4.2, we have to study carefully the τ h,θ defined in (4.16), using the method of performed in [He] and [Ou]. L 2 (W n ;C 2 ) W n ([P h (x, hD); χ]ϕ h,Ω )(x).(χϕ h,Ω )(x + ω)dx .

We get from the formula of [P h (x, hD); χ] and from the estimate (4.13) that

ρ + ω = - h 2 ϕ h,θ 2 L 2 (W n ;C 2 ) W n
(∇χ(x -ω)∇ϕ h,Ω )(x -ω)).(χϕ h,Ω )(x)dx + (6.2) + O(he -S 0 /h ) and

ρ - ω = - h 2 ϕ h,θ 2 L 2 (W n ;C 2 ) W n
(∇χ(x)∇ϕ h,Ω )(x)).(χϕ h,Ω )(x + ω)dx + O(he -S 0 /h ) . Theorem 6.1 Under the assumption of Theorem 4.2, if for any ω ∈ {±β 1 , . . . , ±β n } such that the Agmon distance in R n between 0 and ω is the least one, (i.e. d(0, ω) = S 0 ), there exists one or a finite number of minimal geodesics joining 0 and ω, then there exists η 0 > 0 and h 0 > 0 such that b h 1 = η 0 h 1/2 e -S 0 /h 1 + O(h 1/2 ) , ∀h ∈]0, h 0 [ . Sketch of the proof. Following the proof of splitting in [He-Sj-1] and [He], in (6. + e d(x)/h χ 0 (ϕ h,Ω -ϕ h B.K.W. ) 2 L 2 (W n ;C 2 ) ≤ h p 0 C p 0 , (6.3) where χ 0 is a cut-off function supported in W n ∩ O and equal to 1 in a neighborhood of the minimal geodesics between 0 and the ±β k . We have assumed that ϕ h,Ω L 2 (W n ;C 2 ) = 1 and then χ 0 ϕ h B.K.W. ) 2 L 2 (W n ;C 2 ) -1 = O(h p ) for any p > 0.

)

  and φ is the real function satisfying the eikonal equation|∇φ(x)| 2 = v(x) -|w(x)| -E 0 ,(5.3) equal to d(x) in a neighbourhood of 0.(r 11 and the a k,11 are defined by (1.1) and (1.3). E 0 and τ 1 are defined by (2.2) and (2.4)).

Using

  W n ;C 2 ) W n ([P h (x, hD); χ]ϕ h,Ω )(x -ω).(χϕ h,Ω )(x

  2) we can change W n into W n ∩ O, where O is any neighbourhood of the minimal geodesics between 0 and the ±β k , such that d(x) = d(0, x) ∈ C ∞ (O). In this case the B.K.W. method is valid in W n ∩ O. If ϕ h B.K.W. is the B.K.W. approximation of ϕ h,Ω in W n ∩ O,then, thanks to (4.1), for any p 0 > 0 there exists C p 0 such thath n k=1 (D x k I -A k )e d(x)/h χ 0 (ϕ h,Ω -ϕ h B.K.W. ) 2 L 2 (W n ;C 2 )

  .,11 (x)) = It is feasible thanks to the existence of a smooth real and Γ periodic function ψ(x) such that ∆ψ(x) = div(a .,11 (x)). Let O be any open set of R n (or one can take also O = T n ). Conjugation of P h (x, hD) by the unitary operator J ψ on L 2 (O; C 2 ) : .,11 (x) for a .,11 (x) -∇ψ(x) and a .,21 (x) for e 2iψ(x) a .,21 (x); the new a .,22 (x) is equal to minus the new a .,11 (x), and the new a .,12 (x) remains the conjugate of the new a .,21 (x). So by (2.8) we have a .,11 = β∇α -α∇β + δ∇ρ -ρ∇δ -∇ψ , a .,21 = e 2iψ [(ρ∇β -β∇ρ + δ∇α -α∇δ) + i(α∇ρ -ρ∇α + δ∇β -β∇δ)] a .,22 = -a .,11 , a .,12 = a .,21 ∆ψ = div(β∇α -α∇β + δ∇ρ -ρ∇δ)

	J ψ =	e iψ 0 e -iψ , 0	(5.5)
	leads to changing a
		n
		∂ x k a k,11 (x) = 0 .	(5.4)
		k=1

  .,11 ) 2 + r 11 (i∇ + a .,11 ).a .,21 + a .,21 .(i∇ -a .,11 ) + r 12 a .,21 .(i∇ + a .,11 ) + (i∇ + a .,11 ).a .,21 + r 21 (i∇ -a .,11 ) 2 + r 22

		-2i∇φ.a .,21	11 )	-2i∇φ.a .,21 ∆φ -2i∇φ.(i∇ -a .,11 )	,
	and
	W 2 (x, D) =	(i∇ + a

with κ 0 (x) = (W 0 (x) -e 0 I)f 0 (x) κ 1 (x) = (W 1 (x, D) -e 1 I)f 0 (x) + (W 0 (x) -e 0 I)f 1 (x) κ 2 (x) = (W 2 (x, D) -e 2 I)f 0 (x) + (W 1 (x, D) -e 1 I)f 1 (x) + (W 0 (x) -e 0 I)f 2 (x) κ j (x) = -e j f 0 (x) -j-3 ℓ=1 e j-ℓ f ℓ (x) + (W 2 (x, D) -

We recall that f j = (f + j , f - j ) and we want that κ j (x) = 0, ∀j. 1) Term of order 0

(5.8)

2) Term of order 1.

As |∇φ(x)| has a simple zero at x 0 = 0, the equation κ + 1 (x) = 0 can be solved only when e 1 = ∆φ(0). In this case there exists a unique function f + 0 (x) such that f + 0 (0) = 1 and κ + 1 (x) = 0. We can conclude that the study of the term of order 1 leads to

The unknown function f + 1 (x) must give κ + 2 (x) = 0 in (5.10). This equation, with the initial condition f + 1 (0) = 0, can be solved only if

As (6.2) remains valid if we change ϕ h,Ω into χ 0 ϕ h,Ω , the estimate (6.3) allows also to change ϕ h,Ω into χ 0 ϕ h B.K.W. . As a consequence, Theorem 6.1 follows easily, if in W n ∩ O, χ(x) = χ 1 (d(x)) for a decreasing function χ 1 on [0, +∞[ with compact support, equal to 1 in a neighborhood of 0. In this case (6.2) becomes

) . We remind that for any y in a minimal geodesic joining 0 to ±β k , if y = 0 and y = ±β k , then the function d(x) + d(x ∓ β k ), when it is restricted to any hypersurface orthogonal to the geodesic through y, has a non degenerate minimum S 0 at y

Acknowledgement

The authors are grateful to Bernard Helffer for many discussions.