Cross-Diffusion Systems with Excluded-Volume Effects and Asymptotic Gradient Flow Structures

In this paper, we discuss the analysis of a cross-diffusion PDE system for a mixture of hard spheres, which was derived in Bruna and Chapman (J Chem Phys 137:204116-1–204116-16, 2012a) 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., Bendahmane et al. 2009;Bruna and Chapman 2012a, b;Burger et al. 2010Burger et al. , 2016Burger et al. , 2012Di Francesco and Fagioli 2016;Painter 2009;Schlake 2011;Simpson et al. 2009). 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 Burger et al. (2010), Simpson et al. (2009). The case when particles are not confined to a regular lattice and undergo instead a Brownian motion with hardcore interactions was considered in Bruna and Chapman (2012a) 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 Ladyzhenskaia et al. (1968) or Amann (1985Amann ( , 1989, 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 function 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., Ambrosio et al. (2008). It has been used to analyze existence and long-time behavior of systems, see for example Carrillo et al. (2014), Jüngel and Zamponi (2014), Liero and Mielke (2013), Zinsl and Matthes (2015). 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 Bruna and Chapman (2012a), 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 (Kipnis and Landim 2013). More recently, the microscopic origin of entropy structures, which connects gradient flows and the large deviation principle, was analyzed in Adams et al. (2011), Liero et al. (2015).
In this paper, we introduce the idea of an asymptotic gradient flow structure as a generalization of a standard or, as we also call it, full gradient flow structure for systems derived as an asymptotic expansion such as that in Bruna and Chapman (2012a). 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 discretization 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 unregularized system in the limit (similar to the deep quench limit for the Cahn Hilliard equation Elliott and Garcke 1996). 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 a certain threshold. It relates to the fact that the model assumptions break down if the maximum density is reached. We discuss this problem in more detail in Sect.

3.
This paper is organized as follows: we introduce the mathematical model in Sect. 2 and discuss the cases for which the system has either a full or an asymptotic gradient flow structure. In Sect. 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 Sect. 4. Finally, we give a global in time existence result in the case of particles of same size and diffusivity in Sect. 5.

The Mathematical Model
In this paper, we analyze a cross-diffusion system for a mixture of hard spheres derived in Bruna and Chapman (2012a), 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Ṽ 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 hardcore 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 ddimensional 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 individual-based model was derived in Bruna and Chapman (2012a) 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 nondimensionalized 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 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. This means that r = N rr and b = N bb for the probability densitiesr ,b. Consequently, meaningful solutions satisfy r ≥ 0 with r dx = N r and b ≥ 0 with b dx = N b . For more detailed information, see Bruna and Chapman (2012a). 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 Bruna and Chapman (2012a) 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 (Bruna and Chapman 2012a) where 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, j = b and vice versa, 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 Bruna and Chapman (2012a) 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 higherorder 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 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).

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 identical 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 setups 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 Eqs. (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 (Bruna and Chapman 2012b).
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 Eq. (11) can be rewritten in the following form where we have used that β = α + γ . It is straightforward 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

Cross-Diffusion System for Particles of Different Size and Diffusivity
In this section, we attempt to write a gradient flow structure for the general crossdiffusion 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 . Note that as br = ( r + b )/2, it holds that a br = ((a 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 with 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 a b , 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.

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. Schlake (2011). 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 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 we will also use in the existence proof presented in Sect. 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.
Theorem 3.1 (Linear stability) The stationary solutions of the system (17) are unique and linearly stable with respect to small perturbations ξ, η ∈ L 2 (0, T ; H 1 ( )) with zero mean.
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 linearization 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 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 firstorder 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.
Note that we assumed ξ, η with zero mean, which corresponds to the mass conservation property of the system. Next, we consider the well posedness close to equilibrium. We shall make use of the following auxiliary lemma: Then, the gradient of the dual entropy functional E * : (16), is continuous.
for κ > 0 sufficiently small. Then, there exists a unique solution to system (17) in where 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 constructed by considering 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 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 Ladyzhenskaia et al. (1968) or Evans (1998) 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 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 .

Asymptotic Gradient Flow Structure
We have seen in Sect. 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: Definition 3.3 Let 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 positive definite 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 Bruna and Chapman (2012a) or Bruna and Chapman (2012b), 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 linearized 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 linearized 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 Schlake (2011). 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 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 linearized problem as 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.

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 Sect. 3.2.
Note that for = 0, the equilibrium solutions are given by (r ∞ , b ∞ ) = (C r e −V r , 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 3.1. 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
where X = H 3 ( ) and R depending on and T > 0 only.
Proof We follow the ideas detailed in Sect. 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. Gilbarg and Trudinger (2015) or Evans (1998), we obtain that the solution (ũ,ṽ) to Eq. (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 Sect. 3.2, we rewrite (21) as The linearization of Eq. (38) around (r * , b * ) is given by the following system for (r ,b):

Using the linearized entropy variables
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 Sect. 3.2, (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. Kato (2013).

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 V 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)  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)-(23)]. We recall that the parameters θ r and θ b are only zero in the symmetric case with identical particles, so that can think of these as a measure of the asymmetry in the system (either in size, diffusivity, or both). In particular, we do  (7) and (40), respectively, in the case with θ r = θ b = 0. The parameter values are d = 2, 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, a one-parameter sweep with θ r , increasing it from 0 (as in Fig. 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 Fig. 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 Fig. 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 Fig. 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). We look for a weak solution (r, b) : × (0, T ) → S to the system Note the special fluxes J r and J b , where the multiplication with (1 − γ ρ) is due to the fact that J r , J b are not well defined for maximal packing. From a variational perspective, one might expect additional Lagrange parameters to change J r , J b if the constraint γρ ≤ 1 is active.
Theorem 5.1 (Global existence in the case of small volume fraction) Let T > 0, let (25), be a measurable function such that E(r 0 , b 0 ) < ∞. Then, there exists a weak solution (r, b) : × (0, T ) → S to system (41) satisfying Moreover, the solution satisfies the following entropy dissipation inequality: where 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 discretization 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 1 τ 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).

Lemma 5.1 The entropy densitỹ
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.
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).
Inequality (54) becomes which provides us the following a priori estimates. Note that from now on K denotes a generic constant.
Lemma 5.3 (A priori estimates) There exists a constant K ∈ R + , such that the following bounds hold: where T = × (0, T ).

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 ( )).
Together with Lemma 5.4, we get a solution to where weakly in L 2 ( T ). In order to identify the limit terms, we multiply Eq. (63) by (1 − γ ρ).
Proof The estimates from Lemma 5.3 and Lemma 5.4 allow us to use Aubin's lemma to deduce the existence of a subsequence (not relabeled) such that, as τ → 0: ρ τ → ρ strongly in L 2 ( T ).
The convergence of (ii) follows from the L ∞ bounds, the a priori estimate (58) as well as from the convergences (65) and (70). The strong convergences of (iii) and (iv) can be shown by applying (70) in (iii) and the generalized Aubin-Lions lemma with f (r τ , b τ ) = r τ b τ in (iv).
Analogous results hold for Eq. (64) which allows us to perform the limit τ → 0 giving a weak solution to system (41). The only thing which remains to verify is the entropy inequality (42). Since E is convex and continuous, it is weakly lower semicontinuous. Because of the weak convergence of (r τ (t), b τ (t)), h (r (t), b(t)) dx ≤ lim inf τ →0 h (r τ (t), b τ (t)) dx for a.e. t > 0.
We cannot expect the identification of the limit of

Conclusion
Gradient flow techniques provide a natural framework to study the behavior of time evolving systems that are driven by an energy. This energy is decreasing along solutions as fast as possible, a property inherent in nature. Hence, many partial differential equation models exhibit this structure. Most of these systems arise in the mean field limit of a particle system, which has a gradient structure itself. Passing from the microscopic level to the macroscopic equations often relies on closure assumptions and approximations, which perturb the original gradient flow structure.
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. We expect that the presented approach would also be extendable to a higher number of species. However, one would have to consider all types of pairwise interactions which would result in much more complicated systems.
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 Sect. 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 toward 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.