Clotilde Fermanian Kammerer

Caroline Lasser
email: classer@ma.tum.de

Didier Robert
email: didier.robert@univ-nantes.fr

Caroline Lasser Didier

Robert To

Ari Laptev

C Fermanian Kammerer

Adiabatic and non-adiabatic evolution of wave packets and applications to initial value representations

est destinée au dépôt

Introduction

We consider semi-classical systems of the form

iε∂ t ψ ε (t) = H(t)ψ ε (t), ψ ε |t=t 0 = ψ ε 0 (1.1)
where (ψ ε 0 ) is a bounded family in L 2 (R d , C N ), ψ ε 0 L 2 (R d ,C N ) = 1 and Ĥ(t) is the ε-Weyl quantization of a smooth Hermitian symbol H(t, x, ξ ) satisfying suitable growth assumptions. We are interested in approximate realizations of the unitary propagator associated with this equation relying on the use of continuous superpositions of Gaussian wave packets. Before explaining these ideas in more detail (in particular the definition of the quantization that the reader will find below), let us first emphasize different types of systems that we have in mind.

The simplest case is the one where H(t) is independent of (x, ξ ), which implies that the operator Ĥ(t) coincides with the matrix H(t). In this case it is known that for N ≥ 2 new phenomena appear by comparison with the scalar case N = 1. When the eigenvalues of H(t) are not crossing, the adiabatic theorem says that if ψ ε 0 is an eigenfunction of H(t 0 ), then at every time ψ ε (t) has an asymptotic expansion for ε → 0, and the leading term is an eigenfunction of H(t). The first complete proof was given by Friedrichs ([9] after first results by Born-Fock [START_REF] Born | Beweis des Adiabatensatzes[END_REF] and Kato [START_REF] Kato | On the adiabatic theorem of quantum mechanics[END_REF]. In [9, Part II], Friedrichs considered a non-adiabatic situation where H(t) is a 2 × 2 Hermitian matrix with two analytic eigenvalues h ± (t) such that h + (t) = h -(t) for t = 0 and

(h + -h -)(0) = 0, d dt (h + -h -)(0) = 0.
In this spirit, we here consider a toy model for space dependent Hamiltonians of the following form: We choose d = 1, N = 2, θ ∈ R + , k ∈ R * , and

H k,θ = ε i d dx I C 2 + kx 0 e iθ x e -iθ x 0 . (1.2)
Its semiclassical symbol H k,θ (x, ξ ) = ξ + kx 0 e iθ x e -iθ x 0 has the eigenvalues and associated eigenvectors (the latter depending only in x)

h ± (x, ξ ) = ξ ± kx, V ± (x) = 1 √ 2 e iθ x ±1 (1.3)
So for k = 0 we have a crossing at x = 0 that we call smooth crossing because there exists smooth eigenvalues and eigenprojectors. We shall use this toy-model all along this article. Notice that one can prove that the operator H k,θ is essentially self-adjoint in L 2 (R, C 2 ) by using a commutator criterion with the standard harmonic oscillator (see [START_REF] Maspero | On time dependent Schrödinger equations: global well-posedness and growth of Sobolev norms[END_REF]Appendix A]). As we shall see later, the solutions of the Schrödinger equation for H k,θ can be computed by solving an ODE asymptotically as ε → 0, see (5.4), and by using non-adiabatic asymptotic results from Friedrichs [START_REF] Friedrichs | On the adiabatic Theorem in Quantum Theory[END_REF] and Hagedorn [START_REF] Hagedorn | Adiabatic Expansions near Eigenvalue Crossings[END_REF].

More generally, of major interest because of their applications to molecular dynamics, are the Schrödinger Hamiltonians with matrix-valued potential considered by Hagedorn in his monograph [START_REF] Hagedorn | Molecular Propagation through Electron Energy Level Crossings[END_REF]Chapter 5],

H S = - ε 2 2 ∆ x I C N +V (x), V ∈ C ∞ (R d , C N×N ), (1.4)
or the models arising in solid state physics in the context of Bloch band decompositions studied by Watson and Weinstein in [START_REF] Watson | Wavepackets in inhomogeneous periodic media: propagation through a one-dimensional band crossing[END_REF].

H A = A(-iε∇ x ) +W (x)I C 2 , A ∈ C ∞ (R d , C N×N ), W ∈ C ∞ (R d , C). (1.5)
The notation H refers to the Weyl quantization that we shall use extensively in this article. For a ∈ C ∞ (R 2d , C N,N ) with adequate control on the growth of derivatives, the operator a is defined by its action on functions f ∈ S (R d , C N ):

op w ε (a) f (x) := a f (x) := (2πε) -d R d a x + y 2
, ξ e i ε ξ •(x-y) f (y)dy dξ .

It turns out that the propagator U ε H (t,t 0 ) associated with H is well defined according to [START_REF] Maspero | On time dependent Schrödinger equations: global well-posedness and growth of Sobolev norms[END_REF]Theorem 5.15] provided that the map (t, z) → H(t, z) is in C ∞ (R × R 2d , C N×N ) valued in the set of self-adjoint matrices and that it has subquadratic growth, i.e.

∀α ∈ N 2d , |α| ≥ 2, ∃C α > 0, sup (t,z)∈R × R 2d ∂ α H(t, z) C N,N ≤ C α . (1.6)
Using the commutator methods of [START_REF] Maspero | On time dependent Schrödinger equations: global well-posedness and growth of Sobolev norms[END_REF], we could extend main of our results to a more general setting. However, the assumptions (1.6) are enough to guarantee the existence of solutions to equation (1.1) in L 2 (R d , C N ). One of our objectives here is to describe different Gaussian-based approximations of the semi-group U ε H (t,t 0 ) in the limit ε → 0.

Let g ε

z denote the Gaussian wave packet centered in z = (q, p) ∈ R 2d with standard deviation √ ε:

g ε z (x) = (πε) -d/4 exp - 1 2ε |x -q| 2 + i ε p • (x -q) . (1.7)
The family of wave packets (g ε z ) z∈R 2d forms a continuous frame and provides for all square integrable functions f ∈ L 2 (R d ) the reconstruction formula

f (x) = (2πε) -d z∈R 2d g ε z , f g ε z (x)dz. (1.8)
The leading idea is then to write the unitary propagation of general, square integrable initial data

ψ ε 0 ∈ L 2 (R d ) as U ε H (t,t 0 )ψ ε 0 = (2πε) -d z∈R 2d g ε z , ψ ε 0 U ε H (t,t 0 )g ε z dz, (1.9)
and to take advantage of the specific properties of the propagation of Gaussian states to obtain an integral representation that allows in particular for an efficient numerical realization of the propagator.

Such a program has been completely accomplished in the scalar case (N = 1). It appeared first in the 80's in theoretical chemistry [START_REF] Herman | A semiclassical justification for the use of non-spreading wavepackets in dynamics calculations[END_REF][START_REF] Kay | Integral expressions for the semi-classical time-dependent propagator[END_REF][START_REF] Kay | The Herman-Kluk approximation: derivation and semiclassical corrections[END_REF] and has led to the socalled Herman-Kluk approximation. The mathematical proof of the convergence of this approximation is more recent [START_REF] Swart | A mathematical justification for the Herman-Kluk Propagator[END_REF][START_REF] Robert | On the Herman-Kluk Semiclassical Approximation[END_REF]. Here, we revisit these results in Section 2 and Section 5, and discuss some extensions to the case of systems (N 1), first for the gapped situation in Section 3 and then for smooth crossings in Section 4.

Propagation of Gaussian states and the Herman-Kluk approximation for scalar equations

In this section, N = 1 and the equation (1.1) is associated with a scalar Hamiltonian H(t) = h(t) of subsquadratic growth (1.6). In that case, the approximate propagation of Gaussian states is described by classical quantities, leading to a simple Herman-Kluk approximation. We set

J = 0 I R d -I R d 0
and for z 0 ∈ R 2d we consider the classical Hamiltonian trajectory z(t) = (q(t), p(t)) defined by the ordinary differential equation

ż(t) = J∂ z h(t, z(t)), z(t 0 ) = z 0 .
The associated flow map is then denoted by

Φ t,t 0 h (z 0 ) = z(t) = (q(t), p(t)), (2.1)
and the blocks of its Jacobian matrix

F(t,t 0 , z 0 ) = ∂ z Φ t,t 0 h (z 0 ) by F(t,t 0 , z 0 ) = A(t,t 0 , z 0 ) B(t,t 0 , z 0 ) C(t,t 0 , z 0 ) D(t,t 0 , z 0 ) .
We note that the Jacobian satisfies the linearized flow equation

∂ t F(t,t 0 , z 0 ) = JHess z h(t, z(t)) F(t,t 0 , z 0 ), F(t 0 ,t 0 , z 0 ) = I R 2d . (2.2)
Thus, F is a smooth in t,t 0 , z with any derivative in z bounded. We will also use the action integral

∂ t S(t,t 0 , z 0 ) = p(t) • q(t) -h(t, z(t)), S(t 0 ,t 0 , z 0 ) = 0. (2.3)
Then, the Herman-Kluk approximation (also called frozen Gaussian approximation in the literature) of the unitary propagator U ε h (t,t 0 ) writes as follows.

Theorem 2.1 ( [START_REF] Swart | A mathematical justification for the Herman-Kluk Propagator[END_REF][START_REF] Robert | On the Herman-Kluk Semiclassical Approximation[END_REF]). Assuming (1.6), the evolution through the scalar equation with N = 1 and H = hI C in (1.1) satisfies for every J ≥ 0,

U ε h (t,t 0 ) = I ε,J h (t,t 0 ) + O(ε J+1 )
in the norm of bounded operators on L 2 (R d ), where the Herman-Kluk propagator is defined for all ψ ∈ L 2 (R d ),

I ε,J h (t,t 0 )ψ = (2πε) -d R 2d g ε z , ψ u ε,J (t,t 0 , z)e i ε S(t,t 0 ,z) g ε Φ t,t 0 h (z) dz, (2.4
)

with u ε,J (t,t 0 , z) = ∑ 0≤ j≤J ε j u j (t,t 0 , z),
where every u j is smooth in t,t 0 , z, with any derivative in z bounded. The function u 0 is the Herman-Kluk prefactor, u 0 (t,t 0 , z) = 2 -d/2 det 1/2 (A(t,t 0 , z) + D(t,t 0 , z) + i(C(t,t 0 , z) -B(t,t 0 , z))) , (2.5)

which has the branch of the square root determined by continuity in time.

Let us remark that the Gaussian wave packets in (2.4) all have the same covariance matrix Γ = I C d , that is, in the terminology put forward by E. Heller [START_REF] Heller | Frozen Gaussians: A very simple semiclassical approximation[END_REF], the Gaussians are frozen. The first statement of the form (2.4) at leading order (J = 0) is due to M. Herman and E. [START_REF] Herman | A semiclassical justification for the use of non-spreading wavepackets in dynamics calculations[END_REF] [START_REF] Herman | A semiclassical justification for the use of non-spreading wavepackets in dynamics calculations[END_REF]. Then, W.H. Miller (2002) [START_REF] Willliam | An alternate derivation of the Herman-Kluk (coherent state) semiclassical initial value representation of the time evolution operator[END_REF] noticed that (2.4) can be deduced from the Van Vleck approximation of the propagator (hence related with a Feynman integral). In more recent work [START_REF] Swart | A mathematical justification for the Herman-Kluk Propagator[END_REF][START_REF] Robert | On the Herman-Kluk Semiclassical Approximation[END_REF], (2.4) was established using Fourier-Integral operator techniques. Here we shall comment in more details on a proof based on the propagation of Gaussian wave packets, via a thawed Gaussian approximation (see Section 5). This proof has the advantage that it can be used in the case of systems of Schrödinger equations, which were not treated in the preceding references.

A Gaussian wave packet is a wave packet with profile function belonging to the class of Gaussian states with variance taken in the Siegel set S + (d) of d × d complexvalued symmetric matrices with positive imaginary part,

S + (d) = Γ ∈ C d×d , Γ = Γ τ , ℑΓ > 0 .
With Γ ∈ S + (d) and z ∈ R, we associate the Gaussian state

g Γ,ε z (x) = c Γ (πε) -d/4 exp i ε p • (x -q) + i 2ε Γ(x -q) • (x -q) ,
where c Γ = det 1/4 (ℑΓ) is a positive normalization constant in L 2 (R d ). We note that the standardized Gaussian defined in (1.7) satisfies g ε z = g i Id,ε z .

We shall use the notation WP ε z (ϕ) to denote the bounded family in L 2 (R d ) associated with ϕ ∈ S (R d ) and z = (q, p) ∈ R 2d by

WP ε z (ϕ)(x) = ε -d/4 e i ε p•(x-q) ϕ x -q √ ε . (2.6)
With wave packet notation, the above Gaussian states satisfy g Γ,ε z = WP ε z (g Γ,1 0 ).

Theorem 2.2 ( [START_REF] Combescure | Coherent states and applications in mathematical physics[END_REF][START_REF] Robert | Propagation of coherent states in quantum mechanics and applications[END_REF]). Assuming (1.6) with N = 1 and H = hI C , there exists a family of time-dependent Schwartz functions (ϕ j (t)) j∈N such that for all N 0 ∈ N, in L 2 (R d ),

U ε h (t,t 0 )g Γ 0 ,ε z 0 = e i ε S(t,t 0 ,z 0 ) g Γ(t,t 0 ,z 0 ),ε Φ t,t 0 (z 0 ) + N 0 ∑ j=1 ε j/2 WP ε Φ t,t 0 (z 0 ) (ϕ j (t)) + O(ε (N 0 +1)/2 ) (2.7) with Γ(t,t 0 , z 0 ) = (C(t,t 0 , z 0 ) + D(t,t 0 , z 0 )Γ 0 )(A(t,t 0 , z 0 ) + B(t,t 0 , z 0 )Γ 0 ) -1 (2.8)
and c Γ(t,t 0 ,z 0 ) = c Γ 0 det -1/2 (A(t,t 0 , z 0 ) + B(t,t 0 , z 0 )Γ 0 ), where the complex square root is continuous in time.

Note that we have Γ(t 0 ,t 0 , z 0 ) = Γ 0 in the statement above. Actually, explicit information is obtained on the profiles (ϕ j (t)) j∈N , in particular about ϕ 1 (t) (see [START_REF] Combescure | Coherent states and applications in mathematical physics[END_REF]Section 4.1.2] or Proposition 2.3 in [START_REF] Fermanian Kammerer | Propagation of wave packets for systems presenting codimension one crossings[END_REF]). As we shall see in Section 5, this result can be a starting point for proving the Herman-Kluk approximation of Theorem 2.1.

The result above also holds for general wave packets as in (2.6) with profiles that are not necessarily Gaussians (see [START_REF] Combescure | Coherent states and applications in mathematical physics[END_REF]Section 4.1.2]). Moreover, also the norm can be generalized to

f Σ ε k = sup |α|+|β |≤k x α (ε∂ x ) β f L 2 , k ∈ N.
Example 2.3. Both Theorem 2.1 and 2.2 are exact when the Hamiltonian h is quadratic. For the Hamiltonians h ± (z) = p ± kq associated with the toy model (see (1.3)), the classical flow is linear Φ t,t 0 h ± (z 0 ) =: (q ± (t), p ± (t)) = z 0 + (tt 0 )(1, ∓k), z 0 = (q 0 , p 0 ), the width of the Gaussian is constant, Γ ± (t,t 0 , z 0 ) = Γ 0 , and the actions only depend on q 0 and are given by S ± (t,t 0 , q 0 ) = ∓kq 0 (t

-t 0 ) ∓ k 2 (t -t 0 ) 2 .
3 Herman-Kluk formula in the adiabatic setting

We need adaptions for generalizing the scalar ideas to systems. The first one replaces the scalar Herman-Kluk prefactor u 0 (t,t 0 , z 0 ) by a vector-valued one U(t,t 0 , z 0 ) that is expanded in a basis of eigenvectors of H(t, z), taken along the classical trajectory Φ t,t 0 h (z) of the corresponding eigenvalue (denoted here by h(t, z)). This is done by parallel transport, and it is sufficient for an order ε approximation as long as the system is gapped.

Parallel transport

The following construction generalizes [3, Proposition 1.9], which was inspired by the work of G. Hagedorn, see [START_REF] Hagedorn | Molecular Propagation through Electron Energy Level Crossings[END_REF]Proposition 3.1]. The details are given in [START_REF] Fermanian Kammerer | Propagation of wave packets for systems presenting codimension one crossings[END_REF]. In the sequel, we denote the complementary orthogonal projector by Π ⊥ (t, z) = I C N -Π(t, z) and assume that

H(t, z) = h(t, z)Π(t, z) + h ⊥ (t, z)Π ⊥ (t, z) (3.1)
with the second eigenvalue given by h ⊥ (t, z) = tr(H(t, z))h(t, z). We introduce the auxiliary matrices

Ω(t, z) = -1 2 h(t, z) -h ⊥ (t, z) Π(t, z){Π, Π}(t, z)Π(t, z), K(t, z) = Π ⊥ (t, z) (∂ t Π(t, z) + {h, Π}(t, z)) Π(t, z), Θ(t, z) = iΩ(t, z) + i(K -K * )(t, z),
that are smooth and satisfy some algebraic properties. In particular, Ω is skew-symmetric and Θ is self-adjoint: Ω = -Ω * and Θ = Θ * .

Proposition 3.1 ([7]

). Let H(t, z) be a smooth Hamiltonian that satisfies (1.6) and has a smooth spectral decomposition (3.1). We assume that both eigenvalues are of subquadratic growth as well. We consider

V 0 ∈ C ∞ 0 (R 2d , C N ) and z 0 ∈ R 2d such that there exits a neighborhood U of z 0 such that for all z ∈ U V 0 (z) = Π(t 0 , z) V 0 (z) and V 0 (z) C N = 1.
Then, there exists a smooth normalized vector-valued function V (t,t 0 ) satisfying

V (t,t 0 , z) = Π(t, z) V (t,t 0 , z) for all z ∈ Φ t,t 0 h (U), such that for all t ∈ R and z ∈ Φ t,t 0 h (U), ∂ t V (t,t 0 , z) + {h, V }(t,t 0 , z) = -iΘ(t, z) V (t,t 0 , z), V (t 0 ,t 0 , z) = V 0 (z). (3.2)
We note that the results of this Proposition are valid as long as smooth eigenprojectors and eigenvalues do exist. It does not require an explicit adiabatic situation and we will use that observation in Section 4. In the case of Schrödinger systems, one refers to (3.2) as parallel transport because the vectors ∂ t V (t) + {h, V }(t) belong to the range of Π ⊥ (t) at any time of the evolution.

Example 3.2. For the toy model H k,θ (x, ξ ), the auxiliary matrices are

Ω ± = 0, K ± = iθ 4 1 ±e iθ x ∓e -iθ x -1 , Θ ± = θ 2 -1 0 0 1 .
Initiating the parallel transport equation by the eigenvectors given in (1.3), we obtain

V ± (t,t 0 , q) = 1 √ 2 e iθ /2(t-t 0 ) e iθ (q-t+t 0 ) ±e -iθ /2(t-t 0 ) = 1 √ 2 e -i 2 θ (t-t 0 ) e iθ q ±e -i 2 θ (t-t 0 )
. Note that they do not depend on p.

Herman-Kluk approximation in an adiabatic setting

In the context of the preceding section, and in presence of an eigenvalue gap,

∃δ > 0 : dist(h(t, z), h ⊥ (t, z)}) ≥ δ for all (t, z),
it is well-known that one has adiabatic decoupling: If the initial data are scalar multiples of an eigenvector associated with the eigenvalue, then the solution keeps this property to leading order in ε. For an approximation to higher order in ε, one has to consider perturbations of the eigenprojector, the so-called superadiabatic projectors (see [START_REF] Bily | Propagation d'états cohérents et applications[END_REF][START_REF] Teufel | Adiabatic perturbation theory in quantum dynamics Lecture Notes in Mathematics 1821[END_REF] and the new edition of [START_REF] Combescure | Coherent states and applications in mathematical physics[END_REF] that should appear soon). Adiabatic theory implies a Herman-Kluk approximation of the propagator that we state next. Theorem 3.3. [START_REF] Fermanian Kammerer | Propagation of wave packets for systems presenting codimension one crossings[END_REF] In the situation of Proposition 3.1, we assume the existence of an eigenvalue gap for the eigenvalue h of H and consider initial data of the form

ψ ε 0 = V 0 v ε 0 + O(ε) in L 2 (R d ),
where V 0 is a smooth eigenvector and

(v ε 0 ) ε>0 a family of functions uniformly bounded in L 2 (R d , C) Then, in L 2 (R d ), U ε H (t,t 0 )ψ ε 0 = (2πε) -d R 2d g ε z , v ε 0 U(t,t 0 , z) e i ε S(t,t 0 ,z) g ε Φ t,t 0 h (z) dz + O(ε) with U(t,t 0 , z) = u 0 (t,t 0 , z) V (t,t 0 , Φ t,t 0 h (z))
where u 0 (t,t 0 , z) is given by (2.5) and the eigenvector V (t,t 0 , z) by (3.2).

We describe in Section 5 a proof of this result that crucially uses the approximate evolution of a Gaussian state, that is [25, Section 3],

U ε H (t,t 0 )( V 0 g Γ,ε z ) = e i ε S(t,t 0 ,z) V (t,t 0 )g Γ(t,t 0 ,z),ε Φ t,t 0 h (z) 1 + √ ε a(t) • (x -q(t)) √ ε + O(ε) (3.3)
with a(t) = a(t,t 0 , z) a smooth and bounded map. Note that using superadiabatic projectors [START_REF] Robert | Propagation of coherent states in quantum mechanics and applications[END_REF], one can push the asymptotics further and exhibit O(ε) contributions that will have components on the other mode. Here again, a formula similar to (3.3) can be proved for general wave packets [START_REF] Combescure | Coherent states and applications in mathematical physics[END_REF].

Algorithmic realization of the propagator

Numerical realizations of the Herman-Kluk approximation have been first developed in [START_REF] Kluk | Comparison of the propagation of semiclassical frozen Gaussian wave functions with quantum propagation for a highly excited anharmonic oscillator[END_REF] and are still in practical use in theoretical chemistry [START_REF] Werther | Coherent state based solutions of the time-dependent Schrödinger equation: hierarchy of approximations to the variational principle[END_REF]. Below, we follow the more recent account given in [START_REF] Lasser | Discretising the Herman-Kluk Propagator[END_REF][START_REF] Lasser | Computing quantum dynamics in the semiclassical regime[END_REF]. We stay with the assumptions of Theorem 3.3 and consider initial data associated with the (gapped) mode h.

Step 1: Initial sampling. Choose a set of numerical quadrature points

z 1 = (q 1 , p 1 ), . . . , z N = (q N , p N )
with associated weights w 1 , . . . , w N > 0, and evaluate the initial transform g ε z , v ε 0 in these points. This provides an approximation to the initial data,

ψ ε 0 (x) ∼ (2πε) -d ∑ 1≤ j≤N g ε z j , v ε 0 V 0 (z j )g ε z j (x)w j ,
which is of order ε in L 2 (R d ), as long as the chosen quadrature rule is sufficiently accurate.

Step 2: Transport. For each of the points z j , j = 1, . . . , N, compute

(1) the classical trajectories z j (t) defined by (2.1), (2) the eigenvectors V (t,t 0 , z j (t)) along the flow by use of equation (??), (3) the Jacobian matrices F(t,t 0 , z j ) using equation (2.2), (4) the action integrals S(t,t 0 , z j ) using equation (2.3), ( 5) the Herman-Kluk prefactor u 0 (t,t 0 , z j ) using the Jacobians and equation (2.5).

If the time-discretization is symplectic and sufficiently accurate, then the overall accuracy of order ε is not harmed, see [START_REF] Lasser | Discretising the Herman-Kluk Propagator[END_REF]Theorem 2].

Step 3: Conclusion. At the end of these thwo steps, we are left with the Hermann-Kluk quadrature formula:

ψ ε (t, x) ∼ (2πε) -d ∑ 1≤ j≤N g ε z j , v ε 0 U(t,t 0 , z)e i ε S(t,t 0 ,z j ) g ε z j (t) w j .
Of course, the algorithm generalizes to initial data which has several components on separated eigenspaces. In higher dimensional applications, often Monte-Carlo quadrature is used. Then, the initial transform could be written as

(2πε) -d g ε z , v ε 0 dz = r ε 0 (z)µ ε 0 (dz),
where

µ ε 0 (dz) = | g ε z , v ε 0 |dz -1 | g ε z , v ε 0 |dz
is a probability measure and r ε 0 (z) a complex-valued function of L 1 (R 2d , dµ ε 0 ) according to the Radon-Nikodym Theorem. The quadrature nodes z 1 , . . . , z N are then chosen independently and identically distributed according to the measure µ ε 0 , while the weights are all the same, w j = 1/N for all j.

What about smooth crossings ?

For simplicity we assume here that the Hermitian matrix H(t, z) is 2 × 2 (N = 2) (the same results are also proved for two eigenvalues with arbitrary multiplicity [START_REF] Fermanian Kammerer | Propagation of wave packets for systems presenting codimension one crossings[END_REF]). We assume that it has smooth eigenvalues h 1 (t, z) and h 2 (t, z) and smooth eigenprojectors Π 1 (t, z) and Π 2 (t, z). To ensure control on the derivatives of the eigenprojectors, we suppose a non crossing assumption at infinity: there exist c 0 , n 0 , r 0 > 0 such that

|h 1 (t, z) -h 2 (t, z)| c 0 z -n 0 for all (t, z) with |z| r 0 ,
where we denote z = (1 + |z| 2 ) 1/2 . We assume the following form of the matrix:

Asumption (SC). There exist scalar functions v, f ∈ C ∞ (R 2d+1 , R) and a vector- valued function u ∈ C ∞ (R 2d+1 , R 3 ) with |u(t, z)| = 1 for all (t, z) such that H(t, z) = v(t, z)Id + f (t, z) u 1 (t, z) u 2 (t, z) + iu 3 (t, z) u 2 (t, z) -iu 3 (t, z) -u 1 (t, z) .
Besides, the crossing is generic in ϒ in the sense that

(∂ t f + {v, f })(t , z ) = 0, ∀(t , z ) ∈ ϒ.
In that case, the crossing set ϒ is a submanifold of codimension one, and all the classical trajectories that reach ϒ are transverse to it. In contrast to the previous adiabatic situation, initial data associated with one eigenspace generates a component on the other eigenspace, that is larger than the adiabatic O(ε), namely O( √ ε). Starting at time t 0 with a Gaussian wave packet, that is associated with the eigenvalue h 1 and localized far from the crossing set ϒ, an approximation of the form (3.3) holds as long as the trajectory does not reach ϒ. The apparition of a √ ε contribution on the other mode then occurs exactly on ϒ. One can interpret this phenomenon in terms of hops: starting at time t 0 from some point z 0 / ∈ ϒ for which the Hamiltonian trajectory z 1 (t,t 0 ) = Φ t,t 0 1 (z 0 ) passes through the crossing at time t = t (z 0 ) and point z = z 1 (t ,t 0 ), we generate a new trajectory Φ t,t 2 (z ) associated with the mode h 2 . This results in the construction of a hopping trajectory that hops from one mode to the other one at time t .

Such an interpretation is crucial for describing the dynamics of systems with eigenvalue crossings. It has been introduced around the 70s in the chemical literature for avoided and conical crossings of eigenvalues and has been widely used since then (see [START_REF] Tully | Perspective: Nonadiabatic dynamics theory[END_REF]). The first mathematical results on the subject are more recent and analyse the propagation of Wigner functions through singular crossings (see [START_REF] Kammerer | Propagation through generic level crossings: a surface hopping semigroup[END_REF][START_REF] Kammerer | An Egorov Theorem for avoided crossings of eigenvalue surfaces[END_REF]).

We now aim at a precise description of these new contributions in the case of smooth crossing, first for initial data that are Gaussian wave packets, then we lift this result to a Herman-Kluk formula for smooth crossings in the case of the toy model.

Wave packet propagation

Let us now give a precise statement for the propagation of wave packets through a smooth crossing. We start with initial data of the form

ψ ε 0 = V 0 v ε 0 , v ε 0 = g Γ 0 ,ε z 0 with H(t 0 , z) V 0 (z) = h 1 (t 0 , z) V 0 (z)
in a neighborhood of z 0 . We use Proposition 3.1 to construct two families of time-dependent eigenvectors: ( V 1 (t, z)) t≥t 0 is associated with the eigenvalue h 1 (t, z) and initial data at time t 0 given by V 1 (t 0 , z) = V 0 (z), while ( V 2 (t, z)) t≥t is constructed for t t (the crossing time introduced in the previous paragraph) for the eigenvalue h 2 (t, z) with initial data at time t = t (z 0 ) satisfying

V 2 (t , z) = -γ(t , z) -1 Π 2 (∂ t Π 2 + {v, Π 2 }) V 1 (t , z) with γ(t , z) = (∂ t Π 2 + {v, Π 2 }) V 1 (t , z) C N .
We next introduce a family of transformations, which describes the non-adiabatic effects for a wave packet that passes the crossing. For parameters (µ, α, β

) ∈ R × R 2d and ϕ ∈ S (R d ), we set T µ,α,β ϕ(y) = +∞ -∞ e i(µ-α•β /2)s 2 e isβ •y ϕ(y -sα)ds.
This operator maps S (R d ) into itself if and only if µ = 0. Moreover, for µ = 0, it is a metaplectic transformation of the Hilbert space L 2 (R d ), multiplied by a complex number. In particular, for any Gaussian function g Γ , the function T µ,α,β g Γ is a Gaussian:

T µ,α,β g Γ = c µ,α,β ,Γ g Γ µ,α,β ,Γ ,
where Γ µ,α,β ,Γ ∈ S + (d) and c µ,α,β ,Γ ∈ C can be computed explicitly (see [7, Appendix E]).

Combining the parallel transport for the eigenvector and the metaplectic transformation for the non-adiabatic transitions, we obtain the following result. Theorem 4.2 (Propagation through a smooth crossing). Let Assumption (SC) on the Hamiltonian matrix H(t) hold and that the crossing is generic. Assume that the initial data (ψ ε 0 ) ε>0 are wave packets as above. Let T > 0 be such that the interval [t 0 ,t ] is strictly included in the interval [t 0 ,t 0 + T ]. Then, for all k ∈ N there exists a constant C > 0 such that for all t ∈ [t 0 ,t ) ∪ (t ,t 0 + T ] and for all ε ≤ |tt | 9/2 ,

ψ ε (t) -V 1 (t)v ε 1 (t) - √ ε1 t>t V 2 (t)v ε 2 (t) L 2 (R d ) ≤ C ε m ,
with an exponent m 5/9. The components of the approximate solution are

v ε 1 (t) = U ε h 1 (t,t 0 )g Γ 0 ,ε z 0 and v ε 2 (t) = U ε h 2 (t,t )v ε 2 (t ) with v ε 2 (t ) = γ e iS /ε WP ε z (T ϕ 0 (t )),
where ϕ 0 (t) = g Γ 1 (t,t 0 ,z),1 0

is the leading order profile of the coherent state v ε 1 (t) given by Theorem 2.2, and

γ = γ(t , z ) = ({v, Π 2 } + ∂ t Π 2 ) V 1 (t , z ) C N .
The transition operator T = T µ ,α ,β is defined by the parameters

µ = 1 2 (∂ t f + {v, f }) (t , z ) and (α , β ) = Jd z f (t , z ).
The constant C = C(T, k, z 0 , Γ 0 ) > 0 is ε-independent but depends on the Hamiltonian H(t, z), the final time T , and on the initial wave packet's center z 0 and width Γ 0 .

Note that by the transversality assumption we have µ = 0, which guarantees that T ϕ 0 (t ) is Schwartz class. The coefficient γ quantitatively describes the distortion of the eigenprojector Π 1 (t) during its evolution along the flow generated by h 1 (t). If the matrix H(t, z) is diagonal (or diagonalizes in a fixed orthonormal basis that is (t, z)-independent), then γ = 0: the equations are decoupled (or can be decoupled), and one can then apply the result for a system of two independent equations with a scalar Hamiltonian and, of course, there is no interaction of order √ ε between the modes.

The previous theorem extends to more general wave packets as defined in (2.6) and also holds with respect to Σ ε k -norms for k ∈ N (see [START_REF] Fermanian Kammerer | Propagation of wave packets for systems presenting codimension one crossings[END_REF]Theorem 3.8]). As mentioned alongside the proof [START_REF] Fermanian Kammerer | Propagation of wave packets for systems presenting codimension one crossings[END_REF], the argument contains the germs for a full asymptotic expansion in powers of √ ε (with log ε corrections).

Example 4.3. Notice that for the toy model H k,θ , we have at any point of

ϒ = {x = 0}, µ = k 2 , α = 0, β = -k, γ = |θ | 2 ,
and T ϕ(y) = 2π ik e k 2i y 2 ϕ(y) for all ϕ ∈ S (R) and all y ∈ R. Besides, if t 0 ≤ 0, the trajectories of the minus mode that reach ϒ are those arising from points z = (q, p) with q < 0. One then has t = t 0q, p = pkq and the trajectory on the plus mode issued from z = (0, p ) is

Φ t,t + (0, p ) = (t -t , p -k(t -t ) = (t -t 0 + q, p -2kq -k(t -t 0 )).

Towards a Herman-Kluk approximation

The preceding result implies that the leading order of the propagation is still driven by the modes in which the initial data had been taken and we have an Herman-Kluk formula similar to the one obtained in the adiabatic regime, however, with a remainder which is worse. One can conjecture that a more accurate Herman-Kluk approximation holds in a weaker sense (see [START_REF] Fermanian Kammerer | Herman-Kluk approximation for systems presenting codimension one crossings[END_REF]). We define the operator I ε sc (t) by its actions on functions of the form

ψ ε 0 = V 0 v ε 0 + O( √ ε) in L 2 (R d ) with v ε 0 ∈ L 2 (R d ) as I ε sc (t)ψ ε 0 (x) = (2πε) -d R 2d g ε z , v ε 0 U 1 (t,t 0 , z) e i ε S 1 (t,t 0 ,z) g ε Φ t,t 0 h 1 (z) dz + √ ε (2πε) -d × R 2d 1 t>t (z) g ε z , v ε 0 U 2 (t,t (z), z)e i ε S 1 (t (z),t 0 ,z)+ i ε S 2 (t,t (z),z (z)) g ε Φ t,t (z) h 2 (z ) dz with z (z) = Φ t (z),t 0 h 1
(z) and with some adequate formula (taking into account the trans-

fer coefficients γ (z)) for the prefactor U 2 (t,t (z), z) = v 2 (t,t (z), z) V 2 (t,t (z), z (z)).
The conjecture is, that if Assumptions (SC) are satisfied, then, for all

χ ∈ C ∞ 0 (R), in L 2 (R d ), one has R χ(t) (I ε sc (t)ψ ε 0 -U ε H (t,t 0 )ψ ε 0 ) dt = o( √ ε) (4.1)
Estimates that are "averaged in time" have been previously obtained for systems (see [START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF][START_REF] Kammerer | Propagation through generic level crossings: a surface hopping semigroup[END_REF][START_REF] Kammerer | An Egorov Theorem for avoided crossings of eigenvalue surfaces[END_REF] for example). They correspond to an observation period that is a short, but non negligible, time interval. An approximation as (4.1) can be easily proved for the toy-model (1.2) (see Section 5.5) below. The proof in the general case is work in progress [START_REF] Fermanian Kammerer | Herman-Kluk approximation for systems presenting codimension one crossings[END_REF]. It involves more refined estimates than those of Theorem 3.3, that we present in the next section.

The authors believe that the technics used for treating the apparition of a new contribution when trajectories reach a hypersurface (the crossing one in this special case) will be useful for developing Herman-Kluk approximations for avoided crossings and conical ones, using the hopping trajectories of [START_REF] Kammerer | An Egorov Theorem for avoided crossings of eigenvalue surfaces[END_REF] and [START_REF] Kammerer | Propagation through generic level crossings: a surface hopping semigroup[END_REF] respectively. However, we point out, that conical crossings pose the additional difficulty that Gaussians states do not remain Gaussian, even at leading order, as emphasized in [START_REF] Hagedorn | Molecular Propagation through Electron Energy Level Crossings[END_REF].

A sketchy proof for Herman-Kluk approximations

We consider here the scalar and the adiabatic case and we discuss the proof of Theorems 2.1 and 3.3.

The proof strategy

As mentioned in the introduction, the underlying idea for constructing Gaussian based approximations for unitary propagators is to start form equation (1.9). Let us develop a proof strategy based on this idea.

Step 1. For each z ∈ R 2d , we build a thawed wave packet approximation ψ ε th,z (t) to the Schrödinger system (1.1) with initial data

ψ ε |t=t 0 = ψ ε 0 = g ε z for N = 1, V 0 g ε z for N > 1.
We prove that ψ ε th,z (t) satisfies an evolution equation of the form

iε∂ t ψ ε th,z (t) = H(t)ψ ε th,z (t) + ε 2 R ε z (t), ψ ε th,z (t 0 ) = ψ ε 0 (5.1)
with source term R ε z (t).

Step 2. For the thawed Gaussian propagation of general initial data

ψ ε |t=t 0 = ψ ε 0 = v ε 0 for N = 1, V 0 v ε 0 for N > 1, with v ε 0 ∈ L 2 (R d , C) we consider I ε th (t,t 0 )ψ ε 0 = (2πε) -d z∈R 2d g ε z , v ε 0 ψ ε th,z (t) dz.
The error e ε th (t

) = U ε H (t,t 0 )ψ ε 0 -I ε th (t,t 0 )ψ ε 0 satisfies the evolution equation iε∂ t e ε th (t) = H(t)e ε th (t) + ε 2 Σ ε th (t), e ε th (t 0 ) = 0 with source term Σ ε th (t) = (2πε) -d R 2d g ε z , ψ ε 0 R ε z (t)dz.
Since H(t) is self-adjoint, the usual energy argument provides

e ε th (t) ε t t 0 Σ ε th (s) ds.
Step 3. For analysing the source term Σ ε th (t), we consider the integral operator

ψ → (2πε) -d R 2d g ε z , ψ R ε z (t)dz
and its Bargmann kernel

k B (t; X,Y ) = (2πε) -2d R 2d g ε z , g ε Y g ε X , R ε z (t) dz, X,Y ∈ R 2d .
We aim at establishing constants C 1 (t),C 2 (t) > 0, that do do not depend on the semiclassical parameter ε, such that sup

X R 2d |k B (t; X,Y )|dY C 1 (t), sup Y R 2d |k B (t; X,Y )|dX C 2 (t), (5.2)
since then, by the Schur test,

Σ ε th (t) C 1 (t)C 2 (t) ψ ε 0 .
Step 4. For systems, that is, for N > 1, we use the additional observation that

I ε th (t,t 0 )ψ ε 0 = (2πε) -d z∈R 2d g ε z , v ε 0 V (t,t 0 , Φ t,t 0 h (z))e i ε S(t,t 0 ,z) g Γ(t,t 0 ,z),ε Φ t,t 0 h (z)
dz + O(ε).

Step 5. We turn the thawed propagation in a frozen one, in proving that

I ε th (t,t 0 ) = I ε fr (t,t 0 ) + O(ε) in the norm of bounded operators on L 2 (R d ),
where the frozen propagator is defined by the Herman-Kluk formula

I ε fr (t,t 0 )ψ ε 0 = (2πε) -d z∈R 2d g ε z , v ε 0 U(t,t 0 , z) e i ε S(t,t 0 ,z) g ε Φ t,t 0 h (z) dz with U(t,t 0 , z) = u 0 (t,t 0 , z) for N = 1, u 0 (t,t 0 , z) V (t,t 0 , Φ t,t 0 h (z)) for N > 1.
Once the previous steps have been carried out, we have proven the basic J = 0 version of the scalar Herman-Kluk formula of Theorem 2.1 and its generalization to the adiabatic situation given in Theorem 3.3.

The thawed remainder

We first consider scalar wave packet propagation as described in Theorem 2.2 with an accuracy of order ε, that is, for N 0 = 1. The corresponding thawed Gaussian wave packet ψ ε th,z (t) that is defined by the right hand side of (2.7) satisfies an evolution equation of the form (5.1) with a source term

R ε z (t) = e i ε S(t,t 0 ,z) op w ε (L z (t,t 0 ))g Γ(t,t 0 ,z),ε Φ t,t 0 h (z)
, where w → L z (t,t 0 , w) is a smooth function, that is polynomially bounded. It depends on the remainder of Taylor expansions of h(t, •) around the point Φ

t,t 0 h (z), see [4, Sec- tion 4.3.1].
For adiabatic propagation by systems with eigenvalue gaps, as presented in Theorem 3.3, we work with

ψ ε th,z (t) = e i ε S(t,t 0 ,z) V (t,t 0 ) 1 + √ ε a(t,t, 0 , z) • x -q(t,t 0 , z) √ ε g Γ(t,t 0 ,z),ε Φ t,t 0 h (z)
, where the vector a(t,t 0 , z; x) can be constructed explicitly. This wave packet satisfies an evolution equation of the form (5.1) with source term

R ε z (t) = e i ε S(t,t 0 ,z) op w ε ( L z (t,t 0 ))g Γ(t,t 0 ,z),ε Φ t,t 0 h (z)
,

where the vector-valued function w → L z (t,t 0 , w) contains remainder terms of Taylor expansions of h(t, •) around the classical trajectory. We note that L z (t,t 0 , •) has contributions both in the range of V (t,t 0 , •) but also in the orthogonal complement.

The Schur estimate

We now analyse the Bargmann kernel of the source term Σ ε th (t). Since

| g ε z , g ε Y | = e -|Y -z| 2 4ε , we have |k B (t; X,Y )| (2πε) -2d R 2d e -|Y -z| 2 4ε | g ε X , R ε z (t) | dz.
Hence, the crucial estimate that is required concerns the Bargmann transform of the remainder R ε z (t). We write the remainder as

R ε z (t) = e i ε S(t,t 0 ,z) op w ε (L z (t,t 0 ))g Γ z ,ε
Φ z , where w → L z (t,t 0 , w) is a smooth function on phase space, that grows at most polynomially, and Γ z = Γ(t,t 0 , z), Φ z = Φ t,t 0 (z) are short-hand notations for the classical quantities defined in (2.1) and (2.8), respectively. We express the Bargmann transform as a phase space integral

g ε X , R ε z (t) = e i ε S(t,t 0 ,z) R 2d L z (t,t 0 , w) Wig(g ε X , g Γ z ,ε Φ z )(w) dw
with respect to the cross-Wigner function of two Gaussian wave packets with different centers and different width. One can prove (see [START_REF] Lasser | Computing quantum dynamics in the semiclassical regime[END_REF]Lemma 5.20]) that Wig(g ε X , g

Γ z ,ε Φ z )(w) = (πε) -d γ X,z exp i ε J(X -Φ z ) • w + i 2ε G z (t)(w -m X,z ) • (w -m X,z ) ,
where m X,z = 1 2 (X +Φ z ) is the mean of the two centres, while γ X,z is a complex number with |γ X,z | 1 and G z (t) ∈ S + (2d). Repeated integration by parts, see [START_REF] Lasser | Computing quantum dynamics in the semiclassical regime[END_REF]Proposition 5.21], yields an upper bound

| g ε X , R ε z (t) | c z (t) X -Φ t,t 0 (z) √ ε -(d+1)
,

where the constant c z (t) > 0 depends on bounds of the function L z (t,t 0 , w) and is inversely proportional to the smallest eigenvalue of ℑG z (t). A combination of arguments given in the proofs of [START_REF] Lasser | Computing quantum dynamics in the semiclassical regime[END_REF]Lemma 5.18] and [START_REF] Robert | On the Herman-Kluk Semiclassical Approximation[END_REF]Lemma 3.2], reveals that the spectrum of ℑG z (t) is bounded away from zero uniformly in z, which implies the existence of constant c(t) > 0 such that

| g ε X , R ε z (t) | c(t) X -Φ t,t 0 (z) √ ε - (d+1)
.

This gives us enough decay to deduce the existence of constants C 1 (t),C 2 (t) > 0 such that the Schur estimate (5.2) holds.

Passing from thawed to frozen approximation

Here we present a slight variant of [START_REF] Robert | On the Herman-Kluk Semiclassical Approximation[END_REF]Proposition 4.1] for passing from a thawed to a frozen Gaussian approximation by a linear deformation argument.

Proposition 5.1. Let U(t,t 0 , z) be a smooth function with values in C N , N 1, whose derivatives are at most of polynomial growth. Then,

(2πε) -d R 2d g ε z , ψ U(t,t 0 , z) e i ε S(t,t 0 ,z) g Γ(t,t 0 ,z),ε Φ t,t 0 h (z) dz = (2πε) -d R 2d g ε z , ψ u 0 (t,t 0 , z) U(t,t 0 , z) e i ε S(t,t 0 ,z) g ε Φ t,t 0 h (z) dz + O(ε)
uniformly for all ψ ∈ L 2 (R d , C) with norm one.

Proof. For notational simplicity, we omit the time-dependance in S = S(t,t 0 , z), Φ = Φ t,t 0 h (z), Γ = Γ(t,t 0 , z), and U = U(t,t 0 , z). We consider both operators, the thawed and the frozen one, as special members of a class of linear operators of the form

I ψ = (2πε) -d R 2d g ε z , ψ (x -Φ q (z)) α W (z) e i ε S(z) g G (z),ε Φ(z) dz ,
that are defined by two smooth functions G : R 2d → S + (d) and W : R 2d → C N . The Siegel half-space S + (d) is invariant under inversion in the sense that any G ∈ S + (d) is invertible with -G -1 ∈ S + (d). We require that the smallest eigenvalue of ℑ(G (z)) and ℑ(-G -1 (z)) are bounded away from zero uniformly in z. The monomial powers (x -Φ q (z)) α with α ∈ N d 0 are included for technical reasons, that will become clear soon. These operators are bounded on L 2 (R d ) and satisfy

I C ε |α|/2 , (5.3)
where the constant C > 0 independent of ε, and |α|/2 denotes the smallest integer |α|/2, see [START_REF] Lasser | Computing quantum dynamics in the semiclassical regime[END_REF]Proposition 5.12]. We linearly link the thawed matrix function z → Γ(z) and the frozen z → i Id by setting

G (z, s) = (1 -s)Γ(z) + is Id ∈ S + (d), s ∈ [0, 1],
and consider the corresponding Gaussian function, that is only partially normalised,

g G (z,s),ε Φ(z) (x) = (πε) -d/4 e i ε Φ p (z)•(x-Φ q (z))+ i 2ε G (z,s)(x-Φ q (z))•(x-Φ q (z)) .
We now aim at constructing a smooth function W (z, s) with two properties.

1. Firstly, we require that W (z, 0) = det -1/2 (A(z) + iB(z)) U(z), ensuring that the deformation value s = 0 corresponds to the thawed approximation.

2. Secondly, we hope to acchieve

∂ ∂ s (2πε) -d R 2d g ε z , ψ W (z, s) e i ε S(z) g G (z,s),ε Φ(z) dz = O(ε),
uniformly for all ψ ∈ L 2 (R d ) of norm one.

Once this construction has been carried out, we will verify that the deformation yields the frozen approximation for s = 1, that is, W (z, 1) = u 0 (z) U(z). As a first step, we open the inner product involving the standard Gaussian g ε z and examine the multivariate exponential function

g ε z (y) e i ε S(z) g G (z,s),ε Φ(z)
.

Using the derivative properties of the action S(z), one obtains the identity

(x -Φ q (z))g ε z (y) e i ε S(z) g G (z,s),ε Φ(z) = ε i M -T G (z,s) (z)(i∂ q + ∂ p ) -f (x, z, s) g ε z (y) e i ε S(z) g G (z,s),ε Φ(z)
where

M G (z,s) (z) = -iG (z, s)A(z) -G (z, s)B(z) + iC(z) + D(z) is an invertible complex d × d and f (x, z, s) = M -T G (z,s) (z) 1 2 (i∂ q j + ∂ p j )G (z, s) (x -Φ q (z)) • (x -Φ q (z)) d j=1
is a vector-valued function, that is quadratic in x-Φ q (z). Since i 2ε (x-Φ q (z))• f (x, z, s) is cubic in x -Φ q (z), we use (5.3) and recognize its contribution as a term of order ε. Thus, we have

(2πε) -d R 2d g ε z , ψ W (z, s) e i ε S(z) ∂ s g G (z,s),ε Φ(z) dz = (2πε) -d R 2d W (z, s) L(z, s)(x -Φ q (z)) • (i∂ q + ∂ p ) g ε z , ψ e i ε S(z) g G (z,s),ε Φ(z) dz + O(ε), with L(z, s) = 1 2 M -1 G (z,s) (z)∂ s G (z, s).
We perform an integration by parts and arrive at

-(2πε) -d R 2d d ∑ k=1 (i∂ q k + ∂ p k ) W (z, s) (L(z, s)(x -Φ q (z))) k × g ε z , ψ e i ε S(z) g G (z,s),ε Φ(z) dz + O(ε).
Computing the derivative we obtain several terms that are linear in x -Φ q (z), and thus of order ε. The contributions we have to keep are

-W (z, s) d ∑ k, =1 L(z, s) k (i∂ q k + ∂ p k )(x -Φ q (z)) = W (z, s) d ∑ k, =1 L(z, s) k (iA(z) + B(z)) k = W (z, s) tr(L(z, s)(iA(z) + B(z))) .
We observe that

L(z, s)(iA(z) + B(z)) = -1 2 M -1 G (z,s) (z)∂ s M G (z,s) (z).
This suggests that W (z, s) should solve the differential equation

∂ s W (z, s) -1 2 tr M -1 G (z,s) (z)∂ s M G (z,s) (z) W (z, s) = 0, that is, W (z, s) = 2 -d/2 det 1/2 (M G (z,s) (z)) U(z),
by using Liouville's formula for the differentiation of determinants. Checking for the initial condition at s = 0, we observe that

M G (z,0) (z) = -i (Γ(z)(A(z) -iB(z)) -(C(z) -iD(z))) = -i Γ(z) -Γ(z) (A(z) -iB(z)) = 2 ℑΓ(z) (A(z) -iB(z)) = 2 (A(z) + iB(z)) -T
which implies W (z, 0) = det -1/2 (A(z) + iB(z)) U(z), indeed.

Herman-Kluk approximation for the toy model

The Schrödinger equation associated with (1.2) writes as a transport equation and, integrating along curves s → (s, x + s) it reads

iε d ds ψ ε (s, x + s) = k(x + s)V θ (x + s)ψ ε (s, x + s) (5.4) with V θ (x) = 0 e iθ x e -iθ x 0 .
The equation reduces to the system of ODEs

iε d dσ η ε (σ ) = kσV θ (σ )η ε (σ ). (5.5)
This problem was solved first in [START_REF] Friedrichs | On the adiabatic Theorem in Quantum Theory[END_REF] then in a more general setting in [START_REF] Hagedorn | Adiabatic Expansions near Eigenvalue Crossings[END_REF] where an asymptotic expansion in power of ε 1/2 (with power of log ε corrections), at any order, is established for R ε k,θ (σ , σ 0 ), the propagator (or resolvent matrix) of the linear differential equation (5.5); we shall use this result below . It is then possible to prove the Herman-Kluk approximation of the conjecture (4.1). Proposition 5.2. Consider t 0 < 0 and an initial data of the form

ψ ε 0 = V -v ε 0 + O(ε) in L 2 (R d ), (5.6)
where V -is the smooth eigenvector for the minus mode (see (1.3)). Set

I ε sc (v ε 0 ) =(2πε) -d z∈R 2d v ε 0 , g ε z e i ε S -(t,t 0 ,q) V -(t,t 0 , q + t -t 0 )g ε Φ t,t 0 -(z) dz + √ ε κ -(2πε) -d z∈R 2d 1 q>0 1 t>t (z) v ε 0 , g ε z e i ε (S -(t 0 ,t (z),q)+S + (t,t (z),0)) × V + (t,t (z),t -t (z))g ε,1-ik Φ t,t (z) + (0,p (z)) dz (5.7)
where κ := 2π ik θ 2 and for z = (q, p), t (z) = t 0q, and p (z) = pkq have been

computed in Example 4.3. Then, for all χ ∈ L 1 (R), in L 2 (R d ) R χ(t) I ε sc (v ε 0 ) -U ε H k,θ (t,t 0 )ψ ε 0 (x) dt = o( √ ε).
The scalar Herman-Kluk prefactor is 1 because the width of the Gaussians wave packet stays constant along the propagation (see Ex. 2.3). Note also (see Ex. 4.3) T g 1 (y) = e -i 2ε ky 2 g 1 = g 1,(1-ik/2) ; in the statement above, we have chosen not to froze the Gaussian after crossing time. The value of the coefficient κ arises from Theorem 4.2 and Ex. 4.3.

The proof relies on the analysis of the propagator R ε k,θ (σ , σ 0 ) as performed in [START_REF] Hagedorn | Adiabatic Expansions near Eigenvalue Crossings[END_REF] (see p. 280 therein). We denote by Y ± (σ , σ 0 ) the time-dependent eigenprovectors of V θ (σ ) for the eigenvalues E ± (σ ) = ±kσ (see Remark 3.2): Y ± (σ , σ 0 ) = V ± (σ , σ 0 , σ ) = V ± (σ + r, σ 0 + r, σ ), ∀r ∈ R .

(5.8)

Then, the solutions η(σ ) with initial data at time σ 0 < 0 of the form η(σ 0 ) = η 0 Y -(σ 0 ): We observe that, with the notations of Ex. 4.

Similarly, using S + (t,t (z), q) = -k 2 (tt 0 + q) 2 , S -(t ,t 0 , q) = 3k 2 q 2 and Φ t,t (z) + (0, p (z)) = (tt 0 + q, p -2kqk(tt 0 )), , we obtain V -(t,t 0 , q -(t))g ε Φ t,t 0

g ε z (x -t + t 0 ) = e i ε k(x-t+t 0 )(
-(z) dz + O(ε)
in L 2 (R d ), which allows to identify the first term of (5.7) since q -(t) = q + tt 0 .

In the second term, one can treat similarly the term V + (x, 0, x) that turns into (using also (5.8))

V (tt 0 + q, 0,tt 0 + q) = V (t,t 0q,tt 0 + q) = V (t,t (z),tt (z)

It remains to turn the discontinuous function x → 1 0<x<t-t 0 into 1 0<q + (t)<t-t 0 = 1 q<0 1 t>t 0 -q = 1 q<0 1 t>t (z) .

One then takes advantage of the averaging in time to use the estimate of 5.3 despite the discontinuity. One regularizes the discontinuous function which will differ from its regularization on a set of small Lebesgue measure in the variable t. One can then use the preceding argument on the regularized term and gets rid of the correction ones by estimating them thanks to the estimate (5.3) and using the smallness of integrals on χ on sets of small Lebesgue measures.

Example 4 . 1 .

 41 For the toy model (1.2), u(x) = (0, cos(θ x), sin(θ x)), v(ξ ) = ξ and f (ξ ) = kξ . Hence we have Assumption (SC) : ϒ = {x = 0} and {v, f } = k = 0.

1 .e i ε k σ σ 0

 10 If σ < 0, then η(σ ) = τdτ Y -(σ , σ 0 )η 0 + O(ε).

  Y -(σ , σ 0 )η 0 -√ ε(1i) √ π (-k) -1/2 h(0)e i ε λ e -i ε k σ 0 τdτ Y + (σ , 0)η 0 + O(ε) with h(0) = Y + (σ , 0), d dσ Y -(σ , 0) σ =0= -iθ /2 and λ = k 0 σ 0 τdτ.

  Proof of Proposition 5.2 . By(5.4),ψ ε (t, x) := U ε H k,θ (t,t 0 )ψ ε 0 (x) = R ε k,θ (x, xt + t 0 )ψ ε 0 (xt + t 0 ),and, in view of Friedrichs' description, we deduceψ ε (t, x) =e ik ε x(t-t 0 )-ik 2ε (t-t 0 ) 2 V -(t,t 0 , x)v ε 0 (xt + t 0 ) + √ εκ 1 0<x<t-t 0 e ik ε x(t-t 0 )-ik 2ε (t-t 0 ) 2 -i ε x 2 v ε 0 (xt + t 0 ) V + (x, 0, x) + o(

	3,
	-(1 -i) √ π (-k) -1/2 h(0) =	2π ik	θ 2	=	2π ik	γ = κ.
						√
						ε)
	where we have used
	Y z∈R 2d	v ε 0 , g ε z g ε

-(x, xt + t 0 ) = V -(t,t 0 , x), Y + (σ , 0) = V + (x, 0, x) and λ = -k 2 (xt + t 0 ) 2 .

Using (1.8), we write

v ε 0 (xt + t 0 ) = (2πε) -d z (xt + t 0 )dz

and we observe that, in view of S -(t,t 0 , z) = kq(tt 0 ) + k 2 (tt 0 ) 2 , we have

g ε z (xt + t 0 ) = e -ik ε x(t-t 0 )+ i 2ε (t-t 0 ) 2 e i ε S -(t,t 0 ,z) g ε Φ t,t 0 -(z)

  2q+(t-t 0 ) g ε -(t ,t 0 ,q)+ i ε S + (t,t ,0)+ k iε (x-q-t+t 0 ) 2Using Taylor expansion and Lemma 5.3, we can transform the first part of the righthand side of the preceding equation:

		(2πε) -d	z∈R 2d	v ε 0 , g ε z e	i ε S -(t,z) V -(t,t 0 , x)g ε Φ t,t 0 -(z)	dz
			= (2πε) -d	z∈R 2d	v ε 0 , g ε z e	i ε S -(t,z)
						Φ + t,t (z)	(0,p (z))	(x)
					= e	i ε S g ε Φ + t,t (z)	(0,p (z))	(x)
					= e	i ε S -(t ,t 0 ,q)+ i ε S + (t,t ,0) WP ε Φ + t,t (z)	(0,p (z))	(T g 1 )
	Putting these elements together, we are left with
	ψ ε (t, x) = o( √ ε) + (2πε) -d	z∈R 2d	v ε 0 , g ε z e	Φ t,t 0 -(z)	dz
	+	√ εκ(2πε) -d	z∈R 2d	1 0<x<t-t 0 v ε 0 , g ε z e

i ε S -(t,z) V -(t,t 0 , x)g ε i ε S -(t ,t 0 ,q)+ i ε S + (t,t ,0) V + (x, 0, x)g ε Φ t,t + (z)

dz.


Margot Selosse
email: margot.selosse@gmail.com

Isobel Claire Gormley
email: claire.gormley@ucd.ie

Julien Jacques
email: julien.jacques@univ-lyon2.fr

Christophe Biernacki
email: christophe.biernacki@inria.fr

A bumpy journey: exploring deep Gaussian mixture models

The deep Gaussian mixture model (DGMM) is a framework directly inspired by the finite mixture of factor analysers model (MFA) and the deep learning architecture composed of multiple layers. The MFA is a generative model that considers a data point as arising from a latent variable (termed the score) which is sampled from a standard multivariate Gaussian distribution and then transformed linearly. The linear transformation matrix (termed the loading matrix) is specific to a component in the finite mixture. The DGMM consists of stacking MFA layers, in the sense that the latent scores are no longer assumed to be drawn from a standard Gaussian, but rather are drawn from a mixture of factor analysers model. Thus the latent scores are at one point considered to be the input of an MFA and also to have latent scores themselves. The latent scores of the DGMM's last layer only are considered to be drawn from a standard multivariate Gaussian distribution. In recent years, the DGMM has gained prominence in the literature: intuitively, this model should be able to capture complex distributions more precisely than a simple Gaussian mixture model. We show in this work that while the DGMM is an original and novel idea, in certain cases it is challenging to infer its parameters. In addition, we give some insights to the probable reasons of this difficulty. Experimental results are provided on github: https://github.com/ansubmissions/ICBINB, alongside an R package that implements the algorithm and a number of ready-to-run R scripts.

Mixture of Factor Analysers

 is an extension of the factor analysis model to a finite mixture of G factor analysis models. We assume that an observation x i ∈ R J can be expressed through a latent random vector z i ∈ R R such that R ≤ J. In addition, the latent labels {v 1 , . . . , v n } are assumed to be independent unobserved realisations of a categorical G-dimensional random vector,

such that v ig = 1 indicates that x i is generated by the gth factor analyser. The generative process of this model is as follows:

v i ∼ M(1, (π 1 , . . . , π G )), z i ∼ N (0, I R ), u i ∼ N (0, ψ g ), x i = η g + Λ g z i + u i ,
where g is the index of the element of v i that is equal to 1, Λ g is a J × R matrix and η g ∈ R J is the mean vector of the gth factor analysis model. In addition, ψ g is a diagonal covariance matrix of dimension J, and π = (π 1 , . . . , π G ) are the mixing proportions. The latent variable v i = (v i1 , . . . , v iG ) represents the partitions for observation i and p(v ig = 1) = π g . The marginal density of x i is then a Gaussian mixture model (GMM):

f (x i ; θ) = G g=1 π g f g (x i ; η g , Σ g ),
where θ = (η g , Λ g , ψ g , π g ) g∈{1,...,G} denotes all parameters and:

• f g (.; η g , Σ g ) is the Gaussian probability density function for component g with mean η g and covariance matrix Σ g and

• Σ g = Λ g Λ T g + ψ g .
Hence, the observed-data log-likelihood of N samples x = (x 1 , . . . , x N ) is given by:

l o (θ; x) = N i=1 log G g=1 π g f g (x i ; η g , Σ g ).
2 Deep Gaussian Mixture Models

Coupling Deep Learning and Gaussian Mixture Models

One of the earliest papers to define a model that stacks MFA layers in a deep learning fashion is that of [START_REF] Tang | Deep mixtures of factor analysers[END_REF] where the authors refer to the model as the "Deep Mixture of Factor Analysers". In this model, the layers are not fully-connected and the architecture is a tree: a neuron receives an input from only one neuron of the previous layer. In [START_REF] Van Den Oord | Factoring variations in natural images with deep Gaussian mixture models[END_REF], the authors define a DGMM model with fully-connected layers, but the matrices of weights Λ are assumed to be square. In other words, there is no dimension reduction between each MFA layer, which leads to important unidentifiability issues. In the model proposed by [START_REF] Yang | Deep mixtures of factor analyzers with common loadings: A novel deep generative approach to clustering[END_REF], instead of using the MFA model as a layer, the authors suggest using a mixture of factor analysers model with common factor loadings (MCFA) [START_REF] Baek | Mixtures of factor analyzers with common factor loadings: Applications to the clustering and visualization of high-dimensional data[END_REF]: the matrices Λ g of the same layer are therefore assumed to be equal for all g. Finally, [START_REF] Viroli | Deep Gaussian mixture models[END_REF] defines the deep Gaussian mixture model (DGMM) with dimension reduction at each layer and provides an EM algorithm for the inference of the parameters. This present work aims at investigating the statistical efficiency of this algorithm.

Definition of the Deep Gaussian Mixture Model

This section relies on [START_REF] Viroli | Deep Gaussian mixture models[END_REF]. Let us assume that there are L layers and that the data x = (x i ) i∈{1,...,N } are generated as follows:

x i = z (0) i = η (1) g1 + Λ (1) g1 z (1) i + u (1)
i with prob. π (1) g1 , g 1 ∈ {1, . . . , G (1) } z

(1) i = η (2) g2 + Λ (2) g2 z (2) i + u (2) i with prob. π (2) g2 , g 2 ∈ {1, . . . , G (2) } . . . z (l) i = η (l+1) g l+1 + Λ (l+1) g l+1 z (l+1) i + u (l+1) i with prob. π (l+1) g l+1 , g l+1 ∈ {1, . . . , G (l+1) } (1)
. . .

z (L-1) i = η (L) g L + Λ (L) g L z (L) i + u (L) i
with prob. π (L) g L , g L ∈ {1, . . . , G (L) }, where:

• z (L) i
is a vector of size R (L) and is assumed to be drawn from the Gaussian distribution

N (0, I R (L) ), • u (l)
i is a specific random error that follows a Gaussian distribution with expectation 0 and covariance matrix ψ (l)

g l , • η (l) g l is a vector of length R (l) , • Λ (l) g l is a matrix of dimension R (l-1) × R (l) , • J > . . . > R (l-1) > R (l) for all l.
The vectors u (l) i are assumed to be independent of the scores z (l) i . From these considerations, we see that at each layer l, the conditional distribution of z

(l) i given z (l+1) i
is a mixture of Gaussian distributions such that:

f (z (l) i |z (l+1) i ; θ) = G (l+1) g l+1 =1 π g l+1 N (η (l+1) g l+1 + Λ (l+1) g l+1 z (l+1) i
, ψ (l+1) g l+1 ).

(

) 2
Figure 1 represents the graphical model of this DGMM in the case of L = 2. Therefore, the DGMM has parameters θ = (η (l)

g l , Λ (l)
g l , ψ

(l)

g l , π (l)
g l ) l∈{1,...,L};g l ∈{1,...,G (l) } . It also has the latent variables v (l) and z (l) where v (l) corresponds to the partition of the lth layer.

In addition, it is important to notice that the marginal probability distribution f (x i ; θ) is a GMM where each component corresponds to one of the possible paths of the network. Therefore, by noting G as the set of all possible paths through the network, the DGMM can be written:

f (x i ; θ) = g∈G π g N (µ g , Σ g ), (3)
where g = (g 1 , . . . , g L ) is one of the possible paths of the network,

π g = L l=1 π (l) g l , (4)
µ g = η (1) g1 + Λ (1) g1 η (2) g2 + Λ (2) g2 . . . (η (L-1) g L-1 + Λ (L-1) g L-1 η (L) g L ) and (5)
Σ g = Λ (1) g1 Λ (2) g2 . . . (Λ (L) g L Λ (L) g L T + ψ (L) g L ) . . . Λ (2) g2 T + ψ (2) g2 Λ (1) g1 T + ψ (1) g1 . (6)
We refer to the GMM of Equation (3) as the "global GMM" of the DGMM so as not to confuse it with the GMMs of the layers as expressed in Equation ( 2). In addition, the number of possible paths through the network is denoted as G. In fact, Equation ( 3) can be generalised to all layers: the marginal distributions of all the scores z (l)

i are Gaussian mixtures. So, by noting G(l) as all the possible paths from the lth layer, and by integrating out the latent variables of the layers that follow the lth layer, we have:

f (z (l) i ; θ) = g(l) ∈ G(l) πg (l) N ( μ(l+1) g(l) , Σ(l+1) g(l) ), (7)
where g(l) = (g l+1 , . . . , g L ) is a possible path of the network from the lth layer,

πg (l) = L l =l+1 π g l , (8)
μg (l) = η (l+1) g l+1 + Λ (l+1) g l+1 η (l+2) g l+2 + Λ (l+2) g l+2 . . . (η (L-1) g L-1 + Λ (L-1) g L-1 η (L)
g L ) and ( 9)

Σg (l) = Λ (l+1) g l+1 Λ (l+2) g l+2 . . . (Λ (L) g L Λ (L) g L T +ψ (L) g L ) . . . Λ (l+2) g l+2 T +ψ (l+2) g l+2 Λ (l+1) g l+1 T +ψ (l+1) g l+1 . (10)

Inference of the model

Similar to MFA, the DGMM has parameters and latent variables to estimate. Therefore, the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is a candidate to optimise the log-likelihood function with respect to the parameters and latent variables. However, in the MFA model, the computation of the expectation of the complete data log-likelihood does not involve p(z; θ) since z is assumed to be drawn from the standard Gaussian distribution N (0, I). For the DGMM, unless l = L, z

l is assumed to be drawn from an MFA with parameters (π (l+1)

g l+1 , η (l+1) g l+1 , Λ (l+1) g l+1 , ψ (l+1)
g l+1 ) g l+1 ∈{1,...,G (l+1) } . Thus, the parameters of the model are involved in the computation of p(z (l) ; θ) which makes it difficult to optimise the auxiliary function with all the layers at the same time. A solution to this issue is to perform the optimisation layer by layer in iterations of the EM algorithm. In this case, when we optimise the auxiliary function with respect to the parameters of layer l, the fact that p(z (l) ; θ) depends on the parameters of layer l + 1 is not a problem and an EM algorithm can be performed. We refer to [START_REF] Viroli | Deep Gaussian mixture models[END_REF] for further details about the EM algorithm steps. Finally, choosing the number of clusters for each layer can be done with a model selection criterion such as the BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF].

Interpreting the inferred partitions under the DGMM

There are two main ways to use the inferred partitions of the DGMM. The first way focuses on the first layer of the network, and then have a look at the partitions v (1) . In this case, the number of clusters is equal to G (1) and the distribution of the observations of a cluster is assumed not to be Gaussian since the latent scores z (1) are assumed to be drawn through an MFA, as detailed in Equation ( 2). The second way of using the DGMM is to observe the clusters of the global GMM of Equation [START_REF] Peel | Modelling data by mixtures of factor analyzers[END_REF]. In this case, the number of clusters is equal to G (1) × . . . × G (L) , and the observations of a cluster are now assumed to follow a Gaussian distribution.

Experiments with the Deep Gaussian Mixture Model

We simulate data via the generative process of the DGMM with known parameters and partitions. Then, we run the EM algorithm to see if it estimates the parameters and partitions well. Because the number of parameters is large, we focus on the Adjust Rand Index (ARI), which gives an indicator of the performances of the algorithm on the estimation of the partitions and hence on the estimation of the parameters too. In these experiments we start by focusing on the partitions of the first layer, and then we focus on the partitions of the global GMM.

Settings for the simulated data set

For the data set x, we have N = 4, 000 and J = 17. The network was built with three layers (L = 3) and G (1) = 3, G (2) = 3 and G (3) = 2. In addition, the latent scores dimensions were set to R (1) = 11, R (2) = 6 and R (2) = 2 respectively; note this DGMM is identified in a factorial sense.. The parameters π (l) g were set to 1/G (l) for all g and for all l. The parameters (η

(l) g , Λ (l) g , ψ (l)
g ) g,l were sampled from Gaussian distributions whose parameters can be found in the supplementary material files.

Results on the first layer

In this section, we focus on the partitions of the first layer v (1) . We run 20 times the DGMM EM algorithm and the classic GMM by using the R package mclust [START_REF] Scrucca | mclust 5: clustering, classification and density estimation using Gaussian finite mixture models[END_REF]. For the DGMM runs, we fix G (2) and G (3) to their ground-truth values and vary G (1) from 2 to 6. In addition, the algorithm is run with random starting values. For the GMM runs, we vary the number of components G from 2 to 6. Table 1 shows for each run the number of clusters with the best BIC value. We see that mclust was not able to find the true number of clusters. This result is expected since the clusters of the first layer are not supposed to be Gaussian: mclust overfits the data by finding more clusters than there truly are. However, the DGMM seems to yield the same behaviour, which is not expected since its first layer is supposed to capture non-Gaussian distributions. DGMM 6 6 6 5 6 4 4 5 5 6 5 6 5 6 5 6 6 6 6 5 mclust 6 6 6 6 6 6 6 6 6 6 5 6 6 6 6 6 6 5 5 6

Table 1: Number of clusters chosen with BIC for each run of the DGMM and mclust.

Figure 2 represents the ARI of the runs with the best BIC value. We see that the DGMM yields better results than mclust. However, these results are not satisfying because the data was generated through the DGMM generative process: one would hope that at least some of the runs should achieve an ARI equal to 1. Let us note that similar behaviour was observed across a large number of data sets generated with a range of settings of the DGMM architecture but for brevity only one such is presented here.

Results on the global GMM

For every experiment described below, the number of clusters (G (l) ) l are fixed to their true value. The DGMM's inference algorithm is run 20 times with random starting values.

Figure 3 shows the resulting ARI. The DGMM yields poor results since the highest ARI value achieved is 0.67 even though the data set was generated through the DGMM's generative process.

In Figure 4, we plot the BIC values over the iterations of the EM algorithm for a single random initialisation; the BIC reaches a plateau, meaning that the EM algorithm achieves at least a local maximum of the likelihood. If the EM algorithm is able to find local maxima, but fails to find the right partitions, it might mean that there are many local maxima. To verify this, we run the EM algorithm a hundred times and compute the ARI of each resulting clustering with respect to every other resulting clustering. The mutual ARI are plotted in Figure 5 (left and right), where we see that none of the mutual ARI values are 1 (except the ARI of a result with respect to itself). This means that the EM algorithm found 100 different local maxima leading to very different partitions.

The last result we detail regards another unexpected behaviour of the DGMM's inference algorithm.

When we list the resulting partition for each observation on each layer, we can see that all the clusters of every layer are represented. However, when we list the resulting partitions of the global GMM as in Table 2, we see that not all the clusters are represented. This means that even if each neuron of each layer has observations, not every combination of neurons of the three layers does. Therefore, there is a severe problem of empty clusters at the global GMM's level.

Finally, let us note that the package mclust was not able to find a solution for G = l G (l) = 18. From this result, the DGMM can be seen as a candidate when the number of clusters is high, even if the estimation of partitions is not perfect. We see that the latent scores (z (l) ) l do not seem to be estimated well.

Possible reasons for the poor results

There are several possible reasons for the poor results described in the previous section. First of all, the MFA model already suffers from unindentifiability issues that could worsen with the multilayer architecture of the DGMM. The most common solutions to tackle the unidentifiablity issues are:

• To constrain Λ (l)
g so that it is orthonormal and to order the columns by decreasing variance of the corresponding latent scores.

• To constrain Λ (l) g T ψ (l) g Λ (l)
g to be diagonal.

• To constrain Λ (l)
g so that it is a lower triangular full rank matrix with diagonal elements strictly positive.

• To choose an informative rotation matrix. Many heuristic methods try to find rotation matrices P that can be used to modify Λ (l) g (and the latent factors) so as to try to increase the interpretability.

While we already tried to apply these constraints (results not detailed in this document), we did not obtained better results regarding the partitions and parameters estimation.

Second, the DGMM has many latent variables, which can be a severe problem when it comes to estimation. Indeed, in the context of the configuration of Section 3.1, the latent scores alone represent N × (R (1) + R (2) + R (3) ) = 4, 000 × (11 + 6 + 2) = 76, 000 continuous values to estimate, which is very high. This large number of latent variables is for us one of the most important obstacles to correct estimation. To support this point, Table 3 shows the ARI that the EM algorithm achieved at each layer. At the first layer, the ARIs are globally higher than the other layers and can even be equal to 1. However, the ARI for the other layers worsens with the depth: the second layer's higher ARI value is 0.51 and the third layer's ARI value are all near to zero. This probably means that the latent scores are poorly estimated since the DGMM of layer 1 depends on the data x, the DGMM of layer 2 depends on the scores z (1) and the DGMM of layer 3 depends on the scores z (2) . A solution to bypass the latent score estimation is to use another kind of algorithm, such as the Newton-Raphson algorithm and its variants. However, we encountered different problems when using this kind of algorithm such as computational errors due to matrix inversions. Finally, we have also observed that even though the clusters of each layer do not necessarily become empty, the clusters of the global GMM do so, which is also a problem for estimation of the partition.

Conclusion

In theory, the DGMM is a powerful model that can estimate complex distributions by stacking MFA layers in a deep learning fashion. We have seen in this Section 3.3 that this algorithm is able to provide an estimation of partitions and parameters in cases where mclust is not able to do so, for instance when the number of clusters is high. However, this estimation is not correct and we have seen that the EM algorithm suffers from drawbacks and can struggle to find the global maxima.

Figure 1 :

 1 Figure 1: Graphical model of an example of the DGMM.

Figure 2 :

 2 Figure 2: ARI of the DGMM's first layer and of the mclust result.

Figure 3 :

 3 Figure 3: Boxplot of the ARI for the global GMM of the DGMM using random initialisations.

Figure 4 :

 4 Figure 4: BIC values over the iterations of the EM algorithm using a single random initialisation.

Figure 5 :

 5 Figure 5: Mutual ARI of 100 runs of the EM algorithm. None of them gets the same final partition.

Table 2 :

 2 The non-empty clusters that are represented in the global GMM for the 20 runs.

	Cluster numbers 1, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 18 1, 2, 3, 5, 6, 7, 9, 10, 11, 12, 14, 15, 16 1, 3, 4, 5, 6, 7, 8, 9, 10, 13, 14, 15, 16, 17	ARI of layer 1 ARI of layer 2 ARI of layer 3 0.18 0.5 0 0.17 0.17 0.01 0.2 0.25 0.01
	2, 3, 4, 5, 6, 7, 11, 12, 13, 14, 15, 16, 18	0.45	0.19	0
	1, 2, 3, 5, 6, 7, 9, 10, 12, 14, 15, 16, 18	0.32	0.29	0.01
	1, 2, 4, 6, 7, 9, 10, 11, 13, 14, 15, 16, 17, 18	0.38	0.09	0
	1, 2, 5, 6, 8, 10, 11, 12, 14, 15, 17 2, 3, 4, 5, 6, 9, 11, 12, 13, 14, 16, 18 2, 3, 4, 7, 9, 10, 11, 12, 13, 15, 16, 17, 18 1, 2, 3, 4, 5, 8, 9, 10, 11, 12, 13, 17, 18 1, 2, 3, 5, 6, 7, 9, 10, 11, 14, 15, 16, 18	0.43 0.4 0.55 0.11 0.17 0.17	0.07 0.22 0.17 0.24 0.37 0.19	0 0 0.01 0.01 0 0
	2, 3, 5, 6, 7, 8, 9, 11, 12, 14, 15, 16	1	0.08	0.01
	1, 2, 4, 5, 6, 7, 8, 10, 11, 13, 14, 15, 16, 17, 18	0.19	0.42	0.01
	1, 2, 3, 4, 5, 8, 9, 10, 11, 12, 14, 17, 18	0.5	0.51	0
	1, 2, 4, 5, 8, 9, 10, 11, 13, 14, 17, 18 1, 2, 3, 4, 5, 6, 7, 10, 11, 12, 13, 14, 15, 17 1, 2, 3, 4, 6, 7, 9, 10, 11, 12, 15, 16, 17 1, 4, 5, 8, 10, 13, 14, 15, 17, 18 1, 3, 4, 7, 8, 9, 10, 11, 12, 13, 15, 17, 18	0.22 0.4 0.29 0.55 0.52	0.51 0.21 0.09 0.07 0.21	0 0 0 0 0
	1, 5, 7, 8, 9, 10, 14, 16, 17, 18

Table 3 :

 3 ARI results for the GMM of each layer.

Acknowledgments and Disclosure of Funding

We would like to warmly thank Cinzia Viroli and Geoff McLachlan for their helpful discussions from which the work benefitted.


Couro Kane
email: fcouro.kane@polytechnique.fr

Fran Cois Jouve
email: francois.jouve@polytechnique.fr

Marc Schoenauer
email: marc.schoenauerg@polytechnique.fr

STRUCTURAL TOPOLOGY OPTIMIZATION IN LINEAR AND NONLINEAR ELASTICITY USING GENETIC ALGORITHMS

In this paper, structural topology optimization is addressed through Genetic Algorithms. A set of designs is evolved following the Darwinian survival-of-ttest principle. The standard crossover and mutation operators are tailored for the needs of 2D topology optimization. The genetic algorithm based on these operators is experimented on plane stress problems of cantilever plates: the goal is to optimize the weight of the structure under displacement constraints.

The main advantage of this approach is that it can both nd out alternative optimal solutions, as experimentally demonstrated on a problem with multiple solutions, and handle di erent kinds of mechanical model: some results in elasticity with large displacements are presented. In that case, the nonlinear geometrical e ects of the model lead to non viable solutions, unless some constraints are imposed on the stress eld.

INTRODUCTION

Since the seminal work of [START_REF] Holland | Adaptation in natural and articial systems[END_REF] and the comprehensive study of [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF], Genetic Algorithms (GAs) have gradually been recognized as powerful stochastic optimization algorithms. Their strength proceeds from their wide range of applications: GAs can handle non derivable, non continuous and even non analytically de ned functions 1 .

The main interest of stochastic methods in Engineering Sciences is to break the limits of standard deterministic methods in many optimization problems: when the search space involve both discrete and continuous domains (e.g. for the optimal design of truss structures); when the objective function or the constraints lack regularity; or when the objective function admits a huge number of local optima. In counterpart, stochastic methods are computationally expensive: GAs for instance, are slower than classical 1 e.g. functions only available through computation.

optimization methods by about one or two orders of magnitude (when comparison is possible, i.e. when classical methods apply).

This paper focuses on applying GAs to some wellstudied problems in mechanical engineering, namely the structural topology optimization of cantilever plates. Section 2 brie y reviews some related works in topology optimization; their advantages and their limitations regarding the speci cities of constrained optimization in mechanical engineering are discussed. Section 3 gives the broad lines of GAs; it then describes genetic operators tailored for topology optimization. Experimentations in the linear elasticity case are presented in section 4: The rst problem deals with minimizing the weight of a cantilever plate under displacement constraints. Further, the standard cantilever plate is modi ed to highlight the ability of the GA approach to handle problems having multiple optimal or quasi-optimal solutions. A compliance-based optimization is nally achieved; it permits to compare GA-based and homogenizationbased approaches. Section 5 presents the rst results obtained in linear elasticity with large displacements for the cantilever plate problem. The nonlinear geometrical e ects clearly show the need to take into account the stress eld in the tness function to obtain realistic solutions.

RELATED WORKS

The main trends in structural optimization can be sketched as follows.

A rst approach consists in continuously varying the domain at hand. This approach is e ective as far as a good solution can actually be obtained by a continuous transformation of the initial domain. In contrast, it does not allow to nd an optimal shape that includes holes (except if the number of holes is known beforehand, of course).

An approach to topology design, introduced by [START_REF] Bendsoe | Generating Optimal Topologies in Structural Design Using a Homogenization Method[END_REF], is that of homogeniza-tion; it consists in dealing with a continuous density of material. In the end of this deterministic optimization, the current density is forced toward value 1 or 0, that respectively stands for material present or absent. However, this approach requires the design of the homogenized operator, as thoroughly described in [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructures[END_REF], and is insofar limited to the linear elasticity case. Moreover, it cannot address loadings that apply on the actual boundary of the shape to be determined (e.g. pressure load), and hardly handles optimization for multiple loadings.

Another approach to topology design is that of stochastic optimization, such as involved in simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]) and genetic algorithms. Both methods have been applied to structural optimization: in the framework of discrete truss structures, for cross-section sizing by [START_REF] Goldberg | Engineering Optimization via Genetic Algorithms[END_REF], [START_REF] Hajela | Genetic Algorithms in Automated Structural Synthesis[END_REF], [START_REF] Lin | Genetic Search Strategies in Large Scale Optimization[END_REF], [START_REF] Schoenauer | Discrete optimal design of structures by genetic algorithms[END_REF] among others, as well as for topological optimization by [START_REF] Hajela | Genetic Algorithms in Structural Topology Optimization[END_REF], [START_REF] Grierson | Discrete Optimal Design using a Genetic Algorithm[END_REF], Wu (1995); for the optimization of composite materials by [START_REF] Leriche | Optimization of laminate stacking sequence for Buckling Load Maximization by Genetic Algorithms[END_REF]; and more recently for structural components optimization by [START_REF] Anagnostou | A computational procedure for part design[END_REF], [START_REF] Ghaddar | Analysis of a Part Design Procedure[END_REF] using simulated annealing, and by [START_REF] Jensen | Topological Structural Design using Genetic Algorithms[END_REF], [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF] with genetic algorithms.

This paper continues and extends the GA-based approaches of [START_REF] Jensen | Topological Structural Design using Genetic Algorithms[END_REF] and [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF] in several respects. First, speci c genetic operators (crossover, mutation) are introduced. Second, a problem exhibiting multiple solutions is studied, demonstrating the ability of the GA approach to successfully nd out di erent quasi-optimal solutions. Last but not least, this method is not limited to linear elasticity, as demonstrated by considering the large displacements model.

GA: RECALLS AND ADAPTATIONS

This section gives the broad lines of basic GAs; then it describes genetic operators devised for the needs of topology optimization.

Principle

Given a search space E and a tness function F dened from E onto IR + , GAs evolve a set of p individuals (points of E), termed population. This evolution crudely mimics the Darwinian evolution: according to the Darwinian survival-of-the-ttest principle, the ttest individuals, i.e. the near-optimal points of tness function F will appear in population P i for some i.

The basic step in GAs, called generation, is the transformation used for population P i to give birth to population P i+1 . This transformation involves three steps :

Selection builds population P 0 i by copying elements of P i ; the number of copies of an individ-ual increases with its tness, the total number of elements in P 0 i being same as in P i . Crossover applies on population P 0 i to build population P" i . From two individuals x and y in P 0 i , crossover builds two o springs x 0 and y 0 with probability p c (p c usually varies from .2 to 1.). When considering a bitstring representation (E = f0; 1g N ), a crossover c can be represented 2 as a bitstring itself, c = (c 1 ; : : : c N ):

x 1 :: x N y 1 :: y N ! x 1 0 :: x N 0 y 1 0 :: y N 0 with x i 0 = x i if c i = 1 y i if c i = 0 and y i 0 = y i if c i = 1

x i if c i = 0

Mutation applies on population P" i to build population P i+1 . Mutation transforms an individual x in P" i into an o spring x 0 ; when considering a bitstring representation, a mutation can similarly be represented by a bitstring m = (m 1 ; : : : m N ):

x 1 :: x N ! x 1 0 :: x N 0 with x i 0 = 1 x i if m i = 1

x i if m i = 0

The probability for m i to take value 1 is noted p m (p m usually varies from 102 to 10 4 ).

3.2 Adaptation for structural topology optimization Representation. We restrict ourselves to considering 2D domains (see sections 4 and5). An individual is here a shape embedded into the problem domain. This domain is discretized according to a grid of quadrangular elements. An individual can then be represented by an array of bits whose dimension is the number of such elements: an element is void (respectively material) whether the corresponding bit takes value 0 (resp. 1).

Note that such arrays do not all stand for viable shapes: e.g. a sound shape must be connected to the load point(s) involved in the mechanical problem. Non-viable shapes are handled through the tness calculation (they will be given a null tness; see 3.3).

Further, the usual one-dimensional bitstring representation is not faithful with respect to topology: On the one hand, adjacent bits in the bitstring representation do not necessarily correspond to neighbor elements of the domain; And on the other hand, vertical adjacent positions are quite far from one another in the bitstring representation. This remark leads to design new genetic operators, in order to take into account the two-dimensional topology.

Crossover. The classical 1-point, 2-points and uniform crossovers were the rst choice for crossover operators. As expected, both 1-point and 2-points crossovers gave poor results. And they were constantly outperformed by the uniform crossover: This can be explained by the fact that uniform crossover involves a lesser bias than 1-point and 2-point crossovers on the two-dimensional representation of 2D shapes.

Two di erent crossover operators have been purposely devised for a better transmission of topological properties between parents and o spring:

the diagonal crossover (Figure 1) involves (a) the random selection of line D (position and slope) and

(b) the exchange of the lower parts of the structures at hand; the block crossover (Figure 2), rst introduced by Jensen (1992), involves here (a) the random selection of points X 1 ; X 2 ; Y 1 ; Y 2 , de ning 9 blocks in both structures, (b) the random choice of the number of blocks to exchange (2 or 3), (c) the random selection of the blocks among the 9 blocks and (c) the actual exchange of the blocks.

1: The diagonal crossover Y Y 1 2 Y Y 1 2 X X 1 2 X X 1 2
Figure 2: The block crossover A thorough study of the in uence of the crossover operator on the standard cantilever plate problem by [START_REF] Kane | Algorithmes g en etiques et Optimisation topologique de formes[END_REF] demonstrates a slight improvement of both the diagonal operator and the block operator over the uniform operator, while both signi cantly outperform any standard 1D crossover (e.g. 1 point crossover used by [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF]). All experiments presented in this paper use the block crossover.

Mutation. Mutation classically proceeds by ipping randomly selected bits. Two other kinds of mutation are proposed, where the bits altered are chosen after either the topology of the current individual, or the whole current population. In the rst mode, termed boundary mutation, the bits ipped are preferably selected among the elements near from the boundary of the current structure (the computation of this boundary is described in section 3.3). This mutation performs slight modi cations of the individual at hand. In particular, it hardly modies the topology of the structure, and could be related to the classical "domain variation" method of shape optimization.

In the second mode, termed epistatic mutation, the goal is to counteract the genetic drift: It is well known that genetic population classically tends to become homogeneous over generations, i.e. the set of bits with constant value over the population increases. This loss of diversity is prejudicial to optimization, since it arbitrarily restricts the exploration. Epistatic mutation thus selects and ips bits whose value is almost constant over the population, in order to reintroduce the disappearing values.

The epistatic mutation performs statistically better than both the classical uniform mutation and the boundary mutation, as experimented thoroughly by [START_REF] Kane | Algorithmes g en etiques et Optimisation topologique de formes[END_REF], and all experiments presented in this paper use the epistatic mutation. However, when the population has converged, the boundary mutation should allow to re ne the solution(s) found so far. Such an iterated schema (epistatic mutation in the beginning, boundary mutation afterwards) remains to be tested.

Fitness computation

A given individual is evaluated in a 2-step process. Some seed material is imposed at the point(s) where dynamical loading is applied. The connected component containing that seed is computed. Grid elements are connected if and only if they share an edge. Note that no seed material is prescribed on the part of the boundary where the plate is xed, as it is in the approach of [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF]. The optimization process does choose where to hang the structure on the xed vertical boundary. This allows for a greater exibility in solving the optimization problem, as witnessed in the range of alternative solutions proposed for the modi ed cantilever plate problem (see section 4.3). On the other hand, more structures are likely to be disconnected from that xed boundary, leading to an ill-posed problem: such structures Once the component connecting the seed and the xed boundary has been computed (if any), the FEM analysis is performed on that component, on the same regular quadrangular mesh used to represent the individuals; the FEM tool is detailed in [START_REF] Jouve | Mod elisation math ematique de l' il en elasticit e non-lin eaire Recherches en Math ematiques Appliqu ees[END_REF]. The actual material boundary is used, in contrast with both [START_REF] Jensen | Topological Structural Design using Genetic Algorithms[END_REF] and [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF] in which the FEM analysis was done on the whole design domain, assigning a very low Young modulus to void elements. Though the results of the analysis does not di er signi cantly from one method to the other, computing the actual boundary allows to take into account loading applied on this boundary (e.g. pressure loading or heat exchanges).

Last, the tness of the structure is computed (the analytic expression of the tness is discussed in next sections), involving, in the case of elasticity: the weight of the structure; the maximum value of the displacement, or the stress eld or the compliance of the structure, computed from the above FEM analysis. Note that this approach allows as well to consider multiple loadings (with as many FEM analyses as loading cases), though no such example is given here.

LINEAR ELASTICITY

All experiments in this section consider the standard plane stress linear elastic model. The e ects of gravity are neglected.

The general framework is that of the cantilever plate: a rectangular plate is xed on the left vertical part of its boundary (both displacements are set to 0), and the loading is made of a single force applied on the middle of its right vertical boundary. Figure 4 shows the 2 1 cantilever plate.

Penalized tness

As said before, unstable structures (the connected component containing the load point does not contain any xed edge) are assigned zero tness. Following [START_REF] Jensen | Topological Structural Design using Genetic Algorithms[END_REF], the tness associated to stable structures is rst de ned by:

Fitness penalized = 1 Area + (D Max D Lim ) +
where Area is the area of the actual structure (i.e. of the connected component used in the FEM analysis), D max is the displacement of the load point (computed from the FEM analysis) and D Lim is a prescribed constraint on that displacement. Positive parameter is the penalization factor, and x + denotes the positive part of x, max(x,0).

As pointed out by [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF], the quality of the solution greatly depends on : When is large, the area is not that important and heavy structures are retained. When is small, the structures optimal according to Fitness penalized do not satisfy the displacement constraints. Nevertheless, the above tness function, for large values of , allows to precisely respect the constraint on the maximum displacement of the structure, which is not the case, neither with the compliance criterion (see [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructures[END_REF], or section 4.5) nor with the sti ness criterion used by [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF].

Iterated penalized tness

The solution proposed to overcome the critical tuning of parameter is a 2-step optimization: A medium value of is rst considered and a population of light structures is thereby determined. Then, during the evolution, the value of is increased by a factor of 10 to ensure the satisfaction of the constraints. Such ideas about iterated GAs has been demonstrated powerful in truss structure optimization, by [START_REF] Schoenauer | Constrained GA optimization[END_REF], as well as in other domains of application (see de [START_REF] Garis | Genetic Programming: building arti cial nervous systems using genetically programmed neural networks modules[END_REF], [START_REF] Schoenauer | Iterated Genetic Algorithms[END_REF]). All experiments presented in this paper are based on the penalized tness, and involve this 2-step optimization.

First results

Figure 3 show the rst results on the 2 1 cantilever plate, discretized according to a 32 22 regular mesh. The population size for all runs is 125, and the number of generations arbitrarily xed to 1000. One run thus require about 100000 FEM analyses, taking approximately 6 hours of a powerful HP workstation for the 32 22 discretization. The genetic operators are described in section 3.2: The block crossover is applied at a rate of 0.6, and the epistatic mutation at a rate of 0.1 (the probability to ip a given pixel varies from 0 to 0.01, depending of the diversity of the values for this particular pixel in the whole population, following a parabolic rule). All these parameters were adjusted after exhaustive tests, and details can be found in [START_REF] Kane | Algorithmes g en etiques et Optimisation topologique de formes[END_REF].

Each one of the structures on Figure 3 corresponds to a given value of the constraint on the maximum displacement D Lim . And each one is the most significant result for that value of D Lim obtained out of 5 runs with di erent random initial populations. Due to the stochastic aspect of the GA, the results differ from one run to the other. Nevertheless, on the 2 1 cantilever plate, such di erences are not very important. However, note that the second structure of Figure 3, though optimized using a stronger constraint on the displacement, is lighter and more rigid than the third one. This suggested a modi cation of the geometry and boundary conditions to obtain a problem exhibiting known multiple solutions.

Multiple solutions

The problem considered in that section is the 1 4 cantilever plate, discretized according to a 10 40 mesh and for which both left and bottom boundaries are xed (Figure 6-a). Depending on the height of the point where the loading is applied and on the constraint on the displacement, the same problem can have multiple solutions. A simple example of such situation is given in Figure 6 (b) and (c): the loading is applied at height 10, and, provided that the displacement constraint is large enough, both structures (a) and (b) are optimal solutions. Moreover, if the height of the loading point is xed, and the displacement constraint D Lim is gradually relaxed, di erent quasi-optimal solutions exist for some ranges of D Lim . And the GA method is then able to nd out such multiple solutions, as demonstrated in Figure 5. Here again only the most signi cant results are shown, for di erent values of D Lim , while the loading is applied at height 15. The GA parameters are those of section 4.3, except for the population size (100) and the number of generations (500).

Compliance tness

In this section, a compliance-based tness is used in order to compare our results with those of the homogenization method:

Fitness compliance = 1 Area + C where C = R Fu is the energy of the external load, or compliance of the structure under the load F.

The comparisons are made on the 1 2 cantilever plate using a 10 20 discretization. Figure 7 shows typical solutions of the GA approach for taking values 1, 0.1 and 0.01. The population size is 75, runs are limited to 500 generations, other parameters are described in section 4.3. = 1 = 0:1 = 0:01 Figure 7: The GA solutions, for di erent values of , using the compliance based tness, on a coarse mesh.

The results of Figure 7 strongly di er from those of the homogenization method (the perfect ">" shape): Figure 8-a shows the results of the latter method, as described in [START_REF] Allaire | Optimal design for minimum weight and compliance in plane stress using extremal microstructures[END_REF], for = 1 on the same coarse mesh. The strength of the homogenization method lies in its convergence property, when the element size of the mesh vanishes, and whatever the parameter , toward the optimal homogenized solution, i.e. a solution with continuous density of material between 0 (void) and 1 (pure material). However, its Achilles'heel sometimes appears when it comes to design an actual part, in which the density has to be boolean (0 or 1). The usual method of using a relaxed cost function favoring a 0/1 density at the end of homogenization might still give results far from a 0/1 material, as attested by Figure 8-a (the density of grey is proportional to the density of material of the solution). Of course, things get better when rening the mesh, and Figures 8 (b) and (c) witness that a more feasible solution can be reached easily for the 2 1 cantilever plate. But the results of Figure 7 demonstrate the ability of the GA method to provide good solutions for a given mesh, though totally di erent from those of the homogenization method. Anyway, the aim of GA-based topological optimization is not to compete with the homogenization method, as GAs weakness is well known to be its computational cost. Next section presents results that cannot be obtained by homogenization, in a more complex mechanical background.

NONLINEAR GEOMETRICAL EFFECTS IN TOPOLOGY OPTIMIZATION

This section considers standard plane stress problems in the context of large displacements. The material still obeys a linear law (the extension to any other law is straightforward), but the nonlinear geometric e ects due to the large displacement hypothesis are taken into account. A thorough description of the theoretical model together with the numerical algorithm can be found in [START_REF] Ciarlet | Mathematical Elasticity, Vol I : Three-Dimensional Elasticity[END_REF]. Details on the FEM model and implementation used here can be found in [START_REF] Jouve | Mod elisation math ematique de l' il en elasticit e non-lin eaire Recherches en Math ematiques Appliqu ees[END_REF].

The rst experiments use the penalized tness function described in section 3.3. The idea was to use di erent loadings F, with a xed ratio F=D Lim (D Lim is the constraint on the displacement). In the purely linear case, all such problems are of course equivalent.

A typical (best) result of the GA based optimization in the large displacement context is shown in Figure 9-b, together with the result on the same problem in the pure linear model (Figure 9-a). And di erent values of F and D Lim do produce similar results. But a closer look at the maximal stress of the solutions Max gives some hint on what is hapenning: First, the large displacement model can have di erent solutions for the same loading, some with higher stress than others. Second, the stress eld itself, on a rough domain like the structures obtained for the standard cantilever plate problem, and with such coarse discretization, presents some singularities. The value of Max in Figure 9 Moreover, Table 1 gives some idea of the nonlinear geometrical e ects in the case of two simple structures, the perfect ">" shape and the straight beam; both the displacement and the maximum stress are given for varying loads. This suggests to incorporate a constraint on the maximal stress in the tness function. Of course, some nice structures might be missed, like the ">" shape for F = 910 3 . But, hopefully, another structure very similar to it can arise, with the same mechanical properties, and with the FEM model giving a reasonable solution among the possible solutions. Following these ideas, the tness becomes

1 Area + (D Max D Lim ) + + ( Max Lim ) +
where Lim is the constraint on the stress, Max is the maximum of the stress on the structure, and is some positive penalty parameter which has to be adjusted. Figure 10 shows the optimal designs obtained with that modi ed tness: Reasonable solutions are found, as long as the stress is imposed a strong constraint.

CONCLUSION AND PERSPECTIVES

The feasibility of GA-based optimal design has already been witnessed by [START_REF] Jensen | Topological Structural Design using Genetic Algorithms[END_REF] and [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF]. This paper intends to emphasize on its exibility: When multiple quasi-optimal solutions exist, all can be obtained by this method. Actually, these were found by successive runs of the GA starting from di erent random populations. But it is possible, using niching techniques like the sharing scheme of [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF], to locate in a single run more than one near-optimal solutions.

In the context of large displacements, the experiments presented in this paper acknowledge for the usual "the extension to other mechanical models is straightforward", though the tness function had to be adjusted in order to take into account the maximal stress of the structure: there is no free lunch for nonlinearity.

Many problems remain to be addressed. The tuning of the penalization parameter, even in the context of iterated tness functions presented in this paper, remains to be done for every new problem. The idea of evolving this parameter together with the current population will be investigated.

The comparison with the results of the homogenization method shows the need to use ner meshes. The separation of the design domain into sub-domain proposed by [START_REF] Chapman | Genetic Algorithms as an approach to Conguration and Topology Design[END_REF] is a promising direction. Some iterated GA on the same population of designs, the mesh being gradually re ned, will also be tested.

Further research also include straightforward (!) modi cations of the tness function, to handle some multiple loading problems, and to extend this work to any nonlinear material.

Figure

  Figure 1: The diagonal crossover

Figure 4 :

 4 Figure 4: The standard 2 1 cantilever plate.

  The homogenization solutions, for = 1, using di erent meshes.

  -b witnesses that phenomenon.

Figure 9 :

 9 Figure 9: Optimal designs with displacement constraints. F = 0:009 and D Lim = 0:02285.

F

  Optimal cantilever plates for di erent values of the displacement constraint D Lim . D Max is the actual maximum displacement of the structure.

	D Lim D Max Area	2.463 2.439 0.2075	0.1293 0.129286 0.4675	0.133 0.133 0.49
	Figure 3:

are arbitrarily assigned zero tness, and are therefore eliminated by the next Darwinian selection.

  Optimal structures on problems with multiple solutions. F is applied at height 15.

	D Lim D Max Area	0.32 0.294 0.38	0.3927 0.33	0.4	0.3919 0.34	0.42 0.3415 0.26	7.66 0.15	7.66	6.743 0.19	46.6 46.64 0.1
		Figure 5:

Table 1 :

 1 The displacement is what was expected, but the stress does present weird values. Nonlinear e ects on two reference shapes. D Max and Max are respectively the maximum displacement and the maximum stress.

	">" shape D Max Max 910 6 0.0002360 0.0005968 0.04161 0.03035 straight beam Load D Max Max 910 5 0.00235 0.0059 0.34499 0.28344 910 4 0.02286 0.05312 0.80763 1.272 910 3 1.06731 10.99 1.08208 1.11683 910 2 1.249 2.0763 1.6871 4.3787 910 1 2.9574 break 3.177 break

Most authors only consider one-point crossovers, corresponding to masks (1::10::0), or two-point crossovers, corresponding to masks (1::10::01::1). The general case represented here is called uniform crossover, introduced by[START_REF] Syswerda | Uniform crossover in genetic algorithms[END_REF].

Acknowledgments

The authors are indebted to Mich ele Sebag and to Jo el Frelat (LMS -Ecole Polytechnique) for many useful discussions and thorough proofreading. Thanks too to the anonymous reviewer whose critical comments were quite helpful to improve this paper.