# Predators-Prey models with competition Part.1:
Existence, Bifurcation and qualitative properties

Henri Berestycki, Alessandro Zilio

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

► To cite this version:

Henri Berestycki, Alessandro Zilio. Predators-Prey models with competition Part.1: Existence, Bifurcation and qualitative properties. Communications in Contemporary Mathematics, 2018, 20 (07), pp.1850010. 10.1142/S0219199718500104. hal-03884837

HAL Id: hal-03884837

https://cnrs.hal.science/hal-03884837

Submitted on 5 Dec 2022

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>
# PREDATORS-PREY MODELS WITH COMPETITION
## PART I: EXISTENCE, BIFURCATION AND QUALITATIVE PROPERTIES

HENRI BERESTYCKI AND ALESSANDRO ZILIO

**ABSTRACT.** We study a mathematical model of environments populated by both preys and predators, with the possibility for predators to actively compete for the territory. For this model we study existence and uniqueness of solutions, and their asymptotic properties in time, showing that the solutions have different behavior depending on the choice of the parameters. We also construct heterogeneous stationary solutions and study the limits of strong competition and abundant resources. We then use these information to study some properties such as the existence of solutions that maximize the total population of predators. We prove that in some regimes the optimal solution for the size of the total population contains two or more groups of competing predators.

1. INTRODUCTION

Systems of reaction-diffusion equations are ubiquitous in mathematical biology, as they serve as a basic framework for modeling a diversity of biological and ecological mechanisms. In particular, the study of population dynamics often involves such systems. In a recent paper [4], we have introduced a new reaction-diffusion system to describe the emergence of territoriality in predatory animals. As a matter of fact, this question, to a large extent, remains puzzling. Specifically, our aim was to understand whether selfish and non organized behaviors of predators are sufficient mechanisms to explain the emergence of territoriality.

In this model, we only consider basic mechanisms that characterize an environment inhabited by predators and preys. Given a region $\Omega \subset \mathbb{R}^n$ with $n \le 2$ in practice (the restriction on the dimension is purely for modeling reasons) occupied by $N + 1$ groups of animals. The density of the first one, which we denote as $u$, of preys while the remaining $N$ densities, denoted

<footnote n="2010">2010 Mathematics Subject Classification. Primary: 35Q92; secondary: 35A01, 35B25, 35B32, 92D50.

Key words and phrases. systems of parabolic equations, asymptotic analysis, stability of solutions, bifurcation analysis, non-constant solutions, competition, optimization.</footnote>
by $w_1, \dots, w_N$, of predators. Each of the densities evolves in time following a typical law of Lotka-Volterra type. The model proposed in [4] is synthesized in the form of the system

$$ (1.1) \qquad \begin{cases} w_{i,t} - d_i \Delta w_i = (-\omega_i + k_i u - \mu_i w_i - \beta \sum_{j \neq i} a_{ij} w_j) w_i \\ u_t - D \Delta u = (\lambda - \mu u - \sum_{i=1}^N k_i w_i) u \end{cases} $$

for $(x,t) \in \Omega \times (0, +\infty)$, completed by homogeneous Neumann boundary and smooth initial conditions. The parameters of the model are easily explained: some terms in the equations model internal mechanism in the populations

* $D, d_1, \dots, d_N$ are the diffusion coefficients of the different populations, and thus are always considered positive in the following;

* $\lambda > 0$ is the effective reproduction rate of the preys, and $\mu \ge 0$ stands for the possible saturability of the environment due to an excess of preys;

* $\omega_i$ is the mortality coefficient of the group $i$ of predators that takes into account the starvation caused by the absence of the prey $u$, and $\mu_i$ takes into account possible saturation phenomena in the predator populations (for instance, an internal-low level-competition between member of the same density).

The other terms are, on the other hand, responsible for the interaction between different densities

* $k_1, \dots, k_N$ govern the predation rates. That is, $k_i$ is the success of predator $w_i$ in catching prey $u$ as a factor of the probability of an encounter;

* the elements $a_{ij} > 0$ represent how the presence of the density $w_j$ affects the density $w_i$, and the particular choice of the sign suggest that we only consider competing interactions. The parameter $\beta \ge 0$ on the other hand expresses the strength of the interaction: the higher the value of $\beta$, the more aggressive is the behavior of predators between different groups.

Similar models have already been introduced in the ecological and mathematical literature, starting from the seminal paper by Volterra [25] on predator-preys interactions, to the more recent contribution [13] that started the study of strongly interaction systems of elliptic equations (in that case, modeling populations of competing predators without preys), the study
on the evolution of dispersal by means of systems of many interacting predators [15], and the papers on the qualitative properties of the solutions to such systems [7, 6, 14]. The novelty in our model, that complicates the analysis but allows for more profound results, is the inclusion in the system of the equation for the preys and of the competition between the predators. A more in depth comparison with the results in the scientific literature can be found in [4], to which we refer the interested reader for more detailed biological considerations.

The main results in this paper regarding model (1.1) are summarized as follows. In the following sections we will give more general statements. For sufficiently smooth and positive initial data, the system (1.1) admits a unique, bounded and smooth solution, defined for all $t \ge 0$ (see Lemma 2.1). The competition is a driving force in the heterogeneity of the set of solutions, indeed the set of stationary solutions of the system collapse to the set of constant solutions if $\beta$ is small (see Proposition 2.2), and is very rich for $\beta$ large (see for instance Theorems 3.9, 3.15 and 5.1), but under some assumptions the asymptotic behavior $\beta \to +\infty$ can be described accurately (see Lemma 2.1 and Proposition 2.2). The solutions of the stationary system are regular, independently of the strength of the competition, and they converge to segregated configurations when $\beta \to +\infty$ (see Proposition 4.1) in which territories of different competing groups do not overlap. This allows us to define also limit solutions to (1.1) in the case $\beta = +\infty$. Aggressiveness ($\beta \gg 1$) may help to resist an invasion by a foreign group (see Propositions 2.4 and 2.8).

On the other hand, the strong competition limits the number of different groups that can survive in a given domain. Indeed we have

**Theorem.** For a given smooth domain $\Omega \subset \mathbb{R}^n$, there exist $\tilde{N} \in \mathbb{N}$ and $\bar{\beta} > 0$ such that if $\beta > \bar{\beta}$ and $w_\beta = (w_{1,\beta}, \dots, w_{N,\beta}, u_\beta)$ is a solution to (1.1) then

* either at most $\tilde{N}$ components of $(w_{1,\beta}, \dots, w_{N,\beta})$ are strictly positive and the others are zero;

* or the solution is such that

$$ \|(w_{1,\beta}, \dots, w_{N,\beta})\|_{C^{0,\alpha}(\Omega)} + \|u_{\beta} - \lambda/\mu\|_{C^{2,\alpha}(\Omega)} = o_{\beta}(1) $$

for every $\alpha \in (0, 1)$.
Furthermore, the threshold value $\tilde{N}$ has the following upper-bound:

$$ \tilde{N} \lesssim \frac{|\Omega|}{4\pi} \max_{i=1,...,k} \frac{\lambda k_i - \mu \omega_i}{d_i \mu} $$

if $n = 2$, and similar estimates hold in any dimension.

In the theorem $A \lesssim B$ stands for $A \le B + o(B)$ where $o(B)/B \to 0$ as $B \to +\infty$. This theorem states that, when $\beta$ is sufficiently large, solutions to (1.1) are either close to constant (and small) solutions, or they have at most $\tilde{N}$ (+1) non trivial components. This result has important repercussions on the biological interpretations of the model, as it imposes an upper-bound on the total number of hostile groups of predators that can survive in a given environment. Moreover, the upper-bound itself has important ecological consequences: We have explored them in [4].

Finally, under some assumptions on the coefficients, there exist a number of densities of predators $N \in \mathbb{N}_0$ and a solution $(w_1, ..., w_N, u)$ of (1.1) that maximize the total population of predators. We show furthermore that, in many cases, this optimal configuration is given by two or more densities of predators that have very aggressive behavior between each other, rather than by a simple homogeneous group that displays no aggressiveness between its components.

**Theorem** (see Theorems 6.5 and 7.4). For any given smooth domain $\Omega$, there exist a number $\tilde{N} \in \mathbb{N}$ and a solution $(w_1, ..., w_N, u)$ of (1.1) at most $\tilde{N} + 1$ non trivial components (possibly with $\beta = +\infty$) that maximizes the functional

$$ P(w_1, ..., w_N, u) = \int_{\Omega} \sum_{i=1}^{N} w_i $$

among the set of all non negative solutions of (1.1). Moreover, if $\Omega$ is a rectangular domain (in any dimension) and $\mu$ is sufficiently small, the maximum is attained by a solution with two or more densities of predators, that is, $\tilde{N} \ge 2$ and $\beta > 0$.

An interesting open problem remains for the second conclusion of the theorem, concerning the case when $\tilde{N} \ge 2$. Indeed, we believe that the result holds true for rather general domains, and numerical simulations sustain our claim, but the proof is so far elusive.

A consequence of the theorem is that competition between predators can be beneficial not only for the preys, but for the total population of predators as well.
**Structure of the paper.** The paper is organized as follows: in Section 2 we consider some basic properties of the system, such as existence and regularity of solutions, together with some asymptotic properties of the system, focusing on the stability properties of specific solutions. In Section 3, thanks to a bifurcation analysis, we show that the set of stationary solutions is very rich. In Section 4, we present some uniform estimates that we later use to give a precise description of the solutions for large competition. In Section 5, we show a more precise description of the bifurcation diagram in dimension one. In Section 6, we investigate some properties of the system with a large number of components, and finally, in Section 7, we show configurations in which the maximizers of the integral of the densities $w_i$ are spatially heterogeneous. ■

The interested reader will find a companion paper [4] to this one, where we investigate the biological and ecological interpretations of the mathematical results herein contained. The present paper is the first of a series in which we investigate properties of the system (1.1). In a second part [1] we give deeper a priori estimates of the solution to the elliptic counterpart of system (1.1). We exploit these properties in [2] to give a more precise description of the set of solutions. Finally, in [3] we prove results regarding the parabolic version of the model.

## 2. BASIC PROPERTIES OF THE SOLUTIONS

In this section we investigate some basic properties of the system. First, we establish existence and uniqueness results for the solutions. Then we analyze the long time behavior of the set of solutions. We also consider stability properties of a special class of solutions, namely those with only one predator and one prey, that is where all the $w_i$ are zero but one.

We recall that the system reads:

$$ (2.1a) \qquad \begin{cases} w_{i,t} - d_i \Delta w_i = (-\omega_i + k_i u - \mu_i w_i - \beta \sum_{j \neq i} a_{ij} w_j) w_i \\ u_t - D \Delta u = (\lambda - \mu u - \sum_{i=1}^N k_i w_i) u \end{cases} $$
in a domain $Q := \Omega \times (0, \infty)$, with $\Omega \in \mathbb{R}^n$ open, smooth, bounded and connected. It is completed by boundary and smooth initial conditions

$$ (2.1b) \qquad \begin{cases} \partial_{\nu} w_i = \partial_{\nu} u = 0 & \text{on } \partial\Omega \times (0, +\infty) \\ w_i(x, 0) = w_i^0(x) \ge 0, u(x, 0) = u^0(x) \ge 0 & \text{on } \partial\Omega \times \{0\}, \end{cases} $$

where $\nu$ denotes the unit outward normal vector field on $\partial\Omega$. We start with the following existence result

**Lemma 2.1.** Let $(w_1^0, \dots, w_N^0, u^0) \in C^{0,\alpha}(\Omega)$ be a non-negative initial condition for the system (2.1). There exists a unique solution $(w_1, \dots, w_N, u) \in C_x^{2,\alpha}C_t^{1,\alpha/2}(Q)$ for all $\alpha \in (0,1)$ which is defined globally for all $t > 0$. Moreover the solution is bounded in $L^\infty(Q)$ and for any $\varepsilon > 0$ there exists $T_\varepsilon > 0$ such that

$$ \sup_{(x,t) \in \Omega \times [T_\varepsilon + \infty)} w_i(x,t) \le \frac{\lambda k_i - \mu \omega_i}{\mu \mu_i} + \varepsilon $$

$$ \sup_{(x,t) \in \Omega \times [T_\varepsilon + \infty)} u(x,t) \le \frac{\lambda}{\mu} + \varepsilon. $$

Consequently, if there exists and index $i \in \{1, \dots, N\}$ such that $\lambda k_i \le \mu \omega_i$, then

$$ \lim_{t \to +\infty} \sup_{x \in \Omega} w_i(x, t) = 0. $$

In view of the last property, we shall also assume in the following that

(H) the relation $\lambda k_i > \mu \omega_i$ holds for all $i = 1, \dots, N$.

*Proof.* The existence of solution for $t \in [0, t_0]$ with $t_0 > 0$ small follows from standard arguments, since the semi-linear terms of the system are locally Lipschitz continuous: in order to extend the existence result for all time $t > 0$, it suffices to show an a priori $L^\infty$ uniform bound on the solutions.

First of all, we can observe that each single equation of the system (2.1) is satisfied by the trivial solution ($u = 0$ and $w_i = 0$ for some $i$). Consequently, the comparison principle applied to each equation implies that the solutions, when defined, are strictly positive for positive $t$. Using this information, we focus our attention on the equation satisfied by the density $u$, that
is

$$ (2.2) \qquad \begin{cases} u_t - D\Delta u = (\lambda - \mu u - \sum_{i=1}^N k_i w_i) u \\ u(x, 0) = u^0(x). \end{cases} $$

Let $U \in C^1(\mathbb{R}^+)$ be the solution of the initial value problem

$$ \left\{ \begin{aligned} \ddot{U} &= \lambda U - \mu U^2 && \text{for } t > 0 \\ U(0) &= \max\{\lambda/\mu, \sup_{x \in \Omega} u^0(x)\} > 0. \end{aligned} \right. $$

The family of solution $U$ is decreasing in $t > 0$ and $U(t) \to \lambda/\mu$ as $t \to +\infty$: as a result, for any $\varepsilon > 0$ there exists $T_\varepsilon \ge 0$ finite such that $U(t) \le \lambda/\mu + \varepsilon$ for any $t \ge T_\varepsilon$. Clearly, since each $w_i$ is non-negative, $u(x,t) \le U(t)$ for all $x \in \Omega$. Therefore, $u(x,t)$ is then bounded uniformly. Taking into account this information, we see that each $w_i$ satisfies

$$ (2.3) \qquad \begin{cases} w_{i,t} - d_i \Delta w_i \le (-\omega_i + k_i U - \mu_i w_i - \beta \sum_{j \ne i} a_{ij} w_j) w_i \\ w_i(x, 0) = w_i^0(x). \end{cases} $$

Using a similar reasoning as before, we can introduce the auxiliary function $W_i \in C^1(\mathbb{R}^+)$ solution to the initial value problem

$$ \left\{ \begin{aligned} \dot{W}_i &= (-\omega_i + k_i U - \mu_i W_i) W_i && \text{for } t > 0 \\ W_i(0) &= \max \left\{ \frac{\lambda k_i - \mu \omega_i}{\mu \mu_i}, \sup_{x \in \Omega} w_i^0(x) \right\} > 0. \end{aligned} \right. $$

Clearly $W_i$ is uniformly bounded in $t$ and moreover, $W_i(t) \to (\lambda k_i - \mu \omega_i)/(\mu \mu_i)$ as $t \to +\infty$. Hence, we deduce that for any $\varepsilon > 0$ there exists $T_\varepsilon \ge 0$ finite such that $W_i(t) \le (\lambda k_i - \mu \omega_i)/(\mu \mu_i) + \varepsilon$ for any $t \ge T_\varepsilon$. Again, since each $w_i$ is non negative and $u \le U$, we see that $W_i$ is a super-solution for (2.3) and thus $w_i$ is bounded uniformly.

The previous uniform upper bounds are enough to ensure that the solution can be extended for all time $t > 0$ and also yield the asymptotic estimates. $\square$

Before going further we recall a result in [9] about the asymptotic behavior in time of solutions to systems of reaction diffusion equations. We let $L$ be the Lipschitz constant of the
semi-linear term in (1.1) on the invariant region of Lemma 2.1. That is, letting

$$F(s_1, \dots, s_N, S) = \begin{pmatrix} (-\omega_i + k_i S - \mu_i s_i - \beta \sum_{j \neq i} a_{ij} s_j) s_i \\ (\lambda - \mu S - \sum_{i=1}^N k_i s_i) S \end{pmatrix}$$

we define

$$L = \sup \left\{ |\nabla F(s_1, \dots, s_N, S)| ; 0 < s_i < \frac{\lambda k_i - \mu \omega_i}{\mu \mu_i}, 0 < S < \frac{\lambda}{\mu} \right\}.$$

We observe that, thanks to the assumptions, $L$ is finite and positive. We also let

$$d = \min\{d_1, \dots, d_N, D\}$$

and define $\gamma_1$ to be the first non trivial (that is, positive) eigenvalue of the Laplace operator $-\Delta$ in $\Omega$ with homogeneous boundary conditions. Finally, for any solution $(w_1, \dots, w_N, u)$ of (1.1), we let

$$\bar{w}_i(t) = \frac{1}{|\Omega|} \int_{\Omega} w_i(x, t) dx, \quad \bar{u}(t) = \frac{1}{|\Omega|} \int_{\Omega} u(x, t) dx.$$

Applying [9, Theorem 3.1] to our system (1.1) we have the following result on the asymptotic behavior of the solutions for large time.

**Proposition 2.2.** Let

$$\sigma = d\gamma_1 - L.$$

If $\sigma > 0$, then for any non negative initial condition $(w_1^0, \dots, w_N^0, u^0) \in C^{0,\alpha}(\Omega)$, the corresponding unique solution of the system (2.1) converges exponential towards spatially homogeneous solutions, that is, for any $0 < \sigma' < \sigma$ there exists a constant $C > 0$ such that

$$\sum_{i=1}^{N} \|\nabla w_i\|_{L^2(\Omega)} + \|\nabla u\|_{L^2(\Omega)} \leq C e^{-\sigma't}$$

$$\sum_{i=1}^{N} \|w_i(\cdot, t) - \bar{w}_i(t)\|_{L^{\infty}(\Omega)} + \|u(\cdot, t) - \bar{u}(t)\|_{L^{\infty}(\Omega)} \leq C e^{-\sigma't/n}.$$

Moreover, the vector $(\bar{w}_1, \dots, \bar{w}_M, \bar{u})$ is a solution of the system of ordinary differential equations of the form

$$\begin{cases} \bar{w}'_i = (-\omega_i + k_i \bar{u} - \mu_i \bar{w}_i - \beta \sum_{j \neq i} a_{ij} \bar{w}_j) \bar{w}_i + g_i(t) \\ \bar{u}'_t = (\lambda - \mu \bar{u} - \sum_{i=1}^N k_i \bar{w}_i) \bar{u} + g(t) \end{cases}$$
with

$$ \bar{w}_i(0) = |\Omega|^{-1} \int_{\Omega} w_i(x)^0 dx, \quad \bar{u}(0) = |\Omega|^{-1} \int_{\Omega} u(x)^0 dx. $$

and

$$ \sum_{i=1}^{N} |g_i(t)| + |g(t)| \leq C e^{-\sigma't}. $$

*Proof.* The proof is a straightforward application of [9, Theorem 3.1]. We only observe that by Lemma 2.1 we know that from any positive initial data and any $\varepsilon > 0$ there exists $T_\varepsilon > 0$ such that the corresponding unique solution is contained in the region

$$ \{0 < w_i(x,t) < \frac{\lambda k_i - \mu \omega_i}{\mu \mu_i} + \varepsilon, 0 < u(x,t) < \frac{\lambda}{\mu} + \varepsilon, \forall x \in \Omega\} $$

for all $t \geq T_\varepsilon$. Now, if $\sigma > 0$, by regularity of $F$ for any $\varepsilon > 0$ sufficiently small

$$ \sigma' = d\gamma_1 - \sup \left\{ |\nabla F(s_1, \dots, s_N, S)|; 0 < s_i < \frac{\lambda k_i - \mu \omega_i}{\mu \mu_i} + \varepsilon, 0 < S < \frac{\lambda}{\mu} + \varepsilon \right\} > 0 $$

and we can apply [9, Theorem 3.1] to obtain the stated exponential estimates. $\square$

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

The important consequence of the previous proposition is that the behavior of the solutions, in the regime $\sigma > 0$ is well described by the corresponding system of ordinary differential equations. It also gives us a complete characterization of the set of stationary solutions of (2.1), which is then given by the (spatially constant) solutions of $F(w_1, \dots, w_N, u) = 0$. For instance, by studying the stability of the stationary and homogeneous solutions (see Proposition 2.4 and Lemma 3.3 below) we will see that in this case, when $\beta > 0$ the only stable stationary solutions are those that have $u > 0$ and only one component of $(w_1, \dots, w_N)$ non trivial (and positive).

Finally, we observe that the condition $\sigma > 0$ can be violated in three different ways:

(i) lowering the diffusion coefficients,

(ii) enlarging the domain or

(iii) augmenting the Lipschitz constant $L$.

This last possibility, which is the one that we mainly explore later, can be result for instance from taking $\mu$ small and $\beta$ large enough.

We now start investigating the equilibria of the system, in particular we want to analyze what is the impact of the competition parameter on the possible heterogeneity of the solutions of the system.
We first recall the well known result by Dockery et al. [15] on a related simpler model

$$ (2.4) \qquad \begin{cases} w_{i,t} - d_i \Delta w_i = (a(x) - \sum_{j=1}^N w_j) w_i & \text{in } \Omega \times (0, +\infty) \\ \partial_v w_i = 0 & \text{on } \partial\Omega \times (0, +\infty). \end{cases} $$

This system describes N populations that share the same spatially distributed resource a. These populations do not compete actively against each other, but they do suffer from overpopulation, which is model by the logistic term in the equations. Here a is a smooth non constant function such that the principal eigenvalue of each of the elliptic operators

$$ \begin{cases} -d_i \Delta w = aw + \lambda w & \text{in } \Omega \\ \partial_v w = 0 & \text{on } \partial\Omega, \end{cases} $$

denoted by $\lambda(d_i, a)$, is strictly negative (implying, in particular, the instability of the zero solution). Exploiting the particular symmetric structure of the interaction/competition term, Dockery et al. were able to show that if a is not constant, the only asymptotically stable equilibrium of the system is the stationary solution that has all the components $w_i$ zero except for the one with the smallest diffusion coefficient $d_i$. Moreover, the same result holds if we introduce a small mutation term in the system, which in terms imply also an evolutionary advantage for small diffusion rates. The classic interpretation of this result is that, since the densities $w_i$ in (2.4) are equivalent if not for the diffusion rates, the density which can concentrate more on favorable zones (maxima of a) will benefit more than the others and will end up eliminating them.

In what follows, we shall show that this is not the case for the model we are considering, and in particular we prove that for $\beta$ sufficiently large, all the solutions that have only one nontrivial density of predators are asymptotically stable.

**Remark 2.3.** In order to justify the link between the model (2.4) and our model (2.1), let us consider the limit case of (2.1) in which the density $u$ has a very fast dynamic with respect to the other components, that is, let us assume that for each $t > 0$, the density $u$ reaches
instantaneously its non-trivial inviscid equilibrium state,

$$ \lambda u - \mu u^2 - u \sum_{i=1}^{N} k_i w_i = 0 \implies u = \frac{1}{\mu} \left( \lambda - \sum_{i=1}^{N} k_i w_i \right). $$

Substituting the previous identity in the equations satisfied by $w_i$ we obtain

$$ w_{i,t} - d_i \Delta w_i = \left( \frac{k_i \lambda}{\mu} - \omega_i - \frac{k_i}{\mu} w_i - \sum_{j \neq i} \left( \beta a_{ij} + \frac{k_i}{\mu} \right) w_j \right) w_i $$

In the simplified case $k_i = \mu$, $\omega_i = \omega$ and $\beta = 0$, we obtain finally

$$ w_{i,t} - d_i \Delta w_i = \left( \lambda - \omega - \sum_{j=1}^{N} w_j \right) w_i = \left( a - \sum_{j=1}^{N} w_j \right) w_i. $$

We thus obtain the model of Dockery et al. [15] with $a = \lambda - \omega$. Notice that we could consider that $\lambda$ and $\omega$ depend on the location in space (certain locations being more favorable than others). More details can be found in [4].

We have the following

**Proposition 2.4.** For a fixed $i \in \{1, \dots, N\}$, let W be the stationary solution of (1.1) which has only the i-th densities of predator which is nontrivial. Then **v** is constant and **v** = (0, ..., $\tilde{w}_i$, ..., 0, $\tilde{u}$) with

$$ \tilde{w}_i = \frac{\lambda k_i - \mu \omega_i}{k_i^2}, \quad \tilde{u} = \frac{\omega_i}{k_i}. $$

There exists $\bar{\beta} \ge 0$ such that if $\beta \ge \bar{\beta}$ then W is asymptotically stable with respect to perturbations in $C^{2,\alpha}(\bar{\Omega})$. More explicitly, this holds whenever $\beta$ satisfies the system of inequalities

$$ \beta \ge \frac{k_j}{a_{ji} \tilde{w}_i} \left( \frac{\omega_i}{k_i} - \frac{\omega_j}{k_j} \right) \quad \forall j \neq i. $$

*Proof.* First of all, by [17, Theorem 1] and [10], we have that the only solution of the system with only the i-th density of predator non trivial is the constant solution W. The study of the stability of this solution is based on a simple analysis of the linearized system around it: we consider the Gâteaux differential around **v** of the operator describing the system, which is given by

$$ L(v)[w_1, \dots, w_N, u] = \begin{cases} -d_i \Delta w_i - k_i \tilde{w}_i u + \beta \tilde{w}_i \sum_{j \neq i} a_{ij} w_j & \text{for } j \neq i \\ -d_j \Delta w_j + \left[ k_j \left( \frac{\omega_i}{k_j} - \frac{\omega_i}{k_i} \right) + \beta \tilde{w}_i a_{ji} \right] w_j & \text{for } j \neq i \\ -D \Delta u + \mu \frac{\omega_i}{k_i} u + \omega_i w_i & \text{for } j = i \end{cases} $$
for all $(w_1, \dots, w_N, u) \in C^{2,\alpha}(\bar{\Omega})$ with homogenous Neumann boundary conditions. To ensure the stability of the solution we need to show that the spectrum of $L$ is contained in $C^+ = \{z \in \mathbb{C} : \mathrm{Re}(z) > 0\}$, that is for any $(w_1, \dots, w_N, u) \neq 0$ and $\gamma \in \mathbb{C}$

$$L(\mathbf{v})[w_1, \dots, w_N, u] = \gamma(w_1, \dots, w_N, u) \implies \mathrm{Re}(\gamma) > 0.$$

In the previous system, the components corresponding to $j \neq i$ are decoupled from the others, and thus their presence does not influence the stability of $W$. This solution $W$ is stable if and only if

$$k_j \left( \frac{\omega_j}{k_j} - \frac{\omega_i}{k_i} \right) + \beta \tilde{w}_i a_{ji} > 0 \quad \forall j \neq i$$

which gives the condition established by the proposition; indeed, under this assumption the components $w_j$ with $j \neq i$ are necessarily trivial. Let us show that this condition is enough to ensure the stability: we suppose that the previous system of inequalities holds but there exist $(w_1, \dots, w_N, u) \neq 0$ and $\gamma \in \mathbb{C}$ with $\mathrm{Re}(\gamma) \leq 0$ solution to

$$L(\mathbf{v})[w_1, \dots, w_N, u] = \gamma(w_1, \dots, w_N, u).$$

Then necessarily $w_j = 0$ for all $j \neq i$, and the system is reduced to

$$(2.5) \qquad \begin{cases} -d_i \Delta w_i = \gamma w_i + k_i \tilde{w}_i u \\ -D \Delta u = -\omega_i w_i + \left( \gamma - \mu \frac{\omega_i}{k_i} \right) u \\ \partial_\nu w_i = \partial_\nu u = 0 \end{cases} \quad \text{on } \partial \Omega$$

Since any weak solution to the previous system is regular, the stability in $C^{2,\alpha}(\Omega)$ can be deduced from the solvability of the system in $H^1(\Omega)$. To analyze it, let $\{(\gamma_h, \psi_h)\}_{h \in \mathbb{N}}$ be the spectral resolution of the Laplace operator with homogeneous Neumann boundary condition in $\Omega$ (let us recall that $\gamma_0 = 0$ and $\gamma_h > 0$ for $h > 0$); since $\{\psi_h\}_{h \in \mathbb{N}}$ is a complete basis of $L^2(\Omega)$, we can write

$$w_i = \sum_{h=0}^{\infty} a_h \psi_h \quad \text{and} \quad u = \sum_{h=0}^{\infty} b_h \psi_h$$

as series converging in $L^2(\Omega)$. Inserting these relations in (2.5) and using the orthogonality of the eigenfunctions, we see that the linear system (2.5) is equivalent to the sequence of algebraic
eigenvalue problems

$$ egin{cases} d_i \gamma_h a_h - k_i \tilde{w}_i b_h = \gamma a_h \\ (D \gamma_h + \mu \frac{\omega_i}{k_i}) b_h + \omega_i a_h = \gamma b_h \end{cases} \quad \text{for } h \in \mathbb{N}. $$

By direct inspection, we can observe that $a_h = 0$ if and only if $b_h = 0$. Thus, solving the first equation in $b_h$ and substituting the result in the second, we find that $\gamma$ must be a solution to

$$ (D\gamma_h + \mu\frac{\omega_i}{k_i} - \gamma)(d_i\gamma_h - \gamma) + k_i\tilde{w}_i\omega_i = 0 $$

that is

$$ \gamma = \frac{1}{2} \left[ \left( (D+d_i)\gamma_h + \mu \frac{\omega_i}{k_i} \right) \pm \sqrt{ \left( (D+d_i)\gamma_h + \mu \frac{\omega_i}{k_i} \right)^2 - 4k_i \tilde{w}_i \omega_i } \right] $$

and in particular $\text{Re}(\gamma) > 0$. $\square$

We observe that the diffusion rates do not play any role in the stability of the solutions, while a crucial role is played by the ratio $\omega_i/k_i$. In particular if $i$ is such that

$$ \frac{\omega_i}{k_i} < \frac{\omega_j}{k_j} \quad \forall j \neq i $$

then the solution $W$ is asymptotically stable also in a slightly cooperative environment, that is for $\beta < 0$ and small in absolute value. This is a consequence of the fact that the semi-trivial solutions are constant and the different diffusion rates do not play a direct role in the stability of the solution (that is, advantage of low/high diffusion rate). In this setting, the quantity $\omega_i/k_i$ can be interpreted as the fitness of the $i$-th population.

One could then wonder whether the previous stability result is a spurious consequence either of the fact that the simple solutions are constant or of another specific feature of this particular formulation of the system. To clarify this issue, we shall now adapt the proof to a very general framework. Let us consider the following operator

$$ S_{\beta}(\mathbf{v}) := \begin{cases} \mathcal{L}_i w_i - [f_i(x, u, w_i) - \beta \sum_{j \neq i} g_{ij}(x, w_i, w_j)] w_i & \text{for all } i \in \{1, \dots, N\} \\ \mathcal{L} u - f(x, u, w_1, \dots, w_N) u & \end{cases} $$

defined for $\mathbf{v} = (w_1, \dots, w_N, u)$ in the set

$$ \mathcal{F}(\Omega) = \{ \mathbf{v} \in \mathcal{C}^{2,\alpha}(\bar{\Omega}; \mathbb{R}^{N+1}) : \partial_{\nu_i} w_i = \partial_{\nu} u = 0 \text{ on } \partial\Omega \} $$
where the respective operators $\mathcal{L}_i$ and $\mathcal{L}$ stand for linear strongly elliptic operators of the form

$$ \mathcal{L}_i w_i = -\operatorname{div}(A_i(x) \nabla w_i), \quad \mathcal{L}u = -\operatorname{div}(A(x) \nabla u) $$

associated with some smooth and uniformly elliptic symmetric matrices $A_i$ and $A$, and the $v_i$ and $\nu$ denote the co-normal vector fields associated to the corresponding elliptic operators. We assume in the following that all the terms in the operator $S_\beta$ are sufficiently smooth to justify the following computations, and moreover we suppose that there exists positive constants $C$ such that for any $\mathbf{v} \in \mathcal{F}(\Omega)$ of non negative components we have

$$ \left\{ \begin{array}{l} f_i(x, u, w_i) \le C(1+u-w_i) \\ f(x, u, w_1, \dots, w_N) \le C(1-u) \\ g_{ij}(x, w_i, w_j) \ge 0 \end{array} \right. $$

Based on the previous notation, a function $\mathbf{v} \in \mathcal{F}(\Omega)$ is a solution of the generalized model if

$$ S_\beta(\mathbf{v}) = 0 $$

while a function $\mathbf{v} \in C^1(\mathbb{R}^+; \mathcal{F}(\Omega)) \cap C(\overline{\mathbb{R}}^+; \mathcal{F}(\Omega))$ is a solution to the parabolic model if

$$ \left\{ \begin{array}{ll} \partial_t \mathbf{v} + S_\beta(\mathbf{v}) = 0 & t > 0 \\ \mathbf{v}(0) = \mathbf{v}_0 & \mathbf{v}_0 \in \mathcal{F}(\Omega). \end{array} \right. $$

Using the previous assumptions, we have

**Lemma 2.5.** For any non-negative initial datum $\mathbf{v}_0 \in \mathcal{F}(\Omega)$ there exists a unique solution $\mathbf{v}$ of the previous parabolic problem. Moreover, there exists $T > 0$ and $M > 0$, independent of $\beta$, such that

$$ 0 \le w_1(t,x), \dots, w_N(t,x), u(t,x) \le M \quad \text{for all } t \ge T, x \in \overline{\Omega}. $$

If there exist $i \in \{1, \dots, N\}$, $t > 0$ and $x_0 \in \overline{\Omega}$ such that $w_i(t, x_0) = 0$ (respectively, $u(t, x_0) = 0$), then $w_i \equiv 0$ (respectively, $u \equiv 0$).

*Proof.* The proof follows directly from the maximum principle, and thus we omit it (see Lemma 2.1 for reasoning of this kind). □
In an analogous fashion, we have a corresponding result for the stationary model. Among the class of all possible solutions, we are interested in the particular case of solutions that have only one component among the first N which is non-trivial.

**Definition 2.6.** For a given $i \in \{1, \dots, N\}$, a solution $\mathbf{v} \in \mathcal{F}(\Omega)$ is said to be *i*-simple if $w_j \equiv 0$ for all $j \neq i$ and the other components are positive.

Let us observe that if $\mathbf{v} \in \mathcal{F}(\Omega)$ is an *i*-simple solution for $S_\beta$, then it is an *i*-simple solution for any value of $\beta$.

For a given solution $\mathbf{v} \in \mathcal{F}(\Omega)$, let $L(\mathbf{v})$ be the Gâteaux derivatives of $S_\beta$ in $\mathcal{F}(\Omega)$, that is for any $\boldsymbol{\varphi} \in \mathcal{F}(\Omega)$:

$$L(\mathbf{v})[\boldsymbol{\varphi}] = \lim_{\varepsilon \to 0} \frac{S_{\beta}(\mathbf{v} + \varepsilon\boldsymbol{\varphi}) - S_{\beta}(\mathbf{v})}{\varepsilon} = \left\{ \begin{array}{ll} \mathcal{L}_i \varphi_i - [f_i(x, u, w_i) - \beta \sum_{j \neq i} g_{ij}(x, w_i, w_j)] \varphi_i & \\ -f_{i,u}(x, u, w_i)w_i\varphi - f_{i,w_i}(x, u, w_i)w_i\varphi_i & \\ +\beta \sum_{j \neq i} g_{ij,w_i}(x, w_i, w_j)w_i\varphi_i + \beta \sum_{j \neq i} g_{ij,w_j}(x, w_i, w_j)w_i\varphi_j & \\ \mathcal{L}\varphi - f(x, u, w_1, \dots, w_N)\varphi - f_{,u}(x, u, w_1, \dots, w_N)\varphi & \\ - \sum_{i=1}^N f_{,w_i}(x, u, w_1, \dots, w_N)\varphi_i & \end{array} \right.$$

Analogously, for any fixed $i \in \{1, \dots, N\}$ we define the *i*-th partial derivatives $L_i(\mathbf{v})$ as the Gâteaux derivatives of $S_\beta$ in $\mathcal{F}(\Omega)$ with respect to the direction $\boldsymbol{\varphi} \in \mathcal{F}(\Omega)$ such that $\boldsymbol{\varphi} = (0, \dots, \varphi_i, \dots, 0, \varphi)$, that is

$$L_i(\mathbf{v})[\boldsymbol{\varphi}] = \left\{ \begin{array}{ll} \mathcal{L}_i \varphi_i - [f_i(x, u, w_i) - \beta \sum_{j \neq i} g_{ij}(x, w_i, w_j)] \varphi_i & \\ -f_{i,u}(x, u, w_i)w_i\varphi - f_{i,w_i}(x, u, w_i)w_i\varphi_i & \\ +\beta \sum_{j \neq i} g_{ij,w_i}(x, w_i, w_j)w_i\varphi_i & \\ 0 & \text{for } j \neq i \\ \\ \mathcal{L}\varphi - f(x, u, w_1, \dots, w_N)\varphi - f_{,u}(x, u, w_1, \dots, w_N)\varphi & \\ -f_{,w_i}(x, u, w_1, \dots, w_N)\varphi_i & \end{array} \right.$$
Accordingly, we recall that a solution $\mathbf{v} \in \mathcal{F}(\Omega)$ is (strongly) stable if any non-trivial solution $(\gamma, \varphi)$ of

$$L(\mathbf{v})[\varphi] = \gamma\varphi$$

has necessarily $\mathrm{Re}(\gamma) > 0$. For $i$-simple solutions we have

**Definition 2.7.** For a given $i \in \{1, \dots, N\}$, an $i$-simple solution $\mathbf{v} \in \mathcal{F}(\Omega)$ is one-predator stable if any non-trivial solution $(\gamma, \varphi)$ of

$$L_i(\mathbf{v})[\varphi] = \gamma\varphi$$

with $\varphi = (0, \dots, \varphi_i, \dots, 0, \varphi)$ has necessarily $\mathrm{Re}(\gamma) > 0$.

An $i$-simple solution is thus one-predator stable if it is stable with respect to all the admissible perturbations that leave unchanged the zero components $w_j$ for $j \neq i$. Clearly, if an $i$-simple solution is stable it is also one-predator stable: under suitable conditions, the inverse is true.

**Proposition 2.8.** For a given $i \in \{1, \dots, N\}$, let us assume that

$$\inf_{x \in \Omega} g_{ji}(x, 0, s) > 0 \quad \text{for all } s > 0 \text{ and } j \neq i.$$

If $\mathbf{v} \in \mathcal{F}(\Omega)$ is an $i$-simple one-predator stable solution $\mathbf{v} \in \mathcal{F}(\Omega)$, then there exists $\bar{\beta}$ such that $\mathbf{v}$ is a stable solution for all $\beta > \bar{\beta}$.

*Proof.* The $i$-simple solution $\mathbf{v} = (0, \dots, w_i, \dots, 0, u)$ is stable if

$$\left\{
\begin{aligned}
& \mathcal{L}_i \varphi_i - [f_i(x, u, w_i) - \beta \sum_{j \neq i} g_{ij}(x, w_i, 0)] \varphi_i \\
& \quad - f_{i,u}(x, u, w_i) w_i \varphi - f_{i,w_i}(x, u, w_i) w_i \varphi_i \\
& \quad + \beta \sum_{j \neq i} g_{ij,w_i}(x, w_i, 0) w_i \varphi_i + \beta \sum_{j \neq i} g_{ij,w_j}(x, w_i, 0) w_i \varphi_j = \lambda \varphi_i \\
& \mathcal{L}_j \varphi_j - [f_j(x, u, 0) - \beta g_{ji}(x, 0, w_i)] \varphi_j = \lambda \varphi_j \\
& \mathcal{L} \varphi - f(x, u, 0, \dots, w_i, \dots, 0) \varphi - f_{,u}(x, u, 0, \dots, w_i, \dots, 0) \varphi \\
& \quad - \sum_{i=1}^N f_{,w_i}(x, u, 0, \dots, w_i, \dots, 0) \varphi_i = \lambda \varphi
\end{aligned}
\right.$$
has a nontrivial solution $\varphi \in \mathcal{F}(\Omega)$ if and only if $\mathrm{Re}(\lambda) > 0$. Let us consider the equations of index $j \neq i$, which are decoupled from the other equations in the system. They read:

$$ egin{cases} 
\mathcal{L}_j \varphi_j = [f_j(x, u, 0) - \beta g_{jh}(x, 0, w_i) + \lambda] \varphi_j & \text{in } \Omega \\ 
\partial_v^{\mathcal{L}_j} \varphi_j = 0 & \text{on } \partial\Omega 
\end{cases} $$

Since the operator $\mathcal{L}_j$ is self-adjoint<fnref n="1" />, any non-trivial solution of the system must have $\lambda \in \mathbb{R}$.
The solution $v$ being an $i$-simple solution, by the maximum principle it follows that

$$ \inf_{x \in \Omega} w_i(x) = c > 0. $$

As a result, thanks to our assumptions, there exists $\bar{\beta} \geq 0$ such that

$$ \bar{\beta} \geq \sup_{x \in \Omega} \frac{f_j(x, u, 0)}{g_{ji}(x, 0, w_i)} \quad \text{for all } j \neq i. $$

Choosing $\beta > \bar{\beta}$ and testing the equation in $\varphi_j$ by $\varphi_j$ itself, we obtain

$$ \int_{\Omega} A_j(x) \nabla w_j \cdot \nabla w_j = \int_{\Omega} [f_j(x, u, 0) - \beta g_{ji}(x, 0, w_i) + \lambda] \varphi_j^2 < \lambda \int_{\Omega} \varphi_j^2 $$

thus either $\lambda > 0$ or the component $\varphi_j = 0$. On the other hand, assuming that $\mathrm{Re}(\lambda) \leq 0$, we find a contradiction with the internal stability of the solution $v$. □

### 3. EXISTENCE OF NON HOMOGENEOUS SOLUTIONS: A BIFURCATION ANALYSIS

We continue the investigation of the asymptotic properties of the system (1.1), by now studying the set of solutions of the corresponding (stationary) elliptic problem. We consider here the model (1.1) under the assumption that the domain $\Omega$ is occupied by only two groups of predators, having the same parameters. This system reads:

$$ egin{cases} 
-d\Delta w_1 = (-\omega + ku - \beta w_2) w_1 & \text{in } \Omega \\ 
-d\Delta w_2 = (-\omega + ku - \beta w_1) w_2 & \text{in } \Omega \\ 
-D\Delta u = (\lambda - \mu u - k(w_1 + w_2)) u & \text{in } \Omega \\ 
\partial_v w_i = \partial_v u = 0 & \text{on } \partial\Omega 
\end{cases} $$

<footnote n="1">¹More precisely, the operator is self-adjoint if seen as an operator acting on $H^1(\Omega)$ functions, and the conclusion can be reached using the regularity assumptions on its coefficients.</footnote>
for which we look for solutions $(w_1, w_2, u) \in C^{2,\alpha}(\bar{\Omega})$. Let us point out that here we take $\mu_1 = \mu_2 = 0$. Alternatively, we can easily generalize the results that will we show in the following to the case of positive saturation coefficients (though the computations are inevitably more involved). Since we are looking for stationary solutions, the system can be simplified by some linear substitutions. Indeed, letting

$$u \mapsto \frac{d}{D}u, \lambda \mapsto \lambda D, \mu \mapsto \mu \frac{D^2}{d}, k \mapsto kD, \omega \mapsto \omega d, \beta \mapsto \beta d$$

we can reformulate the system as

$$ (3.1) \qquad \left\{ \begin{array}{ll} -\Delta w_1 = (-\omega + ku - \beta w_2) w_1 & \text{in } \Omega \\ -\Delta w_2 = (-\omega + ku - \beta w_1) w_2 & \text{in } \Omega \\ -\Delta u = (\lambda - \mu u - k(w_1 + w_2)) u & \text{in } \Omega \\ \partial_v w_i = \partial_v u = 0 & \text{on } \partial\Omega \end{array} \right. $$

We recall the definition of the set

$$ \mathcal{F}(\Omega) := \{ (w_1, w_2, u) \in C^{2,\alpha}(\bar{\Omega}) : \partial_v w_1 = \partial_v w_2 = \partial_v u = 0 \text{ on } \partial\Omega \}. $$

We are interested in non negative solutions of the system. Letting all the other parameters of the model fixed, we shall study the set of the solutions of (3.1) by varying the competition strength $\beta$. Let us recall that the assumption (H) holds, that is, in this context, $\lambda k > \mu\omega$.

We start by recalling a result concerning the regularity of solutions of system (3.1). This result follows from Lemma 2.1, but we report it here for the reader's convenience.

**Lemma 3.1.** Let $(w_1, w_2, u) \in H^1(\Omega)$ be a non negative weak solution to (3.1). Then

* the solutions are classical. More precisely, $(w_1, w_2, u) \in C^\infty(\Omega) \cap C^{2,\alpha}(\bar{\Omega})$ for any $\alpha < 1$ and the regularity is limited only by that of $\Omega$;

* $(w_1, w_2, u)$ are non negative and bounded uniformly in $\beta$, that is

$$ (3.2) \qquad \left\{ \begin{array}{l} w_1 \ge 0, w_2 \ge 0, 0 \le u \le \lambda/\mu \\ u + w_1 + w_2 \le \frac{(\lambda + \omega)\lambda}{\mu\omega} \end{array} \right. $$

and either all the inequalities are strict or the solution is constant;
*Proof.* All the assertions in the statement are rather straightforward consequences of the maximum principle and the classical regularity theory of elliptic equations. We only observe that the last inequality follows by summing the three equations together. This yields to

$$-\Delta(u + w_1 + w_2) \leq (\lambda + \omega)u - \omega(u + w_1 + w_2) \leq (\lambda + \omega)\frac{\lambda}{\mu} - \omega(u + w_1 + w_2).$$

We conclude again by virtue of the maximum principle. □

Lemma 3.1 gives a description of the solutions of the system (3.1), but it contains no information about the existence of such solutions. In the following, our aim is to complete this aspect, showing that the set of solutions is rich. Before doing so, we need to introduce some notation.

For a given solution $(w_1, w_2, u) \in \mathcal{F}(\Omega)$ of the system (3.1), the Gâteaux derivate in $\mathcal{F}(\Omega)$ associated to (3.1) computed at $(w_1, w_2, u)$ is given by

$$L_{\beta}\varphi = -\Delta\varphi - A_{\beta}\varphi, \quad \text{for any } \varphi \in \mathcal{F}(\Omega)$$

where $A_{\beta} = A_{\beta}(w_1, w_2, u) \in C^{2,\alpha}(\bar{\Omega}, \mathbb{R}^{3\times3})$ is

$$A = A_{\beta} = \begin{pmatrix} -\omega + ku - \beta w_2 & -\beta w_1 & kw_1 \\ -\beta w_2 & -\omega + ku - \beta w_1 & kw_2 \\ -ku & -ku & \lambda - 2\mu u - kw_1 - kw_2 \end{pmatrix}.$$

The solution $(w_1, w_2, u)$ is said to be (strongly linearly) stable if any non-trivial solution $(\gamma, \varphi)$ of the linearized equation

$$L_{\beta}\varphi = \gamma\varphi$$

has necessarily $\text{Re}(\gamma) > 0$ and weakly stable we can only infer that $\text{Re}(\gamma) \geq 0$. It is said to be (strongly linearly) unstable if, on the contrary, there exists a non-trivial solution with $\text{Re}(\gamma) < 0$.

If the solution $(w_1, w_2, u)$ in the previous definition is constant, its stability can be directly deduced from the spectrum of the matrix $A_{\beta}$ or, more explicitly, from that of $-A$. We start with the simplest scenario, that is the limit case $\beta = 0$. Under this assumption, since the densities of predators do not interact directly with each other, we can simplify drastically the system and give a complete description of the set of solutions of the system.
**Lemma 3.2.** The unique non negative solutions $(w_1, w_2, u)$ to the system (3.1) with $\beta = 0$ are the two unstable constant solutions

$$ (0,0,0), \left(0,0,\frac{\lambda}{\mu}\right) $$

and the one-parameter family of (weakly) stable ones

$$ s \in [0,1] \mapsto \left( \frac{\lambda k - \mu \omega}{k^2} s, \frac{\lambda k - \mu \omega}{k^2} (1-s), \frac{\omega}{k} \right). $$

*Proof.* In this proof, we shall only classify the solutions; the study of the stability will be postponed until Lemma 3.3, where we shall address the question about stability of constant solutions for $\beta \ge 0$ more generally.

Since for $\beta = 0$ the densities of predators do not interact directly with each other, we can simplify the system introducing the new variable $V = w_1 + w_2$, which, together with $u$ is a solution of the classical (i.e. one predator) Lotka-Volterra system

$$ \left\{ \begin{array}{ll} -\Delta V = -\omega V + kVu & \text{in } \Omega \\ -\Delta u = \lambda u - \mu u^2 - kVu & \text{in } \Omega \\ \partial_{\nu} V = \partial_{\nu} u = 0 & \text{on } \partial\Omega. \end{array} \right. $$

By a classical result of Mimura [17, Theorem 1] it follows that the previous system has only constant solutions, that is solutions of the algebraic system

$$ \left\{ \begin{array}{l} (ku - \omega)V = 0 \\ (\lambda - \mu u - kV)u = 0 \end{array} \right. $$

When $V = 0$, we have the solutions $u = 0$ or $u = \lambda/\mu$ which correspond to the first two solutions in the statement (recall that $w_1$ and $w_2$ are non negative, that is, in this case, $w_1 = w_2 = 0$). On the other hand, if $u = \omega/k$, we obtain the solution $V = w_1 + w_2 = (\lambda k - \mu \omega)/k^2$. Substituting this information in (3.1) we obtain that both $w_1$ and $w_2$ are harmonic functions, hence constants. $\square$

As we shall see later, the value $\beta = 0$ corresponds to a bifurcation point of multiplicity one for the system (3.1) around the solution

$$ (w_1, w_2, u) = \left( \frac{\lambda k - \mu \omega}{2k^2}, \frac{\lambda k - \mu \omega}{2k^2}, \frac{\omega}{k} \right), $$
so that the one-parameter family of solutions of Lemma 3.2 is nothing but the branch of solutions emanating from it.

**Lemma 3.3.** When $\beta > 0$, system (3.1) admits four different types of constant solutions:

(a) the solution $(0,0,0)$, which is strongly unstable;

(b) the solution

$$w_1 = 0, w_2 = 0, u = \frac{\lambda}{k}$$

which is strongly unstable;

(c) the solutions

$$w_1 = \frac{\lambda k - \mu \omega}{k^2}, w_2 = 0, u = \frac{\omega}{k} \quad \text{and} \quad w_1 = 0, w_2 = \frac{\lambda k - \mu \omega}{k^2}, u = \frac{\omega}{k}$$

which are strongly stable;

(d) the family of solutions

$$w_1 = w_2 = \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2}, u = \frac{\lambda \beta + 2k \omega}{\mu \beta + 2k^2}$$

which are unstable for $\beta > 0$. In particular, in this latter case,

$$\sigma(A_{\beta}) = \left\{ \beta \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2}, \gamma_{1,\beta}, \gamma_{2,\beta} \right\}$$

where $\gamma_{1,\beta}$ and $\gamma_{2,\beta}$ are two, possibly complex conjugate, eigenvalues with negative real part.

Later we will prove that the solutions in Lemma 3.3 are the only solution of (3.1) when $\beta > 0$ is sufficiently small (compare with Proposition 3.14).

*Proof.* The proof is a rather straightforward computation, but we include it here in order to glean from it an interpretation of the results.

The solution $(0,0,0)$ corresponds to the matrix

$$A_{\beta} = \begin{pmatrix} -\omega & 0 & 0 \\ 0 & -\omega & 0 \\ 0 & 0 & \lambda \end{pmatrix}$$

which is already in a diagonal form. The instability of this solution is caused by the eigenvalues $-\lambda < 0$ of $-A$, which corresponds to the constant eigenfunction $(0,0,1)$. As a result, in
complete accordance with other biological models, it implies that a logistic growth law in the prey population is responsible for an (initial) exponential growth, uniform in all the domain $\Omega$, at least when the population is small. Let us observe that none of the spectral and stability properties of the trivial solution depends on the competition $\beta$.

Similar computations hold for the solution $(0, 0, \lambda/\mu)$, whose associated matrix is

$$A_{\beta} = \begin{pmatrix} \frac{\lambda k - \mu \omega}{\mu} & 0 & 0 \\ 0 & \frac{\lambda k - \mu \omega}{\mu} & 0 \\ -\lambda k / \mu & -\lambda k / \mu & -\lambda \end{pmatrix}.$$

Here the eigenvalues are $\frac{\lambda k - \mu \omega}{\mu}$ (of multiplicity 2) and $-\lambda$. The eigenspaces of the matrix is generated by the vectors

$$\begin{pmatrix} 1 \\ 0 \\ -\frac{\lambda k}{\lambda k - \mu\omega + \lambda\mu} \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \\ -\frac{\lambda k}{\lambda k - \mu\omega + \lambda\mu} \end{pmatrix} \quad \text{and} \quad \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}$$

respectively.

To discuss the solutions of type (c), consider for example the solution $w_1 = (\lambda k - \mu \omega)/k^2$, $w_2 = 0$ and $u = \omega/k$. In this case the matrix $A_{\beta}$ is

$$A_{\beta} = \begin{pmatrix} 0 & -\beta w_1 & kw_1 \\ 0 & -\beta w_1 & 0 \\ -\omega & -\omega & -\mu\omega/k \end{pmatrix}.$$

Thus, the spectrum of $A_{\beta}$ consists of

$$-\beta w_1, -\frac{\mu\omega/k \pm \sqrt{(\mu\omega/k)^2 - 4k\omega w_1}}{2},$$

which implies strong stability of these solutions. As already observed in the previous section, this result in unchanged even when the two populations of predators have different parameters, as long as $\beta > 0$ (compare Proposition 2.4).
In the case of the constant solutions of type (d), recalling that here $w_1 = w_2$, the matrix $A_β$ reduces to

$$A_β = \begin{pmatrix} 0 & -βw_1 & kw_1 \\ -βw_1 & 0 & kw_1 \\ -ku & -ku & -μu \end{pmatrix}.$$

By direct inspection, we see that $βw_1 ≥ 0$ is an eigenvalue, implying in particular that these solutions are unstable for $β > 0$. Using this information, we can factorize the characteristic polynomial of $A_β$, yielding

$$\det(A - γ\mathrm{Id}) = (γ - βw_1) [γ^2 + (βw_1 + μu)γ + (2k^2uw_1 + βμuw_1)]$$

and the spectrum of $A_β$ consists of

$$βw_1, -\frac{(βw_1 + μu) \pm \sqrt{(βw_1 + μu)^2 - (2k^2uw_1 + βμuw_1)}}{2},$$

and this concludes the proof. □

The set of non-trivial constant solutions undergoes a transformation as $β$ changes from $β = 0$ to $β > 0$, see Lemma 3.2. Moreover, the spectrum of the matrix $A_0$, computed on the linear set of solutions is given by

$$0, -\frac{\mu\omega/k \pm \sqrt{(\mu\omega/k)^2 - 4k\omega(w_1 + w_2)}}{2}.$$

The zero eigenvalue underlines the degeneracy of the constant solutions, as they form a linear subspace, while the other two strictly negative eigenvalues confirm that this set of solutions is stable with respect to perturbations that move away from this configuration, i.e. non homogeneous perturbation (see Proposition 2.2).

**Remark 3.4.** The stability of the solutions belonging to the classes (a), (b) and (c) does not depend on $β$. More precisely, in the classes (a) and (b) the spectrum of $A_β$ is independent of $β$, while in the third case (c) the spectrum is also contained in $C^- := \{z \in C : \mathrm{Re}(z) < 0\}$.

We can say more about constant solution, and in particular we have that if a component is constant, so are the other.
**Lemma 3.5.** For a solution $(w_1, w_2, u)$ of (3.1), if one component is constant, then also all the other components are constant.

*Proof.* The case for $\beta = 0$ is already considered in Lemma 3.2. Thus we can assume $\beta > 0$. We start by assuming that $u$ is constant. If $u$ is zero, we conclude directly by the maximum principle. Assuming that $u$ is a positive constant, from the equation in $u$, we find that necessarily

$$w_1 + w_2 = \frac{\lambda - \mu u}{k}$$

is a non-negative constant. This yields in the equation for $w_i$, $i = 1,2$:

$$\begin{cases} -\Delta w_i = (-\omega + ku + \beta \frac{u}{k} u - \beta \frac{\lambda}{k} + \beta w_i) w_i \\ \partial_v w_i = 0 \end{cases} \quad \text{on } \partial \Omega.$$

Summing up the two equations, we obtain moreover

$$w_1^2 + w_2^2 = \frac{\lambda - \mu u}{k} \left( \frac{\lambda - \mu u}{k} - \frac{ku - \omega}{\beta} \right) \geq 0$$

As a result, we have obtained the identities

$$w_1 + w_2 = a, \quad w_1^2 + w_2^2 = b$$

for some non-negative constant $a$ and $b$, which directly implies that $w_1$ and $w_2$ are constant. Using this information, it is also possible to compute explicitly the solutions, and in particular we find $u = \omega/k$.

We now assume that $w_1$ is constant. From the equation in $w_1$ we find that $w_1 = 0$ or $\beta w_2 = ku - \omega$. The former case is equivalent to assuming $\beta = 0$. We then need only to address the latter. Substituting the previous identity in the equation in $u$ we find

$$\begin{cases} -\Delta u = (A - Bu) u & \text{in } \Omega \\ \partial_v u = 0 & \text{on } \partial \Omega. \end{cases}$$

where $A$ is a real constant and $B$ is a strictly positive constant. If $A \le 0$ then $u$ is zero. On the other hand, if $A > 0$ by the maximum principle (see Lemma 3.6 below), we find that $u$ is again a constant, and we can conclude as above. $\square$
Note that the arguments of the proof are only valid in the case of two predators $w_1, w_2$.
In the previous result, we made use of the following classical consequence of the maximum
principle.

**Lemma 3.6.** Let $\Omega \subset \mathbb{R}^n$ a bounded smooth domain, A and B positive constants. If $u \in H^1(\Omega)$ is a non negative solution to

$$
\left\{
\begin{array}{ll}
-\Delta u = (A - Bu)u & \text{in } \Omega \\
\partial_{\nu} u = 0 & \text{on } \partial\Omega
\end{array}
\right.
$$

then $u \equiv 0$ or $u \equiv A/B$.

Here, we are mostly interested in solutions which are not homogeneous in space (i.e., non
constant solutions). We will derive their existence through bifurcation arguments.

First, we introduce some notation.

**Definition 3.7.** We denote with $\mathcal{P} \subset \mathbb{R}^+ \times \mathcal{F}(\Omega)$ the set of all solutions $(\beta, w_{1,\beta}, w_{2,\beta}, u_\beta)$ of (3.1) with competition parameter $\beta > 0$ such that all of its components are strictly positive. Let also $\mathcal{S}_0$ stand for the set of constant solutions $(\beta, w_1, w_2, u)$ of the form

$$
w_1 = w_2 = \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2}, u = \frac{\lambda \beta + 2k \omega}{\mu \beta + 2k^2}
$$

for all values of $\beta > 0$. Lastly, we let

$$
\mathcal{S}_1 = \mathcal{P} \setminus \mathcal{S}_0
$$

and $\mathcal{S}_* = \overline{\mathcal{F}}_1$, where the closure is taken in the $\mathbb{R} \times \mathcal{F}(\Omega)$ topology.

Observe that the solutions in $\mathcal{S}_0$ are parameterized in $\beta$. We start with the asymptotic
analysis when $\beta \to \infty$ of the spectrum associated to solutions of type $\mathcal{S}_0$.

**Lemma 3.8.** Let $(w_1, w_2, u) \in \mathcal{S}_0$. The eigenvalues of $A_\beta$ behave like

$$
\beta \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2} \sim \frac{\lambda k - \mu \omega}{\mu}, \quad \gamma_{1,\beta} \sim -\left(\lambda + \frac{\lambda k - \mu \omega}{\mu}\right), \quad \gamma_{2,\beta} \to 0^{-} \quad \text{as } \beta \to \infty.
$$

As a consequence, the supremum of the spectrum of the matrix $A_\beta$ is described, in terms of $\beta$, by the
curve

$$
\beta \mapsto \beta \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2}.
$$
Moreover, this supremum of the spectrum is monotone increasing in $\beta$ and its limit as $\beta \to \infty$ can be made arbitrarily large by taking $\mu$ small accordingly. In particular, in the limit case $\mu = 0$, the spectrum is unbounded.

Lastly, the unstable direction of $A_\beta$ is spanned by the eigenvector $(1, -1, 0)$.

We can now derive a result concerning the existence of non constant solutions: our construction is implicit and uses the topological degree argument through a bifurcation analysis of the set of constant solutions. Let $0 = \gamma_0 < \gamma_1 \le \gamma_2 \le \dots$ denoted the unbounded sequence of eigenvalues of the Laplace operator with homogeneous Neumann boundary conditions and let $\{\psi_i\}$ be the corresponding eigenfunctions:

$$ (3.3) \qquad \begin{cases} -\Delta\psi_i = \gamma_i\psi_i & \text{in } \Omega \\ \partial_\nu\psi_i = 0 & \text{on } \partial\Omega. \end{cases} $$

We define $n^* \in \mathbb{N}$ to be the largest index corresponding to an eigenvalue $\gamma_{n^*}$ such that

$$ \gamma_{n^*} < \frac{\lambda k - \mu\omega}{\mu}. $$

We assume in the following that $n^* \ge 1$. Let also $\beta_n > 0$ be defined by

$$ \beta_n \frac{\lambda k - \mu\omega}{\mu\beta_n + 2k^2} = \gamma_n. $$

Observe that $n^*$ can be made as large as desired by taking $\mu$ small accordingly.

**Theorem 3.9.** For any $1 \le n \le n^*$, if $\gamma_n$ is an eigenvalue of odd multiplicity, then $(\beta_n, w_{1,n}, w_{2,n}, u_n)$ is a bifurcation point from the branch of solutions $\mathcal{S}_0$ into non constant solutions $\mathcal{S}_*$. More precisely, there exists a maximal closed and connected subset $\mathcal{C}_n \subset \mathcal{S}_*$ of solutions of (3.1) such that $\mathcal{C}_n$ contains the point $(\beta_n, w_{1,n}, w_{2,n}, u_n)$ and either

* $\mathcal{C}_n$ is unbounded in $\beta$, or

* $\mathcal{C}_n$ contains another point $(\beta_m, w_{1,m}, w_{2,m}, u_m)$ for a different value of $m$.

Aside from these bifurcation points emanating from the branch of solutions $\mathcal{S}_0$, the set $\mathcal{C}_n$ consists of solutions $(w_1, w_2, u)$ which are non constant.
**Remark 3.10.** One could wonder what happens for the eigenvalue $\gamma_0 = 0$. Actually, this is already contained in the previous remarks: indeed, $\gamma_0$ corresponds to the value $\beta = 0$, and we have already observed in Lemma 3.2 that in this situation there exists a one dimensional subspace of constant solutions emanating from this point. Thus it is also a bifurcation point from $\mathcal{S}_0$ (in such case, the branch is explicit and the solutions are constant). This point is particular and indeed we do not take it into account in the statement of Theorem 3.9.

**Remark 3.11.** If one assumes some symmetry properties for the domain $\Omega$, one can then also give a more detailed description of the branches in Theorem 3.9. In particular we can show that the symmetries of the eigenfunctions are preserved along a global branch of solutions, see for instance [20, 21].

*Proof.* The theorem follows from the classical bifurcation theorem of Rabinowitz, see [18, 19]. For $\beta > 0$ and a corresponding nontrivial constant solution $(w_1, w_2, u)$ with $w_1 = w_2$, we look for a new solution of the form $(w_1 + \varphi_1, w_2 + \varphi_2, u + \varphi)$, for small perturbations $\varphi = (\varphi_1, \varphi_2, \varphi) \in \mathcal{F}(\Omega)$. Inserting this ansatz in the system (3.1) we obtain

$$ (3.4) \qquad - \Delta \varphi = A_{\beta} \varphi + \begin{pmatrix} k\varphi_1\varphi - \beta\varphi_1\varphi_2 \\ k\varphi_2\varphi - \beta\varphi_1\varphi_2 \\ -\mu\varphi^2 - k(\varphi_1 + \varphi_2)\varphi \end{pmatrix} = A_{\beta} \varphi + H(\beta, \varphi) \quad \text{in } \Omega $$

completed by homogeneous Neumann boundary conditions. Here the nonlinear functional $H : (\mathbb{R}, \mathcal{F}(\Omega)) \to \mathcal{F}(\Omega)$ is continuous and $\|H(\beta, \varphi)\|_{\mathcal{F}(\Omega)} \le C\|\varphi\|_{\mathcal{F}(\Omega)}^2$ for a constant $C > 0$ that can be chosen uniformly on compact sets of $\beta \in [0, +\infty)$. Let us now introduce the operator $L \in \mathcal{K}(\mathcal{F}(\Omega); \mathcal{F}(\Omega))$ defined as the linear map such that for any $\mathbf{u}, \mathbf{f} \in \mathcal{F}(\Omega)$

$$ \mathbf{u} = L\mathbf{f} \iff \begin{cases} -\Delta\mathbf{u} + \mathbf{u} = \mathbf{f} & \text{in } \Omega \\ \partial_{\nu}\mathbf{u} = 0 & \text{on } \Omega. \end{cases} $$

We can rewrite the perturbed system as

$$ (3.5) \qquad \varphi = (A_{\beta} + \text{Id})L\varphi + LH(\beta, \varphi) = (A_{\beta} + \text{Id})L\varphi + h(\beta, \varphi) $$

where now $h : (\mathbb{R}, \mathcal{F}(\Omega)) \to \mathcal{F}(\Omega)$ is a compact operator. Furthermore, it is such that $\|h(\beta, \varphi)\|_{\mathcal{F}(\Omega)} \le C\|\varphi\|_{\mathcal{F}(\Omega)}^2$ with a constant $C > 0$ that again can be chosen uniformly on
compact sets of $\beta$. We are now in a position to apply the global bifurcation theorem of Rabinowitz [18, 19]. Indeed, as $\beta$ varies, $(A_\beta + \text{Id})L$ is a homotopy of compact operators. It is known that a sufficient condition for a value $\bar{\beta}$ to be a bifurcation point for the equation (3.5) is that the set of solutions to the linear equation

$$ \varphi = (A_{\bar{\beta}} + \text{Id})L\varphi $$

has odd dimension. This equation translates into the 3-component system

$$ \begin{cases} -\Delta\varphi = A_{\bar{\beta}}\varphi & \text{in } \Omega \\ \partial_v\varphi = 0 & \text{on } \Omega. \end{cases} $$

We have already studied the spectral properties of the matrix $A_{\bar{\beta}}$ in Lemma 3.3. The matrix has a unique positive eigenvalue $\bar{\beta}(\lambda k - \mu\omega)/(\mu\bar{\beta} + 2k^2)$ that correspond to the eigenvector $(-1, 1, 0)$. As a consequence, $(\gamma_i, \psi_i)$ is an eigenvalue-eigenvector couple of (3.3) and $\bar{\beta}(\lambda k - \mu\omega)/(\mu\bar{\beta} + 2k^2) = \gamma_i$ if and only if $\varphi = (\psi_i, -\psi_i, 0)$ solves the previous system for the prescribed value of $\bar{\beta}$. In particular if $\gamma_i$ has odd multiplicity, then $\bar{\beta}$ is a bifurcation point in the sense of the theorem.

It remains to show that the continua $\mathcal{C}_n$ are either unbounded in $\beta$ or meet the set $\mathcal{S}_0$ in another bifurcation point: recalling the global bifurcation theorem of Rabinowitz [18, Theorem 1.3], we already know that each continuum is either unbounded in $\mathbb{R} \times \mathcal{F}(\Omega)$ or touches the set $\mathcal{S}_0$ at another bifurcation point. Hence the statement in the theorem is reduced to showing that if the continuum $\mathcal{C}_n$ is unbounded, then it must be unbounded in $\beta$. We recall that, by Lemma 3.1, the non-negative solutions satisfy the inequalities

$$ w_1 \ge 0, \quad w_2 \ge 0, \quad 0 \le u \le \lambda/\mu, \quad u + w_1 + w_2 \le \frac{(\lambda + \omega)^2}{4\mu\omega} $$

and either all the inequalities are strict or the solution is constant. It follows that if $\beta$ is bounded on $\mathcal{C}_n$, there must exists on $\mathcal{C}_n$ a solution which is constant. In view of Lemma 3.3, discarding the solutions on $\mathcal{S}_0$, the only possibilities are solutions which are either case (a) and (b) strongly unstable or case (c) strongly stable. Recall that these properties do not depend on $\beta > 0$. We shall exclude these possibilities in the following results. $\square$

**Lemma 3.12.** The set of solutions $\mathbb{R} \times (0, 0, 0)$ is isolated in $\mathcal{P}$ for $\beta$ bounded.
**Proof.** We assume that there exists a sequence $(\beta_n, w_{1,n}, w_{2,n}, u_n) \in \mathcal{P} \setminus \mathbb{R} \times (0,0,0)$ of solutions of (3.1) such that $\beta_n > 0$ and, as $n \to +\infty$, we have $\beta_n \to \bar{\beta} \in [0, +\infty)$ and $\mathbf{v}_n := (w_{1,n}, w_{2,n}, u_n) \to (0,0,0)$ in $\mathcal{F}(\Omega)$. We can rewrite (3.1) as

$$ -\Delta \mathbf{v}_n = \begin{pmatrix} -\omega & 0 & 0 \\ 0 & -\omega & 0 \\ 0 & 0 & \lambda \end{pmatrix} \mathbf{v}_n + \begin{pmatrix} k u_n w_{1,n} - \beta_n w_{1,n} w_{2,n} \\ k w_{2,n} u_n - \beta_n w_{1,n} w_{2,n} \\ -\mu u_n^2 - k(w_{1,n} + w_{2,n}) u_n \end{pmatrix} = A_{0,\beta} \mathbf{v}_n + H_0(\beta_n, \mathbf{v}_n) $$

in $\Omega$, where $H_0 : (\mathbb{R}, \mathcal{F}(\Omega)) \to \mathcal{F}(\Omega)$ is continuous and $\|H_0(\beta, \mathbf{v})\|_{\mathcal{F}(\Omega)} \le C \|\mathbf{v}\|_{\mathcal{F}(\Omega)}^2$ locally at $\beta = \bar{\beta}$. We follow a reasoning similar to that of Theorem 3.9. We recall that the eigenvalues of the Laplacian with Neumann boundary conditions are non negative. Thus, by the stability analysis of Lemma 3.3, for $n \to +\infty$ we have that there exists $\varepsilon_n \to 0$ and an eigenpair $(\gamma_i, \psi_i)$ of (3.3) such that

$$ \lambda = \gamma_i \quad \text{and} \quad \mathbf{v}_n = \varepsilon_n \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \psi_i + o(\varepsilon_n). $$

Since $\lambda > 0$, it must be that the index $i$ is strictly positive. But then the eigenfunction $\psi_i$ changes sign in $\Omega$, and for $n$ sufficiently large, so does $u_n$. We reach the desired contradiction, as we are considering only solutions that are non negative in $\Omega$. $\square$

In Lemma 4.5 we will show that the same conclusion holds for $\beta_n \to +\infty$, that is, the sets $\mathbb{R} \times (0,0,0)$ and $\mathcal{P} \setminus \mathbb{R} \times (0,0,0)$ are at a positive distance.

We now turn to the other line of constant solutions.

**Lemma 3.13.** The set of solutions $\mathbb{R} \times (0,0,\lambda/\mu)$ is isolated in $\mathcal{P}$ for $\beta$ bounded.

**Proof.** The proof is rather similar to that of Lemma 3.12. We omit the details. We only point out that this time the conclusion is reached exploiting the expansion of the solutions as

$$ \frac{\lambda k - \mu \omega}{\mu} = \gamma_i \quad \text{and} \quad \mathbf{v}_n = \varepsilon_n \begin{pmatrix} 1 \\ 0 \\ -\frac{\lambda k}{\lambda k - \mu \omega + \lambda \mu} \end{pmatrix} \psi_i + o(\varepsilon_n) $$

and the assumption (H), that implies again $\gamma_i > 0$. $\square$
Observe that here, in contrast with the case of $\mathbb{R} \times (0,0,0)$, the sets $\mathbb{R} \times (0,0,\lambda/\mu)$ and $\mathcal{P} \setminus \mathbb{R} \times (0,0,\lambda/\mu)$ are at distance 0. This is due to the presence of the set $\mathcal{S}_0$.

We can strengthen the result in Theorem 3.9, by showing that on all branches of non constant solutions, $\beta$ is bounded away from zero. More precisely, we have

**Proposition 3.14.** There exists $\bar{\beta} > 0$ such that the set of solution of (3.1) consists only of constant solutions if $\beta \in [0, \bar{\beta})$.

*Proof.* Let us assume that there exists a sequence $(\beta_n, w_{1,n}, w_{2,n}, u_n)$ of non constant solutions of (3.1) such that $\beta_n > 0$ and $\beta_n \to 0$ as $n \to +\infty$. Up to striking out a subsequence, $v_n := (w_{1,n}, w_{2,n}, u_n)$ converges to a constant solution $\bar{v} := (\bar{w}_1, \bar{w}_2, \bar{u})$ with $\beta = 0$. We have already classified these solutions in Lemma 3.2. By the results in Lemmas 3.12 and 3.13, we know that the sequence $(w_{1,n}, w_{2,n}, u_n)$ must converge to a solution in the linear space of solution of Lemma 3.2, that is the segment

$$s \in [0, 1] \mapsto v_s = \left( \frac{\lambda k - \mu \omega}{k^2} s, \frac{\lambda k - \mu \omega}{k^2} (1-s), \frac{\omega}{k} \right).$$

We now prove that they must converge to the solution $v_{1/2}$. Indeed, assume that there exists $s^* \in [0, 1] \setminus \{1/2\}$ such that $v_n \to v_{s^*}$. Without loss of generality, we can assume that $s^* > 1/2$. Thus, for $n$ sufficiently large, we have $w_{1,n} > w_{2,n}$ in $\Omega$. We define

$$g_n = -\beta_n w_{1,n} w_{2,n} < 0 \quad \text{in } \Omega.$$

For $n$ large, we deduce from (3.1) that $w_{1,n}$ and $w_{2,n}$ are both distinct solutions of the linear equation

$$\begin{cases} -\Delta w_{i,n} + (\omega - k u_n) w_{i,n} = g_n & \text{in } \Omega \\ \partial_v w_{i,n} = 0 & \text{on } \partial \Omega. \end{cases}$$

But then, by Fredholm's alternative it must be that the difference of any two distinct solutions is orthogonal to the zero order term $g_n$, that is

$$0 = \int_{\Omega} g_n (w_{1,n} - w_{2,n}) = -\beta_n \int_{\Omega} w_{1,n} w_{2,n} (w_{1,n} - w_{2,n}) < 0$$

an obvious contradiction.
As a result, we have $\mathbf{v}_n \to \mathbf{v}_{1/2}$ in $\mathcal{F}(\Omega)$. Thus, up to a subsequence, we can write,

$$ \mathbf{v}_n = \left( \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2}, \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2}, \frac{\lambda \beta_n + 2k \omega}{\mu \beta_n + 2k^2} \right) + \boldsymbol{\varphi}_n $$

where $\boldsymbol{\varphi}_n = (\varphi_{1,n}, \varphi_{2,n}, \varphi_n)$ is such that $\boldsymbol{\varphi}_n \to 0$ in $\mathcal{F}(\Omega)$ as $\beta_n \to 0$. We let

$$ \boldsymbol{\varphi}_n = \bar{\boldsymbol{\varphi}}_n + \hat{\boldsymbol{\varphi}}_n \quad \text{where} \quad \bar{\boldsymbol{\varphi}}_n = \frac{1}{|\Omega|} \int_{\Omega} \boldsymbol{\varphi}_n $$

Plugging these relations in system (3.1), we find

$$ \begin{cases} -\Delta \hat{\boldsymbol{\varphi}}_n = A_{\beta_n} \hat{\boldsymbol{\varphi}}_n + A_{\beta} \bar{\boldsymbol{\varphi}}_n + H(\beta_n, \boldsymbol{\varphi}_n) & \text{in } \Omega \\ \partial \hat{\boldsymbol{\varphi}}_n = 0 & \text{on } \partial \Omega \end{cases} $$

where we used notations similar to those in the proof of Theorem 3.9. In particular, we have $\|H(\beta_n, \boldsymbol{\varphi}_n)\|_{\mathcal{F}(\Omega)} \le C \|\boldsymbol{\varphi}_n\|_{\mathcal{F}(\Omega)}^2$ for a positive constant C. Let us now test the equation against $\hat{\boldsymbol{\varphi}}_n$. We find

$$ (3.6) \qquad \int_{\Omega} |\nabla \hat{\boldsymbol{\varphi}}_n|^2 - \int_{\Omega} \langle A_{\beta_n} \hat{\boldsymbol{\varphi}}_n, \hat{\boldsymbol{\varphi}}_n \rangle = \int_{\Omega} H(\beta_n, \boldsymbol{\varphi}_n) \hat{\boldsymbol{\varphi}}_n $$

We now derive estimates for the two sides of the previous equation. First, recalling the definition of the function $H$ in (3.4), we see that the right hand side of (3.6) is bounded from above by

$$ \left| \int_{\Omega} H(\beta_n, \boldsymbol{\varphi}_n) \hat{\boldsymbol{\varphi}}_n \right| \le C \|\boldsymbol{\varphi}_n\|_{\mathcal{F}(\Omega)} \int_{\Omega} |\hat{\boldsymbol{\varphi}}_n|^2 $$

for some positive constant C > 0. On the other hand, by the results in Lemma 3.3 and 3.8 we find

$$ \begin{aligned} \int_{\Omega} \langle A_{\beta_n} \hat{\boldsymbol{\varphi}}_n, \hat{\boldsymbol{\varphi}}_n \rangle = & \beta_n \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2} \int_{\Omega} |(1, -1, 0) \cdot \hat{\boldsymbol{\varphi}}_n|^2 \\ & + \operatorname{Re}(\gamma_{1,n}) \int_{\Omega} |e_{1,n} \cdot \hat{\boldsymbol{\varphi}}_n|^2 + \operatorname{Re}(\gamma_{2,n}) \int_{\Omega} |e_{2,n} \cdot \hat{\boldsymbol{\varphi}}_n|^2 \end{aligned} $$

where $\gamma_{1,n}$ and $\gamma_{2,n}$ are two complex conjugate eigenvalues of $A_{\beta_n}$ with strictly negative real part and $e_{1,n}$ and $e_{2,n}$ are the corresponding eigenvectors. More precisely the following holds

$$ \beta_n \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2} = o_n(1) \quad \text{and} \quad \operatorname{Re}(\gamma_{i,n}) = -\frac{\mu \omega}{2k} + o_n(1) < 0. $$
Since by construction the $\hat{\phi}_n$ have zero average on $\Omega$, by Poincaré's inequality, we find the following lower bound for the left hand side of (3.6)

$$ \int_{\Omega} |\nabla \hat{\phi}_n|^2 - \int_{\Omega} \langle A_{\beta_n} \hat{\phi}_n, \hat{\phi}_n \rangle \geq C \int_{\Omega} |\hat{\phi}_n|^2. $$

Here, the constant $C > 0$ can be chosen independent on $n$ for $n$ large enough.

Combining the two bounds, we infer from (3.6) that

$$ \int_{\Omega} |\hat{\phi}_n|^2 \leq C \|\varphi_n\|_{\mathcal{F}(\Omega)} \int_{\Omega} |\hat{\phi}_n|^2 $$

for $C > 0$, a contradiction for $n$ sufficiently large.

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

Thus we also conclude the proof of Theorem 3.9.

We can go further and make the conclusion of Theorem 3.9 more precise by noticing that we are in a position to apply the analytic bifurcation theory developed by Dancer in [11, 12] (see also [5, Theorem 9.1.1]). This approach yields the following.

**Theorem 3.15.** Under the assumptions of Theorem 3.9, for any continua of solutions $\mathcal{C}_n$ there exists a curve $\mathcal{C}_n := \{(\beta(s), v(s)) : \mathbb{R} \mapsto \mathbb{R}^+ \times \mathcal{F}(\Omega)\} \subset \mathcal{C}_n$, which contains the bifurcation point from which $\mathcal{C}_n$ emanate, such that

* at any point, the curve $\mathcal{C}_n$ can be locally reparametrized as an analytic curve;

* the set of possible secondary bifurcation points on $\mathcal{C}_n$ has no accumulation points.

Moreover

* either $\mathcal{C}_n$ is a closed loop, and meets the set $\mathcal{S}_0$ in two distinct bifurcation points;

* or the set $\mathcal{C}_n$ is unbounded in $(\bar{\beta}, +\infty) \times \mathcal{F}(\Omega)$, and more specifically

$$ \beta(s) \to +\infty \quad \text{as } s \to \infty. $$

## 4. A PRIORI ESTIMATES AND STRONG COMPETITION SINGULAR LIMITS

To shed more light on the solutions, we now derive new regularity estimates on the solutions. We are chiefly interested in estimates that are uniform in the competition parameter $\beta \geq 0$. Here we consider non negative solutions, that is solutions $v = (w_1, \dots, w_N, u)$ with $w_i \geq 0$ for all $i$ and $u \geq 0$. Let us recall that, if $\beta$ is bounded, the regularity result of Lemma 3.1 applies also two the system of $N$ (+1) components.
**Proposition 4.1.** Let $\Omega \subset \mathbb{R}^N$ be a smooth domain, and let $\beta$, $D$, $d_i$, $\omega_i$, $k_i$, $a_{ij} = a_{ji}$ for $1 \le i, j \le k$ be positive parameters. We consider a non negative solution $\mathbf{v} = (w_1, \dots, w_N, u) \in \mathcal{F}(\Omega)$ of the system

$$ (4.1) \qquad \left\{ \begin{array}{ll} -d_i \Delta w_i = (-\omega_i + k_i u - \beta \sum_{j \neq i} a_{ij} w_j) w_i & \\ -D \Delta u = (\lambda - \mu u - \sum_{i=1}^N k_i w_i) u & \\ \partial_v w_i = \partial_v u = 0 & \text{on } \partial \Omega. \end{array} \right. $$

Then all components of $\mathbf{v}$ are uniformly bounded in $L^\infty(\Omega)$ with respect to $\beta > 0$ and moreover there exists $C$ (independent of $\beta$ and $N$) such that

$$ \|(w_1, \dots, w_N)\|_{Lip(\bar{\Omega})} + \|u\|_{C^{2,\alpha}(\bar{\Omega})} \le C $$

If $\{\mathbf{v}_\beta\}_\beta$ is a family of non negative solutions as above, defined for $\beta \to +\infty$, then, up to subsequences, there exists $\mathbf{v} = (w_1, \dots, w_N, u)$ with $(w_1, \dots, w_N) \in Lip(\bar{\Omega})$ and $u \in C^{2,\alpha}(\bar{\Omega})$ and

$$ (4.2) \qquad (w_{1,\beta}, \dots, w_{N,\beta}) \to (w_1, \dots, w_N) \text{ in } C^{0,\alpha} \cap H^1(\bar{\Omega}), u_\beta \to u \text{ in } C^{2,\alpha}(\bar{\Omega}) $$

for any $\alpha \in (0,1)$. Any such limit satisfies in the sense of measures the following system of inequalities

$$ (4.3) \qquad \left\{ \begin{array}{ll} -d_i \Delta w_i \le (-\omega_i + k_i u - \mu_i) w_i & \\ -\Delta \left( d_i w_i - \sum_{j \neq i} d_j w_j \right) \ge (-\omega_i + k_i u) w_i - \sum_{j \neq i} (-\omega_j + k_j u) w_j & \text{in } \Omega \\ -D \Delta u = \left( \lambda - \mu u - \sum_{i=1}^k k_i w_i \right) u & \\ \partial_v w_i = \partial_v u = 0 & \text{on } \partial \Omega. \end{array} \right. $$

Lastly, if the limit has two or more non zero components of $(w_1, \dots, w_N)$, then the subset $\{x \in \Omega : \sum_{i=1}^N w_i = 0\}$ is a rectifiable set of co-dimension 1, made of the union of a finite number of $C^{1,\alpha}$ smooth sub-manifolds.

**Remark 4.2.** It can be shown that the limit in (4.2) does not hold in the Lipschitz norm (the case $\alpha = 1$), although the sequence is bounded in the Lipschitz norm. This follows from Hopf's Lemma applied to the limit functions on a regular portion of the free boundary $\{x \in \Omega : \sum_{i=1}^N w_i = 0\}$ in conjunction with the obvious fact that $C^1(\Omega)$ is a closed subset of $Lip(\Omega)$.
**Remark 4.3.** From the theorem we can deduce the following. Let $M$ be the number of non zero components of $(w_1, \dots, w_N)$. If $M = 0$, then either $u = 0$ or $u = \lambda/\mu$. If $M = 1$, by [17, Theorem 1] the function $u$ and the unique non zero component of $(w_1, \dots, w_N)$ are positive constants and $\{x \in \Omega : \sum_{i=1}^N w_i = 0\} = \emptyset$.

We shall not prove the result in all of its details, since it follows from already known ones. See for instance [7] adn [6] for the uniform estimates in Hölder spaces, [22] for the uniform estimate in the optimal Lipschitz norm, [8], [6] and [23] for the study of a closely related free-boundary problem. In any case, a complete proof for a much strong result can be found in [1].

We observe that, in order to simplify the exposition, we have assumed $a_{ij} = 1$ for all $1 \le i, j \le k$: the results that follow can be generalized without difficulty. A much harder case arises when the competition matrix is not symmetric, that is when $a_{ij} \neq a_{ji}$ for some $i \neq j$: even though most of the results are also valid in this case, we will not consider it here, since we can only obtain a less complete description of the solutions. We refer the reader to [24] to understand the new difficulties of this case.

We first use the uniform estimates to study more closely those bifurcation branches $\mathcal{C}_n$ which are unbounded in $\mathbb{R}^+ \times \mathcal{F}(\Omega)$. More generally, we will look here at any sequence of solutions $(\beta_n, w_{1,n}, w_{2,n}, u_n)$ such that $\beta_n \to +\infty$. We recall that assumptions (H) holds, in particular $\lambda k > \mu\omega$. In the analysis of the singular limit, we use a blow-up technique first introduced in [13] to study a similar situation.

We start with a key property.

**Lemma 4.4.** There exists $M > 0$ such that

$$ \frac{1}{M} \|w_{2,\beta}\|_{L^{\infty}(\Omega)} \le \|w_{1,\beta}\|_{L^{\infty}(\Omega)} \le M \|w_{2,\beta}\|_{L^{\infty}(\Omega)} $$

for all $v_\beta \in \mathcal{P}$ and $\beta$ sufficiently large.

*Proof.* We argue by contradiction, assuming that there exists a sequence of solutions in $\mathcal{P}$ that invalidates the conclusion. Without loss of generality, let us assume that $\|w_{1,n}\|_{L^{\infty}(\Omega)} \le \|w_{2,n}\|_{L^{\infty}(\Omega)}$ and that the ratio $\|w_{1,n}\|_{L^{\infty}(\Omega)} / \|w_{2,n}\|_{L^{\infty}(\Omega)} \to 0$ as $\beta_n \to +\infty$. We introduce the
renormalized functions

$$ \bar{w}_{i,n} = \frac{w_{i,n}}{\|w_{i,n}\|_{L^{\infty}(\Omega)}} \quad \text{for } i = 1,2 $$

which are solutions to

$$ \left\{ \begin{array}{ll} -\Delta \bar{w}_{1,n} = -\omega \bar{w}_{1,n} + k \bar{w}_{1,n} u - \beta_n \|w_{2,n}\|_{L^{\infty}(\Omega)} \bar{w}_{1,n} \bar{w}_{2,n} & \text{in } \Omega \\ -\Delta \bar{w}_{2,n} = -\omega \bar{w}_{2,n} + k \bar{w}_{2,n} u - \beta_n \|w_{1,n}\|_{L^{\infty}(\Omega)} \bar{w}_{1,n} \bar{w}_{2,n} & \text{in } \Omega \\ -\Delta u_n = \lambda u_n - \mu u_n^2 - k(w_{1,n} + w_{2,n}) u_n & \text{in } \Omega \\ \partial_v \bar{w}_{i,n} = \partial_v u = 0 & \text{on } \partial \Omega \end{array} \right. $$

We distinguish between two different cases.

1) $\beta_n \|w_{2,n}\|_{L^\infty(\Omega)}$ is bounded. In this case, all the terms in the equations are bounded uniformly with respect to $n$, and thus it is easy to see that the sequence $\bar{w}_{i,n}$, $u_n$ and also $w_{i,n}$, are uniformly bounded in $W^{2,p}(\Omega)$ for any $p < \infty$. Up to striking out a subsequence, we derive the strong convergence of the renormalized densities to some limit profile $(\bar{w}_{1,\infty}, \bar{w}_{2,\infty}, u_\infty)$ with both $\bar{w}_{1,\infty}$ and $\bar{w}_{2,\infty}$ positive, while by assumption $w_{i,n} \to 0$ uniformly in $\Omega$. Moreover, by assumption we have that

$$ \beta_n \|w_{2,n}\|_{L^\infty(\Omega)} \to C \geq 0 \quad \text{while} \quad \beta_n \|w_{1,n}\|_{L^\infty(\Omega)} \to 0. $$

As a result, the limit profiles solve

$$ \left\{ \begin{array}{ll} -\Delta \bar{w}_{1,\infty} = -\omega \bar{w}_{1,\infty} + k \bar{w}_{1,\infty} u - C \bar{w}_{1,\infty} \bar{w}_{2,\infty} & \text{in } \Omega \\ -\Delta \bar{w}_{2,\infty} = -\omega \bar{w}_{2,\infty} + k \bar{w}_{2,\infty} u & \text{in } \Omega \\ -\Delta u_\infty = \lambda u_\infty - \mu u_\infty^2 & \text{in } \Omega \\ \partial_v \bar{w}_{i,\infty} = \partial_v u = 0 & \text{on } \partial \Omega. \end{array} \right. $$

By the maximum principle (see also Lemma 3.6), we know that the equation for $u_\infty$ has only the constant solutions $u_\infty = 0$ or $\lambda/\mu$. Inserting this information in the equation satisfied by $\bar{w}_{2,\infty}$ we see that

$$ \left\{ \begin{array}{ll} -\Delta \bar{w}_{2,\infty} = -\omega \bar{w}_{2,\infty} + k \bar{w}_{2,\infty} u_\infty = C' w_{2,\infty} & \text{in } \Omega \\ \partial_v \bar{w}_{2,\infty} = 0 & \text{on } \partial \Omega. \end{array} \right. $$
where the constant $C'$ is non zero by the assumption (H). It follows that necessarily $\bar{w}_{2,\infty} \equiv 0$, in contradiction with $\|\bar{w}_{2,\infty}\|_{L^\infty(\Omega)} = 1$.

2) $\beta_n \|w_{2,n}\|_{L^\infty(\Omega)} \to +\infty$, along a subsequence. We test the equation in $\bar{w}_{i,n}$ by $\bar{w}_{i,n}$ itself. Recalling that $\bar{w}_{i,n} \ge 0$ and that $u_n \le \lambda/\mu$ we have

$$ \int_{\Omega} |\nabla \bar{w}_{i,n}|^2 \le k \frac{\lambda}{\mu} |\Omega|, $$

where $|\Omega|$ is the measure of the set $\Omega$. Consequently, the $\bar{w}_{i,n}$ are bounded uniformly in $H^1(\Omega)$ and thus converge weakly to some limit $\bar{w}_{i,\infty} \in H^1(\Omega)$. Moreover, the compact embedding of $H^1(\Omega)$ in $L^2(\Omega)$ yields $\bar{w}_{i,n} \to \bar{w}_{i,\infty}$ strongly in $L^2(\Omega)$ and furthermore, since by construction $\|w_{i,n}\|_{L^\infty(\Omega)} = 1$, we have that $\bar{w}_{i,n} \to \bar{w}_{i,\infty}$ strongly in $L^p(\Omega)$ for any $p \ge 2$. Recalling that the equation for $u_n$ contains only uniformly bounded terms, up to a subsequence we have $u_n \to u_\infty$ in $W^{2,p}(\Omega)$ for any $p < \infty$. Let us show that each component of the limit configuration $(\bar{w}_{1,\infty}, \bar{w}_{2,\infty}, u_\infty)$ is non zero. From the equations satisfied by $\bar{w}_{i,n}$ we know that

$$ \begin{cases} -\Delta \bar{w}_{i,n} + \omega \bar{w}_{i,n} \le k \bar{w}_{i,n} u_n & \text{in } \Omega \\ \partial_v \bar{w}_{i,n} = 0 & \text{on } \partial \Omega. \end{cases} $$

Now, letting $g_{i,n} \in H^1(\Omega)$ be the solution of

$$ \begin{cases} -\Delta g_{i,n} + \omega g_{i,n} = k \bar{w}_{i,n} u_n & \text{in } \Omega \\ \partial_v g_{i,n} = 0 & \text{on } \partial \Omega \end{cases} $$

we have, from the previous discussion, that the sequence $\{g_{i,n}\}_n$ is compact in $W^{2,p}(\Omega)$ for any $p > 1$ and, in particular, in $C^{0,\alpha}(\Omega)$ for some $\alpha > 0$. On the other hand, the maximum principle yields $0 \le \bar{w}_{i,n} \le g_{i,n}$. We then assume, by way of contradiction, that either $\bar{w}_{i,\infty} = 0$ or $u_\infty = 0$. Then it follows that $g_{i,n} \to 0$ uniformly, that is, $\bar{w}_{i,n} \to 0$ uniformly. This is in contradiction with $\|w_{i,n}\|_{L^\infty(\Omega)} = 1$.

Testing the equations in $\bar{w}_{1,n}$ by $\varphi \in H^1(\Omega)$, we get

$$ (4.4) \quad \beta_n \|w_{2,n}\|_{L^\infty(\Omega)} \int_{\Omega} \bar{w}_{1,n} \bar{w}_{2,n} \varphi = \int_{\Omega} (ku_n - \omega) \bar{w}_{1,n} \varphi - \int_{\Omega} \nabla w_{1,n} \cdot \nabla \varphi \le C $$
so that, using our assumption

$$ \beta_n \|w_{1,n}\|_{L^\infty(\Omega)} \int_{\Omega} \bar{w}_{1,n} \bar{w}_{2,n} \varphi = \frac{\|w_{1,n}\|_{L^\infty(\Omega)}}{\|w_{2,n}\|_{L^\infty(\Omega)}} \cdot \beta_n \|w_{2,n}\|_{L^\infty(\Omega)} \int_{\Omega} \bar{w}_{1,n} \bar{w}_{2,n} \varphi \to 0. $$

As a result, $\bar{w}_{2,\infty}$ is a weak solution of the equation

$$ \begin{cases} -\Delta \bar{w}_{2,\infty} = -\omega \bar{w}_{2,\infty} + k \bar{w}_{2,\infty} u_{\infty} & \text{in } \Omega \\ \partial_{\nu} \bar{w}_{2,\infty} = 0 & \text{on } \partial \Omega \end{cases} $$

where $0 \le u_{\infty} \le \lambda/\mu$. By the maximum it follows that either $\bar{w}_{2,\infty} \equiv 0$ or $\bar{w}_{2,\infty}$ is bounded away from 0. The former case was already excluded, thus the latter holds. But then equation (4.4), with $\varphi = 1$, yields

$$ \beta_n \|w_{2,n}\|_{L^\infty(\Omega)} \cdot \int_{\Omega} \bar{w}_{1,n} \bar{w}_{2,n} \le C \implies \int_{\Omega} \bar{w}_{1,\infty} \bar{w}_{2,\infty} = 0 $$

which implies $\bar{w}_{1,\infty} = 0$, in contradiction with the previous discussion. $\square$

To push forward our analysis, in the remaining of this section we now impose a strengthening of assumption (H). We require that $(\lambda k - \mu\omega)/k$ is not an eigenvalue of the Laplace operator with Neumann boundary conditions.

**Lemma 4.5.** The set $\mathcal{P}$ is a pre-compact subset of $C^{0,\alpha} \times C^{0,\alpha} \times C^{2,\alpha}(\bar{\Omega})$ for any $\alpha \in (0,1)$. Moreover any converging subsequence $(w_{1,n}, w_{2,n}, u_n) \to (w_{1,\infty}, w_{2,\infty}, u_\infty)$ with $\beta_n \to +\infty$ is such that

* either $(w_{1,\infty}, w_{2,\infty}, u_\infty)$ has all non zero components and, letting $V = w_{1,\infty} - w_{2,\infty}$, $V$ changes sign and $(V, u_\infty) \in C^{2,\alpha}(\bar{\Omega})$ is a non-trivial solution of

$$ \begin{cases} -\Delta V = -\omega V + kVu & \text{in } \Omega \\ -\Delta u = \lambda u - \mu u^2 - k|V|u & \text{in } \Omega \\ \partial_{\nu} V = \partial_{\nu} u = 0 & \text{on } \partial \Omega, \end{cases} $$

* or

$$ (\beta_n w_{1,n}, \beta_n w_{2,n}, u_n) \sim \left( \frac{\lambda k - \mu\omega}{\mu}, \frac{\lambda k - \mu\omega}{\mu}, \frac{\lambda}{\mu} \right) $$

as $n \to +\infty$ in $L^p(\Omega)$ for any $p < \infty$ and weakly in $H^1(\Omega)$.
Proof. The compactness in strong topology of the sequence of solutions was already established in Proposition 4.1. We are left with the study of the asymptotic profiles. First of all we exclude the case $u_n \to 0$ (which would hold uniformly in $\Omega$ by the compactness properties). Indeed, in this situation we would have

$$ egin{cases} -\Delta w_{i,n} = -\omega w_{i,n} + k w_{i,n} u_n - \beta_n w_{i,n} w_{j,n} \le -\frac{\omega}{2} w_{i,n} & \text{in } \Omega \\ \partial_v w_{i,n} = 0 & \text{on } \partial \Omega \end{cases} $$

for $n$ sufficiently large, which implies that necessarily $w_{i,n} \equiv 0$ for $n$ large, in contradiction with the assumptions.

Let us now assume that

$$ w_{1,n}, w_{2,n} \to 0 \quad \text{uniformly in } \Omega. $$

Passing to the limit in the equation in $u_n$, we see that $u_\infty$ satisfies

$$ egin{cases} -\Delta u_\infty = \lambda u_\infty - \mu u_\infty^2 & \text{in } \Omega \\ \partial_v u_\infty = 0 & \text{on } \partial \Omega \end{cases} $$

which implies that $u_n \to \lambda/\mu$ in $C^{2,\alpha}(\Omega)$ (recall that we have already excluded the case $u_n \to 0$). We introduce the renormalized functions

$$ \bar{w}_{i,n} := \frac{w_{i,n}}{\|w_{1,n}\|_{L^\infty(\Omega)}} $$

which are solutions to

$$ egin{cases} -\Delta \bar{w}_{1,n} = -\omega \bar{w}_{1,n} + k \bar{w}_{1,n} u_n - \beta_n \|w_{1,n}\|_{L^\infty(\Omega)} \bar{w}_{1,n} \bar{w}_{2,n} & \text{in } \Omega \\ -\Delta \bar{w}_{2,n} = -\omega \bar{w}_{2,n} + k \bar{w}_{2,n} u_n - \beta_n \|w_{1,n}\|_{L^\infty(\Omega)} \bar{w}_{1,n} \bar{w}_{2,n} & \text{in } \Omega \\ \partial_v \bar{w}_{i,n} = 0 & \text{on } \partial \Omega. \end{cases} $$

Let us observe that, thanks to Lemma 4.4, we know that the $\bar{w}_{i,n}$ are bounded by some positive constant $M > 0$. By redefining $\beta_n$ as $\beta_n \|w_{1,n}\|_{L^\infty(\Omega)}$ and exploiting the uniform estimates of Proposition 4.1, we have that the sequence $(\bar{w}_{1,n}, \bar{w}_{2,n})$ are precompact in $C^{0,\alpha}(\bar{\Omega})$. Moreover, using the same initial steps as in Case 2) of Lemma 4.4, we see that $\bar{w}_{i,n} \to \bar{w}_{i,\infty}$ in $L^p(\Omega)$ for
any $p < \infty$ and weakly in $H^1(\Omega)$ and also $\bar{w}_{i,\infty} \neq 0$. Letting $V_n = \bar{w}_{1,n} - \bar{w}_{2,n}$, we see that
$$ \|V_n\|_{L^\infty(\Omega)} \leq M + 1 \text{ and} $$

$$
\begin{cases}
-\Delta V_n = (-\omega + ku_n)V_n & \text{in } \Omega \\
\partial_v V_n = 0 & \text{on } \partial\Omega.
\end{cases}
$$

As a consequence of the strong convergence $u_n \to \lambda/\mu$, we see that $V_n \to V_\infty$ in $C^{2,\alpha}(\Omega)$. This function $V_\infty$ is a solution of the following limit equation

$$
\begin{cases}
-\Delta V_{\infty} = \frac{\lambda k - \mu \omega}{\mu} V_{\infty} & \text{in } \Omega \\
\partial_{\nu} V_{\infty} = 0 & \text{on } \partial \Omega
\end{cases}
$$

where, by assumption, $(\lambda k - \mu \omega)/\mu \neq \gamma_i$, the eigenvalues of the Laplacian with Neumann boundary conditions. Consequently, $V_\infty \equiv 0$ and thus $\bar{w}_{1,\infty} = \bar{w}_{2,\infty} \neq 0$. Testing the equation in $\bar{w}_{i,n}$ by $\bar{w}_{i,n}$ itself, we find

$$
\int_{\Omega} |\nabla \bar{w}_{i,n}|^2 + \beta_n \|w_{1,n}\|_{L^{\infty}(\Omega)} \int_{\Omega} \bar{w}_{i,n}^2 \bar{w}_{j,n} = \int_{\Omega} (-\omega + ku_n) \bar{w}_{i,n}^2 \leq C
$$

which implies, in particular, that $\beta_n \|w_{1,n}\|_{L^\infty(\Omega)} \to C$ for some constant $C \ge 0$. Finally, passing to the limit in the equation in $\bar{w}_{i,n}$ we find

$$
\begin{cases}
-\Delta \bar{w}_{\infty} = \left(-\omega + k \frac{\lambda}{\mu}\right) \bar{w}_{\infty} - C \bar{w}_{\infty}^2 & \text{in } \Omega \\
\partial_{\nu} \bar{w}_{\infty} = 0 & \text{on } \partial \Omega.
\end{cases}
$$

If $C = 0$, since $(\lambda k - \mu \omega)/\mu > 0$ and $\bar{w}_\infty$ is non negative by the maximum principle, it must be the case that $\bar{w}_\infty \equiv 0$, in contradiction with the renormalization. Thus $C > 0$ and, from a direct application of the maximum principle (see Lemma 3.6), we have that the only non negative solution to the previous equation are the constant. In particular, it must be the case that $\bar{w}_\infty \equiv 1$, thus

$$
C = \frac{\lambda k - \mu \omega}{\mu}. \quad \square
$$

5. ONE DIMENSIONAL CASE

In the one-dimensional case the description of the bifurcation branches is more complete.
**Theorem 5.1.** Under the assumptions of Theorem 3.9, let us moreover suppose that $\Omega \subset \mathbb{R}$ is an open and bounded interval. Then any eigenvalue $\gamma$ of (3.3) is of multiplicity one, and the corresponding continuum of solutions $\mathcal{C}_n$ (and $\mathcal{C}_n$) generated from the set $\mathcal{S}_0$ at the value

$$ \beta_n \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2} = \gamma_n $$

is unbounded and it intersects the set $\mathcal{S}_0$ only once.

*Proof.* The proof follows again the main ideas presented in [18]. Under the assumptions, there exist $a < b \in \mathbb{R}$ such that $\Omega = (a, b) \subset \mathbb{R}$. We can explicitly compute the eigenvalues of (3.3), which are given by

$$ \gamma_n := \left( \frac{\pi}{b-a} n \right)^2 \quad \text{for any } n \in \mathbb{N}. $$

Based on the discussion in Theorem 3.9, any value of $\gamma_n$ corresponds to a bifurcation point. This is true even for the set $\mathcal{C}_0$ generating from $\gamma_0 = 0$, which is given by a trivial linear subspace of constant solutions, subject of Lemma 3.2.

Let us consider, for a fixed $\gamma_n$ with $n \ge 1$ as before, the continuum of solutions $\mathcal{C}_n$ that emerges from the set $\mathcal{S}_0$. By the perturbations analysis conducted in Theorem 3.9, we know that the solutions are of the form

$$ (w_{1,\beta}, w_{2,\beta}, u_{\beta}) = \left( \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2}, \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2}, \frac{\lambda \beta_n + 2k \omega}{\mu \beta_n + 2k^2} \right) + \varepsilon(\psi_n, -\psi_n, 0) + o(\varepsilon) $$

where $\varepsilon$ is a parameter such that $\varepsilon \to 0$ when $\beta \to \beta_n$, $\psi_n$ is a normalized eigenfunction of 3.3 in $\Omega = (a, b)$ and $o(\varepsilon)$ is a perturbation in $C^{2,\alpha}([a, b])$ of order less than $\varepsilon$. In particular, letting

$$ v_{\beta,n} = w_{1,\beta} - w_{2,\beta} = 2\varepsilon\psi_n + o(\varepsilon) $$

(where we have explicitly stated the index $n$ of the eigenfunction which spans $v_{\beta,n}$) we have that $v_{\beta,n}$ solves

$$ (5.1) \qquad \begin{cases} -v''_{\beta,n} = (-\omega + ku_{\beta})v_{\beta,n} & \text{in } (a,b), \\ v'_{\beta,n}(a) = v'_{\beta,n}(b) = 0. \end{cases} $$

As a result, when $\varepsilon$ is small, $v_\beta$ has exactly $n$ distinct simple zeroes in $(a, b)$, located closely to the zeroes of the eigenfunction $\psi_n$. We recall that the solutions of the system (3.1) are bounded
in Lip([a, b]) uniformly with respect to $\beta$ and in particular the last component $u_\beta$ is bounded in $C^{2,\alpha}([a, b])$ for all $\alpha < 1$: it follows that there exists a parametrization of the continuum $\mathscr{C}_n$ with respect to which the functions $v_{\beta,n}$ vary smoothly and they also are uniformly bounded in $C^{2,\alpha}([a, b])$ for all $\alpha < 1$.

We claim that on each continuum of solutions $\mathscr{C}_n$, the number of zeroes of the function $v_{\beta,n}$ does not change. To this end, we adapt a classical argument for scalar equations to the present situation of a system. To prove the claim, we first observe that the solutions depend smoothly in $C^{2,\alpha}([a, b])$ on the parametrization. Thus, if $v_{\beta,n}$ changes the number of zeros, there exists a solution $v$ of (5.1) inside $\mathscr{C}_n$ that has a zeros of multiplicity at least two. This point could be located at the interior or at the boundary of the interval $(a, b)$, thanks to the homogeneous Neumann condition. The uniqueness theorem for ordinary differential equations with smooth coefficients applied to (5.1) implies that the function $v$ must be identically 0. As a result, the corresponding solution $(w_1, w_2, u)$ has first and second components that are equal: reasoning as in Lemma 3.2, by letting $V = w_1 + w_2$, we obtain a solution to

$$ egin{cases} -V'' = -\omega V + kVu - \beta V^2 & \text{in } (a, b) \\ -u'' = \lambda u - \mu u^2 - kVu & \text{in } (a, b) \\ V' = u' = 0 & \text{on } \{a, b\}. \end{cases} $$

and again thanks to the results in [17], the solution $(w_1, w_2, u)$ must be a constant solution. By Proposition 3.14, we already know that $\beta > 0$. Thus $(w_1, w_2, u)$ belongs necessarily to the set $\mathscr{S}_0$ (we recall that if $\beta > 0$ the only bifurcation points belong to $\mathscr{S}_0$, see Remark 3.4). From the previous discussion it must be that the point corresponds to a different eigenvalue $\gamma_m$, $m \neq n$, and locally the solutions can be written as a perturbation along the line spanned by the eigenfunction $\psi_m$. In the same vein of the previous argument, the difference of the first two components can be asymptotically expanded as

$$ v_{\beta,m} = w_{1,\beta} - w_{2,\beta} = 2\epsilon\psi_m + o(\epsilon) $$
and, again, for $\varepsilon \to 0$, the solution $v_{\beta,m}$ has $m$ distinct simple zeroes on $(a,b)$. In particular, the solution must have $m \neq n$ zeros in a neighborhood of the bifurcation point, leading us to a contradiction. □

We can show the following stronger version of Lemma 4.5, which completes the analysis of the bifurcation diagram in dimension one.

**Lemma 5.2.** Let $\Omega \subset \mathbb{R}$ be an open and bounded interval. Then same conclusion of Lemma 4.5. Moreover, let $(w_{1,n}, w_{2,n}, u_n) \in \mathcal{P}$ be any converging sequence for $\beta_n \to +\infty$ and let $(w_{1,\infty}, w_{2,\infty}, u_\infty)$ be its limit. Then either $(w_{1,\infty}, w_{2,\infty}, u_\infty)$ has all its components non zero or for $n$ large it holds

$$ (w_{1,n}, w_{2,n}, u_n) = \left( \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2}, \frac{\lambda k - \mu \omega}{\mu \beta_n + 2k^2}, \frac{\lambda \beta_n + 2k \omega}{\mu \beta_n + 2k^2} \right). $$

*Proof.* As in Theorem 5.1, we can use the auxiliary function $v_\beta = w_{1,\beta} - w_{2,\beta}$ to study more accurately the second case of the lemma. The conclusion is reached once again by counting the number of zeroes of $v_\beta$ and observing that this must be constant on each bifurcation branch. □

As a direct consequence, we have that there exists $\delta > 0$ such that, for $\beta'$ sufficiently large,

$$ \mathrm{dist}(\mathcal{P} \setminus \mathcal{S}_0, \mathcal{S}_0) > \delta \quad \text{for all } \beta \geq \beta' $$

where the distance is taken in the sense of the $C^{0,\alpha} \times C^{0,\alpha} \times C^{2,\alpha}(\bar{\Omega})$ norm for any $\alpha \in (0,1)$. Moreover each branch of solutions constructed in Theorem 5.1 converge (up to a subsequence) to a disjoint set of solutions for the limit problem, characterized by the different number of zeroes for the function $V$ and do not re-collapse on the set $\mathcal{S}_0$ for $\beta = +\infty$.

## 6. MAXIMUM NUMBER OF DENSITIES

We continue the investigation of the model by addressing a fundamental question: can the model be used to predict the optimal repartition of the domain $\Omega$ in hunting territories, that is, the optimal number of packs?

To answer this question, we first focus on the limit stationary system satisfied by the densities in the case of segregation. We shall prove two complementary results in this direction.

(1) We show that each bounded domain $\Omega \subset \mathbb{R}^n$ can sustain a maximum number of densities of predators (see Lemma 6.1 and Theorem 6.4). As we will see, this implies
that there exists a number $k \ge 1$ of packs that, for a given configuration of parameters, maximizes the total population of predators.

(2) Then, we show that, for a particular choices of the parameters, the total population of predators in the case of two packs is strictly higher than that of one pack only, implying that in these cases the optimal configuration is given by a finite number of packs strictly greater than two.

We start with the following result, which states that for each environment $\Omega$ there is a maximal number of densities of predators $\tilde{N}$ that can be sustained.

**Lemma 6.1.** For a given smooth domain $\Omega \subset \mathbb{R}^n$, there exists $\tilde{N} \in \mathbb{N}$ such that any non negative solution $(w_1, \dots, w_N, u) \in H^1(\Omega)$ of (4.3) has at most $\tilde{N}$ non zero components among $(w_1, \dots, w_N)$.

In order to prove the previous result, we need to recall the notion of optimal partition (see for instance [16] for a general survey and some fundamental results). For consistency with the theory of optimal partitions, in the next two results, eigenvalues will be counted starting from the index 1. For any $1 \le h \in \mathbb{N}$ we say that a family $\mathcal{D} = \{D_1, \dots, D_h\}$ of open subsets of $\Omega$ is a $h$ (open) partition of $\Omega$ if

$$D_i \cap D_j = \emptyset \quad \forall i \neq j \text{ and } \overline{\cup_{i=1}^h D_i} = \overline{\Omega}.$$

For each $D_i$, we define the generalized first eigenvalue as

$$\gamma_1(D_i) := \inf \left\{ \int_{\Omega} |\nabla \varphi|^2 / \int_{\Omega} \varphi^2 : \varphi \in H^1(\Omega), \varphi = 0 \text{ in } \Omega \setminus D_i \right\}$$

and for the partition $\mathcal{D}$ we assign the total value

$$\Lambda(\mathcal{D}) = \max_i \gamma_1(D_i).$$

A partition $\mathcal{D}$ is optimal if it minimizes the value of $\Lambda(\mathcal{D})$ among all $N$-partitions. We recall the following result (see [16, Corollary 5.6]), which follows from the Courant-Fischer-Weyl characterization of the eigenvalues of compact hermitian operators.
**Theorem 6.2.** Let $\gamma_h(\Omega)$ be the h-th eigenvalue (counted with multiplicity) of

$$ \begin{cases} -\Delta\varphi = \gamma_h\varphi & \text{in } \Omega \\ \partial_v\varphi = 0 & \text{on } \partial\Omega. \end{cases} $$

Then

$$ \Lambda(\mathcal{D}) \geq \gamma_h(\Omega) $$

for all h-partitions $\mathcal{D}$ of $\Omega$.

*Proof of Lemma 6.1.* Consider the solution $(w_1, \dots, w_N, u)$ of (4.3). If the component $u$ is zero, all the remaining components of the solution are zero. This can be derived by testing the inequalities in (4.3) with $w_i \in Lip(\bar{\Omega})$. Indeed, we thus obtain

$$ \int_{\Omega} d_i |\nabla w_i|^2 + (\omega_i - k_i u) w_i^2 = 0 $$

which yields the claim taking into account that $u \equiv 0$. As a result, we can assume $u \ge 0$ and $u \ne 0$. By the maximum principle applied to the equation in $u$, we find that $u$ is strictly positive. Moreover, since $w_i \ge 0$ for $i = 1, \dots, k$, we have

$$ -D\Delta u = \left( \lambda - \mu u - \sum_{i=1}^{N} k_i w_i \right) u \le (\lambda - \mu u) u \quad \text{thus } u \le \frac{\lambda}{\mu}. $$

On the other hand, we have

$$ -d_i \Delta w_i = (-\omega_i + k_i u) w_i \le \left( -\omega_i + k_i \frac{\lambda}{\mu} \right) w_i. $$

That is, letting $\Omega_i := \{w_i > 0\}$, $w_i$ satisfies

$$ \left\{ \begin{array}{ll} -d_i \Delta w_i \le (-\omega_i + k_i \frac{\lambda}{\mu}) w_i & \text{in } \Omega_i \\ w_i = 0 & \text{on } \partial \Omega_i \cap \Omega \\ \partial_v w_i = 0 & \text{on } \partial \Omega_i \cap \partial \Omega \end{array} \right. $$
Considering the first eigenvalue $\gamma_1(\Omega_i)$ of $\Omega_i$, that is the minimal value of the following mixed problem in $\Omega_i$:

$$ egin{cases} -\Delta\varphi_i = \gamma_1(\Omega_i)\varphi_i & \text{in } \Omega_i \\ \varphi_i = 0 & \text{on } \partial\Omega_i \cap \Omega \\ \partial_v\varphi_i = 0 & \text{on } \partial\Omega_i \cap \partial\Omega \end{cases} $$

with $\varphi_i \neq 0$. It is known that $\gamma_1(\Omega_i) \ge 0$ and $\varphi_i > 0$ in $\Omega_i$. From the comparison principle, it follows that

if $\frac{\lambda k_i - \mu \omega_i}{d_i \mu} < \gamma_1(\Omega_i)$ then $w_i \equiv 0$.

In particular we see that if all the components $w_1, \dots, w_N$ are $\neq 0$ then necessarily

$$ \max_{i=1,\dots,k} \gamma_1(\Omega_i) < \max_{i=1,\dots,k} \frac{\lambda k_i - \mu \omega_i}{d_i \mu} = \bar{\gamma}. $$

Since $\Omega_1, \dots, \Omega_N$ is a $N$-partition of the set $\Omega$, we evince by Theorem 6.2 that necessarily

$$ \gamma_N(\Omega) < \bar{\gamma}. $$

As a result, we reach the desired conclusion recalling that the sequence of eigenvalues $\gamma_1 < \gamma_2 \le \gamma_3 \le \dots$ is unbounded. $\square$

Using Weyl's asymptotic law for the Laplacian with Neumann boundary conditions, we can obtain a more explicit bound on the constant $\bar{N}$. This is stated in the next theorem.

**Theorem 6.3.** The number $\bar{N}$ of Lemma 6.1 admits the following asymptotic upper bound

$$ \bar{N} \lesssim \frac{\omega_n}{(2\pi)^n} |\Omega| \left( \max_{i=1,\dots,k} \frac{\lambda k_i - \mu \omega_i}{d_i \mu} \right)^{n/2} \quad \text{for } \max_{i=1,\dots,k} \frac{\lambda k_i - \mu \omega_i}{d_i \mu} \to +\infty. $$

where $\omega_n$ is the volume of the unit sphere in $\mathbb{R}^n$.

We recall that $A \lesssim B$ for $B \to +\infty$ means that $A \le B + o(B)$ and $o(B)/B \to 0$ as $B \to +\infty$. This estimate agrees with the intuition that doubling the size of the domain $\Omega$ would also double the number of groups of predators it can sustain.
*Proof.* For a fixed smooth domain $\Omega \subset \mathbb{R}^n$, if we let $N(\gamma)$ stand for the number of eigenvalues for the Laplace operator with homogeneous Neumann boundary conditions in $\Omega$ which are less than $\gamma$, it can be shown that

$$N(\gamma) = \frac{\omega_n}{(2\pi)^n} |\Omega|^{\gamma^{n/2}} + o(\gamma^{n/2}).$$

Substituting this expression in the bound found in Lemma 6.1 proves the statement. □

We can now extend the result in Lemma 6.1 to the original competitive system. We show in particular

**Theorem 6.4.** There $\bar{\beta} > 0$ such that if $\beta > \bar{\beta}$ and $v_\beta = (w_{1,\beta}, \dots, w_{N,\beta}, u_\beta)$ is a solution to (1.1) then

* either at most $\bar{N}$ components of $w_\beta = (w_{1,\beta}, \dots, w_{N,\beta})$ are strictly positive and the others are zero;

* or the solution is such that

$$\|(w_{1,\beta}, \dots, w_{N,\beta})\|_{C^{0,\alpha}(\Omega)} + \|u_\beta - \lambda/\mu\|_{C^{2,\alpha}(\Omega)} = o_\beta(1)$$

for every $\alpha \in (0,1)$.

We recall that, by Proposition 4.1, the set of solutions of (3.1) is pre-compact in $H^1(\Omega) \cap C^{0,\alpha}(\bar{\Omega})$.

*Proof.* The statement will follow from some approximation results in combination with Lemma 6.1. First of all, we want to show that if $h > \bar{N}$ components of $w_\beta$ are non zero and $\beta$ sufficiently large, then the solutions is close to the solutions $(0, \dots, 0, 0)$ or $(0, \dots, 0, \lambda/\mu)$. We have

**Claim.** Let $v_\beta = (w_\beta, u_\beta) \in \mathcal{F}(\Omega)$ be a family of solutions to (4.1) and let us assume that there exists a solution $v = (w, u) = (w_1, \dots, w_N, u)$ with $(w_1, \dots, w_N) \in \text{Lip}(\bar{\Omega})$ and $u \in C^{2,\alpha}(\bar{\Omega})$ to (4.3) with $h$ components of $w$ non zero (with $1 \le h \le N$) such that $v_\beta \to v$ as $\beta \to +\infty$ in $H^1(\Omega) \cap C^{0,\alpha}(\bar{\Omega})$. Then there exists $\bar{\beta} > 0$ sufficiently large such that the solution $v_\beta$ has exactly $h$ components of $w_\beta$ that are non zero for $\beta \ge \bar{\beta}$.
Let us first show how to use the claim in order to derive the first part of the theorem. Consider a family of solutions $(w_β, u_β) \in \mathcal{F}(\Omega)$ of the system (4.1), with $\beta \to +\infty$. By Lemma 6.1 we already know that any solution of (4.3) has at most $\tilde{N}$ non trivial components among those of $w_β$. Let us assume that $w_β$ contains a sub-family (which we shall not relabel) with $\beta \to +\infty$ that has more than $\tilde{N}$ non zero components. Proposition 4.1 implies that, up to a subsequence, $v_β \to v$ in $H^1(\Omega) \cap C^{0,\alpha}(\bar{\Omega})$, where $v$ solves the limit system (4.3), and thus has at most $\tilde{N}$ components among those of $w$, in contradiction with our claim.

We now prove the claim, arguing by contradiction and adopting the scheme of Lemma 4.4. Let $v_n$ be any sequence satisfying the assumptions of the claim and let $v$ be its limit for $\beta_n \to +\infty$. By the maximum principle, if the corresponding $w$ has $h$ non trivial components, then necessarily $u$ is strictly positive. Up to a relabelling, we can assume that the first $h$ components $(w_1, \dots, w_h)$ are also non zero, while $w_{h+1} = \dots = w_N = 0$. As a result, the sub-vector $(w_1, \dots, w_h, u)$ satisfies the conclusions of Proposition 4.1, and, in particular, the set

$$ \mathcal{N} := \left\{ x \in \Omega : \sum_{i=1}^{h} w_i(x) = 0 \right\} $$

is a rectifiable set of co-dimension 1, made of the union of a finite number of $C^{1,\alpha}$ smooth sub-manifolds (points if $\Omega \subset \mathbb{R}$, curves if $\Omega \subset \mathbb{R}^2$, and in general embedded surfaces if $\Omega \subset \mathbb{R}^N$). For any $n \in \mathbb{N}$, we introduce the renormalized solution

$$ \bar{v}_n := \left( \frac{w_{1,\beta_n}}{\|w_{1,\beta_n}\|_{L^\infty(\Omega)}}, \dots, \frac{w_{N,\beta_n}}{\|w_{N,\beta_n}\|_{L^\infty(\Omega)}}, u_{\beta_n} \right) $$

which is well defined since, by assumption, for any $n$, and $1 \le i \le N$, $w_{i,n} > 0$. Let us observe that, since the first $h$ components of $v_n$ converge to non zero functions as $\beta \to +\infty$, the $L^\infty$ norm of these components is bound from above and away from zero uniformly in $n$. That is, for the first $h$ components of $\bar{v}_n$ we have

$$ \bar{w}_{i,\beta_n} = C_{i,n} w_{i,\beta_n} \quad \text{with} \quad 0 < C < C_{i,n} < C^{-1} \quad \forall 1 \le i \le k, n \in \mathbb{N}. $$

By the convergence of $w_{\beta_n}$ we deduce that the first $h$ components of $\bar{w}_{\beta_n}$ converge in $H^1(\Omega) \cap C^{0,\alpha}(\bar{\Omega})$ to a non zero scaling of the first $h$ components of $\bar{w}$. Consequently $(\bar{w}_1, \dots, \bar{w}_h)$ and $(w_1, \dots, w_h)$ share the zero set $\mathcal{N}$.
Including the scaling in the system, it follows that $\bar{v}_n$ is a solution of

$$ (6.1) \quad \left\{ \begin{array}{ll} -d_i \Delta \bar{w}_{i, \beta_n} = (-\omega_i + k_i u_{\beta_n} - \beta \sum_{j \neq i} a_{ij} \|w_{j, \beta_n}\|_{L^\infty(\Omega)} w_{j, \beta_n}) \bar{w}_{i, \beta_n} & \\ -D \Delta u_{\beta_n} = (\lambda - \mu u_{\beta_n} - \sum_{i=1}^N k_i w_{i, \beta_n}) u_{\beta_n} & \\ \partial_v \bar{w}_{i, \beta_n} = \partial_v u_{\beta_n} = 0 & \text{on } \partial \Omega \end{array} \right. $$

We are chiefly interested in the equations satisfied by the densities $\bar{w}_{i, \beta_n}$ for $h + 1 \le i \le N$. The maximum principle implies that $u_n \le \lambda/\mu$: testing the $i$-th equation with the density $\bar{w}_{i, \beta_n}$ itself and using its positivity, we find

$$ \int_{\Omega} |\nabla \bar{w}_{i, \beta_n}|^2 \le \frac{k_i}{d_i} \frac{\lambda}{\mu} |\Omega|. $$

Since by definition $\|\bar{w}_{1, \beta_n}\|_{L^\infty(\Omega)} = 1$, $\bar{w}_{i, \beta_n}$ is uniformly bounded in $H^1(\Omega)$ and it admits a weak limit $\bar{w}_i \in H^1(\Omega)$. The compact embedding in $L^2(\Omega)$ and the boundedness in $L^\infty(\Omega)$ also imply that $\bar{w}_{i, \beta_n} \to \bar{w}_i$ strongly in $L^p(\Omega)$ for any $p < \infty$. Let us show that $\bar{w}_i$ is not zero: for each $h + 1 \le i \le k$ (the other components are non trivial by assumption) we have that

$$ \left\{ \begin{array}{ll} -d_i \Delta \bar{w}_{i, \beta_n} \le (-\omega_i + k_i u_{\beta_n}) \bar{w}_{i, \beta_n} & \\ \partial_v \bar{w}_{i, \beta_n} = 0 & \text{on } \partial \Omega. \end{array} \right. $$

Let $g_{i,n} \in H^1(\Omega)$ be a solution to

$$ \left\{ \begin{array}{ll} -d_i \Delta g_{i,n} + \omega_i g_{i,n} = k_i u_{\beta_n} \bar{w}_{i, \beta_n} & \\ \partial_v g_{i,n} = 0 & \text{on } \partial \Omega. \end{array} \right. $$

By standard arguments we know that $0 \le \bar{w}_{i, \beta_n} \le g_{i,n}$ and that

$$ \|g_{i,n}\|_{C^{0,a}(\Omega)} \le C \|g_{i,n}\|_{W^{2,p}(\Omega)} \le C \frac{k_i}{d_i} \frac{\lambda}{\mu} \|\bar{w}_{i, \beta_n}\|_{L^p(\Omega)} $$

for any $N/2 < p < \infty$ and suitable $C$ and $\alpha > 0$. As a result, using the order relationship between $w_{i, \beta_n}$ and $g_{i,n}$, we have

$$ 1 = \|\bar{w}_{i, \beta_n}\|_{L^\infty(\Omega)} \le C \frac{k_i}{d_i} \frac{\lambda}{\mu} \|\bar{w}_{i, \beta_n}\|_{L^p(\Omega)}. $$

Thus, the strong limit $\bar{v}$ in $L^p(\Omega)$ of $\bar{v}_n$ has all of its components which are non trivial. We are now in a position to reach the aimed contradiction. Let us consider the equation satisfied by
$\bar{w}_{i,\beta}$ for $h+1 \le i \le k$: for simplicity, we scale back the first $h$ densities and find that

$$ egin{cases} -d_i \Delta \bar{w}_{i, \beta_n} \le (-\omega_i + k_i u_{\beta_n} - \beta_n \sum_{1 \le j \le h} w_{j, \beta_n}) \bar{w}_{i, \beta_n} \\ \partial_v \bar{w}_{i, \beta_n} = 0 \end{cases} \quad \text{on } \partial \Omega. $$

From the previous discussion, $\bar{w}_{i,\beta_n} \to \bar{w}_i$ in $L^2(\Omega)$ and the limit is non trivial. We let

$$ c_i := \int_{\Omega} \bar{w}_i^2 > 0. $$

For any $\varepsilon > 0$ and $n \in \mathbb{N}$, we consider the sets

$$ \Omega_{\varepsilon,n} := \left\{ x \in \Omega : \mathrm{dist}(x, \partial\Omega) \ge \varepsilon \text{ and } \inf_{m \ge n} \left( \sum_{i=1}^{h} w_{i, \beta_m}(x) \right) \ge \varepsilon \right\}. $$

By the uniform convergence of $v_{\beta_n}$ and the properties of its limit configuration we know that each $\Omega_{\varepsilon,n}$ is closed and $\Omega_{\varepsilon_1,n_1} \subseteq \Omega_{\varepsilon_2,n_2}$ whenever $\varepsilon_1 > \varepsilon_2$ and $n_1 < n_2$. Finally,

$$ \bigcup_{\varepsilon>0, n \in \mathbb{N}} \Omega_{\varepsilon,n} = \Omega \setminus \mathcal{N}. $$

As we have already recalled, $\mathcal{L}^N(\mathcal{N}) = 0$. Hence, it follows that for any $\delta > 0$, there exist $\bar{\varepsilon} > 0$ and $n^* \in \mathbb{N}$ such that $\mathcal{L}^N(\Omega \Delta \Omega_{\varepsilon,n}) \le \delta$ for $0 < \varepsilon < \bar{\varepsilon}$ and $n \ge n^*$ (here $A \Delta B$ is the symmetric difference of the sets $A$ and $B$). By the absolute continuity of the Lebesgue integral and the uniform integrability of converging sequences, there exists $\bar{\delta} > 0$ (and consequently $\bar{\varepsilon}$ and $n^*$) such that

$$ \int_{\Omega_{\varepsilon,n^*}} \bar{w}_{i, \beta_m}^2 \ge \frac{c_i}{2} > 0 \quad \text{for } m \in \mathbb{N} \text{ sufficiently large.} $$

On the other hand, testing the equation in $\bar{w}_{i,\beta_m}$ by $\bar{w}_{i,\beta_m}$ itself, we obtain

$$ \int_{\Omega} \left[ d_i |\nabla \bar{w}_{i, \beta_m}|^2 + \omega_i \bar{w}_{i, \beta_m}^2 + \beta_m \left( \sum_{1 \le j \le h} w_{j, \beta_m} \right) \bar{w}_{i, \beta_m}^2 \right] \le k_i \int_{\Omega} u_{\beta_m} \bar{w}_{i, \beta_m}^2 \le k_i \frac{\lambda}{\mu} |\Omega|. $$

Since the terms of the left hand side are positive, we can localize the integral on the sets $\Omega_{\varepsilon,n^*}$ and find

$$ \beta_m \inf_{\Omega_{\varepsilon,n^*}} \left( \sum_{1 \le j \le h} w_{j, \beta_m} \right) \int_{\Omega_{\varepsilon,n^*}} w_{i, \beta_m}^2 \le k_i \frac{\lambda}{\mu} |\Omega| $$

that is

$$ 0 < \frac{c_i}{2} < \int_{\Omega_{\varepsilon,n^*}} \bar{w}_{i, \beta_m}^2 \le \frac{1}{\beta_m \bar{\varepsilon}} \cdot k_i \frac{\lambda}{\mu} |\Omega| $$

a contradiction when $\beta_m$ is sufficiently large, and this proves the first claim.
At this point, we have established that positive solutions must converge to either one of the solutions $(0, \dots, 0, 0)$ or $(0, \dots, 0, \lambda/\mu)$. To conclude the proof, we show that they can only converge to the latter. For this, we can follow the same reasoning as for Lemma 4.5. That is, we suppose that $v_n \to (0, \dots, 0, 0)$ in $H^1(\Omega) \cap C^{0,\alpha}(\bar{\Omega})$ and find that for $n$ large enough

$$ \begin{cases} -d_i \Delta w_{i,\beta_n} = (-\omega_i + k_i u_{\beta_n} - \beta_n \sum_{j \neq i} w_{j,\beta_n}) w_{i,\beta_n} \leq -\frac{\omega_i}{2} w_{i,\beta_n} \\ \partial_v w_{i,\beta_n} = 0 \end{cases} \quad \text{on } \partial\Omega $$

which implies that each $w_{i,\beta_n}$ must be identically zero, in contradiction with our positivity assumption. This concludes the proof of Theorem 6.4. □

Combining the previous results, we can state the following

**Theorem 6.5.** Let $\delta > 0$ and, for arbitrary $N \ge 1$, let us consider an arbitrary family of coefficients such that

$$ \delta < d_1, \dots, d_N, \omega_1, \dots, \omega_N, k_1, \dots, k_N < \frac{1}{\delta}. $$

For $\beta \ge 0$, let $\mathcal{S}$ be the set of solutions $(w_1, \dots, w_N, u) \in \mathcal{F}(\Omega)$ to (4.1) with any number of components and coefficients as above. For any $(w_1, \dots, w_N, u) \in \mathcal{S}$ we associated

$$ P(w_1, \dots, w_N, u) = \int_{\Omega} \sum_{i=1}^{N} w_i. $$

Then $\bar{N} \in \mathbb{N}$ for which we have two alternatives

* either there exists $(\bar{w}_1, \dots, \bar{w}_{\bar{N}}, \bar{u}) \in \mathcal{S}$ such that

$$ P(\bar{w}_1, \dots, \bar{w}_{\bar{N}}, \bar{u}) = \max_{N \ge 1, (w_1, \dots, w_N, u) \in \mathcal{S}} P(w_1, \dots, w_N, u); $$

* there exist a sequence $(\bar{w}_{1,n}, \dots, \bar{w}_{\bar{N},n}, \bar{u}_n) \in \mathcal{S}$ and functions $(w_1, \dots, w_N) \in \text{Lip}(\bar{\Omega})$ and $u \in C^{2,\alpha}(\bar{\Omega})$ such that

- $(\bar{w}_{1,n}, \dots, \bar{w}_{\bar{N},n}, \bar{u}_n)$ are solutions of (4.1) for $\beta_n \to +\infty$;

- $(w_{1,\beta}, \dots, w_{N,\beta}) \to (w_1, \dots, w_N)$ in $C^{0,\alpha} \cap H^1(\bar{\Omega})$, $u_\beta \to u$ in $C^{2,\alpha}(\bar{\Omega})$;

- $(w_1, \dots, w_N, u)$ solves (4.3) and

$$ P(\bar{w}_1, \dots, \bar{w}_{\bar{N}}, \bar{u}) = \sup_{N \ge 1, (w_1, \dots, w_N, u) \in \mathcal{S}} P(w_1, \dots, w_N, u). $$
We stress the fact that we have imposed no conditions on $\beta > 0$ and $\mu > 0$. The proof of this theorem in contained in the previous results.

## 7. THE QUESTION OF OPTIMAL REPARTITION OF DENSITIES

We now show that, in some settings, the configuration that maximizes the total population of predators (that is, the solution mentioned in Theorem 6.5) contains more than one non zero component of $w_i$.

To this end, we first consider a simplified version of the system (3.1), viz.

$$ (7.1) \qquad \begin{cases} -\Delta w_1 = (-\omega + ku - \beta w_2) w_1 \\ -\Delta w_2 = (-\omega + ku - \beta w_1) w_2 \\ -\Delta u = (\lambda - kw_1 - kw_2) u \\ \partial_v w_i = \partial_v u = 0 \end{cases} \quad \text{on } \partial\Omega. $$

We observe that we are considering here a system with indistinguishable densities of predators (all the characterizing parameters in the equations are the same for all groups). Furthermore, we assume here that the parameter $\mu$ vanishes. An extension of the result in the case $\mu > 0$ will be presented later. We recall (see Proposition 2.4 and Definition 2.6) that (7.1) only admits as simple solutions (that is, solutions with only one of $w_1$ and $w_2$ non zero) the constant solution:

$$ (W, U) = \left(\frac{\lambda}{k}, \frac{\omega}{k}\right) \implies \int_{\Omega} W = \frac{\lambda}{k}|\Omega|, \int_{\Omega} U = \frac{\omega}{k}|\Omega| $$

Under these assumptions, we have

**Lemma 7.1.** Let $(w_1, w_2, u)$ be an arbitrary solution of (7.1) with $\beta > 0$. If all of its components are non zero and non constant, then

$$ \int_{\Omega} \sum_{i}^{2} w_{i} = \frac{\lambda}{k} |\Omega| + \frac{1}{k} \int |\nabla \log u|^2 > \frac{\lambda}{k} |\Omega| $$

and

$$ \int_{\Omega} u = \frac{\omega}{k} |\Omega| + \frac{\beta}{\lambda} \int w_1 w_2 + \frac{\omega}{k\lambda} \int |\nabla \log u|^2 > \frac{\omega}{k} |\Omega|. $$

**Remark 7.2.** Equivalently, we could have compared the solutions of (7.1) with $\beta > 0$ with any non zero solution in the case $\beta = 0$ (see Lemma 3.2).
*Proof.* We recall that, owing to Lemma 3.5, if one component is constant, so are the others. Thus, we can assume that $u$ is not constant. The existence of non constant solutions is already known thanks to Theorems 3.9 and 3.15. We consider the equation in $u$. By the maximum principle, $u > 0$ and thus we can divide the two sides of the equation by $u$ and integrate over $\Omega$. This yields

$$ (7.2) \qquad \int_{\Omega} \sum_{i}^{2} w_i = \frac{1}{k} \int \left( \lambda + \frac{\Delta u}{u} \right) = \frac{\lambda}{k} |\Omega| + \frac{1}{k} \int |\nabla \log u|^2 > \frac{\lambda}{k} |\Omega|. $$

Here, the strict inequality follows by the fact that $u$ is not a constant. Similarly, integrating directly the equations of the system, summing them and using the previous identity, we get

$$ \int_{\Omega} u = \frac{\omega}{k} |\Omega| + \frac{\beta}{\lambda} \int w_1 w_2 + \frac{\omega}{k\lambda} \int |\nabla \log u|^2 > \frac{\omega}{k} |\Omega| $$

and this concludes the proof.

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

As a result, according to system (7.1), for any $\beta > 0$, competition is always advantageous both for the predators and for the preys. Indeed, the total population of predators (and preys) is greater in the case of two groups competing for the same territory, than in the case of only one group.

We now extend this result to the model (3.1), that is when $\mu > 0$. To this end, let us first observe that the same computations as above yield the identity

$$ \int_{\Omega} \sum_{i}^{2} w_i = \frac{\lambda}{k} |\Omega| - \frac{\mu}{k} \int u + \frac{1}{k} \int |\nabla \log u|^2 \\ = \frac{\lambda k - \mu \omega}{k^2} |\Omega| + \frac{\mu}{k} \left( \frac{\omega}{k} |\Omega| - \int u \right) + \frac{1}{k} \int |\nabla \log u|^2 $$

which, by uniform convergence of the densities as $\beta \to +\infty$, is valid also in the limit case of segregation $\beta = +\infty$. Unlike the case $\mu = 0$, a direct comparison of the previous formula with the case of only one population of predators is not immediate. Indeed, in general, we can show that the second term in the last expression is negative, that is $\int u > \frac{\omega}{k} |\Omega|$.

As a result, we need to carefully estimate the various contributions in this identity, in order to show that, when $\mu$ is sufficiently small, the total population of predators is large in the case of two non trivial components $w_1$ and $w_2$.
We wish to emphasize that this is a rather delicate task. Indeed, an a priori estimate on the solutions which is independent of $\mu$ may not be true, for two main reasons.

* From the equation in $u$, we can only say a priori that $u \le \lambda/\mu$. If $\mu \to 0$, we have no reason to conclude that the solutions of (3.1) converge to solutions of (7.1).

* One may wonder whether the previous bound is not sharp and that it may be achieved only by "spurious" solutions such as $(0,0,\lambda/\mu)$. But this assertion is not true in general. To see this we can recall that, by Theorem 3.9, non constant (as it were "genuine") solutions bifurcate from

$$ \left( \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2}, \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2}, \frac{\lambda \beta + 2k \omega}{\mu \beta + 2k^2} \right) \quad \text{for } \beta \frac{\lambda k - \mu \omega}{\mu \beta + 2k^2} = \gamma_n $$

where $\gamma_n$ is the n-th eigenvalue of the Laplace operator with Neumann boundary conditions. For $\mu$ sufficiently small and $\gamma_n$ large (and, consequently, $\beta$ large), we have non constant solutions for which $u$ is close (at least) in the uniform topology to the upper bound $\lambda/\mu$.

We thus focus on the one-dimensional case, for which (see Theorem 5.1) we have already established the existence of segregated solutions and pointed out their symmetries (Remark 3.11). As a result, for $\Omega = (-a,a)$, $a > 0$ and $\mu > 0$ sufficiently small, we have a continuum of solutions such that $w_1 - w_2$ vanishes only for $x = 0$. Sending the competition parameter to infinity $\beta \to +\infty$, by Lemma 5.2 we can start by considering the system

$$ (7.3) \qquad \begin{cases} -w'' = (-\omega + ku)w \\ -u'' = (\lambda - \mu u - kw)u \\ w(0) = w'(a) = u'(0) = u'(a) = 0 \end{cases} \quad \text{in } (0,a) $$

for $\mu > 0$, for which the identity (7.2) reduces to

$$ \int_0^a w = \frac{\lambda}{k}a - \frac{\mu}{k}\int_0^a u + \frac{1}{k}\int_0^a |(\log u)'|^2 $$

**Proposition 7.3.** Let $(w, u)$ be any classical solution of (7.3) with both components non negative and nontrivial. For $\mu > 0$ sufficiently small

$$ \int_0^a w > \frac{\lambda}{k}a. $$
Let us observe that, since the solutions of (3.1) converge for $\beta \to +\infty$ to segregated solutions, the previous result implies an improvement of Theorem 6.5, and in particular we have

**Theorem 7.4.** Under the assumptions of Theorem 6.5, let us assume moreover that the coefficients in (4.1) do not depend on the index $i$ and that $\Omega = (a, b) \subset \mathbb{R}$. If $\mu > 0$ sufficiently small the solution of (4.1) that maximizes

$$P(w_1, w_2, \dots, w_N, u) = \int_{\Omega} \sum_{i=1}^{N} w_i$$

has at least $N \ge 2$ non trivial components and $\beta > 0$.

We divide the proof of Theorem 7.4 in two separate results. Letting all the parameters in (7.3) fixed a part from $\mu > 0$, we have

**Lemma 7.5.** Let $(w, u)$ be an arbitrary classical solution of (7.3) with both components non negative and non zero. For any $\varepsilon > 0$ there exists $\bar{\mu} > 0$ such that

$$\mu \int_{0}^{a} u \le \varepsilon \quad \text{if } \mu \in (0, \bar{\mu}).$$

*Proof.* Let $(w_n, u_n)$ be a sequence of positive solutions to (7.3) with $\mu = \mu_n \to 0$, and let us assume, by contradiction, that

$$\mu_n \int_0^a u_n > C > 0.$$

The maximum principle, as already observed, implies that $u_n \le \lambda/\mu_n$, hence

$$\mu_n \int_0^a u_n \le a\lambda.$$

Thus, we can assume that, up to a subsequence,

$$\lim_{n \to +\infty} \mu_n \int_0^a u_n = \gamma \in (0, a\lambda].$$

for some constant $\gamma > 0$. We can introduce the scaled functions $(\bar{w}_n, \bar{u}_n)$ as

$$\bar{w}_n := \left( \int_0^a (w'_n)^2 \right)^{-1/2} w_n, \quad \bar{u}_n := \left( \frac{1}{a} \int_0^a u_n \right)^{-1} u_n,$$
which are solutions to

$$
(7.4) \qquad \begin{cases}
-\bar{w}_n'' = (-\omega + k_n \bar{u}_n) \bar{w}_n \\
-\bar{u}_n'' = (\lambda - \bar{u}_n \bar{u}_n - k'_n \bar{w}_n) \bar{u}_n & \text{in } (0, a) \\
\bar{w}_n(0) = \bar{w}_n'(a) = \bar{u}_n'(0) = \bar{u}_n'(a) = 0
\end{cases}
$$

where we have defined

$$
k_n := k \frac{1}{a} \int_0^a u_n, \quad \bar{\mu}_n := \mu_n \frac{1}{a} \int_0^a u_n, \quad k'_n := k \left( \int_0^a (w'_n)^2 \right)^{1/2}.
$$

We observe that, from the starting assumption,

$k_n \to +\infty \quad \text{and} \quad \bar{\mu}_n \to \bar{\mu}_\infty \in (0, \lambda),$

while we have no information on $k'_n$. Moreover, by definition and the Dirichlet boundary con-
dition at zero, the sequence $\{\bar{w}_n\}_{n\in\mathbb{N}}$ is bounded in $H^1(0,a)$, and by positivity, also $\{\bar{u}_n\}_{n\in\mathbb{N}}$
is bounded in $H^1(0,a)$; Indeed, by testing the equation in $\bar{u}_n$ with $\bar{u}_n$ itself, we obtain

$$
\int_0^a (\bar{u}'_n)^2 + \mu_n \bar{u}_n^3 \le \lambda \int_0^a \bar{u}_n^2
$$

and the claim follows from the assumption $\bar{\mu}_n \to C > 0$. By the embedding theorems we have
that, up to a subsequence, both $\{\bar{w}_n\}_{n \in \mathbb{N}}$ and $\{\bar{u}_n\}_{n \in \mathbb{N}}$ converge uniformly in $(0, a)$ to their
respective weak limits in the space $H^1(0, a)$, $\bar{w}_\infty$ and $\bar{u}_\infty$. Moreover, by renormalization and
strong convergence, we have

$$
\int_0^a \bar{u}_\infty = a
$$

and thus, in particular, $\bar{u}_\infty$ is non zero. Finally, from the equation in $\bar{w}_n$ we see that

$$
k_n \int_0^a \bar{u}_n \bar{w}_n^2 = \int_0^a (\bar{w}_n')^2 + \omega \bar{w}_n^2 \le C'.
$$

Since $k_n \to +\infty$, by the uniform convergence we have that

$$
\bar{u}_n \bar{w}_n \to \bar{u}_\infty \bar{w}_\infty \equiv 0 \quad \text{uniformly in } (0, a).
$$

**Step 1)** We now proceed and exclude the possibility that the sequence $k'_n$ is bounded. By the uniform convergence we have

$$
k'_n \bar{u}_n \bar{w}_n \to 0 \quad \text{uniformly in } (0, a).
$$
Passing to the limit in the equation satisfied by $\bar{u}_n$, it follows that $\bar{u}_\infty$ is a non trivial solution of

$$ (7.5) \qquad \begin{cases} -\bar{u}_\infty'' = (\lambda - \bar{\mu}_\infty \bar{u}_\infty) \bar{u}_\infty & \text{in } (0, a) \\ \bar{u}_\infty'(0) = \bar{u}_\infty'(a) = 0. \end{cases} $$

From Lemma 3.6, $\bar{u}_\infty$ can thus be only the constant $\lambda/\bar{\mu}_\infty$, and finally, thanks to renormalization, $\bar{\mu}_\infty = \lambda$ and $\bar{u}_\infty \equiv 1$. By the uniform convergence of $\bar{u}_n$ to its limit, for $n$ sufficiently large we find

$$ \bar{u}_n > \frac{1}{2} \quad \text{in } (0, a). $$

Using this estimate in the equation for $\bar{w}_n$, we get that

$$ \begin{cases} -\bar{w}_n'' = (-\omega + k_n \bar{u}_n) \bar{w}_n > (-\omega + k_n/2) \bar{w}_n \\ \bar{w}_n(0) = \bar{w}_n'(a) = 0. \end{cases} $$

Since $k_n \to +\infty$, for $n$ large enough we find

$$ (-\omega + k_n/2) > \left(\frac{\pi}{2a}\right)^2. $$

The right hand side is the principal eigenvalue of the operator. By the comparison principle we see that $\bar{w}_n \equiv 0$, against the assumption that the functions $(w_n, u_n)$ are positive in $(0, a)$. This possibility is thus ruled out.

**Step 2)** Hence it must be the case that $k'_n \to +\infty$: similarly we can show that $k'_n \ge Ck_n$. Indeed, let us assume by contradiction that $k'_n/k_n \to 0$. Testing the equation in $\bar{w}_n$ with $\varphi$ smooth and compactly supported, we find

$$ k_n \int_0^a \bar{u}_n \bar{w}_n \varphi = \int_a^b (\bar{w}_n' \varphi + \omega \bar{w}_n \varphi) $$

and the right-hand side is bounded in $n$, so that

$$ k'_n \int_0^a \bar{u}_n \bar{w}_n \varphi = \frac{k'_n}{k_n} k_n \int_0^a \bar{u}_n \bar{w}_n \varphi \to 0 $$

for all test function $\varphi$. It follows that the weak and uniform limit $\bar{u}_\infty$ of $\bar{u}_n$ solves again (7.5), and thus

$$ \bar{u}_n \to 1 \quad \text{uniformly in } (0, a). $$
We can then reach a contradiction as in Step 1).

**Step 3)** We now show that $k_n \ge Ck'_n$. Again by contradiction, let us assume that $k_n/k'_n \to 0$. Integrating the equation in $\bar{u}_n$ we obtain

$$k'_n \int_0^a \bar{u}_n \bar{w}_n = \int_0^a \bar{u}_n (\lambda - \mu_n \bar{u}_n)$$

and the right-hand side is bounded uniformly in $n$. It follows that

$$k_n \int_0^a \bar{u}_n \bar{w}_n = \frac{k_n}{k'_n} k'_n \int_0^a \bar{u}_n \bar{w}_n \to 0.$$

But then, testing the equation in $\bar{w}_n$ with $\bar{w}_n$ itself, we obtain, thanks to the renormalization of $\bar{w}_n$ and the uniform convergence,

$$0 < C' < \int_0^a (\bar{w}'_n)^2 + \omega \bar{w}_n^2 = k_n \int_0^a \bar{u}_n \bar{w}_n^2 \le \|\bar{w}_n\|_{L^\infty} \cdot k_n \int_0^a \bar{u}_n \bar{w}_n \to 0$$

a contradiction.

**Step 4)** In summary, so far we have shown that

$$k'_n \sim k_n \to +\infty.$$

We already know that, up to a subsequence, the sequence $\{\bar{w}_n\}_{n \in \mathbb{N}}$ converges uniformly to continuous function $\bar{w}_\infty$. We now show that $\bar{w}_\infty \equiv 0$. We argue by contradiction and assume that this is not the case. Then there exist $0 \le x_0 < x_1 \le a$ such that

$$\inf_{x \in [x_0, x_1]} \bar{w}_n > C > 0 \quad \text{for all } n \text{ sufficiently large}$$

Then for some positive constants $A, B$ and $n$ sufficiently large, we have

$$\begin{cases} \bar{u}_n < A & \text{in } (0, a) \\ -\bar{u}_n'' = (\lambda - \mu_n \bar{u}_n - k'_n \bar{w}_n) \bar{u}_n < -Bk'_n \bar{u}_n & \text{in } (x_0, x_1). \end{cases}$$

By comparison with the super-solution

$$x \mapsto A \cosh \left[ (Bk'_n)^{1/2} \left( x - \frac{x_0 + x_1}{2} \right) \right] / \cosh \left[ (Bk'_n)^{1/2} \left( \frac{x_1 - x_0}{2} \right) \right]$$

we have that, for $\varepsilon > 0$ small

$$\sup_{x \in [x_0 + \varepsilon, x_1 - \varepsilon]} \bar{u}_n \le A \frac{\cosh \left[ (Bk'_n)^{1/2} \left( \frac{x_1 - x_0}{2} - \varepsilon \right) \right]}{\cosh \left[ (Bk'_n)^{1/2} \left( \frac{x_1 - x_0}{2} \right) \right]}.$$
Turning back to the equation in $\bar{w}_n$, and recalling that $k_n = O(k'_n)$, we can pass to the limit and obtain

$$-\bar{w}_{\infty}'' = -\omega \bar{w}_{\infty} \quad \text{in } [x_0 + \epsilon, x_1 - \epsilon].$$

We now observe that the previous reasoning holds is true for any $\epsilon$ and any interval of positivity $[x_0, x_1]$ of $\bar{w}_{\infty}$. As a result, in any interval of positivity $[x_0, x_1] \subset [0, a]$ we have

$$\text{either} \quad \begin{cases} -\bar{w}_{\infty}'' = -\omega \bar{w}_{\infty} \\ \bar{w}_{\infty}(x_0) = \bar{w}_{\infty}(x_1) = 0 \end{cases} \quad \text{or} \quad \begin{cases} -\bar{w}_{\infty}'' = -\omega \bar{w}_{\infty} \\ \bar{w}_{\infty}(x_0) = \bar{w}_{\infty}'(x_1) = 0 \end{cases}$$

In both case, arguing as in Step 1), we see that $\bar{w}_{\infty} \equiv 0$ in $[x_0, x_1]$, meaning that there are no intervals in $(0, a)$ where $\bar{w}_{\infty}$ is positive, and thus

$$\bar{w}_n \to 0 \quad \text{uniformly in } (0, a).$$

Now we can repeat the reasoning of Step 3): integrating the equation in $\bar{u}_n$ we obtain

$$k_n \int_0^a \bar{u}_n \bar{w}_n \sim k'_n \int_0^a \bar{u}_n \bar{w}_n = \int_0^a \bar{u}_n (\lambda - \mu_n \bar{u}_n).$$

From the equation in $\bar{w}_n$ and the uniform limit proved before, we have

$$0 < C' < \int_0^a (\bar{w}_n')^2 + \omega \bar{w}_n^2 = k_n \int_0^a \bar{u}_n \bar{w}_n^2 \le \|\bar{w}_n\|_{L^\infty} \cdot k_n \int_0^a \bar{u}_n \bar{w}_n \to 0,$$

which is a contradiction. □

**Lemma 7.6.** Let $(w, u)$ be an arbitrary classical solution of (7.3) with both components non negative and non zero. There exist two constants $C > 0$ and $\bar{\mu} > 0$ such that

$$\int_0^a |(\log u)'|^2 > C \quad \text{if } \mu \in (0, \bar{\mu}).$$

*Proof.* Let us consider a sequence $(w_n, u_n)$ of positive solutions to (7.3) for $\mu = \mu_n$ such that

$$\lim_{n \to +\infty} \int_0^a |(\log u_n)'|^2 = 0$$

By the embedding theorems we know that

$$\lim_{n \to +\infty} \frac{\sup_{x,y \in (0,a)} |u_n(x) - u_n(y)|}{\|u_n\|_{L^\infty(0,a)}} = 0.$$
Moreover, since $u_n$ is positive, we have

$$ (7.6) \qquad \lim_{n \to +\infty} \frac{\inf_{(0,a)} u_n}{\sup_{(0,a)} u_n} = 1. $$

Let us now show that (7.6) implies that $u_n$ converges to a positive and finite constant. For this it suffices to show that $\{u_n\}_{n \in \mathbb{N}}$ is uniformly bounded from above and away from zero. Indeed, if $\inf u_n \to +\infty$, then for $n$ sufficiently large

$$ \left\{ \begin{array}{l} -w_n'' = (-\omega + ku_n) w_n > \left(\frac{\pi}{2a}\right)^2 w_n \\ w_n(0) = w_n'(a) = 0 \end{array} \right. \quad \implies \quad w_n \equiv 0. $$

On the other hand, if $\sup u_n \to 0$, then for $n$ sufficiently large

$$ \left\{ \begin{array}{l} -w_n'' = (-\omega + ku_n) w_n < \left(\frac{\pi}{2a}\right)^2 w_n \\ w_n(0) = w_n'(a) = 0 \end{array} \right. \quad \implies \quad w_n \equiv 0 $$

and in both cases we reach a contradiction with the positivity of $w_n$. This, together with (7.6) implies that $u_n$ converges uniformly to a positive constant $C$. Again by the equation in $w_n$, we see that necessarily

$$ u_n \to \frac{\omega}{k} + \frac{1}{k} \left(\frac{\pi}{2a}\right)^2 \quad \text{uniformly in } (0, a). $$

Up to a renormalization, we also infer that

$$ (\|w_n\|_{L^\infty(0,a)})^{-1} w_n \to \sin\left(\frac{\pi}{2a}x\right) $$

strongly in $H^1(0, a)$ and also uniformly. Integrating the equation in $u_n$, we have

$$ k \int_0^a u_n w_n = \lambda \int_0^a u_n - \mu_n \int_0^a u_n^2. $$

Thus, we have that $\{w_n\}_{n \in \mathbb{N}}$ is uniformly bounded in $(0, a)$, so that

$$ w_n \to C \sin \left( \frac{\pi}{2a} x \right) $$

for a non negative constant $C$. Consequently, using these information in the equation satisfied by $u_n$, we have that $u_n$ is bounded uniformly in $H^1(0, a)$, and we can thus take the weak limit of the equation to derive

$$ 0 = \left[ \lambda - kC \sin \left( \frac{\pi}{2a} x \right) \right] \left[ \frac{\omega}{k} + \frac{1}{k} \left( \frac{\pi}{2a} \right)^2 \right], $$
which is impossible to enforce and gives us the desired contradiction. The proof of Theorem
7.4 is thereby complete.
□

**Acknowledgements:** This work has been supported by the ERC Advanced Grant 2013 n. 321186 "ReaDi – Reaction-Diffusion Equations, Propagation and Modelling" held by Henri Berestycki while Alessandro Zilio was post-doc fellow of the research grant. This work was also partially supported by the French National Research Agency (ANR), within project NON-LOCAL ANR-14-CE25-0013.

REFERENCES

[1] Henri Berestycki and Alessandro Zilio. Predators-prey models with competition Part II: uniform regularity estimates. *In preparation*.

[2] Henri Berestycki and Alessandro Zilio. Predators-prey models with competition Part III: classification of solutions and rigidity of the model. *In preparation*.

[3] Henri Berestycki and Alessandro Zilio. Predators-prey models with competition Part IV: the properties of the parabolic system. *In preparation*.

[4] Henri Berestycki and Alessandro Zilio. Predators-prey models with competition: the emergence of packs and territoriality. *Preprint*.

[5] Boris Buffoni and John Toland. *Analytic theory of global bifurcation*. Princeton Series in Applied Mathematics. Princeton University Press, Princeton, NJ, 2003. An introduction.

[6] L. A. Caffarelli, A. L. Karakhanyan, and Fang-Hua Lin. The geometry of solutions to a segregation problem for nondivergence systems. *J. Fixed Point Theory Appl.*, 5(2):319–351, 2009.

[7] Monica Conti, Susanna Terracini, and G. Verzini. Asymptotic estimates for the spatial segregation of competitive systems. *Adv. Math.*, 195(2):524–560, 2005.

[8] Monica Conti, Susanna Terracini, and Gianmaria Verzini. A variational problem for the spatial segregation of reaction-diffusion systems. *Indiana Univ. Math. J.*, 54(3):779–815, 2005.

[9] Edward Conway, David Hoff, and Joel Smoller. Large time behavior of solutions of systems of nonlinear reaction-diffusion equations. *SIAM J. Appl. Math.*, 35(1):1–16, 1978.

[10] Edward D. Conway and Joel A. Smoller. Diffusion and the predator-prey interaction. *SIAM J. Appl. Math.*, 33(4):673–686, 1977.

[11] E. N. Dancer. Bifurcation theory in real Banach space. *Proc. London Math. Soc.* (3), 23:699–734, 1971.

[12] E. N. Dancer. Bifurcation theory for analytic operators. *Proc. London Math. Soc.* (3), 26:359–384, 1973.
[13] E. N. Dancer and Yi Hong Du. Competing species equations with diffusion, large interactions, and jumping nonlinearities. *J. Differential Equations*, 114(2):434–475, 1994.

[14] E. N. Dancer, Kelei Wang, and Zhitao Zhang. Dynamics of strongly competing systems with many species. *Trans. Amer. Math. Soc.*, 364(2):961–1005, 2012.

[15] Jack Dockery, Vivian Hutson, Konstantin Mischaikow, and Mark Pernarowski. The evolution of slow dispersal rates: a reaction diffusion model. *J. Math. Biol.*, 37(1):61–83, 1998.

[16] B. Helffer, T. Hoffmann-Ostenhof, and S. Terracini. Nodal domains and spectral minimal partitions. *Ann. Inst. H. Poincaré Anal. Non Linéaire*, 26(1):101–138, 2009.

[17] Masayasu Mimura. Asymptotic behaviors of a parabolic system related to a planktonic prey and predator model. *SIAM J. Appl. Math.*, 37(3):499–512, 1979.

[18] Paul H. Rabinowitz. Some global results for nonlinear eigenvalue problems. *J. Functional Analysis*, 7:487–513, 1971.

[19] P.H. Rabinowitz and H. Berestycki. *Théorie du degré topologique et applications à des problèmes aux limites non linéaires*. Université Paris VI, Laboratoire analyse numérique, 1975.

[20] D. Ruelle. Bifurcations in the presence of a symmetry group. *Archive for Rational Mechanics and Analysis*, 51:136–152, January 1973.

[21] D. Sattinger. *Branching in the Presence of Symmetry*. Society for Industrial and Applied Mathematics, 1983.

[22] Nicola Soave and Alessandro Zilio. Uniform bounds for strongly competing systems: the optimal Lipschitz case. *Arch. Ration. Mech. Anal.*, 218(2):647–697, 2015.

[23] Hugo Tavares and Susanna Terracini. Regularity of the nodal set of segregated critical configurations under a weak reflection law. *Calc. Var. Partial Differential Equations*, 45(3-4):273–317, 2012.

[24] Susanna Terracini, Gianmaria Verzini, and Alessandro Zilio. Spiraling asymptotic profiles of competition-diffusion systems. *Preprint*.

[25] Vito Volterra. Variations and fluctuations of the number of individuals in animal species living together. *Journal du Cons. Int. Explor. Mer*, 3:5–51, 1928.

E-mail address: hb@ehess.fr

ÉCOLE DES HAUTES ÉTUDES EN SCIENCES SOCIALES, PSL RESEARCH UNIVERSITY PARIS, CENTRE D'ANALYSE ET DE MATHÉMATIQUE SOCIALES (CAMS), CNRS, 54 BOUVELARD RASPAIL, 75006, PARIS

E-mail address: azilio@math.univ-paris-diderot.fr

UNIVERSITÉ PARIS DIDEROT - PARIS 7, LABORATOIRE J.-L. LIONS (CNRS UMR 7598), PARIS, FRANCE, 8 PLACE AURÉLIE NEMOURS, 75205, PARIS CEDEX 13