A nonnegativity preserving scheme for the relaxed Cahn-Hilliard equation with single-well potential and degenerate mobility

We propose and analyze a finite element approximation of the relaxed Cahn-Hilliard equation with singular single-well potential of Lennard-Jones type and degenerate mobility that is energy stable and nonnegativity preserving. The Cahn-Hilliard model has recently been applied to model evolution and growth for living tissues: although the choices of degenerate mobility and singular potential are biologically relevant, they induce difficulties regarding the design of a numerical scheme. We propose a finite element scheme and we show that it preserves the physical bounds of the solutions thanks to an upwind approach adapted to the finite element method. Moreover, we show well-posedness, energy stability properties, and convergence of solutions of the numerical scheme. Finally, we validate our scheme by presenting numerical simulations in one and two dimensions.


Introduction
Being of fourth order, the Cahn-Hilliard equation does not fit usual softwares for finite elements. To circumvent this difficulty a relaxed version has been proposed in [32] and the presentation of a finite element numerical scheme that preserves the physical properties of the solutions is the purpose of the present work. The relaxed version of the Cahn-Hilliard (RDCH in short) equation reads t > 0, x ∈ Ω, (1.1) and is set in a regular bounded domain Ω ⊂ R d with d = 1, 2, 3. It describes the evolution in time of the (relative) volume fraction n ≡ n(t, x) of one of the two components in a binary mixture. The system is equipped with nonnegative initial data n(0, x) = n 0 (x) ∈ H 1 (Ω), 0 ≤ n 0 (x) < 1, x ∈ Ω, and with zero-flux boundary conditions on the boundary ∂Ω of Ω where ν is the unit normal vector pointing outward ∂Ω.
As pointed at the end of [32], System (1.1) admits a rewriting that can proved to be useful for numerical simulations. Indeed, in the following of this article we use the fact that defining w = n− σ γ ϕ, System (1.1) can be rewritten as System (1.1) was proposed in [32] as an approximation, in the asymptotic regime whereby the relaxation parameter σ vanishes (i.e., σ → 0), of the fourth order Cahn-Hilliard equation [12,13]. The Cahn-Hilliard (CH) equation describes spinodal decomposition phenomena occurring in binary alloys after quenching: an initially uniform mixed distribution of the alloy undergoes phase separation and a two-phase inhomogeneous structure arises. In its original form, the Cahn-Hilliard equation is written in the form of an evolution equation for n:  The homogeneous free energy ψ describes repulsive and attractive interactions between the two components of the mixture while the regularizing term γ 2 |∇n| 2 accounts for partial mixing between the pure phases, leading to a diffuse interface separating the states n ≡ 0 and n ≡ 1, of thickness proportional to √ γ. The parameter γ > 0 is related to the surface tension at the interface (see, e.g., [30]) and the function b is called mobility.
In most of the literature, ψ is a double-well logarithmic potential, often approximated by a smooth polynomial function, with minimums located at the two attraction points that represent pure phases (see, e.g., [18,23,21]). The mobility can be either constant [23,21] or degenerate at the pure phases [6,22]. We refer to the introductory chapters [20,31] and to the recent review [30] for an overview of the derivation of the Cahn-Hilliard equation, its analytical properties and its variants.
Recently, the Cahn-Hilliard equation has been considered as a phenomenological model for the description of cancer growth; see, for instance, [3,17,37]. In this context, n represents the volume fraction of the tumor in a two-phase mixture containing cancerous cells and a liquid phase, such as water and other nutrients. In biological contexts, a double-well potential appears to be nonphysical. In fact, as suggested by Byrne and Preziosi in [11], a single-well potential of Lennard-Jones type allows for a more suitable description of attractive and repulsive forces acting in the mixture. Following this intuition and building upon previous works [2,17], in this paper we consider a single-well homogeneous free energy ψ : [0, 1) → R, that can be decomposed in a convex and a concave part ψ ± such that ψ(n) = ψ + (n) + ψ − (n), ±ψ ± (n) ≥ 0, 0 ≤ n < 1. (1.4) Additionally, we consider a singularity at n = 1 to represent saturation by one phase (see [11]). Hence, the potential is called single-well logarithmic and the singularity is contained in the convex part of the potential. Furthermore, we assume that and we extend the smooth concave part defined on [0, 1] to R with The latter assumption is necessary to bound the energy of the system from below. In particular, the example of potential that we have in mind is, for n ∈ [0, 1), where in this case ψ + is convex if n * ≤ 0.7. In the above form, ψ models cell-cell attraction at small densities (ψ (·) < 0 for 0 < n ≤ n and ψ (0) = 0) and repulsion in overcrowded zones (ψ (·) > 0 for n ≥ n ); cf. Figure 1. The quantity n > 0 represents the value of the cellular density at which repulsive and attractive forces are at equilibrium. With a potential of the form (1.7), the pure phases are represented by the states n = 0 and n = 1, where n = 1 is a singularity for ψ in such a way to avoid overcrowding. In this work, we consider a degenerate mobility b ∈ C 1 ([0, 1]; R + ), such that b(0) = b(1) = 0, b(n) > 0 for 0 < n < 1. (1.8) Furthermore, the admissible mobility functions that we consider admit the decomposition b(n) = b 1 (n)b 2 (n), with the extension on R defined by b 1 (n) > 0, for n > 0, b 1 (n) = 0, for n ≤ 0, (1.9) and b 2 (n) > 0, for n < 1, b 2 (n) = 0, for n ≥ 1. (1.10) The typical expression in the applications we have in mind is b(n) = n(1 − n) 2 .
(1.11) Therefore, we can easily see that, considered as transport equations, System (1.1) imposes formally the property 0 ≤ n ≤ 1. We also assume an additional property that there is some cancellation at 1 such that b(·)ψ + (·) ∈ C(R; R + ). (1.12) For the examples of mobility and potential we described above this assumption is satisfied.
The Cahn-Hilliard equation (1.3) with the logarithmic single-well potential defined in (1.7) and a mobility given by (1.11) has been studied by Agosti et al. in [2], where the authors prove wellposedness of the equation for d ≤ 3. For the relaxed version of the Cahn-Hilliard equation (1.1), well-posedness and convergence to the original Cahn-Hilliard equation as σ → 0 have been proved in [32].
We also mention here some important variants of the Cahn-Hilliard equation that have been used in the context of the modelling of tumors and in which the potential of interaction considered is of double-well form. In particular, Garcke et al. [29] proposed a Cahn-Hilliard-Darcy system that models tumor growth and in which the cells can move due to chemotaxis, attractive and repulsive forces. This model also comprises the effect that cells are crawling in a porous medium (i.e. the extracellular matrix). This effect is represented by a velocity term defined by Darcy's law. The tumor cells can also proliferate and die depending on the availability of nutrients. Garcke et al. [28] extended the model to the multiphase case to study the effect of necrosis. In order to take into account the viscosity of the tissue, and since it is not a valid approximation to consider tissues as a porous medium, Ebenbeck et al. [19] proposed a Cahn-Hilliard-Brinkman model in which Brinkman's law gives the flow velocity this time.
Summary of previous results and specific difficulties. Numerous numerical methods have been developed to solve the Cahn-Hilliard equation (1.3) with smooth and/or logarithmic double-well potential as well as with constant or degenerate mobility. Generally, a numerical scheme for the Cahn-Hilliard equation is evaluated by several aspects: i) its capacity to keep the energy dissipation (energy stability) and the physical bounds of the solutions; ii) if it is convergent, and if error bounds can be established; iii) its efficiency; iv) its implementation simplicity. To meet the first point concerning the energy stability, several implicit schemes have been proposed. The main drawback of these methods is the necessity to use an iterative method to solve the resulting nonlinear system. To circumvent this issue, unconditionally energy-stable schemes have been proposed based on the splitting of the potential in a convex and a non-convex part. This idea comes from Eyre [26] and leads to unconditionally energy-stable explicit-implicit (i.e. semi-implicit) approximations of the model. For references on all the previous numerical methods discussed above, we refer the reader to the review paper [35].
For finite element approximations, most of these results are based on the second-order splitting where, µ is called chemical potential; see, e.g., [2,6,21]. In [23], Elliot and Songmu propose a finite element Galerkin approximation for the resolution of (1.3) with a smooth double-well potential and constant mobility. The more challenging case of a degenerate mobility and singular potentials has been considered by Barrett et al. in [6], where the authors propose a finite element approximation which employs the second-order splitting (1.13 and a degenerate mobility of the form (1.11). As the authors remark, the main issues arising when considering a single-well logarithmic potential is that the positivity of the solution is not ensured at the discrete level, since the mobility degeneracy set {0; 1} does not coincide with the singularity set of the potential, i.e., n = 1. Therefore, the absence of cells represents an unstable equilibrium of the potential. In [2], the authors design a finite element scheme which preserves positivity by the means of a discrete variational inequality, as also suggested in [6]. More recently, in [ x ψ(x) Figure 1: Single-well potential of Lennard-Jones type as in (1.7) with n = 0.6. a discontinuous Galerkin finite element discretization of the equation where, again, the positivity of the discrete solution is ensured thanks to a discrete variational inequality.
In this paper, we take advantage of the equivalent rewriting (1.2) of the RDCH system and use previous results obtained on the analysis of finite element schemes derived for degenerate parabolic equations. Namely, we use results obtained by Cancès and Guichard [14] for a Control-Volume-Finite-Element (CVFE) used to simulate the anisotropic porous medium equation. In this work, the nonnegativity of the discrete solution is achieved by a suitable definition of the convective mobility that is discretized on control volumes. The compactness of time and space translates are achieved using discrete energy estimates, and the convergence analysis is obtained using the Frechet-Kolmogorov theorem. This work is the basis of a subsequent paper by Cancès et al. [15] in which a CVFE scheme for an anisotropic Keller-Segel system is analyzed. Therefore, using the fact that the RDCH system is close to the Keller-Segel model (i.e. they are both composed of one degenerate parabolic equation and one elliptic equation), we combine the results of the previously presented works on CVFE schemes to perform our analysis.
Contents of the paper. The aim of this paper is to present and analyze two finite element approximations of the relaxed Cahn-Hilliard equation with single-well potential (1.7) and degenerate mobility (1.11) in dimensions d = 1, 2, 3. More in detail, we propose two different time discretizations leading to one nonlinear scheme on which we can prove analytically some important properties and one efficient linear scheme that retrieve these important properties during numerical simulations (but we can not prove them analytically). However, for the nonlinear scheme, we prove analytically: (i) well-posedness of the numerical approximation; (ii) nonnegativity of discrete solutions ensured by adapting the upwind technique to the finite element approximation method; (iii) dissipation of a discrete energy; (iv) convergence of discrete solutions.
The main novelty of our work is to propose an alternative to the second-order splitting (1.13) by replacing the chemical potential µ by its relaxed approximation ϕ, solution to a second order elliptic equation with diffusivity 0 < σ γ. The relaxed system is based on the analysis performed in [32], where the authors prove well-posedness of the system as well as the convergence, as σ → 0, of weak solutions of (1.1) to the ones of the original Cahn-Hilliard equation (1.3). For the analysis that follows, it is worth noticing that System (1.1) admits the energy functional that, as proved in [32], is decreasing in time, i.e., We also notice that the convex/concave splitting of ψ is different from the one employed, e.g., in [2] and is motivated by the need to retrieve energy dissipation as well as by the fact that we can take advantage of the linearity of ψ − to achieve regularity results on w. Furthermore, we observe that the relaxed Cahn-Hilliard system bears some similarities with the Keller-Segel model with additional cross diffusion, proposed and analyzed in [7,16].
In this work, we use another definition of the continuous solutions for the equivalent system (1.2). To do so, we define the non-decreasing continuous function η : R → R by Similarly, we also define the non-decreasing continuous function ζ : R → R, using the properties of b(s)ψ + (s) stated in (1.12), such that With the definition of these new functions, we define the solutions of the RDCH system (1.2).
are the weak solutions of the relaxed-degenerate Cahn-Hilliard model (1.2) in the following sense: For all test functions χ 1 , χ 2 ∈ C ∞ c (Ω T , R + ) with χ 1 (T, ·) = χ 2 (T, ·) = 0, we have This paper is organized as follows. Section 2 presents the details of the Control-Volume-Finite-Element framework and the assumptions on the mesh. In Section 3, we introduce a nonlinear semi-implicit finite element approximation of the relaxed Cahn-Hilliard equation (1.1). The definition of the upwind mobility coefficient is given using a constant approximation of the mobility term on specific control volumes. Subsection 3.2 and Subsection 3.3 are dedicated to the presentation of the properties of this nonlinear scheme, such as the nonnegativity of the discrete solutions, mass conservation, and apriori estimates. The existence of discrete solutions to this problem is given in Subsection 3.4. In Subsection 3.5 and Subsection 3.6, we derive compactness estimates and prove the convergence of the discrete solutions to the weak solutions of the continuous relaxed Cahn-Hilliard model defined in Definition 1.1. Then, Section 4 is devoted to the description of an efficient linear semi-implicit scheme. The existence and the nonnegativity of a unique solution are given for a suitable upwind approximation of the mobility coefficient. Finally, in Section 5, we present numerical simulations using our linear semi-implicit scheme in one and two dimensions. These numerical simulations are in good agreement with previous numerical results obtained for the degenerate Cahn-Hilliard equation with single-well logarithmic potential.

Notations and assumptions
We first set up the notations we will use in the numerical discretization and recall some well-known properties we employ in the analysis.
Geometric and functional setting. Let Ω ⊂ R d with d = 1, 2, 3 be a polyhedral domain. We indicate the usual Lebesgue and Sobolev spaces by respectively L p (Ω), W m,p (Ω) with H m (Ω) := W m,2 (Ω), where 1 ≤ p ≤ +∞ and m ∈ N. We denote the corresponding norms by || · || m,p,Ω , || · || m,Ω and seminorms by | · | m,p,Ω , | · | m,Ω . The standard L 2 inner product will be denoted by (·, ·) Ω and the duality pairing between (H 1 (Ω)) and H 1 (Ω) by < ·, · > Ω . Let T h , h > 0 be a conformal mesh on the domain Ω which is defined by N el disjoint piecewise linear mesh elements, denoted by K ∈ T h , such that Ω = K∈T h K. The elements are triangles for d = 2 and tetrahedra for d = 3. We let h := max K h K refers to the level of refinement of the mesh, where h K := diam(K) for K ∈ T h . We define by κ K the minimal perpendicular length of K and κ h = min K∈T h κ K . We assume that the mesh is quasi-uniform, i.e., it is shape-regular and there exists a constant C > 0 such that Moreover, we assume that the mesh is acute, i.e., for d = 2 the angles of the triangles do not exceed π 2 and for d = 3 the angle between two faces of the same tetrahedron do not exceed π 2 . We define the set of nodal points J h = {x j } j=1,...,N h , where the number of nodes is N h := card(J h ), and we assume that each x j is a vertex of a simplex K ∈ T h . We denote by Λ i the set of nodes connected to the node x i by an edge and G h the maximal number of connected nodes, i.e.
For two nodes x i , x j , if they are connected by an edge, we denote this latter e ij .
We also define the barycentric dual mesh associated to T h . On each element K ∈ T h , the barycentric coordinates of an arbitrary point x ∈ K are defined by the real numbers λ i with i = 1, . . . , n K such that where n K is the number of nodes of the element K. We define the barycentric subdomains associated to the vertex P i ∈ K k (where K k refers to the k-th element of T h and k = 1, . . . , N el ), by Therefore, for each node of the mesh T h , we define the barycentric cell of the dual mesh We also add another subdivision of the mesh that relies on the definition of the barycentric cells and that will be useful to define the upwind approximation of the mobility. Using the barycentric coordinates, we subdivide each element K ∈ T h in d+1 subdomains, such that, for i = 1, 2, 3, j = 2, 3, and i = j, we have the subdomaiñ In the following of the article, we use the terminology Diamond cell for the twoD K ij that share the edge e ij . To illustrate what are these subdomains, we represent graphically what isD T ij on Figure 2 as well as what is a diamond cell for d = 2.
We introduce the set of piecewise linear functions χ j ∈ C(Ω) associated with the nodal point x j ∈ J h , that satisfies χ j (x i ) = δ ij , where δ ij is the Kronecker's delta. We introduce the P 1 conformal finite element space V h associated with T h , where P 1 (K) denotes the space of polynomials of order 1 on K: Furthermore, we let K h be the subset containing the nonnegative elements of V h , namely We denote by π h : C(Ω) → V h the Lagrange interpolation operator corresponding to the discretized domain T h , defined as We also use the lumped spaceV h defined bŷ On C(Ω) we define the approximate scalar product as where x i ∈ K denotes the vertices of the element K. We denote the corresponding discrete semi-norm Inequalities. We summarize important inequalities that will be used later on for the analysis of the numerical schemes. We will use the following inequalities (see, e.g., [33]): (2.1) Concerning the interpolation operator, the following result holds [9]: and we have [34], For the L 2 projection operator, the following inequality holds [9] v − P 0 Finite element matrices. We define M and Q respectively the mass and stiffness matrix. M l is the lumped mass matrix, that is the diagonal matrix where each coefficient is the sum of the associated We recall some important properties of the stiffness matrix. We start by Furthermore, from the fact that the triangulation is of acute type, we know from [27] that 3 Nonlinear CVFE scheme

Description of the scheme
Given N T ∈ N, let ∆t := T /N T be the constant time step size and t k := k∆t, for k = 0, . . . , N T − 1.
We consider a partitioning of the time Using the finite element approximations of n and w, where {n k j } j=1,...,N h and {w k j } j=1,...,N h are the unknown degrees of freedom, we introduce the following finite element approximation of System (1.1).
whereB(·) is the discrete upwind approximation of the continuous mobility b(·) for the convective term whileB(·) and Bψ + (·) are the constant approximations of the mobility and second derivative of the convex part of the potential multiplied by the mobility. These latter quantities are defined in Equations (3.11), (3.12), and (3.13) respectively. The initial condition for the discrete problem is chosen such that for n 0 (x) ∈ H 1 (Ω), we have and w 0 h is the unique solution of To preserve the nonnegativity of the discrete solutions, we compute the mobility coefficient employing an implicit upwind method adapted to the finite element setting. The explanation on how to adapt the upwind method requires us to define the matrix formulation of the problem (3.1)-(3.2).
Matrix form. For k = 0, . . . , N T − 1, let n k and w k be the vectors We can then rewrite System (3.1)-(3.2) in its matrix form where ψ − is the vector containing the values at the nodal values In the above matrix form, we have used the definitions of the finite element matrices associated with and From (2.6) and (2.7), we can interpret the first equation of System (3.1)-(3.2) as using the finite volume method i.e. , for every . Thus, the scheme reads Let us now describe in detail how the coefficientsB k+1 ij ,B k+1 ij and Bψ + k+1 ij are computed.
Constant approximation of mobility and second derivative of potential. We define the upwind mobility coefficient associated with the convective term bỹ The mobility function and second derivative of the convex part of the potential associated to the diffusion term are defined byB k+1 ij and The approximation (3.11) is similar to the one used in [4] for the one-dimensional finite volume discretization of the Keller-Segel system. Furthermore, our adaptation of the upwind method in the finite element context is also close in spirit to the one proposed by Baba and Tabata in [5], where the authors also used barycentric coordinates to define the basis functions.
We remark that the coefficientsB k+1 ij ,B k+1 ij and Bψ + k+1 ij are constant and uniquely defined along each edge of the mesh (i.e. in the diamond cells). Yet, this method is well suited for a standard assembling procedure and, as a result, is simpler to implement in already existing finite element software since it requires only the adaptation of the calculation of a non-constant matrix. This method can also be adapted for the simulation of other advection-diffusion equations to preserve the nonnegativity of solutions.
From this approximation of the mobility we have the following properties concerning the convective flux.
Proposition 3.1 (Properties of the convective flux). For any (a, b, c) ∈ R 3 , we have: 5. there is a modulus of continuity ω : 3.2 Discrete maximum principle, non-negativity, and mass conservation Proof.
Step 1: Non-negativity. We prove the non-negativity by induction. We assume that From the non-negativity ofB k+1 ij , and Bψ + k+1 ij as well as (2.7), we know that , and from the extension of b 1 (s) by zero for s < 0, Step 2: Upper bound. To prove the upper bound we repeat the previous argument but this time we assume ∀x i ∈ J h , 0 ≤ n k i ≤ 1, and n k+1 ) and from the non-negativity ofB k+1 ij and Bψ + k+1 ij as well as (2.7), we have  Proof. To prove mass conservation, we use the identity, for each (3.14) Summing over the nodes in Equation (3.1), we obtain Using the symmetry of the matrices U k+1 and R k+1 , the property (3.14), we obtain which implies mass conservation.

Energy stability and a priori estimates
The finite element numerical scheme (3.1)-(3.2) preserves the dissipation of the energy at the discrete level.
Proof. Multiplying Equation (3.9) by n k+1 i and summing over the nodes x i ∈ J h leads to The term on the right hand side can be bounded using Young's Inequality, and the boundedness of the mobility to obtain, for any κ > 0, ∆t Therefore, it leads to Therefore, from the boundedness of bothB k+1 ij andB k+1 ij , as well as the fact thatB k+1 ij ≥B k+1 ij by definition, and since κ can be chosen arbitrarily ,we know that it exists a finite constant (we recall that the term γ σ does not induce any difficulty since σ < γ). Using the property (a − b)a ≥ 1 2 (a 2 − b 2 ), Proposition 3.5 as well as the property (2.2), and summing for k = 0 → N T − 1, we obtain the result. Proof. From Proposition 3.6 and the non-negativity of Bψ + k+1 ij , we know that Furthermore, from the definition (3.12) and the definition of η provided by (1.15), we have The same is obtained for N T −1 k=0 ∆t ζ k+1 h 2 1 using similar arguments.

Existence of discrete solution
We prove well-posedness of our problem.
Proof. The proof of the existence of a solution for the discrete problem relies on the use of Brouwer's fixed point theorem. A necessary step before applying the theorem is to change the unknown of our problem, we define m k+1 h = n k+1 h − α, where α = 1 |Ω| Ω n 0 dx. Therefore, the discrete problem (3.1)- where F : S → S is an application defined on the space S with S = {m | M l m · (1, . . . , 1) = 0 and − α ≤ m ≤ 1 − α}.
The first constraint in the definition of the space S reflects the conservation of the initial mass while the second comes from Inequality (3.17). As a result, S is a convex and compact subspace of R N h . Using the matrix formulation (3.5)-(3.6) for the problem P , the application F is defined by In the previous definition of F , we precise that the matrix U k+1 associated to the nonlinear convection term is given by and for the diffusion term To check if F is continuous, we compute for m 1 , m 2 ∈ S, F (m 1 ) − F (m 2 ) , and obtain For the first term on the right-hand side of the previous inequality, we have from the properties of the mesh, the upper bound presented by Varah [36] for the inverse of M-matrices, and our assumptions on n k h and w k h as well as the function ψ − (·). Lastly, from the definitions of the matrix U (·) and of the convective flux G (see Definition 3.1), we have The same applies using the properties of the matrix R(·) to obtain Altogether, using the properties of the standard finite element matrices and the fact that the mobilitỹ B(m k h,ε + α) is bounded, we obtain which proves that the mapping F is Lipschitz continuous. Therefore, applying Brouwer's fixed point theorem, we know that it exists a solution m k+1 ∈ S such that F (m k+1 ) = m k+1 . Therefore, it exists a m k+1 h ∈ V h that gives the existence of a pair {n k+1 h , w k+1 h } ∈ V h × V h solution of the Problem P .

Compactness estimates
We now derive bounds for the time and space translates of the discrete solutions. This results follow the lines of the works [14,15]. We first define the space and time discrete spaces V h ∆t andV h ∆t . These sets are composed of piecewise constant functions in time with values in V h andV h respectively.
To study the convergence as h, ∆t → 0, we use an index m such that, as m → ∞, h m , ∆t m → 0. Therefore, we study the convergence of sequences of functions in the spaces V hm ∆tm andV hm ∆tm .
We define by η hm,∆tm , ζ hm,∆tm and w hm,∆tm the piecewise affine in space and piecewise constant in time approximation of the functions η, ζ and w. We also denote byη hm,∆tm ,ζ hm,∆tm andŵ hm,∆tm , the piecewise constant in space and time corresponding approximations. For each node x i ∈ J h and time t = t k , we have the notation η hm,∆tm (x i , t k ) =η hm,∆tm (x i , t k ) = η k i . The associated vector containing all the value for all nodes is denoted by η k .
Time translate estimates We start by estimates on the time translates. We denote by Q T −τ = Ω × (0, T − τ ), for all τ ∈ (0, T ). We have the following result Lemma 3.9 (Time translate estimates). There are constants C η,t , C ζ,t and C w,t independent of h and τ such that the following inequalities hold Proof. The proof is similar to the proof of Lemma 4.3 in [14]. We give here the details for the sake of clarity since the equation under study is different.
Since the functionsη hm,∆tm are piecewise constant in time, we define υ(t), for t ∈ (0, T ], the unique integer such that t υ(t) < t ≤ t υ(t)+1 . Therefore, we write ∀t ∈ (0, T − τ ), However, from definition (1.15), we know that it exists a constant C such that Then, using the definition of the integer υ(·), we find From the first equation of the scheme (3.1), we obtain Using Young's inequality, we have To handle the first sum in each of these quantities, we introduce some additional notations. We define ρ(k, t, τ ) the characteristic function such that ρ(k, t, τ ) = 1, if t < k∆t m ≤ t + τ, 0, otherwise.
Using the previous argument and Corollary 3.7, we obtain This achieves the proof of Inequality (3.18). The proofs of Inequalities (3.19) and (3.20) are very similar and we do not repeat them.
Space translate estimates We now turn to the space translate estimates. Lemma 3.11 (Space translate estimates). There are three constants C η,s , C ζ,s and C w,s independent of m and y such that Proof. The proof of this result follows Lemmas 4.1 and 4.2 in [14]. For the sake of clarity, we here present the main steps.
First, using Corollary 3.7, we know that ∇η hm,∆tm (L 2 (Ω×(0,T ))) d is bounded. Since the time and space domain considered is of finite measure, Hölder inequality indicates that ∇η hm,∆tm (L 1 (Ω×(0,T ))) d is also bounded. Hence, since η hm,∆tm is bounded as well, its extension by zeros outside Ω × (0, T ) lies in L ∞ BV (R d+1 ). Therefore, we have from which, using the triangular inequality and Remark 3.12 stated below, we find Inequality (3.21). The same is found for Inequality (3.22) and Inequality (3.23) following the same arguments.
Remark 3.12. One of the important tool that we did not present in the previous proof is found from the use of Lemma A.2 in [14].

Convergence analysis
This section is organized as follow, we first apply Fréchet-Kolmogorov theorem to show strong convergence in L 1 (Ω T ), then we show that the limit is a solution of the continuous RDCH system.
To prove the convergence (3.27), we start by stating that the inverse η −1 of the continuous function η exists and is continuous. Furthermore, we have thatn hm,∆tm is bounded in L ∞ from Proposition 3.2. Therefore, applying the dominated convergence theorem ton hm,∆tm = η −1 (η hm,∆tm ), we arrive to the convergence (3.27) in which the limit is defined as the unique function n(t, x) = η −1 (η) (η being the limit function in (3.25)).
Proof. The proof of this theorem is very close to Section 5 in [15]. For the sake of clarity, our proof can be found in Appendix A.
For each k = 0, . . . , We define the following finite elements matrices and We write the matrix form of Equation (4.1a) M l + ∆tL k n k+1 = −∆tU k ϕ k+1 + M l n k , and since U k has zero row sum, we can rewrite the previous equation for each node where Λ i is the set of nodes connected to the node x i by an edge. In the definition of (4.2) we compute the mobility coefficient in function of the direction of ∇ϕ k+1 h . As for the nonlinear case, the mobility coefficient is given byB otherwise. Even though we cannot redo the same analysis as for the nonlinear scheme, we can establish the existence and the nonnegativity of discrete solutions of (4.1a) and (4.1b). Let Ω ⊂ R d , d = 1, 2, 3, and assume that T h is a quasi-uniform acute mesh of Ω, and the condition

5)
(where Λ i is the set of nodes connected to the node x i by an edge) is satisfied. Then, the linear finite element scheme (4.1a)-(4.1b) with initial condition n 0 h ∈ K h admits a unique solution Proof.
Step 1. Existence of a unique solution. Assuming that {n k h , ϕ k h } ∈ K h × V h , from the Lax-Milgram theorem, it exists a unique solution ϕ k+1 h ∈ V h of Equation (4.1b) and Equation (4.1a) admits a unique solution n k+1 h ∈ V h . Therefore, it exists a unique pair of discrete solutions {n k+1 h , ϕ k+1 h } ∈ V h × V h for the system (4.1a)-(4.1b). Next, we need to prove that n k+1 h is nonnegative and bounded from above by 1.
Step 2. Nonnegativity and upper bound on n k+1 h for d = 1, 2, 3. First, from the fact that (M l + ∆tL k ) is a M-matrix, we know that its inverse is non-negative, i.e. (M l +∆tL k ) −1 ≥ 0. Therefore, to preserve the non-negativity of n k+1 h , we need that For every node x i ∈ J h , the previous condition reads where Λ i is the set of nodes connected to the node x i by an edge. From the fact that the mesh is acute, we know that Q ij is negative. Therefore, using the definition of the mobility coefficient (3.11), we need to focus on the case ϕ k+1 j − ϕ k+1 i < 0. In that situation, we have However, from (2.7), we find the following condition to ensure the non-negativity of n k+1 Then, we need to prove that for every node x i ∈ J h we have n k+1 i < 1. We use the upper bound presented by Varah [36] for the inverse of M-matrices, and write Therefore, to retrieve the upper bound on the discrete solution, the condition has to be satisfied. Note in the previous equation that we have considered the case ϕ k+1 j − ϕ k+1 i > 0 since in the other case the bound will be satisfied trivially. Then, subtracting n k i from both sides of the previous inequality, we obtain and we retrieve the same condition as before with a strict inequality.
Altogether, we proved the existence of a unique solution {n k+1 h , ϕ k h } ∈ K h × S h for the system (4.1a)-(4.1b) with 0 ≤ n k+1 h < 1 if the stability condition (4.5) is satisfied.

Numerical simulations
In this section, we use the previously presented linear scheme (4.1a)-(4.1b) for the RDCH system. The numerical simulations are performed using the MATLAB software. At each time step the matrices U k and L k are reassembled. The linear system is solved using the function linsolve of the MATLAB software. This function uses the LU factorization.
Even though we are presenting numerical results obtained using the linear scheme (4.1a)-(4.1b), the evolution of the energy during the simulations is given from the computation of the discrete formulation of the continuous energy First of all, we present test cases in one and two dimensions to validate our method. The physical properties of the solutions such as the shape of the aggregates, the energy decay, the mass preservation and the non-negativity of the solution are the key characteristics we need to observe to validate our method. A comparison with previous results from the literature is also of main importance. The reference used for this study is the work of Agosti et al. [2]. The analysis of the long-time behavior of the solutions of the RDCH equation [32] gives us some insights about what we should observe at the end of the simulations. The solutions should evolve to steady-states that are minimizers of the energy functional. Depending on the initial mass, three regions of the cell density should appear. The first being the region of absence of cells, the second the continuous interface linking the bottom and the top of the aggregates. If the initial mass is sufficiently large, the third expected region is a plateau of the cell density close to n = n . This evolution is related to the clustering of tumor cells in in-vitro biological experiments (see e.g. [10,3]).
The study of the effect of the regularization on the numerical scheme is the purpose of the last subsection.

Numerical results: test cases
1D test cases Table 1 summarizes the parameters used for the one dimensional test cases. The initial cell density is a uniform distributed random perturbation around the values n 0 . Figures 3 show the evolution in time of the solutions n h for the three different initial masses.
We can observe that the solution for each of the three test cases remains nonnegative and the mass is conserved throughout the simulations. From Figure 4, we observe that the energies decrease monotonically for the three simulations but at different speeds. They all display at the end of the computation a stable (or metastable) state that is a global (or respectively a local) minimizer of the discrete energy.
For the initial condition n 0 = 0.3 (Figures 3 b1), b2), b3) and Figure 4 at the middle), the energy decreases rapidly and reaches a plateau showing that the solution evolves rapidly to a steady state. The solution at t = 10 presents aggregates that are not saturated (i.e. the maximum density is below n ). The explanation behind this observation is that the initial mass is not sufficiently large for the system to produce any saturated aggregates. However, the clusters appear to be of similar thickness and are relatively symmetrically distributed in the domain.
For n 0 = 0.36 > n /2 (Figures 3 c1), c2), c3)), the aggregates are thicker. The top of the aggregate located at the center of the domain is flatter than for the other simulation. The maximum density is closer to the value n than for the initial condition n 0 = 0.3. Likewise, the symmetry in the domain is respected. Using Figure 4 on the right, we observe that at different times, the energy evolves through several meta-stable equilibria. This reflects the fact that the solution went to different meta-stable states before reaching a stable equilibrium that better minimizes the energy.  For the initial condition n 0 = 0.05 (Figures 3 a1), a2), a3)), the shape of the final solution is different. The aggregates appear to be thinner and far from each other. The symmetry is not retrieved in the domain. Furthermore, from Figure 4 (on the left), we can observe that the evolution of the solution is slow compared to the other initial conditions. The energy seems to be constant in the first moment of the simulation (i.e. after the spinodal decomposition phase). The slow evolution of the solution is explained from the fact that the mobility is degenerate and the amount of mass available in the domain is small. Using Figure 4, we can also see that the energy continues to decrease even at the end of the simulation. To keep comparable simulation times, we did not reach the complete steady state.
Let us compare qualitatively these results with the ones obtained in [2] for the one dimensional case. For the two test cases n 0 = 0.3 and n 0 = 0.36, there is no differences in the shape the aggregates or the distribution of the mass in the domain. For n 0 = 0.05, some small discrepancies with the final solutions are observed. In particular, the symmetry of the aggregates in the domain is not respected in our case whereas it is in the reference work. We must stress that doing other simulations, the symmetry was sometimes reached at the time t ≈ 100 for the initial condition n 0 = 0.05. The reason is that the system will evolve to respect the symmetry but the time at which this stable-steady state is reached depends on the initial distribution of the cell density.
Altogether, the solutions obtained at the end the three simulations are in accordance with the description of the steady-states made in [32]. The three regions of interest are indeed retrieved at the end of each simulation.

2D test cases
For the two-dimensional test cases, the domain is a square of length L = 1. The initial density is computed in the same way as for the one-dimensional test cases i.e. a random uniformly distributed perturbation around n 0 . The summary of the values of parameters can be found in Table 2. Figure 5 depicts the results of the three test cases with different initial masses. The three simulations respect the nonnegativity of the cell density, the conservation of the initial mass and the monotonic decay of the discrete energy. However, different shapes can be observed for the aggregates. Figures 5 a1),a2),a3) show the evolution of the solution through time for the small initial mass n 0 = 0.05. Starting from a uniform random distribution of the cell density in the domain, the solution evolves into a more organized configuration. Progressively, a separation of the two phases of the mixture occurs. At the end of the simulation, small clusters are formed. They display a circular shape and are of similar width. The organization of the clusters in the domain tries to maximize the distance between each other. Using the Figure 6 (left), we observe a drop of the energy in the first moments of the simulation denoting a fast reorganization of the randomly distributed initial condition. Then, the solution appears to evolve very slowly, i.e. a meta-stable state was reached. A second drop of the energy follows around t ≈ 15, the system enters the "coarsening" phase: the small aggregates become more dense and merge with others. At the end, the evolution is very slow. The system continues to rearrange but due to the degeneracy of the mobility and the small amount of initial mass this process is very slow. Figures 5 b1),b2),b3) show the evolution of the solution for n 0 = 0.3. The two phases that are the spinodal decomposition and the coarsening are retrieved. Between the Figures 5 b1) and 5 b2), we observe that the solution evolves from a random uniform configuration to an organization in small aggregates that are not saturated. Then ( Figure 5 b3)), the cell density is distributed in elongated and saturated aggregates. The separation of the two phases is clear. However, using Figure 6 (middle), we observe that at the end of the simulation the cell density continues to rearrange. Due to the degeneracy of the mobility, this evolution is very slow.
In Figures 5 c1),c2),c3), we can observe the evolution of the solution for n 0 = 0.36. Again, the solution goes through the spinodal decomposition and coarsening phases. The only difference that needs to be highlighted for this simulation is the different shape of the aggregates at the end. Indeed, the initial mass being n 0 = 0.36 > n /2, the aggregates are wider and more connected to each others.
Therefore, depending on the initial mass of cells in the domain, the 2D simulations of the model show very different spatial organizations of the cell density.  Compared to the reference work [2], the organizations of the cells for the different initial cell densities are the same. No clear difference can be established regarding the simulation involving the relaxed model and the original one.
The three regions corresponding to a steady-state described in [32] are retrieved at the end of the simulations for these 2D test cases.

Effect of the relaxation parameter σ
In this section we evaluate the effect of the relaxation parameter σ for the stability of the scheme, and in particular to satisfy the CFL-like condition (4.5). This conditions is necessary to preserve the nonnegativity of the solutions of the linear discrete scheme. To evaluate the effect of this parameter on the choice of the time step ∆t, we compute the amplification matrix H defined by Here, X k is called the state vector. Using the matrix form of the scheme (4.1a)-(4.1b) we can decomposed the amplification matrix by H = H −1 1 H 2 with We denote by λ i , i = 1, . . . , N h , the eigenvalues of the amplification matrix H.
To analyze the stability of the numerical scheme due to the relaxation parameter, we compute the spectral radius of the amplification matrix for a smooth initial conditions. The scheme is stable when the maximum value of the modulus of the eigenvalues is less or equal to 1. Figure 7 represents the spectral radius in function of the time step ∆t for two values of σ (the other parameters are the ones taken from the one dimensional test cases with n 0 = 0.3). We can observe that the scheme remains stable when ∆t is small for the two test cases. However, we see that increasing σ allows to take larger time steps while remaining stable. This result can be explained due to the fact that increasing σ diminishes the value max present in the stability condition (4.5). Therefore, the regularization induced by the relaxation parameter allows for faster simulations, but it has an effect on the accuracy of the solution compared to the solution given by the non-relaxed model. However, at the moment it remains unclear how to compare the solution given by a simulation of the relaxed model and a solution of the original degenerate model (without relaxation). This will be the subject of a further work.

Conclusion
We described and studied a finite element method to solve the relaxed degenerate Cahn-Hilliard equation with single-well logarithmic potential. We considered two time discretization schemes: nonlinear semi-implicit and linear semi-implicit. We showed that the nonlinear scheme preserves the physical properties of the solutions of the continuous model. We proved that this scheme is well-posed and convergent in dimension d = 1, 2, 3. The nonnegativity of the solutions is retrieved thanks to the use of an upwind method adapted within the finite element framework. The linear semi-implicit scheme allows faster simulations and we proved that it is well-posed and preserves the nonnegativity of the solutions as well. We presented some numerical simulations using this linear scheme in one and two dimensions. The numerical simulations validated the nonnegativity-preserving and energy decaying properties of the scheme. The numerical solutions of the finite element approximation of the RDCH model are in good agreement with previous works dealing with the non-relaxed model. We showed that the relaxation parameter σ allows us to take a larger time step in the scheme (as long as the condition for the nonnegativity is preserved and σ < γ). We point out that thanks to the spatial relaxation, our numerical scheme can be easily implemented and simulations of the relaxed degenerate Cahn-Hilliard model can be computed efficiently using standard softwares.
First, we choose a test function χ ∈ C ∞ c (Ω T , R + ) with χ(T, ·) = 0. Then, we multiply Equation Next, we show the limit of each quantity when m → ∞.
For B m , since χ N T h = 0, we have (see [24]) For C m , we aim to show that C m → 0 as m → ∞. To do so, we rewrite From the Cauchy-Schwarz Inequality, we have The first term on the right-hand side is bounded from (3.15). The second term, that we denote by R m , is handled using Lemma A.1 in [14]. Indeed, we denote for each K ∈ T hm , η k+1 K = max and we define the uniformly continuous function √ b • η, defined on the interval [0, η(1)], such that for each K ∈ T hm , B k+1 Therefore, we have From the regularity of χ, we know that N h i=0 x j ∈Λ i |Q ij | (χ k+1 i − χ k+1 j ) 2 is bounded and we obtain From Lemma A.1 in [14], we obtain R m → 0, as m → ∞, hence, C m → 0. Using similar arguments, we have E m → 0, as m → ∞.
However, we need to show that |D m − D m | → 0 as m → ∞. We use similar arguments as to show that C m → 0. Indeed, we rewrite with a k+1 K = Θ k+1 K , we are in position to repeat the arguments presented for the convergence of C m → 0 and we do not repeat them here. Therefore, we obtain, as m → ∞, D m → γ σ Ω T b(n)∇η(n) · ∇χ dx dt.
The same arguments are also applied to the convergence of F m to obtain, as m → ∞, The last term of the first equation to analyze is G m . We repeat again the same arguments. We use the previously defined quantities Θ hm,∆tm and Γ hm,∆tm . Then, we define G m = γ σ Ω T (Θ hm,∆tm ) 2 ∇w hm,∆tm · ∇χ hm,∆tm dx dt, and we have, as m → ∞, We still have to prove that D m − D m → 0 as m → ∞. Since we have we obtain the result using similar arguments as for the convergence of R m . In the limit m → ∞, we obtained the first equation of the weak system defined in (1.1).
For the second equation, we obtain the limit equation using the weak convergence of w hm,∆tm given by (3.26), the strong convergence of w hm,∆tm given by (3.25), the continuity of ψ − (·), and the strong convergence of n hm,∆tm given by (3.27).
This concludes the proof of Theorem 3.14.