{'fraction_non_alphanumeric': 0.05618213844414649, 'fraction_numerical': 0.020427478507768753, 'mean_word_length': 3.9918635170603674, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
TREND TO EQUILIBRIUM FOR SYSTEMS WITH SMALL CROSS-DIFFUSION

This paper presents new analytical results for a class of nonlinear parabolic systems of partial different equations with small cross-diffusion which have been proposed to describe the dynamics of a variety of large systems of interacting particles. Under suitable assumptions, we prove existence of classical solutions and we show exponential convergence in time to the stationary state. Furthermore, we consider the special case of one mobile and one immobile species, for which the system reduces to a nonlinear equation of Fokker-Planck type. In this framework, we improve the equilibration result obtained for the general system and we derive L∞-bounds for the solutions in two spatial dimensions. We conclude by illustrating the behaviour of solutions with numerical experiments in one and two spatial dimensions.


Introduction
1.1. Background and motivation. In this paper we focus on a relatively general class of nonlinear crossdiffusion systems with sufficiently small off-diagonal diffusion terms. The smallness assumption ensures that the system we consider is "close", in a suitable sense, to a linear, decoupled parabolic system. This allows us to adapt techniques developed in the theory of linear, parabolic systems of partial differential equations (PDEs). More specifically, we consider a parabolic system of PDEs which, in its most general form, reads (1) ∂u i ∂t − α,β,j ∂ ∂x α A αβ ij (x, u) ∂u j ∂x β − B α ij (x, u)u j = 0, for 1 ≤ i, j ≤ m, and 1 ≤ α, β ≤ d. Here m ≥ 1 represents the number of components (or species) and d ∈ {1, 2, 3} the number of space dimensions. We denote by u i the i-th species. We will specify the assumptions on the diffusion tensor A and the drift B in Section 2. Cross-diffusion systems arise in multiple contexts in Physics, Life Sciences and Social Sciences; in particular, they have been derived as formal macroscopic limits of several microscopic models describing multispecies systems in the presence of finite volume effects, size exclusion or joint population pressures (see, for example, [12,13,27,28]). Such effects enhance the macroscopic diffusion as each species tends to avoid overcrowding. One classical example is the derivation from a lattice based microscopic description, see, for instance, [8]. A subsequent formal Taylor expansion of the associated master equation, as considered in [15,14], yields systems of cross-diffusing species. Another derivation was performed in [12,11], where the authors derived a cross-diffusion system from an underlying stochastic microscopic representation using the method of matched asymptotic expansions. We highlight that the two approaches yield different continuum models, however, they share the property of having small cross-diffusion terms. The smallness assumption is justified since the cross-diffusion terms are of the same order of magnitude as the microscopic particle size (see, e.g., Examples 1.1-1.3 in [3]).

Gradient flow techniques.
In recent years, gradient flow methods have been successfully employed to study certain families of cross-diffusion systems, see, among others, [23,20,19]. Techniques such as the boundedness-by-entropy principle introduced in [24], provide a mathematical framework to ensure existence and uniqueness of solutions to general nonlinear cross-diffusion systems that exhibit a gradient flow structure. Another prominent tool developed in this context, is the Bakry-Emery strategy (see [6]) which allows to establish convex Sobolev inequalities and compute exponential decay rates in the trend towards equilibrium. The cross-diffusion systems mentioned above lack a full gradient flow structure (in the Wasserstein sense) in certain parameter ranges, even though the underlying microscopic system possesses a natural one. This lack is caused by approximations made due to the finite volume effects. From the perspective of modelling, a first connection between the large deviations of stochastic particle systems and macroscopic Wasserstein gradient flows was established in [2] without finite volume constraints.
Unfortunately, the results and techniques mentioned above cannot be applied at this stage to the systems of PDEs we want to study, nevertheless in this work we present an alternative strategy that relies on the smallness of the off-diagonal cross diffusion terms.
1.3. Numerical methods. The numerical discretisation of non-linear cross-diffusion systems, especially the development of structure preserving numerical schemes received considerable interest in the last years. These methods are designed in such a fashion that they preserve important physical and structural features such as positivity, conservation of mass or the dissipation of the associated entropy. Owing to the fact that Wasserstein gradient flows are posed in the set of probability measures, conservation of mass (or probability) of solutions is an important physical feature and finite volume discretisations are a natural framework to guarantee this. In addition, in recent years, several advances have been made in designing flux approximations that are in agreement with the energy dissipation. Bessemoulin-Chatard and Filbet, see [7], were among the first to present a finite volume method for nonlinear degenerate parabolic equations, which resolved the long-time behaviour correctly. Based on their scheme, different finite volume schemes have been proposed for systems in the last years, see, for example, [17,16]. Other numerical approximations are based on the underlying variational Wasserstein gradient flow structure. These so-called variational schemes are often restricted to one space dimension, as the computational complexity of computing the Wasserstein distance in higher space dimension is significant, see, for instance, [18]. Also, convergence results are, to the best of our knowledge, restricted to one spatial dimension; cf. [26]. 1.4. Summary of the main results. We consider a family of systems of PDEs with small cross-diffusion terms (namely system (1)) for which existence and uniqueness of solutions were investigated in [3], as recalled in Proposition 2.1. We extend the analysis presented in [3] by showing in Theorem 2.2 that (under suitable assumptions) the solutions are global-in-time and classical. Furthermore, we provide insights into the equilibration behaviour of such cross-diffusion systems in Theorem 3.2.
A special instance of system (1), to which we shall return in Section 4, was analysed in [10] using gradient flow techniques; existing results for such system and its stationary states are presented in Proposition 4.1. These results were further improved in [9] for a reduced nonlinear PDE of Fokker-Planck type, to which the system reduces in the special case of one mobile species and one immobile species (also considered in Section 4). The reduced PDE (namely problem (26)) can be interpreted as the diffusion of hard-core interacting particles through a domain with obstacles distributed according to the given immobile function. Exponential convergence to equilibrium for the scalar case is somehow simpler and it is proved in Theorem 4.2. Furthermore, in Theorem 4.3 we establish sharper L ∞ -bounds for the special case of the reduced PDE in two spatial dimensions.
We conclude by illustrating our analytic results with numerical simulations. In doing so we construct a finite volume scheme in one and two spatial dimensions, investigate the exponential rate of convergence to equilibrium numerically and illustrate the dynamics of the system in physically relevant scenarios. The impact of the immobile species is clearly visible in that the immobile species diminishes the mobility of the mobile species; cf. Section 5.
The rest of paper is organised as follows: we start by recalling well-posedness results for general systems with small cross diffusion and prove existence of global classical solutions in Section 2. Section 3 states the main exponential convergence result for the full system. In Section 4 we consider a specific example of an asymptotic gradient flow system, describing a mobile species diffusing in the presence of an immobile given species. In particular we are able to prove L ∞ -bounds for such scalar, nonlinear equations of Fokker-Planck type. Numerical simulations illustrating the analytic results are presented in Section 5.

Analysis of systems with small cross-diffusion
We start by introducing the cross-diffusion system under consideration and by discussing analytic properties of solutions. We recall existing well-posedness results and establish existence of classic solutions under certain conditions on the initial datum.
2.1. Review of existing well-posedness results. We consider the following general initial-boundary value problem involving a parabolic system of PDEs with small cross-diffusion: for 1 ≤ i, j ≤ m, and 1 ≤ α, β ≤ d. We indicate by ν the outward normal of Ω, m represents the number of components and d ∈ {1, 2, 3} corresponds to the number of space dimensions. Throughout the paper we denote the i-th species by u i and the related initial distribution by u 0 i ≥ 0. Note that, due to the no-flux boundary condition, the total "mass" (or number of particles in the microscopic case) of each species is conserved.
Global existence and regularity of solutions for system (2) was proved in [3]. An important assumption for their results was that the system is close to a diagonal, decoupled, linear problem. Such "reference problem" is given by the weak formulation of for 1 ≤ i, j ≤ m, and 1 ≤ α, β ≤ d.
Let us recall the precise assumptions we make (notice that they are analogous to those in [3]).

H1
We assume that D αβ ij ∈ C 1 (Ω) is symmetric in the space indices α and β and that it is "diagonal" in the component indices i and j, i.e., (4) D αβ ij = 0, for all i = j. Furthermore, we suppose that there exist two constants Λ ≥ λ > 0 such that for every x ∈ Ω and ξ ∈ R d , it holds that for any i = 1, . . . , m. H2 We make the following structural assumptions: where δ is a small parameter such that 0 ≤ δ ≤ δ 0 , for a suitable constant δ 0 depending on u 0 , A, B and Ω (see [3] for further details on the smallness assumption). H3 The dependence on u of the nonlinear terms is (at least) of class C 2 : H4 We assume that, for i = 1, . . . , m, the potentials V i : H5 The domain Ω is bounded, connected and sufficiently smooth (e.g. of class C 2 ). The initial datum u 0 in (2) is assumed to be non-negative and to belong to H 2 (Ω) satisfying the following compatibility condition on the boundary, ∂Ω, for i = 1, . . . , m: Before recalling the well-posedness result from [3] we introduce the following notation.
Definition 2.1. We denote by W (Q T ) the Banach space of functions with two spatial weak derivatives taking values in L 2 (Ω) continuously in time, and one time derivative in L 2 (0, T ; H 1 (Ω)), that is, Note that, given T > 0, we denote the parabolic cylinder by Q T = (0, T ) × Ω. Also, we write H 2 (Ω) for H 2 (Ω; R m ), and similarly for other spaces.

2.2.
Existence of global classical solutions. We are now going to show that solutions (in the sense of Definition 2.1) of problem (2) are, indeed, classical solutions if the initial data is sufficiently smooth. The arguments in this section follow closely those presented in [4] for a different class of parabolic systems. We start by showing that classical solutions exist for a short time (see Proposition 2.3) and then extend their maximal time interval of definition to the positive half-line using uniform estimates (see Proposition 2.4). This yields the main result of this section: Theorem 2.2. Let hypotheses H1-H5 hold and assume that u 0 is (at least) of class C 2 (Ω). Then problem (2) admits a unique, bounded, global, classical solution in the space C 0 ([0, ∞); C 2 (Ω)) ∩ C 1 ([0, ∞); C 0 (Ω)).
Remark 1. We denote by J(u 0 ) ⊆ [0, ∞) the maximal interval of definition of classical solutions of (2) given an initial datum u 0 .
Proposition 2.4 (Criterion for global existence, [5]). Let all assumptions of Proposition 2.2 hold and consider a solution u of problem (2). If an exponent ε ∈ (0, 1) (not depending on time) exists, such that for any T > 0, Remark 2 (Notation). Proposition 2.4 is a rephrasing of Theorem 2 in [5] with θ = 0. Notice also that the "convective term", indicated by f in the notation of [5], is affine in the gradient in our case. Additionally, we do not use the notation BU C ε to denote the space of "bounded, uniformly ε-Hölder continuous functions" but simply write that the Hölder exponent does not depend on time.
The following interpolation results provides sufficient conditions for the criterion in Proposition 2.4 to be satisfied.
and, in turn, Proof. Thanks to the higher order extensions for Sobolev functions, we can define f on a larger cylindrical domain R ⊆ R d+1 containing Q T . Introducing a cut-off function, we further extend f to the whole space ensuring sufficiently fast decay at infinity. Let us call g such an extension. We observe that the norm This means that ω + κ 2 ĝ ∈ L 2 (R d+1 ) and we obtain ω r κ s |ĝ| ≤ ω + κ 2 r+ s 2 |ĝ|. Thus we obtain the desired fractional Sobolev regularity provided that r + s 2 ≤ 1. We are also using the following inequality, relating the norms of g and f , g X ≤ 2 f X . Then the Hölder regularity follows from standard embeddings for fractional Sobolev spaces (see, e.g., [21]). In particular, for r, s > 1 2 , we take η = r − 1 2 and θ = s − 1 2 .
Proof of Theorem 2.2. Thanks to Proposition 2.3, we know that classical solutions exist for short times and that, as explained in [1], they may be extended to a maximal interval of existence denoted by J(u 0 ) by standard methods. In order to show that such solutions exist for arbitrarily large time we are going to use the criterion provided by Proposition 2.4. In particular, we need Hölder continuity of u with respect to time, as well as a uniform L ∞ -bound in the space variable. Thanks to Proposition 2.1, we know that u ∈ W with a bound that is uniform in time. This implies that u ∈ X and hence we apply Lemma 2.5 in order to obtain obtain uniform Hölder estimates. Thus Proposition 2.4 allows us to conclude the proof.

Trend to equilibrium
In this section we investigate the equilibration behaviour of solutions to system (2). We show that they converge exponentially fast to the stationary state in the L 2 -norm if assumptions H1-H5 are satisfied. We choose a more compact notation in the following. Let u = (u 1 , . . . u m ), then system (3) can be written as and the corresponding stationary problem: where we used the compact notation Recall that Ω is a bounded domain of class C 2 in R d (with d = 1, 2, 3), and 0 ≤ δ ≤ δ 0 ≤ 1, u 0 i ∈ C 2 (Ω) and V i ∈ C 2 (Ω).
Before stating and proving the main theorem of this section, we introduce some notation and simple inequalities that we will use in the proof.

Remark 3. Given a function
in particular, given two functions w, z : R d → R m such that w, z ≤ M for some M > 0 and for a.e. x ∈ R d , then we have (by Taylor expansion) Remark 4. We shall use the following notation: For completeness, we present the following easy modification of the standard Poincaré inequality. We recall that, if Ω is convex, we may take C P ≤ diam(Ω)π −1 . Lemma 3.1 (A Poincaré-type inequality). Consider a bounded, connected, Lipschitz domain Ω, a function f ∈ H 1 (Ω) and a subset S ⊂ Ω with positive Lebesgue measure. Suppose that S f (x)dx = 0, then there exists a constant C P > 0 (depending on Ω) such that Proof. We argue by contradiction. If the statement does not hold, then for each n > 0 there exists a function f n ∈ H 1 (Ω) such that S f n (x)dx = 0 and f n L 2 (Ω) > n ∇f n L 2 (Ω) .
Without loss of generality, we can assume that f n L 2 (Ω) = 1 (by homogeneity of the norm). This implies that ∇f n L 2 (Ω) < n −1 . Therefore, as n → ∞, the sequence converges (up to a subsequence) strongly in L 2 (Ω) to a limit function f ∞ ∈ H 1 (Ω) which satisfies the following properties The first two conditions above imply that f ∞ = 0, which contradicts the third.

Remark 5.
Under the same assumptions of Lemma 3.1, we have Notice that the condition C P ∇V ∞ < 1 gives a restriction on the size of Ω (since the Poincaré constant is proportional to diam(Ω)) or, alternatively, on the L ∞ -norm of the drift term ∇V . We will use this relation in estimate (18).
We now state the main result of this section: Theorem 3.2 (Exponential convergence). Let assumptions H1-H5 hold and suppose that u ∈ W (Q T ) and u * ∈ W 1,∞ (Ω) are solutions of problem (10) and (11) respectively. Let us denote by M > 0 a constant not depending on T such that where Γ M is given by (19). Then the following inequality holds for any t 1 , t 2 ≥ 0: Proof. Let h := u − u * , then h satisfies the following system (with no-flux boundary conditions) in Ω: Testing system (16) where Next we estimate the order δ terms in (17) and define In particular, making extensive use of Remarks 3, 4, and 5, we derive the following four estimates: as well as (17) and (12) yields where the constant Γ M is given by By assumption (14), we have the bound This bound in conjunction with estimate (18) yields Using Poincaré's inequality for g and integrating in time over the interval (t 1 , Thanks to Gronwall's inequality we obtain the following continuous-dependence estimate Finally, let us express this continuous-dependence estimate in terms of the original variables. To this end we recall that u − u * = h = g exp(−V ).
Substituting this relation into the preceding inequality, we deduce that

Application to an asymptotic gradient flow system
In this section we focus on an asymptotic gradient flow system which is used to describe the dynamics of multiple diffusing species in the presence of volume constraints. Such systems often have small cross-diffusion terms and therefore fall into the class of problems introduced in Section 2. We will discuss the equilibration behaviour of their solutions and establish sharper L ∞ -bounds (recall that, thanks to Proposition 2.1, we already know that solutions are essentially bounded if δ < δ 0 ).

4.1.
A cross diffusion model for two interacting diffusing species. We consider the asymptotic gradient flow system for two species that was introduced by Bruna and Chapman in [11]. It describes the behaviour of two different types of hard spheres and was derived from a stochastic system of hardcoreinteracting Brownian particles using the method of matched asymptotic expansions. Few analytic results are available, and the little results that exist seem to rely on the assumption that the asymptotic gradient flow system is close to its stationary state (see [10]) or that it is close to being a diagonal, decoupled linear problem; cf. [4].
The asymptotic gradient flow system presented in [11] can be written in the general form of (2). In particular: (20) A(x, u) = D 1 0 where δ i,j for i ∈ {1, 2}, j ∈ {1, 2, 3} are asymptotically small parameters with max i,j δ i,j =: δ < δ 0 . The δ i,j depend on the size of the particles and on the corresponding diffusion coefficients D 1 , D 2 ≥ 0. In [10] it was noted that system (2) with the nonlinearities in (20) and (21) has an asymptotic gradient flow structure.
In particular, it is well known that certain cross-diffusion systems possess a formal gradient-flow structure, that is, they can be formulated as where M ∈ R m×m is usually referred to as the mobility matrix and δE/δu is the variational derivative of the entropy functional, E. In order to highlight the connection between gradient flows and asymptotic gradient flows, let us consider the following entropy functional as well as the mobility matrix The cross-diffusion system (2) with diffusion and drift matrices (20) and (21), can be rewritten as where G is given by with δ 1,3 = (d − 1)δ 1,2 + dδ 2,2 and δ 2,3 = (d − 1)δ 2,2 + dδ 1,2 , a relation resulting from the derivation of the model, cf. [11]. The discrepancy between the system above and the gradient-flow induced by (23) is of order δ 2 , and one can show that they actually coincide if and only if both species have the same diffusivities D 1 = D 2 ≥ 0 as well as the same particle sizes, cf. [10]. In particular, having the same size and same diffusivities implies that δ i,1 = (d − 1)(δ 1,2 + δ 2,2 ) for i = 1, 2 as well as δ 1,j = δ 2,j for j = 1, 2, 3, cf. [10] and therefore G ≡ 0. Moreover, one can show that this is the only possibility for G to vanish. We refer to [10] for more details as well as the existence proof of unique stationary solutions in both cases. More specifically, we recall the following result: Proposition 4.1. Let δ = max{δ i,j }, assume 0 < δ < δ 0 and suppose that the potentials in (2) satisfy V i ∈ H 3 (Ω). Then systems (22) and (24) admit unique stationary states in H 3 (Ω), denoted by u ∞ and u * , respectively. Additionally, there exists a constant C > 0 such that Notice that, as a consequence, we deduce that the stationary states are in the space W 1,∞ . This fact will be useful in the next sections.

4.2.
The scalar problem. In the following we focus on a special case of system (24), namely where r = r(t, x) describes the density of the mobile species diffusing in the presence of a given immobile species b(x). Notice that we have simplified the notation by letting δ j = δ i,j and δ = max j δ j for i, j = 1, 2. We chose to use a different notation to emphasise that (26) is a scalar equation. It can be obtained from system (24) by setting r = u 1 , b = u 2 and setting the diffusion coefficient as well as the external potential of the species b to zero. Note that this equation can be derived as a macroscopic limit of a stochastic system with two types of particles, one species diffusing in a domain with fixed obstacles of a certain size interacting via hard-core collisions. Hence, the particles diffuse in a "perforated domain" with obstructions distributed according to the density b(x); cf. [9] for more details. (26) is a special case of (24), namely with immobile blue particles, Proposition 4.1 ensures that equation (26) admits a stationary state r * ∈ H 3 (Ω) provided that that δ is sufficiently small, that b, V ∈ H 3 (Ω) and that for some M 0 > 0 it holds that

Remark 6. As equation
From the discussion in Section 4.1 we know that equation (26) does not exhibit a full gradient flow structure, but only an asymptotic one. Hence, taking the entropy functional and mobility matrix M(r) = r(1 − δ 2 b), we observe that equation (26) can be cast in the form (24) where G(r, δ) = rb(δ 1 δ 2 ∇r + δ 2 δ 3 ∇b).

Remark 7.
As a consequence of Proposition 2.1, we have that Problem 26 admits bounded, regular solutions in the space W , however, in the rest of this section we will need less regularity for the solutions and we will work in a weaker "parabolic" space, Z, that we introduce in Definition 4.1.
Definition 4.1 (Parabolic space). Consider the space Z = C 0 (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω)). We use the so-called "parabolic norm" Similarly to Theorem 3.2, we obtain a result on exponential convergence to equilibrium. Note that in the scalar case, we can avoid the smallness assumption on the potential V (x) as well as the use of the L ∞ -bound on r.
Theorem 4.2 (Exponential convergence to equilibrium). Let d ∈ {1, 2, 3} and assume that b(x), V (x) ∈ H 3 (Ω) are such that V l ≤ V (x) ≤ V u , for some V l , V u ∈ R and for all x ∈ Ω. Moreover, if d = 1, 2, we suppose that r ∈ Z is a weak solution of problem (26) such that r Z ≤ L, for some constant L > 0; if d = 3 we also suppose that 0 ≤ r ≤ M for some constant M > 0 and a.e. (t, x) ∈ [0, T ] × Ω. Let the stationary state r * ∈ W 1,∞ (Ω). Then the following inequality holds for any t 1 , t 2 ≥ 0: for δ sufficiently small, satisfying relation (36).
Proof. In order to transform the terms in front of the potential V to order δ, we perform the change of variables w = r exp(V ) implying ∇r = (∇w − w∇V ) exp(−V ) and Now let us consider a perturbation, h, of the stationary state w * of (28) (which exists by Proposition 4.1), namely where Ω h(x, t) = 0 since the mass is conserved. We observe that h satisfies the following equation: We test equation (29) against h to obtain Integration in time over (t 1 , t 2 ) ⊂ [0, T ] yields Using the fact that ∇V, w * , ∇w * , ∇b ∈ L ∞ (Ω), we deduce that where where, for d = 2, we set K 2 = e −2V l ∇V L ∞ (Ω) C GN L and we used a special case of Gagliardo-Nirenberg interpolation inequality known as Ladyzhenskaya's inequality to bound h L 4 (Ω) , as well as h Z ≤ L; on the other hand, for d = 3, we set K 2 = 2e −2V l ∇V L ∞ (Ω) C S M and we used Sobolev's embedding in order to estimate h L 4 (Ω) . We define and we substitute equations (32) and (33) into equation (31), which, in turn, becomes Recalling that h + w * = w is non-negative, we choose δ sufficiently small so that the following relation is satisfied: which implies Thus we have obtained the following inequality: Using the Poincare's inequality and the fact that Thanks to Gronwall's inequality we obtain the following estimate: Finally, let us switch back to the original variables by recalling the fact that h = w − w * . We obtain Since V is bounded from above and below by assumption, a simple change of variable r = we −V implies the exponential convergence to equilibrium for the variable r as claimed in the statement, which concludes the proof.

4.3.
A maximum principle for the scalar problem (d = 2). In this section we will improve the L ∞bounds for solutions of problems of type (26), in dimension d = 2 (the same result can be obtained in one dimension). We adapt the strategy presented in [25], Chapter V 1 . In the context of elliptic equations, this approach is referred to as De Giorgi's method (see, e.g., [22]).

Theorem 4.3 (Maximum principle).
Let d = 2 and suppose that δ is sufficiently small and let r ∈ Z be a solution of problem (26). Then we have that for a.e. t ∈ [0, ∞), x ∈ Ω and for some σ ∈ ( 1 2 , 1) depending only on d ∈ {1, 2, 3}. Note that the constant C > 0 does not depend on time.
The proof of Theorem 4.3 uses several technical results which we shall recall in the following. Notice that, if x is close to a maximum point of V , estimate (39) is "almost sharp", in the sense that r is estimated by r 0 L ∞ (Ω) with a multiplicative constant that is very close to 1.
The following lemma, which is necessary to prove the following results, was originally stated in [25] Chapter II, and it can be proved proved directly by induction.
Lemma 4.4. Suppose that a sequence (a n ) n∈N ⊂ R + , of non-negative numbers satisfies the recursive relation a n+1 ≤ κ ζ n a 1+ε n and a 0 ≤ κ − 1 ε ζ − 1 ε 2 , for κ, ε > 0 and ζ > 1. Then we have that The following result is essential in order to prove Theorem 4.3.
Lemma 4.5. Let d = 2 and suppose that z ∈ Z is a weak solution of the problem We assume that δ ∈ (0, 1) is sufficiently small (see (44) and the assumptions below) and, furthermore, we suppose that Proof. The proof consists of several steps. We begin with an energy estimate involving Stampacchia's truncation and continue with an estimate for the measure of the superlevel sets S k,t := {x ∈ Ω : z(x, t) > k}, for k ∈ R + and t > 0. The strategy of the proof is the following: via careful use of a priori estimates, we will construct a sequence of the form for n ∈ N and suitable M > 0 as well as q > 0. Subsequently, we will apply Lemma 4.4 to a n in order to deduce that |S M,t | = 0 for a.e. t ∈ [0, T ]. This implies that z is essentially bounded.
Step 0: Stampacchia's truncation LetM ≥M (to be chosen later) and k ∈ R + , k > max(M , 1). It is easily verified that (z − k) + ∈ Z. Thus we may test equation (40) against the truncated function (z − k) + and obtain d dt Ω Rearranging this relation we get Step 1: upper bound for the RHS in terms of |{z ≥ k}| Let us focus on the estimate of the right-hand side: where we set C V = max{ F 1 L ∞ (Ω) , F 2 L ∞ (Ω) } and used the fact that k ≥ 1.
Notice that the integral on the right-hand side of equation (41) may be estimated further by Since z is a Lebesgue-integrable function, it cannot be unbounded on a region of positive measure, therefore, for a sufficiently large k (depending only on Ω z 0 dz thanks to Chebyshev's inequality), the function (z −k) + vanishes on a region of positive measure. This implies that a modified version of Poincaré's inequality given in Lemma 3.1 holds. Since d = 2, we derive a bound for (z − k) + 2 L 4 (Ω) applying Gagliardo-Nirenberg and Poincaré's inequalities (see Lemma 3.1); namely we get: This may be further estimated upon observing that Using this in conjunction with assumption (5) in the statement we obtain (z − k) + L 2 (Ω) ≤ z Z ≤ L.
Then the previous inequality leads to the following estimate: Notice that for the remaining term in estimate (42), we have Hence estimate (42) becomes Altogether we have obtained the inequality: Since ω ≥ µ and k ≥M , after an integration in time, we obtain Supposing that δ is such that where Γ = 7C V (4µ 2 min{µ/2, µλ − δC F }) −1 .
Step 2: lower bound for the LHS in terms of |{z ≥ k}| Notice that, as above, there exists a constant C e > 0 such that the following inequality holds for any 1 ≤ p < ∞ (it is a direct consequence of Sobolev's and Poincaré's inequality): We also know that, for any f ∈ Z, and thus in particular for f = (z − k) + , there holds f L ∞ (0,T ;L 2 (Ω)) < f Z . Thanks to the Riesz-Thorin interpolation theorem, we deduce that, for ν = p −1 ∈ (0, 1), where C e does not depend on time. By Chebyshev's inequality, for any a ≥ 1, we obtain whence, upon integration, we obtain Choosing a = 4 − ν, we get (47) Step 3: combining all previous estimates Combining (45), (46), and (47) we obtain T 0 |{z > h}|dt for any n ∈ N. In order to apply Lemma 4.4 we need to define a recursive relation for the sequence (a n ) n∈N . Let us set k = k n , h = k n+1 , then we observe that we observe that, from (48), Thus the sequence (a n ) n∈N satisfies the following relation Let us set ε = 1 − ν/2, ζ = 4, and κ = 16Γ δ and apply Lemma 4.4. We deduce that a n → 0 as n → ∞, and that the following inequality holds: provided that This will be satisfied choosingM sufficiently large. Thus, we shall now concentrate on giving an explicit lower bound forM . To this end, let θ > 1 be such thatM = θM . From (46) and (47) and upon choosing h = θM and k =M , we also know that, Notice that the quantity (z −M ) + Z is bounded uniformly in time by C 0 = L +M . Therefore we impose that θ satisfies which, in turn, gives such that the upper bound for a 0 is met. We conclude that T 0 |{z > θM }|dt = 0, after passing n → ∞, provided thatM and consequently we have obtained that, for a.e. t, x, With this result at hand the proof of Theorem 4.3 is immediate. Proof of Theorem 4.3. We apply Lemma 4.5 making the following choices: A(·, r) = r − b, L = Γ 0 (see Proposition 2.1), Recall that rewriting equation (26) in terms of the new unknown z = u exp(V ) we obtain with suitable initial and boundary conditions.

Numerical simulations
Let us recall the equation with zero flux conditions at the boundary. We note that the equation can be cast into a transport part and a remainder part, i.e., This formulation is the basis of the finite volume scheme used for the following section. The scheme is based on the schemes studied in [16,17,7]. For simplicity we introduce it here in one space dimension but an extension to two (or more) dimensions is immediate. We note that we could have used a different formulation to define the numerical fluxes, for example, based on (22) with entropy (27). However, the formulation chosen by us, is somewhat more natural, since it clearly highlights the transport part.
Next, we discretise the initial data by computing the cell averages of the continuous initial data on each cell, i.e., for all i ∈ I. Upon integrating the equation over [t n , t n+1 ) × C i , we obtain the following finite volume approximation for i ∈ I. Here the numerical fluxes are given by as well as where 0 < ε = 10 −7 1 is a small constant that is commonly chosen to regularise the log. The numerical no-flux boundary condition F 1/2 = F N +1/2 = 0. As usual, we use (z) ± to denote the positive (resp. negative) part of z, i.e., (z) + := max(z, 0), and (z) − := min(z, 0). 5.1. One Dimensional Explorations. Let us fix Ω = [−5, 5] as computational domain and set δ 1 = 10 −1 in this subsection. The immobile species is given by b(x) = 1/ √ 2σ 2 π exp(−x 2 /(2σ 2 )), with σ = 10 −1 . The active species is initially distributed according to r 0 (x) = 1/ √ 2σ 2 π exp(−(x + 3) 2 /(2σ 2 )), with σ = 1/2. The external potential is chosen as V (x) = (7/2−x) 2 χ {x<7/2} (x)+(x−7/2) 5 χ {x≥7/2} (x); cf. Figure 1. This choice forces the mobile species to penetrate the immobile species as it migrates towards the minimum of the external potential. In Figure 2 we present the evolution of the mobile species for three different choices of δ 2 = δ 3 . It is apparent that an increase in the parameters δ 2 = δ 3 leads to a more and more inhibited migration. This behaviour is perfectly physical. In fact, the parameters δ i are directly linked to the radii of particles of the mobile (resp. immobile) species. Since it is harder for large particles to traverse other large particles the motion is increasingly slowed down. For the second numerical exploration we chooseṼ (x) = x 2 /2 and assume that the mobile species is initially distributed according tor 0 (x) = 1 2 r 0 (x − 3) + 1 2 r 0 (x + 3), with r 0 as above; cf. Figure 1. Figure 3 shows the migration of the two initial bumps of the mobile species towards the minimum of the external potential. Since the immobile species inhabits this region we observe a competition for space leading to an indentation in the density of the mobile species, showcasing the strong effect of finite size exclusion effects. In Figure 4 we show the exponential convergence of the solution towards equilibrium.  . We provide two simulations based on the same set of parameters for two types of initial data (left), external potential (centre), and a given immobile species (right).

Two Dimensional Explorations.
In this subsection we set Ω = (−1.75, 1.75) 2 and δ 2 = δ 3 = 10 −1 . Note that this way both parameters are in the physical range in the sense of [11].
5.2.1. Immobile species as porous medium. In this example we assume the passive species is spread 'heterogeneously' according to the distribution b(x, y) = N −1 (1 + cos(5x) + cos(5y)), where N is chosen such that Ω b(x, y) dxdy = 1. Initially the active species is spread around the origin, i.e., r 0 (x, y) = 1 √ 2πσ 2 exp − r 2 2σ 2 , with σ 2 = 0.5. Moreover, we also choose δ 1 = 10 −1 and we assume the absence of any external forces. Thus the dynamics are dictated only by the internal dispersal of the mobile species and its interaction with the immobile species. As the time evolution continues the volume exclusion effects imposed on the active species Figure 2. The influence of finite size effects upon the evolution of the mobile species are displayed. At x = 0, the maximum of the immobile species, we observe the formation of a bump of the mobile species. As the radii of red and blue particles are increased (left to right) the mobile species peaks increasingly higher as it encounters the inactive species due to the inhibited traversal caused by the size exclusion effects.
(c) δ 2 = δ 3 = 3 × 10 −1 . Figure 3. Again, we can observe the impact of the size exclusion effect. The immobile species is located at x = 0. As we increase the size of the particles of both species (left to right) the mobile species attains its maximum at a lower level. This is in perfect agreement with the physial derivation of the model as the size exclusion affects the density in regions of coexistence, i.e., around x = 0.
by the passive species become visible. In the final state, complementary regions are occupied by the passive species and the active species, respectively; cf. Figure 5. After some time the evolution slows down as the mobile species reaches a stationary state. In Figure 6 we show the numerical long-time asymptotics on a semi-log scale. It can be seen that the dynamic relaxes to the numerically computed stationary state at an exponential rate.

5.2.2.
Immobile species as barrier. In this section we present another stunning example of how the immobile species can inhibit the evolution of the mobile species. We set δ 1 = 10 −2 and consider the initial datum r 0 (x, y) = 1 √ 2πσ 2 exp − (x + 3/4) 2 + (y + 3/4) 2 2σ 2 , Figure 4. The L 2 -distance (blue) between the evolution and the (numerically computed) minimiser u on a semi-log scale. A line of slope −1.4 (green) is superimposed to highlight the exponential convergence rate. The parameter choice corresponds to that of Figure 3 (a).
In Figure 8 it can be observed that the motion of the active species is slowed down and that its density is lower in the region occupied by the passive species compared to the unoccupied regions. Since the potential is relatively strong the active species moves through the immobile one only slightly changing its shape. In Figure 9 we rescaled the potential by a factor of ten. Again, we see that the motion of the mobile species is slowed down and that its density is reduced in the region occupied by the passive species in comparison with unoccupied regions. Since, now, the potential is relatively weak the motion of the active species is inhibited by the passive species which incentivises circumnavigating the immobile species instead of penetrating it. This behaviour is reflected in the numerical simulation; cf. Figure 9. Compared to Figure 8 the active species closes in from the side as it approaches the minimum of V , rather than directly moving into it.    Figure 9. The active species is incentivised to move around the passive species.