Laurence Cherfils
email: laurence.cherfils@univ-lr.fr

Hussein Fakih
email: hussein.fakih@math.univ-poitiers.fr

Alain ||

Miranville
email: alain.miranville@math.univ-poitiers.fr

A Cahn-Hilliard system with a fidelity term for color image inpainting

Keywords: Cahn-Hilliard system, fidelity term, color image inpainting, wellposedness, exponential attractor, simulations. AMS Subject Classification: 35B40, 35K55, 68M07, 80A23

In this paper, we propose a model for multi-color image inpainting composed of n colors. In particular, as in the binary model, i.e., the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation [4], we add a fidelity term to the corresponding Cahn-Hilliard system. We are interested in the study of the asymptotic behavior, in terms of finite-dimensional attractors, of the dynamical system associated with the problem. The main difficulty here is that we no longer have the conservation of mass, i.e., of the spatial average of the order parameter c, as in the Cahn-Hilliard system. Instead, we prove that the spatial average of c is dissipative. We finally give numerical simulations which confirm and extend previous ones on the efficiency of the binary model.

Introduction

Image inpainting involves filling in part of an image or video from the surrounding area. It is essentially some kind of interpolation. Its applications include restoration of old paintings by museum artists [START_REF] Emile-Male | The restorer's handbook of easel painting[END_REF], removing scratches from old photographs [START_REF] Braverman | Photoshop retouching handbook[END_REF], altering scenes in photographs [START_REF] King | The Commissar vanishes[END_REF], and restoration of motion pictures [START_REF] Kokaram | Motion Picture Restoration: Digital Algorithms for Artefact Suppression in Degraded Motion Picture Film and Video[END_REF].

The authors in [START_REF] Bertalmio | Image inpainting[END_REF] introduced image inpainting for digital image processing. Their model is based on nonlinear PDE's. A large number of PDE's and variational approaches have been considered in view of such applications, such as the total variation model [START_REF] Chan | Variationnal restoration of nonflat image features: Models and algorithms[END_REF][START_REF] Chan | Mathematical models of local nontexture inpaintings[END_REF][START_REF] Chan | Variationnal image inpainting[END_REF][START_REF] Dolcetta | Area-preserving curve-shortening flows: From phase separation to image processing[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Rudin | Total variation based image restoration with free local constraints[END_REF][START_REF] Taylor | Linking anisotropic sharp and diffuse surface motion laws via gradient flows[END_REF], the Euler elastica model [START_REF] Chan | Euler's elastica and cuvature-based inpainting[END_REF][START_REF] Chan | Total variation wavelet inpainting[END_REF][START_REF] Masnou | Level lines based disocclusion[END_REF], the active contour model based on the Mumford and Shah segmentation [START_REF] Tsai | Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification[END_REF], the inpainting scheme based on the Mumford-Shah-Euler image model [START_REF] Esedoglu | Digital inpainting based on the Mumford-Shah-Euler image model[END_REF], inpainting with the Navier-Stokes equations [START_REF] Bertalmio | Navier-Stokes, fluid dynamics, and image and video inpainting[END_REF], or wavelet-based inpainting [START_REF] Chan | Total variation wavelet inpainting[END_REF][START_REF] Dobrosotskaya | A Wavelet-Laplace variational technique for image deconvolution and inpainting[END_REF].

The authors of [START_REF] Bertozzi | Analysis of a two-scale Cahn-Hilliard model for binary image inpainting[END_REF] introduced the following Cahn-Hilliard model for binary image inpainting:

∂c ∂t = -∆(ε∆c - 1 ε f (c)) + λ(x)(h -c), (1.1)
where

λ(x) = 0 if x ∈D, λ 0 if x ∈Ω\D, h (
F (s) = s 2 (s -1) 2 .
Equation (1.1) is identical to the standard Cahn-Hilliard equation [START_REF] Bates | The dynamics of nucleation for the Cahn-Hilliard equation[END_REF][START_REF] Novick-Cohen | Nonlinear aspects of the Cahn-Hilliard equation[END_REF], except for the second term on the right-hand side. This term is added to keep the solution constructed close to the given image h(x) in the complement of the inpainting domain, where image information is available.

Well-posedness results for (1.1) have been obtained in [START_REF] Bertozzi | Analysis of a two-scale Cahn-Hilliard model for binary image inpainting[END_REF] (see also [START_REF] Burger | Cahn-Hilliard inpainting and a generalization for grayvalue images[END_REF] for the study of the stationary problem).

Furthermore, the asymptotic behavior, in terms of finite-dimensional attractors, of (1.1) has been studied in [START_REF] Cherfils | Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation in image inpainting[END_REF]. More precisely, such sets provide information on all the possible dynamics of the system and are expected to have a rich geometric structure; they contain in particular all steady states and heteroclinic orbits. Moreover, the finitedimensionality means, roughly speaking, that, even though the initial phase space is infinite-dimensional, the reduced dynamics can be characterized by a finite number of parameters. We refer the reader to, e.g., [START_REF] Efendiev | Exponential attractors for a singularly perturbed Cahn-Hilliard system[END_REF][START_REF] Miranville | Attractors for dissipative partial differential equations in bounded and unbounded domains[END_REF][START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF] for discussions on this subject.

Finally, the existence of local (in time) solutions to (1.1) with logarithmic nonlinear terms has been studied in [START_REF] Cherfils | On the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation with logarithmic nonlinear terms[END_REF] (note indeed that the original Cahn-Hilliard equation was actually proposed with thermodynamically relevant logarithmic nonlinear terms which follow from a mean-field model; regular (and, in particular, cubic) nonlinear terms are approximations of such logarithmic nonlinear terms). In that case, we can obtain better results than those obtained with polynomial nonlinear terms in [START_REF] Cherfils | Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation in image inpainting[END_REF], as far as the convergence time is concerned, in particular. We also note that double obstacle nonlinear terms give better results than those obtained with polynomial nonlinear terms, see [START_REF] Bosch | Fast solvers for Cahn-Hilliard inpainting[END_REF].

While the articles mentioned above only deal with binary images (i.e., black and white; see however [START_REF] Burger | Cahn-Hilliard inpainting and a generalization for grayvalue images[END_REF] for grayvalue images), it is natural and important to also address color images. In particular, one idea in this direction is to extend the approach of [START_REF] Bertozzi | Analysis of a two-scale Cahn-Hilliard model for binary image inpainting[END_REF] to multi-component Cahn-Hilliard systems, each phase corresponding to a given color. Multi-component Cahn-Hilliard systems have been proposed and studied in, e.g., [START_REF] Conti | Multi-component Cahn-Hilliard systems with dynamic boundary conditions[END_REF][START_REF] Elliott | A generalised diffusion equation for phase separation of a multicomponent mixture with interfacial free energy[END_REF][START_REF] Eyre | Systems of Cahn-Hilliard equations[END_REF][START_REF] Garcke | On anisotropic order parameter models for multi-phase systems and their sharp interface limits[END_REF][START_REF] Li | Global attractor for the Cahn-Hilliard system with fast growing nonlinearity[END_REF]. Our main aim in this article is to propose such a model for color image inpainting, based on a fidelity term as in [START_REF] Bertozzi | Analysis of a two-scale Cahn-Hilliard model for binary image inpainting[END_REF].

As far as the mathematical treatment of the problem is concerned, we have to face two essential difficulties. The first one, which we already encountered in the binary model in [START_REF] Cherfils | Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation in image inpainting[END_REF], is that, due to the presence of the fidelity term, we no longer have the conservation of mass, i.e., of the spatial average of the order parameter, contrary to the original Cahn-Hilliard system. This is overcome by deriving a global (in time) and dissipative (i.e., independent of time and bounded sets of initial data, at least for large times) estimate on the order parameter; in particular, a key step is to obtain such an estimate on the spatial average of the order parameter. We can note that this estimate is a crucial one, already in order to prove the existence of a solution. The second difficulty is that the order parameter has to satisfy some constraint, namely, the sum of the phases must be equal to one. As a consequence, we have to be careful in the choice of the functional spaces and in the (weak) formulation of the problem.

This article is organized as follows. In Section 2, we present the model constructed from the multi-phase Cahn-Hilliard system by adding a fidelity term. Furthermore, we give the assumptions on the nonlinear term. Then, in Section 3, we derive a priori estimates which allow us to prove, in Section 4, the existence and uniqueness of solutions, as well as the existence of the global attractor. In Section 5, we construct finite-dimensional attractors and, in Section 6, we prove that the model is algebraically consistent with the diphasic one. Finally, in Section 7, we give some numerical simulations which confirm that the one step algorithm with threshold proposed in [START_REF] Cherfils | Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation in image inpainting[END_REF] (see also [START_REF] Cherfils | On the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation with logarithmic nonlinear terms[END_REF]) is efficient, also in the context of multi-color inpainting.

Remark 1.1. After this work was completed, we discovered the recent preprint [START_REF] Bosch | A fractional inpainting model based on the vector-valued Cahn-Hilliard equation[END_REF] in which a similar model is considered for grayvalue inpainting. Note however that, there, the mathematical analysis of the problem is not addressed.

Proposed model and setting of the problem

Let c = (c 1 , c 2 , ...., c n ) be the phase variable and Ω be a bounded domain of R N (N = 2, 3) with a regular boundary Γ. We assume that n i=1 c i = 1 and 0 ≤ c i ≤ 1, i = 1, ..., n.

We define the hyperplane

S := {c ∈ R n such that n i=1 c i = 1}.
(2.1)

We postulate that the free energy can be written as follows:

F = Ω F (c) + ε 2 2 n i=1 |∇c i | 2 dx, (2.2)
where

F (c) = 1 n n i=1 c 2 i (1 -c i ) 2 .
The time evolution of c is governed by the gradient of the energy with respect to the H -1 -inner product under the additional constraint (2.1). This constraint has to hold everywhere and at every time. In order to ensure this constraint, we use a variable Lagrangian multiplier β(c) [START_REF] Garcke | On anisotropic order parameter models for multi-phase systems and their sharp interface limits[END_REF]. The time dependence of c i is given by the following system of Cahn-Hilliard equations describing each phase:

∂c i ∂t = ∆µ i , i = 1, ..., n, (2.3)
µ i = ∂F (c) ∂c i -ε 2 ∆c i + β(c), i = 1, ..., n. (2.4)
To compute β(c), we write the equation satisfied by C = c 1 + c 2 + ... + c n and we want C = 1 to be solution to this equation,

∂C ∂t = ∆ n i=1 ∂F ∂c i -ε 2 ∆C + nβ(c) . (2.5) Therefore, β(c) = -1 n n i=1
∂F ∂c i and we obtain

∂c i ∂t = ∆µ i , i = 1, ..., n, (2.6
)

µ i = ∂F (c) ∂c i -ε 2 ∆c i - 1 n n i=1 ∂F ∂c i , i = 1, ..., n. (2.7)
Finally, we need to add the fidelity term λ(x)(h i -c i ) to each equation to keep the solution constructed close to the given image h i (x) as in the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard model for binary image inpainting; more precisely, h i (x) is a part of the original image corresponding to the color representing the phase c i . We thus have the following model:

∂c i ∂t = ∆µ i + λ(x)(h i -c i ), i = 1, ..., n, (2.8)
µ i = ∂F (c) ∂c i -ε 2 ∆c i - 1 n n i=1 ∂F ∂c i , i = 1, ..., n, (2.9)
where

λ(x) = 0 if x ∈D, λ 0 if x ∈Ω\D,
D Ω is the inpainting domain, and h = (h 1 , ..., h n ) ∈ (L 2 (Ω)) n . Equations (2.8)-(2.9) are endowed with the Neumann boundary conditions

∂c i ∂ν = ∂µ i ∂ν = 0, on ∂Ω, i = 1, ..., n. (2.10)
We note that the additional fidelity term in (2.8) can be derived from a gradient flow under the L 2 -inner product for the energy

λ 0 2 n i=1 Ω\D (c i -h i ) 2 dx.
(2.11)

Our proposed model (2.8)-(2.9) can thus be thought of as the superposition of a gradient descent with respect to the H -1 -inner product for (2.2) and a gradient descent with respect to the L 2 -inner product for (2.11). However, it is not the gradient descent, neither in the H -1 -nor the L 2 -inner product, for the sum of the energies (2.2) and (2.11).

We again assume that n i=1 c i = 1, so that the additional fidelity term in (2.8) must also satisfy the constraint (2.1) everywhere and at every time,

n i=1 h i = 1, hence λ(x) n i=1 (c i -h i ) = 0, ∀t ≥ 0.
We rewrite (2.8)-(2.10) in the form

∂c ∂t = ∆µ + λ 0 χ Ω\D (x)(h -c), (2.12
)

µ = f (c) -ε 2 ∆c, (2.13
)

∂c ∂ν = ∂µ ∂ν = 0, (2.14)
where c = (c i ) i , µ = (µ i ) i , h = (h i ) i , and f (c) = (f i (c)) i , f i (c) = 2 n [c i (1 -c i ) 2 -c 2 i (1 -c i )] - 1 n n i=1 ∂F ∂c i , i = 1, ..., n.
We define the tangent space of S (defined in (2.1)) as

T S := {φ = (φ i ) i ∈ R n such that n i=1 φ i = 0}. (2.15)
Then, we introduce the projection

P : R n →T S φ →Pφ = φ - 1 n (φ.e)e,
where e = (1, ..., 1). Notice that, if φ ∈ S ∩ (L 2 (Ω)) n , then

P φ = φ.
We denote by φ the spatial average of a function φ in

(L 1 (Ω)) n , φ := 1 Vol(Ω) Ω φ(x)dx
and set c := c -c , c being the mean-free part of c. We further introduce the following spaces:

L = (H -1 (Ω)) n ∩ S, H = (L 2 (Ω)) n ∩ S, V = (H 1 (Ω)) n ∩ S, W = v ∈ (H 2 (Ω)) n ∩ S, ∂v ∂ν = 0 on Γ , Ḣ-1 (Ω) = {φ ∈ (H -1 (Ω)) n ∩ T S, φ, e (H -1 (Ω)) n ,(H 1 (Ω)) n = 0}, Ḣ = {φ ∈ (L 2 (Ω)) n ∩ T S, φ = 0}, V = {φ ∈ (H 1 (Ω)) n ∩ T S, φ = 0}, and Ẇ = {φ ∈ (H 2 (Ω)) n ∩ T S, φ = 0}.
We finally introduce the bilinear form

a(u, v) = ((∇u, ∇v))
which is continuous on V . Here and below, we denote by . the usual L 2 -norm (with associated scalar product ((., .))); in general, . X denotes the norm on the Banach space X. We can associate with a the linear operator A in V with domain D(A) = Ẇ . For u ∈ D(A), Au is defined by

((Au, v)) = a(u, v), ∀v ∈ V ,
which amounts to saying that Au = -∆u. It is easy to show that A is an isomorphism from V onto V (the dual of V ) and from D(A) onto Ḣ. Furthermore, in light of the compact injection of (H 1 (Ω)) n in (L 2 (Ω)) n (and also of V in Ḣ), it follows that A -1 is self-adjoint and compact. There thus exists an orthonormal basis of Ḣ associated with the eigenvalues λ j , j ≥ 1, of A,

         λ 1 = inf v∈ V \{0} v 2 V v 2 , AN j = λ j N j , N j ∈ D(A), j ≥ 1, 0 < λ 1 ≤ λ 2 ...,
The family {N j } j may be assumed to be normalized in the norm of Ḣ, i.e.,

((N i , N j )) = δ ij ,
where

δ ij = 1 if i = j, 0 otherwise.
We now rewrite problem (2.12)-(2.13) in the form

∂c ∂t + APµ + λ 0 χ Ω\D (x)(c -h) = 0, (2.16
)

Pµ = Pf (c) + ε 2 Ac. (2.17)
For simplicity, we set ε = 1 and have

∂c ∂t + APµ + λ 0 χ Ω\D (x)(c -h) = 0, (2.18
)

Pµ = Pf (c) + Ac. (2.19)
Integrating (2.18) over Ω, we obtain

d c dt = -λ 0 χ Ω\D (x)(c -h) (2.20)
We can thus rewrite (2.18)-(2.19) in the form

∂c ∂t + APµ + λ 0 χ Ω\D (x)(c -h) = 0, (2.21
)

Pµ = Pf (c) + Ac. (2.22)
Finally, we can reformulate the problem as follows:

∂A -1 c ∂t + Pµ + λ 0 A -1 (χ Ω\D (x)(c -h)) = 0, (2.23
)

Pµ = Pf (c) + Ac. (2.24)
As far as the function f is concerned, we assume that f ∈ C 2 (R n ) and

((f (c)v, v)) ≥ -ξ 0 v 2 , ξ 0 ≥ 0, ∀c ∈ R n , ∀v ∈ T S, (2.25)
((f (c), c)) ≥ α Ω (c 2 . c2 + c 2 . c2 )dx, α > 0, ∀c ∈ S, (2.26)
where

ϕ 2 = (ϕ 2 i ) i , ϕ = (ϕ i ) i ∈ R n , and f (c) = (∇f i (c)) i denotes the Jacobian matrix of f , i = 1, ..., n.
Remark 2.1. In particular, assumptions (2.25) and (2.26) are satisfied by our proposed model (2.8)-(2.9). Indeed, we recall that

f i (c) = 2 n [c i (1 -c i ) 2 -c 2 i (1 -c i )] - 1 n n i=1
∂F ∂c i , i = 1, ..., n. We thus have

∂f i ∂c i (c) = 2 n [6c 2 i -6c i + 1] - 2 n 2 [6c 2 i -6c i + 1], i = 1, ..., n,
and ∂f i ∂c j (c) = - 2 n 2 [6c 2 j -6c j + 1], i, j = 1, ..., n, j = i. Therefore, the i-th component of f (c)v reads (f (c)v) i = 2 n [6c 2 i -6c i + 1]v i - 2 n 2 n j=1 [6c 2 j -6c j + 1]v j , i = 1, ..., n, hence
((f (c)v, v)) = 2 n Ω n i=1 [6c 2 i -6c i + 1]v 2 i dx - 2 n 2 Ω n j=1 [6c 2 j -6c j + 1]v j n i=1 v i dx = 2 n Ω n i=1 [6c 2 i -6c i + 1]v 2 i dx, since v ∈ T S. It thus follows that ((f (c)v, v)) ≥ - n i=1 ξ i v i 2 ,
where ξ i > 0, i = 1, ..., n. On the other hand, we have

((f (c), c)) = ((f (c), c)) = Ω n i=1 f i (c)c i dx = 2 n Ω n i=1 [2c 3 i -3c 2 i + c i ]c i dx - 1 n Ω n i=1 ∂F ∂c i n i=1 ci = 2 n Ω n i=1 [2c 3 i -3c 2 i + c i ]c i dx = 2 n Ω n i=1 [2c 3 i -3c 2 i + c i ]c i dx,
for all c ∈ S, since c ∈ T S. We then deduce from Remark 2.11 in [START_REF] Cherfils | On a generalized Cahn-Hilliard equation with biological applications[END_REF] that

((f (c), c)) ≥ α Ω n i=1 [c 4 i + c i 2 c2 i ]dx, α > 0. We note that, for v = (v i ) i ∈ S, v → ( v -v 2 -1 + v 2 ) 1 2 , v → ( v -v 2 + v 2 ) 1 2 , v → ( A 1 2 v 2 + v 2 ) 1 2 , and v → ( Av 2 + v 2 ) 1 2
are norms in L, H, V , and W , respectively, which are equivalent to the usual ones.

Here,

• -1 = A -1 2
• . Throughout this article, the same letter ξ (and, sometimes, ζ and η) denotes constants which may vary from line to line, or even in a same line. Similarly, the same letter Q denotes monotone increasing functions which may vary from line to line, or even in a same line.

A priori estimates

We first multiply (2.23) by c and have

1 2 d dt c 2 -1 = -((Pµ, c)) -λ 0 ((A -1 χ Ω\D (x)(c -h), c)) = -A 1 2 c 2 -((f (c), c)) -λ 0 ((χ Ω\D (x)(c -h), A -1 c)) ≤ -A 1 2 c 2 -α Ω (c 2 . c2 + c 2 . c2 )dx + λ 0 ( c + h ) c ≤ -A 1 2 c 2 -α Ω (c 2 . c2 + c 2 . c2 )dx + ζ( c 2 + c 2 ) + η( h 2 + 1) ≤ -A 1 2 c 2 -α Ω (c 2 . c2 + c 2 . c2 )dx + α 2 Ω (c 2 . c2 + c 2 . c2 )dx + η( h 2 + 1) ≤ -A 1 2 c 2 - α 2 Ω (c 2 . c2 + c 2 . c2 )dx + η( h 2 + 1),
owing to (2.24) and (2.26), hence

d dt c 2 -1 + A 1 2 c 2 + α Ω (c 2 . c2 + c 2 . c2 )dx ≤ ξ, (3.1)
which yields

d dt c 2 -1 + ξ c 2 -1 ≤ ζ, ξ > 0. (3.2)
It thus follows from (3.2) and Gronwall's lemma that

c 2 -1 ≤ e -ξt c0 2 -1 + ζ, ξ > 0, t ≥ 0. (3.3)
Let B be a bounded subset of Ḣ-1 (Ω) and t 0 be such that c0 ∈ B and t ≥ t 0 implies c(t) ∈ B 0 , where B 0 = {φ ∈ Ḣ-1 (Ω), φ 2 -1 ≤ 2ζ}, ζ being the constant in (3.3). We then deduce from (3.1) that, for t ≥ t 0 ,

t+r t A 1 2 c 2 ds ≤ ξ(r), t+r t ds Ω (c 2 . c2 + c 2 . c2 )dx ≤ ξ(r), (3.4)
for r > 0 fixed. We then multiply (2.23) by Ac and find 1 2

d dt c 2 = -((Pµ, Ac)) -λ 0 ((χ Ω\D (x)(c -h), c)) = -Ac 2 -((f (c), Ac)) -λ 0 ((χ Ω\D (x)(c -h), c)) ≤ -Ac 2 -((f (c)A 1 2 c, A 1 2 c)) + λ 0 ( c + h ) c ≤ -Ac 2 + ξ 0 A 1 2 c 2 + ζ Ω ( Ω c2 . c2 + c 2 . c2 )dx + η h 2 ≤ -Ac 2 + ξ 0 A 1 2 c 2 + ζ Ω (c 2 . c2 + c 2 . c2 )dx + h 2 + 1 ,
owing to (2.24) and (2.25), hence

d dt c 2 + Ac 2 ≤ ξ 0 A 1 2 c + ζ Ω (c 2 . c2 + c 2 . c2 )dx + h 2 + 1 . (3.5)
It thus follows from (3.4), (3.5), and the uniform Gronwall's lemma that

c 2 ≤ ξ, t ≥ t 0 + r, (3.6)
where the constant ξ is independent of c0 and t, hence, employing (3.6) and integrating (3.1) and (3.5) between 0 and t 0 + r,

c(t) 2 ≤ Q( c0 ), t ≥ 0. (3.7)
Now, writing c = c + c in (2.20), we find

d dt c + λ 0 Vol(Ω) Ω\D ( c + c -h)dx = 0.
Therefore,

d dt c + ξ c = - λ 0 Vol(Ω) Ω\D (c -h)
| c | ≤ Q( c0 ) + | c 0 | e -ξt + ζ, ∀t ≥ 0, (3.8)
where ξ and ζ are two constants which are nonnegative and independent of t and c 0 . We finally multiply (2.18) by A 2 c to obtain, noting that

((APµ, A 2 c)) = ((Pµ, A 3 c)) = ((µ, A 3 c)),
the differential equalities

1 2 d dt Ac 2 = -((µ, A 3 c)) -λ 0 ((χ Ω\D (x)(c -h), A 2 c)) = -A 2 c 2 -((Af (c), A 2 c)) -λ 0 ((χ Ω\D (x)(c -h), A 2 c)).
(3.9)

Using the fact that

∆f (c) = f (c)∆c + f (c)|∇c| 2 ,
we have (see, e.g., [START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF])

|((Af (c), A 2 c))| = |((-∆f (c), A 2 c))| ≤ ξ ∆f (c) 2 + 1 4 A 2 c 2 ≤ 1 2 A 2 c 2 + ξ.
(3.10) Furthermore,

λ 0 |((χ Ω\D (x)(c -h), A 2 c))| ≤ ξ c -h A 2 c ≤ ξ( c 2 + c 2 ) + 1 2 A 2 c 2 + ζ h 2 . (3.11)
We finally deduce from (3.9), (3.10), and (3.11) that

d dt Ac 2 ≤ ξ( c 2 + c 2 ) + ζ. (3.12)
Note that, integrating (3.5) over (t, t + r), we have, for t ≥ t 0 + r,

t+r t Ac 2 dx ≤ ξ(r), (3.13)
owing to (3.6). Furthermore, owing to (2.20), we have

d dt c 2 = 2 c d dt c ≤ 2λ 0 | c | Ω\D (c -h)dx ≤ ξ( c 2 + ( c ) 2 ) + ζ( h 2 + 1) (3
where r > 0 is fixed.

Well-posedness and existence of the global attractor

We have the Theorem 4.1. We assume that c 0 ∈ H. Then, (2.18)-(2.19) possesses a unique solution c such that

c ∈ L ∞ ([0, T ], H) ∩ L 2 ([0, T ], W ) ∩ L 4 ([0, T ], (L 4 (Ω)) n ∩ S), ∀T > 0.
Proof. We first denote by E m the space

E m = span{N 1 , N 2 , ..., N m }
and by P m the orthogonal projection from V onto E m ,

P m h = m j=1 ((h, N j ))N j .
The variational formulation of the problem reads: Find

c m : [0, T ] → c m + E m , c m (t) = c m (t) + c m (t) = c m (t) + m j=1 c m j (t)N j , such that d dt A -1 m i=1 cm (t)N i , N j - m i=1 cm i (t)A 1 2 N i , A 1 2 N j + ((Pf ( c m + c m ), N j )) + λ 0 ((A -1 (χ Ω\D (x)( c m + c m -h)), N j )) = 0, j = 1, ..., m, (4.1)
d dt c m + ξ c m = - λ 0 Vol(Ω) Ω\D (c m -h)dx, (4.2
)

cm (0) = P m c 0 , (4.3)
c m (0) = c 0 . (4.4)
We can rewrite (4.1) and (4.3) as

∂A -1 cm ∂t + Pf ( c m + c m ) + Ac m + λ 0 A -1 (χ Ω\D (x)( c m + c m -h)) = 0, in E m , (4.5
)

cm (0) = P m c 0 , (4.6)
and we finally rewrite equation (4.1) as follows:

M -1 dY dt + M Y + H(Y ) = 0, where M = ((AN i , N j )) i,j=1,...,m , Y =       c m 1 . . . c mm      
, and

H(Y ) =       ((Pf ( c m + c m ) + λ 0 A -1 (χ Ω\D (P( c m + c m -h))), N 1 )) . . . ((Pf ( c m + c m ) + λ 0 A -1 (χ Ω\D (P( c m + c m -h))), N m ))      
, noting that c m can be expressed as a function of c m by solving (4.2). The matrix M is invertible and positive definite and H(Y ) depends continuously on Y . Applying Cauchy's theorem, we find that there exists a time t m ∈ (0, T ) and a solution Y to

dY dt + M 2 Y + M H(Y ) = 0
f (c m ) f (c) weakly in L 4 3 (Ω T ) ∩ T S,
owing to the weak dominated convergence theorem. Finally, we deduce that A -1 ∂cm ∂t A -1 ∂c ∂t weakly in L 4 3 (Ω T ). Thus, passing to the limit in (4.5), we obtain

∂A -1 c ∂t + Pf (c) + Ac + λ 0 A -1 (χ Ω\D (x)(Pc -Ph)) = 0, in L 4 3 (Ω T ) ∩ T S. (4.10)
Finally, we easily pass to the limit in (4.2), at least in a weak sense, i.e., solving explicitly this ODE. Let now c 1 and c 2 be two solutions to (2.18)-(2.19) with initial data c 0,1 and c 0,2 , respectively. We set c = c 1 -c 2 , c 0 = c 0,1 -c 0,2 , and µ = µ 1 -µ 2 and have

∂c ∂t + APµ + λ 0 χ Ω\D (x)c = 0, (4.11
)

Pµ = P(f (c 1 ) -f (c 2 )) + Ac, (4.12
)

c| t=0 = c 0 = c 0,1 -c 0,2 . (4.13)
Integrating (4.11) over Ω, we obtain

d dt c = - λ 0 Vol(Ω) Ω\D cdx. (4.14)
We can thus rewrite the problem as

∂c ∂t + APµ + λ 0 χ Ω\D (x)c = 0, (4.15
)

Pµ = P(f (c 1 ) -f (c 2 )) + Ac, (4.16)
or else

∂A -1 c ∂t + Pµ + λ 0 A -1 χ Ω\D (x)c = 0, (4.17
)

Pµ = P(f (c 1 ) -f (c 2 )) + Ac. (4.18)
We multiply (4.17) by c and find

1 2 d dt c 2 -1 = -((µ, c)) -λ 0 ((A -1 χ Ω\D (x)c, c)) = -((f (c 1 ) -f (c 2 ), c)) + ((f (c 1 ) -f (c 2 ), c )) -A 1 2 c 2 -λ 0 ((χ Ω\D (x)c, A -1 c)) ≤ξ 0 c 2 + ((f (c 1 ) -f (c 2 ), c )) -A 1 2 c 2 -λ 0 ((χ Ω\D (x)c, A -1 c)) ≤ξ 0 c 2 + ((f (c 1 ) -f (c 2 ), c )) -A 1 2 c 2 + ξ c c -1 ≤ξ( c 2 + c 2 ) + ζ c 2 -1 -A 1 2 c 2 + ((f (c 1 ) -f (c 2 ), c )) ,
owing to (2.25) and (4.18). We note that, owing also to (2.25),

((f (c 1 ) -f (c 2 , c )) = c Ω c 1 0 f (c 1 + s(c 2 -c 1 ))dsdx ≤ ξ| c | Ω (|c 1 | 2 + |c 2 | 2 + 1)|c|dx ≤ ξ| c |( c 1 2 (L 4 (Ω)) n + c 2 2 (L 4 (Ω)) n + 1) c ≤ ξ( c 1 2 (L 4 (Ω)) n ∩S + c 2 2 (L 4 (Ω)) n ∩S + 1)(| c | 2 + c 2 ), hence 1 2 d dt c 2 -1 + A 1 2 c 2 ≤ ξ( c 1 2 (L 4 (Ω)) n ∩S + c 2 2 (L 4 (Ω)) n ∩S +1)(| c | 2 + c 2 ) + ζ c 2 -1 . (4.19)
Noting then that 1 2

d dt | c | 2 ≤ | c | d dt c ≤ ξ| c || Ω\D cdx| ≤ ξ| c | c ≤ ξ( c 2 + c 2 ) (4.20)
and recalling the interpolation inequality

c 2 ≤ ξ c -1 A 1 2 c , it follows that d dt ( c 2 -1 + | c | 2 ) + A 1 2 c 2 ≤ ξ( c 1 4 (L 4 (Ω)) n ∩S + c 2 4 (L 4 (Ω)) n ∩S + 1)(| c | 2 + c 2 -1 ). (4.21)
We thus deduce from (4.21) and Gronwall's lemma that

c 1 (t) -c 2 (t) L ≤ Q(T, c 0,1 , c 0,2 ) c 0,1 -c 0,2 L , 0 ≤ t ≤ T, (4.22)
hence the uniqueness, as well as the continuous dependence with respect to the initial data in the L-norm.

It follows from Theorem 4.1 that we have the continuous (with respect to the Lnorm) semigroup

S(t) : H → H, c 0 → c(t), t ≥ 0 (i.e., S(0) = I, S(t + s) = S(t) • S(s), t, s ≥ 0).
It then follows from (3.16) that S(t) possesses a bounded absorbing set B 1 which is compact in H and bounded in W . We thus deduce from standard results (see, e.g., [START_REF] Miranville | Attractors for dissipative partial differential equations in bounded and unbounded domains[END_REF][START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF]) the following theorem.

Theorem 4.2. The semigroup S(t) possesses the connected global attractor A such that A is compact in H and bounded in W .

Existence of exponential attractors

Let c 1 and c 2 be two solutions to (2.18)-(2.19) with initial data c 0,1 and c 0,2 , respectively. We set c = c 1 -c 2 , c 0 = c 0,1 -c 0,2 , and µ = µ 1 -µ 2 and have

∂c ∂t + APµ + λ 0 χ Ω\D (x)c = 0, (5.1)
Pµ = P(f (c 1 ) -f (c 2 )) + Ac, (5.2)
c| t=0 = c 0 . (5.3)
Furthermore, it is sufficient here to take initial data belonging to the bounded absorbing set B 1 defined in the previous section. We multiply (5.1) by tc to find, owing to (2.25) and (5.2),

1 2 d dt (t c 2 ) + λ 0 t χ Ω\D (x)c 2 = 1 2 c 2 -t Ac 2 -t((A 1 2 (f (c 1 ) -f (c 2 )), A 1 2 c)) ≤ 1 2 c 2 + ξ 0 t A 1 2 c 2 , which yields d dt (t c 2 ) ≤ c 2 + ξt A 1 2 c 2 .
(5.4)

Integrating (5.4) between 0 and t, we obtain

c 2 ≤ ξ 1 + t t t 0 c 2 V ds, t > 0. (5.5)
Noting that (4.21) is equivalent to

d dt ( c 2 -1 + | c | 2 ) + c 2 V ≤ ξ( c 1 4 (L 4 (Ω)) n ∩S + c 2 4 (L 4 (Ω)) n ∩S + 1)(| c | 2 + c 2 -1 ), (5.6)
it follows from (3.1), (4.22), and (5.6) that

t 0 c 2 V ds ≤ ξe ζt c 2 L . (5.7)
We then deduce from (5.4) and (5.7) that

t c 2 H ≤ ξe ζt c 2 L , t > 0. (5.8)
Next, we derive a Hölder (both with respect to space and time) estimate. Actually, owing to (4.22), it suffices to prove the Hölder continuity with respect to time. We have

c(t 1 ) -c(t 2 ) L = t 2 t 1 ∂c ∂t dτ L ≤ t 2 t 1 ∂c ∂t L dτ ≤ |t 1 -t 2 | 1 2 t 2 t 1 ∂c ∂t 2 L dτ 1 2
.

(5.9)

Multiplying (2.23) by ∂c ∂t , we find

∂c ∂t 2 -1 = -f (c), ∂c ∂t - 1 2 d dt A 1 2 c 2 + λ 0 χ Ω\D (x)(h -c), A -1 ∂c ∂t ≤ A 1 2 f (c) ∂c ∂t -1 - 1 2 d dt A 1 2 c 2 + λ 0 c -h ∂c ∂t -1 ,
owing to (2.24), hence

d dt A 1 2 c 2 + ∂c ∂t 2 -1 ≤ ξ c 2 V + ζ, (5.10)
owing to (2.25) and the continuous embedding H 2 (Ω) ⊂ C(Ω). Noting that

∂c ∂t 2 ≤ ξ c 2 + ζ, (5.11)
owing to (2.20), we then deduce from (3.16), (5.10), and (5.11) that

t 2 t 1 ∂c ∂t 2 L dτ ≤ ξ, (5.12)
where the constant ξ depends on B 1 and T such that t 1 , t 2 ∈ [0, T ], so that

c(t 1 ) -c(t 2 ) L ≤ ξ|t 1 -t 2 | 1 2 , (5.13)
where the constant ξ depends on B 1 and T such that t 1 , t 2 ∈ [0, T ]. We finally deduce from (4.22), (5.8), and (5.13) the following result (see, e.g., [START_REF] Eden | Expenential Attractors for Dissipative Evolution Equations[END_REF][START_REF] Efendiev | Exponential attractors for a nonlinear reaction-diffusion system in R 3[END_REF]).

Theorem 5.1. The semigroup S(t) possesses an exponential attractor M ⊂ B 1 , i.e., (i) M is compact in L;

(ii) M is positively invariant, S(t)M ⊂ M, ∀t ≥ 0;

(iii) M has finite fractal dimension in L;

(iv) M attracts exponentially fast the bounded subsets of H,

∀B ⊂ H bounded, dist L (S(t)B, M) ≤ Q( B H )e -ξt , ξ > 0, t ≥ 0,
where the constant ξ is independent of B and dist L denotes the Hausdorff semidistance between sets defined by

dist L (A, B) = sup a∈A inf b∈B a -b L .
Remark 5.2. Setting M = S(1)M, we can prove that M is an exponential attractor for S(t), but now in the topology of H (see, e.g., [START_REF] Efendiev | Exponential attractors for a singularly perturbed Cahn-Hilliard system[END_REF]).

Since M (or M) is a compact attracting set, we deduce from Theorem 5.1 the following corollary. Remark 5.4. We can more generally consider a nonlinear term F of the form

F (v) = n i=1 2p+2 k=0 a k v k i , a 2p+2 > 0, v = (v i ) i ∈ S,
but, now, assumption (2.26) is replaced by

((f (c), c)) ≥ α Ω n i=1 p k=0 c2p+2-2k i c 2k dx, α > 0, c = (c i ) i ∈ S.
Indeed, since

F (c) = n i=1 2p+2 k=0 a k c k i , a 2p+2 > 0, we have ∂F (c) ∂c i = 2p+1 k=0 (k + 1)a k+1 c k i , i = 1, ..., n, hence ∂ 2 F (c) ∂c 2 i ≥ -ξ i , ξ i ≥ 0, i = 1, .
.., n, and, owing to Remark 2.1,

((f (c)v, v)) = Ω n i=1 ∂ 2 F (c) ∂c 2 i v 2 i dx ≥ - n i=1 ξ i v i 2 , ξ i ≥ 0,
for all v ∈ T S. Furthermore, it follows from Remark 2.11 in [START_REF] Cherfils | On a generalized Cahn-Hilliard equation with biological applications[END_REF] that

(( ∂F (c) ∂c i , ci )) ≥ α Ω p k=0 c2p+2-2k i c 2k dx, α > 0, c ∈ S, i = 1, .
.., n, and, owing to Remark 2.1,

((f (c), c)) = Ω n i=1 ∂F (c) ∂c i ci dx ≥ α Ω n i=1 p k=0 c2p+2-2k i c 2k dx, α > 0,
for all c ∈ S.

Remark 5.5. We studied the (global) dynamics of the system when λ 0 is constant and ε = 1. More generally, all constants here and in the previous sections depend on λ 0 and ε and grow as these quantities go to +∞ and 0, respectively. In particular, the dependence on λ 0 is an important issue in view of inpainting applications (see [START_REF] Bertozzi | Analysis of a two-scale Cahn-Hilliard model for binary image inpainting[END_REF]) and a natural question is whether one can have an upper bound on the dimension of the global attractor which is independent of this quantity (here, the upper bound that we obtain explodes as λ 0 goes to +∞). A natural (but involved, see, e.g., [START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF]) problem would be to find a lower bound on this dimension (possibly, in terms of λ 0 ). This will be investigated elsewhere.

Algebraic consistency with the two-phase model

In order to ensure the "physical" constraints of our proposed model, the model must (at least, see below) satisfy two properties:

(i) When phases (i 1 , i 2 , i 3 , ..., i n-2 ) ∈ {1, ..., n} are not present, the n-phase free energy is equal to the one of the two-phase model.

(ii) When phases (i 1 , i 2 , i 3 , ..., i n-2 ) ∈ {1, ..., n} are not present at the initial time (i.e., c i j (0) = 0 and h i j = 0, j = 1,...,n -2), these phases do not appear artificially during the evolution of the system.

In that case, we say that the model is algebraically consistent with the two-phase model (see [START_REF] Boyer | Study of a three component Cahn-Hilliard flow model, M2AN[END_REF]). Of course, the model should also be consistent with the k-phase model, k < n, but we have not been able to prove this more general property here.

Theorem 6.1. Model (2.8)-(2.9) is algebraically consistent with the diphasic system.

Proof. We first recall that the total free energy of (2.8)-(2.9) (n-phase model) can be written as follows:

F n ((c i ) i ) = Ω 1 n n i=1 (c 2 i (1 -c i ) 2 ) + ε 2 2 n i=1 |∇c i | 2 dx + λ 0 2 n i=1 Ω\D (c i -h i ) 2 dx. (6.1)
We assume that c i = c, c j = 1 -c, and c k = h k = 0, for i and j fixed in {1, 2, ..., n} and k = 1, ..., n, k = i, j. We set h i = h (i.e., h j = 1 -h) and c = (0, ..., 0, c, 0, .., 0, 1c, 0, .., 0) ∈ R n . Then, the total free energy can be rewritten as

F n (c) = Ω 2 n (c 2 (1 -c) 2 ) + ε 2 |∇c| 2 dx + λ 0 Ω\D (c -h) 2 dx (6.2)
which is equal, up to a multiplicative constant, which, therefore, does not change anything in the description of the problem, to the total free energy of the two-phase model (1.1),

F(c) = Ω 1 n (c 2 (1 -c) 2 ) + ε 2 2 |∇c| 2 dx + λ 0 2 Ω\D (c -h) 2 dx, (6.3)
so that the model satisfies (i). Furthermore, (ii) says that our particular c has to be solution to (2.8)-(2.9), which is equivalent to

∂F ∂c k - 1 n n l=1 ∂F ∂c l | c k =0,c j =1-c i = 0. (6.4)
Moreover, we have

∂F ∂c k - 1 n n l=1 ∂F ∂c l | c k =0,c i =c,c j =1-c i = - 2 n 2 c(1 -c) 2 -c 2 (1 -c) + (1 -c)c 2 -(1 -c) 2 c = 0, (6.5)
which yields that the model satisfies (6.4) and then (ii).

Numerical simulations

As far as the numerical simulations are concerned, we rewrite the problem in the form ∂c i ∂t + ∆µ i + λ 0 χ Ω\D (x)(c i -h i ) = 0, i = 1, ..., n, (7.1)

µ i = 2 ∆c i -f i (c), i = 1, ..., n, (7.2)
∂c i ∂ν = ∂µ i ∂ν = 0 on Γ, i = 1, ..., n, (7.3)

c |t=0 = c 0 , (7.4)
where

f i (c) = 2 n [c i (1 -c i ) 2 -c 2 i (1 -c i )] -1 n n i=1 ∂F (c)
∂c i , i = 1, ..., n. Problem (7.1)-(7.4) has the advantage of splitting the fourth-order (in space) equation into a system of two second-order ones (see [START_REF] Elliott | A second order splitting method for the Cahn-Hilliard equation[END_REF][START_REF] Grasselli | A splitting method for the Cahn-Hilliard equation with inertial term[END_REF][START_REF] Injrou | Stable discretizations of the Cahn-Hilliard-Gurtin equations[END_REF]). Consequently, we use a P1-finite element for the space discretization, together with a semi-implicit Euler time discretization (i.e., implicit for the linear terms and explicit for the nonlinear ones). The numerical simulations are performed with the software Freefem++ (see [33]).

In the numerical results presented below, Ω is a (0, 0.5) × (0, 0.5)-square. The triangulation is obtained by dividing Ω into 100 × 100 rectangles and by dividing each rectangle along the same diagonal.

In order to obtain the final inpainting results, we use a dynamic one step algorithm with threshold involving the diffuse interface thickness ε which allows us to connect regions across large inpainting domains (see [START_REF] Cherfils | Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation in image inpainting[END_REF]). Here, the thresholding consists in replacing max i=1,...,n c i by 1 and c j = max i=1,...,n c i by 0 for all j and at every point (x, y) ∈ Ω (more generally, at every point (x, y, z) ∈ Ω); this means that we replace the dominant phase (color) by 1 at every point of Ω and the other phases (colors) by 0 to obtain the final inpainting result. The green region in Figure 1(a) corresponds to the inpainting region. We run the modified Cahn-Hilliard system (n = 3) with ε = 0.03, ∆t = 0.001, and λ 0 = 500000. We are close to a steady state at t = 0.006 and we replace the dominant color by 1 and the other colors by 0 to obtain the final inpainting in Figure 1(b). In Figure 2(a), the light blue region corresponds to the inpainting region. We run the modified Cahn-Hilliard system (n = 9) with ε = 0.03, ∆t = 0.05, and λ 0 = 900000.

Nine colors inpainting

  x) is a given binary image, and D Ω is the inpainting domain. Here, c satisfies the Neumann boundary conditions ∂c ∂ν = ∂∆c ∂ν = 0, on ∂Ω, and F is the antiderivative of f such that

| 0 ee-ξt t 0 e ξs c ds = ζe -ξt t 0 +r 0 e

 000 dx, where ξ = λ 0 Vol(Ω\D) Vol(Ω) , hence d dt (e ξt c ) = -λ 0 Vol(Ω) e ξt Ω\D (c -h)dx and c = e -ξt c 0c | ≤ e -ξt | c 0 | + ζe -ξt t ξs ( c + h )ds, ∀t ≥ 0, ξs h ds ≤ ζe -ξt h e ξt ≤ ζ h ≤ η. Furthermore, for t ≥ t 0 + r, ζe ξs c ds + ζe -ξt t t 0 +r e ξs c ds ≤ Q( c0 )e -ξt e ξ(t 0 +r) + ζe -ξt t t 0 +r e ξs c ds ≤ Q( c0 )e -ξt + ζe -ξt t t 0 +r e ξs c ds ≤ Q( c0 )e -ξt + ηe -ξt (e ξt -e ξ(t 0 +r) ) ≤ Q( c0 )e -ξt + ζ, where we have used (3.6) (resp., (3.7)) to estimate the first (resp., the second) integral in the right-hand side of the first equality. Finally, we obtain

4 3 4 3 4 3

 444 on the time interval [0, t m [. Having c m , we then deduce c m (from (4.2)). It follows from the a priori estimates derived in the previous section for the solution c m (t) (c(t) being replaced by c m (t)) that any local solution to (2.18)-(2.19) is actually a global solution defined on the whole interval [0, T ]. It then follows from the a priori estimates that, up to a subsequence which we do not relabel, c m c weakly in L 2 (0, T, W ), m → ∞. It follows from (4.7), (4.8), and the Aubin-Lions compactness theorem thatc m → c strongly in L (0, T, (L (Ω)) n ∩ S)(4.9) and c m (t, x) → c(x, t) a.e. (t, x) ∈ [0, T ] × Ω. Moreover, since f is a continuous, f (c m (t, x)) → f (c(t, x)) a.e. and, since f (c m ) is bounded in L (Ω T ), Ω T = (0, T ) × Ω,

Corollary 5 . 3 .

 53 The semigroup S(t) possesses the finite-dimensional global attractor A ⊂ B 1 .

Figure 1 :

 1 Figure 1: (a) Inpainting region in green, random initial datum between 0 and 1 in inpainting region, ε = 0.03. (b) Solution at t = 0.006; replacing the dominant color by 1 and the other colors by 0.

Figure 2 :

 2 Figure 2: (a) Inpainting region in light blue, random initial datum between 0 and 1 in inpainting region, ε = 0.03. (b) Solution at t = 0.2; replacing the dominant color by 1 and the other colors by 0.

Figure 4 :

 4 Figure 4: (a) Inpainting region in light blue, random initial datum between 0 and 1 in inpainting region, ε = 0.03. (b) Solution at t = 0.2; replacing the dominant color by 1 and the other colors by 0.

We are close to a steady state at t = 0.2 and we replace the dominant color by 1 and the other colors by 0 to obtain the final inpainting in Figure 2 The light green region (horizontal bars) in Figure 3(a) corresponds to the inpainting region. We run the modified Cahn-Hilliard system (n = 10) with ε = 0.008, λ 0 = 3000000, and ∆t = 0.01. We are close to a steady state at t = 0.04 and we replace the dominant color by 1 and the other colors by 0 to obtain the final inpainting in Figure 3(b).

Consistency with the two-phase model example

In Figure 4, we take the same example of nine colors inpainting as in Figure 2, but, now, we take h i = c i (0) = 0, i = 1, ..., 7. We run again the modified Cahn-Hilliard system (n = 9) with the same parameters as in Figure 2. We are close to a steady state at t = 0.2 and we replace the dominant color by 1 and the other colors by 0 to obtain the final inpainting in Figure 4(b). Note that |c i (t)| < 10 -5 , ∀t ≥ 0, i = 1, ..., 7.


A-L Basdevant
email: anne-laure.basdevant@u-paris10.fr

N Enriquez
email: nenriquez@u-paris10.fr

L Gerin
email: gerin@cmap.polytechnique.fr

J-B Gouéré

Discrete Hammersley's lines with sources and sinks

Keywords: MSC 2010 Classification: 60K35, 60F15 Hammersley's process, Ulam's problem, longest increasing subsequences

Introduction

In a celebrated paper, J.M.Hammersley used Poissonization to attack the so-called Ulam problem of the typical length (n) of the longest increasing subsequence of a uniform permutation of size n. Namely, he reduced this problem to finding the greatest number of points of a Poisson point process inside a square, an increasing path can go through. He proved ( [START_REF] Hammersley | A few seedlings of research[END_REF], Theorem 4) that (n)/ √ n converges in probability to some constant c, sometimes refered to as the Ulam constant, and conjectured that c = 2.

The proof of c = 2 was achieved independently by Logan and Shepp and by Vershik and Kerov in 1977, using algebraic methods. Various authors were then interested in finding a more probabilistic proof of this result. First, Aldous and Diaconis [START_REF] Aldous | Hammersley's interacting particle process and longest increasing subsequences[END_REF] gave one, using the properties of what they called Hammersley's process, which was implicitly introduced in [START_REF] Hammersley | A few seedlings of research[END_REF] (p.358 and following). Hammersley's process is continuous in time and space and Aldous and Diaconis studied its properties on the infinite line, in particular its stationary distributions. A few years later, Groeneboom [START_REF] Groeneboom | Hydrodynamical methods for analyzing longest increasing subsequences[END_REF] and Cator and Groeneboom [START_REF] Cator | Hammersley's process with sources and sinks[END_REF] studied Hammersley's process on a quarter plane. By adding what they called Poisson sinks and sources on the x and y-axis, they also found a stationary version of this process on the quarter plane. Using this point of view, they were able to recover again the value of c.

In this paper, we study two discrete variants of Ulam's problem. Namely, for all p in [0, 1], n, m ≥ 1 , we replace the original Poisson point process by the following random set ξ of integer points: each integer point of the rectangle [1, n] × [1, m] is chosen independently with probability p. We are interested in the two following quantities.

• The length of the longest increasing subsequence through points of ξ: L [START_REF] Aldous | Hammersley's interacting particle process and longest increasing subsequences[END_REF] (n,m) = max {L; there exists (i 1 , j 1 ), . . . , (i L , j L ) ∈ ξ, i 1 < • • • < i L and j 1 < j 2 < • • • < j L } .

• The length of the longest non-decreasing subsequence through points of ξ:

L (2)
(n,m) = max {L; there exists (i 1 , j 1 ), . . . , (i L , j L ) ∈ ξ, i

1 < • • • < i L and j 1 ≤ j 2 ≤ • • • ≤ j L } .
One aim of the present work is to recover by simple and unified probabilistic arguments the first order asymptotics of L [START_REF] Aldous | Hammersley's interacting particle process and longest increasing subsequences[END_REF] (n,m) and L [START_REF] Cator | Hammersley's process with sources and sinks[END_REF] (n,m) already obtained by T. Seppäläinen in two independent papers, for L [START_REF] Aldous | Hammersley's interacting particle process and longest increasing subsequences[END_REF] (n,m) in [START_REF] Seppäläinen | Increasing sequences of independent points on the planar lattice[END_REF] and for L [START_REF] Cator | Hammersley's process with sources and sinks[END_REF] (n,m) in [START_REF] Seppäläinen | Exact limiting shape for a simplified model of first-passage percolation on the plane[END_REF]. In the following result the variables

L (i) (n,m) (n,m)
are coupled in the obvious way. Moreover, for any a, b > 0, L

(an,bm) stands for L (i) ( an , bn ) .

Theorem 1. For a, b > 0, we have, when n tends to infinity

L (1) (an,bn) n a.s. →      √ p(2 √ ab -(a + b) √ p) 1 -p if p < min{a/b, b/a}, min{a, b} otherwise, (1) L (2) (an,bn) n a.s.
To prove this, Seppäläinen associates to each problem a particle system on the infinite line Z and checks that Bernoulli product measures are stationary. The position of particles at a given time is then characterized as the solution of a discrete optimization problem, and the a.s. limit of L (i) (an,bn) /n is identified using convex analysis arguments. Note also that for the second question, Johansson used later ( [START_REF] Johansson | Discrete orthogonal polynomial ensembles and the Plancherel measure[END_REF], Th.5.3) another description of L

(2) (an,bn) /n and proved that the rescaled fluctuations converge to the Tracy-Widom distribution.

In [START_REF] Rolla | The discrete and continuum broken line process. Markov Process[END_REF], Rolla-Sidoravicius-Surgailis-Vares considered the closely related model of last-passage percolation with geometric weights. Here is a quick description of the model. Let ( ξ(i, j)) 1≤i≤n,1≤j≤m be a family of independent geometric random variables with parameter p. Set:

L (3) (n,m) = max t ξ(γ(t))
where the maximum is taken over all paths γ from (1, 1) to (n, m) with (1, 0) or (0, 1) steps. With respect to

L (1) (n,m) or L (2)
(n,m) , we have simply changed the weight of each point (from Bernoulli to geometric) and the set of paths considered. Rolla-Sidoravicius-Surgailis-Vares gave a new proof of the following result:

L (3) (n,n) n a.s. -----→ n→+∞ 2 √ p 1 - √ p .
(The original proof was given by Jockusch-Propp-Shor in [START_REF] Jockusch | Random domino tilings and the arctic circle theorem[END_REF].) Our proof of Theorem 1 is inspired by the proof of Rolla-Sidoravicius-Surgailis-Vares of the above result. The strategy is to consider a particle system on a bounded domain which turns out to coincide with the restriction of Seppäläinen's particle system on Z. This simplifies the definition of the process (especially, in the case of the second problem). Moreover, it turns out that a local balance property around a single site, see Lemma 1 and 2 below, is enough to check the stationarity of the process. Our proofs are essentially the same for both models. This kind of remarkable local balance property also occurs in last-passage percolation with geometric weights (see Rolla-Sidoravicius-Surgailis-Vares , Section 3.1 in [START_REF] Rolla | The discrete and continuum broken line process. Markov Process[END_REF] or Seppäläinen, Lemma 2.3 in [START_REF] Seppäläinen | Lecture notes on the corner growth model[END_REF]). The particle systems studied here provide two other examples where such a property holds. Theorem 1 is then proven by investigating the behaviour of the models under the different stationary measures which have been exhibited.

A nice by-product of our proof is that we obtain non-asymptotic estimates which provide an elementary proof of c = 2 for the original Ulam problem (see the discussion at the end of Section 4).

In the last section, which can be read independently, we complete the study of the original infinite discrete Hammersley process introduced by Seppäläinen [START_REF] Seppäläinen | Increasing sequences of independent points on the planar lattice[END_REF] for the Problem 1. We prove that the only extremal translation-invariant stationary measures of this process are the Bernoulli product measures.

Discrete Hammersley's processes

Like Hammersley did, we construct two sets of broken lines whose cardinality is respectively the variable L

(n,m) and L

(2) (n,m) . For each case, we first introduce a partial order: Case 1:

Case 2: (x, y) ≺ (x , y ) iff (x, y) ≺ (x , y ) iff x < x and y < y

x < x and y ≤ y Now, Hammersley's lines are paths starting from the top side of the rectangle [1, n] × [1, m], ending at its right side and making only south or east steps. They are constructed recursively. The first line is the highest non-increasing path connecting the minimal points of ξ for ≺. We withdraw these points from ξ and connect the new minima to get the second line, and so on. In the below picture, n = m = 8 and L We claim that the number

H (i) (n,m) of Hammersley's lines is equal to L (i)
(n,m) . Indeed, note that, in both cases, L (i) (n,m) is equal to the length of the longest increasing subsequence trough points of ξ for the partial order ≺. Each Hammersley's line connecting minimal points of a subset of ξ, the L (i) (n,m) points of a maximizing sequence necessarily lie on distinct Hammersley's lines and therefore

H (i) (n,m) ≥ L (i)
(n,m) . To prove the opposite inequality, we observe that Hammersley's lines give a construction of an increasing sequence of points a 1 , . . . a

H (i) (n,m)
of ξ in the following way.

For a

H (i) (n,m)
we take any point of ξ which lies in the top right Hammersley's line. Then, for each , a -1 is taken as one of the points of ξ which lies in the -1-th line and such that a -1 ≺ a . The existence of such a -1 is granted. Indeed, if no such point existed, a would have belonged to the ( -1)-th line. This proves

L (i) (n,m) ≥ H (i) (n,m) .
Let us note that the construction of lines above is consistent in n and m, i.e. for n < n and m < m the restriction of Hammersley's lines to the rectangle [1, n ] × [1, m ] coincides with Hammersley's lines constructed from the points of ξ in this smaller rectangle. In particular, for n ≥ 1, we can define the two discrete time processes (X 1 t ) t≥0 := (X 1 t (x), x ∈ 1, n ) t≥0 and (X 2 t ) t≥0 := (X 2 t (x), x ∈ 1, n ) t≥0 defined by

X i t (x) = 1 if
there is a Hammersley's line on the vertical edge {(x, t), (x, t + 1)} 0 otherwise.

Note that on each vertical edge there is at most one Hammersley's line. Moreover, in Model 1, there is also at most one Hammersley's line on any horizontal edge and Hammersley's lines do not intersect each others. We will say that there is a particle at time t at position x in Model i if X i t (x) = 1. Hence, for both models we start with the initial condition X i 0 ≡ 0 and one can check that both processes are Markovian. Besides, since Hammersley's lines start all from the top side of the rectangle, for both models, the number of particles at time t is equal to

L (i) (n,t) .
In both models, the dynamic of the particle system is quite simple. As the description of the dynamic is not needed in the proof of Thereom 1, we do not give such a description here. However, we refer the interested reader to Appendix A.

Sinks and sources

Since, in both models, the number of particles can only increase as time goes on, both processes converge to their unique stationary measure X i ∞ :≡ 1. In this section, we are going to modify a little bit the problems we are looking at, such that the particle systems associated to the new problems admit less trivial stationary measures. To do so, we must first introduce some notation.

From now on, for any integer x ≥ 1 and t ≥ 1, we define the random variable ξ t (x) by

ξ t (x) = 1 if (x, t) ∈ ξ 0 otherwise. ( 4
)
With this notation, we can write the quantities

L (i) (n,m) as L (i) (n,m) = max s ξ t(s) (x(s)) (5)
where the maximum is taken over all the discrete paths (x(•),

t(•)) in [1, n] × [1, m] starting from (1, 1
) which make only steps in N × N for Model 1 and in N × Z + for Model 2. Given two sequences of integers (ξ t (0)) t∈ 1,m and (ξ 0 (x)) x∈ 1,n , we now consider the quantities

L (i) (n,m) = max s ξ t(s) (x(s)) (6)
where the maximum is taken over all the discrete paths (x(•), t(•)) in [0, n] × [0, m] starting from (0, 0) such that:

• There exists an integer u ∈ Z + such that: either they first make u (0, 1) steps or they first make u (1, 0) steps.

• Then, they only make steps in N × N for Model 1 and in N × Z + for Model 2.

Using the same terminology as in [START_REF] Cator | Hammersley's process with sources and sinks[END_REF], we say that there is k sinks (resp. sources) at site (0, t)

(resp. (x, 0)) if ξ t (0) = k (resp. ξ 0 (x) = k).
The sources will be drawn on the x-axis and the sinks on the y-axis. Hence, in other words, L

(n,m) is equal to the maximum number of points an increasing (resp. non-decreasing) path can go trough when the path can pick either sources or either sinks before going into N × N. The law of the sources and sinks will be specified later, let us just say that the sources will take values in {0, 1} n whereas the sinks will take values in {0, 1} m for Model 1 and in Z m + for Model 2.

(n,n) = 5 Model 2, L (1)
(n,n) = 7

As before, we can construct a set of broken lines whose cardinality is equal to

L (i)
(n,m) . Now, Hammersley's lines are paths starting from the top side or the left side of the rectangle [0, n] × [0, m], ending at its bottom side or its right side and making only south or east steps. They are still constructed recursively. Assume that the sink of minimal height is located at (0, t 1 ) and the leftmost source at (x 1 , 0). Then, the first line is the highest non-increasing path starting from (0, t 1 ), ending at (x 1 , 0) and connecting minimal points (x, t) of ξ for the partial order ≺ introduced in the previous section such that (x 1 , 0) ⊀ (x, t) and (0, t 1 ) ⊀ (x, t). If there is no sink (resp. source) at all, we do the same construction except that the line starts from the top of the box (resp. ends on the right of the box). Then, we withdraw these points from ξ, one of the sink located at (0, t 1 ) and the source located at (x 1 , 0) and we do the same procedure with this new set of sources, sinks, and set of points in 1, n × 1, m to obtain the second line. The algorithm goes on until there is no more points, sinks and sources. Using the same argument as in the previous section, one can easily prove that the number of lines obtained with this procedure is equal to

L (i) (n,m) . With a slight abuse of notation, we still denote (X 1 t ) t≥0 := (X 1 t (x), x ∈ 1, n ) t≥0 and
n = 5 Model 2, L (1)
n = 7 (X 2 t ) t≥0 := (X 2 t (x), x ∈ 1, n ) t≥0 (2)
the two discrete time system of particles defined by

X i t (x) = 1 if
there is a Hammersley line (with sources and sinks) on the edge {(x, t), (x, t + 1)} 0 otherwise.

Hence, the processes now start from the configuration given by the sources i.e X i 0 (x) = ξ 0 (x).

Theorem 2 (Stationarity of Hammersley's processes on a bounded interval, with sources and sinks). For all n,

• Model 1. For all p, α ∈ (0, 1), the process (X 1 t (x), x ∈ 1, n ) t≥0 is stationary if sources are i.i.d. Ber(α) and sinks are i.i.d. Ber(α ) (and sinks independent from sources) with

α = p(1 -α) α + p(1 -α) . (7)
• Model 2. For all p, α ∈ (0, 1) such that α > p, the process (X 2 t (x), x ∈ 1, n ) t≥0 is stationary if sources are i.i.d. Ber(α) and sinks are i.i.d. (Geo(α ) -1) (i.e. P(ξ(0, t) = k) = α (1α ) k for k = 0, 1, . . . ) (and sinks independent from sources) with

α = α -p α(1 -p) . ( 8
)
(Note that α ∈ (0, 1) if α > p.)
Remark 1. In [START_REF] Seppäläinen | Increasing sequences of independent points on the planar lattice[END_REF] (resp. in [START_REF] Seppäläinen | Exact limiting shape for a simplified model of first-passage percolation on the plane[END_REF]), Seppäläinen associates to Model 1 (resp. Model 2) a particle system defined on the whole line Z. One can check that the restrictions of these systems to 1, n have respectively the same transition probabilities that our processes X i t , i = 1, 2 when sinks and sources are distributed as in Theorem 2. Hence,

• For Model 1 (resp. Model 2), Theorem 2 should be compared to Lemma 2.1 in [START_REF] Seppäläinen | Increasing sequences of independent points on the planar lattice[END_REF] (resp.

Proposition 1 in [START_REF] Seppäläinen | Exact limiting shape for a simplified model of first-passage percolation on the plane[END_REF]) which states that Bernoulli product measures are stationary for the analogous model on the infinite line.

• The condition α > p in Model 2 reminds the condition on the density of the initial configuration needed to define the process on Z in [START_REF] Seppäläinen | Exact limiting shape for a simplified model of first-passage percolation on the plane[END_REF].

• The theorem states stationarity in time ( i.e. from bottom to top in our figures). In fact, the proof also shows stationarity from left to right.

Let us collect a few consequences for further use. Let L

L (i),0,0 (n,m) = T (i),0,0 (n,m) d = L (i) (n,m) , L (i),α,β (n,m) = T (i),α,β
(n,m) + # {sinks between 1 and m} .

(

) 9
This relies on the fact that any lines has to exit by the top or by a sink. Note that in the right-hand side of ( 9) the two terms are not independent. Moreover, using the interpretation of

L (i) (n,m) and L (i),α,β
(n,m) as the length of a longest subsequence (cf. ( 5) and ( 6)), we have

L (i) (n,m) ≤ L (i),α,β (n,m) . (10)
Besides, T

m) being by definition equal to the number of particles at time m, Theorem 2 implies that

T (i),α,α (n,m) law = Binom(n, α). (11)
Proof. We aim to prove Theorem 2. We first prove for both models a local balance property, in Lemmas 1 and 2. These lemmas are elementary but do not seem to be written elsewhere, and form the heart of our proof. Model 1.

We first focus on what happens around a single point (x, t). Recall that ξ t (x) is a Bernoulli(p). Denote by X (resp. Y, X , Y ) be the indicator that a Hammersley's line hits (x, t) from the bottom (resp. left,top,right).

Lemma 1 (Local balance for Model 1). The random variables (X , Y ) are measurable with respect to (X, Y, ξ t (x)) and we have

(X, Y, ξ t (x)) (d) = Ber(α) ⊗ Ber(α ) ⊗ Ber(p) ⇒ (X , Y ) (d) = Ber(α) ⊗ Ber(α )
where α is as in [START_REF] Liggett | Interacting particle systems[END_REF].

X ∼Ber(α)

Y ∼Ber(α )

X ∼Ber(α) Y ∼Ber(α ) ind.
ind.

ξ t (x)

Proof of Lemma 1. Recalling that, in Model 1, Hammersley's lines are non increasing lines which do not touch each other and noticing that X = Y = 1 iff X = Y = 0 and ξ t (x) = 1, we get that

(X , Y ) =    (X, Y ) if X = Y (1, 1) if X = Y = 0 and ξ t (x) = 1. (0, 0) otherwise.
Hence, we see that (X , Y ) ∈ σ(X, Y, ξ t (x)) and one can easily check that if (X , Y ) ∼ Ber(α) ⊗ Ber(α ) with α is as in [START_REF] Liggett | Interacting particle systems[END_REF], then (X , Y ) ∼ Ber(α) ⊗ Ber(α ).

We now explain how Lemma 1 proves Theorem 2 for Model 1. For every (x, 1) define the corresponding (X (x), Y (x)) be the output of (x, 1). Independence of sources and sinks and Lemma 1 ensure that the output of (1, 1) is independent of sources at (2, 0), (3, 0), . . . and sinks at (0, 2), (0, 3), . . . . Then a simple induction on x proves that the random variables (X (x), x ∈ 1, n ) are distributed as Ber(α). Model 2. The strategy is similar:

Lemma 2 (Local balance for Model 2). The random variables (X , Y ) are measurable with respect to (X, Y, ξ t (x)) and we have

(X, Y, ξ t (x)) (d) = Ber(α) ⊗ (Geo(α ) -1) ⊗ Ber(p) ⇒ (X , Y ) (d) = Ber(α) ⊗ (Geo(α ) -1)
where α is as in [START_REF] Rolla | The discrete and continuum broken line process. Markov Process[END_REF].

X ∼Ber(α)

Y ∼Geo(α ) -1 X ∼Ber(α) Y ∼Geo(α ) -1 ind.
ind.

ξ t (x)

Proof of Lemma 2. Recall that X, X ∈ {0, 1} whereas Y, Y ∈ Z + . Besides, using that Hammersley's lines are non-increasing, the balance between incoming and outcoming lines at a site yields

X + Y = X + Y . Moreover, X = 1 iff (ξ t (x) = 1 or (X, Y ) = (1, 0)).
We observe that the different cases can be summed up in

(X , Y ) = (1, 0) if (X, Y ) = (1, 0) (ξ t (x), Y -X + ξ t (x)) otherwise.
One can easily check that this equality implies Lemma 2. For instance, for k ≥ 1,

P(X = 1, Y = k) = αα (1 -α ) k should be equal to P(X = 1, Y = k) = P(ξ t (x) = 1, X = 0, Y = k -1) + P(ξ t (x) = 1, X = 1, Y = k) = p(1 -α)α (1 -α ) k-1 + pαα (1 -α ) k = p(1 -α ) k-1 α [(1 -α) + α(1 -α )] ,
which requires (8).

Theorem 2 for Model 2 then follows from Lemma 2 in the same way as before Theorem 2 for Model 1 follows from Lemma 1.

Law of Large Numbers in Hammersley's processes

In this section, we explain how Theorem 2 on the stationary measure of the processes with sources and sinks implies Theorem 1.

Proof of Theorem 1. Trivial case.

We first consider Model 1 in the case p ≥ min{a/b, b/a}. The upper bound is straightforward. The lower bound can be proven as follows. Let us consider for example the case p ≥ ab -1 . We have to prove:

lim inf n→∞ L (1) (an,bn) n ≥ a. ( 12
)
One can build an increasing path as follows. The first point (1, y 1 ) of the path is the lowest point of ξ having first coordinate equal to 1. The second point (2, y 2 ) of the path is the lowest point of ξ having first coordinate equal to 2 and second coordinate strictly larger than y 1 . The other points are defined in a similar fashion. As pb ≥ a, a study of this path easily provides (12). The case p ≥ ba -1 can be handled in a similar way. The same strategy treats the case p ≥ a/(a + b) in Model 2.

Non trivial case

We only need to compute lim

1 n E[L (i)
(an,bn) ] since almost sure convergence follows from superadditivity. Model 1. Upper bound. We assume that p ≥ min{a/b, b/a}. For any α ∈ (0, 1), taking α as in [START_REF] Liggett | Interacting particle systems[END_REF] and using ( 9), ( 10) and ( 11), we get

1 n E[L (1)
(an,bn) ] ≤

1 n E[L (
= aα + bα = aα + b p(1 -α) α + p(1 -α) =: φ (1) (a, b, α). ( 13
)
The latter is minimized for

α(a, b) := √ p b a -p 1 -p , α (a, b) := √ p a b -p 1 -p , (14)
which both are in (0, 1) if p < min{a/b, b/a}. This yields

1 n E[L (1)
(an,bn) ] ≤ φ (1) (a, b, α(a, b)) = √ p(2 √ ab -(a + b) √ p) 1 -p .
Model 2. Upper bound. We assume that p ≥ a/(a + b). Using [START_REF] Seppäläinen | Lecture notes on the corner growth model[END_REF] with α = α-p α(1-p) and that E[Geo(α

) -1] = 1 α -1, we get 1 n E[L (2)
(an,bn) ] ≤

1 n E[T (i),α,α (an,bn) ] + 1 n E[# {sinks between 1 and bn}] = aα + b( 1 α -1) = aα + b p(1 -α) α -p .
We minimize the latter by taking

α = p + b a p(1 -p).
Note that this choice is allowed since for p ∈ (0, a/(a + b)) then 0 < p < α < 1.

Model 1 and 2. Lower bound.

In [START_REF] Seppäläinen | Increasing sequences of independent points on the planar lattice[END_REF][START_REF] Seppäläinen | Exact limiting shape for a simplified model of first-passage percolation on the plane[END_REF] the lower bound was obtained with a convexity argument on the scaling limit. Instead, we adapt a more probabilist argument due to Rolla-Sidoravicius-Surgailis-Vares (proof of Theorem 4.1 in [START_REF] Rolla | The discrete and continuum broken line process. Markov Process[END_REF]) for the closely related model of last passage percolation with geometric weights. Here, we only give the details for Model 1 but the same argument applies for Model 2 (there are only minor modifications due mainly to the fact that sinks have no more a Bernoulli distribution but a geometric distribution).

Let us consider Model 1 on the rectangle [0, an] × [0, bn] with sinks and sources with optimal source intensity α = α(a, b) and α = α (a, b) (see ( 14)). For ε ∈ [0, 1], denote by L

(1) (an,bn) (ε) the length of the largest non-decreasing subsequence defined by

L (1) (an,bn) (ε) = max {L; (i 1 , 0), . . . , (i k , 0), (i k+1 , j k+1 ), . . . , (i L , j L ) ∈ ξ, 0 < i 1 < • • • < i k ≤ anε < i k+1 < . . . < i L ≤ an and 0 < j k+1 < • • • < j L ≤ bn} .
Here is an example where L

(an,bn) (ε) = 4. Note that because of the constraint L

(an,bn) (ε) is smaller than the number of Hammersley's lines L (1),α,α (an,bn) .

i k i k+1 anε

Recall that φ (1) has been defined in (13).

Lemma 3 (Largest subsequence using some sources). There exist positive and non-decreasing functions f , g on (0, 1], which depend on a, b, such that

P L (1) (an,bn) (ε) ≥ n φ (1) (a, b, α(a, b)) -f (ε) ≤ exp(-ng(ε)). Proof of Lemma 3. We write L (1) (an,bn) (ε) = I 1 + I 2 with I 1 := max {k; (i 1 , 0), . . . , (i k , 0) ∈ ξ, i 1 < • • • < i k ≤ anε} I 2 := max {k; (i 1 , j 1 ), . . . , (i k , j k ) ∈ ξ, anε < i 1 < • • • < i k ≤ an and 0 < j 1 < • • • < j k ≤ bn} .
The random variables I 1 and I 2 are independent, I 1 is binomially distributed with parameters ( anε , α(a, b)). The random variable I 2 has the law of L . Using (13), we get:

1 n E L (1)
(an,bn) (ε) ≤ This is due to the fact that α(a, b) = α(a , b ) and that there is a unique α ∈ [0, 1] which minimizes φ(a , b , α). Thus,

1 n E [I 1 + I 3 ] ≤ aεα(a, b) + φ (1) a(1 -ε), b, α(a(1 -ε), b) ,
aεα(a, b) + φ (1) a(1 -ε), b, α(a(1 -ε), b) < aεα(a, b) + φ (1) a(1 -ε), b, α(a, b) = aεα(a, b) + a(1 -ε)α(a, b) + bα (a, b) = aα(a, b) + bα (a, b) = φ (1) (a, b, α(a, b)).
Denote by d(ε) the difference between the right and left hand side of the inequality. The function d is positive on (0, 1] and continuous. We set f (ε) = 1 2 min{d(δ), δ ∈ [ε, 1]} which is non-decreasing and positive. Then, we have

P L (1) (an,bn) (ε) ≥ n φ (1) (a, b, α(a, b)) -f (ε) ≤ P I 1 + I 3 ≥ E(I 1 + I 3 ) + n d(ε) -f (ε) ≤ P (I 1 + I 3 ≥ E(I 1 + I 3 ) + nf (ε)) .
Recalling that I 1 is binomially distributed and I 3 is the sum of two binomial random variables, function g is obtained by applying Hoeffding's inequality to I 1 and I 3 .

We still consider Model 1 on the rectangle [0, an]×[0, bn] with sinks and sources with optimal source intensity α = α(a, b) and α = α (a, b). Let π n be an optimal path for L (1),α,α (an,bn) . If there are several optimal paths, we choose one of them in an arbitrary way. Define D n as the number of sources (i, j) ∈ π n with j = 0. Lemma 4 (Optimal paths do not take many sources). There exists a positive function h on (0, 1] such that for and any δ ∈ (0, 1] and for any n large enough

P (D n > anδ) ≤ exp(-nh(δ)).
In particular D n /n converges to 0 in L 1 .

Proof of Lemma 4. By definition, if D n > anδ then for some ε ≥ δ such that εan ∈ N and ≤ 1

• there are more than anδ sources in {1, . . . , εan}; ,α,α (an,bn) + # {sinks between 1 and bn}.

• L (1) (an,bn) (ε) = L (1),α,α (an,bn) = T (1)
Therefore

P {D n > anδ} ≤ P T (1),α,α (an,bn) ≤ n(aα - 1 2 f (δ)) + P # {sinks} ≤ n(bα - 1 2 f (δ)) + P L (1)
(an,bn) (ε) > n (aα + bαf (δ)) for some ε as above

≤ exp(-ng(δ)) + 1≥ε≥δ; εan∈N P L (1) (an,bn) (ε) > n φ (1) (a, b, α(a, b)) -f (δ) ≤ exp(-ng(δ)) + 1≥ε≥δ; εan∈N P L (1) (an,bn) (ε) > n φ (1) (a, b, α(a, b)) -f (ε) ≤ exp(-ng(δ)) + 1≥ε≥δ; εan∈N exp(-ng(ε))
≤ exp(-ng(δ)) + an exp(-ng(δ))

for large n, where we used Lemma 3 and where g is some positive function. This implies the first part of the lemma. The convergence of D n /n to 0 in probability and in L 1 follows, since the sequence is bounded.

We now conclude the proof of the lower bound noticing that we have

L (1),α,α (an,nb) -D n -D n ≤ L (1) (an,bn)
where D n is the number of sinks (i, j) ∈ π n with i = 0. By the previous lemma we know that D n /n tends to 0 in L 1 . By symmetry, the same result holds for D n /n. Taking expectations of both sides in the previous inequality we get lim inf

n→+∞ 1 n E[L (1)
(an,bn) ] ≥ aα + bα = φ (1) (a, b, α(a, b)).

Back to Ulam's constant

We observe that if we take p = 1/n, Theorem 1 for Model 1 suggests L

n ≈ 2 √ n, which is consistent with the asymptotics of Ulam's problem, since ξ is then close, after renormalization, to a Poisson point process with intensity n.

In fact, one can rigorously recover that c = 2 using our proof of Theorem 1. To do so, consider a Poisson point process Ξ with intensity n in the unit square and denote by (n) the greatest number of points of Ξ an increasing path can go through. To get a lower bound and an upper bound of (n), we divide the square [0, 1] 2 in two different ways.

First, we fix some k ≥ 1 and we divide it into small squares of length side 1/(k √ n). Say that ξ j (i) = 1 if at least one point of Ξ is in the square with top-right corner (i/(k n), j/(k √ n)) and consider the quantity

L (1) k √ n := L (1) (k √ n,k √ n) associated to the family (ξ j (i)) i,j≤k √ n . It is clear that (n) ≥ L (1) k √ n . (15)
Denoting

p k = P(ξ j (i) = 1) = P Poiss(1/k 2 ) ≥ 1 = 1 -e -1/k 2 , Theorem 1 implies L (1) k √ n √ n a.s. -----→ n→+∞ 2k √ p k √ p k + 1 .
Using (15) and letting k tend to infinity, we get c := lim (n)/ √ n ≥ 2. To prove the upper bound, we divide now the square into small squares of length side 1/n 4 . Say that ξ j (i) = 1 if at least one point of Ξ is in the square with top-right corner (i/n 4 , j/n 4 ) and consider the quantity L (1)

n 4 := L (1)
(n 4 ,n 4 ) associated to the family (ξ j (i)) i,j≤n 4 . The parameter of these Bernoulli random variables is now

pn = P Poiss(1/n 7 ) ≥ 1 = 1 -e -1/n 7 .
With probability higher than 1n -2 , all the columns and lines of width 1/n 4 contain at most one point of Ξ. On this event that we denote F n , L

n 4 coincides with (n). We now use some intermediate results of the proof of the upper bound for a = b = 1. Inequality [START_REF] Seppäläinen | Exact limiting shape for a simplified model of first-passage percolation on the plane[END_REF] still holds despite the dependence of pn on n i.e.

L

(1)

n 4 ≤ L (1),αn,α n (n 4 ,n 4 ) , with α n = α n = √ pn-pn 1-pn ∼ n -7/2 . Using (13),we get E[ (n)] ≤ E[|Ξ|1 Fn ] + E[L (1),αn,α n (n 4 ,n 4 ) ] ≤ E[Poiss(n) 2 ]/n 2 + n 4 (α n + α n ),
where we used the Cauchy-Schwarz inequality. Dividing by n and letting n go to infinity, we obtain c ≤ 2.

Stationary measures of Hammersley's process on Z

In this section, which can be read independently from the rest of the article, we study the analogous on the infinite line Z of the process (X 1 t ) t≥0 defined in (3). This infinite Hammersley process was introduced by T.Seppäläinen to prove the Law of Large Numbers for Problem 1.

In order to generalize (X t ) t≥0 to a process taking its values in {0, 1} Z we need a construction that does not rely on Hammersley lines. As before, (ξ t (i)) i∈Z,t∈N denote i.i.d. Bernoulli(p) random variables we say that there is a cross at time t located at x if ξ t (x) = 1. Informally, the infinite Hammersley process is defined as follows:

At time t, if there is a particle at x and if the particle immediately on its left is at y, then the particle at x jumps at time t + 1 at the rightmost cross of ξ t+1 in interval (y, x) (if any, otherwise it stays still).

For our purpose, we need a construction of (X t ) in which we let the crosses act one by one on configurations. We first draw X 0 ∈ {0, 1} Z at random according to some distribution µ such that µ-almost surely, for all i there is j < i such that X 0 (j) = 1.

We now explain how to construct (X t+1 ) from (X t ). A cross will act on a given configuration of particles in the following way: if a cross is located at x, then the leftmost particle of (X t ) in the interval x, +∞ (if any) moves to x; if there is no such particle then a particle is created at x.

With this definition, we can now construct the value of (X t+1 ) as a function of (X t ) and the crosses at time t + 1: We define X t+1 as the result of the successive actions on X t of all crosses at time t + 1 from the right to the left. Note that by definition, there only a finite number of crosses in ξ t+1 that can modify X t (i): those between X t (i) and the first particle on its left.

An example is drawn in this figure where the circles represent the particles and the crosses the locations where ξ t+1 equals 1. Here, we have X 1 t+1 = (1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, ...) (as in the previous sections, time goes from bottom to top in our pictures): This theorem is precisely the discrete counterpart of ( [START_REF] Aldous | Hammersley's interacting particle process and longest increasing subsequences[END_REF], Lemma 7). The fact that Bernoulli product measures are stationary was proved by T.Seppäläinen ([9] Lemma 2.1). The reason for which we only focus on Model 1 on the infinite line is that the analogous of Model 2 is much more complicated to analyze. Indeed, the evolution of a particle in Model 2 depends on the whole configuration on its left. We don't know if a similar statement to Theorem 3 holds for Model 2.

In order to prove that there are no other extremal measures we adapt the proof of Jockush-Propp-Shor ( [START_REF] Jockusch | Random domino tilings and the arctic circle theorem[END_REF] p.17-20) for the discrete-time TASEP (see also [START_REF] Liggett | Interacting particle systems[END_REF] Ch.VIII). The main difference here with Jockush-Propp-Shor is our Lemma 5 below that replaces their block argument.

Proof of extremality. Let α ∈ (0, 1] and µ be an extremal translation-invariant stationary measure on {0, 1} Z with marginals µ(X(i) = 1) = α, we want to prove that µ = Ber(α) ⊗Z .

The proof is made of the following steps:

1. we introduce a measure π on {0, 1} Z × {0, 1} Z which is a minimal coupling between Ber(α) ⊗Z and µ;

2. we prove (Lemma 6) that some patterns that would decrease card {i; x i = y i }, have πprobability zero;

3. we conclude: µ = Ber(α) ⊗Z .

Let M be the set of translation-invariant measures on {0, 1} Z × {0, 1} Z that are a coupling of Ber(α) ⊗Z and µ. The set M is non-empty (it contains Ber(α) ⊗Z ⊗ µ) and compact for the weak topology.

We will prove that in M there is a coupling ((

X(i)) i∈Z , (Y (i)) i∈Z ) such that X = Y a.s. Set D : M → [0, 1] π → π (X(0) = Y (0)) .
The function D is continuous on the compact set M and thus attains its minimum δ ≥ 0 at some coupling π. If δ = 0, then the theorem is proved: by translation-invariance X = Y π-a.s. Let (X 0 , Y 0 ) ∼ π, we use the same random variables (ξ t (i)) i∈Z,t∈N to define the dynamics of X and Y . We denote by (X t , Y t ) the pair of configurations at time t, its joint distribution is denoted by π t . By construction, for every t, π t is translation-invariant. For every k ≥ 1 we denote by ∆ k (X t , Y t ) the number of i ∈ {1, . . . , k} such that X t (i) = Y t (i).

Let n ≥ 2 and t ≥ 0. We say that there is a n-forbidden pattern at location x ∈ Z at time t if the configuration (X t , Y t ) between locations x and x + n -1 is either

  x x + n -1 X t = ( 1 0 . . . 0 0 ) Y t = ( 0 0 . . . 0 1 )   or   x x + n -1 X t = ( 0 0 . . . 0 1 ) Y t = ( 1 0 . . . 0 0 )   . Lemma 5. For n ≥ 2, let E x,n be the event ξ 1 (x) = 1, ξ 1 (x+1) = ξ 1 (x+2) = • • • = ξ 1 (x+n-1) = 0. Let j ≥ 1,
denote by F (n, j) the subset of locations of {0, n, 2n, . . . , (j -1)n} at which is located a forbidden pattern at time 0. Let A(n, j) = card{x ∈ F (n, j), E x,n occurs}.

Then:

∆ jn (X 1 , Y 1 ) ≤ ∆ jn (X 0 , Y 0 ) + 1 -2A(n, j).
Proof of Lemma 5. We let the crosses at time 1 act one by one from the right to the left. We consider the impact of each cross action on the discrepancy ∆ jn between the two configurations. Let us note that, when a cross located at i acts on some configuration of particles, it changes the value of at most two sites of the configuration: i and the leftmost 1 of the configuration in the interval i, ∞). Thus, we deduce the following facts:

• The crosses strictly to the right of jn have no impact on the discrepancy.

• A cross located at a point x in {1, . . . , jn} cannot increase the discrepancy. To check this fact, one can study all the cases. They are all shown below, up to symmetries between X and Y .

  x ( 1 ) ( 1 )   ,   x ( 1 ? • • • ? 1 ) ( 0 0 • • • 0 1 )   ,   x ( 1 ? • • • ? 0 ) ( 0 0 • • • 0 1 )   ,   x ( 0 • • • 0 1 ) ( 0 • • • 0 1 )   ,   x ( 0 • • • 0 0 0 0 0 1 ) ( 0 • • • 0 1 ? • • • ? 1 )   ,   x ( 0 • • • 0 0 0 0 0 1 ) ( 0 • • • 0 1 ? • • • ? 0 )   .
• For a forbidden pattern at some in, if E in,n occurs, the cross at in decreases the discrepancy by 2.

• In the worst case, the right-most cross which is strictly to the left of 1 increases the discrepancy by 1.

• The crosses to the left of the previous one have no impact on the discrepancy.

The result follows.

Lemma 6 (Forbidden patterns). Let n ≥ 2. The π 0 probability of a n-forbidden pattern at any location is 0. Now we can obtain our contradiction. Indeed, there exists an integer n such that π

  1 n X = ( 1 ? . . . ? 0 ) Y = ( 0 ? . . . ? 1 )   > 0 or π   1 n X = ( 0 ? . . . ? 1 ) Y = ( 1 ? . . . ? 0 )   > 0.
Without loss of generality we do the first case. We can assume that X, Y coincide between 2 and n -1, otherwise we can decrease n. Let 1 < k < n be the leftmost position at which

X(k) Y (k) = 1
1 . With positive probability we can turn this 1 1 into 1

0 :     × ------ 1 k n X t = ( 1 0 . . . 0 1 ? 0 ) Y t = ( 0 0 . . . 0 1 ? 1 )     t+1 →     1 k n X t+1 = ( ? 0 . . . 0 1 ? 0 ) Y t+1 = ( ? 0 . . . 0 0 ? 1 )    
We repeat this process between positions k and n and finally for some t, i, j we have

π   i j X t = ( 1 0 . . . 0 0 ) Y t = ( 0 0 . . . 0 1 )   > 0,
which contradicts Lemma 6.

A Dynamic of the underlying interacting particle system

In this appendix we explicit the dynamic of the underlying interacting particle system in Models 1 and 2. Consider the processes (X i t ) t≥0 defined by (3) (without sources and sinks). One can explicitely construct the dynamic of these processes in term of the family of i.i.d. random variables {ξ t (x), x ∈ 1, n , t ∈ N} with law Bernoulli(p) (see ( 4)). To do so, let say that there is a cross at time t located at x if ξ t (x) = 1. A cross will act on a given configuration of particles in the following way : if a cross is located at x, then the left-most particle in the interval x, n (if any) moves to x ; if there is no such particle then a particle is created at x.

With this definition, we can now construct the value of X i t+1 as a function of X i t and the crosses at time t + 1: Model 1. We define X 1 t+1 as the result of the successive actions on X 1 t of all crosses at time t + 1 from the right to the left.

An example is drawn in this figure where the circles represent the particles and the crosses the locations where ξ t+1 equals 1. Here, we have X 1 t+1 = (1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0) (as it is usual in the literature on Hammersley's processes, time goes from bottom to top in our pictures):

1 n

X 1 t+1 X 1 t
Let us note that there is an alternative way to construct the model 1 : At time t, if there is a particle at x and if the particle immediately on its left is at y, then the particle at x jumps at time t + 1 at the rightmost cross of ξ t+1 in interval (y, x) (if any, otherwise it stays still). If there are some points of ξ t+1 at the right of the rightmost particle of X 1 t then a new particle appears at the leftmost point among them. Model 2. We define X 2 t+1 as the result of the successive actions on X 2 t of all crosses at time t + 1 from the left to the right. We let the reader convince themselves that this construction coincides with the definition given in [START_REF] Groeneboom | Hydrodynamical methods for analyzing longest increasing subsequences[END_REF].

Although the definition of both models seems, at first glance, very close, the nature of the two processes is in fact quite different. Indeed, in the first model, to find the location at time t + 1 of a particle located at time t at x, one just need to know the location y of the particle immediately on his left at time t and the position of the crosses of ξ t+1 in the interval (y, x). In particular, with no difficulty, one can define a process with similar transitions on the whole line Z as it is done in [START_REF] Seppäläinen | Increasing sequences of independent points on the planar lattice[END_REF]. For Model 2, the dependances are more intricate. Indeed, to determine the location at time t + 1 of a particle located at time t at x, one need to know the whole configuration of X 2 t and ξ t+1 on the interval 1, x . In particular, the definition on the whole line Z of a similar process is more delicate and requires a condition between the density of crosses and particles (see [START_REF] Seppäläinen | Exact limiting shape for a simplified model of first-passage percolation on the plane[END_REF]).

Finally, let us note that the construction can be extended to the models with sources and sinks. The set of sources simply give the initial configuration of the particle system. The sinks acts as the crosses. In Model 1 sinks act after the other crosses while in Model 2 sinks act before the other crosses.

Aknowledgements. The authors are glad to aknowledge David Aldous who was at the origin of their interest in this problem.

→ 2 abp( 1 -

 1 p) + (ab)p if p < a/(a + b), aotherwise.

  m) = 5 (crosses denote points of ξ, Hammersley's lines are in blue, one of the longest subsequences is shown in a red dashed line):

  (i),α,β (n,m) be the number of Hammersley's lines in Model i with parameters α (sources), β (sinks) in the box [0, n] × [0, m] and let T (i),α,β (n,m) be the number of Hammersley's lines that leave the same box from the top. Then, for i = 1, 2 and every α, β ∈ [0, 1]

  {sinks between 1 and bn}]

  ε),bn) , thus is dominated by a random variableI 3 with law L (1),α(a(1-ε),b),α (a(1-ε),b) (an(1-ε),bn)

  where, by convention, α(a , b ) = 1 when p ≥ a /b and α(a , b ) = 0 when p ≥ b /a . One can check that for p < min{a/b, b/a} and a , b > 0 such that a/b = a /b we have φ (1) (a , b , α(a, b)) > φ (1) (a , b , α(a , b )).

Theorem 3 .

 3 The only extremal translation-invariant stationary measures of Hammersley's process on Z are measures Ber(α) ⊗Z for all α ∈ (0, 1].

Proof of Lemma 6. Let j ≥ 1. By Lemma 5 we have: ∆ jn (X 1 , Y 1 ) ≤ ∆ jn (X 0 , Y 0 ) + 1 -2A(n, j).

Note that the probability that, at some location in, there is a n-forbidden pattern and that E in,n occurs is ζp(1p) n-1 where ζ is the probability of a n-forbidden pattern at a given location. We aim to prove ζ = 0. Taking expectation in the previous display we get

Dividing by nj and letting j → ∞ we get:

But the minimality of π yields D(π 1 ) ≥ D(π 0 ) and we get ζ = 0.

Lemma 7 (There is a forbidden pattern somewhere). Assume that δ > 0, there exists n such that

The two terms on the right side are equal:

Then 0 < δ/2 = π(X(0) = 1, Y (0) = 0). Assume that the lemma is false (for every n), we have π(A ∪ B) = 1, where

As seen above, δ > 0 implies then that π(A), π(B\A) > 0. Besides, it is easy to check that A and B\A are preserved by the dynamics. Hence, considering the restriction of π on these two subsets, π A := (ν A , µ A ) and π B\A := (ν B\A , µ B\A ), we see that ν A (resp. µ A ) and ν B\A (resp. µ B\A ) define two translation-invariant stationary measures such that Ber(α) ⊗Z

The extremality of Ber(α) ⊗Z and µ implies that Ber(α) ⊗Z = ν A = ν B\A and µ = µ A = µ B\A . But, by definition of A and B, ν A µ A and µ B\A ν B\A . Thus, necessarily, µ = Ber(α) ⊗Z .