The boundedness-by-entropy method for cross-diffusion systems

The global-in-time existence of bounded weak solutions to a large class of physically relevant, strongly coupled parabolic systems exhibiting a formal gradient-flow structure is proved. The main feature of these systems is that the diffusion matrix may be generally neither symmetric nor positive semi-definite. The key idea is to employ a transformation of variables, determined by the entropy density, which is defined by the gradient-flow formulation. The transformation yields at the same time a positive semi-definite diffusion matrix, suitable gradient estimates as well as lower and/or upper bounds of the solutions. These bounds are a consequence of the transformation of variables and are obtained without the use of a maximum principle. Several classes of cross-diffusion systems are identified which can be solved by this technique. The systems are formally derived from continuous-time random walks on a lattice modeling, for instance, the motion of ions, cells, or fluid particles. The key conditions for this approach are identified and previous results in the literature are unified and generalized. New existence results are obtained for the population model with or without volume filling.


Introduction
Many applications in physics, chemistry and biology can be modeled by reaction-diffusion systems with cross diffusion, which describe the temporal evolution of the densities or mass fractions of a multicomponent system. Physically, we expect that the concentrations are nonnegative or even bounded (examples are given in section 2). Since generally no maximum principle holds for parabolic systems, the proof of these bounds is a challenging problem. A second difficulty arises from the fact that in many applications, the diffusion matrix is neither symmetric nor positive semi-definite.
In this paper, we present a general technique which allows us, under certain assumptions, to prove simultaneously the global existence of a weak solution as well as its boundedness from below and/or above. The key idea is to exploit the so-called entropy structure of the parabolic system, which is assumed to exist, leading at the same time to gradient estimates and lower/upper bounds. More specifically, we consider reaction-diffusion systems of the form subject to the boundary and initial conditions Here, u(t) = (u 1 , . . . , u n )(·, t) : → R n is a vector-valued function (n 1), representing the densities or mass fractions u i of the components of the system, A(u) = (a ij (u)) ∈ R n×n is the diffusion matrix and the reactions are modeled by the components of the function f : R n → R n . Furthermore, ⊂ R d (d 1) is a bounded domain with Lipschitz boundary and ν is the exterior unit normal vector to ∂ . The divergence div(A(u)∇u) and the expression (A(u)∇u) · ν consist of the components 3), we consider more general cross-diffusion systems satisfying the main assumptions stated below. By identifying the key elements of the existence analysis, we are able to extend significantly known results and to state a unified method. Before summarizing the state of the art of cross-diffusion systems and detailing our main results, we explain the key idea of our method and provide an illustrating example.

Idea of the method
The main assumption in this paper is that system (1) possesses a formal gradient-flow structure, i.e. (1) can be formulated as where B is a positive semi-definite matrix and δH/δu is the variational derivative of the entropy (or free energy) functional H[u] = h(u)dx. The function h : D → [0, ∞) is called the entropy density and it is assumed to be convex. Identifying δH/δu with its Riesz representative Dh(u) (the derivative of h) and introducing the entropy variable w = Dh(u), the above formulation can be understood as where B = B(w) = A(u)(D 2 h(u)) −1 ∈ R n×n and D 2 h(u) ∈ R n×n is the Hessian of h. For transforming back from the w-to the u-variable, we need to assume that Dh : D → R n is invertible such that u = (Dh) −1 (w). Here and in the following, we identify the vectors u, w ∈ R n , linked by w = Dh(u) and the solution w to (3) with u = u(w) := (Dh) −1 (w). Equation (3) is parabolic in the sense of [3], since the diffusion matrix B(w) is positive semidefinite and the matrix Du(w) = [D 2 h((Dh) −1 (w))] −1 in ∂ t u(w) = Du(w)∂ t w is positive definite.
The gradient-flow formulation has two important consequences. First, calculating the formal time derivative of H, (3) and integrating by parts yields where A : B = i,j a ij b ij for matrices A = (a ij ) and B = (b ij ). Thus, if f (u) · w 0 and since B(w) is assumed to be positive semi-definite, H is a Lyapunov functional. We refer to H as an entropy and to the integral of ∇w : B∇w as the corresponding entropy dissipation. Under certain conditions, it gives gradient estimates for u needed to prove the global-in-time existence of solutions to (2) and (3). We remark that the positive semi-definiteness of B(w) is in fact a consequence of the existence of an entropy. It was shown in [20,36] that both properties are equivalent and moreover, B(w) may be even symmetric.
Second, supposing that there exists a weak solution w to (3), the invertibility of Dh : D → R n shows that the original variable u = (Dh) −1 (w) satisfies u(·, t) ∈ D for t > 0. Thus, if D is bounded, we automatically obtain L ∞ bounds for u, without using a maximum principle. If D is only a cone, for instance D = (0, ∞) n , we conclude the positivity of u(t).
We call the above technique the boundedness-by-entropy method since it provides lower and/or upper bounds for the solutions to (1)-(2) by the use of the entropy density. Summarizing, the method is based on two main hypotheses: H1: There exists a function h ∈ C 2 (D; [0, ∞)) whose derivative is invertible on R n . This yields the bound u(·, t) ∈ D. H2: The matrix D 2 h(u) A(u) is positive semi-definite for all u ∈ D. This condition is necessary to derive a priori estimates for u.
Note that the positive semi-definiteness of D 2 h(u) A(u) is equivalent to that of B(w) = A(u)(D 2 h(u)) −1 since for all z ∈ R n , z D 2 h(u) A(u)

z = (D 2 h(u)z) B(w)(D 2 h(u)z). Hypothesis H2 avoids the inversion of D 2 h(u).
In fact, we need a stronger hypothesis than H2 since it does not allow us to infer gradient estimates. We need to suppose that D 2 h(u) A(u) is positive definite in such a way that we obtain L 2 gradient estimates for u m i , where m > 0 is some number. Moreover, we need an estimate for the time derivative of u i which makes it necessary to impose some growth conditions on the coefficients of A(u). We explain these requirements with the help of the following example.

An illustrative example
We consider a multicomponent fluid consisting of three components with mass fractions u 1 , u 2 and 1 − u 1 − u 2 and equal molar masses under isobaric, isothermal conditions. The model equals (1) with the diffusion matrix A(u) = 1 2 + 4u 1 + u 2 1 + 2u 1 u 1 2u 2 2 + u 2 (5) where we have chosen particular diffusivities to simplify the presentation (see section 2.1). This example is already treated in [34] but well illustrates the key elements of our method. Notice that the nonnegativity of u 1 and u 2 can be proved easily by a maximum principle argument but the proof of upper bounds is less clear. The logarithmic entropy density for u = (u 1 , u 2 ) ∈ D, where D = {(u 1 , u 2 ) ∈ (0, 1) 2 : u 1 + u 2 < 1} (7) satisfies Hypothesis H1, since the inverse transformation of variables gives u = (Dh) −1 (w) = e w 1 1 + e w 1 + e w 2 , e w 2 1 + e w 1 + e w 2 ∈ D, (8) where w = (w 1 , w 2 ) ∈ R 2 . Thus, once the existence of weak solutions w to the transformed system (3) is shown, we conclude the bounds 0 < u 1 , u 2 < 1 automatically by transforming back to the original variable. In particular, no maximum principle is used. In fact, the inverse transformation even shows that 1−u 1 −u 2 > 0, as required from the application. Furthermore, the matrix with δ(u) = u 1 u 2 (1 − u 1 − u 2 )(2 + 4u 1 + u 2 ) is symmetric and positive definite, thus fulfilling Hypothesis H2. However, the matrix B = A(u)(D 2 h(u)) −1 degenerates at u 1 = 0 or u 2 = 0 and we cannot expect to conclude gradient estimates for w from (4). Instead, it is more appropriate to derive these estimates for the original variable u. Indeed, calculating the time derivative of the entropy according to (4) and using ∇w = D 2 h(u)∇u, we find that dH dt + ∇u : and the entropy dissipation can be estimated according to Thus, assuming that the integral involving the reaction terms can be bounded uniformly in u, we obtain H 1 ( ) estimates for u 1/2 1 and u 1/2 2 . Using the boundedness of u i , a priori estimates for ∂ t u 1 and ∂ t u 2 can be proven, taking into account the particular structure of A(u).
This example shows that we need additional assumptions on the nonlinearities of system (1) in order to derive suitable a priori estimates, detailed in section 1.4.

Remark 1 (Notion of entropy).
There exists an intimate relation between the boundednessby-entropy method and non-equilibrium thermodynamics. In particular, the entropy variable w = Dh(u) is related to the chemical potentials of a mixture of gases and the special transformation (8) is connected with a special choice of thermodynamic activities; see appendix A for details. The entropy density defined above equals the negative thermodynamic entropy. Since the physical entropy is increasing and we wish to investigate nonincreasing functionals, we have reversed the sign as usual in entropy methods. Moreover, the logarithmic entropy (6) is motivated by Boltzmann's entropy for kinetic equations. For these reasons, we refer to the functional H[u] = h(u)dx as a (mathematical) entropy. We note that in some applications, free energy may be a more appropriate notion. Furthermore, we remark that the notion of entropy for parabolic systems, as used in this paper, is closely related to the entropy in hyperbolic equations. For details, we refer to, e.g. [36].

State of the art
We have already mentioned that the analysis of cross-diffusion systems is delicate since standard tools like maximum principles and regularity results generally do not apply. For instance, there exist Hölder continuous solutions to certain cross-diffusion systems which are not bounded and there exist bounded weak solutions which develop singularities in finite time [59]. In view of these counterexamples, it is not surprising that additional conditions are required to prove that (local in time) weak solutions are bounded and that they can be continued globally in time.
Ladyženskaya et al [39, chapter 7] reduced the problem of finding a priori estimates of local-in-time solutions u to quasilinear parabolic systems to the problem of deriving L ∞ bounds for u and ∇u. Under some growth conditions on the nonlinearities, the global-intime existence of classical solutions was shown. A fundamental theory of strongly coupled systems was developed by Amann [2]. He formulated the concept of W 1,p weak solutions and their local existence and proved that the solutions exist globally if their L ∞ and Hölder norms can be controlled. The above mentioned counterexamples show that the control on both norms is necessary. Le and Nguyen [42] proved that bounded weak solutions are Hölder continuous if certain structural assumptions on the diffusion matrix are imposed. The regularity of the solutions to systems with diagonal or full diffusion matrix was investigated in, for instance, [3,27,45,53]. The global existence of classical solutions was proved under certain conditions avoiding blow-up in some higher norms (which could be possible in case of weak solutions).
The boundedness of weak solutions to strongly coupled systems has been proved using various methods. Invariance principles were employed by Küfner [38] and Redlinger [54], requiring severe restrictions on the initial data. Truncated test functions, which are nonlinear in the solutions to 2 × 2 systems, were suggested by Le [40] who proved the boundedness under some structural assumptions. In the work of Lepoutre et al [43], the existence of bounded solutions to strongly coupled systems with spatially regularized arguments in the nonlinearities was shown using Hölder theory for nondivergence parabolic operators. Other methods are based on the derivation of L p bounds uniform in p and the passage to the limit p → ∞ (Moser-type or Alikakos-type iterations [1]).
The idea of proving the boundedness of weak solutions using the entropy density (6) was, to our best knowledge, first employed by Burger et al [9] in a size-exclusion model for two species (see section 2.1). It was applied to a tumor-growth model in [33] and extended to Maxwell-Stefan systems for fluids with arbitrary many components [34]. A different entropy density was suggested in [29] to prove L ∞ bounds in one space dimension. In fact, the idea of using entropy variables already appeared in the analysis of parabolic systems from non-equilibrium thermodynamics [20], originally used to remove the influence from the electric potential and goes back to the use of so-called Slotboom variables in semiconductor modeling [47, section 3.2].
In this paper, we identify the key elements of this idea and provide a general global existence result for bounded weak solutions to certain systems. Furthermore, the technique is applied to a novel claess of cross-diffusion systems derived from a random-walk master equation on a lattice, underlying the strength and flexibility of the method. The novelty of the results is detailed in section 1.5.

Main results
Our first main result concerns the existence of bounded weak solutions to (1)-(2) under some general structural assumptions. Motivated by the comments in sections 1.1 and 1.2, we impose the following hypotheses: . . , n) be such that for all z = (z 1 , . . . , z n ) ∈ R n and u = (u 1 , . . . , u n ) ∈ D,

H2 :
There exists a * > 0 such that for all u ∈ D and i, j = 1, . . . , n for which m j > 1, it holds that |a ij (u)| a * |α j (u j )|. H3: It holds A ∈ C 0 (D; R n×n ) and there exists C f > 0 such that for all u ∈ D, Hypothesis H1 shows that the inverse transformation u = (Dh) −1 (w) is well defined. If w(t) is a weak solution to (3), we conclude that u(t) = (Dh) −1 (w(t)) ∈ D, yielding the desired L ∞ bounds on u(t), since D is assumed to be bounded in Assumption H2 . Assumption H2 , which implies H2, is employed to prove a priori estimates for ∇u m i i . Hypothesis H2 is required to show a bound on ∂ t u i . The condition on f (u) in Hypothesis H3 is needed to derive an a priori estimate for the solution. For instance, if the entropy density is given by (6), we may choose f i (u) = u α i i (1 − u 1 − u 2 ) β i for α i , β i > 0 and i = 1, 2. It is related to the quasi-positivity assumption f i (u) 0 for all u ∈ (0, ∞) 2 with u i = 0 [5, section 6]. In contrast to the structural assumptions of [40], Hypotheses H1-H3 are easy to verify as soon as an entropy density is found (often motivated from the application at hand). Under the above conditions, the following result holds. Theorem 2 (General global existence result). Let Hypotheses H1, H2 , H2 and H3 hold and let u 0 ∈ L 1 ( ; R n ) be such that u 0 (x) ∈ D for x ∈ . Then there exists a weak solution u to (1)-(2) satisfying u(x, t) ∈ D for x ∈ , t > 0 and Note that since D is bounded by Hypothesis H2 , the theorem yields global L ∞ bounds on u. We remark that we may also assume the slightly weaker condition u 0 (x) ∈ D. Indeed, we may approximate u 0 by u 0 η (x) ∈ D satisfying u 0 η → u 0 a.e. apply the theorem to u 0 η and perform the limit η → 0. We refer to [18] for details.
For the proof of theorem 2, we first semi-discretize (3) in time with step size τ > 0 and regularize this equation by the expression ε( |α|=m (−1) m D 2α w + w), where ε > 0, D 2α is a partial derivative of order 2|α| and m ∈ N is such that H m ( ) → L ∞ ( ). This regularization ensures that the approximate solution w (τ ) is bounded. For the proof of approximate solutions, we only need Hypotheses H1, H2 and H3. The discrete version of the entropy identity (4) and Hypotheses H2 and H2 yield uniform estimates for ∇u (τ ) = ∇((Dh) −1 (w (τ ) )) and the discrete time derivative of u (τ ) . Then, by the discrete Aubin lemma from [26], the limit (τ, ε) → 0 can be performed yielding the existence of a weak solution to (1)- (2).
Assumptions H2 and H2 may be too restrictive in certain applications. For instance, it may be impossible to find a bounded set D ⊂ R n satisfying H1 or the inequality in H2 is not satisfied (see the examples in section 2). However, we show that our technique can be adapted to situations in which variants of Hypotheses H2 and H2 hold. To illustrate this idea, we choose a class of cross-diffusion systems modeling the time evolution of two population species. These models are derived from a random-walk master equation in the diffusive limit (see appendix B).
We consider two situations. In the first case, we assume that volume limitations lead to a limitation of the population densities (volume-filling case). Then the densities are nonnegative and bounded by a threshold value which is normalized to one. This situation occurs, for instance, in the volume-filling Keller-Segel model [51] and in ion-channel modeling [30]. The diffusion matrix in (1) reads as (see appendix B) where u 1 and u 2 are the densities of the species with bounded total density, u 1 + u 2 1 and The function q is related to the transition probability of a species to move from one cell to a neighboring cell and β > 0 is the ratio between the transition rates of both species. Biologically, q vanishes when the cells are fully packed, i.e. if u 1 + u 2 = 1, so q(0) = 0 and q is nondecreasing. When only one species is considered, the diffusion equation corresponds to the equation for the cell density in the volume-filling chemotaxis model [63]. The special case q(u 3 ) = u 3 was analyzed by Burger et al [9]. The novelty of the following theorem is the extension to nonlinear functions q. For the sake of presentation, we state the result for power functions q(u 3 ) = u s 3 (s 1) but also more general functions are allowed (see theorems 6 and 9).  (2) with diffusion matrix (9) and f = 0 satisfying 0 u 1 , u 2 1 and u 3 : where β 1 = 1, β 2 = β and ·, · is the dual product between H 1 ( ; R 2 ) and H 1 ( ; R 2 ). The initial datum is satisfied in the sense of H 1 ( ; R 2 ) .
Since we can write the expression in (10) for i = 1 as and similarly for i = 2, we see that (10) is the weak formulation of (1) with matrix (9).
We may also consider reaction terms f (u) satisfying a particular structure; see remark 8.
The key idea of the proof is to introduce the entropy density where u = (u 1 , u 2 ) ∈ D = {(u 1 , u 2 ) ∈ (0, 1) 2 : u 1 + u 2 < 1} and 0 < c < 1. The use of this entropy functional is new. A computation shows that Hypotheses H1 and H2 are satisfied but not H2 . Indeed, we will prove in section 4 (see (32)) that and the factors degenerate at u 1 + u 2 = 1 since q(0) = 0. From this, we are able to conclude a gradient estimate for u 3 but not for u 1 or u 2 . However, this is sufficient to infer the existence of solutions, employing an extension of Aubin's compactness lemma, which was first used in [9] (see appendix C).
In the second case, volume-filling effects are not taken into account. Then the diffusion matrix reads as (see appendix B) where p 1 and p 2 are related to the transition probabilities of the two species and ∂ i p j = ∂p j /∂u i . Note that each row of A(u) is the gradient of a function such that which allows for additional L 2 estimates for u 1 and u 2 , using the duality estimates of Pierre and Schmitt [52]. This observation was exploited in [23]. When the transition probabilities depend linearly on the densities, p i (u) = α i0 + α i1 u 1 + α i2 u 2 (i = 1, 2), we obtain the well-known population model of Shigesada, Kawasaki and Teramoto [57] with the diffusion matrix A(u) = α 10 + 2α 11 u 1 + α 12 u 2 α 12 u 1 α 21 u 2 α 20 + α 21 u 1 + 2α 22 The maximum principle implies that u 1 and u 2 are nonnegative. Less results are known concerning upper bounds. In fact, in one space dimension and with coefficients α 10 = α 20 , Shim [58] proved uniform upper bounds. Moreover, if cross-diffusion is weaker than selfdiffusion (i.e. α 12 < α 22 , α 21 < α 11 ), weak solutions are bounded and Hölder continuous [40]. Uniform in time L ∞ bounds on the solution were derived in [35], imposing some conditions on the diffusion coefficients. Upper bounds without any restrictions are not known so far (at least to our best knowledge). This model has received a lot of attention in the mathematical literature. One of the first existence results is due to Kim [37] who neglected self-diffusion (α 11 = α 22 = 0) and assumed equal coefficients (α ij = 1). The tridiagonal case α 21 = 0 was investigated by Amann [2], Le [41] and more recently, by Desvillettes and Trescaces [25]. Yagi [64] proved an existence theorem under the assumption that the diffusion matrix is positive definite (α 12 < 8α 11 , α 21 < 8α 22 , α 12 = α 21 ). The first global existence result without any restriction on the diffusion coefficients (except positivity) was achieved in [29] in one space dimension and in [15,16] in several space dimensions. The case of concave functions p 1 and p 2 , for instance, and i = 1, 2, u = (u 1 , u 2 ), 0 < s < 1, was analyzed by Desvillettes et al [23]. Global existence results for small data in one space dimension for s > 1 were proved in [4]. We are able to generalize these results to the case 1 < s < 4 and for 'large' initial data but we need to restrict the size of the cross-diffusion coefficients α 12 and α 21 . More general functions p 1 and p 2 are possible; see section 5. The extension to n species is currently under investigation. (15) below. Then there exists a weak solution u = (u 1 , u 2 ) to (1)- (2) with diffusion matrix (12) and

Theorem 4 (No volume-filling case). Let (14) hold with 1 < s < 4 and
Nonvanishing reaction terms f (u) can be treated if they satisfy an appropriate growth condition such that f (u (τ ) ) is bounded in some L p space with p > 1, where u (τ ) is a solution to an approximate problem. We leave the details to the reader.
The idea of the proof is to employ the entropy density where c > 0. This functional was also employed in [23] but under different conditions on the diffusivities. Hypothesis H2 is only satisfied if the restriction (1 − 1/s)α 12 α 21 α 11 α 22 holds. It is an open problem whether there exists another entropy density fulfilling Hypothesis H2 without any restriction on α ij (except positivity). In order to satisfy Hypothesis H1, we need to regularize the entropy density (15) as in [23]: This regularization is motivated from the population model with diffusion matrix (13), see section 2.2. The range of Dh ε equals R 2 , as required in Hypothesis H1. This regularization makes necessary to regularize also the diffusion matrix (in a different way than in [23]), Then the regularized product D 2 h ε (u)A ε (u) is positive semi-definite (Hypothesis H2) only if s < 4. This restriction may be improved by developing a better regularization procedure.

Novelty of the results
As already mentioned, the transformation to exponential entropy variables has been employed already in certain diffusion problems to prove the global existence of weak solutions, but all these results concern particular diffusion systems. The key novelty of this paper is the unification of the entropy approach and the generalization of these particular results. We are able to treat rather general classes of diffusion matrices (being not diagonal and not triangular) and possible degeneracies, not considered in the literature up to now. Compared to the results in the literature (in particular [9,15,16,20,33]), we treat not only much more general situations but also simplify the strategy of the existence proof. In fact, the proof relies on the algebraic Hypotheses H1-H3 and the analytic approximation procedure, thus separating in a clear way the algebraic and analytic parts of the existence analysis.
The assumptions of theorem 2 are satisfied for the tumor-growth model [33] and the Maxwell-Stefan equations [34]-detailed in section 2-, including the example in section 1.2. Furthermore, the diffusion matrix with entropy density (6) and domain (7) also satisfies Hypotheses H1, H2 and H2 . The model with this diffusion matrix describes the aggregation of two population species with cross-diffusion terms related to the drift term of the Keller-Segel system. We stress the fact that theorem 2 provides a global existence results for all these models, thus unifying known results and providing a general existence theorem. In fact, up to our knowledge, it is the first existence result for bounded weak solutions to cross-diffusion systems with an arbitrary number of equations under physically motivated hypotheses. A special case of theorem 3 is given by n = 2 and q(s) = s, modeling the ion transport through nanopores [10]. We have been able to generalize the existence result of [9] to a class of general nonlinear functions q (which model transition rates of the underlying lattice model). Theorem 3 can be also seen as a generalization of the volume-filling Keller-Segel model for one cell density [63] to two cell densities (but ignoring the equation for the chemoattractant density). Finally, we mention that the compactness result in lemma 13 slightly generalizes the corresponding result contained in [9] to the time-discrete case.
Theorem 4 provides a global existence result to generalizations of the population model (13) of Shegasada, Kawasaki and Teramoto. Existence results for functions (14) with s = 1 [15,16] or 0 < s < 1 [23] are known. We have been able to extend these results to the case 1 < s < 4. In fact, more general choices of these functions are possible; see (48). Moreover, we provide the first (formal) derivation of the population model of [57] and its extensions from a random-walk lattice model (see appendix B).
Note. Since this work was submitted, further improvements of entropy methods for crossdiffusion systems were made. We mention the paper [24] in which the population model with diffusion matrix (12) and functions (14) was investigated and the global existence of weak solutions for general s > 1 was shown. Furthermore, diffusion systems including both nonlinear functions p i and q were analyzed in [65].
The paper is organized as follows. In section 2, we consider some examples of crossdiffusion systems, studied in the literature and discuss the validity of our hypotheses. The general existence theorem 2 is shown in section 3. The proofs of theorems 3 and 4 are presented in sections 4 and 5, respectively. Some further results and open problems are discussed in section 6. The appendix is concerned with some relations of our method to non-equilibrium thermodynamics; the derivation of a general population diffusion model, containing the diffusion matrices (9) and (12) as special cases; and the proof of a variant of the Aubin compactness lemma, needed in the proof of theorem 3. To simplify the presentation, we write in the following L p (0, T ; H k ( )) instead of L p (0, T ; H k ( ; R n )).

Examples
We present some diffusion problems from physical and biological applications and discuss the validity of Hypotheses H1-H3. The examples illustrate that these hypotheses are not just abstract conditions but that they are satisfied by a number of systems arising from applications.

Examples with volume filling
Volume-filling model of Burger. The transport of ions in biological cells or in multicomponent fluid mixtures may be modeled by the Poisson-Nernst-Planck equations for the ion concentrations and the electric potential. They can be derived from microscopic particle models [48]. The rigorous derivation, however, is unclear when volume filling effects are included, an issue which naturally arises in ion transport through channels and pores. In [10], a modification of the classical Poisson-Nernst-Planck equations was formally derived from a lattice-based hopping model. Neglecting electric effects (in order to concentrate on cross-diffusive phenomena), the system consists of equation (1) with the diffusion matrix and D i > 0 are some constants. It was shown in [10] that the entropy with density , is a Lyapunov functional along solutions to (1) with f = 0. This entropy density satisfies Hypothesis H1 since The case of n = 2 species was investigated in [9]. Then the diffusion matrix becomes which is a special case of (9) with q(u 3 ) = u 3 and β = D 2 /D 1 . A computation shows that for Consequently, assuming without loss of generality that D 2 > D 1 (otherwise, change the indices 1 and 2), Hypothesis H2 is satisfied but not Hypothesis H2 since the quadratic form degenerates at ρ = 1. Burger et al have shown in [9] that a global existence analysis is still possible. We generalize this result to general functions q in theorem 3.
In the general case n > 2 and the case of equal diffusivity constants D 1 = D i for all i = 2, . . . , n, summing up all equations yields ∂ t ρ − ρ = 0 in with homogeneous Neumann boundary conditions. Thus, by the classical maximum principle, ρ is bounded from below and above and this yields L ∞ bounds for the components u i 0. The analysis of different diffusivity constants and n > 2 is much more delicate and an open problem.
Tumor-growth model. Jackson and Byrne [32] derived a continuous mechanical model for the growth of symmetric avascular tumors in one space dimension. They assumed that the tumorhost environment consists of the tumor cells, the extracellular matrix (ECM) and interstitial fluid (water). The mixture is supposed to be saturated, i.e. the volume fractions of tumor cells u 1 , the ECM u 2 and water u 3 sum up to one, u 1 + u 2 + u 3 = 1. The model is derived from the mass balance equations and the momentum balance equations for each phase, neglecting inertial and external forces. The pressures for the cell and ECM phase are assumed to be proportional to their respective volume fraction. A nonlinearity is introduced by supposing that the tumor cells increase the ECM pressure but not inversely. For a small cell-induced pressure coefficient (i.e. θ = β = 1 in [32]), the resulting diffusion matrix reads as .
With the entropy density (6), defined in the domain D which is given by (7), we compute for Thus, since u 2 1, Hypotheses H1, H2 and H2 are satisfied and the global existence result follows from theorem 3 if the reaction terms satisfy Hypothesis H3. This result was first proven in [33].
Maxwell-Stefan equations. The dynamics of a multicomponent gaseous mixture with vanishing barycentric velocity can be described by the Maxwell-Stefan equations, which model the diffusive transport of the components of the mixture. Applications arise in many fields, e.g. in sedimentation, dialysis, respiratory airways, electrolysis and chemical reactors [62]. Under some simplifying assumptions (ideal gas, isobaric and isothermal conditions, same molar masses for all components), the equations for the molar concentrations u 1 , . . . , u n are given by the mass balance equations and the reduced force balances [6], where D ij > 0 are the binary diffusion coefficients. The variables u i satisfy n j =1 u j = 1. The inversion of the flux-gradient relations is not straightforward since the linear system in J i has a singular matrix; see [5,46]. The idea of [34] was to replace the last component u n by the remaining ones by u n = 1 − n−1 j =1 u i and to analyze the remaining n − 1 equations. For n = 3 components, the inversion can be made explicit, leading to system (1) for the remaining n − 1 = 2 equations with the diffusion matrix where we abbreviated Employing the entropy density (6) as in the previous example, we compute Therefore, Hypothesis H1 and H2 are fulfilled with the domain D given by (7). Moreover, Hypothesis H2 is satisfied with m 1 = m 2 = 1 2 and also Hypothesis H2 holds. Theorem 2 yields the existence of global bounded weak solutions to this problem. It can be shown that this result holds true in the case of n > 3 components; see [34].
Multispecies chemotaxis equations. In order to explore the movement of different cell types in both crowded (dense) and uncrowded (dispersed) tissues, Painter assumed a heterogeneous tissue composed of two cell types defined on a discrete lattice [50]. The transition rates are chosen to describe a movement which takes place either through movement unimpeded into neighboring unoccupied space or by an active cell, which pulls forward, displacing a passive cell, that is pulled back, in neighboring occupied space via 'location-swapping'. The diffusion limit in the master equations leads to a cross-diffusion system for the cell densities u 1 and u 2 with a (simplified) matrix of the form where α 1 , α 2 > 0 are some parameters. In the model, we have neglected the chemotactic signal to simplify the presentation. The diffusion system possesses the entropy density (6) and it holds In particular Hypotheses H1 and H2 are satisfied. Moreover, we obtain an L 2 bound for ∇ √ u 1 u 2 . However, it is unclear how to derive a gradient estimate for u 1 and u 2 in the spirit of H2 . The global existence of weak solutions to (1), (2) with the above diffusion matrix is an open problem.

Models without volume filling
Population model of Shigesada, Kawasaki and Teramoto. Shigesada, Kawasaki and Teramoto [57] suggested a population model to describe the spatial segregation of interacting populations and to study the coexistence of two similar species. The authors assumed that the populations of two species prefer the same environment and that they are influenced by intra-specific and inter-specific population pressures. Formally, the equations can be derived from a randomwalk lattice model with transition rates which depend linearly on the population densities u 1 and u 2 (see appendix B). The diffusion system has the diffusion matrix where the coefficients α ij are nonnegative. Introducing the entropy density a computation shows that Hypothesis H2 is satisfied, i.e. for all Thus, we see that Hypothesis H2 with m i = 1 2 or m i = 1 and Hypothesis H2 hold true. Since D = (0, ∞) 2 is not bounded, we cannot apply theorem 2. However, as the coefficients of A(u) depend only linearly on u, the proof of theorem 4 can be adapted to this case and we conclude the existence of global nonnegative weak solutions to this problem. This statement was first proved in [15] using a different approximation procedure. In theorem 4 we generalize this result.
Semiconductor model with electron-hole scattering. In high-injection situations, the carrier transport through semiconductors is strongly affected by electron-hole scattering, i.e. by collisions between electrons and so-called defect electrons. This scattering needs to be taken into account on the mesoscopic level by electron-hole collision operators in the semiconductor Boltzmann equations [55]. Then, performing the diffusion limit in the Boltzmann system for electrons and holes, drift-diffusion equations can be derived which include cross-diffusion terms arising from the interaction of the electrons and holes. The equations for the electron density u 1 and the hole density u 2 possess the diffusion matrix where µ 1 and µ 2 denote the (positive) mobility constants for the electrons and holes, respectively. The global existence of weak solutions to this model was shown in [17].
Introducing the entropy density Thus, Hypothesis H2 is satisfied but not H2 since the quadratic form is not uniformly positive. However, as shown in [17], bounds on the electron and hole masses (i.e. the integrals of u 1 and u 2 ) together with estimates from the entropy dissipation yield an H 1 bound for u 1/2 1 and u 1/2 2 . Then, with the entropy variables w i = ∂h/∂u i = log u i , the existence of global weak solutions was proved in [17].

Proof of theorem 2
The proof is based on the solution of a time-discrete and regularized problem, for which only Hypotheses H1, H2 and H3 are needed.

Lemma 5 (Existence for the regularized system (18)). Let Assumptions
Proof. Let y ∈ L ∞ ( ; R n ) and δ ∈ [0, 1] be given. We solve first the following linear problem: where The forms a and F are bounded on H m ( ; R n ). The matrix B(y) = A(u(y))(D 2 h) −1 (u(y)) is positive semi-definite since, by Assumption H2, Hence, the bilinear form a is coercive: It holds S(y, 0) = 0 for all y ∈ L ∞ ( ; R n ). We claim that the operator S is continuous. Indeed, let δ η → δ in R and y η → y strongly in L ∞ ( ; R n ) as η → 0 and set w η = S(y η , δ η ) ∈ H m ( ; R n ). Then, by continuity, u(y η ) → u(y), B(y η ) → B(y) and f (u(y η )) → f (u(y)) strongly in L ∞ ( ; R n ) and the above coercivity estimate gives a uniform bound for (w η ) in H m ( ; R n ). This implies that, for a subsequence which is not relabeled, w η w weakly in H m ( ; R n ). Then, performing the limit η → 0 in (20), it follows that w = S(y, δ). In view of the compact embedding H m ( ; R n ) → L ∞ ( ; R n ), we infer that for a subsequence, S(y η , δ η ) = w η → w = S(y, δ) strongly in L ∞ ( ; R n ). Since the limit w is unique, the whole sequence (w η ) is converging which shows the continuity. Furthermore, the compact embedding H m ( ; R n ) → L ∞ ( ; R n ) yields the compactness of S.
It remains to prove a uniform bound for all fixed points of S(·, δ) in L ∞ ( ; R n ). Let w ∈ L ∞ ( ; R n ) be such a fixed point. Then w solves (20) with y replaced by w. With the test function φ = w, we find that The convexity of h implies that h(x) − h(y) Dh(x) · (x − y) for all x, y ∈ D. Choosing x = u(w) and y = u(w k−1 ) and using Dh(u(w)) = w, this gives Taking into account the positive semi-definiteness of B(w) and Assumption H3, (21) can be estimated as follows: Choosing τ < 1/C f , this yields an H m bound for w uniform in δ (not uniform in τ or ε). The Leray-Schauder fixed-point theorem shows that there exists a solution w ∈ H m ( ; R n ) to (20) with y replaced by w and δ = 1.

Proof of theorem 3
We impose more general assumptions on q than in theorem 3: Let q ∈ C 1 ([0, 1]) be positive and nondecreasing on (0, 1) such that For instance, the functions q(y) = y/(1 + y) s with 0 s 1 and q(y) = y s with 1 s 2 fulfill (28). The case q(y) = y s with any s 1 is also possible, see theorem 9.
Proof. First we observe that we may assume that β 1 since otherwise we rescale the equations by t → βt, which removes this factor in the second equation but yields the factor 1/β > 1 in the first equation.
To verify Assumption H2, we show the following lemma.

Lemma 7 (Positive definiteness of (D 2 h)A). The matrix (D 2 h)A is positive semi-definite.
Moreover, if yq (y) 2q(y) holds for all 0 y 1 (see (28)), it holds A straightforward computation shows that from which we infer the result.
Step 2: A priori estimates. By lemma 5, there exists a sequence (w k ) of weak solutions to (18) satisfying the entropy-dissipation inequality (19) and (29), the discrete entropy inequality (19) can be written as where u k = (u k 1 , u k 2 ) and u k Resolving the recursion as in the second step of the proof of theorem 2, we infer that The assumption q(y) κq (y) 2 (see (28)) implies an H 1 estimate for u k 3 : The boundedness of (u k i ) in L ∞ yields an L 2 estimate for (u k i ) and hence, also for (u k 3 ). Therefore, We cannot perform the simultaneous limit (τ, ε) → 0 as in the proof of theorem 2 since we need a compactness result like in lemma 13 (see appendix C) in which the discrete time derivative is uniformly bounded in H 1 ( ) and not in the larger space H m ( ) with m > d/2. Therefore, we pass to the limit in two steps.
Functions q(y) = y s with s > 2 can also be considered.

Theorem 9.
Let β > 1 and q(y) = y s with s 1. Let u 0 ∈ L 1 ( ) with 0 < u 0 < 1 in . Then there exists a weak solution to (1)-(2) satisfying the weak formulation (10), Proof. It is sufficient to consider the case s > 2 since the case s 2 is contained in theorem 6. Assumptions H1 and H2 hold since this requires only q(0) = 0 and q (s) 0. By lemma 5, there exists a sequence (w k ) of weak solutions to (18) satisfying the entropydissipation inequality (31), which reads as Observing that we infer, after resolving the above recursion, the uniform estimates Compared to the proof of theorem 2, we cannot conclude an H 1 bound for u k 3 but only for The function z → z 2/s for z 0 is Hölder continuous since s > 2. Therefore, by the lemma of Chavent and Jaffre [14, p. 141], Since the embedding W 2/s,s ( ) → L s ( ) is compact, we conclude the existence of a subsequence (not relabeled) such that u (ε) 3 → u 3 strongly in L s ( ). At this point, we can proceed as in Step 3 of the proof of theorem 2.
For the limit τ → 0, we need another compactness argument. Let u (τ ) i be defined as in the proof of theorem 2. We have the uniform estimates Then, by the generalization of the Aubin lemma in [19, theorem 3], we obtain, up to a subsequence, as τ → 0, u (τ ) 3 → u 3 strongly in L s (0, T ; L s ( )). The remaining proof is exactly as Step 4 of the proof of theorem 2.
Since the entropy density h(u), defined in (15), may not fulfill Hypothesis H1, we need to regularize: Then Dh ε (u) = Dh(u) + ε(log u 1 , log u 2 ) , and the range of Dh ε is R 2 . The Hessian of h ε equals showing that each component of Dh ε is strictly increasing and thus, (Dh ε ) −1 : R 2 → D = (0, ∞) 2 is well defined. Hence, h ε fulfills Hypothesis H1. We also regularize the diffusion matrix by setting The entropy density is regularized similarly in [23] but the diffusion matrix is regularized here in a different way.
Step 1: Verification of Hypotheses H1 and H2. Set H = D 2 h(u). First, we show that H A and H ε A ε are positive definite under additional conditions on the functions a ij . H = D 2 h(u), where h is defined by (15) and let A = A(u) be given by (12) with p i as in (47). If for some s 1, sa 12 (u 2 ) u 2 a 12 (u 2 ), sa 21 (u 1 ) u 1 a 21 (u 1 ) and a 11 (u 1 )a 22 (u 2 ) (1 − 1/s)a 12 (u 2 )a 21 (u 1 ) for all u 1 , u 2 0, then H A is positive definite and for all z = (z 1 , z 2 ) ∈ R 2 ,

Lemma 11 (Positive definiteness of H A and H ε A ε ). Let
Furthermore, if additionally 4a 21 (u 1 ) u 1 a 21 (u 1 ) 0 and 4a 12 (u 2 ) u 2 a 12 (u 2 ) 0 for all u 1 , u 2 0, then H ε A ε is positive definite and for all z = (z 1 , If p 1 and p 2 are given by (14), H A is positive definite if s 1, α 11 α 22 (1 − 1/s)α 12 α 21 and H ε A ε is positive definite if additionally s 4. The proof also works in the case s < 1; in this situation the restriction α 11 α 22 (1 − 1/s)α 12 α 21 is not needed. The last three terms are nonnegative for all z 1 , z 2 ∈ R if and only if a 11 a 21 0 and a 12 a 22 (1 − 1/s)a 21 a 12 . This shows the first statement of the lemma.

Further results and open problems
In this section, we discuss some further results and open problems related to cross-diffusion systems and entropy methods.
(a) Dirichlet boundary conditions: Dirichlet boundary conditions may be treated under some conditions on the nonlinearites. To see this, we assume that u i = u D,i on ∂ for i = 1, . . . , n and we set w D = Dh(u D ). For simplicity, we assume that u D = (u D,1 , . . . , u D,n ) depends on the spatial variable only. Then, using the test function If h(u) grows superlinearly (such as u i log u i ) and w D is smooth, we may estimate, after integration over time, the first integral on the left-hand side as (h(u) − w D u)dx 1 2 h(u)dx − C for some constant C > 0. Under Hypotheses H2 and H3, the second integral on the left-hand side and the second integral on the right-hand side are estimated as follows: Then, assuming that A(u) and f (u) are such that, for instance, we can estimate the remaining integrals on the right-hand side of (56). In a similar way, we may treat Robin boundary conditions by combining the ideas for Dirichlet and homogeneous Neumann conditions. A uniqueness result for bounded weak solutions to an n-species system related to (9) was recently presented in [65], but assuming that the transition rate function q is the same for all species. Uniqueness results under weaker conditions are an open problem (at least to our best knowledge). In fact, it is well known that this question is delicate since non-uniqueness even occurs for some scalar equations; see, e.g. [31]. (c) Quadratic reaction terms: In this paper, our focus was rather on the diffusion part than on the reaction part. The question is whether the results can be extended to diffusion systems with, say, quadratic reaction rates as they arise in reversible chemistry. For instance, a global existence result for the following problem is generally unknown: for i = 1, 2, 3, 4, with boundary and initial conditions (2).
, the global existence of weak solutions was shown in [22]. This result is based on the observation that h(u) = 4 i=1 u i (log u i − 1) possesses an L 2 estimate. This is proved by using the duality method of M. Pierre [52] and the proof is based on the diagonal structure of A. A global existence result for a system modeling the more general reaction A 1 + A 2 A 3 + A 4 was given in [11] in the two-dimensional case, the three-dimensional case being an open problem. Reaction-diffusion systems with diffusivities d i depending on u and quadratic rate functions were recently analyzed in [7], but still only in the diagonal case. (d) Long-time behavior of solutions: We expect that the long-time behavior of weak solutions to (1)-(2) with f = 0 can be proven under some additional conditions. For specific diffusion matrices, the exponential decay has been already shown; see [33,34]. The idea is to estimate the entropy dissipation in terms of the entropy. For instance, let u ∞ = (u ∞,1 , . . . , u ∞,n ) be a constant steady-state to (1)- (2). To simplify, we assume that the entropy is given by and that the entropy dissipation can be estimated from below according to for some λ > 0. Then the entropy-dissipation inequality (4) becomes dH dt + λH 0, and Gronwall's lemma implies exponential convergence in terms of the relative entropy H * . By the Csiszár-Kullback inequality, we conclude the exponential decay in the L 1 norm, u(t) − u ∞ L 1 ( ) C exp(−λt/2) for t > 0. For details, we refer to [12]. The main task is to derive the entropy-dissipation relation (57). According to Hypothesis H2 , we need to prove If α i (u i ) = u 1/2 i (for instance), this inequality is a consequence of the logarithmic Sobolev inequality. For more general functions α i and power-type entropy densities, one may employ Beckner-type inequalities; see [13, lemma 2]. The general case, however, is an open problem. (e) Entropies for population models: It is an open problem to find an entropy for the general population model with diffusion matrix A = (a ij ) defined in (B.5). Just a simple combination of the entropies for the systems analyzed in theorems 3 and 4, i.e. summing (11) and (15), seems to be not sufficient. A partial result was shown recently in [65]. It is also an open problem if an entropy for the population system with matrix (12), where p 1 and p 2 are given by (14), exists without restrictions on the coefficients α ij (except positivity). (f) Gradient systems and geodesic convexity: There might to be a relation between our entropy formulation and the gradient structure for reaction-diffusion systems developed by Liero and Mielke [44] but no results are known so far. Using purely differential methods, they have shown the geodesic λ-convexity for particular reaction-cross-diffusion systems. A generalized McCann condition for geodesic convexity of the internal energy of certain diffusion systems was given in [66].
It is an open problem if such geometric properties also hold for the examples presented in section 2. (g) Coupled PDE-ODE problems: Coupled reaction-diffusion-ODE problems occur, for instance, in chemotaxis-haptotaxis systems modeling cancer invasion [60]. Denoting by u 1 , u 2 , u 3 the densities of the cancer cells, matrix-degrading enzymes, extracellular matrix, respectively, the evolution is governed by (1) with the exemplary diffusion matrix It is an open problem whether these systems possess an entropy structure and whether our method can be extended to such problems. One idea could be to reduce the problem to the subsystem (u 1 , u 2 ) and to estimate terms involving ∇u 3 directly by differentiating the ODE for u 3 . (h) Bounded weak solutions: A challenging task is to determine all diffusion systems whose weak solutions are bounded. As a first step in this direction, we consider n = 2 and diffusion matrices A(u) = (a ij ) depending linearly on the variables, where α ij , β ij , γ ij ∈ R. We wish to find conditions on the coefficients for which the logarithmic entropy density (6) satisfies Hypothesis H2. Requiring that B = A(D 2 h) −1 is symmetric, we can fix some coefficients, The matrix (c ij ) = (D 2 h)A is positive semi-definite if and only if Other examples can be easily constructed. For instance, for α 11 = α 22 = β 11 = β 12 = 1 we have and (c ij ) = (D 2 h)A is positive semi-definite if and only if which is the case if 0 c 22 < ∞. However, it seems to be difficult to solve inequalities (58)- (59) in the general situation. Possibly, techniques from quadratic optimization with inequality constraints and quantifier elimination may help. Some results using different ideas were presented in [35].

Acknowledgments
The author acknowledges partial support from the Austrian Science Fund (FWF), grants P22108, P24304 and W1245 and the Austrian-French Program of the Austrian Exchange Service (ÖAD). Part of this manuscript was written during the stay of the author at the King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi-Arabia. The author thanks Peter Markowich for his kind invitation and support

Appendix A. Relations to non-equilibrium thermodynamics
We show that the entropy variable w = Dh(u), defined in section 1, is strongly related to the chemical potentials of a fluid mixture and that the particular change of unknowns associated with the logarithmic entropy density (6) is related to a special choice of the thermodynamic activities.
We introduce first the thermodynamic setting. Consider a fluid consisting of N components with the same molar mass under isobaric and isothermal conditions. We write ρ i instead of u i to denote the mass density of the ith component. The evolution of the mass densities is governed by the mass balance equations where J i are the diffusion fluxes. We have assumed that the barycentric velocity vanishes and that there are no chemical reactions. Furthermore, we assume for simplicity that the total mass density is constant, N j =1 ρ j = 1. Let s(ρ) = s(ρ 1 , . . . , ρ N ) be the thermodynamic entropy of the system. Then the chemical potentials µ i are defined (in the isothermal case) by where µ 0 i is the Gibbs energy which generally depends on temperature and pressure. Since we have supposed an isobaric, isothermal situation, µ 0 i is constant and, for simplicity, we set µ 0 i = 0 for i = 1, . . . , N. In order to model non-ideal gases, it is usual in thermodynamics to introduce the thermodynamic activity a i and the activity coefficient γ i by µ i = µ 0 i + log a i = log a i , where a i = γ i ρ i . If γ i = 1, we recover the ideal-gas case. In the volume-filling case, Fuhrmann [28] has chosen γ i = 1 + N −1 j =1 a j for numerical purposes. Then

Appendix B. Formal derivation of two-species population models
We derive formally cross-diffusion systems from a master equation for a continuous-time, discrete-space random walk in the macroscopic limit. We consider only random walks on onedimensional lattices but the derivation extends to higher dimensions in a straightforward way. The lattice is given by cells x i (i ∈ Z) with the uniform cell distance h > 0. The densities of the populations in the ith cell at time t > 0 are denoted by u 1 (x i , t) and u 2 (x i , t), respectively. We assume that the population species u 1 and u 2 move from the ith cell into the neighboring (i ± 1)th cells with transition rates S ± i and T ± i , respectively. Then the master equations can be formulated as follows [49]: where i ∈ Z. We further suppose that the transition rates S ± i and T ± i depend on the departure cell i and the arrival cells i ± 1: where σ 0 > 0 is some number. Abbreviating p j (x i ) = p j (u 1 (x i ), u 2 (x i )) and q j (x i ) = q j (1 − u 1 (x i ) − u 2 (x i )), the master equations (B.1)-(B.2) become The functions p 1 (x i ) and p 2 (x i ) model the tendency of the species to leave the cell i, whereas q 1 (x i ) and q 2 (x i ) describe the probability to move into the cell i. The latter functions allow us to model the so-called volume-filling effect. Indeed, we may interpret u 1 and u 2 as volume fractions satisfying u 1 + u 2 1. Then 1 − u 1 − u 2 describes the volume fraction not occupied by the two species. If the ith cell is fully occupied, i.e. u 1 (x i ) + u 2 (x i ) = 1 and q j (1 − u 1 (x i ) − u 2 (x i )) = q j (0) = 0, the probability to move into the ith cell is zero.
It seems very difficult-if not impossible-to explore the entropy structure of this equation in full generality. Therefore, we investigated two special cases in this paper. First, we assumed that p 1 = p 2 = 1, q 1 = q and q 2 = βq, where β > 0. This corresponds to the volume-filling case (9). Second, we have set q 1 = q 2 = 1. This gives the matrix (12).
Because of the relative compactness of (y τ ) in L 2 , the right-hand side converges to zero uniformly in τ > 0. This finishes the proof.