Jérémy Le Deunff
email: jeremy.ledeunff@lmpt.univ-tours.fr

Amaury Mouchet
email: mouchet@lmpt.phys.univ-tours.fr

Instantons revisited: dynamical tunnelling and resonant tunnelling

Keywords: numbers: 05.45.Mt, 03.65.Sq, 03.65.Xp, 05.60.Gg, Tunnelling, dynamical tunnelling, resonant tunnelling, instantons, semiclassical methods, Wick rotation, complex time

Starting from trace formulae for the tunnelling splittings (or decay rates) analytically continued in the complex time domain, we obtain explicit semiclassical expansions in terms of complex trajectories that are selected with appropriate complex-time paths. We show how this instanton-like approach, which takes advantage of an incomplete Wick rotation, accurately reproduces tunnelling effects not only in the usual double-well potential but also in situations where a pure Wick rotation is insufficient, for instance dynamical tunnelling or resonant tunnelling. Even though only onedimensional autonomous Hamiltonian systems are quantitatively studied, we discuss the relevance of our method for multidimensional and/or chaotic tunnelling.

I. INTRODUCTION: EVENTS OCCUR IN COMPLEX TIME

Instantons generally refer to solutions of classical equations in the Euclidean space-time, i.e. once a Wick rotation t → -it has been performed on time t in the Minkowskian space-time. Since the mid-seventies, they have been extensively used in gauge field theories to describe tunnelling between degenerate vacua [START_REF]Advanced Series in Mathematical Physics[END_REF]. In introductory texts [2, 3, for instance], they are first presented within the framework of quantum mechanics: when the classical Hamiltonian of a system with one degree of freedom has the usual form

H(p, q) = 1 2 p 2 + V (q) (1)
(p and q denote the canonically conjugate variables), the Wick rotation induces an inversion of the potential and, then, some classical real solutions driven by the transformed Hamiltonian p 2 /2 -V (q) can be exploited to quantitatively describe a tunnelling transition. As far as we know, in this context, only the simplest situations have been considered, namely the tunnelling decay from an isolated minimum of V to a continuum and the tunnelling oscillations between N degenerated minima of V that are related by an N -fold symmetry. In those cases, what can be captured is tunnelling at the lowest energy only. However, not to speak of the highly non-trivial cases of tunnelling in non-autonomous and/or non-separable multidimensional systems, there are many situations that cannot be straightforwardly treated with a simple inversion of the 1d, time-independent, potential. First, tunnelling -i.e. any quantum phenomenon that cannot be described by real classical solutions of the original (non Wick-rotated) Hamilton's equationsmay manifest itself through a transition that is not necessarily a classically forbidden jump in position [4]. For instance, the reflection above an energy barrier, as a forbidden jump in momentum, is indeed a tunnelling process [5,6]. In the following, we will consider the case of a simple pendulum whose dynamics is governed by the potential V (q) = γ cos q; for energies larger than the strength γ > 0 of the potential, one can observe a quantum transition between states rotating in opposite direction while two distinct rotational classical solutions, obtained one from each other by the reflection symmetry, are always disconnected in real phase-space.

A second example is provided by a typical situation of resonant tunnelling: when, for instance, V has a third, deeper, well which lies in between two symmetric wells (see Fig. 6a), the oscillation frequency between the latter can be affected by several orders of magnitude, when two eigenenergies get nearly degenerate with a third one, corresponding to a state localised in the central well. Then, we lose the customary exponential weakness of tunnelling and it is worth stressing, coming back for a second to quantum field theory, that having a nearly full tunnelling transmission through a double barrier may have drastic consequences in some cosmological models. In this example, one can immediately see (Fig. 6b) that -V also has an energy barrier and working with a complete Wick rotation only remains insufficient in that case.

In order to describe the tunnelling transmission at an energy E below the top of an energy barrier, which may be crucial in some chemical reactions, the pioneer works by Freed [7], George and Miller [8,9], have shown that the computation of a Green's function G(q f , q i ; E) (or a scattering matrix element) requires taking into consideration classical trajectories with a complex time. These complex times come out when looking for the saddlepoint main contributions to the Fourier transform of the time propagator, G(q f , q i

; E) = 1 √ 2π ∞ 0 G(q f , q i ; t) e iEt/ dt , (2)
which is, up to now, the common step shared by all the approaches involving complex time [10-14, for instance]. Though well-suited for the study of scattering, indirect computations are required to extract from the poles of the energy Green's function (2) some spectral signatures of tunnelling in bounded systems. Generically, these signatures appear as small splittings between two quasi-degenerate energy levels and can be seen as a narrow avoided crossing of the two levels when a classical parameter is varied [15][16][17][18]. In the present article, we propose a unified treatment that provides a direct computation of these splittings (formulae 9 or 10); as shown in section II, it takes full advantage of the possibility of working not necessarily with purely imaginary time, but with a general parametrisation of complex time as first suggested in [19]. The semiclassical approach naturally follows (sec. III) and some general asymptotic expansions can be written (equations (C6) and ( 40)) and simplified (equation ( 73)); they constitute the main results of this paper. To understand where these formulae come from and how they work, we will start with the paradigmatic case of the double-well potential (sec. IV) and the simple pendulum (sec.V); we defer some general and technical justifications in the appendices. Then we will treat the resonant case in detail in sec. VI, where an appropriate incomplete Wick rotation t → e -iθ t provides the key to showing how interference effects à la Fabry-Pérot between several complex trajectories reproduce the nonexponential behaviour of resonant tunnelling, already at work in open systems with a double barrier [20,[START_REF] Bohm | Quantum theory[END_REF]. After having shown how to adapt our method to the computations of escape rates from a stable island in phase-space (sec. VII), we will conclude with more long-term considerations by explaining how our approach provides a natural and new starting point for studying tunnelling in multidimensional systems.

II. TUNNELLING SPLITTINGS

A particularly simple signature of tunnelling can be identified when the Hamiltonian has a two fold symmetry and, therefore, we will consider quantum systems whose time-independent Hamiltonian Ĥ = H(p, q) commute with an operator Ŝ such that Ŝ2 = 1 (the ˆallows to distinguish the quantum operators from the classical phase-space functions or maps). In most cases, Ŝ stands for the parity operator: H(-p, -q) = H(p, q) .

(

) 3
The spectrum of Ĥ can be classified according to Ŝ and, for simplicity, we will always consider a bounded system whose discrete energy spectrum and the associated orthonormal eigenbasis are defined by

Ĥ |φ ± n = E ± n |φ ± n ; Ŝ |φ ± n = ± |φ ± n ( 4
)
where n is a natural integer. When the Planck constant is small compared to the typical classical actions, standard semi-classical analysis [START_REF] Takahashi | [END_REF][START_REF] Heller | Chaos et Physique Quantique-Chaos and Quantum Physics, Les Houches, école d'été de physique théorique 1989, session LII[END_REF][START_REF] Torres-Vega | [END_REF] shows that one can associate some classical regions in phase-space to each eigenstate |φ ± n . This can be done by constructing a phase-space representation of |φ ± n , typically the Wigner or the Husimi representation, and look where the corresponding phase-space function φ ± n (p, q) is mainly localised (all the more sharply than is small). For a Hamiltonian of the form (1) where V has local minima, some of the eigenstates remain localised in the neighbourhood of the stable equilibrium points. For instance, for a double-well potential whose shape is shown in Fig. 1a), the symmetric state |φ + n and the antisymmetric state |φ - n at energies below the local maximum of the barrier have their Husimi representation localised around both stable equilibrium points (p, q) = (0, ±a), more precisely along the lines

H(p, q) = E (1d tori) at energy E ≃ E + n ≃ E - n .
Because of the non-exact degeneracy of E + n and E - n , any linear combination of |φ + n and |φ - n constructed in order to be localised in one well only will not be stationary anymore and will oscillate back and forth between -a and a at a tunnelling frequency |E - n -E + n |/(2π ). At low energy this process is classically forbidden by the energy barrier. This is very general, even for multidimensional non-integrable systems: the splitting ∆E n = |E - n -E + n | between some nearly degenerated doublets provides a quantitative manifestation of the tunnelling between the phase-space regions where the corresponding eigenstates are localised.

Rather than computing the poles of E → G(a, -a, E), another systematic strategy to obtain one individual splitting [START_REF]Statistical approaches have also been proposed[END_REF] is to start with Herring's formula [25][START_REF] Creagh | Tunneling in complex systems[END_REF][START_REF] Garg | [END_REF] that relies on the knowledge of the eigenfunctions outside the classically allowed regions in phase-space. Here, we will propose alternative formulae [28] that involve traces of a product of operators, among them the evolution operator,

Û (T ) def = e -i ĤT / = ∞ n=0 |φ + n φ + n | e -iE + n T / + |φ - n φ - n | e -iE - n T / , (5)
analytically continued in some sector of the complex time domain. This approach, which privileges the timedomain, provides, of course, a natural starting point for a semiclassical analysis in terms of classical complex orbits. The simplest situation occurs when the tunnelling doublet is made of the two lowest energies E + 0 E - 0 = E + 0 + ∆E 0 of the spectrum. If ω denotes the energy difference between the nearest excited states and E ± 0 (for the double well potential E ±

1 ≃ E ± 0 + ω where ω is the classical frequency around the stable equilibrium points), by giving a sufficiently large imaginary part to -T ,

-ω Im(T ) ≫ 1 , (6)
we can safely retain the n = 0 terms only, which exponentially dominate the trace of (5). To be valid, this approximation requires that we remain away from a quantum resonance where the definition of the tunnelling doublet is made ambiguous by the presence of a third energy level in the neighborhood of E + 0 and E - 0 . Then we have immediately

tan T ∆E 0 2 ≃ -i tr Ŝ Û (T ) tr Û (T ) . (7)
When the tunnelling splitting is smaller than ω by several orders of magnitude, we can work with a complex time such that

|T |∆E 0 2 ≪ 1 (8)
remains compatible with condition (6) and, therefore,

∆E 0 ≃ ∆ 0 (T ) def = 2 iT tr Ŝ Û (T ) tr Û (T ) (9)
will provide a good approximation of the tunnelling splitting for a wide range of T . Though condition ( 8) is widely fulfilled in many situations this is not an essential condition since one can keep working with tan -1 . Numerically, the estimation (9) has also the advantage of obviating a diagonalization.

When we want to compute a splitting, due to tunnelling, between an arbitrary doublet, the selection of the corresponding terms in the right hand side of (5) can be made with an operator Πn that will mimic the (a priori unknown) projector |φ

+ n φ + n | + |φ - n φ - n |.
It will be chosen such that its matrix elements are localised in the regions of phase-space where φ ± n (p, q) are dominant. Under the soft condition |T |∆E n /(2 ) ≪ 1, we will therefore take

∆E n ≃ ∆ n (T ) def = 2 iT tr Ŝ Πn Û (T ) tr Πn Û (T ) . ( 10
)
The localization condition on the matrix elements of Πn is a selection tool that replaces (6). In that case, there is a battle of exponentials between the exponentially small matrix elements φ ± m | Πn |φ ± m and the time dependent terms e -i(Em-En)T / . The last term would eventually dominate for the lower energy states (E m < E n ) if -Im T could be increased arbitrarily. But once a Im T < 0 is given, one expects to recover a good approximation of the excited splitting by increasing Re T since, when Re T ≫ | Im T |, we actually recover the real time case.

The next step consists in computing ∆ n (T ) by semiclassical techniques and, then, we will add some more specific prescriptions on the choice of T , in order to improve the accuracy of ∆E n . We will illustrate how this works in the examples of secs. IV, V and VI.

Let us mention another way to select excited states, that we did not exploit further. With the help of a positive smooth function F (u) that has a deep, isolated minimum at u = 0, say F (u) = u 2N with N strictly positive, we can freeze the dynamics around any energy E by considering a new Hamiltonian H ′ (p, q) def =F H(p, q) -E) . Classically, the phase-space portrait is the same as the original one obtained with Hamiltonian H except that the set of points H(p, q) = E now consists of equilibrium points. The quantum Hamiltonian Ĥ′ def =F ( Ĥ -E) has the same eigenfunctions as Ĥ but the corresponding spectrum is now

F (E ± n -E). By choosing E = E + n , the doublet E ±
n yields to the ground-state doublet ∆E ′ n = F (∆E n ) and we can use the approximation (9) with Û ′ = exp(-i Ĥ′ T ). With the F given above, we have

∆E n ≃ [∆ ′ n (T )] 1/(2N ) def = 2 iT tr Ŝ Û ′ (T ) tr Û ′ (T ) 1/(2N ) . (11)

III. SEMICLASSICAL EXPRESSIONS A. Hamiltonian dynamics with complex time

Formally, the numerator and the denominator of the right hand side of (9) can be written as a phase-space path integral of the form

P e iS[p,q;t]/ D[p]D[q] ( 12
)
where the continuous action is the functional

S[p, q; t] def = s f si p(s) dq ds (s) -H p(s), q(s) dt ds (s) ds.
(13) The subset P of phase-space paths, the measure D[p]D[q] and the action (13) appear as a continuous limit (τ → 0) of a discretized expression whose precise definition depends on the choice of the basis for computing the traces but, in any cases, involves a typical, finite, complex time step τ (see also appendix A). The complex continuous time path s → t(s) is given with fixed ends t(s i ) = t i = 0, t(s f ) = t f = T . Because the slicing of T in small complex time steps of modulus of order τ is arbitrary, the integrals of the form (12) remain independent of the choice of t(s) for s i < s < s f as long as Im t(s) is non-increasing in order to keep the evolution operators well-defined for any slice of time [19]. In the following we will denote by [t] such an admissible time-path. In a semiclassical limit (keeping the order lim →0 lim τ →0 [49]), the dominant contributions to integrals (12)

with appropriate boundary conditions imposed on some canonical variables at s i and/or s f . When H is an analytic function of the phase-space coordinates (p, q), we can take the real and imaginary parts of equations ( 14), use the Cauchy-Riemann equations that render explicit the entanglement between the real part and the imaginary part of any analytic function f (z) : Re(df

/dz) = ∂(Re f )/∂(Re z) = ∂(Im f )/∂(Im z) and Im(df /dz) = ∂(Im f )/∂(Re z) = -∂(Re f )/∂(Im
Therefore, the dynamics described in terms of complex canonical variables is equivalent to a dynamics that remains Hamiltonian -though not autonomous with respect to the parametrisation s whenever dt/ds varies with s -, involving twice as many degrees of freedom as the original system, namely (Re q, Im q) and their respectively conjugated momenta (Re p, -Im p). The new Hamiltonian function is then Re(Hdt/ds) but it would have been equivalent (though not canonically equivalent) to choose the other constant of motion Im(Hdt/ds) as a Hamiltonian.

What will be important in what follows is that [t] will not be given a priori. Unlike in the standard instanton approach where t is forced to remain on the imaginary axis, we will see that for describing tunnelling it is a more efficient strategy to look for some complex paths p(s), q(s) that naturally connect two phasespace regions and then deduce [t] from one of the equations [START_REF] Creagh | [END_REF]. It happens that in the two usual textbook examples that we mentioned in the first paragraph of our Introduction, the complex time [t] has a vanishing real part, but in more general cases where tunnelling between excited states is studied, this is no longer true.

B. Trace formulae

Let us privilege the q-representation and consider the analytic continuation for complex time T of the well-known Van Vleck approximation for the propagator:

G(q f , q i ; T ) ∼ →0 o (-1) κo det i 2π ∂ 2 S o ∂q i ∂q f e iSo/ .
(16) The sum involves (complex) classical trajectories o in phase-space i.e. solutions of equations (15) with q(s i ) = q i , q(s f ) = q f for a given [t] such that t(s i ) = 0, t(s f ) = T . The action S o is computed along o with definition [START_REF] Maitra | Classical, semiclassical and quantum dynamics in atoms[END_REF] and is considered as a function of (q f , q i ; T ). The integer κ o encapsulates the choice of the Riemann sheet where the square root is computed; it keeps a record of the number of points on o where the semiclassical approximation ( 16) fails. As far as we do not cross a bifurcation of classical trajectories when smoothly deforming [t], the number of o's, the value of S o and κ o do not depend of the choice of [t]. The numerator (resp. denominator) of ∆ 0 (T ) are given by the integral G(ηq, q; T )dq where η = -1 (resp. η = +1). Within the semiclassical approximation, when o contributes to the propagator G(ηq, q; T ), we have to evaluate (-1) κo dq det i 2π

∂ 2 S o ∂q i ∂q f (ηq,q;T )
e iSo(ηq,q,T )/ .

(17) The steepest descent method requires the determination of the critical points of S o (ηq, q, T ). Since the momenta at the end points of o are given by

p i = p(s i ) = -∂ qi S o (q f , q i ; T ) ; (18a) p f = p(s f ) = ∂ q f S o (q f , q i ; T ) , (18b)
the dominant contributions to tr Û (T ) come from periodic orbits, i.e. when (p f , q f ) = (p i , q i ), whereas the dominant contributions to tr Ŝ Û (T ) come from half symmetric periodic orbits, i.e. when (p f , q f ) = (-p i , -q i ). As explained in detail in the appendix A, one must distinguish the contributions of the zero length orbits e (the equilibrium points) from the non-zero length periodic orbits o. For one degree of freedom, their respective contributions are given, up to sign, by e -iH(pe,qe)T / e λeT /2 -ηe -λeT /2 , (

with ±λ e the two Lyapunov exponents of the equilibrium point e = (p e , q e ), and

(-1) µo β T β √ -2ηiπ dE o dT e iSo/ , (19b)
where o is a periodic orbit of period T o = T if η = +1 and half a symmetric periodic orbit of half period

T o = T if η = -1.
The sum runs over all the branches β that compose the geometrical set of points belonging to o. T β is the characteristic time (A12) on the branch β. The energy E o is implicitly defined by (A10) and µ o essentially counts the number of turning points on o.

Expressions (19) are purely geometric; their classical ingredients do not rely on a specific choice of canonical coordinates and they are independent of the choice of the basis to evaluate the traces.

In the general case, the most difficult part consists in determining which periodic orbits contribute to the semiclassical approximation of the traces. Condition T o = T is necessary but far from being sufficient; the structure of the complex paths (keeping a real time) may appear to be very subtle [29] and have been the subject of many recent delicate works [30]. Even for a simple oscillating integral, the determination of the complex critical points of the phase that do contribute is a highly non-trivial problem because it requires a global analysis: one must know how to deform the whole initial contour of integration to reach a steepest descent paths.

Our strategy consists in retaining the terms (19b) for which we can choose the complex time [t] to constrain the (half) periodic orbit o to keep one of the canonical coordinate (say q) real. This o surely contributes because we do not have to deform the q part of the integration domain of ( 12) in the complex plane. But we will see in the next section that two different periodic orbits with real q correspond to two different choices of [t]: For a chosen [t], only some isolated points, if any, in phase space, will provide a starting point (p i , q i ) of a trajectory with q(s) real all along o. When sliding slightly the initial real coordinate q i in the integral (17), it requires to change [t] as well to maintain q(s) real on the whole o. This is not a problem since all the quantities involved in expression ( 17) are [t]-independent if the deformation of [t] is small enough not to provoke a bifurcation of o, that is, whenever the initial point does dot cross one of the turning points, which are the boundaries of the branches β. The integral (17) may be computed with a fixed timepath or an adaptive one for each separated branch. This computation, presented in the appendix A, generates the sum over all branches that appears in formula (19b). To sum up, to keep a contribution to the traces, we must know if we can choose a shape of [t] in order to pick a real-q periodic orbit with a specific q i . This construction is certainly not unique (one may choose other constraints) so we will use the intuitive principle that the periodic orbits we choose will connect the two regions of phase-space that are concerned by tunnelling i.e. where φ ± n are dominant; this is justified by the presence of the operator Πn in the right hand side of (10). The prefactor

Π n (q, q ′ ) def = q| Πn |q ′ will select in the integral dq dq ′ Π n (q, q ′ ) G(ηq ′ , q; T ) , (20)
a domain around the projection onto the q-space of the appropriate tori. In the specific case of the ground-state splitting, condition (6) does the job of Π0 : the large value of Im T requires that the orbit approach at least one equilibrium point and then follow a separatrix line.

In general we will not be able to prove that other complex periodic orbits give sub-dominant contributions but the examples given in the following sections are rather convincing. Moreover, our criterion of selection allows us to justify the four rules presented in [31, sec. IIA] for computing the contributions of complex orbits to the semiclassical expansions of the energy Green's function.

IV. APPLICATION TO THE DOUBLE-WELL POTENTIAL

Let us show first that our strategy leads to the usual instanton results for a Hamiltonian of the form (1) with an even V having two stable symmetric equilibrium points at q = ±a (Fig. 1a); in their neighbourhood, the frequency of the small vibrations is ω. We are interested in the ground-state splitting ∆E 0 for small enough in order to have ∆E 0 ≪ ω.

Before we make any semiclassical approximations, we check in Fig. 2 the validity of the estimation ( 9) on the quartic potential by a) verifying that ∆ 0 (T ) is almost real and independent of the choice of the complex T provided that conditions ( 6) and ( 8) are fulfilled, and b) by checking that this constant gives a good approximation of the "exact" ∆E 0 computed by direct numerical diagonalization of the Hamiltonian [32].

As explained at the end of section sec. III, in phasespace we will try to find some time-path [t] that allows the existence a (half) symmetric periodic orbit o that connects two tori at (real) energy E 0 in the neighbourhood of (p, q) = (0, ±a) while q remains real. If we impose t(s i ) = 0 and t(s f ) = T , then E = E o (q f , q i , T ) is implicitly given by relation (A10). We will denote by q r (E) (resp. q ′ r (E)) the position of the turning point at energy E > 0 that lies in between q = 0 and q = a (resp. that is larger that a). The two branches p ± (q, E) = ± 2 E -V (q) are either purely real when V (q) E or purely imaginary when V (q) E. Then, from equation (14b),

dt ds = 1 p dq ds (21)
is purely real in the classically allowed region, while purely imaginary in the forbidden region. Therefore the complex time path [t] must have the shape of a descending staircase whose steps are made of pure real or pure imaginary variations of time (see Fig. 3). The complex orbit with real q can be represented in phase-space as a continuous concatenation of paths, following the lines E = (Re p) 2 /2 + V q) in the allowed region and -E = (Im p) 2 /2 -V (q) in the forbidden region. It is natural to represent o in the three dimensional section Im q = 0 of the complex phase-space with axes given by (Re q, Re p, Im p): the junctions at the turning points lie necessarily on the (Re p = 0, Im p = 0) axis (see Figs. 1c andd). A periodic orbit is made of a succession of repetitions of 14) with a purely real q evolving from one well to the other. In phase-space, this trajectory appears to be a concatenation of two types of curves that join on the (Re q)-axis: 1) a trajectory that lies in the phase-space plane (Re q, Re p) with real variations of t at energy E with the potential V and 2) a trajectory that lies in the phase-space plane (Re q, Im p) with imaginary variations of t at energy -E with the potential -V shown in b). In c) a family of constant energy curves is shown (horizontal green (light gray) for the first type, vertical red (dark gray) for the second type). In d) for a given energy, we show how the three curves l, c, r glue together at the turning points (0, ±qr). S denotes the separatrix Im p = ± 2V (q).

(i) primitive real periodic orbits r with energy E in the right region, that is, such that q r (E) q(s) q ′ r (E) and Im p(s) = 0; (ii) primitive complex periodic orbits c in the central region with purely imaginary p and real q, that is, such that -q r (E) q(s) q r (E) and Re p(s) = 0; (iii) primitive real periodic orbits l with energy E in the left region. They are obtained from the periodic orbits r by the symmetry S.

By denoting T r (E) and T c (E) the (real, positive) peri- 2)+ln( ω/π) in order to emphasize the contribution of the prefactor. The blue solid line corresponds to ln |∆0| -Λ( ) for T ≃ -4i and become indistinguishable from the previous one for 1 / 3. With the same substraction, our first proposal (39) gives the constant (ln √ π) ≃ 0.6 and our second one (44), which is the same as the formula in Landau and Lifshitz [33, sec. 50, problem 3], corresponds to 0 (dotted line). Garg's formula [27, eq.1.1] gives the constant (ln π/e ≃ 0.07) (dash-dotted line).

ods of the primitive periodic orbits r and c at energy E respectively, we have

T o (E) = w r T r (E) -iw c T c (E) -i 1 -η 4 T c (E) (22)
with the winding numbers w r and w c being non-negative integers. The trajectory o may contribute to the denominator (η = +1) or to the numerator (η = -1) of the right hand side of ( 9) provided that T o = T . For η = +1, o is periodic whereas, for η = -1, o is half a symmetric periodic orbit. We will keep the contributions of all orbits for which a staircase [t] can be constructed. We can understand from figure 3 that orbits differing by a small sliding of their initial q i can be obtained by a small modification of the length or height of the first step. All these contributions are summed up when performing the integral (17) on one branch β and correspond to one term in the sum β in formula (19b). When at least one of the winding number is strictly larger than one, several staircase time-paths can be constructed while keeping relation [START_REF] Takahashi | [END_REF]: they differ one from each other by a different partition into steps of the length and/or heights of the staircase. The corresponding orbits o can be obtained one from each other by a continuous smooth sliding of the steps of the staircase but during this process one cannot avoid an orbit starting at a turning point where a bifurcation occurs.

a) a') c') c) b') l q i -q i r q ′ i -q ′ i q ′ i -q ′ i q i -q i -q i q i q ′ i -q ′ i T c /2 τ T r -τ T r T c /2 τ T r -τ T r T c /2 T r τ ′ T r -τ ′ b) 0 0 0 T T T c Im p
Re q Re p Figure 3: (colour online) For the double well in Fig. 1, at a given energy 0 < E < Vmax, when T = 2Tr(E) -iTc(E)/2, for each one-step time path, there exists a unique half symmetric periodic orbit with energy E keeping Im q = 0 starting from the left allowed region. Its initial point is uniquely defined: For the time-path depicted in the left column, (pi, qi) is such that it takes exactly the real time τ to reach the turning point (0, -qr) following the primitive orbit l. Then it follows half of the orbit c and ends in the right well winding once along r (a). When starting at q ′ i < qi, there still exists a half periodic orbit joining (p ′ i , q ′ i ) to (-p ′ i , -q ′ i ) in time T for some p ′ i but unlike the previous one it gets outside the three dimensional section Im q = 0 (its projection is schematically represented by the dashed line in a'). However an appropriate shift of the time step to τ ′ > τ (central column) turns the latter smoothly into a trajectory with Im q = 0 (b'). During this process, the first trajectory starting at qi cannot keep a purely real q anymore (b). If we shift the step by exactly Tr (right column), we recover a trajectory starting at qi with purely real q (c) but with a different topology since it is now winding once along l. The change of topology occurs when the initial position reaches a turning point while sliding the step.

In the right hand side of ( 16), the sum involves several trajectories differing one from each other by the sequence of turning points that are successively encountered along o. Therefore, to compute the dominant contribution of the non-zero length orbits to the numerator and to the denominator of (9), we will add all the contributions of the topological classes of orbits, each of them uniquely characterised by an ordered sequence of turning points [ρ 1 , ρ 2 , . . . ], in other words by a partition of w r and w c into integers and by the branch β where its starting point lies. We can therefore express our result in a way that can be applied to cases more general than the double well: the total contribution of the non-zero length paths to the numerator (η = -1) and to the denominator (η = +1) of ( 9) is

s [ρ1,ρ2,... ] (-1) µo T β √ -2ηiπ dE o dT e iSo/ . (23)
s denotes a section of an energy surface in complex phase space corresponding to one purely real canonical variable.

[ρ 1 , ρ 2 , . . . ] is an ordered sequence of (not necessarily distinct) turning points that belong to a section s. In the case of the double-well, for an energy below V max , a section s for real q has four turning points (0, ±q r ) and (0, ±q ′ r ). Only the points on s such that -q ′ r < q < q ′ r can provide starting points of a periodic orbit. They belong to one of the three closed loops r,c, l that connect on the axis (p = 0, Im q = 0) at the turning points (0, ±q r ). Once T is given, the condition T o = T will select a finite set of energies for which we have [START_REF] Takahashi | [END_REF] with integers w r and w c . Any o with winding numbers w r and w c will have an action given by

S o = w r S r + iw c S c + i 1 -η 4 S c (24)
and an index

µ o = w r + w c + 1 -η 2 (25)
(in the real case, 2µ o is computed in the same way as the Maslov index: it counts the number of turning points encountered along o). These quantities are independent on the choice of the six possible starting branch β (each r, c, l is made of two branches). When the orbit o starts on r or on l, we will have T β = T r /2 and when o starts on c, we have

T β = -iT c /2.
For -Im T larger than the oscillation period in the central well of -V , w c = 1 is the minimum value of the winding number when η = +1 (just one back and forth trip around c) while it is w c = 0 when η = -1 (just one half of c is concerned). For these orbits, condition (6) forces c to stay near the separatrix S defined by Im p = ± 2V (q) and, thus, r must lie in the immediate neighbourhood of the equilibrium point (p, q) = (0, a). These orbits will give the dominant contribution because they have the smallest Im S o among all the other possible orbits involving repetitions of c. Indeed, for a fixed T , all the orbits o that may contribute semiclassically are such that Im S o = -2 Im(T )S c /T c and S c /T c is a decreasing function of the energy

Ẽ def = -E when c is inside the separatrix since d d Ẽ S c ( Ẽ) T c ( Ẽ) = d d Ẽ 4 T c ( Ẽ) qr (-Ẽ) 0 2 Ẽ + V (q) dq -Ẽ = - 4 T c ( Ẽ) 2 dT c ( Ẽ) d Ẽ qr (E) 0 |p(q, E)|dq < 0 . ( 26
)
Therefore Im S o reaches its minimum when Ẽ is at its maximum, that is E → 0 + . The only possible equilibrium point contributing to tr Ŝ Û (T ) is the origin (p f , q f ) = (-p i , q i ) = (0, 0); it is also subdominant because H(0, 0) = V (0) = V max (which can be seen as the limit of

S c ( Ẽ)/T c ( Ẽ) when Ẽ → -V max ) is strictly larger than S c ( Ẽ)/T c ( Ẽ) for Ẽ > -V max .
Assuming that only orbits with real q do contribute to tr Ŝ Û (T ) , we have proven that the dominant contribution is given by the half symmetric orbits o at energy E such that T c (E) = -2 Im T (w c = 0). However, for such an orbit to exist we cannot choose Re T arbitrarily since it must be an integer multiple of T r (E). To put it in another way, for a given E, we will choose T = w r T r (E) -iT c (E)/2 such that an orbit o with a real q exists. Condition (6) will be fulfilled if E is sufficiently small. We must now enumerate all the topological classes concerned by the sum ( 23): When o starts on the upper branch Im p > 0 of c, it will reach the turning point ρ 1 = (0, q r ) then wind w r times around r alternatively crossing (0, q ′ r ) and (0, q r ) before turning back on the lower branch Im p < 0 of c. Starting on the lower branch of c corresponds to the symmetric trajectory and provide the same contribution with T β = -iT c /2. When starting on the upper branch of r, o crosses ρ 1 = (0, q ′ r ) first, then reaches ρ 2 = (0, q r ). Then it may go on winding r times around r then take the lower branch of c up to the turning point (0, -q r ), wind w r -1 -r times around l and eventually join the symmetric of its starting point on the lower branch of l. There are exactly w r such topological classes because we can take r = 0, ..., w r -1. If we start from the lower branch of r or on one of the two branches of l, we obtain the same contribution and exhaust the possible topological classes. The sum of T β on all β's and classes is then 2(-iT c /2) + 4w r T r /2 = 2T and keeping only the w c = 0 solutions, the sum [START_REF] Heller | Chaos et Physique Quantique-Chaos and Quantum Physics, Les Houches, école d'été de physique théorique 1989, session LII[END_REF] reduces to

tr Ŝ Û (T ) ∼ →0 2T (-1) wr+1 √ 2iπ dE o dT e iSo/ (27)
with o being one half symmetric orbit of energy E o defined implicitly by

T = w r T r (E) -iT c (E)/2 (28)
with

T r (E) = 2 q ′ r (E) qr (E) dq 2 E -V (q) , (29a)
T c (E) = 4 qr(E) 0 dq 2 V (q) -E . ( 29b
)
We have

S o (E) = w r S r (E) + iS c (E)/2 with S r (E) + ET r (E) = Sr (E) def = 2 q ′ r (E) qr (E) 2 E -V (q) dq, (30a)
S c (E) -ET c (E) = Sc (E) def = 4 qr (E) 0 2 V (q) -E dq. ( 30b
)
The dominant contributions to tr Û (T ) comes from the two stable equilibrium points e = (0, ±a) for which λ e = iω. The contribution of (0, 0) is sub-dominant as well as the contribution of any periodic orbit o which necessarily turns around about e during Re T then follow an orbit c near S during -Im T before coming back to its initial point. Together with ( 6), η = 1, the two stable equilibrium points give two identical contributions (19a) and we have

tr Û (T ) ∼ →0 2 e -iωT /2 (31)
which of course could have been deduced directly from E ± 0 ≃ ω/2. Collecting all these results in the right hand side of (9), we obtain

∆ 0 ∼ →0 2i π dT o dE -1/2 × exp - 1 2 Sc (E) + (E -ω/2)T c (E) × exp iw r Sr (E) -(E -ω/2)T r (E) -π . ( 32
)
When E → 0 + , we have the following asymptotic expansions (see appendix B)

Sc (E) = Sc (0) + 4E ω ln √ 2E 2aω - 2(2A + 1) ω E + o(E) ; (33a) Sr (E) = 2πE ω + BE 2 + o(E 2 ) ; (33b)
with

A def = a 0 ω 2V (q) - 1 a -q dq ; (34)
B def = π 24ω 7 5V (3) (a) 2 -3ω 2 V (4) (a) (35)
(the superscript in parenthesis indicates the order of the derivative of V ). The differentiation of expressions [START_REF] Landau | Quantum Mechanics (non relativistic theory[END_REF] with respect to E leads to the asymptotic expansions for -T c and T r . From the first one we can extract the exponential sensitivity of E on Im T :

E ≃ 2a 2 ω 2 e 2A e ω Im T . (36)
From relation (28) we can see that 2dT o /dE = 2w r dT r /dEi dT c /dE is dominated by the last term if E is small.

dT o dE ≃ - i 2 dT c dE ≃ i ωE . (37)
Inserting all these asymptotic expansions in the right hand side of (32), we get

∆ 0 ∼ →0 4 a 2 ω 3 π e A e E/( ω) e -Sc(0)/(2 )
× e -iwrBE(E-ω)/ . [START_REF] Khuat-Duy | [END_REF] This expression can be turned into the usual jwkb expansion exp a 0 (E)/ + a 1 (E) ln + a 2 (E) + o(1) : As soon as we have condition ( 6), from expression [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] we see that E is exponentially small and we obtain a 0 (0) = -Sc (0)/2 and a 1 (0) = 1/2. To obtain the correct value of a 2 , we must proceed to fine-tune the choice of T . A criterion is to impose on ∆ 0 to have a vanishing imaginary part at any order in consistent with the jwkb expansions used so far. From ( 32), we will choose E such that Sr (E) = π, which is exactly the Einstein-Brillouin-Keller quantization condition for the ground state in one well. This leads to E = ω/2 + o( ).

Then T c = -2 ln /(4a 2 ω) /ω + 4A/ω + o(1) and

∆ 0 ∼ →0 ω √ π e -Sc( ω/2)/(2 ) = 2aω
e ω π e A e -Sc(0)/(2 ) , (39) which differs from [27, eqs. (1.1) and (1.4)] by a reasonable factor √ e ≃ 1.6. This discrepancy, already noticed in [34, sec. V], which appears in the third order term in the -expansion, comes from the different kind of approximations involved in our approach on the one hand and in Herring's formula on the other hand.

We are also able to obtain a formula for the splitting of the excited states that is consistent with the result given in [27, eq. (B1)]. Using a semiclassical approximation for the matrix element of Πn , we explain in detail in appendix C how to obtain ∆ n (T ). For one dimensional systems whose energy surface E n is made of two branches (two Riemann sheets in the complex plane), we can insert (C8) into (C6) and get one of the main result of this paper:

∆ n (T ) ∼ →0 2T [ρ1,ρ2,... ] (-1) µo+1 e i So(bβ′ ,b β ,En)/ . ( 40
)
To see how formula (40) works in the case of the doublewell potential, we choose the quasi-mode |Φ n localised on the right torus r at energy E n . This torus is made of two branches labelled by the sign of p and we can choose a common base point for these two branches, namely b ± = q r (E). On the symmetric torus l, the two base points will be b ± = -q r (E). Then So (b β ′ , b β , E n ) and the index µ o do not depend on the choice of the initial and final branch. The orbits o that go from r to l must correspond to a T of the form

T = τ + w r T r (E n ) -i w c + 1 2 T c (E n ) (41)
for non-negative integer (w c , w r ) and a fraction of time τ strictly smaller that T c (E n ) that depends on the initial and final conditions (those are not necessarily symmetric). Then we have

So (b β ′ , b β , E n ) = w r Sr (E n ) + i w c + 1 2 Sc (E n ) (42)
and we take

µ o = w r + w c + 1 . (43)
For the same reason as previously explained the dominant contributions come from those orbits where w c = 0.

In order to mimic a real T , we will choose large winding numbers w r such that Re

T = w r T r (E n ) ≫ Im T = T c (E n )/2.
Because of the quantization condition in the right well Sr (E n ) = (n + 1/2)2π , the rapid oscillations exp(i Re S o / ) disappear (or inversely if we want to maintain ∆ n real to first order, we recover the usual quantization condition). Then we obtain

∆ n ∼ →0 2 T r (E n ) e -Sc(En)/(2 ) . ( 44
)
The classical frequency 1/T r (E n ) attached to r is essentially of order ω/(2π). ∆ n (T ) becomes independent of T for large w r : the behaviour of T is mainly governed by w r T r , then the 1/T prefactor in ( 40) is compensated by the increasing number of identical terms in the sum since we have seen that the number of topological classes of orbits increases linearly with 4w r (the factor 4 comes from the two possible initial branches β = ± and the two possible final branches β ′ = ±; in other words from the sequences of [ρ 1 , ρ 2 , . . . ] beginning either by q r or q ′ r and ending either by -q r or -q ′ r ). The discrepancy between (44) and Garg's formula is just the factor g n given by [27, eq. ( B2)] (see also [35, eq. (3.41)]) that tends to 1 when n increases: g 0 = π/e ≃ 1.075, g 1 ≃ 1.028, . . . There is also a ratio of order one, more precisely 2/π ≃ 0.8, between estimations (39) and ( 44) taken for n = 0; the second is slightly better and coincides with the formula given in Landau and Lifshitz [33, sec. 50, problem 3]. Again, these discrepancies come from the different nature of the approximations that are involved.

Let us end this section by a short comment on the connection with the usual instanton theory where Re T = 0 and where Im T → -∞. This regime, that allows to select the ground-state doublet only, is included in our approach because the instanton trajectories appear to be the limit of c getting closer to the separatrix S whereas the classical real oscillations in the wells shrink to the equilibrium points. All along this paper we emphasize that the phase-space representation is particularly appropriate and it is straightforward to recover the usual picture of instantons (for instance q(it) versus it) from our figure 1 d).

V. DYNAMICAL TUNNELLING FOR THE SIMPLE PENDULUM

The simple pendulum corresponds to V (q) = -γ cos q with γ > 0 and strictly periodic boundary conditions that identify q = -π and q = π. At energy E > γ, the classical rotation with p = 2(E -V (q)) can never switch to the inverse rotation with -p = -2(E -V (q)). At the quantum level, the Schrödinger's equation for the stationary wave-function φ

- 2 2 d 2 dq 2 -γ cos q φ(q) = Eφ(q) (45)
leads to the Mathieu equation [START_REF] Abramowitz | Handbook of mathematical functions[END_REF],

y ′′ (x) + a -2g cos(2x) y(x) = 0 (46)
with x = q/2, y(x) = φ(2x), a = 8E/ 2 and g = -4γ/ 2 . The 2π-periodicity of φ forces y to be πperiodic. The eigenfunctions can be classified according to the parity operator: Even π-periodic solutions exist only for a countably infinite set of characteristic values of a denoted by {a 2n } with n = 0, 1, . . . Odd solutions correspond to another set, {b 2n }, with n = 1, 2 . . . (the {a 2n+1 , b 2n+1 } correspond to π-antiperiodic solutions and will be rejected). The discrete energy spectrum corresponding to even and odd solutions of ( 45) is then {E + n = 2 a 2n /8, n = 0, 1 . . . } and {E - n = 2 b 2n /8, n = 1, 2 . . . } respectively. Any eigenstate |φ ± n with energy E ± n ≃ E n > γ has its Husimi distribution spread symmetrically between the two half phase-space of positive and negative p, near the lines ± 2(E n -V (q)) that define two disconnected tori in the cylindrical phasespace. If we prepare a wave-packet localised on the line p = 2(E n -V (q)), since it is no longer a stationary state, its average momentum will oscillate between two opposite values, with a tunnelling frequency equal to ∆E n / where

∆E n = E + n -E - n = 2 8 (a 2n -b 2n ) ( 47
)
is the splitting between the two quasi-degenerate eigenenergies. To compute ∆E n , we will use (10) with the operator Πn very much like the exact projector |φ

+ n φ + n | + |φ - n φ -
n |: Its matrix element p ′ | Πn |p will vanish rapidly as soon as p or p ′ lie outside the region of the two tori. The main contribution to the semiclassical expansion of tr Ŝ Πn Û (T ) will come from classical trajectories that connect two symmetric tori. The trace will be semiclassically computed in the momentum basis and we will choose the complex time path to maintain p real.

To construct one trajectory at energy E > γ connecting the two tori requires to have purely imaginary q whenever -p r < p < p r with p r def = 2(E -γ) being the classical turning point in momentum. More precisely, q is given by cos q = -1 -(p 2 r -p 2 )/(2γ). From (14a) dt ds = -1 γ sin q dp ds [START_REF]Statistical approaches have also been proposed[END_REF] we see immediately that sin q(s) and dt/ds will be real when p(s) > p r or p(s) < -p r and purely imaginary otherwise.

In the latter case, since we keep p 2 (s) def = Im p(s) = 0, equations (15b,c) lead to two possible families of solutions (i) q 1 (s) def = Re q(s) = 0 and (ii) q 1 (s) = ±π. Then, with p 1 def = Re p and q 2 def = -Im q, equations (15a,d) become

dp 1 ds = ±γ sh q 2 i dt ds ; (49a) dq 2 ds = p 1 i dt ds (49b)
and are associated with a real time dynamics governed by the Hamiltonian

H(p 1 , q 2 ) = p 2 1 /2 ∓ γ ch q 2 . ( 50
)
with "-" corresponding to case (i) and "+" corresponding to case (ii). Then, instantons correspond to trajectories evolving in the transformed potential Ṽ (q) def = ∓γ ch q rather than the usual inverted potential -V (q). We will choose a one step complex time path as in Fig. 3 and we will represent the orbits in a three dimensional space (Re q, Re p, Im q) and the connection between the allowed and the forbidden trajectories occurs on the axes (Re q = 0, Im q = 0) and (Re q = ±π, Im q = 0) (Fig. 4a). Trajectories in family (i) escape from the unstable point at Re p = 0, Im q = 0 without coming back to the plane Re p = 0, Re q = 0. Only orbits in family (ii) can be used to produce periodic orbits. A typical periodic orbit o connecting two symmetric rotations of the pendulum is given in Fig. 4b). We will choose Im T to be precisely the half period of the periodic orbit c of family (ii) at energy E n and Re T to be (an integer multiple of) the period of rotation of the pendulum at the same energy. This choice exhibits the dominant contribution to tr Ŝ Πn Û (T ) . We can reproduce the same reasoning that led to (44) with the rôle of p and q being exchanged. Now T r = d Sr /dE is the typical frequency on the torus at energy E n .

Keeping only the contributions that provide a positive All these expressions can be written in terms of the complete elliptic integrals [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF] (defined for |u| < 1)

K(u) def = π/2 0 dx 1 -u 2 sin 2 x ; (54a) E(u) def = π/2 0 1 -u 2 sin 2 x dx . ( 54b
)
Namely,

Sr (E) = 4 2(E + γ) E 2γ E + γ ; (55a) Sc (E) = 8 2(E + γ) K E -γ E + γ -E E -γ E + γ . ( 55b
)
Then we obtain

∆ n ∼ 2 T r (E n ) e -Sc(En)/(2 ) (56) with T r (E) = d Sr dE = 2 √ 2 √ E + γ K 2γ E + γ . ( 57
)
The energies of the highly excited states are approximately given by the free rotations:

E n ≃ n 2 2 /2 ≫ γ. Then T r (E n ) ≃ 2π/(n ).
The last expression corresponds exactly to equation (3.44) of [35] obtained with standard uniform semiclassical analysis. We see on figure 5 that (56) is a very good approximation even when the energies E n get close to γ the energy of the separatrix.

VI. RESONANT TUNNELLING AND FABRY-P ÉROT EFFECT

We are now ready the see how formula (40) allows us to reproduce the resonant tunnelling between two wells related by parity when the potential in (1) has also a 47)) versus 1/ for n = 1, . . . 9 (black lines). The dots are given by the semiclassical approximation ln(∆n) (see eq. ( 58b)) up to the maximum value of 1/ (vertical dotted lines) where E + n become lower that the separatrix energy γ. The semiclassical limit is obtained when the states become more and more localised in the region of phase space corresponding to rotations: for a fixed n this requires to increase En ≃ n 2 2 /2 that is to increase ; or for a fixed energy E > γ, we must increase n ∝ √ 2E/ . To control the validity of the prefactors, we can check in b) that ln(∆En/∆n) approaches zero in the semiclassical limit.

deeper central well (Fig. 6 a). The minimum of the right and left wells is fixed at zero, the minimum of the central well is V (0) = V min < 0 and the local maximum between the wells is denoted V max . When 0 < E < V max , we will denote by q r (E) < q ′ r (E) < q ′′ r (E) the three positive solutions of V (q) = E. As explained in section (IV), we will try to construct appropriate time paths [t] to exhibit complex trajectories o with purely real q that connect the two symmetric tori from r to l at some We still maintain the position q to be real and can use a three dimensional section of phase-space (see text). Here we have chosen V given by (65) with a = 7/4 and b = 1/2. energy 0 < E < V max . These tori are delimited by the two turning points ±ρ ′ = (0, ±q ′ r ) and ±ρ ′′ = (0, ±q ′′ r ). What is new of course is the existence of a central real torus m delimited by ρ = (0, q r ) and -ρ. Using the three dimensional. representation of the section Im q = 0 of phase space (Fig. 6 c ande), we see that the orbits o must be a series of concatenation of five trajectories connecting at the turning points ±ρ, ±ρ ′ . First we start with one portion living on r with dt/ds > 0 and Im p = 0, then connect at ρ ′ to a trajectory with Re p = 0 where idt/ds > 0. It follows the energy curve c whose equation is (Im p) 2 /2 -V (q) = -E. Then o can connect at ρ to a real trajectory on m with Im p = 0 and dt/ds > 0, then can cross -c from -ρ to -ρ ′ before reaching l = -r. The corresponding time path will necessary have at least two steps (Fig. 6 d) each of them having a height which is a half-integer multiple of T c , the real period of the primitive periodic orbit c. Among all the possible o's, we will keep only the exponentially dom-inant contributions, when o remains as shortly as possible with complex p. Then, for such orbits to exist, we must choose T of the form

T = τ + w r T r (E) + w m + 1 2 T m (E) -iT c (E) (59)
involving the periods of the primitive orbits and the corresponding winding numbers w r , w m which are nonnegative integers. τ denotes a positive fraction of time smaller than T r (E). The base points for the two branches β = ± defined by Re p ≷ 0 on r coincide with the turning point b ± = q ′ r , the action So (b β ′ , b β , E n ) and the index µ o are independent of the choice of the branch where o starts:

So (b β ′ , b β , E n ) = w r Sr (E) + w m + 1 2 Sm (E) + i Sc (E) . (60)
µ o = w r + w m + 3 . (61)
Explicitly, we have

Sr (E) def = 2 q ′′ r (E) q ′ r (E) 2 E -V (q) dq ; (62a) Sm (E) def = 4 qr (E) 0 2 E -V (q) dq ; (62b) Sc (E) def = 2 q ′ r (E) qr (E) 2 V (q) -E dq ; (62c)
and the corresponding periods T r , T m , T c are obtained by deriving with respect to E.

As we have seen, formula (40) will provide an approximation of the exact splitting ∆E n that becomes better, the better the condition Re T (E n ) ≫ Im T (E n ) is satisfied. Not only this condition render the precise value of τ irrelevant, it also requires large w m and/or w r , especially if we work with E n 0 for which T c (E n ) ≫ T m (E) and T c (E n ) ≫ T r (E). For a given pair of w r and w m , there are 4(w r + 1) topological classes of orbits corresponding to two possible initial branches, two possible final branches and r = 0, . . . , w r possible windings on r for w r -r windings on l. Then all different w r and w m , such that (59) holds, give a contribution to (40

): If we define R(T ) def = Re T -τ -T m /2 ≃ Re T then, ∆ n (T ) ∼ →0 2 T e -Sc(En)/
× {wr,wm} pos. int. such that wrTr(En)+wmTm(En)=R(T )

(w r +1) e iwr[ Sr(En)/ -π]+iwm[ Sm(En)/ -π] . (63)
We immediately see the resonance at work since the sum reaches a maximum when both ν r 30. The (blue) dots correspond to estimation (A13) with T = |T |e -iθ = (K + 1/2)Tm(En) -iTc(En) in the limit K → ∞ (θ → 0). The upper thin (red) curve provides ln |∆En/∆n|. The inset provides a magnification of the third spike around 1/ = 6.6 for which Tm(En) = 2Tr(En) ≃ 1.58 and Tc(En) ≃ 1.50. Here, we have also plotted ln |∆n(T )| for K = 100 (thin orange (light gray), θ ≃ 0.01) K = 30 (thin green θ ≃ 0.03), K = 5 (dashed thin magenta θ ≃ 0.17).

E ± n ≃ ω/2 with ω = 2a 2(a 2 -b 2 ) ≃ 6.
degenerate (up to 2 terms) with the energy doublet in the lateral wells. Then the contributions of the repetitions of m interfere constructively like the optical rays in a Fabry-Pérot interferometer[21, secs. 12.14-12.17]. To estimate the sum in the right hand side of (63), let us take a rational approximation of the ratio T m /T r , namely

T m T r ≃ m r (64)
with m and r being coprimes positive integers. For the polynomial potential

V (q) = (q 2 -a 2 ) 2 (q 2 -b 2 ) (65)
with a > b > 0, the argument presented in [START_REF] Khuat-Duy | [END_REF] can be generalised to show that (64) is actually exact with r = 1 and m = 2 for any energy 0 < E < V max . If K denotes the integer part of R(T )/(rT m ), we can compute and approximate for K ≫ 1 the right hand side of (63) and obtain

|∆ n (T )| ∼ →0 T r
e -Sc(En)/ sin π(mν r -rν m ) .

(66) Fig. 7 shows that this latter expression provides a good approximation for ∆E n even in the immediate neighbourhood of a resonance where estimation (10) is not justified any more. If we had continued working with a finite K, the sum (63) would have involved a finite number of terms and the singularities due to the vanishing denominators in (66) would have been smoothed down (inset in Fig. 7). In other words, for a fixed |T |, rotating down T in the lower half plane Im T < 0, destroys very quickly (θ ≃ 0.2) the large resonant fluctuations of tunnelling. This effect has already been shown in the case of a kicked system [28].

VII. ESCAPE RATES

So far we have focused our analysis of tunnelling in bounded systems only, but the philosophy we presented here can be extended to more general situations. For instance, let us show how we can compute the escape rate from a metastable state localised in a confining potential V whose shape has the form given in Fig. 8 a). The potential has a local minimum at q = 0 (say V min = 0) and an energy barrier for 0 < q < q max whose height is V max . For q > q max , the potential remains non-positive and therefore in real phase-space (Re p, Re q), V defines around the origin an island of stability made of tori with positive energy. One state whose Husimi distribution is initially localised in the island, say a quasi-mode |Φ n at energy 0 < E n < V max , will progressively decay outside the well. The decay rate Γ n is then defined from the overlap:

Φ n | Û (T ) |Φ n = e -ΓnT /2-iEnT / .
(67)

If we choose a complex T such that

|Γ n T | ≪ 1 , (68)
then we obtain a trace formula for Γ n with the help of the projector-like operator Πn

Γ n ≃ - 2 Im T Im e iEnT / tr Πn Û (T ) (69)
which allows an explicit semiclassical expansion in terms of classical solutions with a complex time path. The dominant contributions will be provided by periodic orbits in time T with real q starting on the torus at energy E n , then, while idt/ds > 0, going forth outside the well before coming back to its initial starting point (Fig. 8 c,d)). Then, up to a dimensionless factor f n of order one, we get

Γ n ≃ f n T c (E n )
e -Sc(En)/ (70) The escape rate at energy E from an island created by a confining potential V whose shape is shown in a) can be semiclassically computed in terms of complex classical orbits where one canonical coordinate, say q, is kept real. In phase-space, this trajectory appears to be a concatenation of two types of curves that join on the (Re q)-axis: 1) a trajectory that lies in the phase-space plane (Re q, Re p) with real variations of t at energy E with the potential V and 2) a trajectory that lies in the phase-space plane (Re q, Im p) with imaginary variations of t at energy -E with the potential -V shown in b). In c) a family of constant energy curves is shown (horizontal green for the first type, vertical red for the second type). In d) for a given energy, we show how three curves l, c, r glue together at the turning points.

with

Sc (E) def = 2 q ′ r (E) qr (E) 2 V (q) -E dq ( 71
)
where 0 < q r (E) < q ′ r (E) are the two right turning points at energy E. Here we have supposed that the action between the two left turning points is larger than (71) and gives a subdominant contribution. If the potential is symmetric, two symmetric complex orbits would contribute with the same weight and therefore the escape rate should twice as large. In the case of an island with sharp boundaries around q max , Sc (E) which is the area enclosed by the primitive orbit c (Fig. 8 d) is mainly given by the portion of c where we can keep an harmonic approximation for the potential: V (q) ∼ ω 2 q 2 /2, q r (E) ≃ √ 2E/ω, q ′ r (E) ≃ A/(πω) where A is the area of the island in the real phase space. Then, with a def =ωA/(2πE), we have

Sc (E) ≃ 2E ω √ a √ a -1 -ln √ a + √ a -1 . ( 72
)
Inserting this expression in (69) with E n ≃ (n + 1/2) ω, we exactly recover the expression (5) used in [39] with an elegant and simple interpretation. The chaotic sea that surrounds the integrable island in the mixed system considered by Bäcker et al. acts as a sharp effective potential barrier as the one draw in picture 8; our complex trajectory that allows to escape from the regular region has its main features governed by the integrable (and even harmonic) approximation of the dynamics about the island, following precisely the general philosophy of [39] and [40]. This computation of the "direct" tunnelling (by opposition to resonant tunnelling where the model of a pure quadratic kinetic energy fails) can also be reproduced within the standard one-dimensional JWKB theory used for computing transmission coefficients.

Here again, we can check easily that the traditionnal instanton method is included in our approach : The regime ℜT → 0 selects the instanton solution and provides the escape rate from the equilibrium point as given by equation (2.47) of [2, chap. 7].

VIII. CONCLUSIONS

The explicit semiclassical expansions of trace formulae for tunnelling splittings (or escape rates) in terms of classical orbits constructed with complex-time paths provide an interesting alternative approach to Herring formulae essentially because they do not require to analytically continue the wave functions in the complex plane. In multidimensional tunnelling and/or for a nonautonomous Hamiltonian system, the generic lack of constants of motion isolates the stable islands (if any) from each other by chaotic seas. The analytic continuation of the kam tori that build the islands is prevented by the existence of natural boundaries (see for instance the recent discussion in [30, part II, appendix B] and references therein). The approach we have presented here seems to circumvent these difficulties but, of course, the problem of how to select, with an appropriate [t], the relevant trajectories among an a priori exponentially growing number of classical complex solutions remains open. We expect the tunnelling splittings between two states at energy around E localised in two symmetric islands to be approximated by the expansion of the form

f T [ρ1,ρ2,... ] e i N i=1 wi[ Si(En)/ -π] e -M j=1 w ′ j [ S′ j (En)/ ] (73)
where the sum runs over all possible sequences of turning points [ρ 1 , ρ 2 , . . . ] at energy E such that we can choose a complex-time path that leads to a trajectory, made of primitive orbits, that connects in time T the two (real) tori. The winding {w i } and the actions {S i } (resp. {w ′ j } and { S′ j }) refer to the primitive orbits obtained when the variations dt/ds are purely real (resp. purely imaginary). For dimensions larger that one, the dimensionless prefactor f may appear as a power law in . Here, inspired by the study in section VI, we can qualitatively see how the constructive interferences between repeated paths emerge in a speckle-like forest because of the presence of resonances. As shown in [28], a progressive complex rotation of time provides a natural way to select the main resonance effects. If we want to expand the splittings (or the escape rates) according to elementary process as proposed in [40], our approach offers a promising tool to interpret and compute semiclassically all of the ingredients of such an expansion. up to a global sign. Using (A5) again, the contribution (A6) can be written as

(-1) νo (-η) d det(1 -ηM o ) e iSo/ (A7)
where (-1) νo fixes the sign of the square root and d demotes the number of degrees of freedom. Recall that in the definition of the path-integrals ( 12), a time-slice τ is always implicit. Let us keep for a moment an explicit discretization, for instance with T being an integer multiple of τ , o referring to a discrete set of points, Hamilton's equations ( 14) being discretized into a phase-space map, path integral (13) being turned into a discrete (Riemann) sum, etc. Then, the contribution of each of the points of a given o, such that (p f , q f ) = (ηp i , ηq i ), is the same and remains given by (A7); but this is correct only if, for a given τ , is small enough, because when performing the steepest-descent method on (17), one must be able to split the oscillating integral into separate contributions coming from two distinct points of o. In the continuous time limit, i.e. when the limit τ → 0 is taken before the semiclassical limit → 0 (see [??]), the orbits o of non-zero length appear as a one-dimensional continuum none of whose points can be considered separately anymore. The only isolated critical points q c are given by the position of some equilibrium points. Under the symmetric condition (q i , p i ) = (-p i , -q i ), only the origin must be examined. Linearising the Hamiltonian flow about a non-degenerated fixed point e = (p e , q e ) leads to a monodromy matrix whose eigenvalues can be collected by pairs exp(±λ e,α T ) where {λ e,α } α=1,...d are the Lyapunov exponents. Then (A7) becomes (-1) νe e -iH(pe,qe)T / d α=1 (e λe,αT /2 -ηe -λe,αT /2 ) (A8) with a possible adjustment of the sign. For a generic choice of T , the denominator does not vanish. For a non-zero length path o, the critical q's are degenerate along the trajectory; for a system with several degrees of freedom, one must treat separately the (gaussian) integrals on the transverse coordinates q ⊥ along which S o (ηq, q, T ) varies (quadratically) from the longitudinal coordinates q along which S o (ηq, q, T ) is constant. Of course, the dimensions of q ⊥ and q depend crucially on the presence of kam tori. However, multidimensional tunnelling is beyond the scope of this paper and, the quantitative studies presented here concern onedimensional systems only. The contribution to the trace of such a path is then, up to a global sign,

e iSo/ √ -2iπ dq ∂ 2 S o ∂q i ∂q f (ηq,q;T ) . ( A9
)
The conservation of energy along o, H(p, q) = E, implicitly defines a function p(q, E) in the neighbourhood of any point where ∂ p H = 0. Globally along the trajectory o, we may encounter several possible branches p β (q, E) for the graph of these functions (the two possible signs of a square root when the Hamiltonian has the form (1)) which become singular but pairwise connect smoothly at the turning points, defined by ∂ p H = 0. A relation between E, T , q f and q i can be obtained by integrating s f si (dt/ds) ds and using (14b):

T = T o (q f , q i ; E) def = s f si dq/ds ∂ p H p q(s), E , q(s)
ds .

(A10) Everywhere but at the turning points, the value of s dictates the choice of the branch used for the integrand. This relation implicitly defines E o (q f , q i ; T ). The usual expression for the derivative of implicit functions leads to the relations:

∂ E p = 1/∂ p H and ∂ qi E o = -∂ qi T o /∂ E T o = ∂ T E o /∂ p H(p i , q i ).
If we differentiate (18b) with respect to q i when p f is given by p q f , E o (q f , q i , T ) , we obtain

∂ 2 S o ∂q i ∂q f (q f ,qi;T ) = 1 ∂ p H(p f , q f ) 1 ∂ p H(p i , q i ) ∂E o ∂T (q f ,qi;T )
.

(A11) With property (3), we get ∂ p H(ηp, ηq) = η∂ p H(p, q) and then, by differentiating (A10), we have ∂ q T o (ηq, q, E) = 0. Therefore, the energy E o of the path o depends only on T , not on its starting point q. The square root of ∂ T E o (ηq, q, T ) = dE o /dT can be got out from the integral in (A9). When adding the contribution of each path whose starting point lie on the branch β, we obtain

T β (E) def = dq ∂ p H p β q, E , q . (A12)
This is not exactly the right hand side of (A10) because the domain of integration in (A12) is the domain of the branch where the starting point of o lives. Each branch β is delimited by two turning points and T β is the time spent to go from one point to the other. The integral (A9) involves all the possible starting points for a trajectory o and therefore we must add all the branches that patchwork smoothly in phase-space to form the geometrical set of points crossed by o. Referring to the purely geometrical quantities (i.e. independent of the choice of the parametrization), we have the contribution

(-1) µo β T β √ -2ηiπ dE o dT e iSo/ (A13)
only if T o = T ; for η = +1 the path o is a periodic orbit and, for η = -1, the path o half a symmetric periodic orbit (the whole periodic orbit being of period equal to 2T ). The sum concerns all the geometrical branches β crossed by o (even if o passes several times by the same points, each branch is only counted once). As before, the conversion of a product of two square roots of complex numbers to the square root of the product may introduce a sign that can be absorbed in the definition of µ o ; the exact computation of the index µ o is difficult but since it may change at the bifurcation points only, where the semiclassical approximation fails, it is sufficient to know that it depends on the nature and the number of the turning points encountered on o. Therefore it is an additive quantity when several primitive orbits are repeated or concatenated together.

Appendix B

In this appendix we explain how to obtain the asymptotic expansions [START_REF] Landau | Quantum Mechanics (non relativistic theory[END_REF] as E → 0 + .

First consider Sr and split it in two parts S+ + Swhere

S+ def = q ′ r (E) a 2 2 E -V (q) dq ; (B1a) S- def = a qr (E)
2 2 E -V (q) dq .

Setting ǫ def = a -q r (E), rewrite Sas

S-= 2 √ 2ǫ 1 0 V (a -ǫ) -V (a -sǫ) ds, (B2)
Now expand the integrand as a power series in ǫ up to the fourth order, compute the integrals that appear in each coefficient and insert the expansion of ǫ in E obtained from the implicit equation V (a -ǫ) = E:

ǫ = √ 2 ω √ E + V (3) (a) 3ω 4 E - √ 2 3ω 2 V (4) (a) -5 V (3) (a) 2 36ω 7 E 3/2 + O(E 2 ) . (B3)
Proceed in an analogous way for the computation of the first three terms of the asymptotic expansion in √ E for S+ . When summing S + and S -, expression (33b) is obtained with (35).

The expansion of Sc is more subtle since it is not differentiable at E = 0. Its derivative is given by where we denote ǫ def = a -q r (E) and define L(s, ǫ)

def = 2 √ 2 V (a -s) -V (a -ǫ) . ( B5
)
The function L(s, ǫ) is not continuous but we can extract the discontinuous part from

L(s, ǫ) = 4 ω √ s 2 -ǫ 2   1 + n 3 (-1) n 2V (n) (a) n! s n -ǫ n s 2 -ǫ 2   - 1 2
(B6) by expanding the last factor: 3) (a) 3ω 3 s 3 -ǫ 3 (s 2 -ǫ 2 ) 3/2 + M (s, ǫ) (B7) where now M (s, ǫ) is a continuous function of its two variables. Then, a standard theorem in analysis assures that ǫ → a ǫ M (s, ǫ) ds is continuous and its limit when ǫ → 0 is with A given by (34). The two other integrals obtained by inserting (B7) in the right hand side of (B4) can be computed exactly and expanded as ǫ → 0 + up to order o(1). Then, inserting (B3), we obtain

L(s, ǫ) = 4 ω √ s 2 -ǫ 2 + 2V (
d Sc dE = 4 ω ln √ 2E 2aω - 4A ω + o(1) . (B9)
Its integration leads directly to (33a).

when we insert the semiclassical expressions ( 16), we obtain a sum of integrals of the form dq dq ′ e i q b β p β (x, E n ) dx

∂ p H p β (q, E n ), q e - i q ′ b β ′ p β ′ (x, E n ) dx
∂ p H p β ′ (ηq ′ , E n ), q ′ × ∂ 2 S o ∂q i ∂q f (ηq ′ ,q;T ) e iSo(ηq ′ ,q;T )/ . (C3)

The stationary conditions p β (q, E n ) = -∂ qi S o (ηq ′ , q; T ) ; (C4a) ηp β ′ (q ′ , E n ) = ∂ q f S o (ηq ′ , q; T ) (C4b) select the classical trajectories o with energy E n that go from (p i , q i ) = (p β (q, E n ), q) at t(s i ) = 0 to (p f , q f ) = (ηp β ′ (q ′ , E n ), ηq ′ ) at time t(s f ) = T . Then the value of the exponent

S o,β,β ′ (E n , T ) = -E n T + So (b β ′ , b β , E n ) (C5)
depends only on the branches where the starting and ending points lie and not on the precise location of these points on the branches. Since q and q ′ correspond to the same torus, such a trajectory must connect the two symmetric tori for η = -1. At a given E n and T , for a fixed q = q i on the branch β, q f and β ′ are uniquely given and we can make the stationary phase approximation for the integral on q ′ . Then if we insert (A10) at energy E n into (C4) and differentiate it with respect to q or q ′ , we obtain some identities that, with (A11), allow us to simplify the combination of the prefactors and the remaining integral in q turns out to be precisely of the form of the right hand side of (A12). A priori, the domain of integration is included in the domain of the branch β but is not necessarily equal to it because when sliding the starting point on the whole branch β, the endpoint may cross a turning point and correspond to a jump of β ′ . However, we obtain characteristic times that depend only on the geometry of the orbit, not on the number of times the considered branch may be repeated as s goes from s i to s f . As discussed in the case of the double well, if there exist different topological classes of o, each of them being characterised by an ordered sequence of turning points [ρ 1 , ρ 2 , . . . ], we must add such contributions. Then, using directly tr Πn Û (T ) ≃ 2e -iEnT / , we have proven that ∆ n (T ) ∼ →0 T

[ρ1,ρ2,... ]

(-1) µo A β,β ′ (E n , T ) e i So(bβ′ ,b β ,En)/ (C6) where the sum runs over all the sequences of turning points on the section s at energy E n where one canonical variable is maintained real. There must exist for such a sequence, one half symmetric orbit o starting on the branch β of the torus at energy E n , crossing successively all the sequences [ρ 1 , ρ 2 , . . . ] and ending on the branch β ′ at time T . The dimensionless coefficients A β,β ′ (E n , T ) have a -independent modulus of order one and depend only on the geometrical properties of the branches. If some parts of the trajectory are repeated, their repetition numbers do not appear in A

+ q(si) b β p β (x, E n ) dx - q(s f ) b β ′ p β ′ (x, E n ) dx . (C7)
In the case of two branches, the computation of the coefficient can be done exactly using the appropriate choice of phase conventions for the base points b β and A n,β ; we obtain:

A β,β ′ (E n , T ) = -1/2 . (C8)

Figure 1 :

 1 Figure1: (colour online) In the case of a double well potential shown in a), tunnelling can be described by a solution of (14) with a purely real q evolving from one well to the other. In phase-space, this trajectory appears to be a concatenation of two types of curves that join on the (Re q)-axis: 1) a trajectory that lies in the phase-space plane (Re q, Re p) with real variations of t at energy E with the potential V and 2) a trajectory that lies in the phase-space plane (Re q, Im p) with imaginary variations of t at energy -E with the potential -V shown in b). In c) a family of constant energy curves is shown (horizontal green (light gray) for the first type, vertical red (dark gray) for the second type). In d) for a given energy, we show how the three curves l, c, r glue together at the turning points (0, ±qr). S denotes the separatrix Im p = ± 2V (q).

Figure 2 :

 2 Figure 2: (colour online) For V (q) = (q 2a 2 ) 2 with a = 1 and ≃ 1/12, we have plotted in the upper graph a) the real part of ∆0 defined by (9) as a function of the complex T . It becomes constant for Im T large enough for having (6) (1/ω = 1/(2 √ 2a) ≃ .35) regardless of Re T . In the same range, Im ∆0(T )) is negligible compared to ∆E0 ∼ 4.4 10 -10 and condition (8) is largely fulfilled. In the lower graph b), we have computed the "exact" value ∆E0 by direct diagonalization. The black thick solid line corresponds to ln(∆E0) from which we have substracted Λ( ) def = -Sc( ω/2)/(2)+ln( ω/π) in order to emphasize the contribution of the prefactor. The blue solid line corresponds to ln |∆0| -Λ( ) for T ≃ -4i and become indistinguishable from the previous one for 1/ 3. With the same substraction, our first proposal(39) gives the constant (ln √ π) ≃ 0.6 and our second one(44), which is the same as the formula in Landau and Lifshitz [33, sec. 50, problem 3], corresponds to 0 (dotted line). Garg's formula[27, eq.1.1] gives the constant (ln π/e ≃ 0.07) (dash-dotted line).

Figure 4 :

 4 Figure4: (colour online) With the choice of a complextime path given in 3) (one-step path), the contributions to tr Ŝ Πn Û (T ) come from one half symmetric periodic orbits o (b) that are a composite of 1) (repetitions of) a periodic orbit r of period Tr in the phase-space plane (Re q, Re p) at energy E > γ (recall that the planes q = ±π are identified); 2) one half of a periodic orbit c of period Tc in the phase-space plane (Im q, Re p) at energy E.

√

  The asymptotic expansion of Sc (E) for large E leads to ∆ n ∼ 2 π 2E n e 2(ln(γ/En)+2-3 ln 2)

Figure 5 :

 5 Figure 5: For the pendulum with γ = 1, we have plotted in a) the exact value of ln(∆En) (see eq. (47)) versus 1/ for n = 1, . . . 9 (black lines). The dots are given by the semiclassical approximation ln(∆n) (see eq. (58b)) up to the maximum value of 1/ (vertical dotted lines) where E + n become lower that the separatrix energy γ. The semiclassical limit is obtained when the states become more and more localised in the region of phase space corresponding to rotations: for a fixed n this requires to increase En ≃ n 2 2 /2 that is to increase ; or for a fixed energy E > γ, we must increase n ∝ √ 2E/ . To control the validity of the prefactors, we can check in b) that ln(∆En/∆n) approaches zero in the semiclassical limit.

Figure 6 :

 6 Figure 6: (colour online) To describe resonant tunnelling between two symmetric wells centered at ±a and separated by a third central well (a), we use orbits made of five primitive trajectories with a time path given in (d). We still maintain the position q to be real and can use a three dimensional section of phase-space (see text). Here we have chosen V given by (65) with a = 7/4 and b = 1/2.

def=Figure 7 :

 7 Figure 7: (colour online) For V given by (65) with a = 7/4 and b = 1/2, we plot ln |∆En| (thick black line) for the lowest doublet in the symmetric wellsE ± n ≃ ω/2 with ω = 2a 2(a 2b 2 ) ≃ 6.30. The (blue) dots correspond to estimation (A13) with T = |T |e -iθ = (K + 1/2)Tm(En) -iTc(En) in the limit K → ∞ (θ → 0). The upper thin (red) curve provides ln |∆En/∆n|. The inset provides a magnification of the third spike around 1/ = 6.6 for which Tm(En) = 2Tr(En) ≃ 1.58 and Tc(En) ≃ 1.50. Here, we have also plotted ln |∆n(T )| for K = 100 (thin orange (light gray), θ ≃ 0.01) K = 30 (thin green θ ≃ 0.03), K = 5 (dashed thin magenta θ ≃ 0.17).

Figure 8 :

 8 Figure 8: (colour online)The escape rate at energy E from an island created by a confining potential V whose shape is shown in a) can be semiclassically computed in terms of complex classical orbits where one canonical coordinate, say q, is kept real. In phase-space, this trajectory appears to be a concatenation of two types of curves that join on the (Re q)-axis: 1) a trajectory that lies in the phase-space plane (Re q, Re p) with real variations of t at energy E with the potential V and 2) a trajectory that lies in the phase-space plane (Re q, Im p) with imaginary variations of t at energy -E with the potential -V shown in b). In c) a family of constant energy curves is shown (horizontal green for the first type, vertical red for the second type). In d) for a given energy, we show how three curves l, c, r glue together at the turning points.

  β,β ′ (E n , T ) but only in the cumulative quantities: the index µ o and the action So (bβ ′ , b β , E n ) given by So (b β ′ , b β , E n ) =

  The sum concerns all s (different energies may be possible) and [ρ 1 , ρ 2 , . . . ] such that we can construct on s, with an appropriate choice of [t], a periodic orbit o if η = +1 (a half symmetric periodic orbit o if η = -1) of period T o = T . The branch β is the one where o starts, the sequence of turning points that are successively crossed by o is exactly [ρ 1 , ρ 2 , . . . ].

Acknowledgments

We have received a lot of benefits from discussions with Olivier Brodier, Dominique Delande, Akira Shudo, Denis Ullmo and Jean Zinn-Justin, with special thanks to Stephen Creagh, who has had a good intuition about these issues for a long time and has shared his insights with us. We acknowledge Stam Nicolis for his careful reading of the manuscript and Olivier Thibault for his efficient skills in terms of computer maintenance. One of us (A.M.) is grateful to the Laboratoire Kastler Brossel for its hospitality.

We illustrate in the main body of this article, how to compute the sum in the right hand side of (C6).

Appendix A

In this appendix, we derive the dominant contributions (19) to the traces involved in (9). Both cases tr Û (T ) and tr Ŝ Û (T ) will be treated simultaneously by defining the sign η to be +1 in the first case and -1 in the second case. The semiclassical arguments underpinning the derivation are relatively standard and may be found in one way or another in the literature. For instance, the contribution (A13) for η = 1 can be found in [41, eq. (2.12)] within a more restricted context (see also [42, chap. 8]). Nevertheless we found it useful to provide all the steps in the precise context of this work, not only to render the presentation self-contained, but also because we are working in the time domain with general Hamiltonians that have not necessarily the form [START_REF]Advanced Series in Mathematical Physics[END_REF].

Given a [t], for any classical phase-space path o, we can consider the final coordinates (p f , q f ) at time t(s f ) = T as smooth functions of the initial coordinates (p i , q i ) at time t(s i ) = 0. The monodromy matrix M o is defined as the differential of these functions: to first order, an initial small perturbation implies the final perturbation

The sub-matrices can be expressed from the second derivatives of S o (q f , q i ; T ) by differentiating equations (18):

which is nothing but the expression that det M o = 1 once we use the identity

where (A, B, C, D) are square matrices of the same size and C is invertible.

Coming back to the oscillating integral (17), if the integration path can be deformed in order to pass through an isolated critical point q c of q → S o (ηq, q, T ), we will have a contribution of the form det

where all the functions are evaluated at (ηq c , q c , T ). Physically q c is interpreted as the initial position of the phase space paths o such that (p f , q f ) = (ηp i , ηq i ). With the help of (A3), the prefactor of the exponential can be written as det

The quasi-mode

2 is localised on one torus at energy E + n ≃ E - n . Standard jwkb techniques [START_REF] Keller | [END_REF]44] provide a semiclassical approximation to its wave function:

(β labels the possible several branches of the torus, A n,β are dimensionless coefficients of unit modulus, b β is a base point of the branch β and T β (E n ) the characteristic time (A12) spent on the branch β). Within the semiclassical approximation, it is, therefore, consistent to construct Π n (q ′ , q) by substituting (C1) in the matrix elements of the projector operator |φ (20) 2 dq dq ′ Φ s.c. n (q) Φ s.c. n (q ′ ) * G(ηq ′ , q; T ) , (C2)