Cross-diffusion systems with excluded volume effects and asymptotic gradient flows

In this paper we discuss the analysis of a cross-diffusion PDE system for a mixture of hard spheres, which was derived by Bruna and Chapman from a stochastic system of interacting Brownian particles using the method of matched asymptotic expansions. The resulting cross-diffusion system is valid in the limit of small volume fraction of particles. While the system has a gradient flow structure in the symmetric case of all particles having the same size and diffusivity, this is not valid in general. We discuss local stability and global existence for the symmetric case using the gradient flow structure and entropy variable techniques. For the general case, we introduce the concept of an asymptotic gradient flow structure and show how it can be used to study the behavior close to equilibrium. Finally we illustrate the behavior of the model with various numerical simulations.


Introduction
Systems of interacting particles can be observed in biology (e.g. cell populations), physics or social sciences (e.g. animal swarms or large pedestrian crowds). Macroscopic models describing the individual interactions of these particles among themselves as well as their environment lead to complex systems of differential equations (cf. e.g. [5,6,7,8,9,10,12,22,23,24]). In microscopic models the dynamics of each particle is accounted for explicitly, while the macroscopic models typically consist of partial differential equations for the population density. Passing from the microscopic model to the macroscopic equations in a systematic way is, in general, very challenging, and often one relies on closure assumptions, which can be made rigorous under certain scaling assumptions on the number and size of particles. In particular, when crowding due to the finite-size of particles is included in the model, the limiting process is quite subtle and, using different assumptions and closure relations, a variety of macroscopic equations have been derived. For instance, the macroscopic equations of a two-species system where particles undergo a simple exclusion process on a lattice can be derived using formal Taylor expansions, see for example [8,24]. The case when particles are not confined to a regular lattice and undergo instead a Brownian motion with hard-core interactions was considered in [6] using matched asymptotic expansions. Cross-diffusion is a common feature of all these models and poses a particular challenge for the analysis since maximum principles do not hold. Classical examples of cross-diffusion systems are reaction diffusion systems or systems describing multicomponent gas mixtures. These quasi-linear parabolic systems were analyzed by Ladyzenskaya [19] or Amann [2,3], which however rely on strong parabolicity assumptions that break down for the degenerate cross-diffusion systems derived from the interacting particle systems mentioned above.
The canonical form for a two-species system of interacting particles (called red and blue in the following) is where D = D(r, b) is the diffusion matrix and F = F (r, b) is the drift matrix due to a convective flux.
Systems like (1) often have a gradient flow structure where M is a mobility matrix and ∂ r E and ∂ b E denote the functional derivative of an entropy functions E with respect to r and b, respectively. The gradient flow formulation provides a natural framework to study the analytic behavior of such systems, cf. e.g. [4]. It has been used to analyze existence and long-time behavior of systems, see for example [11,16,20,26]. As a result, being able to express a PDE system as gradient flows of an entropy is a very desirable feature; yet, this is not possible in general. The lack of the gradient flow structure on the PDE level can result from the approximations made when passing from the microscopic description to the macroscopic equations. This is the case of the cross-diffusion system derived in [6], which was derived using the method of matched asymptotics. There has been a lot of research on the passage from microscopic models to the continuum equations, for example in the hydrodynamic limit [18]. More recently the microscopic origin of entropy structures, which connects gradient flows and the large deviation principle was analyzed in [1,21].
In this paper we introduce the idea of an asymptotic gradient flow structure as a generalization of a standard or full gradient flow for systems derived as an asymptotic expansion such as that in [6]. In this paper we provide several analytic results for these cross-diffusion systems and introduce the notion of asymptotic gradient flows. We discuss how the closeness of these asymptotic gradient flow structures can be used to analyze the behavior of the system close to equilibrium. Furthermore we present a global in time existence result in the case of particles of same size and diffusivity (in which the system has a full gradient flow structure). The existence proof is based on an implicit Euler discretisation and Schauder's fixed point theorem. We study the linearized system with an additional regularization term in the entropy to ensure boundedness of the solutions and deduce existence results for the unregularised system in the limit (similar to the deep quench limit for the Cahn Hilliard equation [13]). This is, to the authors' knowledge, the first global in time existence result for this system so far. We note however that it is only valid if the total density stays strictly below the maximum density.
This paper is organized as follows: we introduce the mathematical model in Section 2 and discuss the cases for which the system has either a full or an asymptotic gradient flow structure. In Section 3 we define the notion of asymptotic gradient flows formally and discuss how they can be used to analyze the behavior of stationary solutions close to equilibrium. Several numerical examples illustrating the deviation of stationary solutions from the equilibrium solutions for asymptotic gradient flows are presented in Section 4. Finally, we give a global in time existence result in the case of particles of same size and diffusivity in Section 5.

The mathematical model
In this paper we analyze a cross-diffusion system for a mixture of hard spheres derived in [6], which we present below. The system is obtained as the continuum limit of a stochastic system with two types of interacting Brownian particles, referred to as red and blue particles. In particular, we consider N r red particles of diameter ǫ r , constant diffusion coefficient D r and external potential V r , and N b blue particles of diameter ǫ b , diffusion coefficient D b , and external potentialṼ b . Each particle evolves according to a stochastic differential equation (SDE) with independent Brownian dynamics, and interacts with the other particles in the system via hard-core collisions. This means that the centers of two particles with position X i and X j in space are not allowed to get closer than the sum of their radii, that is, X i − X j ≥ (ǫ i + ǫ j )/2, where ǫ i denotes the radius of the ith particle. We define the total number of particles in the system by N = N r + N b , and the distance at contact between a red and blue particles by ǫ br = (ǫ r + ǫ b )/2. The situation detailed above can be described by the overdamped Langevin SDEs where X i ∈ Ω ⊂ R d , d = 2, 3, is the position of the ith particle and W i a d-dimensional standard Brownian motion. We assume that Ω is a bounded domain. The boundary conditions due to collisions between particles and with the domain walls are where n denotes the outward unit normal. The continuum-level model associated to this individualbased model was derived in [6] using the method of matched-asymptotic expansions in the limit of low but finite volume fraction. If v d (ǫ) is the volume of a d-dimensional ball of diameter ǫ, then the volume fraction in the system is assuming that the problem is nondimensionalised such that the domain Ω has unit volume, |Ω| = 1.
Because particles cannot overlap each other, in addition to the global constraint Φ ≪ 1 there is also a local constraint on the total volume density, defined as In particular, the local volume density cannot exceed the theoretical maximum allowed volume fraction, given by the Kepler conjecture. We note that Φ and φ are related via Φ = Ω φ dx.
The cross-diffusion model in [6] is valid for any number of blue and red particles, N b and N r . However, here we will consider the case that the number of both particles is large, such that N r −1 ≈ N r ,N b − 1 ≈ N b , as it simplifies the model slightly. In this case, the model reads [6] where r = r(x, t) and b = b(x, t) are the number densities of the red and blue species, respectively, depending on space and time. Consequently, meaningful solutions satisfy r ≥ 0 with Ω r dx = N r and b ≥ 0 with Ω b dx = N b . In (7), V i =Ṽ i /D i are the rescaled potentials, and the parameters α, β i and γ i depend on the geometry of the particles. For balls, they are given by for i = r or b, and space dimension d = 2 or 3. This system is an asymptotic expansion in ǫ r , ǫ b (assuming that both small parameters are of the same asymptotic order, ǫ r ∼ ǫ b ∼ ǫ), valid up to order ǫ d . The nonlinear terms in (7) correspond to the leading-order contribution of the pairwise particle interactions. The asymptotic method used in [6] could be extended if desired to evaluate higher-order terms coming from three or more particle interactions, as well as higher-order corrections in the pairwise interaction. This would result in higher-order terms in ǫ i in (7) (of order ǫ (d+1) i and higher) with quite some effort. However, it seems impossible to derive the full infinite series expansion.
We will consider the system (7) in Ω × (0, T ) with no-flux boundary conditions on ∂Ω × (0, T ) and initial values (10) r In order to analyze the cross-diffusion system (7), it is convenient to consider its associated gradient flow structure of the form (2). However, only the system in the symmetric case where red and blue particles have same size and diffusivity can be rewritten in that form. For the general case we introduce a generalization of a gradient flow, namely an asymptotic gradient flow, motivated by the underlying structure of the general system (7).
2.1. Cross-diffusion system for particles of the same size and diffusivity. In this section we suppose that red and blue particles are of the same size, that is ǫ r = ǫ b := ǫ, and have the same diffusion coefficient, D r = D b . Without loss of generality, we take the diffusion coefficient equal one (this can be achieved by rescaling time). In this case, the cross-diffusion system (7) can be written as where β i and γ i , for i = r, b, are now equal and simplify to γ = π/d and β = 2 d−1 γ, respectively.
This cross-diffusion system can be used to describe a mixture of particles that are physically identically but that are driven by different potentials V r and V b (for example cells that are attracted to different food sources, or pedestrians that want to move in different directions). Moreover, it can also be used to model the scenario where the red and the blue particles are in fact identical, but one has knowledge about the initial distributions of each sub-population, r 0 and b 0 . This is the scenario in many experimental set-ups that use noninvasive fluorescent tagging. On the other hand, if the red and blue particles are identical and initially indistinguishable, then one has that r/N r = b/N b := p for all times. In this case, both equations (11a) and (11b) reduce to the same equation, which coincides with the equation for the evolution of a single population of hard spheres as expected [7].
In the following we defineᾱ = ǫ d α,γ = ǫ d γ, and the total number density When particles have the same size and diffusivity we find that where φ is the total volume density given in (6). Using ρ, the equations (11) can be rewritten in the following form where we have used that β = α + γ. It is straight-forward to see that the system (14) has a formal gradient flow structure, with an entropy functional given by Using the corresponding entropy variables the system can be written in the form with the symmetric mobility matrix 2.2. Cross-diffusion system for particles of different size and diffusivity. In this section we attempt to write a gradient flow structure for the general cross-diffusion system (7) guided by the symmetric case in the previous subsection, (15) and (18). We will see that this requires a generalization of the definition of gradient flow structure. We define the following entropy and mobility matrix As mentioned earlier, we suppose that the red and blue particle sizes are of the same asymptotic order, namely ǫ r ∼ ǫ b . It is then convenient to introduce a single small parameter ǫ and the order one parameters a r , a b and a br such that ǫ d i = a i ǫ d . Then the entropy and mobility can be expressed as Using (19), the general cross-diffusion system (7) can be rewritten as where G(r, b) is the vector Then it is easy to see that the gradient flow structure induced by (19) and our system (7) (or (21)) agree up to order ǫ d , which is the order of the asymptotic expansion that produced (7) in the first place. In other words, the discrepancy between the system (7) and the gradient-flow induced by (19) is of order ǫ 2d . Therefore, up to order ǫ d , we can see (21) as a gradient flow structure of our system. We will call this an asymptotic gradient flow structure; the precise definition will be made clear in the following section. Finally, we note that the system (19) coincides with the gradient-flow structure in the case D r = D b and ǫ r = ǫ b , see (15) and (18). Note that G ≡ 0 for the parameter values of the simpler system (11), as expected. Specifically, we find that if D r = D b and ǫ r = ǫ b , then θ r = θ b = 0. A natural question to ask is whether there are other parameter values for which G(r, b) ≡ 0 for all r, b. Imposing that θ r = θ b = 0 leads to the condition a 2 br = a r ab, which in turn leads to ǫ r = ǫ b , and thus that D r = D b . Therefore, the only case for which (19) is an exact gradient flow for the system is the case which we have already studied, that is when the particle sizes and diffusivities are equal.

Gradient Flows and Asymptotic Gradient Flows close to Equilibrium
In the following we provide a more detailed discussion on gradient flow structures and implications for the behavior close to equilibrium.
3.1. Full gradient flow structure case. In this subsection we analyze the behavior of system (17) close to equilibrium. We follow the strategy outlined in the previous subsection, by proving uniqueness of equilibrium solutions and studying the stability and well-posedness of the system close to this equilibrium solution.
We have seen in the previous subsection that the linear stability analysis for gradient flow structures reduces to showing that the mobility matrix M is positive definite in the case of a strictly convex entropy functional E, cf. [23]. We assume from now on: We recall that in case of assumption (AI) an equilibrium solution (r ∞ , b ∞ ) exists and that the corresponding entropy variables u ∞ and v ∞ are constant. The determinant of the mobility matrix M defined in (18) is given by (24) det M = rb(1 −γρ).
Together with the positivity of diagonal entries we see that M is positive definite if ρ < 1/γ. This constraint gives a local bound on the total local volume density (using (13)), namely 2φ < 1. This is consistent with the asymptotic assumption that φ ≪ 1. Hence, we define the set which will also use in the existence proof presented in Section 5. For stability and uniqueness it will be crucial to have solutions staying strictly in the interior of S, due to the degeneracy of the mobility matrix on the boundary of S.
Proof. Due to the gradient flow structure, any stationary solution of (17) is a minimizer of the entropy subject to the constraints of given mass and (r(x), b(x)) ∈ S almost everywhere. Due to the strict convexity of the entropy and the convexity of the constraint set, the minimizer is unique. Let us consider the linearisation of system (17) around the unique equilibrium, which corresponds to the constant entropy variables (u ∞ , v ∞ ). As we have seen before this is equivalent to have a linear expansion in (r, b) and in the entropy variables (u, v) , i.e. u = u ∞ + ξ, v = v ∞ + η. In the latter setting we obtain the following first order approximation where E ⋆′′ denotes the Hessian of the dual entropy functional. Note that for the first order approximation, we also have no flux boundary conditions. A simple calculation shows that E ⋆′′ (u ∞ , v ∞ ) as well as M (r, b) are positive definite for (r, b) in the interior of S, which is guaranteed everywhere for the stationary solution (r ∞ , b ∞ ). Stability of this linear system is equivalent to nonpositivity of all the real parts of eigenvalues λ in Note that due to the symmetry of the eigenvalue problem, all eigenvalues are real. Moreover, we find Since E ⋆′′ and M (r, b) are positive definite, we conclude that λ < 0, which implies linear stability.
Theorem 3.2 (Well-posedness). Consider system (17) with initial data u 0 , v 0 ∈ H 2 (Ω) and poten- for κ > 0 sufficiently small. Then, there exists a unique solution to system (17) in and R is a constant depending on κ and T > 0 only.
Proof. The proof is based on Banach's fixed point theorem. The corresponding fixed point operator is defined by the evolution of u − u ∞ and v − v ∞ , which can be written as where we used that (r, b) = E ⋆′ (u, v). Note that by using a similar argumentation as in the proof of Lemma 3.1, we can show that the stationary solutions r ∞ , b ∞ are in H 3 (Ω) assuming that the potentials V r , V b are in H 3 (Ω). Consider (u, v) ∈ X × X with the corresponding function r = r(u, v), b = b(u, v) and let L denote the solution of (30) for a given right-hand side. Then the fixed point operator is given by the concatenation of L and F , that is Standard results for linear parabolic equations, see [19] or [14], ensure that the solution (ũ − u ∞ ,ṽ − v ∞ ) to equation (30) lie in X × X.
To apply Banach's fixed point theorem, it remains to show that the operator J is self-mapping into the ball B R and contractive. The selfmapping property follows from the fact that For the contractivity we consider (u 1 , v 1 ) ∈ X × X and (u 2 , v 2 ) ∈ X × X and deduce that: for some constant C 1 > 0. Hence, we have that for some C > 0. Choosing κ and R such that R < 1 C , we can apply Banach's fixed point theorem which guarantees the existence of unique solutions (u, v) ∈ B R .

3.2.
Asymptotic Gradient Flow Structure. We have seen in Section 2.1 that system (7) with particles of same size satisfies a gradient flow structure, which is not valid for the general system due to terms of higher order in ǫ. However, we want to interpret the latter as an asymptotic gradient flow structure, motivated by the fact that it was derived from an asymptotic expansion in ǫ. For further motivation, consider a gradient flow structure for the density w of the form where both the mobility M and the entropy E depend on a small parameter δ > 0. With an expansion of M and E in terms of δ as Truncating the expansion on the right-hand side at a finite k does not yield a gradient flow structure in general, but up to terms of order δ k it coincides with the gradient flow structure with mobility k j=0 δ j M j (w) and entropy k j=0 δ j E j (w). In our case we deal with the example k = 1 (with δ = ǫ d ), where we have . Adding a term of order δ 2 , namely δ 2 ∇ · (M 1 (w)∇E ′ 1 (w)), this equation becomes a gradient flow. This motivates a more general definition: F (.; δ) be a densely defined operator on some Hilbert space for δ ∈ (0, δ * ). Then the dynamical system is called an asymptotic gradient flow structure of order k if there exist densely defined operators G j , j = k + 1, . . . , 2k such that for δ ∈ (0, δ * ) for some (parametric) energy functional E(·; δ), and M(w; δ) is a densely defined formally positivedefinite operator for each w.
If an expansion of mobility and entropy up to order k are available, it seems natural to perform a separate expansion to derive a lower order model that is a gradient flow as well. For complicated models and types of expansions as in [6] or [7] it seems not suitable to derive such however. Hence, we shall work with the asymptotic gradient flow concept below. Note that with the above notations we can rewrite (32) as which opens the door to perturbation arguments in the analysis of (32) for δ sufficiently small. In the remainder of this section we will highlight in particular the use of asymptotic gradient flow structures close to equilibrium. Let w δ ∞ denote the equilibrium solution, which is a minimizer of the energy functional on the manifold defined by M. Hence w δ ∞ solves M(w; δ)E ′ (w δ ∞ ; δ) = 0 for any w. In the case of (31) it typically means that E ′ (w δ ∞ ; δ) is constant. In order to prove the existence of a stationary solution of (32) one can then try the following strategy: first of all compute w δ ∞ (or prove at least its existence and uniqueness by variational principles) and then use the equation as the basis of a fixed-point argument, freezing w on the right-hand side. Since the terms on the right-hand side are of high order in δ or of second order in terms of w − w δ ∞ , there is some hope of contractivity of the fixed-point operator close to equilibrium w δ ∞ . Such an approach can also yield some structural insight into the stationary solution, since it will be a higher order perturbation of w δ ∞ . The same idea can be employed to analyze transient solutions of (32), since If M(w δ ∞ ; δ) is invertible and E(·; δ) is strictly convex on its domain, one can directly apply variational techniques to analyze the fixed point operator. In particular it can be rather beneficial to set up the fixed-point operator in dual (or entropy) variables z = E ′ (w; δ) instead.
Finally let us comment on the linear stability analysis around a stationary solution w δ * . Using a similar way of expanding the equation around w δ ∞ , the linearised problem for a variablew is given by , where we denote by E ′ and M ′ the derivatives with respect to w at fixed δ. Due to positive definiteness of E ′′ (w δ * ; δ), this system can be interpreted as a linear equation for the linearised entropy variablez = E ′′ (w δ * ; δ)w, which is equivalent to considering linear stability directly in the transformed equation for the entropy variable z as performed in [23]. Using the simplified notation A = E ′′ (w δ * ; δ) −1 and B = M(w δ ∞ ; δ), we obtain ). In the case of a gradient flow (G j ≡ 0, w δ * = w δ ∞ ) this reduces to A∂ tz + Bz = 0, which is stable if A and B are positive definite. In the asymptotic gradient flow case, with w δ * = w δ ∞ + O(δ k+1 ), we can formally write the linearised problem as (35) A∂ tz + (B + δ k+1 C)z = 0, and hence expect linear stability also for w δ * if δ is sufficiently small. The application of the above strategies to prove existence of solutions and linear stability to a concrete model obviously depends on an appropriate choice of topologies. In the remaining part of this section we focus on the analysis of the asymptotic gradient flow of the general model.

3.3.
Asymptotic gradient flow structure case. First we study the existence of stationary solutions to (21). Then we discuss stability of stationary states following the ideas presented in subsection 3.2.
Note that for ǫ = 0, the equilibrium solutions are given by (r ∞ , b ∞ ) = (C r e −Vr , C b e −V b ), with constants C r and C b depending on the initial masses only. Hence (r ∞ , b ∞ ) are bounded for V r and V b satisfying assumption (AI). For ǫ > 0, the equilibrium solutions are a O(ǫ d ) perturbation in L ∞ and therefore also uniformly bounded.
Theorem 3.4 (Existence of stationary solutions). Consider system (21) with potentials V r , V b ∈ H 3 (Ω). Then there exists a unique stationary state (u * , v * ) to system (21) in where X = H 3 (Ω) and R depending on ǫ and T > 0 only. Proof. We follow the ideas detailed in Subsection 3.2 and define a fixed point operator close to equilibrium. Denote by (r ∞ , b ∞ ) the minimizer of the entropy functional E ǫ (r, b), which exists as the entropy functional is strictly convex. Then any stationary solution to system (21) exists has to satisfy Similar arguments as in Lemma 3.1 ensure that for (u, v) ∈ X × X the functions r = r(u, v) and b = b(u, u) lie in X × X. Let L denote the solution operator to (36) for a given right-hand side F (u, v). Then the fixed point operator is constructed by: Hence, we can conclude that F maps from Employing results about the elliptic operator, cf. [15] or [14], we obtain that the solution (ũ,ṽ) to equation (36) is in X × X.
To apply Banach's fixed point theorem, it remains to show that the operator J is self-mapping into the ball B R and contractive. The self-mapping property follows from the fact that For the contractivity we consider (u 1 , v 1 ) ∈ X × X and (u 2 , v 2 ) ∈ X × X. Then for some constants C 1 , C 2 , C 3 > 0 and therefore for some C > 0. Choosing R and ǫ such that we can apply Banach's fixed point theorem which guarantees the existence of unique solutions (u * , v * ) ∈ B R .
A direct consequence of the proof is the closeness of the stationary solution (u * , v * ) to the gradient flow solution (u ∞ , v ∞ ): Corollary 3.1. Let the assumptions of Theorem 3.4 be satisfied. Then there exists a constant C > 0 such that for ǫ sufficiently small Proof. We use (36) rewritten as and the properties of the operators used above immediately imply the assertion.
We conclude this section by discussing linear stability of system (21) close to its stationary states (u * , v * ). Following the ideas presented in Section 3.2 we rewrite (21) as The linearisation of equation (38) around (r * , b * ) is given by the following system for (r,b): Using the linearised entropy variables (ũ,ṽ) = E ′′ (r * , b * )(r,b) we obtain where A = E ′′ −1 (r * , b * ) is a positive and B = M(r ∞ , b ∞ ) are negative semidefinite operator. Note that with the usual settings for elliptic systems, B is elliptic and hence invertible on the space of function pairs in H 1 (Ω) with zero means.
As already mentioned in Section 3.2, (r * , b * ) = (r ∞ , b ∞ ) + O(ǫ 2d ) and (39) can be written as for some bounded operator C on H 1 (Ω) 2 . As B is symmetric and negative definite except on the two-dimensional space of constant functions also annihilated by C, the nonzero eigenvalues of B+ǫ 2d C stay negative for ǫ sufficiently small, yielding linear stability for (r * , b * ), cf. [17].

Numerical investigations of steady states
In this section we compute the stationary solutions of (7). For the symmetric system (11), the solutions can be computed exactly as the minimizers of the entropy E in (15). If the mobility matrix (18) is positive definite (which it is under the assumptions), the equilibrium states can be computed by finding constants χ r ∈ R and χ b ∈ R such that ∂ r E = χ r and ∂ b E = χ b subject to normalization constraints. In the case of system (11) we have System (40) defines a nonlinear operator equation F (r ∞ , b ∞ , χ r , χ b ) = 0, which can be solved via Newton's method. Note that the no-flux boundary conditions are automatically satisfied by assuming that ∂ r E and ∂ b E are constant.
For the general case (7) we only obtain an asymptotic gradient flow structure with the entropy E ǫ ; if we use (40) to solve for the stationary solutions we will be committing an order ǫ 2d error. Instead, we compute the exact stationary states (r * , b * ) of the general system by solving the time-dependent problem (7) for long-times, until the system has equilibrated. To solve (7), we use a second-order accurate finite-difference scheme in space and the method of lines with the inbuilt Matlab ode solver ode15s in time.
We set d = 2 and consider one-dimensional external potentialsṼ r =Ṽ r (x) andṼ b =Ṽ b (x) so that the stationary states will be also one-dimensional. In particular, we take linear potentialsṼ r = v r x andṼ b = v b x and solve for the full system (7) and for the minimizers (40) in [−1/2, 1/2], which is split into 200 intervals. The Newton solver is initialized with the stationary state solution in the case of point particles and terminated if F (r, b, χ r , χ b ) L 2 (0,1) ≤ 10 −8 . Example 1. First we consider the case: ǫ r = ǫ b and D r = D b , that is particles of the same size and diffusivity. In this case, system (11) has a full gradient flow structure and hence we expect that the stationary states computed with the two approaches to be the same. We plot the two pairs, (r * , b * ) computed as the long-time limit of (11), and (r ∞ , b ∞ ), computed from (40) in Figure 1. The parameters are D r = D b = 1, ǫ r = ǫ b = 0.01, N b = N r = 200 and v r = 2, v b = 1. As expected, the solutions are identical. Example 2. From Corollary 3.1 we expect the stationary solutions corresponding to the case of an asymptotic and a full gradient flow equation agree up to order O(ǫ d ). To investigate this, we again compare the solutions (r * , b * ) and (r ∞ , b ∞ ) as we move away from the case with an exact gradient-flow structure (which corresponds to θ r = θ b = 0, see (21) and (23)).
In particular, we do a one-parameter sweep with θ r , increasing it from 0 (as in Figure 1) to 9 · 10 −5 , while keeping ǫ r = ǫ b = 0.01 and D b = 1 fixed. This ensures that when θ r = 0 then θ b = 0. The reds diffusivity D r is varied according to (23). We plot the result for θ r = 8 · 10 −5 in Figure 2. As expected, the error between the stationary solutions is apparent.
The absolute error and the relative error between the solutions, r ∞ − r * and b ∞ − b * and r ∞ − r * / r ∞ and b ∞ − b * / b ∞ , respectively, as a function of θ r is shown in Figure 3.
To conclude this section, we compute the stationary solutions of the (exact) full system and that approximated by the asymptotic gradient flow system as we vary ǫ, where ǫ = ǫ b = ǫ r , while keeping all the other parameters fixed. We plot the results in Figure 4. As expected from Corollary 3.1, the errors scale with ǫ 2d = ǫ 4 .

Global existence for the full gradient flow system
In this section we present a global in time existence result for the system with particles of same size and diffusivity (14).
Moreover, the solution satisfies the following entropy dissipation inequality: where and C ≥ 0 is a constant. PSfrag replacements Figure 2. Stationary solutions (r * , b * ) and (r ∞ , b ∞ ) from solving the long-time limit of (7) and (40), respectively, in a case with θ r = 8 · 10 −5 . The parameter values are d = 2, D r = 0.2,  We recall that system (14) can be written as a gradient flow: where Note that if r, b and ρ ∈ S • , then the matrix M is positive definite. We perform a time discretisation of system (43) using the implicit Euler scheme. The resulting recursive sequence of elliptic problems is then regularized. Let N ∈ N and let τ = T /N be the time step size. We split the time interval into the subintervals Then for given functions (r k−1 , b k−1 ) ∈ S, which approximate (r, b) at time τ (k − 1), we want to find (r k , b k ) ∈ S solving the regularized time discrete problem where we use the modified entropỹ with associated entropy variables The additional term in the entropy provides upper bounds on the solutions and the higher order regularization terms guarantee coercivity of the elliptic system in H 1 (Ω), which is needed to show existence of weak solutions to a linearized version of the problem (44) using Lax-Milgram. The existence result of the corresponding nonlinear problem is concluded by applying Schauder fixed point theorem.
Finally uniform a priori estimates in τ and the use of a generalized version of the Aubin-Lions lemma allow to pass to the limit τ → 0 leading to the existence of (41). Note that the compactness results are sufficient for 1 −γρ > 0 to pass to the correct limit in the flux terms J r and J b , i.e leading to the global existence of weak solutions to system (14).
is strictly convex and belongs to C 2 (S • ). Its gradienth ′ : S • → R 2 is invertible and the inverse of the Hessianh ′′ : S • → R 2×2 is uniformly bounded.
Proof. Note thath The matrixh ′′ is positive definite on the set S • , soh is strictly convex. We can easily deduce that the inverse ofh ′′ exists and is bounded on S • .
To apply Schauer's fixed point theorem, we need to show that the map S: (i) maps a convex, closed set onto itself, (ii) is compact, (iii) is continuous. Since S is convex and closed, property (i) is satisfied; (ii) follows from the compact embedding H 1 (Ω, R 2 ) ֒→ L 2 (Ω, R 2 ). Continuity (iii): let (r k ,b k ) be a sequence in S converging strongly to (r,b) in L 2 (Ω, R 2 ) and let (ũ k ,ṽ k ) be the corresponding unique solution to (48) in H 1 (Ω; R 2 ). As the matrix M only contains sums and products of r and b, we have that M (r k ,b k ) → M (r,b) strongly in L 2 (Ω, R 2 ). The positive semidefiniteness of the matrix M for (r, b) ∈ S provides a uniform bound for (ũ k ,ṽ k ) in H 1 (Ω; R 2 ). Hence, there exists a subsequence with (ũ k ,ṽ k ) ⇀ (ũ,ṽ) weakly in H 1 (Ω; R 2 ). The L ∞ bounds of M (r k ,b k ) and the application of a density argument allow us to pass from test functions (Φ 1 , Φ 2 ) ∈ W 1,∞ (Ω, R 2 ) to test functions (Φ 1 , Φ 2 ) ∈ H 1 (Ω, R 2 ). So, the limit (ũ,ṽ) as the solution of problem (48) with coefficients (r,b) is well defined. Due to the compact embedding H 1 (Ω, R 2 ) ֒→ L 2 (Ω, R 2 ), we have a strongly converging subsequence of (ũ k ,ṽ k ) in L 2 (Ω, R 2 ). Since the limit is unique, the whole sequence converges. From Lemma 5.1 we know that (r, b) = h ′−1 (ũ,ṽ) is Lipschitz continuous, which yields continuity of F .
Hence, we can apply Schauder's fixed point theorem, which assures the existence of a solution (r, b) ∈ S to (48) with (r,b) replaced by (r, b).
Lemma 5.4. The discrete time derivatives of r τ and b τ are uniformly bounded, i.e.
Proof. Let Φ ∈ L 2 (0, T ; H 1 (Ω)). Using the a priori estimates from Lemma 5.3 gives A similar estimate can be deduced for b which concludes the proof.
Even though the a priori estimates from Lemma 5.3 are enough to get boundedness for all terms in (55) in L 2 (Ω T ), the compactness results are not enough to identify the correct limits for τ → 0. From Lemma 5.3 we get that, as τ → 0 τũ τ , τṽ τ → 0 strongly in L 2 (0, T ; H 1 (Ω)).
In order to identify the limit terms, we multiply equation (63) by (1 − γρ).
Note that the L ∞ bounds for b τ and r τ imply that, up to a subsequence, (67) r τ ⇀ r, b τ ⇀ b weakly * in L ∞ (Ω T ).
The convergence of (ii) follows from the L ∞ bounds, the a priori estimate (58) as well as from the convergences (65) and (70).
In this paper we studied a mean-field model for two species of interacting particles which was derived using the method of matched asymptotics in the case of low volume fraction. This asymptotic expansions results in a cross-diffusion system which has a gradient flow structure up to a certain order. We therefore introduce the notion of asymptotic gradient flows for systems whose gradient flow structure is perturbed by higher order terms. We show that this 'closeness' to a classic gradient flow structure allows us to deduce existence and stability results for the perturbed or as we call them asymptotic gradient flow system.
While the presented results on linear stability (Theorem 3.1), well-posedness (Theorem 3.2) and existence of stationary solutions (Theorem 3.4) also hold on unbounded domains, the proof of the global existence result in Section 5 uses embeddings which do not hold on unbounded domains in general, e.g. H 2 (Ω) is compactly embedded in L 2 (Ω).
The presented work is a first step towards the development of a more general framework for asymptotic gradient flows. It provides the necessary tools to understand the impact of high order perturbations on the energy dissipation as well as the behavior of solutions and opens interesting directions for future research.