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

Marie Doumic, Sophie Hecht, Benoît Perthame, Diane Peurichard

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marie Doumic, Sophie Hecht, Benoît Perthame, Diane Peurichard. Multispecies cross-diffusions: from a nonlocal mean-field to a porous medium system without self-diffusion. Journal of Differential Equations, inPress, 389, pp.228-256. 10.1016/j.jde.2024.01.017. hal-04108050v2

HAL Id: hal-04108050

https://hal.science/hal-04108050v2

Submitted on 9 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Multispecies cross-diffusions: from a nonlocal mean-field to a porous medium system without self-diffusion

Marie Doumic* Sophie Hecht† Benoît Perthame‡ Diane Peurichard§

June 9, 2023

## Abstract

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.

2010 Mathematics Subject Classification. 35Q92, 35K55, 35R09, 35B25.

Keywords and phrases. Aggregation equation; Cross-diffusion; Nonlocal interactions; Localisation limit; Mathematical biology; multispecies models

## 1 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

$$ \partial_t u - \text{div}[u \nabla K * u] = 0, \quad t \ge 0, x \in \mathbb{R}^d, \quad d \ge 2, \qquad (1) $$

where $u(t, x)$ models the evolution of a single population as a function of time $t \in \mathbb{R}_+$ and space $x \in \mathbb{R}^d$, and $K(x)$ is a self-interaction potential. However, a lot of biological systems are composed of many

<footnote>*Inria and CMAP, IP Paris, Ecole polytechnique, CNRS, 91128 Palaiseau cedex. Email: Marie.Doumic@inria.fr</footnote>

<footnote>†Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris. Email: sophie.hecht@sorbonne-universite.fr</footnote>

<footnote>‡Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris. Email: Benoit.Perthame@sorbonne-universite.fr</footnote>

<footnote>§Inria Paris, team MAMBA, Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris. Email: diane.a.peurichard@inria.fr</footnote>
interacting species and a generalisation of Eq. (1) to multiple species $N \ge 1$ is written as follows:

$$ \partial_t u^i - \operatorname{div}[u^i \nabla \sum_{j=1}^{N} K^{ij} * u^j] = 0, \quad t \ge 0, x \in \mathbb{R}^d, \qquad (2) $$

with $u^i(t, x)$ for $i \in [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 aggressive 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 $\varepsilon > 0$,

$$ K_{\varepsilon}(x) = \frac{1}{\varepsilon^d} K\left(\frac{x}{\varepsilon}\right), $$

The localisation limit consists in considering $\varepsilon \to 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 $\infty$. 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 $\rho$ such that

$$ K(x) = \tilde{\rho} * \rho(x), \qquad (3) $$

with $\tilde{\rho}(\cdot) = \rho(-\cdot)$. This hypothesis allows to recover a priori estimates from the classical entropy $\mathcal{E}(u_{\varepsilon}) = \int u_{\varepsilon} \ln u_{\varepsilon}$. 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$,

$$ \int_{\mathbb{R}^d} pK * p \ge 0. \qquad (4) $$

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 $\rho$ a non-negative compactly supported function. A key point of the approach is to prove the compactness of the family $(u_{\varepsilon} * \rho_{\varepsilon})_{\varepsilon>0}$ and therefore its convergence. In Sec. 3, we extend this localisation limit to the multiple-species 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 $A := (\alpha^{kj}) \in \mathcal{M}_{p,N}(\mathbb{R})$ with $p \ge N$ and $(\rho_{\varepsilon}^i)_{i \in [1,N]} \in (L^1(\mathbb{R}^d))^N$ under the form

$$K^{ij}(x) = \sum_{k=1}^{p} \alpha^{ki} \alpha^{kj} \tilde{\rho}^{i} * \rho^{j}(x), \quad \forall (i, j) \in [1, N]^2. \qquad (5)$$

Given this decomposition, assuming in addition that the matrix $A$ is of rank $N$ allows to obtain compactness for the sequence $(u_{\varepsilon}^i * \rho_{\varepsilon}^i)_{\varepsilon>0}$ for all $i \in [1, 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 $\rho^i$ have a bounded moment of order $\frac{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 \ge 2$, and let the simpler case $d = 1$ to the reader.

## 2 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

$$\partial_t u_{\varepsilon} - \operatorname{div}[u_{\varepsilon} \nabla \tilde{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}] = 0, \quad t \ge 0, x \in \mathbb{R}^d, \qquad (6)$$

where the convolution kernel $\rho_{\varepsilon}$ satisfies

$$\rho_{\varepsilon}(x) = \frac{1}{\varepsilon^d} \rho(\frac{x}{\varepsilon}), \quad \rho(x) \ge 0, \quad \int_{\mathbb{R}^d} \rho(x)dx = 1. \qquad (7)$$

We are interested in the localisation limit of this equation, i.e., the case when $\varepsilon \to 0$. We aim to establish that when $\varepsilon \to 0$, the solution of (6) converges toward a solution of

$$\partial_t u_0 - \operatorname{div}[u_0 \nabla u_0] = 0, \quad t \ge 0, x \in \mathbb{R}^d. \qquad (8)$$

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 $\int |x|^{\frac{d+2}{2}} \rho(x)dx < \infty$, see Sec. 3. We assume that the initial data satisfies

$$u_{\varepsilon}^{0} \ge 0, \quad \int_{\mathbb{R}^d} (1 + u_{\varepsilon}^{0} + |x|^2 + |\ln u_{\varepsilon}^{0}|)u_{\varepsilon}^{0}(x)dx \le C, \qquad (9)$$

for some constant $C$ independent of $0 < \varepsilon \le 1$.
**Theorem 1** For $\varepsilon \in (0, 1)$, let $\rho_\varepsilon$ be defined by (7) and $R > 0$ such that $\text{Supp}(\rho) \subset B(0, R)$. Assuming that $u_\varepsilon^0$ verifies (9), letting $(u_\varepsilon)_\varepsilon$ be a weak solution of (6), then there exists a converging subsequence (that is not relabelled) such that when $\varepsilon \to 0$, we have

$$ \rho_{\varepsilon} * u_{\varepsilon} \to u_0 \quad \text{strongly in} \quad L^p(0, T; L_{\text{loc}}^q(\mathbb{R}^d)) \quad \forall 1 \le p < 2, 1 \le q < \frac{2d}{d-2}, \quad (10) $$

$$ u_{\varepsilon} \to u_0 \quad \text{weakly in} \quad L_{\text{loc}}^{1}((0, T) \times \mathbb{R}^d), \quad (11) $$

and $u_0$ is a solution of (8) in the distributional sense.

In Sec. 2.1 we present a priori estimates, in Sec. 2.2 we show compactness of $(\rho_\varepsilon * u_\varepsilon)_\varepsilon$, and finally in Sec. 2.3 we show the convergence in Eq. (6).

**Remark 2** We notice that for any functions $f$ and $g$ we have the following - frequently used - identity

$$ \int (\rho * f)(x)g(x)dx = \iint \rho(x - y)f(y)g(x)dxdy = \iint \tilde{\rho}(y - x)g(x)dxf(y)dy = \int (\tilde{\rho} * g)(y)f(y)dy. \quad (12) $$

## 2.1 A priori estimates

**Proposition 3** For $\varepsilon \in (0, 1)$, we assume (7) and (9). Let $u_\varepsilon$ be solutions of (6), then for all $t \in [0, T]$, $u_\varepsilon \ge 0$ and, for some constants $C_1(T)$, $C_2(T)$ independent of $\varepsilon$

$$ \int_{\mathbb{R}^d} |x|^2 u_\varepsilon(t) \le C_1(T) \quad \text{and} \quad \int_{\mathbb{R}^d} u_\varepsilon(t) |\ln u_\varepsilon(t)| \le C_2(T). $$

Moreover, we have the following estimates, uniform with respect to $\varepsilon$,

$$ u_{\varepsilon} \in L^{\infty}(0, T; L^{1}(\mathbb{R}^{d})), \quad (13) $$

$$ \rho_{\varepsilon} * u_{\varepsilon} \in L^{\infty}(0, T; L^{1}(\mathbb{R}^{d})), \quad (14) $$

$$ \rho_{\varepsilon} * u_{\varepsilon} \in L^{2}(0, T; H^{1}(\mathbb{R}^{d})), \quad (15) $$

$$ \rho_{\varepsilon} * u_{\varepsilon} \in L^{2}(0, T; L^{\frac{2d}{d-2}}(\mathbb{R}^{d})), \quad \text{for } d \ge 3, \quad (16) $$

$$ \sqrt{u_{\varepsilon}} \nabla K_{\varepsilon} * u_{\varepsilon} \in L^{2}((0, T) \times \mathbb{R}^{d}). \quad (17) $$

For $d = 2$, the estimates are the same except (16) which is replaced by $\rho_\varepsilon * u_\varepsilon \in L^2(0, T; L^p(\mathbb{R}^d))$ for any $1 \le p < \infty$, see also Remarks 4 and 5.

**Proof.**

**Step 1.** Positivity, $L^1$ bound and $L^2$ bound. Let us first prove that $u_\varepsilon \ge 0$: we multiply the equation by $-\mathbb{1}_{\{u_\varepsilon \ge 0\}}$ to get

$$ \partial_t |u_\varepsilon|_- = \text{div} [|u_\varepsilon|_- \nabla_x (\tilde{\rho}_\varepsilon * \rho_\varepsilon * u_\varepsilon)], $$

so that integrating in space we get

$$ \frac{d}{dt} \int_{\mathbb{R}^d} |u_\varepsilon|_- dx = 0, $$
hence $|u_\epsilon(t, \cdot)|_- = |u_\epsilon(0, \cdot)| = 0$. Integrating the equation in space we get that $\int u_\epsilon(t, x)dx = \int u_\epsilon^0(x)dx$ which gives (14) and (13). In addition, we can remark that

$$ \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} |\rho_\epsilon * u_\epsilon(t, x)|^2 dx = \int_{\mathbb{R}^d} (\partial_t u_\epsilon)(\check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon) dx. $$

Thus multiplying (6) by $\check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon$ and integrating in space, the above equation gives

$$ \begin{aligned} \frac{1}{2} \frac{d}{dt} \int_{\mathbb{R}^d} |\rho_\epsilon * u_\epsilon(t, x)|^2 dx &= \int_{\mathbb{R}^d} (\check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon) \mathrm{div}[u_\epsilon \nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon] dx \\ &= - \int_{\mathbb{R}^d} \nabla (\check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon) \cdot (u_\epsilon \nabla (\check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon)) dx \\ &= - \int_{\mathbb{R}^d} u_\epsilon |\nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon|^2 dx, \end{aligned} $$

using (12) and the Green's formula. After integrating in time we have

$$ \frac{1}{2} \int_{\mathbb{R}^d} |\rho_\epsilon * u_\epsilon(t, x)|^2 dx + \int_0^T \int_{\mathbb{R}^d} u_\epsilon |\nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon|^2 dx \le \frac{1}{2} \int_{\mathbb{R}^d} |\rho_\epsilon * u_\epsilon^0(x)|^2 dx \le \frac{1}{2} \int_{\mathbb{R}^d} |u_\epsilon^0(x)|^2 dx. $$

Hence the estimate (17) is verified and $\rho_\epsilon * u_\epsilon \in L^\infty(0, T; L^2(\mathbb{R}^d))$.

**Step 2. Second moment control.** We integrate the equation multiplied by the weight $|x|^2$ and find

$$ \begin{aligned} \frac{d}{dt} \int_{\mathbb{R}^d} |x|^2 u_\epsilon dx &= - \int_{\mathbb{R}^d} 2[u_\epsilon x \cdot \nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon] dx \\ &\le 2 \left( \int_{\mathbb{R}^d} |x|^2 u_\epsilon(x) dx \right)^{1/2} \left( \int_{\mathbb{R}^d} u_\epsilon(x) |\nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon|^2 dx \right)^{1/2}. \end{aligned} $$

Since $W(t) := \int_{\mathbb{R}^d} |x|^2 u_\epsilon(x) dx$ remains positive, we may write

$$ \frac{d}{dt} W^{1/2} = \frac{1}{2W^{1/2}} \frac{d}{dt} W \le \left( \int_{\mathbb{R}^d} u_\epsilon(x) |\nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon|^2 dx \right)^{1/2}. $$

After integration in time we conclude

$$ W^{1/2}(t) \le W^{1/2}(0) + \int_0^t \left( \int_{\mathbb{R}^d} u_\epsilon(x) |\nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon|^2 dx \right)^{1/2} ds. $$

Applying the Cauchy-Schwarz inequality to the second term of the right-hand side and taking the square, we finally obtain

$$ \int_{\mathbb{R}^d} |x|^2 u_\epsilon(t, x) dx \le 2 \int_{\mathbb{R}^d} |x|^2 u_\epsilon(0, x) dx + 2t \int_0^t \int_{\mathbb{R}^d} u_\epsilon(s, x) |\nabla \check{\rho}_\epsilon * \rho_\epsilon * u_\epsilon|^2(s, x) dx ds, $$

which is finite thanks to Step 1.

**Step 3. Lower bound for the entropy.** Let us decompose $\int_{\mathbb{R}^d} u_\epsilon |\ln(u_\epsilon)|_- dx$ as follows:

$$ \int_{\mathbb{R}^d} u_\epsilon |\ln(u_\epsilon)|_- dx = \int_{\mathbb{R}^d} u_\epsilon |\ln(u_\epsilon)|_{-\mathbb{H}} u_\epsilon \ge e^{-|x|^2} dx + \int_{\mathbb{R}^d} u_\epsilon |\ln(u_\epsilon)|_{-\mathbb{H}} u_\epsilon \le e^{-|x|^2} dx = A + B $$
We then bound each term, noticing first that if $u_{\varepsilon} \ge e^{-|x|^2}$ then $|\ln(u_{\varepsilon})|_{-} = -\ln(u_{\varepsilon})\mathbb{1}_{u_{\varepsilon} \le 1} \le |x|^2$:

$$A = \int_{\mathbb{R}^d} u_{\varepsilon} |\ln(u_{\varepsilon})|_{-} \mathbb{1}_{u_{\varepsilon} \ge e^{-|x|^2}} dx \le \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon} dx < \infty,$$

and, noticing now that if $x \ge 1$ and if $y \le e^{-|x|^2}$ then $y|\ln(y)| \le e^{-|x|^2}|\ln(e^{-|x|^2})|$ because $s \mapsto s|\ln(s)|$ is increasing on $(0, e^{-1})$ (with a maximum on $e^{-1}$)

$$B = \int_{\mathbb{R}^d} u_{\varepsilon} |\ln(u_{\varepsilon})|_{-} \mathbb{1}_{u_{\varepsilon} \le e^{-|x|^2}} dx \le \int_{|x| \le 1} u_{\varepsilon} |\ln(u_{\varepsilon})|_{-} \mathbb{1}_{u_{\varepsilon} \le 1} dx + \int_{|x| \ge 1} |x|^2 e^{-|x|^2} dx < \infty.$$

**Step 4. Space compactness.** We consider the classical entropy $\int_{\mathbb{R}^d} u_{\varepsilon} \ln(u_{\varepsilon})$ and compute its derivative

$$\begin{align*}
\frac{d}{dt} \int_{\mathbb{R}^d} u_{\varepsilon} \ln(u_{\varepsilon}) dx &= \int_{\mathbb{R}^d} (1 + \ln(u_{\varepsilon})) \operatorname{div}[u_{\varepsilon} \nabla \check{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}] dx \\
&= - \int_{\mathbb{R}^d} \nabla (1 + \ln(u_{\varepsilon})) \cdot [u_{\varepsilon} \nabla \check{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}] dx \\
&= - \int_{\mathbb{R}^d} \nabla u_{\varepsilon} \cdot \nabla (\check{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}) dx \\
&= - \int_{\mathbb{R}^d} (\nabla u_{\varepsilon} * \rho_{\varepsilon}) \cdot \nabla (\rho_{\varepsilon} * u_{\varepsilon}) dx \\
&= - \int_{\mathbb{R}^d} |\nabla u_{\varepsilon} * \rho_{\varepsilon}|^2 dx.
\end{align*}$$

Then integrating in time, we find

$$\int_{\mathbb{R}^d} u_{\varepsilon} |\ln(u_{\varepsilon})| + \int_0^T \int_{\mathbb{R}^d} |\nabla \rho_{\varepsilon} * u_{\varepsilon}|^2 dx dt \le \int_{\mathbb{R}^d} u_{\varepsilon}^0 \ln(u_{\varepsilon}^0) + 2 \int_{\mathbb{R}^d} u_{\varepsilon} |\ln(u_{\varepsilon})|_{-},$$

which gives the entropy bound and estimate (15). Additionally, thanks to the Gargliano-Nirenberg-Sobolev interpolation inequality, see [29], we have for $\theta \in [0, 1]$ and $p$ such that $\frac{1}{p} = \theta(\frac{1}{2} - \frac{1}{d}) + \frac{1-\theta}{q}$ with $q = 2$:

$$\|\rho_{\varepsilon} * u_{\varepsilon}\|_{L^p(\mathbb{R}^d)} \le C \|\nabla \rho_{\varepsilon} * u_{\varepsilon}\|_{L^2(\mathbb{R}^d)}^{\theta} \|\rho_{\varepsilon} * u_{\varepsilon}\|_{L^q(\mathbb{R}^d)}^{1-\theta}.$$

Then if $d \ne 2$, for $\theta = 1$ we have

$$\|\rho_{\varepsilon} * u_{\varepsilon}\|_{L^2(0,T; L^{\frac{2d}{d-2}}(\mathbb{R}^d))} \le \|\nabla \rho_{\varepsilon} * u_{\varepsilon}\|_{L^2(0,T; L^2(\mathbb{R}^d))}$$

which shows the estimate (16). This concludes the proof of Proposition 3.

**Remark 4** Taking $\theta < 1$, we get $\rho_{\varepsilon} * u_{\varepsilon} \in L_t^{\frac{2}{\theta}}(L_x^p)$ with $\frac{1}{p} = \theta(\frac{1}{2} - \frac{1}{d}) + 1 - \theta$, hence $\rho_{\varepsilon} * u_{\varepsilon} \in L^p((0, T) \times \mathbb{R}^d)$ for $p = 2 + \frac{2}{d}$.

**Remark 5** For $d = 2$ we take $\theta < 1$ and $q = 1$ hence $p = \frac{1}{1-\theta} > 1$. We get $\rho_{\varepsilon} * u_{\varepsilon} \in L^{\frac{2}{\theta}}(0, T; L^p(\mathbb{R}^d))$. For $\theta = 2/3$ we have $\rho_{\varepsilon} * u_{\varepsilon} \in L_t^3(L_x^3)$.
## 2.2 Compactness

Our next step is to prove space and time compactness. We first remark that, thanks to the bounds (14) and (15), we have

$$ \rho_{\varepsilon} * u_{\varepsilon} \in L^{\infty}(0, T; L_{\text{loc}}^{1}(\mathbb{R}^{d})), \quad \nabla \rho_{\varepsilon} * u_{\varepsilon} \in L_{\text{loc}}^{1}((0, \infty) \times \mathbb{R}^{d}), \quad (18) $$

which provides us with space compactness. To pass to the limit $\varepsilon \to 0$ we need to obtain some compactness in time. To this aim, we compute the equation of $\rho_{\varepsilon} * u_{\varepsilon}$,

$$ \partial_t \rho_{\varepsilon} * u_{\varepsilon} = \operatorname{div} \underbrace{[\rho_{\varepsilon} * u_{\varepsilon} \nabla (\tilde{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon})]}_{:=Q_{\varepsilon}(t,x)}, \quad (19) $$

with $Q_{\varepsilon}(t, x) = \int \rho_{\varepsilon}(x - y)u_{\varepsilon}(y)\nabla \tilde{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}(t, y)dy$ and

$$ \begin{aligned} \|Q_{\varepsilon}(t, x)\|_{L^1(\mathbb{R}^d)} &\le \|\rho_{\varepsilon}\|_{L^1(\mathbb{R}^d)} \int u_{\varepsilon}(y) |\nabla \tilde{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}(t, y)| dy \\ &\le \frac{1}{2} (\|\sqrt{u_{\varepsilon}}\|_{L^2(\mathbb{R}^d)}^2 + \|\sqrt{u_{\varepsilon}}|\nabla \tilde{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}(t, y)\|_{L^2(\mathbb{R}^d)}^2). \end{aligned} $$

Thus

$$ \|Q_{\varepsilon}(t, x)\|_{L^1(0,T;L^1(\mathbb{R}^d))} \le \frac{1}{2} \left( T \|u_{\varepsilon}\|_{L^{\infty}(0,T;L^2(\mathbb{R}^d))}^2 + \|\sqrt{u_{\varepsilon}}|\nabla \tilde{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}(t, y)\|_{L^2(0,T;L^1(\mathbb{R}^d))}^2 \right) \le CT. $$

Let us show that $\rho_{\varepsilon} * u_{\varepsilon}$ is compact in $L_{\text{loc}}^1((0, T) \times \mathbb{R}^d)$. We first have, for fixed $t$, and on a compact $K \subset \mathbb{R}^d$, thanks to the inequality (15) and to the Taylor-Lagrange inequality

$$ \|\rho_{\varepsilon} * u_{\varepsilon}(t, x+h) - \rho_{\varepsilon} * u_{\varepsilon}(t, x)\|_{L^1((0,T) \times K)} \le C(K,T)h \|\nabla u_{\varepsilon} * \rho_{\varepsilon}\|_{L^2((0,T) \times K)} \le C(K,T)h. $$

Then we take a mollifier sequence $\omega_{\eta}(x) = \frac{1}{\eta^d}\omega(\frac{x}{\eta})$ and write, for $v = u_{\varepsilon} * \rho_{\varepsilon}$,

$$ \|v(t+k, x) - v(t, x)\|_{L^1((0,T) \times K)} \le \| (v - v * \omega_{\eta})(t+k, x) - (v - v * \omega_{\eta})(t, x) \|_{L^1((0,T) \times K)} \\ + \|v * \omega_{\eta}(t+k, x) - v * \omega_{\eta}(t, x)\|_{L^1((0,T) \times K)}. $$

We evaluate the second term on the right-hand side using (19)

$$ \begin{aligned} \|v * \omega_{\eta}(t+k, x) - v * \omega_{\eta}(t, x)\|_{L^1((0,T) \times K)} &= \left\| \int_{t}^{t+k} \partial_t v * \omega_{\eta}(t+s, x) \mathrm{d}s \right\|_{L^1((0,T) \times K)} \\ &= \left\| \int_{t}^{t+k} \mathrm{div}(Q_{\varepsilon}) * \omega_{\eta}(t+s, x) \mathrm{d}s \right\|_{L^1((0,T) \times K)} \\ &= \left\| \int_{t}^{t+k} Q_{\varepsilon} * \mathrm{div}(\omega_{\eta})(t+s, x) \mathrm{d}s \right\|_{L^1((0,T) \times K)} \\ &\le \int_{t}^{t+k} \|Q_{\varepsilon} * \mathrm{div}(\omega_{\eta})(t+s, x)\|_{L^1((0,T) \times K)} \mathrm{d}s \\ &\le \int_{t}^{t+k} \frac{1}{\eta} \|Q_{\varepsilon}(t+s, \cdot)\|_{L^1((0,T) \times K)} \|\omega(t+s, \cdot)\|_{W^{1,1}((0,T) \times K)} \mathrm{d}s \\ &\le C(K)^{\frac{k}{\eta}} \end{aligned} $$

For the first term of the right-hand side, we apply the Cauchy-Schwarz inequality, use (15) and a classical convolution inequality:

$$ \|(v-v*\omega_{\eta})(t+k, x)-(v-v*\omega_{\eta})(t, x)\|_{L^1((0,T) \times K)}dt \le 2\eta C(\alpha) \sqrt{|K|(T+k)} \|\nabla v\|_{L^2((0,T+k) \times K)} \le C(K,T)\eta $$

so finally choosing for instance $\eta = \sqrt{k}$ we have satisfied the assumptions of the Weil-Kolmogorov-Frechet theorem on $L^1((0, T) \times K)$, hence the sequence $\rho_{\varepsilon} * u_{\varepsilon}$ is compact in this space.
## 2.3 Convergence

Thanks to the compactness obtained in the previous section we have for a subsequence

$$ 
\begin{array}{lll}
\rho_{\varepsilon} * u_{\varepsilon} \to u_0 & \text{strongly in} & L_{\text{loc}}^{1}((0, T) \times \mathbb{R}^d), \\
\nabla \rho_{\varepsilon} * u_{\varepsilon} \to \nabla u_0 & \text{weakly in} & L_{\text{loc}}^{1}((0, T) \times \mathbb{R}^d), \\
\partial_t \rho_{\varepsilon} * u_{\varepsilon} \to \partial_t u_0 & \text{weakly in} & L_{\text{loc}}^{1}((0, T) \times \mathbb{R}^d).
\end{array}
$$

In addition, from the uniform bound $\rho_{\varepsilon} * u_{\varepsilon} \in L^2(0, T; L^{\frac{2d}{d-2}}(\mathbb{R}^d))$ thanks to estimates (16), and the strong limit (10), we have

$$ \rho_{\varepsilon} * u_{\varepsilon} \to u_0 \quad \text{strongly in} \quad L^p(0, T; L_{\text{loc}}^q(\mathbb{R}^d)) \quad \forall 1 \le p < 2, 1 \le q < \frac{2d}{d-2}. \quad (20) $$

It remains to pass to the limit in the equation. Let $\phi \in C_c^{\infty}([0, T] \times \mathbb{R}^d)$, then multiplying Eq. (6) by $\phi$ and integrating by parts, we have

$$ \int_0^T \int_{\mathbb{R}^d} \partial_t u_{\varepsilon} \phi = \int_0^T \int_{\mathbb{R}^d} \text{div}[u_{\varepsilon} \nabla \check{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}] \phi, $$

thus

$$ \int_0^T \int_{\mathbb{R}^d} u_{\varepsilon} \partial_t \phi = \int_0^T \int_{\mathbb{R}^d} u_{\varepsilon} \nabla \check{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon} \Phi. $$

where we have used the notation $\Phi = \nabla \phi$. The convergence is immediate for the time derivative but more delicate for the nonlinear drift term. We decompose:

$$ 
\begin{aligned}
\int_{\mathbb{R}^d} \Phi(t, x) u_{\varepsilon} \nabla (\check{\rho}_{\varepsilon} * \rho_{\varepsilon} * u_{\varepsilon}) &= \int_{\mathbb{R}^d} \rho_{\varepsilon} * [\Phi(t, x) u_{\varepsilon}] \cdot \nabla (\rho_{\varepsilon} * u_{\varepsilon}) \\
&= \int_{\mathbb{R}^d} [\Phi(t, x) \rho_{\varepsilon} * u_{\varepsilon} + r_{\varepsilon}] \cdot \nabla (\rho_{\varepsilon} * u_{\varepsilon})
\end{aligned}
$$

with

$$ r_{\varepsilon}(t, x) = \int_{\mathbb{R}^d} [\Phi(t, y) - \Phi(t, x)] u_{\varepsilon}(t, y) \rho_{\varepsilon}(x - y) dy. $$

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_{\Phi}$ the 1-Lipschitz norm of $\Phi$, we write

$$ |r_{\varepsilon}(t, x)| \le L_{\Phi} \int_{\mathbb{R}^d} |x - y| \rho_{\varepsilon}(x - y) u_{\varepsilon}(t, y) dy \le \varepsilon L_{\Phi} \rho_{\varepsilon} * u_{\varepsilon}, $$

since $\rho_{\varepsilon}$ is compactly supported. We apply the Cauchy-Schwarz inequality

$$ \int_0^T \int_{\mathbb{R}^d} |r_{\varepsilon}(t, x)| ||\nabla(\rho_{\varepsilon} * u_{\varepsilon})| dt dx \le \varepsilon L_{\Phi} ||\rho_{\varepsilon} * u_{\varepsilon}||_{L^2(0,T; L^2(\mathbb{R}^d))} ||\nabla(\rho_{\varepsilon} * u_{\varepsilon})||_{L^2(0,T; L^2(\mathbb{R}^d))} \to 0. $$

Hence at the limit $u_0$ is solution of (8).
### 3 The case of systems

We may now extend the above proof to a system of N non-local equations (2), that is

$$ \partial_t u_{\varepsilon}^i - \operatorname{div}[u_{\varepsilon}^i \sum_{j=1}^N \nabla K_{\varepsilon}^{ij} * u_{\varepsilon}^j] = 0, \quad t \ge 0, x \in \mathbb{R}^d, \quad i \in [1, N], \quad (21) $$

where $K_{\varepsilon}^{ij}$ are kernel functions which satisfy

$$ K_{\varepsilon}^{ij}(x) = \frac{1}{\varepsilon^d} K^{ij}(\frac{x}{\varepsilon}), \quad \gamma^{ij} := \int_{\mathbb{R}^d} K^{ij}(x)dx = \int_{\mathbb{R}^d} K_{\varepsilon}^{ij}(x)dx < +\infty. $$

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 $\varepsilon \to 0$ is the system of porous media equations

$$ \partial_t u_0^i - \sum_{j=1}^N \operatorname{div}[\gamma^{ij} u_0^i \nabla u_0^j] = 0, \quad t \ge 0, x \in \mathbb{R}^d, \quad i \in [1, N]. \quad (22) $$

This limit has been studied in the case $N=2$ in [9] using Wasserstein metrics related methods, under the assumption that $\min\{\gamma^{11}, \gamma^{22}\} > \frac{\gamma^{12}+\gamma^{21}}{12} \ge 0$; in [22], the limit is proved with $n=2$ and $\gamma^{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 $\varepsilon$,

$$ u_{\varepsilon}^{i0} \ge 0, \quad |x|^2 u_{\varepsilon}^{i0} \in L^1(\mathbb{R}^d) \quad u_{\varepsilon}^{i0} \in L^1(\mathbb{R}^d) \cap L^2(\mathbb{R}^d), \quad u_{\varepsilon}^{i0} \ln(u_{\varepsilon}^{i0}) \in L^1(\mathbb{R}^d). \quad (23) $$

For the interaction function, we consider kernels $K^{ij}$ that can be decomposed thanks to a matrix $A := (\alpha^{kj}) \in \mathcal{M}_{p,N}(\mathbb{R})$ with $p \ge N$ and $(\rho_{\varepsilon}^i)_{i \in [1,N]} \in (L^1(\mathbb{R}^d))^N$ under the form

$$ K^{ij}(x) = \sum_{k=1}^{p} \alpha^{ki} \alpha^{kj} \tilde{\rho}^{i} * \rho^{j}(x), \quad \forall (i,j) \in [1, N]^2, \quad (24) $$

with the notation $\tilde{\rho}^i(\cdot) = \rho^i(\cdot, \cdot)$. Additionally, we assume that

$$ \operatorname{rank}(A) = N, \quad \int_{\mathbb{R}^d} \rho^i(x)dx = 1, \quad \rho^i(x) \ge 0, \quad \int_{\mathbb{R}^d} |x|^{\frac{d+2}{2}} \rho^i(x)dx < +\infty. \quad (25) $$

Under these assumptions, the matrix $K^{ij}(x)$ is not necessarily symmetric, but we remark that $\gamma^{ij} = \sum_{k=1}^p \alpha^{ki} \alpha^{kj} = (A^T A)_{ij}$ and that the matrix $G := (\gamma^{ij})_{i,j}$ is symmetric positive definite.

It is then natural to define

$$ \rho_{\varepsilon}^{i}(x) = \frac{1}{\varepsilon^{d}} \rho^{i}(\frac{x}{\varepsilon}) $$
which gives $K_{\varepsilon}^{ij} = \sum_{k=1}^{p} \alpha^{ki} \alpha^{kj} \tilde{\rho}_{\varepsilon}^{i} * \rho_{\varepsilon}^{j}$. With these conditions, we may write

$$ 
\begin{aligned}
\sum_{i,j=1}^{N} \int_{\mathbb{R}^d} f_{\varepsilon}^{i} K_{\varepsilon}^{ij} * f_{\varepsilon}^{j} dx &= \sum_{i,j=1}^{N} \sum_{k=1}^{p} \int_{\mathbb{R}^d} \alpha^{ki} \alpha^{kj} f_{\varepsilon}^{i} \tilde{\rho}_{\varepsilon}^{i} * \rho_{\varepsilon}^{j} * f_{\varepsilon}^{j} dx \\
&= \sum_{k=1}^{p} \int_{\mathbb{R}^d} \left( \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^{i} * f_{\varepsilon}^{i} \right) \left( \sum_{j=1}^{N} \alpha^{kj} \rho_{\varepsilon}^{j} * f_{\varepsilon}^{j} \right) dx \quad (26) \\
&= \sum_{k=1}^{p} \int_{\mathbb{R}^d} \left( \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^{i} * f_{\varepsilon}^{i} \right)^2 dx,
\end{aligned}
$$

and for $k \in [1, p]$ we have

$$ \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^{i} * f_{\varepsilon}^{i} = (A(\rho_{\varepsilon}^{i} * f_{\varepsilon}^{i})_{i})_{k} \quad (\text{as a matrix-vector product}). $$

We will often use the following computation: with the notation $A^{-1} = (\beta_{ij})$ the left inverse of $A$, we may write

$$ \rho_{\varepsilon}^{j} * f_{\varepsilon}^{j} = \sum_{k=1}^{p} \beta_{jk} \left( \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^{i} * f_{\varepsilon}^{i} \right), \quad (27) $$

hence

$$ \| \rho_{\varepsilon}^{j} * f_{\varepsilon}^{j} \|_{L^2(\mathbb{R}^d)} \leq \sum_{k=1}^{p} |\beta_{jk}| \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^{i} * f_{\varepsilon}^{i} \|_{L^2(\mathbb{R}^d)}. $$

We are ready to state our main theorem.

**Theorem 6** For $\varepsilon \in (0, 1)$, let $(\rho_{\varepsilon})_{i \in [1, N]}$ be defined by (24)-(25) and $(u_{\varepsilon}^{i0})_{i \in [1, N]}$ be a family of initial data satisfying (23). Let $(u_{\varepsilon}^{i})_{i \in [1, N]}$ be a weak solution of (21), then for all $i \in [1, N]$ there exist subsequences $(u_{\varepsilon}^{i})$ and $(\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i})$ that are not relabelled such that

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \to u_{0}^{i} \quad \text{strongly in } L^{p}(0, T; L_{\text{loc}}^{q}(\mathbb{R}^{d})) \quad \forall 1 \le p < 2, 1 \le q < \frac{2d}{d-2}, \quad (28) $$

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \to u_{0}^{i} \quad \text{weakly in } L^{2}(0, T; H^{1}(\mathbb{R}^{d})), \quad (29) $$

$$ u_{\varepsilon}^{i} \to u_{0}^{i} \quad \text{weakly in } L_{\text{loc}}^{1}((0, T) \times \mathbb{R}^{d}). \quad (30) $$

and $(u_{0}^{i})_{i \in [1, N]}$ is a solution of (22) in the distributional sense.

This result is a generalisation of the result presented in Sec. 1 with two new features: the multi-species aspect of the population, and the fact that the interaction functions $\rho_{\varepsilon}^{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 $(\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i})_{\varepsilon}$, 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]

$$ K_{\varepsilon}^{ij}(x) = \frac{\gamma^{ij}}{(2\pi\varepsilon^2)^{d/2}} \exp\left(-\frac{|x|^2}{2\varepsilon^2}\right) \quad (31) $$
satisfies assumptions (24)-(25). Here, $\gamma^{ij}$ is defined such that there exists $A = (\alpha^{ij})$ invertible such that $\gamma^{ij} = \sum_{k=1}^{p} \alpha^{ki} \alpha^{kj}$. We remark that

$$ 
\begin{aligned}
K_{\varepsilon}^{ij}(x) &= \gamma^{ij} \frac{e^{-\frac{|x|^2}{2\varepsilon^2}}}{(2\pi\varepsilon^2)^{d/2}} = \gamma^{ij} \int_{\mathbb{R}^d} \frac{e^{-\frac{|x-y|^2}{\varepsilon^2}}}{(\pi\varepsilon^2)^{d/2}} \frac{e^{-\frac{|y|^2}{\varepsilon^2}}}{(\pi\varepsilon^2)^{d/2}} dy \\
&= \sum_{k=1}^{p} \alpha^{ki} \alpha^{kj} \tilde{\rho}_{\varepsilon}^{i} * \rho_{\varepsilon}^{j}(x),
\end{aligned}
$$

with $\rho_{\varepsilon}^{i}(x) = \frac{e^{-\frac{|x|^2}{\varepsilon^2}}}{(\pi\varepsilon^2)^{d/2}}$. Note that in this case $\rho_{\varepsilon}^{i} = \tilde{\rho}_{\varepsilon}^{i}$ and $\rho_{\varepsilon}^{i}$ does not depend on $i$.

## 3.1 A priori estimates

**Proposition 8** Assume (23)-(25) and let $(u_{\varepsilon}^{i})_{i \in [1, N]}$ be solutions of (21). Then, for all $t \ge 0$ and $i \in [1, N]$, $u_{\varepsilon}^{i} \ge 0$ and, for some constants $C_1(T)$, $C_2(T)$ independent of $\varepsilon$

$$ \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^{i}(t) \le C_1(T) \quad \text{and} \quad \int_{\mathbb{R}^d} u_{\varepsilon}^{i}(t) |\ln u_{\varepsilon}^{i}(t)| \le C_2(T), \quad \forall t \in [0, T]. $$

Moreover, we have the following estimates, uniform with respect to $\varepsilon$,

$$ u_{\varepsilon}^{i} \in L^{\infty}(0, T; L^{1}(\mathbb{R}^{d})), \qquad (32) $$

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \in L^{\infty}(0, T; L^{1}(\mathbb{R}^{d})), \qquad (33) $$

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \in L^{2}(0, T; H^{1}(\mathbb{R}^{d})), \qquad (34) $$

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \in L^{2}(0, T; L^{\frac{2d}{d-2}}(\mathbb{R}^{d})), \quad \text{for } d \ge 3, \qquad (35) $$

$$ \sqrt{u_{\varepsilon}^{i}} \sum_{j=1}^{N} \nabla K_{\varepsilon}^{ij} * u_{\varepsilon}^{j} \in L^{2}((0, T) \times \mathbb{R}^{d}). \qquad (36) $$

For $d=2$, the estimates are the same except (35) which is replaced by $\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \in L^{2}(0, T; L^{p}(\mathbb{R}^{d}))$ for any $1 \le p < \infty$, see also Remarks 4 and 5 of the one species case.

**Proof.** Step 1. Positivity, $L^1$ bound and $L^2$ bound. Let us first prove that $u_{\varepsilon}^{i} \ge 0$ for $i \in [1, N]$. We multiply Eq. (21) by $-\mathbb{1}_{u_{\varepsilon}^{i} \ge 0}$ to get

$$ \partial_t |u_{\varepsilon}^{i}|_{-} = \text{div} [|u_{\varepsilon}^{i}|_{-} \nabla_x (K_{\varepsilon}^{ij} * u_{\varepsilon}^{i})]. $$

Integrating in space gives

$$ \frac{d}{dt} \int_{\mathbb{R}^d} |u_{\varepsilon}^{i}|_{-} dx = 0. $$

Hence $|u_{\varepsilon}^{i}(t, \cdot)|_{-} = |u_{\varepsilon}^{i0}|_{-} = 0$ for all $i \in [1, N]$, which gives us the positivity. Additionally, for a given $i \in [1, N]$, integrating Eq. (21) in space gives

$$ \frac{d}{dt} \int_{\mathbb{R}^d} u_{\varepsilon}^{i}(t, x) dx = \int_{\mathbb{R}^d} \sum_{j=1}^{N} \text{div}[u_{\varepsilon}^{i} \nabla K_{\varepsilon}^{ij} * u_{\varepsilon}^{j}] dx = 0. $$
Therefore $\int_{\mathbb{R}^d} u_{\varepsilon}^i(t, x)dx = \int_{\mathbb{R}^d} u_{\varepsilon}^{i0}dx$ and we have estimate (32). Since $\rho^i \in L^1(\mathbb{R}^d)$, we deduce easily the estimate (33). The $L^2$ control comes from the following computation: on the one hand we have

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \sum_{i,j=1}^{N} \int_{\mathbb{R}^d} u_{\varepsilon}^i K_{\varepsilon}^{ij} * u_{\varepsilon}^j dx &= \frac{1}{2} \sum_{i,j=1}^{N} \int_{\mathbb{R}^d} \partial_t u_{\varepsilon}^i K_{\varepsilon}^{ij} * u_{\varepsilon}^j dx + \frac{1}{2} \sum_{i,j=1}^{N} \int_{\mathbb{R}^d} u_{\varepsilon}^i K_{\varepsilon}^{ij} * \partial_t u_{\varepsilon}^j dx \\
&= \sum_{i,j=1}^{N} \int_{\mathbb{R}^d} \partial_t u_{\varepsilon}^i K_{\varepsilon}^{ij} * u_{\varepsilon}^j dx \\
&= - \sum_{i=1}^{N} \int_{\mathbb{R}^d} u_{\varepsilon}^i |\sum_{j=1}^{N} \nabla K_{\varepsilon}^{ij} * u_{\varepsilon}^j|^2 dx \le 0.
\end{aligned}
$$

On the other hand, given the formula (26), we have

$$ \sum_{i,j=1}^{N} \int_{\mathbb{R}^d} u_{\varepsilon}^i K_{\varepsilon}^{ij} * u_{\varepsilon}^j = \sum_{k=1}^{p} \int_{\mathbb{R}^d} \left( \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^i * u_{\varepsilon}^i \right)^2 . $$

Thus, inserting it in the previous equation, and after integrating in time, we find

$$ \sum_{k=1}^{p} \int_{\mathbb{R}^d} \left( \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^i * u_{\varepsilon}^i \right)^2 + \sum_{i=1}^{N} \int_{0}^{T} \int_{\mathbb{R}^d} u_{\varepsilon}^i |\sum_{j=1}^{N} \nabla K_{\varepsilon}^{ij} * u_{\varepsilon}^j|^2 \le \sum_{k=1}^{p} \int_{\mathbb{R}^d} \left( \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^i * u_{\varepsilon}^{i0} \right)^2 . $$

This immediately proves (36) and we also have that $\sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^i * u_{\varepsilon}^i \in L^{\infty}(0, T; L^2(\mathbb{R}^d))$ for all $k = 1, ..., p$. Moreover, given Eq. (27) we deduce

$$ \rho_{\varepsilon}^i * u_{\varepsilon}^i = \sum_{k=1}^{p} \beta_{ik} \left( \sum_{i=1}^{N} \alpha^{ki} \rho_{\varepsilon}^i * u_{\varepsilon}^i \right) \in L^{\infty}(0, T; L^2(\mathbb{R}^d)) \quad \forall i = 1, ..., N. $$

**Step 2. Second moment control.** For a given $i \in [1, N]$, we compute

$$ 
\begin{aligned}
\frac{d}{dt} \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^i &= \int_{\mathbb{R}^d} |x|^2 \operatorname{div}[u_{\varepsilon}^i (\sum_{j=1}^N \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^i)] \\
&= - \int_{\mathbb{R}^d} 2x u_{\varepsilon}^i (\sum_{j=1}^N \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^i) \\
&\le 2 \left( \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^i \right)^{1/2} \left( \int_{\mathbb{R}^d} u_{\varepsilon}^i \left| \left( \sum_{j=1}^N \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^i \right) \right|^2 \right)^{1/2}.
\end{aligned}
$$

Hence

$$ \frac{d}{dt} \left( \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^i \right)^{1/2} \le \left( \int_{\mathbb{R}^d} u_{\varepsilon}^i \left| \sum_{j=1}^N \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^i \right|^2 \right)^{1/2} . $$

finally, after integration in time we find

$$ \left( \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^i \right)^{1/2} \le \left( \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^{i0} \right)^{1/2} + \int_0^t \left( \int_{\mathbb{R}^d} u_{\varepsilon}^i \left| \sum_{j=1}^N \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^i \right|^2 \right)^{1/2} . $$
Applying the Cauchy-Schwarz inequality to the second term of the right-hand side of the equation and taking the square, we obtain the announced estimate

$$ \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^i \le 2 \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^{i0} + 2t \int_0^t \int_{\mathbb{R}^d} u_{\varepsilon}^i \left| \left( \sum_{j=1}^N \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^i \right) \right|^2 < +\infty, $$

thanks to estimate (36) and the initial condition (23).

**Step 3. Space compactness.** We consider the classical entropy $\sum_{i=1}^N \int_{\mathbb{R}^d} u_{\varepsilon}^i \ln(u_{\varepsilon}^i)$ and compute its derivative

$$ \begin{aligned} \frac{d}{dt} \sum_{i=1}^N \int_{\mathbb{R}^d} u_{\varepsilon}^i \ln(u_{\varepsilon}^i) &= - \sum_{i,j=1}^N \int_{\mathbb{R}^d} \nabla u_{\varepsilon}^i K_{\varepsilon}^{ij} * \nabla u_{\varepsilon}^j dx \\ &= - \sum_{k=1}^p \int_{\mathbb{R}^d} \left( \sum_{i=1}^N \alpha^{ki} \nabla \rho_{\varepsilon}^i * u_{\varepsilon}^i \right)^2 dx, \end{aligned} $$

where we have used (26). Then, after integrating in time we get

$$ \sum_{i=1}^{N} \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)| + \sum_{k=1}^{p} \int_{0}^{T} \int_{\mathbb{R}^d} \left( \sum_{i=1}^{N} \alpha^{ki} \nabla \rho_{\varepsilon} * u_{\varepsilon}^i \right)^2 dx dt \le \sum_{i=1}^{N} \int_{\mathbb{R}^d} u_{\varepsilon}^{i0} |\ln(u_{\varepsilon}^{i0})| + 2 \sum_{i=1}^{N} \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-}. $$

Moreover for a given $i \in [1, N]$, $\int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} dx$ can be decomposed as follows:

$$ \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} dx = \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} \mathcal{H}_{u_{\varepsilon}^i \ge e^{-|x|^2}} dx + \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} \mathcal{H}_{u_{\varepsilon}^i \le e^{-|x|^2}} dx. $$

We then bound each term, noticing first that when $u_{\varepsilon} \ge e^{-|x|^2}$ then $|\ln(u_{\varepsilon}^i)|_{-} = -\ln(u_{\varepsilon}^i) \mathcal{H}_{u_{\varepsilon}^i \le 1} \le |x|^2$

$$ \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} \mathcal{H}_{u_{\varepsilon}^i \ge e^{-|x|^2}} dx = \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} \mathcal{H}_{u_{\varepsilon}^i \ge e^{-|x|^2}} dx \le \int_{\mathbb{R}^d} |x|^2 u_{\varepsilon}^i dx < +\infty, $$

and, noticing now that when $x \ge 1$ and $y \le e^{-|x|^2}$ then $y|\ln(y)| \le e^{-|x|^2}|\ln(e^{-|x|^2})|$ because $x \mapsto x|\ln(x)|$ is increasing on $(0, e^{-1})$ (with a maximum on $e^{-1}$)

$$ \begin{aligned} \int_{\mathbb{R}^d} u_{\varepsilon} |\ln(u_{\varepsilon}^i)|_{-} \mathcal{H}_{u_{\varepsilon}^i \le e^{-|x|^2}} dx &= \int_{\mathbb{R}^d} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} \mathcal{H}_{u_{\varepsilon}^i \le e^{-|x|^2}} dx \\ &\le \int_{|x| \le 1} u_{\varepsilon}^i |\ln(u_{\varepsilon}^i)|_{-} \mathcal{H}_{u_{\varepsilon}^i \le 1} dx + \int_{|x| \ge 1} |x|^2 e^{-|x|^2} dx < +\infty. \end{aligned} $$

Then $\sum_{i=1}^N \alpha^{ki} \nabla \rho_{\varepsilon}^i * u_{\varepsilon}^i \in L^2(0, T; L^2(\mathbb{R}^d))$ for all $k = 1, ..., p$. Using again the inversion formula (27), we deduce

$$ \nabla \rho_{\varepsilon}^j * u_{\varepsilon}^j \in L^2(0, T; L^2(\mathbb{R}^d)) \quad \forall j = 1, ..., N, $$

which gives the estimates (34) (since we have already proved the corresponding $L^2$ bound). In addition using the Gargliano-Nirenberg-Sobolev inequality, see [29], we have

$$ \| \rho_{\varepsilon}^j * u_{\varepsilon}^j \|_{L^2(0,T; L^{\frac{2d}{d-2}}(\mathbb{R}^d))} \le \| \nabla \rho_{\varepsilon}^j * u_{\varepsilon}^j \|_{L^2(0,T; L^2(\mathbb{R}^d))}, $$

for $d \ge 3$, which shows estimate (35). This concludes the proof of Proposition 8.
## 3.2 Compactness

Our next step is to prove space and time compactness. We first remark that, thanks to the bounds (33) and (34), we have

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \in L^{\infty}(0, T; L_{\text{loc}}^{1}(\mathbb{R}^{d})), \quad \nabla \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \in L_{\text{loc}}^{1}((0, \infty) \times \mathbb{R}^{d}), \quad \forall i \in [1, N], \quad (37) $$

which provides us with space compactness. To pass to the limit $\varepsilon \to 0$ we need to obtain some compactness in time. To this aim, we compute the equation of $\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}$ for a given $i \in [1, N]$,

$$ \partial_t \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} = \operatorname{div}[\rho_{\varepsilon}^{i} * \left( u_{\varepsilon}^{i} \sum_{j=1}^{N} \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^{j} \right)], \\ := Q_{\varepsilon}^{i}(t,x) $$

with $Q_{\varepsilon}^{i}(t, x)$ defined by $Q_{\varepsilon}^{i}(t, x) = \int_{\mathbb{R}^d} \rho_{\varepsilon}^{i}(x-y)u_{\varepsilon}^{i}(y) \sum_{j=1}^{N} \nabla_x K_{\varepsilon}^{ij} * u_{\varepsilon}^{j}(y)dy$. Thus, we may write

$$ \begin{aligned} \|Q_{\varepsilon}^{i}(t, x)\|_{L^{1}(\mathbb{R}^{d})} &\le \|\rho_{\varepsilon}^{i}\|_{L^{1}(\mathbb{R}^{d})} \int_{\mathbb{R}^{d}} u_{\varepsilon}^{i}(y) \sum_{j=1}^{N} \nabla_{x} K_{\varepsilon}^{ij} * u_{\varepsilon}^{j}(y) dy \\ &\le \frac{1}{2} \left( \|\sqrt{u_{\varepsilon}^{i}}\|_{L^{2}(\mathbb{R}^{d})}^{2} + \|\sqrt{u_{\varepsilon}^{i}} \sum_{j=1}^{N} \nabla_{x} K_{\varepsilon}^{ij} * u_{\varepsilon}^{j}\|_{L^{2}(\mathbb{R}^{d})}^{2} \right), \end{aligned} $$

and

$$ \|Q_{\varepsilon}^{i}(t, x)\|_{L^{1}(0,T;L^{1}(\mathbb{R}^{d}))} \le \frac{1}{2} \left( T \|u_{\varepsilon}^{i}\|_{L^{\infty}(0,T;L^{1}(\mathbb{R}^{d}))} + \|\sqrt{u_{\varepsilon}^{i}} \sum_{j=1}^{N} \nabla_{x} K_{\varepsilon}^{ij} * u_{\varepsilon}^{j}\|_{L^{2}(0,T;L^{2}(\mathbb{R}^{d}))}^{2} \right). \quad (38) $$

Therefore $Q_{\varepsilon}^{i}$ is bounded in $L^{1}((0, T) \times \mathbb{R}^{d})$. This, together with (37), allows us to show that for a compact set $K \subset \mathbb{R}^{d}$ and $k > 0$, we have

$$ \|\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t+k, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t, x)\|_{L^{1}((0,T) \times K)} \le C(K,T)k \|\nabla u_{\varepsilon}^{i} * \rho_{\varepsilon}^{i}\|_{L^{2}((0,T) \times K)} \le C(K,T)\sqrt{k}. \quad (39) $$

Indeed considering a mollifier sequence $\omega_{\eta}(x) = \frac{1}{\eta^{d}}\omega(\frac{1}{\eta})$, we have

$$ \begin{aligned} \|\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t+k, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t, x)\|_{L^{1}((0,T) \times K)} &\le \|\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t+k, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t+k, x)\|_{L^{1}((0,T) \times K)} \\ &\quad + \|\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t, x)\|_{L^{1}((0,T) \times K)} \\ &\quad + \|\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t+k, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t, x)\|_{L^{1}((0,T) \times K)}, \end{aligned} $$

and

$$ \begin{aligned} & \|\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t+k, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t+k, x)\|_{L^{1}((0,T) \times K)} \\ &= \int_{0}^{T} \int_{K} \int_{K} (\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t+k, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t+k, x-y)) \omega_{\eta}(y) dy dx dt \\ &\le C(K,T)\eta \int_{0}^{T} \int_{K} \int_{K} \nabla \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t+k, x) \omega_{\eta}(y) dy dx dt \\ &\le \tilde{C}(K,T)\eta \|\nabla u_{\varepsilon} * \rho_{\varepsilon}\|_{L^{2}((0,T) \times K)}. \end{aligned} $$
Similarly, we have

$$ \| \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}(t, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t, x) \|_{L^{1}((0,T) \times K)} \leq \tilde{C}(K, T) \eta \| \nabla u_{\varepsilon} * \rho_{\varepsilon} \|_{L^{2}((0,T) \times K)} $$

Finally, we obtain

$$ \begin{aligned} \| \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t+k, x) - \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(t, x) \|_{L^{1}((0,T) \times K)} &= \| \int_{t}^{t+k} \partial_{t} \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} * \omega_{\eta}(s, x) ds \|_{L^{1}((0,T) \times K)} \\ &\leq \| \int_{t}^{t+k} \mathrm{div}(Q_{\varepsilon}^{i}) * \omega_{\eta}(s, x) ds \|_{L^{1}((0,T) \times K)} \\ &\leq \int_{t}^{t+k} \| Q_{\varepsilon}^{i} * \mathrm{div}(\omega_{\eta})(s, x) \|_{L^{1}((0,T) \times K)} ds \\ &\leq C(K) \frac{k}{\eta} \| Q_{\varepsilon}^{i} \|_{L^{1}((0,T) \times K)}. \end{aligned} $$

With the choice $\eta = \sqrt{k}$ we have the result.

Then, given (37) and (39), for all $i \in [1, N]$, $\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i}$ satisfies the assumptions of the Weil-Kolmogorov-Frechet theorem on $L^{1}((0, T) \times K)$, hence the sequence $(\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i})$ is compact in this space.

### 3.3 Convergence

Given the previous estimates, we can extract a subsequence (we do not relabel it for the sake of simplicity) such that

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \to u_{0}^{i} \quad \text{strongly in} \quad L^{1}(0, T; L_{\text{loc}}^{1}(\mathbb{R}^{d})), \qquad (40) $$

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \to u_{0}^{i} \quad \text{weakly in} \quad L^{2}(0, T; H^{1}(\mathbb{R}^{d})), \qquad (41) $$

$$ u_{\varepsilon}^{i} \to u_{0}^{i} \quad \text{weakly in} \quad L_{\text{loc}}^{1}((0, T) \times \mathbb{R}^{d}). \qquad (42) $$

In addition, from the uniform bound $\rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \in L^{2}(0, T; L^{\frac{2d}{d-2}}(\mathbb{R}^{d}))$ thanks to estimates (35), and the strong limit (40), we have

$$ \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \to u_{0}^{i} \quad \text{strongly in} \quad L^{p}(0, T; L_{\text{loc}}^{q}(\mathbb{R}^{d})) \quad \forall 1 \le p < 2, 1 \le q < \frac{2d}{d-2}. \qquad (43) $$

We are left to show the convergence in Eq. (21). Let us take $\phi^{i} \in C_{c}^{\infty}([0, T] \times \mathbb{R}^{d})$ for $i = 1, ..., N$. Then multiplying Eq. (21) by $\phi^{i}$ and integrating by parts, we have

$$ \begin{aligned} \int_{0}^{T} \int_{\mathbb{R}^{d}} \partial_{t} u_{\varepsilon}^{i} \phi^{i} &= \int_{0}^{T} \int_{\mathbb{R}^{d}} \mathrm{div}[u_{\varepsilon}^{i} \sum_{j=1}^{N} \nabla_{x} K_{\varepsilon}^{ij} * u_{\varepsilon}^{j}] \phi^{i} \\ &= - \int_{0}^{T} \int_{\mathbb{R}^{d}} u_{\varepsilon}^{i} (\sum_{j=1}^{N} \nabla_{x} K_{\varepsilon}^{ij} * u_{\varepsilon}^{j}) \Phi^{i} \end{aligned} $$

where we have used the notation $\Phi^{i} = \nabla \phi^{i}$. For the term on the right-hand side, the weak limit (42) shows that

$$ \int_{0}^{T} \int_{\mathbb{R}^{d}} u_{\varepsilon}^{i} \partial_{t} \phi^{i} \to \int_{0}^{T} \int_{\mathbb{R}^{d}} u_{0}^{i} \partial_{t} \phi^{i}. $$
The term on the right-hand side can be rewritten as

$$
\begin{align*}
\int_{\mathbb{R}^d} u_{\varepsilon}^i \nabla \left( \sum_{j=1}^N K_{\varepsilon}^{ij} * u_{\varepsilon}^j \right) \Phi^i &= \sum_{k=1}^p \sum_{j=1}^N \alpha^{ki} \alpha^{kj} \int_{\mathbb{R}^d} \Phi^i(t, x) u_{\varepsilon}^i \cdot \nabla (\tilde{\rho}_{\varepsilon}^i * \rho_{\varepsilon}^j * u_{\varepsilon}^j) \\
&= \sum_{k=1}^p \sum_{j=1}^N \alpha^{ki} \alpha^{kj} \int_{\mathbb{R}^d} \rho_{\varepsilon}^i * [\Phi^i(t, x) u_{\varepsilon}^i] \cdot \nabla (\rho_{\varepsilon}^j * u_{\varepsilon}^j) \\
&= \sum_{k=1}^p \sum_{j=1}^N \alpha^{ki} \alpha^{kj} \int_{\mathbb{R}^d} [\Phi^i(t, x) \rho_{\varepsilon}^i * u_{\varepsilon}^i + r_{\varepsilon}^i] \cdot \nabla (\rho_{\varepsilon}^j * u_{\varepsilon}^j),
\end{align*}
$$

with

$$
r_{\varepsilon}^{i}(t, x) = \int_{\mathbb{R}^{d}} [\Phi^{i}(t, y) - \Phi^{i}(t, x)] u_{\varepsilon}^{i}(t, y) \rho_{\varepsilon}^{i}(x - y) dy.
$$

Then, we obtain a decomposition in two terms

$$
\int_{\mathbb{R}^d} u_{\varepsilon}^i \nabla \left( \sum_{j=1}^N K_{\varepsilon}^{ij} * u_{\varepsilon}^j \right) \Phi^i = \sum_{k=1}^p \sum_{j=1}^N \alpha^{ki} \alpha^{kj} \int_{\mathbb{R}^d} \rho_{\varepsilon}^i * u_{\varepsilon}^i \Phi^i(t, x) \cdot \nabla (\rho_{\varepsilon}^j * u_{\varepsilon}^j) + \sum_{k=1}^p \sum_{j=1}^N \alpha^{ki} \alpha^{kj} \int_{\mathbb{R}^d} r_{\varepsilon}^i \cdot \nabla (\rho_{\varepsilon}^j * u_{\varepsilon}^j).
$$

The first term passes to the limit because of strong-weak limits. Indeed for $j, k \in [1, N]$, we have successively

$$
\alpha^{ki} \rho_{\varepsilon}^{i} * u_{\varepsilon}^{i} \to \alpha^{ki} u_{0}^{i} \quad \text{strongly in } L^{p}(0, T; L_{\text{loc}}^{q}(\mathbb{R}^{d})) \quad \forall 1 \le p < 2, 1 \le q < \frac{2d}{d-2}, \quad (44)
$$

$$
\alpha^{ki} \nabla \rho_{\varepsilon}^{j} * u_{\varepsilon}^{j} \to \alpha^{kj} \nabla u_{0}^{j} \quad \text{weakly in } L^{2}((0, T) \times \mathbb{R}^{d}), \quad (45)
$$

and thus

$$
\sum_{k=1}^{p} \sum_{j=1}^{N} \alpha^{ki} \alpha^{kj} \int_{\mathbb{R}^d} \rho_{\varepsilon}^i * u_{\varepsilon}^i \Phi^i(t, x) \cdot \nabla (\rho_{\varepsilon}^j * u_{\varepsilon}^j) \\
\to \sum_{k=1}^{p} \sum_{j=1}^{N} \alpha^{ki} \alpha^{kj} \int_{\mathbb{R}^d} u_0^i \Phi^i(t, x) \cdot \nabla u_0^j = \int_{\mathbb{R}^d} u_0^i \Phi^i(t, x) \cdot \sum_{j=1}^{N} \gamma^{ij} \nabla u_0^j
$$

For the second term, since $\nabla(\rho_{\varepsilon}^j * u_{\varepsilon}^j) \in L^2((0, T) \times \mathbb{R}^d)$, we need to prove that $r_{\varepsilon}^{ki}$ converges to 0 in $L^2((0, T) \times \mathbb{R}^d)$. With $L_{\phi}^i$ the Lipschitz constant for the function $\Phi^i$, we have

$$
\begin{align*}
r_{\varepsilon}^{i}(t, x) &\le L_{\phi}^{i} \int_{\mathbb{R}^d} |x-y| u_{\varepsilon}^{i}(t, y) |\rho_{\varepsilon}^{i}(x-y)| dy \\
&\le L_{\phi}^{i} \int_{\mathbb{R}^d} |x-y| u_{\varepsilon}^{i}(t, y) \rho_{\varepsilon}^{i}(x-y) dy \\
&\le L_{\phi}^{i} \int_{\mathbb{R}^d} |x-y| \left( u_{\varepsilon}^{i}(t, y) |\rho_{\varepsilon}^{i}(x-y)| \right)^{\frac{d}{d+2} + \frac{2}{d+2}} dy \\
&\le L_{\phi}^{i} \left[ \int_{\mathbb{R}^d} u_{\varepsilon}^{i}(t, y) \rho_{\varepsilon}^{i}(x-y) dy \right]^{\frac{d}{d+2}} \left[ \int_{\mathbb{R}^d} |x-y|^{\frac{d+2}{2}} u_{\varepsilon}^{i}(t, y) \rho_{\varepsilon}^{i}(x-y) dy \right]^{\frac{2}{d+2}} \\
&\le L_{\phi}^{i} \left[ u_{\varepsilon}^{i}(t) * \rho_{\varepsilon}^{i} \right]^{\frac{d}{d+2}} \left[ \int_{\mathbb{R}^d} |x-y|^{\frac{d+2}{2}} u_{\varepsilon}^{i}(t, y) \rho_{\varepsilon}^{i}(x-y) dy \right]^{\frac{2}{d+2}}.
\end{align*}
$$
Next, we use the Hölder inequality with $p = \frac{d+2}{d-2}$, $p' = \frac{d+2}{4}$ (notice however that the case $d=2$ is slightly different and left to the reader, then we choose $p < \infty$ and $p' > 1$, but close to, and the same arguments work) and we write

$$
\begin{aligned}
\|r_{\varepsilon}^{i}\|_{L^{2}(\mathbb{R}^{d})}^{2} & \leq L_{\phi}^{i 2} \int_{\mathbb{R}^{d}} \left[ u_{\varepsilon}^{i}(t, \cdot) * \rho_{\varepsilon}^{i} \right]^{\frac{2d}{d+2}} \left[ \int_{\mathbb{R}^{d}} |x-y|^{\frac{d+2}{2}} u_{\varepsilon}^{i}(t, y) \rho_{\varepsilon}^{i}(x-y) dy \right]^{\frac{4}{d+2}} dx \\
& \leq L_{\phi}^{i 2} \|u_{\varepsilon}^{i}(t, \cdot) * \rho_{\varepsilon}^{i}\|_{L^{\frac{2d}{d-2}}(\mathbb{R}^{d})}^{\frac{2d}{d+2}} \left[ \int_{\mathbb{R}^{2d}} |x-y|^{\frac{d+2}{2}} u_{\varepsilon}^{i}(t, y) \rho_{\varepsilon}^{i}(x-y) dy dx \right]^{\frac{4}{d+2}} \\
& \leq L_{\phi}^{i 2} \|u_{\varepsilon}^{i}(t, \cdot) * \rho_{\varepsilon}^{i}\|_{L^{\frac{2d}{d-2}}(\mathbb{R}^{d})}^{\frac{2d}{d+2}} \|u_{\varepsilon}^{i}(t, y)\|_{L^{1}(\mathbb{R}^{d})}^{\frac{4}{d+2}} \varepsilon^{2} \left[ \int_{\mathbb{R}^{d}} |z|^{\frac{d+2}{2}} \rho^{i}(z) dz \right]^{\frac{4}{d+2}},
\end{aligned}
$$

where we have used the change of variable $x \to z = \frac{x-y}{\varepsilon}$. Since we have proved that $u_{\varepsilon}^{i}(t) * \rho_{\varepsilon}$ is uniformly bounded in $L^{2}(0, T; L^{\frac{2d}{d-2}}(\mathbb{R}^{d}))$, we conclude that

$$
\begin{aligned}
\|r_{\varepsilon}^{i}\|_{L^{2}(0,T;L^{2}(\mathbb{R}^{d}))}^{2} & \leq \varepsilon^{2} L_{\phi}^{i 2} \|u_{\varepsilon}^{i}\|_{L^{\infty}(0,T;L^{1}(\mathbb{R}^{d}))}^{\frac{4}{d+2}} \|z^{\frac{d+2}{2}} \rho\|_{L^{1}(\mathbb{R}^{d})}^{\frac{4}{d+2}} \int_{0}^{T} \|u_{\varepsilon}^{i}(t,y) * \rho_{\varepsilon}^{i}\|_{L^{\frac{2d}{d-2}}(\mathbb{R}^{d})}^{\frac{2d}{d+2}} dt \\
& \leq \varepsilon^{2} L_{\phi}^{i 2} \|u_{\varepsilon}^{i}\|_{L^{\infty}(0,T;L^{1}(\mathbb{R}^{d}))}^{\frac{4}{d+2}} \|z^{\frac{d+2}{2}} \rho^{i}\|_{L^{1}(\mathbb{R}^{d})}^{\frac{4}{d+2}} T^{\frac{2}{d+2}} \left( \int_{0}^{T} \|u_{\varepsilon}^{i}(t,y) * \rho_{\varepsilon}^{i}\|_{L^{\frac{2d}{d-2}}(\mathbb{R}^{d})}^{2} dt \right)^{\frac{d}{d+2}}
\end{aligned}
$$

Thus

$$
\begin{aligned}
& \sum_{k=1}^{p} \sum_{j=1}^{N} \alpha^{ki} \alpha^{kj} \int_{0}^{T} \int_{\mathbb{R}^{d}} r_{\varepsilon}^{ki}(t, x) \nabla(\rho_{\varepsilon}^{kj} * u_{\varepsilon}^{j}) \\
& \leq \left( \sup_{(i,j) \in [1,N]^2} |\alpha^{ki}| \right)^{2} \sum_{k=1}^{p} \sum_{j=1}^{N} \|r_{\varepsilon}^{ki}\|_{L^2(0,T;L^2(\mathbb{R}^d))}^2 \|\nabla(\rho_{\varepsilon}^{kj} * u_{\varepsilon}^{j})\|_{L^2(0,T;L^2(\mathbb{R}^d))}^2, \to 0.
\end{aligned}
$$

Therefore, we have proved that at the limit the equation holds

$$
\int_{0}^{T} \int_{\mathbb{R}^{d}} u_{0}^{i} \partial_{t} \phi^{i} - \sum_{j,k=1}^{N} \alpha^{ki} \alpha^{kj} \int_{0}^{T} \int_{\mathbb{R}^{d}} u_{0}^{i} \nabla(\sum_{j=1}^{N} u_{0}^{j}) \nabla \phi^{i} = 0,
$$

and thus $u_0^i$ is solution of (22).

# 4 Numerical simulations

We illustrate the localisation limit in 2D thanks to numerical simulations. We consider a spatial square domain $\Omega = [x_{\min}, x_{\max}] \times [y_{\min}, y_{\max}] \in \mathbb{R}^2$ taken large enough so that the solution (spatial) support stays far from the boundaries. The spatial domain is discretized into $N_x \times N_x$ regularly spaced points with space step $\Delta 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 $\rho^i(x)$ for $i = 1...4$ are chosen as Gaussian functions of variance $\sigma_i^2$:

$$
\rho^i(x) = \frac{1}{2\pi\sigma_i^2} \exp\left(-\frac{|x|^2}{2\sigma_i^2}\right).
$$
With this choice, $K^{ij}(x)$ is again a Gaussian function of variance $\sigma_i^2 + \sigma_j^2$. Therefore, $K^{ij}$ are of the form:

$$K^{ij}(x) = \frac{(A^T A)_{ij}}{2\pi(\sigma_i^2 + \sigma_j^2)} \exp^{-\frac{|x|^2}{2(\sigma_i^2 + \sigma_j^2)}}.$$

Note that thus defined, the coefficients $\alpha^{ij}$ of the matrix A control the intensities of the interactions, while the variances $\sigma_i^2$ 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., $\sigma_i = \sigma_j \quad \forall (i, j) = 1...4$. These conditions satisfy hypothesis (25).

* case 2: Species-dependent repulsion intensities: $(A^T A)_{ij} = \min(i, j)$ (the matrix A is an upper triangular matrix) and species-dependent interaction distances, i.e., $\sigma_i \neq \sigma_j \forall i \neq j$. 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., $\sigma_i \neq \sigma_j, \forall i \neq 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:

$$\rho_i(x) = \frac{1}{\pi S^2} \mathbb{1}_{B(0,S)}(x), \quad \forall 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 $\varepsilon \to 0$ in the different cases for the choice of the interaction kernel.

## 4.1 Case 1: same interaction distances for all species $\sigma_i^2 = \sigma_j^2 = 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 $\varepsilon$ ($\varepsilon = 5$, first column, $\varepsilon = 1$, second column and $\varepsilon = 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 $\varepsilon = 5$ (first column), but seems to get closer as $\varepsilon$ decreases. For large $\varepsilon$, 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 $\varepsilon$ 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(\lambda)$, which gives the average density on rings of size $d\lambda$:

$$g_i^\varepsilon(\lambda, t)d\lambda = \int_{B(0,\lambda+d\lambda)\setminus B(0,\lambda)} u_i^\varepsilon(x,t)dx,$$
<figure><img src="image_3.png" /><figcaption>Figure 1: Case 1: Solutions of the nonlocal model (21) for different values of $\epsilon$ ($\epsilon = 5$, first column, $\epsilon = 1$, second column and $\epsilon = 0.5$, third column) and of the local model (22) (fourth column) at time $t = 100$, for each of the 4 species (different lines).</figcaption></figure>

and we denote by $g_i^0(x, t)$ the quantity obtained in the same way with the solution of the local model $u_i^0(x, t)$. In the left subplots of Fig. 2, we show the quantities $g_i^0(\lambda, 100)$ (plain lines) $g_i^\epsilon(\lambda, 100)$ (dotted lines) for different values of $\epsilon$ (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 $\epsilon$ (red curves), while the radial distributions of the nonlocal model get closer to the one of the local model when decreasing $\epsilon$ (compare the green and blue lines). In order to quantify the convergence of the non-local to the local model as $\epsilon \to 0$, we plot in Fig. 2 the discrete $L^2$ distance between $u^\epsilon$ and $u^0$ at time 100 as function of $\epsilon$ for each of the 4 species.

As one can see, the nonlocal model seems to converge to the local one as $\epsilon$ 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 $\epsilon$.

## 4.2 Case 2: different interaction distances as function of the species $\sigma_i \neq \sigma_j$, different interaction intensities $(A^T A)_{ij} = \min(i, j)$

We show in Fig. 3 the simulations obtained in Case 2 with $\sigma_1^2 = 0.1, \sigma_2^2 = 0.2, \sigma_3^2 = 0.3, \sigma_4^2 = 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 $\epsilon \to 0$ of the nonlocal one is the same here as for Case 1.

As one can observe in Fig. 3, for $\epsilon = 5$, a quite different behavior is obtained when the interaction
<figure><img src="image_4.png" /><figcaption>Figure 2: Case 1: Solutions of the nonlocal model (21) for different values of $\epsilon$ ($\epsilon = 5$, first column, $\epsilon = 1$, second column and $\epsilon = 0.5$, third column) and of the local model (22) (fourth column) at time $t = 100$, for each of the 4 species (different lines).</figcaption></figure>

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 $\epsilon = 5$ indeed, species with larger indices do not diffuse faster than species with lower indices as observed previously, although the coefficients $\gamma^{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 $\rho^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 $\epsilon$). For small values of $\epsilon = 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 $\epsilon \to 0$ than the other species. This is because species 4 is the one the most impacted by the interaction distance controlled by $\epsilon$.

## 4.3 Case 3: different interaction distances as function of the species $\sigma_i \neq \sigma_j$, same interaction intensities $(A^T A)_{ij} = 1, \forall (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 \forall (i, j)$.

As before, we show in Fig. 5 the simulations obtained in Case 3 with $\sigma_1^2 = 0.1$, $\sigma_2^2 = 0.2$, $\sigma_3^2 =$
<figure><img src="image_5.png" /><figcaption>Figure 3: Case 2: Solutions of the nonlocal model (21) for different values of $\varepsilon$ ($\varepsilon = 5$, first column, $\varepsilon = 1$, second column and $\varepsilon = 0.5$, third column) and of the local model (22) (fourth column) at time $t = 100$, for each of the 4 species (different lines).</figcaption></figure>

$0.3, \sigma_4^2 = 0.4$ (increasing interaction distance with increasing species index). We adopt the same representation as in Figs. 1 and 3.

In this case, the diffusion is even slower for species with large indices compared to Case 2 for $\varepsilon = 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 $\varepsilon$ 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 $\varepsilon \to 0$.

These numerical results suggest that the nonlocal model converges to the local one when $\varepsilon \to 0$ in a more general framework than the one considered in this paper, i.e., without the full rank condition for A.

# 5 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
<figure><img src="image_6.png" /><figcaption>Figure 4: Case 2: Solutions of the nonlocal model (21) for different values of ε (ε = 5, first column, ε = 1, second column and ε = 0.5, third column) and of the local model (22) (fourth column) at time t = 100, for each of the 4 species (different lines).</figcaption></figure>

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

$$ \partial_t u_0^i - \Delta[u_0^i(1 + \sum_{j=1}^N \gamma^{ij} u_0^j)] = 0, \quad t \ge 0, x \in \mathbb{R}^d, \quad i \in [1, N], $$

which is also an active field of research, [18, 31] for which the nonlocal to local convergence is also studied.

**Acknowledgments.** DP was supported by Sorbonne Alliance University with an Emergence project MATHREGEN, grant number S29-05Z101 and by Agence Nationale de la Recherche (ANR) under the project grant number ANR-22-CE45-0024-01. The authors warmly thank Markus Schmidtchen for his careful reading and suggestions for the bibliography.
<figure><img src="image_7.png" /><figcaption>Figure 5: Case 3: Solutions of the nonlocal model (21) for different values of $\varepsilon$ ($\varepsilon = 5$, first column, $\varepsilon = 1$, second column and $\varepsilon = 0.5$, third column) and of the local model (22) (fourth column) at time $t = 100$, for each of the 4 species (different lines).</figcaption></figure>

References

[1] Luca Alasio, Maria Bruna, Simone Fagioli, and Simon Schulz. Existence and regularity for a system of porous medium equations with small cross-diffusion and nonlocal drifts. *Nonlinear Analysis*, 223:113064, 2022.

[2] Julien Barré, José Carrillo, Pierre Degond, Diane Peurichard, and Ewelina Zatorska. Particle interactions mediated by dynamical networks: Assessment of macroscopic descriptions. *J Nonlinear Sci*, 28:235–268, 2018.

[3] Julien Barré, Pierre Degond, and Ewelina Zatorska. Kinetic theory of particle interactions mediated by dynamical networks. *Multiscale Modeling & Simulation*, 15(3):1294–1323, 2017.

[4] Julien Barré, Pierre Degond, Diane Peurichard, and Ewelina Zatorska. Modelling pattern formation through differential repulsion. *Networks and Heterogeneous Media*, 15(3):307–352, 2020.

[5] Andrea L Bertozzi, José A Carrillo, and Thomas Laurent. Blow-up in multidimensional aggregation equations with mildly singular interaction kernels. *Nonlinearity*, 22(3):683, 2009.

[6] Andrea L. Bertozzi and Thomas Laurent. Finite-Time Blow-up of Solutions of an Aggregation Equation in $R^n$. *Communications in Mathematical Physics*, 274(3):717–735, 2007.
<figure><img src="image_8.png" /><figcaption>Figure 6: Case 3: Solutions of the nonlocal model (21) for different values of $\epsilon$ ($\epsilon = 5$, first column, $\epsilon = 1$, second column and $\epsilon = 0.5$, third column) and of the local model (22) (fourth column) at time $t = 100$, for each of the 4 species (different lines).</figcaption></figure>

[7] Andrea L. Bertozzi, Thomas Laurent, and Jesús Rosado. $L^p$ theory for the multidimensional aggregation equation. *Communications on Pure and Applied Mathematics*, 64(1):45–83, 2023/04/20 2011.

[8] Maria Bruna, S Jonathan Chapman, and Martin Robinson. Diffusion of particles with short-range interactions. *SIAM Journal on Applied Mathematics*, 77(6):2294–2316, 2017.

[9] Martin Burger and Antonio Esposito. Porous medium equation and cross-diffusion systems as limit of nonlocal interaction, 2022. arXiv 2202.05030.

[10] Martin Burger and Marco Di Francesco. Large time behavior of nonlocal aggregation models with nonlinear diffusion. *Networks and Heterogeneous Media*, 3(4):749–785, 2008.

[11] José A. Carrillo, Charles Elbar, and Jakub Skrzeczkowski. Degenerate Cahn-Hilliard systems: From nonlocal to local. *hal-04043503*, 23pp, March 2023.

[12] José A Carrillo, Stephan Martin, and Vladislav Panferov. A new interaction potential for swarming models. *Physica D: Nonlinear Phenomena*, 260:112–126, 2013.

[13] José Antonio Carrillo, Katy Craig, and Francesco S Patacchini. A blob method for diffusion. *Calculus of Variations and Partial Differential Equations*, 58:1–53, 2019.
[14] José Antonio Carrillo, Marco Di Francesco, Antonio Esposito, Simone Fagioli, and Markus Schmidtchen. Measure solutions to a system of continuity equations driven by newtonian nonlocal interactions. *Discrete and Continuous Dynamical Systems*, 40(2):1191–1231, 2020.

[15] José Antonio Carrillo, Antonio Esposito, and Jeremy Sheung-Him Wu. Nonlocal approximation of nonlinear diffusion equations. *arXiv preprint arXiv:2302.08248*, 2023.

[16] José Carrillo, Massimo Diffrancesco, Alessio Figalli, T. Laurent, and D. Slepčev. Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations. *Duke Mathematical Journal*, 156, 02 2011.

[17] José A. Carrillo, Alina Chertock, and Yanghong Huang. A finite-volume method for nonlinear nonlocal equations with a gradient flow structure. *Communications in Computational Physics*, 17(1):233–258, 2015.

[18] Li Chen and Ansgar Jüngel. Analysis of a multidimensional parabolic population model with strong cross-diffusion. *SIAM J. Math. Anal.*, 36(1):301–322, 2004.

[19] Xiuqing Chen, Esther S Daus, and Ansgar Jüngel. Global existence analysis of cross-diffusion population systems for multiple species. *Archive for Rational Mechanics and Analysis*, 227(2):715–747, 2018.

[20] Rinaldo M. Colombo and Magali Lécureux-Mercier. Nonlocal crowd dynamics models for several populations. *Acta Mathematica Scientia*, 32(1):177–196, 2012.

[21] Gianluca Crippa and Magali Lécureux-Mercier. Existence and uniqueness of measure solutions for a system of continuity equations with non-local flow. *Nonlinear Differential Equations and Applications NoDEA*, 20(3):523–537, 2013.

[22] Noemi David, Tomasz Dębiec, Mainak Mandal, and Markus Schmidtchen. A degenerate cross-diffusion system as the inviscid limit of a nonlocal tissue growth model. *arXiv preprint arXiv:2303.10620*, 2023.

[23] Marco Di Francesco, Antonio Esposito, and Simone Fagioli. Nonlinear degenerate cross-diffusion systems with nonlocal interaction. *Nonlinear Analysis*, 169:94–117, 2018.

[24] Marco Di Francesco and Simone Fagioli. Measure solutions for non-local interaction PDEs with two species. *Nonlinearity*, 26:2777, 09 2013.

[25] M. R. D'Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes. Self-propelled particles with soft-core interactions: Patterns, stability, and collapse. *Phys. Rev. Lett.*, 96:104302, Mar 2006.

[26] Charles Elbar and Jakub Skrzeczkowski. Degenerate Cahn-Hilliard equation: From nonlocal to local. *J. Differential Equations*, 364:576–611, 2023.

[27] Valeria Giunta, Thomas Hillen, Mark Lewis, and Jonathan R. Potts. Local and global existence for nonlocal multispecies advection-diffusion models. *SIAM Journal on Applied Dynamical Systems*, 21(3):1686–1708, 2022.

[28] Ansgar Jüngel, Stefan Portisch, and Antoine Zurek. Nonlocal cross-diffusion systems for multi-species populations and networks. *Nonlinear Analysis*, 219:112800, 2022.
[29] Elliott H. Lieb and Michael Loss. Analysis, volume 14 of Graduate Studies in Mathematics. American Mathematical Society, Providence, RI, second edition, 2001.

[30] Pierre-Louis Lions and Sylvie Mas-Gallic. Une méthode particulaire déterministe pour des équations diffusives non linéaires. C. R. Acad. Sci. Paris Sér. I Math., 332(4):369-376, 2001.

[31] Ayman Moussa. From nonlocal to classical Shigesada-Kawasaki-Teramoto systems: triangular case with bounded coefficients. SIAM J. Math. Anal., 52(1):42-64, 2020.

[32] Karl Oelschläger. Large systems of interacting particles and the porous medium equation. Journal of differential equations, 88(2):294-346, 1990.

[33] Robert Philipowski. Interacting diffusions approximating the porous medium equation and propagation of chaos. Stochastic Processes and their Applications, 117(4):526-538, 2007.

[34] Jonathan R. Potts and Mark A. Lewis. Spatial memory and taxis-driven pattern formation in model ecosystems. Bulletin of Mathematical Biology, 81(7):2725-2747, 2019.

[35] Pawel Romanczuk, Markus Bär, Werner Ebeling, Benjamin Lindner, and Lutz Schimansky-Geier. Active brownian particles: From individual to collective stochastic dynamics. The European Physical Journal Special Topics, 202:1-162, 2012.