Global existence analysis of cross-diffusion population systems for multiple species

The existence of global-in-time weak solutions to reaction-cross-diffusion systems for an arbitrary number of competing population species is proved. The equations can be derived from an on-lattice random-walk model with general transition rates. In the case of linear transition rates, it extends the two-species population model of Shigesada, Kawasaki, and Teramoto. The equations are considered in a bounded domain with homogeneous Neumann boundary conditions. The existence proof is based on a refined entropy method and a new approximation scheme. Global existence follows under a detailed balance or weak cross-diffusion condition. The detailed balance condition is related to the symmetry of the mobility matrix, which mirrors Onsager's principle in thermodynamics. Under detailed balance (and without reaction), the entropy is nonincreasing in time, but counter-examples show that the entropy may increase initially if detailed balance does not hold.


Introduction
Shigesada, Kawasaki, and Teramoto suggested in their seminal paper [24] a diffusive Lotka-Volterra system for two competing species, which is able to describe the segregation of the population and to show pattern formation when time increases. Starting from an on-lattice random-walk model, this system was extended to an arbitrary number of species in [30,Appendix]. While the existence analysis of global weak solutions to the two-species model is well understood by now [3,4], only very few results for the n-species model under very restrictive conditions exist (see the discussion below). In this paper, we provide for the first time a global existence analysis for an arbitrary number of population species using the entropy method of [15], and we reveal an astonishing relation between the monotonicity of the entropy and the detailed balance condition of an associated Markov chain.
Here, u i models the density of the ith species, u = (u 1 , . . . , u n ), Ω ⊂ R d (d ≥ 1) is a bounded domain with Lipschitz boundary, and ν is the exterior unit normal vector to ∂Ω. The diffusion coefficients are given by a ik u s k , i, j = 1, . . . , n, where a i0 , a ij ≥ 0 and s > 0. The functions p i are the transition rates of the underlying random-walk model [16,30]. The source terms f i are of Lotka-Volterra type, b ij u j , i = 1, . . . , n, and we suppose that b i0 , b ij ≥ 0 (competition case). Note that (1) can be written more compactly as . . , f n (u)).
State of the art. From a mathematical viewpoint, the analysis of (1)-(2) is highly nontrivial since the diffusion matrix A(u) is neither symmetric nor generally positive definite. Although the maximum principle may be applied to prove the nonnegativity of the densities, it is generally not possible to show upper bounds. Moreover, there is no general regularity theory for diffusion systems, which makes the analysis very delicate. Equations (1) can be written in the form which allows for the proof of an L 2+s estimate by the duality method [8,21], but we will not exploit this method in the paper. The case of n = 2 species and linear transition rates s = 1 corresponds to the original population model of Shigesada, Kawasaki, and Teramota [24], (6) ∂ t u 1 − ∆ u 1 (a 10 + a 11 u 1 + a 12 u 2 ) = f 1 (u), The numbers a i0 are the diffusion coefficients, a ii are the self-diffusion coefficients, and a ij for i = j are called the cross-diffusion coefficients. This model attracted a lot of attention in the mathematical literature. The first global existence result is due to Kim [17] who studied the equations in one space dimension, neglected self-diffusion, and assumed equal coefficients (a ij = 1). His result was extended to higher space dimensions in [11]. Most of the papers made restrictive structural assumptions, for instance supposing that the diffusion matrix is triangular (a 21 = 0), since this allows for the maximum principle in the second equation [1,18,20]. Another restriction is to suppose that the cross-diffusion coefficients are small, since in this situation the diffusion matrix becomes positive definite [11,28]. Significant progress was made by Amann [1] who showed that a priori estimates in the W 1,p norm with p > d are sufficient for the solutions to general quasilinear parabolic systems to exist globally in time, and he applied his result to the triangular case. The first global existence result without any restriction on the diffusion coefficients (except positivity) was achieved in [14] in one space dimension and in [3,4] in several space dimensions. The results were extended to the whole space in [12]. The existence of global classical solutions was proved in, e.g., [19], under suitable conditions on the coefficients.
As already mentioned, there are very few results for more than two species. The existence of positive stationary solutions and the stability of the constant equilibrium was investigated in [2,23]. The existence of global weak solutions in one space dimension assuming a positive definite diffusion matrix was proved in [27], based on Amann's results. Using an entropy approach, the global existence of solutions was shown in [10] for three species under the condition 0 < s < 1/ √ 3 (which guarantees that det(A(u)) > 0). To our knowledge, a global existence theorem under more general conditions seems to be not available in the literature. In this paper, we prove such a result and relate a structural condition on the coefficients a ij with Onsager's principle of thermodynamics.
Key ideas. Before we state the main results, let us explain our strategy. The idea is to find a priori estimates by employing a Lyapunov functional approach with where π i > 0 are some numbers and Because of the connection of our method to nonequilibrium thermodynamics [16,Section 4.3], we refer to H[u] as an entropy and to h(u) as an entropy density. Introducing the so-called entropy variable w = (w 1 , . . . , w n ) (called chemical potential in thermodynamics) by equations (1) can be written as which yields gradient estimates for u α/2 i . This strategy was employed in many papers on cross-diffusion systems; see, e.g., [3,4,9,12,14,15,30]. In this paper, we introduce two new ideas which we explain for the case s = 1 (s = 1 is studied below).
It is known that the entropy (7) with π i = 1 is a Lyapunov functional for the two-species model (6) with f 1 = f 2 = 0. This property is generally not satisfied for the corresponding n-species system. Our first idea is to introduce the numbers π = (π 1 , . . . , π n ) in the entropy (7). It turns out that (7) is a Lyapunov functional and H(u)A(u) is symmetric and positive definite if (10) π i a ij = π j a ji for all i, j = 1, . . . , n.
More precisely, this property is equivalent to the symmetry of H(u)A(u) (see Proposition 19). We recognize (10) as the detailed balance condition for the Markov chain associated to (a ij ). The equivalence of the symmetry and the detailed balance condition is new but not surprising. In fact, the latter condition means that π is a reversible measure, and time-reversibility of a thermodynamic system is equivalent to the symmetry of the so-called Onsager matrix B(w), so symmetry and reversibility are related both from a mathematical and physical viewpoint. We detail these relations in Section 5.1. In Section 2.1, we derive a refined estimate for H(u)A(u) leading to and thus giving an . This is the key estimate for the global existence result. (Below we also take into account the reaction terms (4).) One may ask whether the detailed balance condition is necessary for the monotonicity of the entropy. It is not. We show that if self-diffusion dominates cross-diffusion in the sense (12) η and detailed balance may be not satisfied, then the estimate leading to (11) still holds (with different constants), and global existence follows. (Throughout this paper, we set π i = 1 when detailed balance does not hold.) However, if conditions (10) or (12) are both not satisfied, there exist coefficients a ij and initial data u 0 such that t → H[u(t)] is increasing on [0, t 0 ] for some t 0 > 0; see Section 5.3. Numerical experiments (not shown) indicate that after the initial increase, the entropy decays and, in fact, it stays bounded for all time. We conjecture that the entropy is bounded for all time for all nonnegative coefficients and nonnegative initial data and that global existence of weak solutions holds for any (positive) coefficients a ij . Our results can be extended to nonlinear transition rates of type (3). One may choose more general terms a ij u s j j with different exponents s j but the results are easier to formulate if all exponents are equal. Coefficients with exponents s = 1 were also considered in [9,10,15] but in the two-species case only. We generalize these results to the multispecies case for any n ≥ 2. The entropy method has to be adapted since the inverse of h ′ s (z) = (s/(s − 1))(z s−1 − 1) cannot be defined on R and thus, u(w) = (h ′ ) −1 (w) is not defined for all w ∈ R n . This issue can be overcome by regularization as in [9,15]. In fact, we introduce Then h ′ ε : (0, ∞) n → R n can be inverted and (h ′ ε ) −1 : R n → (0, ∞) n is defined on R n . As a consequence, u i = (h ′ ε ) −1 (w) i is positive for any w ∈ R n and even strongly positive if w varies in a compact subset of R n .
Unfortunately, the product H ε (u)A(u), where H ε (u) = h ′′ ε (u), is generally not positive definite and we need to approximate A(u). In contrast to the approximations suggested in [9,15], we employ a non-diagonal matrix; see (23) below. More specifically, we introduce A ε (u) = A(u) + εA 0 (u) + ε η A 1 (u) with non-diagonal A 0 (u), diagonal A 1 (u), and η ≤ 1/2 such that The choice of the non-diagonal approximation satisfying this inequality is nontrivial, and this construction is our second idea.
Main results. First, we show that global existence of weak solutions holds for linear transition rates (s = 1). In the following, we set Q T = Ω × (0, T ).
The lower bound s > 1 − 2/d can be avoided if the regularity u i ∈ L 2+s (Q T ) holds, which is expected to follow from the duality method [8,21]. Unfortunately, this method is not compatible with our approximation scheme (see (24) below). This issue can possibly be overcome by employing the scheme proposed in [10] which is specialized to diffusion systems like (5). In this paper, however, we prefer to employ scheme (24).
The paper is organized as follows. Section 2 is concerned with the positive definiteness of the matrices H(u)A(u) and H ε (u)A ε (u). The existence theorems are proved in Sections 3 and 4, respectively. In the final Section 5, we detail the connection between the detailed balance condition and the symmetry of H(u)A(u), prove a nonlinear Aubin-Lions compactness lemma needed in the proof of Theorem 2, and show that the entropy may be increasing initially for special initial data.

Positive definiteness of the mobility matrix
We derive sufficient conditions for the positive definiteness of the matrix H(u)A(u). Let The following result is valid for any s > 0.
Lemma 3. Let s > 0. Then, for any z ∈ R n and u ∈ R n + , Proof. The elements of the matrix H(u)A(u) equal Therefore, for z ∈ R n , The sum I 1 is the same as the first term on the right-hand side of (19), and I 2 equals the first part of the last term on this right-hand side. The remaining terms are written as The second term corresponds to the second term on the right-hand side of (19). Thus, it remains to prove that For this, we employ twice the inequality b 2 + c 2 ≥ 2bc: This finishes the proof. Lemma 4 (Detailed balance). Let 0 < s ≤ 1 and π i a ij = π j a ji for all i = j. Then, for all z ∈ R n and u ∈ R n + , Proof. The sum of the terms I 1 and I 2 in (20) is exactly the first term on the right-hand side of (21). Using detailed balance, we find that The sum of the first three terms equal the second term on the right-hand side of (21), and the remaining two terms are nonnegative since s ≤ 1.

Remark 5.
In the existence proof, we will choose z i = ∇u i (with a slight abuse of notation). Then the first term in (21) gives an estimate for ∇u . If a ii = 0, we lose the latter regularity. This loss can be compensated by the last term in (21) giving and consequently a bound for ∇(u i u j ) s/2 in L 2 . This observation is used in Remark 12.
It is possible to show the positive definiteness of H(u)A(u) without any restriction on (a ij ) (except positivity) if we restrict the choice of the parameter s; see the following lemma.
Lemma 7. Let a ij + a ji > 0 for i, j = 1, . . . , n and 0 < s ≤ s 0 , where Then, for all z ∈ R n and u ∈ R n + , Proof. We choose π i = 1 for i = 1, . . . , n. With the notation of the proof of Lemma 3, we only need to show that I 3 + I 4 ≥ 0. Employing the inequality b 2 + c 2 ≥ 2bc, we find that and this expression is nonnegative if s ≤ s 0 .
2.2. Superlinear transition rates. Again, we assume first that detailed balance holds.
Lemma 8 (Detailed balance). Let s > 1 and π i a ij = π j a ji for all i = j. If Proof. It is sufficient to estimate the sum I 3 + I 4 , defined in the proof of Lemma 3: This expression simplifies because of the detailed balance condition: and we end up with from which we conclude the result.
Lemma 10 (Non detailed balance). Let s > 1 and let Proof. We choose π i = 1 for i = 1, . . . , n. Then, as in the previous proof, By definition of η 2 , the result follows.

Approximate matrices.
Our theory requires that the range of the derivative h ′ equals R n . Since this is not the case if s = 1, we need to approximate the entropy density and consequently also the diffusion matrix. The approximate entropy density possesses the property that the range of its derivative is R n . We set H(u) = h ′′ (u) = (δ ij sπ i u s−2 i ) i,j=1,...,n for its Hessian and where η < 1/2 and The approximation ε η A 1 (u) is needed to achieve bounds for ε (η+1)/2 ∇u i in L 2 , which are necessary for the limit ε → 0. The off-diagonal terms in A 0 (u) are needed to preserve the entropy structure in the sense that H ε (u)A ε (u) is still positive definite. This is shown in the following lemma.
Lemma 11. Let s > 0. Then, for all z ∈ R n and u ∈ R n + , Proof. We decompose the product H ε (u)A ε (u) as The ε 2 -term becomes We obtain for z ∈ R n : Next, we consider the ε-terms: Summing these expressions and neglecting some positive contributions, we find that Here we see how we constructed A 0 ij (u): The off-diagonal coefficients are chosen in such a way that the mixed terms in z i z j cancel, and the diagonal elements (namely µ i ) are sufficiently large to obtain positive definiteness of H 0 (u)A 0 (u). Finally, we have (H ε (u)A 1 (u)) ij = δ ij (sπ i u s−1 i + ε) and which proves the lemma.
Step 1: solution of an approximated problem. Given w k−1 ∈ L ∞ (Ω; R n ) for k ∈ N, we wish to find w k ∈ H m (Ω; R n ) such that is a partial derivative of order m. If k = 1, we define w 0 = h ′ (u 0 ). Equation (24) is an implicit Euler discretization of (1) including an H m regularization term.
Remark 12 (Detailed balance and vanishing self-diffusion). In the detailed balance case, we may allow for vanishing self-diffusion. If a ii = 0 but a i0 > 0, Lemma 4 implies that only ∇(u (τ ) i ) 1/2 is bounded in L 2 (Q T ). This situation was considered in [4] for the two-species case, and we sketch the generalization to the n-species case.
Applying the Gagliardo-Nirenberg inequality similarly as in Step 2 of the previous proof, we conclude that (u and thus, (u ). This loss of regularity is problematic for the estimate of the discrete time derivative of u (τ ) i . In order to compensate this, we need the last sum in (21). Indeed, Remark 5 shows that for any . We infer from the Gagliardo-Nirenberg inequality that (u Next we exploit the structure of the equations, Thus, to show that A ij (u (τ ) )∇u (τ ) j is bounded, we only need to verify that ∇(u The estimate for the Lotka-Volterra term is more delicate since we have only the regularity u (τ ) i ∈ L 1+1/d (Q T ). Here, we need to suppose that b ii > 0, since this assumption provides an estimate for (u . By the Aubin-Lions lemma, there exists a subsequence (not relabeled) such that, as (ε, τ ) → 0, The problem now is to show that (a subsequence of) the discrete time derivative of u (τ ) i converges to ∂ t u i since L 1 (0, T ; W m,q (Ω) ′ ) is not reflexive. The idea is to apply a result from [29] which provides a criterium for weak compactness in L 1 (0, T ; X), where X is a reflexive Banach space. For details, we refer to [4].

Nonlinear transition rates: proof of Theorem 2
The strategy of the proof is similar to the proof of Theorem 1 but the nonlinear transition rates complicate the proof significantly. As outlined in Section 2.3, we approximate the entropy density by (22) and the diffusion matrix by (23). Again, we assume without loss of generality that u 0 Step 1: solution of an approximated problem. We employ the transformation w i = ∂h ε /∂u i and define B ε (w) = A ε (u(w))H ε (u(w)) −1 . Given w k−1 ∈ L ∞ (Ω; R n ), we wish to find w k ∈ H m (Ω; R n ) solving The construction of h ε ensures that Hypothesis (H1) of [15] is satisfied. By Lemma 11, Hypothesis (H2) holds as well. Also Hypothesis (H3) holds true since, for some C f > 0, where σ = 1 if s > 1 and 0 ≤ σ ≤ max{0, 2s − 1 + 2/d} if s < 1. We apply Lemma 5 in [15] to deduce the existence of a weak solution w k ∈ H m (Ω; R n ) to the above problem, which satisfies the discrete entropy inequality Setting u k := u(w k ) and employing Lemma 11, the second integral can be estimated as follows: where C s = s −1 (s + 1) min{a 11 π 1 , . . . , a nn π n , η 0 , η 1 π 1 , . . . , η 1 π n , η 2 }.
To finish this step, we wish to write the "very weak" formulation for the solution u (τ ) , which is defined from u k as in the previous section. First, we observe that Next, we choose a test function φ = (φ 1 , . . . , φ n ) ∈ L q (0, T ; W m,q ν (Ω)), where m > max{1, d/2} and q ≥ 2 will be determined below. Recall that W m,q ν (Ω) is defined in (17). Integrating by parts in (35), u (τ ) solves Step 2: uniform estimates. Arguing as in Step 2 of the proof of Theorem 1, we obtain from (36) and (37) for suffiently small τ > 0 the following uniform estimates.
Then, performing the limit (ε, τ ) → 0 in (38) with φ ∈ L ∞ (0, T ; W m,∞ ν (Ω)) , it follows that u solves (18) for such test functions. A density argument shows that, in fact, u solves (18) for φ ∈ L q (0, T ; W m,q ν (Ω)), finishing the proof. Remark 16 (Weak formulation). In the superlinear case s > 1, the solution constructed in the previous proof satisfies (1) even in the weak sense (13) with test functions φ ∈ L q (0, T ; W 1,q (Ω)). In order to see this, it is sufficient to show that Because of the structure of A ij , we only need to verify that Indeed, we have the convergences u (τ ) i → u i strongly in L γ (Q T ) for any 2 < γ < p(s) and (u (τ ) i ) s ⇀ u s i weakly in L 2 (0, T ; H 1 (Ω)) and hence, u choosing q ′ = 2γ/(γ + 2) > 1. For the remaining convergence, we need to integrate by parts. It holds for φ i ∈ L q (0, T ; W 2,q ν (Ω)) that A density argument shows that the weak formulation also holds for φ i ∈ L q (0, T ; W 1,q (Ω)).
Remark 17 (Vanishing self-diffusion). Assume that a i0 > 0 and a ii = 0. The difficulty is that we obtain a uniform bound only for ∇(u . In order to compensate this loss of regularity, we need additional assumptions, namely either s > max{1, d/2} (superlinear rates); or 0 < s < 1, d = 1, and σ < s + 1 (sublinear rates). Under these conditions, the statement of Theorem 2 holds true.

Additional and auxiliary results
5.1. Detailed balance condition. We wish to interpret the detailed balance condition (10) and to explain how the numbers π i can be computed from the coefficients (a ij ). We assume that the coefficients are normalized in the sense that a ij ≥ 0 and k=1, k =j a kj ≤ 1 for all i, j. The idea is to use a probabilistic approach, interpreting the coefficients a ij as the transition rates between two discrete states i and j of the state space S := {1, . . . , n}. Then a ij = P(X k = j|X k−1 = i) is the conditional probability for a random variable X : N → S. This variable represents the Markov chain associated to the stochastic matrix Q = (Q ij ) i,j ∈ R n×n , defined by Q ij = a ij for i = j and Q ii = 1 − i=1, i =j a ij for i = 1, . . . , n. A Markov chain is called reversible if there exists a probability distribution π = (π 1 , . . . , π n ) on S (called an invariant measure) such that (45) π i a ij = π j a ji , i, j = 1, . . . , n.
The Markov chain can be interpreted as a directed graph, where the states i ∈ S are the nodes and the edges are labeled by the probabilities a ij going from state i to state j. The state space S can be partitioned into so-called communicating classes. We write i → j if there exist i 0 , i 1 , . . . , i n+1 ∈ S such that a i 0 ,i 1 a i 1 ,i 2 · · · a in,i n+1 > 0 for i 0 = i and i n+1 = j. We say that i communicates with j if both i → j and j → i. A set of states σ ⊂ S is a communicating class if every pair in σ communicates with each other. This defines an equivalence relation, and communicating classes are the equivalence classes.
Consider the following properties: (A1) For all i, j ∈ S, it holds that either a ij = a ji = 0 or a ij a ji > 0.
The detailed balance condition (45) implies (A1) and (A2). It is shown in [25] that the converse is true and that the invariant measure π can be constructed explicitly.
Proposition 18. Let (A1)-(A2) hold. Then there exists an invariant measure π = (π 1 , . . . , π n ) such that the detailed balance condition (45) is satisfied. Moreover, π can be computed explicitly by choosing an i 0 in each communicating class and defining π j for i 0 and j belonging in the same class by depending only on i 0 and j, where i 1 , i 2 , . . . , i n = j are such that a i k ,i k+1 > 0 for k = 0, . . . , n − 1.
The following result relates the detailed balance condition and the symmetry of the matrix H(u)A(u).
(ii) Detailed balance condition: π i a ij = π j a ji for i = j.
(iii) Symmetry: The matrix H(u)A(u) is symmetric.
Proof. The implication (i) ⇐ (ii) is shown in Proposition 18. The converse can be proved directly using the detailed balance condition. Finally, the equivalence (ii) ⇔ (iii) follows from an explicit calculation of H(u)A(u).
Remark 20. The equivalence of the symmetry of H(u)A(u) and the detailed balance condition is related to the Onsager principle of thermodynamics. Indeed, the diffusion is the vector of entropy variables, is the Onsager matrix which is symmetric, according to Onsager, if and only if the thermodynamic system is time-reversible. Time-reversibility means that the Markov chain associated to the matrix (a ij ) is reversible, and the symmetry of B is equivalent to the symmetry of H(u)A(u). Thus, the equivalence (ii) ⇔ (iii) corresponds to the equivalence of the symmetry of B and the time-reversibility. For details on the detailed balance principle in thermodynamics, we refer to [7].
If there exist uniform estimates for the gradient (∇u (τ ) ) and the discrete time derivative τ −1 (u (τ ) − σ τ u (τ ) ), then, by the Aubin-Lions theorem and under suitable conditions on the spaces, (u (τ ) ) is relatively compact in some L q space. In the case of nonlinear transition rates, we obtain uniform estimates only for (∇(u (τ ) ) s ), where s > 0. Then relative compactness follows from a nonlinear version of the Aubin-Lions theorem [5]. We recall a special case of this result.
Theorem 21 can be extended to the case s < 1/2 if (u (τ ) ) is additionally bounded in L ∞ (0, T ; L 1 (Ω)) which generally follows from the entropy inequality. This result is new.
One may ask if a similar result as above holds if the diffusion coefficients a i0 do not vanish, since they give positive contributions to the entropy production. The next lemma shows that the entropy may be increasing even if a i0 > 0 is chosen arbitrarily.