Multispecies cross-diffusions: from a nonlocal mean-field to a porous medium system without self-diffusion

Systems describing the long-range interaction between individuals have attracted a lot of attention in the last years, in particular in relation with living systems. These systems are quadratic, written under the form of transport equations with a nonlocal self-generated drift. We establish the localisation limit, that is the convergence of nonlocal to local systems, when the range of interaction tends to 0. These theoretical results are sustained by numerical simulations. The major new feature in our analysis is that we do not need diffusion to gain compactness, at odd with the existing literature. The central compactness result is provided by a full rank assumption on the interaction kernels. In turn, we prove existence of weak solutions for the resulting system, a cross-diffusion system of quadratic type.


Introduction
Aggregation models have been used to study a wide range of living systems arising in biology and ecology, such as prey-predator system, movement of animal herds, aggregation of cells, etc.The simplest example is the non-local one species aggregation equation where u(t, x) models the evolution of a single population as a function of time t ∈ R + and space x ∈ R d , and K(x) is a self-interaction potential.However, a lot of biological systems are composed of many interacting species and a generalisation of Eq. ( 1) to multiple species N ≥ 1 is written as follows: with u i (t, x) for i ∈ [ [1, N ]] the local density of the i-th species.Notice that we do not include diffusion, which is the major challenge of our study.The kernel functions K ij describe the self and cross-interaction of the species, which can be either attractive or repulsive.In the literature, typical choices for the kernel function are the attractive-repulsive Morse potentials [25], Gaussian potentials, compactly-supported or yet quadratic potentials [35].These different potentials allow to model aggregative phenomena in population dynamics such as swarming [25,35,12,3,2,4].
The non-local terms in Eqs. ( 1) and ( 2) offer various mathematical challenges.In the case of a single population, the model has been widely studied.Several properties are now established as existence of solutions [7,16], well-posedness of a solution [7], long time behaviour [10], possible blow-up [6,5].The case N > 1, where the model includes cross-diffusion, is a more recent subject.Most of the literature considers the equation with an additional diffusion term [34,8,23,27,1].In the case without diffusion, a mathematical theory has been studied in [20,21] for smooth kernel function, see also [14] for a system with singular Newtonian potentials in dimension one.In addition, the existence and uniqueness of measure solutions for symmetrisable systems, i.e., K ij = CK ji , is proved in [24] for N = 2. Another existence proof is presented in [9], Theorem 6.1.
Our interest here is about the localisation limit for Eqs.(1) and (2).Let us define, for ε > 0, The localisation limit consists in considering ε → 0, meaning that the kernel functions converge toward the Dirac delta distribution.The motivation of this limit comes from the derivation of Eqs. ( 1) and (2) from many-particle systems [32] when the number of particles goes to ∞.The localisation limit is the next natural step to recover the macroscopic system linked to the many-particle system.
For the single species model, a sketch of the proof of the localisation limit is proposed in [30].Additional studies consider the limit in the case where diffusion is present [32,33].In particular [33] presents the successive limits from an interacting particle system with Brownian motion toward a macroscopic equation without diffusion.More recently, proofs using a gradient-flow approach have also been proposed [13,15].A common hypothesis found in these papers is that the kernel function is an auto-convolution, i.e., there exists ρ such that with ρ(•) = ρ(−•).This hypothesis allows to recover a priori estimates from the classical entropy E(u ε ) = u ε ln u ε .In the case of multiple species N > 1, the localisation limit has been proven with additional diffusion [19,28].In [19], the limit is proven in the case of small initial data and in [28], the hypotheses made on the interaction kernel are weaker than (3), and only consider that for a given function p, This hypothesis, together with the diffusion term is enough to obtain estimates on the density.Up to our knowledge, only [9] worked out the limit without diffusion for two species, using tools based on Wasserstein distance, a method very different from the one we propose here, with different conditions on the interaction matrix.The Cahn-Hilliard limit is studied in [11], following the one-species case in [26].The degenerate two-species cross-diffusion system with growth is also obtained as the limit of a nonlocal system with growth and interaction via a velocity potential in [22].
The paper is organised as follows.For the sake of clarity, in Sec. 2, we first detail the sketch of the proof proposed in [30] in the case where the kernel function is of the form (3) with ρ a non-negative compactly supported function.A key point of the approach is to prove the compactness of the family (u ε * ρ ε ) ε>0 and therefore its convergence.In Sec. 3, we extend this localisation limit to the multiplespecies case (2).The hypotheses made on the interaction kernel are then slightly different and consist in considering kernels K ij that can be decomposed thanks to a matrix Given this decomposition, assuming in addition that the matrix A is of rank N allows to obtain compactness for the sequence N ]] and pass to the limit in the equations.While the assumption ( 5) is stronger than the assumption (4), it allows us to obtain the limit in the case at hand, where there is no diffusion in the system.Additionally, considering that the functions ρ i have a bounded moment of order d+2 2 permit to extend the result to non-compactly supported potentials.Finally, in Sec. 4, we illustrate the localisation limit in 2D thanks to numerical simulations.
Note that we assume throughout d ≥ 2, and let the simpler case d = 1 to the reader.

The one-species case
In order to explain the ideas needed for systems, we first focus on the simpler one-species case.With the assumption (3), Eq. ( 1) can be rewritten as where the convolution kernel ρ ε satisfies We are interested in the localisation limit of this equation, i.e., the case when ε → 0. We aim to establish that when ε → 0, the solution of (6) converges toward a solution of Existence, uniqueness as well as the localisation limit have already been stated, and the proofs sketched, in [30]; we revisit this note here and provide full details for these results.To avoid technicalities and keep the proofs simple, we moreover restrict ourselves to compactly supported kernels; this may be relaxed by assuming |x| d+2 2 ρ(x)dx < ∞, see Sec. 3. We assume that the initial data satisfies for some constant C independent of 0 < ε ≤ 1.
Remark 2 We notice that for any functions f and g we have the following -frequently used -identity

A priori estimates
Proposition 3 For ε ∈ (0, 1), we assume (7) and (9).Let u ε be solutions of (6), then for all t ∈ [0, T ], u ε ≥ 0 and, for some constants Moreover, we have the following estimates, uniform with respect to ε, For d = 2, the estimates are the same except (16) which is replaced by for any 1 ≤ p < ∞, see also Remarks 4 and 5. Proof.
Step 1. Positivity, L 1 bound and L 2 bound.Let us first prove that u ε ≥ 0 : we multiply the equation by −⊮ uε≥0 to get so that integrating in space we get Integrating the equation in space we get that u ε (t, x)dx = u 0 ε (x)dx which gives (14) and (13).In addition, we can remark that Thus multiplying (6) by ρε * ρ ε * u ε and integrating in space, the above equation gives using (12) and the Green's formula.After integrating in time we have Hence the estimate ( 17) is verified and Step 2. Second moment control.We integrate the equation multiplied by the weight |x| 2 and find |x| 2 u ε (x)dx remains positive, we may write After integration in time we conclude ds.
Applying the Cauchy-Schwarz inequality to the second term of the right-hand side and taking the square, we finally obtain which is finite thanks to Step 1.
Step 3. Lower bound for the entropy.Let us decompose We then bound each term, noticing first that if Step 4. Space compactness .We consider the classical entropy R d u ε ln(u ε ) and compute its derivative Then integrating in time, we find which gives the entropy bound and estimate (15).Additionally, thanks to the Gargliano-Nirenberg-Sobolev interpolation inequality, see [29], we have for θ ∈ [0, 1] and p such that Then if d ̸ = 2, for θ = 1 we have which shows the estimate (16).This concludes the proof of Proposition 3.
Remark 5 For d = 2 we take θ < 1 and q = 1 hence p

Compactness
Our next step is to prove space and time compactness.We first remark that, thanks to the bounds ( 14) and ( 15), we have which provides us with space compactness.To pass to the limit ε → 0 we need to obtain some compactness in time.To this aim, we compute the equation of with Thus We first have, for fixed t, and on a compact K ⊂ R d , thanks to the inequality (15) and to the Taylor-Lagrange inequality Then we take a mollifier sequence We evaluate the second term on the right-hand side using ( 19) For the first term of the right-hand side, we apply the Cauchy-Schwarz inequality, use (15) and a classical convolution inequality: so finally choosing for instance η = √ k we have satisfied the assumptions of the Weil-Kolmogorov-Frechet theorem on L 1 ((0, T ) × K), hence the sequence ρ ε * u ε is compact in this space.

Convergence
Thanks to the compactness obtained in the previous section we have for a subsequence In addition, from the uniform bound thanks to estimates (16), and the strong limit (10), we have It remains to pass to the limit in the equation.
, then multiplying Eq. ( 6) by ϕ and integrating by parts, we have where we have used the notation Φ = ∇ϕ.The convergence is immediate for the time derivative but more delicate for the nonlinear drift term.We decompose: When we integrate in time, the first term passes to the limit because of strong-weak limits, and we want to prove that the second term vanishes.For L Φ the 1-Lipschitz norm of Φ, we write since ρ ε is compactly supported.We apply the Cauchy-Schwarz inequality Hence at the limit u 0 is solution of (8).

The case of systems
We may now extend the above proof to a system of N non-local equations (2), that is where K ij ε are kernel functions which satisfy This models the evolution of a multi-species population, where each species i can sense another species j through the non-local operator K ij .The aim of this section is to establish that its limit as ε → 0 is the system of porous media equations This limit has been studied in the case N = 2 in [9] using Wasserstein metrics related methods, under the assumption that min{γ 11 , γ 22 } > γ 12 +γ 21

12
≥ 0 ; in [22], the limit is proved with n = 2 and γ ij = 1 for i, j = 1, 2. In this work the method used to prove the convergence is very different from the work in [9] as well as our hypothesis.
We now list a set of assumptions on the interaction function and on the initial data.For the initial data we assume that, with uniform bounds with respect to ε, For the interaction function, we consider kernels K ij that can be decomposed thanks to a matrix with the notation ρi (.) = ρ i (−.).Additionally, we assume that Under these assumptions, the matrix K ij (x) is not necessarily symmetric, but we remark that γ ij = p k=1 α ki α kj = (A T A) ij and that the matrix G := (γ ij ) i,j is symmetric positive definite.
It is then natural to define With these conditions, we may write and We will often use the following computation: with the notation A −1 = (β ij ) the left inverse of A, we may write hence We are ready to state our main theorem.
be a weak solution of (21), then for all i ∈ [[1, N ]] there exist subsequences (u i ε ) and (ρ i ε * u i ε ) that are not relabelled such that and This result is a generalisation of the result presented in Sec. 1 with two new features: the multispecies aspect of the population, and the fact that the interaction functions ρ i ε do not necessarily have a compact support.The proof is structured as follows.In Sec.3.1 we present a priori estimates.In Sec.3.2 we show compactness of (ρ i ε * u i ε ) ε , and finally in Sec.3.3 we show the convergence in Eq. ( 21).
Remark 7 An example of interaction kernel.The Gaussian function used in [28] satisfies assumptions (24)- (25).Here, γ ij is defined such that there exists A = (α ij ) invertible such that γ ij = p k=1 α ki α kj .We remark that . Note than in this case ρ i ε = ρi ε and ρ i ε does not depend on i.

A priori estimates
Proposition 8 Assume (23)- (25) and let ] be solutions of (21).Then, for all t ≥ 0 and i Moreover, we have the following estimates, uniform with respect to ε, For d = 2, the estimates are the same except (35) which is replaced by for any 1 ≤ p < ∞, see also Remarks 4 and 5 of the one species case.
Integrating in space gives , which gives us the positivity.Additionally, for a given integrating Eq. ( 21) in space gives 0 dx and we have estimate (32).Since ρ i ∈ L 1 (R d ), we deduce easily the estimate (33).The L 2 control comes from the following computation: on the one hand we have On the other hand, given the formula (26), we have Thus, inserting it in the previous equation, and after integrating in time, we find . This immediately proves (36) and we also have that N i=1 α ki ρ i ε * u i ε ∈ L ∞ (0, T ; L 2 (R d )) for all k = 1, ..., p.Moreover, given Eq. ( 27) we deduce Step 2. Second moment control.For a given i ∈ [[1, N ]], we compute . finally, after integration in time we find .
Applying the Cauchy-Schwarz inequality to the second term of the right-hand side of the equation and taking the square, we obtain the annouced estimate thanks to estimate (36) and the initial condition (23).
Step 3. Space compactness.We consider the classical entropy N i=1 R d u i ε ln(u i ε ) and compute its derivative where we have used (26).Then, after integrating in time we get

Compactness
Our next step is to prove space and time compactness.We first remark that, thanks to the bounds ( 33) and (34), we have which provides us with space compactness.To pass to the limit ε → 0 we need to obtain some compactness in time.To this aim, we compute the equation of ], with y)dy.Thus, we may write and Therefore This, together with (37), allows us to show that for a compact set K ⊂ R d and k > 0, we have Indeed considering a mollifier sequence η ), we have and Similarly, we have Finally, we obtain With the choice η = √ k we have the result.Then, given (37) and ( 39), for all i ∈ [[1, N ]], ρ i ε * u i ε satisfies the assumptions of the Weil-Kolmogorov-Frechet theorem on L 1 ((0, T ) × K), hence the sequence (ρ i ε * u i ε ) is compact in this space.

Convergence
Given the previous estimates, we can extract a subsequence (we do not relabel it for the sake of simplicity) such that In addition, from the uniform bound thanks to estimates (35), and the strong limit (40), we have We are left to show the convergence in Eq. ( 21).Let us take .., N .Then multiplying Eq. ( 21) by ϕ i and integrating by parts, we have where we have used the notation Φ i = ∇ϕ i .For the term on the right-hand side, the weak limit (42) shows that The term on the right-hand side can be rewritten as Then, we obtain a decomposition in two terms The first term passes to the limit because of strong-weak limits.Indeed for j, k ∈ [[1, N ]], we have successively and thus ϕ the Lipschitz constant for the function Φ i , we have .
Next, we use the Hölder inequality with p = d+2 d−2 , p ′ = d+2 4 (notice however that the case d = 2 is slighly different and left to the reader, then we choose p < ∞ and p ′ > 1, but close to, and the same arguments work) and we write where we have used the change of variable x → z = x−y ε .Since we have proved that Therefore, we have proved that at the limit the equation holds and thus u i 0 is solution of ( 22).

Numerical simulations
We illustrate the localisation limit in 2D thanks to numerical simulations.We consider a spatial square domain Ω = [x min , x max ] × [y min , y max ] ∈ R 2 taken large enough so that the solution (spatial) support stays far from the boundaries.The spatial domain is discretized into N x × N x regularly spaced points with space step ∆x and we consider four species.The numerical schemes used to solve the nonlocal model ( 21) and the limiting model (22) are adapted from [17,2].We consider interaction kernels of the form (24) where the ρ i (x) for i = 1 . . . 4 are chosen as Gaussian functions of variance σ 2 i : With this choice, K ij (x) is again a Gaussian function of variance σ 2 i + σ 2 j .Therefore, K ij are of the form: ) .
Note that thus defined, the coefficients α ij of the matrix A control the intensities of the interactions, while the variances σ 2 i control the interaction distance between the different species.We consider three cases: • case 1: Species-dependent repulsion intensities: (A T A) ij = min(i, j) (the matrix A is an upper triangular matrix with coefficient equal to 1) and same interaction distance for all species, i.e., σ i = σ j ∀(i, j) = 1 . . . 4. These conditions satisfy hypothesis (25).
• case 3: Same repulsion intensities for all species: (A T A) ij = 1 (matrix A is a matrix of rank 1) and species-dependent interaction distances, i.e., σ i ̸ = σ j , ∀i ̸ = j.These conditions violate the full rank condition for matrix A, therefore we do not satisfy hypothesis (25).
All species are initially supposed to be uniformly distributed in a ball of radius S = 2 centered in the center of the domain: ρ i (x) = 1 πS 2 1 B(0,S) (x), ∀ i = 1 . . .4, and we let the system evolve until time t = 100.We aim to study numerically the convergence of the nonlocal system (21) to the local system (22) as ε → 0 in the different cases for the choice of the interaction kernel.
4.1 Case 1: same interaction distances for all species σ 2 i = σ 2 j = 0.1, different interaction intensities (A T A) ij = min(i, j) We show in Fig. 1 the spatial distributions at time t = 100 obtained with the nonlocal model (21) for different values of ε (ε = 5, first column, ε = 1, second column and ε = 0.5, third column) and with the local model ( 22) (fourth column) for each of the 4 species (different lines).
As one can observe in Fig. 1 and as expected, the solution evolves radially in space.Moreover, for this choice of interaction kernel, we observe that the species with larger indices diffuse faster than the first species (compare the support of the solution from top to bottom).Indeed, for interaction intensities of the form (A T A) ij = min(i, j), species with larger indices interact more strongly with themselves than species with smaller indices, resulting in stronger repulsion.Comparing the columns from left to right, we observe that the nonlocal model is quite far from the local model for ε = 5 (first column), but seems to get closer as ε decreases.For large ε, we observe the concentration of the solution in rings close to the boundary of the support due to the Gaussian form of the interaction, while reducing ε diminishes this effect and the nonlocal model gets closer to the local one.
Since the solutions are spreading radially, an efficient way to characterize the dynamics is to compute the radial distribution g(λ), which gives the average density on rings of size dλ: and we denote by g 0 i (x, t) the quantity obtained in the same way with the solution of the local model u 0 i (x, t).In the left subplots of Fig. 2, we show the quantities g 0 i (λ, 100) (plain lines) g ε i (λ, 100) (dotted lines) for different values of ε (different colors) and for each of the 4 species (different subplots).We again observe the concentration of the solutions of the nonlocal model in rings for large ε (red curves), while the radial distributions of the nonlocal model get closer to the one of the local model when decreasing ε (compare the green and blue lines).In order to quantify the convergence of the non-local to the local model as ε → 0, we plot in Fig. 2 the discrete L 2 distance between u ε and u 0 at time 100 as function of ε for each of the 4 species.
As one can see, the nonlocal model seems to converge to the local one as ε decreases for all species, with a faster convergence rate for species 1 than for larger indices.Again, the ring effect is visible for large values of ε.

Case 2: different interaction distances as function of the species σ
We show in Fig. 3 the simulations obtained in Case 2 with σ 2 1 = 0.1, σ 2 2 = 0.2, σ 2 3 = 0.3, σ 2 4 = 0.4 (increasing interaction distance with increasing species indices).We adopt the same representation as in Fig. 1.It is noteworthy that the local model obtained as the limit ε → 0 of the nonlocal one is the same here as for Case 1.
As one can observe in Fig. 3, for ε = 5, a quite different behavior is obtained when the interaction distance is species-dependent compared to when the species interact at the same distance (compare the first column of this figure with the one of Fig. 1).In this case ε = 5 indeed, species with larger indices do not diffuse faster than species with lower indices as observed previously, although the coefficients γ ij are the same.This is due to the fact that the pairs interacting the stronger (large species indices) are also the ones interacting the farer, reducing de facto the interaction intensity (because of the choice of the Gaussian forms for ρ i ).As a result, all species diffuse less than in the previous case, but this effect is damped when localizing the interaction (i.e., decreasing ε).For small values of ε = 0.5 indeed, we recover a profile close to the one obtained with the local model (compare the last two columns of Fig. 3).
Regarding the convergence of the nonlocal to local model (Fig. 4), we observe that contrary to Case 1, species 4 seems to converge faster to the local model as ε → 0 than the other species.This is because species 4 is the one the most impacted by the interaction distance controlled by ε.
4.3 Case 3: different interaction distances as function of the species σ i ̸ = σ j , same interaction intensities (A T A) ij = 1, ∀(i, j) Finally, we aim to study the convergence of the nonlocal to the local model when the hypothesis (25) of our convergence theorem is not satisfied.Namely, we consider here the case where matrix A is of rank 1, i.e., all species interact with the same intensity (A T A) ij = 1∀(i, j).
In this case, the diffusion is even slower for species with large indices compared to Case 2 for ε = 5.It is noteworthy that for this case, all species behave in the same way in the localisation limit, as expected since all coefficients for the interactions are the same and the interaction distances do not play a role in the localisation limit.We still observe a quite good correspondence between the nonlocal and local model as ε decreases, even though our interaction kernel does not satisfy the full rank hypothesis of A. These results are confirmed in Fig. 6, where we observe that the error decreases for all species as ε → 0.
These numerical results suggest that the nonlocal model converges to the local one when ε → 0 in a more general framework than the one considered in this paper, i.e., without the full rank condition for A.

Conclusion and perspectives
We have established the convergence of solutions from non-local to local systems describing aggregation for both single and multi-species population, a subject that have attracted a lot of attention recently.The major new feature in our analysis is that we do not need diffusion to gain compactness, at odd with much of the existing literature including one of the most advanced result in [28].The central compactness result is provided by a full rank assumption on the interaction kernel.To remove this condition is a challenging question which is certainly possible, in some cases at least, in view of the numerical simulations of Sec.4.3.
In turn, we prove existence of weak solutions for the resulting system, a cross-diffusion system of quadratic type, Eq. ( 21).This system is always symmetric due to the use of an energy (gradient flow) structure which is fundamental to provide the necessary a priori estimates.To extend the type of interaction is also an interesting question.
Notice that the form of this system differs from a closely related class, namely which is also an active field of research, [18,31] for which the nonlocal to local convergence is also studied.