Athmane Bakhta
email: athmane.bakhta@cermics.enpc.fr

Virginie Ehrlacher
email: virginie.ehrlacher@enpc.fr

Cross-diffusion systems with non-zero flux and moving boundary conditions

, we prove the existence of a global weak solution to the obtained system. In addition, existence of a solution to an optimization problem defined on the fluxes is established under the assumption that the solution to the considered cross-diffusion system is unique. Lastly, we prove that in the case when the imposed external fluxes are constant and positive and the entropy density is defined as a classical logarithmic entropy, the concentrations of the different species converge in the long-time limit to constant profiles at a rate inversely proportional to time. These theoretical results are illustrated by numerical tests.

Introduction

The aim of this work is to propose and analyze a mathematical model for the description of a Physical Vapor Deposition (PVD) process, the different steps of which are described in details for instance in [START_REF] Donald | Handbook of physical vapor deposition (PVD) processing[END_REF]. Such a technique is used in several contexts, for instance for the fabrication of thin film crystalline solar cells. The procedure works as follows: a wafer is introduced in a hot chamber where several chemical elements are injected under a gaseous form. As the latter deposit on the substrate, an heterogeneous solid layer grows upon it. Two main phenomena have to be taken into account: the first is naturally the evolution of the surface of the film; the second is the diffusion of the various species in the bulk, due to the high temperature conditions. Experimentalists are interested in controlling the external gas fluxes that are injected into the chamber, so that, at the end of the process, the spatial distributions of the concentrations of the diverse components inside the new layer are as close as possible to target profiles.

In this article, a one-dimensional model which takes into account these two factors is studied. We see this work as a preliminary step before tackling more challenging models in higher dimensions, including surfacic diffusion effects for instance. This will be the object of future work. Our main motivation for the study of such a model concerns the optimization of the external fluxes injected in the chamber during a PVD process.

More precisely, let us assume that at a time t ≥ 0, the solid layer is composed of n + 1 different chemical species and occupies a domain of the form (0, e(t)) ⊂ R + , where e(t) > 0 denotes the thickness of the film. The evolution of e(t) is determined by the fluxes of atoms that are absorbed at the surface of the layer. At time t > 0 and point x ∈ (0, e(t)), the local volumic fractions of the different species are denoted respectively by u 0 (t, x), • • • , u n (t, x). Let us point out that if the molar volume of the solid is uniform in the thin film layer and constant during all the process, then u i (t, x) is also equal (up to a constant multiplicative constant) to the local concentration of the i th species at time t > 0 and point 0 ≤ x ≤ e(t). Up to some renormalization condition, it is natural to expect that these functions are non-negative and satisfy a volumic constraint which reads as follows: ∀0 ≤ i ≤ n, u i (t, x) ≥ 0 and n i=0 u i (t, x) = 1.

(1)

Because of the constraint [START_REF] Nicholas D Alikakos | Lp bounds of solutions of reaction-diffusion equations[END_REF], it holds that u 0 (t, x) = 1 -n i=1 u i (t, x) for all t > 0 and x ∈ (0, e(t)). Thus, the knowledge of the n functions u 1 , • • • , u n is enough to determine the dynamics of the whole system. Replacing u 0 by 1 -n i=1 u i , and denoting by u the vector-valued function (u 1 , • • • , u n ), the evolution of the concentrations inside the bulk of the solid layer is modeled through a system of cross-diffusion equations of the form

∂ t u -∂ x (A(u)∂ x u) = 0, for t > 0, x ∈ (0, e(t)), (2)
with approriate boundary and initial conditions, where A : [0, 1] n → R n×n is a matrixvalued function encoding the cross-diffusion properties of the different species.

Such systems have received much attention from the mathematical community in the case when no-flux boundary conditions are imposed on a fixed domain [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF][START_REF] Amann | Dynamic theory of quasilinear parabolic equations. ii. reactiondiffusion systems[END_REF][START_REF] Le | Everywhere regularity of solutions to a class of strongly coupled degenerate parabolic systems[END_REF][START_REF] André | Local existence, uniqueness and smooth dependence for nonsmooth quasilinear parabolic problems[END_REF]. Then, in arbitrary dimension d ∈ N * , the system reads

∂ t u -div x (A(u)∇ x u) = 0, for t > 0, x ∈ Ω,
for some fixed bounded regular domain Ω ⊂ R d and boundary conditions

(A(u)∇ x u) • n = 0 on ∂Ω and t ≥ 0,
where n denotes the outward normal unit vector to Ω. Such systems appear naturally in the study of population's dynamics in biology, and in chemistry, for the study of the evolution of chemical species concentrations in a given environment [START_REF] Painter | Continuous models for cell migration in tissues and applications to cell sorting via differential chemotaxis[END_REF][START_REF] Hillen | A user's guide to pde models for chemotaxis[END_REF]. The analysis of these systems is a challenging task from a mathematical point of view [START_REF] Lepoutre | Global well-posedness of a conservative relaxed cross diffusion system[END_REF][START_REF] Nicholas D Alikakos | Lp bounds of solutions of reaction-diffusion equations[END_REF][START_REF] Horst | Invariant regions for quasilinear reaction-diffusion systems and applications to a two population model[END_REF][START_REF] Redlinger | Invariant sets for strongly coupled reaction-diffusion systems under general boundary conditions[END_REF][START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF][START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF][START_REF] Di | Fully parabolic keller-segel model for chemotaxis with prevention of overcrowding[END_REF]. Indeed, the obtained system of parabolic partial differential equations may be degenerate and the diffusion matrix A is in general not symmetric and/or not positive definite. Besides, in general, no maximum principle can be proved for such systems. Nice counterexamples are given in [START_REF] Stará | Some (new) counterexamples of parabolic systems[END_REF]: there exist Hölder continuous solutions to certain cross-diffusion systems which are not bounded, and there exist bounded weak solutions which develop singularities in finite time.

It appears that some of these cross-diffusion systems have a formal gradient flow structure. Recently, an elegant idea, which consists in introducing an entropy density that appears to be a Lyapunov functional for these systems, has been introduced by Burger et al. in [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF]. This analysis strategy, which was later extended by Jüngel in [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] and named boundedness by entropy technique, enables to obtain the existence of global in time weak solutions satisfying [START_REF] Nicholas D Alikakos | Lp bounds of solutions of reaction-diffusion equations[END_REF] under suitable assumptions on the diffusion matrix A. It was successfully applied in several contexts (see for instance [START_REF] Jungel | Existence analysis of maxwell-stefan systems for multicomponent mixtures[END_REF][START_REF] Juengel | Entropy structure of a cross-diffusion tumorgrowth model[END_REF][START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF][START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF]).

However, there are very few works which focus on the analysis of such cross-diffusion systems with non zero-flux boundary conditions and moving domains. To our knowledge, only systems containing at most two different species have been studied, so that n = 1 and the evolution of the concentrations inside the domain are decoupled and follow independent linear heat equations [START_REF] Jacobus | Well-posedness of a parabolic movingboundary problem in the setting of wasserstein gradient flows[END_REF].

The one-dimensional model [START_REF] Amann | Dynamic theory of quasilinear parabolic equations. ii. reactiondiffusion systems[END_REF] we propose and analyze in this paper describes the evolution of the concentration of n+1 different atomic species, with external flux boundary conditions, in the case when the diffusion matrix A satisfies similar assumptions to those needed in the no-flux boundary conditions case studied in [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF].

The article is organized as follows: the results of [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] in the case of no-flux boundary conditions in arbitrary dimension are recalled in Section 1. We illustrate them on a prototypical example of diffusion matrix A, which is introduced in Section 1.1.

Our results in the case of a one-dimensional moving domain with non-zero flux boundary conditions are gathered in Section 2. We prove the existence of a global in time weak solution to (2) with appropriate boundary conditions and evolution law for e(t) in Section 2.2.1. The long time behaviour of a solution is analyzed in the case of constant external absorbed fluxes in 2.2.2 and an optimization problem is studied in 2.2.3. The proofs of these results are gathered in Section 3.

A numerical scheme used to approximate the solution of such systems is described in Section 4 and our theoretical results will be illustrated by several numerical tests. We refer the reader to [START_REF] Bakhta | Mathematical models and numerical simulation of photovoltaic devices[END_REF] for comparisons between our proposed model and experimental results obtained in the context of thin film solar cells fabrication.

Case of no-flux boundary conditions in arbitrary dimension

In Section 1.1, a particular cross-diffusion model on a fixed domain with no-flux boundary conditions is presented. The latter is a prototyical example of the systems of equations considered in this paper. Its formal gradient flow structure is highlighted in Section 1.1.2. Using slight extensions of results of [START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF][START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF], it can be seen that this system can be analyzed using the theoretical framework developped in [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF][START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF], which is recalled in Section 1.2.

Throughout this section, let us denote by d ∈ N * the space dimension, Ω ⊂ R d the regular bounded domain occupied by the solid. The local concentrations at time t > 0 and position x ∈ Ω of the n + 1 different atomic species entering in the composition of the material are respectively denoted by u 0 (t, x), • • • , u n (t, x). We also denote by n the normal unit vector pointing outwards the domain Ω.

Example of cross-diffusion system

Presentation of the model

As mentioned above, we have one particular example of system of cross-diffusion equations in mind, which is used to illustrate more general theoretical results. This system, with no-flux boundary conditions, reads as follows : for any

0 ≤ i ≤ n,            ∂ t u i -div x 0≤j =i≤n K ij (u j ∇ x u i -u i ∇ x u j ) = 0, for (t, x) ∈ R * + × Ω, 0≤j =i≤n K ij (u j ∇ x u i -u i ∇ x u j ) • n = 0, for (t, x) ∈ R * + × ∂Ω, (3)
where for all 0 ≤ i = j ≤ n, the positive real numbers K ij satisfy K ij = K ji > 0. They represent the cross-diffusion coefficients of atoms of type i with atoms of type j. This set of equations can be formally derived from a discrete stochastic lattice hopping model, which is detailed in the Appendix.

The initial condition

(u 0 0 , • • • , u 0 n ) ∈ L 1 (Ω; R n+1
) of this system is assumed to satisfy:

∀0 ≤ i ≤ n, u 0 i (x) ≥ 0, n i=0 u 0 i (x) = 1 and u i (0, x) = u 0 i (x) a.e. in Ω. (4)
The relationship n i=0 u 0 i (x) = 1 is a natural volumic constraint which encodes the fact that each site of the crystalline lattice of the solid has to be occupied (vacancies being treated as a particular type of atomic species).

Summing up the n + 1 equations of (3), we observe that a solution

(u 0 , • • • , u n ) must necessarily satisfy ∂ t ( n i=0 u i ) = 0.
It is thus expected that the following relationship should hold:

∀0 ≤ i ≤ n, u i (t, x) ≥ 0, n i=0 u i (t, x) = 1, a.e. in R * + × Ω. (5)
Plugging the expression u 0 (t, x) = 1 -n i=1 u i (t, x) in (3), it holds that for all

1 ≤ i ≤ n, 0= ∂ t u i -div x   1≤j =i≤n K ij (u j ∇ x u i -u i ∇ x u j )   -div x   K i0     1 - 1≤j =i≤n u j -u i   ∇ x u i -u i ∇ x   1 - 1≤j =i≤n u j -u i       = ∂ t u i -div x   1≤j =i≤n (K ij -K i0 ) (u j ∇ x u i -u i ∇ x u j ) + K i0 ∇ x u i   .
Thus, the system can be rewritten as a function of

u := (u 1 , • • • , u n ) T as follows    ∂ t u -div x (A(u)∇ x u) = 0, for (t, x) ∈ R * + × Ω, (A(u)∇ x u) • n = 0, for (t, x) ∈ R * + × ∂Ω, u(0, x) = u 0 (x), for x ∈ Ω, (6)
where

u 0 := (u 0 1 , • • • , u 0 n ) T and the matrix-valued application A : [0, 1] n → R n×n u := (u i ) 1≤i≤n → (A ij (u)) 1≤i,j≤n is defined by ∀1 ≤ i ≤ n, A ii (u) = 1≤j =i≤n (K ij -K i0 )u j + K i0 , ∀1 ≤ i = j ≤ n, A ij (u) = -(K ij -K i0 )u i . (7)
Despite their importance in chemistry or biology, it appears that the mathematical analysis of systems of the form [START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF], taking into account constraints [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF], is quite recent [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF][START_REF] André | Local existence, uniqueness and smooth dependence for nonsmooth quasilinear parabolic problems[END_REF][START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF][START_REF] Liero | Gradient structures and geodesic convexity for reaction-diffusion systems[END_REF]. Let us point out here that the non-negativity of the solutions to (6) through time is a mathematical issue, linked to the absence of a maximum principle for such systems.

At least up to our knowledge, the first proof of existence of global weak solutions of (6) satisfying constraints [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF] with non-identical cross-diffusion coefficients is given in [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF] for n = 2 with coefficients K ij such that K i0 > 0 for i = 1, 2 and K 12 = K 21 = 0. These results were later extended in [START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF] to a general number of species n ∈ N * with crossdiffusion coefficients satisfying K i0 > 0 and K ij = 0 for all 1 ≤ i = j ≤ n; the authors of the latter article proved in addition the uniqueness of such weak solutions. In [START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF], the case n = 2 with arbitrary positive coefficients K ij > 0 is covered, though no uniqueness result is provided. The main difficulty of the mathematical analysis of such equations relies in the bounds [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF], which are not obvious since no maximum principle can be proved for these systems in general. In all the articles mentioned above, the analysis framework used by the authors is the so-called boundedness by entropy method. The main idea of this technique is to write the above system of equations as a formal gradient flow and derive estimates on the solutions (u 0 , • • • , u n ) using the decay of some well-chosen entropy functional. We present in Section 1.1.2 the formal gradient flow structure of (6) and recall the results of [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] in Section 1.2.

Remark 1. This model is linked to the so-called Stefan-Maxwell model, studied in [START_REF] Jungel | Existence analysis of maxwell-stefan systems for multicomponent mixtures[END_REF][START_REF] Boudin | A mathematical and numerical analysis of the maxwell-stefan diffusion equations[END_REF]. Indeed, the model considered in the latter paper reads

   ∂ t u -div x A(u) -1 ∇ x u = 0, for (t, x) ∈ (0, T ] × Ω, (A(u)∇ x u) • n = 0, for (t, x) ∈ (0, T ] × ∂Ω, u(0, x) = u 0 (x), for x ∈ Ω, (8)
where A is defined by [START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF].

1.1.2 Formal gradient flow structure of [START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF] We detail in this section the formal gradient flow structure of the system [START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF]. Let D ⊂ R n be defined by

D := (u 1 , • • • , u n ) ∈ (R * + ) n , n i=1 u i < 1 ⊂ (0, 1) n . ( 9
)
Let us introduce the classical entropy density h (see for instance [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF], [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF], [START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF] and [START_REF] Liero | Gradient structures and geodesic convexity for reaction-diffusion systems[END_REF])

h :    D -→ R u := (u i ) 1≤i≤n -→ h(u) = n i=1 u i log u i + (1 -ρ u ) log(1 -ρ u ), (10)
where ρ u := n i=1 u i . Some properties of h can be easily checked:

(P1) the function h belongs to C 0 (D) ∩ C 2 (D); consequently, h is bounded on D;

(P2) the function h is strictly convex on D;

(P3) its derivative Dh :

D -→ R n (u i ) 1≤i≤n → log u i 1-ρu 1≤i≤n
, is invertible and its inverse is given by

(Dh) -1 : R n -→ D (w i ) 1≤i≤n → e w i 1+ n j=1 e w j .
In the following, we denote by D 2 h the Hessian of h. The entropy functional E is defined by

E : L ∞ (Ω; D) -→ R u -→ E(u) := ´Ω h(u(x)) dx. (11)
Throughout the article, for all u ∈ L ∞ (Ω; D), we shall denote by DE(u) the measurable vector-valued function defined by

DE(u) : Ω → R n x → Dh(u(x)).
The system (6) can then be formally rewritten under the following gradient flow structure

   ∂ t u -div x (M (u)∇ x DE(u)) = 0, for (t, x) ∈ R * + × Ω, (M (u)∇ x DE(u)) • n = 0, for (t, x) ∈ R * + × ∂Ω, u(0, x) = u 0 (x), for x ∈ Ω, (12)
where M : D → R n×n is the so-called mobility matrix of the system defined for all u ∈ D by

M (u) := A(u)(D 2 h(u)) -1 .
More precisely, it holds that for all u ∈ D, M (u) = (M ij (u)) 1≤i,j≤n where for all

1 ≤ i = j ≤ n, M ii (u) = K i0 (1 -ρ u )u i + 1≤j =i≤n K ij u i u j and M ij (u) = -K ij u i u j . (13)

Existence of global weak solutions by the boundedness by entropy technique

The formal gradient flow formulation of a system of cross-diffusion equations is a key point in the boundedness by entropy technique. In the example presented in Section 1.1, it implies in particular that E is a Lyapunov functional for the system (6) [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF][START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF]. However, the mobility matrix obtained for these systems is not a concave function of the densities, so that standard gradient flow theory arguments (such as the minimizing movement method) cannot be applied in this context [START_REF] Zinsl | Transport distances and geodesic convexity for systems of degenerate diffusion equations[END_REF][START_REF] Dolbeault | A new class of transport distances between measures[END_REF][START_REF] Jordan | The variational formulation of the fokker-planck equation[END_REF][START_REF] Liero | Gradient structures and geodesic convexity for reaction-diffusion systems[END_REF]. However, the existence of a global weak solution to (6) can still be proved. Let us recall here a simplified version of Theorem 2 of [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] which is adapted to our context.

Theorem 1 (Theorem 2 of [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF]). Let D ⊂ R n be the domain defined by [START_REF] Dolbeault | A new class of transport distances between measures[END_REF]. Let A : u ∈ D → A(u) := (A ij (u)) 1≤i,j≤n ∈ R n×n be a matrix-valued functional defined on D satisfying A ∈ C 0 (D; R n×n ) and the following assumptions:

(H1) There exists a bounded from below convex function h ∈ C 2 (D, R) such that its derivative Dh : D → R n is invertible on R n ;

(H2) There exists α > 0, and for all

1 ≤ i ≤ n, there exist 1 ≥ m i > 0, such that for all z = (z 1 , • • • , z n ) T ∈ R n and u = (u 1 , • • • , u n ) T ∈ D, z T D 2 h(u)A(u)z ≥ α n i=1 u 2m i -2 i z 2 i .
Let u 0 ∈ L 1 (Ω; D) so that w 0 := Dh(u 0 ) ∈ L ∞ (Ω; R n ). Then, there exists a weak solution u with initial condition u 0 to

∂ t u = div x (A(u)∇ x u), for (t, x) ∈ R * + × Ω, (A(u)∇ x u) • n = 0, for (t, x) ∈ R * + × ∂Ω, (14)
such that for almost all

(t, x) ∈ R * + × Ω, u(t, x) ∈ D with u ∈ L 2 loc (R + ; H 1 (Ω, R n )) and ∂ t u ∈ L 2 loc (R + ; (H 1 (Ω; R n )) )
. Lemma 1 states that the prototypical example presented in Section 1.1 falls into the framework of Theorem 1. The proof of the latter is given Section 3.1 for the sake of completeness, and relies on ideas introduced in [START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF].

Lemma 1. Let D ⊂ R n be the domain defined by [START_REF] Dolbeault | A new class of transport distances between measures[END_REF] and A : u ∈ D → A(u) := (A ij (u)) 1≤i,j≤n ∈ R n×n be the matrix-valued function defined by [START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF]. Then, A ∈ C 0 (D; R n×n ) and satisfies assumptions (H1)-(H2) of Theorem 1, with h given by [START_REF] Dreher | Compact families of piecewise constant functions in lp (0, t; b)[END_REF], α = min 1≤i =j≤n K ij and m i = 1 2 for all 1 ≤ i ≤ n. The existence of global weak solutions to (6) is then a direct consequence of Theorem 1 and Lemma 1.

Let us point out that the uniqueness of solutions to general systems of the form [START_REF] Jordan | The variational formulation of the fokker-planck equation[END_REF] remains an open theoretical question, at least up to our knowledge. It can be obtained in some particular cases. When the diffusion matrix A is defined by [START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF] and when all the diffusion coefficients K ij are identically equal to some constant K > 0, the uniqueness of the solution can be trivially obtained since the system boils down to a set of n decoupled heat equation for the evolution of the density of each species. In the sequel, we restrict the study to the case when d = 1. In this section, we propose a model for the description of a PVD process and present theoretical results whose proofs are postponed to Section 3. The global existence of a weak solution is proved. The long-time behaviour of such a solution is studied in the case of constant external fluxes. Lastly, under the assumption that the coefficients K ij are chosen so that there is a unique solution to the system, we prove the existence of a solution to an optimization problem.

Presentation of the model

For the sake of simplicity, we assume that non-zero fluxes are only imposed on the righthand side of the domain occupied by the solid. At some time t > 0, this domain is denoted by Ω t := (0, e(t)) where e(t) > 0 models the thickness of the layer. Initially, we assume that the domain Ω 0 occupied by the solid at time t = 0 is the interval (0, e 0 ) for some initial thickness e 0 > 0.

The evolution of the thickness of the film e(t) is determined by the external fluxes of the atomic species that are absorbed at its surface. More precisely, let us assume that there are n + 1 different chemical species composing the solid layer and let

(φ 0 , • • • , φ n ) belong to L ∞ loc (R + ; R n+1 + ).
For all 0 ≤ i ≤ n, the function φ i (t) represents the flux of the species i absorbed at the surface at time t > 0 and is assumed to be non-negative. In this one-dimensional model, the evolution of the thickness of the solid is assumed to be given by e(t)

:= e 0 + ˆt 0 n i=0 φ i (s) ds. (15)
In the following, we will denote by 1).

ϕ := (φ 1 , • • • , φ n ) T (see Figure
For all t ≥ 0 and 0 ≤ i ≤ n, the local concentration of species i at time t and point x ∈ (0, e(t)) is denoted by u i (t, x). The evolution of the vector u := (u 1 , • • • , u n ) is given by the system of cross-diffusion equations

∂ t u -∂ x (A(u)∂ x u) = 0, for t ∈ R * + , x ∈ (0, e(t)), (16)
where A : D → R n×n is a well-chosen diffusion matrix satisfying (H1)-(H2).

We consider that for every t > 0, the system satisfies the following conditions on the boundary ∂Ω t :

(A(u)∂ x u) (t, 0) = 0 and (A(u)∂ x u) (t, e(t)) + e (t)u(t, e(t)) = ϕ(t). (17)
An easy calculation shows that these boundary conditions, in addition to ( 15) and ( 16), ensure that, for all 0

≤ i ≤ n, d dt ˆΩt u i (t, x) dx = φ i (t).
Indeed, it holds that

d dt ˆΩt u(t, x) dx = ˆe(t) 0 ∂ t u(t, x) dx + e (t)u(t, e(t)), = ˆe(t) 0 ∂ x (A(u)∂ x u) + e (t)u(t, e(t)), = (A(u)∂ x u)(t, e(t)) + e (t)u(t, e(t)) -(A(u)∂ x u)(t, 0), = ϕ(t).
The calculation for the 0 th species reads:

d dt ˆΩt u 0 (t, x) dx = d dt |Ω t | - n i=1 ˆΩt u i (t, x) dx = e (t) - n i=1 d dt ˆΩt u i (t, x) dx = n i=0 φ i (t) - n i=1 φ i (t) = φ 0 (t).
To sum up, the final system of interest reads:

           e(t) = e 0 + ´t 0 n i=0 φ i (s) ds, for t ∈ R * + , ∂ t u -∂ x (A(u)∂ x u) = 0, for t ∈ R * + , x ∈ (0, e(t)), (A(u)∂ x u) (t, 0) = 0, for t ∈ R * + , (A(u)∂ x u) (t, e(t)) + e (t)u(t, e(t)) = ϕ(t), for t ∈ R * + , u(0, x) = u 0 (x), for x ∈ (0, e 0 ), (18)
where u 0 ∈ L 1 (0, e 0 ) is an initial condition satisfying u 0 (x) ∈ D for almost all x ∈ (0, e 0 ). We assume in addition that w 0 := Dh(u 0 ) belongs to L ∞ ((0, e 0 ); R n ).

Rescaled version of the model

We introduce here a rescaled version of system [START_REF] Horst | Invariant regions for quasilinear reaction-diffusion systems and applications to a two population model[END_REF]. For all 0 ≤ i ≤ n, t ≥ 0 and y ∈ (0, 1), let us denote by v i (t, y) := u i (t, e(t)y). It holds that

∂ t v(t, y) = ∂ t u(t,
where v := (v 1 , • • • , v n ).
Thus, u is a solution of (18) if and only if v is a solution to the following system:

                   e(t) = e 0 + ˆt 0 n i=0 φ i (s) ds, for t ∈ R * + , ∂ t v -1 e(t) 2 ∂ y (A(v)∂ y v) -e (t) e(t) y∂ y v = 0, for (t, y) ∈ R * + × (0, 1), 1 e(t) (A(v)∂ y v)(t, 1) + e (t)v(t, 1) = ϕ(t), for (t, y) ∈ R * + × (0, 1), 1 e(t) (A(v)∂ y v)(t, 0) = 0, for (t, y) ∈ R * + × (0, 1) v(0, y) = v 0 (y), for y ∈ (0, 1), (19)
where v 0 (y) := u 0 (e 0 y).

Proving the existence of a global weak solution to ( 18) is equivalent to proving the existence of a global weak solution to [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF].

Actually, it can be seen that the entropy of the system (19) satisfies a formal inequality at the continuous level which is at the heart of the proof of our existence result. Indeed, let us denote by

E(t) := ˆ1 0 h(v(t, y)) dy,
where v is a solution to [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF]. Then, formal calculations yield that

dE dt (t) = ˆ1 0 ∂ t v(t, y) • Dh(v(t, y)) dy = 1 e(t) 2 ˆ1 0 ∂ y (A(v(t, y))∂ y v(t, y)) • Dh(v(t, y)) dy + e (t) e(t) ˆ1 0 y∂ y v(t, y) • Dh(v(t, y)) dy = - 1 e(t) 2 ˆ1 0 ∂ y v(t, y) • D 2 h(v(t, y))A(v(t, y))∂ y v(t, y) dy + 1 e(t) 2 (A(v(t, 1))∂ y v(t, 1)) • Dh(v(t, 1)) + e (t) e(t) ˆ1 0 y∂ y (h(v(t, y))) dy = - 1 e(t) 2 ˆ1 0 ∂ y v(t, y) • D 2 h(v(t, y))A(v(t, y))∂ y v(t, y) dy + 1 e(t) ϕ(t) -e (t)v(t, 1) • Dh(v(t, 1)) + e (t) e(t) h(v(t, 1)) - e (t) e(t) ˆ1 0 h(v(t, y)) dy.
Denoting by f (t) := ϕ(t) e (t) , it holds that f (t) ∈ D for all t > 0. Besides, using assumption (H2), we obtain that

- ˆ1 0 ∂ y v(t, y) • D 2 h(v(t, y))A(v(t, y))∂ y v(t, y) dy ≤ 0, which yields that dE dt (t) ≤ e (t) e(t) h(v(t, 1) + Dh(v(t, 1)) • f (t) -v(t, 1) - ˆ1 0 h(v(t, y)) dy .
Using the convexity of h, we obtain that h(v(t, 1) + Dh(v(t, 1))

• f (t) -v(t, 1) ≤ h(f (t)), so that dE dt (t) ≤ e (t) e(t) h(f (t)) - ˆ1 0 h(v(t, y)) dy . (20)
Inequality [START_REF] Le | Everywhere regularity of solutions to a class of strongly coupled degenerate parabolic systems[END_REF] is not an entropy dissipation inequality in the sense that the quantity E(t) may increase with time. However, using the fact e ∈ L ∞ loc (R + ; R + ) and assumption (H3), it implies that the quantity E(t) cannot blow up in finite time, which is sufficient for our purpose.

Theoretical results

Global in time existence of weak solutions

Our first result deals with the global in time existence of bounded weak solutions to [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] (and thus to [START_REF] Horst | Invariant regions for quasilinear reaction-diffusion systems and applications to a two population model[END_REF]).

Theorem 2. Let D := {(u 1 , • • • , u n ) T ∈ (R * + ) n , n i=1 u i < 1} ⊂ (0, 1) n .
Let A : D → R n×n be a matrix-valued functional satisfying A ∈ C 0 (D; R n×n ) and assumptions (H1)-(H2) of Theorem 1 for some well-chosen entropy density h : D → R. We assume in addition that (H3) h ∈ C 0 (D).

Let e 0 > 0, u 0 ∈ L 1 ((0, e 0 ); D) so that w

0 := (Dh) -1 (u 0 ) ∈ L ∞ ((0, e 0 ); R n ) and (φ 0 , • • • , φ n ) ∈ L ∞ loc (R + ; R n+1 + ).
Let us define for almost all y ∈ (0, 1), v 0 (y) := u 0 (e 0 y) and ϕ := (φ 1 , • • • , φ n ) T . Then, there exists a weak solution v with initial condition v 0 to [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] such that for almost all (t, y) ∈ R * + × (0, 1), v(t, y) ∈ D. Besides,

v ∈ L 2 loc (R + ; H 1 ((0, 1); R n )) and ∂ t v ∈ L 2 loc (R + ; (H 1 ((0, 1); R n )) ).
In particular, v ∈ C 0 (R + ; L 2 ((0, 1); R n )).

Let us point out that the example described in Section 1.1 satisfies all the assumptions of Theorem 2 since the entropy density h defined by [START_REF] Dreher | Compact families of piecewise constant functions in lp (0, t; b)[END_REF] belongs to C 0 (D). Let us also point here that the form of ( 19) is different from the system considered in [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] through i) the boundary conditions and ii) the existence of the drift term e (t) e(t) y∂ y v. The strategy of proof developped in [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF][START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] is still adapted to our case though, because a discrete entropy inequality can still be obtained. The proof of Theorem 2 is given in full details in Section 3.2.

Long-time behaviour for constant fluxes

In the case when the fluxes are constant in time, we obtain long-time asymptotics for the functions v i , provided that the entropy density h is given by [START_REF] Dreher | Compact families of piecewise constant functions in lp (0, t; b)[END_REF]. More precisely, the following result holds: Proposition 1. Under the assumptions of Theorem 2, let us make the following additional hypotheses:

(T1) for all 0 ≤ i ≤ n, there exists φ i > 0 so that φ i (t) = φ i , for all t ∈ R + ;

(T2) for all u ∈ D, the entropy density h can be chosen so that h(u)

= n i=1 u i log u i + (1 -ρ u ) log(1 -ρ u ).
For all 0 ≤ i ≤ n, let us define f i :=

h : D → R u → n i=1 u i ln u i f i + (1 -ρ u ) ln 1-ρu 1-ρ f
the relative entropy associated to h and f . Then, there exists a global weak solution v to [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] and a constant C > 0 such that

ˆ1 0 h (v(t, y)) dy ≤ C t + 1 , (21)
and

∀1 ≤ i ≤ n, v i (t, •) -f i L 1 (0,1) ≤ C √ t + 1 and 1 -ρ v(t,•) -f 0 L 1 (0,1) ≤ C √ t + 1 . ( 22
)
The proof of Proposition 1 is given in Section 3.3. Numerical results presented in Section 4 illustrate the rate of convergence of the rescaled concentrations to constant profiles in O 1 t . Let us comment here on assumption (T2). Actually, in the proof, we use the following properties of the logarithmic entropy density:

• h ≥ 0 and if u ∈ D satisfies h(u) = 0, then necessarily u = f ;
• There exists a constant vector g ∈ R n such that for all u ∈ D, Dh(u) = Dh(u) + g.

• The relative entropy density h satisfies a Csizàr-Kullback type inequality.

Similar long-time asymptotics results can be obtained for general entropy densities satisfying these three properties. For the sake of simplicity, we chose to restrict ourselves to the case of logarithmic entropy density in Proposition 1.

The central ingredient of the proof is the following formal entropy inequality. In the case when h is given by [START_REF] Dreher | Compact families of piecewise constant functions in lp (0, t; b)[END_REF], it can be easily seen that h is also a valid entropy density for the diffusion coefficient A in the sense that h also satisfies assumptions (H1)-(H2)-(H3). Thus, inequality [START_REF] Le | Everywhere regularity of solutions to a class of strongly coupled degenerate parabolic systems[END_REF] holds with h instead of h so that

dE dt (t) ≤ e (t) e(t) h(f ) - ˆ1 0 h(v(t, y)) dy = e (t) e(t) h(f ) -E(t) ,
where for all t > 0, E(t) := ´1 0 h(v(t, y)) dy. Denoting by V := n i=0 φ i , it holds that e (t) = V and e(t) = e 0 + V t for all t ≥ 0. Finally, using the fact that h ≥ 0 and that h(f ) = 0, we obtain that

e 0 V + t dE dt (t) + E(t) = d dt e 0 V + t E(t) ≤ 0.
This inequality implies that there exists a constant C > 0 such that for all t ≥ 0,

E(t) ≤ C t + 1 .
The rates on the L 1 norm of the solutions are then obtained using the Csizàr-Kullback inequality.

Optimization of the fluxes

As mentioned in the introduction, our main motivation for studying this system is the control of the gazeous fluxes injected during a PVD process. It is assumed here that the wafer remains in the hot chamber where the different atomic species are injected during a time T > 0. The cross-diffusion phenomena occur in the bulk of the thin film layer because of the high temperatures that are imposed during the process. Once the wafer is taken out of the chamber, the composition of the film is freezed and no diffusion phenomena take place anymore. The profiles of the local volumic fractions of the different chemical species in the film thus remain unchanged after the time T . It is of practical interest to adapt the fluxes through time so that these final concentration profiles are as close as possible to target functions chosen a priori.

Let e 0 > 0 be the initial thickness of the solid. In practice, the maximal value of the fluxes which can be injected is limited due to device constraints. Let F > 0 and let us then denote by Let us point out here the uniqueness of a solution to (18) (or [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF]) remains an open problem in general. When the diffusion matrix A is defined by [START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF], the only case for which uniqueness of a global solution can be obtained is the trivial case where the cross-diffusion coefficients K ij are identical to some constant K > 0 for all 0 ≤ i = j ≤ n. Indeed, in this case, it can be seen that the system [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] can be written as a set of n independent advectiondiffusion PDEs on each of the rescaled concentration profiles v i (1 ≤ i ≤ n). Thus, we will have to make some assumption on the cross-diffusion coefficients (K ij ) 0≤i =j≤n in the general case.

Ξ := Φ ∈ L ∞ ((0, T ); R n+1 + ), Φ L ∞ ≤ F . For all Φ := (φ 0 , • • • , φ n ) ∈ Ξ,
We make the following assumption on the diffusion matrix A:

(C1) For any Φ ∈ Ξ, there exists a unique global weak solution v Φ to system [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] so that for almost all (t, y) ∈ R * + × (0, 1), v Φ (t, y) ∈ D.

The goal of the optimization problem consists in the identification of optimal timedependent non-negative functions Φ ∈ Ξ so that the final thickness of the film e Φ (T ) and the (rescaled) concentration profiles for the different chemical species v Φ (T, •) at the end of the fabrication process are as close as possible to desired targets denoted by e opt > e 0 and v opt ∈ L 2 ((0, 1); D).

The real-valued functional J : Ξ → R defined by

∀Φ ∈ Ξ, J (Φ) := |e Φ (T ) -e opt | 2 + v Φ (T, •) -v opt 2 L 2 (0,1) , (23)
is the cost function we consider here. More precisely, we have the following result, which is proved in Section 3.4.

Proposition 2. Under the assumptions of Theorem 2, let us make the additional assumption (C1). Then, the functional J is well-defined and there exists a minimizer Φ * ∈ Ξ to the minimization problem

Φ * ∈ argmin Φ∈Ξ J (Φ). (24)
Of course, uniqueness of such a solution Φ * is not expected in general.

Proofs

Proof of Lemma 1

Let us prove that the matrix-valued function A defined in [START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF] satisfies the assumptions of Theorem 1 with the entropy functional h given by [START_REF] Dreher | Compact families of piecewise constant functions in lp (0, t; b)[END_REF].

As mentioned in Section 1.1, the entropy density h belongs to C 0 (D; R)∩C 2 (D; R) (thus is bounded on D), is strictly convex on D, and its derivative Dh : D → R n is invertible. As a consequence, h satisfies assumption (H1) of Theorem 1.

Let us now prove that assumption (H2) of Theorem 1 is satisfied with m i = 1 2 for all 1 ≤ i ≤ n. To this aim, we follow the same strategy of proof as the one used in [START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF]. Let us prove that there exists β > 0 such that for all u ∈ D,

H(u)A(u) ≥ βΛ(u), (25)
where

H(u) := D 2 h(u), Λ(u) := diag 1 u i 1≤i≤n and β := min 0≤i =j≤n K ij .
This inequality implies (H2) with α = β and

m i = 1 2 for all 1 ≤ i ≤ n. Let u ∈ D. We have for all 1 ≤ i, j ≤ n, H ii (u) = 1 u i + 1 1 -ρ u and H ij (u) = 1 1 -ρ u if i = j.
Introducing P (u) := (P ij (u)) 1≤i,j≤n , where for all 1 ≤ i, j ≤ n,

P ii (u) = 1 -u i and P ij (u) = -u i if i = j,
it holds that H(u)P (u) = Λ(u). Thus, H(u)A(u) -βΛ(u) = H(u)(A(u) -βP (u)). It can be easily checked that A(u) -βP (u) = A(u) + βD(u), where A(u) has the same structure as A(u) but with diffusion coefficients K ij -β instead of K ij , and D(u) := (D ij (u)) 1≤i,j≤n where D ij (u) = u i for all 1 ≤ i ≤ n.

On the one hand, H(u)D(u) = 1 1-ρu Z where Z is the n×n matrix whose all coefficients are identically equal to 1. Since the matrix Z is a semi-definite positive matrix, so is H(u)D(u).

On the other hand, since h is strictly convex on D, H(u) A(u) is semi-definite positive if and only if M (u) := A(u)H(u) -1 is semi-definite positive. Indeed, for all z ∈ R n , we have z T H(u) A(u)z = (H(u)z) T A(u)H(u) -1 (H(u)z). It can be observed that M (u) = ( M ij (u)) 1≤i,j≤n , where for all 1 ≤ i, j ≤ n,

M ii (u) = (K i0 -β)(1-ρ u )u i + 1≤j =i≤n (K ij -β)u i u j and M ij (u) = -(K ij -β)u i u j if j = i. For all z = (z 1 , • • • , z n ) T ∈ R n , we have z T M (u)z = n i=1 (K i0 -β)(1 -ρ u )u i z 2 i + n i=1 1≤j =i≤n (K ij -β)u i u j (z 2 i -z i z j ), = n i=1 (K i0 -β)(1 -ρ u )u i z 2 i + 1≤i =j≤n (K ij -β)u i u j 1 2 z 2 i + 1 2 z 2 j -z i z j , ≥ 0.
The matrix M (u) is indeed a semi-definite positive matrix. Hence we have proved inequality [START_REF] Painter | Continuous models for cell migration in tissues and applications to cell sorting via differential chemotaxis[END_REF], which yields the desired result.

Proof of Theorem 2

For the sake of simplicity, we will prove the existence of a solution v on the finite time interval [0, T ] where T > 0 is an arbitrary positive constant. Actually, the proof can be easily adapted to obtain the existence of a global solution for an infinite time horizon.

The proof follows similar lines as the proof of Theorem 2 of [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] and is divided in three main steps. Firstly, an approximate time-discrete problem is introduced for which uniform bounds are proved in a second step. Lastly, passing to the limit in this approximate problem using the obtained bounds enables to obtain the existence of a weak solution.

Step 1 : Approximate time-discrete problem

Let us first assume at this point that φ 0 , • • • , φ n belong to C 0 ([0, T ]).

Let N ∈ N, τ = T N and > 0. For all k ∈ N * so that kτ ≤ T , let us denote by e k := e(kτ ), e k := e (kτ ) and

ϕ k = (φ 1,k , • • • , φ n,k ) T := ϕ(kτ ). Let us also define f k := ϕ k e k if e k > 0, 0 otherwise, ( 26
) so that f k ∈ D and ϕ k = e k f k .
By assumption, w 0 (y) := Dh(v 0 (y)) belongs to L ∞ ((0, 1); R n ). In the rest of the proof, for any w ∈ R n , we denote by v(w) := (Dh) -1 (w) = (v i (w)) 1≤i≤n and by B(w) := M (v(w)).

Let us already mention at this point that the (formal) weak formulation of (19) reads as follows: for all ψ ∈ L 2 ((0, T ); H 1 ((0, 1);

R n )), ˆT 0 ˆ1 0 ∂ t v • ψ + ˆT 0 ˆ1 0 ∂ y 1 e 2 ψ • (A(v)∂ y v) + ˆT 0 ˆ1 0 e e (v • ψ + yv • ∂ y ψ) = ˆT 0 1 e ϕ • ψ(•, 1).
Let us first prove the following lemma.

Lemma 2. Assume that φ 0 , • • • , φ n ∈ C 0 ([0, T ]).
Then, for all k ∈ N * such that kτ ≤ T , there exists w k ∈ H 1 ((0, 1); R n ) solution of

1 τ ˆ1 0 v(w k ) -v(w k-1 ) • ψ + 1 e 2 k ˆ1 0 ∂ y ψ • (B(w k )∂ y w k ) + ˆ1 0 (∂ y w k • ∂ y ψ + w k • ψ) (27)
+ e k e k ˆ1 0 (v(w k ) • ψ + yv(w k ) • ∂ y ψ) = 1 e k ϕ k • ψ(1),
for all ψ ∈ H 1 ((0, 1); R n ). Besides, the following discrete inequality holds for all k ∈ N * such that kτ ≤ T ,

1 τ ˆ1 0 h(v(w k )) + ˆ1 0 |∂ y w k | 2 + |w k | 2 + 1 e 2 k ˆ1 0 ∂ y w k • (B(w k )∂ y w k ) (28) ≤ 1 τ ˆ1 0 h(v(w k-1 )) + e k e k h(f k ) - ˆ1 0 h(v(w k )) .
The proof of this lemma is postponed until Section 3.2.4. Let us point out the following fact: from (28), we obtain

1 τ + e k e k ˆ1 0 h(v(w k )) + ˆ1 0 (|∂ y w k | 2 + |w k | 2 ) + 1 e 2 k ˆ1 0 ∂ y w k • B(w k )∂ y w k (29)
≤ 1 τ ˆ1 0 h(v(w k-1 )) + e k e k h L ∞ (D) ,
which implies

1 τ ˆ1 0 h(v(w k )) + ˆ1 0 (|∂ y w k | 2 + |w k | 2 ) + 1 e 2 k ˆ1 0 ∂ y w k • B(w k )∂ y w k (30) ≤ 1 τ ˆ1 0 h(v(w k-1 )) + 2 e k e k h L ∞ (D) .

Step 2: Uniform bounds

For all 0 ≤ i ≤ n, let (φ i,p ) p∈N be a sequence of non-negative functions of C 0 ([0, T ]) which weakly-* converges to φ i in L ∞ (0, T ) as p goes to infinity, and for all p ∈ N, It holds that (e p ) p∈N * strongly converges to e in L ∞ (0, T ). Indeed, let ε > 0. Since e is continuous on [0, T ], it is uniformly continuous, and there exists η > 0 so that for all

φ i,p L ∞ (0,T ) ≤ φ i L ∞ (0,T ) .

Let us define

t, t ∈ [0, T ] satisfying |t -t | ≤ η, then |e(t) -e(t )| ≤ ε/2. Let M ∈ N * and 0 = s 0 < s 1 < • • • < s M = T so that for all 0 ≤ j ≤ M -1, |s j -s j+1 | ≤ η. Then, it holds that max 0≤j≤M |e p (s j ) -e(s j )| -→ p→+∞ 0,
because of the weak-* convergence in L ∞ [0, T ] of (φ i,p ) p∈N * to φ i for all 0 ≤ i ≤ n. Thus, there exists p 0 ∈ N * large enough such that for all p ≥ p 0 , max 0≤j≤M |e p (s j )-e(s j )| ≤ ε/2. Besides, the non-negativity of the functions φ i and φ i,p implies that e and e p are nondecreasing functions, so that for all 0 ≤ j ≤ M -1 and all p ∈ N * , ∀s ∈ [s j , s j+1 ], e(s j ) ≤ e(s) ≤ e(s j+1 ) and e p (s j ) ≤ e p (s) ≤ e p (s j+1 ).

As a consequence, for all p ≥ p 0 , all 0 ≤ j ≤ M -1 and all s ∈ [s j , s j+1 ], Hence, for all p ≥ p 0 , ee p L ∞ (0,T ) ≤ ε, which yields the strong convergence of the sequence (e p ) p∈N * to e in L ∞ (0, T ).

For all k ∈ N * such that kτ ≤ T , we denote by w k,p a solution to ( 27) associated to the fluxes (φ i,p ) 0≤i≤n . The time-discretized associated quantities are denoted (using obvious notation) by ϕ k,p , e k,p and e k,p .

Let us define the piecewise constant in time functions w ( ,τ,p) (y, t), v ( ,τ,p) (y, t), σ τ v ( ,τ,p) (y, t), e (τ,p) (t) and e d (τ,p) (t) as follows: for all k ≥ 1 such that kτ ≤ T , (k -1)τ < t ≤ kτ and almost all y ∈ (0, 1), w ( ,τ,p) (y, t) := w k,p (y), v ( ,τ,p) (y, t) := Dh(w k,p (y)), σ τ v ( ,τ,p) (y, t) = Dh(w k-1,p (y)), e (τ,p) (t) = e k,p , e d (τ,p) (t) := e k,p , ϕ (τ,p) := ϕ k,p .

Besides, let us set w ( ,τ,p) (0, •) = Dh(v 0 ) and v ( ,τ,p) (0, •) = v 0 . Let us also denote by (v

( ,τ,p) 1 , • • • , v ( ,τ,p) n
) the n components of v ( ,τ,p) .

Then, the following system holds for all piecewise constant in time functions ψ : (0, T ) → H 1 ((0, 1); R n ),

1 τ ˆT 0 ˆ1 0 v ( ,τ,p) -σ τ v ( ,τ,p) ) • ψ dy dt + ˆT 0 1 e 2 (τ,p) ˆ1 0 ∂ y ψ • (B(w ( ,τ,p) )∂ y w ( ,τ,p) ) dy dt (31)
+ ˆT 0 ˆ1 0 (∂ y w ( ,τ,p) • ∂ y ψ + w ( ,τ,p) • ψ) dy dt + ˆT 0 e d (τ,p) e (τ,p) ˆ1 0 v(w ( ,τ,p) ) • ψ + yv(w ( ,τ,p) ) • ∂ y ψ) dy dt = ˆT 0 1 e (τ,p) ϕ (τ,p) • ψ(1) dt.
The set of piecewise constant functions in time ψ : (0, T ) → H 1 ((0, 1); R n ) is dense in L 2 ((0, T ); H 1 ((0, 1); R n )), so that (31) also holds for any ψ ∈ L 2 ((0, T ); H 1 ((0, 1); R n )).

Using the fact that A satisfies assumption (H2) of Theorem 1 and the fact that

∂ y w k,p = D 2 h(v k,p )∂ y v k,p , we obtain for all k ∈ N * such that kτ ≤ T , ˆ1 0 ∂ y w k,p • (B(w k,p )∂ y w k,p ) = ˆ1 0 ∂ y v(w k,p ) • D 2 h(v(w k,p ))A(v(w k,p ))∂ y v(w k,p ) dy ≥ n i=1 ˆ1 0 α v i (w k,p ) 2m i -2 |∂ y v i (w k,p )| 2 dy = n i=1 ˆ1 0 |∂ y G i (v i (w k,p ))| 2 dy = ˆ1 0 |∂ y G(v(w k,p ))| 2 dy,
where G i (s) := √ α m i |s| m i for all s ∈ (0, 1) and G(z) = (G i (z i )) 1≤i≤n for all z := (z i ) 1≤i≤n ∈ (0, 1) n . It follows from [START_REF] Zamponi | Analysis of degenerate cross-diffusion population models with volume filling[END_REF] that for all k ∈ N * such that kτ ≤ T ,

ˆ1 0 h(v(w k,p )) + τ ˆ1 0 |∂ y α(v(w k,p ))| 2 + τ ˆ1 0 |∂ y w k,p | 2 + |w k,p | 2 ≤ 2τ h L ∞ (D) e k,p e k,p + ˆ1 0 h(v(w k-1,p )).
Summing these inequalities yields, for k ∈ N * so that kτ ≤ T ,

ˆ1 0 h(v(w k,p )) + τ k j=1 ˆ1 0 |∂ y G(v(w j,p ))| 2 + τ k j=1 ˆ1 0 (|∂ y w j,p | 2 + |w j,p | 2 ) (32) ≤ 2τ h L ∞ (D) k j=1 e j,p e j,p + ˆ1 0 h(v 0 ), ≤ 2 h L ∞ (D) 1 e 0 k j=1 τ e j,p + ˆ1 0 h(v 0 ), ≤ 2 h L ∞ (D) (n + 1) Φ L ∞ (0,T ) e 0 T + ˆ1 0 h(v 0 ).
In the sequel, C will denote an arbitrary constant, which may change along the calculations, but remains independent on , τ , p and Φ. We are deliberately keeping here the explicit dependence of the constants on Φ L ∞ (0,T ) in view of the proof of Proposition 2. It then holds that e d (τ,p) L ∞ (0,T ) ≤ C Φ L ∞ (0,T ) and 0 < e 0 ≤ e (τ,p) L ∞ (0,T ) ≤ C Φ L ∞ (0,T ) .

We also obtain from (32) and the fact that

G i L ∞ (0,1) ≤ √ α m i for all 1 ≤ i ≤ n that G(v ( ,τ,p) ) L 2 ((0,T );H 1 (0,1) n ) ≤ C 1 + Φ L ∞ (0,T ) (33)
and √ w ( ,τ,p) L 2 ((0,T );H 1 (0,1)

n ) ≤ C 1 + Φ L ∞ (0,T ) . ( 34
)
Since for all 1 ≤ i ≤ n, m i ≤ 1, this implies that

∂ y v ( ,τ,p) i L 2 ((0,T );L 2 (0,1)) = v ( ,τ,p) i 1-m i m i ∂ y |v ( ,τ,p) i | m i L 2 ((0,T );L 2 (0,1)) (35) = v ( ,τ,p) i 1-m i √ α ∂ y G i (v ( ,τ,p) i ) L 2 ((0,T );L 2 (0,1)) ≤ C ∂ y G i (v ( ,τ,p) i ) L 2 ((0,T );L 2 (0,1)) ≤ C 1 + Φ L ∞ (0,T ) .
Besides,

A(v ( ,τ,p) ∂ y v ( ,τ,p) 2 L 2 ((0,T );L 2 (0,1) n ) ≤ A(v ( ,τ,p) ) 2 L ∞ ((0,T );L ∞ (0,1) n×n ) ∂ y v ( ,τ,p) 2 L 2 ((0,T );L 2 (0,1) n ) (36) ≤ C 1 + Φ L ∞ (0,T ) ,
using the fact that A ∈ C 0 (D; R n×n ).

This yields that for all ψ ∈ L 2 ((0, T ); H 1 ((0, 1); R n )),

1 τ ˆT τ ˆ1 0 (v ( ,τ,p) -σ τ v ( ,τ,p) ) • ψ dy dt ≤ 1 e 2 0 A(v ( ,τ,p) ∂ y v ( ,τ,p) L 2 ((0,T );L 2 (0,1) n ) ∂ y ψ L 2 ((0,T );L 2 (0,1) n ) + w ( ,τ,p) L 2 ((0,T );H 1 (0,1) n ) ψ L 2 ((0,T );H 1 (0,1) n ) + 2 e d (τ,p) L ∞ (0,T ) e 0 v ( ,τ,p) L 2 ((0,T );H 1 (0,1) n ) ψ L 2 ((0,T );H 1 (0,1) n ) + 1 e 0 Φ L ∞ (0,T ) ψ L 2 ((0,T );H 1 (0,1) n ) , ≤ C 1 + Φ L ∞ (0,T ) ψ L 2 ((0,T );H 1 (0,1) n ) .
This last inequality shows that These inequalities, together with the fact that (e p ) p∈N * weakly-* converges to e in L ∞ (0, T ) (respectively that (e p ) p∈N * strongly converges to e in L ∞ (0, T )), imply that the sequence e d (τp,p) p∈N * (respectively e (τp,p) p∈N * ) also weakly-* converges to e in L ∞ (0, T ) (respectively strongly converges to e in L ∞ (0, T )).

1 τ v ( ,τ,p) -σ τ v ( ,τ,p) L 2 ((τ,T );(H 1 (0,1) n ) ) ≤ C 1 + Φ L ∞ (0,T ) . ( 37
In the following, we consider such a subsequence (τ p ) p∈N * . The uniform estimates (37) and (35) allow us to apply the Aubin lemma in the version of Theorem 1 of [START_REF] Dreher | Compact families of piecewise constant functions in lp (0, t; b)[END_REF]. Up to extracting a subsequence which is not relabeled, there exists v = (v i ) 1≤i≤n ∈ H 1 ((0, T ); (H 1 ((0, 1); R n )) ) ∩ L 2 ((0, T ); H 1 ((0, 1); R n )) so that as p goes to infinity and goes to 0,

v ( ,τp,p) -→ p→+∞, →0 v,    strongly in L 2 ((0, T ); L 2 ((0, 1); R n )), weakly in L 2 ((0, T ); H 1 ((0, 1); R n )),
and a.e. in (0, T ) × (0, 1),

1 τ p v ( ,τp,p) -σ τp v ( ,τp,p) p→+∞, →0 ∂ t v weakly in L 2 ((0, T ); (H 1 ((0, 1); R n )) ).
Because of the boundedness of v ( ,τp,p) in L ∞ ((0, T ); L ∞ ((0, 1); R n )), the convergence even holds strongly in L q ((0, T ); L q ((0, 1); R n )) for any q < +∞, which is a consequence of the dominated convergence theorem. The latter theorem, together with A ∈ C 0 (D; R n×n ) implies also that the convergence A(v ( ,τp,p) ) -→ A(v) holds strongly in L q ((0, T ); L q ((0, 1); R n×n )). Moreover, using (36) and (34), up to extracting another subsequence, there exists V ∈ L 2 ((0, T ); L 2 ((0, 1); R n )) so that

A(v ( ,τp,p) )∂ y v ( ,τp,p) V weakly in L 2 ((0, T ); L 2 ((0, 1); R n )),
w ( ,τp,p) -→ 0 strongly in L 2 ((0, T ); H 1 ((0, 1); R n )).

The strong convergence of A(v ( ,τp,p) ) in L q ((0, T ); L q ((0, 1); R n )) and the weak convergence of ∂ y v ( ,τp,p) in L 2 ((0, T ); L 2 ((0, 1); R n )) implies necessarily that V = A(v)∂ y v.

We are now in position to pass to the limit → 0 and p → +∞ in [START_REF] Zinsl | Transport distances and geodesic convexity for systems of degenerate diffusion equations[END_REF] with τ = τ p and ψ ∈ L 2 ((0, T ); H 1 ((0, 1); R n )). Let us recall that e (τp,p) p∈N * (respectively e d (τp,p) p∈N * ) converges strongly (respectively weakly-*) to e (respectively e ) in L ∞ (0, T ). We obtain that v is a solution to

ˆT 0 ˆ1 0 ∂ t v • ψ dy dt + ˆT 0 1 e(t) 2 ˆ1 0 ∂ y ψ • (A(v)∂ y v) dy dt (38)
+ ˆT 0 e (t) e(t) ˆ1 0 (v • ψ + yv • ∂ y ψ) dy dt = ˆT 0 1 e(t) ϕ • ψ(1) dt,
yielding the result.

Proof of Lemma 2

Proof of Lemma 2. We prove Lemma 2 by induction using the Leray-Schauder fixed-point theorem. Let z ∈ L ∞ ((0, 1); R n ) and δ ∈ [0, 1]. We consider the following linear problem:

find w ∈ H 1 ((0, 1); R n ) solution of ∀ψ ∈ H 1 ((0, 1); R n ), a z (w, ψ) = l δ,z (ψ), (39)
where

a z (w, ψ) := 1 e 2 k ˆ1 0 ∂ y ψ • B(z)∂ y w + ˆ1 0 (∂ y w • ∂ y ψ + w • ψ) and l δ,z (ψ) := - δ τ ˆ1 0 (v(z) -v(w k-1 )) • ψ + δ e k ϕ k • ψ(1) -δ e k e k ˆ1 0 (v(z) • ψ + yv(z) • ∂ y ψ).
As a consequence of (H2), the matrix B(z) is positive semi-definite for any z ∈ R n . Thus, the bilinear form a z is coercive and continuous on H 1 ((0, 1); R n ), and it holds that

∀ψ ∈ H 1 ((0, 1); R n ), a z (ψ, ψ) ≥ ψ 2 H 1 (0,1) . (40)
Since v(z) ∈ L ∞ ((0, 1); R n ) and v(z) L ∞ (0,1) ≤ 1, the linear form l δ,z is continuous. From the Agmon inequality, there exists C > 0 independent of Φ :

= (φ 0 , • • • , φ n ), or τ such that for all ψ ∈ H 1 ((0, 1); R n ), |l δ,z (ψ)| ≤ 2 τ + C Φ L ∞ (0,T ) ψ H 1 (0,1) , (41)
where

Φ L ∞ (0,T ) = max i=0,••• ,n φ i L ∞ (0,T )
. It immediately follows from the Lax-Milgram theorem that there exists a unique solution w ∈ H 1 ((0, 1); R n ) to (39).

We define the operator S : 39). We are going to prove that there exists a fixed-point w k ∈ H 1 ((0, 1); R n ) of the equation S(1, w k ) = w k using the Leray-Schauder fixed-point theorem (Theorem 3 in the Appendix). This will end the proof of Lemma 2 since such a fixed-point w k is a solution of [START_REF] Redlinger | Invariant sets for strongly coupled reaction-diffusion systems under general boundary conditions[END_REF].

[0, 1] × L ∞ ((0, 1); R n ) → L ∞ ((0, 1); R n ) as follows. For all δ ∈ [0, 1] and χ ∈ L ∞ ((0, 1); R n ), S(δ, χ) is the unique solution w ∈ H 1 ((0, 1); R n ) → L ∞ ((0, 1); R n ) of (
Let us check that all the assumptions of Theorem 3 are satisfied:

(A1) For all χ ∈ L ∞ ((0, 1); R n ), S(0, χ) = 0;

(A2) Let us prove that S is a compact map. To this aim, let us first prove that it is continuous. Let (δ n ) n∈N and (χ n ) n∈N be sequences in [0, 1] and L ∞ ((0, 1); R n ) respectively, δ ∈ [0, 1] and χ ∈ L ∞ ((0, 1); R n ) so that δ n -→ n→+∞ δ and χ n -→ n→+∞ χ strongly in L ∞ ((0, 1); R n ). For all n ∈ N, let w n := S(δ n , χ n ). From assumption (H1) and the global inversion theorem, h :

D → R n is a C 2 -diffeomorphism.
Thus, together with the fact that A ∈ C 0 (D; R n×n ), it holds that the applications

z ∈ R n → v(z) = (Dh) -1 (z) and z ∈ R n → B(z) = A(v(z))D 2 h((Dh) -1 (z)) = A(v(z)D Dh -1 (z) are continuous. Hence, v(χ n ) -→ n→+∞ v(χ) and B(χ n ) -→ n→+∞ B(χ)
strongly in L ∞ ((0, 1); R n ) and L ∞ ((0, 1); R n×n ) respectively.

Besides, the uniform coercivity and continuity estimates (40) and (41) imply that (w n ) n∈N is a bounded sequence in H 1 ((0, 1); R n ). Thus, up to the extraction of a subsequence which is not relabeled, (w n ) n∈N weakly converges to some w in H 1 ((0, 1); R n ).

Passing to the limit n → +∞ in (39) implies that w = S(δ, χ). The uniqueness of the limit yields that the whole sequence (w n ) n∈N weakly converges to S(δ, χ) in H 1 ((0, 1); R n ). The convergence thus holds strongly in L ∞ ((0, 1); R n ) because of the compact embedding H 1 ((0, 1); R n ) → L ∞ ((0, 1); R n ). This proves the continuity of the map S and its compactness follows again from the compact embedding

H 1 ((0, 1); R n ) → L ∞ ((0, 1); R n ).
(A3) Let δ ∈ [0, 1] and w ∈ L ∞ ((0, 1); R n ) so that S(δ, w) = w. It holds that (taking ψ = w as a test function in (39) with χ = w),

1 e 2 k ˆ1 0 ∂ y w • (B(w)∂ y w) + ˆ1 0 (|∂ y w| 2 + |w| 2 ) = (42)
- δ τ ˆ1 0 (v(w) -v(w k-1 )) • w + δ e k ϕ k • w(1) -δ e k e k ˆ1 0 (v(w) • w + yv(w) • ∂ y w). (43)
Let us consider separately the different terms appearing in (43). First, by convexity of h, and using the fact that w = Dh(v(w)), it holds that

δ τ ˆ1 0 (v(w)-v(w k-1 ))•w = δ τ ˆ1 0 (v(w)-v(w k-1 ))•Dh(v(w)) ≥ δ τ ˆ1 0 (h(v(w))-h(v(w k-1 ))).
(44) Besides, using an integration by parts,

δ e k e k ˆ1 0 (v(w) • w + yv(w) • ∂ y w) = δ e k e k v(w)(1) • w(1) - ˆ1 0 yw • ∂ y v(w) , = δ e k e k v(w)(1) • Dh(v(w)(1)) - ˆ1 0 yDh(v(w)) • ∂ y v(w) , = δ e k e k v(w)(1) • Dh(v(w)(1)) - ˆ1 0 y∂ y (h(v(w))) , = δ e k e k v(w)(1) • Dh(v(w)(1)) -h(v(w)(1)) + ˆ1 0 h(v(w)) . (45)
Using [START_REF] Jacobus | Well-posedness of a parabolic movingboundary problem in the setting of wasserstein gradient flows[END_REF], we obtain

δ e k ϕ k • w(1) = δ e k e k f k • Dh(v(w)(1)). (46)
Finally, using (43), ( 44), ( 45) and ( 46), and again the convexity of h, we obtain

δ τ ˆ1 0 h(v(w)) + ˆ1 0 (|∂ y w| 2 + |w| 2 ) + 1 e 2 k ˆ1 0 ∂ y w • (B(w)∂ y w) (47) ≤ δ τ ˆ1 0 h(v(w k-1 )) + δ e k e k (f k -v(w)(1)) • Dh(v(w)(1)) + h(v(w)(1)) - ˆ1 0 h(v(w)) = δ τ ˆ1 0 h(v(w k-1 )) + e k e k h(f k ) - ˆ1 0 h(v(w)) .
This inequality implies that

w 2 H 1 ((0,1);R n ) ≤ 2 τ + C Φ L ∞ (0,T ) h L ∞ (D) ,
for some constant C > 0 independent of , τ of Φ.

All the assumptions of the Leray-Schauder fixed-point theorem are thus satisfied. This yields the existence of a fixed-point solution w k ∈ H 1 ((0, 1); R n ) to S(1, w k ) = w k . Besides, using (47) with δ = 1, we have the discrete entropy inequality (28).

Proof of Proposition 1

Let us define by

V := n i=0 φ i ∈ R * + , ϕ := (φ 1 , • • • , φ n ) T and f := ϕ V .
From (T1), the vector f := f i 1≤i≤n obviously belongs to the set D.

If h defined by ( 10) is an entropy density for which A satisfies assumptions (H1)-(H2)-(H3), then A satisfies the same assumptions with the entropy density

h : D → R u → n i=1 u i log u i f i + (1 -ρ u ) log 1-ρu 1-ρ f .
Indeed, for all u ∈ D, Dh(u) = Dh(u) + g, where g := log

1-ρ f f i 1≤i≤n
is a constant vector in R n and D 2 h(u) = D 2 h(u). Moreover, the entropy density h has the following interesting property: f is the unique minimizer of h on D so that h(u) ≥ h(f ) = 0 for all u ∈ D. In the rest of the proof, for all w ∈ R n , we will denote by v(w) = (v i (w)) 1≤i≤n := (Dh) -1 (w) = Dh -1 (wg).

Let w ,k k∈N be a sequence of solutions to the regularized time-discrete problems (27) defined in Lemma 2 with the constant fluxes (φ 0 , • • • , φ n ) and the entropy density h. The entropy inequality (28) then reads

1 τ ˆ1 0 h(v(w ,k )) + ˆ1 0 (|∂ y w ,k | 2 + |w ,k | 2 ) + 1 e 2 k ˆ1 0 ∂ y w ,k • B(w ,k )∂ y w ,k (48)
≤ 1 τ ˆ1 0 h(v(w ,k-1 )) + e k e k h(f ) - ˆ1 0 h(v(w ,k )) .
In our particular case, for all k ∈ N, e k = V , e k = e 0 + V kτ and h(f ) = 0, so that we obtain

e 0 + V (k + 1)τ τ ˆ1 0 h(v(w ,k )) - e 0 + V kτ τ ˆ1 0 h(v(w ,k-1 )) ≤ 0.
This implies that for all k ∈ N and > 0,

(e 0 + V (k + 1)τ ) ˆ1 0 h(v(w ,k )) ≤ (e 0 + V τ ) ˆ1 0 h(v(w 0 )). ( 49
)
Let us denote by w ( ,τ ) : R * + → H 1 ((0, 1); R n ) the piecewise constant in time function defined by for a.a. y ∈ (0, 1), w ( ,τ ) (t, y) = w ,k (y) if (k -1)τ < t ≤ kτ.

Let T > 0 and ξ ∈ L 1 (0, T ) such that ξ ≥ 0 a.e. in (0, T ). Inequality (49) and Fubini's theorem for integrable functions implies that ˆT 0 ˆ1 0 (e 0 + V (k + 1)τ )h(v(w ( ,τ ) )) -(e 0 + V τ )h(v(w 0 )) ξ(t) dy dt ≤ 0.

From the proof of Theorem 2, we know that up to the extraction of a subsequence which is not relabeled, v(w ( ,τ ) ) ,τ >0 converges strongly in L 2 loc (R * + ; L 2 ((0, 1); R n )) and a.e. in R * + × (0, 1) as and τ go to zero to a global weak solution v to [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF]. Using Lebesgue dominated convergence theorem, and passing to the limit , τ → 0 in the above inequality yields ˆT 0 ˆ1 0 (e 0 + V t)h(v) -e 0 h(v(w 0 )) ξ(t) dy dt ≤ 0, which implies that there exists C > 0 such that for almost all t > 0,

(e 0 + V t) ˆ1 0 h(v) ≤ C, (50)
which yields inequality [START_REF] Lepoutre | Global well-posedness of a conservative relaxed cross diffusion system[END_REF]. In the rest of the proof, C will denote an arbitrary positive constant independent on the time t > 0. Furthermore, since v ∈ H 1 ((0, T ); (H 1 ((0, 1); R n )) )∩ L 2 ((0, T ); H 1 ((0, 1); R n )), it holds that v ∈ C 0 ((0, T ); L 2 ((0, 1); R n )) from [START_REF] Louis | Non-homogeneous boundary value problems and applications[END_REF], and the Lebesgue dominated convergence theorem implies that t ∈ R * + → ´1 0 h(v(t, y)) dy is a continuous function. Inequality (50) then holds for all t > 0.

For all 0 ≤ i ≤ n, let us denote by v i (t) := ´1 0 v i (t, y) dy. By convention, we define v 0 (t, y) := 1 -ρ v(t,y) and f 0 := 1 -ρ f . It can be checked from the weak formulation of (27) that ˆ1 0 v i w ,k = kφ i τ + e 0 ´1 0 v 0 i e 0 + V (k + 1)τ .

Passing to the limit , τ → 0 using the Lebesgue dominated convergence theorem, we obtain that for almost all t > 0,

v i (t) = e 0 ´1 0 v 0 i (y) dy + tφ i e 0 + V t , so that |v i (t) -f i | ≤ C e 0 + V t
. The continuity of v i implies that this equality holds for all t > 0.

The Csizàr-Kullback inequality states that for all t > 0,

v i (t, •) -v i (t) 2 L 1 (0,1) ≤ 2 ˆ1 0 v i (t, y) log v i (t, y) v i (t) dy = 2 ˆ1 0 v i (t, y) log v i (t, y) f i dy+2 ˆ1 0 v i (t, y) log f i v i (t)
dy.

Thus,

n i=0 v i (t, •) -f i L 1 (0,1) ≤ n i=0 v i (t, •) -v i (t) L 1 (0,1) + |f i -v i (t)| ≤ 2 ˆ1 0 h(v) + n i=0 2 log v i (t) f i + |f i -v i (t)| ≤ C e 0 + V t .
Hence inequality ( 22) and the result.

Proof of Proposition 2

Let (Φ m ) m∈N ⊂ Ξ be a minimizing sequence for J i.e such that

lim m→+∞ J (Φ m ) = inf Φ∈Ξ J (Φ).
By definition of the set Ξ, the sequence (Φ m ) m∈N is bounded in L ∞ (0, T ). Thus, up to a non relabeled extraction, it weakly-* converges to some limit Φ * ∈ Ξ in L ∞ (0, T ).

As a consequence, d dt e Φ m m∈N (respectively (e Φ m ) m∈N ) converges weakly-* (respectively strongly) in L ∞ (0, T ) to d dt e Φ * (respectively e Φ * ). For each m ∈ N, let v Φ m be the unique global weak solution to [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] associated to the fluxes Φ m . Its uniqueness is a consequence of assumption (C1). From the bounds obtained in the proof of Theorem 2 and the boundedness of (Φ m ) m∈N in L ∞ (0, T ), it holds that the sequences ∂ t v Φ m L 2 ((0,T );(H 1 (0,1)) , A(v Φ m )∂ y v Φ m L 2 ((0,T );L 2 (0,1)) and ∂ y v Φ m L 2 ((0,T );L 2 (0,1)) are also uniformly bounded in m.

Thus, up to the extraction of a subsequence which is not relabeled, using the compact injection of L 2 ((0, T ); H 1 ((0, 1); R n ))∩H 1 ((0, T ); (H 1 ((0, 1); R n )) ) into C((0, T ); L 2 ((0, 1); R n )) (see [START_REF] Louis | Non-homogeneous boundary value problems and applications[END_REF]), there exists v * ∈ L 2 ((0, T ); H 1 ((0, 1); R n )) ∩ H 1 ((0, T ); (H 1 ((0, 1); R n )) ) and

V * ∈ L 2 ((0, T ); L 2 ((0, 1); R n )) so that v Φ m v * weakly in L 2 ((0, T ); H 1 ((0, 1); R n )) ∩ H 1 ((0, T ); (H 1 ((0, 1); R n )) ), v Φ m -→ v *
strongly in C((0, T ); L 2 ((0, 1); R n )) and a.e. in (0, T ) × (0, 1),

A(v Φ m )∂ y v Φ m V * weakly in L 2 ((0, T ); L 2 ((0, 1); R n )).
Using similar arguments as in the proof of Theorem 2, we also obtain that V * is necessarily equal to A(v * )∂ y v * . Passing to the limit m → +∞, we obtain that for all ψ ∈ L 2 ((0, T ); H 1 ((0, 1); R n )),

ˆT 0 ˆ1 0 ∂ t v * • ψ dt dy + ˆT 0 ˆ1 0 1 e Φ * (t) 2 ∂ y ψ • (A(v * )∂ y v * ) dt dy + ˆT 0 d dt e Φ * (t) e Φ * (t) ˆ1 0 (v * • ψ + yv * • ∂ y ψ) dt dy = ˆT 0 1 e Φ * (t) ϕ * (t) • ψ(1) dt.
Assumption (C1) yields v * = v Φ * . The above convergence results then imply that

J (Φ m ) -→ m→+∞ J (Φ * ),
and hence Φ * is a minimizer of problem [START_REF] Donald | Handbook of physical vapor deposition (PVD) processing[END_REF]. Hence the result.

Numerical tests

In this section, we present some numerical tests illustrating the results of Section 2 on the prototypical example of Section 1.1. In Section 4.1, we present the numerical scheme used in our simulations to compute an approximation of a solution of [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF]. In Section 4.2 and Section 4.3, some numerical tests which illustrate Proposition 1 and Proposition 2 are detailed.

Discretization scheme

In view of the optimization problem [START_REF] Donald | Handbook of physical vapor deposition (PVD) processing[END_REF] we are aiming at, it appears that a fully implicit unconditionally stable scheme is needed to allow the use of reasonably large time steps.

We present here the numerical scheme used for the discretization of ( 19), for the particular model presented in Section 1.1. We do not provide a rigorous numerical analysis for this scheme here.

Let M ∈ N * and ∆t := T M . We define for all 0 ≤ m ≤ M , t m := m∆t. The discrete external fluxes are characterized for every 0 ≤ i ≤ n by vectors φ In addition, let Q ∈ N * and ∆y := 1 Q and y q := (q -0.5)∆y. For all 0 ≤ i ≤ n, 1 ≤ q ≤ Q and 0 ≤ m ≤ M , we denote by v m,q i the finite difference approximation of v i at time t m and point y q ∈ (0, 1). Here again, we use the convention that v 0 = 1 -ρ v .

i := φ m i 1≤m≤M ∈ R M + , where φ m i = ´tm t m-1 φ i (s) ds. For every 1 ≤ m ≤ M ,
We use a centered second-order finite difference scheme for the diffusive part of the equation, and a first-order upwind scheme for the advection part, together with a fully implicit time scheme. Assuming that the approximation v m-1,q i 0≤i≤n,1≤q≤Q

is known, one computes ( v m,q i ) 0≤i≤n, 1≤q≤Q as solutions of the following sets of equations. For all 0 ≤ i ≤ n and 2

≤ q ≤ Q -1, v m,q i -v m-1,q i ∆t = e d m e m y q v m,q+1 i -v m,q i ∆y (51)
+ 0≤j =i≤n K ij e 2 m v m,q j v m,q+1 i + v m,q-1 i -2 v m,q j 2∆y 2 -v m,q i v m,q+1 j + v m,q-1 j -2 v m,q j 2∆y 2
together with boundary conditions which reads for all 0 ≤ i ≤ n,

0≤j =i≤n K ij e m v m,1 j v m,2 i -v m,1 i ∆y -v m,1 i v m,2 j -v m,1 j ∆y = 0, ( 52
) 0≤j =i≤n K ij e m v m,Q j v m,Q-1 i -v m,Q i ∆y -v m,Q i v m,Q-1 j -v m,Q j ∆y = -e d m v m,Q i + φ m i . (53)
The nonlinear system of equations ( 51)-( 52)-( 53), whose unknowns are ( v m,q i ) 0≤i≤n,1≤q≤Q is solved using Newton iterations with initial guess v m-1,q i 0≤i≤n,1≤q≤Q

. The obtained solution does not satisfy in general the desired non-negativeness and volumic constraints. This is the reason why an additional projection step is performed. For all 0 ≤ i ≤ n and

1 ≤ q ≤ Q, we define v m,q i := [ v m,q i ] + n j=0 [ v m,q j ] + , so that v m,q i ≥ 0 and n j=0 v m,q j = 1.
We numerically observe that this scheme is unconditionally stable with respect to the choice of discretization parameters ∆t and ∆y. In itial con centration s

v M 0 v M 1 v M 2 v M 3 3 i=0 v M i (

Optimization of the fluxes

The optimization problem ( 24) is solved in practice using an adjoint formulation associated to the discretization scheme described in Section 4.1. We refer the reader to [START_REF] Bakhta | Mathematical models and numerical simulation of photovoltaic devices[END_REF] for more details and comparisons between our model and experimental results obtained in the context of thin film CIGS solar cell fabrication. To illustrate Proposition 2, we proceed as follows: first, we perform a simulation of [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] with external fluxes Φ sim for a duration T to obtain a final thickness e Φ sim (T ) and final concentrations v Φ sim (T, •), then, we solve the minimization problem [START_REF] Donald | Handbook of physical vapor deposition (PVD) processing[END_REF] to obtain optimal fluxes Φ * where the target concentrations are

v opt (y) = v Φ sim (T, y) ∀ y ∈ (0, 1)
and the target thickness is e opt = e Φ sim (T ).

Lastly, we perform another simulation of [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] with the obtained optimal fluxes Φ * and compare the final concentrations v Φ * and the final thickness e Φ * to the target concentrations v opt and the target thickness e opt . In Figures 4-( We use a quasi-Newton iterative gradient algorithm for the resolution of the problem. At each iteration of the algorithm, the approximate hessian is updated by means of a BFGS procedure and the optimal step size is the solution of a line search subproblem. More details on the numerical optimization algorithms can be found in [START_REF] Lemaréchal | Jean Charles Gilbert and Claudia Sagastizàbal[END_REF]. The initial guess Φ 0 is taken of the form (56) where β i = 1 for all 0 ≤ i ≤ n.

The algorithm is run until one of the following stopping criterion is reached : either (J (Φ) ≤ ε) or ( ∇ Φ J (Φ) L 2 ≤ ν) with ε = 10 -5 and ν = 10 -5 .

In Figure 8-(a) we plot the evolution of the value of the cost J (Φ) with respect to the number of iterations. We refer the reader to [START_REF] Bakhta | Mathematical models and numerical simulation of photovoltaic devices[END_REF] for more details and comparison between different minimization approaches.

We numerically observe that all the concentrations are well reconstructed and that the value of the optimal thickness e Φ * = 483.4022 is close to the target thickness e Φ sim = 483.4. Unlike the external fluxes Φ sim , the optimal fluxesΦ * are not piecewise constant. These tests show that the uniqueness of a solution to the optimization problem [START_REF] Donald | Handbook of physical vapor deposition (PVD) processing[END_REF] can not be expected in general. We refer the reader to [START_REF] Bakhta | Mathematical models and numerical simulation of photovoltaic devices[END_REF] for results obtained in the case of the control of PVD process for the fabrication of this film solar cells.

Conclusion

In this work, we propose and analyze a one-dimensional model for the description of a PVD process. The evolution of the local concentrations of the different chemical species in the bulk of the growing layer is described via a system of cross-diffusion equations similar to the ones studied in [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF][START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF]. The growth of the thickness of the layer is related to the external fluxes of atoms that are absorbed at the surface of the film.

The existence of a global weak solution to the final system using the boundedness by entropy method under assumptions on the diffusion matrix of the system close to those needed in [START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] is established. In addition, the entropy density h is required to be continuous (hence bounded) on the set D = u = (u i ) 1≤i≤n ∈ R n + , n i=1 u i ≤ 1 . We prove the existence of a solution to an optimization problem under the assumption [START_REF] Donald | Handbook of physical vapor deposition (PVD) processing[END_REF]. that there exists a unique global weak solution to the obtained system, whatever the value of the external fluxes.

Lastly, in the case when the entropy density is defined by h(u) = n i=1 u i log u i + (1 -ρ u ) log(1 -ρ u ), we prove in addition that, when the external fluxes are constant and positive, the local concentrations converge in the long time to a constant profile at a rate which scales like O 1 t . A discretization scheme, which is observed to be unconditionnaly stable, is introduced for the discretization of [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF]. This scheme enables to preserve constraints (5) at the discretized level.

We see this work as a preliminary step before tackling related problems in higher dimension, including surfacic diffusion effects. Besides, the proof of assumption (C1) remains an open question in general at least to our knowledge. Lastly, a nice theoretical question which is not tackled in this paper, but will be the object of future research, is the characterization of the set of reachable concentration profiles.

Figure 1 : 2 2

 12 Figure 1: Illustration of the composition of the film layer at time t in the case n = 2

  and by f := (f i ) 1≤i≤n ∈ D. Let us also denote by

  we denote by e Φ : t ∈ [0, T ] → e 0 + ´t 0 n i=0 φ i (s) ds the time-dependent thickness of the film, and by v Φ a solution to (19) associated with the external fluxes Φ.

ϕ

  p := (φ 1,p , • • • , φ n,p ) T, and e p (t) := e 0 + ˆt 0 n i=0 φ i,p (s) ds.

  |e(s) -e p (s)| ≤ max (|e(s j+1 ) -e p (s j )|, |e p (s j+1 ) -e(s j )|) ≤ max (|e(s j+1 ) -e(s j )| + |e(s j ) -e p (s j )|, |e p (s j+1 ) -e(s j+1 )| + |e(s j+1 ) -e(s j )|) ≤ ε.

  the thickness of the thin film and it derivative at time t m are approximated respectively by e m := e 0 + ∆t ≈ e(t m ), and e d m := n i=0 φ m i ≈ e (t m ).

Figure 3 :

 3 Figure 3: Long-time behavior in the case of non negative constant external fluxes.

• 7 •

 7 a), 5-(a), 6-(a) and 7-(a) are plotted the final concentration profiles v Φ sim (T, •) resulting from the simulation of (19) with the following parameters : • T = 120, M = 120, Q = 100, ∆t = 1, ∆y = 0.01, e 0 = 1. • Cross-diffusion coefficients K ij j = 0 j = 1 j = 2 j = 3 i = External fluxes Φ sim of the form (54) with i = 0 i = 1 i = 2 i Initial concentrations v 0 i of the form (55) with w 0 (y) = y, w 1 (y) = 2y, w 2 (y) = √ y and w 3 (y) = 0.

Figure 4 :

 4 Figure 4: Reconstruction of the final concentration of the species i = 0.

Figure 5 :

 5 Figure 5: Reconstruction of the final concentration of the species i = 1.

Figure 6 :

 6 Figure 6: Reconstruction of the final concentration of the species i = 2.

Figure 7 :

 7 Figure 7: Reconstruction of the final concentration of the species i = 3.

  (log scale) Cost fucntion (log scale)

Figure 8 :

 8 Figure 8: Convergence of the BFGS gradient descent algorithm for the minimization problem (24).

Acknowledgements

We are very grateful to Eric Cancès and Tony Lelièvre for very helpful advice and discussions. We would like to thank Jean-François Guillemoles, Marie Jubeault, Torben Klinkert and Sana Laribi from IRDEP, who introduced us to the problem of modeling PVD processes for the fabrication of thin film solar cells. The EDF company is acknowledged for

A standard practice in the production of thin film CIGS (Copper, Indium, Gallium, Selenium) solar cells by means of PVD process is to consider piecewise-constant external fluxes. We refer the reader to [START_REF] Bakhta | Mathematical models and numerical simulation of photovoltaic devices[END_REF] for further details. In the following numerical tests, we consider time-dependent functions of the form

where 0 < τ i 1 < τ i 2 < T and (α i 1 , α i 2 , α i 3 ) ∈ (R + ) 3 are non-negative constants for all 0 ≤ i ≤ n. Besides, we consider initial condition of the form

where w i : [0, 1] → R + are functions which will be precised below. In the whole section, system ( 19) is simulated with four species (i.e. n = 3).

In Figure 2 are plotted the results obtained for the simulation of ( 19) with the following parameters :

Long-time behaviour results

In this section is given a numerical illustration of Proposition 1. We consider timedependent functions of the form

where

In Figure 3 are plotted the results obtained for the the simulation of [START_REF] Olga | Nn ural ceva, linear and quasilinear equations of parabolic type, translated from the russian by s. smith[END_REF] with the following parameters : For all 0 ≤ i ≤ n, let vi := β i / n j=0 β j . We consider the time-dependent quantity

where the relative entropy h is defined in [START_REF] Lepoutre | Global well-posedness of a conservative relaxed cross diffusion system[END_REF]. We also consider the quantities

In Figure 3 The evolution of (η i (t)) 0≤i≤n (respectively η(t) and γ(t)) with respect to t is shown in Figure 3-(c) (respectively 3-(d) and 3-(e)). We numerically observe that these quantities are affine functions of t in the asymptotic regime which illustrates the theoretical result of Proposition 1.

funding. We would also like to thank Daniel Matthes for very helpful discussions on the theoretical part and the reviewers whose comments helped in very significantly improving the quality of the paper.

Appendix

Formal derivation of the diffusion model (3)

We present in this section a simplified formal derivation of the cross-diffusion model ( 3) from a one-dimensional microscopic lattice hopping model with size exclusion, in the same spirit than the one proposed in [START_REF] Burger | Nonlinear cross-diffusion with size exclusion[END_REF].

We consider here a solid occupying the whole space R and discretize the domain using a uniform grid of step size ∆x > 0. At any time t ∈ [0, T ], we denote by u k,t i the number of atoms of type i (0

Let ∆t > 0 denote a small enough time step. We assume that during the time interval ∆t, an atom i located in the k th interval can exchange its position with an atom of type j (j = i) located in one of the two neighbouring intervals with probability p ij = p ji > 0. In average, we obtain the following evolution equation for u k,t i :

This yields that

Choosing ∆t and ∆x so that these quantities satisfy a classical diffusion scaling 2∆x 2 ∆t = α > 0, denoting by K ij := αp ij and letting the time step and grid size go to 0, we formally obtain the following equation for the evolution of u i on the continuous level:

which is identical to the system of equations (3) introduced in the first section. Of course, this formal argument can be easily extended to any arbitrary dimension.

Leray-Schauder fixed-point theorem

Theorem 3 (Leray-Schauder fixed-point theorem). Let B be a Banach space and S : B × [0, 1] → B be a continuous map such that (A1) S(x, 0) = 0 for each x ∈ B;

(A2) S is a compact map;

(A3) there exists a constant M > 0 such that for each pair (x, σ) ∈ B ×[0, 1] which satisfies x = S(x, σ), we have x < M .

Then, there exists a fixed-point y ∈ B satisfying y = S(y, 1).