S Benzoni-Gavage

M D Rosini

Weakly nonlinear surface waves and subsonic phase boundaries

Keywords: 2000 Mathematics Subject Classification: 35C20, 35L50, 35L67, 35R35, 76T10 Amplitude equation, nonlocal Burgers equation, subsonic phase transitions

The aim of this work is twofold. In a first, abstract part, it is shown how to derive an asymptotic equation for the amplitude of weakly nonlinear surface waves associated with neutrally stable undercompressive shocks. The amplitude equation obtained is a nonlocal generalization of Burgers' equation, for which an explicit stability condition is exhibited. This is an extension of earlier results by J. Hunter. The second part is devoted to 'ideal' subsonic phase boundaries, which were shown by the first author to be associated with linear surface waves. The amplitude equation for corresponding weakly nonlinear surface waves is calculated explicitly and the stability condition is investigated analytically and numerically.

Introduction

This work is concerned with the multi-dimensional theory of -possibly nonclassical -shock waves that are neutrally stable, which means that their linearized stability analysis yields neutral normal modes. More specifically, we are interested in cases when these neutral modes are of finite energy, that is, when these modes are (genuine) surface waves. A program initiated by Hunter [START_REF] Hunter | Nonlinear surface waves[END_REF] has shown that surface waves are usually associated, in the weakly nonlinear regime, to amplitude equations that are nonlocal generalizations of Burgers' equation. Our main purpose is to apply this program in the framework of 'shocks', including undercompressive ones, with application to phase boundaries. Indeed, it was shown in [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF] that nondissipative, dynamic and subsonic phase boundaries in van der Waals-like fluids are neutrally stable, with surface waves (also see [START_REF] Benzoni-Gavage | Effects of surface tension on the stability of dynamical liquid-vapor interfaces[END_REF]). The present paper contains two main parts. In the first one we derive the amplitude equation associated with surface waves along neutrally stable shocks in an abstract framework, and give an alternative version of Hunter's stability condition that is easy to check in practice. In the second part we perform the computations in the explicit case of surface waves along dynamic subsonic phase boundaries. It turns out, as our numerical results show, that Hunter's stability condition is not satisfied by the amplitude equation associated with subsonic phase boundaries. This is in contrast with what happens in Elasticity for instance, where the amplitude equation associated with Rayleigh waves is known to satisfy Hunter's condition [START_REF] Hunter | Nonlinear surface waves[END_REF][START_REF] Parker | Waveform evolution for nonlinear surface acoustic waves[END_REF][START_REF] Parker | Analysis and computation for nonlinear elastic surface waves of permanent form[END_REF].

Derivation of the amplitude equation 2.1 General framework

We consider a hyperbolic system of conservation laws

d i=0 ∂ i f i (u) = 0 n , x ∈ R d , t > 0 , (2.1)
where the unknown is u = t (u 1 , . . . , u n ) : (t, x) ∈ [0, ∞)×R d → u(t, x) ∈ R n , ∂ 0 stands for the partial derivative with respect to t and ∂ i denotes the partial derivative with respect to x i , i = 1, . . . , d. Here, f i = t (f i 1 , . . . , f i n ) : u ∈ U → f i (u) ∈ R n , i = 0, . . . , d, are given smooth fluxes (at least C 2 ) on an open subset U of R n . We shall denote by A i := (∂ u k f i j ) 1≤j,k≤n : u ∈ U → A i (u) ∈ R n×n the Jacobian matrix of f i , i = 0, . . . , d, and assume that:

• for all u ∈ U, the matrix A 0 (u) is nonsingular,

• for all u ∈ U and all η ∈ R d \ {0 d }, the matrix A 0 (u) -1 d i=1 η i A i (u), has n real eigenvalues λ 1 (u, η) ≤ λ 2 (u, η) ≤ . . . ≤ λ n (u, η) and n linearly independent corresponding eigenvectors r 1 (u, η), . . . , r n (u, η) ∈ R n , i.e. d i=1 η i A i (u)λ j A 0 (u) r j = 0 n , j = 1, . . . , n .

We are concerned here with special, shock-like weak solutions to (2.1) that are C 1 outside a smooth moving interface. Recall that for a hypersurface Σ :

= (t, x) ∈ [0, ∞) × R d : Φ(t, x) = 0 , (2.2)
where Φ : [0, ∞) × R d → R is a C 1 function, a mapping u : (0, T ) × R d → R n that is C 1 on either side of Σ is a weak solution of (2.1), if and only if,

d i=0 ∂ i f i (u ± ) = 0 n , ±Φ(t, x) > 0 , t ∈ (0, T ) ,
where u ± is the restriction of u to the domain

Ω ± := (t, x) ∈ [0, ∞) × R d : ± Φ(t, x) > 0 ,
and

d i=0 f i (u) ∂ i Φ = 0 n , Φ(t, x) = 0 , (2.3)
where the brackets [ • ] give the 'strength' of the jump across the interface.

It is well-known that the Rankine-Hugoniot jump conditions in (2.3) are not sufficient in general to ensure uniqueness of weak solutions. They must be supplemented with admissibility conditions. For "classical" shocks, standard admissibility conditions are given by the Lax inequalities (see for instance [START_REF] Serre | I. Hyperbolicity, entropies, shock waves[END_REF]), which require that the number of characteristics outgoing the shockfront is less than the number of incoming characteristics. More precisely, for Laxian shocks in a states space of dimension n, the number of outgoing characteristics is n -1 and the number of incoming ones is n + 1. For "nonclassical" shocks, the situation is different, and in particular for undercompressive ones, the number r of incoming characteristics is less or equal to n. Then a number p := n + 1r of additional jump conditions is needed. In what follows, we consider undercompressive shocks for which these additional jump conditions can be written as

d i=0 g i (u) ∂ i Φ = 0 p , Φ(t, x) = 0 , (2.4)
where g i : u ∈ U → g i (u) ∈ R p , i = 0, . . . , d, are smooth (at least C 2 ). For both Laxian and undercompressive shocks, the resulting system is

           d i=0 ∂ i f i (u) = 0 n , Φ(t, x) = 0 , d i=0 f i (u) ∂ i Φ = 0 n+p , Φ(t, x) = 0 , (2.5)
where

• for classical shocks: p = 0, f i (u) := f i (u) ∈ R n ;
• for undercompressive shocks:

p ≥ 1, f i (u) := f i (u) g i (u) ∈ R n+p .
This work is motivated by nondissipative subsonic phase boundaries in van der Waals fluids, which can viewed as undercompressive shocks with p = 1 and (2.4) given by the so-called capillarity criterion [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF]. More precisely, for isothermal phase boundaries, the interior equations are given by the conservation of mass and of momentum -with a non-monotone pressure law ρ → p(ρ) -and the additional jump condition is given by the conservation of total energy (in fact, the free energy plus the kinetic energy).

Our purpose is to describe nontrivial approximate solutions to the fully nonlinear problem (2.5). The starting point will be a planar stationary noncharacteristic shock-like solution.

Assumption 1 There exists u = (u l , u r ) ∈ R n × R n such that    u(t, x) := u l , x d < 0 , u r , x d > 0 , Φ(t, x) := x d ,
is a solution of the nonlinear problem (2.5). In addition, we assume that the matrices A d (u l ) and A d (u r ) are nonsingular.

The linearized problem

We are interested in solutions of (2.5) close to the planar stationary solution u given by Assumption 1. In this respect, we shall concentrate on solutions (v, Ψ) for which the location of the shock front is given by an equation of the form

Ψ(t, x) = 0 , where Ψ(t, x) = x d -χ(t, x 1 , . . . , x d-1 ) for a smooth map χ : (t, x 1 , . . . , x d-1 ) ∈ [0, ∞)×R d-1 → χ(t, x 1 , . . . , x d-1 ) ∈ R. Then the system (2.5) applied to (v, Ψ) instead of (u, Φ) becomes            d i=0 ∂ i f i (v) = 0 n , x d = χ(t, x 1 , . . . , x d-1 ) , d-1 i=0 f i (v) ∂ i χ = f d (v) , x d = χ(t, x 1 , . . . , x d-1 ) , (2.6)
where

[ f i (v)] := f i (v r ) -f i (v l ) ∈ R n+p , being v l (t, x 1 , . . . , x d-1 ) := lim x d րχ(t,x 1 ,...,x d-1 ) v(t, x 1 , . . . , x d ) , v r (t, x 1 , . . . , x d-1 ) := lim x d ցχ(t,x 1 ,...,x d-1 ) v(t, x 1 , . . . , x d ) .
As usual for free boundary value problems, we start by making a change of variables that leads to a problem in a fixed domain. Introducing the new unknowns v

± : [0, ∞) × R d-1 × [0, ∞) → R n , related to v by v ± (y 0 , y 1 , . . . , y d ) := v(y 0 , . . . , y d-1 , χ(y 0 , . . . , y d-1 ) ± y d ) ,
and redefining v as

v := (v -, v + ) : [0, ∞) × R d-1 × [0, ∞) → R 2n ,
we are led to the boundary value problem

L(v, ∇χ) • v = 0 2n , y d > 0 , b(v, ∇χ) = 0 n+p , y d = 0 , (2.7)
with

L(v, ∇χ) := d-1 i=0 A i (v)∂ y i + A d (v, ∇χ)∂ y d b(v, ∇χ) := d-1 i=0 (∂ i χ) f i (v) -f d (v) ∈ R n+p ,
where, for i = 0, . . . , d,

A i (v) := A i (v -) 0 n×n 0 n×n A i (v + ) , Ȋ2n := -I n 0 n×n 0 n×n I n Ȃi (v) := Ȋ2n A i (v) , A d (v, ∇χ) := Ȃd (v) - d-1 i=0 (∂ i χ) Ȃi (v) ∈ R 2n×2n ,
and, for convenience, the brackets [v] now stand for v +v -for any vector v = (v -, v + ) ∈ R 2n (or even C2n).

Using an observation of Métivier [START_REF] Métivier | Stability of multidimensional shocks[END_REF], we may simplify the boundary conditions in (2.7), at least for solutions close to u, provided that the following assumption holds true.

Assumption 2 The jump vectors

[ f 0 (u)], . . . , [ f d-1 (u)] are independent in R n+p .
Under Assumption 2, there exist a neighborhood V ⊆ U × U of u and a map

Q : v ∈ V → Q(v) ∈ GL n+p (R) (the group of non-singular (n + p) × (n + p) matrices) such that Q(v) d-1 i=0 ξ i+1 f i (v) = ξ 0 n+p-d ∈ R n+p for all ξ ∈ R d and all v ∈ V.
Therefore, the boundary value problem (2.7) can be rewritten with "simpler" boundary conditions:

L(v, ∇χ) • v = 0 2n , y d > 0 , J∇χ + h(v) = 0 n+p , y d = 0 , (2.8)
where

J := I d 0 (n+p-d)×d ∈ R (n+p)×d , h(v) := -Q(v) f d (v) ∈ R n+p .
The linearization of the simplified problem (2.8) about its constant solution (v ≡ u, χ ≡ 0) readily gives the equations for the perturbations v and χ of v and χ respectively,

L(u, 0) • v = 0 2n , y d > 0 , J∇ χ + H(u) • v = 0 n+p , y d = 0 , (2.9)
where H(u) ∈ R (n+p)×2n denotes the Jacobian matrix of h at u.

We shall now make a further assumption regarding the solutions of (2.9) that go to zero as y d goes to +∞. First of all, we introduce, for η =

(η 0 , η 1 , . . . , η d-1 ) ∈ R × R d-1 , the operator L(u, iη) := A(u, iη) + Ȃd (u) ∂ y d , with A(u, iη) := d-1 k=0 i η k A k (u) ,
obtained from L(u, 0) by Fourier transform in the tangential variable y = (y 0 , y 1 , . . . , y d-1 ). Observe that by the noncharacteristicity of the shock u (Assumption 1), the (2n × 2n) block-diagonal matrix Ȃd (u) is nonsingular. In what follows, we also use the notation L(u, iη) for vectors η for which η 0 = -iτ ∈ C, Re(τ ) > 0, the operator L(u, τ, iη 1 , • • • , iη d-1 ) arising when we perform a Laplace transform in y 0 instead of a Fourier transform. The hyperbolicity of (2.1) implies, by a classical observation due to Hersh [START_REF] Hersh | Mixed problems in several variables[END_REF], that for all η = (η 0 , η 1 , . . . , η d-1 ) ∈ C × R d-1 with Im(η 0 ) < 0, the matrix

A(u, η) := -Ȃd (u) -1 A(u, iη) (2.10)
is hyperbolic, that is, has no purely imaginary spectrum. It is well known that the well-posedness of the linear problem (2.9) crucially depends on the properties of the invariant subspaces of A(u, η). The following is a natural generalization of the Lopatinskiȋ condition to undercompressive shocks [START_REF] Freistühler | Some results on the stability of non-classical shock waves[END_REF] (regarding Laxian shocks, see the seminal work by [START_REF] Majda | Compressible fluid flow and systems of conservation laws in several space variables[END_REF] ).

Assumption 3 For all η = (η 0 , η 1 , . . . , η d-1 ) ∈ C × R d-1 with Im(η 0 ) < 0, the stable subspace E s (u, η) of A(u, η) is of dimension q := n + p -1, and there is no nontrivial (X, V ) ∈ C × E s (u, η) such that XJη + H(u)V = 0 n+p . (2.11)
Assumption 3 is known to be necessary for the well-posedness of (2.9) associated with suitable initial data. To investigate the actual well-posedness of this initial-boundary-value problem we need to go further and consider the subspace E(u, η) obtained as

E(u, η) := lim b ր 0 E s (u, η 0 + ib, η 1 , . . . , η d-1 )
(in the Grassmannian of q-dimensional subspaces of C 2n ). As the hyperbolicity of the matrix A(u, η) fails in general for real η 0 , the limiting space E(u, η) decomposes as

E(u, η) = E -(u, η) ⊕ E 0 (u, η) ,
where E -(u, η) is the (genuine) stable subspace of A(u, η), of dimension say m ≤ q, and E 0 (u, η) is a subspace of the center subspace of A(u, η).

Assumption 4 There exists η = (η 0 , η 1 , . . . , η d-1 ) ∈ R d , η 0 = 0, and

(X η , V η ) ∈ C × E(u, η) such that { (X, V ) ∈ C × E(u, η) ; XJη + H(u)V = 0 } = C {(X η , V η )} ,
and the vector V η belongs to E -(u, η)\{0 2n }.

Assumption 4 means that (2.5) admits surface waves, that is, solutions that are exponentially decaying in y d and oscillating in y = (y 0 , y 1 , . . . , y d-1 ). As observed in [START_REF] Benzoni-Gavage | Multi-dimensional hyperbolic partial differential equations: First-order systems and applications[END_REF]Chap. 7], even though surface waves signal a failure of the so-called uniform Kreiss-Lopatinskiȋ condition, their existence is still compatible with the well-posedness of constant-coefficients linear homogeneous boundary value problems, such as (2.9). For non-linear problems, the resolution of which relies on non-homogeneous linear problems, surface waves are responsible for a loss of regularity, see in particular the work of Coulombel and Secchi [START_REF] Coulombel | Nonlinear compressible vortex sheets in two space dimensions[END_REF].

Our purpose here is to adapt the method proposed by Hunter [START_REF] Hunter | Nonlinear surface waves[END_REF] to derive an amplitude equation for weakly non-linear surface waves associated with weakly stable shocks -i.e. shocks satisfying in particular Assumption 4.

Finally, we shall assume that frequencies of surface waves do not correspond to 'glancing points'. This is the purpose of the following.

Assumption 5 For all η as in Assumption 4, the matrix A(u, η) is diagonalizable.

In particular, for nondissipative isothermal subsonic phase transitions considered, our five assumptions are satisfied; see Section 3 for more details. The existence of surface waves has also been evidenced by Serre [START_REF] Serre | Second-order initial-boundary value problems of variational type[END_REF] in a general framework, when the evolution equations derive from a variational principle.

We enter now into more technical details. Assumption 5 and the fact that A(u, η) has purely imaginary coefficients implies the existence of eigenvalues

β ± i ∈ C and associated eigenvectors R ± i ∈ C 2n , for i ∈ {1, . . . , q ± } with q -:= q = n + p -1, q + := n -p + 1, ( A(u, η) -β ± i I 2n ) R ± i = 0 2n , or equivalently, ( A(u, iη) + β ± i Ȃd (u) ) R ± i = 0 2n , with Re(β ± i ) ≷ 0 , β + i = -β - i , R + i = R - i , i ∈ {1, . . . , m} , Re(β ± i ) = 0 , R ± i ∈ R 2n , i ∈ {m + 1, . . . , q ± } , and
C 2n = Span{R - 1 , . . . , R - q -, R + 1 , . . . , R + q + } = E -(u, η) ⊕ E c (u, η) ⊕ E + (u, η) ,
where E ± (u, η) := Span{R ± 1 , . . . , R ± m } , (we recall that E -(u, η) is the stable subspace of A(u, η), and similarly, E + (u, η) is its unstable subspace), and

E c (u, η) := Span{R - m+1 , . . . , R - q -, R + m+1 , . . . , R + q + }
is the center subspace of A(u, η). Let L ± i ∈ C 2n be such that (L ± i ) * Ȃd (u) are left eigenvectors of the matrix A(u, η), and more precisely,

(L ± i ) * ( A(u, iη) + β ± i Ȃd (u) ) = 0 * 2n .
Above " * " gives the conjugate of the transpose, i.e. A * = t (A). Like the right eigenvectors, they can be chosen so that L + i = L - i , i = 1, . . . , m, and

L ± i ∈ R 2n , i = m + 1, .
. . , q ± . We make the following further assumption.

Assumption 6

(L ± i ) * Ȃd (u) R ± j = 0 , i, j = 1, . . . , q ± , i = j , (L ± i ) * Ȃd (u) R ∓ j = 0 , i = 1, . . . , q ± , j = 1, . . . , q ∓ .
Observe that Assumption 6 is automatic if all the β ± i are distinct. We rescale the eigenvectors so that

(L ± j ) * Ȃd (u) R ± j = 1 , j = 1, . . . , q ± .
Now, Assumption 4 may be interpreted in terms of the eigenvectors R - 1 , . . . , R - q only. We first make some further reductions. Observing that A(u, kη) = k A(u, η) for any k ∈ R (which is due to scale invariance), we see that the subspace E(u, η) is positively homogeneous degree 0 in η. Therefore, the wave vectors η for which Assumption 4 holds true form a positive cone, and for all k > 0

X kη = 1 k X η , V kη = V η .
Thus, without loss of generality, we may assume that η 0 = 1. Then we observe that (2.11) equivalently reads

X = -H 1 (u) V , C(u, η) V = 0 q , (2.12)
where H 1 (u) = dh 1 (u) is the first row of the Jacobian matrix H(u) = dh(u), and C(u, η) ∈ R q×2n is defined by

C(u, η) := T (η) H(u) , T (η) :=      -η 1 . . . I q -η d-1 0 q-d+1      ∈ R q×(q+1) .
Hence, by Assumption 4, there exists

γ ∈ C m \{0 m } such that m j=1 γ j C(u, η)R - j = 0 q ,
the components γ j of γ merely being the components of V η in the basis

{R - 1 , . . . , R - m } of E -(u, η). Moreover, Assumption 4 means that the q × q matrix [C(u, η)R - 1 , . . . , C(u, η)R - q ] is of rank q -1, so that there exists σ ∈ C q \{0 q } such that σ * C(u, η)R - j = 0 , j ∈ {1, . . . , q} . (2.13)
Since the matrix C(u, η) and the vectors R - m+1 , . . . , R - q have real coefficients, and R - j = R + j for j ∈ {1, . . . , m}, we also have, by conjugation,

σ * C(u, η)R + j = 0 , j ∈ {1, . . . , m} , σ * C(u, η)R - k = 0 , k ∈ {m + 1, . . . q} .
(2.14)

Weakly nonlinear surface waves

We can now turn to the derivation of an amplitude equation for weakly nonlinear surface waves in (2.5). Following Hunter's approach [START_REF] Hunter | Nonlinear surface waves[END_REF], we consider an expansion for v, χ, of the form

v ε (y) = u + εv 1 (η 0 y 0 + η • y, y d , εy 0 ) + ε 2 v 2 (η 0 y 0 + η • y, y d , εy 0 ) + O(ε 3 ) , χ ε (y) = εχ 1 (η 0 y 0 + η • y, εy 0 ) + ε 2 χ 2 (η 0 y 0 + η • y, εy 0 ) + O(ε 3 ) ,
where η and y stand for the (d -1)-dimensional vectors defined by η = (η 1 , . . . , η d-1 ) and y = (y 1 , . . . , y d-1 ), and v 1,2 is supposed to go to zero as y d goes to infinity. The above ansatz for v describes a small amplitude wave that is changing slowly in reference frame moving with the wave.

From now on, we use the notation (ξ, z, τ ) = (η 0 • y 0 + η • y, y d , εy 0 ) for the new independent variables. Using Taylor expansions for f i and h,

A i (v ε ) = A i (u) + ε dA i (u) • v 1 + O(ε 2 ) , h(v ε ) = εH(u) • v 1 + ε 2 (H(u) • v 2 + 1 2 d 2 h(u) • (v 1 , v 1 )) + O(ε 3 ) ,
and equating to zero the coefficients of ε and ε 2 in (2.8), we find

L(u, η) • v 1 = 0 2n , z > 0 Jη ∂ ξ χ 1 + H(u) • v 1 = 0 n+p , z = 0 , (2.15)
and

L(u, η) • v 2 + M(u, η; v 1 , ∂ ξ χ 1 ) • v 1 = 0 2n , z > 0 Jη ∂ ξ χ 2 + H(u) • v 2 + G(u; v 1 , ∂ τ χ 1 ) = 0 n+p , z = 0 , (2.16)
where

L(u, η) := A(u, η) ∂ ξ + Ȃd (u) ∂ z , M(u, η; v 1 , ∂ ξ χ 1 ) := A 0 (u)∂ τ + dA(u, η) • v 1 • ∂ ξ + d Ȃd (u) • v 1 • ∂ z -(∂ ξ χ 1 ) Ȃ(u, η) ∂ z , with Ȃ(u, η) := Ȋ2n A(u, η) , and
G(u; v 1 , ∂ τ χ 1 ) := (∂ τ χ 1 ) e 1 + 1 2 d 2 h(u) • (v 1 , v 1 )
, with e 1 denoting the first vector of the canonical basis in C q+1 . We recall that by definition, [START_REF] Benzoni-Gavage | Multi-dimensional hyperbolic partial differential equations: First-order systems and applications[END_REF], or equivalently, (2.12). Denoting by P j and Q j the real and imaginary parts, respectively, of γ j R - j , we have

V η = m j=1 γ j R - j and X η = -H 1 (u) V η solve (2.
V η = P + i Q , with P := m j=1 P j , Q := m j=1 Q j .
For convenience, we also introduce the notations ̺ j and δ j for the real and imaginary parts, respectively, of β - j (j ∈ {1, . . . , m}). In what follows, H stands for the Hilbert transform, such that for any L 2 function w,

FH[w](k) = -i sgn(k) F[w](k) , ∀k ∈ R ,
where F denotes the Fourier transform, with the convention [START_REF] Alì | Nonlinear surface waves on a tangential discontinuity in magnetohydrodynamics[END_REF]) that are square integrable in ξ and such that v 1 goes to zero as z → +∞ are of the form

F[w](k) = w(k) := +∞ -∞ w(x) e -ikx dx , ∀k ∈ R . Proposition 2.1 The solutions (ξ, z) → (v 1 , χ 1 )(ξ, z) of (2.
v 1 (ξ, z) = (w * ξ r)(ξ, z) , r(ξ, z) := - 1 π m j=1 z ̺ j P j + (z δ j + ξ) Q j (̺ j z) 2 + (δ j z + ξ) 2
or equivalently,

v 1 (k, z) = w(k) r(k, z) , r(k, z) =            m j=1 γ j e β - j k z R - j , k > 0 , z > 0 , m j=1 γ j e β + j k z R + j , k < 0 , z > 0 ,
and

v 1 (ξ, 0) = w(ξ) P -H[w](ξ) Q , χ 1 (ξ) = H 1 (u) +∞ ξ v 1 (ζ, 0) dζ ,
where w is an arbitrary L 2 function of zero mean.

Proof. By Fourier transform in the variable ξ, (2.15) becomes

   k A(u, iη) v 1 (k, z) + Ȃd (u) ∂ z v 1 (k, z) = 0 2n , z > 0 , i k χ 1 (k) Jη + H(u) v 1 (k, 0) = 0 n+p .
(2.17)

Similarly as in (2.12), we may eliminate χ 1 from the boundary condition in (2.17). We thus obtain

χ 1 (k) = i k H 1 (u) v 1 (k, 0) , C(u, η) v 1 (k, 0) = 0 q . (2.18)
Since by Assumption 1 the matrix Ȃd (u) is nonsingular, the first line in (2.17) is a genuine ODE on v 1 , which may equivalently be written as

∂ z v 1 = A(u, kη) v 1 , (2.19)
where A(u, kη) is defined as in (2.10) (note that A is homogeneous degree one in η). Then, the vanishing of v 1 at z = +∞ implies that for k > 0, there exists

W (k) ∈ C such that v 1 (k, 0) = W (k) V η = W (k) m j=1 γ j R - j ,
hence, by the solving (2.19),

v 1 (k, z) = exp(z A(u; kη)) v 1 (k, 0) = W (k) m j=1 γ j e β - j k z R - j .
Observing that (k, z) → v 1 (-k, z) solves the same problem as (k, z) → v 1 (k, z), we find that for k < 0, there exists also

W (k) ∈ C such that v 1 (k, z) = W (k) m j=1 γ j e β + j k z R + j , z ≥ 0 .
The conclusion follows by inverse Fourier transform, with w := F -1

a 0 (k) ∂ τ w(k, τ ) + +∞ -∞ a 1 (k -ℓ, ℓ) w(k -ℓ, τ ) w(ℓ, τ ) dℓ = 0 , (2.20)
where a 0 and a 1 are given by (2.24) and (2.25) below.

Proof. By Fourier transform in ξ, (2.16) becomes

   k A(u, iη) v 2 + Ȃd (u) ∂ z v 2 + m 1 = 0 2n , z > 0 , i k χ 2 Jη + H(u) • v 2 + g 1 = 0 n+p , z = 0 , (2.21)
where

m 1 := F[M(u, η; v 1 , ∂ ξ χ 1 ) • v 1 ] , g 1 := F[G(u; v 1 , ∂ τ χ 1 )] .
A crucial fact in what follows on will be that m 1 (k, z, τ ) decays exponentially fast to zero as z goes to +∞, as v 1 itself. We first eliminate χ 2 from the boundary condition in (2.21), as we have made for χ 1 in (2.17). This yields

χ 2 (k) = i k H 1 (u) v 2 (k, 0, τ ) + g 1 1 (k, τ ) ,
where g 1 1 denotes the first component of g 1 , and

C(u, η) v 2 (k, 0, τ ) + T (η) g 1 (k, τ ) = 0 . (2.22) Now, decomposing v 2 as v 2 (k, z, τ ) = q - j=1 ν - j (k, z, τ ) R - j + q + j=1 ν + j (k, z, τ ) R + j ,
thanks to the normalization of left and right eigenvectors, we see that the first equation in (2.21) is equivalent to

∂ z ν ± j -k β ± j ν ± j + (L ± j ) * m 1 = 0 , j ∈ {1, . . . , q ± } .
Solving these ODEs, taking into account the signs of Re(β ± j ) and the fact that m 1 decays exponentially fast to zero as z goes to +∞, we find that for v 2 to decay to zero as z goes to +∞, necessarily

ν - j (k, 0, τ ) = +∞ 0 e -kβ - j z (L - j ) * m 1 (k, z, τ ) dz
for (k < 0 and j ∈ {1, . . . , q -}) or (k > 0 and j ∈ {m + 1, . . . , q -}), and

ν + j (k, 0, τ ) = +∞ 0 e -kβ + j z (L + j ) * m 1 (k, z, τ ) dz
for (k > 0 and j ∈ {1 . . . q + }) or (k < 0 and j ∈ {m + 1 . . . q + }).

Going back to the boundary condition in (2.22) and multiplying it successively by σ * and σ * , we get, thanks to (2.13) and (2.14),

q + j=1 (σ * C(u, η) R + j ) ν + j (k, 0, τ ) + σ * T (η)g 1 (k, τ ) = 0 , m j=1 (σ * C(u, η) R - j ) ν - j (k, 0, τ ) + q + j=m+1 (σ * C(u, η) R + j ) ν + j (k, 0, τ ) + σ * T (η)g 1 (k, τ ) = 0 .
Substituting the integrals found above for ν ± j (k, 0, τ ), we obtain for all k = 0,

+∞ 0 L(k, z) m 1 (k, z, τ ) dz + σ(k) T (η)g 1 (k, τ ) = 0 , (2.23) with σ(k) := σ * for k < 0 and σ(k) := σ * for k > 0, L(k, z) := q + j=1 (σ * C(u, η) R + j ) e -kβ + j z (L + j ) * , k > 0 , L(k, z) := L(-k, z) = m j=1 (σ * C(u, η) R - j ) e -kβ - j z (L - j ) * + q + j=m+1 (σ * C(u, η) R + j )e -kβ + j z (L + j ) * , k < 0 .
Finally, the compatibility equation (2.23) may be rewritten explicitely in terms of the amplitude function w = w(k, τ ) and of the linear surface wave function r = r(k, z) (given by Proposition 2.1). Indeed, substituting w(k, τ ) r(k, z) for v 1 (k, z, τ ) in the definition of m 1 , we get

m 1 (k, z, τ ) = ∂ τ w(k, τ ) A 0 (u) r(k, z) + +∞ -∞ m(u, η; k-ℓ, ℓ, z) w(k-ℓ, τ ) w(ℓ, τ ) dℓ , 2π • m(u, η; k, ℓ, z) = i ℓ dA(u, η) • r(k, z) • r(ℓ, z) + i ℓ dA d (u) • r(k, z) • F[r ′ ](ℓ, z) + i ℓ H 1 (u) r(k, 0) A(u, η)F[r ′ ](ℓ, z) ,
where we have introduced a new vector-valued function r ′ , defined by

i ℓ F[r ′ ](ℓ, z) := Ȋ2n ∂ z r(ℓ, z) ,
or equivalently,

F[r ′ ](k, z) =            -i Ȋ2n m j=1 γ j β - j e β - j k z R - j , k > 0 , z > 0 , -i Ȋ2n m j=1 γ j β + j e β + j k z R + j , k < 0 , z > 0 .
To find the last term in the kernel m, we have used the expression of χ 1 given by (2.18). This expression is also useful to compute

g 1 (k, τ ) = i k ∂ τ w(k, τ )H 1 (u) r(k, 0) e 1 + +∞ -∞ g(u, η; k -ℓ, ℓ) w(k -ℓ, τ ) w(ℓ, τ ) dℓ , g(u, η; k, ℓ) := 1 4π d 2 h(u) • ( r(k, 0), r(ℓ, 0)) .
We have thus obtained the nonlocal equation (2.20) for w, with

a 0 (k) := +∞ 0 L(k, z) A 0 (u) r(k, z) dz + i k (H 1 (u) r(k, 0)) σ(k) T (η) e 1 , (2.24
)

a 1 (k, ℓ) := +∞ 0 L(k + ℓ, z) m(u, η; k, ℓ, z) dz + σ(k + ℓ) T (η) g(u, η; k, ℓ) .
(2.25) ✷ Remark 2.3 Since L(k, z) and r(k, z) are linear combinations of exponentials e -kβ + j z and e kβ - p z , and, by construction, L(k, z) = L(-k, z), r(k, z) = r(-k, z), σ(k) = σ(-k), we see on (2.24) that a 0 is of the form

a 0 (k) =    α 0 k , k > 0 , -α 0 k , k < 0 , (2.26)
with α 0 ∈ C. More explicitly, this number is given by

α 0 = σ * C(u, η) R + j β + j -β - p (L + j ) * A 0 (u)(γ p R - p ) + i (γ p H 1 (u)R - p ) σ * T (η) e 1 , (2.27)
where we have used the usual summation convention on the repeated indices, with j ∈ {1, . . . , q + }, p ∈ {1, . . . , m}.

Remark 2.4

The kernel a 1 is obviously not symmetric in (k, ℓ). However, it can easily be symmetrized. Indeed, by change of variables ℓ → kℓ, the nonlocal equation (2.20) is equivalent to

a 0 (k) ∂ τ w(k, τ ) + +∞ -∞ a s 1 (k -ℓ, ℓ) w(k -ℓ, τ ) w(ℓ, τ ) dℓ = 0 , with a s 1 (k, ℓ) := +∞ 0 L(k + ℓ, z) m s (u, η; k, ℓ, z) dz + σ(k + ℓ) T (η) g(u, η; k, ℓ) ,
(2.28) 4π m s (u, η; k, ℓ, z)

:= i(k + ℓ) dA(u, η) • r(k, z) • r(ℓ, z) + iℓ dA d (u) • r(k, z) • F[r ′ ](ℓ, z) + ik dA d (u) • r(ℓ, z) • F[r ′ ](k, z) + iℓ H 1 (u) r(k, 0) A(u, η)F[r ′ ](ℓ, z) + ik H 1 (u) r(ℓ, 0) A(u, η)F[r ′ ](k, z) .
(2.29) (The first term is indeed symmetric by the symmetry of dA, as a linear combination of second order differentials d 2 f j . For the same reason, g being defined by means of d 2 h, it is symmetric in (k, ℓ).) In addition, both the integral and the last term in a s 1 (k, ℓ) are positively homogeneous degree zero in (k, ℓ). Theorem 2.5 Under the Assumptions 1, 2, 3, 4, 5, 6, we also assume that the number α 0 defined in (2.27) is nonzero. Then weakly nonlinear surface waves for the nonlinear model (2.5) are governed by the nonlocal amplitude equation

∂ τ w + ∂ ξ Q [w] = 0 , (2.30)
where Q is given by

Q [v] (ξ) = (K * (v ⊗ v))(ξ, ξ) (2.31)
for all v ∈ S (R) (the Schwartz class), the kernel K being the real tempered distribution on

R 2 K := 2πF -1 (Λ) , (2.32)
with Λ ∈ L ∞ (R 2 ) defined as in (2.33).

Proof. With the notations introduced above, we define for k = 0, ℓ = 0,

k + ℓ = 0, Λ(k, ℓ) := a s 1 (k, ℓ) i (k + ℓ) a 0 (k + ℓ) . (2.33)
Then (2.20) can be rewritten as

∂ τ w(k, τ ) + ik +∞ -∞ Λ(k -ℓ, ℓ) w(k -ℓ, τ ) w(ℓ, τ )dℓ = 0 . (2.34)
By inverse Fourier transform this gives (2.30) with, formally,

Q [w] (ξ, τ ) := 1 2π +∞ -∞ +∞ -∞ e i(k+ℓ)ξ Λ(k, ℓ) w(k, τ ) w(ℓ, τ ) dℓ dk , or, Q [w] (ξ, τ ) = 2πF -1 (Λ w(•, τ ) ⊗ w(•, τ ))(ξ, ξ) ,
where F here denotes the Fourier transform on S ′ (R 2 ). Since F -1 ( w⊗ w) = w ⊗ w, we find that

2πF -1 (Λ w(•, τ ) ⊗ w(•, τ )) = K * (w(•, τ ) ⊗ w(•, τ )) ,
with K := 2πF -1 (Λ). To justify the above computations, we first observe that the kernel Λ has some nice properties inherited from the properties of a s 1 and a 0 . It is indeed smooth (analytic) outside the lines k = 0, ℓ = 0, and k + ℓ = 0, symmetric in (k, l), like a s 1 , and positively homogeneous degree zero, like a s 1 and k → ka 0 (k). In addition, since a 0 (-k) = a 0 (k) and

a s 1 (-k, -ℓ) = a s 1 (k, ℓ), we have Λ(-k, -ℓ) = Λ(k, ℓ).
To summarize, we have for all k = 0, ℓ = 0, and

θ > 0, Λ(k, ℓ) = Λ(ℓ, k) , Λ(-k, -ℓ) = Λ(k, ℓ) , Λ(θk, θℓ) = Λ(k, ℓ) .
(2.35)

Using these properties and noting that Λ(1, θ) and Λ(-1, θ) are uniformly bounded for θ ∈ (0, 1), we easily check that Λ is bounded on (R\{0}) 2 . Thus it can be viewed as a tempered distribution, and K is therefore well-defined by (2.32) as a tempered distribution. Furthermore, the second property in (2.35) shows that K is a real distribution.

To conclude, for all v ∈ L 2 (R), (2.36)

k → +∞ -∞ Λ(k -ℓ, ℓ) v(k -ℓ) v(ℓ)dℓ
He had in particular checked it was satisfied in the case of weakly nonlinear surfaces waves in Elasticity [START_REF] Hunter | Nonlinear surface waves[END_REF][START_REF] Parker | Waveform evolution for nonlinear surface acoustic waves[END_REF][START_REF] Parker | Analysis and computation for nonlinear elastic surface waves of permanent form[END_REF]. More recently, he and co-workers derived and investigated a stronger condition [START_REF] Alì | Hamiltonian equations for scale-invariant waves[END_REF][START_REF] Alì | Nonlinear surface waves on a tangential discontinuity in magnetohydrodynamics[END_REF], which ensures that (2.34) has a Hamiltonian structure (see [START_REF] Hunter | Short-time existence for scale-invariant Hamiltonian waves[END_REF] for a local-in-time existence under this condition in a periodic setting). It turns out that (2.36) is in fact exactly what we need to get a priori estimates without loss of derivatives for (2.34), see [START_REF] Benzoni-Gavage | On nonlocal Burgers equations[END_REF]. This is the condition we are going to investigate further in our abstract framework and afterwards in the explicit case of subsonic phase boundaries.

Proposition 2.6 For Λ defined as in (2.33) with a 0 given by (2.26) (2.27), α 0 being assumed to be nonzero, and a s 1 given by (2.28) (2.29), the stability condition (2.36) is equivalent to requiring that a(P ) and a(Q) be real, with a the linear form a : C 2n → C defined by

α 0 a(R) = -iσ * T d 2 h • (R, V ) + q + j=1 m p=1 σ * C R + j β + j -β - p (L + j ) * (dA -iβ - p d Ȃd ) • (γ p R - p ) • R -iβ - p (H 1 R) Ȃ(γ p R - p ) ,
where, for simplicity, underlined letters correspond to quantities evaluated at u and/or η, while, as in Proposition 2.1,

V = m p=1 γ p R - p , P = Re m p=1 γ p R - p , Q = Im m p=1 γ p R - p .
Proof. By direct computation we find that

4π α 0 Λ(1, 0+) = -i σ * T d 2 h • (V, V ) + σ * C R + j β + j -β - p (L + j ) * (dA -iβ - p d Ȃd ) • (γ p R - p ) • V -iβ - p (H 1 V ) Ȃ(γ p R - p ) , 4π α 0 Λ(1, 0-) = -i σ * T d 2 h • (V , V ) + σ * C R + j β + j -β - p (L + j ) * (dA -iβ - p d Ȃd ) • (γ p R - p ) • V -iβ - p (H 1 V ) Ȃ(γ p R - p ) ,
where for simplicity we have used the convention of summation over repeated indices. Observe that (2.36) is equivalent to require that Λ(1, 0+) + Λ(1, 0-) ∈ R and Λ(1, 0+) -Λ(1, 0-) ∈ iR, or that the sum of the above equalities divided by 2α 0 and their difference divided by 2iα 0 must be real. ✷

Application to van der Waals fluids

In this section we apply the method of the previous section to a concrete model for fluids exhibiting phase changes, and obtain an explicit form for the kernel as in Theorem 2.5.

Introduction

The Euler equations governing the motion in R d , d ≥ 1, of a compressible, non-viscous, isothermal fluid of van der Waals are

∂ t ρ + ∇ • (ρv) = 0 , ∂ t (ρv) + ∇ • (ρv ⊗ v) + ∇p = 0 d . (3.37)
Above ρ > 0 denotes the density, v ∈ R d the velocity and p > 0 the pressure of the fluid obeying the pressure law

p(V ) = R T V -b - a V 2 ,
where V := 1/ρ is the specific volume, T is the temperature, R is the perfect gas constant and a, b are positive constants. Below the critical temperature, T c := 8a/(27bR), van der Waals fluids can undergo transitions between two phases, the liquid phase for 1/ρ ∈ (0, V * ) and the vapor phase for 1/ρ ∈ (V * , ∞), for the presence of the nonphysical region (V * , V * ), called the spinodal region. The van der Waals law is considered here for concreteness, but our results do not depend on the actual form of this law. They basically depend on the existence of three zones, namely the intervals (0, V * ) and (V * , ∞) where the pressure is decreasing with 1/ρ and the system (3.37) is hyperbolic, and the interval (V * , V * ) where the pressure is increasing with 1/ρ and the system (3.37) becomes elliptic. In this situation it is natural to consider (weak) solutions to (3.37) that avoid the spinodal region.

The simplest weak solutions in this case are piecewise C 1 functions which satisfy (3.37) outside a moving interface Σ(t), and, at least, the Rankine-Hugoniot jump conditions across the interface. We are interested here in dynamic discontinuities, for which there is some mass tranfer across the interface, and especially the subsonic ones, for which the Mach numbers with respect to the interface are lower than one on both sides. In the terminology of hyperbolic conservation laws these discontinuities are undercompressive, the number of outgoing characteristics being equal to that of incoming ones, and an additional jump condition is thus needed. In the continuation of [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF], we have chosen a simple and explicit additional condition, referred to as the capillarity criterion merely because it is equivalent to the existence of travelling capillarity profiles. It can be understood as the conservation of 'total energy', namely the kinetic energy plus the free energy, across the interface. It amounts to neglecting dissipation due to viscosity, which is reasonable in some physicals contexts (e.g. for water in extreme conditions or for superfluids). The well-posedness of the full nonlinear problem in this situation has been proved by Coulombel and Secchi [START_REF] Coulombel | Nonlinear compressible vortex sheets in two space dimensions[END_REF].

The nonlinear problem and the reference phase boundary

We consider a problem of the form (2.6) with n = d + 1, p = 1, and

u := ρ  , f 0 (u) := u, f i (u) :=  i p(ρ) e i +  i ρ  , f 0 (u) := u  2 2ρ + ρ f (ρ) , f i (u) := f i (u)  2 2ρ 2 + f (ρ) + p(ρ) ρ  i ,
where (e 1 , . . . , e d ) is the canonical basis of R d ,

 = ( 1 , . . . ,  d ) := ρ v ∈ R d
is the momentum, and f = f (ρ) is the free specific energy of the fluid, characterized by

p(ρ) = ρ 2 f ′ (ρ). (3.38)
By definition, a solution of (2.6) with χ ≡ 0 and u constant on either side of the hyperplane {x ∈ R d :

x d = 0} is characterized by [ f d (u)] = 0 d+2 , that is, [ d ] = 0 , p(ρ) e d +  d ρ  = 0 d ,  2 2ρ 2 + f (ρ) + p(ρ) ρ  d = 0 . (3.39)
As is very well-known, the first two equations imply that for a dynamical discontinuity, for which  d = 0, the jump of the tangential velocity must be zero, that is,

[v 1 ] = • • • = [v d-1 ] = 0
. By a change of Galilean frame we may assume without loss of generality that the tangential velocity of the left and right reference states is zero. With this simplification, the jump conditions in (3.39) reduce to

[ d ] = 0 , p(ρ) +  2 d ρ = 0 ,  2 d 2ρ 2 + f (ρ) + p(ρ) ρ = 0 . (3.40)
For later use, it is convenient to introduce the functions

q : (ρ, j) → p(ρ) + j 2 ρ , z : (ρ, j) → j j 2 2ρ 2 + f (ρ) + p(ρ) ρ .
Notice that using these functions the jump conditions in (3.40) equivalently read

[ d ] = 0 , [q(ρ,  d )] = 0 , [z(ρ,  d )] = 0 . (3.41)
It is not difficult to show that for a non-monotone pressure law ρ → p(ρ),

there exist ρ l , ρ r , v l , v r satisfying (3.41) with ρ l v l = ρ r v r =:  d > 0, that is, q(ρ l ,  d ) = q(ρ r ,  d ) and z(ρ l ,  d ) = z(ρ r ,  d ), together with q ρ (ρ l,r ,  d ) = 0, that is p ′ (ρ l,r ) -v 2 l,r = 0; see [2, page 249].
The corresponding reference states u l = t (ρ l , 0, . . . , 0, v l ) and u r = t (ρ r , 0, . . . , 0, v r ) are thus connected by a planar dynamical subsonic phase boundary located at x d = 0, and u = (u l , u r ), satisfies Assumption 1. From now on, we fix u as above, and we introduce the notations c l,r for the sound speeds on each side of the reference phase boundary: c l,r := p ′ (ρ l,r ) .

Linearization

Proceeding as in Section 2.2, we may reformulate the free boundary problem (2.5) with our specific fluxes in terms of ρ ± (y 0 , y 1 , . . . , y d ) := ρ(y 0 , . . . , y d-1 , χ(y 0 , . . . , y d-1 ) ± y d ) ,  ± (y 0 , y 1 , . . . , y d ) := (y 0 , . . . , y d-1 , χ(y 0 , . . . , y d-1 ) ± y d ) ,

as            L(ρ ± ,  ± , ∇χ) •     ρ -  - ρ +  +     = 0 2d+2 , y d > 0 b(ρ ± ,  ± , ∇χ) = 0 d+2 , y d = 0 . (3.42)
Linearizing this problem about (ρ -≡ ρ l ,  -≡ (0, . . . , 0, ρ l v l ), ρ + ≡ ρ r ,  + ≡ (0, . . . , 0, ρ r v r ), χ ≡ 0), we readily get a system of the form (2.9), without having to invoke Assumption 2 for the reduction of the boundary conditions. Indeed, for the specifix fluxes we are considering, the linearized version of the jump conditions in (2.6) turns out to reduce to

                     [ρ] ∂ t χ = [ d ] , [p] ∂ i χ = [v i ] , i ∈ {1, . . . , d -1} , (c 2 -v 2 ) ρ + 2 v d = 0 , 1 2 ρ v 2 + ρ f ∂ t χ = (c 2 -v 2 )v ρ + (f + p ρ + 3 2 v 2 ) d , (3.43)
which is obviously of the form

J(u)∇ χ + H(u) • ( ρ, ) = 0 d+2 , (3.44)
with J(u) a matrix depending only on the reference state, as well as H(u), and

( ρ, ) :=     ρ- - ρ+ +     .
Regarding the linearized version of the interior equation in (2.6), it is given by the block-diagonal operator

L(ρ l,r , 0 2(d-1) , v l,r , 0 d ) = L -(ρ l , v l ) 0 0 L + (ρ r , v r ) , ,
the operators L ± being defined in tangential Fourier variables by

L ± (ρ, v, η) =       iη 0 iη ±∂ y d ip ′ (ρ) t η (iη 0 ± v∂ y d )I d-1 0 d-1 ±(p ′ (ρ) -v 2 )∂ y d iv η iη 0 ± 2v∂ y d      
, where η := (η 1 , . . . , η d-1 ). The subsonicity of the reference phase boundary (c l,r > v l,r ) and the additional assumption

η 2 0 < (c 2 l,r -v 2 l,r ) η 2 (3.45)
imply that Assumption 5 is satisfied (see [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF]). In this case, with the notations of Section 2.2 we have,

n = d + 1, p = 1, q -= q + = d + 1, m = 2 , the eigenvalues β ± 1 being the roots of (c 2 l -v 2 l )β 2 + 2 iη 0 v l β + η 2 0 -c 2 l η 2 = 0 ,
and the eigenvalues β ± 2 being the roots of

(c 2 r -v 2 r )β 2 -2 iη 0 v r β + η 2 0 -c 2 r η 2 = 0 ,
which gives explicitly,

β + 1 = 1 c 2 l -v 2 l (α l -iη 0 v l ) , β - 1 = 1 c 2 l -v 2 l (-α l -iη 0 v l ) , α l := c l (c 2 l -v 2 l ) η 2 -η 2 0 , β + 2 = 1 c 2 r -v 2 r (α r + iη 0 v r ) , β - 2 = 1 c 2 r -v 2 r (-α r + iη 0 v r ) , α r := c r (c 2 r -v 2 r ) η 2 -η 2 0 . (3.46)
The other, purely imaginary eigenvalues are

β + 3 := i η 0 v l , β - 3 := -i η 0 v r ,
of multiplicity d -1. Right eigenvectors may be chosen as follows

R ± 1 = r ± 1 0 d+1 , R ± 2 = 0 d+1 r ± 2 , r - 1 =   -iη 0 + v l β - 1 ic 2 l η α l   , r + 1 =   iη 0 -v l β + 1 -ic 2 l η α l   , r - 2 =   -iη 0 -v r β - 2 ic 2 r η -α r   , r + 2 =   iη 0 + v r β + 2 -ic 2 r η -α r   , R + i = r + i 0 d+1 , R - i = 0 d+1 r - i , r + i =   0 η 0 ěi-2 v l η • ěi-2   , r - i =   0 η 0 ěi-2 v r η • ěi-2   , i = 3, . . . , d + 1 , (3.47)
where (ě 1 , . . . , ěd-1 ) is an arbitrary basis of R d-1 . For left eigenvectors, we may take

L ± 1 = l ± 1 0 d+1 , L ± 2 = 0 d+1 l ± 2 , l - 1 =    iη 0 + 2v l β - 1 -iη -β - 1    , l + 1 =    -iη 0 -2v l β + 1 iη β + 1    , l - 2 =    -iη 0 + 2v r β - 2 iη -β - 2    , l + 2 =    iη 0 -2v r β + 2 -iη β + 2    , L + i = l + i 0 d+1 , L - i = 0 d+1 l - i , l + i =   -v 2 l η • ě′ i-2 η 0 ě′ i-2 v l η • ě′ i-2   , l - i =   -v 2 r η • ě′ i-2 η 0 ě′ i-2 v r η • ě′ i-2   , i = 3, . . . , d + 1 , (3.48)
where

(ě ′ 1 , . . . , ě′ d-1 ) is another arbitrary basis of R d-1 . Recalling that Ȃd (u) = -A d (ρ l , v l ) 0 (d+1)×(d+1) 0 (d+1)×(d+1) A d (ρ r , v r ) , A d (ρ, v) =       0 0 * d-1 1 0 d-1 vI d-1 0 d-1 p ′ (ρ) -v 2 0 * d-1 2v       , we easily compute ∓ (L ± i+2 ) * Ȃd (u)R ± k+2 = v 3 l,r (η • ě′ i ) (η • ěk ) + v l,r η 2 0 (ě ′ i • ěk )
for i, k ∈ {1, . . . , d -1}. So, even though the eigenvalues β ± 3 are non-simple, it is possible to choose the bases (ě 1 , . . . , ěd-1 ) and (ě ′ 1 , . . . , ě′ d-1 ) to have

(L ± i ) * Ȃd (u)R ± k = 0 , i, k ∈ {3, . . . , d + 1} , i = k .
For instance, we can take ě1 = ě′ 1 = η , and choose 'dual' bases (ě 2 , . . . , ěd-1 ) and (ě ′ 2 , . . . , ě′ d-1 ) of η⊥ . The left and right eigenvectors above are not normalized to have (

L ± i ) * Ȃd (u)R ± i = 1.
from (3.43). Indeed, subtracting (g + 1 2 v 2 ) times the first equation to the (d + 2)th equation in (3.43), we may replace the latter by

-[p] ∂ t χ = (c 2 -v 2 )v ρ + v 2 d .
Then, substituting X e i(η 0 t + η•y) for χ in (3.43), we can complete the elimination of χ. Since η 0 = 0 (and therefore also η = 0 by (3.45)), we have

X = [ d ] iη 0 [ρ] e -i(η 0 t + η•y) ,
and we are left with the following algebraic system for ( ρ, ),

           η 0 [v ] • η -η 2 v l v r [ d ] = 0 , (c 2 -v 2 ) ρ + 2 v d = 0 , (c 2 -v 2 )v ρ + (v 2 + v l v r ) d = 0 ,
with the additional condition that [v ] be colinear to η. The rest of the proof is a matter of elementary algebra and is left to the reader. ✷ Remark 3.2 As was observed in [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF], the linear surface waves found in Lemma 3.1 are slow, in that for (η 0 , η) ∈ R d solution of (3.50) with (3.45) and (3.46), we have the inequality

η 2 0 < v l v r η 2 . (3.52)
This inequality will be used later on.

Notice that, in terms of the abstract form (3.44) of (3.43) for χ = χ(η 0 t + η • y), the method of elimination used in Lemma 3.1 above leads to the system

(∂ ξ χ) J(u) η + H(u, η) • ( ρ, ) = 0 d+2 , with J(u) η =     -η 0 [ρ] 0 d-1 0 0     , H(u, η) • ( ρ, ) =       [ d ] Π η [v ] η 0 [v ] • η -η 2 v l v r [ d ] (c 2 -v 2 ) ρ + 2 v d (c 2 -v 2 )v ρ + (v 2 + v l v r ) d      
, where the brackets stand as usual for "jumps" (e.g. [ d ] = ( d ) + -( d ) -), and Π η denotes the (d-2)×(d-1) matrix whose rows are t ěi for i = 2, . . . , d-1. (Recall that (ě 2 , . . . , ěd-1 ) has been chosen to be a basis of η⊥ .) We have H(u, η) = E(u, η)H(u, η) and J(u) = E(u, η)J(u), with

E(u, η) =       1 0 * d-1 0 0 0 d-2 Π η 0 d-2 0 d-2 -η 2 v l v r η 0 t η 0 0 0 0 * d-1 1 0 v l v r -(g + 1 2 v 2 ) 0 * d-1 0 1       .
Denoting by C(u, η) the linear mapping

( ρ, ) →     Π η [v ] η 0 [v ] • η -η 2 v l v r [ d ] (c 2 -v 2 ) ρ + 2 v d (c 2 -v 2 )v ρ + (v 2 + v l v r ) d     ,
(that is, we retain all but the first row in H(u, η)•( ρ, )), Lemma 3.1 says that the matrix made up with the column vectors (C(u, η)R - 1 , . . . , C(u, η)R - d+1 ) is of rank d and

γ 1 C(u, η)R - 1 + γ 2 C(u, η)R - 2 = 0 d+1 .
By definition (see (2.13)), the vector σ = (σ -d+2 , . . . , σ -1 , σ 1 , σ 2 , σ 3 ) must be orthogonal (in C d+1 equipped with the standard hermitian product) to

C(u, η)R - 1 =     0 d-2 η 2 v l (v r α l -iη 0 c 2 l ) -v l α l + iη 0 c 2 l v l (-v r α l + iη 0 c 2 l )     , C(u, η)R - 2 =     0 d-2 η 2 v r (v l α r + iη 0 c 2 r ) -v r α r -iη 0 c 2 r -v r (v l α r + iη 0 c 2 r )     , C(u, η)R - i =     η 0 v r Π η ěi-2 (η 2 0 -v l v r η 2 ) v r (η • ěi-2 ) 2v 2 r (η • ěi-2 ) v 2 r (v l + v r )(η • ěi-2 )     , i = 3, . . . , d + 1,
or, with the choice of the vectors ěj made above,

C(u, η)R - 3 =     0 d-2 (η 2 0 -v l v r η 2 ) v r η 2 2v 2 r η 2 v 2 r (v l + v r ) η 2     , C(u, η)R - i =     η 0 v r Π η ěi-2 0 0 0    
for i = 4, . . . , d + 1. Since the (d -2) × (d -2) matrix (Π η ě2 , . . . , Π η ěd-1 ) is nonsingular, we easily see that σ must be of the form σ = (0, . . . , 0, σ 1 , σ 2 , σ 3 ). Furthermore, taking for instance

σ 2 = -v l α r + iη 0 c 2 r ,
we find that

σ 3 -σ 1 η 2 = α r -iη 0 c 2 r /v r = - γ 1 v r ,
and

σ 1 = -(v r -v l ) v r α r + iη 0 c 2 r η 2 0 + η 2 v 2 r = γ 1 v r -v l η 2 0 + η 2 v 2 r .
Knowing that σ * C(u, η)R - 1,2 = 0, that C(u, η) has real coefficients and that R + 1,2 = R - 1,2 , we readily get

σ * C(u, η)R + 1 = 2 σ * Re C(u, η)R - 1 = 2 α l v l (v r σ 1 η 2 -σ 2 -v r σ 3 ) = -2 α l α r v l (v r -v l ) , σ * C(u, η)R + 2 = 2 σ * Re C(u, η)R - 2 = 2 α r v r (v l σ 1 η 2 -σ 2 -v l σ 3 ) = 2 i η 0 c 2 r α r (v r -v l ) . Finally, since C(u, η)R + 3 =     0 d-2 -(η 2 0 -v l v r η 2 ) v l η 2 -2v 2 l η 2 -v 2 l (v l + v r ) η 2     resembles C(u, η)R - 3 , it is not difficult to evaluate σ * C(u, η)R + 3 = γ 1 [v] 2 η 2 0 + η 2 v 2 r η 2 v l v r (v l v r η 2 -η 2 0 ) .
(Observe that by (3.52), the last factor here above is positive.) To summarize we have

                       σ * C(u, η)R + 1 = -2 α l α r v l [v] , σ * C(u, η)R + 2 = 2 i η 0 c 2 r α r [v] , σ * C(u, η)R + 3 = γ 1 η 2 v l v r v l v r η 2 -η 2 0 η 2 0 + η 2 v 2 r [v] 2 ,
σ * C(u, η)R + j = 0 , j ∈ {4, . . . , d + 1} .

(3.53)

Weakly nonlinear surface waves for phase boundaries

We are now going to derive the explicit form of the first and second order systems associated with our specific free boundary problem (3.42). Since the boundary condition in (3.42) is not as decoupled as in (2.8), the resulting second order system will look slightly more complicated than (2.16).

To avoid multiple indices we prefer using here notations with dots instead of subscripts 1 and 2 for the first order and second orders of the expansion. Then the first and second order systems associated with (3.42)

are of the form    L(u, η) • ( ρ, ) = 0 2d+2 , z > 0 , (∂ ξ χ) J(u)η + H(u) • ( ρ, ) = 0 d+2 , z = 0 , (3.54)    L(u, η) • (ρ,) + M(u, η; ρ, , ∂ ξ χ) • ( ρ, ) = 0 2d+2 , z > 0 , (∂ ξ χ) J(u)η + H(u) • (ρ,) + G(u, η; ρ, , ∂ τ χ, ∂ ξ χ) = 0 d+2 , z = 0 , (3.55) with ( ρ, ) :=     ρ- - ρ+ +     , (ρ,) :=     ρ- - ρ+ +     .
The linear terms in (3.54) and (3.55) are of course reminiscent of the linearized problem considered in Section 3.3. The operator L(u, η) is blockdiagonal and defined by

L(u, η) = L -(ρ l , v l , η) 0 (d+1)×(d+1) 0 (d+1)×(d+1) L + (ρ r , v r , η) , L ± (ρ, v, η) =       η 0 ∂ ξ η∂ ξ ±∂ z p ′ (ρ) t η∂ ξ (η 0 ∂ ξ ± v∂ z )I d-1 0 d-1 ±(p ′ (ρ) -v 2 )∂ z v η∂ ξ η 0 ∂ ξ ± 2v∂ z       .
In the boundary condition we have

J(u)η :=     -η 0 [ρ] -[p]η 0 -η 0 [ρ(g + 1 2 v 2 ) -p]     , H(u) • (ρ,) =     [ d ] [ v  ] (c 2 -v 2 ) ρ + 2 v d (c 2 -v 2 )v ρ + (g + 3 2 v 2 ) d     .
We recall that c denotes the sound speed (c(ρ) 2 = p ′ (ρ)) and g denotes the chemical potential (g(ρ) = f (ρ) + p(ρ) ρ ) of the fluid. The other terms in (3.55) are of the form

M(u, η; ρ, , ∂ ξ χ) • ( ρ, ) = ∂ τ ( ρ, ) + ∂ ξ (Q(u, η)( ρ, )) + ∂ z (P(u)( ρ, )) -(∂ ξ χ) ∂ z (N (u, η)( ρ, )) ,
where Q(u, η) and P(u) are quadratic mappings and N (u, η) is a linear mapping, and

G(u, η; ρ, , ∂ τ χ, ∂ ξ χ) = (∂ τ χ) e -(∂ ξ χ) N (u, η)( ρ, ) + P (u)( ρ, ) ,
where P (u) is another quadratic map (involving P(u)), N (u, η) is another linear map (involving N (u, η)), and

e(u) :=   -[ρ] 0 d -[ρ(g + 1 2 v 2 ) -p]   .
Explicit formulas are

N (u, η)( ρ, ) = -N 0 (ρ l , v l , η)( ρ-, -) N 0 (ρ r , v r , η)( ρ+ , + ) , N 0 (ρ, v, η)( ρ, ) =   η 0 ρ + η •  η 0  + p ′ (ρ) ρη η 0 d + v η •    , Q(u, η)( ρ, ) = Q 0 (ρ l , v l , η)( ρ-, -) Q 0 (ρ r , v r , η)( ρ+ , + ) , Q 0 (ρ, v, η)( ρ, ) =   0 1 ρ (η • )  + 1 2 p ′′ (ρ)( ρ) 2 η 1 ρ (η • )( d -v ρ)   , P(u)( ρ, ) = -P 0 (ρ l , v l )( ρ-, -) P 0 (ρ r , v r )( ρ+ , + ) , P 0 (ρ, v)( ρ, ) =    0 1 ρ ( d -v ρ)  ( 1 2 p ′′ (ρ) + v 2 ρ )( ρ) 2 + 1 ρ d ( d -2v ρ)    , P (u)( ρ, ) = 2 P 0 (ρ r , v r )( ρ+ , + ) -P 0 (ρ l , v l )( ρ-, -) π(ρ r , v r )( ρ+ , + ) -π(ρ l , v l )( ρ-, -) , π(ρ, v)( ρ, ) = v ρ ( 1 2  2 + ( d ) 2 ) + 1 2 v(p ′′ (ρ) -p ′ (ρ)-3v 2 ρ )( ρ) 2 + p ′ (ρ)-3v 2 ρ ρ d , N (u, η)( ρ, ) = 2 N 0 (ρ r , v r , η)( ρ+ , + ) -N 0 (ρ l , v l , η)( ρ-, -) ν(ρ r , v r , η)( ρ+ , + ) -ν(ρ l , v l , η)( ρ-, -) , ν(ρ, v, η)( ρ, ) = η 0 (g(ρ) -1 2 v 2 ) ρ + (g(ρ) + 1 2 v 2 ) (η • ) + η 0 v d . From Lemma 3.
1 the resolution of (3.54) is given by the following analogue of Proposition 2.1.

Proposition 3.3

The solutions (ξ, z) → ( ρ, , χ)(ξ, z) of (3.54) that are square integrable in ξ and such that ρ,  go to zero as z → +∞ are of the form

( ρ, )(ξ, z) = (w * ξ r)(ξ, z) , χ(ξ) = (w * ξ s)(ξ) , r(k, z) = γ 1 e β - 1 k z R - 1 + γ 2 e β - 2 k z R - 2 , k > 0 , z > 0 , γ 1 e β + 1 k z R + 1 + γ 2 e β + 2 k z R + 2 , k < 0 , z > 0 , s(k) =          - γ 2 α r + γ 1 α l ikη 0 [ρ] , k > 0 , - γ 2 α r + γ 1 α l ikη 0 [ρ] , k < 0 ,
where w is an arbitrary L 2 function.

Now, since the triangular matrix E(u, η) is nonsingular (for η = 0), the boundary condition in the second order system (3.55) is equivalent to

(∂ ξ χ) J(u)η + H(u, η) • (ρ,) + G(u, η; ρ, , ∂ τ χ, ∂ ξ χ) = 0 d+2 , with G(u, η; •) = E(u, η)G(u, η; •) (and as before J(u) = E(u, η)J(u), H(u, η) = E(u, η)H(u, η)), or, isolating the first row,    -η 0 [ρ](∂ ξ χ) + [ d ] -η 0 [ρ](∂ τ χ) -2 [η 0 ρ + η • ] (∂ ξ χ) = 0 , C(u, η) • (ρ,) + g(u, η; ρ, , ∂ τ χ, ∂ ξ χ) = 0 d+1 , (3.56)
where

g(u, η; ρ, , ∂ τ χ, ∂ ξ χ) := e(u, η) G(u, η; ρ, , ∂ τ χ, ∂ ξ χ) . e(u, η) :=     0 d-2 Π η 0 d-2 0 d-2 -η 2 v l v r η 0 t η 0 0 0 0 * d-1 1 0 v l v r -(g + 1 2 v 2 ) 0 * d-1 0 1     .
It turns out that the factor of ∂ τ χ in g reduces to

e(u, η) e(u) =     0 d-2 η 2 v l v r [ρ] 0 0     . Theorem
:= 2πF -1 (Λ) ∈ S ′ (R 2 ), in which Λ is defined as in (2.33) by Λ(k, ℓ) := a s 1 (k, ℓ) i (k + ℓ) a 0 (k + ℓ) , a 0 (k) =    α 0 k , k > 0 , α 0
k , k < 0 , with α 0 and a s 1 defined in (3.59) (3.60) (3.61) (3.62) below. This kernel Λ is well defined because α 0 = 0. In addition, it satisfies the reality-symmetryhomogeneity properties in (2.35), and the stability condition (2.36) is equivalent to requiring that a(Re(γ

1 R - 1 + γ 2 R - 2 )) and a(Im(γ 1 R - 1 + γ 2 R - 2
)) be real, with the linear form defined by (3.63) below.

Proof. Similarly as in the abstract framework of Section 2.3, using the reformulation (3.56) of the boundary condition in (3.55), we find that for (3.55) to have a L 2 solution, the first order solution ( ρ, , χ) of (3.54) must satisfy .

Concerning the last term in (3.62) we find the value

i η 0 [ρ] (γ 2 α r + γ 1 α l ) σ 1 v l v r [ρ] η 2 = i [v]v l v r η 2 η 0 (η 2 0 + η 2 v 2 r ) (α l |γ 1 | 2 + α r γ 2 γ 1 ) .
And finally, after computing that

(L + 3 ) * R - 1 = c 2 l η 2 c 2 l -v 2 l γ 2 and β + 3 -β - 1 = γ 2 v l (c 2 l -v 2 l )
, we find that the third term in (3.62) equals to

-[v] 2 |γ 1 | 2 c 2 l η 2 v l (v l v r η 2 -η 2 0 ) v r µ l µ r ,
with µ l,r := η 2 0 + η 2 v 2 l,r . (It can easily be checked that each of these terms has the dimension of c 8 β 2 , or equivalently x 6 t -8 in the physical space-time variables.) Observing that, thanks to the dispersion relation (3.50), iη 0 γ 2 γ 1 =η 2 0 (α l v l c 2 r + α r v r c 2 l ) ∈ (-∞, 0) , we readily find that Re(α 0 ) = -[v]α r η 2 0 (α l v l c 2 r + α r v r c 2 l )

2θrc 2 r |γ 1 | 2 (β + 2 -β - 2 ) 2 + v l vr η 2 µrη 2 0 -[v] 2 |γ 1 | 2 c 2 l η 2 v l (v l vr η 2 -η 2 0 ) vr µ l µr , Im(α 0 ) = 2 [v] α l η 0 η 2 0 α r |γ 2 | 2 (α l v l c 2 r + α r v r c 2 l ) θ l v l (β + 1 -β - 1 ) 2 + v l v r η 2 2µ r |γ 1 | 2 .
Since [v], α l,r , θ l,r , and µ l,r are all positive real numbers, we see that Im(α 0 ) is nonzero (it is of the same sign as [v]). Concerning Re(α 0 ), it is always nonzero for [v] > 0, which corresponds to an expansive phase transition (typically, vaporization), and it is also nonzero for -[v] > 0 and not too big. To evaluate We may observe that γ 2 α r + γ 1 α l = SV , where S = (-S 0 , S 0 ) with S 0 = (0, 0 * d-1 , -1). Therefore, Hunter's stability condition (2.36) for phase boundaries is equivalent to a(Re(V )) and a(Im(V )) being real, with a(R) :=

a s 1 (1, 0±) = lim
× (L + j ) * β + j -β - p (Q 2 -iβ - p P 2 )(R, γ p R - p ) -iβ - p SR η 0 [ρ] N (γ p R - p ) - i α 0 σ * e P 2 (V, R) + SV 2η 0 [ρ] N (R) + SR 2η 0 [ρ]
N (V ) .

(3.63)

✷

The condition (3.63) can be tested numerically. We present on Figure 3.4 numerical results in a realistic situation, which show that (3.63) is not satisfied. This might explain why surface waves are hardly ever observed in liquid-vapor flows. )/Re(a(Re(V ))) in terms of the mass transfer flux j for phase transitions in water at T = 600K (thermodynamic coefficients taken from [START_REF] Weast | CRC Handbook of Chemistry and Physics[END_REF]).

  defines an L ∞ function by the Cauchy-Schwarz inequality, whose inverse Fourier transform, Q[v], is a tempered distribution. If moreover v belongs to the Schwartz class, Q[v] is a function, explicitly given in terms of K by(2.31).✷In[START_REF] Hunter | Nonlinear surface waves[END_REF], Hunter had pointed out the following stability condition for equations of the form (2.34) with Λ satisfying (2.35), Λ(1, 0+) = Λ(1, 0-) .

+∞ 0 L 3 j=1σγ 2 θ l α l v l (β + 1 -β - 1 ) 2 + γ 2 γ 1 iη 0 θ r c 2 r

 03211212 (k, z) m 1 (k, z, τ ) dz + σ(k) g 1 (k, τ ) = 0 , (3.57) with m 1 := F [M(u, η; ρ, , ∂ ξ χ) • ( ρ, )] , g 1 := F[ g(u, η; ρ, , ∂ τ χ, ∂ ξ χ)] , L(k, z) := * C(u, η) R + j (L + j ) * Ȃd R + j e -kβ + j z (L + j ) * , k > 0 , L(k, z) := L(-k, z) , k < 0 , σ(k) := σ * , k > 0 , and σ(k) := σ(k) , k < 0 .(The fact that the sum is limited to j ≤ 3 comes from σ * C(u, η)R + j = 0 for j ≥ 4, see (3.53).)Substituting w(k, τ ) r(k, z) for F[( ρ, )](k, z, τ ), and w(k, τ ) s(k) for F[ χ](k) in the definition of m 1 and g 1 , we can rewrite (3.57) asa 0 (k) ∂ τ w(k, τ ) + +∞ -∞ a s 1 (kℓ, ℓ) w(kℓ, τ ) w(ℓ, τ ) dℓ = 0 , z) r(k, z) dz + s(k) σ(k) e(u, η)e(u) , + ℓ, z) m s (u, η; k, ℓ, z) dz + σ(k + ℓ) e(u, η) γ(u, η; k, ℓ) ,(3.59)is a negative real number (because c 2 r > v 2 r ). Therefore, we find that the first two terms in (3.62) equal to 2 [v] α r -γ 1

ℓ→0± +∞ 0 L 0 L( 1 L( 1

 0011 (1+ℓ, z) m s (u, η; 1, ℓ, z) dz + σ * e(u, η) γ(u, η; 1, 0±) , we go back to the definitions (3.60) (3.61) of m s and γ, and also to the definition of s (seeProposition 3.3). We thus see that limℓ→0+ +∞ + ℓ, z) m s (u, η; 1, ℓ, z) dz = σ * CR + j 4π (L + j ) * Ȃd R + j × (L + j ) * β + jβ - p iQ 2 (V, γ p R - p ) + P 2 (V, γ p β - p R - p ) + i γ 2 α r + γ 1 α l η 0 [ρ] N (γ p β - p R - p ) + ℓ, z) m s (u, η; 1, ℓ, z) dz = σ * CR + j 4π (L + j ) * Ȃd R + j × (L + j ) * β + jβ - p iQ 2 (V , γ p R - p ) + P 2 (V , γ p β - p R - p ) + i γ 2 α r + γ 1 α l η 0 [ρ] N (γ p β - p R - p ) ,with, as before,V = γ 1 R - 1 + γ 2 R - 2 , and 2π γ(u, η; 1, 0+) = P 2 (V, V ) + γ 2 α r + γ 1 α l η 0 [ρ] N (V ) , 2π γ(u, η; 1, 0-) = P 2 (V, V ) + γ 2 α r + γ 1 α l 2η 0 [ρ] N (V ) + γ 2 α r + γ 1 α l 2η 0 [ρ] N (V ) ,

3 j=1 2 p=1σ

 32 * CR + j 2α 0 (L + j ) * Ȃd R + j

Figure 1 :

 1 Figure1: Ratio Im(a(Re(V )))/Re(a(Re(V ))) in terms of the mass transfer flux j for phase transitions in water at T = 600K (thermodynamic coefficients taken from[START_REF] Weast | CRC Handbook of Chemistry and Physics[END_REF]).

  3.4 Under the Assumptions of Lemma 3.1, weakly nonlinear surface waves for the nonlinear model (3.37) (3.39) are governed by a nonlocal amplitude equation of the form (2.30), where Q is related by (2.31) to K

Acknowledgment: The authors thank the anonymous referee for his/her very careful reading of the manuscript. The second author was supported by University of L'Aquila, by INdAM, and by the European network HYKE under contract HPRN-CT-2002-00282.

Instead, we have

(3.49) In addition, we do have

Lemma 3.1 We assume that, as described above, the states t (ρ l , 0, . . . , 0, v l ) and t (ρ r , 0, . . . , 0, v r ) are connected by a planar dynamical subsonic phase boundary located at x d = 0. Without loss of generality we assume that the velocities v l and v r are positive. The associated right eigenvectors are defined as in (3.47). Then a linear combination of the form

solves the linearized jump conditions in (3.43) for some

such that the frequency η 0 = 0 and the wave vector η satisfy (3.45), if and only if,

where α l,r are defined as in (3.46), γ i = 0 for i ≥ 3, and

Proof. This is part of the main result in [START_REF] Benzoni-Gavage | Stability of multi-dimensional phase transitions in a van der Waals fluid[END_REF], in different variables though.

Let us give a sketch of computations. First recall that the Rankine-Hugoniot conditions (3.41) imply that

where g(ρ) := f (ρ) + p(ρ) ρ (which corresponds to the chemical potential of the fluid). These jump relations will enable us to eliminate χ = X e i(η 0 t + η•y) 4π m s (u, η; k, ℓ, z) := i(k + ℓ) Q 2 (u, η)( r(k, z), r(ℓ, z))

where Q 2 (u, η), P 2 (u), and P 2 (u) denote the symmetric bilinear mappings associated with the quadratic mappings Q(u, η), P(u), and P (u) respectively, and

(unlike what we did in the abstract framework of Section 2.3 we do not insert the matrix Ȋ2n here). By computation we find that a 0 (k) = α 0 /k for k > 0, and a 0 (k) = α 0 /k for k > 0, with

σ * e(u, η)e(u) .

(3.62) (We have used here the observation that (L + 2 ) * R - 1 = 0 and (L + j ) * R - 2 = 0 for j = 1 or j ≥ 3, which is an obvious consequence of the 'block form' of these vectors.) To check whether the number α 0 is nonzero we recall from (3.46), (3.49), (3.51), and (3.53) the values of β ± j , γ p , (L + j ) * Ȃd R + j and σ * C(u, η) respectively. In particular, we observe that

, and thus

In addition, going back to the definitions (3.47) and (3.48), we easily compute that

which is a positive real number (because c 2 l > v 2 l ), and similarly,


Afaf Bouharguane

Angelo Iollo

Lisl Weynans

Numerical solution of the Monge-Kantorovich problem by Picard iterations

Keywords: Optimal transport, Monge-Kantorovich problem, Numerical methods, Picard iterations Transport Optimal, Problème de Monge-Kantorovich, Méthodes numériques, Itérations de Picard

We present an iterative method to numerically solve the L 2 Monge-Kantorovich problem. The method is based on a Picard fixed point iteration of the linearized problem. Examples relative to the transport of two-dimensional densities show that the present method can significantly reduce the computational time over existing methods, especially when the Wasserstein distance between the densities is small.

Introduction

Optimal transportation is increasingly used to model problems in mechanics, physics, image analysis and other fields, see e.g. [START_REF] Villani | Optimal Transport, old and new[END_REF] and references therein. Because of all these applications, this old topic first introduced by Monge in 1781 [START_REF] Monge | Memoire sur la Theorie des Déblais et des Remblais[END_REF], has attracted a lot of attention these last years especially from a numerical point of view [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF], [START_REF] Angenent | Minimizing flows for the Monge-Kantorovich problem[END_REF], [START_REF] Dean | Numerical methods for fully nonlinear elliptic equations of the Monge Ampere type[END_REF], [START_REF] Iollo | A lagrangian scheme for the solution of the optimal mass transfer problem[END_REF], [START_REF] Benamou | Numerical solution of the second boundary value problem for the Elliptic Monge-Ampère equation[END_REF], [START_REF] Papadakis | Optimal Transport with Proximal Splitting[END_REF]. Indeed, compared to the theoretical results already obtained, the discrete solution of this problem still poses challenging problems in terms of computational burden and accuracy.

In this work we focus on the numerical solution of the L 2 Monge-Kantorovich problem (MKP) defined as follows. Let ρ 0 (ξ), ρ 1 (x) be two non-negative scalar density functions with compact support Ω 0 and Ω 1 , where ξ, x ∈ R d and d is the space dimension. We assume that Ω0 ρ 0 (ξ)dξ = Ω1 ρ 1 (x)dx.

Let X : Ω 0 → Ω 1 a smooth one-to-one map such that X(ξ) realizes the transfer of ρ 0 onto ρ 1 , i.e., a map that satisfies the following Jacobian equation:

ρ 0 (ξ) = det(∇X(ξ))ρ 1 (X(ξ)).
This equation is underdetermined with respect to X(ξ) and a solution is selected among all possible maps by introducing the following L 2 Kantorovich-Wasserstein distance:

inf X Ω0 ρ 0 (ξ)|X(ξ) -ξ| 2 dξ.
The L 2 MKP corresponds to finding a map X * such that this infimum is achieved. It has been proved that this problem admits a unique solution [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal Transport, old and new[END_REF], which is the gradient of a.e. convex function Ψ : Ω 0 → R:

X * (ξ) = ∇Ψ(ξ).

In the scientific computing literature, there exist two main approaches to approximate this problem. The first one is based on a direct solution of the Monge-Ampere equation (MAE):

ρ 0 (ξ) = det(∇ 2 Ψ(x))ρ 1 (∇Ψ(x)).

The difficulty in approaching the problem by this way is that the boundary conditions of this equation are not known a priori. Instead, the solution must verify the constraint X(Ω 0 ) = Ω 1 .

For the solution of the MAE with Dirichlet b.c. a recent numerical study in two dimensions is discussed in [START_REF] Dean | Numerical methods for fully nonlinear elliptic equations of the Monge Ampere type[END_REF], where the solution is obtained via a least-square formulation. The full MKP solution via a MAE was considered in [START_REF] Benamou | Numerical solution of the second boundary value problem for the Elliptic Monge-Ampère equation[END_REF]. The numerical method employed is based on the solution of MAE with boundary conditions that are iteratively updated to converge to the MKP solution.

Another class of methods relies on ideas from continuum mechanics. Benamou and Brenier (BB) numerically solved the MKP by using an augmented Lagrangian method [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]. In their formulation a temporal dimension is introduced so that, given Π:

[0, 1]×Ω 0 → R d , with Π(0, ξ) = ξ, Π(1, ξ) = X(ξ), x = Π(t, ξ) and ∂ t Π = v(t, x), the MKP amounts to the solution of inf ρ,v R d ρ(t, x)|v(t, x)| 2 dx,
where the infimum is taken among all densities ρ(t, x) ≥ 0 and velocity fields v(t, x) ∈ R d satisfying the continuity equation

∂ t ρ + ∇ • (ρv) = 0,
and the initial and final conditions:

ρ(0, •) = ρ 0 , ρ 1 (1, •) = ρ 1 .
The BB method results in a robust and viable discrete minimization problem under constraints which admits a unique solution. However, since it is a gradient method in space-time, the computational cost may be relevant. Also, numerical diffusion of the transported densities is observed in the simulation of the transport for intermediate times. A recent improvement of the minimization method at the base of the BB algorithm is proposed in [START_REF] Papadakis | Optimal Transport with Proximal Splitting[END_REF]. Using proximal splitting schemes the authors were able to solve difficult transport problems in presence of geometric constraints.

A yet different approach is devised in [START_REF] Angenent | Minimizing flows for the Monge-Kantorovich problem[END_REF]. The idea is to first consider a mass preserving mapping, not necessary optimal, between the initial and final distributions and then to solve a PDE up to steady state in order to rearrange the non-optimal mapping into an optimal one. It is shown that also this approach leads to a gradient-based minimization problem for which many gradient steps are usually needed to converge.

Here, we investigate two computationally efficient algorithms, easily implemented, to solve the optimal mass transfer problem. The paper is organised as follows. In the next section we explain how to get a relevant initialisation of the mapping and in section 3, we present the main idea of the numerical methods. Sections 4 and 5 are respectively devoted to the implementation of the Lagrangian and Eulerian methods and to numerical tests. Comparisons with the BB algorithm are presented in section 6, where a final discussion is provided.

The initial mapping

The initial mapping is found by a modification of the approach proposed in [START_REF] Angenent | Minimizing flows for the Monge-Kantorovich problem[END_REF]. As a first step, we compute a mapping that transports ρ 0 onto ρ 1 . This mapping will not necessarily be the gradient of a convex potential. Without loss of generality we assume that Ω = [-Λ, Λ] × [-Λ, Λ], Λ ∈ R. We assume also that the initial and final densities have compact support included in Ω. Let a : [-Λ, Λ] → R be the solution of

a ′ (ξ 1 ) Λ -Λ ρ 1 (a(ξ 1 ), x 2 ) dx 2 = Λ -Λ ρ 0 (ξ 1 , ξ 2 ) dξ 2 (1)
with x = X(ξ), ξ = (ξ 1 , ξ 2 ), x = (x 1 , x 2 ) and a(-Λ) = -Λ. Also let b : Ω → R be solution of the parametric ordinary differential equation

a ′ (ξ 1 ) ∂b(ξ 1 , ξ 2 ) ∂ξ 2 ρ 1 (a(ξ 1 ), b(ξ 1 , ξ 2 )) = ρ 0 (ξ 1 , ξ 2 ) (2) with b(ξ 1 , -Λ) = -Λ, ∀ξ 1 ∈ [-Λ, Λ]. The mapping X ω = (a(ξ 1 ), b(ξ 1 , ξ 2
)) takes ρ 0 onto ρ 1 by construction. Other choices of a(ξ) respecting monotonicity, regularity and appropriate initial conditions are possible. The mapping X ω = a(ξ 1 ) b(ξ 1 , ξ 2 ) is not in general the gradient of a complex potential. This requirement is a necessary condition for optimality. However, this mapping can always be decomposed as the sum of an irrotational field and a solenoidal field:

X ω = ∇ ξ Ψ i + ∇ × A. (3)
Hence, the actual initial mapping ∇ ξ Ψ i is found by solving ∇×∇×A = ∇×X ω with homogeneous Neumann boundary conditions, and computing ∇ ξ Ψ i = X ω -∇ × A. In two space dimensions, this amounts to the solution of one Laplace equation.

Inria 3 Picard iteration

The main idea of solution relies on the assumption that we dispose of an initial mapping that is a perturbation of the optimal mapping. Then, we derive a linear equation for the perturbation that is used to iteratively improve the initial guess.

Linearization

Let us assume that the optimal mapping is

X o (ξ) = ∇ ξ Ψ o
and that the mapping obtained at iteration n is

X n (ξ) = ∇ ξ Ψ i = ∇ ξ Ψ o + ∇ ξ Ψ n ǫ ,
where the error Ψ n ǫ satisfies Ψ n ǫ ≈ ǫ. We define ρ n 0 (ξ) as the initial density that mapped by X n (ξ) gives the exact final density ρ 1 (x). Then, taking a first-order Taylor expansion, we have

ρ n 0 (ξ) := ρ 1 (X n (ξ)) det (∇ ξ X n (ξ)) = ρ 1 (X o ) det (∇ ξ X o ) =ρ0(ξ) + det (∇ ξ X o )ρ 1 (X o )Tr (∇ ξ X o ) -1 ∇ 2 ξ Ψ n ǫ + det (∇ ξ X o )∇ x ρ 1 (X o ) • ∇ ξ Ψ n ǫ + o(ǫ),
where Tr designs the matrix trace operator. At first order in ǫ, we have

ρ n 0 (ξ) -ρ 0 (ξ) det (∇ ξ X o ) = ρ 1 (X o )Tr (∇ ξ X o ) -1 ∇ 2 ξ Ψ n ǫ + ∇ x ρ 1 (X o ) • ∇ ξ Ψ n ǫ
and finally, using the assumption that Ψ n ǫ ≈ ǫ, we replace X o by X n so that

ρ n 0 (ξ) -ρ 0 (ξ) det (∇ ξ X n ) = ρ 1 (X n )Tr (∇ ξ X n ) -1 ∇ 2 ξ Ψ n ǫ + ∇ x ρ 1 (X n ) • ∇ ξ Ψ n ǫ . (4)
The mapping update is then found by the iteration

X n+1 = X n -α ∇Ψ n ǫ ,
that converges to X o as a geometric series, for α ∈ [0, 1], in the linearized approximation.

In particular, if the optimal mapping is a perturbation of identity, i.e.,

X o (ξ) = ξ + ∇ ξ Ψ o with ∇ ξ Ψ o = o(ǫ)
, then at first order in ǫ equation ( 4) reduces to

ρ n 0 (ξ) -ρ 0 (ξ) = ∇ ξ • (ρ 1 (X n (ξ))∇ ξ Ψ n ǫ ) (5)
that is equivalent to a semidiscretization in time of the continuity equation written in Eulerian form.

Regularization and continuation

The above equations are well defined only on the support of the initial and final densities. Therefore, in actual computations the initial and final densities are regularized by adding a small constant ζ so that the leading order differential operator is well defined on the entire domain:

(ρ 1 (X n ) + ζ) Tr (∇ ξ X n ) -1 ∇ 2 ξ Ψ n ǫ with ζ ≈ 10 -3 in practice.
The initial mapping determined as explained in section 2 may in general be far from the optimal mapping. This implies that the perturbation equation ( 4) is not a good approximation of the original problem. In those cases it possible to solve the problem by continuation adding to the initial and final densities a constant. This constant is then iteratively brought to O(10 -3 ). Taking ζ k ∈ R + large enough, we define

ρ k 0 = ρ 0 + ζ k and ρ k 1 = ρ 1 + ζ k so that the condition Ψ ǫ ≈ ǫ is verified. Once the optimal mapping Ψ k o for such a ζ k is determined, we define ζ k+1 = (1 -β)ζ k with β ∈ [0, 1]
, and Ψ k o as the initial mapping for the next optimal transport problem. In practice, we choose β = 0.1. In order to give a first estimate of ζ, it is natural to assume that it is proportional to the Wasserstein distance between the two densities. However, this distance in unknown a priori. Therefore, we approximate the Wasserstein distance using the initial mapping computed as explained in 2. We take

ζ ∝ ρ 0 |∇ ξ Ψ i -ξ| 2 dx.
For example, for test case TC3 (see section 5) we took ζ = 1 as the Wasserstein distance approximation above is of the same.

Numerical solution

We detail two solution approaches. A Lagrangian approach, that can readily be coded as it implies the solution of a rather simple elliptic PDE and an Eulerian approach, that can more easily be applied in more general geometries.

Lagrangian approach

According to previous sections, the following iterative method (Algorithm 1) is considered. In the limits of the first order perturbation analysis of this section, Algorithm 1 converges as a geometric series, i.e.,

Ψ n ǫ ≈ Ψ o (1 -(1 -α) n ) . (6)
Therefore, defining the normalized residual as r = Ψ n ǫ -Ψ o / Ψ o , we have:

log r = n log (1 -α), (7)
where log(1 -α) is the rate of convergence of the iterative scheme.

To approximate the problem, we discretize equation ( 4) by a standard second-order finite difference scheme on a Cartesian grid. At the discrete level, we impose homogeneous Dirichlet Inria Algorithm 1 Lagrangian iterative algorithm 1. n = 0; 2. compute X 0 = ∇Ψ i as explained in section 2;

3. ρ n 0 (ξ) = ρ 1 (X n (ξ)) det ∇ ξ X n (ξ);
4. compute Ψ n ǫ by solving equation (4); 5. X n+1 = X n -α ∇Ψ n ǫ ; 6. n = n + 1;

7. goto 3 if convergence is not attained; boundary conditions to equation [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF]. In practice, since the initial and final densities have compact support, the boundary conditions imposed at the borders of the computational domain do not significantly affect the solution as ζ → 0. Mass transportation is performed by the continuity equation in Lagrangian coordinates.

In the next sub-section we will provide an alternative approach where the continuity equation is integrated in Eulerian coordinates. For numerical examples, see section 5, we will assume that the convergence is attained when a threshold criterion on ρ n 0 -ρ 0 ∞ is satisfied.

Eulerian approach

In the Eulerian framework, the solution of the L 2 MKP consists in determining the initial velocity u(0, .) = u 0 satisfying:

     ∂ t (ρu) + ∇ • (ρu) = 0 (8)
∂ t (ρu) + ∇ • (ρu ⊗ u) = 0 (9) ρ(0, •) = ρ 0 , ρ(1, •) = ρ 1 (10)
The optimal initial velocity u 0 is written under the form u 0 = ∇Ψ, where the potential Ψ is the Lagrange multiplier of constraints ( 8) and [START_REF] Villani | Optimal Transport, old and new[END_REF]. We propose an alternative solution method for u 0 by adapting Algorithm 1 to the Eulerian framework.

Starting from an initial guess for u 0 , we numerically solve equations ( 8)-( 9) using ρ(0, .) = ρ 0 as the initial condition for density. Solving system (8)-( 9) as detailed hereafter, we get a final density ρ1 and a final velocity field that we denote ũ1 . Since in general the initial velocity guess does not correspond to the optimal solution, the density ρ1 is different from ρ 1 . The next step is to find a potential velocity field transporting the data ρ1 on ρ 1 . This velocity field is considered as a correction to the initial mapping.

To find the additional velocity field, we approach the continuity equation by the following elliptic problem:

∇ • ρ 1 + ρ1 2 ∇ Ψ = ρ1 -ρ 1 .
This problem is solved using classical second-order centered finite-differences. From the potential Ψ we compute a velocity field ũ = ∇ Ψ, which can be seen as the velocity field necessary to advect ρ1 to ρ 1 . As we want to use it as a correction to u 0 , we have to advect it backward in time to make it match with the initial density ρ 0 . We therefore perform an integration backward in time of the system (8)-( 9) with final conditions ρ1 and ũ1 , obtained by the forward numerical computation. The result of the backward advection of ũ, that we denote ũn 0 , is added to u 0 as a corrective term. The full method is detailed in Algorithm 2.

The numerical resolution of equations ( 8)-( 9) is performed using a remeshed particle method, similar to the one developed in [START_REF] Weynans | Consistency, accuracy and entropy behaviour of remeshed particle methods[END_REF] for compressible Euler equations. In this class of methods, the fluid is discretized on small masses concentrated on points. These points, the particles, are displaced in a Lagrangian way. New particles uniformly distributed are created at regular time intervals by an interpolation of the values of the existing particles, what is usually called remeshing the particles. The remeshing step creates new particles in a conservative way, by distributing the quantities carried by the particles at the nodes of an underlying grid.

Because the particles themselves are moved in a Lagrangian way, remeshed particle methods are submitted to less restrictive stability conditions, in the context of advection problems, than more classical grid-based methods such as finite-differences or finite-volume methods. Indeed their stability condition is typically proportional to the inverse of the velocity gradient (meaning that the particles trajectories do not cross), instead of a classical CFL condition. In the context of the method that we present here for optimal mass transport, we need to compare the numerical result computed at the final time to the final exact density ρ 1 . Therefore, the use of larger time steps is an advantage because it means that the final state is computed with less numerical dissipation due to the temporal integration than in the case of more classical grid-based methods.

Algorithm 2 Eulerian iterative algorithm 1. n = 0; 2. initialize u 0 ; 3. for initial data ρ 0 and u 0 = X n , compute ρ n 1 and u n 1 by numerically solving system (8)-( 9); 4. compute Ψ by solving equation:

∇ • ρ 1 + ρ1 2 ∇ Ψ = ρ1 -ρ 1 5
. integrate backward in time system ( 8)-( 9) with final conditions ρ1 = ∇ Ψ and ũ1 to get at initial time the velocity correction ũn 0 ;

6. X n+1 = X n + α ũn 0 ;

7. n = n + 1;

8. goto 3 if convergence is not attained;

Numerical illustrations

In this section, we provide several numerical validations of our methods. We perform the following test cases (TC):

• TC1. We consider two density distributions having a non-negligible support intersection.

Inria

We take two Gaussians of unit mass, same variance and displaced of 0.2 in the vertical direction. The final density has a cross correlation of 0.5, the initial is isotropic. We have:

ρ 0 (ξ) = 3.97887e -12.5(ξ 2 1 +(ξ2+0.1) 2 ) ρ 1 (x) = 4
.59441 e -16.66666(x 2 1 +x1(x2-0.1)+(x2-0.1) 2 ) .

• TC2. This is similar to TC1, but the centroid of the initial density is further displaced. We take the initial and final density distributions:

ρ 0 (ξ) = 1
.98944 e -12.5((ξ1-0.2) 2 +(ξ2+0.1) 2 ) + e -12.5((ξ1+0.2) 2 +(ξ2+0.1) 2 )

ρ 1 (x) = 3
.97887 e -12.5(x 2 +y 2 ) .

• TC3. In this case there is mass separation and virtually no intersection between the initial and final density support. The initial and final density are ρ 0 (ξ) = 1.98944 e -12.5((ξ1-0.5) 2 +(ξ2+0.3) 2 ) + e -12.5((ξ1+0.5) 2 +(ξ2+0.3) 2 )

ρ 1 (x) = 3
.97887 e -12.5(x 2 +y 2 ) .

• TC4. A final case showing the robustness of the Picard iteration is presented. The mapping between two scans of 168 2 pixels relative to the abdomen of a breathing patient, see figure 9, is determined.

Test case TC1, TC2, TC3 are performed using both the Lagrangian and the Eulerian schemes. We consider a square domain with Λ = 2 discretized using a uniform 200 × 200 grid. Compared to the Lagangian case, the corresponding Eulerian TCs are rotated of π/4 in order to show the accuracy of the particle remeshing in a transverse direction with respect to the grid.

For TC1, with α = 0.2, the expected rate of convergence is of 0.09691, see [START_REF] Monge | Memoire sur la Theorie des Déblais et des Remblais[END_REF]. The fitting of the convergence curve with a straight line gives a convergence rate of 0.0958315 for the algorithm 1, see figure 1. The relative error in the max norm is of the order of 10 -3 after 23 iterations. The Eulerian solution is given in figure 2.

In the next test case (TC2), the expected rate of convergence is of .0457575 (α = 0.1), the fitting of the convergence curve with a straight line gives a convergence rate of 0.444000, see figure 3. The Eulerian solution is given in figure 4.

Finally, in TC3 the Lagrangian scheme was initialized with a regularization ζ = 1 and in 40 steps it was reduced to 10 -3 . The initial condition along with the optimal mapping stream curves are shown in figure 5. The intermediate and final solution and error are shown in figures 6 and 7. The Eulerian solution is given in figure 8.

The last case presented, TC4, is relative to front abdomen sections of a breathing patient, figure 9 (courtesy of B.D. de Senneville). The figure shows the initial velocity field relative to the optimal transport between the two subsequent scans. In certain focalized beam therapies, like focalized ultrasounds, it is crucial to accurately predict the movement of a patient in order to calibrate the displacement of the targeted region. In applications, real-time optical flow techniques based on heuristic arguments are employed. The methods presented here are not real time, although the solution of the MKP took a few seconds in this case. However, the solution of the MKP to determine the displacement between two images offers the advantage over optical flow techniques of being objective in the sense that it is based on a clearly identified model. The Picard iteration approach presented in this paper is suitable for transports where the distance travelled by the elementary masses is small compared to the characteristic length of the density distributions. The reason is that the unknowns of this approach are the initial transport velocities. If the distance travelled becomes large, small deviations of the initial condition lead to large errors. The BB approach has the advantage of solving for a space-time saddle point. Hence the trajectories are not identically straight line, as in the case of the present method, and the problem remains well conditioned also for large mass displacements as shown in [START_REF] Papadakis | Optimal Transport with Proximal Splitting[END_REF]. This is of course at the price of a costly computational problem in d + 1 dimensions.

In order to show the computational advantage of the Picard iteration for suitable cases, we compare the CPU time to solve the MKP by the BB method, to the CPU time of algorithms presented here. Since the convergence criteria are different for the two methods, we will determine the CPU time so that the initial criteria will be divided by 30, 60 and 100.

We recall that the convergence criterion of the BB method is based on the following residual of the Hamilton-Jacobi equation (see [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]):

res n = ∂ t φ n + |∇φ n | 2 ,
where φ is the Lagrange multiplier of constraints ( 8) and [START_REF] Villani | Optimal Transport, old and new[END_REF]. The convergence criterion is given by:

crit n b = 1 0 Ω ρ n |res n | 1 0 Ω ρ n |∇φ n | 2 .
For the Picard iteration method presented here, we consider the error between the given final distribution ρ 1 and the final distribution ρ n 1 at iteration n:

crit n = ||ρ 1 -ρ n 1 || ∞
in the Eulerian case. For the Lagrangian case we consider the equivalent error based on the initial density. We consider TC1. For the BB algorithm, we discretize the time domain using 32 nodes and space with 200 × 200 grid points. We take optimal numerical parameters in the Uzawa iteration in order to converge as quickly as possible. 1 shows the evolution of the CPU time for different convergence thresholds. We can see that CPU time is very small for the present methods compared to the BB method. Note also that it is impossible to reduce the initial criteria by 100 in a reasonable time for BB. This is the reason we consider a convergence IC/60 < IC ≤ IC/100 in the table above. Usually, in order to get to convergence of the BB algorithm (IC/100), we have about 3000 Uzawa iterations with a computing time larger than three hours.

BB method

Cases where the Wasserstein distance between the densities is small are particularly favorable for the Picard fixed point solution because there is no need to regularize the solution by adding an initial constant ζ to both the initial and final distributions. In cases where the Wasserstein distance is larger, like TC3, the computational advantage of the Picard iteration with respect to the BB scheme is reduced (initially, ζ = 1). For example, in TC3 the solution is obtained with the Lagrangian scheme in about 5 min CPU time, where the BB scheme takes 20 min to get to IC/60 < IC ≤ IC/100. Still , the Picard iteration is more advantageous.

In conclusion, the Picard iteration that we presented can lead to a significant improvement in terms of convergence rate over the BB scheme, in particular when the Wasserstein distance is small. This can be a significant advantage for large three-dimensional problems in imagery and computational mechanics.

Figure 1 :

 1 Figure 1: TC1. Error after 23 iterations and convergence for algorithm 1

Figure 2 :

 2 Figure 2: TC1. Top picture row: plot of density isolines ρ(t, x) for t = 0, 0.25, 0.5, 0.75, 1 along with the optimal path computed by algorithm 2 . Bottom: red (resp. blue) line represents the error in L 1 norm (resp. L ∞ norm).

Figure 3 :

 3 Figure 3: TC2: initial distribution, final distribution and convergence. Isolines from 0 to the max spaced of 0.2.

Figure 4 :

 4 Figure 4: TC2. Top picture row: plot of the isolines of the density ρ(t, x) for t = 0, 0.25, 0.5, 0.75, 1. Bottom: convergence of algorithm 2. Red (resp. blue) data represent the log 10 (||ρ n 1 -ρ 1 || 1 ) (resp. log 10 (||ρ n 1 -ρ 1 || ∞ )) as a function of the number of iterations

Figure 5 :

 5 Figure 5: TC3. Initial and final densities and stream curves corresponding to ∇ ξ Ψ o , the optimal mapping.

Figure 6 :InriaFigure 7 :

 67 Figure 6: TC3. Intermediate densities for X α = ξ + χ (∇ ξ Ψ o -ξ) with χ = 0, 1 4 , 1 2 , 3 4 from left to right, top to bottom. Isolines from 0 to the max spaced of 0.2.

Figure 8 :Figure 9 :

 89 Figure 8: TC3. Plot of the isolines of the density ρ(t, x) along the optimal path computed with algorithm 2

Table 1 :

 1 CPU time: BB method vs. Our methods Table

				Eulerian method	Lagrangian method
		crit n b	CPU time	crit n	CPU time	crit n	CPU time
	Initial Criteria (IC)	0.707107	0m0s	2.063231	0m0s	2.063231	0m0s
	≈ IC/30	0.023547	2m36s	0.06004	0m4.13s	0.060026	0m2.6s
	≈ IC/60	0.011784	5m12s	0.030851	0m5.32s	0.034312	0m3.2s
	IC/60 < IC ≤ IC/100 0.009841	19m31s	0.018421	0m6.22s	0.027887	0m3.74s

RR n°8477

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399