Isabelle Gallagher
email: isabelle.gallagher@math.jussieu.fr

AND Thierry Paul
email: thierry.paul@math.polytechnique.fr

Laure Saint-Raymond
email: laure.saint-raymond@ens.fr

ON THE PROPAGATION OF OCEANIC WAVES DRIVEN BY A MACROSCOPIC CURRENT

Keywords: 2010 Mathematics Subject Classification. 35Q86, 76M45, 35S30, 81Q20 Semiclassical analysis, microlocal analysis, Mourre estimates, Geophysical flows

, where the profile was assumed to be zonal only and where explicit calculations were made possible due to the 1D setting.

Here the diagonalization of the system, which allows to identify Rossby and Poincaré waves, is proved by an abstract semi-classical approach based on normal forms. The dispersion of Poincaré waves is also obtained by a more abstract and more robust method using Mourre estimates. Only some partial results however are obtained concerning the Rossby propagation, as the two dimensional setting complicates very much the study of the dynamical system.

Introduction

This paper is a continuation of [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] so before discussing the matter of this paper (in Section 2) let us review the contents of that work. We shall start by recalling briefly the model, then we shall explain the methods and results obtained in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] and discuss their limitations.

1.1. The model. The goal of [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] is to understand, through the study of a toy model, the persistence of oceanic eddies observed long past by physicists among which [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF][START_REF] Gill | Resonant interactions between planetary waves[END_REF][START_REF] Greenspan | The theory of rotating fluids[END_REF][START_REF] Majda | Introduction to PDEs and waves for the atmosphere and ocean[END_REF][START_REF] Pedlosky | Geophysical fluid dynamics[END_REF][START_REF] Pedlosky | Ocean Circulation Theory[END_REF], who gave heuristic arguments to explain their formation due both to wind forcing and to convection by a macroscopic current.

The ocean is considered in this toy model as an incompressible, inviscid fluid with free surface submitted to gravitation and wind forcing, and we further make the following classical assumptions: we assume that the density of the fluid is homogeneous ρ = ρ 0 = constant, that the pressure law is given by the hydrostatic approximation p = ρ 0 gz, and that the motion is essentially horizontal and does not depend on the vertical coordinate. This leads to the so-called shallow water approximation.

For the sake of simplicity, the effects of the interaction with the boundaries are not discussed and the model is purely horizontal with the longitude x 1 and the latitude x 2 both in R.

The evolution of the water height h and velocity v is then governed by the shallow-water equations with Coriolis force (1.1) ∂ t (ρ 0 h) + ∇ • (ρ 0 hv) = 0

∂ t (ρ 0 hv) + ∇ • (ρ 0 hv ⊗ v) + ω(ρ 0 hv) ⊥ + ρ 0 gh∇h = ρ 0 hτ
where ω denotes the vertical component of the Earth rotation vector Ω, v ⊥ := (-v 2 , v 1 ), g is the gravity and τ is the -stationary -forcing responsible for the macroscopic flow. The vertical component of the Earth rotation is therefore Ω sin(x 2 /R), where R is the radius of the Earth; note that it is classical in the physical literature to consider the linearization of ω (known as the betaplane approximation) ω(x 2 ) = Ωx 2 /R. We consider here general functions ω in the sequel, with some restrictions that are be made precise later.

We consider small fluctuations around the stationary solution ( h, v) satisfying

h = constant, ∇ • (v ⊗ v) + ωv ⊥ = τ, div v = 0.
In [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] the study is restricted to the case of a shear flow, in the sense that v(x) = (v 1 (x 2 ), 0), with v1 a smooth, compactly supported function. Moreover some orders of magnitude and scalings allow to transform the previous system into the following one:

(1.2)

∂ t η + 1 ε ∇ • u + ū • ∇η + ε 2 ∇ • (ηu) = 0 , ∂ t u + 1 ε 2 bu ⊥ + 1 ε ∇η + ū • ∇u + u • ∇ū + ε 2 u • ∇u = 0 ,
where (u, η) are the fluctuations of the velocity and height of the water, b := ω/|Ω| and where ε is a small parameter (of the order of Fr 2 , where Fr is the Froude number, and of Ro

1 2
where Ro is the Rossby number).

1.2. Methods and results in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF]. Most of the analysis in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] concerns the linear version of (1.2), namely the following system:

(1.3)

ε 2 i∂ t v + A(x 2 , εD, ε)v = 0 v = (v 0 , v 1 , v 2 ),
where D := 1 i ∂, and the linear propagator is given by

A(x 2 , εD, ε) := i   εū 1 ε∂ 1 ε∂ 1 ε∂ 2 ε∂ 1 εū 1 ε∂ 1 -b(x 2 ) + ε 2 ū 1 ε∂ 2 b(x 2 ) εū 1 ε∂ 1   ,
The first step of the analysis consists in diagonalizing (approximately) the system (1.3). The computation of a kind of characteristic polynomial associated with (1.3), in symbolic form, allows to construct three symbols the quantization of which provides three scalar propagators.

Two of those propagators, called Poincaré propagators, are then proved to satisfy dispersive estimates; that result relies on a spectral analysis (usual semi-classical theory does not operate here due to the very large time scales at play) using Bohr-Sommerfeld quantization, which requires that b has only one, non degenerate critical value and which also uses very much the fact that the motion is translation-invariant in x 1 . A stationary phase argument on the spectral decomposition of any solution to the Poincaré propagation gives the result: Poincaré modes exit any compact set in finite time.

The last propagator is the Rossby one, which is one order of magnitude (in ε) smaller than the Poincaré modes. This allows to analyse the propagation by semi-classical analysis tools.

In particular the precise study of the dynamical system associated with those waves, which is an integrable system due to translation invariance in x 1 , allows to derive a condition on the initial microlocalization of the solution which guarantees that the Rossby waves are trapped for all times in a compact set.

Those results on the linear system (1.3) can finally be transposed to the original system (1.2) due to the high power of ε in front of the nonlinearity, and due to the semi-classical setting, which allows to exhibit vector fields which almost-commute with the linear operator A(x 2 , εD, ε).

1.3.

Limitations of the methods of [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF]. The restriction which is the most used in the analysis described briefly in the previous paragraph is the fact that the stationary flow ū is a shear flow of the type ū = (ū 1 (x 2 ), 0). Indeed

• It allows to Fourier-transform in the direction x 1 , which makes the diagonalization procedure much easier; • It simplifies the spectral analysis of Poincaré waves, again due to the Fourier transform (in particular the dual variable ξ 1 is fixed during the propagation, and there is a wavelike behaviour in x 1 ); • It allows the Rossby dynamical system to be integrable, which is a tremendous help in the analysis.

However it is important for physical reasons to consider truly 2D convection flows.

An additional restriction in the previous arguments is that in order to prove the dispersion of Poincaré waves, the rotation amplitude b should have only one, non degenerate critical value: this allows to use a Bohr-Sommerfeld quantization argument to compute the eigenvalues of the Poincaré operator. This assumption on b is not really restrictive from the physical point of view.

1.4. On the nonlinear term. As explained above, most of the analysis in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] is concerned with the linear system (1.3). In order to transpose the linear results to the nonlinear setting, one uses the following arguments (along with the fact that the coupling is vanishing when ε goes to zero):

• Uniform existence which is obtained via an almost-commutation result;

• Bilinear estimates in anisotropic semi-classical spaces;

• A Gronwall lemma, which requires an L ∞ (R 2 ) bound on the linear solution. This is not known in general, due to the bad Sobolev embeddings in semi-classical settings, so the nonlinear result is proved for vanishing couplings only.

It is important to notice that none of those three steps require that ū is a shear flow. In the whole of the present paper we shall therefore only focus on the linear equation, and leave to the reader the transposition to the nonlinear equation, using the above steps.

2. Main result of this paper and strategy of the proof 2.1. The model. In this paper we shall be concerned with the linear system (2.1)

ε 2 i∂ t v + A(x, εD, ε)v = 0 v = (v 0 , v 1 , v 2 ),
where the linear propagator is given by

(2.2) A(x, εD, ε) := i   εū • ε∇ ε∂ 1 ε∂ 2 ε∂ 1 εū • ε∇ + ε 2 ∂ 1 ū1 -b(x 2 ) + ε 2 ∂ 2 ū1 ε∂ 2 b(x 2 ) + ε 2 ∂ 1 ū2 εū • ε∇ + ε 2 ∂ 2 ū2   .
We shall assume throughout the paper that b is smooth, with a symbol-like behaviour: for all α ∈ N, there is a constant C α such that for all x 2 ∈ R,

(2.3) |b (α) (x 2 )| ≤ C α 1 + b 2 (x 2 ) 1 2 .
We shall further assume that lim

|x 2 |→∞ b 2 (x 2 ) = ∞,
and for the sake of simplicity, that b 2 has only non degenerate critical points.

We shall also suppose that the initial data is microlocalized in some compact set C of T * R 2 satisfying (2.4)

C ∩ {ξ 2 1 + ξ 2 2 + b 2 (x 2 ) = 0} = ∅
or actually rather

(2.5) C ∩ {ξ 1 = 0} = ∅ .
We shall prove that assumption (2.4) is propagated by the flow, while (2.5) is propagated only by the Poincaré component. We recall (see for instance [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF], Appendix B) that a function f is microlocalized in a compact set C of T * R 2 if for any (x 0 , ξ 0 ) in the complement of C in R 4 (we shall identify T * R 2 to R 4 in the following), there is a smooth function χ 0 , bounded as well as all its derivatives and equal to one at (x 0 , ξ 0 ), satisfying

(2.6) Op W ε (χ 0 )u 0 L 2 (R 2 ) = O(ε ∞ ),
where Op W ε denotes the Weyl quantization:

(2.7)

Op W ε (χ 0 )u 0 (x) := 1 (2πε) 4 e i(x-y)•ξ/ε χ 0 ( x + y 2 , ξ)u 0 (y) dydξ.
We also recall that (2.6) means that for any N ∈ N, there are ε 0 and C such that

∀ε ∈]0, ε 0 ], Op W ε (χ 0 )u 0 L 2 (R 2 ) ≤ Cε N .
In the following, to simplify some formulations, we shall denote by (µ)Supp f the projection of the (micro)support of f onto the = 0 axis.

2.2.

Statement of the main result and organization of the paper. Let us state the main theorem proved in this paper.

Theorem 1. Let v ε,0 be a family of initial data, microlocalized in a compact set C satisfying Assumptions (2.4)-(2.5). For any parameter ε > 0, denote by v ε the associate solution to (2.1). Then for all t ≥ 0 one can write v ε (t) as the sum of a "Rossby" vector field and a "Poincaré" vector field: v ε (t) = v R ε (t) + v P ε (t), satisfying the following properties:

(1) µSupp v R ε (t) and µSupp v P ε (t) satisfy (2.4) for all times.

(2) For any compact set Ω in R 2 , one has

∀t > 0, v P ε (t) L 2 (Ω) = O(ε ∞ ). (3) µSupp x 2 v R ε (t) lies in a bounded subset of R uniformly in time.
Compared to [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF], the main difficulties are due to the presence of a x 1 -dependent underlying flow ū. The diagonalization of the system (exhibiting Rossby and Poincaré-type waves, with very different qualitative features) must be revised, and obtained in a less explicit way. Moreover the proof of (2) in Theorem 1, namely the dispersion of Poincaré waves can also not be proved in the same way (note that it is not assumed here that b 2 has only one non degenerate critical value). Finally the trapping of Rossby waves seems much harder to obtain since the underlying dynamical system no more decouples; the behaviour of the Rossby waves is therefore much less precise than in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF].

Let us explain our strategy here, compared with that in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] described above.

2.2.1. The diagonalization. The construction of the Rossby and Poincaré modes is not as direct as in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF] due to the lack of translation invariance in x 1 . We choose therefore to follow a more abstract way to recover those modes in Section 3, which relies on semi-classical analysis, and normal forms (instead of explicit computations as in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF]). Finding the propagators associated with those modes requires a microlocalisation assumption of the type (2.4), in order for the eigenvalues of the matrix of principal symbols to be well separated. The diagonalization result is therefore in this paragraph conditional to the fact that the solution to the propagation equation is correctly microlocalized (that corresponds to Point 1 of Theorem 1).

Dispersion of Poincaré waves and propagation of the nondegeneracy assumption (2.5).

In order to prove (2) in Theorem 1 we again rely on a more abstract, and more efficient method than that followed in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF]. It is based on Mourre estimates and the assumption (2.5) on the initial data: we start by proving, by a semi-classical argument, that after a very short time (of the order of ε) the support in x 1 of the solution escapes the support of ū. Then we use Mourre estimates to prove that the solution remains outside the support of ū for all times, and actually escapes any compact set in x 1 in finite time (to prove this last point we use the fact that the equation reduces to a translation-invariant equation in x 1 since the support of the solution is outside the support in x 1 of ū). This allows finally to check that the nondegeneracy assumption (2.5) does hold for all times. This analysis is achieved in Section 4.

Study of Rossby waves and propagation of the nondegeneracy assumption (2.4).

In Section 5 we first prove that the nondegeneracy assumption (2.4) does hold during the propagation of Rossby waves. That is due to semi-classical analysis, by the study of the dynamical system associated with those waves. The study of that system is also the key to the proof of Point (3), which is also proved in Section 5. This fact indicates that the only possible escape is along x 1 .

Reduction to scalar propagators

In this section we shall construct the three operators T + , T -and T R diagonalizing A(x, εD, ε).

We shall start by proving a general diagonalization result, and at the end we shall apply the general result to our context.

Before stating the general result, let us give some notation. A semi-classical symbol is a function a = a(x, ξ; ε) defined on R 2d ×]0, ε 0 ] for some ε 0 > 0, which depends smoothly on (x, ξ) and such that for any α ∈ N 2d and any compact K ⊂ R 2d , there is a constant C such that for any ((x, ξ), ε)

∈ K×]0, ε 0 ], |∂ α a((x, ξ), ε)| ≤ C.
We shall consider the Weyl quantization of such symbols, as recalled in (2.7). We shall say that a pseudodifferential operator Op W ε (a) is supported in a set K if for any smooth function χ equal to one in a neighborhood of K one has aχ = a.

Finally we shall say that a matrix is pseudodifferential if each of its entries is a pseudodifferential operator.

Let us first prove the following general result. Theorem 2. Let K be a compact subset of R 2d , and consider a N × N hermitian pseudodifferential matrix A ε = A(x, εD, ε), supported in K. Assume that • the (matrix) principal symbol of A(x, εD, 0), denoted by A 0 , is diagonalizable, in the sense that there are some unitary and diagonal matrices of symbols, U and D, such that

U -1 A 0 U = D, • the eigenvalues (δ 1 (x, ξ), . . . , δ N (x, ξ)) satisfy (3.1) ∀i = j, inf (x,ξ)∈K |δ i (x, ξ) -δ j (x, ξ)| ≥ C > 0.
Then there exists a family of unitary and diagonal pseudodifferential operators V ε and D ε supported in K, such that:

(3.2) V * ε A ε V ε = D ε + O(ε ∞ ), V * ε V ε = I + O(ε ∞ ). Moreover one has (3.3) D ε = Op W ε (D) + εD 1 + O(ε 2 )
, where the principal symbol of D 1 is given by

D 1 = diag ∆ 1 - D 0 I 1 + I 1 D 0 2
with the notations

(3.4) ∆ 1 = 1 ε Op W ε (U * )A ε Op W ε (U) -D 0 , I 1 = 1 ε Op W ε (U * )Op W ε (U) -I
More explicitly, let us denote by a ij (x, ξ) the matrix elements of A 1 , subsymbol of A(x, εD, ε) defined by:

A 1 (x, ξ) := ∂A ∂ε (x, ξ, 0)
and by u nj (x, ξ), i = 1 . . . N, the coordinates of any unit eigenvector of A 0 (x, ξ) of eigenvalue δ n (x, ξ). We have

(3.5) (D 1 ) nn = jk=1...N (u jn {a jk , u kn }) + a jk {u jn , u kn } 2i + (U * A 1 U) nn ,
where {f, g}

:= ∇ ξ f ∇ x g -∇ x f ∇ ξ g is the Poisson bracket on T * R n .
Here and in all the sequel, we say that a pseudo-differential operator V is unitary if it satisfies

V * V = I + O(ε ∞ ).
The proof is divided into two parts: in Section 3.1 we present the formal construction and in Section 3.2 we show that the series of symbols formally constructed are indeed symbols. Finally Section 3.3 is devoted to the case of the matrix given by (2.2).

3.1. The formal construction. The proof of Theorem 2 is a combination of semiclassical and perturbation methods. Let us start by defining

U 0 = Op W ε (U).

Elementary properties of the Weyl quantization imply then that U

* 0 A ε U 0 = D 0 + O(ε).
The following proposition shows that one can construct a unitary pseudodifferential operator

U ∞ such that U * ∞ A ε U ∞ = D 0 + O(ε). Lemma 3.1. Let U be a pseudodifferential matrix such that U * U = I + εI 1 , where I is the identity. Then one can find V ∼ ∞ k=0 ε k V k such that (3.6) (U + εV ) * (U + εV ) = I + O(ε ∞ ). Proof. Let us denote V 0 := - 1 2 I 1 U . On easily check that (U + εV 0 ) * (U + εV 0 ) = I + O(ε 2 ). Indeed (U + εV 0 ) * (U + εV 0 ) = U * U + ε 2 (I 1 U * U + U * U I 1 ) + O(ε 2 ) = I + εI 1 -εI 1 + O(ε 2 ).
Then one concludes by iteration.

That lemma allows to define the pseudo-differential operator of (semiclassical) order 0

∆ 1 = 1 ε (U * ∞ A ε U ∞ -D 0 ) ,
where U ∞ is a unitary operator. Now our aim is to find a unitary operator

V ∞ (up to O(ε ∞ )) such that (U ∞ V ∞ ) * A ε (U ∞ V ∞ ) = D ∞ + O(ε ∞ ),
where D ∞ = D 0 + εD 1 + . . . is a diagonal matrix satisfying the conclusions of the theorem.

We shall write V ∞ = e iεW , with W selfadjoint (so V ∞ thus constructed is automatically unitary). We recall that if W is a pseudodifferential operator, then so is e iεW (simply by writing

e iεW ∼ ∞ 0 (iε) k k! W k ).
We look for W under the form W ∼ ∞ 0 ε k W k , and compute the W k recursively. Since

V * ∞ (D 0 + ε∆ 1 )V ∞ = (D 0 + ε∆ 1 ) + iε[(D 0 + ε∆ 1 ), W ] + (iε) 2 2 [[(D 0 + ε∆ 1 ), W ], W ] + . . . we see that, if W 1 satisfies (3.7) [D 0 , W 1 ] + ∆ 1 = D 1 + O(ε 2 ), D 1 diagonal,
then we have that

(3.8) e -iεW 1 (D 0 + ε∆ 1 )e iεW 1 = D 0 + εD 1 + ε 2 ∆ 2 ,
where ∆ 2 is a zero order pseudodifferential operator. The following lemma is a typical normal form type result, and is crucial for the following.

Lemma 3.2. Let D 0 be a diagonal pseudodifferential matrix whose principal symbol D 0 has a spectrum satisfying (3.1) and let ∆ 1 be a pseudodifferential matrix.

Then there exist two pseudodifferential matrices W and D 1 , with D 1 diagonal, such that:

(3.9) [D 0 , W ] + ∆ 1 = D 1 + ε∆ 2 ,
where ∆ 2 is a pseudodifferential matrix of order 0.

Moreover the principal symbol of D 1 is the diagonal part of the principal symbol of ∆ 1 .

Proof. By the non degeneracy condition of the spectrum of D 0 we know, by standard arguments (see [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF] for instance), that there exists a matrix W 0 and a diagonal one

D 1 such that [D 0 , W 0 ] + D 1,0 = D 1 , where D 1,0 is the principal symbol of ∆ 1 .
Indeed it is enough to take D 1 as the diagonal part of D 1,0 and

(3.10) (W 0 (x, ξ)) i,j = (D 1,0 (x, ξ)) i,j δ i (x, ξ) -δ j (x, ξ)
and notice that the Weyl quantization of W 0 satisfies (3.9). By Lemma 3.2 we know that there exists W 1 satisfying (3.7). Developing (3.8) as

e -iεW 1 (D 0 + ε∆ 1 )e iεW 1 = D 0 + ε(∆ 1 + [D 0 , W 1 ]) + ε 2 ∆ 2 ,
we get immediately (3.8).

It is easy to get convinced that all the W k will satisfy recursively an equation of the form

[D 0 , W k ] + ∆ k = D k + O(ε),
which can be solved by Lemma 3.2.

The expression for the principal symbol of D 1 follows by construction and the following well known lemma (see [START_REF] Martinez | An introduction to semiclassical and microlocal analysis[END_REF] for instance): In order to derive (3.5) we have to compute the subprincipal symbol of the diagonal part of the right-hand side of (3.4), that is, for each n = 1 . . . N ,

jk Op W ε (U jn )Op W ε ((A 0 + εA 1 ) jk )Op W ε (U kn ),
since U is unitary.

The term εA 1 is obviously responsible for the second term in the right-hand side of (3.5). Using Lemma 3.3 and the distributivity of the Poisson bracket, we get the following expression for the first one:

jk 1 2i {U jn , (A 0 ) jk U kn } + U jn {(A 0 ) jk , U kn } = jk 1 2i U jn {(A 0 ) jk , U kn } + (A 0 ) jk {U jn , U kn } + U kn {U jn , (A 0 ) jk }
Interverting j and k in half of the terms and noticing that, since A 0 is Hermitian, (A 0 ) jk = (A 0 ) kj , we get easily (3.5).

Symbolic properties.

With the hypothesis that both D and U are pseudodifferential matrices it is quite obvious that V ε and D ε are pseudodifferential matrices as well. Indeed the formal construction in the preceding section shows that the iterative process uses only three things: multiplications of pseudodifferential operators, computation of subprincipal symbols and solving equation (3.9).

For (3.9), the formula (3.10) used in the proof of Lemma 3.2, together with the non-degeneracy condition (3.1) which shows clearly that (δ i (x, ξ) -δ j (x, ξ)) -1 is a symbol, implies that W 0 is a pseudodifferential operator.

Note that the microlocalization assumption is crucial in order that the expansions obtained by this iterative construction do define symbols. We have indeed no uniform control on the growth at infinity.

The Rossby-Poincaré case.

In the case of oceanic waves A(x, εD, ε) is given by (2.2):

A(x, εD) := i   εū • ε∇ ε∂ 1 ε∂ 2 ε∂ 1 εū • ε∇ + ε 2 ∂ 1 ū1 -b(x 2 ) + ε 2 ∂ 2 ū1 ε∂ 2 b(x 2 ) + ε 2 ∂ 1 ū2 εū • ε∇ + ε 2 ∂ 2 ū2   . Therefore A 0 (x, ξ) :=   0 ξ 1 ξ 2 ξ 1 0 -ib(x 2 ) ξ 2 ib(x 2 ) 0   ,
and

A 1 (x, ξ) :=   ū • ξ 0 0 0 ū • ξ 0 0 0 ū • ξ   = ū • ξ   1 0 0 0 1 0 0 0 1   .
A straightforward computation shows that the spectrum of

A 0 is 0, ξ 2 1 + ξ 2 2 + b 2 (x 2 ), -ξ 2 1 + ξ 2 2 + b 2 (x 2 ) .
3.3.1. Microlocalization. The three eigenvalues of A 0 are separated if and only if

ξ 2 1 + ξ 2 2 + b 2 (x 2 ) = 0. Therefore, considering a compact subset K of R 2d such that K ∩ {(x 1 , x 2 , ξ 1 , ξ 2 ) / ξ 2 1 + ξ 2 2 + b 2 (x 2 ) = 0} = ∅ ensures that
• the eigenvalues do not cross, so that it is possible to get a unitary diagonalizing matrix with regular entries; • the non degeneracy condition (3.1) is satisfied.

In other words, A(x, εD, ε) satisfies the assumptions of Theorem 2 provided that one considers only its action on vector fields which are suitably microlocalized.

We assume of course that this microlocalization condition is satisfied by the initial datum, which is the condition (2.4). Furthermore, we will prove in the next two sections that the propagation by the scalar operators T ± and T R (to be defined now) preserves this suitable microlocalization, thus justifying a posteriori the diagonalization procedure for all times.

Computation of the Poincaré and Rossby Hamiltonians. The above computations show that one can define the two Poincaré Hamiltonians as follows:

τ ± := ± ξ 2 1 + ξ 2 2 + b 2 (x 2
) and we shall denote the associate operator constructed via Theorem 2 by T ± . Now let us consider the Rossby Hamiltonian. In all this paragraph, for the sake of readability, we will denote

ξ b = ξ 2 1 + ξ 2 2 + b 2 (x 2
) . An easy computation shows that a (normalized) eigenvector of A 0 (x, ξ) of zero eigenvalue is

u 0 = 1 ξ b   b iξ 2 -iξ 1  
By Theorem 2, the Rossby Hamiltonian is then given by the formula

(3.11) τ R = j,k=1...3 (u j0 {a jk , u k0 }) + a jk {u j0 , u k0 } 2i + j,k=1...3 (A 1 ) jk u j0 u k0 .
In order to compute the different Lie brackets, we start with a couple of simple remarks :

{ξ j , f } = ∂ x j f and {b(x 2 ), f } = -b (x 2 )∂ ξ 2 f .
In particular, if f does not depend on x 1 , then {ξ 1 , f } = 0.

The contribution of the first term in the parenthesis in (3.11) is

j,k=1...3 (u j0 {a jk , u k0 }) = b ξ b ξ 2 , -iξ 1 ξ b + iξ 2 ξ b ib, -iξ 1 ξ b + iξ 1 ξ b ξ 2 , b ξ b + ib, iξ 2 ξ b = -ibξ 1 ξ b ∂ x 2 1 ξ b - iξ 2 ξ 1 b ξ b ∂ ξ 2 1 ξ b + iξ 1 ξ b ∂ x 2 b ξ b + iξ 1 b ξ b ∂ ξ 2 ξ 2 ξ b = 2iξ 1 b ξ 2 b •
Using the distributivity of the Poisson brackets, we get the contribution of the second term in a very similar way j,k=1...3

a jk {u j0 , u k0 } 2 = ξ 1 b ξ b , iξ 2 ξ b -ξ 2 b ξ b , iξ 1 ξ b + ib iξ 1 ξ b , iξ 2 ξ b = ξ 1 ib ξ b 1 ξ b , ξ 2 + iξ 2 ξ b b, 1 ξ b + i ξ 2 b {b, ξ 2 } - iξ 2 ξ 1 ξ b b, 1 ξ b - ibξ 1 ξ b 1 ξ b , ξ 2 = -ib ξ 1 ξ 2 b •
The computation of the second term of the right hand side of (3.11) is trivial since A 1 is a multiple of the identity. Adding with the two previous expressions we get finally

τ R = ξ 1 b ξ 2 1 + ξ 2 2 + b(x 2 ) 2 + ū
• ξ and the associate operator will be denoted by T R . Remark 3.4. Since the elementary steps of the diagonalization process use only multiplications, computations of subprincipal symbols and solving normal forms equations, all the subsymbols of T R and T ± depend on x 1 only through ū and its derivatives.

Study of the Poincaré waves

In Section 3 we constructed two linear operators, called T ± , whose principal symbols are

τ ± = ± ξ 2 1 + ξ 2 2 + b 2 (x 2 ).
We now want to study the propagation equation associated to those operators, namely the linear equation in

R × R 2 (4.1) iε 2 ∂ t ϕ ± = T ± ϕ ± , ϕ ±|t=0 = ϕ 0 ±
where ϕ 0 ± are microlocalized in a compact set C satisfying Assumption (2.4). Before studying that equation we need to check that it makes sense, since a priori T ± is only defined on vector fields microlocalized on such a compact set. This is achieved in the coming section, where we check that the separation of eigenvalues (3.1) required in the statement of Theorem 2 holds because ξ 2 1 + ξ 2 2 + b 2 (x 2 ) remains bounded away from zero during the propagation.

Then we shall show that the solutions to these equations exit any compact set in finite time (Point (2) of Theorem 1).

4.1. Microlocalization. Let us prove the following result, which provides the first part of Point (1) in Theorem 1 and allows to make sense of Equation (4.1) for all times.

Proposition 4.1. Under the assumptions of Theorem 1, the operators T ± are self-adjoint, and the function ϕ(t) = e i t ε 2 T ± ϕ 0 ± are such that µSuppϕ ± (t) satisfies (2.4) for all times.

Proof. The proof of that result relies on a spectral argument. Due to the form of the principal symbols of T ± recalled above, the operators T ± are self-adjoint. We can therefore define two families (ψ ± n ) n∈N of (pseudo)-eigenvectors of T ± in L 2 (R 2 ) and two sequences of eigenvalues λ ± n such that if the initial data writes

ϕ 0 ± (x) = n c ±,0 n ψ ± n (x), then ϕ ± (t, x) = n e i λ ± n t ε 2 c ±,0 n ψ ± n (x).
Since the eigenfunctions ψ ± n are microlocalized on the energy surfaces of the Poincaré Hamiltonians, the result follows.

Dispersion.

In this paragraph we shall prove Point (2) of Theorem 1. The strategy is the following. In Section 4.2.1 we prove using semi-classical analysis that for a very short time, the solution to (4.1) remains microlocalized in a compact set satisfying assumption (2.5), and such that µSupp x 1 ϕ ± becomes disjoint from Supp x 1 ū. Section 4.2.2 is then devoted to the long-time behaviour of the solution, and Mourre estimates allow to prove that the solution exits any compact set after some time, and that it remains microlocalized far from ξ 1 = 0.

The result of the analysis carried out in this paragraph is that the behaviour of µSupp x 1 ϕ ± is as depicted in the following figure .   4.2.1. Short time behaviour. The aim of this paragraph is to prove the following result. It shows that the solutions of (4.1) exit the support of ū after a time t exit ε, for |t exit | large enough (independent of ε). We only state the forward in time result: the backwards result is identical, up to changing the sign of time. We shall further restrict the analysis to T + since the argument for T -is identical, up to some sign changes. Proposition 4.2. Let ϕ 0 be a function, microlocalized in a compact set C satisfying Assumption (2.4), and let ϕ be the associate solution of (4.1). Let [u -, u + ] be a closed interval of R containing Supp x 1 ū. There exists a constant t exit > 0 such that for any ε ∈]0, 1[, the function ϕ(εt exit , •) is microlocalized in a compact set K such that the projection of K onto the x 1 -axis does not interesect 3)

[u -, u + ]. Moreover µSupp ξ 1 ϕ is unchanged. More precisely, if µSupp ξ 1 ϕ 0 ⊂ R + \{0}, then µSupp x 1 ϕ(εt exit , •) ⊂]u + , +∞[, and if µSupp ξ 1 ϕ 0 ⊂ R -\ {0}, then µSupp x 1 ϕ(εt exit , •) ⊂] -∞, u -[.
E ∆ i[H, A]E ∆ ≥ θE ∆ .
Note that Assumptions (1 -3) can be replaced by the fact that [f (H), A] is bounded for any smooth, compactly supported function f (see [START_REF] Hunziker | Minimal escape velocities[END_REF]).

Under those assumptions, for any integer m ∈ N and for any θ ∈]0, θ[, there is a constant C such that

χ -(A -a -θ t)e -iHt g(H)χ + (A -a) ≤ Ct -m
where χ ± is the characteristic function of R ± , g is any smooth compactly supported function in ∆, and the above bound is uniform in a ∈ R. As pointed out in [START_REF] Hunziker | Minimal escape velocities[END_REF], this implies in particular that for any ϕ 0 in the image of E ∆ , the function e -iHt ϕ 0 has spectral support in [a + tθ , ∞[ with respect to A, up to t -∞ .

Let us apply this theory to our situation. We consider equation (4.1) with data e i t exit ε T + ϕ 0 + , and let us define the operator T 0 + as the operator T + where ū has been chosen identically zero. We shall start by studying the equation

(4.4) iε 2 ∂ t ϕ = T 0 + ϕ, ϕ |t=εt exit = e i t exit ε T + ϕ 0 + ,
for which we shall prove Proposition 4.3. Then we shall prove that the solution ϕ actually solves the original equation (4.1) with the same data e i t exit ε T + ϕ 0 + at t = εt exit up to O(ε ∞ ), because its support in x 1 lies outside the support of ū and because the symbolic expansion of T + depends on x 1 only through ū and its derivatives (see Remark 3.4). So let us start by applying Mourre's theory to (4.4). Let us write the projection of K onto the ξ 1 -axis as included in [d 0 , d 1 ] with 0 < d 0 < d 1 < ∞. We recall that on the support of e i t exit ε T + ϕ 0 + , x 1 remains to the right of the support of ū. Then we apply the theory to H = T 0 + and to A = x 1 (the pointwise multiplication). Assumptions (1) to (3) are easy to check, in particular because this is a semiclassical setting, so only the principal symbols need to be computed. Similarly finding a lower bound for

E ∆ i[T 0 + , x 1 ]E ∆ boils down to computing the Poisson bracket {τ + , x 1 } where {f, g} = ∇ ξ f • ∇ x g -∇ x f • ∇ ξ g,
and one finds (4.5)

{τ + , x 1 } = ξ 1 ξ 2 1 + ξ 2 2 + b 2 (x 2 )
• Since T 0 + has constant coefficients in x 1 , one can take the Fourier transform of (4.4) and ξ 1 is preserved, so in particular for all times one has µSupp ξ 1 ϕ(t) ⊂ [d 0 , d 1 ]. One can furthermore choose for ∆ an interval of R of the type ]D 0 , D 1 [ where the constants D 0 and D 1 are chosen so that for any (x, ξ) ∈ K, one has (4.6)

D 0 < ξ 2 1 + ξ 2 2 + b 2 (x 2 ) < D 1 .
As the microlocal supports of the eigenfunctions of T 0 + lie on energy surfaces, we know that the solution to (4.4) will remain in E ∆ for all times. Now let us apply the results of [START_REF] Jensen | Multiple commutator estimates and resolvent smoothness in quantum scattering theory[END_REF] and [START_REF] Hunziker | Minimal escape velocities[END_REF]. By Lemma 3.3, (4.5), (4.6) and the assumption on ξ 1 written above, we have that

E ∆ i[H, A]E ∆ ≥ ε d 0 D 0 E ∆ , so (
ε 2 T 0 + e i t exit ε T + ϕ 0 + = e i (t-εt exit ) ε 2 T + e i t exit ε T + ϕ 0 + in L 2
locally uniformly in t (see the comparison result in the Appendix).

The proposition follows.

5. Propagation of the Rossby waves 5.1. Semiclassical transport equations and microlocalization.

Because of the scaling of the Rossby hamiltonian (which is smaller than the Poincaré hamiltonians by one order of magnitude), on the times scales considered here, the propagation of energy by Rossby waves is described by the hamiltonian dynamics

dx i dt = ∂τ R ∂ξ i , dξ i dt = - ∂τ R ∂x i ,
which can be written explicitly (5.1)

dx 1 dt = b (x 2 ) ξ 2 b -2ξ 2 1 ξ 4 b + u 1 (x), dx 2 dt = -2b (x 2 ) ξ 1 ξ 2 ξ 4 b + u 2 (x), dξ 1 dt = -∂ 1 u 1 (x)ξ 1 -∂ 1 u 2 (x)ξ 2 , dξ 2 dt = ξ 1 2b(b ) 2 -b ξ b ξ 4 b -∂ 2 u 1 (x)ξ 1 -∂ 2 u 2 (x)ξ 2
where we recall that ξ b = ξ 2 1 + ξ 2 2 + b 2 (x 2 ) . In order for the dynamics to be well defined and also in order to justify the diagonalization process, we need the quantity ξ b to remain bounded from below for all times. Proposition 5.1. Let C be some compact subset of R 4 such that

C ∩ {(x 1 , x 2 , ξ 1 , ξ 2 ) / ξ 2 1 + ξ 2 2 + b 2 (x 2 ) = 0} = ∅ .
Then the bicharacteristics t → (x(t), ξ(t)) of the Rossby Hamiltonian starting from any point (x 0 1 , x 0 2 , ξ 0 1 , ξ 0 2 ) of C are defined globally in time, and ∀t ∈ R, inf

(x 0 1 ,x 0 2 ,ξ 0 1 ,ξ 0 2 )∈C (ξ 1 (t) 2 + ξ 2 (t) 2 + b 2 (x 2 (t)) > 0.
Proof. As b , b , u and Du are Lipschitz, by the Cauchy-Lipschitz theorem the system of ODEs (5.1) has a unique maximal solution. In order to prove that this solution is defined globally, it is enough to prove that the time derivative of this solution is uniformly bounded. This comes from assumption (2.3) giving an upper bound on b / 1 + b 2 (x 2 )

1 2 and b / 1 + b 2 (x 2 ) 1 2
, and from the lower bound on ξ b to be established now.

The crucial assumption here is the fact that b and b do not vanish simultaneously, and more precisely the existence of η, β > 0 such that for all x 2 ∈ R.

|b(x 2 )| < η ⇒ |b (x 2 )| ≥ β .
Along a trajectory of the Rossby Hamiltonian, τ R is conserved. Therefore,

• either b(x 2 ) ≥ η, and ξ 2 b > η 2 , • or b(x 2 ) < η. In this last case, if |ξ| ≥ η, ξ 2 b > η 2 . Else, |τ | ≥ η ξ b -ū ∞ η
from which we deduce that

ξ b ≥ η |τ | + ū ∞ η •
In any case, we have the lower bound Proof. • Let us start by describing the dynamics in the absence of zonal flow: ξ 1 is then an invariant of the motion, so that the dynamics in (x 2 , ξ 2 ) can be decoupled. Furthermore, as the energy surfaces are compact τ = ξ 1 b (x 2 ) ξ 2 1 + ξ 2 2 + b 2 (x 2 ) the motion along x 2 is periodic (with infinite period for homoclinic and heteroclinic orbits).

ξ b ≥ η min 1, 1 |τ | + ū ∞ η ,
The motion along x 1 is then determined by the equation

dx 1 dt = b (x 2 ) ξ 2 b -2ξ 2 1 ξ 4 b •
It is trapped if and only if the average of the right-hand side over one period is zero. Outside from saddle points, this quantity depends continuously on ξ 1 , so that we expect the initial data leading to trapped trajectories to belong to a manifold of codimension 1. This can be proved rigorously if b 2 has only one non degenerate critical points (see [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF]).

• Let us now turn to the influence of the zonal flow. We will first check that the only possible escape direction is again x 1 . Indeed the energy surfaces corresponding to τ = 0 are bounded in the x 2 direction : as

x 2 → ±∞, b (x 2 )ξ 1 ξ 2 b + ū(x) • ξ → 0 .
Consider now a trajectory on the energy level τ R = 0, and some point of this trajectory (y 1 , y 2 , ξ 1 , ξ 2 ) such that y 2 / ∈ Supp x 2 ū. One has b (y 2 )ξ 1 = 0 .

-If b (y 2 ) = 0, then

dx 1 dt = dx 2 dt = dξ 1 dt = dξ 2 dt = 0 .
The uniqueness in Cauchy-Lipschitz theorem implies then that the trajectory is nothing else than a fixed point, and therefore in particular is bounded.

-If ξ 1 = 0, then

dx 2 dt = dξ 1 dt = dξ 2 dt = 0 and dx 1 dt = b (y 2 ) ξ 2 2 + b 2 (y 2 ) -ξ 2 1 ξ 4 b ,
meaning that the trajectory is a uniform straight motion along x 1 . In particular, it is bounded in the x 2 -direction.

Finally, we conclude that trajectories on the energy level τ R = 0 are either trapped in the support Supp x 2 ū, or trivial in the x 2 -direction.

Perspectives.

As recalled in the introduction, it is generally believed that in the situation depicted in this paper (a flow around a large macroscopic current) Rossby waves are trapped. However due to the 2-dimensional setting (compared to the work in [START_REF] Cheverry | Semiclassical and spectral analysis of oceanic waves[END_REF]) the trapping in the x 1 direction seems difficult to prove.

One way to be convinced of the trapping phenomenon should be by implementing the dynamical system numerically.

It should be pointed out however that actually in order to get physically relevant predictions for the oceanic eddies, one should consider 3D models, or at least 2D models involving the influence of stratification. The methods presented here seem to be robust and should be extended to such complex models, up to again the study of the hamiltonian system describing the Rossby dynamics.

Lemma 3 . 3 .

 33 Let a and b two symbols. Then the principal symbol of Op W ε (a)Op W ε (b) is ab and its subprincipal symbol is 1 2i {a, b}.

Proof. 3 )

 3 Define the function ψ(s) := ϕ(εs). Then (4.1) reads (4.2) iε∂ s ψ = T + ψ, ψ |s=0 = ϕ 0 , and any result proved on ψ on [0, T ] will yield the same result for ϕ on [0, T ε]. Notice that (4.2) is written in a semi-classical setting, so by the propagation of the microsupport (The operator i[H, A] is bounded from below and closable, and the domain D(B 1 ) where iB 1 is its closure, contains D(H). More generally for all n ∈ N the operator i[iB n , A] is bounded from below and closable and the domain D(B n+1 ) of its closure iB n+1 contains D(H), and finally B n+1 extends to a bounded operator from D(H) to its dual. (4) There exists θ > 0 and an open interval ∆ of R such that if E ∆ is the corresponding spectral projection of H, then (4.

5 . 2 .

 52 which can be made uniform for initial data in C replacing |τ | by max C |τ |. Dynamics outside from the support of ū.Using the fact that ū has compact support, and simple properties of the Rossby dynamics in the absence of zonal flow, we can prove the following Proposition 5.2. Let C be some compact subset of R 4 such thatC ∩ {(x 1 , x 2 , ξ 1 , ξ 2 ) / ξ 2 1 + ξ 2 2 + b 2 (x 2 ) = 0} = ∅ .Then the bicharacteristics of the Rossby Hamiltonian starting from any point of C are bounded in x 2 : sup(x 0 1 ,x 0 2 ,ξ 0 1 ,ξ 0 2 )∈C |x 2 (t)| < ∞.

  [START_REF] Gill | Atmosphere-Ocean Dynamics[END_REF].3) holds with θ = εd 0 /D 0 . It follows that the solution e i (t-εt exit )

				ε 2	T 0 + e i t exit ε T + ϕ 0 + to (4.4)
	has a support in x 1 such that
		x 1 > u + +	d 0 D 0	t ε
	which proves the result for (4.4).
	Since µSupp x 1 e i (t-εt exit ) ε 2	T 0

+ e i t exit ε T + ϕ 0 + does not cross Supp x 1 ū, one has actually e i (t-εt exit )

Acknowledgements. The authors are grateful to J.-F. Bony and N. Burq for introducing them to Mourre estimates, and for interesting discussions. I. Gallagher and L. Saint-Raymond are partially supported by the French Ministry of Research grant ANR-08-BLAN-0301-01.

theorem (see for instance [START_REF] Martinez | An introduction to semiclassical and microlocal analysis[END_REF] and the references therein), the microsupport of ψ is propagated by the bicharacteristics, which are the integral curves of the principal symbol. Recall that the principal symbol of T + is

) and the bicharacteristics are given by the following set of ODEs:

Notice that τ + is independent of x 1 , so ξt 1 is identically zero and therefore ξ t 1 ≡ ξ 0 1 . So for all s ≥ 0, the microlocal support in ξ 1 of ψ(s) remains unchanged, and in particular is far from ξ 1 = 0. Moreover one has

• Now we recall that the bicharacteristic curves lie on energy surfaces, meaning that on each bicharacteristic, τ

2 ) is a constant on each bicharacteristic, so that for all times,

• If ξ 0 1 > 0, then x 1 is propagated to the right and eventually escapes to the right of the support in x 1 of ū, whereas if ξ 0 1 < 0, the converse (to the left) occurs. Note that since the support of the data is compact, the exit time is uniform. Proposition 4.2 is proved. 4.2.2. Long time behaviour. The aim of this paragraph is to prove the following result, which again is only proved for positive times for simplicity. Proposition 4.3. Under the assumptions of Proposition 4.2, let ϕ + be the solution of (4.1) associated with the data ϕ(εt exit , •). Then µSupp x 1 ϕ + (t) does not intersect Supp x 1 ū for t ≥ εt exit , and µSupp ξ 1 ϕ + (t) remains unchanged for t ≥ εt exit . Finally µSupp x 1 ϕ + (t) exits any compact set in x 1 in finite time.

Proof. Before going into the proof, we shall simplify the analysis by only studying the case of T + (the case T -is obtained by identical arguments), and we shall only deal with the case when the support in ξ 1 of the data lies in the positive half space. The other case is obtained similarly.

The proof is based on Mourre's theory which we shall now briefly recall, and we refer to [START_REF] Jensen | Multiple commutator estimates and resolvent smoothness in quantum scattering theory[END_REF] and [START_REF] Hunziker | Minimal escape velocities[END_REF] for all details. Let us consider two self-adjoint operators H and A on a Hilbert space H. We make the following assumptions:

(1) the intersection of the domains of A and H is dense in the domain D(H) of H .

(2) t → e itA maps D(H) to itself, and for all ϕ 0 ∈ D(H),

He itA ϕ 0 < ∞.

Appendix : a comparison result

For the sake of completeness, we state here the result which shows the stability of the propagation under a O(ε ∞ ) error on the propagator. This result has been implicitly used in the proof of the diagonalization when comparing A and T ± , T R , and in the proof of dispersion when comparing T ± and T 0 ± .

Proposition 5.3. Let A ε and Ãε be two pseudo-differential operators such that

Let φ be a solution to i∂ t φ + Ãε φ = 0 microlocalized in Ω, and ϕ be the solution to

with the same initial data. Then, for all N ∈ N,

Proof. The proof is based on a simple energy inequality and is completely straightforward. We have

This leads to ϕ(t) -φ(t) 2 L 2 (R d ) = O(ε ∞ )t , which concludes the proof.


Hélène Henon

Fabrice Anizon

Bruno Pfeiffer

Michelle Prudhomme

published or not. The documents may come

Synthesis of dipyrrolo [3,4-a:3,4-c]carbazole-1,3,4,6-tetraones bearing a sugar moiety.


Marek Bucki
email: marek.bucki@imag.fr

Mohammad Ali Nazari

Yohan Payan

Finite Element Speaker-Specific Face Model Generation for the Study of Speech Production

In situations where automatic mesh generation is unsuitable, the Finite Element (FE) mesh registration technique known as Mesh-Match-and-Repair (MMRep) is an interesting option to quickly create a subject-specific FE model by fitting a predefined template mesh onto the target organ. The irregular or poor quality elements produced by the elastic deformation are corrected by a "mesh reparation" procedure ensuring that the desired regularity and quality standards are met. Here, we further extend the MMRep capabilities and demonstrate the possibility of taking into account additional relevant anatomical features. We illustrate this approach with an example of speaker's face biomechanical model generation comprising face muscles insertions. While taking advantage of the a priori knowledge about tissues conveyed by the template model, this novel fast and automatic mesh registration technique makes it possible to achieve greater modeling realism by accurately representing the organ surface as well as inner anatomical or functional structures of interest.

Introduction

The physical speech signal is the consequence of the interaction between the mechanical structure of the articulatory apparatus (jaw, tongue, face, velum, vocal folds) and the forces that act on it. Part of these forces is the result of a central control while the other part depends on interactions between articulators and on peripheral feedback. Thus, to study the control of speech production from the articulatory or acoustic signals, some researchers have proposed to model separately the biomechanics of the articulators (with the anatomical, mechanical and dynamical properties of their tissues) and the neurophysiology of the control, [START_REF] Laboissière | The control of multi-muscle systems: Human jaw and hyoid movements[END_REF]], [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: a modeling study[END_REF]]. Some clues were then provided in terms of biomechanical models of the jaw ( [START_REF] Laboissière | The control of multi-muscle systems: Human jaw and hyoid movements[END_REF]], [START_REF] Hannam | A dynamic model of jaw and hyoid biomechanics during chewing[END_REF]]), tongue [START_REF] Payan | Synthesis of V-V sequences with a 2D biomechanical tongue model controlled by the Equilibrium Point Hypothesis[END_REF], [START_REF] Dang | Construction and control of a physiological articulatory model[END_REF]], [START_REF] Sanguineti | A control model of human tongue movements in speech[END_REF]], [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: a modeling study[END_REF]], [START_REF] Buchaillard | Muscle activations and the impact of gravity on tongue positioning during cardinal vowel production: a modeling study[END_REF]]), velum [START_REF] Payan | Physical modelling of the airflow-walls interactions to understand the sleep apnea syndrome[END_REF]]) or vocal folds ( [START_REF] Hunter | A three-dimensional model of vocal fold abduction/adduction[END_REF]], [START_REF] Titze | A two-dimensional biomechanical model of vocal fold posturing[END_REF]]).

As regards the face, excepting the pioneer works in computer graphics for facial animation, [START_REF] Lee | Realistic modeling for facial animationIn: Proceedings of SIGGRAPH[END_REF]], [START_REF] Sifakis | Automatic Determination of Facial Muscle Activations from Sparse Motion Capture Marker DataIn[END_REF]], and computer aided surgery, [START_REF] Chabanas | Patient specific finite element model of the face soft tissues for computer-assisted maxillofacial surgery[END_REF]], [START_REF] Gladilin | A framework for biomechanical simulation of cranio-maxillofacial surgery interventionsIn[END_REF]], few works were provided, to our knowledge, in terms of accurate biomechanical modeling of the facial soft tissues for speech production study, i.e. with a focus on facial muscles activations and their control for a specific speech task. Some modeling works were provided with a focus on lips modeling, [START_REF] Gomi | A physiologically based model of perioral dynamics for various lip deformations in speech articulation[END_REF]], [START_REF] Kim | Model-based investigation of control and dynamics in human articulatory motion[END_REF]], while [START_REF] Nazari | A continuous biomechanical model of the face: a study of muscle coordinations for speech lip gesturesIn[END_REF]] started to model facial tissues deformations due to muscles activations for speech gestures.

Once models such as the one manually developed by [START_REF] Nazari | A continuous biomechanical model of the face: a study of muscle coordinations for speech lip gesturesIn[END_REF]] on a generic single example are able to simulate speech production, they will have to be geometrically adapted to morphologies of various speakers in order to compare measured data (speakers face deformations and acoustic signals) with numerical simulations provided by the models. This paper addresses such a perspective.

In terms of modeling, the mesh adaptation step requires a geometrical fit of the generic model onto each speaker anatomy. The generic face model proposed by [START_REF] Nazari | A continuous biomechanical model of the face: a study of muscle coordinations for speech lip gesturesIn[END_REF]] was built from data (CT exam) collected on a maxillo-facial patient. A three layers Finite Element (FE) mesh, symmetrical with respect to the sagittal plane, was manually designed to fit the patient skin (external mesh surface) and skull (internal mesh surface) geometries reconstructed from CT images. Inside the 3D mesh, Nazari and colleagues have manually defined a set of muscular macrofibres which are implemented in reference to anatomical key points, independently of the mesh structure, by cable (tension only) elements, as shown in Fig. 1.

Our goal is to provide a robust automatic procedure that will match the Nazari et al. generic FE mesh onto data collected on given speakers, while fulfilling the following constraints:

1. Each generated mesh must be as close as possible to the speaker's anatomy reconstructed from CT and/or MRI data; to this end, high accuracy must be guaranteed between original data and mesh geometry, with specific focus on the lip region which plays a crucial role in acoustical outcome; 2. The topology of the mesh has to be maintained so that the three anatomical layers are still represented (epidermis, dermis and hypodermis, from the external topmost layer to the internal one respectively); 3. The symmetry of the mesh along the sagittal plane has to be maintained; 4. The location of the speaker muscular structures must be inferred from the registration, assuming that some clues can be manually provided such as muscles insertion points onto the skull.

Our group has already provided partial responses to points 1) and 2) by proposing the Mesh-Match-and-Repair (MMRep) algorithm. This procedure automatically computes an elastic 3D transform that matches a generic mesh onto the anatomy of a new subject, while maintaining the topology and the FE regularity of the mesh. It was developed in the context of computer aided surgery and evaluated on data acquired on maxillo-facial patients. This algorithm therefore did not focus specifically on the lips region. Section 2 describes the MMRep algorithm. Section 3 explains how we suggest to extend the MMRep method and proposes a "MMRep-for-Speech" algorithm compatible with the above mentioned constraints 3) and 4). Finally, Section 4 demonstrates our new technique on image data collected from a speaker.

2 The Mesh-Match-and-Repair Algorithm Finite Element mesh generation remains an important issue for subject specific biomechanical modeling. While some techniques, such as the widely used software Tetgen4 , make automatic mesh generation possible, in most cases, manual mesh generation is preferred for better control over the elements layout that it provides. Yet, this option is time consuming and requires interaction with an expert operator.

The Mesh-Match-and-Repair (MMRep) algorithm discussed in this paper is a twofold procedure. First, the subject data and the generic (the term "Atlas" is often used) mesh surface nodes are registered using the elastic deformation procedure described by [Bucki(2008)]. This deformation is then applied to the inner Atlas nodes yielding a FE mesh that represents the modeled domain with sufficient accuracy.

As a consequence of this deformation, the Atlas elements may suffer distortions and become either "irregular elements" which make FE analysis impossible, or "poor quality elements" in which case the computation, although feasible, can exhibit numerical instabilities. To recover mesh regularity and reach an acceptable quality level a fast and robust mesh reparation procedure has been proposed by [Lobos(2009)]. The goal of the mesh reparation procedure is to optimize mesh regularity and quality criteria formulated as relaxation energy terms.

The MMRep procedure has been tested on 10 clinical cases in the field of orthopedics and 50 cases in the context of orthognathic surgery where it proved capable of automatically generating patient specific hexahedral dominant quality meshes within a couple of minutes time, [START_REF] Bucki | A Fast and Robust Patient Specific Finite Element Mesh Registration Technique: Application to 60 Clinical Cases[END_REF]]. In this study, the patients' skin and skull surfaces were fully recovered from manually processed and segmented Computed Tomography (CT) scans, [Tilotta(2008)], and the face meshes registered so as to accurately represent both anatomical surfaces, as shown in Fig. 2.

The MMRep-for-Speech Algorithm

The MMRep-for-Speech algorithm proposed in this article provides the following functions: 1. In order to improve the feature correspondence computed by the elastic registration, the ability to constrain the deformation with a set of control points which positions are known both in Atlas and speaker space, has been added to the mesh registration procedure ( §3.1); 2. The elastic registration procedure which initially produced asymmetrical face meshes, as it can be seen on Fig. 2-b, has been further adapted to suit the prerequisites of speech production modeling by enforcing a symmetrical deformation field with reference to the central sagittal subject plane, while maintaining the regularity which makes it suitable for our mesh deformation procedure ( §3.2); 3. As speech production simulation requires guaranteed mesh registration accuracy in specific anatomical areas such as the speaker's lips region, a local registration refinement procedure has been implemented which strongly reduces the residual registration error ( §3.3); 4. The face muscles that were integrated into the generic Atlas mesh proposed by [START_REF] Nazari | A continuous biomechanical model of the face: a study of muscle coordinations for speech lip gesturesIn[END_REF]] were automatically deformed by the computed elastic registration, therefore inferring the subject's face muscles location ( §3.4).

Control points driven registration

We define an elastic registration function as a mapping R : R 3 → R 3 that superimposes a set of Atlas mesh nodes A onto a speaker specific target data set S, which can either be a point cloud or a surface mesh. The computed elastic registration procedure complies with continuum mechanics conditions of motion, [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]], as R defines a C 1 -diffeomorphic, one-toone and non-folding correspondence between geometric data sets, [Bucki(2008)], [START_REF] Bucki | A Fast and Robust Patient Specific Finite Element Mesh Registration Technique: Application to 60 Clinical Cases[END_REF]].

The iterative assembly of R is driven by the minimization of a "registration energy" E reg which measures the similarity between the deformed Atlas point set and the target data set S. As geometrical shape similarity is sought, E reg is defined as the sum of Euclidean distances between deformed A and S :

E reg (R ) = a∈A d(R(a), S ) (1)
To increase registration performance, a distance map ( [START_REF] Saito | New algorithms for Euclidean distance transformations of an n-dimensional digitized picture with applications[END_REF]]) is computed based on the structure defined by S, and used for the evaluation of E reg at each registration iteration.

In situations where the modeled organ geometry does not provide enough constraints for the registration to converge properly, for subjects with soft features, or where the face mesh is poorly reconstructed, the mesh deformation algorithm must be provided with additional data. This information can be supplied as a set of "control points" automatically or manually identified within both Atlas, āi , i ∈ [1, C], and subject's space, si , i ∈ [1, C]. Any quantity of control points can be defined, but in our experience, a small number (C = 5) of key features seems sufficient.

To be effective, the influence of the small set of control points within the overall registration energy must be magnified by multiplicative "weight" factors

W i , i ∈ [1, C],
representing the strength of the attraction between the deformed control points R(ā i ) and their desired positions si in subject space.

As the elastic deformation is being iteratively computed, the value of the weight factors for control points presenting inaccurate registration is increased until the desired accuracy standards are met. Within this new framework, the minimized registration energy defined in Eq. 1 is rewritten as:

E reg (R ) = a∈A d(R(a), S ) + i∈[1,C] W i d(R(ā i ), si ) (2)
In order to provide greater control over features positioning during mesh registration, 5 skin control points were defined on half of Atlas and subject face prior to the symmetrical registration procedure described in §3.2. Fig. 3 displays an example of segmented half speaker's face along with the 5 chosen control points defined in Atlas space (dark gray disks) and subject space (light gray disks). The dark dots represent the Atlas FE mesh nodes that need to be registered onto the subject's skin and skull surfaces. The skull is not represented on the figure, for sake of clarity.

Symmetry enforcement

Speech production simulation run times can be greatly reduced by assuming sagittal symmetry and performing FE analysis on half of the face mesh. To this end, the registration must be constrained on the sagittal mid-plane so that the extension to the right half of the model of the deformation computed on the left half retains the continuity and smoothness properties mentioned above. This aspect differs from our previous application of the MMRep procedure in the context of orthognathic surgery where pathological patient morphology was essentially asymmetrical. The first step in the symmetrical registration process is the proper identification of the pseudo-symmetry5 plane. The medical image volume sagittal plane is not necessarily a good candidate as the speaker may have tilted his or her head before the CT volume acquisition. A practical solution to this problem is proposed by [START_REF] Chang | Reverse engineering of a symmetric object[END_REF]] and relies on rigid registration of the segmented skin surface onto itself after application of a rough planar symmetry (image sagittal mid-plane) and subsequent analysis of the correspondences between the surface points established by the rigid registration.

Once the sagittal pseudo-symmetry plane has been estimated, the elastic deformation R is computed as follows. The initial Atlas point set is embedded within a deformable "virtual elastic grid". We arbitrarily set the shape of the grid to be the bounding box of the points, extended by a 10% margin. The considered deformation R is then formed by combinations of elementary deformations of the grid, r n , n ∈ [1, N ], all having the desired regularity properties, much in the same way as proposed by [START_REF] Rueckert | Diffeomorphic registration using B-splinesIn: Proceedings of Medical Image Computing and Computer-Assisted Intervention MICCAI[END_REF]].

Each elementary deformation step n is performed by moving a single grid node and propagating the local deformation to the neighboring grid cells. The nodal displacement is transferred to the Atlas points located in the affected cells using a C 1 -diffeomorphic shape function ϕ : R 3 → [0, 1] taking the value 1 at the displaced node and 0 at the boundary and outside the cells surrounding this node. The computation of each r n is driven by the minimization of the registration energy defined in Eq. 2. As the registration progresses, the virtual grid is subdivided into smaller cells so that correspondence for finer anatomical details can be computed, [START_REF] Bucki | A Fast and Robust Patient Specific Finite Element Mesh Registration Technique: Application to 60 Clinical Cases[END_REF]]. Thus, the expression of the overall registration function is:

R = r N • . . . • r 1
Fig. 4 illustrates in 2D our multi-grid iterative registration technique. The Atlas points A are the gray dots and S is not represented. In (a) the initial set A 0 := A is embedded within the virtual grid (black square) discretized at the coarsest level 1 and the gradients of the energy, as defined in Eq. 2, are computed at the 4 mesh nodes. In (b) the nodal displacement yielding the highest energy decrease is applied to the grid and the embedded Atlas points move accordingly giving a new Atlas point set A 1 . This operation is repeated i times until no significant energy decrease can be generated at level 1. The grid mesh is then refined (c) and A i is embedded within it. The energy gradients are computed at the 9 grid nodes, the best nodal displacement is applied (d) and the resulting Atlas point set A i+1 is computed. The registration iterations stop when a user defined highest grid refinement level is reached.

The enforcement of symmetrical deformation in this framework is quite straightforward. At each virtual grid refinement level, the grid nodes located on the face sagittal mid-plane are identified and their movements are constrained within the sagittal plane. As a consequence, the sagittal symmetry plane remains invariant for all considered elementary registration functions r n , n ∈ [1, N ]. Furthermore, thanks to the diffeomorphic nature of ϕ all continuity and smoothness properties are maintained throughout the sagittal mid-plane after the left half deformation has been transferred to the right side of the model.

Local registration refinement

The registration accuracy of the MMRep algorithm is strongly related to the resolution of the distance map used to speed-up the evaluations of E reg at each step, which enables FE mesh generation within seconds. This resolution, limited by practical considerations such as distance map computation time and memory occupation, constrains the size of the smallest virtual grid cell reached during the top-down grid refinement process. Indeed, performing registration on a finer grid than the represented target set S is meaningless and leads to unnecessary local space distortions.

Yet, in the context of a speech production study, a submillimetric domain representation accuracy must be guaranteed in strategic areas such as the speaker's lips region. To reduce the residual surface representation error in these regions after the initial deformation R has been applied, we propose to compute a finer local registration function Q based on more detailed lips geometry information. This local registration adjustment relies on a finer distance map and hence, enables smaller virtual grid cells to be used, while requiring acceptable computational overhead.

To maintain the consistency of the resulting elastic registration Q • R, continuity and smoothness must be preserved not only on the local sub-domain subject to refinement but also at the overall face scale. Again, this is accomplished by controlling the displacements of the virtual grid nodes at each grid refinement level. Fig. 5 illustrates our local elastic registration refinement strategy. First, the overall registration R is computed over the raw virtual grid (black). Then, a finer distance map is used to increase the surface representation accuracy in the region surrounded by the finer grid (dashed) while the points located outside remain unchanged. The continuity and smoothness of Q•R through the boundary between both grids is preserved by considering a registration function Q that leaves all its underlying grid boundary nodes fixed. The definition of the shape functions ϕ used to transfer the virtual grid nodes displacements to the enclosed Atlas mesh nodes ensures that the overall registration function Q • R retains the smoothness properties discussed in §3.1.

Speaker muscles location inferred by the elastic registration

Muscle courses are identified from CT data and anatomical facts. Muscle insertion points can be recognized from the landmarks on the skull surface. Muscles meet points can also be used to define muscle tracks, especially for those which do not have skull insertion points. Since muscles appear as fiber bundles, the course of these fibers can be modeled by a few linear piecewise continuous tracks.

The number of tracks per muscle depends on the extent of each muscle. For example in orbicularis oris two tracks are attributed for peripheralis and marginalis parts. Each track is characterized by a series of points which are picked up from the CT data based on anatomical facts and relying on surgeon's expertise, [START_REF] Nazari | A continuous biomechanical model of the face: a study of muscle coordinations for speech lip gesturesIn[END_REF]].

In our approach, the muscle tracks identified in the Atlas mesh by an expert operator are deformed by the elastic registration function computed as described above. The muscle tracks structure is preserved by the deformation while the inferred subject specific insertions points match the considered speaker's anatomy.

Results

In order to check its applicability in situations where less attention is paid to the regularity of the reconstructed subject anatomical structures, the MMRep-for-Speech algorithm has been tested on raw clinical data, as opposed to the manually processed skull and face meshes used in our previous study, [START_REF] Bucki | A Fast and Robust Patient Specific Finite Element Mesh Registration Technique: Application to 60 Clinical Cases[END_REF]]. Popular surface segmentation and reconstruction techniques often yield uneven surface meshes, possibly exhibiting segmentation artifacts. Automatic clean-up is not always satisfactory and a manual time consuming post-processing is required in most cases. As it is usually difficult to achieve smooth surface reconstruction, FE mesh registration procedures must be able to deal with this limitation.

The new algorithm was tested on a sample speaker data set. Bone and skin surfaces were automatically extracted from the CT volume using the Marching Cubes method, [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithmIn[END_REF]]. The segmented and reconstructed skull surface comprises 263,994 vertices and 421,410 triangles, while the skin model has 128,011 vertices and 253,654 triangles. Both skin and skull surfaces are taken into account in the subject specific FE mesh produced by the MMRep algorithm, as illustrated in Fig. 6. After overall mesh registration R and subsequent reparation, the surface representation errors measured as Atlas mesh nodes-to-surface distances are given in table 1. Mesh registration and reparation time was less than one minute. The registration was computed using a point-to-point distance map with a resolution of 1.5mm but the errors reported below are measured using the actual triangular surface meshes. An overall submillimetric accuracy was achieved and the maximal errors are mainly due, for the skull, to an approximative representation of the temporal bones in the Atlas, and, for the skin, to the absence of ears in the generic mesh, as can be seen in Fig. 6-b. 1. Mesh registration errors, in millimeters, measured as distances between the Atlas mesh skull or skin nodes and their corresponding anatomical surface (mean, maximal error and standard deviation).

As for the face control points registration accuracy, the errors measured on the set of 5 features presented in Fig. 3 are: mean err. = 0.7 mm, max. err. = 1.0 mm and σ = 0.2 mm. These figures could be improved by using a distance map with a greater resolution or stronger control points weights. However this seems unnecessary as the purpose of control points is merely to ensure that proper features geometry is reconstructed by the registration process. Although the Atlas control points near the lips do not perfectly overlap their subject specific counterparts, it is clear that the lips geometry has been perfectly captured by the registration procedure, as Fig. 6-c andd show.

Within the Atlas mesh, a subset of 350 skin nodes has been marked as "lips nodes". These nodes were used, first, to measure the mesh registration accuracy in this region and, second, to define a sub-region where registration refinement would take place, should the achieved accuracy be unsatisfactory.

Table 2 shows the lips registration accuracy after overall half face registration (i.e. without local registration refinement), after local registration refinement using a 0.5mm resolution distance map and finally after reparation of the locally registered mesh.

The lips surface representation accuracy is higher after local registration refinement, with a maximal error of 0.15 mm. Nevertheless, as the resulting mesh needs to undergo the reparation procedure, i.e. regularity recovery and quality optimization, some lips surface nodes positions are altered resulting in a slight accuracy loss, shown in the bottom line of table 2. Nevertheless, the resulting FE mesh presents an overall submillimetric mean accuracy and, in the particular lips region, a maximal surface representation error of 0.29 mm which is an acceptable accuracy standard for speech production modeling.

As stated initially, the Atlas registration procedure was able to recover an approximate localization of the speaker's face muscles, as illustrated in Fig. 7-b. Unfortunately, in the absence of anatomical information about the position of the actual muscle insertions, no quantitative comparison could be made between the positions of these anatomical landmarks and those inferred by the elastic registration algorithm. 2. Lips registration errors, in millimeters, after overall half face registration, refined lips registration (before mesh reparation) and reparation of the refined mesh registration. Finally the produced mesh could be used to simulate the contraction of two muscles on the subject's face. The effects of the contractions of zygomaticus and orbicularis oris peripheralis are respectively shown in Figures 8 and9. As can be observed, the activation of zygomaticus draws the corner of the mouth upward and laterally and activation of orbicularis oris peripheralis shows lips protrusion which qualitatively fits the anatomical facts even if the magnitude of protrusion seems too small.

Conclusions and Discussion

A series of extensions of MMRep, our fast and robust mesh generation technique, has been proposed and successfully tested on an illustrative case reproducing the constraints likely to be encountered in a speech production study.

The proposed algorithm offers an interesting balance between manual FE mesh generation with fully controlled element features and automatic "from scratch" mesh generation, by taking as starting point a user defined, tailor made generic mesh. Our procedure is capable of automatically generating within minutes a subject specific mesh directly exploitable by a commercial FE software such as ANSYS Workbench (ANSYS Inc., USA). Among the proposed extensions, only the definition of the control points in both Atlas and speaker space required user interaction. Yet the 5 control points used here could easily be identified by manual selection in less than a minute. If necessary, for these simple face features an automatic extraction algorithm may also be envisaged.

Future work includes enhancement of the mesh reparation procedure in order to take into account the surface representation constraint and avoid accuracy loss such as that shown in table 2. The current algorithm considers all nodes indistinctly and the search for an acceptable mesh configuration is performed by moving in priority the nodes presenting the worst regularity or quality conditions. A first enhancement would be to begin this search by moving the inner mesh nodes whose locations do not affect the surface representation accuracy. Furthermore, the displacements of surface nodes could also be considered, but should be constrained by the previously recovered subject surface or its vicinity.

The elastic registration algorithm used here preserves a certain regularity of space but can also reproduce finer details when the subject data are irregular, for morphological reasons or due to the presence of noise. This capability explains the jagged aspect of some muscles fibers in Fig. 7-b as opposed to the regular shapes featured in Fig. 7-a. Indeed, for the poorly reconstructed speaker data presented here, a point-to-point distance measure was used to evaluate the registration energy, mainly to reduce computation times. A point-to-surface measure would probably define a smoother energy landscape and contribute to preserving muscle regularity throughout the registration process.

The muscle deformation regularity could be further enhanced by considering during the registration computation a mechanical regularization term, as suggested by [START_REF] Bucki | A Fast and Robust Patient Specific Finite Element Mesh Registration Technique: Application to 60 Clinical Cases[END_REF]], which would specifically reflect the bending of the muscle fibers. As a consequence, the muscles network embedded within the face mesh would behave like a wire-frame structure preserving the surrounding tissues as well as its own smoothness.

Finally, from an anatomical point of view, a validation of the position of muscle insertions inferred from the registration procedure must be carried out based on actual muscle insertion coordinates in speaker space. If it appeared to be insufficiently accurate, then these insertion points could be used as registration control points, which should lead, as suggested above, to realistic muscle positioning within the reconstructed subject specific mesh and more accurate biomechanical modeling of the face's soft tissues.

Fig. 1 .

 1 Fig. 1. (a) Generic FE mesh of the face soft tissues. (b) Interactive segmentation of muscles fibers (arrows) on CT data. (c) Location of some muscles involved in facial mimics (represented for the left side of the skull only).

Fig. 2 .

 2 Fig. 2. MMRep face registration. (a) Segmented skull and skin surfaces from a sample CT volume. (b) Registered FE mesh.

Fig. 3 .

 3 Fig. 3. Face registration control points in Atlas space (dark gray disks) and subject space (light gray disks): inner and outer corner of the left eye, left corner of the lips, upper and lower center of the lips. The dark dots represent the Atlas FE mesh nodes prior to registration. The light gray surface is the CT segmented subject's skin, the reconstructed skull surface is not represented here. (a) Anterior view. (b) Sagittal view.

Fig. 4 .

 4 Fig. 4. Elastic registration overview. (a) A0 at refinement level 1. (b) A1 after deformation at level 1. (c) Ai at refinement level 2. (d) Ai+1 after deformation at level 2.

Fig. 5 .

 5 Fig. 5. Registration refinement. Black grid: overall half face registration R. Dashed grid: accurate lips registration Q. The boundary between the overall deformation domain and the lips sub-domain is fixed.

Fig. 6 .

 6 Fig. 6. Sample speaker data. (a) Bone and skin surfaces reconstructed from CT. (b) FE mesh produced by the MMRep-for-Speech procedure. (c) Sagittal close-up view on the lips region. (d) Anterior close-up view on the lips region.

Fig. 7 .

 7 Fig. 7. (a) Muscles layout within the generic Atlas mesh, shown in transparency. (b) Muscles location inferred by the elastic deformation computed between the Atlas mesh and the considered subject's face.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Effects of zygomaticus muscle (white arrow) contraction on the face mesh. (a) Rest position. (b) Deformed face.

http://tetgen.berlios.de

As subject data is not strictly speaking symmetrical.

Acknowledgments

The authors wish to thank Mr. Badin for his remarkable anatomy.