C L Ément

Canc Ès

AND Daniel Matthes

Flore Nabet

A TWO-PHASE TWO-FLUXES DEGENERATE CAHN-HILLIARD MODEL AS CONSTRAINED WASSERSTEIN GRADIENT FLOW

Keywords: 1991 Mathematics Subject Classification. 35K65, 35K41, 49J40, 76T99 Multiphase flow, Cahn-Hilliard type system, constrained Wasserstein gradient flow

We study a non-local version of the Cahn-Hilliard dynamics for phase separation in a two-component incompressible and immiscible mixture with linear mobilities. In difference to the celebrated local model with nonlinear mobility, it is only assumed that the divergences of the two fluxes -but not necessarily the fluxes themselves -annihilate each other. Our main result is a rigorous proof of existence of weak solutions. The starting point is the formal representation of the dynamics as a constrained gradient flow in the Wasserstein metric. We then show that time-discrete approximations by means of the incremental minimizing movement scheme converge to a weak solution in the limit. Further, we compare the non-local model to the classical Cahn-Hilliard model in numerical experiments. Our results illustrate the significant speed-up in the decay of the free energy due to the higher degree of freedom for the velocity fields.

1. Motivation and presentation of the model 1.1. Introduction. We are interested in a non-local Cahn-Hilliard system

∂ t c i -∇ • (m i c i ∇ (µ i + Ψ i )) = m i θ i ∆c i for i ∈ {1, 2}, (1a)
c 1 + c 2 = 1, (1b) µ 1 -µ 2 = -α∆c 1 + χ(1 -2c 1 ), (1c)
modelling the flow of an incompressible and immiscible mixture in a bounded convex domain Ω. Above, m 1 and m 2 are positive mobility constants, χ > 0 is a parameter for the chemical activity, and the non-negative numbers θ 1 and θ 2 control the degree of thermal agitation. The system is complemented with initial conditions [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] (c i ) |t=0 = c 0 i ∈ H 1 (Ω) with c 0 i ≥ 0 and c 0 1 + c 0 2 = 1 in Ω, and boundary conditions [START_REF] Barrett | Finite element approximation of the Cahn-Hilliard equation with degenerate mobility[END_REF]

c i ∇(µ i + Ψ i ) • n = ∇c i • n = 0 on ∂Ω × (0, ∞).
Equation (1a) can be rewritten as a conservation law (4)

∂ t c i + ∇ • J i = 0,
where

J i = -m i c i ∇ (µ i + Ψ i + θ i log(c i ))
. Summing (4) over i ∈ {1, 2} and using (1b) yields ∇ • J tot = 0 where J tot = J 1 + J 2 denotes the total flux. In particular, we do not require that J tot = 0 as for the classical (or local) degenerate Cahn-Hilliard model [START_REF] De Gennes | Dynamics of fluctuations and spinodal decomposition in polymer blends[END_REF][START_REF] Gurtin | Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance[END_REF][START_REF] Elliott | On the Cahn-Hilliard equation with degenerate mobility[END_REF]. The relaxation of the constraint from vanishing total flux to divergence free total flux was initially proposed by E and Palffy-Muhoray in [START_REF] Palffy-Muhoray | Phase separation in incompressible systems[END_REF] and studied formally by Otto and E in [START_REF] Otto | Thermodynamically driven incompressible fluid mixtures[END_REF]. It was in particular noticed in [START_REF] Otto | Thermodynamically driven incompressible fluid mixtures[END_REF] that the system (1) can be interpreted as the constrained Wasserstein gradient flow of some Ginzburg-Landau energy where the velocity field V = (v 1 , v 2 ) transporting the concentrations c = (c 1 , c 2 ) has to preserve the constraint (1b).

The nonlocal model is derived formally in Section 1.2, then compared to the classical (or local) degenerate Cahn-Hilliard model in Section 1.3. Numerical illustrations of its behavior are given in Section 1. [START_REF] Barrett | On fully practical finite element approximations of degenerate Cahn-Hilliard systems[END_REF]. In Section 2, we introduce the necessary material to prove our main result, that is the global existence of a weak solution to the model ( 1)- [START_REF] Barrett | Finite element approximation of the Cahn-Hilliard equation with degenerate mobility[END_REF]. This existence result is obtained by showing the convergence of a minimizing movement scheme à la Jordan, Kinderlehrer, and Otto [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF]. Section 3 is devoted to the proof of the convergence of the minimizing movement scheme. 1.2. Derivation of the model. We consider an incompressible mixture composed of two phases flowing within a bounded open convex subset Ω of R d with d ≤ 3. The fluid is incompressible, so its composition at time t ≥ 0 is fully described by the saturations c i (x, t) ∈ [0, 1], i ∈ {1, 2}, i.e., the volume ratio of the phase i in the fluid. This leads to the constraint [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] c 1 + c 2 = 1 in Ω × (0, ∞).

We assume that this constraint is already satisfied at the initial time t = 0, where the composition of the mixture is given by c 0 = (c 0 1 , c 0 2 ) : Ω → [0, 1] 2 , i.e. ( 6)

c 0 1 + c 0 2 = 1 in Ω.
We further assume that both phases have positive mass. The motions of each phase is governed by a linear transport equation ( 7)

∂ t c i + ∇ • (c i v i ) = 0 in Ω × (0, ∞),
where v i : Ω × R + → R d denotes the speed of the phase i, and V = (v 1 , v 2 ). To enforce mass conservation, the boundary ∂Ω of Ω is assumed to be impervious, hence

c 1 v 1 • n = c 2 v 2 • n = 0 on ∂Ω, (8)
where n denotes the outward normal to ∂Ω, so that Gauss' theorem implies [START_REF] Cahn | The Cahn-Hilliard equation with a concentration dependent mobility: motion by minus the Laplacian of the mean curvature[END_REF])

Ω c i (x, t)dx = Ω c 0 i (x)dx = |Ω|c i , t ≥ 0, i ∈ {1, 2}.
Consequently, at each time t ≥ 0, the saturations c(t) belong to the set A = A 1 × A 2 , where

A i = c i : Ω → [0, 1] measurable Ω c i (x)dx = |Ω|c i , i ∈ {1, 2}.
To each configuration c ∈ A, we associate the energy

E(c) = E Dir (c) + E chem (c) + E therm (c) + E cons (c) + E ext (c),
whose components are as follows.

• Dirichlet energy: for given α 1 , α 2 > 0,

E Dir (c) :=      i∈{1,2} α i 2 Ω |∇c i | 2 dx if c ∈ H 1 (Ω) 2 , +∞ otherwise,
penalizes the spatial variations of the saturation profiles. • Chemical energy: for a given χ > 0,

E chem (c) = χ Ω c 1 c 2 dx
measures the impurity of the mixture.

• Thermal energy: for given θ 1 , θ 2 ≥ 0,

E therm (c) = i∈{1,2} θ i Ω H(c i )dx, with H(c) = c log(c) -c + 1 ≥ 0,
is the free thermal energy. The case θ 1 = θ 2 = 0 is called the deep quench limit [START_REF] Cahn | The Cahn-Hilliard equation with a concentration dependent mobility: motion by minus the Laplacian of the mean curvature[END_REF].

• Constraint:

E cons (c) = Ω E cons (c)dx, where E cons (c) = 0 if c 1 + c 2 = 1, +∞ otherwise, realizes the constraint (5). • Exterior potential: for given external potentials Ψ 1 , Ψ 2 ∈ H 1 (Ω), E ext (c) = i∈{1,2} Ω c i (x)Ψ i (x)dx
is the potential energy related to exterior forces like gravity or electrostatic forces. All of these components are convex in c except for the chemical energy E chem , which, however, is smooth. Denote by

X = c ∈ H 1 (Ω; [0, 1]) 2 c 1 + c 2 = 1 a.e.
in Ω the domain of E, i.e., [START_REF] Cancès | Simulation of multiphase porous media flows with minimizing movement and finite volume schemes[END_REF]

E(c) < ∞ ⇔ c ∈ X .
Before entering into the rigorous derivation of the PDEs that govern the gradient flow of the energy E(c) with respect to a tensorized Wasserstein distance, we provide formal calculations based on the framework of generalized gradient flows of [START_REF] Mielke | A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems[END_REF][START_REF] Peletier | Variational modelling: Energies, gradient flows, and large deviations[END_REF][START_REF] Cancès | The gradient flow structure of immiscible incompressible two-phase flows in porous media[END_REF] in order to identify the underlying PDEs.

The motion of the phases induces the viscous dissipation

(11) D(c, V ) = i∈{1,2} 1 2m i Ω c i |v i | 2 dx, ∀c ∈ A, ∀V ∈ Z(c),
where m i is the mobility coefficient of the ith phase, and where

Z(c) = V = (v 1 , v 2 ) : Ω → R d 2 c i v i • n = 0 on ∂Ω
denotes the space of the admissible (but unconstrained) vector fields. This quantity is closely related to the tensorized Wasserstein distance through Benamou-Brenier formula [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF], see [START_REF] Junge | A fully discrete variational scheme for solving nonlinear Fokker-Planck equations in multiple space dimensions[END_REF] later on. We suppose as in [START_REF] Cancès | The gradient flow structure of immiscible incompressible two-phase flows in porous media[END_REF] that at each time t ≥ 0, the phase speeds V = (v 1 , v 2 ) is selected by the following steepest descent condition:

(12) V ∈ argmin V =( v1, v2)∈Z(c)   D(c, V ) + sup w∈∂E(c) i∈{1,2} Ω c i v i • ∇w i dx   , where E's subdifferential ∂E(c) = w ∈ L 2 (Ω; R d ) 2 E( c) -E(c) - Ω w • ( c -c)dx ≥ o c -c L 2 (Ω)
is non-empty at each c ∈ X with c 1 ∈ H 2 (Ω) and f (c 1 ) ∈ L 2 (Ω), and there, its elements w = (w 1 , w 2 ) are characterized by (c.f. [START_REF] Cancès | The gradient flow structure of immiscible incompressible two-phase flows in porous media[END_REF])

(13) w 1 -w 2 = -α∆c 1 + χ(1 -2c 1 ) + f (c 1 ) + (Ψ 1 -Ψ 2 ).
In formula [START_REF] Cancès | Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations[END_REF], we have set α = α 1 + α 2 and ( 14)

f (c 1 ) = log c θ1 1 (1 -c 1 ) θ2 = θ 1 log(c 1 ) -θ 2 log(c 2 ).
The integral in [START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF] has to be understood as the evaluation of the element w = (w 1 , w 2 ) in E's subdifferential at c on the infinitesimal changes ∂ t c i = -∇ • (c i v i ) induced by V . Since condition [START_REF] Cancès | Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations[END_REF] only determines the difference w 1w 2 , but not the w i individually, is easily seen that the maximum above is +∞ unless ∇ • (c 1 v 1 + c 2 v 2 ) = 0 holds. That implies that any minimizer V is such that the constraint ( 5) is preserved.

In order to identify V , we swap minimization and maximization in the variational problem [START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF]. The inner minimization is a quadratic problem in V , with solution ( 15)

v i = -m i ∇w i on {c i > 0} for i ∈ {1, 2}.
The outer maximization in then a quadratic problem for w ∈ ∂E(c), which amounts to

-∇ •   i∈{1,2} m i c i ∇w i   = 0. ( 16
)
Together with the condition [START_REF] Cancès | Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations[END_REF], this elliptic equation determines w uniquely.

To sum up, the system of partial differential equations implied on the solution to the variational problem [START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF] is

∂ t c i -∇ • (m i c i ∇w i ) = 0 for i ∈ {1, 2}, (17a) c 1 + c 2 = 1, (17b) w 1 -w 2 = -α∆c 1 + χ(1 -2c 1 ) + f (c 1 ) + (Ψ 1 -Ψ 2 ), (17c)
to be satisfied in Ω × (0, ∞). The system relation [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] is complemented by homogeneous Neumann boundary conditions [START_REF] Dolbeault | A new class of transport distances between measures[END_REF] -

∇c 1 • n = 0 on ∂Ω × R +
and the no-flux conditions [START_REF] Blowey | The Cahn-Hilliard gradient theory for phase separation with nonsmooth free energy. II. Numerical analysis[END_REF]. Notice that with these boundary conditions, the set (17) of equations is (formally) closed in the sense that if c is known at some instance of time, then the w i are uniquely determined (up to irrelevant global additive constants) by (17c) and by the elliptic equation ( 16) that follows from adding (17a) for i = 1, 2 in combination with the conservation

∂ t (c 1 + c 2 ) = 0 implied by (17b
µ i = w i -θ i log(c i ) -Ψ i , i ∈ {1, 2}, so that (17c) turns to (20) µ 1 -µ 2 = -α∆c 1 + χ(1 -2c 1 ).
So far, µ= (µ 1 , µ 2 ) is only defined up to an additive constant. This degree of freedom is eliminated by imposing for almost all t > 0 that ( 21)

Ω µ(x, t)dx = 0, where µ = c 1 µ 1 + c 2 µ 2 .
Whilst the phase chemical potential µ i have a clear physical sense only on {c i > 0}, the global chemical potential µ remains meaningful in the whole Ω. In particular, its spatial variations can be controlled, and thus µ itself too thanks to a Poincaré-Wirtinger inequality.

With the help of the chemical potentials, the system (17) turns into (1), to be complemented with boundary conditions (3) and initial conditions (2).

In the next section, we highlight some differences between the non-local model (1) and the local degenerate Cahn-Hilliard model that has been studied for instance in [START_REF] Elliott | On the Cahn-Hilliard equation with degenerate mobility[END_REF].

1.3. Comparison with the classical degenerate Cahn-Hilliard model. Even in the simple situation where the external potentials Ψ i are equal to 0 and where θ 1 = θ 2 = 0, the system (17) differs as soon as d ≥ 2 from the local degenerate Cahn-Hilliard model that can be written as ( 22)

∂ t c -∇ • (η(c)∇µ) = 0, µ = -α∆c + χ(1 -2c), η(c) = m 1 m 2 c(1 -c) m 1 c + m 2 (1 -c) .
We refer to [START_REF] Elliott | On the Cahn-Hilliard equation with degenerate mobility[END_REF] for the existence of weak solutions to [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] (complemented with suitable boundary conditions) and to [START_REF] Barrett | On fully practical finite element approximations of degenerate Cahn-Hilliard systems[END_REF] for the extension of the model to the case of N phases (N ≥ 3). Here, µ is the generalized chemical potential that is defined as the difference of the phase chemical potentials.

The energy E associated to [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] is similar to the one of our problem, i.e.,

E(c) = α 2 Ω |∇c| 2 dx + χ Ω c(1 -c)dx.
But both the equation governing the motion [START_REF] Blanchet | Convergence of the mass-transport steepest descent scheme for the subcritical Patlak-Keller-Segel model[END_REF] and the dissipation [START_REF] Cancès | The gradient flow structure of immiscible incompressible two-phase flows in porous media[END_REF] have to be modified. More precisely, the continuity equation ( 7) must be replaced by its nonlinear counterpart

∂ t c + ∇ • (η(c)v) = 0,
while the dissipation is now given by

D(c, v) = Ω η(c)|v| 2 dx.
Therefore, the PDEs [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] can still be interpreted as the gradient flow of the energy E, but the geometry is different: rather than considering some classical quadratic Wasserstein distance for each phase and to constrain the sum of the concentrations to be equal to 1 (as it will be the case for our approach), the set

A = c ∈ L 1 (Ω; R + ) Ω c(x, t)dx = Ω c 0 1 (x)dx
has to be equipped with the weighted Wasserstein metric corresponding to the concave mobility η. We refer to [START_REF] Dolbeault | A new class of transport distances between measures[END_REF] for the description of the corresponding metric and to [START_REF] Lisini | Cahn-Hilliard and thin film equations with nonlinear mobility as gradient flows in weighted-Wasserstein metrics[END_REF] for the rigorous recovery of ( 22) by a gradient flow approach.

The difference between the non-local model [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] and the local one [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] can also be seen as follows. Summing the first equation of [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] for i ∈ {1, 2} yields [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF] ∇ • J tot = 0, where

J tot = i∈{1,2} c i v i = - i∈{1,2} m i c i ∇µ i .
The equation for c 1 can then be rewritten (24)

∂ t c 1 + ∇ • (ρ(c 1 )J tot -η(c 1 )∇(µ 1 -µ 2 )) = 0, with ρ(c) = m 1 c m 1 c + m 2 (1 -c) .
Thus our model [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] boils down to the local Cahn-Hilliard equation as soon as J tot ≡ 0. This is the case when d = 1 because of ( 23), but no longer if d ≥ 2. Since our non-local model does not impose that J 1 = -J 2 , it allows for additional motions. These motions -corresponding to the transport term ∇ • (ρ(c 1 )J tot ) in ( 24)-contribute to the dissipation as shows the formula

D(c, -∇µ) = Ω |J tot | 2 m 1 c 1 + m 2 c 2 dx + D(c 1 , -∇(µ 1 -µ 2 )).
Therefore, and as already noticed by Otto and E in [START_REF] Otto | Thermodynamically driven incompressible fluid mixtures[END_REF], the instantaneous dissipation corresponding to a phase configuration c is greater for the non-local model [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] than for the local model ( 22) and the energy decreases faster.

1.4. Numerical illustration. The goal of this section is to illustrate the behavior of the model [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] and to compare it with the classical degenerate Cahn-Hilliard problem [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. In order to solve numerically [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] we use an implicit in time finite volume scheme with upstream mobility described in [START_REF] Cancès | Finite volume approximation of a degenerate immiscible two-phase flow model of cahnhilliard type[END_REF] and inspired from the oil engineering context [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. The mesh is triangular and assumed to fulfill the so-called orthogonality condition [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusionconvection problem on a triangular mesh[END_REF][START_REF] Eymard | Finite volume methods[END_REF] (this amounts to requiring the mesh to be Delaunay) so that the diffusive fluxes can be approximated thanks to a two-point flux approximation in a consistent way. As it is exposed in [START_REF] Cancès | Finite volume approximation of a degenerate immiscible two-phase flow model of cahnhilliard type[END_REF], the scheme is positivity preserving (i.e., 0 ≤ c i,h ≤ 1), it is energy diminishing (the discrete counterpart of the energy is decreasing) and entropy stable. It leads to a nonlinear system of algebraic equations to be solved at each time step. It is shown in [START_REF] Cancès | Finite volume approximation of a degenerate immiscible two-phase flow model of cahnhilliard type[END_REF] that this system admits (at least) one solution that is computed thanks to the Newton-Raphson method.

Concerning the problem ( 22), we use a similar approach, but since the mobility function η is no longer monotone, we have to use an implicit Godunov scheme to discretize it as a generalization of the upstream mobility (see for instance [START_REF] Cancès | Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations[END_REF]). Here again, the discrete solution remains bounded between 0 and 1, the energy is decreasing and the entropy remains bounded. Here again, the resulting nonlinear system is solved at each time step by the mean of the Newton-Raphson method.

Remark 1.1. Alternative numerical methods have been proposed in order to solve degenerate Cahn-Hilliard problems. We won't perform here an exhaustive list, but let us mention the contributions of Barrett et al. based on conformal finite elements [START_REF] Barrett | Finite element approximation of the Cahn-Hilliard equation with degenerate mobility[END_REF][START_REF] Barrett | On fully practical finite element approximations of degenerate Cahn-Hilliard systems[END_REF]. Even though very efficient, these methods have the drawback of requiring a small stabilization to be tuned following the mesh size. Unless one considers non-smooth energies as in [START_REF] Blowey | The Cahn-Hilliard gradient theory for phase separation with nonsmooth free energy. II. Numerical analysis[END_REF], the scheme does not preserve the bounds 0 ≤ c i,h ≤ 1. These difficulties are overpassed in our approach by using some entropy stable hyperbolic fluxes to discretize the mobilities.

Since our model has a Wasserstein gradient flow structure, it would be natural to use a Lagrangian method as for instance in [START_REF] Blanchet | Convergence of the mass-transport steepest descent scheme for the subcritical Patlak-Keller-Segel model[END_REF][START_REF] Matthes | A convergent Lagrangian discretization for a nonlinear fourth-order equation[END_REF][START_REF] Junge | A fully discrete variational scheme for solving nonlinear Fokker-Planck equations in multiple space dimensions[END_REF][START_REF] Carrillo | A lagrangian scheme for the solution of nonlinear diffusion equations using moving simplex meshes[END_REF]. The main problem with this approach is that both phase move with their own speed, therefore such an approach would impose to move two meshes simultaneously. It is then rather unclear how to manage the constraint (5) in this case. For this reason, it seems more suitable to stick to an Eulerian description. An alternative approach to solve numerically our problem would therefore be to adapt the ALG2-JKO algorithm of Benamou et al. [START_REF] Benamou | An augmented Lagrangian approach to Wasserstein gradient flows and applications[END_REF] to our setting , see [START_REF] Cancès | Simulation of multiphase porous media flows with minimizing movement and finite volume schemes[END_REF].

We propose two different test cases that will allow to illustrate the difference between the local model ( 22) and the non-local model [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF]. For both of them, we do not consider any exterior potential, i.e., Ψ i = 0, and we neglect the thermal diffusion, i.e., θ i = 0. Both phase mobilities m i are assumed to be equal to 1.

1.4.1. Test case 1: from a cross to a circle. We start from an initial data that is the characteristic function of a cross and we choose α = (3.6).10 -4 and χ = 0.8. Since α χ, it follows from the Modica and Mortola' result [START_REF] Modica | Un esempio di γ-convergenza[END_REF] that the free energy is close to the perimeter of a characteristic set (up to a multiplicative constant). This means that (up to a small regularization) both the local and the non-local Cahn-Hilliard models aim at minimizing the perimeter of the sets {c 1 = 0} and {c 1 = 1} corresponding to pure phases. Since the non-local model allows for more movements (cf. Section 1.3), the energy (thus the perimeter) should decay faster for the nonlocal model. This is indeed what we observe on Figures 1 and2 Similarly to the local model, the non-local model is able to reproduce the spinodal decomposition for mixtures. In order to illustrate this fact, we start from an initial data which consists in a constant concentration plus a small random perturbation:

c 0 1 (x) = 0.5 + r(x), r 1.
Since c 1 = 0.5 is very unfavorable from an energetic point of view, both phase will separate very rapidly, letting areas with pure phase appear. Then these area will cluster in order to minimize their perimeter. We choose α = 3.10 -4 and χ = 0.96 and we plot on Figure 3 some snapshots to illustrate the spinodal decomposition corresponding to models [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF] and [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. On Figure 4, we compare the evolution of the energy along time for spinodal decomposition corresponding to both models. As expected, the energy decay is faster for the non-local model than for the local one. But contrarily to Test case 1, the solutions seem to converge towards different steady states. 17) and the local one [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] for the spinodal decomposition test case.

2.1. Wasserstein distance. As a preliminary to the introduction of the minimizing movement scheme, we introduce some necessary material related to Wasserstein (or Monge-Kantorovich) distances between nonnegative measures of prescribed mass that are absolutely continuous w.r.t. to the Lebesgue measure. We refer to Santambrogio's monograph [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling[END_REF] for an introduction to optimal transportation and to the Wasserstein distances, and to Villani's big book [START_REF] Villani | Optimal transport, volume 338 of Grundlehren der Mathematischen Wissenschaften[END_REF] for a more complete presentation.

Given two elements c i and či of

A i (i ∈ {1, 2}), a map t : Ω → Ω is said to send c i on či (we write či = t#c i ) if A či (x)dx = t -1 (A) c i (x)dx,
for all Borel subset A of Ω.

The Wasserstein distance W i (c i , či ) with quadratic cost function between c i and či is then defined by

(25) W i (c i , či ) = inf t s.t. či=t#ci 1 m i Ω |x -t(x)| 2 c i (x)dx 1/2 , i ∈ {1, 2}.
In [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusionconvection problem on a triangular mesh[END_REF], the infimum is in fact a minimum, and t is the gradient of a convex function. In our context of fluid flows, the cost for moving the mass of the phase i from a configuration c i to another configuration či regardless to the other phase is equal to W 2 i (c i , či ). The multiplying factor 1/m i is natural since the more mobile is the phase, the less expensive are its displacements. We can then define the tensorized Wasserstein distance W on A by

W 2 (c, č) = W 2 1 (c 1 , č1 ) + W 2 2 (c 2 , č2 ), ∀c, č ∈ A.
In the core of the proof, we will make an extensive use of the Kantorovich dual problem. More precisely, we will use the fact that

(26) 1 2 W 2 i (c i , či ) = sup ϕi∈L 1 (ci),ψi∈L 1 (či) ϕi(x)+ψi(y)≤ |x-y| 2 2m i Ω c i (x)ϕ i (x)dx + Ω či (y)ψ i (y)dy.
Here, L 1 (ρ) denotes the sets of integrable functions for the measure with density ρ. Here again, the supremum is in fact a maximum, and the Kantorovich potentials (ϕ i , ψ i ) achieving the sup in ( 26) are dc i ⊗ dč i unique up to an additive constant. The optimal transportation t i sending c i on či achieving the inf in [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusionconvection problem on a triangular mesh[END_REF] is related to the Kantorovich potential by

t i (x) = x -m i ∇ϕ i (x), ∀x ∈ Ω
with ϕ i achieving the sup in [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF]. As a consequence, the formula ( 25) provides ( 27)

W 2 i (c i , či ) = m i Ω c i |∇ϕ i | 2 dx, i ∈ {1, 2}
to be used in the sequel.

2.2. The JKO scheme and the approximate solution.

We have now all the necessary material at hand to define the minimizing movement scheme. Let τ > 0 and c n-1 ∈ A ∩ X , then define the functional F n τ : A → R ∪ {+∞} by setting

F n τ (c) = E(c) + 1 2τ W 2 (c, c n-1 ), ∀c ∈ A.
The functional F n τ is bounded from below since all its components are. Then we define (28)

c n ∈ argmin c∈A F n τ (c).
In case that F n τ possesses several mimimizers, a selection must be made: if c n-1 1 and c n-1 2 are both positive a.e. in Ω, then any minimizer is fine. However, if one of the two densities vanishes on a set of positive measure, then we need to select c n such that it is an accumulation point in the weak H 2 (Ω) 2 -topology for δ → 0 of the set of minimizers of the functionals

F n,δ τ (c) = E(c) + 1 2τ W 2 (c, c n-1,δ ), with c n-1,δ i = (1 -δ)c n-1 i + δc i (29)
where c i is defined by [START_REF] Cahn | The Cahn-Hilliard equation with a concentration dependent mobility: motion by minus the Laplacian of the mean curvature[END_REF]. Under certain conditions, for instance in the thermally agitated situation where θ 1 , θ 2 > 0, it is known a priori that both c n 1 and c n 2 are strictly positive, and hence no selection is necessary, see e.g. Lemma 3.1. The proof of the existence of c n solution to ( 28) is given below. The proof that one can indeed select the minimizer in the way described above is postponed to Corollary 3.3 in Section 3.4, where we show compactness properties of the set of minimizers. Proposition 2.1. For any c n-1 ∈ A ∩ X , there exists (at least) one solution c n ∈ A ∩ X to the minimization scheme [START_REF] Lieb | Analysis[END_REF].

Proof. Let c n,k k≥0 be a minimizing sequence with c n,0 = c n-1 and such that ( 30)

E(c n,k ) ≤ E(c n-1 ) < ∞ for all k ≥ 0.
We infer from (10) that c n,k ∈ X , hence 0 ≤ c n,k i ≤ 1 for all k ≥ 0 and i ∈ {1, 2}. Moreover, in view of the presence of the Dirichlet energy in E, it follows by means of the Poincaré-Wirtinger inequality that the c n,k i are k-uniformly bounded in H 1 (Ω). Passing to a subsequence if necessary, the c n,k i converge to limits c n i , weakly in H 1 (Ω) thanks to Alaoglu's theorem, and also strongly in L 1 (Ω) by Rellich's compactness lemma. All the components of E are sequentielly lower semi-continuous with respect to this convergence, and so

E(c n ) ≤ lim inf k→∞ E(c n,k ).
Finally, also the Wasserstein distance is lower semi-continuous with respect to strong convergence in L 1 (Ω),

W (c n , c n-1 ) ≤ lim inf k→∞ W (c n,k , c n-1 ).
As a consequence, c n is a minimizer of F n τ . Each sequence (c n ) n≥1 of iterated solutions to the scheme ( 28) is accompanied by approximate phase potentials (µ n ) n≥1 : these are introduced in such a way that in the time-continuous limit, the continuity equations

∂ t c i = ∇ • [m i c i ∇(µ i + Ψ i )]
hold. The suitable quantities µ n i are identified (somewhat a posteriori ) by comparing the optimality conditions for (28) with the persued limit PDE system [START_REF] Ambrosio | A user's guide to optimal transport[END_REF].

In order to justify the formal definition of µ n that we give below in [START_REF] Mielke | A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems[END_REF], we anticipate an auxiliary result from Section 3.4, that can be understood as a formulation of the time-discrete Euler-Lagrange equations. It involves the (backward) Kantorovich potentials ϕ n = (ϕ n 1 , ϕ n 2 ) sending c n on c n-1 . A subtle point is to overcome the inherent non-uniqueness of ϕ n 1 and ϕ n 2 -particularly if one of the c n i 's vanishes on a set of positive measure -by making a suitable selection, as will be explained in Section 3.4. We note that the intricate selection procedure for the minimizer c n in ( 28) enters precisely at this point. Lemma 2.2. At each step n ≥ 1, there exist Kantorovich potentials ϕ n = (ϕ n 1 , ϕ n 2 ) for sending c n to c n-1 such that F n : Ω → R, given by (31)

F n = ϕ n 1 τ - ϕ n 2 τ -α∆c n 1 + f (c n 1 ) + χ(1 -2c n 1 ) + Ψ 1 -Ψ 2 , satisfies F n ≤ 0 a.e. in {c n 1 > 0}, F n ≥ 0 a.e. in {c n 2 > 0}. ( 32
)
With the particular choice of ϕ n from the lemma, we define

µ n = (µ n 1 , µ n 2 ) : Ω → R 2 by µ n 1 := - ϕ n 1 τ -Ψ 1 -θ 1 log(c n 1 ) + (F n ) + , and µ n 2 := - ϕ n 2 τ -Ψ 2 -θ 2 log(c n 2 ) + (F n ) -. (33)
Actually, since F n is invariant under simultaneous addition of a global constant to both ϕ n 1 and ϕ n 2 , we can assume that their normalization is chosen such that

Ω µ n dx = 0, where µ n = c n 1 µ n 1 + c n 2 µ n 2 . ( 34
)
Note that it follows directly from the definition of µ n that

µ n 1 -µ n 2 = -α∆c n 1 + χ(1 -2c n 1 )
a.e. in Ω. [START_REF] Otto | Thermodynamically driven incompressible fluid mixtures[END_REF] Moreover, thanks to (32), we have that

µ n i = - ϕ n i τ -Ψ i -θ i log(c n i ) a.e. in {c n i > 0}. ( 36
)
With these quantities at hand, we define further the piecewise constant interpolants c τ : R + → A ∩ X and µ τ : R + → R 2 by (37)

c τ (t) = (c n 1 , c n 2 ), µ τ (t) = (µ n 1 , µ n 2 ) for t ∈ ((n -1)τ, nτ ].
Any pair (c τ , µ τ ) that has been obtained in this way from an iterated minimizer (c n ) n≥1 will be referred to as τ -approximate solution below.

2.3. Weak solutions. The goal of this section is to state our main result, that is the convergence of the JKO scheme. It requires the introduction of the notion of weak solution that will be obtained at the limit when the approximation parameter τ tends to 0. • Regularity of concentrations:

c i ∈ L ∞ (R + ; H 1 (Ω)) ∩ L 2 loc (R + ; H 2 (Ω)) ∩ C(R + ; L 2 (Ω)); • Regularity of potentials: µ i ∈ L 2 loc (R + ; L q d (Ω)) with q d = 2 if d ∈ {1, 2}
, and q d = 3 2 if d = 3, as well as c i ∇µ i ∈ L 2 (R + , L 2 (Ω)) d , and (normalization)

Ω (c 1 (x, t)µ 1 (x, t) + c 2 (x, t)µ 2 (x, t)) dx = 0,
for a.e. t ≥ 0;

• Initial and boundary conditions:

c i (•, 0) = c 0 i on Ω, and
∇c i • n = 0 on ∂Ω × R + ; • Volume constraint: c 1 + c 2 = 1 a.e. in Ω × R + ; • Continuity equations: for all ξ ∈ C 2 (Ω) and all t 1 , t 2 ∈ R + with t 2 ≥ t 1 , (38)
Ω (c i (x, t 2 ) -c i (x, t 1 ))ξ(x)dx + m i t2 t1 Ω (c i ∇(µ i + Ψ i ) + θ i ∇c i ) • ∇ξdxdt = 0;
• Steepest descent: (20) holds almost everywhere in Ω × R + , i.e.

(39)

µ 1 -µ 2 = -α∆c 1 + χ(1 -2c 1 ).
Here is the convergence theorem for the minimization scheme. The global existence of a weak solution is a by-product. All along the paper, T denotes an arbitrary finite time horizon, and we make use of the shorten notation Q T for the space-time cylinder Ω × (0, T ). Theorem 2.4 (Convergence of the minimizing movement scheme). Assume that parameters α 1 , α 2 , χ > 0 and θ 1 , θ 2 ≥ 0 as well as external potentials Ψ 1 , Ψ 2 ∈ H 1 (Ω) are prescribed. Let c 0 = (c 0 1 , c 0 2 ) be an initial condition satisfying [START_REF] Benamou | An augmented Lagrangian approach to Wasserstein gradient flows and applications[END_REF].

Then, for any sequence (τ n ) n≥1 ⊂ (0, 1) with τ n → 0 and a corresponding sequence of τ napproximate solutions c τn , µ τn n≥1 , one can select a sub-sequence (not relabeled) such that,

c i,τn -→ n→∞ c i in the L ∞ ((0, T ); H 1 (Ω))-weak-sense, c i,τn (•, t) -c i (•, t) L 2 (Ω) -→ n→∞ 0 for all t ∈ [0, T ], c i,τn -→ n→∞ c i in L 2 ((0, T ); W 1,d (Ω)), c i,τn -→ n→∞ c i weakly in L 2 ((0, T ); H 2 (Ω)), µ i,τn -→ n→∞ µ i for the weak topology of L 2 ((0, T ); L q d (Ω)), c i,τn ∇µ i,τn -→ n→∞ c i ∇µ i weakly in L 2 (Q T ) d ,
and the limit (c, µ) is a weak solution in the sense of Definition 2.3.

We recall that the definition of τ -approximate solutions involves the selection of particular minimizers in (28) unless all densities have full support.

Proof of Theorem 2.4

We first establish some estimates on the approximate solution c τ . The very classical energy estimate and some straightforward consequences are exposed in Section 3.1. In Section 3.2, we show that the approximate solution remains bounded away from 0 and 1 if the thermal diffusion coefficients θ i are positive. Then in Section 3.3, we make use of the flow interchange technique initially introduced in [START_REF] Matthes | A family of nonlinear fourth order equations of gradient flow type[END_REF] to get enhanced regularity estimates on the approximate solutions. The Euler-Lagrange equation are then obtained in Section 3.4 thanks to a linearization technique inspired from the work of Maury et al. [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF]. The convergence of the approximate solution towards a weak solution is finally established in Section 3.5.

3.1.

Energy and distance estimate. By definition of c n in ( 28), one has

F n τ (c n ) ≤ F n τ (c n-1 ), i.e., (40) E(c n ) + 1 2τ W 2 (c n , c n-1 ) ≤ E(c n-1 ), ∀n ≥ 1.
Summing (40) over n, and using that E(c) ≥ E > -∞ for all c ∈ A, we obtain the square distance estimate

(41) 1 τ n≥1 W 2 (c n , c n-1 ) ≤ 2 E(c 0 ) -E < +∞.
This readily gives the approximate 1/2-Hölder estimate

(42) W (c n2 , c n1 ) ≤ C |n 2 -n 1 |τ , ∀n 1 , n 2 ≥ 0.
Bearing in mind the definition [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling[END_REF] of the approximate solution c τ , we get that

(43) W (c τ (t), c τ (s)) ≤ C |t -s| + τ , ∀s, t ≥ 0.
We also deduce from the energy estimate (40) that ( 44)

E(c n ) ≤ E(c 0 ) < ∞.
We deduce that ( 45)

Ω |∇c n 1 | 2 dx ≤ 2 α   E(c 0 ) + i∈{1,2} Ψ i L 1   < ∞, ∀n ≥ 0.
3.2. Positivity of the discrete solution in presence of thermal agitation. The formula [START_REF] Peletier | Variational modelling: Energies, gradient flows, and large deviations[END_REF] suggests to give a proper sense to the quantity θ i log(c n i ). This is the purpose of the following lemma, which is an adaptation to our framework of [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling[END_REF]Lemma 8.6].

Lemma 3.1. Let c n be a minimizer of F n τ as in [START_REF] Lieb | Analysis[END_REF]. Assume that θ i > 0, then c n i > 0 a.e. in Ω. Moreover, θ i log(c n i ) ∈ L 1 (Ω). Proof. With a slight abuse of notation, we denote by c = (c 1 , c 2 ) the constant element of X ∩ A given by ( 9). Fix n. Given ∈ (0, 1), we introduce c ∈ X ∩ A by (46)

c i = c i + (1 -)c n i = c n i + (c i -c n i ).
Note that c i > 0 everywhere in Ω, and that

c i -c i = (1 -)(c n i -c i ) (47)
By optimality of c n , the inequality

E therm (c n ) -E therm (c ) ≤ E Dir (c ) -E Dir (c n ) + E chem (c ) -E chem (c n ) + E ext (c ) -E ext (c n ) + 1 2τ W 2 (c , c n-1 ) -W 2 (c n , c n-1 )
holds for all ∈ (0, 1). Then -with generic constants C that may change from line to line, but are independent of -it follows directly from the definition of c in ( 46) that

E Dir (c ) ≤ E Dir (c n ), that E chem (c ) ≤ E chem (c n ) + Ω (c 1 -c n 1 )(1 -2c n 1 )dx ≤ E chem (c n ) + C ,
and that

E ext (c ) = E ext (c n ) + i∈{1,2} c n i -c i , Ψ i L ∞ (Ω),L 1 (Ω) ≤ E ext (c n ) + C .
Moreover, the convexity of W 2 yields

W 2 (c , c n-1 ) ≤ W 2 (c n , c n-1 ) + W 2 (c, c n-1 ) -W 2 (c n , c n-1 ) ≤ W 2 (c n , c n-1 ) + C .
Combining the above inequalities, we obtain that

(48) E therm (c n ) -E therm (c ) = i∈{1,2} θ i Ω (H(c n i ) -H(c i ))dx ≤ C τ .
Let i ∈ {1, 2} be such that θ i > 0, and denote by

A = {x ∈ Ω | c n i (x) > 0}. Then convexity of H implies that (49) H(c n i ) -H(c i ) ≥ (c n i -c i ) log(c i ) = (c n i -c i ) log(c i ) a.
e. in A. On the other hand, since c i = c i a.e. in A c , (50)

H(c n i ) -H(c i ) = 1 -H( c i ) = -c i log( c i ) + c i a.
e. in A c . Integrating (49) over A and (50) over A c and using (48), we obtain after division by > 0 that (51)

-c i log( c i )|A c | + A (c n i -c i ) log(c i ) ≤ C τ .
In view of (47), we have

(c n i -c i ) log(c i ) ≥ (c n i -c i ) log(c i ), (52)
and therefore

-c i log( c i )|A c | + A (c n i -c i ) log(c i ) ≤ C τ .
Letting tend to 0 produces a contradiction unless |A c | = 0. Thus we have proved that A = Ω (up to a negligible set). Moreover, in view of (52), and since (c n ic i ) log(c i ) ∈ L 1 (Ω), Fatou's Lemma applies and leads to

Ω (c n i -c i ) log(c n i )dx ≤ lim inf →0 Ω (c n i -c i ) log(c i )dx ≤ C τ .
This latter inequality imposes that log(c n i ) belongs to L 1 (Ω) for fixed τ > 0.

3.3. Flow interchange and entropy estimate. In the next lemma, our goal is to get an improved regularity estimate on c by the mean of the flow interchange technique.

Lemma 3.2.

There exists C depending only on α, χ, Ψ i , c 0 , m i , T , such that

(53) i∈{1,2} θ i T 0 Ω |∇ √ c i,τ | 2 dxdt + T 0 Ω |∆c τ | 2 dxdt ≤ C, ∀T > 0.
Since Ω is convex, it implies

(54) c τ L 2 ((0,T );H 2 (Ω)) ≤ C.
Moreover, there holds

(55) ∇c n i • n = 0 on ∂Ω.
Proof. Below, s ≥ 0 denotes an auxiliary time variable. Let či (i ∈ {1, 2}) be the unique solution to (56)

     ∂ s či = ∆č i in Ω × (0, ∞), ∇č i • n = 0 on ∂Ω × (0, ∞), (č i ) |s=0 = c n i in Ω. It is easy to check that Ω či (x, s)dx = Ω c 0 i (x)dx, hence each č(s) = (č 1 (s), č2 (s)
) is an admissible competitor in [START_REF] Lieb | Analysis[END_REF], i.e., F n τ (c n ) ≤ F n τ (č(s)). After rearraging terms in F n τ and dividing by s > 0, the passage to the limit s → 0 in this inequality produces

-lim sup s→0 d ds E(č(s)) ≤ d ds s=0 1 2τ W 2 (č(s), c n-1 ) . (57)
To estimate the derivative on the right hand side above, we use that the heat equation ( 56) is the gradient flow of the Boltzmann entropy functional H(c) = Ω H(c)dx, which is displacement convex. Therefore, solutions č to (56) satisfy the Evolution Variational Inequality [1, Definition 4.5] centered at c n-1 i , that is

1 2 d ds W 2 i (č i (s), c n-1 i ) ≤ H(c n-1 i ) -H(č i (s)) m i .
Division by τ and summation over i leads to

(58) d ds 1 2τ W 2 (č(s), c n-1 ) ≤ i∈{1,2} H(c n-1 i ) -H(č i (s)) m i τ , ∀s > 0.
The estimate for the left-hand side of (57) is more combersome. We estimate the s-derivate of the various parts of E individually. To begin with, notice that there is no contribution from E cons , since the constraint č1 (s) + č2 (s) = 1 is preserved by (56). For working on the other parts of E, we use that the solution či of (56) belongs to C ∞ ((0, ∞); H 2 (Ω)), is positive on Ω for s > 0, and satisfies homogeneous Neumann boundary conditions. That makes či regular enough to justify the following calculations, and in particular the integration by parts. For the Dirichlet energy, we obtain

(59) d ds E Dir (č(s)) = d ds α 2 Ω |∇č 1 (s)| 2 dx = -α Ω ∆č 1 (s)∂ s č1 (s)dx = -α Ω |∆č 1 (s)| 2 dx.
The derivative of the chemical energy amounts to

d ds E chem (č(s)) = d ds χ Ω č1 (s)(1 -č1 (s))dx = χ Ω (1 -2č 1 (s))∂ s č1 (s)dx = 2χ Ω |∇č 1 (s)| 2 dx,
and since s → E Dir (č(s)) is non-increasing by the previous calculation, we can further conclude that

(60) d ds E chem (č(s)) ≤ 2χ Ω |∇c n 1 | 2 dx ≤ C,
where the second inequality follows from (45). Next, for the thermal energy,

d ds E therm (č(s)) = i∈{1,2} θ i d ds Ω H(č i (s))dx = i∈{1,2} θ i Ω log(č i (s))∂ s či (s)dx = i∈{1,2} θ i Ω log(č i (s))∆č i (s)dx = - i∈{1,2} θ i Ω ∇ log(č i (s)) • ∇č i (s)dx,
and since či (t) is smooth on Ω and bounded away from 0, we conclude that

(61) d ds E therm (č(s)) = -4 i∈{1,2} θ i Ω |∇ či (s)| 2 dx.
Finally, we obtain for the derivative of the potential

d ds E ext (č(s)) = d ds i∈{1,2} Ω Ψ i či (s)dx = Ω (Ψ 1 -Ψ 2 )∂ s č1 (s)dx = Ω (Ψ 1 -Ψ 2 )∆č 1 (s)dx ≤ α 2 ∆č 1 (s) 2 L 2 + 1 2α Ψ 1 -Ψ 2 2 L 2 .
In combination with (59), this implies

(62) d ds E Dir (č(s)) + E ext (č(s)) ≤ - α 2 ∆č 1 (s) 2 L 2 + C, ∀s > 0.
Since the initial condition c n i in (56) belongs to H 1 (Ω), we have that či (s) → c n i in H 1 (Ω) as s → 0. In particular, H(č i (s)) → H(c n i ) as s → 0. Now we substitute (58) and ( 60), (61), (62) into (57), obtaining

lim sup s→0    Ω |∆č 1 (s)| 2 dx + i∈{1,2} θ i Ω ∇ či (t) 2 dx    ≤ C   1 + i∈{1,2} H(c n-1 i ) -H(c n i ) m i τ   . (63)
Since the initial condition c n i in (56) belongs to H 1 (Ω), we have that či (s) → c n i in H 1 (Ω) as s → 0. Convexity of the the functionals c → |∆c| 2 dx and c → |∇ √ c| 2 dx implies lower semi-continuity with respect to convergence in H 1 (Ω), and therefore

Ω |∆c n 1 | 2 dx + i∈{1,2} θ i Ω ∇ c n i 2 dx ≤ C   1 + i∈{1,2} H(c n-1 i ) -H(c n i ) m i τ   , ∀n ≥ 1.
Multiplying by τ and summing over n ∈ 1, . . . , T τ leads to (53); here we use that 0 ≤ i∈{1,2} H(c i ) ≤ 2|Ω| for all measurable c 1 , c 2 : Ω → [0, 1], and that |∆c n 2 | = |∆c n 1 |. Concerning the boundary condition (55): observe that the solution či to (56) satisfies in particular ∇č i (s) • n = 0 at each s > 0. Since Ω is convex, one can show (see [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF]Chapter 3]) that

Ω ∇ 2 či (s) 2 dx ≤ Ω |∆č i (s)| 2 dx.
In combination with the fact that či (s) has values in [0, 1] only, it follows that the L 2 -norm of ∆č i (s) controls the full H 2 -norm of či (s), i.e., či (s

) 2 H 2 (Ω) ≤ C 1 + Ω |∆č i (s)| 2 dx . ( 64
)
On the other hand, (63) implies in particular that ∆č i (s) remains bounded in L 2 as s → 0. By Alaoglu, and since we already know that či (s) tends to c n i in H 1 (Ω), it follows that či (s) converges to c n i weakly in H 2 (Ω) as s → 0. Now, since the trace operator mapping c ∈ H 2 (Ω) to ∇c • n ∈ L 2 (∂Ω) is weakly continuous, we conclude that ∇c n i • n = 0 as well. Finally, we obtain (54): since we have just shown that c n i satisfies homogeneous Neumann boundary conditions, the same argument as above shows that (64) holds also with c n i in place of či (s). Summation of (53) over n directly produces (54).

Corollary 3.3. At each step n ≥ 1, there exists among the minimizers in (28) at least one which is a weak H 2 (Ω) 2 -limit point for δ → 0 of the set of minimizers of the functionals F n,δ τ defined in [START_REF] Lisini | Cahn-Hilliard and thin film equations with nonlinear mobility as gradient flows in weighted-Wasserstein metrics[END_REF].

Proof. Since c n-1,δ → c n-1 in L 1 (Ω) 2 as δ → 0, the F n,δ τ Γ-converge to F n
τ , for instance in the topology induced by W on A. In complete analogy to the previous proof, we obtain for each minimizer c n,δ of F n,δ τ that

c n,δ 2 H 2 ≤ C   1 + i∈{1,2} H(c n-1,δ i ) -H(c n,δ i ) m i τ   ≤ C τ .
Therefore, there is a sequence c n,δ k that converges weakly in H 2 (Ω) 2 to a limit c n that minimizes F n τ . 3.4. Euler Lagrange equations. The goal of this section is to characterize the minimizer c n of (28). We have anticipated the main result already in Lemma 2.2. It is a consequence of the following linearized optimality condition for c n . Lemma 3.4. In each step n ≥ 1, there exist Kantorovich potentials ϕ n = (ϕ n 1 , ϕ n 2 ) such that the quantity F n given in [START_REF] Matthes | A convergent Lagrangian discretization for a nonlinear fourth-order equation[END_REF] satisfies the following linearized optimality condition:

(65) Ω F n c n 1 dx ≤ Ω F n c 1 dx, ∀c = (c 1 , c 2 ) ∈ A.
Proof. The proof is inspired from [32, Lemma 3.1], see also [START_REF] Cancès | Incompressible immiscible multiphase flows in porous media: a variational approach[END_REF]Lemma 3.2]. Assume first that c n-1 i > 0 almost everywhere in Ω, so that the Kantorovich potentials ϕ n 1 , ψ n 1 from c n 1 to c n-1 1 are unique up to addition of a global constant, that we fix at an arbitrary value for this proof. For any given c ∈ A ∩ X , we perturb c n into (66) c = (1 -)c n + c, ∀ ∈ (0, 1).

Clearly, c belongs to A ∩ X . Defining (ϕ i , ψ i ) the Kantorovich potential from c i to c n-1 i , we infer from [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF] that

     1 2 W 2 i (c i , c n-1 i ) = Ω ϕ i (x)c i (x)dx + Ω ψ i (y)c n-1 i (y)dy, 1 2 W 2 i (c n i , c n-1 i ) ≥ Ω ϕ i (x)c n i (x)dx + Ω ψ i (y)c n-1 i (y)dy.
Subtracting the two above relations and using the definition (66) of c , one gets

1 2τ W 2 i (c i , c n-1 i ) -W 2 i (c n i , c n-1 i ) ≤ Ω ϕ i τ (c i -c n i ) dx.
Hence, using c n , c ∈ X , one gets that

(67) 1 2τ W 2 (c , c n-1 ) -W 2 (c n , c n-1 ) ≤ Ω ϕ 1 τ - ϕ 2 τ (c 1 -c n 1 ) dx.
On the other hand, the convexity of E Dir and E therm , the linearity of E ext , and the concavity of

E chem yield (68) E(c ) -E(c n ) ≤ Ω α∇c 1 • ∇(c 1 -c n 1 )dx + Ω χ(1 -2c n 1 )(c 1 -c n 1 )dx + Ω f (c 1 )(c 1 -c n 1 )dx + Ω (Ψ 1 -Ψ 2 ) (c 1 -c n 1 )dx.
Bearing in mind that c n 1 is a minimizer, the combination of ( 67) with (68) leads to

0 ≤ F n τ (c ) -F n τ (c n ) ≤ Ω ϕ 1 τ - ϕ 2 τ (c 1 -c n 1 ) dx + Ω α∇c 1 • ∇(c 1 -c n 1 )dx + Ω χ(1 -2c n 1 )(c 1 -c n 1 )dx + Ω f (c 1 )(c 1 -c n 1 )dx + Ω (Ψ 1 -Ψ 2 ) (c 1 -c n 1 )dx, ∀ > 0.
We can consider the limit → 0 in the right-hand side of the above expression. From the definition (66) of c , it is clear that c 1 converges in H 1 (Ω) towards c n 1 and that f (c 1 ) converges in L 1 (Ω) towards f (c n 1 ), while ϕ i converges uniformly towards ϕ n i (see for instance [37, Theorem 1.52]), so that (65) holds thanks to Lemma 3.2.

Assume now that c n-1 i = 0 on some part of Ω. By definition of c n in and after [START_REF] Lieb | Analysis[END_REF], there exists a sequence of minimizers c n,δ of the respective functionals F n,δ τ given in (29) that converge weakly in H 2 (Ω) 2 to c n . Defining F n,δ as F n in (31) upon replacing c n i by c n,δ i and ϕ n i by the Kantorovich potential ϕ n,δ i sending c n,δ i to c n-1,δ i , then we get from the reasoning above that (69)

Ω F n,δ (c 1 -c n,δ 1 ) ≥ 0, ∀c = (c 1 , c 2 ) ∈ A ∩ X .
The weak convergence of the c n,δ i to c n i in H 2 (Ω), and the induced uniform convergence of the Kantorovich potentials ϕ n,δ i to Kantorovich potentials ϕ n i sending c n i to c n-1 i (cf. [37, Theorem 1.52]) is sufficient to deduce (65) from (69) in the limit δ → 0.

We are now in the position to prove Lemma 2.2 which has been essential for the definition of the potentials µ.

we deduce from (72) the desired L 2 ((0, T ); L q d (Ω)) estimates on the phase potentials µ n i with q d defined in Definition 2.3. Finally, the combination of the relations (36), ( 27), ( 53) and (41) yields

Q T c i,τ |∇µ i,τ | 2 dxdt ≤ C.
The following lemma is a first step towards the recovery of the weak formulation [START_REF] Villani | Optimal transport, volume 338 of Grundlehren der Mathematischen Wissenschaften[END_REF].

Lemma 3.6. For any ξ ∈ C 2 (Ω), there holds

(75) Ω c n i -c n-1 i ξdx + τ m i Ω (c n i ∇ (µ n i + Ψ i ) + θ i ∇c n i ) • ∇ξdx ≤ 1 2 W 2 i (c n i , c n-1 i ) D 2 ξ ∞ .
Proof. The optimal transport map

t n i (x) = x -m i ∇ϕ n i (x), ∀x ∈ Ω sending c n i to c n-1 i maps Ω into itself because Ω is convex. Therefore, since c n-1 i = t n i #c n i and thanks to (36), one gets that Ω c n i -c n-1 i ξdx + τ m i Ω (c n i ∇ (µ n i + Ψ i ) + θ i ∇c n i ) • ∇ξdx = Ω (ξ(x) -ξ(t n i (x)) -m i ∇ξ(x) • ∇ϕ n i (x)) c n i (x)dx
for all ξ ∈ C 2 (Ω). The Taylor expansion of ξ at point x provides

|ξ(t n i (x)) -ξ(x) + m i ∇ξ(x) • ∇ϕ i (x)| ≤ 1 2 D 2 ξ ∞ |t n i (x) -x| 2 , ∀x ∈ Ω, so that Ω (c n i -c n-1 i )ξdx + τ m i Ω (c n i ∇ (µ n i + Ψ i ) + θ i ∇c n i ) • ∇ξdx ≤ 1 2 D 2 ξ ∞ Ω |t n i (x) -x| 2 c n i (x)dx,
which is exactly the desired result.

3.5. Convergence towards a weak solution. The goal of this section is to consider the limit τ → 0. This requires some compactness on the approximate phase field c τ and on the approximate potential µ τ . In what follows, A is equipped with the topology corresponding to the distance W .

Proposition 3.7.

There exist c ∈ C([0, T ]; L 2 (Ω)) ∩ L 2 ((0, T ); H 2 (Ω)) ∩ L ∞ ((0, T ); H 1 (Ω)) with c(t) ∈ A ∩ X for a.e. t ∈ [0, T ],
and µ ∈ L 2 ((0, T ); L q d (Ω)) such that, up to the extraction of a subsequence, the following convergence properties hold:

c i,τ -→ τ →0 c i in the L ∞ ((0, T ); H 1 (Ω))-weak-sense, (76a) c i,τ (•, t) -c i (•, t) L 2 (Ω) -→ τ →0 0 for all t ∈ [0, T ], (76b) c i,τ -→ τ →0 c i in L 2 ((0, T ); W 1,d (Ω)), (76c) c i,τ -→ τ →0 c i weakly in L 2 ((0, T ); H 2 (Ω)), (76d) µ i,τ -→ τ →0
µ i for the weak topology of L 2 ((0, T ); L q d (Ω)), (76e)

c i,τ ∇µ i,τ -→ τ →0 c i ∇µ i weakly in L 2 (Q T ) d . (76f)
Proof. All the convergence properties stated below occur up to the extraction of a subsequence when τ tends to 0. We deduce from Estimate (45) that the family (c τ ) τ >0 is bounded in L ∞ ((0, T ); H 1 (Ω)). Hence we can assume that c τ tends to some c in the L ∞ ((0, T ); H 1 (Ω))-weak-sense. Moreover, since 0 ≤ c i,τ ≤ 1, we also have that 0 ≤ c i ≤ 1 a.e. in Q T . We also infer from Estimate (54) that c τ converges weakly in L 2 ((0, T ); H 2 (Ω)) towards c.

As a consequence of the L ∞ bound on c i,τ and of the Benamou-Brenier formula, we get that c

ic

(2) i H 1 (Ω) ≤ 1 m i W i (c (1)
i , c

i ), ∀c

i , c

(2) i ∈ A i , (see more precisely [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling[END_REF]Lemma 3.4]). Therefore, we infer from (43) that This implies in particular some strong convergence in L 2 (Q T ), which can be combined with the weak convergence in L 2 ((0, T ); H 2 (Ω)) thanks to interpolation arguments to derive some strong convergence in L 2 ((0, T ); H s (Ω)) for any s < 2. The continuous embedding of H s (Ω) into W 1,d (Ω) when s ≥ 1 + max(0, d-2 2 ) (see for instance [START_REF] Di Nezza | Hitchhikers guide to the fractional Sobolev spaces[END_REF]Theorem 6.7]) ensures that (77)

c τ -→ τ →0 c in L 2 ((0, T ); W 1,d (Ω)) 2 .
Let us switch to the phase potentials µ. Thanks to Lemma 3.5, we have (uniform w.r.t. τ ) L 2 ((0, T ); L q d (Ω)) estimates on µ i,τ . Hence there exists µ in L 2 ((0, T ); L q d (Ω)) 2 such that (78) µ i,τ -→ τ →0

µ i weakly in L 2 ((0, T ); L q d (Ω)).

In Lemma 3.5, we also established a (uniform w.r.t. τ ) L 2 (Q T ) d estimate on √ c i,τ ∇µ i,τ τ >0 .

Since 0 ≤ c i,τ ≤ 1, it implies a uniform L 2 (Q T ) d estimate on (c i,τ ∇µ i,τ ) τ >0 . Therefore, there exists ϑ i ∈ L 2 (Q T ) d such that c i,τ ∇µ i,τ converges weakly in L 2 (Q T ) d to ϑ i as τ tends to 0. It remains to show that ϑ i = c i ∇µ i . First, the distributions c i and ∇µ i can be multiplied since c i belongs to L 2 ((0, T ); W 1,d (Ω)) and ∇µ i belongs to L 2 ((0, T ); W -1,q d (Ω)) with q d ≤ d d-1 . Moreover, for all φ ∈ C ∞ c (Q T ) d , one has (79)

Q T c i,τ ∇µ i,τ • φdxdt = - Q T µ i,τ (∇c i,τ • φ + c i,τ ∇ • φ) dxdt.
Thanks to (77) and (78), we can pass in the limit in the right-hand side of the above expression. This leads to (80)

Q T c i,τ ∇µ i,τ • φdxdt -→ τ →0 - Q T µ i (∇c i • φ + c i ∇ • φ) dxdt = c i ∇µ i , φ D ,D .
As a consequence, ϑ i = c i ∇µ i in the distributional sense, thus also in L 2 (Q T ). Since c τ (t) converges in L 2 (Ω) 2 towards c(t) for all t ∈ [0, 1], there holds

Ω c i (t)dx = Ω c 0 i dx, c 1 (x, t) + c 2 (x, t) = 1 for all t ∈ [0, T ].
Moreover, since c(t) belongs to H 1 (Ω) for a.e. t ∈ (0, T ), c(t) belongs to A ∩ X for a.e. t ∈ (0, T ). This concludes the proof of Proposition 3.7.

We have all the necessary convergence properties to pass to the limit τ → 0 and to identify the limit (c, µ) exhibited in Proposition 3.7 as a weak solution in the sense of Definition 2.3. Proof. Since c i,τ and µ i,τ tend weakly in L 2 ((0, T ); H 2 (Ω)) and L 2 ((0, T ); L q d (Ω)) towards c i and µ i respectively, we can pass to the limit in [START_REF] Otto | Thermodynamically driven incompressible fluid mixtures[END_REF]. Moreover, one can also pass to the limit in the relation 0 = ∇c 1,τ • n established in Lemma 3.2, leading to ∇c 1 • n = 0 on ∂Ω.

It only remains to recover the weak formulation [START_REF] Villani | Optimal transport, volume 338 of Grundlehren der Mathematischen Wissenschaften[END_REF] τ Ω (c i,τ ∇ (µ i,τ + Ψ i ) + θ i ∇c i,τ ) • ∇ξdxdt

≤ 1 2 D 2 ξ ∞ t 2 τ n= t 1 τ W 2 i (c n i , c n-1 i ) ≤ Cτ,
the last inequality being the consequence of the squared distance estimate (41). We can pass to the limit τ → 0 in the above relation thanks to Proposition 3.7.

We have finally proved Theorem 2.4 that is a combination of Propositions 3.7 and 3.8.

  .

Figure 1 . 1 . 4 . 2 .

 1142 Figure 1. Evolution of the energy E(c) along time for the non-local model (17) and the local one[START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF]. The decay of the energy is faster for the non-local model, as shown in Section 1.3.

Figure 2 .

 2 Figure 2. Evolution along time of the numerical solution to the non-local problem (17) (top) and to the local problem (22) (bottom). Snapshots at time t = 10 -2 (left), t = 2.10 -2 (middle), and t = 10 -1 (right).

Figure 3 .

 3 Figure 3. Snapshots at times t = 6.10 -3 (left), t = 5.10 -2 (middle), and t = 1 (right) illustrating the spinodal decomposition governed by model (17) (top) and model (22) (bottom).

2 .Figure 4 .

 24 Figure 4. Evolution of the energy E(c) along time for the non-local model (17) and the local one[START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] for the spinodal decomposition test case.

Definition 2 . 3 .

 23 A pair (c, µ) of phase concentrations c = (c 1 , c 2 ) : Ω → [0, 1] and phase potentials µ = (µ 1 , µ 2 ) : Ω → R is said to be a weak solution to the problem (17), (2), and (3) if

c

  i,τ (t)c i,τ (s) H 1 (Ω) ≤ C |t -s| + τ , ∀s, t ∈ [0, T ]. Let ∆t > 0 and let t ∈ [0, T -∆t], then c i,τ (t + ∆t)c i,τ (t) 2 L 2 (Ω) ≤ c i,τ (t + ∆t)c i,τ (t) H 1 c i,τ (t + ∆t)c i,τ (t) (H 1 ) ≤ C √ ∆t + τ .Bearing in mind the L ∞ ((0, T ); H 1 (Ω)) estimate on c i,τ , we can apply a refined version of the Arzelà-Ascoli theorem [2, Prop. 3.3.1] to obtain that c i ∈ C([0, T ]; L 2 (Ω)) and thatc i,τ (t) -→ τ →0 c i (t) in L 2 (Ω) for all t ∈ [0, T ].Together with the estimate 0 ≤ c i,τ ≤ 1, we deduce that c ∈ C([0, T ]; L p (Ω)) 2 and thatc τ -→ τ →0 c in L p (Q T ) 2 , ∀p ∈ [1, +∞).

Proposition 3 . 8 .

 38 Let (c, µ) be as in Proposition 3.7, then (c, µ) is a weak solution to the problem (17), (2)-(3) in the sense of Definition 2.3.

  ). Hence c's time derivatives ∂ t c i are determined as well. At this point, it is natural to define the chemical potential µ i of the phase i by

	(19)

  . Let t 1 , t 2 ∈ [0, T ] with t 2 ≥ t 1 , then summing (75) over n ∈ t1 τ + 1, . . . , t2

	t 2 τ	τ
	t 1
	τ

τ yields Ω (c i,τ (t 2 )c i,τ (t 1 )) ξdx + m i

Acknowledgements. The authors warmly thank the referees for their valuable remarks and suggestions. This research was supported by the DFG Collaborative Research Center TRR 109, "Discretization in Geometry and Dynamics" and by the French National Research Agency (ANR) through grants ANR-13-JS01-0007-01 (project GEOPOR) and ANR-11-LABX-0007-01 (Labex CEMPI).

Proof of Lemma 2.2. Since 0 ≤ c n 1 ≤ 1, the minimization problem (65) can be solved thanks to the bathtub principle [START_REF] Lieb | Analysis[END_REF]Theorem 1.14]. It amounts to saturate the sublevel sets of F n with c n 1 = 1 until all the mass Ω c 0 i has been allocated. This implies the existence of some ∈ R such that (70)

Given F n , the solution c n 1 to the minimization problem (65) is in general not unique since a prescribed amount of mass can be distributed on the level set {F n = } in different ways without changing the value of the functional. But this lack of uniqueness does not affect the relations (70).

Recall that we still have the freedom to change the Kantorovich potential ϕ n 1 in the definition (31) of F n by addition of a global constant. We choose that constant to enforce = 0 in (70) above. With that normalization, (70) implies [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF].

Lemma 3.5. The chemical potentials µ satisfy the following τ -uniform estimates:

Proof. Recall the definition of µ in [START_REF] Mielke | A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems[END_REF]. It follows from ( 53) that (72)

The definition of µ n implies that

Using the triangle inequality, Cauchy-Schwarz inequality, and 0 ≤ c n i ≤ 1, we get

Since c n ∈ A, one has Ω c n i dx = Ω c 0 i dx. Thus it follows from ( 27), (41), ( 45) and ( 72) that

Bearing [START_REF] Modica | Un esempio di γ-convergenza[END_REF] in mind, we can use the Poincaré-Sobolev inequality and get that (74) µ τ L 2 ((0,T );L d/(d-1) (Ω)) ≤ C.


Jyotirmoy Chatterjee

Eric Gautier

Marc Veillerot

Ricardo C Sousa

Stéphane Auffret

Bernard Dieny

Physicochemical origin of improvement of magnetic and transport properties of STT-MRAM cells using Tungsten on FeCoB storage layer

We investigated and compared the structural and magnetic properties of MgO/FeCoB based out-of-plane magnetized tunnel junctions at thin film level as well as the magneto-transport properties of corresponding patterned STT-MRAM cells comprising either Ta 1 nm or W2 /Ta1 nm cap layers for different annealing temperatures up to 455°C. W material in the cap was found to improve the structural stiffness of the perpendicular magnetic tunnel junctions and most importantly prohibits Fe diffusion from the FeCoB storage layer to the cap layer, remarkably improving the thermal robustness and magneto-transport properties of the stacks and of the corresponding patterned memory cells. As a result, the interfacial anisotropy constant of the MgO/FeCoB interfaces is improved by 17-29% compared to Ta cap. The STT-MRAM cells fabricated from the pMTJ stacks with W/Ta cap reveal a significant improvement of tunneling magnetoresistance and thermal stability factor, which are respectively 120% and 52 as compared to 70% and 35 for the stack with Ta cap. This improvements are ascribed to the enhancement of MgO crystallinity upon higher temperature annealing (425°C) as well as prohibition Fe out-diffusion.

Perpendicular spin transfer torque random access memory (pSTT-MRAM) based on perpendicular magnetic tunnel junction (pMTJ) stack has attracted considerable interest due to their combination of assets of non-volatility, high thermal stability, low critical current for current-induced spin transfer torque (STT) magnetization switching (few tens of microamps at sub-40nm), high speed (write access time typically in the range 10ns-30ns) and high density memory array (4Gbit capacity and 1.5F pitch demonstrated). [1][2][3][4][START_REF] Apalkov | Proc. IEEE[END_REF][START_REF] Nguyen | [END_REF] The core component of the stack, where spin transport phenomena such as tunneling magnetoresistance (TMR) and magnetization reversal of storage layer by STT occurs is an FeCoB/MgO/FeCoB based tunnel junction. [7][8][9] Here one of the FeCoB layer is storage layer, while other has a fixed magnetization, called as reference layer. The reference layer is commonly pinned by exchange interaction with Co/Pt or Co/Pd based synthetic antiferromagnet (SAF). 10,11 The perpendicular magnetic anisotropy (PMA) required to provide a long enough retention of the memory cell originates from the MgO/FeCoB interface due to Fe-O orbitals hybridization. [12][13][14] However, a cap layer used on top of storage layer, can also indirectly influence the magnetic and transport properties of the pMTJ by modifying the interfacial and bulk physicochemical and electronic properties. [15][16][17][18] The main physicochemical modifications are Boron absorption from the storage layer 17 or interdiffusion during sputtering or post deposition anneals. 15 In an earlier paper, we reported significant improvement of annealing tolerance using thick W/Ta cap as compared to Ta because of structural stiffening of the stack and much weaker interdiffusion with W than with Ta. In this report, we show that the interfacial PMA of MgO/FeCoB interface with W2/Ta1 nm cap is larger than with Ta 1nm cap after annealing at all investigated annealing temperatures (315°C-455°C) and explain the physicochemical reasons behind this improvement. In addition, we investigated the magneto-transport properties of sub-100nm diameter memory cells fabricated from pMTJ stacks with these two different cap layers and with different annealing up to 425°C. Note that this annealing temperature is higher than back-end-of-line (BEOL) process temperature. A significant improvement in TMR (from 70% up to 120%) and thermal stability factor (∆) (from 35 up to 52) is reported for 80 nm devices when a Ta cap layer is replaced by a W/Ta cap.

The samples were deposited on Si wafer by magnetron sputtering under an Ar pressure of 2×10 -3 mbar. The MgO tunnel barrier was obtained by naturally oxidizing a metallic Mg layer under an oxygen pressure of 3×10 -2 mbar with a flow rate of 100 sccm. On top of this oxidized layer, a second Mg layer 0.5 nm thick was deposited. All the samples were annealed for 10 min at different temperatures under high vacuum (5×10 -6 mbar). The effective perpendicular anisotropy was calculated from the area between perpendicular and in-plane M(H) loops measured by vibrating sample magnetometer (VSM). At first, the magnetic parameters, such as magnetic dead layer thickness (td), saturation magnetization (Ms) and interfacial anisotropy constant (Ki) of the storage electrodes consisting of Ta 3/FeCoB 0.4/Mg 0.8/Oxidation 30s/Mg 0.5/FeCoB 1.2/cap/Pt 3nm were derived as a function of annealing temperature and are listed in table 1. The values of Ms and td were calculated from the slope and x-intercept of the m/A (magnetic moment per unit area) versus t plots, where t is the nominal thickness determined from the deposition rate. The dead layer thicknesses are observed to be similar for both caps. However, the values of Ms of FeCoB at annealing temperatures ≥ 400 °C are higher for W/Ta than with Ta cap. Interfacial anisotropy constants (Ki) were evaluated from the y-intercepts after fitting the linear zone of Keff.(t-td) vs. (t-td) graphs using the equation

= - -2 (1)
Figure 1 shows Keff.(t-td) vs. (t-td) graphs at various annealing stages for W/Ta and Ta caps. Ki for W/Ta cap is 1.05 erg/cm 2 , which is about 24% higher than that of Ta cap.

At this point, a question is raised, why W in cap layer yields higher interfacial anisotropy?

In order to understand the reason of interfacial anisotropy enhancement for storage layer caped with W instead of Ta, chemical profiling with dynamic secondary ion mass spectroscopy (DSIMS) were carried out for two stacks with simpler structure Si/ Mg 0.7/Oxidation/Mg 0.5/FeCoB 1.6/(W2 or Ta2)/Pt 3 nm. We discuss the comparisons between the SIMS profiles for Ta and W caps depicted in Fig. 2, focusing in particular at the MgO top interface, as this interface plays a determinant role in the PMA amplitude. Note that, the MgO-barrier/FeCoB interface is nominally identical in both structures. However, the comparisons of the Fe elemental profiles with Ta cap (Fig. 2a) and with W cap (Fig. 2d) in the as-deposited samples already reveal a clear difference between the two samples. In the sample with Ta cap, an enrichment of Fe content is observed next to the Ta interface. A pronounced peak of Fe concentration is indeed observed next to this interface and a weaker one next to the MgO interface implying preferred Fe migration towards the Ta cap during deposition. In contrast, in the sample with W cap, the Fe concentration is higher next to the MgO interface as shown in W cap is ascribed to the much larger negative enthalpy of formation of FeTa (-3468 meV/atom) than of FeW (-554 meV/atom). 19 Comparatively, the enthalpy of formation of CoTa (-253meV/atom) and CoW(-84meV/atom) are relatively weak which explains the absence of significant Co migration towards the cap layer. As a result of using W cap, a higher Fe concentration along the MgO interface can be maintained upon annealing than in the case of Ta cap. Therefore, larger interfacial anisotropy with W cap than with Ta cap is produced as the interfacial anisotropy at Fe/MgO interface is larger than at Co/MgO interface. 13,14 As a consequence, the effective magnetic critical thicknesses (tc), listed in table 1, above which the FeCoB layer becomes in-plane magnetized are higher for W cap than for Ta cap. In addition to Fe, there is also a significant difference in B profile between the two caps observed from the SIMS profile. Both ass deposited samples have higher Boron concentration near the MgO interface as well as near the interface with cap layer. The peak near MgO suggests formation of B-oxides, which is in agreement with previous reports. 15,17 Although both cap layers absorb B out from the FeCoB layer upon annealing, Ta serves as a better B getter than W due to the decreasing absolute value of enthalpy of formation of metal borides along the series (Hf, Ta, W, Re, Os) . 20 It is clear from the figure that for both caps, Boron is not completely absorbed by the cap layer. After annealing at 400°C or 450°C, some Boron still remains along the MgO interface. For W cap, higher amount of B seems to be left over near the MgO interface as compared to the Ta cap case. Nevertheless, W cap yields higher PMA suggesting that the prohibition of preferential Fe migration with W cap is the dominant mechanism.

After having focused on the magnetic properties of the top electrode storage layer and investigated the role of W in the cap layer, pMTJ stacks using Ta 1nm and W2/Ta 1nm caps were patterned to characterize the magneto-electric properties of memory cells as a function of annealing temperature. For this purpose, three samples were deposited on 100 mm Si wafers. antiferromagnetically coupled by 0.9 nm Ru layer. 16,[START_REF] Cuchet | Magnetic and Transport Properties of Single and Double Perpendicular Magnetic Tunnel Junctions[END_REF] Co/Pt MLs based SAF are widely used in the pMTJ stacks [START_REF] Cuchet | [END_REF]23 because of their large PMA [START_REF] Weller | Proc. Mater. Res. Soc. Conf. (Materials Res. Soc[END_REF][START_REF] Chatterjee | [END_REF][26] , which makes them good candidates for realizing magnetically hard layers.

One pMTJ with

TMR versus FeCoB thicknesses for the three samples are shown in Fig. 3. TMR increases as a function of FeCoB thickness due to the improvement of the crystalline quality of the FeCoB layer and to the reduced amount of thermal fluctuations at room temperature in thicker magnetic layer. Above the critical thickness, TMR reduces again as the PMA of the storage layer decreases and correlatively the degree of magnetic antiparallel alignment between storage layer and reference layer. For W/Ta cap, the maximum TMR (120%) is larger than the one with Ta cap (70%). Moreover, the stack with W/Ta after annealing 425°C yields the highest TMR (120%) since both the storage layer and the MgO barrier then exhibit the best crystallinity among the three wafers thanks to the high temperature annealing. 27 High-resolution transmission electron microscopy imaging of pMTJ stacks were performed to investigate the crystallinity upon annealing at 340°C and 425°C, which are shown in Fig. 4(a) and (b) respectively. An improved MgO crystallinity after annealing at 425°C as compared with 340°C annealing can be observed, which is the one of the possible reasons of obtaining larger TMR for W cap. Moreover, the reduced Fe out-diffusion when using W cap likely improves the spin-polarization of ∆1 tunneling states resulting in higher TMR than with Ta cap. This is in agreement with the report of F. Bonell et al. where they demonstrated that TMR decreases for higher Co content in FeCo, due to onset of minority ∆1 states at Fermi energy. 28 Lastly, ab initio calculations of FeCo/MgO/FeCo tunnel junctions with different cap layers also demonstrated increased TMR with W cap compared with Ta due to the improvement of spinpolarization of ∆1 tunneling states. 29 The W cap layer on storage layer shows nanocrystalline or even amorphous structure (see Fig. 4). In our previous work 16 , the improvement of thermal tolerance was explained in terms of overall improvement of the mechanical stiffness due to the incorporation of W as its melting temperature is very high, 3422°C. The observation that W is nanocrystalline or even amorphous further supports this explanation ruling out the possibility of crystallinity of W being responsible for high thermal tolerance as stated before by G. An et al. 30 The thermal stability factors of patterned memory cells were calculated from statistical measurements of resistance versus magnetic field (R(H)) loops using the switching field distribution model. 31,32 According to this model, the thermal stability factor ∆ can be derived by fitting the switching probability of the cells using the equation below :

= 1 -exp [-2 Δ !√∆ 1 - $] (2)
In this equation f0, r and Hk are the attempt frequency (1 GHz), field sweep rate (typically 40 kOe/s) and anisotropy field respectively. For each memory cells, 300 R(H) loops were measured. Then the probability of switching was calculated and fitted with the equation to extract ∆. Following this method, ∆ of 80 nm patterned cells were calculated and plotted as a function of effective storage layer thickness in Fig. 5. The sample with Ta 1 nm cap layer shows an average ∆ about 35 at the effective thickness of 0.85 nm, where the maximum TMR was obtained. For the samples with W2/Ta1 nm cap, the average value of ∆ obtained is 52 at the effective thickness 0.95 nm (shown in Fig. 5 (b) and (c)), which is about 50% higher than in the sample with Ta cap. One must note that the error bars are large due to process variation and selection of 80 nm cells with Mg wedge thickness in the range of 1.5 Å. The thermal stability factor increases as a function of thickness, which can be fitted assuming the energy barrier of magnetization reversal as EB= KeffVn where, Keff is effective perpendicular anisotropy of storage layer and Vn is the nucleation volume. The nucleation diameter is comparable to the domain wall width (δw) and therefore, ∆ can be expressed by the equation below. 33 Annealing temperature (°C)

Fig. 2d .

 2d Fig. 2d. In the annealed sample with Ta cap, the Fe concentration enhancement near the Ta

  Ta 1nm cap, annealed at 340 °C and the other two with W2/Ta1 nm annealed at 400 °C and 425 °C for 30mins. The layer configurations of these pMTJ stacks was Si/bottom electrode/Pt 5/(Co/Pt) ML based SAF/Ta 0.3/FeCoB 1.2/Mg-wedge/Oxidation-240s/Mg 0.5/FeCoB-wedge/cap ( Ta 1/Pt 3 or W 2/Ta 1/Pt 3 nm). One must note that the directions of wedges are orthogonal to each other. The SAF layer consists of hard layer (HL: [Co 0.5/Pt 0.25]6/Co 0.5 nm) and reference layer (RL: [Co 0.5/Pt 0.25]3/Co 0.5/Ta 0.3/FeCoB 1.2 nm)

Figure 1 :Figure 2 :FIG. 3 : 16 FIG. 4 :

 123164 Figure 1: Keff.(t-td) vs. (t-td) at different annealing temperatures of FeCoB storage electrodes

FIG. 5 :

 5 FIG. 5: Thermal stability factor (∆) of 80 nm memory cells as a function of storage layer

Table 1 :

 1

			Ta 1nm cap			W2/Ta1 nm cap
		Ms	td	tc	Ki	Ms	td	tc	Ki
		(emu/cm 3 )	(Å)	(Å)	(erg/cm 2 )	(emu/cm 3 )	(Å)	(Å)	(erg/cm 2 )
	315	1216±26 5.5 9.6 0.78±0.02		-	-	-
	340	1262±14 5.6 9.6 0.87±0.03	1210±26	5.5 10.5 1.02±0.037
	400	1252±54 6.1 8.1 0.85±0.025 1368±49 6.25 9.4 1.06±0.028
	425	-	-			1340±50	6.2 9.4 1.05±0.034
	455	-	-			1375±41	6.5 9.4 1.13±0.038

Saturation magnetization (Ms), dead layer thickness (td), effective magnetic critical thickness (tc) and interfacial anisotropy constant (Ki) of FeCoB storage electrode with Ta 1 and W2/Ta1 nm cap at various annealing temperatures.

Acknowledgement

This work was funded under the ERC Adv grant MAGICAL N° 669204 and J.C. acknowledges LabEx Minos ANR-10-LABX-55-01 for fellowship.

 34In summary, we demonstrated that W on top of FeCoB prohibits Fe out-diffusion, which occurs in case of Ta cap. The resulting interfacial compositional integrity yields a 17-29 % higher interfacial anisotropy of the storage layer with W2/Ta1 nm cap than with Ta 1 nm cap.

Transport properties of 80 nm memory cells also exhibit 50 % enhancement in thermal stability factor (∆) as well as a significant improvement in TMR (70% for Ta cap, 120% for W/Ta) with W/Ta cap compared with Ta cap. This improvement is attributed to the higher Fe concentration at the interface with MgO as well as improvement of crystallinity of the tunnel junction upon high temperature annealing up to 425°C, which is higher than BEOL thermal budget. By fitting the variation of ∆ vs. effective FeCoB thickness, the exchange stiffness constants (Aex) of FeCoB with W/Ta cap were evaluated to be 29.3 and 26.6 pJ/m after 400°C and 425°C

annealing respectively.