Jean-Christophe Sibel
email: jean-christophe.sibel@ensea.fr

Sylvain Reynal
email: reynal@ensea.fr

David Declercq
email: declercq@ensea.fr

Experimental results about the dynamics of the Generalized Belief Propagation used on LDPC codes

Keywords: iterative decoder, LDPC, region-graph, chaos

In the context of channel coding, the Generalized Belief Propagation (GBP) is an iterative algorithm used to recover the transmission bits sent through a noisy channel. To ensure a reliable transmission, we apply a map on the bits, that is called a code. This code induces artificial correlations between the bits to send, and it can be modeled by a graph whose nodes are the bits and the edges are the correlations. This graph, called Tanner graph, is used for most of the decoding algorithms like Belief Propagation or Gallager-B. The GBP is based on a non unic transformation of the Tanner graph into a so called region-graph. A clear advantage of the GBP over the other algorithms is the freedom in the construction of this graph. In this article, we explain a particular construction for specific graph topologies that involves relevant performance of the GBP. Moreover, we investigate the behavior of the GBP considered as a dynamic system in order to understand the way it evolves in terms of the time and in terms of the noise power of the channel. To this end we make use of classical measures and we introduce a new measure called the hyperspheres method that enables to know the size of the attractors.

I. INTRODUCTION

T HE channel coding is a research field whose purpose is to protect an information to transmit from environmental disturbances. The first step is the encoding of the information, a procedure in which the information, modeled as a sequence of k bits u 1 , . . . , u k , is mapped to a larger sequence of N bits x 1 , . . . , x N . The map consists in artificial correlations called constraints or parity-check equations. In [START_REF] Gallager | Low-Density Parity-Check Codes[END_REF] are introduced the Low-Density Parity-Check (LDPC) codes which are a widespread technique to encode the information. Such a code can be represented by a Tanner graph [START_REF] Kschischang | Factor Graphs and the Sum-Product Algorithm[END_REF], a graphical representation which turns out to be very useful in the second step, the decoding. In this part, the bits transmitted though a random noisy channel are iteratively handled by a decoding algorithm to create an associated output sequence of N bits that verify the whole set of parity-check equations and that must be as close as possible to the input sequence. One of the most famous decoding algorithm is the Belief Propagation (BP) introduced in [START_REF] Pearl | Probablistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF], extensively studied in [START_REF] Richardson | Error floors of LDPC codes[END_REF], [START_REF] Chung | Analysis of Sum-Product Decoding of Low-Density Parity-Check Codes Using a Gaussian Approximation[END_REF], [START_REF] Richardson | The Capacity of Low-Density Parity-Check Codes Under Message-Passing Decoding[END_REF], which is deemed to be the optimal message-passing algorithm in the case the Tanner graph of the LDPC code is loopfree. However, in most cases the Tanner graph is not loopfree [START_REF] Murphy | Loopy belef propagation for approximate inference: an empirical study[END_REF] that involves that the BP becomes suboptimal. To circumvent this problem has been proposed the Generalized Belief Propagation (GBP) [START_REF] Yedidia | Constructing free energy approximations and Generalized Belief Propagation algorithms[END_REF], [START_REF] Pakzad | Estimation and marginalization using Kikuchi approximation methods[END_REF] which is an adaptation of the Kikuchi approximation [START_REF] Bethe | Statistical Theory of Superlattices[END_REF], [START_REF] Kikuchi | A Theory of Cooperative Phenomena[END_REF] used in statistical physics. This algorithm is a message-passing algorithm which maps the Tanner graph to a non unic region-graph whose edges are media messages. The non uniqueness of the region-graph involves a large degree of freedom, thus it implies a possible way to go beyond the BP.

Along the whole paper, we focus on a particular LDPC code, the Tanner code [START_REF] Tanner | A Class of Group-Structured LDPC Codes[END_REF] whose main property in our study is the fact that it can be entirely described by a set of particular combinations of loops called the Trapping-Sets [START_REF] Richardson | Error floors of LDPC codes[END_REF], [START_REF] Sankaranarayanan | Failures of the Gallager B Decoder: Analysis and Applications[END_REF] that we use for the region-graph construction. The second section first deals with the preliminaries about the Tanner graph and the Belief Propagation. In the third section is explained the region-graph construction rules and the equations of the GBP with a novel construction of the region-graph for the Tanner code. The last section is dedicated to the exposure of some measures of dynamics to better understand the GBP and to highlight its relevant properties. Also we introduce a new measure called the hyperspheres method that enables to get the size of the attractors, that helps to evaluate the divergence of the algorithm.

II. PRELIMINARIES A. The Tanner graph

Let us consider a set of N binary random variables X = {X 1 , . . . , X N } whose global state is denoted by x = [x 1 , . . . , x N ]. In the following of the article, we use the notation x = [x 1 , . . . , x N ] to denote both the variables and their states. An LDPC code is built by a set of M constraints, or parity-check equations, C = {c 1 , . . . , c M } such that for each check c j :

c j = xi∈Nj x i
where N j is a subset of {x 1 , . . . , x N } called the neighborhood of c j depending on the LDPC code, and the sum is computed over the Galois field GF [START_REF] Kschischang | Factor Graphs and the Sum-Product Algorithm[END_REF]. The neighborhood of a variable x i is denoted by N i and it is built by the set of checks {c j } j such that x i ∈ N j . We consider that a check c j and a variable x i , such that x i ∈ N j , form an edge e ij between two nodes inside an undirected bipartite graph called the Tanner graph G = (X ∪ C, {e ij }). This graph is used as the media for the propagation of messages for the BP algorithm detailed in [START_REF] Pearl | Probablistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF], [START_REF] Kschischang | Factor Graphs and the Sum-Product Algorithm[END_REF]. An example of a Tanner graph is displayed on the figure 1. By applying the constraints C on the variables X we obtain a Markov random field whose joint probability distribution is p(X = x) written simply p(x). Each variable X i has its own marginal probability distribution p i (x i ). We denote respectively by b(x) and {b i (x i )} 1≤i≤N the estimate, by any decoding algorithm, of the joint probability distribution and the marginal probability distributions, called the beliefs. The ouput sequence x that is the estimate of the input sequence in the transmission channel is the most likely sequence given by:

x 0 x 1 x 2 x 3 x 4 x 5 x 6 c 0 c 1 c 2
x = arg max x b(x)
The noisy transmission channel that we use in all of our simulations is an additive white Gaussian noise channel whose power is σ 2 . The observations Y = {y 1 , . . . , y N } we get at the output of the channel enable the computation of the likelihoods

L(x) = {l 1 (x 1 ) = p(y 1 |x 1 ), . . . , l N (x N ) = p(y N |x N )}.

B. Belief Propagation

The BP is an iterative algorithm that consists in passing messages between the variables and the constraints on the edges of a given Tanner graph. An iteration k ≥ 1 of this algorithm is classically done in two half iterations:

• for each edge e ij , we compute the messages n (k) ji (x i ) from the the constraint C j to the variable x i :

n (k) ji (0) = 1 2 + 1 2 X i ∈Nj \Xi 2m (k-1) i j (0) -1 (1) n (k) ji (1) = 1 -n (k) ji (0)
• for each edge e ij , we compute the messages m

(k) ij (x i ) from the variable x i to the constraint C j : m (k) ij (x i ) = l i (x i ) Z ij C j ∈Ni\Cj n (k) j i (x i ) (2)
where Z ij is a normalization factor such that:

m (k) ij (0) + m (k) ij (1) = 1,
To initialize the algorithm, we use the likelihoods:

∀e ij , m (0) ij = l i (x i ) (3)
The beliefs on a variable x i are computed as the geometric averages on the incoming messages:

b (k) i (x i ) = l i (x i ) Z i j∈Ni n (k) ji (x i ) (4)
where Z i is a normalization factor such that: b i (0)+b i (1) = 1.

III. THE GENERALIZED BELIEF PROPAGATION A. The region-graph construction

The region-graph R is a directed graph of depth D built level by level. We decompose this construction in two steps.

The first step is the construction of the first level R 0 . The principle is to gather the nodes of the Tanner graph in order to absorb some harmful topological structures, as loops or combinations of loops. These gatherings are the nodes of R 0 if and only if each check c j is included in at least one gathering accompanied by its neighborhood N j . We call a gathering a region.

The second step is the construction of the next levels R 1 , . . . , R D-1 . To build a level R l we search for the intersections between the regions of the previous level R l-1 . To this end, we first define C r the set of constraints and X r the set of variables that both form the region r. A set of variables and checks r = X r ∪ C r is a region of R l if and only if there is a set of n ≥ 2 regions {r 1 , . . . , r n } ∈ R n l-1 such that:

• ∀c j ∈ C r , ∀r k ∈ {r 1 , . . . , r n }, c j ∈ C k , • ∀x i ∈ X r , ∀r k ∈ {r 1 , . . . , r n }, x i ∈ X k .
A region-graph of the Hamming code presented previously is displayed on the figure 2. The nodes of the first level are constructed considering only one check per region.

c 0 x 0 , x 1 , x 2 , x 4 c 1 x 0 , x 1 , x 3 , x 5 c 2 x 0 , x 2 , x 3 , x 6 x 0 , x 1
x 0 , x 2 x 0 , x 3

x 0

Fig. 2. A region-graph of the Hamming code

The efficiency of the GBP detailed below is fully dependent on the region-graph that is built from the Tanner graph. The construction of the first regions R 0 determines the performance of the algorithm. Therefore, it is relevant to adapt the construction to the topology of the considered Tanner graph. The difficulty is to find a construction that both implies relevant performance and fair complexity. In the case of the previous example where we consider only one check per region of R 0 , the construction is easy to implement and the complexity is low enough to compare it with other algorithms though the performance are not improved enough or even worse.

B. The message-passing algorithm

The region-graph is a Bayesian network whose probability distribution can be computed by the following algorithm. The principle of the GBP is to convey messages on the edges of the region-graph where a message m (k) rq from a region r to a region q is the a posteriori probability distribution of the region q given by the region r at iteration k. An iteration of the algorithm consists in the computation, or update, of all the messages {m (k) rq } r,q and the computation of all the beliefs {b (k) r } r and {b (k) i } i whose equations are detailed just below. We denote by K the maximum value of k.

1) The update equations: The family F r of a region r is defined as the set of regions such that for any region q ⊆ F r :

c j ∈ C q ⇒ c j ∈ C r x i ∈ X q ⇒ x i ∈ X r
We define the descendants D r of a region r as F r \r. The children C r of r are the regions q of D r such that there is an edge from r to q. According to [START_REF] Yedidia | Constructing free energy approximations and Generalized Belief Propagation algorithms[END_REF], [START_REF] Pakzad | Estimation and marginalization using Kikuchi approximation methods[END_REF] the equations of the messages from a region r to a region q at iteration k are:

m (k) rq (x q ) = x r\q L r\q (x r\q )c r\q (x r ) p⊂R\Fr s⊂Fr \Fq m (k-1) ps (x s ) Z rq p⊂Dr \Fq s⊂Dq m (k) ps (x s )
(5) where:

• x r\q is the state of X r\q = X r \X q , • L r\q (x r\q ) = xi⊂X r\q l i (x i ), • c r\q (x r ) = cj ∈Cr\Cq c j (x Nj ),
• Z rq is a normalization factor to ensure that:

xq m (k) rq (x q ) = 1.
For any region r at iteration k the beliefs are:

b (k) r (x r ) = 1 Z r L r (x r )c r (x r ) p⊂R\Fr s⊂Fr m (k) ps (x s ) (6)
where Z r is a normalization factor to ensure that

xr b (k) r (x r ) = 1.
To get the beliefs for the single variables we only need to compute the marginal probability distributions of the regions beliefs. Thus for any region r:

∀x i ∈ X r , b (k) i (x i ) = 1 Z i x r\i b (k) r (x r ) (7)
where Z i is a normalization factor to ensure that

xi b (k) i (x i ) = 1.
2) The relaxation coefficient: The update equations ( 5) can be summarized into an implicit equation:

∀(r, q) ∈ R 2 s.t. q ⊂ C r , m (k) rq = F rq {m (k-1) ps } p,s , {l i } i (8)
In the following, we lighten the notations by the use of F rq alone. An important point for an iterative algorithm is the convergence, that we define here as:

∀(r, q) ∈ R 2 s.t. q ⊂ C r , m (k+1) rq = m (k) rq
In [START_REF] Yedidia | Constructing free energy approximations and Generalized Belief Propagation algorithms[END_REF] is introduced the fact that the GBP suffers from a poor convergence. To circumvent this phenomenon is included a decreasing term of relaxation w k such that:

m (k) rq = w k F rq + (1 -w k )m (k-1) rq (9)
The main issue is to find the most suitable function w k which provides both convergence and relevant performance.

We propose three functions for w k that are displayed on the figure 3. The balance between the memory m (k-1) rq and the update F rq is different from case to case. In the case of a parabolic function, the memory is quickly highlighted at the expense of the update, whereas in the case of an inverse parabolic function the memory is practically ignored for a long time. The affine coefficient enables a progressive oversight of the update up to a complete use of the memory.

On the figure 4 are displayed the bit-error-rates (BER) on the Hamming code for the different relaxation coefficients. It appears that the BER are not the same, which indicates that the choice of w k is not trivial. We observe that the inverse parabolic coefficient involves the best performance whatever w 0 , unlike the affine coefficient or the parabolic coefficient which is the worst one. Therefore, we can assert that the best choice for w k is the inverse parabolic one. Another idea would be to create a function w k whose values are computed at each iteration k according to the difference between the update and the memory, that can be considered as an adaptative relaxation coefficient. Yet, the results are really poorer than what we can get from the function proposed previously so we dropped this function. Finally, extensive simulations let us conclude that the choice of the relaxation coefficient is completely dependent on the code we use in the transmission. Therefore, we need to test different functions w k for each code to select the best one.

C. A novel construction

We present here a particular construction of the region-graph for a given LDPC code. We consider the Tanner code [START_REF] Tanner | A Class of Group-Structured LDPC Codes[END_REF] of length N = 155 with M = 93 parity-check equations. The Tanner graph of this code can be established by a sophisticated concatenation of topological structures called Trapping-Sets (TS). A TS(a,b) introduced in [START_REF] Richardson | Error floors of LDPC codes[END_REF] and studied in [START_REF] Sankaranarayanan | Failures of the Gallager B Decoder: Analysis and Applications[END_REF] among other, is a Tanner graph which contains a variables and b unsatisfied parity-check equations. In the studied case the TS(5, 3), whose Tanner graph is represented on the figure 5, are sufficient to cover the whole Tanner graph of the code.

: variable : constraint : unsatisfied constraint Such a structure is known to damage the decoding [START_REF] Sankaranarayanan | Failures of the Gallager B Decoder: Analysis and Applications[END_REF] because only the total null state enables to verify all the parity-check equations. Thus to build the first level R 0 of the region-graph, it would be relevant to absorb them into regions. However, the complexity of the message-passing would soar because these regions are too large. Therefore, we had better break these structures into several smaller regions as on figure 6 with:

x 0 x 1 x 2 x 3 x 4 c 0 c 1 c 2 c 3
• r = {X r = {x 0 , x 1 , x 3 }, C r = {c 0 , c 2 , c 6 }}, • p = {X p = {x 1 , x 2 , x 3 }, C p = {c 1 , c 4 , c 7 }}, • q = {X q = {x 1 , x 3 , x 4 }, C q = {c 3 , c 5 , c 8 }}. X r , C r X p , C p X q , C q X r∩p∩q , C r∩p∩q

D. Results

The region-graph of the TS(5, 3) is loopfree which makes the decoding optimal. When we apply this construction to all the TS(5, 3) of the whole code, unfortunately we do not get a loopfree region-graph because all the variables belong to several TS [START_REF] Chung | Analysis of Sum-Product Decoding of Low-Density Parity-Check Codes Using a Gaussian Approximation[END_REF][START_REF] Pearl | Probablistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]. Nevertheless, the performance are still relevant to make the GBP a good candidate for the decoding. We present on the figures 7 the BER in terms of the iteration of the GBP on the Tanner code for error events made by the noisy channel that are truly harmful for the Belief Propagation particularly because of the TS [START_REF] Chung | Analysis of Sum-Product Decoding of Low-Density Parity-Check Codes Using a Gaussian Approximation[END_REF][START_REF] Pearl | Probablistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]. On the figure (b) is displayed the BER of the BP which is oscillating, making the decoding divergent and wrong. On the same figure is displayed the BER of the GBP when the memory is not taken into account. It clearly appears that such an algorithm does not bring any improvement noticing that it is often worse than the BP. When w k is set to 0.5 then the BER is dramatically reduced. On the figure (a) is presented the fact that to keep a smooth evolution of the BER with relevant values the relaxation coefficient must not imply a total oversight at the first iteration. Extensive simulations make us convinced that for the Tanner code, the initial value w 0 = 0.85 is the most relevant to get the best results, whatever the relaxation coefficient decreasing rule. The most interesting point is the comparison between figures (c) and (d) where n is the order of the parabola and of the inverse parabola:

• parabolic:

w k = w 0 K-k K n • inverse parabolic: w k = w 0 1 -k K n
The relaxation coefficients have dual behaviors as it was explained previously. This reveals that considering the update rule for too many iterations damage the performance of the GBP. Therefore it is more efficient to use this update for some iterations at the beginning and quickly to foster the memory. The parabolic relaxation coefficient offers the most accurate estimate of the input sequence with the most stable evolution.

In the following of the article, we deeply study this last assumptions to bring out the dynamics of the GBP with a parabolic relaxation coefficient.

IV. DYNAMICS

In this section, we expose some measures to evaluate the dynamical behaviors of the GBP compared with the BP. To this end, we first define the mathematical environment in which these measures are relevant.

A. State space definition and properties

To compare the dynamics of the BP and GBP, we need to define two state spaces that are similar. Using the fact that both are message-passing algorithms, we should consider the messages as the state variables. This choice fits perfectly with the convergence conditions written in the previous section for the GBP and written down here for the BP:

∀e ij , n (k) ji = n (k-1) ji (10)
The update equations of the BP can be implicitly noted:

∀e ij , m (k) ij = f ij {n (k-1) yx } (x,y) , l i (11)
n (k) ji = g ji {m (k) xy } (x,y) (12)
The messages {n

(k)
ji } are then computed by a composed function g • f that maps the messages {n (k-1) ji } and the likelihoods. Thus the BP can be described by a unic equation:

∀e ij , n (k) ji = G ij ({n (k-1) nm } (m,n) , {l i } i ) (13)
We consider respectively {G ij } (i,j) and {F rq } (r,q) as two sets of iterated maps on the state variables {n

(k)
ji } (i,j) and {m (k) rq } (r,q) which are called trajectories in the associated state spaces, denoted E BP and E GBP , whose dimensions are the number of messages to compute at each iteration. We denote by U (k) = {n (k) ji } (i,j) and by V (k) = {m (k) rq } (r,q) the points of the trajectories of the BP and the GBP at iteration k.

B. Parameters and scaling

In [START_REF] Zheng | Study of bifurcation behavior of LDPC decoders[END_REF] the value of σ 2 , or the corresponding Signal to Noise Ratio (SNR), is used as a parameter such that different values imply different motions of the BP. However, most of their simulations are done for particular noise realizations and scaled on the SNR, that prevents from evaluating a statistical behavior. A reason is that the noise realizations that lead the BP not to converge or to converge to a wrong estimate are rare events, essentially because the LDPC codes and the iterative algorithms are created to this end. A way to have some statistical evaluations of the behavior of the BP and the GBP is the following:

1) find some of these noise realizations, 2) store the corresponding initializations on the state variables, 3) average the quantities to measure for a sufficient set of initializations that are close in the state space in the sense of the Euclidean distance. By this way, we can target the guenine critical values. For all the measures presented in the following of the paper, we use this method to get statistical results which are relevant enough to describe the behavior of the BP and the GBP. We call critical values the SNR that correspond to a radical change in the behavior of the motion of the algorithms, which are the bifurcations.

C. Bifurcation diagram

A relevant method to extract the critical values of the SNR is the bifurcation diagram computation. For a given noise realization, it consists in evaluating the value of a function E that is computed from the state variables at their steady state for J values of the SNR. We get a sequence [E σ1 , . . . , E σ J ] that represents the behavior of the dynamic system in terms of the SNR.

We consider the following function exposed in [START_REF] Zheng | Study of bifurcation behavior of LDPC decoders[END_REF] called the mean square beliefs:

∀σ, E σ = 1 N N -1 i=0 b 2 i (0) ( 14
)
where the input sequence in the channel is the null sequence and the beliefs are computed at the last iteration K of the considered algorithm, BP or GBP. Obviously, there is no reason that the associated dynamic system has reached any steady state at this iteration but we need to suppose it for computation time's sake. This function presents three important values:

• E σ = 1: all the beliefs indicate that the ouput sequence is the null sequence which is a good decoding,

• E σ = 1
4 : all the beliefs are uniform distributions thus there is no information about the true state of the transmitted bits, which is a missed decoding,

• E σ = 0: all the beliefs indicate that the output sequence is the complementary of the sent sequence, which is a completely wrong decoding. The display of [E σ1 , . . . , E σ J ] enables to know two properties of the used algorithm: the amplitudes provides information about the decoding performance, and the variation between successive values gives us the critical values of the SNR. We display on the figure 8 the mean bifurcation diagrams of the BP and the GBP for a given noise realization computed as we exposed at the beginning of the current subsection.

D. Reduced trajectory

Another use of the mean square beliefs function is the representation of the trajectory in a 3-dimensional state space. To this end, we use the phase space reconstruction introduced in [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF]. The method is first to compute this function at each iteration to get the following sequence

E σ = [E σ (k)] 0≤k≤K .
Afterthat we share this one dimensional sequence in a three dimensional sequence as follows:

Ẽσ =    E σ (0) E σ (1) E σ (2) . . . . . . . . . E σ (K -2) E σ (K -1) E σ (K)    (15)
On the figures 9 and 10 are displayed some reduced trajectories of the BP for typical values of the SNR deduced from the previous bifurcation diagram. It Fig. 10. Reduced trajectory for the BP on the Tanner code with SNR = 2.70 dB and 3.00 dB typical behaviors that match with the four intervals exposed in the previous paragraph. We obtain a very small sized attractor for SNR = 2.10 dB that can be considered as a fixed point, whereas the reduced trajectory transforms to a limit cycle when the SNR is between 2.19 dB and 2.49 dB. A crucial point is that the thickness of the trajectory along this limit cycle increases as the SNR is getting greater up to 2.50 dB. At the same time this limit cycle interleaves with other limit cycles, that can be understood as a sequence of period doubling bifurcations in terms of dynamic system, as is displayed on the figure 11 with two interleft cycles. Such a phenomenon is a typical route to chaos, that is observable from 2.51 dB. A chaotic motion means that there is not any periodic motion or fixed point convergence anymore, as it is displayed for 2.70 dB. When the SNR reaches 2.99 dB the trajectory collapses to a single point that is a true fixed point. Concerning the GBP, whose reduced trajectories are displayed on the figures 12 and 13, we cannot share the SNR values so accurately because the reduced trajectory does not transform blatantly. However, the reduced trajectory makes reveal also four different behaviors that follow the same order than that of the BP: small attractor, limit cycle, chaos, fixed point. The corresponding SNR intervals match with what have been revealed previously. Reduced trajectory for the GBP on the Tanner code with SNR = 2.80 dB and 3.00 dB We have to be cautious because Ẽσ is not a true trajectory, it does not respect the Cauchy-Lipschitz condition [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF] due to the non bijection between the messages and the beliefs. Thus, this sequence only has the role of giving clues about the true behavior of the considered algorithm as the possible shape of the actual trajectory in E BP or E GBP that are: convergence to a fixed point, convergence to a limit cycle, convergence to a chaotic attractor. To distinguish these shapes, we need a criterion that reflects the behavior by its own value. A good candidate is the Lyapunov exponent.

E. Lyapunov exponents

A common measure is the Lyapunov exponent λ [START_REF] Rosenstein | A practical method for calculating largest Lyapunov exponents from small data sets[END_REF], [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF], [START_REF] Wolf | Determining Lyapunov Exponents from a Time Series[END_REF]. Its computation consists in evaluating at each iteration k ≤ K the Euclidean distance d k between two initially close trajectories, and computing the log-ratio:

λ = ln d K d 0 (16)
The sign of the Lyapunov exponent λ reveals the behavior of the system around the corresponding initialization of the trajectories.

• if λ is positive then the trajectories have moved away one from the other, which is an evidence of a chaotic behavior, • if λ is negative then the trajectories have got closer, which is an evidence of a convergent behavior to a small sized volume of the state space. This volume is reduced to a fixed point if and only if λ → ∞. When λ crosses the x-axis the system suffers from a bifurcation meaning that the algorithm has changed of motion. The corresponding SNR are the critical values. We display on the figure 14 the Lyapunov exponents of the BP and the GBP on the Tanner code. As we have observed on the bifurcation diagrams, the evolutions are really different as soon as the SNR is greater than 2.07 dB. The BP curves is perfectly consistent with the associated bifurcation diagram in the sense that the critical values we extract are the same and the behaviors we could imagine by the bifurcation diagram are also revealed by the Lyapunov exponent. Concerning the GBP, we observe that the Lyapunov exponent does not follow blatantly different motions as we noticed about the reduced trajectory. Moreover, it reaches greater values than the BP, which is not in its favour.

However, a crucial point is that the SNR that corresponds to the fall of the Lyapunov exponent is clearly lower for the GBP than for the BP. In association with the results about the error correction power in the previous section, it appears that from this SNR of 2.88 dB the GBP presents a stable evolution and better skills to correct the errors on the transmitted bits.

A relevant analyze we need to effect is the comparison with the reduced trajectory we exposed previously so as to associate accurately with each reduced motion a particular evolution of λ. Here are these associations for the BP:

• SNR ∈ [0dB; 2.19db]: the reduced trajectory converges to a very small sized volume of E BP that we can consider as a fixed point whereas λ is close to the null value, • SNR ∈ [2.20db; 2.49db]: the reduced trajectory is trapped into a limit cycle whereas λ has gone over a stair, • SNR ∈ [2.50db; 2.98db]: the reduced trajectory does not converge to any fixed point, limit cycle or quasi-limit cycle but to a chaotic attractor whereas λ soars to high values, • SNR ∈ [2.99db; +∞db]: the reduced trajectory converges to a fixed point corresponding to a good decoding. As we mentionned previously, the GBP does not present behaviors that are easy to distinguish therefore it would not be relevant to associate any possible critical value between the bifurcation diagram, the reduced trajectory and the Lyapunov exponent.

F. Hyperspheres method

We propose here a novel method to evaluate the unstability of the BP and the GBP, based on their own trajectory in E BP and E GBP . This method is complementary to the Lyapunov exponent measure because it reveals the size of the attractor that the trajectory falls into and some other properties about the limit cycles.

This method consists in computing the rays R k of the hyperspheres circumscribe to the trajectory inside a given temporal window centered around each point U (k) (or V (k) ) of the trajectory. On the figures 15 and 16 are displayed the evolution of two rays that correspond to two initially close trajectories in the Euclidean sense. The motion we observe for SNR = 2.30 dB is consistent with the limit cycle we observed on the reduced trajectory. The curve of the ray enable to estimate the period of the trajectory around 23 iterations.

Moreover we can assert that this limit cycle is stable because the two rays cannot be distinguished. For 2.70 dB the rays moved away one from the other as it was predicted by the Lyapunov exponent observations. More accurately we can see that the rays has different oscillations step. This is due to the period doubling bifurcations explained previously. During 92 iterations in average for k ≤ 500, the trajectory is trapped in a given limit cycle and for the next 92 iterations the trajectory falls into another limit cycle of different ray. For both it is possible to measure the period or pseudo-period that is the same as the period of the first limit cycle, that is 23 iterations.

For k ≥ 500 we cannot really distinguish these different phases of evolution because the period doubling has led to the chaos. Such an observation makes our method relevant to bring out crucial information by a one dimensional function.

Another important aspect of the hyperspheres method is the raising of the behaviors difference between two initially close trajectories: we easily observe that the evolution of the rays cannot be distinguished while k ≤ 200 but as k is getting greater, the evolution of the rays move away one from the other but they follow the same kind of motion. For both of them the hypervolumes of the state space in which they are locking in are quite of the same size. When the SNR reaches the last critical value we observe that one of the rays decreases to the null value because the BP has converged to a fixed point.

The other ray has not collapsed yet because the SNR is just at the critical value, if it was increased a little we would see the two rays going to zero. Concerning the GBP, for low SNR values we did not display the evolutions of R k because it follows the same shape as the BP one. On the figure 17, at SNR = 2.70 dB we observe that the rays of the couple of trajectories are quite smaller than that of the BP. It appears that the maximum value of the rays is still lower than the minimum value of the BP ray.

The comparison between both rays on the GBP enables to deduce that even though each trajectory eventually reaches a particular attractor of the state space, their size are not very different and they collapse at the end, which is not true about the BP, making the GBP more stable than the BP for the suited construction we have proposed.

V. CONCLUSION

In this paper, we have presented a novel method of construction of the region-graph for the GBP algorithm. The results have shown firstly that the error correction power was improved compared with the BP one, and secondly that the GBP was more stable than the BP. The measure of the hyperspheres we have introduced enable to bring out this property and also the fact that the GBP converges towards small sized attractors contrary to the BP. Finally, this work leads to conclude that a suited region-graph construction enables the GBP to be a good candidate to surpass the famous BP. In future works, we propose to find a mathematical criterion to evaluate the relevancy of any region-graph according to the associated GBP performance, and also we propose other handlings of the TS(5, 3) to build region-graphs that could present better correction power and improved stability.

Fig. 1 .

 1 Fig. 1. Tanner graph of the Hamming code

Fig. 3 .

 3 Fig. 3. Relaxation coefficient : (A) affine, (P) parabolic, (IP) inverse parabolic

  parabolic: w 0 = 0.9 : w 0 = 0.7 . : w 0 = 0.5+ : w 0 = 0.3 -: w 0 = 0.1

Fig. 4 .

 4 Fig. 4. Binary Error Rate on the Hamming code for different relaxation coefficents

8 Fig. 5 .

 85 Fig. 5. TS(5, 3) : 5 variables and 3 unsatisfied constraints

Fig. 6 .

 6 Fig.6. Region-graph of the TS[START_REF] Chung | Analysis of Sum-Product Decoding of Low-Density Parity-Check Codes Using a Gaussian Approximation[END_REF][START_REF] Pearl | Probablistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]

Fig. 7 .

 7 Fig. 7. BER of the BP and the GBP for different relaxation coefficients on the Tanner code with the suited region-graph construction

Fig. 8 .

 8 Fig. 8. Bifurcation diagrams of the BP and the GBP

Fig. 11 .

 11 Fig. 11. Reduced trajectory for the BP on the Tanner code with SNR = 2.40 dB

Fig. 12 .

 12 Fig. 12.Reduced trajectory for the GBP on the Tanner code with SNR = 2.00 dB and 2.20 dB

Fig. 15 .

 15 Fig. 15. Evolution of two hyperspheres rays corresponding to two initially close trajectories of the BP on the Tanner code at SNR = 2.10 dB, SNR = 2.30 dB

Fig. 16 .

 16 Fig. 16. Evolution of two hyperspheres rays corresponding to two initially close trajectories of the BP on the Tanner code at SNR = 2.70 dB, SNR = 3.00 dB

Fig. 17 .

 17 Fig. 17. Evolution of two hyperspheres rays corresponding to two initially close trajectories of the GBP on the Tanner code at SNR = 2.70 dB, SNR = 3.00 dB

  while the SNR is lower than 2.49 dB, for values in [2.5 dB, 2.98 dB] however it does not appear any known evolution which is an indication of chaos. Concerning the GBP, it appears globally three intervals: in [2.08 dB, 2.43 dB] the shape of E σ indicates some irregular oscillations whereas for SNRs in [2.44 dB,2.69 dB] we cannot assert anything except that the chaos would appear. From 2.70 dB up to 3.02 dB, the GBP tends to the right decoding state whereas the BP does not present this behavior before 2.99 dB. Moreover the shapes of the two whole signals indicate that the GBP function E σ is globally beyond the BP one, that shows that the GBP tends faster than the BP to the right state.

We observe that for SNRs lower than 2.19 dB, the BP follows a regular increasing steady motion. Such a behavior is analoguous to that of the GBP except that the critical value is 2.07 dB. When the SNR is greater than these two critical values, the algorithms follow two distinct motions. The BP seems to oscillate

  appears for the BP four

	0.94	0.94
	0.88	0.88
	0.82	0.82
	0.94	0.94
	0.88	0.88
	0.82	0.82
	0.82	0.82
	0.88	0.88
	0.94		0.94
	Fig. 9. Reduced trajectory for the BP on the Tanner code with SNR = 2.15
	dB and 2.30 dB
	0.94	1
		0.94
	0.88
		0.88
	0.82	0.82
	0.94	1
		0.94
	0.88
		0.88
	0.82	0.82
	0.82	0.82
	0.88	0.88	0.94
	0.94		1

ACKNOWLEDGMENT

Supported by the French ANR Defis program under contract ANR-08-EMER-003 (COCQ project).


J C Sibel

S Reynal

On the Belief Propagation and its Dynamics

 as an inference method to estimate the joint probability distribution of a given Bayesian network. Many works have raised up the use of BP in numerous applications as in image processing, computer vision, neural network, statistical physics and channel coding.

Among examples in the channel coding realm, we focus on the Low-Density

Parity-Check (LDPC) codes used to protect a sequence of bits sent through a noisy channel. Any LDPC code can be represented by a graph whose edges (resp. the vertices) stand for the parity-check equations (resp. the bits) which protect the information against the noise of the channel. This graph is equivalent to a Markov random field. It is used to run the algorithm by iteratively passing messages along the edges in such a way that a message is the probability distribution of the destination conditionned to the emitter. LDPC codes are sufficiently robust against the channel noise for the algorithm to converge after some iterations to the most likely state of the estimate of the joint probability distribution. LDPC codes provide a lower bit-error rate than a no coded transmission for any values of the Signal-to-Noise Ratio (SNR).

However, most LDPC codes exhibit loop-like topological structures that ruin the BP performance in terms of the number of corrected errors. The short loops and the combinations thereof are especially harmful. These defects drive the algorithm toward states that are not the most likely ones, or even toward complex attractors that need to be defined.

A few studies have been carried out that propose an analytic approach to the BP convergence [START_REF] Chung | Analysis of sum-product decoding of low-density parity-check codes using a gaussian approximation[END_REF][START_REF] Ruffer | On robust stability of the belief propagation algorithm for LDPC decoding[END_REF]. Other authors have focused on experimental results to extract properties about the behavior of the BP algorithm [START_REF] Zheng | Study of bifurcation behavior of LDPC decoders[END_REF]. In our work, we make use of well-known estimators and new ones to describe the dynamics of the algorithm in terms of the SNR.

One is the bifurcation diagram which we use to extract the critical values of the SNR for which the algorithm blatantly changes its behavior. Experiments we conducted on many noise realizations that lead to complex attractors provide three critical values that divide the SNR values into four intervals. Each interval Email address: jean-christophe.sibel@ensea.fr (J.C. Sibel) Preprint submitted to Satellite Meeting Complex Behavior in Discrete Dynamical SystemsJuly 9, 2012 corresponds to particular dynamics.

To visualize these dynamics, we create a three dimensional reduced state space by the phase space reconstruction method. It brings out four different behaviors announced by the bifurcation diagram. Whatever the harmful noise realization we get back these four patterns as the SNR is increased: small sized attractor, limit cycle, chaos and fixed point.

To assert that the algorithm suffers from chaos and to evaluate its intensity we compute the largest Lyapunov exponent λ. The simulations provide the confirmation of these four behaviors and assert that the third interval matches with chaos of the BP.

We use the Jacobian matrix to underline the state parameters responsible for chaos. Though, the loop-like structure propagates the information along most edges, whatever its nature. Thereby we obtain from this estimator that all the messages suffer from a chaotic dynamic.

To complete the investigation, we propose a new estimator. It consists in computing the radius of the circumscribed hypersphere of the true trajectory in a given temporal window. By sliding this window, we get a sequence of radii. This sequence has the advantage to be one dimensional, whereas the trajectory is multi dimensional. Furthermore, it describes accurately the shape of the attractor and its diameter. This helps to understand the nature of the considered attractor. An extra use of these method is to compute the iterative evolution of the mutual distance between the hyperspheres of two initially close trajectories. Thus this estimator highlights the stability brought out by the largest Lyapunov exponent and it also provides the attractor nature.

To conclude, this study helps to classify four distinguished dynamical behaviors of the BP along the SNR values. It confirms that the BP is widely infuenced by the topology of the considered LDPC code. The estimators, either well-known or new, provide qualitative and quantitative descriptions of the algorithm, especially for the complex dynamics.

Keywords: Inference, message-passing, loops, stability, chaos.


Yohann Le Floch
email: yohann.lefloch@univ-rennes1.fr

Singular Bohr-Sommerfeld conditions for 1D Toeplitz operators: elliptic case

In this article, we state the Bohr-Sommerfeld conditions around a global minimum of the principal symbol of a self-adjoint semiclassical Toeplitz operator on a compact connected Kähler surface, using an argument of normal form which is obtained thanks to Fourier integral operators. These conditions give an asymptotic expansion of the eigenvalues of the operator in a neighbourhood of fixed size of the singularity. We also recover the usual Bohr-Sommerfeld conditions away from the critical point. We end by investigating an example on the two-dimensional torus.

Introduction

Let M be a compact, connected Kähler manifold of complex dimension 1, with fundamental 2-form ω. Assume M is endowed with a prequantum line bundle L. Let K be another holomorphic line bundle and define the quantum Hilbert space H k as the space of holomorphic sections of L k ⊗ K, for every positive integer k. The operators acting on H k that we consider are Berezin-Toeplitz operators ( [START_REF] Boutet De Monvel | The spectral theory of Toeplitz operators[END_REF][START_REF] Borthwick | Semiclassical spectral estimates for Toeplitz operators[END_REF][START_REF] Charles | Berezin-Toeplitz operators, a semi-classical approach[END_REF][START_REF] Ma | Toeplitz operators on symplectic manifolds[END_REF], and many others). The semiclassical parameter is k, and the semiclassical limit is k → +∞. Formally, k is the inverse of Planck's constant .

Our aim is to understand the spectrum of a given self-adjoint Toeplitz operator, in the semiclassical limit. In the setting of ( -)pseudodifferential operators, the similar study was done by Colin de Verdière in [START_REF] De Verdière | Spectre conjoint d'opérateurs pseudo-différentiels qui commutent[END_REF]. In his article [START_REF] Charles | Symbolic calculus for Toeplitz operators with half-form[END_REF], Charles obtained the description of the intersection of the spectrum of a self-adjoint Toeplitz operator with an interval of regular values of its principal symbol, in the semiclassical limit: the eigenvalues are selected by an integrality condition for some geometric quantities (actions) associated to the symbol of the operator (these are the Bohr-Sommerfeld conditions).

In this article, we extend these conditions around a global minimum of the principal symbol; since we work with only one degree of liberty, we expect to have a very precise description of the eigenvalues near the critical point.

Main theorem

Let A k be a self-adjoint Toeplitz operator on M ; its normalized symbol a 0 + a 1 + . . . is real-valued. Assume that its principal symbol a 0 admits a global minimum at m 0 ∈ M , with a 0 (m 0 ) = 0. Denote by λ [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform[END_REF] k ≤ λ [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform. Part II. A family of related function spaces. Application to distribution theory[END_REF] k ≤ . . . ≤ λ (j) k ≤ . . . the eigenvalues of A k . Our main result is the following theorem.

Theorem (Theorem 6.2). There exist E 0 > 0, a sequence g(., k) of functions of C ∞ (R, R) which admits an asymptotic expansion of the form g(., k) = ≥0 k -g in the C ∞ topology, and a positive integer k 0 such that for every integer N ≥ 1 and for every E ≤ E 0 , there exists a constant C N > 0 such that for k ≥ k 0 :

λ (j) k ≤ E or E (j) k ≤ E ⇒ λ (j) k -E (j) k ≤ C N k -N where E (j) k = g k -1 j + 1 2 , k , j ∈ N.
This allows to compute asymptotic expansions to all order for eigenvalues of A k lower than E 0 , except that so far, we do not know who are the g , ≥ 0, or how to compute them. In fact, g(., k) is constructed as the local inverse of a sequence f (., k) which also admits an asymptotic expansion f (., k) = ≥0 k -f in the C ∞ topology, and the first terms f 0 , f 1 are related to geometric quantities (actions) associated to A k .

Link with the usual Bohr-Sommerfeld conditions

Let I be a set of regular values of the principal symbol a 0 ; for every E in I, the level set Γ E := a -1 0 (E) is diffeomorphic to S 1 . Fix an orientation on Γ E depending continuously on E. Define the principal action c 0 ∈ C ∞ (I) in such a way that the parallel transport in L along Γ E is the multiplication by exp(ic 0 (E)). Of course, c 0 (E) is defined up to an integer multiple of 2π, but we can always choose a determination of c 0 that is smooth on I.

Let (δ, ϕ) be a half-form bundle, that is a line bundle δ → M together with an isomorphism of line bundles ϕ : δ 2 → Λ 2,0 T * M . It is known that for any connected compact Kähler manifold of complex dimension 1, such a couple exists. Introduce the hermitian holomorphic line bundle L 1 such that K = L 1 ⊗ δ. For E in I, define the subprincipal form κ E as the 1-form on Γ E such that κ E (X a 0 ) = -a 1

where X a 0 stands for the Hamiltonian vector field associated to a 0 . Denote by j E the embedding Γ E → M , and introduce the connection ∇ E on j * E L 1 → Γ E defined by

∇ E = ∇ j * E L 1 + 1 i κ E ,
with ∇ j * E L 1 the connection induced by the Chern connection of L 1 on j * E L 1 . Define the subprincipal action c 1 ∈ C ∞ (I) like the principal action, replacing L by L 1 endowed with this connection.

Finally, define an index from the half-form bundle δ as follows: the map ϕ E : δ 2 E → T * Γ E ⊗ C, u → j * E ϕ(u) is an isomorphism of line bundles. The set u ∈ δ E ; ϕ E (u ⊗2 ) > 0 has one or two connected components. In the first case, we set E = 1, and in the second case E = 0. In fact, E is a constant E = for E in I.

If we define carefully c 0 and c 1 , the following result holds.

Proposition (Proposition 6.4). Set I =]0, E 0 [. Then

f 0 = 1 2π c 0 , f 1 = 1 2π c 1 (1)
on I.

Thus, we recover the regular Bohr-Sommerfeld conditions away from the minimum.

Structure of the article

The paper is organized as follows: we start by recalling some properties of Toeplitz operators on a compact manifold. Then, we briefly explain how to adapt the theory in the case where the phase space is the whole complex plane. The fourth section is devoted to the construction of Fourier integral operators that we use to construct our microlocal normal form in the following part. In section 6, we state the Bohr-Sommerfeld conditions and some consequences. In the last section, we investigate an example to give some numerical evidence of our results.

Preliminaries and notations

First, we introduce the notations and conventions that we will adopt through this whole article. They are already written in [START_REF] Charles | Symbolic calculus for Toeplitz operators with half-form[END_REF] for instance, but we recall them here for the sake of completeness.

Quantum spaces

Let M be a connected compact Kähler manifold, with fundamental 2-form ω ∈ Ω 2 (M, R). Assume M is endowed with a prequantum bundle L → M , that is a Hermitian holomorphic line bundle whose Chern connection ∇ has curvature 1 i ω. Let K → M be a Hermitian holomorphic line bundle. For every positive integer k, define the quantum space H k as:

H k = H 0 (M, L k ⊗ K) = holomorphic sections of L k ⊗ K .
The space H k is a subspace of the space L 2 (M, L k ⊗ K) of sections of finite L 2 -norm, where the scalar product is given by

ϕ, ψ = M h k (ϕ, ψ)µ M
with h k the hermitian product on L k ⊗ K induced by those of L and K, and µ M the Liouville measure on M . Since M is compact, H k is finite dimensional, and is thus given a Hilbert space structure with this scalar product.

Geometric notations

Unless otherwise mentioned, "smooth" will always mean C ∞ , and a section of a line bundle will always be assumed to be smooth. The space of sections of a bundle E → M will be denoted by Γ(M, E).

Let L P → P and L N → N be two prequantum bundles over Kähler manifolds, whose fundamental 2-forms are denoted by ω P and ω N . Denote by p 1 and p 2 the projections of P ×N on each factor, and

L P L N = p * 1 L P ⊗ p * 2 L N ; then, if P × N is endowed with the symplectic form p * 1 ω P + p * 2 ω N , L P L N → P × N is a prequantum bundle.
Let P op be the manifold P endowed with the symplectic form -ω P and the (almost) complex structure opposed to the one of P , and let L -1 P be the inverse (dual) bundle of L P with induced hermitian and holomorphic structure and connection; then L -1 P → P op is a prequantum bundle. If k is a positive integer, we can identify the Schwartz kernel of an operator T : Γ(P,

L k P ) → Γ(N, L k N ) to a section of L k N L -k P → N × P op via the following formula: T s(x) = P T (x, y).s(y)µ P (y),
where µ P is the Liouville measure on the manifold P .

Admissible and negligible sequences

Let M be a compact connected Kähler manifold. Let (s k ) k≥1 be a sequence such that for each k, s k belongs to Γ(M, L k ⊗ K). We say that (s k ) k≥1 is

• admissible if for every positive integer , for every vector fields X 1 , . . . , X on M , and for every compact set C ⊂ M , there exist a constant c > 0 and an integer N such that

∀m ∈ C ∇ X 1 . . . ∇ X s k (m) ≤ ck N ,
• negligible if for every positive integers and N , for every vector fields X 1 , . . . , X on M , and for every compact set C ⊂ M , there exists a constant c > 0 such that

∀m ∈ C ∇ X 1 . . . ∇ X s k (m) ≤ ck -N .
We say that (s k ) k≥1 is negligible over an open set U ⊂ M if the previous estimates hold for every compact subset of U . We denote by O(k

Toeplitz operators

Let Π k be the orthogonal projector of

L 2 (M, L k ⊗ K) onto H k . A Toeplitz operator is any sequence (T k : H k → H k ) k≥1 of operators of the form T k = Π k M f (.,k) + R k (2) where f (., k) is a sequence of C ∞ (M ) with an asymptotic expansion f (., k) = ≥0 k -f for the C ∞ topology, M f (.,k)
is the operator of multiplication by f (., k) and R k is a smoothing operator. They are the semiclassical analogue of the Toeplitz operators studied by Boutet de Monvel and Guillemin in [START_REF] Boutet De Monvel | The spectral theory of Toeplitz operators[END_REF].

We recall the following essential theorem about Toeplitz operators, which is a consequence of the works of Boutet de Monvel and Guillemin [START_REF] Boutet De Monvel | The spectral theory of Toeplitz operators[END_REF] (see also [START_REF] Bordemann | Toeplitz quantization of Kähler manifolds and gl(N ), N → ∞ limits[END_REF][START_REF] Charles | Berezin-Toeplitz operators, a semi-classical approach[END_REF][START_REF] Guillemin | Star products on compact pre-quantizable symplectic manifolds[END_REF]).

Theorem 2.1. The set T of Toeplitz operators is a star algebra whose identity is (Π k ) k≥1 . The contravariant symbol map

σ cont : T → C ∞ (M )[[ ]]
sending T k into the formal series ≥0 f is well defined, onto, and its kernel is the ideal consisting of O(k -∞ ) Toeplitz operators. More precisely, for any integer ,

T k = O(k -) if and only if σ cont (T k ) = O( ).
We will mainly work with the normalized symbol

σ norm = Id + 2 ∆ σ cont
where ∆ is the holomorphic Laplacian acting on C ∞ (M ); unless otherwise mentioned, when we talk about a subprincipal symbol, this refers to the normalized symbol. This symbol has the good property that, if T k and S k are Toeplitz operators with respective principal symbols t 0 and s 0 , then

σ norm (T k S k ) = t 0 s 0 + 2i {t 0 , s 0 } + O( 2 ).
Finally, we will need to apply functional calculus to Toeplitz operators.

Proposition 2.2 ([6]

). Let T k be a self-adjoint Toeplitz operator with symbol ≥0 t and g be a function of C ∞ (R, C). Then g(T k ) is a Toeplitz operator with principal symbol g(t 0 ).

Toeplitz operators on the complex plane

Bargmann spaces

We consider the Kähler manifold C R 2 with coordinates (x, ξ), standard complex structure and symplectic form ω 0 = dξ ∧dx. Let L 0 = R 2 ×C → R 2 be the trivial fiber bundle with standard hermitian metric h 0 and connection ∇ 0 with 1-form 1 i α, where α u (v) = 1 2 ω 0 (u, v); endow L 0 with the unique holomorphic structure compatible with h 0 and ∇ 0 . For every positive integer k, the quantum space that we consider is

H 0 k = H 0 (R 2 , L k 0 ) ∩ L 2 (R 2 , L k 0 );
this means that in this case, we make the arbitrary choice that the auxiliary line bundle K is the trivial bundle with flat connection. These spaces coincide with Bargmann spaces [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform[END_REF][START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform. Part II. A family of related function spaces. Application to distribution theory[END_REF], which are spaces of square integrable functions with respect to a Gaussian weight. More precisely, we choose the holomorphic coordinate z = x-iξ √ 2 and note

B k = f ψ k ; f : C → C holomorphic, R 2 |f (z)| 2 exp(-k|z| 2 ) dλ(z) < +∞ with ψ : C → C, z → exp -1 2 |z| 2 , ψ k its k-th
(z) = k n+1 2πn! z n ψ k , is an orthonormal basis of B k .
We denote by Π 0 k the orthogonal projector from L 2 (R 2 , L k 0 ) onto B k .

Admissible and negligible sequences

Since we will only deal with C ∞ sections, we can adopt the same definitions for admissible and negligible sequences as in the previous section.

Toeplitz operators

To consider Toeplitz operators acting on Bargmann spaces without raising technical issues, we could only work with operators with compactly supported kernels. However, we would miss the simple case of the harmonic oscillator. So we need to introduce symbol classes, very similar to the ones used when dealing with -pseudodifferential operators (see for instance [START_REF] Dimassi | Spectral asymptotics in the semi-classical limit[END_REF]). The proofs of the results of this part are collected in the appendix. Let d be a positive integer. For u in C d , set m(u) = 1 + u 2 1 2 . For every integer j, we define the symbol class S d j as the set of sequences of functions of C ∞ (C d ) which admit an asymptotic expansion of the form a(., k) = ≥0 k -a in the sense that

• ∀ ∈ N ∀α, β ∈ N 2d ∃ C ,α,β > 0 |∂ α z ∂ β z a | ≤ C ,α,β m j , • ∀L ∈ N * ∀α, β ∈ N 2d ∃ C L,α > 0 ∂ α z ∂ β z a -L-1 =0 k -a ≤ C L,α,β k -L m j .
Set S d = j∈Z S d j . Now, let a(., k) be a symbol in S 1 j , and consider the operator

A k = Op(a(., k)) = Π 0 k M a(.,k) Π 0 k (3)
acting on the subspace

S k = ϕ ∈ B k ; ∀j ∈ N sup z∈C |ϕ(z)|(1 + |z| 2 ) j/2 < +∞ of B k .
As shown in [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform. Part II. A family of related function spaces. Application to distribution theory[END_REF], S k corresponds to the Schwartz space via a particular unitary mapping between L 2 (R) and B k , the Bargmann transform. It is easily seen that

A k sends S k into S k ; it is even continuous S k → S k . Note that if j = 0, then A k is bounded B k → B k
, and its norm is lower than sup |a(., k)|.

Let t be the section of L 0 → R 2 with constant value 1. Let F 0 be the section of L 0 L -1 0 given by

F 0 (z 1 , z 2 ) = exp - 1 2 |z 1 | 2 + |z 2 | 2 -2z 1 z2 t(z 1 ) ⊗ t -1 (z 2 ), or equivalently, if u = (x, ξ) where z = 1 √ 2 (x -iξ), F 0 (u, v) = exp - 1 4 u -v 2 - i 2 ω 0 (u, v) t(u) ⊗ t -1 (v).
Adapting the result of section 1.c of [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform[END_REF], with the good normalization for the weight defining our Bargmann spaces, we have the following:

Proposition 3.2. Π 0 k admits a Schwartz kernel given by k 2π F k 0 .
In the rest of the paper, we will use the same letter to designate an operator and its Schwartz kernel. This proposition allows us to compute the Schwartz kernel of any Toeplitz operator.

Lemma 3.3. Let a(., k) be a symbol in S 1 j ; then A k = Op(a(.

, k)) admits a Schwartz kernel given by

A k (z 1 , z 2 ) = k 2π exp - k 2 |z 1 | 2 + |z 2 | 2 -2z 1 z2 ã(z 1 , z 2 , k) +R k exp -Ck|z 1 -z 2 | 2 , ( 4
)
where ã(., ., k) belongs to S 2 j , R k is negligible and C is some positive constant. Moreover, one has

ã(z, z, k) = exp k -1 ∆ a (z, k) ( 5
)
where

∆ = ∂ ∂z ∂ ∂z is the holomorphic Laplacian acting on C ∞ (C 2 )
, in the sense that the asymptotic expansion of ã(., ., k) is obtained by applying the formal asymptotic expansion of the operator exp k -1 ∆ to the asymptotic expansion of a(., k). This leads us to the following definition.

Definition 3.4. A Toeplitz operator is an operator from

S k to S k of the form Π 0 k M a(.,k) Π 0 k + S k , (6)
where a(., k) is a symbol in S 1 and the kernel of S k satisfies

S k (z 1 , z 2 ) = R k (z 1 , z 2 ) exp -Ck|z 1 -z 2 | 2 (7)
with R k negligible and C some positive constant. As in the compact case,

σ cont (A k ) = ≥0 a is called the contravariant symbol of A k .
We denote by T j the set of Toeplitz operators with contravariant symbol belonging to S 1 j .

In fact, lemma 3.3 defines the covariant symbol of A k :

σ cov (A k )(z) = ≥0 ã (z, z).
The following lemma gives an important property of the latter.

Lemma 3.5. If the covariant symbol of A k vanishes, then the Schwartz kernel of A k is of the form (7).

As a corollary of lemmas 3.3 and 3.5, we obtain the stability under composition of the set of Toeplitz operators. Corollary 3.6. Let A k ∈ T j and B k ∈ T j be two Toeplitz operators. Then C k = A k B k belongs to T j+j ; more precisely, its contravariant symbol is given by

σ cont (C k )(z) = exp - ∂ ∂z 1 ∂ ∂ z2 σ cont (A k )(z 1 )σ cont (B k )(z 2 ) |z 1 =z 2 =z , ( 8
) in the same sense as in lemma 3.3.

Define the normalized symbol as in the compact case:

σ norm = Id + 2 ∆ σ cont .
From formula (8), we find that

σ norm (A k B k ) = a 0 b 0 + 2i {a 0 , b 0 } + O( 2 ), as expected.
Definition 3.7. A Toeplitz operator A k ∈ T j is said to be elliptic at infinity if there exists some c > 0 such that for z in C,

|σ cont (A k )(z)| ≥ c(1 + |z| 2 ) j 2 .
Adapting proposition 12 of [START_REF] Charles | Berezin-Toeplitz operators, a semi-classical approach[END_REF] and theorem 39 of [START_REF] De Verdière | Méthodes semi-classiques et théorie spectrale[END_REF], one can show the following: Proposition 3.8 (Functional calculus). If A k belongs to T j for some j ≥ 1, is essentially self-adjoint and elliptic at infinity and if η : R → R is a compactly supported C ∞ function, then η(A k ) belongs to T j for every j < 0.

Fourier integral operators

The aim of this section is to construct microlocally unitary operators between H k and B k , given a local symplectomorphism χ from M to R 2 . In [START_REF] Boutet De Monvel | The spectral theory of Toeplitz operators[END_REF], Boutet de Monvel and Guillemin introduced Fourier integral operators in the homogeneous Toeplitz setting. In the semiclassical Toeplitz theory, such operators between compact manifolds have been used by Charles [START_REF] Charles | Quasimodes and Bohr-Sommerfeld conditions for the Toeplitz operators[END_REF][START_REF] Charles | Semi-classical properties of geometric quantization with metaplectic correction[END_REF], but some difficulties arise when dealing with a non compact manifold. Nevertheless, the ideas, based on Lagrangian sections, are very similar.

Let

χ : Ω 1 ⊂ M → Ω 2 ⊂ R 2 be a symplectomorphism between the open sets Ω 1 and Ω 2 . Then the graph Λ χ = {(u, χ(u)); u ∈ Ω 1 } ⊂ Ω 1 × Ω 2
of χ is a Lagrangian submanifold of the product M ×C op . As in the previous section, let t be the section of L 0 → R 2 with constant value 1. By definition of the connection on L 0 , we have ∇ 0 t = 1 i α ⊗ t, where α is the primitive of ω 0 given by α u (v) = 1 2 ω 0 (u, v). The following lemma is elementary.

Lemma 4.1.

Taking Ω 1 smaller if necessary, we can find a local gauge s of

L → Ω 1 such that ∇s = 1 i χ * α ⊗ s.
We consider the section t Λχ of L L -1 0 over Λ χ given by

t Λχ (u, χ(u)) = s(u) ⊗ t -1 (χ(u)).
Thanks to proposition 2.1 of [START_REF] Charles | Quasimodes and Bohr-Sommerfeld conditions for the Toeplitz operators[END_REF], we can build a local section

E of L L -1 0 → Ω 1 × Ω op 2 such that • E is equal to t Λχ on Λ χ ,
• for every holomorphic vector field Z on Ω 1 × Ω op 2 , the covariant derivative of E with respect to Z is zero modulo a section vanishing to infinite order along Λ χ , and this section is unique modulo a section vanishing to infinite order along Λ χ . Furthermore, we can always assume that E < 1 outside Λ χ , and we will make this assumption until the end of this article.

We consider a sequence of functions of C ∞ (Ω 1 × Ω op 2 ), which admits an asymptotic expansion ≥0 k -a for the C ∞ topology whose coefficients are all supported in a fixed (independent of k)

compact set C ⊂ Ω 1 × Ω 2 . Let S k be the local section of (L k ⊗ K) L -k 0 given by S k (u, v) = k 2π E k (u, v)a(u, v, k),
and consider the operator S k defined on Γ(C, L k 0 ) by

(S k φ)(u) = R 2 S k (u, v).φ(v) dλ(v),
which makes sense since S k (., .)

vanishes outside Ω 1 × Ω 2 . Proposition 4.2. The operator R k = S k Π 0 k maps B k into Γ(M, L k ⊗ K).
Proof. We must show that if ϕ k is a smooth square integrable section of

L k 0 → C, then S k Π 0 k ϕ k is a smooth section of L k ⊗ K → M . It is enough to show that the Schwartz kernel of S k Π 0
k and its derivatives with respect to the first variable are rapidly decreasing in the second variable. Let R k be this kernel; one has

R k (u, v) = p 2 (C) S k (u, w).Π 0 k (w, v)dw, with p 2 the projection from M × C to C. So R k (u, v) = k 2π 2 p 2 (C) f (u, v, w, k)E k (u, w).t k (w) ⊗ t -k (v)dw with f (u, v, w, k) = a(u, w, k) exp -k 4 w -v 2 -ik 2 ω 0 (w, v) . This implies the estimate R k (u, v) ≤ k 2π 2 p 2 (C) |a(u, w, k)| exp - k 4 w -v 2 E k (u, w) dw.
Since E ≤ 1 and a(., ., k) is bounded by some constant c k > 0, this yields

R k (u, v) ≤ k 2π 2 c k p 2 (C) exp - k 4 w -v 2 dw. Using w -v 2 ≥ v 2 -2 w v , we obtain R k (u, v) ≤ k 2π 2 c k exp - k 4 v 2 p 2 (C) exp k 2 w v dw.
Finally, an upper bound for the integral that appears in this inequality is πr 2 exp k 2 r v where p 2 (C) is included in the closed ball of radius r and centered at the origin. This allows us to conclude that

R k (u, v) ≤ k 2π 2 c k πr 2 exp - k 4 v 2 + kr 2 v .
The same kind of estimates hold for the successive derivatives of R k with respect to u; we prove them by differentiating under the integral sign.

Unfortunately, R k has no reason to map holomorphic sections to holomorphic sections; to fix this problem, we set

T k = Π k R k ; defined in this way, T k is an operator from B k to H k . Proposition 4.3. The Schwartz kernel of T k reads T k (u, v) = k 2π E k (u, v)b(u, v, k) + O(k -∞ ) (9)
with b(., ., k) a sequence of smooth functions which admits an asymptotic expansion b(., ., k

) = ≥0 k -b for the C ∞ topology satisfying b 0 (u, χ(u), k) = µ(u)a 0 (u, χ(u), k)
where µ is a smooth, nowhere vanishing function which depends only on the section E.

The proof is the same as the proof of proposition 4.2 of [START_REF] Charles | Quasimodes and Bohr-Sommerfeld conditions for the Toeplitz operators[END_REF]; it is based on an application of the stationary phase lemma.

An operator V k : B k → H k admitting a Schwartz kernel of the form of equation ( 9) and satisfying Π k V k Π 0 k = V k will be called a Fourier integral operator associated to the sequence b(., ., k); let us denote by FI(χ) the set of such operators. We define the full symbol map

σ : FI(χ) → C ∞ (M )[[ ]], V k → ≥0 b (u, χ(u)).
One can show that its kernel consists of smoothing operators. In the same way, we define FI(χ -1 ) : H k → B k . The following property is another application of the stationary phase lemma. Proposition 4.4. Let R k ∈ FI(χ) and S k ∈ FI(χ -1 ) with respective principal symbols r 0 (u, χ(u)) and s 0 (v, χ -1 (v)). Then there exists a smooth nowhere vanishing function ν : R 2 → R such that for every Toeplitz operator

T k on M with principal symbol t 0 , S k T k R k is a Toeplitz operator on R 2 with principal symbol equal to ν(v)s 0 (v, χ -1 (v))t 0 (χ -1 (v))r 0 (χ -1 (v), v) on Ω 2 .
To conclude this section, we prove that we can find microlocally unitary operators mapping B k to H k in the following sense.

Proposition 4.5. There exists a Fourier integral operator U

k : B k → H k such that • U * k U k ∼ Π 0 k on Ω 2 , • U k U * k ∼ Π k on Ω 1 ,
reducing

Ω
k = S k P k ; then U (0) * k U (0) k is a Toeplitz operator on R 2 , with principal symbol ν(v) |p 0 (v)| 2 s 0 χ -1 (v), v 2 . Since ν(v) and s 0 (s, χ -1 (v)) vanish in no point
v, one can choose p 0 such that this principal symbol is equal to 1. Doing so,

U (0) * k U (0) k
has the same principal symbol as Π 0 k , so there exists a Toeplitz operator R (0)

k such that U (0) * k U (0) k ∼ Π 0 k + k -1 R (0) k
on Ω 2 . From now on, when there is no ambiguity, the equality between operators will mean microlocal equality on Ω 2 . Let n ∈ N and assume that there exists an operator

U (n) k : B k → H k and a Toeplitz operator R (n) k (with principal symbol r n ) such that U (n) * k U (n) k = Π 0 k + k -(n+1) R (n) k .
Let T k be a Toeplitz operator on R 2 with principal symbol t 0 , and set

U (n+1) k = U (n) k Π 0 k + k -(n+1) T k . One has U (n+1) * k U (n+1) k = Π 0 k + k -(n+1) T * k + R (n) k + T k + k -(n+2) R (n+1) k with R (n+1) k
a Toeplitz operator. This implies that if we choose t 0 such that 2 (t 0 ) = -r n , then

U (n+1) * k U (n+1) k = Π 0 k + k -(n+2) R (n+1) k .
So we can construct the operators

U (n) k
by induction; it remains to apply Borel's summation lemma to find the desired operator U k .

Composing such a U k by U k on the left and U * k on the right gives:

(U k U * k ) 2 ∼ U k U * k on Ω 1 .
Since U k U * k is an elliptic Toeplitz operator on Ω 1 (its principal symbol vanishes nowhere), it has a microlocal inverse at each point of Ω 1 ; so the preceding equation yields

U k U * k ∼ Π k on Ω 1 .
Of course, such operators satisfy the analogue of Egorov's theorem:

Proposition 4.6. If U k is as above, then, for every Toeplitz operator

T k on M with principal symbol t 0 , S k = U * k T k U k is a Toeplitz operator on R 2 with principal symbol equal to t 0 • χ -1 on Ω 2 .
For a proof, see [12, theorem 47]. The action of a Fourier integral operator at the subprincipal level is much more complicated to compute. Denote by σ 0 (U k ) the principal symbol of U k , and by γ the 1-form on Λ χ such that

∇ Hom(C,K) σ 0 (U k ) = - 1 i γ ⊗ σ 0 (U k )
endowing C with the trivial connection and K with the one inherited from L 1 and δ. Now, notice that the symplectomorphism χ brings the complex structure of C op to a positive complex structure j on M . Introduce the section Ψ of Hom(Ω 1,0 (C), Ω 1,0

j (M )) |Λχ → Λ χ such that for all α ∈ Λ 1,0 (T Ω op 2 ) * and β ∈ Λ 1,0 (T Ω 1 ) * , Ψ(α) ∧ β = (χ * α) ∧ β.
This map is well-defined because the sesquilinear pairing Λ 1,0

j (T m Ω 1 ) * × Λ 1,0 (T m Ω 1 ) * → C, (α, β) → (α ∧ β)/ω m is non-degenerate. Let δ be the 1-form on Λ χ such that ∇Ψ = δ ⊗ Ψ
where ∇ is the connection induced by the Chern connections of Ω 2,0 (C) and Ω 2,0 (M ). In [9, Theorem 3.3], Charles derived the following formula.

Theorem 4.7. With the same notations as in the previous proposition and denoting by t 1 the subprincipal symbol of T k , the subprincipal symbol s 1 of S k is given by:

s 1 (u) = t 1 (m) + γ (m,u) - 1 2 δ (m,u) , X t 0 (m), ((χ -1 ) * X t 0 )(u)
for u ∈ R 2 and m = χ -1 (u).

Microlocal normal form

The local model

Our local model will be the realization of the quantum harmonic oscillator in the Bargmann representation:

Q k = 1 k z ∂ ∂z + 1 2 , with domain C[z] which is dense is B k .
The following lemma is easily shown. Lemma 5.1. Q k is an essentially self-adjoint Toeplitz operator with normalized symbol q 0 . Moreover, the spectrum of

Q k is Sp(Q k ) = k -1 n + 1 2
, n ∈ N .

A symplectic Morse lemma

Let (N, ω) be a two-dimensional symplectic manifold and f a function of C ∞ (N, R). Assume f admits an elliptic critical point at n 0 ∈ N , with f (n 0 ) = 0. Replacing f by -f if necessary, we can assume that this critical point is a local minimum for f . Define

q 0 : R → R 2 , (x, ξ) → 1 2 (x 2 + ξ 2 ).
The following theorem is well-known.

Theorem 5.2.

There exist a local symplectomorphism χ : (N, n 0 ) → (R 2 , 0) and a function g in C ∞ (R, R) satisfying g(0) = 0 and g (0) > 0, such that

f • χ -1 = g • q 0
where χ -1 is defined.

It can be viewed as a consequence of the isochore Morse lemma [START_REF] De Verdière | Le lemme de Morse isochore[END_REF] or of Eliasson's symplectic normal form theorem [START_REF] Eliasson | Normal forms for Hamiltonian systems with Poisson commuting integrals-elliptic case[END_REF], but this case is in fact easier than these two results.

Semiclassical normal form

We consider a self-adjoint Toeplitz operator A k on M ; its normalized symbol a(., ) = a 0 + a 1 + . . . is real-valued. Assume that the principal symbol a 0 admits a non-degenerate local minimum at m 0 ∈ M . Assume also that a 0 (m 0 ) = 0, so that a 0 takes positive values on a neighbourhood of m 0 . Hence, thanks to theorem 5.2, we get a neighbourhood Ω 1 of m 0 in M , a neighbourhood Ω 2 of 0 in R 2 , a local symplectomorphism χ : Ω 1 → Ω 2 and a function g 0 of C ∞ (R, R) with g 0 (0) = 0 and g 0 (0) > 0, such that:

a 0 • χ -1 = g 0 • q 0
on Ω 2 . We denote by f 0 the local inverse of g 0 . Our goal is to show: Theorem 5.3. There exist a Fourier integral operator U k : B k → H k and a sequence f (., k) of functions of C ∞ (R, R) which admits an asymptotic expansion in the C ∞ topology of the form f (., k) = ≥0 k -f , such that:

• U * k U k ∼ Π 0 k on Ω 2 , • U k U * k ∼ Π k on Ω 1 , • U * k f (A k , k)U k ∼ Q k on Ω 2 .
Proof. We consider an operator

U (0)
k satisfying the two first points (see the previous section). We will construct the operator that we seek by successive perturbations by unitary Toeplitz operators on B k . More precisely, we show by induction that for every positive integer n, there exist an operator

U (n) k : B k → H k satisfying the two first points, a sequence f (n) (., k) of functions of C ∞ (R, R) of the form f (n) (., k) = n =0 k -f , with f smooth, and a Toeplitz operator R (n) k acting on B k such that U (n) * k f (n) (A k , k)U (n) k = Q k + k -(n+1) R (n) k on Ω 2 .
The first step is as follows: by the results of the previous section, the operator

U (0) * k f 0 (A k )U (0)
k is a Toeplitz operator on B k , whose principal symbol is equal to f 0 • a 0 • χ -1 = q 0 on Ω 2 . Hence, there exists a Toeplitz operator R (0)

k on B k such that U (0) * k f 0 (A k )U (0) k = Q k + k -1 R (0) k .
We look for U

(1) k of the form U (0) k P k with P k a unitary Toeplitz operator on B k . Moreover, we choose f (1) (., k) = f 0 + k -1 θ 1 • f 0 with θ 1 a smooth function that remains to determine. Expanding, we get

U (1) * k f (1) (A k , k)U (1) k = P * k U (0) * k f 0 (A k )U (0) k P k +k -1 P * k U (0) * k (θ 1 •f 0 )(A k )U (0) k P k which yields U (1) * k f (1) (A k , k)U (1) k = P * k Q k + k -1 R (0) k P k +k -1 P * k U (0) * k (θ 1 •f 0 )(A k )U (0) k P k .
Consequently, we wish to have

P * k Q k + k -1 R (0) k P k + k -1 P * k U (0) * k (θ 1 • f 0 )(A k )U (0) k P k = Q k + k -2 R (1) k where R (1)
k is a Toeplitz operator; this amounts, remembering that P k is unitary, to

Q k P k + k -1 R (0) k P k + U (0) * k (θ 1 • f 0 )(A k )U (0) k P k = P k Q k + k -2 P k R (1) k which we can rewrite [Q k , P k ] + k -1 R (0) k P k + U (0) * k (θ 1 • f 0 )(A k )U (0) k P k = k -2 P k R (1)
k .

This will be true if and only if the subprincipal symbol of the operator on the left of the equality vanishes, that is to say

1 i {q 0 , p 0 } + p 0 (r 0 + θ 1 • q 0 ) = 0
where p 0 and r 0 stand for the respective principal symbols of P k and R (0) k . Set p 0 = exp(iϕ 0 ) with ϕ 0 a smooth, real-valued function (since P k is unitary). The previous equation then becomes

{ϕ 0 , q 0 } = r 0 + θ 1 • q 0 .
This equation is standard and it is well-known that it can be solved. We recall a method from [START_REF] Eliasson | Normal forms for Hamiltonian systems with Poisson commuting integrals-elliptic case[END_REF] to find θ 1 and ϕ 0 smooth such that it is satisfied, since we will need to know how to construct these in part 6.2. Consider the functions:

F (x, ξ) = - 1 2π 2π 0 r 0 (φ t q 0 (x, ξ)) dt and ϕ 0 (x, ξ) = - 1 2π 2π 0 t r 0 (φ t q 0 (x, ξ)) dt,
where φ t q 0 stands for the Hamiltonian flow of q 0 taken at time t: φ t q 0 (x, ξ) = (x cos t + ξ sin t, -x sin t + ξ cos t).

Since we integrate on a compact set and the flow φ t is smooth with respect to (x, ξ), both F and ϕ 0 are smooth. By construction, we have {F, q 0 } = 0. But we have the easy lemma Lemma 5.4. Let f be a function of C ∞ (R 2 , R) such that {f, q 0 } = 0. Then the function g such that

f = g • q 0 belongs to C ∞ (R, R).
So there exists a function θ 1 of C ∞ (R, R) such that F = θ 1 • q 0 . Integrating by parts, it is easy to show that {ϕ 0 , q 0 } = θ 1 • q 0 + r 0

The next steps are practically the same; indeed, let n ≥ 1 and assume that we have found U (n) k and f (n) (., k) satisfying the desired properties. We now look for

U (n+1) k of the form U (n) k (Π 0 k + k -n V k ) with V k a Toeplitz operator on B k such that Π 0 k + k -n V k is unitary. Furthermore, we write f (n+1) (., k) = f (n) (., k) + k -(n+1) θ n+1
• f 0 with θ n+1 an unknown smooth function. We want the existence of a Toeplitz operator R

(n+1) k such that U (n) * k f (n+1) (A k , k)U (n) k Π 0 k + k -n V k = (Π 0 k +k -n V k ) Q k + k -(n+2) R (n+1) k
which gives, expanding,

Q k + k -n Q k V k + k -(n+1) R (n) k + U (n) * k (θ n+1 • f 0 )(A k )U (n) k = Q k + k -n V k Q k + k -(n+2) R (n+1) k + k -(2n+1) V k R (n+1) k
.

Thus, we wish that

[Q k , V k ] + k -1 R (n) k + U (n) * k (θ n+1 • f 0 )(A k )U (n) k = 0;
this will be verified if and only if

1 i {q 0 , v 0 } + r n + θ n+1 • q 0 = 0
which is treated as before.

We conclude thanks to Borel's summation lemma (applied to both f (., k) and U k ).

Bohr-Sommerfeld conditions

Let A k be a self-adjoint Toeplitz operator on M as in the previous section. Moreover, assume that a 0 (m 0 ) is a global minimum of the principal symbol a 0 , and that m 0 is the unique point of M with this property. This implies that there exists E 0 > 0 such that for every E ≤ E 0 , the level set a -1 0 (E) is connected and contained in Ω 1 .

The maximum norm A k ∞ of A k tends to the L ∞ -norm a 0 ∞ of a 0 as k goes to infinity [START_REF] Bordemann | Toeplitz quantization of Kähler manifolds and gl(N ), N → ∞ limits[END_REF]; hence, for k large enough, the spectrum of A k is included in the set -E 1 , E 1 , where E 1 = a 0 ∞ + 1.

Statement of the result

Before stating the Bohr-Sommerfeld conditions, it is convenient to show that the sequence f (., k) can be inverted, and that its inverse still has a good asymptotic expansion. Lemma 6.1. For k large enough, the function f (., k) that appears in theorem 6.2 is a bijection from -E 1 , E 0 to its image; more precisely, it is strictly increasing. Moreover, the inverse sequence g(., k) admits an asymptotic expansion in the C ∞ topology of the form g(., k)

= ≥0 k -g + O(k -∞ ), uniformly on -E 1 , E 0 .
Proof. The first assertion follows from the mean value inequality

∀k ≥ 1 ∀E, Ẽ ∈ -E 1 , E 0 |f (E, k)-f ( Ẽ, k)| ≥ inf [-E 1 ,E 0 ] |f (., k)||E-Ẽ|
and the fact that f (., k) is bounded below by some positive constant. This implies that for k sufficiently large, f (., k) is strictly monotone on -E 1 , E 0 ; since f 0 (0) > 0, f (., k) is in fact strictly increasing. For the second part, the proof is once again based on Borel's summation lemma; it is done by induction thanks to Taylor's formula with integral remainder.

We can therefore introduce the sequences

E (j) k = g k -1 j + 1 2 , k , j ∈ N ( 10
)
for k large enough and for j such that k -1 j + 1 2 belongs to the set

f (-E 1 , k), f (E 0 , k) . Since g(., k) is also strictly increasing, the E (j) k are ordered: ∀j ∈ N, E (j) k < E (j+1) k .
We can be more precise; fix j ∈ N and write

E (j) k = g 0 k -1 j + 1 2 + k -1 g 1 k -1 j + 1 2 + O(k -2 ).
Then, applying Taylor's formula with integral remainder, we get

E (j) k = g 0 (0) =0 +k -1 g 1 (0) + j + 1 2 g 0 (0) + O(k -2 ). ( 11
)
One must be careful with this estimate: the O(k -2 ) remainder is no longer uniform with respect to j. Denote by λ

(1)

k ≤ λ (2) k ≤ . . . ≤ λ (j)
k ≤ . . . the eigenvalues of A k . The main result of this section is the following theorem. Theorem 6.2. There exists a positive integer k 0 ≥ 1 such that for every integer N ≥ 1 and for every E ≤ E 0 , there exist a constant C N > 0 such that for k ≥ k 0 :

λ (j) k ≤ E or E (j) k ≤ E ⇒ λ (j) k -E (j) k ≤ C N k -N . ( 12
)
Moreover, for k large enough, all the eigenvalues of A k smaller than E 0 are simple. In particular, we obtain an asymptotic expansion to every order for the eigenvalues of A k smaller than E 0 .

We will need the following lemma, based on the min-max principle.

Lemma 6.3 ([11, lemma 3.3]).

Let A and B be two self-adjoint operators acting, respectively, on the Hilbert spaces H and H, both bounded from below. Denote by Π A I the spectral projection of A on I and by λ A 1 ≤ λ A 2 ≤ . . . ≤ λ A j ≤ . . . the increasing sequence of eigenvalues below the essential spectrum of A; if there is a finite number j max of such eigenvalues, extend the sequence for j > j max by setting λ A j = λ A ess , where λ A ess is the infimum of the essential spectrum of A. Introduce the same notations for B. Suppose that there exist a bounded operator U : H → H , an interval I = (-∞, E], and constants C > 0, c ∈ (0, 1) such that U Π B I (H) ⊂ Dom(A) and

(U * AU -B)Π B I ≤ C and U * U Π B I -Π B I ≤ c.
Then, for all j such that λ B j ≤ E, one has

λ A j ≤ (λ B j + C) 1 + c 1 -c . Proof of theorem 6.2. Fix E in -E 1 , E 0 . Let J be an open neighbourhood of -E 1 , E such that the open set a -1 0 (J ) is contained in Ω 1
, and let η : R → R be a smooth function equal to 1 on -E 1 , E and 0 outside J . Consider the Toeplitz operator R k = η(A k ) and set

B k = (f (A k , k) -U k Q k U * k ) R k . By the choice of R k , the microsupport of B k is a subset of Ω 1 . Moreover, f (A k , k) is microlocally equal to U k Q k U * k on Ω 1 .
These two facts imply that B k is negligible; since M is compact, this yields that for every N ≥ 1, there exists a positive constant C N such that

B k ≤ C N k -N . Now, let Π f (A k ,k)
≤f (E,k) be the spectral projection associated to f (A k , k) and corresponding to the eigenvalues smaller than f (E, k).

If (λ, ϕ) is an eigencouple for A k with λ ≤ E, then R k ϕ = η(λ)ϕ = ϕ. This implies that for every φ in Π f (A k ,k) ≤f (E,k) (H k ), R k φ = φ, and consequently (f (A k , k) -U k Q k U * k ) Π f (A k ,k) ≤f (E,k) ≤ C N k -N .
Similarly, there exists c N > 0 such that

(Π k -U k U * k ) Π f (A k ,k) ≤f (E,k) ≤ c N k -N .
So lemma 6.3 shows that if f (λ

(j) k , k) ≤ f (E, k), the inequality k -1 j + 1 2 ≤ 1 + c N k -N 1 -c N k -N f (λ (j) k , k) + C N k -N
holds. So for k large enough (independently of E), we have k

-1 j + 1 2 ≤ f (E, k).
Now, let ρ be a smooth function equal to 1 on f 0 (-E 1 ), f 0 (E) and vanishing outside an open neighbourhood K of f 0 (-E 1 ), f 0 (E) such that q -1 0 (K) ⊂ Ω 2 . Thanks to proposition 3.8, we can consider the Toeplitz operator S k = ρ(Q k ), and set

C k = (U * k f (A k , k)U k -Q k ) S k .
Since S k belongs to every T j , j < 0, C k belongs to T 0 and is thus a bounded operator B k → B k . Moreover, by construction, it is negligible. Hence there exists a positive constant CN such that

C k ≤ CN k -N ; modifying C N if necessary, we can assume that CN is equal to C N . So, in- troducing the spectral projection Π Q k ≤f (E,k) corresponding to the eigenvalues of Q k smaller than f (E, k), the inequality (U * k f (A k , k)U k -Q k ) Π Q k ≤f (E,k) ≤ C N k -N
holds. Similarly, we have

U * k U k -Π 0 k Π Q k ≤f (E,k) ≤ c N k -N .
Hence, applying again lemma 6.3, we obtain that f (λ

(j) k , k) ≤ 1 + c N k -N 1 -c N k -N k -1 j + 1 2 + C N k -N
as soon as k -1 j + 1 2 ≤ f (E, k). This shows that if f (λ

(j) k , k) ≤ f (E, k), then f (λ (j) k , k) -k -1 j + 1 2 ≤ C N k -N
for some positive constant C N . Exchanging the roles of Q k and A k , and using lemma 6.1, this gives formula [START_REF] De Verdière | Méthodes semi-classiques et théorie spectrale[END_REF]. Using this result and the fact that there exists c > 0 such that for j ∈ N,

E (j+1) k -E (j)
k is equivalent to ck -1 , we obtain that the λ (j) k are simple for k large enough.

Computation of the principal and subprincipal terms

In order to exploit these results, it remains to compute a few first terms in the asymptotic expansion of the sequence f (., k). What we can do is relate the principal and subprincipal terms to the actions introduced in section 1.2. Proposition 6.4. Set I =]0, E 0 [. Then

f 0 = 1 2π c 0 , f 1 = 1 2π c 1 (13)
on I.

Proof. Let us first compute f 0 . Fix a level E in I. Let 1 i β be the 1-form describing locally the Chern connection on L; then c 0 (E) is given by

c 0 (E) = Γ E β.
Using the relation a 0 • χ -1 = g 0 • q 0 , we can then write

c 0 (E) = C E (χ -1 ) * β
where C E is the circle centered at the origin and with radius 2f 0 (E). Using Stokes' formula and the fact that χ is a symplectomorphism, this yields that c 0 (E) is the area of the disk bounded by C E , that is, if the orientation that we chose is the one giving the positive area (and this is what we will assume in the rest of this section)

c 0 (E) = 2πf 0 (E). (14)
Now, turn back to the proof of our normal form theorem 5.3, where f 1 is constructed from the subprincipal symbol r 0 of U

(0) * k f (A k , k)U (0)
k . By uniqueness of f 1 , instead of starting from any operator U (0) k , we can choose one with symbol u ⊗ v, where u is constant and v is a square root of Ψ. Doing so, we can compute r 0 thanks to theorem 4.7:

r 0 = (a 1 • χ -1 )( f 0 • a 0 • χ -1 ) -ν χ -1 (.) X f 0 •a 0 • χ -1
where ν is the local connection 1-form associated to ∇ L 1 . We have

f 1 = θ 1 • f 0 with θ 1 such that for all (x, ξ) in R 2 (θ 1 • q 0 ) (x, ξ) = - 1 2π 2π 0 r 0 φ t q 0 (x, ξ) dt
where φ t q 0 stands for the Hamiltonian flow of q 0 . Since

q 0 = f 0 • a 0 • χ -1 , this implies that for (x, ξ) in R 2 f 1 • a 0 • χ -1 (x, ξ) = - 1 2π 2π 0 (a 1 f 0 • a 0 ) χ -1 φ t q 0 (x, ξ) dt + 1 2π 2π 0 ν χ -1 (φ t q 0 (x,ξ)) X f 0 •a 0 χ -1 (φ t q 0 (x, ξ)) dt.
So, for m = m 0 in Ω 1 , we have

(f 1 • a 0 )(m) = - 1 2π 2π 0 (a 1 f 0 • a 0 ) χ -1 φ t q 0 (χ(m)) dt + 1 2π 2π 0 ν χ -1 (φ t q 0 (χ(m))) X f 0 •a 0 χ -1 φ t q 0 (χ(m)) dt; thus, if E = a 0 (m), f 1 (E) = 1 2π 2π 0 ν χ -1 (φ t q 0 (χ(m))) X f 0 •a 0 χ -1 φ t q 0 (χ(m)) dt - f 0 (E) 2π 2π 0 a 1 χ -1 φ t q 0 (χ(m)) dt. But χ -1 • φ t q 0 • χ is the Hamiltonian flow of q 0 • χ = f 0 • a 0 , so f 1 (E) = 1 2π 2π 0 ν φ t f 0 •a 0 (m) X f 0 •a 0 (φ t f 0 •a 0 (m)) dt - f 0 (E) 2π 2π 0 a 1 φ t f 0 •a 0 (m) dt.
Therefore, if T E is the period of the flow φ t a 0 along Γ E , we have T E = 2πf 0 (E) for E close to 0, and a change of variable gives

f 1 (E) = 1 2πf 0 (E) 2π 0 ν φ tf 0 (E) f 0 •a 0 (m) X f 0 •a 0 φ tf 0 (E) f 0 •a 0 (m) dt - 1 2π T E 0 a 1 φ tf 0 (E) f 0 •a 0 (m) dt
which yields, since the Hamiltonian vector field associated to f 0 • a 0 is X f 0 •a 0 = (f 0 • a 0 )X a 0 , and hence φ tf 0 (E)

f 0 •a 0 (m) = φ t a 0 (m): f 1 (E) = 1 2πf 0 (E) 2π 0 ν φ t a 0 (m) f 0 (E)X a 0 (φ t a 0 (m)) dt - 1 2π T E 0 a 1 φ t a 0 (m) dt =-Γ E κ E
and by linearity of ν

f 1 (E) = 1 2π 2π 0 ν φ t a 0 (m) X a 0 (φ t a 0 (m)) dt = Γ E ν + 1 2π Γ E κ E .
The right term of this equality is precisely equal to c 1 (E); so we have on I

c 1 = 2πf 1 . ( 15
)
To conclude, we can show that = 1 on I, so the result of theorem 6.2 matches the usual Bohr-Sommerfeld conditions on the set I of regular values.

First terms of the asymptotic expansion of the eigenvalues

Theorem 6.2 and formula [START_REF] Charles | Spectral asymptotics via the semiclassical Birkhoff normal form[END_REF] give an asymptotic expansion for the eigenvalues of A k smaller than E 0 . Fix j ∈ N; for k large enough, one has

λ (j) k = k -1 g 1 (0) + j + 1 2 g 0 (0) + O(k -2 ).
We can be more precise, since we know the value of g 1 (0): by definition of g(., k), we have g 1 = -(f 1 •g 0 )g 0 , and the computation made in the previous part leads to f 1 (0) = -a 1 (0)f 0 (0); consequently, g 1 (0) = a 1 (0) and

λ (j) k = k -1 a 1 (0) + j + 1 2 g 0 (0) + O(k -2 ). ( 16
)
But g 0 (0) = 1 f 0 (0) ; moreover, it is standard that the principal action c 0 is smooth even at the critical value E = 0. Hence, thanks to formula [START_REF] De Verdière | Le lemme de Morse isochore[END_REF], one has

λ (j) k = k -1   a 1 (0) + 2π j + 1 2 c 0 (0)   + O(k -2 ). ( 17
)
In particular, the gap between two consecutive eigenvalues is given by

λ (j+1) k -λ (j) k = 2πk -1 c 0 (0) + O(k -2 ). ( 18
)

An example on the torus

The aim of this section is to give numerical evidence for our results by investigating the case of a particular Toeplitz operator on the torus T of real dimension 2. One can find the details of the quantization of T in [START_REF] Charles | Knot state asymptotics I, AJ conjecture and abelian representations[END_REF], where the authors investigate some conjectures on knot states; let us briefly recall its main ingredients.

The setting

Endow R 2 with the linear symplectic form ω 0 and consider a lattice Λ with symplectic volume 4π. The Heisenberg group H = R 2 × U (1) with product (x, u).(y, v) = x + y, uv exp i 2 ω 0 (x, y) acts on the trivial bundle L 0 → R 2 , with action given by the same formula. This action preserves the prequantum data, and the lattice Λ injects into H; therefore, the fiber bundle L 0 reduces to a prequantum bundle L over T = R 2 /Λ. The action extends to the fiber bundle L k 0 by

(x, u).(y, v) = x + y, u k v exp ik 2 ω 0 (x, y)
and naturally induces an action

T * : Λ → End Γ R 2 , L k 0 , u → T * u .
The Hilbert space H k = H 0 (M, L k ) can naturally be identified to the space H Λ,k of holomorphic sections of L k 0 → R 2 which are invariant under the action of Λ, endowed with the hermitian product

ϕ, ψ = D ϕ ψ |ω 0 |
where D is the fundamental domain of the lattice. Furthermore, Λ/2k acts on H Λ,k . Let e and f be generators of Λ satisfying ω 0 (e, f ) = 4π; one can show that there exists an orthonormal basis (ψ ) ∈Z/2kZ of H Λ,k such that

∀ ∈ Z/2kZ T * e/2k ψ = w ψ T * f /2k ψ = ψ +1 (19)
with w = exp iπ k . The sections ψ can be expressed in terms of Θ functions. Set M k = T * e/2k and L k = T * f /2k . Let (q, p) be coordinates on R 2 associated to the basis (e, f ) and [q, p] be the equivalence class of (q, p). Both M k and L k are Toeplitz operators, with respective principal symbols [q, p] → exp(2iπp) and [q, p] → exp(2iπq), and vanishing subprincipal symbols. Consequently

A k = M k + M -1 k + L k + L -1
k is a Toeplitz operator on T with principal symbol

a 0 (q, p) = 2 (cos(2πp) + cos(2πq)) (20)
and vanishing subprincipal symbol. Its matrix in the basis (ψ

) ∈Z/2kZ is              2α 0 1 0 . . . 0 1
             where α = cos π k .
The principal symbol a 0 is known as Harper's Hamiltonian. It admits a global minimum at m 0 = [1/2, 1/2], with a 0 (m 0 ) = -4. Figure 1 is a contour plot of this function on the fundamental domain. In figure 2, we derived numerically the spectrum of A k . Figure 3 shows the modulus of the eigenfunction associated to the eigenvalue closest to a prescribed level E.

In this situation, we can express c 0 (0), and so, by equation ( 17), the first eigenvalues of A k . Indeed, the symplectic form on T is ω = 4πdp ∧ dq and the hessian of a 0 at m 0 = [1/2, 1/2] is given by:

d 2 a 0 (m 0 ) = 8π 2 I 2 so it is easy to obtain c 0 (0) = 1.
Consequently, the first eigenvalues are given by:

λ (j) k = -4 + 2πk -1 j + 1 2 + O(k -2 ). (21)
This is exactly what our simulations show; we plotted the eigenvalues located in a window of length 20πk -1 around the minimum (so we expect to see about ten eigenvalues) and the result can be seen in figure 4.

Figures

A k (z 1 , z 2 ) = C Π k (z 1 , z 3 )a(z 3 , k)Π k (z 3 , z 2 ) dλ(z 3 ) which is equal to k 2π 2 C exp - k 2 |z 1 | 2 + |z 2 | 2 + 2|z 3 | 2 -2z 1 z3 -2z 3 z2 a(z 3 , k) dλ(z 3 ).
(22) This can be written

A k (z 1 , z 2 ) = k 2π 2 C exp(ikφ(z 1 , z 2 , z 3 ))a(z 3 , k) dλ(z 3 )
with the phase φ given by

φ(z 1 , z 2 , z 3 ) = i 2 |z 1 | 2 + |z 2 | 2 + 2|z 3 | 2 -2z 1 z3 -2z 3 z2 .
It is more convenient to use real variables: let u j = (x j , ξ j ) be the point of

R 2 corresponding to z j = 1 √ 2 (x j -iξ j ). The phase φ reads φ(u 1 , u 2 , u 3 ) = i 4 u 1 -u 3 2 + u 3 -u 2 2 + 2iω 0 (u 1 -u 2 , u 3 ) .
Using the identity

u 1 -u 3 2 + u 3 -u 2 2 = 1 2 u 1 -u 2 2 + 2u 3 -u 1 -u 2 2 ,
we can rewrite φ as

φ(u 1 , u 2 , u 3 ) = i 8 u 1 -u 2 2 + ϕ(u 1 , u 2 , u 3 ) with ϕ(u 1 , u 2 , u 3 ) = i 4 1 2 2u 3 -u 1 -u 2 2 + 2iω 0 (u 1 -u 2 , u 3 ) . So we have A k (u 1 , u 2 ) = exp -k 8 u 1 -u 2 2 I k (u 1 , u 2 ) with I k (u 1 , u 2 ) = k 2π 2 R 2 exp(ikϕ(u 1 , u 2 , u 3 ))a(u 3 , k) dλ(u 3 );
a change of variable finally gives

I k (u 1 , u 2 ) = k 2π 2 R 2 exp(ikψ(u 1 , u 2 , u 3 )) a u 3 + u 1 + u 2 2 , k dλ(u 3 ), with ψ(u 1 , u 2 , u 3 ) = i 2 u 3 2 + iω 0 (u 1 -u 2 , u 3 ) + iω 0 (u 1 , u 2 ) .
To evaluate this integral, we cannot directly use the stationary phase lemma because a(., k) may not be compactly supported; we have to adapt it. In order to do so, we start by computing the critical locus

C ψ = (u 1 , u 2 , u 3 ) ∈ C 3 ; d u 3 ψ(u 1 , u 2 , u 3 ) = 0 and (ψ)(u 1 , u 2 , u 3 ) = 0 of ψ. It is clear that (ψ)(u 1 , u 2 , u 3 ) = 0 if
and only if u 3 = 0; moreover, the derivative of ψ with respect to u 3 is given by

d u 3 ψ(u 1 , u 2 , u 3 ) = i u 3 , . - 1 2 ω 0 (u 1 -u 2 , .), hence C ψ is the set of (u, u, 0), u ∈ R 2 . Now, consider χ ∈ C ∞ (R 2 , R + ) equal
to 1 in the set { x 3 ≤ δ} for some δ > 0 and with compact support, and

decompose I k as I k = k 2π 2 (J k + K k ), with J k (u 1 , u 2 ) = R 2 exp(ikψ(u 1 , u 2 , u 3 )) a u 3 + u 1 + u 2 2 , k χ(u 3 ) dλ(u 3 )
and K k (u 1 , u 2 ) equal to the same integral replacing χ by 1 -χ. First, we show that K k is negligible. Choose R > 0 and consider the points (u 1 , u 2 ) belonging to the ball of R 4 centered at the origin and of radius R. Writing the integral in K k in polar coordinates, we have Performing successive integration by parts, it is easy to prove that for every positive integer f . Furthermore, using the facts that ∂Ψ ∂ρ (ρ, θ) = iρ -1 2 (x 1 -x 2 ) sin θ + 1 2 (y 1 -y 2 ) cos θ, D N A is equal to a linear combination of terms of the form ∂ p Ψ ∂ρ p ∂ q A ∂ρ q and a(., k) belongs to S 1 j , we have the estimate This shows that K k ≤ c N k -N on the ball of radius R for some c N > 0. We treat the derivatives of K k in the same way. It remains to estimate J k . Since the second derivative of ψ with respect to u 3 is equal to iId, and since

K k (u 1 , u 2 ) = 2π
N K k (u 1 , u 2 ) = k -N
ψ(u 1 , u 2 , u 3 ) = i 8 u 1 -u 2 2 - 1 2 ω 0 (u 1 , u 2 ) + 1 4 (i(ξ 1 -ξ 2 ) + 2x 3 ) ∂ x 4 ψ(u 1 , u 2 , u 3 ) + 1 4
(i(x 2 -x 1 ) + 2ξ 3 ) ∂ ξ 4 ψ(u 1 , u 2 , u 3 ), the stationary phase lemma [18, section 7.7] gives

J k (u 1 , u 2 ) = 2π k exp i 8 u 1 -u 2 2 - 1 2 ω 0 (u 1 , u 2 ) ã(u, u, k) + O(k -∞ ),
where ã(., ., k) belongs to S 2 j ; the coefficients ã (u 1 , u 2 , k) of its asymptotic expansion, which we do not write here, are linear combinations of derivatives of the a m , m ≥ 0, evaluated at u 1 +u 2 2 . However, the values of a(., ., k) along the diagonal of C 2 can be easily computed, because a number of terms vanish: for z in C, we have

ã(z, z, k) =   ≥0 k - ! ∆ a   (z, k) = exp k -1 ∆a (z, k).
Putting this together with the fact that K k is negligible, we obtain the result.

Proof of lemma 3.5. Formula (22) shows that the kernel A k is a holomorphic section of L k 0 L -k 0 ; differentiating equation ( 4), this implies that the Thanks to these holomorphy conditions, the fact that ã vanishes on the diagonal implies that it vanishes to all order along the diagonal. We can easily adapt lemma 1 of [START_REF] Charles | Berezin-Toeplitz operators, a semi-classical approach[END_REF] to show that this yields the negligibility of expk 4 |z 1 -z 2 | 2 ã(z 1 , z 2 , k). Injecting this in formula (4) gives the result.

Proof of corollary 3.6. The kernel of C k reads

C k (z 1 , z 2 ) = C A k (z 1 , z 3 )B k (z 3 , z 2 ) dλ(z 3 ).
Using the representations of A k and B k given by lemma 3.3, this yields Using the same technique as in the proof of lemma 3.3, we show that

C k (z 1 , z 2 ) = k 2π
C k (z 1 , z 2 ) = k 2π exp - k 2 |z 1 | 2 + |z 2 | 2 -2z 1 z2 c(z 1 , z 2 , k) +R k exp -C k|z 1 -z 2 | 2
with C > 0, R k negligible, and c(., ., k) ∈ S 2 j+j . Now, consider the function č(., k) defined by č(z, k) = c(z, z, k) for z in C, and put c(., k) = exp -k -1 ∆ č (., k). Then c(., k) belongs to S 1 j+j and, by lemma 3.3, the Toeplitz operator D k = Op(c(., k)) admits a Schwartz kernel of the form

D k (z 1 , z 2 ) = k 2π exp - k 2 |z 1 | 2 + |z 2 | 2 -2z 1 z2 d(z 1 , z 2 , k) +R k exp -C k|z 1 -z 2 | 2 ,
where d(., ., k) belongs to S 2 j , R k is negligible, C is a positive constant and for every z in C, d(z, z, k) = c(z, z, k). |u=v=z up to a negligible term, which was to be proved.

Figure 1 :Figure 2 :

 12 Figure 1: A few level sets of a 0

Figure 3 :

 3 Figure 3: E = a 0 (0.7, 0.6), k = 500

Figure 4 : 3

 43 Figure 4: Eigenvalues in [-4, -4 + 20πk -1 ]; in red squares, the eigenvalues of A k obtained numerically; in blue diamonds, the theoretical eigenvalues up to order k -2 (i.e. forgetting the O(k -2 ) in formula (21))

  (ρ, θ)) A(ρ, θ, k) dρ dθ where Ψ(ρ, θ) = ψ(u 1 , u 2 , (ρ cos θ, ρ sin θ)) andA(ρ, θ, k) = ρ a (ρ cos θ, ρ sin θ) + u 1 + u 2 2 , k (1 -χ(ρ cos θ, ρ sin θ)).

  (ρ, θ)) D N A(ρ, θ, k)dρ dθ where D is the differential operator acting on C ∞ ([δ, +∞[×[0, 2π]) given by Df = i ∂ ∂ρ ∂Ψ ∂ρ -1

+ u 2 2 2 j 2 + u 2 2 2 j 2 ≤ 1 +

 22221 +∞ δ exp (ikΨ(ρ, θ)) D N A(ρ, θ, k)dρ ≤ C N +∞ δ expk 2 ρ 2 w(ρ, θ) dρ with w(ρ, θ) = ρ j 1 1 + ρ(cos θ, sin θ) + u 1for some C N > 0 and j 1 , j 2 ∈ Z. If j 2 < 0, this last integral can be bounded by+∞ δ ρ j 1 expk 2 ρ 2 dρ = O(k -1/2 ); if j 2 > 0, we have 1 + ρ(cos θ, sin θ) + u 1 |ρ| 2 + |R| 2 j 2and hence+∞ δ expk 2 ρ 2 w(ρ, θ) dρ ≤ CN +∞ δ expk 2 ρ 2 ρ j 3 dρ = O(k -1/2 ).

  sequences of functions z ∈ C → ∂ã ∂ z1 (z, z, k) and z ∈ C → ∂ã ∂z 2 (z, z, k) are negligible. Hence, we have for ≥ 0 and z ∈ C ∂ã ∂ z1 (z, z) = 0 = ∂ã ∂z 2 (z, z).

2 C 2 |z 1 | 2 +

 2212 exp (ikφ(z 1 , z 2 , z 3 )) ã(z 1 , z 3 , k) b(z 3 , z 2 , k) dλ(z 3 ) +R k exp -Ck|z 1 -z 2 | 2 , with C > 0, R k negligible and φ(z 1 , z 2 , z 3 ) = i |z 2 | 2 + 2|z 3 | 2 -2z 1 z3 -2z 3 z2 .

  Lemma 3.5 yields that C k = D k + R k exp -C k|z 1 -z 2 | 2 for some C > 0 and R k negligible. It remains to compute the contravariant symbol of C k . For z in C, put ȃ(z, w, k) = ã(z, z + w, k) and b(z, w, k) = b(z + w, z, k). One has č = exp k -1 ∆ w ȃb |w=0with ∆ w the holomorphic laplacian with respect to w; using lemma 3.5, we find č(z, k)= exp k -1 ∂ ∂ ū ∂ ∂v ǎ(u, k)) b(v, k) |u=v=z up to a negligible term. Now, since c(., k) = exp -k -1 ∆ č (., k), this yields c(z, k) = exp -k -1 ∂ ∂u ∂ ∂v a(u, k)b(v, k)

Proposition 3.1. The

  tensor power, and λ the Lebesgue measure on R 2 . It is easily shown that for k ≥ 1, H 0 k is precisely B k . Sometimes, we will use the identification of the section f ψ to the function f in abusive notations, such as talking about the operator ∂ ∂z action on B k , etc. It is standard that B k is closed in L 2 R 2 , exp(-k|z| 2 )dλ(z) , and is thus a Hilbert space; moreover, we know an orthonormal basis of B k . family (ϕ n,k ) n∈N , where ϕ n,k

  [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform[END_REF] and Ω 2 if necessary, where we recall that the symbol ∼ stands for microlocal equality.

Proof. Start from a Fourier integral operator S k associated to s(., ., k) with principal symbol s 0 (v, χ -1 (v)) never vanishing on Ω 1 × Ω 2 . The first step is to construct an operator with the first property; we do it by induction, correcting S k by Toeplitz operators. More precisely, let P k be a Toeplitz operator on R 2 , and denote by p 0 its principal symbol. Set U (0)

Acknowledgements

I would like to thank San Vũ Ngo . c and Laurent Charles for their precious help and numerous readings of the manuscript.