A Ifa

M Rouleux
email: rouleux@univ-tln.fr

The one dimensional semi-classical Bogoliubov-de Gennes Hamiltonian with PT symmetry: generalized Bohr-Sommerfeld quantization rules

We present a method for computing first order asymptotics of semiclassical spectra for 1-D Bogoliubov-de Gennes (BdG) Hamiltonian from Supraconductivity, which models the electron/hole scattering through two SNS junctions. This involves: 1) reducing the system to Weber equation near the branching point at the junctions; 2) constructing local sections of the fibre bundle of microlocal solutions; 3) normalizing these solutions for the "flux norm" associated to the microlocal Wronskians; 4) finding the relative monodromy matrices in the gauge group that leaves invariant the flux norm; 5) from this we deduce Bohr-Sommerfeld (BS) quantization rules that hold precisely when the fibre bundle of microlocal solutions (depending on the energy parameter E) has trivial holonomy. Such a semi-classical treatement reveals interesting continuous symetries related to monodromy. Details will appear elsewhere.

Bogoliubov-de Gennes Hamiltonian

BdG Hamiltonian describes the dynamics of a pair of quasi-particles electron/hole in the Theory of Supraconductivity [START_REF] Bardeen | [END_REF]. We consider a narrow metallic 1-D wire (Normal Metal N) connected to Supraconducting bulks S through a SNS junction, and compute the excitation spectrum in the normal contact region as a function of gate voltage, when electronic levels transform into phase sensitive Andreev levels. The wire, or lead, is identified with a 1-D structure, the interval x ∈ [-L, L] (case of a perfect junction) or x ∈ [-L + ℓ/2, L -ℓ/2] ("dirty junction"), where ℓ ≪ L. The reference energy in the lead is Fermi level E F . The pair electron/hole is acted upon by two kinds of potentials:

(1) the "order parameter" ∆(x) times a phase function e iφ(x)/2 , which is the potential due to Cooper pairs in the supraconducting bulk. This potential, subject to self-consistency relations, is priori unknown. Namely, inside S, ∆(x)e iφ(x)/2 is a solution of Ginzburg-Landau (or Pitaevskiy) equations, and shows typically a vortex profile (in 2-D). In BdG Hamiltonian it is assumed, however, that ∆(x)e iφ(x)/2 is an "effective" potential. Inside N, superconducting gap ∆(x) ≡ 0: quasi-particles live in the "clean metal". For |x| ≥ L + ℓ, ∆(x) = ∆ 0 > 0.

We assume that the phase function φ(x) is constant near the junction, and gauge the interaction by φ -= -φ + = -φ in the superconducting banks, so that φ(x) = sgn(x)φ. We assume further that this equality holds everywhere: since ∆(x) = 0 inside N, the discontinuity of x → φ(x) is irrelevant.

(2) a smooth chemical potential µ(x): typically µ(x) is flat in N and drops smoothly to the band bottom in the superconducting banks S. In our model we assume again µ(x) to be constant in the superconducting bank, i.e. µ(x) = µ 0 when |x| ≥ L + ℓ. Andreev currents at energy E occur only if µ(x) ≥ E in [-L, L].

The case of a perfect junction (∆ "hard-wall potential") has been considered in [5], see also [START_REF] Cayssol | [END_REF] for a SFS junction, and makes use scattering matrix techniques. In this work, justifying semiclassical techniques as in [8] (also in the multi-dimensional case) we rather consider an imperfect (or "dirty") junction: ∆(x)e iφ(x)/2 is a smooth function. In a neighborhood of [-L, L], say x ∈ [-L -ℓ, L + ℓ], the system is described at the classical level by BdG Hamiltonian

P(x, ξ) = ξ 2 -µ(x) ∆(x)e iφ(x)/2 ∆(x)e -iφ(x)/2 -ξ 2 + µ(x) (1)
The energy surface:

Σ E = {det(P -E) = -(ξ 2 -µ(x)) 2 -∆(x) 2 + E 2 = 0} = Λ < E ∪ Λ > E
splits into 2 branches separated in momentum space, so consists of two microlocal wells. Interaction between these wells gives the imaginary parts of the resonances for the electron/hole scattering, and will be ignored in this paper. Because of smoothness of x → ∆(x), the reflections occur inside [-L, L], we denote by (±x E , ξ E ) ∈ Λ > E , the one-parameter family of "branching points" defined by ∆(±x

E ) = E with x E near x 0 ∈ [L -ℓ 2 , L + ℓ 2 ], ∆(x 0 ) > 0.
We do not consider the problem of "clustering" of eigenvalues as E → 0 = E F (Fermi level). In the "hard wall potential" limit for x near x 0 , the potential ∆(x) can be safely approximated by a linear function such that ∆(x 0 ) = E 0 , and µ(x) by a constant µ. So near x 0 we assume that

φ(x) = φ, µ(x) = µ > E, ∆(x) = E + α(x -x E ) for large α > 0. Condition a E = (x E , ξ E ) ∈ Σ E gives ξ 2 E = µ > E, ∆(x E ) = E.
The physical mechanism goes roughly as follows (see [5] for a detailed exposition): An electron e -moving in the metallic lead, say, to the right, with energy 0 < E ≤ ∆ below the gap and kinetic energy K + (x) = µ(x) + E 2 -∆(x) 2 is reflected back as a hole e + from the supraconductor, injecting a Cooper pair into the superconducting contact. The hole has kinetic energy 2 , and a momentum of the same sign as this of the electron. When inf [-L,L] K -(x) > 0 it bounces along the lead to the left and picks up a Cooper pair in the supraconductor, transforming again to the original electron state, a process known as Andreev reflection. This works also the other way in Λ < E , since Hamiltonian system conserves both charge and energy. Actually, the hole can propagate throughout the lead only if inf [-L,L] µ(x) ≥ E. Otherwise, it is reflected from the potential µ(x) in the junction, and Andreev levels are quenched at higher energies, i.e. transform into localized electronic states.

K -(x) = µ(x) -E 2 -∆(x)
For a rescaled "Planck constant" h so that h ≪ ℓ, we consider Weyl h-quantization of BdG Hamiltonian P(x, hD x ) on L 2 (I) ⊗ C 2 , I = [-(L + ℓ), L + ℓ], which is self-adjoint when imposing Dirichlet boundary conditions at ∂I. Phase-sensitive Andreev states carry supercurrents that turn out to be proportional to the φ-derivative of the eigen-energies of P(x, hD x ).

We have σ y P(φ)σ y = -P(-φ), with σ y = 0 -i i 0 , accounting for "negative energies". We shall assume here E > 0. When potentials are even functions (typical for metals), P(x, hD x ) verifies PT symmetry ∨ IP(x, hD x ) = P(x, hD x )I ∨ which is essential for our approach to work. At least formally, since BdG is only defined locally near N, removing boundary conditions leads to "resonances" (i.e. metastable states or quasi-particles with a finite life-time). Thus for simplicity we have assumed that (1), together with its semi-classical quantization, describes the system not only in I, but on the whole real line, provided h ≪ ℓ ≪ L. Thus P(x, hD x ) extends to L 2 (R) ⊗ C 2 , Our general goal is to give a precise mathematical meaning to these "resonances". Here we content to compute their real parts through Bohr-Sommerfeld quantization rules.

2. Monodromy operator, scattering matrix: an outlook a) Schrödinger operator on the real line.

We first recall from [START_REF] Arnold | Geometrical methods in the theory of ordinary differential equations[END_REF] basic facts for a 1-D Schrödinger operator with a compactly supported potential V . The generalized wave-functions u with energy E = k 2 > 0 satisfy

-h 2 u ′′ (x) + V (x)u(x) = Eu(x) (2)
and outside supp V ,

-h 2 u ′′ (x) = k 2 u(x) (3)
defines the state space Z ≈ C 2 of the "free particle", spanned by f 1 (x) = e ikx/h , f 2 (x) = e -ikx/h . The monodromy operator M (k) :

f 1 + Bf 2 → Af 1 is such that M (k) = 1/A -B/A -B/A 1/A ∈ SU(1, 1)
In particular, |A| 2 + |B| 2 = 1. We call |A| 2 the transmission coefficient and |B| 2 the reflection coefficient. Along with the passage from the left to the right of the support of V , consider the passage from the right to the left. The corresponding solution v of ( 2) is e -ikx/h + B 2 e ikx/h to the right of suppV , and A 2 e -ikx/h to the left. The scattering matrix is defined as

S(k) = A B -BA/A A ∈ U(2)
S(k) remains unitary and symmetric for complex values of k. Resonances of (2) are then defined as

E = k 2 ∈ C,
where k is a pole of S, and physical resonances those with Imk > 0. Thus E is a resonance iff the solution of ( 3) is purely outgoing as x → +∞ and x → -∞. The poles coincide with the poles of meromorphic extension of the resolvent (P -k 2 ) -1 from the physical half-plane ImE < 0 to the second sheet ImE > 0.

b) Monodromy matrix for BdG equation: heuristics. Now we discuss BdG equation (P(x, hD x ) -E)U = 0 for large |x|, i.e. (within our approximation above) when |x| ≥ L + ℓ, so ∆(x) = ∆ 0 , µ(x) = µ 0 > E. Solutions are of the form

U (x; h) = a b c d e ikx/h e iℓx/h µ 0 + E ± i∆ 0 ∈ {k 2 , ℓ 2 }, so eigenfrequencies are (±k, ±k), k = √ µ 0 + E + i∆ 0 ,

and the corresponding solutions as follows:

Let φ(x) = sgn(x)φ, Z be the 2-D complex line bundle spanned by

F ± 1 (x) = e iφ(x)/2
-i e ±ikx/h (associated with the scattering process e + → e -), and Z the 2-D complex line bundle spanned by F ± 2 (x) = e iφ(x)/2 i e ±ikx/h (associated with the scattering process e -→ e + ). The space of solutions of exponential type for BdG is Z ⊕ Z, and Z, Z are orthogonal for the usual pointwise Hermitian product in C 2 . Declare that E ∈ C is a Z-resonance iff the Z-component of the wave function solving BdG equation is outgoing and evanescent ("physical solution") at infinity, i.e.

U (x, h)

= A e iφ/2 -i e ikx/h , x → +∞ U (x, h) = B e -iφ/2 -i e -ikx/h , x → -∞
Similarly we say that E is a Z-resonance iff the Z-component of the wave function is outgoing (and evanescent) at infinity, i.e.

U (x, h) = A e iφ/2 i e -ikx/h , x → +∞ U (x, h) = B e -iφ/2 i e ikx/h , x → -∞
So for both sets of resonances, the corresponding solution is simultaneously decaying, and outgoing at ±∞. These sets of resonances need not coincide (although they come up in pairs), but their real parts are given by Bohr-Sommerfeld quantization rules. Namely, define the monodromy operator M Z (k) acting on Z according to the formula

e -iφ/2 -i e ikx/h + B e -iφ/2 -i e -ikx/h → A e iφ/2 -i e ikx/h
and similarly for

M Z (k). It is plausible to expect that M Z (k), M Z (k) ∈ U(1, 1)
, and that the corresponding scattering matrices S Z (k), S Z (k) have a meromorphic extension to the complex plane, their poles defining the resonances E Z and E Z . Actually, we shall construct "relative monodromy operators" in the "classically allowed region". In particular the relative monodromy operators are in U(1,1) for some specific Lorenzian form which is constructed below.

Bohr-Sommerfeld quantization rules

In this work, we content to determine the real parts of the resonances, extending to this setting the method of positive commutators elaborated in [START_REF] Sjöstrand | Density of states oscillations for magnetic Schrödinger operators[END_REF], [9] and [10]. Imaginary parts may be determined as in [START_REF] Rouleux | Tunneling effects for h-Pseudodifferential Operators, Feshbach resonances and the Born-Oppenheimer approximation[END_REF]. We obtain Bohr-Sommerfeld quantization rules for the quasi-particle, alternating even and odd quantum numbers associated with the electron and the hole. In the sequel we will sketch a proof of the following result:

Theorem 1: Let x 0
-x 0 η ρ (y; h) dy be the semi-classical actions (see Proposition 8 below) ρ = 1 for the electron, ρ = -1 for the hole. Bohr-Sommerfeld quantization conditions near E 0 are given at first order by:

γ E η ρ (y; h) dy -hφ + hπ + O(h 2 ) = 2πnh; n ∈ Z
Here γ E denotes integral over the loop γ E obtained by gluing together Λ > E and Λ < E , if we ignore tunneling in momentum space.

Microlocal solutions in Fourier representation near the branching points a) Reduction of the system.

In h-Fourier representation, F h u(ξ) = (2πh) -1/2 e -ixξ/h u(x) dx the local Hamiltonian near a = a E = (x E , ξ E ), P a takes the form :

P a (-hD ξ , ξ) = ξ 2 -µ e iφ/2 (E -αhD ξ -αx E ) e -iφ/2 (E -αhD ξ -αx E ) -ξ 2 + µ (4)
By PT symmetry P a ′ = IP a I near a ′ = a ′ E = (-x E , ξ E ). Solving the system P a (-hD ξ , ξ) U = 0, U = ϕ 1 ϕ 2 gives second order ODE for u(ξ) = exp[-i ξ g(s)ds/h] ϕ 2 (ξ),

P a (-hD ξ , ξ, h)u(ξ) = E 2 α 2 u(ξ) (5)
P a (-hD ξ , ξ, h) = (hD ξ ) 2 + α -2 (ξ 2 -µ) 2 + h 2 (ξ 2 -µ -E) -2 (2ξ 2 + µ + E) After E-dependent scalings β = √ α(2ξ E ) -3/2 > 0, E 1 = (2ξ E ) -2 E, ξ = 2ξ E βωξ ′ + ξ E , ω = ±1 (ξ ′ is "local momentum") we obtain P a ω (-hD ξ ′ , ξ ′ , h)u ω (ξ ′ ) = E 1 β 2 u ω (ξ ′ )
, where

P a ω (-hD ξ ′ , ξ ′ ; h) = (-hD ξ ′ ) 2 + (ξ ′ + βωξ ′2 ) 2 + h 2 β 2 f (ωβξ ′ )
is an anharmonic Schrödinger operator. The lower order term f (z) = (2z 2 + 2z + 3 4 + E 1 )(z 2 + z -E 1 ) -2 has a pole on Λ > E where the linear approximation of ∆(x) breaks down. The linear approximation only holds for small ξ ′ . Consider the map

ι a : ω=±1 Ker h (P a ω - E 1 ω β 2 ) → Ker h (P a -E) (6)
where Ker h denotes the microlocal kernel. The index ω is to be chosen carefully with the complex germ of solutions having the right decay beyond the branching points ±x E . We shall endow the RHS of ( 6) with a Lorenzian structure and "diagonalize" ι a in some orthogonal subspaces.

b) The normal form of Helffer-Sjöstrand When E 1 < 1 4 , we take P a ω microlocally to its normal form, namely:

Proposition 2 [9]: There exists an analytic diffeomorphism t → F 0 (t) defined in a neighborhood of 0, F 0 (0) = 0, with inverse G 0 , and a real analytic phase function φ β (ξ ′ , θ), defined in a neighborhood of (0,0), of the form

φ β (ξ ′ , θ) = ξ ′ θ+g β (ξ ′ , θ), g β (ξ ′ , θ) = O(|ξ ′ , θ| 3 ), parametrizing the canonical transformation κ β : (∂ θ φ β , θ) → (ξ ′ , ∂ ξ ′ φ β ), such that F 0 • p β • κ β = p 0 .
At the semi-classical level, there is a (formally) unitary FIO operator A defined microlocally near (0,0)

Av(ξ ′ , h) = (2πh) -1 e iϕ(ξ ′ ,η,θ)/h c(ξ ′ , η, θ, h)e ib(ξ ′ ,η,θ,h) v(η, h) dηdθ
and a real valued analytic symbol

F (t, β, h) = F 0 (t, β) + hF 1 (t, β) + h 2 F 2 (t, β) + • • • with F 1 (t, β) = -1 2 such that A * F (P ω , β, h)A = P 0 (η, hD η ) = 1 2 (hD η ) 2 + η 2 -h , A * A ≡ Id
The function F 0 , taking the period T (E) of Hamilton vector flow for P a ω at energy (E 1 /β) 2 to 2π, involves an elliptic integral, which requires sometimes the use of formal calculus.

c) Weber equation and parabolic cylinder functions

Weber equation

P 0 v = νhv, through change of variables η = (h/2) 1/2 ζ, v(ζ) = v(η) scales to -v ′′ + 1 4 ζ 2 v = ν + 1 2 v
Fundamental solutions express as parabolic cylinder functions D ν , entire in C. The systems D ν (±ζ), D -ν-1 (±iζ) are fundamental solutions for any choice of ±. Integral representations give asymptotic solutions of (P 0 -νh)u(η) = 0 by stationary phase for real ν, We apply asymptotic stationary phase to AD j , j ∈ {ν, -ν -1}. With h ′ = β 2 h as a "rescaled" Planck constant, we get:

E ′ 2 1 = 2β 2 F (β -2 E 2 1 , β, h) = 2β 2 (ν + 1)h. D ν ε(h/2) -1/2 η = Γ(ν+1) -2iπ √ h h E 2 /4h (0 + ) ∞ exp iΦ ν ε (s; η)/h ds D -ν-1 iε(h/2) -1/2 η = Γ(-ν) 2iπ h -E 2 /4h (0 + ) ∞ exp iΦ -ν-1 ε (s; η)/h ds s with ε = ±1, E = 2(ν + 1)h,
Proposition 3: In Fourier representation, the image K a h (E) = Ker h (P a (-hD ξ , ξ) -E) of ι a is a 2-D vector space spanned by the spinors U j ε,ω = ϕ 1 ϕ 2 j ε,ω , (j, ε, ω) ∈ {ν, -ν -1} × {-1, 1} 2 , of the form: U ν ε,ω = C ν h ′ θω=± θω(ξ 1 ) e iφ/2 (ξ 2 -µ-E) -1/2 X ν ε,ω (ξ 2 -µ-E) 1/2 | a ν ε,ω | exp[i(Φ ν ε,ω + h ′ R ν ω )/h ′ ] + O(h ′ ) U -ν-1 ε,ω = C -ν-1 h ′ θω=± θω(ξ 1 ) ε sgn(θ ω ) e iφ/2 (ξ 2 -µ-E) -1/2 X -ν-1 ε,ω (ξ 2 -µ-E) 1/2 | a -ν-1 ε,ω | exp[i(Φ -ν-1 ε,ω + h ′ R -ν-1 ω )/h ′ ] + O(h ′ )
Here θ ω (ξ 1 ) is a critical point (from stationary phase), Φ j ε,ω + h ′ R j ω ) the h ′ -dependent phase functions, and X j ε,ω , | a j ε,ω | some positive amplitudes. Spinors U j ε,ω verify the symmetry † U j -ε,-ω = U j ε,ω for the "local time" reversal operator † u(ξ 1 ) = u(-ξ 1 ), and the constants

C j h ′ (from Whittaker normalization of D ν , D -ν-1 ) are related by C ν h ′ C -ν-1 h ′ = (2 √ h ′ ) 3 π 2 sin πν -1 .

Normalization

a) The microlocal Wronskian.

We extend to BdG Hamiltonian the classical "positive commutator method" using conservation of some quantity called a "quantum flux' ( [START_REF] Sjöstrand | Density of states oscillations for magnetic Schrödinger operators[END_REF], [9], [START_REF] Rouleux | Tunneling effects for h-Pseudodifferential Operators, Feshbach resonances and the Born-Oppenheimer approximation[END_REF], [10]).

Definition 4: Let P be (formally) self-adjoint, and U a , V a ∈ K h (E) be supported on Λ > E . We call the sesquilinear form

W a ρ (U a , V a ) = i h [P, χ 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 microlocally on ω a ρ (a small neighborhood of supp[P, χ a ] ∩ Λ E near ρ).

A crucial property of the microlocal Wronskian is to be invariant by Fourier transformation:

W a ρ (U a , V a ) = W a ρ ( U a , V a ). The relation W a + (U a , V a ) + W a -(U a , V a )
= 0 doesn't readily follow as in the scalar case [10], the microlocal solutions being neither smooth in spatial of Fourier representation near the branching point, but from a careful inspection, involving also formal calculus. This is used essentially in Propositions 5 and 8 below. Choosing ε, ω such that εω = 1 we define a Lorenzian metric W ρ on the space of microlocal solutions near a. In the basis U j ε,ω , j ∈ {ν, -ν -1} we have, up to a constant factor:

ρW ρ = |C ν h ′ | 2 O(h ′ ) C ν h ′ C -ν-1 h ′ exp[-iπE ′ 1 2 /4h ′ ] 1 + O(h ′ ) C ν h ′ C -ν-1 h ′ exp[iπE ′ 1 2 /4h ′ ] 1 + O(h ′ ) |C -ν-1 h ′ | 2 O(h ′ )
Changing Whittaker normalization for the D ν , D -ν-1 functions, and the microlocal solutions by some constant phase factors, we can reduce to ρW ρ = 0 1 1 0 + O(h ′ ), and prove:

Proposition 5: Under PT symmetry above the microlocal Wronskians W a ρ endow K a h (E) (mod h ′ ) with a Lorenzian form W a = 1 2 (W a + -W a -).
The same holds at a ′ , and the corresponding structures on

K a h × K a * h and K a ′ h × K a ′ * h
are anti-isomorphic. The group of automorphisms preserving W a and W a ′ mod O(h ′ ) is therefore U(1,1).

Spinors in the spatial representation

We compute U a,j ε,ω , U a ′ ,j ε,ω in spatial representation, then extend along the branches ρ = ±1 of Λ > E with WKB solutions.

a) Spinors near the branching points.

Near a, a ′ we apply inverse h-Fourier transform and get:

Proposition 6: Up to a constant phase factor

U ν ε,ω (x, h) = 2ωβξ E e ixξ E /h ρ=± e iφ/2 (ξ 2 -µ-E) -1/2 X ν ε,ω (ξ 2 -µ-E) 1/2 |a ν ε,ω | θ 1 =θω(ξ 1 ),ξ 1 =ξ ρ ω (x) × L ρ ω (x) i -1/2 exp[i Ψ ν,ρ ε,ω (x) + h ′ R ν,ρ ε,ω (x) /h ′ ](1 + O(h ′ )) U -ν-1 ε,ω (x, h) = 2ωβξ E e ixξ E /h ρ=± ε sgn(θ 1 ) e iφ/2 (ξ 2 -µ-E) -1/2 X -ν-1 ε,ω (ξ 2 -µ-E) 1/2 | a -ν-1 ε,ω | θ 1 =θω(ξ 1 ),ξ 1 =ξ ρ ω (x) × L ρ ω (x) i -1/2 exp[i Ψ ν,ρ ε,ω (x) + h ′ R ν,ρ ε,ω (x) /h ′ ](1 + O(h ′ ))
Here L ρ ω (x) -1/2 is a real density (singular at x = x E ), and ρ labels the branch of the Lagrangian manifold. The phases Ψ j,ρ ε,ω (x) + h ′ R j,ρ ε,ω (x), j ∈ {ν, -ν -1} differ only by a constant.

b) WKB spinors away from the branching points

The Lagrangian manifold Λ > E consists of 2 branches Λ >,ρ E (or simply ρ) ρ = ±1 so that ρ = +1 belongs to the electronic state (ξ 1 > 0 in the local coordinates near a above), resp. ρ = -1 to the hole state (ξ 1 < 0). These states mix up when ∆(x) = 0, but we can sort them out semiclassically, outside a, a ′ . Call the vector space of C 2 generated by 1 0 the space of (pure) electronic states, or electronic spinors, and this by 0 1 the space of (pure) hole states, or hole spinors.

The principal symbol P(x, ξ) has eigenvalues λ ρ = ρλ(x, ξ) = ρ ∆(x) 2 + (ξ 2 -µ(x)) 2 . By diagonalizing, we obtain a line bundle Λ ρ E with fiber Looking at the electronic state, we choose ρ = +1 so that λ ρ (x ρ , ξ ρ )-E = 0, while λ -ρ (x ρ , ξ ρ )-E is elliptic. and similarly when looking at the hole state.

Y ρ (x, ξ) = (∆ 2 + (-ξ 2 + µ + ρ ∆ 2 + (ξ 2 -µ) 2 ) 2 ) -1/2 ∆e iφ/2 -ξ 2 + µ + ρ ∆ 2 + (ξ 2 -µ) 2

Proposition 7

The microlocal kernel Ker h (P -E) on Λ >,ρ E is one-dimensional space spanned by W ρ (x, h) = e iSρ(x,h)/h w ρ 0 (x, h)Y ρ (x, ∂ x S ρ ) + O(h) = e iSρ(x,h)/h W ρ (x, h) where w ρ 0 (x)|dx| 1/2 is a smooth half-density. By the uniqueness property of WKB solutions along simple bicharacteristics, the h (or h ′ )-dependent phase function S ρ (x, h) should coincide, up to a constant (in a punctured neighborhood of a) with either one of Ψ j,ρ ε,ω (x) + h ′ R j,ρ ε,ω (x) above, j ∈ {ν, -ν -1}, and similarly for the half-densities.

Relative monodromy matrices

Now we look for connexion formulas. For each ε, ω, ρ = ±1, j ∈ {ν, -ν -1}, the normalized microlocal solutions U a ′ ,j,ρ ε,ω are related to the extension U a,k,ρ -ε,-ω,ext of the normalized microlocal solutions U a,k,ρ ε,ω along the bicharacteristics by a monodromy matrix

M a,a ′ ,ρ = d ρ 11 d ρ 12 d ρ 21 d ρ 22 ∈ U (1, 1)
(defined at least mod O(h ′ )) which we call a relative monodromy matrix. Since there is a pair of particles, the symmetry between the M a,a ′ ,ρ and M a ′ ,a,ρ is order 4; M a ′ ,a,ρ ∈ U (1, 1) is obtained by extending from the left to the right, and applying symmetry

ρM a ′ ,a,ρ = I(M a,a ′ ,ρ ) -1 I =, ρ = ±1 (7)
where I denotes complex conjugation. We compute the coefficients d ij = d ρ ij . Considering behavior of U a ′ ,j,ρ e,ω in the classically forbidden region (according to scattering process e + → e - or e -→ e + ) we obtain

M a,a ′ ,ρ = 0 d ρ 12 d ρ 21 0 , d ρ 12 d ρ 21 = 1
Note that if we do not look too closely at the relevant complex branches, as is the case when computing BS, it makes no difference to choose instead M a,a ′ ,ρ = d ρ 11 0 0 d ρ

22

, with d 11 ρ d ρ 22 = 1. As in [START_REF] Sjöstrand | Density of states oscillations for magnetic Schrödinger operators[END_REF], [9], [START_REF] Rouleux | Tunneling effects for h-Pseudodifferential Operators, Feshbach resonances and the Born-Oppenheimer approximation[END_REF], [10], the argument consists now in extending microlocal solutions obtained above from a to a ′ , and computing the resulting semi-classical action. So take first U 1 equal to U a 1 = U ν,a ε,ω near a, extend it along to a ′ along the bicharacteristics ρ = ±1 by WKB. Evaluating on ρ near a ′ we find U a ′ ,ρ

1 = U ν,a,ρ ε,ω,ext = d ρ 21 U -ν-1,a ′ ,ρ ε,ω
. Similarly, take U 2 starting at a ′ and with

-ν -1 instead of ν, we get U a,ρ 1 = U -ν-1,a ′ ,ρ ε,ω,ext = e ρ 12 U ν,a,ρ ε,ω
, where e ρ 12 = ρ d ρ

21

-1 is the matrix element of M a ′ ,a,ρ given in [START_REF] Gérard | [END_REF]. We compute d ρ 21 in two different ways and compare the result. (1) Using time-reversal and PT symmetries in the microlocal Wronskians, we get

i h [P a ′ , χ a ′ ] ρ U 1 |U ν ε,ω = d ρ 21 i h [P a ′ , χ a ′ ] ρ U -ν-1 ε,ω |U ν ε,ω = = d ρ 21 W a ′ ρ U -ν-1 ε,ω , U ν ε,ω = d ρ 21 W a ′ ρ U -ν-1 ε,ω , U ν ε,ω = = -d ρ 21 W a ρ U -ν-1 -ε,-ω , U ν -ε,-ω = -d ρ 21 W a ρ U ν -ε,-ω , U -ν-1 -ε,-ω = -d ρ

21

(2) Using the extensions described in Proposition 7. Near a ′ we have

U ρ 1,ext = e iφ/2 W ρ (x, h) = d ρ 21 U -ν-1,a ′ ,ρ ε,ω
(by solving transport equation along ρ the amplitude picks up the phase factor e iφ/2 ), so we need to compute

i h [P a ′ , χ a ′ ] ρ W ρ (x, h)|U ν ε,ω . The amplitude W ρ (x, h) is actually defined up to a real, constant factor C ρ . Proposition 8: Let Ψ ν,a ′ ,ρ ε,ω (x) = x ξ E + (2 ξ E ) 3 α Ψ ν,a ′ ,ρ ε,ω (x). We have i h [P a ′ , χ a ′ ] ρ W ρ |U ν,a ′ ,ρ ε,ω = 2 C ρ e iπ/4 exp i S ρ (x; h)/h β(x, h) (χ a ′ 1 ) ′ (x) dx (8)
where the amplitude β(x, h), real mod O(h), is computed from the WKB solutions in Proposition 7, and

S ρ (x, h) = S ρ (x; h) -xξ E + Ψ ν,a ′ ,ρ) ε,ω (x) -h R ν -ω θ -ω (ξ ρ -ω (-x)) = (2 ξ E ) 3 α
Ψ ν,a,ρ -ε,-ω (x 0 ) -x 0 -x 0 η ρ (y; h) dy + h R ν -ω θ -ω (0) Moreover, β(x, h) is also independent of x, so that, comparing the former expression (1) and ( 8) for a suitable choice of C ρ , we get d ρ 21 = -e iτ ρ (h)/h (χ a ′ 1 ) ′ (x) dx = e iτ ρ (h)/h (9)

Here τ ρ (h) = h φ 2 + h π 4 -

x 0 -x 0 η ρ (y; h) dy + Const., where Const. is evaluated at the boundaries x = ±x E , and depends only on E ′ 1 . It will eventually disappear from the final formula, by adding to BS the contribution of the lower branch Λ <,ρ E . Note that x 0 -x 0 η ρ (y; h) dy, η ρ (y; h) being the derivative of the h ′ -depending phase function, is the semi-classical action.

Bohr-Sommerfeld quantization rules

We set F j,a,ρ ε,ω = i h [P a , χ a ] ρ U j,a,ρ) ε,ω , and similarly with a ′ . The set {G j,♭ ε,ω = F j,♭,+ ε,ω -F j,♭,- ε,ω : j ∈ {ν, -ν -1}, ♭ ∈ {a, a ′ }} (or their h-Fourier transform) can be interpreted as a basis of the microlocal co-kernel of P near a, a ′ . Following [10], we introduce Gram matrix G ρ of vectors U ρ 1 and U ρ 2 in this basis, namely

G = U 1 | G -ν-1,a ε,ω U 2 | G -ν-1,a ε,ω U 1 | G ν,a ′ ε,ω Û2 | G ν,a ′ ε,ω
. Using symmetries we get

G = G ρ = 2 1 e ρ 12 -d ρ 21 -1
The condition det(G (ρ) ) = 0 means that U 1 is colinear to U 2 , i.e. there is a global section of Ker h (P -E). Recall e ρ 12 = ρ d ρ

21

-1 ; for ρ = +1 (electronic state) we get Imd + 21 = 0, that is sin τ (+) (h) h = 0. We eventually obtain BS by "surgery": namely (ignoring tunneling) we cut and paste the half-bicharacteristic Λ >,+ E in the upper-half plane ξ > 0 with its symmetric part Λ <,- E in ξ < 0 and add together the contributions. By symmetry, the constant term Const. in τ + (h) drops out, while the other terms h φ 2 + h π 4 -

x 0 -x 0 η ρ (y; h) dy add up, which yields BS for the electronic state. We argue similarly for the hole state. This eventually gives Theorem 1.

Figure 1 .

 1 Figure 1. Phase-space picture

  see[START_REF] Whittaker | A Course of Modern Analysis[END_REF]. This normalization is called Whittaker normalization. Classically forbidden regions |η| > E lie on Stokes lines, classically allowed region |η| < E in between, and 3 Stokes lines stem from each "turning point" η = ±E.

d) Microlocal solutions.

Acknowledgements: We thank Timur Tudorovskiy for having introduced us to the problem. This work has been partially supported by the grant PRC CNRS/RFBR 2017-2019 No.1556 "Multidimensional semi-classical problems of Condensed Matter Physics and Quantum Dynamics".


A Anikin
email: anikin83@inbox.ru

S Dobrokhotov

V Nazaikinskii

& M Rouleux
email: rouleux@univ-tln.fr

SEMI-CLASSICAL GREEN FUNCTIONS

Let H(x, p) ∼ H 0 (x, p) + hH 1 (x, p) + • • • be a semi-classical Hamiltonian on T * R n , and Σ E = {H 0 (x, p) = E} a non critical energy surface. Consider f h a semi-classical distribution (the "source") microlocalized on a Lagrangian manifold Λ which intersects cleanly the flow-out Λ + of the Hamilton vector field X H 0 in Σ E . Using Maslov canonical operator, we look for a semi-classical distribution u h satisfying the limiting absorption principle and H w (x, hD x )u h = f h (semi-classical Green function). In this report, we elaborate (still at an early stage) on some results announced in [AnDoNaRo1] and provide some examples, in particular from the theory of wave beams.

Introduction

Our general motivation is to solve equations like (H -E)u h = f h , with f h (x) = h -n f x h , f ∈ L 2 (R n ), a "localized function" at x 0 = 0. In case H = -h 2 ∆, we apply h-Fourier transform

F h u h (ξ) = e -ixξ/h u h (x) dx to get F h u h (ξ) = F 1 f (ξ)
ξ 2 -E . If E < 0, this is an elliptic equation, and u h has the same form as f h . If E > 0 however, (1) u + (x; E, h) = (2πh) -n e ixξ/h F 1 f (ξ) ξ 2 -E -i0 dξ is defined through regularisation, so to satisfy the limiting absorption principle. Actually,

u + (x; E, h) = E + f h , E + = i h ∞ 0
e -it(H-E)/h dt is the forward parametrix. When f h (x) is replaced by δ(xx 0 ), u + is the Green function. Such a distribution is very singular since there are infinitely many classical trajectories linking the "source"

x 0 to x in time t; so it is suitable to consider instead F 1 f with compact support.

This problem has of course received considerable attention [Ke], [Ba], [Ku]; our special purpose is to reformulate it in terms of Maslov canonical operator, and describe the solution globally, including unfolding of Lagrangian singularities; this is of special importance in the context of wave propagation.

More generally, let f h = f (•; h) be a semi-classical Lagrangian distribution (or oscillatory integral) defined on the manifold M (for all practical purposes, we shall assume here M = R n ), i.e. locally f (x; h) = (2πh) -n/2 e iϕ(x,θ)/h a(x, θ; h) dθ where ϕ(x, θ) is a non-degenerate phase function in the sense of [Hö], and a(x, θ; h) = a 0 (x, θ) +

= | det ϕ ′′ | -1/2 |dξ| 1/2
. The (oscillating) principal symbol of f in Λ ϕ is then defined (up to the factor dµ ϕ ) as e iφ(ξ)/h A 0 (ξ) = e iφ(ξ)/h e iπ sgn ϕ ′′ /4 a 0 (x(ξ), θ(ξ))

where φ is a "reduced phase function". Conversely, assume ι : Λ → T * M is a smooth Lagrangian immersion, with a smooth positive half-density √ dµ, it can be parametrized locally by phase functions ϕ in canonical charts ι ϕ : U ϕ → T * M . These phases can be chosen coherently, and define a class of "reduced phase functions" φ, parametrizing ι locally. This gives the fibre bundle of phases L h , including Maslov indices, equipped with transition functions. We are also given local smooth halfdensities |dµ ϕ | 1/2 on Λ, defining the fibre bundle of half-densities Ω 1/2 , equipped with transition functions. The collection of these objects make a fibre bundle Ω 1/2 ⊗ L h over Λ. A section of Ω 1/2 ⊗ L h will be written as

f (x; h) = K h (Λ,µ) A (x; h) where K h (Λ,µ
) is called Maslov canonical operator. The "reduced phase function" and the "principal symbol" of f are defined invariantly. See [M], [Du], [Iv], [DoZh], [DNS] for details.

So our general aim is to construct, in term of Maslov canonical operators, a representation of E + f . Here are some examples of f (expressed in a single chart):

(1) Λ = {(x, ∂ x ϕ) : x ∈ M }. WKB functions f (x, h) = e iϕ(x)/h a(x, h) or in Fourier representation (2)
f (x; h) = e iπn/4 (2πh) n/2 e i(xp+S(p))/h A(p; h) dp (such an integral conveniently normalized including a phase factor, will be written

* (• • •); here sgn x • p) ′′ = -n).
(2) Gaussian functions f (x; h) = 1 h n exp(-ω 2 •x 2 2h ) or more general (superposition of) coherent states.

(3) A conormal distribution with Λ = T

* N R n , N = {x n = 0} (3) f (x; h) = * e ix n p n /h A(x ′ , p n ) dp n (4) Λ = {x = X(ϕ, ψ) = ϕω(ψ), p = P (ϕ, ψ) = ω(ψ), ϕ ∈ R}, ω ∈ S n-1 , f identifies with a
"Bessel beam", see Sect. 4.

Hypotheses and main result

Let Λ 0 → T * M , Λ 1 → T * M be smooth embedded Lagrangian manifolds, Λ 1 with a smooth boundary ∂Λ 1 (isotropic manifold). Following [MelUh], we say that (Λ 0 , Λ 1 ) is an intersecting pair iff Λ 0 ∩ Λ 1 = ∂Λ 1 and the intersection is clean, i.e.

∀ρ ∈ ∂Λ 1 T ρ Λ 0 ∩ T ρ Λ 1 = T ρ ∂Λ 1
In particular, near each ρ ∈ T * M , there is a canonical transform mapping Λ 0 to T * 0 R n = {x = 0} and Λ 1 to

Λ 0 + = {(x, ξ) : x = (0, x n ), ξ = (ξ ′ , 0), x n ≥ 0}
i.e. the flow-out of T * 0 R n by the Hamilton vector field X ξ n = ((0, • • • , 0, 1), 0) of ξ n , passing through x = 0, ξ 0 = (ξ ′ , 0).

Let now H(x, p; h) = H 0 (x, p) + hH 1 (x, p) + • • • be a symbol of class S 0 (m) and H w (x, hD x ; h) its Weyl quantization, as a bounded operator on L 2 (M ). Assume (A.1) The energy surface

Σ E = H -1 0 (E) is non critical. (A.2) The set L = Σ E ∩ Λ is compact. (A.
3) The Hamiltonian vector field X H 0 is transverse to Λ along L.

For t ≥ 0, let g t H 0 = exp tX H 0 , Λ t = g t H 0 (Λ), and Λ + = t≥0 g t H 0 (L). The pair (Λ, Λ + ) define a Lagrangian intersection, and L = ∂Λ + . We assume also (A.4) The Hamilton flow g t H 0 | L extends for all t > 0. (A.5) |π x (g t H 0 | L )| → ∞ as t → ∞ (non trapping condition). Let ψ be a system of coordinates on L, completed (locally near L) to coordinates (ψ, τ ) on Λ and to coordinates (ψ, t) on Λ + , so that the positive measure on Λ (square of half-density) takes the form dµ = |dτ ∧ dψ|. This implies the measures dµ t = (g t * H 0 ) -1 (dµ) on Λ t and dµ + = |dt ∧ dψ| on Λ + . We shall consider a partition of unity subordinated to a covering of the pair (Λ, Λ + ), with cut-off functions (of small enough support) χ 0 ∈ C ∞ 0 (Λ), that depends only on τ (the "radial variable" on Λ), with χ 0 (τ ) ≡ 1 near L = {τ = 0}, and χ 0 ∈ C ∞ 0 (Λ + ), that depends only on t, χ 0 (t) ≡ 1 near t = 0. Let also χ T ∈ C ∞ (R + ), χ T ≡ 1 for t ≤ T /2, χ T ≡ 0 for t ≥ T /2 where T is so large that x / ∈ π x (g t H 0 | L ) if t ≥ T (using the non-trapping condition (A.5)). For simplicity, we state our main result [START_REF] Anikin | Maslov's canonical operator on a pair of Lagrangian manifolds and asymptotic solutions of stationary equations with localized right-hand sides[END_REF] in the case Λ = T * 0 M , so that all solutions of Hamilton equations start from some (0, p), with p ∈ supp A. Representing E + as a sum of Maslov canonical operators associated with the pair (Λ, Λ + ) in particular we shall retrieve the theorem on propagation of singularities for wave-front sets WF h u ⊂ WF h (Hu) ∪ Char H.

Theorem 1: Let Ω ⊂ M be a bounded domain. Under hypotheses (A.1)-(A.5) above, there is u = u h , solution of H w (x, hD x ; h)u h = f h , satisfying the limiting absorption principle, which can generically be cast in the following form

(4) u(x, h) = * (1 -χ 0 (τ (p))A(p) H(0, p) -E e ipx/h dp + ε 0 0 dt * e iΦ(t,x,η)/h e -iΘ(t,x,η) | det J (t, x, η)| -1/2 χ T (t)χ 0 (t) χ 0 τ (η) A(η) dη dt + K h (Λ + ,dµ + ) χ T (•)(1 -χ 0 (•))e -iΘ A +
up to O(h) terms; the summands on the right hand side are boundary, transient and wave part, respectively. In Eq.( 4), Φ(t, x; η) solves Hamilton-Jacobi equation (see Sect.3), and

Θ(t, x, η) = t 0 H 1 • exp(sX H 0 (0, η))
ds is the integral of the sub-principal 1-form, where we recall H 1 is the sub-principal symbol of H. Moreover J (t, x, η) is a Jacobian computed from Hamilton-Jacobi equation.

The case when Λ is slightly "tilted" with respect to the vertical plane, i.e. has generating function φ(x, ξ) = xξ + S(ξ) can be formulated in a very similar way (see Sect.3). For more general f 's, the boundary part should be written as K h (Λ,µ)

(1-χ 0 )A H (x; h), and the transient part as the integral over t of a phase factor times K h (Λ t ,µ t ) χ T (•) χ 0 (•)χ 0 (•)A t (x; h), A t being the solution of the transport equation along the projection of exp tX H 0 with initial data A(η).

The boundary part is microlocalized on Λ, the wave part on Λ + . Let us make first a few comments:

1) (4) can be interpreted as "integrated" Van Vleck Formula [CdV], which expresses the semiclassical propagator e -itH/h acting on a localized function.

2) The wave part contributes generally at x ∈ Ω only if t → Φ(t, x, η) has a non-degenerate critical point for all η ∈ supp f : let hD x n be the "model" operator, and f (x; h) = * e ixp/h A(p) dp.

Then [MelUh]

u h (x) = i h ∞ 0 χ T (t) dt * e i(x ′ ξ ′ +(x n -t)ξ n )/h A(ξ) dξ
verifies the limiting absorption principle and hD

x n u(x, h) = f (x, h) + O(h ∞ ) for
x n ≤ T /2; still u h has no "wave part".

3) Formula ( 4) is also valid near a focal point x ∈ Ω, or more generally when x is linked to x 0 through a trajectory containing several focal points. The wave-part simplifies outside the focal points to a WKB form, involving non trivial Maslov indices passing the first focal point.

4) The phase function Φ(t, x, η) can be also replaced by a Lagrangian action, in the spirit of [DNS]. This will be discussed in detail in a future work.

We illustrate Theorem 1 by computing u h explicitely in the 2-D case for Helmholtz operator with constant coefficient as in (1), but f with compact support. Let f also be radially symmetric; its Fourier transform g = F 1 f is again of the form g(p) = g(|p|) = g(r) and extends holomorphically to

C 2 . For E = k 2 , k > 0, we rewrite (1) as u h (x) = u(x) = u 0 (x) + u 1 (x) with u 0 (x) = k + i ε (2πh) 2 2π 0 dθ ∞ 0 exp[i|x|r cos θ/h] g(r) r 2 -(k + i ε) 2 dr u 1 (x) = 1 (2πh) 2 2π 0 dθ ∞ 0 exp[i|x|r cos θ/h] g(r) r + k + i ε dr
To compute u 0 we use contour integrals. When θ ∈] -π 2 , π 2 [, we shift the contour of integration to the positive imaginary axis and get by the residues formula (5)

∞ 0 exp[i|x|r cos θ/h] g(r) r 2 -(k + i ε) 2 dr + ∞ 0 exp[-|x|r cos θ/h] g(ir) r 2 + (k + i ε) 2 idr = 2iπ g(k + i ε) 2(k + i ε) exp[i|x|(k + i ε) cos θ/h] while for θ ∈] π 2 , 3π 2 [, (6)
∞ 0 exp[i|x|r cos θ/h] g(r) r 2 -(k + i ε) 2 dr - ∞ 0 exp[|x|r cos θ/h] g(-ir) r 2 + (k + i ε) 2 idr = 0
Summing up ( 5) and ( 6), integrating over θ ∈]0, 2π[ and letting ε → 0, we obtain

u 0 (x) = iπg(k) (2πh) 2 π/2 -π/2 exp[i|x|k cos θ/h] dθ+ ∞ 0 dr r 2 + k 2 π/2 -π/2 g(ir) - 3π/2 π/2 g(-ir) exp[-|x|| cos θ|/h] dθ
Since g(ir) = g(-ir), the latter integral vanishes, so we end up with

u 0 (x) = iπg(k) (2πh) 2 π/2 -π/2 exp[i|x|k cos θ/h] dθ
It is readily seen that

WF h u 0 ⊂ {x = 0} ∪ {(x, k x |x| ), x = 0} = Λ ∪ Λ + Consider now u 1 . We let ε → 0 and set g(r) = g(r) r(r+k) . Since g(r) √ r ∈ L 1 (R + ), we have u 1 (x) = H 0 ( g)( |x| h )
, where H 0 denotes Hankel transform of order 0. Let χ ∈ C ∞ 0 (R 2 ) be radially symmetric, and equal to 1 near 0, since WF h f h = {x = 0}, we have

g = F h (χf h ) + O(h ∞ ) = (2πh) -2 F h (χ) * g + O(h ∞ )
so in the expression for u 1 we may replace mod O(h ∞ ), g(r) by a constant times g(r) = (F h (χ) * g)(r) r(r+m)

(see [Bad] for 2-D convolution and Fourier transform in polar coordinates). To estimate WF h u 1 , we compute again the Fourier transform of (1χ) g where χ is a cut-off equal to 1 near 0, and we find it is again O(h ∞ ) if χ ≡ 1 on supp χ. This shows that WF h u 1 ⊂ {x = 0}. See also [MelUh], Prop.

2.3. Note that the decomposition u 0 + u 1 is directly related with the corresponding one in Theorem 1 as the sum of the boundary (u 1 ), and the wave part (u 0 ). We conjecture that the transient part in Theorem 1 can be removed (taking a limit supp χ 0 → ∅) when the Hamiltonian flow X H 0 enjoys some non-degeneracy properties (see Sect.3), as in the case of a geodesic flow.

Maslov canonical operators associated with (Λ, Λ + )

Proof of Theorem 1 consists first in looking at the propagator e -itH/h acting on Lagrangian distribution f h as above, i.e. the solution of the Cauchy problem hD t +H w (x, hD x ) v = 0, v| t=0 = f h .

Next step is to integrate with respect to t after introducing the partition of unity above.

In this report, we shall content to construct the phase functions by solving Hamilton-Jacobi equation. First we recall from [Hö],Thm 6.4.5 the following:

Φ(y ′ , 0; η) = φ(y ′ ), ∂Φ ∂y (0; η) = η
We consider here the case of a "maximally singular" chart U for Λ where Λ = Λ φ has generating function φ(x, ξ) = S(ξ) + xξ. Even if Λ is a plane, the outgoing manifold Λ + may be very complicated far away from x 0 , but changing the canonical charts, we can proceed step by step.

Without loss of generality we can also assume here E = 0, S(0) = 0, ∂ ξ S(0) = 0. Let η, τ such that τ + H(0, η) = 0, by Theorem 2 (after slightly changing the notations), there exists Φ(x, t; η, τ ) such that ( 7)

∂Φ ∂t + H(x, ∂Φ ∂x ) = 0 Φ(x, 0; η, τ ) = S(η) + xη ∂Φ ∂x (0, 0; η, τ ) = η, ∂Φ ∂t (0, 0; η, τ ) = τ
Moreover, when τ is small enough, the intersecting pair (Λ, Λ + ) for energy level H = 0 extends to a smooth family of intersecting pairs (Λ, Λ + (τ )) for energy levels H = -τ

Λ + (τ ) = {(x, ξ) ∈ T * M, ∃t ≥ 0, ∃(y, η) ∈ Λ, (x, ξ) = exp tX H (y, η), τ + H(y, η) = 0}
intersecting along L(τ ) given by η n = η n (η ′ , τ ), with η n (ξ ′ 0 , 0) = 0. So near x = x 0 we can assume, possibly after permuting the ξ-coordinates, that H(0; ξ ′ 0 , 0) = 0, ∂H ∂ξ n (0; ξ ′ 0 , 0) = 0, and for τ small enough, the equation τ +H(x, ξ) = 0 is equivalent to ξ n = ξ n (x, ξ ′ , τ ), with ξ n (0, ξ ′ 0 , 0) = 0. Now we want to set τ = 0 and therefore, given (x, η), solve the equation ∂Φ ∂t = 0 for t ≥ 0. As we have seen, uniqueness of solutions doesn't always hold, as shows the "model case" hD x n . Namely, the phase function given by Hamilton-Jacobi theory is Φ(x, t; η) = φ(x ′ , x nt, η), so ∂ t Φ = 0 iff η n = 0, for all t. However the phase parametrizing Λ + is given for small t by Taylor expansion

Φ(x, t; η) = xη + S(η) -tH(x, η) + t 2 2 ∂ ξ H(x, η)∂ x H(x, η) + • • •
The following assumption in turn ensures existence of a finite number of solutions t = t(x; η) > 0 of ∂Φ ∂t = 0. Denote by Exp t x 0 η = π x (exp tX H 0 (x 0 , η)) the projection of the bicharacteristic of H starting from (x 0 , η) near (x 0 , ξ 0 ) at time 0. Assume:

(A ′ 0 ) ∃Ω ξ ⊂ R n open (small) neighborhood of ξ 0 , such that if x = Exp t
x 0 η for some t > 0, and η ∈ Ω ξ , then the map ξ → Exp t

x 0 ξ is a local diffeomorphism near η; in other terms, x 0 and x are not conjugated along any trajectory that links them together within time t, with initial momentum ξ.

The set of such x is an open set Ω x .

Proposition 1: Under hypothesis (A ′ 0 ), for all (x, η) ∈ Ω x × Ω ξ there is a finite number of t j = t j (x, η) > 0 (1 ≤ j ≤ N ) solutions of ∂ t Φ(x, t; η) = 0 and t j are non degenerate critical points.

Moreover (A ′

0 ) is generically fulfilled before occurence of the first focal point, as show the following examples.

Let H 0 be a geodesic flow, E > 0, f conormal to N = {x 0 } ; Ω x is a small enough ball centered at x 0 , (geodesic ball), or Ω x is a neighborhood of a minimal geodesic {y = Exp s x 0 η, 0 ≤ s ≤ t 0 } for η in some small neighborhood of η 0 . This applies when f is as in (2) and supp A is localized near η 0 .

The same holds with f as in (3) conormal to the hypersurface N = {x n = 0} (or more generally to a surface N of positive codimension in M ). For instance when N = {x n = 0}, (A ′ 0 ) holds with Ω ξ an open (small) neighborhood of η = (0, • • • , 1), and Ω x a neighborhood of a minimal geodesic from N to some x 1 . This follows from a well-known property of minimal geodesics (see e.g. [HeSj], Proposition 6.3). In case of a Schrödinger operator with principal symbol H 0 (x, ξ) = ξ 2 + V (x) -E (where E is a scattering energy) we use Maupertuis-Jacobi principle to reduce again to a Riemannian metric d E given by ds = (E -V (x)) 1/2 |dx|. So we arrive at the same conclusions, the metric d E being conformal to the standard metric. We proved the following:

Proposition 2: Under the minimality assumptions above, (A ′ 0 ) holds. For η ∈ Ω ξ and x = Exp t

x 0 (η) ∈ Ω x , 0 < t < t 0 , the Lagrangian manifold Λ t = exp tX H (Λ φ ) has same rank as Λ, and is of the form Λ Φ(t,•) . In particular

exp tX H 0 (Λ φ ) = {(x, ∂ x Φ(t, x; η)) : ∂ η Φ(t, x; η) = 0} ⊂ {exp tX H 0 (x, η), (x, η) ∈ T * N ⊥ }
Moreover for all (x, η) ∈ Ω x × Ω ξ there is a unique t(x, η) > 0 solution of ∂ t Φ(x, t; η) = 0 and t(x, η) is a non degenerate crtical point.

Extending the geodesic γ beyond the first focal point occuring in some Λ t for some t = t 0 . we only need another representation of the phase function. Using the mixed representation for Lagrangian manifolds, we know that for any k = 0, • • • , n, (k = 0 corresponds to a maximally singular chart, k = n to a regular chart, they are mapped onto each other by Fourier transform), there exists a partition of variables x = (x ′ , x ′′ ) ∈ R k ×R n-k and ξ = (ξ ′ , ξ ′′ ), such that if π : R 2n → R n , (x, ξ) → (x ′ , ξ ′′ ), then rank d π = n. In such a chart, the generating function for Λ t takes the form φ(x, η) = x ′′ η ′′ +S(x ′ , η ′′ ).

We can reformulate Hamilton-Jacobi equations as in (7) in these coordinates, which has again a unique solution for small tt 0 . Generically (i.e. under an assumption (A ′ k ) modeled after (A ′ 0 ) in the new (x, ξ) coordinates), one still obtains a non degenerate phase function, which contributes to the fibre bundle L h over Λ + together with Maslov index of γ. More explicit formulae will be given elswewhere.

Of course, γ may contribute a finite number of times in the expression for the Green function at x, and a finite number of γ contribute to the wave part of the phase function.

Using eikonal coordinates.

The computations above can be simplified using special coordinates adapted to Λ, called eikonal coordinates.

Let ι : Λ → T * M be a smooth embedded Lagrangian manifold. The 1-form p dx is closed on Λ, so locally exact, and p dx = dS on any simply connected domain U . Such a S is called an eikonal (or action) and is defined up to a constant. Assume dS = 0 on Λ, then S can be chosen as a coordinate on U , which we complete by smooth functions ψ ∈ R n-1 .

We use here eikonal coordinates to construct a phase function solving Hamilton-Jacobi equation in case of a positively homogeneous Hamiltonian of degree m with respect to p. To fix ideas, we take P (s, ψ), Ẋ(s, ψ) = P (s, ψ), ∂ p H(X, P = mH(X, P ) = m and S(t, ψ) = S 0 + mt is the action on Λ + . So m dt = P (t, ψ), dX(t, ψ) = P (t, ψ), Ẋ(t, ψ) dt + P (t, ψ), X ψ (t, ψ) dψ

It follows that P (t, ψ), Ẋ(t, ψ) = m and P (t, ψ), ∂ ψ X(t, ψ) = 0. Now we complete the coordinate system ψ on Λ by a smooth function r such that L is given by r = 1, and set Φ(x, (t, ψ, r)) = mt + r P (t, ψ), x -X(t, ψ)

where r can be interpreted as a Lagrange multiplier. Let us check that Φ satisfies Hamilton-Jacobi equation. We have (10)

∂ t Φ = Φ = m + r Ṗ , x -X(t, ψ)r P, Ẋ = m(1r) + r Ṗ (t, ψ), x -X(t, ψ)

∂ r Φ = P (t, ψ), x -X(t, ψ)

∂ ψ Φ = r ∂ ψ P (t, ψ), x -X(t, ψ)

Last 2 equations in ( Φ, ∂ r Φ, ∂ ψ Φ) = 0 give an homogeneous linear system with determinant det(P, P ψ ).

On x = 0 we get |p| = n(0) > 0, and in dimension n = 2, p = n(0) t (cos ψ, sin ψ). It follows that det (P, P ψ ) = |n(0)| 2 so for small t, we get x -X(t, ψ) 0, so the phase is critical with respect to (ψ, r) for x = X(t, ψ). Substituting into the last equation ( 10) we get Φ = 0 when r = 1. We complete the proof that Φ is a generating function for Λ + by checking d∂ t φ, d∂ ψ Φ, d∂ r Φ are linearly independent on the set x = X(t, ψ). This can be done by examining the variational system associated with Hamilton equations.

Moreover we can reduce this generating function by eliminating t by stationary phase. Again, of course, this holds only for small t > 0, before unfolding of Lagrangian singularities.

  ha 1 (x, θ) + • • • an amplitude. With f we associate the critical set C ϕ = {(x, θ) ∈ M × R N : ∂ θ ϕ = 0} and ι ϕ : C ϕ → T * M with image the Lagrangian submanifold Λ ϕ = {(x, ∂ x ϕ) : (x, θ) ∈ C ϕ }. Choose local coordinates ξ ∈ R d on Λ ϕ and define the half-density in the local chart (C ϕ , ι ϕ ) by dµ ϕ

Theorem 2 :

 2 Denote the variable in T * R d by (y, η) = (y ′ , y d ; η ′ , η d ). Let H be a real valued, smooth Hamiltonian near (0, η) such that H(0, η) = 0, ∂ η d H(0, η) = 0, and let φ be a real valued, smooth function on R d-1 such that ∂ y ′ φ = η ′ . Then there exists in a neighborhood of 0 ∈ R d a unique real valued solution Φ(y; η) of Hamilton-Jacobi equation H(y, ∂ y Φ) = 0 satisfying the boundary condition

P

  Λ = {x = 0} intersects the energy surface H = 1 along L. Let us compute the eikonal S on Λ + . Integrating Hamilton equations we have x = X(t, ψ), p = P (t, ψ) where ψ ∈ R 2 , hence dS = p dx| Λ + = P (t, ψ), dX(t, ψ) . Since dx = 0 on Λ, we have S(0, ψ) = Const. = S 0 (s, ψ), Ẋ(s, ψ) ds By Hamilton equations and Euler identity, we have on H = 1.

Acknowledgements: This work was supported by Grant PRC No. 1556 CNRS-RFBR 2017-2019 "Multi-dimensional semi-classical problems of Condensed Matter Physics and Quantum Mechanics", and RFBR Grant No. 17-51-150006.

Example 2: We take H as in (8) with n(x) = n(|x|), n = 2 for simplicity and Λ = {x = X(ϕ, ψ) = ϕω(ψ), p = P (ϕ, ψ) = ω(ψ), ϕ ∈ R, ω ∈ S n-1 } ψ being the usual angles parametrizing ω ∈ S n-1 . When n = 2 this is the wave-front set of Bessel function f h (x; h) = J 0 ( |x| h ); such functions arise in the wave beam theory (see [Ki], [DoMaNa] and references therein), so we call f h a "Bessel beam".

Computing the action we find p dx| Λ = dϕ so coordinate ϕ will play the role of x in he previous Example. We have n(|x|) = n(ϕ), and ϕ = Const. on L = Λ ∩ {H = 1}. The argument above extends readily to this setting, in particular as in ( 9)

Now we can apply the results of [DNS], Sect.2.2 on eikonal coordinates for a Lagrangian manifold in a general position, with X ψ ( ϕ, ψ) of rank k. Here ϕ, ψ are coordinates in a local chart of Λ + in the extended phase space. More specifically, we take ψ = (t, ψ), ϕ = ϕ, try to make a change of variables ϕ = ϕ(x, t, ψ), and seek for a generating function of Λ + in the form Φ(x, (t, ψ, r)) as in the previous Example. Details will be given elsewhere.

More examples and perspectives.

The methods above apply in a number of situations as:

(1) The water-wave Hamiltonian H(x, hD x ; h) with H 0 (x, p) = |p| tanh |p|D(x) -E, together with Λ = {x = x 0 } has been discussed in [DoNa], [START_REF] Anikin | Asymptotics of Green function for the linear waves equations in a domain with a non-uniform bottom[END_REF], in relationship with (11) or Helmhotz operator H(x, hD x ) = -h 2 ∆n(x)

2 . The localized function f can be a Gaussian (even in x) or a Gaussian times a linear function (odd in x), or can be of antenna type, i.e. its Fourier transform localized in a cone in p.

(2) The kinetic part of Hamiltonian is of Lorenzian type (as -p 2 0 + p 2 1 ), and f a localized function (Gaussian) supported on Λ = {x = x 0 }, so that the semi-classical Green function u h is the linear response to f localized on Kelvin angle (or Mach cone).

It should also be possible to construct semi-classical Green functions in case the pair (Λ, Λ + ) is no longer intersecting cleanly, but glancing. We then need second-microlocalization, and introducing so called 2-phases, see e.g. [LaWi] in the standard (polyhomogeneous) calculus.


Brieux Durand

Aurélie Lecestre

Philippe Menini

Guilhem Larrieu

Miniaturized 3D gas sensors based on silicon nanowires for ppb range detection

I. INTRODUCTION

Current trends go towards embedded multifunction systems. The most notable case is the evolution of mobile phones which resulted the smartphone: a "tool" with many uses where new sensors are integrated to add new functionalities. In that context, gas sensors would be implemented in order to monitor, during all day long, the air quality in polluted urban environments as the concentration of nitrogen dioxide (NO2), or detect domestic gas leaks. However this requires a reorientation of the sensor industry towards components more and more miniaturized to reduce weight, size, power consumption and price. Currently commercial sensors use a metal oxide as sensitive layer on a micro hotplate technology. These devices are relatively large (1-10 mm²) and require high power consumption for supplying the heater (mW range). In addition, they demonstrate low selectivity and sensitivity to high concentrations of gas (>1 ppm), out of step with the recommendations of the various environmental agencies to long period exposure thresholds (50 ppb NO2 [START_REF]Our Nation's Air[END_REF]). In that context, the integration of 1D nanostructures, thanks to their high ratio surface/volume, have been proposed to improve the sensitivity and detection limit [START_REF] Cui | Functional nanoscale electronic devices assembled usingsilicon nanowire building blocks[END_REF], [START_REF] Chen | Nanowire-based gas sensors[END_REF]. However, despite numerous studies of nanowire devices, any demonstration presents satisfactory operation [START_REF] Mcalpine | Highly ordered nanowire arrays on plastic substrates for ultrasensitive flexible chemical sensors[END_REF], [START_REF] Zhang | Detection of NO 2 down to ppb Levels Using Individual and Multiple In 2 O 3 Nanowire Devices[END_REF] at very low concentrations. Here, we propose a 3D sensor architecture based on silicon nanowires networks of (Figure 1). The 3D configuration improves the free surface for gas adsorption in order to enhence the sensitivity, and decrease the detection limit.

II. DEVICE FABRICATION

This approach is compatible with silicon processing using conventional microelectronics tools, allowing the potential co-integration of these sensors with CMOS devices. The transduction is based on resistor configuration where the nanowires bridge the positive and negative electrodes to allow current flow. The device conductance change is monitored. The fabrication process is large scale, reproducible and straightforward, and can be summarized in 3 steps: (i) firstly the nanowire networks patterning, (ii) secondly the bottom electrical contacts fabrication, and (iii) finally the top contacts structuration. A top-down approach has been developed (Figure 2.a-e) to perfectly control all the parameters of NW networks (diameter, height, spacing and the numbers of nanowires…). An approach combining photolithography, plasma etching and sacrificial oxidation is developed to pattern the vertical nanowire networks. The first step consists of substrate patterning by photolithography (i-line stepper), to achieve array of resist nanodots (Figure 2.a-c, f). Then the patterns are transferred into the bulk silicon by deep reactive ionic etching (DRIE) (Figure 2.d). Finally remaining resist is stripped and a sacrificial oxide layer is grown to tune precisely the nanowire diameter (Figure 2.e, g).

The device is composed of two symmetrical aluminum contacts at each end of the NWs. The top contact is obtained by an "air bridge" approach of allowing the definition of a three-dimensional connection to the NW top (Figure 3). For the top and bottom contacts, a standard etch-back process is used. However for the top contact, a sacrificial material is employed to connect only the top of the nanowires, leaving the rest of the surface available to react with gases. Under controlled atmosphere, the device demonstrates excellent gas sensing performances. A typical gas response curve is presented in Figure 4. When NO2 is injected, the resistance of the sensor decreases, but when the gas injection stops, the resistance returns to the base-line. The sensor recovery occurs naturally at room temperature, and does not require specific conditions such as UV illumination [START_REF] Mcalpine | Highly ordered nanowire arrays on plastic substrates for ultrasensitive flexible chemical sensors[END_REF] or flushing [START_REF] Zhang | Detection of NO 2 down to ppb Levels Using Individual and Multiple In 2 O 3 Nanowire Devices[END_REF]. For ammonia, an opposite response was observed: the resistance increases, contrary to NO2 case. Moreover the response of carbon monoxide (CO) and propane (C3H8) have been investigated: in both case there are no response of the sensor. This demonstrates the selectivity of the device, by discriminating NO2 and NH3, without being interfered with CO and C3H8.

Due to the enhanced adsorption surface conferring by the 3D configuration, extremely low sensibility level has been reached. A very high response (30%) is obtained at 50 ppb of NO2 (the dilution limit of our test bench), in comparison of the 25% reached for 200 ppb in [START_REF] Mcalpine | Highly ordered nanowire arrays on plastic substrates for ultrasensitive flexible chemical sensors[END_REF], indicating that the limit of NO2 detection by our device is greatly below 20 ppb. In addition, the sensors retain their outstanding performance in real working condition. In particular, we demonstrated that they are not disturbed by the presence of moisture (tested up to 70% of relative humidity), one of the major drawbacks of metal oxide-based gas sensors. The device is sensitive to humidity and NH3, but in presence of NO2 (even for low concentration) the detection mechanism of NO2 dominates the others mechanisms (moisture and interfering gas).

IV. CONCLUSION

In summary, we demonstrated record detection of NO2 down to ppb level, using resistor gas sensors based on 3D SiNW networks. This device combines for the first time high sensitivity, selectivity, reversibility, low-power consumption, no poisoning by humidity, reliability, miniaturization, cost effective, large scale and CMOS compatible fabrication, and open the way of new possibilities for integration of unidimensional nanostructures for very sensitive embedded sensors.

Figure 1 :

 1 Figure 1: Illustration of the 3D sensor architectures implemented on vertical NW arrays.

Figure 2 :

 2 Figure 2: a)-e) Top-down fabrication of vertical nanowires network, SEM pictures of f) a resist nanodot and g) nanowires

Figure 3 :

 3 Figure 3: SEM picture of the sensor : a side view of one array of vertical silicon nanowires, contacted electrically by an air-bridge structure

Figure 4 :

 4 Figure 4: Response of the device to 50 ppb of NO2 (blue) and 5 ppm to NH3 (orange)


N Desmars

Y Perignon

G Ducrozet

C.-A Guérin

S T Grilli

P Ferrant

PHASE-RESOLVED PREDICTION OF NONLINEAR OCEAN WAVES FROM REMOTE OPTICAL MEASUREMENTS PRÉDICTION À PHASE RÉSOLUE DE CHAMP DE VAGUES NON-LINÉAIRE PAR TÉLÉDÉTECTION OPTIQUE

We investigate the performance of a phase-resolved algorithm for the prediction of nonlinear ocean wave fields based on spatio-temporal remote optical measurements of free surface elevations. To quantify the prediction accuracy, error estimates are calculated for the case of a unidirectional wave field by ensemble averaging a large number of synthetic data sets. For several characteristic wave steepnesses, we compare results based on linear and weakly nonlinear approaches. It is shown that our nonlinear wave model, based on a Lagrangian description of the free-surface dynamics, allows to catch relevant nonlinear effects that play an increasingly large role as the waves get steep, i.e., wave shape asymmetry and phase shift. Experimental data are also used to illustrate the performance of the method applied on real field measurements.

Résumé

Les performances d'une méthode de prédiction de champs de vagues non-linéaires à partir de mesures optiques spatio-temporelles d'élévation de surface sont étudiées. Afin de quantifier la qualité des prédictions, des indicateurs d'erreur sont estimés dans le cas d'une houle unidirectionnelle en moyennant les résultats sur un grand nombre de données synthétiques. Pour des quelques cambrures caractéristiques, nous comparons les résultats issus d'approches linéaire et faiblement non-linéaire. Il est montré que le modèle non-linéaire utilisé, basé sur une description Lagrangienne de la dynamique de surface libre, permet de capter des effets non-linéaires jouant un rôle important et grandissant avec la cambrure, i.e., asymmétrie de forme et déphasage. Des données expérimentales sont également utilisées pour illustrer les performances de la méthode appliquées à des mesures réelles.

I -Introduction: Formulation of the Prediction Problem

To many offshore applications, such as the optimal maneuvering and operational guidance of surface vessels or active control of ocean renewable energy harvesting systems, the capacity to predict in real-time the incoming phase-resolved wave field is an crucial issue. This problem was studied in some earlier work on the basis of free surface elevation time series measured at fixed wave probes [START_REF] Blondel | Deterministic non-linear wave prediction using probe data[END_REF][START_REF] Simanesew | Surface wave predictions in weakly nonlinear directional seas[END_REF][START_REF] Takagi | Prediction of wave time-history using multipoint measurements[END_REF]. However, in situ measurements, typically made from a moving vessel or vehicle, are more challenging since wave reconstruction must rely on data acquired at constantly updated locations surrounding the path of the vehicle, which leads to practical limitations. Solutions to this problem was proposed based on X-band radar measurements made from an on-board sensor, combined with a 3D-FT to reconstruct a large patch of free surface elevation surrounding the sensor [START_REF] Dankert | Wind-and wave-field measurements using marine x-band radar-image sequences[END_REF][START_REF] Hilmer | Deterministic wave predictions from the wamos ii[END_REF][START_REF] Qi | Phase-resolved wave field simulation calibration of sea surface reconstruction using noncoherent marine radar[END_REF][START_REF] Xiao | Large-scale deterministic predictions of nonlinear ocean wave-fields[END_REF]. This initial estimate was then used in a direct numerical simulation of future sea states.

Here, we present recent developments of a method for reconstructing and predicting ocean waves in the short term, ahead of a moving or fixed ocean structure, based on optical measurements such as made with a LIDAR camera mounted on the structure.

The general wave prediction problem is formulated in the form of explicit relationships between the main parameters of the problem, i.e., the incoming wave field, the targeted prediction horizon (time span during which the wave field dynamics can be estimated), and the kinematics and size of the structure of interest (e.g., a moving ship, a marine energy device).

Since typical wave fields are bounded in frequency and direction of propagation (i.e., the energy is not homogeneously spread out over all frequencies/directions), a basic statistical knowledge of the incoming sea state can help define relevant bandwidths for accurately representing the underlying deterministic fields. In the following, the wave field frequency and directional energy density spectrum S (k, θ) is assumed to be known, where k(ω) is the frequency-dependent wavenumber, with ω and θ the wave angular frequency and direction of propagation, respectively. The prediction horizon corresponds to the time constraint of the prediction problem, i.e., the duration for the structure control system or operator to adequately react to changes in incoming waves. The ocean structure motion defines its path (assumed to be known), which together with the structure size allow calculating the spatio-temporal area over which wave dynamics is to be predicted.

A central element in the wave prediction problem, relating the aforementioned parameters, is the concept of accessible prediction zone, that is the spatial area denoted here by P (t), over which sufficient information on the wave field is available to adequately predict its dynamics. This is detailed in the next section.

Accordingly, the wave prediction problem is formulated as follows: based on adequately located real-time ocean surface measurements, extract and process sufficient data to predict the wave dynamics at the structure's future position, fast enough to satisfy the prediction horizon constraint.

II -Application to Optical Measurements

II -1 Accessible Prediction Zone from Nonuniform Observations

Earlier work has shown that the spatio-temporal region over which wave dynamics can be predicted based on a set of free surface measurements, is bounded [START_REF] Naaijen | Limits to the extent of the spatio-temporal domain for deterministic wave prediction[END_REF][START_REF] Qi | Predictable zone for phase-resolved reconstruction and forecast of irregular waves[END_REF][START_REF] Wu | Direct Simulation and Deterministic Prediction of Large-scale Nonlinear Ocean Wave-field[END_REF]. When measurements are made using an optical system, at a specific sampling rate and over a given observation zone, this limits the amount of data that can be assimilated and used in the wave reconstruction process, yielding a reconstructed surface in space/time defined with finite frequency and direction bandwidths. In light of this, the sea-state prediction obtained by propagating the assimilated information is similarly limited to a spatio-temporal region referred to as prediction zone. In the following, we show how the latter can be defined for a set of nonuniform observations made in a unidirectional wave field (i.e., irregularly distributed in space). Both theoretical and experimental studies have shown that directional components in a dispersive wave field, of given amplitude, frequency, and phase, travel at the associated group velocity along their direction of propagation [START_REF] Naaijen | Limits to the extent of the spatio-temporal domain for deterministic wave prediction[END_REF][START_REF] Qi | Predictable zone for phase-resolved reconstruction and forecast of irregular waves[END_REF]. The intersections of the slowest and fastest components thus determine the boundary of the spatio-temporal region over which information is available and a prediction can be issue. Consequently, as time increases, the accessible prediction zone shrinks, to eventually disappear when the assimilated information is completely dispersed in space. Fig. 1 illustrates this phenomenon for a unidirectional wave field propagating in the x direction. The last time used for the assimilation corresponds to the reconstruction time t r . If only spatial data is used in the assimilation (Fig. 1a), the initial prediction zone P (t r ) is the spatial area where observations were made. However, when spatio-temporal data sets are used (over an assimilation time T a ; Fig. 1b), the initial prediction zone expands due to the advection of wave information during T a . A point (x, t ≥ t r ) is included in the prediction zone if,

x b o + c g1 (t -t r ) ≤ x ≤ x e o + c g2 (t -t r ) , (1)
where c g1 and c g2 are the fastest and slowest group velocities, respectively, and x b o and x e o define the beginning and the end of the initial prediction zone as,

   x b o = min k x min ok + c g1 (t r -t k ) , x e o = max k {x max ok + c g2 (t r -t k )} , (2)
where k ∈ {1, ..., K} with K the number of observation times.

II -2 Bandwidths of the Reconstructed Wave Field

The dynamics of a wave field is governed by the propagation of its most energetic components. Since ocean wave generation processes eventually yield energy spectra with fairly narrow frequency bandwidth in their main direction of propagation, the dynamics of the wave field can be accurately represented by finite frequency and direction bandwidths having appropriate cutoff limits k min,max , θ min,max . As mentioned above, the spatio-temporal characteristics of the observation grid constraint the wave information that is accessible. The smallest wavenumber measurable in this grid

k min = 2π/L o is function of the largest distance L o = x e o -x b
o between two observation points at reconstruction time t r (see Fig. 1b). However, L o itself depends on the minimum and maximum group velocities, thus on the cutoff frequencies. This paradox is resolved by conservatively calculating

k min = 2π/L c , with L c = x e o -x min oK ≤ L o .
L c only depends on the spatio-temporal location of the observations and on the minimum group velocity, which is related to the high cutoff wavenumber k max . When reconstructing a signal over a regular observation grid (i.e., with a constant spatial sampling), the maximum high cutoff frequency must satisfy Shannon's condition k max ≤ 2π/ (2 o ) where o is the distance between two observation points. Since the observation grid is highly irregular when using an optical method [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF], k max is set so that the spectral energy truncated at higher frequencies is negligible for the dynamic description of the wave field. Also, since optical measurements provide observations at many spatial locations, there is no constraints on the cutoff directions from the observation grid. Hence, likewise k max , θ min,max are calculated based on the truncation of a negligible amount of energy. If an estimate of the underlying wave field is known, one way to choose these parameters is to evaluate their influence on the quality of the reconstructed or predicted wave fields, since a convergence of the quality of the reconstruction/prediction is expected as k max increases and θ min , θ max gets larger. In operational conditions, such information is not typically available and one must rely on an estimate of the wave spectrum, for instance based on measurements of the motion of the structure of interest, or from earlier measurements of free surface elevations. In this case, relevant values of k max and θ min,max can be selected by ensuring that,

+∞ k max θ max θ min S (k, θ) dθ dk ≤ µ +∞ 0 2π 0 S (k, θ) dθ dk, (3)
where S (k, θ) is the directional wave energy density spectrum and µ 1 is the fraction of total energy that can be considered as negligible on free surface dynamics.

II -3 Group Velocities for the Determination of the Prediction Zone

In practice, the cutoff frequencies calculated for the wave field reconstruction can be too restrictive for the determination of the evolution of the prediction zone. Hence, the wave group velocities c g1,2 that govern the evolution of the prediction zone boundaries are instead defined based on two wavenumbers k 1,2 , such that,

F (k 1 ) = F (k 2 ) = µF (k p ) , (4)
where

F (k) = 2π 0 S (k, θ
) dθ is the one-sided wave spectrum, k p is the wavenumber based on its peak wavelength, i.e., that of maximum energy, and µ is a small fraction, here, of the peak spectral energy. k 1 (k 2 ) is the smallest (largest) wavenumber to respect condition (4). In the following, we will use the linear deep water dispersion relationship to estimate the group velocities from k 1,2 .

III -Wave Models

Computing real-time wave predictions using complex, highly nonlinear, wave models, such as based on High-Order Spectral method, is challenging and may be prohibitive. Hence, similar to earlier work [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF], we seek to apply a model that properly accounts for nonlinear effects, which may significantly affect wave propagation in strong sea-states, while being suitable for the real-time constraint of the prediction problem. In this work, we consider and compare two wave models: a simple model based on linear wave theory (LWT), and a weakly nonlinear wave model, denoted as Improved Choppy Wave Model (ICWM) [START_REF] Guérin | An improved lagrangian model for the time evolution of nonlinear surface waves. Under review for publication[END_REF], derived from a Lagrangian approach of the description of the fluid motion, as an improved version of the earlier CWM [START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF][START_REF] Nouguier | Scattering from nonlinear gravity waves: the "choppy wave" model[END_REF]. When applied to a set of observations, both models provide an analytic expression of the reconstructed free surface elevations, which allows for a very efficient numerical prediction of wave propagation, by simply increasing time in the formulation.

III -1 Linear Wave Theory

Linear Wave Theory (LWT) refers to the Stokes wave model derived based on the classical Eulerian approach and linearized with respect to wave steepness. With this model, the ocean surface is represented in a Cartesian coordinate system (x, y, z) (with x and y axes located on the mean water surface and the z-axis being vertical and positive upward) as the superposition of N individual periodic waves of amplitudes A n , wavelengths λ n , and propagating in directions θ n with respect to the x-axis, as,

η (r, t) = N n=1 A n cos (k n • r -ω n t + ϕ n ) , (5)
with r = (x, y) the horizontal position vector, t the time, ϕ n a set of random phases uniformly distributed in [0, 2π], and k n = k n kn = (k n cos θ n , k n sin θ n ) the wavenumber vectors, with k n = 2π/λ n . The latter are related to angular frequencies ω n through the deep water linear dispersion relationship, ω 2 n = gk n , where g denotes the acceleration of gravity. To simplify the mathematical formalism of wave field reconstruction methods developed hereafter, Eq. ( 5) is recast as,

η (r, t) = N n=1 a n cos ψ n + b n sin ψ n , (6)
where ψ n = k n • r -ω n t are spatio-temporal phases, and (a n , b n ) = (A n cos ϕ n , A n sin ϕ n ) are wave parameters describing the ocean surface.

III -2 Improved Choppy Wave Model

At each data acquisition time, optical measurements of ocean waves provide surface elevations measured at a set of irregularly distributed points of known coordinates in the reference frame. Earlier ocean surface reconstruction and prediction algorithms were implemented to use such irregular data sets [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF][START_REF] Wu | Direct Simulation and Deterministic Prediction of Large-scale Nonlinear Ocean Wave-field[END_REF]. Moreover, Nouguier et al's [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF] algorithm, which is extended here, was based on the Lagrangian CWM [START_REF] Gjosund | A lagrangian model for irregular waves and wave kinematics[END_REF][START_REF] Nouguier | choppy wave" model for nonlinear gravity waves[END_REF][START_REF] Nouguier | Scattering from nonlinear gravity waves: the "choppy wave" model[END_REF], which can simulate higher-order nonlinear wave properties, in terms of steepness, than Eulerian models of the same order. In the following the CWM equations are recast as the solution of Lagrangian dynamical equations in an Eulerian system.

Considering a water particle located on the ocean surface at coordinates ξ in the horizontal plane, its horizontal and vertical displacements are found as a Lagrangian perturbation expansion in wave steepness as [START_REF] Nouguier | Second-order lagrangian description of tridimensional gravity wave interactions[END_REF][START_REF] Pierson | Models of random seas based on the lagrangian equations of motion[END_REF][START_REF] Pierson | Perturbation analysis of the navier-stokes equations in lagrangian form with selected linear solutions[END_REF],

r (ξ, t) = ξ + i D i (ξ, t) , z (ξ, t) = i Z i (ξ, t) , (7)
respectively, where D i and Z i are the horizontal and vertical particle displacements with respect to the particle position at rest ξ associated with the i th -order of expansion, respectively. Noting an inconsistency in the Stokes drift predicted by the CWM, Guérin et al. [START_REF] Guérin | An improved lagrangian model for the time evolution of nonlinear surface waves. Under review for publication[END_REF] added a correction to these equations, yielding the Improved Choppy Wave Model (ICWM),

   r (ξ, t) = ξ + N n=1 kn -a n sin ψn + b n cos ψn + U s t, z (ξ, t) = N n=1 a n cos ψn + b n sin ψn + 1 2 a 2 n + b 2 n k n , (8)
where

ψn = k n • r -ωn t, ωn = ω n -1/2k n • U s and U s = N n=1 a 2 n + b 2 n ω n k n .
To use ICWM in practical simulations, its Lagrangian formulation [START_REF] Grilli | Ocean wave reconstruction algorithms based on spatio-temporal data acquired by a flash lidar camera[END_REF] is first transformed into an equivalent Eulerian model, by deriving an approximate explicit relationship between r and z. First, the horizontal particle shift U s t is implicitly incorporated into a modified angular frequency, which yields,

R (ξ, t) = r (ξ -U s t, t) = ξ + D (ξ) = ξ + N n=1 kn (-a n sin ψ Ln + b n cos ψ Ln ) , Z (ξ, t) = z (ξ -U s t, t) = N n=1 a n cos ψ Ln + b n sin ψ Ln + 1 2 a 2 n + b 2 n k n , (9)
where ψ Ln = k n • r -ω Ln t and ω Ln = ω n + 1/2k n • U s . Then, the particle vertical displacement D is calculated at the current particle location R rather than at its reference location ξ, as,

Z (ξ) = Z (R -D (ξ)) Z (R -D (R)) = η (R) . (10)
Accordingly, the free surface elevation at any location r reads,

η (r, t) = N n=1 a n cos Ψ n + b n sin Ψ n + 1 2 a 2 n + b 2 n k n , with, (11)
Ψ n = k n • r - N n=1 kn (-a n sin ψ Ln + b n cos ψ Ln ) -ω Ln t.
Earlier work has shown [START_REF] Grilli | Ocean wave reconstruction algorithms based on spatio-temporal data acquired by a flash lidar camera[END_REF][START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF] that errors due to the approximation made in Eq. ( 10) are on the order of the mean square slope of the ocean surface, which is expected to be small compared to other sources of error in the reconstruction process.

III -3 Ocean Surface Reconstruction

Similar to earlier work [START_REF] Grilli | Ocean wave reconstruction algorithms based on spatio-temporal data acquired by a flash lidar camera[END_REF][START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF], the ocean surface is reconstructed, using either LWT or ICWM, based on determining model parameters that minimize the mean square error between spatiotemporal ocean observations and their representation in the model, referred to as the Cost Function, i.e.,

C (p) = K k=1 J j=1 η jk (p) -η jk 2 = L =1 (η (p) -η ) 2 , (12)
in which p = {a n , b n } (n ∈ {1, ..., N }) is the control vector containing the 2N unknown model parameters, J and K are the number of spatial observations made at each observation time and the number of observation times, respectively (hence, the total number of spatio-temporal observations assimilated in the model is, J × K = L), η jk are free surface elevations measured at spatial locations r j (j ∈ {1, ..., J}) and times t k (k ∈ {1, ..., K}), η jk are estimates of these made using the wave model, i.e., with Eq. ( 6) for LWT or Eq. ( 11) for ICWM.

III -3.1 Linear Wave Field

The ocean surface is reconstructed through a least square minimization of the cost function, as,

∂C ∂a m = 0, ∂C ∂b m = 0 ⇐⇒ A mn p n = b m , (13)
for n, m ∈ {1, ..., N } (note that index summation is implied for repeated indices), with,

p n = a n , p N +n = b n , (14)
the control vector of 2N unknown wave parameters. According to LWT, the right hand side vector, which contains observations, reads,

b m = L =1 η cos ψ m , b N +m = L =1 η sin ψ m , (15)
and the 2N × 2N matrix A mn reads,

A mn = L =1 cos ψ n cos ψ m , A m,N +n = L =1 sin ψ n cos ψ m , A N +m,n = L =1 cos ψ n sin ψ m , A N +m,N +n = L =1 sin ψ n sin ψ m . ( 16
)
The linear system of Eqs. ( 13) is solved for the model parameters (a n , b n ), with the relevant number of frequency components N being determined through a convergence study.

III -3.2 Nonlinear Wave Field

Using the weakly nonlinear ICWM formulation [START_REF] Hilmer | Deterministic wave predictions from the wamos ii[END_REF], we have,

b m = L =1 η P m , b N +m = L =1 η Q m , (17)
and

A mn = L =1 cos Ψ n + 1 2 a n k n P m , A m,N +n = L =1 sin Ψ n + 1 2 b n k n P m , A N +m,n = L =1 cos Ψ n + 1 2 a n k n Q m , A N +m,N +n = L =1 sin Ψ n + 1 2 b n k n Q m , ( 18
)
in which,

           P m = cos Ψ m -k m (a m sin Ψ m -b m cos Ψ m ) × {sin ψ Lm -[k m (a m cos ψ Lm + b m sin ψ Lm ) + 1] a m ω m k m t } + a m k m , Q m = sin Ψ m -k m (a m sin Ψ m -b m cos Ψ m ) × {-cos ψ Lm -[k m (a m cos ψ Lm + b m sin ψ Lm ) + 1] b m ω m k m t } + b m k m . (19)
Since both A mn and b m now depend on wave parameters (a n , b n ), the system of equations ( 13) is solved iteratively [START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF]. Thus, when solving for p (q+1) n at iteration q + 1, A

mn and b

(q)
m are based on wave parameters obtained at the previous iteration q. The solution is initialized at q = 0 using A (0) mn and b (0) m computed for the linear reconstruction with Eqs. ( 15) and [START_REF] Nouguier | Scattering from nonlinear gravity waves: the "choppy wave" model[END_REF]. Convergence is typically achieved within 5 to 20 iterations depending on wave steepness.

III -3.3 Regularization of the Solution

In applications, the ocean reconstruction problem may become ill-conditioned due to practical constraints, such as the heterogeneous distribution of spatial observation points, the limited ocean area observed by the optical sensor, and the frequency and direction bandwidth cutoffs in the reconstructed wave field. Nevertheless, consistent results can be achieved, independently of the conditioning of the system matrix to invert (i.e., A mn ), through a Tikhonov regularization procedure, in which the matrix inversion is replaced by a minimization problem as,

min ||A mn p n -b m || 2 -r 2 ||p n || 2 , ( 20
)
where r is the regularization parameter. The optimal value of the regularization parameter is found using the "L-curve" method, which consists in finding the r value corresponding to the point of maximal curvature (i.e., corner) of the function (log

||A mn p n -b m || , log ||p n ||).
This method thus provides an optimal compromise between minimizing the residual error and ensuring that the norm of the solution does not become too large, i.e., ill-posed. The L-curve corner can be determined analytically through solving a singular value decomposition problem [START_REF] Calvetti | L-curve and curvature bounds for tikhonov regularization[END_REF][START_REF] Hansen | Computational inverse problems in electrocardiology, Chap. The[END_REF]. Note that this procedure is equivalent to adding a constraint to the minimization problem, physically representing the total energy of the reconstructed wave spectrum, since the latter is proportional to the squared norm of p n ; in this case, -r 2 can simply be interpreted as a Lagrangian multiplier.

IV -Synthetic Application Case

In this application, similar to [START_REF] Grilli | Ocean wave reconstruction algorithms based on spatio-temporal data acquired by a flash lidar camera[END_REF][START_REF] Nouguier | Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera[END_REF], one-dimensional synthetic optical data is geometrically created for irregular sea states described by their wave energy spectrum and generated using a wave model, here a HOS model [START_REF] Bonnefoy | Advances in Numerical Simulation of Nonlinear Water Waves[END_REF][START_REF] Ducrozet | 3-d hos simulations of extreme waves in open seas[END_REF]. The ocean surface is then reconstructed and predicted based on this data set and compared to the original surface η HOS . This is repeated for realistic sea states of various characteristic steepness. Sea state parameters and optical sensor configuration are similar to those used in earlier work [START_REF] Desmars | Phaseresolved reconstruction algorithm and deterministic prediction of nonlinear ocean waves from spatio-temporal optical measurements[END_REF], as briefly recapped below.

IV -1 Setup and Error Definitions

A fixed optical sensor, located at abscissa x c and at an elevation of 30 m observes the ocean surface for x < x c , with a mean viewing angle of 76 • , an aperture of ±10 • , and along J = 64 rays, which are homogeneously distributed over the aperture angle. The sea state is represented by a JONSWAP spectrum of peak period T p = 10 s and peakedness γ = 3.3. The significant wave height H s varies between 1.5 and 6 m. Converged results were achieved using N = 50, a high cutoff wavenumber k max = 10k p , and a number of observation times K = 10 with a sampling rate of 1 Hz. The spatio-temporal target zone for wave prediction T (t) is 100 m long, centered on the sensor location x c , spans a 10 s time, from 20 to 30 s after the first waves enter the region. Hence, a point (x, t) is in the target zone if x c -50 ≤ x ≤ x c + 50 and t r + 20 ≤ t ≤ t r + 30.

Two error indicators are used to quantify the prediction accuracy. The first one is a local "point to point" error,

ε (x, t) = |η (x, t) -η HOS (x, t)| H s . (21)
In the second one, statistics are calculated by ensemble averaging 1,000 independent sea state realizations, generated for different sets of random phases. For each sea state, the average root mean square (RMS) error is calculated as,

ε RMS (t) = T (t) (η (x, t) -η HOS (x, t)) 2 dx T (t) η 2 HOS dx 1/2 , ( 22
)
with η is the predicted elevation using LWT or ICWM. This error is finally averaged over the 1,000 realizations.

IV -2 Results

The ocean surface predicted using the LWT or the ICWM models are compared with each other and with the reference solution based on the two previously defined error indicators. Thus, Fig. 2a shows ε(x, t) computed over a large spatio-temporal domain using ICWM, for a specific realization with H s = 3 m. Boundaries of the prediction zone are calculated according to Eq. ( 1) and group velocities c g1,2 based on Eq. ( 4), for k 1,2 , with µ = 5%. For this realization, the reconstruction error is very low in the spatio-temporal target zone T (t), within the prediction zone P (t). Figure 2b shows the time evolution of the same error averaged over 1,000 realizations at the sensor location x c , for 4 different characteristic steepnesses (i.e., significant wave heights), using LWT or the ICWM model. For each steepness, the error decreases from a relatively large value at t = t r , as waves enter the prediction zone (located between the vertical dashed lines in the figure) where it reaches its lowest value, and then increases to reach 0.2 after 8T p of propagation. Within the prediction zone, ICWM does reduce ε in all cases as compared to LWT, but most significantly for the largest steepnesses (i.e., nonlionearity). Also, due to nonlinear phase shifts for steep waves, the location of the minimum error shifts towards lower times as steepness increases.

Figure 3a further quantifies the prediction algorithm performance over the target zone (box in Fig. 2a), by plotting the time evolution of ε RMS integrated within it, for LWT or the ICWM model. As observed earlier, this error significantly reduces using ICWM for average to high steepness. Moreover, ε RMS only slowly varies with time and, hence, its time-average shown in Fig. 3b can provide a global estimate of the reconstruction model performance, as a function of the characteristic steepness H s /λ p . In all cases this average error reduces using ICWM rather than LWT, with the largest reduction, 13% (from 0.52 to 0.39), occurring for the steepest wave field. This confirms the importance of accurately representing nonlinear effects, with a model such as ICWM, when performing an ocean surface prediction based on observations, and in particular wave asymmetry and phase shifts, which increase as the sea state becomes more severe.

V -Experimental Data

Next, the prediction algorithm performance is assessed for an actual wave field using experimental data acquired in the wave tank of École Centrale de Nantes.

V -1 Experimental Setup

Laboratory experiments were performed in a wave tank, in which irregular waves were generated with a wavemaker, based on a JONSWAP spectrum, with T p = 10 s and γ = 3.3, and a few values of H s (or characteristic wave steepness). Time series of surface elevation were measured using 20 fixed wave gauges (Fig. 4), irregularly spaced according to a geometric law. Measurements made with these gauges are used as a proxy for data acquired at their location using an optical sensor. The equivalent optical sensor used in the model is located at an elevation of 20 m, with a mean viewing anlge of 76 • , a aperture of ±10 • , J = 20 rays, and a 1/50 th scale factor was used to account for wave tank limitations. In results reported here, only the first 16 probes are used, which are closest to the reference gauge (Fig. 4).

The number of observation times is set to K = 100, with a sampling frequency of 1.4 Hz, leading to an horizon of prediction L o = 2.78λ p , i.e., similar to that of the synthetic case (2.55). Unlike synthetic observations, only a finite and smaller number of independent realizations can be performed experimentally. However, an average can be made over dependent realizations by using a large enough time span between two reconstructions of the same wave field. Note that while convergence of result is slower using experimental data, statistical estimates can be approached with the method proposed by [START_REF] Naaijen | Limits to the extent of the spatio-temporal domain for deterministic wave prediction[END_REF]. Note, however, since the purpose of the experimental results presented here is only to illustrate the method using actual field data, such methods are not used here, but will be the object of future work.

V -2 Results

Figure 5 shows time series of surface elevation for four different wave steepness, at the location of the reference probe (red triangle in Fig. 4). Surface elevations reconstructed/predicted by the LWT and the ICWM models are shown, and compared to the laboratory measurements. As expected, the prediction is closets to measurements within the boundaries of the prediction zone. As the sea state characteristic steepness increases, the ICWM solution increasingly differs from the LWT solution, while being closer to the measured surface. Comparing the local surface shapes in the largest steepness case, one clearly sees that nonlinear wave properties are better accounted for when using the ICWM model. In particular, both the amplitude and phase of the steepest waves are better predicted.

VI -Conclusions

The characterization of the proposed prediction algorithm based on optical measurements of free surface elevation was extended through the investigation of the influence of wave steepness on the prediction performance. Using synthetic data, we focused on comparing the gain in accuracy using the ICWM model rather than the LWT model, showing significant improvements for moderate and high steepness. The method was also applied on experimental data to illustrate the performance with real field measurements. Synthetic and experimental results are in agreement, showing that the nonlinear properties modeled by ICWM, i.e., wave asymmetry and phase shifts, are of importance for short-term predictions.

Future efforts will be made on a comprehensive experimental validation of the presented algorithm, which will help further investigate the performance of ICWM for a directional wave

Figure 1 :

 1 Figure 1: Evolution of the prediction zone in time and space, for: (a) spatial data; and (b) spatiotemporal data (dash lines are prediction zones boundaries at time t k ; the increase in prediction zone relative to that of spatial observations is highlighted in red.

Figure 2 :

 2 Figure 2: Numerical error of ocean reconstruction: (a) ε(x, t) using ICWM for one realization with H s = 3 m; white box marks the spatio-temporal target zone T (t), within the prediction zone P (t) marked by oblique lines; dashed line at x = x c marks sensor location. (b) ε (averaged over 1,000 realizations) at x = x c , as a function of time, for different steepness using LWT ( ) or ICWM ( ); vertical dashed lines mark the prediction zone boundaries.

Figure 3 :

 3 Figure 3: Numerical error of ocean reconstruction ε RMS : (a) Instantaneous; and (b) timeaveraged over the target zone, for different steepnesses, using LWT ( ) or ICWM ( ).

Figure 4 :

 4 Figure 4: Locations of 16 (out of 20) irregularly spaced wave gauges (•) used in laboratory experiments; ( ) reference wave gauge.

Figure 5 :

 5 Figure 5: Time series of surface elevation in: ( ) laboratory measurements, ( ) LWT model, and ( ) ICWM model, for H s /λ p = (a) 0.006, (b) 0.018, (c) 0.029 and (d) 0.041. Vertical continuous line corresponds to the reconstruction time t r , and vertical dashed lines are boundaries of the prediction zone.

Acknowledgment

This work is supported by Naval Group and Centrale Nantes through their Joint Laboratory for Marine Technology, and by the PREDEMO-NAV project, grant ANR-15-ASTR-0006.