Abdelwaheb Ifa
email: abdelwaheb.ifa@fsm.rnu.tn

Hanen Louati
email: louatihanen42@yahoo.fr

Michel Rouleux
email: rouleux@univ-tln.fr

BOHR-SOMMERFELD QUANTIZATION RULES REVISITED: THE METHOD OF POSITIVE COMMUTATORS

We revisit the well known Bohr-Sommerfeld quantization rule (BS) of order 2 for a selfadjoint 1-D h-Pseudo-differential operator within the algebraic and microlocal framework of Helffer and Sjöstrand; BS holds precisely when the Gram matrix consisting of scalar products of some WKB solutions with respect to the "flux norm" is not invertible. The interest of this procedure lies in its possible generalization to matrix-valued Hamiltonians, like Bogoliubov-de Gennes Hamiltonian. It is simplified in the scalar case by using action-angle variables.

Introduction.

Let p(x, ξ; h) be a smooth real classical Hamiltonian on T * R ; we will assume that p belongs to the space of symbols S 0 (m) for some order function m with (0.1)

S N (m) = {p ∈ C ∞ (T * R) : ∀α ∈ N 2 , ∃C α > 0, |∂ α (x,ξ) p(x, ξ; h)| ≤ C α h N m(x, ξ)}
and has the semi-classical expansion (0.2) p(x, ξ; h) ∼ p 0 (x, ξ) + hp 1 (x, ξ)

+ • • • , h → 0
We call as usual p 0 the principal symbol, and p 1 the sub-principal symbol. We also assume that p + i is elliptic. This allows to take Weyl quantization of p (0.3) P (x, hD x ; h)u(x; h) = p w (x, hD x ; h)u(x; h) = (2πh) -1 e i(x-y)η/h p( x + y 2 , η; h)u(y) dy dη so that P (x, hD x ; h) is essentially self-adjoint on L 2 (R). In case of Schrödinger operator P (x, hD x ) = (hD x ) 2 + V (x), p(x, ξ; h) = p 0 (x, ξ) = ξ 2 + V (x). We make the geometrical hypothesis of [CdV1],

namely:

Fix some compact interval I = [E -, E + ], E -< E + , and assume that there exists a topological ring A ⊂ p -1 0 (I) such that ∂A = A -∪ A + with A ± a connected component of p -1 0 (E ± ). Assume also that p 0 has no critical point in A, and A -is included in the disk bounded by A + (if it is not the case, we can always change p to -p.) That hypothesis will be referred in the sequel as Hypothesis (H).

We define the microlocal well W as the disk bounded by A + . For E ∈ I, let γ E ⊂ W be a periodic orbit in the energy surface {p 0 (x, ξ) = E}, so that γ E is an embedded Lagrangian manifold.

1

Then if E + < E 0 = lim inf |x,ξ|→∞ p 0 (x, ξ), all eigenvalues of P in I are indeed given by Bohr-Sommerfeld quantization condition (BS) that we recall here, when computed at second order: Theorem 0.1: With the notations and hypotheses stated above, for h > 0 small enough there exists a smooth function S h : I → R, called the semi-classical action, with asymptotic expansion S h (E) ∼ S 0 (E) + hS 1 (E) + h 2 S 2 (E) + • • • such that E ∈ I is an eigenvalue of P iff it satisfies the implicit equation (Bohr-Sommerfeld quantization condition) S h (E) = 2πnh, n ∈ Z. The semi-classical action consists of :

d dE γ E ∆ dt - γ E p 2 dt - 1 2 d dE γ E p 2 1 dt
where

∆(x, ξ) = ∂ 2 p 0 ∂x 2 ∂ 2 p 0 ∂ξ 2 - ∂ 2 p 0 ∂x ∂ξ 2
We recall that S 3 (E) = 0. In constrast with the convention of [CdV], our integrals are oriented integrals, t denoting the variable in Hamilton's equations. This explains why, in our expressions for S 2 (E), derivatives with respect to E (the conjugate variable to t) of such integrals have the opposite sign to the corresponding ones in [CdV]. See also [IfaM'haRo].

There are lots of ways to derive BS: the method of matching of WKB solutions [BenOrz], known also as Liouville-Green method [Ol], which has received many improvements, see e.g. [Ya]; the method of the monodromy operator, see [HeRo] and references therein; the method of quantization deformation based on Functional Calculus and Trace Formulas [Li], [CdV1], [CaGra-SazLiReiRios], [Gra-Saz], [Arg]. Note that the method of quantization deformation already assumes BS, it gives only a very convenient way to derive it. In the real analytic case, BS rule, and also tunneling expansions, can be obtained using the so-called "exact WKB method" see e.g. [Fe], [DePh], [DeDiPh] when P is Schrödinger operator.

Here we present still another derivation of BS, based on the construction of a Hermitian vector bundle of quasi-modes as in [Sj2], [HeSj]. Let K N h (E) be the microlocal kernel of P -E of order N , i.e. the space of microlocal solutions of (P -E)u = O(h N +1 ) along the covering of γ E (see Appendix for a precise definition). The problem is to find the set of E = E(h) such that K N h (E) contains a global section, i.e. to construct a sequence of quasi-modes (QM) (u n (h), E n (h)) of a given order N (practically N = 2). As usual we denote by K h (E) the microlocal kernel of P -E mod O(h ∞ ) ; since the distinction between K N h (E) and K h (E) plays no important rôle here, we shall be content to write K h (E).

Actually the method of [Sj2], [HeSj] was elaborated in case of a separatrix, and extends easily to mode crossing in Born-Oppenheimer type Hamiltonians as in [B], [Ro], but somewhat surprisingly it turns out to be harder to set up in case of a regular orbit, due to "translation invariance" of the Hamiltonian flow. In the present scalar case, when carried to second order, our method is also more intricated than [Li], [CdV1] and its refinements [Gra-Saz] for higher order N ; nevertheless it shows most useful for matrix valued operators with double characteristics such as Bogoliubov-de Gennes Hamiltonian [DuGy] (see [BenIfaRo], [BenMhaRo]). This method also extends to the scalar case in higher dimensions for a periodic orbit (see [SjZw], [FaLoRo], [LoRo]).

The paper is organized as follows:

In Sect.1 we present the main idea of the argument on a simple example, and recall from [HeSj],

[Sj2] the definition of the microlocal Wronskian.

In Sect.2 we compute BS at lowest order in the special case of Schrödinger operator by means of microlocal Wronskian and Gram matrix.

In Sect.3 we proceed to more general constructions in the case of h-Pseudodifferential operator (0.1) so to recover BS at order 2.

In Sect.4 we use a simpler formalism based on action-angle variables, but which would not extend to systems such as Bogoliubov-de Gennes Hamiltonian.

In Sect.5, following [SjZw], we recall briefly the well-posedness of Grushin problem, which shows in particular that there is no other spectrum in I than this given by BS.

At last, the Appendix accounts for a short introduction to microlocal and semi-classical Analysis used in the main text.

Main strategy of the proof.

The best algebraic and microlocal framework for computing quantization rules in the self-adjoint case, cast in the fundamental works [Sj2], [HeSj], is based on Fredholm theory, and the classical "positive commutator method" using conservation of some quantity called a "quantum flux".

a) A simple example

As a first warm-up, consider P = hD x acting on L 2 (S 1 ) with periodic boundary condition u

(x) = u(x + 2π). It is well-known that P has discrete spectrum E k (h) = kh, k ∈ Z, with eigenfunctions u k (x) = (2π) -1/2 e ikx = (2π) -1/2 e iE k (h)x/h . Thus BS quantization rule can be written as γ E ξ dx = 2πkh, where γ E = {x ∈ S 1 ; ξ = E}.
We are going to derive this result using the monodromy properties of the solutions of (hD x -E)u = 0. For notational convenience, we change energy variable E into z. Solving for (Pz)u(x) = 0, we get two solutions with the same expression but defined on different charts

(1.1) u a (x) = e izx/h , -π < x < π, u a ′ (x) = e izx/h , 0 < x < 2π
indexed by angles a = 0 and a ′ = π on S 1 . In the following we take advantage of the fact that these functions differ but when z belongs to the spectrum of P .

Let also χ a ∈ C ∞ 0 (S 1 ) be equal to 1 near a, χ a ′ = 1χ a . We set F a ± = i h [P, χ a ] ± u a , where ± denotes the part of the commutator supported in the half circles 0 < x < π and -π < x < 0 mod 2π.

Similarly F a ′ ± = i h [P, χ a ′ ] ± u a ′ . We compute (u a |F a + ) = u a |(χ a ) ′ u a = π 0 (χ a ) ′ (x) dx = χ a (π) -χ a (0) = -1
Similarly (u a |F a -) = 1, and also replacing a by a ′ so that

(1.2) (u a |F a + -F a -) = -2, (u a ′ |F a ′ + -F a ′ -) = 2
We evaluate next the crossed terms (u

a ′ |F a + -F a -) and (u a |F a ′ + -F a ′ -). Since u a ′ (x) = u a (x) = e izx/h
on the upper-half circle (once embedded into the complex plane), and u a (x) = e izx/h , u a ′ (x) = e iz(x+2π)/h on the lower-half circle we have

(u a ′ |F a + -F a -) = π 0 e izx/h (χ a ) ′ e -izx/h dx - 0 -π e iz(x+2π)/h (χ a ) ′ e -izx/h dx We argue similarly for (u a |F a ′ + -F a ′ -), using also that (χ a ′ ) ′ = -(χ a ) ′ . So we have (1.3) (u a ′ |F a + -F a -) = -1 -e 2iπz/h , (u a |F a ′ + -F a ′ -) = 1 + e -2iπz/h
It is convenient to view F a + -F a -and F a ′ + -F a ′ -as belonging to co-kernel of Pz in the sense they are not annihilated by Pz. So we form Gram matrix

(1.4) G (a,a ′ ) (z) = (u a |F a + -F a -) (u a ′ |F a + -F a -) (u a |F a ′ + -F a ′ -) (u a ′ |F a ′ + -F a ′ -)
and an elementary computation using (1.2) and (1.3) shows that det G (a,a ′ ) (z) = -4 sin 2 (πz/h) so the condition that u a coincides with u a ′ is precisely that z = kh, with k ∈ Z.

Next we investigate Fredholm properties of P as in [SjZw], recovering the fact that hZ is the only spectrum of P .

Notice that (1.4) is not affected when multiplying u a ′ by a phase factor, so we can replace u a ′ by e -izπ/h u a ′ . Starting from the point a = 0 we associate with u a the multiplication operator

v + → I a (z)v + = u a (x)v + on C, i.e. Poisson operator with "Cauchy data" u(0) = v + ∈ C. Define the "trace operator" R + (z)u = u(0).
Similarly multiplication by u a ′ defines Poisson operator I a ′ (z)v + = u a ′ (x)v + , which has the same "Cauchy data" v + at a ′ = π as I a (z) at a = 0.

Consider the multiplication operators

E + (z) = χ a I a (z) + (1 -χ a )e iπz/h I a ′ (z), R -(z) = i h [P, χ a ] -I a ′ (z), E -+ (z) = 2h sin(πz/h) We claim that (1.5) (P -z)E + (z) + R -(z)E -+ (z) = 0
Namely as before (but after we have replaced u a ′ by e -izπ/h u a ′ ) evaluating on 0 < x < π, we have I a (z) = e ixz/h , I a ′ (z) = e -iπz/h e ixz/h , while evaluating on -π < x < 0, I a (z) = e ixz/h , I a ′ (z) = e -iπz/h e i(x+2π)z/h . Now (Pz)E + (z) = [P, χ a ] I a (z)e iπz/h I a ′ (z) vanishes on 0 < x < π, while is precisely equal to 2h sin(πz/h) i h [P, χ a ]I a ′ (z) on -π < x < 0. So (1.5) follows. Hence the Grushin problem

(1.6) P(z; h) u u - = P -z R -(z) R + (z) 0 u u - = v v + with v = 0 has a solution u = E + (z)v + , u -= E -+ (z)v + , with E -+ (z) the effective Hamiltonian.
Following [SjZw] one can show that with this choice of R ± (z), problem (1.6) is well posed, P(z) is invertible, and

(1.7)

P(z) -1 = E(z) E + (z) E -(z) E -+ (z) with (1.8) (P -z) -1 = E(z) -E + (z)E -+ (z) -1 E -(z)
Hence z is an eigenvalue of P iff E -+ (z) = 0, which gives the spectrum z = kh as expected.

These Fredholm properties have been further generalized to a periodic orbit in higher dimensions in several ways [SjZw], [NoSjZw], [FaLoRo] where E -+ (z) is defined by means of the monodromy operator as E -+ (z) = Id -M (z) (in this example M (z) = e 2iπz/h ). In fact our argument here differs essentially from the corresponding one in [SjZw] by the choice of cutt-off χ a . We have considered functions on S 1 , but in Sect.4, we work on the covering of S 1 instead, using a single Poisson operator.

b) The microlocal Wronskian.

We now consider Bohr-Sommerfeld on the real line. Contrary to the periodic case that we have just investigated, where Maslov index is m = 0, we get in general m = 2 for BS on the real line, as is the case for the harmonic oscillator P = (hD x ) 2 + x 2 on L 2 (R). Otherwise, the argument is pretty much the same.

Bohr-Sommerfeld quantization rules result in constructing quasi-modes by WKB approximation along a closed Lagrangian manifold Λ E ⊂ {p 0 = E}, i.e. a periodic orbit of Hamilton vector field H p 0

with energy E. This can be done locally according to the rank of the projection Λ E → R x .

Thus the set K h (E) of asymptotic solutions to (P -E)u = 0 along (the covering of) Λ E can be considered as a bundle over R with a compact base, corresponding to the "classically allowed region" at energy E. The sequence of eigenvalues E = E n (h) is determined by the condition that the resulting quasi-mode, gluing together asymptotic solutions from different coordinates patches along Λ E , be single-valued, i.e. K h (E) have trivial holonomy.

Assuming Here the notation χ a holds for χ a (x, hD x ) as in (0.3), and we shall write P (x, hD x ) (spatial representation) as well as P (-hD ξ , ξ) (Fourier representation). Recall that unitary h-Fourier transform for a semiclassical distribution u(x; h) is given by u(ξ; h) = (2πh) -1/2 e -ixξ/h u(x; h) dx (see Appendix for a review of semi-classical asymptotics).

Definition 1.1: Let P be self-adjoint, and u a , v a ∈ K h (E) be supported microlocally on Λ E . We call

(1.10) W a ρ (u a , v a ) = i h [P, χ a ] ρ u a |v a
the microlocal Wronskian of (u a , v a ) in ω a ρ . Here i h [P, χ a ] ρ denotes the part of the commutator supported on ω a ρ . To understand that terminology, let P = -h 2 ∆+V , x E = 0 and change χ to Heaviside unit stepfunction χ(x), depending on x alone. Then in distributional sense, we have i h [P, χ] = -ihδ ′ + 2δhD x , where δ denotes the Dirac measure at 0, and δ ′ its derivative, so that i h [P, χ]u|u = -ih u ′ (0)u(0)u(0)u ′ (0) is the usual Wronskian of (u, u).

Proposition 1.2: Let u a , v a ∈ K h (E) be as above, and denote by u the h-Fourier (unitary) transform of u.

Then i h [P, χ a ]u a |v a = i h [P, χ a ] u a | v a = 0 (1.11) i h [P, χ a ] + u a |v a = - i h [P, χ a ] -u a |v a (1.12) all equalities being understood mod O(h ∞ ), (resp. O(h N +1 )) when considering u a , v a ∈ K N h (E) instead. Moreover, W a ρ (u a , v a ) does not depend mod O(h ∞ ) (resp. O(h N +1
)) on the choice of χ a as above.

Proof: Since u a , v a ∈ K h (E) are distributions in L 2 , the equality (1.11) follows from Plancherel formula and the regularity of microlocal solutions in L 2 , p + i being elliptic. If a is not a focal point, u a , v a are smooth WKB solutions near a, so we can expand the commutator in w = i h [P, χ a ]u a |v a and use that P is self-adjoint to show that w = O(h ∞ ). If a is a focal point, u a , v a are smooth WKB solutions in Fourier representation, so again w = O(h ∞ ). Then (1.12) follows from Definition 1.1. ♣

We can find a linear combination of W a ± , (depending on a) which defines a sesquilinear form on K h (E), so that this Hermitean form makes K h (E) a metric bundle, endowed with the gauge group U (1). This linear combinaison is prescribed as the construction of Maslov index : namely we take W a (u a , u a ) = W a + (u a , u a ) -W a -(u a , u a ) > 0 when the critical point a of π Λ E is traversed in the -ξ direction to the right of the fiber (or equivalently W a (u a , u a ) = -W a + (u a , u a ) + W a -(u a , u a ) > 0 while traversing a in the +ξ direction to the left of the fiber). Otherwise, just exchange the signs. When Λ E is a convex curve, there are only 2 focal points. In general there may be many focal points a, but each jump of Maslov index is compensated at the next focal point while traversing to the other side of the fiber (Maslov index is computed mod 4), see [START_REF] Bates | Lectures on the geometry of quantization[END_REF]Example 4.13].

As before our method consists in constructing Gram matrix of a generating system of K h (E) in a suitable dual basis; its determinant vanishes precisely at the eigenvalues E = E n (h).

Note that when energy surface p 0 = E is singular, and Λ E is a separatrix ("figure eight", or homoclinic case), equality (1.12) does not hold near the "branching point", see [Sj2] and its generalization to multi-dimensional case [BoFuRaZe].

Bohr-Sommerfeld quantization rules in the case of a Schrödinger operator.

As a second warm-up, we derive the well known BS quantization rule using microlocal Wronskians in case of a potential well, i.e. Λ E has only 2 focal points. Consider the spectrum of Schrödinger

operator P (x, hD x ) = (hD x ) 2 + V (x) near the energy level E 0 < lim inf |x|→∞ V (x), when {V ≤ E} = [x ′ E , x E ] and x ′ E , x E are simple turning points, V (x ′ E ) = V (x E ) = E, V ′ (x ′ E ) < 0, V ′ (x E ) > 0.
For a survey of WKB theory, see e.g. [Dui], [BaWe] or [CdV]. It is convenient to start the construction from the focal points a or a ′ . We identify a focal point a = a E = (x E , 0) with its projection x E . We know that microlocal solutions u of (P -E)u = 0 in a (punctured) neighborhood of a are of the form

(2.1) u a (x, h) = C √ 2 e iπ/4 (E -V ) -1/4 e iS(a,x)/h + e -iπ/4 (E -V ) -1/4 e -iS(a,x)/h + O(h) , C ∈ C
where S(y, x) =

x y ξ + (t) dt, and ξ + (t) is the positive root of ξ 2 + V (t) = E. In the same way, the microlocal solutions of (P -E)u = 0 in a (punctured) neighborhood of a ′ have the form

(2.2) u a ′ (x, h) = C ′ √ 2 e -iπ/4 (E -V ) -1/4 e iS(a ′ ,x)/h + e iπ/4 (E -V ) -1/4 e -iS(a ′ ,x)/h + O(h) , C ′ ∈ C
These expressions result in computing by the method of stationary phase the oscillatory integral that gives the solution of (P (-hD ξ , ξ) -E) u = 0 in Fourier representation. The change of phase factor e ±iπ/4 accounts for Maslov index. For later purposes, we recall here from [START_REF] Hörmander | The Analysis of Partial Differential Operators[END_REF]Thm 7.7.5] that if

f : R d → C, with Im f ≥ 0 has a non-degenerate critical point at x 0 , then (2.3) R d e i h f (x) u(x) dx ∼ e i h f (x 0 ) det( f ′′ (x 0 ) 2iπh ) -1/2 j h j L j (u)(x 0 )
where L j are linear forms, L 0 u(x 0 ) = u(x 0 ), and

(2.4) L 1 u(x 0 ) = 2 n=0 2 -(n+1) in!(n + 1)! (f ′′ (x 0 )) -1 D x , D x n+1 (Φ x 0 ) n u)(x 0 ) where Φ x 0 (x) = f (x) -f (x 0 ) -1 2 f ′′ (x 0 )(x -x 0 ),
xx 0 vanishes of order 3 at x 0 . For the sake of simplicity, we omit henceforth O(h) terms, but the computations below extend to all order in h (practically, at least for N = 2), thus giving the asymptotics of BS. This will be elaborated in Section 3.

The semi-classical distributions u a , u a ′ span the microlocal kernel K h of P -E in (x, ξ) ∈ ]a ′ , a[×R ; they are normalized using microlocal Wronskians as follows.

Let χ a ∈ C ∞ 0 (R 2 ) as in the Introduction be a smooth cut-off equal to 1 near a. Without loss of generality, we can take

χ a (x, ξ) = χ a 1 (x)χ 2 (ξ), so that χ 2 ≡ 1 on small neighborhoods ω a ± , of supp[P, χ a ] ∩ {ξ 2 + V = E} in ±ξ > 0. We define χ a ′ similarly. Since i h [P, χ a ] = 2(χ a ) ′ (x)hD x - ih(χ a ) ′′ , by (2.1) and (2.2) we have, mod O(h): i h [P, χ a ]u a (x, h) = √ 2C(χ a 1 ) ′ (x) e iπ/4 (E -V ) 1/4 e iS(a,x)/h -e -iπ/4 (E -V ) 1/4 e -iS(a,x)/h i h [P, χ a ′ ]u a ′ (x, h) = √ 2C ′ (χ a ′ 1 ) ′ (x) e -iπ/4 (E -V ) 1/4 e iS(a ′ ,x)/h -e iπ/4 (E -V ) 1/4 e -iS(a ′ ,x)/h Let (2.5) F a ± (x, h) = i h [P, χ a ] ± u a (x, h) = ± √ 2C(χ a 1 ) ′ (x)e ±iπ/4 (E -V ) 1/4 e ±iS(a,x)/h
so that:

(u a |F a + -F a -) = |C| 2 e iπ/4 (E -V ) -1/4 e iS(a,x)/h |(χ a 1 ) ′ e iπ/4 (E -V ) 1/4 e iS(a,x)/h ) + |C| 2 (e -iπ/4 (E -V ) -1/4 e -iS(a,x)/h |(χ a 1 ) ′ e -iπ/4 (E -V ) 1/4 e -iS(a,x)/h ) + O(h) = |C| 2 ( a -∞ (χ a 1 ) ′ (x)dx + a -∞ (χ a 1 ) ′ (x)dx) + O(h) = 2|C| 2 + O(h)
(the mixed terms such as e iπ/4 (E -V ) -1/4 e iS(a,x)/h |(χ a 1 ) ′ e -iπ/4 (E -V ) 1/4 e -iS(a,x)/h ) are O(h ∞ ) because the phase is non stationary), thus u a is normalized mod O(h) if we choose C = 2 -1/2 . In the same way, with

(2.6) F a ′ ± (x, h) = i h [P, χ a ′ ] ± u a ′ (x, h) = ± √ 2C ′ (χ a ′ 1 ) ′ (x)e ∓iπ/4 (E -V ) 1/4 e ±iS(a ′ ,x)/h we get (u a ′ |F a ′ + -F a ′ -) = |C ′ | 2 ( ∞ a ′ (χ a ′ 1 ) ′ (x)dx + ∞ a ′ (χ a ′ 1 ) ′ (x)dx) + O(h) = -2|C ′ | 2 + O(h)
and we choose again C ′ = C which normalizes u a ′ mod O(h). Normalization carries to higher order, as is shown in Sect.3 for a more general Hamiltonian.

So there is a natural duality product between K h (E) and the span of functions F a + -F a -and

F a ′ + -F a ′ -in L 2 . As in [Sj2]
, [HeSj] we can show that this space is microlocally transverse to Im(P -E) on (x, ξ) ∈]a ′ , a[×R, and thus identifies with the microlocal co-kernel

K * h (E) of P -E; in general dim K h (E) = dim K * h (E) = 2, unless E is an eigenvalue, in which case dim K h = dim K * h = 1 (showing that P -E is of index 0 when Fredholm, which is indeed the case. )
Microlocal solutions u a and u a ′ extend as smooth solutions on the whole interval ]a ′ , a[; we denote them by u 1 and u 2 . Since there are no other focal points between a and a ′ , they are expressed by the same formulae (which makes the analysis particularly simple) and satisfy :

(u 1 |F a + -F a -) = 1, (u 2 |F a ′ + -F a ′ -) = -1 Next we compute (still modulo O(h)) (u 1 |F a ′ + -F a ′ -) = 1 2 (e iπ/4 (E -V ) -1/4 e iS(a,x)/h |(χ a ′ 1 ) ′ e -iπ/4 (E -V ) 1/4 e iS(a ′ ,x)/h ) + 1 2 (e -iπ/4 (E -V ) -1/4 e -iS(a,x)/h |(χ a ′ 1 ) ′ e iπ/4 (E -V ) 1/4 e -iS(a ′ ,x)/h ) = i 2 e -iS(a ′ ,a)/h ∞ a ′ (χ a ′ 1 ) ′ (x)dx - i 2 e iS(a ′ ,a)/h ∞ a ′ (χ a ′ 1 ) ′ (x)dx = -sin(S(a ′ , a)/h) (taking again into account that the mixed terms are O(h ∞ )). Similarly (u 2 |F a + -F a -) = sin(S(a ′ , a)/h). Now we define Gram matrix (2.7) G (a,a ′ ) (E) = (u 1 |F a + -F a -) (u 2 |F a + -F a -) (u 1 |F a ′ + -F a ′ -) (u 2 |F a ′ + -F a ′ -)
whose determinant -1 + sin 2 (S(a ′ , a)/h) =cos 2 (S(a ′ , a)/h) vanishes precisely on eigenvalues of P in I, so we recover the well known BS quantization condition

(2.8) ξ(x) dx = 2 a a ′ (E -V ) 1/2 dx = 2πh(k + 1 2 ) + O(h)
and det G (a,a ′ ) (E) is nothing but Jost function which is computed e.g. in [DePh], [DeDiPh] by another method.

The general case

By the discussion after Proposition 1.1, it clearly suffices to consider the case when γ E contains only 2 focal points which contribute to Maslov index. We shall content throughout to BS mod O(h 2 ).

; h) = a 0 (x, ξ) + ha 1 (x, ξ) + • • •, such that hD ξ e i(xξ+ψ(ξ))/h a(x, ξ; h) = P (x, D x ; h) e i(xξ+ψ(ξ))/h b(ξ; h)
Expanding the RHS by stationary phase (2.3), we find

hD ξ e i(xξ+ψ(ξ))/h a(x, ξ; h) = e i(xξ+ψ(ξ))/h b(ξ; h) p 0 (x, ξ) -E + h p 1 (x, ξ) + h 2 p 2 (x, ξ) + O(h 3 )
p 0 being the principal symbol of P ,

p 1 (x, ξ) = p 1 (x, ξ) + 1 2i ∂ 2 p 0 ∂x∂ξ (x, ξ), p 2 (x, ξ) = p 2 (x, ξ) + 1 2i ∂ 2 p 1 ∂x∂ξ (x, ξ) - 1 8 ∂ 4 p 0 ∂x 2 ∂ξ 2 (x, ξ)
Collecting the coefficients of ascending powers of h, we get

(p 0 -E)b 0 = (x + ψ ′ (ξ))a 0 (3.1) 0 (p 0 -E)b 1 + p 1 b 0 = (x + ψ ′ (ξ))a 1 + 1 i ∂a 0 ∂ξ (3.1) 1 (p 0 -E)b 2 + p 1 b 1 + p 2 b 0 = (x + ψ ′ (ξ))a 2 + 1 i ∂a 1 ∂ξ (3.1) 2 and so on. Define λ(x, ξ) by p 0 (x, ξ) -E = λ(x, ξ)(x + ψ ′ (ξ)), we have (3.2) λ(-ψ ′ (ξ), ξ) = ∂ x p 0 (-ψ ′ (ξ), ξ) = α(ξ)
This gives a 0 (x, ξ) = λ(x, ξ)b 0 (ξ) for (3.1) 0 . We look for b 0 by noticing that (3.1) 1 is solvable iff

( p 1 b 0 )| x=-ψ ′ (ξ) = 1 i ∂a 0 ∂ξ | x=-ψ ′ (ξ)
which yields the first order ODE L(ξ, D ξ )b 0 = 0, with

L(ξ, D ξ ) = α(ξ)D ξ + 1 2i α ′ (ξ) -p 1 (-ψ ′ (ξ), ξ). We find b 0 (ξ) = C 0 |α(ξ)| -1/2 e i p 1
α with an arbitrary constant C 0 . This gives in turn

(3.3) a 1 (x, ξ) = λ(x, ξ)b 1 (ξ) + λ 0 (x, ξ) with λ 0 (x, ξ) = b 0 (ξ) p 1 + i ∂a 0 ∂ξ x + ∂ ξ ψ
which is smooth near a E . At the next step, we look for b 1 by noticing that (3.1) 2 is solvable iff

( p 1 b 1 + p 2 b 0 )| x=-ψ ′ (ξ) = 1 i ∂a 1 ∂ξ | x=-ψ ′ (ξ) Differentiating (3.3) gives L(ξ, D ξ )b 1 = p 2 b 0 + i∂ ξ λ 0 | x=-ψ ′ (ξ)
, which we solve for b 1 . We eventually

get, mod O(h 2 ) (3.4) u a (ξ; h) = (C 0 + hC 1 + hD 1 (ξ))|α(ξ)| -1/2 exp i h ψ(ξ) + h p 1 (-ψ ′ (ζ), ζ) α(ζ) dζ]
where we have set (for ξ close enough to ξ E so that α(ξ) = 0)

(3.5) D 1 (ξ) = sgn(α(ξ E )) ξ ξ E exp[-i ζ ξ E p 1 α ] i p 2 b 0 -∂ ξ λ 0 | x=-ψ ′ (ζ) |α(ζ)| -1/2 dζ
The integration constants C 0 , C 1 will be determined by normalizing the microlocal Wronskians as follows. We postpone to Sect.3.c the proof of this Proposition making us of the spatial representation of u a .

Proposition 3.1: With the hypotheses above, the microlocal Wronskian near a focal point a E is given by

W a (u a , u a ) = W a + (u a , u a ) -W a -(u a , u a ) = 2 sgn(α(ξ E )) |C 0 | 2 + h 2 Re(C 0 C 1 ) + |C 0 | 2 ∂ x p 1 ∂ x p 0 (ξ E ) + O(h 2 )
The condition that u a be normalized mod O(h 2 ) (once we have chosen C 0 to be real), is then

(3.6) C 1 (E) = - 1 2 C 0 ∂ x p 1 ∂ x p 0 (a E ) so that now W a (u a , u a ) = 2 sgn(α(ξ E ))C 2 0 1 + O(h 2
) . We say that u a is well-normalized mod O(h 2 ). This can be formalized by considering {a E } as a Poincaré section (see Sect.4), and Poisson operator the operator that assigns, in a unique way, to the initial condition C 0 on {a E } the well-normalized (forward) solution u a to (P -E)u a = 0: namely, C 1 (E) and D 1 (ξ), hence also u a , depend linearly on C 0 . Using the approximation

C 0 + hC 1 (E) + hD 1 (ξ) = C 0 + hC 1 (E) + h Re(D 1 (ξ)) exp ih C 0 Im(D 1 (ξ)) + O(h 2 )
the normalized WKB solution near a E now writes, by (3.4)

(3.7) u a (ξ; h) = C 0 + hC 1 (E) + h Re(D 1 (ξ)) |α(ξ)| -1 2 exp i S(ξ, ξ E ; h)/h (1 + O(h 2 ))
with the h-dependent phase function

S(ξ, ξ E ; h) = ψ(ξ) + h ξ ξ E p 1 (-ψ ′ (ζ), ζ) α(ζ) dζ + h 2 C 0 Im(D 1 (ξ))
The modulus of u a (ξ; h) can further be simplified using (3.6) and formula (3.10) below:

C 0 +hC 1 (E)+h Re(D 1 (ξ)) = C 0 1- h 2 ∂ x p 1 ∂ x p 0 | x=-ψ ′ (ξ) = C 0 exp h∂ x p 1 ∂ x p 0 | x=-ψ ′ (ξ) -1/2 +O(h 2 )
which altogether, recalling α(ξ) = ∂ x p 0 (-ψ ′ (ξ), ξ) near ξ E (and assuming α(ξ E ) > 0 to fix the ideas), gives

(3.8) u a (ξ; h) = 1 √ 2 (∂ x p 0 ) exp h∂ x p 1 ∂ x p 0 -1/2 exp i S(ξ, ξ E ; h)/h (1 + O(h 2 ))
b) The homology class of the generalized action: Fourier representation.

Here we identify the various terms in (3.8), which are responsible for the holonomy of u a . First on γ E (i.e. Λ E ) we have ψ(ξ) = -x dξ + Const., and ϕ(x) = ξ dx + Const. By Hamilton equations

ξ(t) = -∂ x p 0 (x(t), ξ(t)), ẋ(t) = ∂ ξ p 0 (x(t), ξ(t)) so p 1 ∂ x p 0 dξ = - p 1 ∂ ξ p 0 dx = -γ E p 1 dt.
The form p 1 dt is called the subprincipal 1-form. Next we consider D 1 (ξ) as the integral over γ E of the 1-form, defined near a in Fourier representation as

(3.9) Ω 1 = T 1 dξ = sgn(α(ξ)) i p 2 b 0 -∂ ξ λ 0 |α| -1/2 e -i p 1 α dξ
Since γ E is Lagrangian, Ω 1 is a closed form that we are going to compute modulo exact forms. Using integration by parts, the integral of Ω 1 (ξ) in Fourier representation simplifies to

√ 2 Re D 1 (ξ) = - 1 2 ∂ x p 1 ∂ x p 0 ] ξ ξ E = - 1 2 ∂ x p 1 ∂ x p 0 (ξ) - C 1 (E) C 0 (3.10) √ 2 ImD 1 (ζ) = ξ ξ E T 1 (ζ) dζ + ψ ′′ 6α ∂ 3 x p 0 + α ′ 4α 2 ∂ 2 x p 0 ξ ξ E (3.11) T 1 = 1 α p 2 - 1 8 ∂ 2 x ∂ 2 ξ p 0 + ψ ′′ 12 ∂ 3 x ∂ ξ p 0 + (ψ ′′ ) 2 24 ∂ 4 x p 0 + 1 8 (α ′ ) 2 α 3 ∂ 2 x p 0 + 1 6 ψ ′′ α ′ α 2 ∂ 3 x p 0 - p 1 α 2 ∂ x p 1 - p 1 2α ∂ 2
x p 0 (3.12) There follows: Lemma 3.2: Modulo the integral of an exact form in A, with T 1 as in (3.12) we have:

(3.13) Re D 1 (ξ) ≡ 0 √ 2 Im D 1 (ξ) ≡ ξ ξ E T 1 (ζ) dζ
Passing from Fourier to spatial representation, we can carry the integration in x-variable between the focal points a E and a ′ E , and in ξ-variable again near a ′ E . Since γ E is smoothly embedded, the microlocal solution u a extends uniquely along γ E .

If f (x, ξ), g(x, ξ) are any smooth functions on A we set Ω(x, ξ) = f (x, ξ) dx + g(x, ξ) dξ. By Stokes formula

γ E Ω(x, ξ) = p 0 ≤E (∂ x g -∂ ξ f ) dx ∧ dξ
where, following [CdV], we have extended p 0 in the disk bounded by A -so that it coincides with a harmonic oscillator in a neighborhood of a point inside, say p 0 (0, 0) = 0. Making the symplectic change of coordinates (x, ξ) → (t, E) in T * R:

p 0 ≤E (∂ x g -∂ ξ f ) dx ∧ dξ = E 0 T (E ′ ) 0 (∂ x g -∂ ξ f ) dt ∧ dE ′
where T (E ′ ) is the period of the flow of Hamilton vector field H p 0 at energy E ′ (T (E ′ ) being a constant near (0,0)). Taking derivative with respect to E, we find

(3.14) d dE γ E Ω(x, ξ) = T (E) 0 (∂ x g -∂ ξ f ) dt
We compute ξ ξ E T 1 (ζ) dζ with T 1 as in (3.12), and start to simplify J 1 = ω 1 , with ω 1 the last term on the RHS of (3.12). Let g 1 (x, ξ) = p 2 1 (x,ξ) ξ) , by (3.14) we get (3.15)

∂ x p 0 (x,
J 1 = 1 2 γ E ∂ x g 1 (x, ξ) ∂ x p 0 (x, ξ) dξ = - 1 2 T (E) 0 ∂ x g 1 (x(t), ξ(t)) dt = - 1 2 d dE γ E g 1 (x, ξ) dξ = - 1 2 d dE γ E p 2 1 (x, ξ) ∂ x p 0 (x, ξ) dξ = 1 2 d dE T (E) 0 p 2 1 (x(t), ξ(t)) dt
which is the contribution of p 1 to the second term S 2 of generalized action in [START_REF] De Verdière | 1. Bohr Sommerfeld rules to all orders[END_REF]Thm2]. Here

T (E) is the period on γ E . We also have

(3.16) ξ ξ E 1 α(ξ) p 2 (-ψ ′ (ξ), ξ) dξ = γ E p 2 (x, ξ) ∂ x p 0 (x, ξ) dξ = - T (E) 0 p 2 (x(t), ξ(t)) dt
To compute T 1 modulo exact forms we are left to simplify in (3.12) the expression ξ) , where we have set according to [CdV] ∆

J 2 = ξ ξ E 1 α - 1 8 ∂ 4 p 0 ∂x 2 ∂ξ 2 + ψ ′′ 12 ∂ 4 p 0 ∂x 3 ∂ξ + (ψ ′′ ) 2 24 ∂ 4 p 0 ∂x 4 dζ + 1 8 ξ ξ E (α ′ ) 2 α 3 ∂ 2 p 0 ∂x 2 dζ + 1 6 ξ ξ E ψ ′′ α ′ α 2 ∂ 3 p 0 ∂x 3 dζ Let g 0 (x, ξ) = ∆(x,ξ) ∂ x p 0 (x,
(x, ξ) = ∂ 2 p 0 ∂x 2 ∂ 2 p 0 ∂ξ 2 - ∂ 2 p 0 ∂x ∂ξ 2
Taking second derivative of eikonal equation p 0 (-ψ ′ (ξ), ξ) = E, we get

(∂ x g 0 )(-ψ ′ (ξ), ξ) α(ξ) = ψ ′′′ α ∂ 3 p 0 ∂x 3 + 2ψ ′′ α ′ α 2 ∂ 3 p 0 ∂x 3 + α ′′ α 2 ∂ 2 p 0 ∂x 2 -2 α ′ α 2 ∂ 3 p 0 ∂x 2 ∂ξ + (α ′ ) 2 α 3 ∂ 2 p 0 ∂x 2
Integration by parts of the first and third term on the RHS gives altogether

ξ ξ E (∂ x g 0 )(-ψ ′ (ξ), ξ) α(ξ) dξ = -3 ξ ξ E 1 α ∂ 4 p 0 ∂x 2 ∂ξ 2 dζ + 2 ξ ξ E ψ ′′ α ∂ 4 p 0 ∂x 3 ∂ξ dζ + ξ ξ E (ψ ′′ ) 2 α ∂ 4 p 0 ∂x 4 dζ + 3 ξ ξ E (α ′ ) 2 α 3 ∂ 2 p 0 ∂x 2 dζ + 4 ξ ξ E ψ ′′ α ′ α 2 ∂ 3 p 0 ∂x 3 dζ + ψ ′′ α ∂ 3 p 0 ∂x 3 ξ ξ(E) + α ′ α 2 ∂ 2 p 0 ∂x 2 ξ ξ E + 3 1 α ∂ 3 p 0 ∂x 2 ∂ξ ξ ξ E
and modulo the integral of an exact form in A

J 2 ≡ 1 24 ξ ξ E (∂ x g 0 )(-ψ ′ (ζ), ζ) α(ζ) dζ = - 1 24 T (E) 0 ∂ x g 0 (x(t), ξ(t)) dt = - 1 24 d dE γ E g 0 (x, ξ) dξ = - 1 24 d dE γ E ∆(x, ξ) ∂ x p 0 (x, ξ) dξ = 1 24 d dE T (E) 0 ∆(x(t), ξ(t)) dt
Using these expressions, we recover the well known action integrals (see e.g. [CdV]):

Proposition 3.3: Let Γ dt be the restriction to γ E of the 1-form

ω 0 (x, ξ) = (∂ 2 x p 0 )(∂ ξ p 0 ) -(∂ x ∂ ξ p 0 )(∂ x p 0 ) dx + (∂ ξ p 0 )(∂ ξ ∂ x p 0 ) -(∂ 2 ξ p 0 )(∂ x p 0 ) dξ We have Re γ E Ω 1 = 0, whereas Im γ E Ω 1 = 1 48 d dE 2 γ E Γ dt - γ E p 2 dt - 1 2 d dE γ E p 2 1 dt c) Well normalized QM mod O(h 2 )
in the spatial representation.

The next task consists in extending the solutions away from a E in the spatial representation.

First we expand u a (x) = (2πh) -1/2 e ixξ/h u a (ξ; h) dξ = (2πh) -1/2 e i(xξ+ψ(ξ))/h b(ξ; h) dξ near x E by stationary phase (2.4) mod O(h 2 ), selecting the 2 critical points ξ ± (x) near x E . The phase functions take the form ϕ ± (x) = xξ ± (x) + ψ(ξ ± (x)).

Lemma 3.4: In a neighborhood of the focal point a E and for x < x E , the microlocal solution of (P (x, hD x ; h) -E)u(x; h) = 0 is given by (with ±∂ ξ p 0 (x, ξ ± (x)) > 0)

(3.17)

u a (x; h) = 1 √ 2 ± e ±iπ/4 ±∂ ξ p 0 (x, ξ ± (x)) -1/2 exp i h ϕ ± (x) -h x x E p 1 (y, ξ ± (y)) ∂ ξ p 0 (y, ξ ± (y) dy 1 + h √ 2 C 1 + D 1 (ξ ± (x)) + hD 2 (ξ ± (x)) + O(h 2 ) with (3.18) D 2 (ξ) = - 1 2i (ψ ′′ (ξ)) -1 b ′′ 0 (ξ) b 0 (ξ) + 1 8i (ψ ′′ (ξ)) -2 ψ (4) (ξ) + 4ψ (3) (ξ) b ′ 0 (ξ) b 0 (ξ) - 5 24i (ψ ′′ (ξ)) -3 (ψ (3) (ξ)) 2
The quantity

√ 2(C 1 + D 1 (ξ))
) has been computed before; with the particular choice of C 1 = C 1 (E) in (3.6) we have:

√ 2(C 1 + D 1 (ξ))) = - 1 2 ∂ x p 1 ∂ x p 0 (-ψ ′ (ξ), ξ) + i √ 2 Im D 1 (ξ) Moreover b ′ 0 (ξ) b 0 (ξ) = - α ′ (ξ) 2 α(ξ) + ip 1 (-ψ ′ (ξ), ξ) α(ξ) b ′′ 0 (ξ) b 0 (ξ) = - α ′ (ξ) 2 α(ξ) + i p 1 (-ψ ′ (ξ), ξ) α(ξ) 2 + d dξ - α ′ (ξ) 2 α(ξ) + ip 1 (-ψ ′ (ξ), ξ) α(ξ)
First, we observe that D 2 (ξ ± (x)) does not contribute to the homology class of the semi-classical forms defining the action, since it contains no integral. Thus the phase in (3.17) can be replaced, mod O(h 3 ) by

(3.19) S ± (x E , x; h) = x E ξ E + x x E ξ ± (y) dy -h x x E p 1 (y, ξ ρ (y)) ∂ ξ p 0 (y, ξ ρ (y) dy + √ 2h 2 Im D 1 (ξ ± (x))
with the residue of √ 2 Im D 1 (ξ ± (x)) , mod the integral of an exact form, computed as in Lemma 3.3.

Proof of Proposition 3.1. We proceed by using Proposition 1.2, and checking directly from (3.17) that normalization relations (u a |F a + ) = 1 2 and (u a |F a -) = -1 2 hold mod O(h 2 ) in the spatial representation, provided C 1 (E) takes the value (3.6). So let us compute F a ± (x) by stationary phase as in (3.17). In Fourier representation we have

(3.20) i h [P, χ a ] u(ξ) = (2πh) -1 e i -(ξ-η)y+ψ(η) /h c(y, ξ + η 2 ; h)(b 0 + hb 1 )(η) dy dη with Weyl symbol (3.21) c(x, ξ; h) ≡ c 0 (x, ξ) + hc 1 (x, ξ) = ∂ ξ p 0 (x, ξ) + h∂ ξ p 1 (x, ξ) χ ′ 1 (x) mod O(h 2 ) Let u ± x (y, η; h) = c( x + y 2 , η; h) ± ∂ ξ p 0 (y, ξ ± (y)) -1/2 exp -i y x E p 1 (z, ξ ± (z)) ∂ ξ p 0 (z, ξ ± (z)) dz × 1 + h √ 2 C 1 + D 1 (ξ ± (x) + hD 2 (ξ ± (x)) + O(h 2 )
with leading order term u (0,±) x (y, η). Applying stationary phase (2.3) gives

F a ± (x; h) = 1 √ 2 e ±iπ/4 e i h ϕ ± (x) u ± x x, ξ ± (x); h + h L 1 u (0,±) x (x, ξ ± (x)) + O(h 2 )
which simplifies as

F a ± (x; h) = ± 1 √ 2 e ±iπ/4 exp i h ϕ ± (x) -h x x E p 1 y, ξ ± (y) ∂ ξ p 0 y, ξ ± (y) dy ± ∂ ξ p 0 (x, ξ ± (x)) 1/2 1 + hZ(ξ ± (x)) + h c 1 (x, ξ ± (x)) c 0 (x, ξ ± (x)) + h 2s ± (x) θ ± (x) + s ′ ± (x) 2ic 0 (x, ξ ± (x)) χ ′ 1 (x) mod O(h 2 )
, where we recall c 0 , c 1 from (3.21). Here we have set

Z(ξ ± (x)) = √ 2 C 1 (E) + D 1 (ξ ± (x)) + D 2 (ξ ± (x) s ± (x) = ( ∂ 2 p 0 ∂ξ 2 )(x, ξ ± (x)) χ ′ 1 (x) = ω ± (x) χ ′ 1 (x) θ ± (x) = - 1 ψ ′′ (ξ ± (x)) α(ξ ± (x)) i p 1 x, ξ ± (x) - ψ ′′′ (ξ ± (x)) α(ξ ± (x)) + ψ ′′ (ξ ± (x)) α ′ (ξ ± (x)) 2 ψ ′′ (ξ ± (x))
and used the fact that

c 0 x, ξ ± (x) ± ∂ ξ p 0 (x, ξ ± (x)) -1/2 = ± ± ∂ ξ p 0 (x, ξ ± (x)) 1/2 χ ′ 1 (x) Since ∂ ξ p 0 (x, ξ ± (x)) = ψ ′′ (ξ ± (x)) α(ξ ± (x)) we obtain (3.22) F a ± (x; h) = ± 1 √ 2 e ±iπ/4 exp i h ϕ ± (x) -h x x E p 1 (y, ξ ± (y)) ∂ ξ p 0 (y, ξ ± (y)) dy ± ∂ ξ p 0 (x, ξ ± (x)) 1/2 χ ′ 1 (x) 1 + h Re Z(ξ ± (x)) + h ∂ ξ p 1 (x, ξ ± (x)) ∂ ξ p 0 (x, ξ ± (x)) -ih ω ± (x) θ ± (x) ∂ ξ p 0 (x, ξ ± (x)) - ih 2 d dx ω ± (x) χ ′ 1 (x) ∂ ξ p 0 (x, ξ ± (x)) χ ′ 1 (x) + O(h 2 )
Taking the scalar product with u a ± gives in particular (3.23)

(u a + |F a + ) = 1 2 +∞ x E χ ′ 1 (x) dx+ h 2 +∞ x E 2 Re Z(ξ ± (x)) + ∂ ξ p 1 (x, ξ + (x)) ψ ′′ (ξ + (x)) α(ξ + (x)) + iω + (x)θ + (x)ψ ′′ (ξ + (x)) α(ξ + (x)) χ ′ 1 (x) dx + ih 4 +∞ x E 1 ψ ′′ (ξ + (x)) α(ξ + (x)) d dx ω + (x) χ ′ 1 (x) dx + O(h 2 ) = 1 2 + h 2 K 1 + ih 4 K 2 + O(h 2 )
There remains to relate K 1 with K 2 . We have

(3.24) 2 Re Z(ξ ± (x)) + ∂ ξ p 1 (x, ξ + (x)) ψ ′′ (ξ + (x)) α(ξ + (x)) + i ω + (x) θ + (x) ψ ′′ (ξ + (x)) α(ξ + (x)) = ω + (x) ψ ′′ (ξ + (x))α(ξ + (x)) i θ + (x) + p 1 (x, ξ + (x)) ψ ′′ (ξ + (x)) α(ξ + (x)) = i ω + (x) 2 ψ ′′ (ξ + (x)) 3 α(ξ + (x)) 2 ψ ′′′ (ξ + (x)) α(ξ + (x)) + ψ ′′ (ξ + (x)) α ′ (ξ + (x))
whence

K 1 = i 2 +∞ x E ω + (x) ψ ′′ (ξ + (x)) 3 α(ξ + (x)) 2 ψ ′′′ (ξ + (x)) α(ξ + (x)) + ψ ′′ (ξ + (x)) α ′ (ξ + (x)) χ ′ 1 (x) dx
Here we have used that

2 Re Z(ξ + (x)) = -∂ x p 1 ∂ x p 0 (-ψ ′ (ξ), ξ) + 2 Re D 2 (ξ + (x)) ω + (x) = ψ ′′′ (ξ + (x)) α(ξ + (x)) + 2ψ ′′ (ξ + (x)) α ′ (ξ + (x)) + ψ ′′ (ξ + (x)) 2 ∂ 2 p 0 ∂x 2 )(x, ξ + (x))
On the other hand, integrating by parts gives

K 2 = ω + (x) χ ′ 1 (x) ψ ′′ (ξ + (x)) α(ξ + (x)) +∞ x E - +∞ x E d dx 1 ψ ′′ (ξ + (x)) α(ξ + (x)) ω + (x) χ ′ 1 (x) dx = - +∞ x E ω + (x) ψ ′′ (ξ + (x)) 3 α(ξ + (x)) 2 ψ ′′′ (ξ + (x)) α(ξ + (x)) + ψ ′′ (ξ + (x)) α ′ (ξ + (x)) χ ′ 1 (x) dx = 2iK 1 This shows (u a + |F a + ) = 1 2 + O(h 2
), and we argue similarly for (u a -|F a -), and Proposition 3.1 is proved.

Away from x E , we use standard WKB theory extending (3.17), with Ansatz (which we review in the Appendix)

(3.25) u a ± (x) = a ± (x; h)e iϕ ± (x)/h
Omitting indices ± and a, we find a(x; h) = a 0 (x) + ha 1 (x) + • • •; the usual half-density is

a 0 (x) = C 0 C 0 |ψ ′′ (ξ(x))| -1/2 b 0 (ξ(x))
with a new constant C 0 ∈ R ; the next term is

a 1 (x) = ( C 1 + D 1 (x))|β 0 (x)| -1/2 exp -i p 1 (x, ϕ ′ (x)) β 0 (x) dx
and D 1 (x) a complex function with

(3.26) Re D 1 (x) = - 1 2 C 0 β 1 (x) β 0 (x) + Const. Im D 1 (x) = C 0 β 1 (x) β 2 0 (x) p 1 (x, ϕ ′ (x)) dx - p 2 (x, ϕ ′ (x)) β 0 (x) dx and β 0 (x) = ∂ ξ p 0 (x, ϕ ′ (x)) = -α(ξ(x)) ξ ′ (x) , β 1 (x) = ∂ ξ p 1 (x, ϕ ′ (x))
. The homology class of the 1-form defining D 1 (x) can be determined as in Lemma 3.2 and coincides of course with this of T 1 dξ (see (3.9)) on their common chart. In particular, Im D 1 (x) = Im D 1 (ξ(x)) (where ξ(x) stands for ξ ± (x)).

We stress that (3.17) and (3.25) are equal mod O(h 2 ), though they involve different expressions.

Normalization with respect to the "flux norm" as above yields

C 0 = C 0 = 1/ √ 2, and C 1 is determined as in Proposition 3.1. As a result (3.27) u(x; h) = 2∂ ξ p 0 exp h∂ x p 1 ∂ ξ p 0 -1 2 exp iS(x E , x; h)/h (1 + O(h 2 ))
This, together with (3.8), provides a covariant representation of microlocal solutions relative to the choice of coordinate charts, x and ξ being related on their intersection by

-x = ψ ′ (ξ) ⇐⇒ ξ = ϕ ′ (x).
d) Bohr-Sommerfeld quantization rule.

Recall from (3.19) the modified phase function of the microlocal solutions u a ± mod O(h 2 ) from the focal point a E ; similarly this of the other asymptotic solution from the other focal point a ′ E takes the form

(3.28) S ± (x ′ E , x; h) = x ′ E ξ ′ E + x x ′ E ξ ± (y) -h x x ′ E p 1 (y, ξ ± (y)) ∂ ξ p 0 (y, ξ ± (y) dy + h 2 x x ′ E T 1 (ξ ± (y))ξ ′ ± (y) dy
Consider now F a ± (x, h) with asymptotics (3.22), and similarly F a ′ ± (x, h). The normalized microlocal solutions u a and u a ′ , uniquely extended along γ E , are now called u 1 and u 2 . Arguing as for (3.23), but taking now into account the variation of the semi-classical action between a E and a ′ E we get

(3.29) (u 1 |F a ′ + -F a ′ + ) ≡ i 2 e iA -(x E ,x ′ E ;h)/h -e iA + (x E ,x ′ E ;h)/h (u 2 |F a + -F a + ) ≡ i 2 e -iA -(x E ,x ′ E ;h)/h -e -iA + (x E ,x ′ E ;h)/h mod O(h 2 )
, where the generalized actions are given by (3.30)

A ρ (x E , x ′ E ; h) = S ρ (x E , x; h) -S ρ (x ′ E , x; h) = x E ξ E -x ′ E ξ ′ E + x ′ E x E ξ ρ (y) dy -h x ′ E x E p 1 (y, ξ ρ (y)) ∂ ξ p 0 (y, ξ ρ (y) dy + h 2 x ′ E x E T 1 (ξ ρ (y))ξ ′ ρ (y) dy
We have

x E x ′ E ξ + (y) -ξ -(y) dy = γ E ξ(y) dy x E x ′ E p 1 (y, ξ + (y)) ∂ ξ p 0 (y, ξ + (y)) - p 1 (y, ξ -(y)) ∂ ξ p 0 (y, ξ -(y)) dy = γ E p 1 dt x E x ′ E T 1 (ξ + (y))ξ ′ + (y) -T 1 (ξ -(y))ξ ′ -(y) dy = Im γ E Ω 1 (ξ(y)) dy
On the other hand, Gram matrix as in (2.7) has determinant

-cos 2 (A -(x E , x ′ E ; h) -A + (x E , x ′ E ; h))/2h)
which vanishes precisely when BS holds. This brings our alternative proof of Theorem 0.1 to an end.

Bohr-Sommerfeld and action-angle variables.

We present here a simpler approach based on Birkhoff normal form and the monodromy operator [LoRo], which reminds of [HeRo]. Let P be self-adjoint as in (0.1) with Weyl symbol p ∈ S 0 (m), and such that there exists a topological ring A where p 0 verifies the hypothesis (H) in the Introduction.

Without loss of generality, we can assume that p 0 has a periodic orbit γ 0 ⊂ A with period 2π

and energy E = E 0 . Recall from Hamilton-Jacobi theory that there exists a smooth canonical transformation (t, τ ) → κ(t, τ ) = (x, ξ), t ∈ [0, 2π], defined in a neighborhood of γ 0 and a smooth

function τ → f 0 (τ ), f 0 (0) = 0, f ′ 0 (0) = 1 such that (4.1) p 0 • κ(t, τ ) = f 0 (τ )
It is given by its generating function S(τ, x) =

x x 0 ξ dx, ξ = ∂ x S, ϕ = ∂ τ S, and (4.2) p 0 (x, ∂S ∂x (τ, x)) = f 0 (τ )
Energy E and momentum τ are related by the 1-to-1 transformation E = f 0 (τ ), and f ′ 0 (E 0 ) = 1. This map can be quantized semi-classically, which is known as the semi-classical Birkhoff normal form (BNF), see e.g. [GuPa] and its proof. Here we take advantage of the fact (see [CdV], Prop.2) that we can deform smoothly p in the interior of annulus A, without changing its semi-classical spectrum in I, in such a way that the "new" p 0 has a non-degenerate minimum, say at (x 0 , ξ 0 ) = 0, with p 0 (0, 0) = 0, while all energies E ∈]0, E + ] are regular. Then BNF can be achieved by introducing the so-called "harmonic oscillator" coordinates (y, η) so that (4.1) takes the form

(4.3) p 0 • κ(y, η) = f 0 ( 1 2 (η 2 + y 2 ))
and

U * P U = f ( 1 2 y ) 2 + y 2 ; h), has full Weyl symbol f (τ ; h) = f 0 (τ ) + hf 1 (τ ) + • • •.
Here f 1 includes Maslov correction 1/2, and U is a microlocally unitary h-FIO operator associated with κ ( [CdVV], [HeSj]). In A, τ = 0, so we can make the smooth symplectic change of coordinates y = √ 2τ cos t, η = √ 2τ sin t, and take 1 2 (hD y ) 2 + y 2 back to hD t . We do not intend to provide an explicit expression for f j (τ ), j ≥ 1 in term of the p j , but only point out that f j depends linearly on p 0 , p 1 , • • • p j and their derivatives. Of course, BNF allows to get rid of focal points. The section t = 0 in f -1 0 (E) (Poincaré section) reduces to a point, say Σ = {a(E)}. Recall from [LoRo] that Poisson operator K(t, E) here solves (globally near γ 0 ) (4.4) (f (hD t ; h) -E)K(t, E) = 0 and is given in the special 1-D case by the multiplication operator on L 2 (Σ) ≈ C K(t, E) = e iS(t;E)/h a(t; E, h)

where S(t, E) verifies the eikonal equation f 0 (∂ t S) = E, S(0, E) = 0, i.e. S(t, E) = f -1 0 (E)t, and a(t, E; h) = a 0 (t, E) + ha 1 (t, E) + • • • satisfies transport equations to any order in h.

Applying (3.25) in the special case where P has constant coefficients, one has

(4.5) a 0 (t, E) = C 0 (f -1 0 ) ′ (E) 1/2 e it S 1 (E) a 1 (t, E) = C 1 (E) + C 0 (β(E) + it S 2 (E)) (f -1 0 ) ′ (E) 1/2 e it S 1 (E)
with C 0 ∈ R a normalization constant as above to be determined as above (4.6)

S 1 (E) = -f 1 (τ )(f -1 0 ) ′ (E) β(E) = - 1 2 (f -1 0 ) ′ (E)f ′ 1 (τ ) S 2 (E) = (f -1 0 ) ′ (E) 1 2 df 2 1 dE -f 2 (τ )
where we recall τ = f -1 0 (E), so that (4.7)

K(t, E) = e iS(t;E)/h (f -1 0 ) ′ (E) 1/2 e it S 1 (E) C 0 + hC 1 (E) + hC 0 β(E) + ithC 0 S 2 (E)
Together with K(t, E) we define K * (t, E) = e -iS(t,E)/h a(t, E; h), and

K * (E) = K * (t, E) dt
The "flux norm" on C 2 is defined by

(4.8) (u|v) χ = i h [f (hD t ; h), χ(t)]K(t; h)u|K(t, h)v
with the scalar product of L 2 (R t ) on the RHS, and χ ∈ C ∞ (R) is a smooth step-function, equal to 0 for t ≤ 0 and to 1 for t ≥ 2π. To normalize K(t, E) we start from

K * (E) i h [f (hD t ; h), χ(t)]K(t, E) = Id L 2 (R) Since i h [f (hD t ; h), χ(t)] has Weyl symbol (f ′ 0 (τ ))+hf ′ 1 (τ ))χ ′ (t)+O(h 2 ) we are led to compute I(t, E) = i h [f (hD t ; h), χ(t)]K(t, E) where we have set Q(τ ; h) = f ′ 0 (τ ) +hf ′ 1 (τ ).
Again by stationary phase (2.3)

I(t, E) = e iS(t,E)/h Q(τ ; h))χ ′ (t)a(t, E; h) -ih∂ τ Q(τ ; h) 1 2 χ ′′ (t)a(t, E; h) + χ ′ (t)∂ t a(t, E; h) + O(h 2 )
Integrating I(t, E) against e -iS(t,E)/h a(t, E; h), we get (4.9) (u|v

) χ = uv Q(τ ; h) χ ′ (t)|a(t, E; h)| 2 - ih 2 ∂ τ Q(τ ; h) χ ′′ (t)|a(t, E; h)| 2 dt -ih∂ τ Q(τ ; h) ∂ t a(t, E; h)a(t, E; h)χ ′ (t) dt + O(h 2 ) Now |a(t, E; h)| 2 = (f -1 0 ) ′ (E) C 2 0 + 2hC 0 C 1 (E) + 2hC 2 0 β(E) + O(h 2 ) is independent of t mod O(h 2 )
, and

(u|v) χ = uv C 2 0 + 2C 0 C 1 (E)h -C 2 0 α(E)(f -1 0 ) ′ (E)f ′′ 0 (τ ) + O(h 2 )
so that, choosing C 0 = 1 and

C 1 (E) = 1 2 (f -1 0 ) ′ (E) 2 f 1 (τ )f ′′ 0 (τ )
we end up with (u|v) χ = uv(1 + O(h 2 ), which normalizes K(t, E) to order 2.

We define K 0 (t, E) = K(t, E) (Poisson operator with data at t = 0), K 2π (t, E) = K(t -2π, E) (Poisson operator with data at t = 2π), and recall from [LoRo] that E is an eigenvalue of f (hD t ; h) iff 1 is an eigenvalue of the monodromy operator

M (E) = K * 2π (E) i h [f (hD t ; h), χ]K 0 (•, E)
, which in the 1-D case reduces again to a multiplication operator. A short computation shows that

M (E) = exp[2iπτ /h] exp[2iπ S 1 (E)] 1 + 2iπh S 2 (E) + O(h 2 )
so again BS quantization rule writes with an h 2 accuracy as

f -1 0 (E) + h S 1 (E) + h 2 S 2 (E) ≡ nh, n ∈ Z Let S 1 (E) = 2π S 1 (E), and S 2 (E) = 2π S 2 (E). Since f -1 0 (E) = τ (E) = 1
2π γ E ξ dx, and we know that S 3 (E) = 0, we eventuelly get

S 0 (E) + hS 1 (E) + h 2 S 2 (E) + O(h 4 ) = 2πnh, n ∈ Z
Note that the proof above readily extends to the periodic case, where there is no Maslov correction in f 1 .

5. The discrete spectrum of P in I.

Here we recover the fact that BS determines asymptotically all eigenvalues of P in I. As in Sect.1 we adapt the argument of [SjZw], and content ourselves with the computations below with an accuracy O(h). It is convenient to think of {a E } and {a ′ E } as zero-dimensional "Poincaré sections" of γ E . Let K a (E) be the operator (Poisson operator) that assigns to its "initial value" C 0 ∈ L 2 ({a E }) ≈ R the well normalized solution u(x; h) = e i(xξ+ψ(ξ))/h b(ξ; h) dξ to (P -E)u = 0 near {a E }. By construction, we have:

(5.1) ±K a (E) * i h [P, χ a ] ± K a (E) = Id a E = 1
We define objects "connecting" a to a ′ along γ E as follows: let T = T (E) > 0 such that exp T H p 0 (a) = a ′ (in case p 0 is invariant by time reversal, i.e. p 0 (x, ξ) = p 0 (x, -ξ) we take T (E) = T (E)/2). Choose χ a f (f for "forward") be a cut-off function supported microlocally near γ E , equal to 0 along exp tH p 0 (a) for t ≤ ε, equal to 1 along γ E for t ∈ [2 ε, T + ε], and back to 0 next to a ′ , e.g. for t ≥ T + 2 ε. Let similarly χ a b (b for "backward") be a cut-off function supported microlocally near γ E , equal to 1 along exp tH p 0 (a) for t ∈ [-ε, T -2 ε], and equal to 0 next to a ′ , e.g. for t ≥ T -ε. By (5.1) we have

K a (E) * i h [P, χ a ] + K a (E) = K a (E) * i h [P, χ a f ]K a (E) = 1 (5.2) -K a (E) * i h [P, χ a ] -K a (E) = -K a (E) * i h [P, χ a b ]K a (E) = 1 (5.3) which define a left inverse R a + (E) = K a (E) * i h [P, χ a f ] to K a (E) and a right inverse R a -(E) = - i h [P, χ a b ]K a (E) to K a (E) * . We define similar objects connecting a ′ to a, T ′ = T ′ (E) > 0 such that exp T ′ H p 0 (a) = a ′ ( T = T ′ if p 0 is invariant by time reversal), in particular a left inverse R a ′ + (E) = K a ′ (E) * i h [P, χ a ′ f ] + to K a ′ (E) and a right inverse R a ′ -(E) = -i h [P, χ a ′ b ]K a ′ (E) to K a ′ (E) * , with the additional requirement (5.4) χ a b + χ a ′ b = 1 near γ E . Define now the pair R + (E)u = (R a + (E)u, R a ′ + (E)u), u ∈ L 2 (R) and R -(E) by R -(E)u -= R a -(E)u a -+ R a ′ -(E)u a ′ -, u -= (u a -, u a ′ -) ∈ C 2 ,
we call Grushin operator P(z) the operator defined by the linear system (5.5)

i h (P -z)u + R -(z)u -= v, R + (z)u = v +
From [SjZw], we know that the problem (5.5) is well posed, and as in (1.7)-(1.8)

P(z) -1 = E(z) E + (z) E -(z) E -+ (z)
with choices of E(z), E + (z), E -+ (z), E -(z) similar to those in Sect.1. Actually one can show that the effective Hamiltonian E -+ (z) is singular precisely when 1 belongs to the spectrum of the monodromy operator, or when the microlocal solutions u 1 , u 2 ∈ K h (E) computed in (3.29) are colinear, which amounts to say that Gram matrix (2.7) is singular. There follows that the spectrum of P in I is precisely the set of z we have determined by BS quantization rule.

Note that the argument used in Sect.4 would need a slightly different justification, since we made use of a single "Poincaré section".

Appendix: Essentials on 1-D semi-classical spectral asymptotics.

Following essentially [BaWe] [CdV2], we recall here some useful notions of 1-D Microlocal Analysis, providing a consistent framework for WKB expansions in different representations.

a) h-Pseudo-differential Calculus

Semi-classical analysis, or h-Pseudodifferential calculus, is based on asymptotics with respect to the small parameter h. This is a (almost straightforward) generalization of the Pseudo-differential calculus of [Hö], based on asymptotics with respect to smoothness, that we refer henceforth as the "Standard Calculus".

The growth at infinity of an Hamiltonian is controlled by an order function, i.e. m ∈ C ∞ (T * R), m ≥ 1, of temperate growth at infinity, that verifies m ∈ S(m); for instance we take m(x, ξ) = 1 + |ξ| 2 for Schrödinger or Helmholtz Hamiltonians with long range potential, m(x, ξ) = 1 + |x, ξ| 2 for Hamiltonians of the type of a harmonic oscillator (with compact resolvant), or simply m = 1 for a phase-space "cut-off".

Consider a real valued symbol p ∈ S(m) as in (0.1), and define a self-adjoint h-PDO p w (x, hD x ; h) on L 2 (R) as in (0.3).

As in the Standard Calculus, h-PDO's compose in a natural way. It is convenient to work with symbols having asymptotic expansions (0.2). A h-PDO P w (x, hD x ; h) is called elliptic if its principal symbol p 0 verifies |p 0 (x, ξ)| ≥ const. m(x, ξ). If P w (x, hD x ; h) is elliptic then it has an inverse Q w (x, hD x ; h) with q ∈ S(1/m). Ellipticity can be restricted in the microlocal sense, i.e. we say that p is elliptic at ρ 0 = (x 0 , ξ 0 ) ∈ T * R if p 0 (ρ 0 ) = 0, so that P w (x, hD x ; h) has also a microlocal inverse Q w (x, hD x ; h) near ρ 0 .

b) Admissible semi-classical distributions and microlocalization

These h-PDO extend naturally by acting on spaces of distributions of finite regularity H s (R) (Sobolev spaces).

It is convenient to view h-PDO's as acting on a family (u h ) of L 2 -functions, or distributions on R, rather than on individual functions. We call u h admissible iff for any compact set K ⊂ R we have u h H s (K) = O(h -N 0 ) for some s and N 0 . We shall be working with some particular admissible distributions, called Lagrangian distributions, or oscillating integrals.

A Lagrangian distribution takes the form

(A.1) u h (x) = (2πh) -N/2
R N e iϕ(x,θ)/h a(x, θ; h) dθ where a is a symbol (i.e. belongs to some S(m)) and ϕ is a non-degenerate phase function, i.e. d x,θ ϕ(x 0 , θ 0 ) = 0, and d∂ θ 1 ϕ, • • • , d∂ θ N ϕ are linearly independent on the critical set

(A.2) C ϕ = {(x, θ) : ∂ϕ ∂θ (x, θ) = 0}
Such a distribution is said to be negligible iff for any compact set K ⊂ R, and any s ∈ R we have

u h H s (K) = O(h ∞ ).
Remark: Negligible Lagrangian distributions up to finite order, as those constructed in this paper, can be defined similarly. Including more general admissible distributions requires to modify the concept of negligible distributions, as well as the frequency set below, in order to take additional regularity into account. The way to do it is to compactify the usual phase-space T * R by "adding a sphere" at infinity [CdV2]. For simplicity, we shall be content with microlocalizing in T * R, let us only mention that microlocalization in case of Standard Calculus is carried in T * R \ 0, where the zero-section has been removed, and the phase functions enjoy certain homogeneity properties in the phase variables.

Microlocal Analysis specifies further the "directions" in T * R where u h is "negligible". To this end, we introduce, following Guillemin and Sternberg, the frequency set FS u h ⊂ T * R by saying that ρ 0 = (x 0 , ξ 0 ) / ∈ F Su h iff there exists a h-PDO A with symbol a ∈ S 0 (m) elliptic at ρ 0 and such that Au h is negligible. Since this definition doesn't depend of the choice of A, and we can take A = χ w (x, hD x ) where χ ∈ C ∞ 0 (T * R) is a microlocal cut-off equal to 1 near ρ 0 . On the set of admissible distributions, we define an equivalence relation at (x 0 , ξ 0 ) ∈ T * R by u h ∼ v h iff (x 0 , ξ 0 ) / ∈ FS(u hv h ), and we say that u h = v h microlocally near (x 0 , ξ 0 ).

As in Standard Calculus, if P ∈ S(m) we have (A.3) FS P u h ⊂ FS u h ⊂ FS P u h ∪ Char P

where Char P = {(x, ξ) ∈ T * R : p 0 (x, ξ) = 0} is the bicharacteristic strip.

For instance, eigenfunctions of P w (x, hD x ; h) with energy E (as admissible distributions) or more generally, solutions, in the microlocal sense, of (P w (x, hD x ; h) -E)u h ∼ 0 are "concentrated" microlocally in the energy shell p 0 (x, ξ) = E, in the sense that FS u h ⊂ Char(P -E). It follows that FS u h is invariant under the flow t → Φ t of Hamilton vector field H p 0 . Assume now that P -E is of principal type (i.e. H p 0 = 0 on p 0 = E), the microlocal kernel of P -E is (at most) one-dimensional, i.e. if u h , v h are microlocal solutions and u h ∼ v h at one point (x 0 , ξ 0 ), then u h ∼ v h everywhere.

The existence of WKB solutions (see below) ensures that the microlocal kernel of P -E is indeed one-dimensional. This fails of course to be true in case of multiple caracteristics, e.g. at a separatrix.

It is convenient to characterize the frequency set in terms of h-Fourier transform (A.4) F h u h (ξ) = (2πh) -1/2 e -ixξ/h u h (x) dx Namely ρ 0 / ∈ FS h (u h ) iff there exists χ ∈ C ∞ 0 (R), χ(x 0 ) = 0, and a compact neighborhood V of ξ 0 such that F h (χ u h )(ξ) = O(h ∞ ) uniformly on V .

Note as above that the frequency set may include the zero section ξ = 0, contrary to the standard wave-front WF, see also [Iv].

Examples:

1) "WKB functions" of the form u h (x) = a(x) e iS(x)/h with a, S ∈ C ∞ , S real valued. We have FS h (u h ) = (x, S ′ (x)) : x ∈ supp(a) . More generally, if u h is as in (A.1) then FS h (u h ) is contained in the Lagrangian manifold Λ ϕ = {(x, ∂ x ϕ(x, θ)) : ∂ θ ϕ(x, θ) = 0}, with equality if a(x, θ; h) = 0 on the critical set C ϕ was defined in (A.2).

2) If u(x) is independent of h, then FS h (u) = WF u ∪ (supp(u) × {0}).

Fourier inversion formula then shows that if U ⊂ R n is an open set, an h-admissible family (u h ) is negligible in U iff π x (FS(u h )) ∩ U = ∅, where π x denotes the projection T * R → R x . So FS u h = ∅ iff u h are smooth and small (with respect to h) in Sobolev norm.

c) WKB method

When P -E is of principal type, and H p 0 is transverse to the fiber in T * R, we seek for microlocal solutions of WKB type, of the form u h (x) = e iS(x)/h a(x; h), where a(x; h) ∼ ∞ j=0 h i j a j (x). Applying P -E, we get an asymptotic sum, with leading term p 0 (x, S ′ (x)) = E, which is the eikonal equation, that we solve by prescribing the initial condition S ′ (x 0 ) = ξ 0 , where p 0 (x 0 , ξ 0 ) = E. The lower order terms are given by (in-)homogeneous transport equations, the first transport equation takes the invariant form L H p 0 a 0 = 0, where L H p 0 denote Lie derivative along H p 0 . Hence e iS(x)/h a 0 (x) gives the Lagrangian manifold Λ S together with the half density a 0 (x) √ dx on it. The right hand side of higher order (non-homegeneous) transport equations or order j involve combinations of previous

a 0 , • • • , a j-1 .
When H p 0 turns vertical, we switch to Fourier representation as in Sect.3. Matching of solutions in such different charts can be done using Gram matrix since, P -E being of principal type, there is only one degree of freedom for choosing the microlocal solution.

  Maslov correction and the integral of the sub-principal 1-formp 1 dt S 1 (E) = πp 1 (x(t), ξ(t))| γ E dt (iii) the second order term S 2 (E) = 1 24

  a) Quasi-modes mod O(h 2 ) in Fourier representation. Let a = a E = (x E , ξ E ) be such a focal point. Following a well known procedure we can trace back to [Sj1], we first seek for WKB solutions in Fourier representation near a of the form u(ξ) = e iψ(ξ)/h b(ξ; h), see e.g. [CdV2] and Appendix below. Here the phase ψ = ψ E solves Hamilton-Jacobi equation p 0 (-ψ ′ (ξ), ξ) = E, and can be normalized by ψ(ξ E ) = 0; the amplitude b(ξ; h) = b 0 (ξ) + hb 1 (ξ) + • • • has to be found recursively together with a(x, ξ

  Λ E is smoothly embedded in T * R, it can always be parametrized by a non degenerate

	phase function. Of particular interest are the critical points of the phase functions, or focal points
	which are responsible for the change in Maslov index. Recall that a(E) = (x E , ξ E ) ∈ Λ E is called a focal point if Λ E "turns vertical" at a(E), i.e. T a(E) Λ E is no longer transverse to the fibers x = Const.
	in T * R. In any case, however, Λ E can be parametrized locally either by a phase S = S(x) (spatial
	representation) or a phase S = S(ξ) (Fourier representation). Choose an orientation on Λ E and for a ∈ 0 (R 2 ) Λ E (not necessarily a focal point), denote by ρ = ±1 its oriented segments near a. Let χ a ∈ C ∞ be a smooth cut-off equal to 1 near a, and ω a ρ a small neighborhood of supp[P, χ a ] ∩ Λ E near ρ.

ξ ξ E

Acknowledgements: We thank a referee for his constructive remarks. This work has been partially supported by the grant PRC CNRS/RFBR 2017-2019 No.1556 "Multi-dimensional semi-classical problems of Condensed Matter Physics and Quantum Dynamics".