Multiple patterns formation for an aggregation/diffusion predator-prey system

1.   Introduction

The mathematical modelling of the collective motion through aggregation/diffusion phenomena has raised a lot of interest in the recent years and it has been deeply studied for its application in several areas, such as biology [9, 39, 47, 48], ecology [35, 42, 43], animal swarming [3, 4, 40, 46] sociology and economics, [49, 50, 51]. One of the common idea in this modelling approach is that a certain population composed by agents evolves according to long-range attraction and short-range repulsion forces between agents. We are interested in modelling the problem of predator-prey interactions, namely we consider two populations that attract (prey) and repel (predators) each others. The pioneering works for this problem are the ones by Lotka, [36] and Volterra[54], which describe the predator-prey interaction via reaction terms in a set of differential equations, possibly combined with diffusion terms, see [41] and the references therein.

As in [21], in this paper we model predator-prey interactions via a transport terms rather than a reaction ones as follows: consider $N$ predators located at $X_1,\ldots,X_N\in {\mathbb{R}^{n}}$, and $M$ prey at $Y_1,\ldots,Y_M\in {\mathbb{R}^{n}}$ with masses $m_X^i>0$ and $m_Y^i>0$ respectively. We assume that each agent of the same species interacts under the effect of a radial non-local force that is attractive in the long range and repulsive in the short range. Moreover, predators are attracted by the prey, while the latter are subject to a repulsive force from the predators, that is proportional to the previous one. This set of modelling assumptions leads to the following system of ODEs:

with $i = 1,\ldots,N$ and $j = 1,\ldots,M$. The potentials $S_1^a$ and $S_2^a$ are called self-interaction and model the long-range attraction among agents of the same species. The potential $K$ is responsible for the predator-prey interaction, and it is called cross-interaction potential. The coefficient $\alpha>0$ models the escape propensity of prey from the predators. The short-range repulsion among particles of the same species is modelled by the non-local forces $S_1^r$ and $S_2^r$, and we can assume that it scales with the number of particles, $S_i^r(z) = N^\beta S(N^{\beta/n}z)$ for a smooth functional $S$, see [40].

The formal limit when the number of particles tends to infinity leads to the following system of partial differential equations

where $\rho$ and $\eta$ are the densities of predators and prey respectively. Through this limit the (non-local) short-range repulsion formally turns to a (local) nonlinear diffusion terms, being $d_1$ and $d_2$ positive constants modelling the spreading velocity, while the long-range attraction takes into account the non-local self-interactions. We can therefore lighten the notation by calling $S_1^a = S_\rho$ and $S_2^a = S_\eta$.

The goal of this paper is to show that the model above catches one of the main features that occur in nature, namely the formation of spatial patterns where the predators are surrounded of empty zones and the prey aggregates around, that is usually observed in fish schools or in flock of sheeps, see [32, 38]. In the fully aggregative case, namely system (2) with $d_1 = d_2 = 0$, the formation of these types of patterns has been studied in several papers, see [15, 27, 21, 46] and references therein.

Existence theory for solutions to system of the form (2) can be performed in the spirit of [10, 20]. In particular, system (2) should be framed in the context of non symmetrizable systems, for which the Wasserstein gradient flow theory developed in [1] and adapted to systems in [22] does not work. In [10, 20, 22], the authors consider different choices for the diffusion part (no diffusion in [22], diagonal nonlinear diffusion in [10] and cross-diffusion with dominant diagonal in [20]), and the existence of solutions is proved via an implicit-explicit version of the JKO scheme [33].

In the following, we reduce our analysis to the one-dimensional setting

We are interested in the existence of stationary solutions to (3), which are solutions to the following system

as well as their properties, e.g., symmetry, compact support, etc. We shall assume throughout the paper that $d_1 = d_2 = \theta$ for simplicity. Note that we can always assume this by a simple scaling argument. Indeed, it is enough to multiply the first equation in (4) by $d_{2}/d_1$, and setting $d_2 = \theta$, $\widetilde{S}_\rho = \frac{d_2}{d_1}S_\rho$, $\widetilde{K} = \frac{d_2}{d_1}K$ and $\widetilde{\alpha} = \frac{d_1}{d_2}\alpha$ to get

The stationary equation for the one species case, i.e.,

is studied several papers, see [2, 7, 11, 16] and therein references. In [7] the Krein-Rutman theorem is used in order to characterise the steady states as eigenvectors of a certain non-local operator. The authors prove that a unique steady state with given mass and centre of mass exists provided that $\theta < \|K\|_{L^1}$, and it exhibits a shape similar to a Barenblatt profile for the porous medium equation; see [52] and [24] for the local stability analysis. Similar techniques are used in [8] in order to partly extend the result to more general nonlinear diffusion, see also [34]. This approach is used in [6] in order to explore the formation of segregated stationary states for a system similar to (3) but in presence of cross-diffusion. Unfortunately, when dealing with systems, it is not possible to reproduce one of the major issues solved in [7], namely the one-to-one correspondence between the diffusion constant (eigenvalue) and the support of the steady state. A support-independent existence result for small diffusion coefficient $\theta$ is obtained in [6] by using the generalised version of the implicit function theorem, see also [5] where this approach is used in the one species case.

In this paper we apply the aforementioned approach in order to show that stationary solutions to (3) are composed of multiple Barenblatt profiles. Let us introduce, for fixed $z_\rho, z_\eta >0$, the following space

Definition 1.1. We say that a pair $(\rho,\eta)\in \mathcal{M}$ is a multiple bumps steady state to (3) if the pair $(\rho,\eta)$ solves (4) weakly and there exist two numbers $N_\rho,N_\eta\in \mathbb{N}$, and two families of intervals $I_\rho^i = \left[l_\rho^i,r_\rho^i\right]$, for $i = 1,...,N_\rho$, and $I_\eta^h = \left[l_\eta^h,r_\eta^h\right]$, for $h = 1,...,N_\eta$ such that

● $I_\rho^i\cap I_\rho^j = \emptyset$, for $i,j = 1,...,N_\rho$, $i\neq j$ and $I_\eta^h\cap I_\eta^k = \emptyset$, for $h,k = 1,...,N_\eta$, $h\neq k$,

● $\rho$ and $\eta$ are supported on

respectively and

where, for $i = 1,...,N_\rho$ and $h = 1,...,N_\eta$, $\rho^i$ and $\eta^h$ are even w.r.t the centres of $I_\rho^i$ and $I_\eta^h$ respectively, non-negative and $C^1$ functions supported on that intervals.

Example 1.1. A possible example of steady states as defined above is a steady state $(\rho,\eta)$ consisting of three bumps for each one of the pair $(\rho,\eta)$ ($N_\rho = N_\eta = 3$), namely, $\rho^1,\rho^2,\rho^3$ and $\eta^1,\eta^2,\eta^3$ respectively, with centers of masses $\{cm_\rho^i: = r_\rho^i-l_\rho^i\}_{i = 1}^{N_\rho}$ and $\{cm_\eta^i: = r_\eta^i-l_\eta^i\}_{i = 1}^{N_\eta}$ as solutions to system (17). A plot is carried out in Figure 1 which shows all the properties stated in Definition 1.1.

Remark 1.1. We remark that one should be careful with finding a shape of steady state. More precisely, one should choose the numbers $N_\rho,N_\eta$ and the centers of masses $\{cm_\rho^i\}_{i = 1}^{N_\rho}$ and $\{cm_\eta^i\}_{i = 1}^{N_\eta}$ such that all the conditions required for the existence of steady states are satisfied, see Theorem 1.1.

In order to simplify the notations, in the following we will denote with $l\in\{\rho,\eta\}$ a generic index that recognise one of the two families. Throughout the paper we shall assume that the kernels satisfy the following:

(A1) $S_\rho$, $S_\eta$ and $K$ are $C^2( \mathbb{R})$ functions.

(A2) $S_\rho$, $S_\eta$ and $K$ are radially symmetric and decreasing w.r.t. the radial variable.

(A3) $S_\rho$, $S\eta$ and $K$ are non-negative, with finite $L^1$-norm supported on the whole real line $\mathbb{R}$.

Note that assumption (A2) together with the sign in front of the nonlocal terms $S_\rho$ and $K$ (in the first equation) and $S_\eta$ in system (3), give the effect of an attractive potential, i.e. a radial interaction potential $G(x) = g(|x|)$ for some $g:[0,+\infty)\rightarrow \mathbb{R}$, such that $g'(r)>0$ for $r>0$. For $K$ in the second equation we obtain the effect of a repulsive potential, i.e. $g'(r)<0$ for $r>0$.

The main result of the paper is the following

Theorem 1.1. Assume that $S_\rho$, $S_\eta$ and $K$ are interaction kernels are under the assumptions (A1), (A2) and (A3). Consider $N_\rho,N_\eta\in\mathbb{N}$ and let $z_l^i$ be fixed positive numbers for $i = 1,2,\cdots,N_l,$ and $l \in \{\rho,\eta\}$. Consider two families of real numbers $\{cm_\rho^i\}_{i = 1}^{N_\rho}$ and $\{cm_\eta^i\}_{i = 1}^{N_\eta}$ such that

(i) $\{cm_\rho^i\}_{i = 1}^{N_\rho}$ and $\{cm_\eta^i\}_{i = 1}^{N_\eta}$ are stationary solutions of the purely non-local particle system, that is, for $i = 1,2,\cdots,N_l,$ for $l,h \in \{\rho,\eta\}$ and $l\neq h$,

(ii) the following quantities

are strictly positive, for all $i = 1,2,\cdots,N_l$, $l,h \in \{\rho,\eta\}$ and $l\neq h$.

where $\alpha_{\rho} = 1$ and $\alpha_{\eta} = -\alpha$. Then, there exists a constant $\theta_0$ such that for all $\theta\in (0, \theta_0)$ the stationary equation (4) admits a unique solution in the sense of Definition 1.1 of the form

where

● each interval $I_l^i$ is symmetric around $cm_l^i$ for all $i = 1,2,\cdots,N_l$, $l \in \{\rho,\eta\}$,

● $\rho^i$ and $\eta^j$ are $C^1$, non-negative and even w.r.t the centres of $I_\rho^i$ and $I_\eta^j$ respectively, with masses $z_\rho^i$ and $z_\eta^j$, for $i = 1,...,N_\rho$ and $j = 1,...,N_\eta$,

● the solutions $\rho$ and $\eta$ have fixed masses

respectively.

The paper is structured as follows. In Section 2 we recall the basics notions on optimal transport and we introduce the $p$-Wasserstein distances in spaces of probability measures. Then, we recall the strategy for proving existence to systems of the form (2). The remaining part of the Section is devoted to a preliminary and partial existence analysis of steady states via the Krein-Rutman theorem of a particular type of stationary solutions that we call mixed steady state. Section 3 is devoted to the proof of Theorem 1.1 in which existence and uniqueness results for multiple bumps stationary solutions are proved in case of small diffusion coefficient using the Implicit Function Theorem. We conclude the paper with Section 4, complementing our results with numerical simulations that also show interesting stability issues of the stationary states, namely transitions between states and others effects such as traveling waves profiles.

2.   Preliminary results

2.1.   Tools in optimal transport

We start collecting preliminary concepts on the Wasserstein distance. Let ${\mathscr{P}}( \mathbb{R}^n)$ be the space of probability measures on $\mathbb{R}^n$ and fix $p \in [1,+\infty)$. The space of probability measures with finite $p$-moment is defined by

For a measure $\mu\in {\mathscr{P}}( \mathbb{R}^n)$ and a Borel map $T: \mathbb{R}^n\rightarrow \mathbb{R}^k$, denote with $T_{\#}\mu\in {\mathscr{P}}( \mathbb{R}^n)$ the push-forward of $\mu$ through $T$, defined by

We endow the space ${\mathscr{P}}_{p}( \mathbb{R}^n)$ with the Wasserstein distance, see for instance [1, 45, 53]

where $\Gamma(\mu_{1},\mu_{2})$ is the class of transport plans between $\mu$ and $\nu$, that is the class of measures $\gamma \in {\mathscr{P}}( \mathbb{R}^n\times \mathbb{R}^n)$ such that, denoting by $\pi^{i}$ the projection operator on the $i$-th component of the product space, the marginality condition $\pi^{i}_{\#}\gamma = \mu_{i}$ $i = 1,2$ is satisfied.

Since we are working in a 'multi-species' structure, we consider the product space ${\mathscr{P}}_{p}( \mathbb{R}^n)\times {\mathscr{P}}_{p}( \mathbb{R}^n)$ endowed with a product structure. In the following we shall use bold symbols to denote elements in a product space. For a $p\in [1,+\infty]$, we use the notation

with $\mathit{\boldsymbol{\mu}} = (\mu_1,\mu_2),\mathit{\boldsymbol{\nu}} = (\nu_1,\nu_2) \in {\mathscr{P}}_{p}( \mathbb{R}^n)\times {\mathscr{P}}_{p}( \mathbb{R}^n)$. In the one-dimensional case, given $\mu\in {\mathscr{P}} ( \mathbb{R})$, we introduce the pseudo-inverse variable $u_\mu \in L^1([0,1]; \mathbb{R})$ as

see [14].

2.2.   Weak solutions for the time-dependent system

In the Section 1 we mentioned that the well-posedness of (3) can be stated according to the results in [10, 20] in an arbitrary space dimension $n$. In these papers, the existence of weak solutions is provided using an implicit-explicit version of the Jordan-Kinderlehrer-Otto (JKO) scheme [33, 22], that we will sketch in the following. A key point in this approach is to associate to (3) a relative energy functional

for a fixed reference couple of measures $(\mu,\nu)$. We state our definition of weak measure solution for (3), in the space ${\mathscr{P}}_2( {\mathbb{R}^{n}})^2: = {\mathscr{P}}_{p}( \mathbb{R}^n)\times {\mathscr{P}}_{p}( \mathbb{R}^n)$.

Definition 2.1. A curve $\mathit{\boldsymbol{\mu}} = (\rho(\cdot),\eta(\cdot)):[0,+\infty)\longrightarrow {\mathscr{P}_{2}({\mathbb{R}^{n}})}^2$ is a weak solution to (3) if

(i) $\rho,\,\eta\in L^{2}([0,T]\times \mathbb{R}^n)$ for all $T>0$, and $\nabla \rho,\,\nabla \eta \in L^2([0,+\infty)\times \mathbb{R}^n)$ for $i = 1,2$,

(ii) for almost every $t\in [0,+\infty)$ and for all $\phi,\varphi\in C_c^{\infty}( {\mathbb{R}^{n}})$, we have

Theorem 2.1. Assume that (A1)-(A3) are satisfied. Let $\mathit{\boldsymbol{\mu}}_0 = (\rho_{1,0},\rho_{2,0})\in {\mathscr{P}_{2}({\mathbb{R}^{n}})}^2$ such that

Then, there exists a weak solution to (3) in the sense of Definition 2.1 with $\mathit{\boldsymbol{\mu}}(0) = \mathit{\boldsymbol{\mu}}_0$.

As already mentioned the proof of Theorem 2.1 is a special case of the results in [10, 20] and consists in the following main steps:

1. Let $\tau>0$ be a fixed time step and consider the initial datum $\mathit{\boldsymbol{\mu}}_0\in {\mathscr{P}_{2}({\mathbb{R}^{n}})}^2$, such that ${\mathcal{F}}_{[\mathit{\boldsymbol{\mu}}_0]}\left(\mathit{\boldsymbol{\mu}}_0\right)<+\infty$. Define a sequence $\left\{\mathit{\boldsymbol{\mu}}_\tau^k\right\}_{k\in\mathbb{N}}$ recursively: $\mathit{\boldsymbol{\mu}}_\tau^0 = \mathit{\boldsymbol{\mu}}_0$ and, for a given $\mathit{\boldsymbol{\mu}}_\tau^k\in {\mathscr{P}_{2}({\mathbb{R}^{n}})}^2$ with $k\geq 0$, we choose $\mathit{\boldsymbol{\mu}}_\tau^{k+1}$ as follows:

For a fixed $T>0$, set $N: = \left[\frac{T}{\tau}\right]$, and

with $\mathit{\boldsymbol{\mu}}_\tau^k$ defined in (9).

2. There exists an absolutely continuous curve $\tilde{\mathit{\boldsymbol{\mu}}}: [0,T]\rightarrow {\mathscr{P}_{2}({\mathbb{R}^{n}})}^2$ such that the piecewise constant interpolation $\mathit{\boldsymbol{\mu}}_\tau$ admits a sub-sequence $\mathit{\boldsymbol{\mu}}_{\tau_h}$ narrowly converging to $\tilde{\mathit{\boldsymbol{\mu}}}$ uniformly in $t\in[0,T]$ as $h\rightarrow +\infty$. This is a standard result coming from the minimising condition (9).

3. There exist a constant $C>0$ such that

and the sequence $\mathit{\boldsymbol{\mu}}_{\tau_{h}}:[0,+\infty)\longrightarrow {\mathscr{P}_{2}({\mathbb{R}^{n}})}^2$ converges to $\tilde{\mathit{\boldsymbol{\mu}}}$ strongly in

The estimate in (10) can be deduced by using the so called Flow-interchange Lemma introduced in [37], see also [25]. In order to deduce the strong convergence we use the extended version of the Aubin-Lions Lemma in [44].

4. The approximating sequence $\mathit{\boldsymbol{\mu}}_{\tau_h}$ converges to a weak solution $\tilde{\mathit{\boldsymbol{\mu}}}$ to (3). This can be showed considering two consecutive steps in the semi-implicit JKO scheme (9), i.e. $\mathit{\boldsymbol{\mu}}_\tau^k$, $\mathit{\boldsymbol{\mu}}_\tau^{k+1}$, and perturbing in the following way

where $P^\epsilon = \mathrm{id}+\epsilon\zeta$, for some $\zeta\in C_c^\infty( {\mathbb{R}^{n}}; {\mathbb{R}^{n}})$ and $\epsilon\ge0$. From the minimizing property of $\mathit{\boldsymbol{\mu}}_\tau^{k+1}$ we have

After some manipulations, sending first $\epsilon\to 0$ and then $\tau\to 0$ the inequality (12) leads to the first weak formulation in Definition 2.1. Perturbing now on $\eta$ and repeating the same procedure we get the required convergence.

2.3.   Stationary states for purely non-local systems

The existence of weak solutions to the purely non-local systems, i.e.,

with generic cross-interaction kernels $K_\rho$ and $K_\eta$ is investigated in [22], whereas studies on the shape of stationary states can be found in [17, 27]. Concerning the predator-prey modelling and patterns formation, in [15, 46] a minimal version of (1) has been considered with only one predator and arbitrarily many prey subject to (different) singular potentials. This model induces the formation of nontrivial patterns in some way to prevent the action of the predators. In [21] the authors study existence and stability of stationary states for the purely aggregative version of system (3), namely equation (3) with $\theta = 0$,

and its relation with the particle system

It is proved that stationary states of system (14) are linear combinations of Dirac's deltas, namely $\bar{\rho},\bar{\eta}\in {\mathscr{P}}( \mathbb{R}^n)$, with

where $\left\{\bar{X}_{k}\right\}_k$, $\left\{\bar{Y}_{h}\right\}_h$ are stationary solutions of system (15), i.e.,

for $i = 1,...,N$ and $j = 1,...,M$, see also [29, 30] for a symilar result in the one-species case. As pointed out in [21], system (17) is not enough to determine a unique steady state, since the linear combination of the first $N$ equations, weighted with $\alpha m_X^i$, and the final $M$ equations weighted with coefficients $-m_Y^j$ get the trivial identity $0 = 0$. System (17) should be coupled with the quantity

that is a conserved quantities, and therefore one would like to produce a unique steady state once the quantity $C_\alpha$ is prescribed. Solutions to system (17) will play a crucial role in the proof of the main Theorem 1.1.

2.4.   Existence of mixed steady state via Krein-Rutman theorem

As a preliminary result, we now prove the existence of one possible shape of steady state, that will be a prototype example for the general case. The steady state is what we can call a mixed steady state, that identifies the case in which the predators can catch the prey, see Figure 2.

The proof of the existence of such steady state follows by using the strong version of the Krein-Rutman theorem, see [26].

Theorem 2.2 (Krein-Rutman). Let $X$ be a Banach space, $K\subset X$ be a solid cone, such that $\lambda K\subset K$ for all $\lambda\geq0$ and $K$ has a nonempty interior $K^{o}$. Let $T$ be a compact linear operator on X, which is strongly positive with respect to $K$, i.e. $T[u]\in K^{o}$ if $u\in K\setminus\{0\}$. Then,

(i) the spectral radius $r(T)$ is strictly positive and $r(T)$ is a simple eigenvalue with an eigenvector $v\in K^{o}$. There is no other eigenvalue with a corresponding eigenvector $v\in K$.

(ii) $|\lambda|<r(T)$ for all other eigenvalues $\lambda\neq r(T)$.

As pointed out in [6], using this strategy we can only obtain existence of stationary states for a diffusion coefficient that depends on the support, without providing any one-to-one correspondence between the diffusion constant (eigenvalue) and the support. Even if non completely satisfactory, the following results give a useful insight on the possible conditions we can expect in order to get existence of steady states, see Remark 3.1.

Let us first introduce a proper definition for mixed steady states as in Figure 2.

Definition 2.2. Let $0<R_\rho<R_\eta$ be fixed. We call a pair $(\rho,\eta)$ a mixed steady state a solution to system (3) with $\rho$ and $\eta$ in $L^1\cap L^{\infty}( \mathbb{R})$, non-negative, symmetric and radially decreasing functions with supports

Let us now assume that $(\rho,\eta)$ is a steady state to system (3) as in Definition 2.2, then (4) can be rephrased as

where $C_{\rho},C_{\eta}>0$ are two constants. Differentiating the two equations in (19) w.r.t. $x\in{\rm{supp}}(\rho)$ and $x\in{\rm{supp}}(\eta)$ respectively, we obtain

By symmetry properties of the kernels $S_\rho$, $S_\eta$ and $K$ and the steady states $\rho$ and $\eta$, for $x>0$, we get

In order to simplify the notations, we set

Notice that $\tilde{G}$, under assumptions (A1)-(A3), is a nonnegative function for $x,y>0$. We also set $p(x) = -\rho_{x}(x)$ for $x\in (-R_\rho,R_\rho)$ and $q(x) = -\eta_{x}(x)$ for $x\in (-R_\eta,R_\eta)$. Hence, (21) is rewritten simply as

Proposition 2.1. Assume that $S_\rho,S_\eta,K$ satisfy (A1), (A2) and (A3) and fix $0<R_\rho<R_\eta$ and $0<z_\rho$, $z_\eta$. Assume that $S_\eta$ and $K$ are strictly concave on $[-R_\eta,R_\eta]$ and $[-R_\rho,R_\rho]$ respectively. Assume that there exists a constant $C$ such that

for any $(k,h)\in C^{1}([0,R_\rho])\times C^{1}([0,R_\eta])$ with $k(0) = h(0) = 0$ and $k^{\prime}(0)>0$, $k (x)>0$ for all $x\in (0,R_\rho)$, and $h^{\prime}(0)>0$, $h(x)>0$ for all $x\in (0,R_\eta)$. Then, there exists a unique mixed steady state $(\rho,\eta)$ in the sense of Definition 2.2 to system (3) with $\theta = \theta(R_\rho,R_\eta)>0$, provided that

where $z_\rho$ and $z_\eta$ are masses of $\rho$ and $\eta$ respectively.

Proof. Let us first introduce the following Banach space

endowed with the $W^{1,\infty}$-norm for the two components $p$ and $q$. Define the operator $T_{R_\rho,R_\eta}[P]$ on the Banach space $\mathcal{X}_{R_\rho,R_\eta}$ as

where $P$ denotes the elements $P = (p,q)\in \mathcal{X}_{R_\rho,R_\eta}$, and $(f,g)$ are given by

The assumptions on the kernels ensure that the operator $T_{R_\rho,R_\eta}$ is compact on the Banach space $\mathcal{X}_{R_\rho,R_\eta}$. Indeed, as the operator $T_{R_\rho,R_\eta}$ is defined on a space of $C^1$ functions defined on compact intervals, and since the kernels $S_\rho,S_\eta,K$ are all in $C^2$, which is from assumption (B1), defined in the operator on compact intervals, then using Arzelá's theorem, it is easy to prove that $T_{R_\rho,R_\eta}$ maps bounded sequences in $\mathcal{X}_{R_\rho,R_\eta}$ into pre-compact ones. Now, consider the subset $\mathcal{K}_{R_\rho,R_\eta}\subseteq \mathcal{X}_{R_\rho,R_\eta}$ defined as

It can be shown that this set is indeed a solid cone in $\mathcal{K}_{R_\rho,R_\eta}$. Moreover, we have that

where $\overset{\circ}{ \mathcal{K}_{R_\rho,R_\eta}}$ denotes the interior of $\mathcal{K}_{R_\rho,R_\eta}$. We now show that the operator $T_{R_\rho,R_\eta}$ is strongly positive on the solid cone $\mathcal{K}_{R_\rho,R_\eta}$ in the sense of Theorem 2.2. Let $(p,q)\in \mathcal{K}_{R_\rho,R_\eta}$ with $p,q\neq 0$, then by the definition of the operator $T_{R_\rho,R_\eta}$, it is easy to see that the first component is non-negative. Concerning the second component, we have

if and only if $\alpha <C$ with $C$ as in (23). Next, it is easy to show that the derivative at $x = 0$ of the first component is strictly positive. The derivative of the second component is given by

Now, we need to find the condition on $\alpha$ such that $A>0$. Chebyshev's inequality in the first integral of $A$ and the concavity assumption of $S_\eta$ on the interval $[-R_\eta,R_\eta]$ yields the bound

Using the concavity assumption of $K$ on the interval $[-R_\rho,R_\rho]$, the other integral can be easily bounded by

Thus, $A>0$ holds under the condition

As a consequence, $T_{R_\rho,R_\eta}[P]$ belongs to $\mathcal{H}_{R_\rho,R_\eta}$, which implies that the operator $T_{R_\rho,R_\eta}$ is strongly positive on the solid cone $\mathcal{K}_{R_\rho,R_\eta}$. Then, the Krein-Rutman theorem applies and guarantees the existence of an eigenvalue $\theta = \theta(R_\rho,R_\eta)$ such that

with an eigenspace generated by one given nontrivial element $(\bar{p},\bar{q})$ in the interior of the solid cone $\mathcal{K}_{R_\rho,R_\eta}$. Moreover, by (i) of Theorem 2.2, there exists no other eigenvalues to $T_{R_\rho,R_\eta}$ with corresponding eigenvectors in $\mathcal{K}_{R_\rho,R_\eta}$ besides the one with eigenfunction $(\bar{p},\bar{q})$, and by (ii) of Theorem 2.2 all other eigenvalues $\tilde{\theta}$ with eigenfunctions in $\mathcal{X}_{R_\rho,R_\eta}$ satisfy $|\tilde{\theta}|<\theta$.

Remark 2.1. Motivated from the purely aggregative case (13), we actually expect a more rich behavior for the possible steady states configurations, such as the separated stationary state in Figure 3. This is expected as a possible winning strategy for the prey, since it corresponds to prey finding a safe distance from the predators. Unfortunately, the symmetrization argument used in the previous proof seems to fail, since in this case we need to require the symmetry around $cm_\eta^2$ for the convolutions. In the next section we prove the existence and uniqueness for the multiple bumps stationary state in the sense of Definition 1.1, that includes the mixed and the separated one, using a completely different approach.

3.   Existence for Multiple Bumps Steady States

In this Section we prove the existence and uniqueness of a multiple bumps steady state in the sense of Definition 1 fixing masses and a small diffusion coefficient. Following the approach in [5, 6], we first formulate the problem in terms of the pseudo-inverse functions and then we use the Implicit Function Theorem (cf. [19, Theorem 15.1]).

We start rewriting our stationary system in terms of pseudo-inverse functions. Let $(\rho,\eta)$ be a solution to the stationary system

where $\alpha_{\rho} = 1$ and $\alpha_{\eta} = -\alpha$. Assume that $(\rho,\eta)$ have masses $z_\rho$ and $z_\eta$ respectively and denote by $cm_{l},\; l \in \{\rho,\eta \},$ the centres of masses

Remember that the only conserved quantity in the evolution, together with the masses, is the joint centre of mass

that we can consider fixed. Define the cumulative distribution functions of $\rho$ and $\eta$ as

Let $u_{l}:[0,z_l]\rightarrow \mathbb{R},\; l \in \{\rho,\eta\},$ be the pseudo-inverse of $F_l,$ namely

supported on

For $\rho$ and $\eta$ multiple bumps in the sense of Definition 1.1 we can denote the mass of each bump as

and the centres of masses accordingly,

and we can always assume that the centres of masses are ordered species by species, i.e. $cm_l^i\geq cm_l^j$ if $i\geq j$. Let us consider the case of centres of masses that are stationary solutions of the purely non-local particle system (17), that we recall for the reader convenience,

coupled with the conservation of the joint centre of mass $CM_\alpha$ in (27), see the discussion in Section 2. For such a density the pseudo-inverse $u_l$ reads as

where

with $z_l^0 = 0$ and $z_l = \sum_{k = 1}^{N_l} z_l^k.$ We are now in the position of reformulating (26) in terms of the pseudo-inverse functions as follows:

The restriction to $z\in J_{l}^i$, $i = 1,2,\cdots,N_l,$ and $l \in \{\rho,\eta\},$ allow us to rephrase (29) in the compact form

Similar to [5, 6], we suggest the linearization

with $v_l^i$, being odd functions defined on $J_l^i$. Using this ansatz in (30), with the scaling $\theta = \delta^{3}$ we have

Multiplying (31) by $\delta\partial_{z}v_l^i$, and taking the primitives w.r.t. $z$, we obtain

where $A_l^i$ are the integration constants. In order to recover the constants $A_l^i$, we substitute $\tilde{z}_l^i$ into equation (32). Denoting by $v_l^i(\tilde{z}_l^i) = \lambda_l^i,$ we obtain

Next, we set $G_{l}$ and $H$ such that $G_l^{\prime} = S_l$ and $H^{\prime} = K$, with $G_l,\; H$ to be odd and satisfy $G_l(0) = H(0) = 0$. Then, multiplying (32) again by $\delta\partial_{z}v_l^i$ and taking the primitives w.r.t. $z\in J_l^i$, we obtain

Let us denote with $\bar{z}_l^i$ the middle point of each interval $J_l^i.$ Then, in order to recover the integration constants $\beta_l^i$, we substitute $\bar{z}_l^i$ in (34) which yields

As a consequence of all above computations, we get the following relation for $z\in J_l^i,$

If we define, for $p = (v_{\rho}^1,\ldots,v_{\rho}^{N_\rho},v_{\eta}^1,\ldots,v_{\eta}^{N_\eta})$

we have that (30) reduces to the equation $\mathcal{F}_l^i[p;\delta](z) = 0$. In the following we compute the Taylor expansion of $\mathcal{F}_l^i[p;\delta](z)$ around $\delta = 0$. Let us begin with the first integral on the r.h.s. of (37), i.e.,

where we used the fact that $\int_{J_l^i}v_l^i(\xi)\,d\xi = 0$ and $R(S'''_l,\delta^4)$ is a remainder term. For the second integral we have

Summing up the contributions in (38) to (39), we get that the self-interaction part in (37) reduces to

Similarly, for the cross-interaction part we obtain

Putting together the contributions of (40) and (41) in the functional equation (37), we get

where we used the notations introduced in (6) and (5), namely

and

Note that since the values $cm_l^i$ satisfy (28) we have that $B_{l}^i = 0$. After the manipulations above, equation $\mathcal{F}_l^i[p;0](z) = 0$ reads as

that gives a unique solution once the value of $\lambda_l^i$ is determined. In order to do that, we evaluate the functional $\mathcal{F}_l^i$ in the end point $\tilde{z}_l^i$ of the corresponding interval

Performing Taylor expansions similar to the ones in (38) and (39) we get that

and we are now in the position to solve

The second equation in (46) admits a solution once we have that $D_l^i > 0$, and $\lambda_l^i$ is uniquely determined by

By construction $\lambda_l^i\geq \lambda_l^j$ if $i\geq j$, and this implies that equation (46) admits the unique solution $\bar{v}_l^i$ which can be recovered as the pseudo inverse of the following Barenblatt type profiles

where the intervals $I_{\rho}^i = \left[l_\rho^i,r_\rho^i\right]$ and $I_{\eta}^h = \left[l_\eta^h,r_\eta^h\right]$ are determined imposing

We are now ready to reformulate (36) as a functional equation on a proper Banach space. Consider the spaces

endowed with the $L^\infty$ norm and take the product spaces

We now introduce the space $\Omega$ defined by

with elements $\omega = (v_\rho^1,\ldots,v_\rho^{N_\rho},\lambda_\rho^1,\ldots,\lambda_\rho^{N_\rho},v_\eta^1,\ldots,v_\eta^{N_\eta},\lambda_\eta^1,\ldots,\lambda_\eta^{N_\eta})$ endowed with the norm

For $\gamma>0$, calling $\tilde{J}_l^i = \left[\bar{z}_l^i,\tilde{z}_l^i\right)$, we consider the norm

and set $\Omega_\gamma : = \{ \omega \in\Omega \; :\; |||\omega|||_\gamma < +\infty\}$. For a given $\omega\in\Omega_{\frac{1}{2}}$, we define the operator $\mathcal{T}:\Omega_{\frac{1}{2}}\rightarrow \Omega_1$

where for $l\in \{ \rho,\eta \}$ we have shorted the notation introducing

The operator $\mathcal{T}$ is a bounded operator for any fixed $\delta\geq 0$ and can be continuously extended at $\delta = 0$ to (43) and (45). In order to prove existence of stationary solutions for small $\delta > 0$ using the Implicit Function Theorem, we need to prove that the Jacobian matrix of $\mathcal{T}$ is a bounded linear operator form $\Omega_{1/2}$ to $\Omega_1$ with bounded inverse. The Jacobian of $\mathcal{T}$ has the following structure

where the components are actually matrices defined by

where $\nu_h^j$ and $a_h^j$ are generic directions. We first compute the diagonal terms in the matrix $D_{v_l} {\mathcal{F}}_l(\delta)$. We have

A Taylor expansion around $\delta = 0$ similar to the ones in (38) - (41) easily gives that in the limit $\delta\to 0$ we obtain

Concerning the other terms in $D_{v_l} {\mathcal{F}}_l(\delta)$ we get

that all vanish in the limit $\delta\to 0$. Let us now focus on the matrix $D_{\lambda_l} {\mathcal{F}}_l(\delta)$. By (37) it is easy to see that the only non-zero terms in $D_{\lambda_l} {\mathcal{F}}_l(\delta)$ are the diagonal ones that are given by

Then, Taylor expansion w.r.t. $\delta$ yields

Since all the entrances in the matrix $D_{\lambda_h} {\mathcal{F}}_l(\delta)$ are zero, the last matrix that concerns ${\mathcal{F}}_l^i$ is $D_{v_h} {\mathcal{F}}_l(\delta)$. The elements of this matrix are given by

that vanish in the limit $\delta\to 0$. We now start in computing the functional derivatives for $\Lambda_l^i$ in (44). Again we should consider the four matrices in (55), and we start from $D_{v_l}\Lambda_l(\delta)$. Note that the terms in the diagonal are zero in this case and the others are given by

The terms in $D_{v_h}\Lambda_l(\delta)$ are

and the usual Taylor expansions around $\delta = 0$, show that both the matrices $D_{v_l}\Lambda_l(0) = D_{v_h}\Lambda_l(0) = 0$. Since $D_{\lambda_h}\Lambda_l(\delta)$ is trivially a zero matrix, only remains to compute the diagonal terms in $D_{\lambda_l}\Lambda_l(\delta)$. We have

The last Taylor expansion gives

We have proved that

where

and $\mathrm{dg}(A_i)$ is a diagonal matrix with elements $A_i$. Let us denote by $\omega_0$ the unique solution to (46), we have the following lemma:

Lemma 3.1. For $\delta>0$ small enough, the operator $D\mathcal{T}[\omega_0;\delta]$ is a bounded linear operator from $\Omega_{1/2}$ to $\Omega_1$.

Proof. Thanks to the previous computations it is easy to see that $D\mathcal{T}$ is a bounded linear operator from $\Omega$ into itself and it is continuous at $\delta = 0$. The definition of the norm in (52) implies that for $z\in \tilde{J}_l^i$ we need to control only that

as $\delta\searrow 0$. We start estimating the third row in (57),

where in the first equality we did a Taylor expansion around the point $x_0 = cm_l^i - cm_h^j + \delta(\lambda_l^i - v_h^j(\xi))$ for the kernel $K\big(cm_l^i - cm_h^j + \delta( v_{l}^i(z)-v_{h}^j(\xi)) \big)$, while in the second equality we did a Taylor expansion around the point $x_0 = cm_l^i - cm_h^j$ for the kernel $K''\big(cm_l^i - cm_h^j + \delta(\lambda_l^i - v_h^j(\xi))\big)$. Similarly, we can show that

The first two rows in (57) can be treated as follows,

Since the functions $v_l^i$ are components of a vector $\omega$ belonging to $\Omega_{1/2},$ the quantities

are uniformly bounded in $\tilde{J}_l^i$, that gives (57).

Lemma 3.2. For any $\delta>0$ small enough, $D\mathcal{T}[\omega_0; 0]:\Omega_{1/2}\rightarrow \Omega_{1}$ is a linear isomorphism.

Proof. Given $w\in\Omega_1$, we have to prove that

admits a unique solution $\omega\in\Omega_{1/2}$ with the property

The determinant of the matrix in (56) is given by

that is always different from zero under the condition $D_l^i > 0$ and since $\big( v_l^i(z)-\lambda_l^i\big) < 0$ on $z\in[\bar{z}_l^i, \tilde{z}_l^i)$. Thanks to the structure in (56) and denoting with $\nu_l^i,a_l^i$ and $\sigma_l^i,k_l^i$ the generic entrances in $\omega$ and $w$ respectively, we easily get that

that implies $||\nu_l^i||_{1/2}\leq C ||\sigma_l^i||_{\infty}$ for $i = 1,2,\cdots,N_l,$ and $l \in \{\rho,\eta\}$. In order to close the argument, it is enough to note that the ratio

is uniformly bounded since $(\lambda_l^i - v_l^i(z))/(\tilde{z}_l^i-z)$ is uniformly bounded, see [5, Lemma 4.4].

We are now in the position of proving the main result of the paper, namely Theorem 1.1, that we recall below for convenience.

Theorem 3.1. Assume that $S_\rho$, $S_\eta$ and $K$ are interaction kernels are under the assumptions (A1), (A2) and (A3). Consider $N_\rho,N_\eta\in\mathbb{N}$ and let $z_l^i$ be fixed positive numbers for $i = 1,2,\cdots,N_l,$ and $l \in \{\rho,\eta\}$. Consider two families of real numbers $\{cm_\rho^i\}_{i = 1}^{N_\rho}$ and $\{cm_\eta^i\}_{i = 1}^{N_\eta}$ such that

(i) $\{cm_\rho^i\}_{i = 1}^{N_\rho}$ and $\{cm_\eta^i\}_{i = 1}^{N_\eta}$ are stationary solutions of the purely non-local particle system, that is, for $i = 1,2,\cdots,N_l,$ for $l,h \in \{\rho,\eta\}$ and $l\neq h$,

(ii) the following quantities

are strictly positive, for all $i = 1,2,\cdots,N_l$, $l,h \in \{\rho,\eta\}$ and $l\neq h$.

where $\alpha_{\rho} = 1$ and $\alpha_{\eta} = -\alpha$. Then, there exists a constant $\theta_0$ such that for all $\theta\in (0, \theta_0)$ the stationary equation (4) admits a unique solution in the sense of Definition 1.1 of the form

where

● each interval $I_l^i$ is symmetric around $cm_l^i$ for all $i = 1,2,\cdots,N_l$, $l \in \{\rho,\eta\}$,

● $\rho^i$ and $\eta^j$ are $C^1$, non-negative and even w.r.t the centres of $I_\rho^i$ and $I_\eta^j$ respectively, with masses $z_\rho^i$ and $z_\eta^j$, for $i = 1,...,N_\rho$ and $j = 1,...,N_\eta$,

● the solutions $\rho$ and $\eta$ have fixed masses

respectively.

Proof. Consider $z_\rho$ and $z_\eta$ fixed masses and a set of points $cm_l^i$ for $i = 1,2,\cdots,N_l,$ and $l \in \{\rho,\eta\}$ that satisfy (i) and (ii). The results in Lemma 3.1 and Lemma 3.2 imply that given $\mathcal{T}$ defined in (53), the functional equation

admits a unique solution

for $\delta >0$ small enough. The entrances $v_l^i(z)$ are solutions to (36) for $z\in J_l^i$. Consider now $u_l^i$ defined for $z\in J_l^i$ as $u_l^i (z) = cm_l^i + \delta v_l^i(z)$. Differentiating (36) twice w.r.t $z$ we get that $v_l^i$ is differentiable and strictly increasing for $z\in J_l^i$ and that $u_l^i$ is a solution to (30). Moreover $u_l^i$ is also strictly increasing and we can define the inverse $F_l^i$ that and its spatial derivative $\rho_l^i = \partial_x F_l^i$ is a solution to (4).

Remark 3.1 Note that conditions (ii) in Theorem 3.1 turn to be conditions on the positions of the centres of masses and on the value of $\alpha$. Indeed, as a sufficient condition for $D_\rho^i> 0$ we can assume that all the differences between the centres of masses are in the range of concavity of the kernels. Moreover, $D_\eta^i> 0$ is satisfied if

Note that the above conditions are comparable to the one we got in the proof of Section 2 using the Krein-Rutman approach.

4.   Numerics and perspectives

In this section, we study numerically solutions to system (3) using two different methods, the finite volume method introduced in [12, 13] and the particles method studied in [23, 28]. We validate the results about the existence of the mixed steady state and the multiple bumps steady states. Moreover, we perform some examples to show the variation in the behaviour of the solution to system (3) under different choices of initial data and the parameter $\alpha$, which, in turn, suggests future work on the stability of the solutions to system (3). Finally, travelling waves are detected under a special choice of initial data and value for the parameter $\alpha$. We begin by sketching the particles method. This method essentially consists in a finite difference discretization in space to the pseudo inverse version of system (3)

as the following: Let $N\in\mathbb{N}$, let $\{z^i\}_{i = 1}^N$ be a sequence of points that partition the interval $[0,1]$ uniformly. Denote by $X_l^i(t): = u_l(t,z^i)$ the approximating particles of the pseudo inverse at each point $z^i$ of the partition. Assuming that the densities $\rho$ and $\eta$ are of unit masses, then we have the following approximating system of ODEs

where the densities are reconstructed as

and the same for the $\eta^i(t)$. To this end, we are ready to solve this particle system by applying the Runge-Kutta MATLAB solver ODE23, with initial positions $X_l(0) = X_{l,0} = \{ X_{l,0}^i \}_{i = 1}^N,\; l = \rho,\eta$ determined by solving

The second method we use is the finite volume method which introduced in [12] and extended to systems in [13], that consists in a $1$D positive preserving finite-volume method for system (3). To do so, we first partition the computational domain into finite-volume cells $U_{i} = [x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ of a uniform size $\Delta x$ with $x_{i} = i\Delta x$, $i \in\{-s,\; \dots,s \}$. Define

the averages of the solutions $\rho,\eta$ computed at each cell $U_{i}$. Then we integrate each equation in system (3) over each cell $U_{i}$, and so we obtain a semi-discrete finite-volume scheme described by the following system of ODEs for $\overline{\rho}^{i}$ and $\overline{\eta}^{i}$

where the numerical flux $F_{i+\frac{1}{2}}^l$, $l = \rho,\eta$, is considered as an approximation for the continuous fluxes $-\rho(\theta\rho + S_{\rho}*\rho + \alpha_{\rho} K*\eta )_{x}$ and $-\eta(\theta\eta + S_{\eta}*\eta + \alpha_{\eta} K*\rho )_{x}$ respectively. More precisely, the expression for $F_{i+\frac{1}{2}}^{\rho}$ is given by

where

and

The minmod limiter in (68) has the following definition

We have the same as the above expresions for $\eta$. Finally, we integrate the semi-discrete scheme (65) numerically using the third-order strong preserving Runge-Kutta (SSP-RK) ODE solver used in [31].

In all the simulations below, we will fix the kernels as a normalised Gaussian potentials

that are under the assumptions on the kernels (A1), (A2) and (A3). This choice helps us in better understanding the variation in the behavior of the solutions w.r.t. the change in the initial data and the parameter $\alpha$. In the first five examples, see Figures 4, 5, 6, 7 and 8 we show:

● in the first row steady states are plotted at the level of density, on the l.h.s. we compare the two methods illustrated above, while on the r.h.s. we show the evolution by the finite volume method,

● in the second row we plot the particles paths for both species obtained with the particles method,

● in the last row we show the pseudo inverse functions corresponding to the steady state densities.

The last example we present shows an interesting traveling waves-type evolution.

The first example is devoted to validating existence of mixed steady state and separated steady state. By choosing the initial data $(\rho_0,\eta_0)$ as

and fixing $\alpha = 0.1$ and $\theta = 0.4$, we obtain a mixed steady state as plotted in Figure 4. Note that, the small value of $\alpha$ allows the predators to dominate the prey which results in the shape shown in Figure 4. Next, we take the initial data as

with the same $\theta$ as above and $\alpha = 0.2$. This choice of the initial data, actually, introduce two equal attractive forces on the right and left hand sides of the predators which, in turn, fix the predators at the centre and gives the required shape of the separated steady state as shown in Figure 5. Finally, a sort of separated steady state can also be obtained starting from the same initial data as in (70) and same diffusion parameter $\theta$. If we choose $\alpha = 6$, the prey $\eta$ will have enough speed to get out of the predators region, producing a transition between the mixed state to the separated one, this is illustrated in Figure 6.

We then, test two cases where we validate the existence of steady states of more bumps. Let us fix $\theta = 0.3$. Then, a four-bumps steady state is performed and plotted in Figure 7, where we consider

as initial data and we take $\alpha = 0.05$. This way of choosing initial data produces a balanced attractive forces which in turn form the required steady state. Similarly, we can obtain a steady state of five bumps by using initial data

and $\alpha = 1$, see Figure 8.

Finally, in this example, we detect existence of traveling waves. Indeed, by choosing initial data as

$\alpha = 1$ and $\theta = 0.2$, we obtain a traveling wave that is shown in Figure 9. Once the initial data and the value of $\theta$ are fixed, if $\alpha$ is taken as a small value, then we will come out with a mixed steady state. If we take a larger value for $\alpha$ then the prey will be fast escaping from the predators. Therefore, the proper value for $\alpha$ will produce a situation where the attack speed of the predators is equal to the escape speed of the prey, which results in a traveling wave of both the densities $\rho$ and $\eta$. All the simulations above motivate further studies on stability vs. instability for such a system, together with the possible existence on traveling wave type solutions.