Jean Cauvin-Vila

Virginie Ehrlacher

AND Greta Marino

Jan-Frederik Pietschmann

STATIONARY SOLUTIONS AND LARGE TIME ASYMPTOTICS TO A CROSS-DIFFUSION-CAHN-HILLIARD SYSTEM

Keywords: 2020 Mathematics Subject Classification. 35D30, 35G31, 35G50 Cahn-Hilliard, cross-diffusion, weak solutions, global existence, degenerate Ginzburg-Landau

We study some properties of a multi-species degenerate Ginzburg-Landau energy and its relation to a cross-diffusion Cahn-Hilliard system. The model is motivated by multicomponent mixtures where crossdiffusion effects between the different species are taken into account, and where only one species does separate from the others. Using a comparison argument, we obtain strict bounds on the minimizers from which we can derive first-order optimality conditions, revealing a link with the single-species energy, and providing enough regularity to qualify the minimizers as stationary solutions of the evolution system. We also discuss convexity properties of the energy as well as long time asymptotics of the time-dependent problem. Lastly, we introduce a structure-preserving finite volume scheme for the time-dependent problem and present several numerical experiments in one and two spatial dimensions.

Introduction

The aim of this work is the study of a multi-species degenerate Ginzburg-Landau energy and its relation to a system of cross-diffusion Cahn-Hilliard equations which was recently studied in [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF]. The latter model describes the evolution of a multicomponent mixture where cross-diffusion effects between the different species are taken into account, and where only one species does separate from the others. This is motivated by multiphase systems where miscible entities may coexist in one single phase, see [START_REF] Klinkert | Comprehension and optimisation of the co-evaporation deposition of Cu (In, Ga) Se2 absorber layers for very high efficiency thin film solar cells[END_REF][START_REF] Wenisch | Nickel-enhanced graphitic ordering of carbon ad-atoms during physical vapor deposition[END_REF] for examples. Within this phase, cross-diffusion between the different species is taken into account in order to correctly account for finite size effects that may occur at high concentrations. We assume that the mixture occupies an open, smooth and bounded domain Ω ⊂ R d with d = 1, 2, 3 and that there are n + 1 species in the mixture. We denote by u i (x, t), i = 0, . . . , n, the volumic fraction of the i th species at point x ∈ Ω and time t ≥ 0 and set u = (u 0 , . . . , u n ). The dynamics of the system is governed by the free energy functional

E(u) := Ω n i=0 (u i ln u i -u i + 1) + ε 2 |∇u 0 | 2 + βu 0 (1 -u 0 ) dx, (1.1)
where ε and β are positive constants. Denoting by µ = D u E(u) the chemical potential, the corresponding evolution system formally reads as

∂ t u = div (M (u)∇µ) in Ω × (0, +∞), (1.2)
where M : R n+1 → R (n+1)×(n+1) is a degenerate mobility matrix. More precisely, for every i ̸ = j = 0, . . . , n, let K ij be positive real numbers satisfying K ij = K ji , then for u ∈ R n+1 + , it has entries M ij (u) := -K ij u i u j for all i ̸ = j = 0, . . . , n, M ii (u) := 0≤k̸ =i≤n K ik u i u k for all i = 0, . . . , n.

(1.3)

As expected, due to their interpretation as volumic fractions, the quantities u i must satisfy 0 ≤ u i (x, t) ≤ 1 for all i = 0, . . . , n and n i=0 u i (x, t) = 1 for a.e. x ∈ Ω, t ∈ (0, +∞), (1.4) and the constraint on the sum is referred to as the volume-filling constraint. The evolution system is supplemented with no-flux boundary conditions as well as initial conditions consistent with the constraints. The main result from [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF] is the existence of a solution to a suitable weak formulation of this problem. The aim of this paper is twofold. First, we study some solutions to the stationary problem 0 = div (M (u)∇µ) in Ω. (1.5) In general, the analysis of this system of coupled, degenerate elliptic equations is by no means straight forward. In this work, motivated by the gradient flow structure of the time-dependent equation highlighted above, we focus our study on the set of local minimizers of the energy functional (1.1). The latter are natural candidates for solutions to (1.5), in the sense that one naturally expects that solutions of the time-dependent system should converge in the long time limit to one of these local minimizers. We acknowledge here that other stationary solutions may exist, but stress on the fact that local energy minimizers are of particular physical relevance for the present system. When the parameters are chosen such that the energy functional is convex, the unique minimizers are constants and we show that solutions to the evolution problem (1.2) converge to them exponentially fast.

In the non-convex case, the dynamics is much more complex which leads us to the second aim of the paper: we introduce a finite volume scheme that preserves the structure of the continuous time-dependent system. The simulations demonstrate the capability of the scheme and allow to explore the dynamics for arbitrary parameter regimes.

Let us briefly review previous contributions on the respective components of our model.

Cross-diffusion systems with size exclusion. Systems of partial differential equations with cross-diffusion have gained a lot of interest in recent years [START_REF] Küfner | Invariant regions for quasilinear reaction-diffusion systems and applications to a two population model[END_REF][START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF][START_REF] Chen | Analysis of a parabolic cross-diffusion population model without self-diffusion[END_REF][START_REF] Lepoutre | Global well-posedness of a conservative relaxed cross diffusion system[END_REF][START_REF] Jüngel | The boundedness-by-entropy method for cross-diffusion systems[END_REF] and appear in many applications, for instance the modelling of population dynamics of multiple species [START_REF] Burger | Lane formation by side-stepping[END_REF] or cell sorting as well as chemotaxis-like applications [START_REF] Painter | Volume-filling and quorum-sensing in models for chemosensitive movement[END_REF][START_REF] Painter | Continuous models for cell migration in tissues and applications to cell sorting via differential chemotaxis[END_REF].

Ginzburg-Landau Energy. In the case n = 1, which implies u 0 = 1 -u 1 , (1.1) reduces to the classical Ginzburg-Landau energy with singular potential as introduced in [START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF]. The works [START_REF] Gelantalis | Energy barrier and Γ-convergence in the d-dimensional Cahn-Hilliard equation[END_REF][START_REF] Gelantalis | Existence and properties of certain critical points of the Cahn-Hilliard energy[END_REF] study the structure of energy minimizers to the functional

E GL (v) = Ω 1 2 |∇v| 2 + 1 4 (1 -v 2 ) 2 dx,
when the system size is large and the mean value of the phase parameter v is close to -1. This is motivated by nucleation phenomena, see [START_REF] Loganayaki | Growth and investigations on the nucleation kinetics of zinc succinate NLO single crystals[END_REF][START_REF] Li | Nucleation rate calculation for the phase transition of diblock copolymers under stochastic Cahn-Hilliard dynamics[END_REF], which correspond to non-constant minimizers. The authors study the case when constant stationary states are local but not global minimizers and estimate the size of the energy barrier, i.e. the difference of the energy at the respective states. In particular, the authors prove bounds on the minimizers using suitable competitors which inspired part of the construction in the proof of Theorem 2.3.

Cahn-Hilliard equation. The scalar Cahn-Hilliard equation with constant mobility was introduced in [START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF] as a model for phase separation. It is indeed the H -1 -gradient flow to (1.1) for two species. Existence of weak solutions was first shown in e.g. [START_REF] Elliott | On the Cahn-Hilliard equation[END_REF][START_REF] Caffarelli | An L ∞ bound for solutions of the Cahn-Hilliard equation[END_REF] in the case of constant mobility, and later extended to degenerate, concentration dependent mobilities [START_REF] Elliott | On the Cahn-Hilliard equation with degenerate mobility[END_REF]. Regarding the long-time behaviour, for a constant mobility and in one spatial dimension, the authors in [START_REF] Otto | Relaxation to equilibrium in the one-dimensional Cahn-Hilliard equation[END_REF] show that for initial data with bounded distance to a so-called kink state, algebraic convergence to equilibrium holds. This was further improved in [START_REF] Otto | Optimal l 1 -type relaxation rates for the Cahn-Hilliard equation on the line[END_REF]. We also refer to [START_REF] Abels | Convergence to equilibrium for the Cahn-Hilliard equation with a logarithmic free energy[END_REF][START_REF] Schimperna | Global attractors for Cahn-Hilliard equations with nonconstant mobility[END_REF][START_REF] Schimperna | Existence of solutions and separation from singularities for a class of fourth order degenerate parabolic equations[END_REF] for long-time analysis in the case of logarithmic nonlinearity. More details can be found in the review [START_REF] Novick-Cohen | The Cahn-Hilliard equation[END_REF] and the monograph [START_REF] Miranville | The Cahn-Hilliard Equation: Recent Advances and Applications[END_REF].

Multi-species Cahn-Hilliard systems have been studied in several earlier works and usually consider an energy functional of the form

E(u) := Ω Ψ(u) + 1 2 ∇u • Γ∇u dx, (1.6)
for some symmetric positive semi-definite matrix Γ ∈ R (n+1)×(n+1) and bulk free-energy functional Ψ. In [START_REF] Elliott | A generalized equation for phase separation of a multi-component mixture with interfacial free energy[END_REF],

Elliott and Luckhaus proved a global existence result for such a multiphase Cahn-Hilliard system with constant mobility and Γ = γI for some γ > 0. In [START_REF] Elliott | Diffusional phase transitions in multicomponent systems with a concentration dependent mobility matrix[END_REF], the authors generalized their result to the case of a degenerate concentration-dependent mobility matrix with a positive definite matrix Γ while [START_REF] Brunk | On existence, uniqueness and stability of solutions to Cahn-Hilliard/Allen-Cahn systems with cross-kinetic coupling[END_REF] study a system of Cahn-Hilliard/Allen-Cahn equations with cross-kinetic coupling. Recently, in [START_REF] Boyer | Hierarchy of consistent n-component Cahn-Hilliard systems[END_REF], the authors proposed a novel hierarchy of multi-species Cahn-Hilliard systems which are consistent with the standard two-species Cahn-Hilliard system, and which read as the model introduced above with Γ positive definite, a particular Ψ and for a constant mobility matrix. Numerical methods for such systems were proposed and analyzed in several contributions, see e.g. [START_REF] Eyre | Unconditionally gradient stable time marching the Cahn-Hilliard equation[END_REF][START_REF] Wu | Multiphase Allen-Cahn and Cahn-Hilliard models and their discretizations with the effect of pairwise surface tensions[END_REF][START_REF] Chen | Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential[END_REF][START_REF] Bailo | Unconditional bound-preserving and energy-dissipating finite-volume schemes for the Cahn-Hilliard equation[END_REF].

Concerning coupled cross-diffusion Cahn-Hilliard systems, other than [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF], the only work we are aware of is [START_REF] Huo | Existence and weak-strong uniqueness for Maxwell-Stefan-Cahn-Hilliard systems[END_REF], which treats the case where all species aim to separate, i.e. the case when Γ in (1.6) is positive definite.

Structure-preserving finite volume schemes. The finite volume method is a classical discretization method to approximate conservation laws, see a pedagogical introduction in [START_REF] Eymard | Finite volume methods[END_REF]. It is a natural physical requirement for a discretization method to preserve as much as possible of the structure of the continuous problem such as conservation laws, nonnegativity or dissipation. In addition, such properties can be mathematically useful, since they enable to "transfer" the mathematical analysis to the discrete level. A method that preserves the dissipation of an energy (resp. entropy) is often called "energy-stable" (resp. entropy-stable). Following the success of the entropy method, there has been considerable effort in order to preserve the entropic structure of scalar parabolic equations [START_REF] Bessemoulin-Chatard | A finite volume scheme for nonlinear degenerate parabolic equations[END_REF][START_REF] Chainais-Hillairet | Entropy method and asymptotic behaviours of finite volume schemes[END_REF][START_REF] Cancès | Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure[END_REF] and parabolic systems [START_REF] Cancès | Finite-volume scheme for a degenerate cross-diffusion model motivated from ion transport[END_REF][START_REF] Cancès | A convergent entropy diminishing finite volume scheme for a cross-diffusion system[END_REF][START_REF] Cancès | Finite volumes for the Stefan-Maxwell cross-diffusion system[END_REF][START_REF] Carrillo | Convergence of a finite volume scheme for a system of interacting species with cross-diffusion[END_REF][START_REF] Jüngel | A finite-volume scheme for a cross-diffusion model arising from interacting many-particle population systems[END_REF][START_REF] Zurek | A convergent structure-preserving finite-volume scheme for the Shigesada-Kawasaki-Teramoto population system[END_REF][START_REF] Daus | Convergence of a finite-volume scheme for a degenerate-singular cross-diffusion system for biofilms[END_REF][START_REF] Jüngel | A discrete boundedness-by-entropy method for finite-volume approximations of cross-diffusion systems[END_REF][START_REF] Cancès | A convergent finite volume scheme for dissipation driven models with volume filling constraint[END_REF][START_REF] Herda | Study of a structure preserving finite volume scheme for a nonlocal cross-diffusion system[END_REF][START_REF] Jüngel | A convergent finite-volume scheme for nonlocal cross-diffusion systems for multi-species populations[END_REF][START_REF] Cancès | Structure preserving finite volume approximation of cross-diffusion systems coupled by a free interface[END_REF] at the discrete level. See also a review of energy-stable schemes for the Cahn-Hilliard equation in [START_REF] Brachet | Convergence to equilibrium for time and space discretizations of the Cahn-Hilliard equation[END_REF].

1.1. Contributions and outline. Our work makes the following contributions • Proving existence and uniform lower and upper bounds for the local minimizers of (1.1) in the L ∞ (Ω) n+1 topology. We emphasise that the latter, in contrast to the results of [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF], requires a construction which has to preserve not only the constraints (1.4) but also the mass of the competitor, which significantly complicates the argument. • Gaining regularity of the minimizers from the Euler-Lagrange system, we show that they qualify as classical solutions to the stationary system. We also show that the Euler-Lagrange equation for the void species decouples, revealing a strong link with the single-species energy. • We study the convexity properties of (1.1) and are able to give explicit quantitative bounds. In a particular parameters regime, we show that the minimizers are constant and that solutions to the dynamical system converge exponentially fast to them, for arbitrary initial data with finite energy. We give explicit rates of convergence. • We introduce a two-point finite volume scheme that approximates the evolution problem (1.2), preserving the constraints (1.4). The discrete free energy is shown to be nonincreasing, adapting the convex-concave splitting of [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF] to the discrete case. We provide numerical simulations to illustrate the behaviour of the scheme and to investigate the variety of stationary solutions in the long-time limit.

Remark 1.1 (Nonlocal and potential contributions to the energy). We remark that most of the results of this work remain valid, after minor modifications, if potential or non-local interaction terms of the form

Ω V i (x)u i (x) dx or c ij Ω u i L * u j dx
are added to the energy. Here, for all 0 ≤ i ≤ n, V i : Ω → R is a given potential and L : Ω → R is an interaction kernel. All these functions must be sufficiently smooth. With these additions, existence of minimizers (lemma 2.1), strict bounds (theorem 2.3) hold without any changes. First order optimality conditions have to be adapted and the regularity of solutions (theorem 2.4) is limited by the regularity of L and V := (V i ) 0≤i≤n . Under suitable assumptions on the matrix C = (c ij ) 0≤i,j≤n the numerical scheme can be adapted and still preserves the structure.

The paper is organised as follows: Section 2 contains an analysis of properties of the energy functional and establishes the link with stationary solutions. Section 3 is dedicated to the large-time asymptotics in a globally stable regime. Section 4 is devoted to the introduction of a structure preserving finite volume scheme and some numerical results are presented in Section 5.

1.2. The model. We now present the system under consideration in full detail. For ε > 0 and β > 0 we consider the energy functional given by (1.1). We define formally the chemical potentials as variational derivatives of the energy by

µ i := D ui E(u) = ln u i for all i = 1, . . . , n, (1.7)
as well as

µ 0 := D u0 E(u) = ln u 0 -ε∆u 0 + β(1 -2u 0 ), (1.8)
so that µ := (µ 0 , µ 1 , . . . , µ n ) = D u E(u). Furthermore, we introduce the auxiliary variables

w i = ln u i -ln u 0 , i = 1, . . . , n, (1.9)
as well as

w 0 := -ε∆u 0 + β(1 -2u 0 ). (1.10)
With these definitions, (1.2) can be rewritten as

∂ t u i = div   0≤j̸ =i≤n K ij u i u j ∇(µ i -µ j )   = div   0≤j̸ =i≤n K ij u i u j ∇(w i -w j )   = div   0≤j̸ =i≤n K ij (u j ∇u i -u i ∇u j ) -K i0 u i u 0 ∇w 0   , (1.11)
for i = 1, . . . , n and

∂ t u 0 = div   n j=1 K 0j u 0 u j ∇(µ 0 -µ j )   = div   n j=1 K 0j u 0 u j ∇(w 0 -w j )   = div   n j=1 K 0j (u j ∇u 0 -u 0 ∇u j + u 0 u j ∇w 0 )   .
(1.12)

The conservative form (1.2) together with the zero-flux boundary conditions suggest that the mass of each species is conserved along the evolution. Therefore, given fixed masses m 0 , . . . , m n > 0 such that

n j=0 m j = |Ω|,
we will look for solutions to (1.5) in the admissible set

A m := u := (u 0 , . . . , u n ) ∈ (L ∞ (Ω)) n+1 : u i ≥ 0, Ω u i dx = m i , i = 0, . . . , n, n j=0 u j = 1 a.e. in Ω and u 0 ∈ H 1 (Ω) .
Note that A m is non-empty, convex, and that for any u ∈ A m , it holds 0 ≤ u i ≤ 1 for all i = 0, . . . , n.

Minimizers of the energy functional

In this section we use the direct method of the calculus of variations to prove the existence of minimizers to the energy (1.1) over the set A m : min u∈Am E(u).

(2.1)

Arguing by means of competitors, we further obtain strict bounds which then allow for higher regularity by making use of the optimality conditions. In consequence, minimizers are solutions to the stationary problem (1.5).

Lemma 2.1. Let E : A m → R be defined by (1.1). Then, E has at least one minimizer.

Proof. We apply the direct method of calculus of variations. First, using the non-negativity of the function [0, 1] ∋ x → x ln x -x + 1, together with the fact that |∇v 0 | 2 ≥ 0 and that v 0 (1 -v 0 ) ≥ 0 for any v = (v 0 , . . . , v n ) ∈ A m , we obtain that E is nonnegative on A m . Moreover, it is clear that A m contains constant solutions with finite energy. Thus, there exists a minimizing sequence

v (p) p∈N ⊂ A m such that E v (p) is bounded and lim p→∞ E v (p) = inf Am E.
In particular, we have that ∇v (p) 0 L 2 (Ω) p∈N is bounded as well. Therefore, without relabelling, up to the extraction of a subsequence, there exists u 0 ∈ H 1 (Ω) such that ∇v (p) 0 ⇀ ∇u 0 weakly in L 2 (Ω), and thanks to the uniform L ∞ -bound we have v (p) 0 → u 0 strongly in L q (Ω) for every 1 ≤ q < ∞ and a.e. in Ω. Furthermore, since v (p) i is bounded in L 2 (Ω) by construction, it follows that, up to the extraction of a subsequence, for all i = 1, . . . , n, there exists u i ∈ L 2 (Ω) such that v (p) i ⇀ u i weakly in L 2 (Ω). We also easily obtain that 0 ≤ u i ≤ 1 almost everywhere on Ω. Then the convexity of the integrands together with the strong continuity of the functional (dominated convergence) imply the lower-semicontinuity

Ω u i ln u i dx ≤ lim inf p→∞ Ω v (p) i ln v (p) i dx
as well as

Ω |∇u 0 | 2 dx ≤ lim inf p→∞ Ω ∇v (p) 0 2 dx.
Furthermore, the weak convergence in L 2 (Ω) yields

Ω -v (p) i + 1 dx → Ω -u i + 1 dx for all i = 0, . . . , n,
while the strong convergence gives

Ω v (p) 0 1 -v (p) 0 dx → Ω u 0 (1 -u 0 ) dx.
This implies

E(u) ≤ lim inf p→∞ E(v (p) ) = inf Am E,
and that Ω u i dx = m i for i = 0, . . . , n. Finally, the weak convergence in L 2 (Ω) also yields that

n i=0 u i = 1
almost everywhere in Ω so that u ∈ A m . The conclusion follows. □ Remark 2.2. We point out that the uniqueness of the minimizer is neither guaranteed nor expected, due to the non-convexity of the energy functional.

A remarkable property is that the minimizers are in the interior of the set A m , i.e. they strictly satisfy the box constraints in (1.4). This is shown by constructing suitable competitors in the following theorem.

Theorem 2.3. Let E be the energy functional given by (1.1). Then, there exists a constant δ > 0 such that for every local minimizer u ∈ A m of E for the L ∞ (Ω) n+1 topology, it holds that δ ≤ u i a.e, for all i = 0, . . . , n, which, together with the volume-filling constraint in (1.4), implies the upper bound u i ≤ 1 -nδ a.e, for all i = 0, . . . , n.

Proof. Let u := (u 0 , . . . , u n ) ∈ A m be a local minimizer of E for the L ∞ (Ω) n+1 topology, i.e. there exists ϵ > 0 such that for all v := (v 0 , . . . , v n ) ∈ A m with ∥v -u∥ L ∞ (Ω) := max i=0,...,n

∥v i -u i ∥ L ∞ (Ω) ≤ ϵ, necessarily E(v) ≥ E(u) holds.
In order to prove the assertion, we proceed as follows: first we show that there exists δ 0 > 0 such that δ 0 ≤ u 0 ≤ 1 -δ 0 . Then we proceed to show that there exists δ i > 0 such that δ i ≤ u i for i = 1, . . . , n. Finally, δ := min(δ 0 , δ 1 , . . . , δ n ) is the constant that appears in the statement.

Step 1: δ 0 ≤ u 0 . Arguing by contradiction, we assume that for all min ϵ, m 0 2|Ω| ≥ δ > 0 the set

C δ = {x ∈ Ω : u 0 (x) < δ} is such that |C δ | > 0.
We further define the set

C avg := x ∈ Ω : u 0 (x) > m 0 |Ω| ,
i.e. the part of Ω on which u 0 strictly exceeds its average. Note that

C m 0 2|Ω| ∩ C avg = ∅ and C m 0 2|Ω| > 0 imply |C avg | > 0.
We then define, for all 0 ≤ λ ≤ 1, the perturbed function

u δ,λ 0 =      δ in C δ (1 -λ)u 0 + λ m0 |Ω| in C avg u 0 otherwise, which satisfies 0 ≤ u δ,λ 0 ≤ 1 and u δ,λ 0 ∈ H 1 (Ω). Moreover, for any 0 ≤ λ ≤ ϵ/2, it holds that ∥u 0 -u δ,λ 0 ∥ L ∞ (Ω) ≤ ϵ. Observe that Ω u δ,0 0 dx = C δ δ dx + Ω\C δ u 0 dx > C δ u 0 dx + Ω\C δ u 0 dx = m 0 ,
while on the other hand

Ω u δ,λ 0 dx = C δ δ dx + Ω\(C δ ∪Cavg) u 0 dx + λ Cavg m 0 |Ω| dx + (1 -λ) Cavg u 0 dx = δ|C δ | + Ω\C δ u 0 dx + λ Cavg m 0 |Ω| -u 0 dx.
In order to estimate the last integral of the previous equality we observe that

0 = Ω m 0 |Ω| -u 0 dx = C m 0 2|Ω| m 0 |Ω| -u 0 dx + Ω\(C m 0 2|Ω| ∪Cavg) m 0 |Ω| -u 0 dx + Cavg m 0 |Ω| -u 0 dx ≥ m 0 2|Ω| |C m 0 2|Ω| | + Cavg m 0 |Ω| -u 0 dx, which implies Cavg m 0 |Ω| -u 0 dx ≤ - m 0 2|Ω| |C m 0 2|Ω| |.
Then, from the previous calculations, we conclude that for any 0 < δ < min ϵ, m0

2|Ω|

and any λ ∈ [0, ϵ/2],

Ω u δ,λ 0 dx ≤ δ -λ m 0 2|Ω| |C m 0 2|Ω| | + m 0 . (2.2)
Let us now assume that δ is chosen so that δ < ϵm0 4|Ω| . Then, it holds that

Ω u δ,ϵ/2 0 dx ≤ δ -ϵ m 0 4|Ω| |C m 0 2|Ω| | + m 0 < m 0 .
Thus, for all 0 < δ < ϵm0 4|Ω| , there exists λ * δ ∈ (0, ϵ/2) such that the function

u δ 0 := u δ,λ * δ 0 satisfies Ω u δ 0 dx = m 0 . (2.3)
Furthermore, it holds from (2.2) that λ * δ necessarily satisfies λ * δ ≤ 2|Ω|δ m0 , so that lim δ→0 λ * δ = 0. While the constructed u δ 0 preserves the mass constraint, the volume-filling constraint is no longer valid. To recover them we have to modify at least one of the other species. To this end, we make the following observation: in the set C δ , it holds that

n i=1 u i = 1 -u 0 ≥ 1 -δ.
Therefore, denoting by u := max i=1,...,n u i , we necessarily have

u ≥ 1 -δ n in C δ and u ≤ n i=1 u i = 1 -u 0 ≤ 1 - m 0 |Ω| in C avg . (2.4)
Let us now define for almost all x ∈ Ω k(x) := min k = 1, . . . , n, u k (x) = max i=1,...,n u i (x) so that u(x) = uk (x) (x) almost everywhere on Ω. Let us then denote for all k = 1, . . . , n,

C k δ := {x ∈ C δ , k(x) = k} so that C δ = n k=1 C k δ and C k δ ∩ C k ′ δ = ∅ as soon as k ̸ = k ′ .
By definition, it then holds that

u k = u in C k δ . (2.5)
On the one hand, let us introduce

m δ 0,k := C k δ (u δ 0 -u 0 ) dx. Then, it holds that 0 ≤ m δ 0,k ≤ δ|C k δ | ≤ δ|Ω|.
On the other hand, from the definition of u δ 0 and from the calculations above, it holds that

Cavg (u 0 -u δ 0 ) dx = C δ (u δ 0 -u 0 ) dx = n k=1 m δ 0,k .
Therefore, using Lemma 6.1, there always exist measurable subsets C δ,k avg for k = 1, . . . , n such that

• n k=1 C δ,k avg = C avg ; • C δ,k avg ∩ C δ,k ′ avg = ∅ as soon as k ̸ = k ′ ; • C δ,k avg (u 0 -u δ 0 ) dx = m δ 0,k = -C k δ (u 0 -u δ 0 ) dx.
We can then define for all k = 1, . . . , n the perturbed function u δ k in the following way:

u δ k = u k + (u 0 -u δ 0 ) in C k δ ∪ C δ,k avg u k otherwise. (2.6)
It can then be easily checked that ∥u δ k -u k ∥ L ∞ (Ω) ≤ ϵ for δ arbitrarily small. Moreover, as a consequence of (2.4) and (2.5), for δ small enough, it holds for all k = 1, . . . , n,

0 ≤ 1 -δ n -δ ≤ u δ k = u k + (u 0 -u δ 0 ) ≤ u k ≤ 1 in C k δ .
On the other hand, using again (2.4) and (2.5), we can estimate

0 ≤ u k ≤ u δ k = u k + (u 0 -u δ 0 ) = u k + λ * δ u 0 - m 0 |Ω| ≤ 1 - m 0 |Ω| + λ * δ 1 - m 0 |Ω| in C δ,k avg . Since λ * δ -→ δ→0 0, choosing δ small enough implies that λ * δ (1 -m 0 /|Ω|) < m 0 /|Ω| and therefore u δ k ≤ 1 in C δ,k avg . Furthermore, Ω u δ k dx = Ω u k dx + C δ,k avg ∪C k δ (u 0 -u δ 0 ) dx = Ω u k dx. (2.7)
Finally, we observe that, almost everywhere in Ω,

n i=1 u δ i = n i=1 u i + (u 0 -u δ 0 ), holds so that n i=0 u δ i = n i=0 u i = 1.
From this, we conclude that u δ := (u δ 0 , . . . , u δ n ) lies in the set A m and satisfies ∥u -u δ ∥ L ∞ (Ω) n+1 ≤ ϵ for δ small enough.

We now show that for δ ≪ 1, it holds E(u δ ) < E(u) strictly, which gives us the desired contradiction. Firstly, let us note that |∇u δ 0 | ≤ |∇u 0 |, which gives

E(u δ ) -E(u) ≤ Ω n i=0 (u δ i ln u δ i -u i ln u i ) -(u δ i -u i ) + β u δ 0 (1 -u δ 0 ) -u 0 (1 -u 0 ) dx
To estimate this difference, we first observe that (2.3)-(2.7) imply

Ω (u δ i -u i ) dx = 0 for i = 0, . . . , n.
Using the convexity of x → x ln x and the concavity of x → x(1 -x) allows to further estimate

E(u δ ) -E(u) ≤ n k=1 Ω ln u δ k + 1 (u δ k -u k ) dx + Ω ln u δ 0 + 1 (u δ 0 -u 0 ) dx + β Ω (1 -2u 0 )(u δ 0 -u 0 ) dx =: ∆E 1 + ∆E 2 + ∆E 3 .
(2.8)

To estimate ∆E 1 we first observe that

∆E 1 = n k=1 Ω ln u δ k + 1 (u δ k -u k ) dx = n k=1 C k δ ∪C δ,k avg ln u δ k + 1 (u 0 -u δ 0 ) dx = n k=1 C k δ ∪C δ,k avg -ln u δ k (u δ 0 -u 0 ) dx.
(2.9)

The quantity -ln

u δ k is nonnegative in C k δ ∪ C δ,k avg while (u δ 0 -u 0 ) is nonnegative in C k δ and nonpositive in C δ,k
avg , and furthermore (2.4) gives

-ln u k = -ln(u) ≤ -ln 1 -δ n ≤ -ln 1 n 1 - m 0 2|Ω| in C k δ .
Therefore (2.9) reduces to

∆E 1 ≤ -ln 1 n 1 - m 0 2|Ω| C δ (u δ 0 -u 0 ) dx.
In order to estimate the second term on the right-hand side of (2.8) we first observe that in the set C avg it holds ln u δ 0 ≥ ln(m 0 /|Ω|), while due to the mass conservation we have

C δ (u δ 0 -u 0 ) dx = - Cavg (u δ 0 -u 0 ) dx. (2.10)
It follows that

∆E 2 = Ω ln u δ 0 (u δ 0 -u 0 ) dx ≤ ln δ C δ (u δ 0 -u 0 ) dx + ln(m 0 /|Ω)| Cavg (u δ 0 -u 0 ) dx = (ln δ -ln(m 0 /|Ω|)) C δ (u δ 0 -u 0 ) dx.
Finally, for the last term in (2.8) we use the fact that -1 ≤ 1 -2u 0 ≤ 1 and (2.10) again to have

∆E 3 ≤ β Ω (1 -2u 0 )(u δ 0 -u 0 ) dx ≤ β Ω |u δ 0 -u 0 | dx ≤ β C δ (u δ 0 -u 0 ) dx -β Cavg (u δ 0 -u 0 ) dx = 2β C δ (u δ 0 -u 0 ) dx. (2.11)
Summarising we have

E(u δ ) -E(u) ≤ ln δ -ln(m 0 /|Ω|) -ln 1 n 1 - m 0 2|Ω| + 2β C δ (u δ 0 -u 0 ) dx.

Thus, since

C δ (u δ 0 -u 0 ) dx > 0 for all δ > 0, taking δ sufficiently small so that the constant in front of this integral becomes negative yields the desired contradiction.

Step 2: u 0 ≤ 1 -δ 0 . We argue again by contradiction and assume that for all min ϵ, m 0 2|Ω| ≥ δ > 0 the set

C δ = {x ∈ Ω : u 0 (x) > 1 -δ} is such that | C δ | > 0.
We further define the set

C avg := x ∈ Ω : u 0 (x) < m 0 |Ω| ,
i.e. the part of Ω on which u 0 is strictly below its average. As before we argue that | C δ | > 0 for all δ > 0 arbitrarily small implies | C avg | > 0. We then define, for 0 ≤ λ ≤ 1, the perturbed function

ũδ,λ 0 =      1 -δ in C δ λu 0 + (1 -λ) m0 |Ω| in C avg u 0 otherwise,
and arguing as in the previous step we show that for all ϵm0 4|Ω| > δ > 0 there exists λ * δ ∈ (0, ϵ/2) such that the function ũδ 0 := ũδ,λ * δ 0 has the same mass as u 0 and satisfies ∥u 0 -ũδ 0 ∥ L ∞ (Ω) ≤ ϵ. We then observe that it holds

n i=1 u i = 1 -u 0 ≥ 1 - m 0 |Ω| in C avg ,
which implies that the function u := max i=1,...,n

u i must satisfy u ≥ 1 n 1 - m 0 |Ω| in C avg and u ≤ n i=1 u i = 1 -u 0 ≤ δ in C δ .
For almost all x ∈ Ω, we denote by k(x) := min k = 1, . . . , n, u k (x) = max i=1,...,n u i (x) . Moreover, for all k = 1, . . . , n we denote by

C k δ := x ∈ C δ , k(x) = k so that C δ = n k=1 C k δ and C k δ ∩ C k ′ δ = ∅ as soon as k ̸ = k ′ .
By definition, it then holds that

u k = u in C k δ .

By the mass conservation property

Cavg (u 0 -ũδ 0 ) dx = C δ (ũ δ 0 -u 0 ) dx, (2.12)
for all δ > 0 sufficiently small, using again Lemma 6.1, there exist measurable subsets C δ,k avg ⊂ C avg for all k = 1, . . . , n such that

• C avg = n k=1 C δ,k avg ; • C δ,k avg ∩ C δ,k ′ avg = ∅ as soon as k ̸ = k ′ ; • C δ,k avg (u 0 -ũδ 0 ) dx = C k δ (ũ δ 0 -u 0 ) dx.
Thus, for δ sufficiently small, for all k = 1, . . . , n, we define the function

ũδ k = u k + (u 0 -u δ 0 ) in C k δ ∪ C δ,k avg u k otherwise which is such that 0 ≤ ũδ k ≤ 1 and Ω ũδ k dx = Ω u k dx.
Using similar arguments as in Step 1, we obtain that for δ sufficiently small, ũδ := (ũ δ 0 , . . . , ũδ n ) ∈ A m and ∥ũ δ -u∥ L ∞ (Ω) n+1 ≤ ϵ. Moreover, it holds that

E(ũ δ ) -E(u) ≤ n k=1 Ω ln ũδ k (ũ δ k -u k ) dx + Ω ln ũδ 0 (ũ δ 0 -u 0 ) dx + β Ω (1 -2u 0 )(ũ δ 0 -u 0 ) dx =: ∆E 1 + ∆E 2 + ∆E 3 .
To estimate ∆E 1 we first observe that

-ln u k ≤ -ln 1 n (1 -m0 |Ω| ) in C δ,k avg ≥ -ln δ in C k δ . Thus we can estimate ∆E 1 = n k=1 Ω ln ũδ k (ũ δ k -uk) dx = n k=1 C δ,k avg ∪ C k δ ln ũδ k (ũ δ k -uk) dx = - n k=1 C k δ ln ũδ k (ũ δ 0 -u 0 ) dx - C δ,k avg ln ũδ k + 1 (ũ δ 0 -u 0 ) dx ≤ -ln δ + ln 1 n (1 - m 0 |Ω| ) C δ (ũ δ 0 -u 0 ) dx,
where we used (2.12) to obtain the last inequality.

The second integral ∆E 2 in (2.8) can be estimated via

∆E 2 = Ω ln ũδ 0 (ũ δ 0 -u 0 ) dx ≤ ln(1 -δ) C δ (ũ δ 0 -u 0 ) dx + ln(m 0 /|Ω)| Cavg (ũ δ 0 -u 0 ) dx,
where we used ln(ũ δ 0 ) ≤ ln(m 0 /|Ω|) in C avg . Using again (2.12), we obtain

∆E 2 ≤ (ln(1 -δ) -ln(m 0 /|Ω|)) C δ (ũ δ 0 -u 0 ) dx.
Estimating the concave part ∆E 3 as in (2.11) and collecting all terms eventually yields

E(ũ δ ) -E(u) ≤ -ln δ + ln 1 n 1 - m 0 |Ω| + ln(1 -δ) -ln(m 0 /|Ω|) + 2β C δ (ũ δ 0 -u 0 ) dx.
As the integral on the right-hand side is strictly negative but the coefficient in front of it becomes positive for δ sufficiently small, we again reach the desired contradiction.

Step 3: δ ≤ u i , i = 1, . . . n. To fix the ideas let us assume that i = 1. While the proof uses the same construction as in Step 1, it is crucial to make sure that the index k(x) used to construct the sets C k δ and the functions such as in (2.6) is such that k(x) is never equal to 0, as applying (2.6) to define u δ 0 with k = 0 would yield a function u δ 0 which does not belong to H 1 (Ω) and thus renders the value of the energy to be infinity. To this end, we may use the upper bound on u 0 established at Step 2 of the proof to calculate j̸ =0,1

u j = 1 -u 0 -u 1 ≥ (δ 0 -δ) in the set C δ,1 := {x ∈ Ω : u 1 (x) < δ}.
As δ 0 from Steps 1 and 2 is fixed at this point, choosing δ < δ 0 we can go on from here to ensure that, defining k(x) := min{k = 2, . . . , n, u k (x) = max i=2,...,n u i (x)}, we have that uk (x) (x) > (δ 0 -δ)/n almost everywhere in C δ,1 . Then, arguing as in Step 1 gives the existence of δ 1 > 0 such that δ 1 ≤ u 1 . As the choice i = 1 was arbitrary it follows that the argument can then be applied to any other u j , using the same construction as just done for u 1 .

Step 4: Conclusion. We then observe that the parameter δ := min(δ 0 , δ 1 , . . . , δ n ) satisfies all the properties in the statement, and therefore the conclusion follows. □

Thanks to these uniform bounds and arguing by elliptic regularity, we can derive first order optimality conditions as given by the following theorem.

Theorem 2.4. Let u ∈ A m be a local minimizer of the energy (1.1) in the L ∞ (Ω) n+1 topology. Then u 0 ∈ H 2 (Ω) and is solution to

-ε∆u 0 = ln 1 -u 0 u 0 -β(1 -2u 0 ) - 1 |Ω| Ω ln 1 -u 0 u 0 -β(1 -2u 0 ) dx, in Ω, ∂u 0 ∂n = 0, on ∂Ω.
(2.13) Moreover, it holds

u i = m i |Ω| -m 0 (1 -u 0 ), i = 1, . . . , n. (2.14)
In addition, u is an element of (C ∞ (Ω)) n+1 and a classical solution to (1.5).

Proof.

Step 1: establishing (2.13) and (2.14). Given ψ := (ψ 1 , . . . , ψ n ) ∈ (C ∞ (Ω)) n , set

ν i := ψ i - 1 |Ω| Ω ψ i dx for all i = 1, . . . , n. (2.15)
Fix now r > 0, consider the perturbations u r,i := u i + rν i for all i = 1, . . . , n, and set

u r,0 := 1 - n i=1 u r,i = 1 - n i=1 u i -r n i=1 ν i = u 0 + rν 0 , with ν 0 := - n i=1 ν i .
We then set ν := (ν 0 , . . . , ν n ), from which u r = u + rν and in turn u r ∈ A m for r sufficiently small, due to the strict lower and upper bounds of minimizers shown in Theorem 2.3, so that E(u) ≤ E(u r ). We now want to calculate the variation of E, that is,

lim r→0 E(u r ) -E(u) r .
Note that, thanks again to the strict bounds of Theorem 2.3, we obtain uniform L ∞ -bounds on ln u i , ln u r,i , for any i = 0, . . . , n and for r sufficiently small. Therefore, we can use dominated convergence to calculate

0 ≤ lim r→0 E(u r ) -E(u) r = Ω n i=1 ln u i -ln u 0 -β(1 -2u 0 ) ν i -ε∇u 0 • ∇ν i dx,
for all such ν i . We repeat the same argument for -ν and then use (2.15) to eventually infer, using Fubini theorem,

0 = Ω n i=1 ln u i -ln u 0 -β(1 -2u 0 ) ψ i -ε∇u 0 • ∇ψ i dx - Ω 1 |Ω| Ω ln u i -ln u 0 -β(1 -2u 0 ) dy ψ i dx for every ψ = (ψ 1 , . . . , ψ n ) ∈ (C ∞ (Ω)) n .
Setting

λ i = 1 |Ω| Ω ln u i -ln u 0 -β(1 -2u 0 ) dy,
we obtain that for any i = 1, . . . , n, for any φ ∈ H 1 (Ω),

Ω ln u i -ln u 0 -β(1 -2u 0 ) φ dx -ε Ω ∇u 0 • ∇φ dx = λ i Ω φ dx. (2.16)
Thanks to the uniform bounds of Theorem 2.3, we know that ln u i -λ i -ln u 0 -β(1 -2u 0 ) ∈ L 2 (Ω) for all i = 1, . . . , n, therefore, using standard elliptic regularity theory (see, e.g., [START_REF] Troianiello | Elliptic differential equations and obstacle problems[END_REF]), we get from (2.16) that u 0 ∈ H 2 (Ω). This in turn implies that the Euler-Lagrange system is satisfied in strong form as follows

-ε∆u 0 = ln u i -λ i -ln u 0 -β(1 -2u 0 ) a.e.
in Ω, for all i = 1, . . . , n.

(2.17)

Let us now proceed to eliminate the Lagrange multipliers from the equations. We write (2.17) for any pair of indexes i ̸ = k ∈ {1, . . . , n} and take the difference of the corresponding expressions. This gives

ln u i -ln u k = λ i -λ k , (2.18)
that is, u k = u i e λ k -λi . Solving this equation for u i by using the volume-filling constraint gives

u i = e λi n k=1 e λ k (1 -u 0 ). (2.19)
Inserting (2.19) in (2.17), the Euler-Lagrange system reduces to the following PDE on u 0 together with n Lagrange multipliers

-ε∆u 0 = ln 1 -u 0 u 0 -β(1 -2u 0 ) -ln n k=1 e λ k ,
e λ k = 1 |Ω| Ω ln 1 -u 0 u 0 -β(1 -2u 0 ) dx =: λ 0 .
Then integrating (2.19) over Ω we obtain that

e λi = m i |Ω| -m 0 e λ0 , i = 1, . . . , n,
and in turn (2.14).

Step 2: establishing (1.5). First, it follows from (2.13) combined with the uniforms bounds on u 0 that, by elliptic regularity, u 0 is smooth in Ω, and as a consequence of (2.14), u is smooth as well. Secondly, it is enough to show that u i satisfies (1.11) for i = 1, . . . , n, since then (1.12) is automatically satisfied thanks to the volume-filling constraint. Remark that, rewriting (2.17) and (2.18) with the notations (1.9)-(1.10) leads to

w i -w 0 = λ i , i = 1, . . . , n, w i -w j = λ i -λ j , i ̸ = j = 1, . . . , n.
Taking into account that u is independent of time and using the previous relations, we obtain from (1.11) that div

  0≤j̸ =i≤n K ij u i u j ∇(w i -w j )   = div   1≤j̸ =i≤n K ij u i u j ∇(w i -w j ) + K i0 u i u 0 ∇(w i -w 0 )   = div   1≤j̸ =i≤n K ij u i u j ∇(λ i -λ j ) + K i0 u i u 0 ∇λ i   = 0,
which concludes the proof. □

Let us make a few comments on the result that we have established here. First, we have proved that any local minimizer of the energy functional E for the L ∞ (Ω) n+1 topology is indeed a solution to the stationary problem (1.5). Of course, we do not expect the converse to be true, in the sense that there may exist other solutions to (1.5) which are not local minimizers of the energy functional E. However, because of the gradient flow structure of the time-dependent system (1.2), it is natural to conjecture that solutions to this system will converge in the long-time limit to some stationary states that are local minimizers of E. We are not able to prove this claim here, but give some numerical evidence in Section 5.

Second, thanks to Theorem 2.4, we can study the properties of the local minimizers of E by studying the scalar equation (2.13), which is now the first-order Euler-Lagrange equation of the single-species Cahn-Hilliard energy

E 0 (u 0 ) = Ω u 0 ln u 0 + (1 -u 0 ) ln(1 -u 0 ) + ε 2 |∇u 0 | 2 + βu 0 (1 -u 0 ) dx.

Convexity properties and long-time behaviour

In this section we present some results obtained on the behaviour of solutions to the time-dependent system (1.2) in the case when the parameters ε and β are chosen such that the functional E is convex. While this range of parameter values may not be practically relevant for some physical applications where separation effects dominate over diffusion, we nevertheless believe that the present analysis is instructive and might be seen as a useful preliminary step towards the study of the long-time behaviour of solutions to (1.2) in the general case. We first give explicit conditions on the parameters ε and β for E to be convex. We then prove that, in a stable regime, solutions of the time-dependent system converge exponentially fast to the minimizer of E, which is proved to be unique in this setting.

Convexity properties.

Let us now study convexity properties of the energy (1.1). We begin by recalling the famous Poincaré-Wirtinger inequality: there exists a constant C p > 0 such that for all v ∈ H 1 (Ω),

Ω v - 1 |Ω| Ω v dy 2 dx ≤ C p Ω |∇v| 2 dx, (3.1)
as well as the logarithmic-Sobolev inequality on bounded domains, [22, Lemma 1],

f log f dx ≤ C sob Ω |∇ f | 2 dx, for every f ∈ H 1 (Ω; R + ) s.t. Ω f dx = 1, (3.2)
where the constant C sob depends on d and Ω, only. We then have the following lemma:

Lemma 3.1 (Convexity of the free energy). Let n ≥ 1 and C p be the Poincaré constant in (3.1). Assume that

1 2m 0 + 1 |Ω| ε 2C p -β ≥ 0. (3.3)
Then, the energy functional (1.1) is convex in the set A m . If n = 1, the result holds under the weaker condition

|Ω| 2m 0 (|Ω| -m 0 ) + 1 |Ω| ε 2C p -β ≥ 0. (3.4)
Moreover, whenever the inequalities are strict, the energy is strictly convex.

Proof. Let u := (u 0 , . . . , u n ), v := (v 0 , . . . , v n ) ∈ A m . This in particular implies that

Ω u i dx = Ω v i dx for all i = 0, . . . , n.
We want to estimate from below the quantity

E(u) -E(v) -E ′ (v)(u -v) = n i=0 Ω u i ln u i v i + ε 2 |∇(u 0 -v 0 )| 2 -β(u 0 -v 0 ) 2 dx.
The first terms can be estimated thanks to the Csiszár-Kullback-Pinsker inequality as

Ω u i ln u i v i dx ≥ 1 2m i Ω |u i -v i | dx 2 , i = 0, . . . , n,
while the gradient term is estimated using the Poincaré inequality (3.1). We obtain

E(u) -E(v) -E ′ (v)(u -v) ≥ n i=0 1 2m i Ω |u i -v i |dx 2 + ε 2C p -β Ω |u 0 -v 0 | 2 dx ≥ n i=0 1 2m i Ω |u i -v i |dx 2 + 1 |Ω| ε 2C p -β Ω |u 0 -v 0 |dx 2 , (3.5)
where we applied Jensen's inequality. If n = 1, using

u 1 = 1 -u 0 , m 1 = |Ω| -m 0 , we obtain E(u) -E(v) -E ′ (v)(u -v) ≥ |Ω| 2m 0 (|Ω| -m 0 ) + 1 |Ω| ε 2C p -β Ω |u 0 -v 0 | dx 2 ,
hence condition (3.4) ensures convexity. However, when n ≥ 2, one cannot easily take advantage of the terms

n i=1 1 2m i Ω |u i -v i |dx 2
. Therefore, using only the fact that they are nonnegative, one obtains from (3.5) that

E(u) -E(v) -E ′ (v)(u -v) ≥ 1 2m 0 + 1 |Ω| ε 2C p -β Ω |u 0 -v 0 | dx 2 ,
which gives again the convexity of E under condition (3.3). □

When applicable, it follows from strict convexity that both the minimization problem and the Euler-Lagrange equation (2.13) have a unique solution. Moreover, the constant state u 0 = m0 |Ω| always solves (2.13). This leads to the following corollary. Corollary 3.2 (Uniqueness of the minimizer). Whenever condition

1 2m 0 + 1 |Ω| ε 2C p -β > 0
holds, the minimization problem (2.1) has a unique solution, given by the constant states

u ∞ i = m i |Ω| for all i = 0, . . . , n. (3.6)
3.2. Large-time asymptotics in the stable regime. We discuss now the large-time asymptotics of the evolution system (1.2) in the framework of Corollary 3.2, i.e. when the energy admits a unique (constant) minimizer. We aim to show that u converges exponentially fast to this unique minimizer, for arbitrary large initial data with finite relative energy.

It follows from [23, Theorem 2.1] that, for any initial condition in A m , the system admits a weak solution u = (u 0 , u 1 , . . . , u n ) that satisfies the constraints (1.4) and such that (see [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF]Definition 2.1] for more details),

• u i ∈ L 2 ((0, T ); H 1 (Ω)) for all i = 1, . . . , n;

• u 0 ∈ L 2 ((0, T ); H 2 (Ω)); • ∂ t u i ∈ L 2 ((0, T ); (H 1 (Ω)) ′ ) for all i = 0, . . . , n; Let u ∞ := (u ∞ 0 , . . . , u ∞ n
) be defined by (3.6). The relative energy functional is defined via the Bregman divergence as, for all u := (u 0 , . .

. , u n ) ∈ A m , RE[u | u ∞ ] = E(u) -E(u ∞ ) -E ′ (u ∞ )(u -u ∞ ) = Ω n i=0 u i ln u i u ∞ i + ε 2 |∇(u 0 -u ∞ 0 )| 2 -β(u 0 -u ∞ 0 ) 2 dx. (3.7)
Our result is stated in the next theorem.

Theorem 3.3 (Global exponential stability of the constant steady state). Let u ∞ be given by (3.6) and let u 0 ∈ A m be an initial condition such that u 0 ∈ H 2 (Ω) n+1 and

RE[u 0 | u ∞ ] < ∞.
Let u be a weak solution to (1.2) as constructed in [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF]. Then, u satisfies

RE[u(t) | u ∞ ] ≤ e -λt RE[u 0 | u ∞ ], t > 0, (3.8)
with the rate

λ = 4k min 1 C sob , 1 C p - 2β ε (3.9)
with k := min 0≤i̸ =j≤n K ij > 0, C p and C sob being the constants in (3.1) and (3.2), respectively. Furthermore, under the condition

ε 2C p -β > 0, (3.10) it holds ∥u i (t) -u ∞ i ∥ 2 L 1 (Ω) ≤ RE[u(t) | u ∞ ] ≤ e -λt RE[u 0 | u ∞ ] → 0 as t → +∞
, for all i = 0, 1, . . . , n. Remark 3.4. The global stability condition (3.10) is stronger than the convexity condition on the energy (3.3) as it does not take the mass into account.

We recall here some results proved in [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF] and in particular the construction of a weak solution to the timedependent system. The method consists in introducing a time-discrete regularised version of the evolution system, together with a suitable weak formulation. More precisely, the following theorem holds: Theorem 3.5 (Theorem 3.1 from [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF]). Define the sets

A := u := (u i ) 1≤i≤n ∈ (L ∞ (Ω)) n : u i ≥ 0, i = 1, . . . , n, u 0 := 1 - n i=1 u i ≥ 0 and B := ϕ := (ϕ i ) 1≤i≤n ∈ (L ∞ (Ω)) n : ϕ 0 := - n i=1 ϕ i ∈ H 1 (Ω) ,
and let τ > 0 be a discrete time step, let p ∈ N, and let u p ∈ A ∩ (H 2 (Ω)) n . Then, there exists a solution (u p+1 , wp+1 ) ∈ (A ∩ (H 2 (Ω)) n ) × (H 2 (Ω)) n to the following coupled system: for all 1 ≤ i ≤ n, for all

ϕ i ∈ H 2 (Ω), Ω u p+1 i -u p i τ ϕ i dx = - Ω 1≤j̸ =i≤n K ij u p+1 i u p+1 j ∇( wp+1 i -wp+1 j ) + K i0 u p+1 i u p+1 0 ∇ wp+1 i • ∇ϕ i dx -τ ⟨ wp+1 i , ϕ i ⟩ H 2 (Ω) , (3.11
)

and for all ψ = (ψ i ) 1≤i≤n ∈ B ∩ (L ∞ (Ω)) n , n i=1 Ω (ln u p+1 i -ln u p+1 0 )ψ i + ε∇u p+1 0 • ∇ψ 0 dx = n i=1 Ω wp+1 i + β(1 -2u p 0 ) ψ i dx, (3.12)
where u p 0 = 1 -n i=1 u p i . Moreover, the function u p+1 satisfies the following property: there exists δ p > 0 such that

u p+1 i ≥ δ p ,
for all 1 ≤ i ≤ n, and u p+1 0

:= 1 - n i=1 u p+1 i ≥ δ p , a.e. in (0, T ) × Ω.
Proof of Theorem 3.3. Let u 0 := (u 0 0 , . . . , u 0 n ) ∈ A m ∩ H 2 (Ω) n+1 . For a given value of time step τ > 0, and starting from u 0 := (u 0 1 , . . . , u 0 n ), we consider the sequence (u p ) p∈N of discrete iterates given by Theorem 3.5. For all p ∈ N, we then define u p = (u p 0 , u p ) where u p 0 := 1 -n i=1 u p i . We first note that they enjoy enough regularity to define the following quantities for all p ≥ 0, using w p+1 := (w p+1 1 , . . . , w p+1 n ),

w p+1/2 0 := -ε∆u p+1 0 + β(1 -2u p 0 ), w p+1 i := w p+1 i + w p+1/2 0 , for i = 1, . . . , n, w p+1 := (w p+1 i ) 1≤i≤n .
Then, one can define the piecewise constant interpolation u (τ ) as well as its time-shifted version σ τ u (τ ) : for all p ∈ N \ {0}, define the discrete time t p = pτ and, for all t ∈ (t p , t p+1 ],

u (τ ) (t) = u p+1 , σ τ u (τ ) (t) = u p-1 , w (τ ) (t) = w p+1 ,
together with u (τ ) (0) = u 0 . For any t ∈ (0, ∞), we also define P (τ ) (t) ∈ N \ {0} to be the lowest integer such that t P (τ ) ≥ t.

Choosing ϕ i = wp+1 i = w p+1 i -w p+1/2 0 in (3.11
) one obtains a discrete (relative) energy-energy dissipation inequality and the authors showed in the proofs of [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF]Lemmas 4.2,4.3,4.4] that the terms in the dissipation can be estimated, which yields the following inequality, for all t > 0,

RE[u (τ ) (t) | u ∞ ] -RE[u 0 | u ∞ ] ≤ -k n i=0 t 0 Ω |∇u (τ ) i | 2 u (τ ) i dxds + P (τ ) (t)-1 p=0 4kβτ ∥∇u p+1 0 ∥ L 2 (Ω) ∥∇u p 0 ∥ L 2 (Ω) -2kε t 0 Ω |∆u (τ ) 0 | 2 dxds -k t 0 Ω u (τ ) 0 (1 -u (τ ) 0 )|∇w (τ ) 0 | 2 dxds -τ t 0 n i=1 ∥w (τ ) i -w (τ ) 0 ∥ 2 H 2 (Ω) ds.
Considering the nonnegativity of the two last terms, and using the inequality

P (τ ) (t)-1 p=0 4kβτ ∥∇u p+1 0 ∥ L 2 (Ω) ∥∇u p 0 ∥ L 2 (Ω) ≤ 2kβ t 0 ∥∇u (τ ) 0 ∥ 2 L 2 (Ω) + ∥∇σ τ u (τ ) 0 ∥ 2 L 2 (Ω) ds,
the energy inequality rewrites

RE[u (τ ) (t) | u ∞ ] -RE[u 0 | u ∞ ] ≤ -4k n i=0 t 0 Ω |∇ u (τ ) i | 2 dxds -2kε t 0 Ω |∆u (τ ) 0 | 2 dxds + 2kβ t 0 ∥∇u (τ ) 0 ∥ 2 L 2 (Ω) + ∥σ τ ∇u (τ ) 0 ∥ 2 L 2 (Ω) ds.
We use the logarithmic Sobolev inequality (3.2), together with the following estimate, which we prove for

any v ∈ H 2 (Ω) with ∂ n v = 0 on ∂Ω: denoting by m = 1 |Ω| Ω v dx, it holds that Ω |∇v| 2 dx = Ω |∇(v -m)| 2 dx = - Ω (u -m)∆(u -m) ≤ Ω (u -m) 2 1/2 Ω (∆(u -m)) 2 1/2 ≤ C p Ω |∇(u -m)| 2 1/2 Ω (∆(u -m)) 2 1/2 ≤ 1 2 Ω |∇u| 2 dx + C p 2 Ω (∆u) 2 i.e. ∀v ∈ H 2 (Ω) with ∂ n v = 0 on ∂Ω, Ω |∇v| 2 dx ≤ C p Ω (∆v) 2 dx.
We thus obtain, using the fact that for all p ∈ N,

u p+1 0 ∈ H 2 (Ω) is such that ∂ n u p+1 0 = 0 (
this is a consequence of the weak variational fomulation (3.12)):

RE[u (τ ) (t) | u ∞ ] -RE[u 0 | u ∞ ] ≤ - 4k C sob t 0 Ω n i=0 u (τ ) i ln u (τ ) i u
∞ i dxds - 2kε C P t 0 ∥∇u (τ ) 0 ∥ 2 L 2 (Ω) ds + 2kβ t 0 ∥∇u (τ ) 0 ∥ 2 L 2 (Ω) + ∥σ τ ∇u (τ ) 0 ∥ 2 L 2 (Ω) ds.
We need to pass to the limit as τ → 0 + before estimating further. On the one hand, using the analysis of [START_REF] Ehrlacher | Existence of weak solutions to a Cross-Diffusion Cahn-Hilliard type system[END_REF], it holds that (u (τ ) 0 ) τ >0 and (σ τ u (τ ) 0 ) τ >0 converge strongly to u 0 in L 2 ((0, t); H 1 (Ω)), so we can pass to the limit in the three gradient terms. On the other hand, (u (τ ) i ) τ >0 converges strongly to u i in L 2 ((0, t); L 2 (Ω)) for any i = 1, . . . , n, so by the continuity of the function [0, 1] ∋ x → x ln x and Lebesgue dominated convergence theorem, we can pass to the limit in the mixing terms, and therefore in the entire energy itself as well. We obtain

RE[u(t) | u ∞ ] -RE[u 0 | u ∞ ] ≤ - 4k C sob t 0 Ω n i=0 u i ln u i u ∞ i dxds -4k( ε 2C p -β) t 0 ∥∇u 0 ∥ 2 L 2 (Ω) ds ≤ -4k min 1 C sob , 1 C p - 2β ε t 0 RE[u(t) | u ∞ ] ds,
so that applying an integral version of Gronwall inequality gives (3.8)-(3.9). Moreover, under the same condition, the Poincaré and the Csiszár-Kullback inequalities (see [START_REF] Jüngel | Entropy methods for diffusive partial differential equations[END_REF]) show that the relative energy (3.7) dominates the square of the L 1 -norm, hence the conclusion. □ Remark 3.6. Note that the scope of the previous result is restricted to a regime where the system evolution is mainly diffusion-driven. In particular, the system does not enjoy phase separation. The dynamics outside of this regime is more relevant, but also much more difficult to study. Some constant states may become unstable, leading to the concepts of spinodal region and spinodal decomposition (see Figure 5). This phenomenon has been studied in the context of the single-species Cahn-Hilliard equation [START_REF] Novick-Cohen | The Cahn-Hilliard equation[END_REF] and for Cahn-Hilliard systems [START_REF] Eyre | Systems of Cahn-Hilliard equations[END_REF]. It is an interesting perspective to study it for cross-diffusion-Cahn-Hilliard systems, which we postpone to future work.

Finite volume scheme

We introduce a two-point finite volume scheme to solve (1.2) and study its large-time asymptotics. In the next subsection, we introduce a suitable discretization of the space-time domain and useful notations. Then we present our scheme and prove several properties related to the preservation of the structure of the continuous system. 4.1. Mesh and notations. An admissible mesh of Ω is a triplet (T , E, (x K ) K∈T ) such that the following conditions are fulfilled.

(i) Each control volume (or cell) K ∈ T is non-empty, open, polyhedral and convex. We assume that

K ∩ L = ∅ if K, L ∈ T with K ̸ = L, while K∈T K = Ω. (ii) Each face σ ∈ E is closed and is contained in a hyperplane of R d , with positive (d -1)-dimensional
Hausdorff (or Lebesgue) measure denoted by

m σ = H d-1 (σ) > 0. We assume that H d-1 (σ ∩ σ ′ ) = 0 for σ, σ ′ ∈ E unless σ = σ ′ .
For all K ∈ T , we assume that there exists a subset E K of E such that ∂K = σ∈E K σ. Moreover, we suppose that K∈T E K = E. Given two distinct control volumes K, L ∈ T , the intersection K ∩ L either reduces to a single face σ ∈ E denoted by K|L, or its (d -1)-dimensional Hausdorff measure is 0. (iii) The cell-centers (x K ) K∈T satisfy x K ∈ K, and are such that, if K, L ∈ T share a face K|L, then the vector x L -x K is orthogonal to K|L.

We denote by m K the d-dimensional Lebesgue measure of the control volume K. The set of the faces is partitioned into two subsets: the set E int of the interior faces defined by

E int = {σ ∈ E | σ = K|L for some K, L ∈ T },
and the set E ext = E \E int of the exterior faces defined by E ext = {σ ∈ E | σ ⊂ ∂Ω}. For a given control volume K ∈ T , we also define

E K,int = E K ∩ E int (respectively E K,ext = E K ∩ E ext
) the set of its faces that belong to E int (respectively E ext ). For such a face σ ∈ E K,int , we may write σ = K|L, meaning that σ = K ∩ L, where L ∈ T . Given σ ∈ E, we let

d σ := |x K -x L | if σ = K|L ∈ E int , |x K -x σ | if σ ∈ E K,ext , and
τ σ = m σ d σ .
In what follows, we use boldface notations for any vector-valued quantity, typically for elements of R n+1 ,

R |T | , R |E| , (R |T | ) n+1 and (R |E| ) n+1
. Moreover, we use uppercase letters to denote discrete quantities, in contrast to lowercase letters used in the previous sections for functions. Given any discrete scalar field V = (V K ) K∈T ∈ R |T | , we define for all cell K ∈ T and interface σ ∈ E K the mirror value V Kσ of V K across σ by setting:

V Kσ = V L if σ = K|L ∈ E int , V K if σ ∈ E ext .
We also define the oriented and absolute jumps of V across any edge by

D Kσ V = V Kσ -V K , and D σ V = |D Kσ V |, ∀K ∈ T , ∀σ ∈ E K .
Note that in the above definition, for all σ ∈ E, the definition of D σ V does not depend on the choice of the element K ∈ T such that σ ∈ E K . Therefore, it can be checked that the following discrete integration by parts formula holds, for any V , W ∈ R |T | :

K∈T σ∈E K D Kσ W .V K = - σ∈Eint D Kσ W .D Kσ V . (4.1)
Concerning the time discretization of (0, T ), we consider P T ∈ N * and an increasing infinite family of times 0 < t 0 < t 1 < • • • < t P T = T and we set ∆t p = t p -t p-1 for p ∈ {1, • • • , P T }.

4.2. Numerical Scheme. Let u 0 = (u 0 0 , . . . , u 0 n ) ∈ A m be an initial condition satisfying the constraints (1.4). It is discretized on the mesh T as

U 0 = U 0 i,K K∈T ,0≤i≤n
by setting

U 0 i,K = 1 m K K u 0 i (x) dx, ∀K ∈ T ; i = 0, . . . , n. (4.2)
Assume that U p = U p i,K K∈T ,0≤i≤n

is given for some p ∈ N, then we have to define how to compute the discrete volume fractions U p+1 = U p+1 i,K K∈T ,0≤i≤n

. For q = p, p + 1 and i = 0, . . . , n, we introduce the notation U q i := U q i,K K∈T

. We also introduce discrete fluxes

J p+1 E = J p+1 i,Kσ K∈T ,σ∈E K ,0≤i≤n
, which are based on edge values U p+1 i,σ for all σ ∈ E, i = 0, . . . , n. For any K ∈ T such that σ ∈ E K , the definition of U p+1 i,σ makes use of the values U p+1 i,K and U p+1 i,Kσ but is independent of the choice of K. The edge volume fractions are then defined through a logarithmic mean as follows

U p+1 i,σ =        0 if min(U p+1 i,K , U p+1 i,Kσ ) ≤ 0, U p+1 i,K if 0 < U p+1 i,K = U p+1 i,Kσ , U p+1 i,K -U p+1 i,Kσ ln(U p+1 i,K )-ln(U p+1 i,Kσ )
otherwise.

(4.3a)

This choice is motivated by a discrete chain rule property: for any i = 0, . . . , n, K ∈ T , σ ∈ E K , if

U p+1 i,K , U p+1 i,Kσ > 0 then D Kσ U p+1 i = U p+1 i,σ D Kσ ln(U p+1 i ). (4.3b)
We employ a time discretization relying on the backward Euler scheme:

m K U p+1 i,K -U p i,K ∆t p + σ∈E K,int J p+1 i,Kσ = 0, ∀K ∈ T , i = 0, . . . , n, (4.3c)
and the discrete fluxes are adapted from formulas (1.11)-(1.12) as

J p+1 i,Kσ = -τ σ 0≤j̸ =i≤n K ij U p+1 j,σ D Kσ U p+1 i -U p+1 i,σ D Kσ U p+1 j + τ σ K i0 U p+1 i,σ U p+1 0,σ D Kσ W p+ 1 2 0 , i = 1, . . . , n, J p+1 0,Kσ = -τ σ n i=1 K i0 U p+1 i,σ D Kσ U p+1 0 -U p+1 0,σ D Kσ U p+1 i -τ σ n i=1 K i0 U p+1 i,σ U p+1 0,σ D Kσ W p+ 1 2 0 , (4.3d)
where the auxiliary variable w 0 is discretized from (1.10) as

W p+1/2 0 = W p+1/2 0,K K∈T
, where for any K ∈ T ,

W p+ 1 2 0,K = - ε m K σ∈E K,int τ σ D Kσ U p+1 0 + β(1 -2U p 0,K ). (4.3e)
Note that, in the latter formula, we apply the same convex-concave splitting as in [START_REF] Eyre | Systems of Cahn-Hilliard equations[END_REF]: the convex part of the energy is discretized implicitly while the concave part is discretized explicitly. This well-known technique is crucial in order to recover free energy dissipation at the discrete level [START_REF] Eyre | Unconditionally gradient stable time marching the Cahn-Hilliard equation[END_REF], see the proof of Proposition 4.5.

Remark 4.1. Rather than defining the system (4.3c) for all i = 0, . . . , n and then check that it holds n i=0 U p+1 i,K = 1 for any K ∈ T , p ∈ N , one could have written only n equations and define U p+1 0,K = 1 -n i=1 U p+1 i,K (it is the chosen approach to define weak continuous solutions). However, it is important not to define the value on the edges U p+1 0,σ as 1 -n i=1 U p+1 i,σ but rather as the logarithmic mean according to (4.3a). In doing so, the saturation constraint does not necessarily hold on the edges anymore, but one can recover free energy dissipation. Proof. Summing (4.3c) over K ∈ T leads to: satisfies the volume-filling constraint n j=0 U p j,K = 1, for any K ∈ T . Then any solution U p+1 to the scheme (4.3) satisfies it as well.

K∈T m K (U p+1 i,K -U p i,K ) = -∆t p K∈T σ∈E K,int J p+1 i,Kσ , i = 0, . . . ,
Proof. Fix K ∈ T and sum (4.3c) over i = 0, . . . , n. Then on the one hand we have the sum of the crossdiffusion contributions

τ σ n i=0 n j=0 K ij U p+1 j,σ D Kσ U p+1 i -U p+1 i,σ D Kσ U p+1 j
that cancels thanks to the symmetry of the coefficients K ij . On the other hand, it is clear by construction that the Cahn-Hilliard term in J p+1 0,Kσ exactly compensates the sum of the Cahn-Hilliard terms in the n other fluxes. Therefore, one obtains

m K 1 ∆t p n i=0 U p+1 i,K -U p i,K = 0.
□ Lemma 4.4 (Weak positivity property). Let p ∈ N and assume that U p is nonnegative. Then any solution U p+1 to the scheme is nonnegative. If moreover U p is positive, then any solution U p+1 is positive as well.

Proof. We begin by proving nonnegativity. Fix i ∈ {0, . . . , n}. Reason by contradiction and assume that U p+1 i has a (strictly) negative minimum U p+1 i,K in the cell K ∈ T . The conservation scheme (4.3c) in the cell K gives that

m K U p+1 i,K -U p i,K ∆t p = - σ∈E K,int
J p+1 i,Kσ , and using that U p i is nonnegative we get that σ∈E K,int

J p+1 i,Kσ > 0. (4.4)
But using that, for any σ ∈ E K,int , U p+1 i,σ = 0 by definition (4.3a), several terms cancel in the fluxes formula (4.3d), and we get

σ∈E K,int   -τ σ 0≤j̸ =i≤n K ij U p+1 j,σ D Kσ U p+1 i   > 0, but since K ij , U p+1 j,σ , D Kσ U p+1 i , τ σ ≥ 0, we obtain a contradiction. Therefore U p+1 i ≥ 0. If U p
i is furthermore (strictly) positive then we can assume the minimum U p+1 i,K to be only nonnegative and still obtain (4.4), so the same reasoning gives U p+1 i > 0. □ Thanks to the previous lemma, the chain rule (4.3b) is always valid provided the initial condition is positive everywhere. We make this assumption in the following. As a consequence, we can define the discrete chemical potentials (see (1.7)-(1.8)) as follows: for any K ∈ T ,

µ p+1 i,K = ln(U p+1 i,K ), i = 1, . . . , n, (4.5)
and

µ p+1 0,K = ln(U p+1 0,K ) + W p+ 1 2 0,K . (4.6)
For all i = 0, . . . , n, we introduce µ p+1 i := (µ p+1 i,K ) K∈T . The discrete fluxes (4.3d) can thus be rewritten in the following entropic form (independently of the discretization formula (4.3e)):

J p+1 i,Kσ = -τ σ 0≤j̸ =i≤n K ij U p+1 i,σ U p+1 j,σ D Kσ µ p+1 i -µ p+1 j , i = 0, . . . , n. (4.7)
This latter formulation of the fluxes is at the core of the discrete free energy dissipation inequality. Let us define, according to (1.1), the following discrete free energy functionals: for all V = (V i,K ) 0≤i≤n,K∈T ,

E T (V ) := E conv,T (V ) + E conc,T (V ), E conv,T (V ) := n i=0 K∈T m K (V i,K ln(V i,K ) -V i,K + 1) + ε 2 σ∈Eint, σ=K|L τ σ |D Kσ V 0 | 2 , E conc,T (V ) := β K∈T m K V 0,K (1 -V 0,K ), (4.8)
with V 0 = (V 0,K ) K∈T . Their differentials at U = (U i,K ) 0≤i≤n,K∈T are given by

DE conv,T (U ) • V = n i=0 K∈T m K ln(U i,K )V i,K + ε σ∈Eint, σ=K|L τ σ D Kσ U 0 D Kσ V 0 , DE conc,T (U ) • V = β K∈T m K (1 -2U 0,K )V 0,K . (4.9)
Remark that, by convexity (resp. concavity), it holds, for p ∈ N:

E conv,T (U p+1 ) -E conv,T (U p ) ≤ DE conv,T (U p+1 ) • (U p+1 -U p ), E conc,T (U p+1 ) -E conc,T (U p ) ≤ DE conc,T (U p ) • (U p+1 -U p ). ( 4

.10)

We establish a discrete free energy dissipation inequality, as stated in the next proposition. Recall the definition (1.3) of the mobility matrix.

Proposition 4.5. Let p ∈ N and U p be (strictly) positive. Then any solution U p+1 to the scheme satisfies

E T (U p+1 ) -E T (U p ) + ∆t p σ∈Eint τ σ (D Kσ µ p+1 ) T M (U p+1 σ )D Kσ µ p+1 ≤ 0, (4.11)
where

U p+1 σ = U p+1 i,σ 0≤i≤n
. In particular, since M (U p+1 σ ) is always a positive semi-definite matrix, it holds

E T (U p+1 ) ≤ E T (U p ) ≤ E T (U 0 ).
Proof. Multiply equations (4.3c) by ∆t p µ p+1 i,K and sum over all species i = 0, ..., n and all cells K ∈ T :

n i=0 K∈T m K U p+1 i,K -U p i,K µ p+1 i,K = -∆t p n i=0 K∈T σ∈E K J p+1 i,Kσ µ p+1 i,K . (4.12)
On the left-hand side, we obtain using the definitions (4.5)-(4.6) of the discrete chemical potentials:

n i=0 K∈T m K U p+1 i,K -U p i,K µ p+1 i,K = K∈T n i=0 m K U p+1 i,K -U p i,K ln U p+1 i,K + K∈T m K U p+1 0,K -U p 0,K W p+ 1 2 0,K ,
where we obtain two different contributions that we want to identify to derivatives of the discrete energies (4.9). The first term identifies to the Boltzmann part, taken at U p+1 , against (U p+1 -U p ). Using the definition (4.3e) of W p+ 1 2 0,K and using the discrete integration by part formula (4.1), the second term can be expressed as:

K∈T m K U p+1 0,K -U p 0,K W p+ 1 2 0,K = -ε K∈T (U p+1 0,K -U p 0,K ) σ∈E K,int τ σ D Kσ U p+1 0 + β K∈T m K (U p+1 0,K -U p 0,K )(1 -2U p 0,K ), = ε σ∈Eint τ σ D Kσ U p+1 0 D Kσ U p+1 0 -U p 0 + β K∈T m K (U p+1 0,K -U p 0,K )(1 -2U p 0,K ),
where we identified the Cahn-Hilliard contributions of (4.9), respectively taken at U p+1 and U p , against U p+1 -U p . Putting everything together, we can identify the total derivative of the energy:

n i=0 K∈T m K U p+1 i,K -U p i,K µ p+1 i,K = DE conv,T (U p+1 ) • U p+1 -U p + DE conc,T (U p ) • U p+1 -U p ,
and we conclude from (4.12), using the convexity (resp. concavity) inequalities (4.10), that:

E T (U p+1 ) -E T (U p ) ≤ -∆t p n i=0 K∈T σ∈E K J p+1 i,Kσ µ p+1 i,K .
On the other hand, using the entropic form of the fluxes (4.7), the right-hand side reads:

-

n i=0 K∈T σ∈E K J p+1 i,Kσ µ p+1 i,K = n i=0 K∈T σ∈E K,int τ σ   0≤j̸ =i≤n K ij U p+1 i,σ U p+1 j,σ D Kσ (µ p+1 j -µ p+1 i )   µ p+1 i,K .
Using once again integration by parts, this is equal to

- σ∈Eint τ σ n i=0 0≤j̸ =i≤n K ij U p+1 i,σ U p+1 j,σ D Kσ (µ p+1 i -µ p+1 j )D Kσ µ p+1 i ,
and having in mind the expression of the mobility matrix (1.3), we finally obtain

- n i=0 K∈T σ∈E K J p+1 i,Kσ µ p+1 i,K = - σ∈Eint τ σ (D Kσ µ p+1 ) T M (U p+1 σ )D Kσ µ p+1 .
We have obtained (4.11). □

Numerical Simulations

The numerical scheme has been implemented in the Julia language. At each time step, the nonlinear system is solved using Newton's method with stopping criterion ∥U p,k+1 -U p,k ∥ ∞ ≤ 10 -10 and adaptive time stepping. We always consider the case n = 2 of three species and the cross-diffusion coefficients are chosen to be K 01 = K 10 = 0.2, K 12 = K 21 = 0.1, K 02 = K 20 = 1 (diagonal coefficients do not play any role). We study the dynamics for different values of ε and β and different initial conditions, with a particular focus on the stationary solutions obtained in the long-time asymptotics and the shape of the free energy over time. 5.1. One-dimensional simulations. We consider the domain (0, 1), a uniform mesh of 100 cells, a maximal time step ∆t 1 = 10 -3 and three initial profiles defined as smooth perturbations of a constant state: for x ∈ (0, 1), u 0 0 (x) = u 0 1 (x) = 1 4

(1 + κ cos(kπx)) , u 0 2 (x) = 1 -u 0 0 (x) -u 0 1 (x), where the perturbation is parametrized by its amplitude κ ∈ (0, 1) and frequency k ∈ N * , making sure that the box constraints in (1.4) are respected. Note that the mass is preserved by the perturbation, so that it holds m 0 = Ω u 0 0 dx = 0.25. Moreover, the Poincaré constant of the domain is C p = 1. We begin by illustrating Theorem 3.3 and we denote by U ∞ the discrete constant state defined according to (3.6). We first consider ε = 4, β = 1 so that condition (3.10) is satisfied. Starting from various initial conditions, we always observe a diffusive behaviour with exponential convergence to U ∞ , which confirms the globally stable behaviour. In Figure 1, we give some snapshots of the evolution and plot the discrete relative free energy RE T [U p | U ∞ ] = E T (U p ) -E T (U ∞ ) (see (4.8)) over the discrete time, starting from a specific initial condition and with time horizon T 1 = 10. For ε = 0.5, β = 2, (3.10) is not anymore satisfied, but the simulations behave similarly. Note that the convexity condition (3.3) on the energy reads

2 + ε 2 -β ≥ 0,
and is verified in this case. It seems that this condition is more determinant for the dynamics than (3.10).

We are now interested in the unstable regime and consider ε = 0.1, β = 10, so that both (3.10) and (3.3) are widely violated. In Figure 2, we display the results of the dynamics starting from the same initial condition as before and with a time horizon T 2 = 8. We observe convergence to a non-constant stationary solution with a diffuse segregation interface and exponential decrease to 0 of the quantity E T (U p ) -E T (U P T 2 ) over the discrete time. Let us make two remarks about this quantity: first, in contrast to the stable situation, we can only measure the difference with respect to the final solution U P T 2 , since we do not know a priori the limit of the dynamics. Second, since the limit is not homogeneous in space, the quantity E T (U p ) -E T (U P T 2 ) differs by a linear term from an approximation of the relative energy (3.7). In Figure 3, we display the results of the dynamics starting from an initial condition with higher frequency k = 2 over a time horizon T 3 = 2. We observe exponential convergence to another stationary solution, for which the free energy is smaller. Note that, in accordance with the results of Theorem 2.3, diffusion prevents the profiles from being exactly zero somewhere, and we always observe a δ > 0 such that all profiles are uniformly bounded between δ and 1 -δ. Finally, we provide numerical evidence that, in the general regime, the solution should converge to a local minimizer of the energy. Although the property of being a local minimizer cannot be easily verified numerically, one can use Theorem 2.4 and verify that the obtained numerical stationary solutions satisfy a discrete version of the optimality conditions (2.13)-(2.14). In Figure 4, we plot, for the two previous simulations, the l ∞ norm of the residual of this system over time and observe exponential convergence to 0, which indicates that the solution converges to a critical point of the energy. 5.2. Two-dimensional simulations. We consider the domain (0, 1) 2 , a uniform mesh of 150 2 squares and a maximal time step ∆t 2 = 5 × 10 -3 . We want to observe spinodal decomposition, so we pick initial profiles defined by random perturbation of constant states: for (x, y) ∈ (0, 1) 2 , u 0 0 (x, y) = 0.5 + 2κ η 0 (x, y) -1 2 , u 0 1 (x, y) = 0.4 + 2κ η 1 (x, y) -1 2 , u 0 2 (x, y) = 1 -u 0 0 (x, y) -u 0 1 (x, y), where, for any (x, y) ∈ (0, 1) 2 , η 0 (x, y) and η 1 (x, y) are independent noises drawn uniformly in (0, 1) and κ = 10 -2 . We choose the parameters ε = 10 -3 , β = 5. The results of the simulation are given in Figure 5. As expected, u 0 quickly separates from the two other species. Then on a slower time scale, the effect of cross-diffusion homogenizes u 1 and u 2 to the constant states. Finally, the coarsening process happens on a much slower time scale, minimizing the interface energy.

4. 3 .

 3 Elements of numerical analysis. Lemma 4.2 (Mass conservation). For all p ∈ N, it holdsK∈T m K U p+1 i,K = Ω u 0 i (x)dx, i = 0, . . . , n.

Figure 1 .

 1 Figure 1. Profiles along time in the globally stable case. ε = 4, β = 1, κ = k = 1.

  (a) t = 0 (b) t = 0.4 (c) T2 = 8 (d) E(u(t)) -E(u(T2))

Figure 2 .

 2 Figure 2. Profiles along time in the non-convex case. ε = 0.1, β = 10, κ = k = 1.

Figure 3 .

 3 Figure 3. Profiles along time in the non-convex case. ε = 0.1, β = 10, κ = 1, k = 2.

Figure 4 .

 4 Figure 4. Residual of the Euler-Lagrange system in the case k = 1 (left) and k = 2 (right).

Figure 5 .

 5 Figure 5. Spinodal decomposition successively at times t = 0, 0.06, 0.13, 0.49, 1.5.

  n, and this quantity is null because of cancellations on both sides of σ ∈ E int . The conclusion follows by induction and using (4.2). □ Lemma 4.3 (Volume-filling constraint). Let p ∈ N and assume that U p = U p

	i,K K∈T ,0≤i≤n

Acknowledgements. The authors acknowledge support from project COMODO (ANR-19-CE46-0002). GM acknowledges the support of Deutsche Forschungsgemeinschaft (DFG) via grant GZ: MA 10100/1-1 project number 496629752. GM is member of the Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). JFP thanks the DFG for support via the Research Unit FOR 5387 POPULAR, Project No. 461909888.

Appendix

Lemma 6.1. Let Ω ⊂ R d be a measurable bounded domain and let u ∈ L 1 (Ω) be such that u ≥ 0 almost everywhere on Ω. Let M := Ω u dx, n ∈ N * and m 1 , . . . , m n ∈ R + such that

Then there exist n measurable subsets Ω k ⊂ Ω for k = 1, . . . , n such that

Proof. Let x 1 ∈ Ω and consider the function f : R + → R + defined as:

Then, it can be easily seen that f is continuous using the Lebesgue convergence theorem, non-decreasing, such that f (0) = 0 and that there exists R > 0 such that for all r ≥ R, f (r) = M . This implies that there exists r 1 ≥ 0 such that f (r 1 ) = m 1 , and we defined Ω 1 = Ω ∩ B r1 . The other sets Ω 2 , . . . , Ω n can be constructed using exactly the same procedure by induction on the set Ω \ Ω 1 . □ Email address: jan-f.pietschmann@uni-a.de


Espen Gaarder
email: espenhaug@mac.com

Stéphane Wojnow
email: wojnow.stephane@gmail.com

The Blackbody CMB Temperature, Luminosity, and Their Relation to Black Hole Cosmology, Compared to Bekenstein-Hawking Luminosity

Keywords:

We use the Stefan-Boltzmann law to derive and investigate what we will call the black hole CMB temperature luminosity. Interestingly, the CMB luminosity is always the same for every black hole at a luminosity distance equal to the Schwarzschild radius, which is fully valid in black hole cosmology, a topic that is actively discussed to this day. The CMB black hole luminosity is always L cmb = T 4 p l 2 p 256⇡ 3 , which can also be written in the form

, regardless of the mass of the black hole. This is in contrast to the Bekenstein-Hawking luminosity, which is linked to the Hawking temperature, that is mass dependent and equal to L Hw = c 6 15360⇡ G2 M 2 .

Introduction

Black holes are considered near-perfect black bodies because they absorb all the radiation that falls on them. It has been suggested that black holes emit blackbody radiation known as Hawking radiation, which depends on the mass of the black hole. However, Hawking radiation has likely never been observed, possibly due to the di culty in detecting it. Nevertheless, we will briefly examine it here, as the equation and derivation for Bekenstein-Hawking luminosity share many similarities with what we will coin as: "CMB luminosity" or alternatively "internal black hole luminosity". We will mathematically explore Hawking luminosity in section 2 and then CMB luminosity in section 3 and 4.

Additionally to Hawking temperature that is assumed to be a black body temperature, there is another black body temperature associated with black holes and cosmology in general including the ⇤-CDM model, even if the ⇤-CDM model possibly will fail short for what we will 1 discuss in this paper, so called R H = ct cosmological models and particular growing blackhole R H = ct will much better fit what we will describe, see [START_REF] Tatum | How a realistic linear r h = ct model of cosmology could present the illusion of late cosmic acceleration[END_REF][START_REF] John | r h = ct and the eternal coasting cosmological model[END_REF][START_REF] John | Generalized Chen-Wu type cosmological model[END_REF][START_REF] John | Comparison of cosmological models using bayesian theory[END_REF][START_REF] Melia | The r h = ct universe[END_REF][START_REF] Melia | The r h = ct universe without inflation[END_REF][START_REF] Melia | The linear growth of structure in the r h = ct universe[END_REF]. In various cosmological models, it is assumed that the Hubble sphere could be a black hole. This idea dates back to at least 1972, as suggested by Pathria [START_REF] Pathria | The universe as a black hole[END_REF], and was later revisited by Stuckey [START_REF] Stuckey | The observable universe inside a black hole[END_REF] in 1994. Although the ⇤-CDM model community may overlook the possibility that the Hubble sphere is a black hole, many researchers continue to work on or discuss black hole cosmology to this day, as shown by [START_REF] Pop | The universe in a black hole in Einstein-Cartan gravity[END_REF][START_REF] Tatum | A Heuristic Model of the Evolving Universe Inspired by Hawking and Penrose in New Ideas Concerning Black Holes and the Universe[END_REF][START_REF] Akhavan | The universe creation by electron quantum black holes[END_REF][START_REF] Lineweaver | All objects and some questions[END_REF]. Assuming, for a moment, that the Hubble sphere is a black hole, we then have a significant amount of information about the interior of a black hole. Furthermore, there is another type of black body radiation associated with the interior of black holes. The CMB temperature is often linked with a black body temperature; for example, Muller et al. in [START_REF]A precise and accurate determination of the cosmic microwave background temperature at z = 0.89[END_REF] state (see also [START_REF] Dhal | Investigation on CMB monopole and dipole using blackbody radiation inversion[END_REF]):

"Observations with the COBE satellite have demonstrated that the CMB corresponds to a nearly perfect black body characterized by a temperature T 0 at z = 0, which is measured with very high accuracy, T 0 = 2.72548 ± 0.00057k."

In recent years, there has been tremendous progress in understanding the CMB (Cosmic Microwave Background) temperature from a practical perspective, particularly in terms of more accurate measurements, as cited in [START_REF] Fixsen | The temperature of the cosmic microwave background at 10 GHz[END_REF][START_REF] Fixsen | The temperature of the cosmic microwave bacground[END_REF][START_REF] Noterdaeme | The evolution of the cosmic microwave background temperature[END_REF][START_REF] Dhal | Calculation of cosmic microwave background radiation parameters using cobe/firas dataset[END_REF]. Additionally, there has been a deeper theoretical understanding of the CMB temperature and its relationship to black hole cosmology. In section 3, we will utilize this knowledge to derive what we will refer to as the CMB luminosity. This luminosity should, in theory, apply to the interiors of all or at least most black holes, even though there are likely to be modifications for di↵erent metrics.

Bekenstein-Hawking luminosity of a black hole

The luminosity of of a black body is given from the Stefan-Boltzman law [START_REF] Stefan | Über die beziehung zwischen der wärmestrahlung und der temperatur[END_REF][START_REF] Boltzmann | Ableitung des stefanschen gesetzes, betre↵end die abhängigkeit der wärmestrahlung von der temperatur aus der electromagnetischen lichttheori[END_REF] and is:

L = 4⇡R 2 T 4 (1)
where is the Stefan-Boltzman constant, T is the black body temperature, and R is the luminosity distance.

For luminosity directly related to a black hole, we will set the luminosity radius equal to the black hole's event horizon radius and the temperature T to the black hole's temperature. There are actually two temperatures associated with black holes. One is the Bekenstein-Hawking temperature, as given by [START_REF] Hawking | Black hole explosions[END_REF], which is described as follows:

T Hw = ~g ck b 2⇡ = ~GM r 2 s ck b 2⇡ = ~c k b 4⇡R s (2)
Next we input this in the Stefan-Boltzman luminosity equation:

L Hw = 4⇡R 2 s T 4 Hw L Hw = 4⇡R 2 s ✓ ~c k b 4⇡R s ◆ 4 L Hw = ~4c 4 k 4 b 64⇡ 3 R 2 s ( 3
)
and since the Stefan-Boltzman constant is =

2⇡ 5 k 4 b 15c 2 h 3
, if we insert this in above equation we get:

L Hw = 2⇡ 5 k 4 b 15c 2 h 3 ~4c 4 k 4 b 64⇡ 3 R 2 s L Hw = ~c2 3840⇡R 2 s . (4)
Further since R s = 2GM c 2 we can also re-write this as

L Hw = ~c6 15360⇡G 2 M 2 .
(

) 5
This is the well-known Bekenstein-Hawking luminosity of a black hole. Please note that the Bekenstein-Hawking luminosity scales with the inverse square law of both the black hole's radius and its corresponding mass.

CMB luminosity of a black hole is the same for every black hole

Haug and Wojnow [START_REF] Haug | How to predict the temperature of the CMB directly using the Hubble parameter and the Planck scale using the Stefan-Boltzman law[END_REF] has derived the CMB temperature from the Stapfan-Boltzman law and got the formula:

T cmb = T p 8⇡ s 2l p R H ⇡ 2.72k. ( 6
)
where T p is the Planck [START_REF] Planck | Der Königlich Preussischen Akademie Der Wissenschaften[END_REF][START_REF] Planck | Vorlesungen über die Theorie der Wärmestrahlung[END_REF] temperature

T p = 1 k b q ~c5 G = mpc 2
k b (see also [START_REF] Unnikrishnan | Standard and derived Planck quantities: selected analysis and derivations[END_REF]). This formula holds for any Schwarschild black hole and in more general forms is given by

T cmb = T p 8⇡ s 2l p R s ⇡ 2.72k. ( 7
)
This correspond also to the formula given by Tatum et al. [START_REF] Tatum | The basics of flat space cosmology[END_REF], who likely is the first to obtained a formula for predicting the CMB temperature using heuristic methods that we soon also will discuss. Furthermore, Haug and Tatum [START_REF] Haug | The hawking Hubble temperature as a minimum temperature, the Planck temperature as a maximum temperature and the CMB temperature as their geometric mean temperature[END_REF] derived the CMB temperature from a very general geometric mean of the shortest and longest possible wavelengths in a black hole. In addition, Haug [START_REF] Haug | Planck quantized bending of light leads to the CMB temperature[END_REF] has derived the same formula from quantized light bending related to black holes.

This can be generalized for any Schwarzschild [START_REF] Schwarzschild | über das gravitationsfeld einer kugel aus inkompressibler flussigkeit nach der einsteinschen theorie[END_REF] black hole:

T cmb = T Hw s R s 2l p = T p 8⇡ s 2l p R s ⇡ 2.72k. ( 8
)
where T cmb is the black body temperature associated with the black hole, the temperature on the inside of the black hole. Therefore, the term 'cosmic microwave background' might be misleading since smaller black holes will have a wavelength spectrum for their temperature that consists of much shorter wavelengths than microwaves, so alternatively we could just have called it the inside temperature of a black hole. The luminosity for any black hole, derived from the internal black body temperature, is then given by:

L cmb = 4⇡R 2 s T p 8⇡ s 2l p R s ! 4 L cmb = 4⇡R 2 s T 4 p 8 4 ⇡ 4 4l 2 p R 2 s L cmb = 4⇡ T 4 p 4l 2 p 8 4 ⇡ 3 L cmb = T 4 p l 2 p 256⇡ 3 . (9)
What is remarkable here is that we observe the black body luminosity as constant no matter the size of the black hole, and furthermore, it is only a function of the Planck scale. The black hole mass is not relevant. We can conclude that all black holes have the same CMB luminosity, despite the CMB temperature scaling with 1 p Rs for di↵erent black holes. Therefore, the temperature inside each black hole varies between black holes with di↵erent size, but the CMB luminosity at a distance of the Schwarzschild radius for the black hole of interest remains the same for all black holes.

Furthermore, since the Stefan-Boltzmann constant is given by =

L cmb = ⇡ 2 k 4 b l 2 p 1920c 2 h 3 T 4 p L cmb = T 4 p k 4 b l 2 p 15360c 2 ⇡~3 . ( 10
)
Alternatively we can derive the black hole CMB luminosity by using the first formula presented by Tatum et. al [START_REF] Tatum | The basics of flat space cosmology[END_REF] and later derived by Haug and Tatum [START_REF] Haug | CMB temperature related to cosmological red-shift[END_REF] from a geometric mean approach, we get:

L cmb = 4⇡R 2 s ~c k b 4⇡ p R s 2l p ! 4 L cmb = 4⇡R 2 s ~c k 4 b 4 4 ⇡ 4 R 2 s 4l 2 p L cmb = ~4c 4 k 4 b 256⇡ 3 l 2 p ( 11
)
as the Stefan-Boltzman constant =

L cmb = 2⇡ 5 k 4 b 15c 2 h 3 ~4c 4 k 4 b 256⇡ 3 l 2 p L cmb = ~c2 15360⇡l 2 p L cmb = ~c6 15360⇡G 2 m 2 p ⇡ 7.50338 ⇥ 10 47 W ( 12
)
where m p is the Planck mass. Equations 9, 10 and 12 are ultimately the same, it is just di↵erent ways to express the same formula. What is important is that the CMB luminosity of any black hole is the same. This means every black hole has a luminosity of 7.50338⇥10 47 W . This is also fully consistent with the observed relation between the cosmological red-shift and the CMB temperature that is assume to follow of T t = T 0 (1 + z) and z = Tt T 0 1 (see for example [START_REF] Riechers | Microwave background temperature at a redshift of 6.34 from h 2 o absorption[END_REF]). This mean we can replace this back into the formula and we get

L cmb = 4⇡R 2 t (T 0 (1 + z)) 4 L cmb = 4⇡R 2 t T 4 0 ✓ 1 + T t T 0 1 ◆ 4 L cmb = 4⇡R 2 t T 4 0 ✓ T t T 0 ◆ 4 L cmb = 4⇡R 2 t T 4 t . (13)
As R t = R s for a Schwarzschild black hole in a R s = R H = ct type cosmological growing black hole models then

T t = ~c k b 4⇡ p Rt2lp
as pointed out by Tatum et. al [START_REF] Tatum | The basics of flat space cosmology[END_REF] and Haug and Tatum [START_REF] Haug | The hawking Hubble temperature as a minimum temperature, the Planck temperature as a maximum temperature and the CMB temperature as their geometric mean temperature[END_REF] this leads to exactly the luminosity equations as we have derived above. That is for the first time we have demonstrated that any black hole has a the same CMB luminosity if measured from the Schwarzschild radius distance of the black body. This naturally when the luminosity distance is set equal to the Schwarzschild radius.

The CMB temperature in yet another perspective

We can naturally also solve the standard Stefan-Boltzman luminosity equation with respect to the temperature T , this gives the well known formula:

T = 4 r L 4⇡R 2 . ( 14
)
As we have found that the black hole luminosity is the same for every Schwarzschild black hole and is given by:

L cmb = ~c6 15360⇡G 2 m 2 p
, this means we for any Schwarzschild black hole we must have a CMB temperature (on the inside) of:

T cmb = 4 s L cmb 4⇡R 2 s T cmb = 4 v u u t ~c6 15360⇡G 2 m 2 p 4⇡R 2 s T cmb = 4 s ~c6 61440⇡ 2 G 2 m 2 p R 2 s T cmb = 1 p ⇡Gm p R s 4 r ~c6 61440 . ( 15
)
In the case we set the Schwarzschild radius to the Hubble radius we get R s = R H which is the case in the critical Friedmann [START_REF] Friedmann | Über die krüng des raumes[END_REF] universe we get:

T cmb = 1 p ⇡Gm p R H 4 r ~c6 61440 ⇡ 2.72 +0.082 0.069 K (16)
we have here used a Hubble parameter in calculation of R H = c H 0 given by the recent published study by Kelly and et. al [33] of 66.6 +4. 1 3.3 (km/s)/M pc. However di↵erent Hubble parameter studies give considerably di↵erent Hubble constant values, so the exact value will naturally vary with input, see for example [START_REF] Mukherjee | First measurement of the Hubble parameter from bright binary black hole gw190521[END_REF][START_REF] Sedgwick | The e↵ects of peculiar velocities in SN Ia environments on the local h0 measurement[END_REF][START_REF] Riess | A comprehensive measurement of the local value of the Hubble constant with 1 km s 1 mpc 1 uncertainty from the Hubble space telescope and the sh0es team[END_REF].

We can naturally also re-write equation 16, by replacing R s with R s = 2GM c 2 this gives:

T cmb = 4 s ~c6 61440⇡ 2 G 2 m 2 p 2GM c 2 2 T cmb = 4 s ~c10 245760⇡ 2 G 4 m 2 p M 2 T cmb = 1 G p ⇡m p M 4 r ~c10 245760 . (17)
If one now input the critical Friedmann mass of the universe M = M c = c 3 2GH 0 , we get again T cmb ⇡ 2.72K. Let us next input the for the Stefan-Boltzman constant: =

T cmb = 1 G p ⇡m p M 4 v u u t ~c10 245760 2⇡ 5 k 4 b 15c 2 h 3 T cmb = 1 G p ⇡m p M 4 v u u t ~c12 4096 ⇡ 2 k 4 b ~3 T cmb = ~c3 k b 8⇡G p Mm p = ~c k b 4⇡ p R s 2l p (18)
This is the same formula as heuristically suggested by Tatum et. al [START_REF] Tatum | The basics of flat space cosmology[END_REF], which has recently been derived from the Stefan-Boltzmann law by Haug and Wojnow [START_REF] Haug | How to predict the temperature of the CMB directly using the Hubble parameter and the Planck scale using the Stefan-Boltzman law[END_REF], and by the geometric mean approach by Haug and Tatum [START_REF] Haug | The hawking Hubble temperature as a minimum temperature, the Planck temperature as a maximum temperature and the CMB temperature as their geometric mean temperature[END_REF], and also by quantization of light bending by Haug [START_REF] Haug | Planck quantized bending of light leads to the CMB temperature[END_REF]. Despite possibly, and even likely, not fitting into the ⇤-CDM model, we think it can be no coincidence that the CMB formula first heuristically suggested by Tatum et. al can just be a coincidence; it predicts the right CMB temperature very accurately compared to measurments and have now been derived from a series of di↵erent complementary angles. We think cosmologies would likely benefit from taking notice of it. That said we naturally do not claim to have all the answers. There is certainly more to investigate here in relation to both theoretical CMB temperature and CMB luminosity; for example, there are multiple metrics that can be applied to black holes that should be investigated, such as the Reissner-Nordström [START_REF] Reissner | Über die eigengravitation des elektrischen feldes nach der Einsteinschen theorie[END_REF][START_REF] Nordström | On the energy of the gravitation field in Einstein's theory[END_REF] metric, the Kerr [START_REF] Kerr | Gravitational field of a spinning mass as an example of algebraically special metrics[END_REF] metric, the Kerr-Newman [START_REF] Newman | Metric of a rotating, charged mass[END_REF][START_REF] Newman | Note on the Kerr spinning-particle metric[END_REF] metric, and the recent Haug-Spavieri [START_REF] Haug | Mass-charge metric in curved spacetime[END_REF] metric, as well as other metrics used in cosmology and for black holes.

Tatum et. al [START_REF] Tatum | High precision Hubble constant determinations based upon a new theoretical relationship between CMB temperature and H 0 . Hal archive[END_REF] has clearly demonstrated that one can also use the CMB temperature formula, and when solved for H 0 , it can dramatically reduce the standard deviation in uncertainty from other methods. This is because we now, for the first time, have a much deeper theoretical understanding of the relation between CMB temperature and the Hubble parameter. The ⇤-CDM model has, on the other hand, not developed any method to theoretically predict the CMB temperature now. All they have is the formula to predict the CMB temperature back in time based on the formula T t = T 0 (1 + z), which is also fully consistent with this work.

There are also other challenges that could be related to this and that could possibly be resolved with a deeper understanding of CMB, namely the Hubble tension [START_REF] Valentino | In the realm of the Hubble tension -a review of solutions[END_REF][START_REF] Krishnan | Does Hubble tension signal a breakdown in FLRW cosmology? Classical and Quantum Gravity[END_REF]. In our model, we must have z = q R H Rt 1 = T 0 Tt 1, as recently pointed out by Haug and Tatum [START_REF] Haug | CMB temperature related to cosmological red-shift[END_REF].

Black hole internatal radiation pressure

The total energy per unit of time per unit of surface area that is radiated by a black body maintained at a temperature T is according to Stefan-Boltzman law the radiation pressure given by:

P = 4 T 4 3c = L 3⇡R 2 c (19)
Next we will derive the black hole internal radiation pressure that also can be called the CMB radiation pressure as it can be derived from what we called the CMB luminosity:

P cmb = = L 3⇡R 2 s c P cmb = = T 4 p l 2 p 768⇡ 4 R 2 s c (20)
So the black hole radiation pressure follow the inverse square law rule, unlike the CMB luminosity that is the same for every black hole. Equation 20 can also be re-written as:

P cmb = T 4 p l 2 p c 3 768⇡ 4 G 2 M 2 (21)
In the special case where M is the critical mass of the Friedmann universe M = M c = c 2 2GH 0 we get a radiation pressure of:

P cmb = T 4 p l 2 p c 3 3072⇡ 4 G 2 M 2 c ⇡ 1.3909 ⇥ 10 14 J • m 3 (22)
If we replace the Stefan-Boltzmann constant with its components: =

Conclusion

We have demonstrated that there is another form of black body radiation or luminosity associated with every black hole. We have named this the CMB luminosity, even though it applies to any black hole. The CMB luminosity is identical for every black hole and is always given by L cmb = T 4 p l 2 p 256⇡ 4 , which can also be expressed as L cmb = c 6 15360⇡ G2 m 2 p . It depends solely on the Planck scale and is not a↵ected by the mass of the black hole. This discovery can have significant implications for both cosmology and other aspects of black hole research.

2⇡ 5 k 4 b

 4 15c 2 h 3 , we can also rewrite equation 9 as:

2⇡ 5 k 4 b 15c 2 h 3

 43 and we can now insert this into the equation above and we get

2⇡ 5 k 4 b 15c 2 h 3

 43 and we get:

2⇡ 5 k 4 b 15c 2 h 3 , we get: P cmb = 2⇡ 5 k 4 b 15c 2 h 3 T 4 p l 2 p c 3 3072⇡ 4 G 2 M 2 P cmb = 2⇡ 5 k 4 b 15c 2 h 3 ~4c 4 l 4 p k 4 b l 2 p c 3

 434342243 3072⇡ 4 G 2 M 2 P cmb = ~c5 184320⇡ 2 G 2 M 2 l 2 p (23)

Declarations

Conflict of interest

The authors declares no conflict of interest.