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

HAL
open science

# Adaptation in a heterogeneous environment II: To be three or not to be

Matthieu Alfaro, François Hamel, Florian Patout, Lionel Roques

► To cite this version:

Matthieu Alfaro, François Hamel, Florian Patout, Lionel Roques. Adaptation in a heterogeneous environment II: To be three or not to be. Journal of Mathematical Biology, 2023, 87 (5), 10.1007/s00285-023-01996-4. doi:10.1007/s00285-023-01996-4

HAL Id: hal-03793309

https://hal.inrae.fr/hal-03793309

Submitted on 30 Sep 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Adaptation in a heterogeneous environment
## II: To be three or not to be *

M. Alfaro<sup>a,b</sup>, F. Hamel<sup>c</sup>, F. Patout<sup>b</sup> and L. Roques<sup>b</sup>

### Abstract

We propose a model to describe the adaptation of a phenotypically structured population in a *H*-patch environment connected by migration, with each patch associated with a different phenotypic optimum, and we perform a rigorous mathematical analysis of this model. We show that the large-time behaviour of the solution (persistence or extinction) depends on the sign of a principal eigenvalue, $λ_H$, and we study the dependency of $λ_H$ with respect to *H*. This analysis sheds new light on the effect of increasing the number of patches on the persistence of a population, which has implications in agroecology and for understanding zoonoses; in such cases we consider a pathogenic population and the patches correspond to different host species. The occurrence of a *springboard* effect, where the addition of a patch contributes to persistence, or on the contrary the emergence of a detrimental effect by increasing the number of patches on the persistence, depends in a rather complex way on the respective positions in the phenotypic space of the optimal phenotypes associated with each patch. From a mathematical point of view, an important part of the difficulty in dealing with $H \ge 3$, compared to $H = 1$ or $H = 2$, comes from the lack of symmetry. Our results, which are based on a fixed point theorem, comparison principles, integral estimates, variational arguments, rearrangement techniques, and numerical simulations, provide a better understanding of these dependencies. In particular, we propose a precise characterisation of the situations where the addition of a third patch increases or decreases the chances of persistence, compared to a situation with only two patches.

**Keywords** Adaptation · Migration · Heterogeneous environment · Spillover · Eigenvalues

## 1 Introduction

In [12], we analysed PDE systems describing the dynamics of adaptation of a phenotypically structured population, under the effects of mutation, selection and migration in a two-patch environment, each patch being associated with a different phenotypic optimum. Consistently with current literature [8, 17, 20] in evolutionary biology, our analysis showed that migration between the two patches leads to a locally reduced fitness. This reduction

<footnote><sup>a</sup> Univ Rouen Normandie, LMRS, CNRS, Rouen, France</footnote>

<footnote><sup>b</sup> INRAE, BioSP, 84914, Avignon, France</footnote>

<footnote><sup>c</sup> Aix Marseille Univ, CNRS, I2M, Marseille, France</footnote>

<footnote>*This work has received funding from the French ANR RESISTE (ANR-18-CE45-0019) and DEEV (ANR-20-CE40-0011-01) projects, from Excellence Initiative of Aix-Marseille Université - A*MIDEX, a French "Investissements d'Avenir" program, and from the région Normandie BIOMA-NORMAN (21E04343) project. The authors acknowledge support of the Institut Henri Poincaré (UAR 839 CNRS-Sorbonne Université) and LabEx CARMIN (ANR-10-LABX-59-01).</footnote>
in fitness is known as a migration load [11] and implies decreased chances of persistence of the global population in a two-patch system compared to a single patch environment.

From an epidemiological viewpoint, the two patches can be interpreted as two different types of hosts (different species, or different genetic variants). Thus, the above result means that it is more difficult for a pathogen to adapt and establish in a two-patch environment connected by migration than in a single patch environment. This observation is in agreement with one of the fundamental principles of agroecology [4, 10], which is that host species diversification should lead to a higher resilience of agroecosystems.

Nevertheless, as discussed in [15], what is right for two patches may not be necessarily right for three patches or more. The presence of a third host may indeed cause a springboard effect, leading to higher chances of persistence of the pathogen, compared to an environment with two hosts. This is a common pattern in zoonoses. For instance the main reservoirs of influenza A virus are the aquatic birds, but it is widely accepted that these viruses need to adapt in an intermediate host (such as pig or poultry) before they lead to an outbreak in human populations [21, 25]. Coronaviruses, including SARS-CoV-1 and SARS-CoV-2, are also zoonotic, with bats as presumed main reservoir. Intermediate hosts are also suspected to have played an important role in the 2003 and 2019 outbreaks [13, 14].

Up to our knowledge, there is no rigorous mathematical framework to study the effect of host diversification on the adaptation of a pathogen, when there are three hosts or more. Here, using the same assumptions as in [12], we study a system for a phenotypically structured population, under the effects of mutation, selection and migration between several hosts. Our main goal is to check whether the introduction of a third host leads to the above-mentioned springboard effect, or whether, on the contrary, it reduces the chances of persistence of the global population compared to a situation with two hosts. Thus, we mostly focus on the case of three hosts, described by the following system:

$$ 
\left\{ 
\begin{aligned}
\partial_t u_1(t, \mathbf{x}) &= \frac{\mu^2}{2} \Delta u_1(t, \mathbf{x}) + f_1(\mathbf{x}, u_1(t, \cdot)) + \delta \left[ \frac{u_2(t, \mathbf{x}) + u_3(t, \mathbf{x})}{2} - u_1(t, \mathbf{x}) \right], \\
\partial_t u_2(t, \mathbf{x}) &= \frac{\mu^2}{2} \Delta u_2(t, \mathbf{x}) + f_2(\mathbf{x}, u_2(t, \cdot)) + \delta \left[ \frac{u_1(t, \mathbf{x}) + u_3(t, \mathbf{x})}{2} - u_2(t, \mathbf{x}) \right], \\
\partial_t u_3(t, \mathbf{x}) &= \frac{\mu^2}{2} \Delta u_3(t, \mathbf{x}) + f_3(\mathbf{x}, u_3(t, \cdot)) + \delta \left[ \frac{u_1(t, \mathbf{x}) + u_2(t, \mathbf{x})}{2} - u_3(t, \mathbf{x}) \right],
\end{aligned}
\right. 
$$

for $t > 0$, and $\mathbf{x} = (x_1, \dots, x_n) \in \mathbb{R}^n$. Here, and as in [12], $\mathbf{x}$ is a breeding value for phenotype (for short, we simply write “phenotype” in the sequel), and corresponds to a set of $n \ge 1$ traits. The unknowns $u_i$ are the phenotype densities in the hosts $i \in \{1, 2, 3\}$, the Laplace operator describes the mutation effects on the phenotype, $\mu > 0$ is a mutational parameter, $\delta > 0$ the migration rate, and the functions $f_i(\mathbf{x}, u_i(t, \cdot))$ describe the growth of the phenotype $\mathbf{x}$ in the host $i$ (the precise assumptions on the functions $f_i$ are given below). Note that the migration and mutation parameters are assumed to be identical over the three hosts. In particular, each host sends migrants to the other hosts, at a rate $\delta$, and the amount of migrants is evenly split between the other hosts, hence the factor $1/2$.

Some of our results also deal with the general case of $H \ge 2$ hosts. Thus, we also consider the following system:

$$ (1) \quad \partial_t u_i(t, \mathbf{x}) = \frac{\mu^2}{2} \Delta u_i(t, \mathbf{x}) + f_i(\mathbf{x}, u_i(t, \cdot)) + \delta \left( \sum_{\substack{k=1 \\ k \neq i}}^{H} \frac{u_k(t, \mathbf{x})}{H-1} - u_i(t, \mathbf{x}) \right), \quad 1 \le i \le H, $$
for $t > 0$ and $\mathbf{x} \in \mathbb{R}^n$. Again, we assume that each host sends migrants to the other hosts, at a rate $\delta$, and that the amount of migrants is evenly split between the other hosts, hence the factor $1/(H - 1)$. With these assumptions, the migration rate from a given host $i$ towards the pool made up of all the $H - 1$ other hosts does not depend on $H$. In the particular case $H = 1$, there is no migration and (1) reduces to the scalar equation $\partial_t u_1(t, \mathbf{x}) = (\mu^2/2)\Delta u_1(t, \mathbf{x}) + f_1(\mathbf{x}, u_1(t, \cdot))$.

Each host is characterised by a phenotype optimum $\mathbf{O}_i \in \mathbb{R}^n$. Namely, the fitness (reproductive success) of a phenotype $\mathbf{x}$ in the host $i$ is described by a function $r_i(\mathbf{x})$ that decreases away from the optimum $\mathbf{O}_i$. As in [12] we use Fisher's geometrical phenotype-to-fitness model (FGM) [16, 23], which assumes that:

$$ (2) \qquad r_i(\mathbf{x}) = r_{\max} - \alpha \frac{\|\mathbf{x} - \mathbf{O}_i\|^2}{2}, $$

where $r_{\max} \in \mathbb{R}$ is the fitness of the optimal phenotype ($\mathbf{O}_i$) in the host $i$, $\alpha > 0$ is a measure of the intensity of selection, and $\| \cdot \|$ denotes the Euclidean norm in $\mathbb{R}^n$ (we use the same notation for any dimension $n \ge 1$, we also denote $\cdot$ the Euclidean inner product in $\mathbb{R}^n$).

In [12], we considered two types of growth functions, either linear (corresponding to a Malthusian growth) or logistic-like. Here, we only consider the logistic-like growth term, which describes a nonlocal competition between the phenotypes within each host. This is the most biologically relevant and mathematically involved case. It corresponds to $f_i(\mathbf{x}, \phi) = \phi(\mathbf{x}) (r_i(\mathbf{x}) - \int_{\mathbb{R}^n} \phi(\mathbf{y})d\mathbf{y})$ for any $\mathbf{x} \in \mathbb{R}^n$ and any continuous $L^1(\mathbb{R}^n)$ function $\phi: \mathbb{R}^n \to \mathbb{R}$, that is,

$$ (3) \qquad f_i(\mathbf{x}, u_i(t, \cdot)) = u_i(t, \mathbf{x}) \left( r_i(\mathbf{x}) - \int_{\mathbb{R}^n} u_i(t, \mathbf{y}) d\mathbf{y} \right). $$

When $H = 2$, we established in [12] the existence and uniqueness of the solution of the Cauchy problem associated with (1) with growth functions (3), under some symmetry assumptions on the initial conditions. Then, still in the case $H = 2$, we obtained a characterisation of the large-time behaviour of the solution (persistence or extinction) based on the sign of a principal eigenvalue, here denoted $\lambda_2$.

In this work, we first extend these results to the general case $H \ge 2$. Let us emphasise that, when $H \ge 3$, symmetry arguments are no longer applicable (except in some very particular configurations). In Section 2, we thus state the existence and uniqueness of the solution of the Cauchy problem associated with (1), and we establish a condition for the persistence of the population, based on the sign of a principal eigenvalue denoted $\lambda_H$, and we describe the large-time behaviour of the total population size. Then, we propose in Section 3 a mathematical analysis of the effects of the parameters on the value of $\lambda_H$, including the effect of the mutation, selection and migration parameters. As explained above, our main objective is to investigate the effect of adding a third host, compared to a baseline situation with two hosts. In that respect, we fix the position of the two optima $\mathbf{O}_1$ and $\mathbf{O}_2$, leading to a given value of $\lambda_2$. Then, in Section 4, we study the sign of $\lambda_2 - \lambda_3$, depending on the position $\mathbf{O}_3$ of the optimum of the third host. This sign determines whether the presence of a third host increases ($\lambda_2 - \lambda_3 > 0$) or decreases ($\lambda_2 - \lambda_3 < 0$) the chances of persistence. We propose a mathematical analysis, completed by some numerical simulations. The outcomes are various, and sometimes surprising, see the discussion in Section 5. We gather all the proofs in Section 6.
# 2 Existence, uniqueness and persistence results

## 2.1 The Cauchy problem

For $H \ge 2$, we consider the nonlinear and nonlocal system (1), where the functions $x \mapsto r_i(x)$ are as in (2), while the functions $f_i(x, u_i(t, \cdot))$ are as in (3). System (1) is supplemented with an initial datum

$$ (4) \quad \mathbf{u}^0 = (u_1^0, \dots, u_H^0), \quad u_i^0 \in C(\mathbb{R}^n) \cap L^\infty(\mathbb{R}^n) \cap L^1(\mathbb{R}^n) \text{ and } u_i^0 \ge 0 \text{ for all } 1 \le i \le H, $$

and we focus on the well-posedness of the Cauchy problem (1)-(4), under some additional decay properties on $\mathbf{u}^0$, see (5) below. The proofs are very different from that for the case $H = 2$ in [12], where symmetry arguments were used to reduce the study of the nonlinear system (1) to that of a linear scalar equation. Such arguments cannot be applied in the general case $H > 2$.

**Theorem 2.1** (Well-posedness). Assume that there exist positive constants $K$ and $\theta$ such that the initial condition $\mathbf{u}^0$ in (4) satisfies

$$ (5) \quad \forall 1 \le i \le H, \quad \forall \mathbf{x} \in \mathbb{R}^n, \quad 0 \le u_i^0(\mathbf{x}) \le K e^{-\theta \|\mathbf{x}\|}. $$

Then, there is a unique solution $\mathbf{u} = (u_1, \dots, u_H)$ of (1)-(4) in $C([0, +\infty) \times \mathbb{R}^n, \mathbb{R}^H) \cap C_{t;\mathbf{x}}^{1;2}((0, +\infty) \times \mathbb{R}^n, \mathbb{R}^H)$ such that

$$ (6) \quad \forall 1 \le i \le H, \quad \forall t \ge 0, \quad \forall \mathbf{x} \in \mathbb{R}^n, \quad 0 \le u_i(t, \mathbf{x}) \le K e^{(r_{\max}+1)t-\min(\theta, 1/\mu)\|\mathbf{x}\|}, $$

and the maps $t \mapsto \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, d\mathbf{x}$ are locally Lipschitz-continuous in $[0, +\infty)$.

## 2.2 The principal eigenvalue

Here, we present some linear material, namely the principal eigenvalue $\lambda_H$ and the principal eigenvector $\Phi = (\varphi_1, \dots, \varphi_H)$ solving $\mathcal{A}\Phi^T = \lambda_H\Phi^T$, where the operator

$$ (7) \quad \mathcal{A} := -\frac{\mu^2}{2} \Delta - \begin{pmatrix} r_1(\mathbf{x}) - \delta & \frac{\delta}{H-1} & \cdots & \frac{\delta}{H-1} \\ \frac{\delta}{H-1} & r_2(\mathbf{x}) - \delta & \cdots & \vdots \\ \vdots & \vdots & \ddots & \frac{\delta}{H-1} \\ \frac{\delta}{H-1} & \cdots & \frac{\delta}{H-1} & r_H(\mathbf{x}) - \delta \end{pmatrix} $$

is obtained by linearizing system (1) around the trivial solution $(0, \dots, 0)$. Since $\delta > 0$ and since the fitness functions, defined in (2), satisfy $r_i(x) \to -\infty$ as $\|x\| \to +\infty$, $\mathcal{A}$ can be seen as a cooperative Schrödinger operator with confining potentials, of which the linear analysis is classical, see [1, 7] in a bounded domain, or [5, 6] in a slightly different setting. Another approach, used in [12], consists in defining $\lambda_H$ as the limit, as $R \to +\infty$, of the Dirichlet principal eigenvalue $\lambda_H^R$ in the open Euclidean ball $B(\mathcal{O}, R)$ of center the origin

$$ \mathcal{O} := (0, \dots, 0) $$

and radius $R > 0$. Similarly, $\Phi = (\varphi_1, \dots, \varphi_H)$ can be defined as in [12] as the locally uniform limit of the functions $\Phi^R$ (after suitable multiplication, say, with $\Phi_1^R(\mathcal{O}) = 1$). The functions $\Phi^R$ belong to $C_0^\infty(B(\mathcal{O}, R))^H$ from standard elliptic estimates,
where $C_0^\infty(\overline{B(\mathcal{O}, R)})$ is the space of $C^\infty(\overline{B(\mathcal{O}, R)})$ functions vanishing on $\partial B(\mathcal{O}, R)$. It is also known that the functions $\varphi_i$'s are positive and decay exponentially as $\|\mathbf{x}\| \to +\infty$, thanks to the confining property of the potentials $r_i$. We thus consider the principal eigenvalue $\lambda_H$, the normalised<fnref n="1" /> and positive (that is, positive componentwise) principal eigenvector $\Phi := (\varphi_1, \dots, \varphi_H) \in (C_0^\infty(\mathbb{R}^n) \cap L^1(\mathbb{R}^n))^H$ (with $C_0^\infty(\mathbb{R}^n)$ being the space of $C^\infty(\mathbb{R}^n)$ functions converging to 0 as $\|\mathbf{x}\| \to +\infty$), satisfying

$$ (8) \quad -\frac{\mu^2}{2} \Delta \varphi_i(\mathbf{x}) - r_i(\mathbf{x}) \varphi_i(\mathbf{x}) - \delta \left( \sum_{\substack{k=1 \\ k \neq i}}^{H} \frac{\varphi_k(\mathbf{x})}{H-1} - \varphi_i(\mathbf{x}) \right) = \lambda_H \varphi_i(\mathbf{x}) \quad \text{in } \mathbb{R}^n, \quad 1 \le i \le H. $$

We also know that $\Phi \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$, with

$$ L_w^2(\mathbb{R}^n) := \{ \psi : \mathbb{R}^n \to \mathbb{R} \text{ such that } \mathbf{x} \mapsto \|\mathbf{x}\| \psi(\mathbf{x}) \in L^2(\mathbb{R}^n) \}, $$

and that the following Rayleigh formula is available:

$$ (9) \quad \lambda_H = Q_H(\Phi) = \min \left\{ Q_H(\Psi) : \Psi \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H, \int_{\mathbb{R}^n} \|\Psi(\mathbf{x})\|^2 d\mathbf{x} = 1 \right\}, $$

where

$$ (10) \quad Q_H(\Psi) = Q_H(\psi_1, \dots, \psi_H) := \sum_{i=1}^{H} \left( \frac{\mu^2}{2} \int_{\mathbb{R}^n} \|\nabla \psi_i(\mathbf{x})\|^2 d\mathbf{x} - \int_{\mathbb{R}^n} r_i(\mathbf{x}) (\psi_i(\mathbf{x}))^2 d\mathbf{x} \right) \\ + \delta \left( 1 - \sum_{1 \le i < j \le H} \frac{2}{H-1} \int_{\mathbb{R}^n} \psi_i(\mathbf{x}) \psi_j(\mathbf{x}) d\mathbf{x} \right). $$

Furthermore, the principal eigenvector $\Phi$ and its opposite $-\Phi$ are the unique normalised minima of $Q_H$ in $(H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$, and $\Phi$ is the unique nonnegative (componentwise) eigenvector of $\mathcal{A}$ in $(H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$.

The principal eigenvalue $\lambda_H$ depends on the parameters $\delta$, $\alpha$, $\mu$, $r_{\max}$ (on a trivial additional manner, since $\lambda_H + r_{\max}$ is independent of $r_{\max}$), as well as on the optima $(\mathbf{O}_i)_{1 \le i \le H}$. In the following sections, depending on the context, we will also use the notations $\lambda_H(\delta)$, $\lambda_H(\mathbf{O}_1, \dots, \mathbf{O}_H)$, or $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ in order to emphasise the effect of the various parameters on $\lambda_H$.

For the particular case $H = 1$, we refer to Section 3.1.

## 2.3 Extinction vs persistence

Let $\mathbf{u} = (u_1, \dots, u_H)$ denote the solution of (1)-(6) given in Theorem 2.1. Remember that $u_i \ge 0$ in $[0, +\infty) \times \mathbb{R}^n$, for each $1 \le i \le H$. Assume without loss of generality that $\mathbf{u}^0 \ne (0, \dots, 0)$ in $\mathbb{R}^n$, that is, there is $1 \le j \le H$ such that $u_j \ne 0$ in $\mathbb{R}^n$. The strong parabolic maximum principle applied to $u_j$ implies that $u_j > 0$ in $(0, +\infty) \times \mathbb{R}^n$, and then $u_i > 0$ in $(0, +\infty) \times \mathbb{R}^n$ from the strong parabolic maximum principle applied to each $u_i$ with $i \ne j$. For $t \ge 0$, we define the population size within each host $1 \le i \le H$ as

$$ N_i(t) := \int_{\mathbb{R}^n} u_i(t, \mathbf{y}) \, d\mathbf{y}, $$

<footnote n="1">¹In the sequel, we say that $\Psi \in L^2(\mathbb{R}^n)^H$ is normalised whenever $\int_{\mathbb{R}^n} \|\Psi(\mathbf{x})\|^2 d\mathbf{x} = 1$.</footnote>
which is a positive real number for $t > 0$, by (6) and the previous observations. Furthermore, the Lebesgue dominated convergence theorem implies that each function $N_i$ is continuous in $[0, +\infty)$. The total population size in the system is defined as

$$N(t) := N_1(t) + \dots + N_H(t),$$

and the mean growth rate within each host as

$$(11) \qquad \bar{r}_i(t) := \frac{\int_{\mathbb{R}^n} r_i(\mathbf{y}) u_i(t, \mathbf{y}) \, d\mathbf{y}}{\int_{\mathbb{R}^n} u_i(t, \mathbf{y}) \, d\mathbf{y}} = \frac{\int_{\mathbb{R}^n} r_i(\mathbf{y}) u_i(t, \mathbf{y}) \, d\mathbf{y}}{N_i(t)}.$$

By integrating (11) over $[t_1, t_2] \times B(\mathcal{O}, R)$ for any $0 < t_1 < t_2 < +\infty$ and $R > 0$, by passing to the limit as $R \to +\infty$ and then as $t_2 - t_1 \to 0$, and by using (6) as in the proof of Theorem 2.1 in Section 6.1 (see especially (38) below), it follows that the population sizes are differentiable in $(0, +\infty)$ and satisfy

$$N_i'(t) = \bar{r}_i(t)N_i(t) - N_i(t)^2 + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} (N_k(t) - N_i(t))$$

for every $1 \le i \le H$ and $t > 0$ (and even at $t = 0$ by continuity as $t \to 0$ in the above formula).

The next result shows that the sign of the principal eigenvalue $\lambda_H$ fully determines the fate of the population at large times.

**Proposition 2.2 (Extinction vs persistence).** Let $\lambda_H$ be given by (9), and let $\mathbf{u}$ be the solution of (1)-(6) given by Theorem 2.1, with a non-trivial initial condition $\mathbf{u}^0$, let $N_i(t)$ be the corresponding population sizes in each habitat, and let $N(t)$ the total population size in the system.

(i) If $\lambda_H > 0$ and $\mathbf{u}^0$ is compactly supported, then

$$N(t) \to 0 \text{ as } t \to +\infty \text{ (extinction of the population).}$$

(ii) If $\lambda_H = 0$ and $\mathbf{u}^0$ is compactly supported, then

$$\liminf_{t \to +\infty} \left( \min_{1 \le i \le H} N_i(t) \right) = 0 \text{ (partial extinction of the population).}$$

(iii) If $\lambda_H < 0$, then

$$(12) \quad \limsup_{t \to +\infty} N(t) \ge \limsup_{t \to +\infty} \left( \max_{1 \le i \le H} N_i(t) \right) \ge -\lambda_H > 0 \text{ (persistence of the population).}$$

## 2.4 Stationary states

In this subsection, we are interested in positive bounded stationary states of (1). By positive, we mean positive componentwise in $\mathbb{R}^n$ since, from the strong elliptic maximum principle applied to this cooperative system, if a bounded stationary state is nonnegative componentwise in $\mathbb{R}^n$, it is either positive componentwise, or identically 0 in $\mathbb{R}^n$. From the confining properties of the fitnesses $r_i$, any positive bounded stationary state of (1) necessarily decays to 0 exponentially as $\|\mathbf{x}\| \to +\infty$. Therefore, it follows from Theorem 2.1
and the proof of Proposition 2.2 that (1) has no positive bounded stationary state if $\lambda_H \ge 0$ (see Remark 6.2 below for further details).

Furthermore, it turns out that, when persistence occurs, that is, $\lambda_H < 0$, the nature of the stationary states deeply depends on the number of hosts. Namely, when $H = 1, 2$, the stationary states are proportional to the principal eigenfunctions (and the study of their “shapes” has recently received a lot of attention, see e.g. [2, 9, 18, 19]). This is no longer true in general in the case $H > 2$, due to a possible symmetry breaking. More precisely, we have the following results.

**Proposition 2.3 (Stationary states).** (i) Assume that $H = 1$ and $\lambda_1 < 0$. The positive bounded stationary state of (1) is unique and is given by:

$$p_1(\mathbf{x}) = -\lambda_1 \frac{\varphi_1(\mathbf{x})}{\int_{\mathbb{R}^n} \varphi_1(\mathbf{y}) \, d\mathbf{y}}, \quad \mathbf{x} \in \mathbb{R}^n.$$

(ii) Assume that $H = 2$ and $\lambda_2 < 0$. Then the principal eigenvector $(\varphi_1, \varphi_2)$ of $\mathcal{A}$ in $(H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^2$ satisfies $\int_{\mathbb{R}^n} \varphi_1(\mathbf{y}) \, d\mathbf{y} = \int_{\mathbb{R}^n} \varphi_2(\mathbf{y}) \, d\mathbf{y}$, and

$$(13) \qquad \mathbf{x} \mapsto (p_1(\mathbf{x}), p_2(\mathbf{x})) := \left( -\lambda_2 \frac{\varphi_1(\mathbf{x})}{\int_{\mathbb{R}^n} \varphi_1(\mathbf{y}) \, d\mathbf{y}}, -\lambda_2 \frac{\varphi_2(\mathbf{x})}{\int_{\mathbb{R}^n} \varphi_2(\mathbf{y}) \, d\mathbf{y}} \right)$$

is a positive bounded stationary state of (1) proportional to $(\varphi_1, \varphi_2)$, and satisfying $\int_{\mathbb{R}^n} p_1(\mathbf{x}) \, d\mathbf{x} = \int_{\mathbb{R}^n} p_2(\mathbf{x}) \, d\mathbf{x} = -\lambda_2$.

(iii) Assume that $H > 2$ and $\lambda_H < 0$. Then positive bounded stationary states of (1), if any, are in general not proportional to the principal eigenvector $(\varphi_1, \dots, \varphi_H)$ of $\mathcal{A}$.

# 3 Effect of the main parameters on persistence

In this section, we investigate the dependence of $\lambda_H$ with respect to the parameters $\delta > 0$, $\alpha > 0$, $\mu > 0$, and the optima $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$. It is useful to start with the “reference case”, namely $H = 1$.

## 3.1 The reference case with a single host

In that case of a single host at position $\mathbf{O}_1$, there is no transfer rate, that is, (8) reduces to a single equation with $\delta = 0$. The principal eigenpair $(\lambda_1, G_1)$ solving

$$(14) \qquad -\frac{\mu^2}{2} \Delta G_1(\mathbf{x}) - r_1(\mathbf{x}) G_1(\mathbf{x}) = \lambda_1 G_1(\mathbf{x}) \quad \text{in } \mathbb{R}^n,$$

is explicitly given by

$$(15) \quad \lambda_1 = \lambda_1(\alpha, \mu) = \frac{\mu n \sqrt{\alpha}}{2} - r_{\max}, \quad G_1(\mathbf{x}) = G(\mathbf{x} - \mathbf{O}_1), \quad G(\mathbf{y}) := C_G e^{-\sqrt{\alpha} \|\mathbf{y}\|^2/(2\mu)},$$

where $C_G > 0$ is a normalising positive constant that guarantees that $\int_{\mathbb{R}^n} G^2(\mathbf{y}) d\mathbf{y} = 1$. Notice that, as the problem (14) is left invariant by any translation of the phenotypic space, the eigenvalue $\lambda_1$ does not depend on the position of the optimum $\mathbf{O}_1$.
## 3.2 The case of H hosts with $H \ge 2$

In this subsection, we assume $H \ge 2$ and denote $\lambda_H = \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ the principal eigenvalue of $\mathcal{A}$ in $(H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$, depending on $\delta > 0$, $\alpha > 0$, $\mu > 0$, and $(\mathbf{O}_i)_{1 \le i \le H}$. As above, the corresponding normalised eigenvector is denoted $\Phi = (\varphi_1, \dots, \varphi_H)$. Observe that, if we consider (8) with $\delta = 0$, the H equations are decoupled and we then define

$$ (16) \qquad \lambda_H(0, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) := \lambda_1(\alpha, \mu) = \frac{\mu n \sqrt{\alpha}}{2} - r_{\max}. $$

For later use, we define, for $1 \le i \le H$,

$$ G_i(\mathbf{x}) := G(\mathbf{x} - \mathbf{O}_i) = C_G e^{-\sqrt{\alpha} \|\mathbf{x} - \mathbf{O}_i\|^2 / (2\mu)} $$

and it is convenient to straightforwardly compute

$$ (17) \qquad \int_{\mathbb{R}^n} G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x} = e^{-\sqrt{\alpha} \|\mathbf{O}_i - \mathbf{O}_j\|^2 / (4\mu)}. $$

The following first result of this section asserts that the principal eigenvalue depends continuously on the parameters, and that an increase of the migration rate $\delta > 0$, or of the intensity of selection $\alpha > 0$, or of the mutational parameter $\mu > 0$ reduces the chances of persistence of the population.

**Proposition 3.1** (Continuity and monotonicity w.r.t. parameters). Let $H \ge 2$ be given. The function $(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \mapsto \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is continuous in $[0, +\infty) \times (0, +\infty)^2 \times (\mathbb{R}^n)^H$, and concave with respect to $(\delta, \alpha, \mu) \in [0, +\infty) \times (0, +\infty)^2$. Furthermore, $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is increasing with respect to $\alpha > 0$ and to $\mu > 0$. Lastly, either the optima $\mathbf{O}_i$ are all identical and $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) = \lambda_1(\alpha, \mu) = \mu n \sqrt{\alpha}/2 - r_{\max}$ for all $\delta \ge 0$, or $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is increasing with respect to $\delta \ge 0$.

The following bounds on $\lambda_H$ shows, in particular, that the case with many hosts is always less favourable for the population than the case with a single host.

**Proposition 3.2** (Bounds on $\lambda_H$ w.r.t. $\lambda_1$). Let $H \ge 2$ be given. Then, for all $\delta > 0$, $\alpha > 0$, $\mu > 0$, and $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$,

$$ \begin{align} \frac{\mu n \sqrt{\alpha}}{2} - r_{\max} = \lambda_1(\alpha, \mu) &\le \lambda_1(\alpha, \mu) + \delta \left( 1 - \sum_{1 \le i < j \le H} \frac{2}{H-1} \int_{\mathbb{R}^n} \varphi_i(\mathbf{x}) \varphi_j(\mathbf{x}) d\mathbf{x} \right) \nonumber \\ &\le \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \nonumber \\ (18) \qquad \le \lambda_1(\alpha, \mu) + \delta \left( 1 - \sum_{1 \le i < j \le H} \frac{2}{H(H-1)} \int_{\mathbb{R}^n} G_i(\mathbf{x}) G_j(\mathbf{x}) d\mathbf{x} \right) \nonumber \\ &\le \lambda_1(\alpha, \mu) + \delta = \frac{\mu n \sqrt{\alpha}}{2} - r_{\max} + \delta. \end{align} $$

Moreover, for every $\mu > 0$, $\alpha > 0$ and $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$, $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is bounded independently of $\delta$.

In the sequel, and without loss of generality, we fix the optima $\mathbf{O}_1$ and $\mathbf{O}_2$ at

$$ (19) \qquad \mathbf{O}_1 := (-\beta, 0, \dots, 0), \quad \mathbf{O}_2 := (\beta, 0, \dots, 0), $$
for some $\beta \ge 0$. In the case $H = 2$, we now denote $\lambda_2(\delta, \alpha, \mu, \beta)$ the principal eigenvalue of (8). In [12], we defined the habitat difference by $m_D := \|\mathbf{O}_1 - \mathbf{O}_2\|^2/2 = 2\beta^2$, and we proved that $\lambda_2(\delta, \alpha, \mu, \beta)$ was an increasing function of $m_D$ (and therefore of $\beta$), see [12, Proposition 4]. We also proved that $\lambda_2(\delta, \alpha, \mu, \beta) \to \lambda_1(\alpha, \mu) = \mu n \sqrt{\alpha}/2 - r_{\max}$ as $\delta \to 0^+$. The result below, which is a direct consequence of Proposition 3.2 (together with (17)), extends this last result to the case of $H \ge 2$ hosts.

**Corollary 3.3.** Let $H \ge 2$ be given. Then, for every $\delta > 0$, $\alpha > 0$, $\mu > 0$, $\beta \ge 0$, and $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$ with $\mathbf{O}_1$ and $\mathbf{O}_2$ as in (19), there holds

$$ \lambda_1(\alpha, \mu) \le \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \le \lambda_1(\alpha, \mu) + \delta \left( 1 - \frac{2}{H(H-1)} e^{-\sqrt{\alpha}\beta^2/\mu} \right). $$

Incidentally, for $H = 2$, $\lambda_2(\delta, \alpha, \mu, \beta) - \lambda_1(\alpha, \mu) \to 0$ when $\alpha \to 0^+$, $\mu \to +\infty$ or $\beta \to 0$.

Having understood, by Propositions 3.1-3.2, the small migration rate regime $\delta \to 0$, we now turn to the asymptotics $\delta \to +\infty$. The case $H = 2$ was investigated in [12, Proposition 4], where we obtained that $\lambda_2(\delta, 1, \mu, \beta) \to \lambda_1(1, \mu) + \beta^2/2 = \mu n/2 - r_{\max} + \beta^2/2$ as $\delta \to +\infty$ (in the case $\alpha = 1$). We generalise here this result to any number $H \ge 2$ of hosts.

**Proposition 3.4 (Large migration rate).** Let $H \ge 2$, $\alpha > 0$, $\mu > 0$, and $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$ be given. Then $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \nearrow \lambda_{H,\infty}(\alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ as $\delta \to +\infty$, where $\lambda_{H,\infty}(\alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is the principal eigenvalue for a population living in a single host and with “effective” fitness equal to the mean of the $H$ original fitness functions, that is,

$$ (20) \qquad -\frac{\mu^2}{2} \Delta \varphi(\mathbf{x}) - R_H(\mathbf{x}) \varphi(\mathbf{x}) = \lambda_{H,\infty}(\alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \varphi(\mathbf{x}), \quad R_H(\mathbf{x}) := \frac{1}{H} \sum_{i=1}^{H} r_i(\mathbf{x}), $$

with principal eigenfunction $\varphi \in C_0^\infty(\mathbb{R}^n)$. As a consequence, given the center of mass

$$ \mathbf{M} := \frac{1}{H} \sum_{i=1}^{H} \mathbf{O}_i, $$

we have

$$ (21) \qquad \lambda_{H,\infty}(\alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) = \lambda_1(\alpha, \mu) - R_H(\mathbf{M}) + r_{\max}. $$

# 4 Effect of a third host

As seen in Proposition 3.2, the inequality $\lambda_1(\alpha, \mu) = \lambda_1 \le \lambda_2 = \lambda_2(\delta, \alpha, \mu, \mathbf{O}_1, \mathbf{O}_2)$ is always true (as is the inequality $\lambda_1 \le \lambda_H$ for any $H \ge 2$), meaning that the presence of a second host always penalises the population. In this section, given $\delta > 0$, $\alpha > 0$, $\mu > 0$, we investigate how the introduction of a third host affects the chances of survival of the population, compared to the case of only two hosts. To do so, we denote $\lambda_2 = \lambda_2(\mathbf{O}_1, \mathbf{O}_2)$ the principal eigenvalue for two hosts with optima $\mathbf{O}_1$ and $\mathbf{O}_2$ as in (19), and $\lambda_3 = \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$ the principal eigenvalue for three hosts with optima $\mathbf{O}_1, \mathbf{O}_2$ and $\mathbf{O}_3$. Our goal is to compare $\lambda_2(\mathbf{O}_1, \mathbf{O}_2)$ and $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$. As revealed in the following, the situation is very rich and the outcomes dramatically depend not only on the different parameters but also the position of the third optimum $\mathbf{O}_3$. In the following, we denote

$$ (22) \qquad \Lambda_{\delta} := \{ \mathbf{O}_3 \in \mathbb{R}^n : \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3) < \lambda_2(\mathbf{O}_1, \mathbf{O}_2) \} $$
the region in the phenotypical space where the presence of a third host corresponding to an optimum $O_3$ causes a springboard effect, leading to higher chances of persistence of the pathogen, compared to an environment with the two hosts with optima $O_1$ and $O_2$.

We first deal in Section 4.1 with the case of small or large migration regimes ($\delta$ small or large), while Section 4.2 provides some comparisons which hold for any given $\delta > 0$. Some complementing numerical simulations are presented in Section 4.3.

## 4.1 Large and small migration regimes

### Large migration regime

Let us first consider the regime corresponding to the limit case $\delta \to +\infty$. We denote

$$ (23) \qquad \Lambda_{\infty} := \{O_3 \in \mathbb{R}^n : \lambda_{3,\infty}(O_1, O_2, O_3) < \lambda_{2,\infty}(O_1, O_2)\}, $$

where the quantities $\lambda_{H,\infty}$ are defined in Proposition 3.4 (there, the dependence of $\lambda_{H,\infty}$ on $\alpha$ and $\mu$ was also emphasised; here $\alpha > 0$ and $\mu > 0$ are fixed, and we are mainly interested in the dependence with respect to the position of the optima, this is why we use the notations $\lambda_{3,\infty}(O_1, O_2, O_3)$ and $\lambda_{2,\infty}(O_1, O_2)$). We also recall that $\lambda_1 = \mu n \sqrt{\alpha}/2 - r_{\max}$ is independent of $\delta$.

For two hosts, the center of mass is $M_2 = (O_1 + O_2)/2 = (0, \dots, 0) = O$, and the effective fitness defined in Proposition 3.4 is $R_2(x) = (r_1(x) + r_2(x))/2 = r_{\max} - \alpha \|x\|^2/2 - \alpha \beta^2/2$. As a result

$$ \lambda_{2,\infty}(O_1, O_2) = \lambda_1 - R_2(M_2) + r_{\max} = \lambda_1 + \frac{\alpha \beta^2}{2}. $$

For three hosts, the third optimum is denoted $O_3 = (a_1, \dots, a_n)$, and the effective fitness of Proposition 3.4 is

$$ R_3(x) = \frac{r_1(x) + r_2(x) + r_3(x)}{3} = r_{\max} - \frac{\alpha}{2} \|x - M_3\|^2 - \frac{\alpha}{2} \left( \frac{2}{3}\beta^2 + \frac{2}{9}(a_1^2 + \dots + a_n^2) \right), $$

where $M_3 := \frac{1}{3}O_3$ is the center of mass of the triangle $O_1O_2O_3$. As a result

$$ \lambda_{3,\infty}(O_1, O_2, O_3) = \lambda_1 - R_3(M_3) + r_{\max} = \lambda_1 + \frac{\alpha}{2} \left( \frac{2}{3}\beta^2 + \frac{2}{9}(a_1^2 + \dots + a_n^2) \right). $$

From this, we immediately identify the asymptotic region (as $\delta \to +\infty$) where the third host increases the chances of persistence (see right column of Figure 1).

**Proposition 4.1** (Region leading to higher chances of persistence for large migration rate). In the phenotypical space $\mathbb{R}^n$, the region $\Lambda_{\infty}$ defined by (23) is the open Euclidean ball $B(O, \sqrt{3/2}\beta)$ centered at the origin $O$ and of radius $\sqrt{3/2}\beta$.

Notice that both $O_1$ and $O_2$ belong to $\Lambda_{\infty}$. On the other hand, $O_3 = (0, \pm\sqrt{3}\beta, 0, \dots, 0)$, which makes $O_1O_2O_3$ an equilateral triangle, does not belong to $\Lambda_{\infty}$. This, as revealed below (see Proposition 4.4), is actually true for all $\delta > 0$.

**Remark 4.2.** Consider the “host in the middle” situation, that is $O_3 = O = (0, \dots, 0)$ so that $\lambda_{3,\infty}(O_1, O_2, O) = \lambda_1 + \alpha\beta^2/3$. We now introduce a fourth host whose optimum is denoted $O_4 = (b_1, \dots, b_n)$. As above we compute

$$ \lambda_{4,\infty}(O_1, O_2, O, O_4) = \lambda_1 + \frac{\alpha}{2} \left( \frac{1}{2}\beta^2 + \frac{3}{16}(b_1^2 + \dots + b_n^2) \right). $$

Then the region where $\lambda_{4,\infty}(O_1, O_2, O, O_4) < \lambda_{3,\infty}(O_1, O_2, O)$ is the ball $B(O, \sqrt{8/9}\beta)$ and, this time, $O_1$ and $O_2$ do not belong to this ball.
Small migration regime

Let us fix $\alpha > 0$, $\mu > 0$, $\beta > 0$, $\mathbf{O}_1$ and $\mathbf{O}_2$ as in (19), and let us now investigate $\Lambda_\delta$ given in (22) for $0 < \delta \ll 1$. For given positions $\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3$, using the notations $\lambda_3 = \lambda_3(\delta)$ and $\lambda_2 = \lambda_2(\delta)$, we focus on the difference $\lambda_3(\delta) - \lambda_2(\delta)$. As mentioned above, by (16), $\lambda_3(0) - \lambda_2(0) = \lambda_1 - \lambda_1 = 0$, and $\lambda_3(\delta) - \lambda_2(\delta) \to \lambda_1 - \lambda_1 = 0$ by Proposition 3.2. Hence, if we define $\Lambda_0$ as in (22) with $\delta = 0$, we have $\Lambda_0 = \emptyset$. Thus, the quantity $\lambda'_3(0) - \lambda'_2(0)$, if it exists, would be an important information. For a general $H \ge 2$, we denote $\lambda_H = \lambda_H(\delta)$, all other parameters $\alpha > 0$, $\mu > 0$ and $(\mathbf{O}_i)_{1\le i\le H} \in (\mathbb{R}^n)^H$ being fixed.

**Proposition 4.3** (Region leading to higher chances of persistence for small migration rate). For any $H \ge 2$, the function $\delta \mapsto \lambda_H(\delta)$ is differentiable in $[0, +\infty)$, with $\lambda'_H(\delta) < 1$ for all $\delta \ge 0$. Furthermore, either the optima $\mathbf{O}_i$ are all identical and $\lambda_H(\delta) = \lambda_1 = \mu n\sqrt{\alpha}/2 - r_{\max}$ for all $\delta \ge 0$, or $\lambda'_H(\delta) > 0$ for all $\delta \ge 0$.<fnref n="2" /> Lastly,

$$ (24) \qquad \lambda'_H(0) = \frac{H - \mu_A}{H - 1}, $$

where

$$ \mu_A := \max_{\mathbf{q}=(q_1, \dots, q_H) \in \mathbb{R}^H, \|\mathbf{q}\|=1} \sum_{1 \le i,j \le H} a_{ij} q_i q_j $$

is the largest eigenvalue of the $H \times H$ symmetric matrix $A = (a_{ij})_{1\le i,j\le H}$ defined by

$$ a_{ij} := e^{-\sqrt{\alpha} \|\mathbf{O}_i - \mathbf{O}_j\|^2 / (4\mu)} > 0. $$

In particular, with $\mathbf{O}_1$ and $\mathbf{O}_2$ as in (19), there holds

$$ (25) \qquad \begin{aligned} \mathcal{P} &:= \{\mathbf{O}_3 \in \mathbb{R}^n : \lambda'_3(0) < \lambda'_2(0)\} \\ &= \{\mathbf{O}_3 \in \mathbb{R}^n : k(\mathbf{O}_3) := g(\mathbf{O}_3) \cos\left(\frac{\arccos(h(\mathbf{O}_3) g(\mathbf{O}_3)^{-3})}{3}\right) > 1\} \end{aligned} $$

with

$$ g(\mathbf{O}_3) := \sqrt{\frac{1 + 2h(\mathbf{O}_3) \cosh(\sqrt{\alpha} a_1 \beta / \mu)}{3}}, \quad h(\mathbf{O}_3) := e^{\sqrt{\alpha}(3\beta^2 - \|\mathbf{O}_3\|^2)/(2\mu)}, $$

and $\mathbf{O}_3 = (a_1, \dots, a_n)$. Thus, if $\mathbf{O}_3 \in \mathcal{P}$, then there is $\delta_0 > 0$ such that $\mathbf{O}_3 \in \Lambda_\delta$ for every $0 < \delta < \delta_0$, with $\Lambda_\delta$ as in (22). On the other hand, if $\mathbf{O}_3 \notin \bar{\mathcal{P}}$, then there is $\delta_0 > 0$ such that $\mathbf{O}_3 \notin \Lambda_\delta$ for every $0 < \delta < \delta_0$.

Let us now briefly comment on the properties of $\mathcal{P}$ given in (25). First of all, since the function $k$ is continuous in $\mathbb{R}^n$ and converges to $\cos(\pi/6)/\sqrt{3} = 1/2 < 1$ as $\|\mathbf{O}_3\| \to +\infty$, the set $\mathcal{P}$ is open and bounded. Furthermore, it is symmetric with respect to the hyperplane $\{0\} \times \mathbb{R}^{n-1}$ and axisymmetric with respect to the axis $\mathbb{R} \times \{0\}^{n-1}$. When $\mathbf{O}_3 = (0, \sqrt{3}\beta, 0, \dots, 0)$, then $\mathbf{O}_1\mathbf{O}_2\mathbf{O}_3$ is an equilateral triangle, and $k(\mathbf{O}_3) = 1$, hence $\mathbf{O}_3 \in \partial\mathcal{P}$. We will actually prove in Proposition 4.4 below that, for any $\delta > 0$, $\mathbf{O}_3 = (0, \sqrt{3}\beta, 0, \dots, 0) \notin \Lambda_\delta$. But the previous observation shows that, for small $\delta$, it was so close!

On the other hand, for an arbitrary point $\mathbf{O}_3 = (a_1, \dots, a_n) \in \mathbb{R}^n$, formula (24) implies that $\lambda'_2(0) = 1 - e^{-\sqrt{\alpha}\beta^2/\mu}$ and, with $\mathbf{q} := (1/\sqrt{3}, 1/\sqrt{3}, 1/\sqrt{3})$,

$$ \lambda'_3(0) \le \frac{3 - \mathbf{q}A\mathbf{q}}{2} = 1 - \frac{e^{-\sqrt{\alpha}\beta^2/\mu} + 2e^{-\sqrt{\alpha}(\beta^2 + \|\mathbf{O}_3\|^2)/(4\mu)} \cosh(\sqrt{\alpha}a_1\beta/(2\mu))}{3}. $$

<footnote n="2">²Notice that the inequalities $0 \le \lambda'_H(0) < 1$ are coherent with Corollary 3.3.</footnote>
Therefore,

$$ \mathcal{P} \supset \{ \mathbf{O}_3 \in \mathbb{R}^n : e^{-3\sqrt{\alpha\beta^2}/(4\mu)} < e^{-\sqrt{\alpha}a_1^2/(4\mu)} \cosh\left(\frac{\sqrt{\alpha}a_1\beta}{2\mu}\right) e^{-\sqrt{\alpha}(a_2^2+\cdots+a_n^2)/(4\mu)} \} $$

By studying the variations of the function $a_1 \mapsto e^{-\sqrt{\alpha}a_1^2/(4\mu)} \cosh(\sqrt{\alpha}a_1\beta/(2\mu))$, one infers that

$$ \mathcal{P} \supset [-\sqrt{3}\beta, \sqrt{3}\beta] \times \{0\}^{n-1}, $$

hence $\mathbf{O}_1 \in \mathcal{P}$ and $\mathbf{O}_2 \in \mathcal{P}$. Lastly,

$$ \mathcal{P} \supset \{0\} \times B'(O', \sqrt{3}\beta), $$

where $B'(O', \sqrt{3}\beta)$ is the open Euclidean ball of radius $\sqrt{3}\beta$ and centered at the origin $O' := (0, \dots, 0)$ in $\mathbb{R}^{n-1}$.

## 4.2 The general case $\delta > 0$

In this subsection, we fix the parameters $\delta > 0$, $\alpha > 0$, and $\mu > 0$. We recall that, without loss of generality, $\mathbf{O}_1$ and $\mathbf{O}_2$ are as in (19) with $\beta \ge 0$, and that $\Lambda_\delta$ denotes the phenotypical region where the presence of $\mathbf{O}_3$ leads to higher chances of persistence, as defined in (22).

### Some fitness loss situations

We show that if the third optimum $\mathbf{O}_3$ is such that $\mathbf{O}_1\mathbf{O}_2\mathbf{O}_3$ forms an equilateral triangle, then $\lambda_2 - \lambda_3$ is always negative (reduced chances of persistence). Since, by uniqueness of the principal eigenpair of (8), the map $\mathbf{O}_3 \mapsto \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$ is invariant by rotation around the line $(\mathbf{O}_1\mathbf{O}_2)$ containing $\mathbf{O}_1$ and $\mathbf{O}_2$, we can assume without loss of generality that $\mathbf{O}_3$ lies in $\mathbb{R}^2 \times \{0\}^{n-2}$, that is, $\mathbf{O}_3 = (0, \pm\sqrt{3}\beta, 0, \dots, 0)$.

**Proposition 4.4 (Equilateral triangle configuration).** Assume that $\delta > 0$, $\alpha > 0$, $\mu > 0$ are fixed, together with $\mathbf{O}_1$ and $\mathbf{O}_2$ as in (19) with $\beta > 0$. If $\mathbf{O}_3 = (0, \sqrt{3}\beta, 0, \dots, 0)$ or $\mathbf{O}_3 = (0, -\sqrt{3}\beta, 0, \dots, 0)$, then $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3) > \lambda_2(\mathbf{O}_1, \mathbf{O}_2)$, hence $\mathbf{O}_3 \notin \Lambda_\delta$.

Now, still for any fixed $\delta > 0$, $\alpha > 0$, $\mu > 0$, and $\mathbf{O}_1, \mathbf{O}_2$ as in (19) with $\beta \ge 0$, we aim at describing the behaviour of $\lambda_3$ when the third optimum $\mathbf{O}_3$ is far away from the two other optima. To do so, we introduce an auxiliary eigenvalue problem with two hosts, corresponding to a situation with migration loss: half of the individuals are lost during the migration process. Namely, we define the principal eigenvalue $\tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2)$ satisfying

$$ (26) \qquad -\frac{\mu^2}{2} \Delta \tilde{\varphi}_i(\mathbf{x}) - r_i(\mathbf{x}) \tilde{\varphi}_i - \delta \left( \frac{1}{2} \tilde{\varphi}_j(\mathbf{x}) - \tilde{\varphi}_i(\mathbf{x}) \right) = \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2) \tilde{\varphi}_i, $$

for $(i, j) = (1, 2)$ and $(i, j) = (2, 1)$ with $(\tilde{\varphi}_1, \tilde{\varphi}_2) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n) \cap C_0^\infty(\mathbb{R}^n) \cap L^1(\mathbb{R}^n))^2$ the corresponding normalised principal eigenvector. Note the factor $1/2$ in front of $\tilde{\varphi}_j$ which describes the migration loss.

Before going further on, we note that an evaluation of the Rayleigh quotient associated with $\tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2)$ at $(G_1, 0)$ (with $G_1 = G(\cdot - \mathbf{O}_1)$ the principal eigenfunction associated with $\lambda_1$) implies that $\lambda_2(\mathbf{O}_1, \mathbf{O}_2) \le \lambda_1 + \delta$. Additionally, multiplying (26) by $\tilde{\varphi}_1$ for $(i, j) = (1, 2)$ (resp. by $\tilde{\varphi}_2$ for $(i, j) = (2, 1)$), integrating by parts and adding the two equations, and using this time the Rayleigh formula associated with $\lambda_2(\mathbf{O}_1, \mathbf{O}_2)$, we observe
that $\tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2) \ge \lambda_2(\mathbf{O}_1, \mathbf{O}_2) + \delta \int_{\mathbb{R}^n} \tilde{\varphi}_1 \tilde{\varphi}_2 > \lambda_2(\mathbf{O}_1, \mathbf{O}_2)$ since $\delta > 0$ and $\tilde{\varphi}_1 > 0, \tilde{\varphi}_2 > 0$ in $\mathbb{R}^n$. Finally,

$$ (27) \qquad \lambda_2(\mathbf{O}_1, \mathbf{O}_2) < \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2) \le \lambda_1 + \delta. $$

When the third optimum “tends to infinity”, we expect the system to resemble this case of two hosts with migration loss. Indeed, roughly speaking, the third host disappears from the system taking with him migrants arriving from hosts 1 and 2, thus acting as a well or a black hole. Precisely, the following holds.

**Proposition 4.5** (Far away third optimum). Assume that $\delta > 0, \alpha > 0, \mu > 0$ are fixed, together with $\mathbf{O}_1$ and $\mathbf{O}_2$ as in (19) with $\beta \ge 0$. For $\|\mathbf{O}_3\|$ large enough, there holds

$$ (28) \qquad \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2) \ge \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3) \ge \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2) - \frac{\delta \sqrt{48 (\tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2) + r_{\max})}}{\sqrt{\alpha} \min(\|\mathbf{O}_3 - \mathbf{O}_1\|, \|\mathbf{O}_3 - \mathbf{O}_2\|)}, $$

so that $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3) \to \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2)$ as $\|\mathbf{O}_3\| \to +\infty$, and $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3) > \lambda_2(\mathbf{O}_1, \mathbf{O}_2)$ for all $\|\mathbf{O}_3\|$ large enough.

**Some situations where the presence of a third host leads to higher chances of persistence**

In Section 4.1, the set $\Lambda_\delta$ where the presence of $\mathbf{O}_3$ leads to higher chances of persistence, defined in (22), was captured for both regimes $\delta \to +\infty$ and $\delta \to 0^+$. For arbitrary $\delta > 0$, a sharp characterisation is very involved (see Proposition 4.7 and the numerical simulations below). Nevertheless, we can provide the following information, saying that a third host with phenotypic optimum close to that of one of the other two hosts is globally beneficial for the chances of persistence.

**Proposition 4.6** (When the third host resembles one of the two others). Assume that $\delta > 0, \alpha > 0, \mu > 0$ are fixed, together with $\mathbf{O}_1$ and $\mathbf{O}_2$ as in (19) with $\beta > 0$. Then there is $\rho > 0$ such that

$$ B(\mathbf{O}_1, \rho) \cup B(\mathbf{O}_2, \rho) \subset \Lambda_\delta. $$

**A host in the middle**

Here we consider a different approach: we fix the parameters $\delta > 0, \alpha > 0, \mu > 0$ and the third optimum $\mathbf{O}_3$ as $(\mathbf{O}_1 + \mathbf{O}_2)/2$. In other words, $\mathbf{O}_3 = (0, \dots, 0) = \mathcal{O}$, with our hypothesis (19). We show that, when the two optima $\mathbf{O}_1$ and $\mathbf{O}_2$ are far enough from each other (meaning $\beta > 0$ large), this “host in the middle” configuration does not maximise the difference $\lambda_2(\mathbf{O}_1, \mathbf{O}_2) - \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$, that is, it does not minimise $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$. This result, consistent with the simulations in Figure 1, is slightly surprising (at least at first glance) and shows the richness of the outcomes.

**Proposition 4.7** (Host in the middle). Assume that $\delta > 0, \alpha > 0$ and $\mu > 0$ are fixed. For $\beta > 0$ large enough in (19), we have

$$ \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathcal{O}) > \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_1). $$

More precisely, $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathcal{O}) \to \lambda_1 + \delta$ while $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_1) \to \lambda_1 + \delta/2$ as $\beta \to +\infty$.

We point out that, since $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_1) = \lambda_3(\mathbf{O}_1, \mathbf{O}_1, \mathbf{O}_2)$ by symmetry, the limit $\lambda_1 + \delta/2 = \lim_{\beta \to +\infty} \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_1) = \lim_{\beta \to +\infty} \lambda_3(\mathbf{O}_1, \mathbf{O}_1, \mathbf{O}_2)$ is coherent with the limit $\lim_{\beta \to +\infty} \lambda_3(\mathbf{O}_1, \mathbf{O}_1, \mathbf{O}_2) = \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_1) = \lambda_1 + \delta/2$, with $\tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_1)$ being the principal eigenvalue of (26) in the case of identical optima (see the proof of Proposition 4.7 in Section 6.3 for further details).
Decreasing $\lambda_3$ by projecting

Next, still with $O_1$ and $O_2$ as in (19) without loss of generality, we show that replacing $O_3$ by its projection on the line $\mathbb{R} \times \{0\}^{n-1}$ containing $O_1$ and $O_2$ always decreases the principal eigenvalue $\lambda_3$, and thus increases the difference $\lambda_2 - \lambda_3$, leading to higher chances of persistence.

**Proposition 4.8** (Better on the line). Let $O_3 \in \mathbb{R}^n$, and let $O_3^\sharp$ be its projection on the line $\mathbb{R} \times \{0\}^{n-1}$, with $O_1$ and $O_2$ as in (19). Then

$$ \lambda_3(O_1, O_2, O_3) \geq \lambda_3(O_1, O_2, O_3^\sharp). $$

Existence of a best third optimum

With $O_1$ and $O_2$ being fixed, the following corollary asserts the existence of a minimum for the function $O_3 \mapsto \lambda_3(O_1, O_2, O_3)$ defined in $\mathbb{R}^n$, thus ensuring the existence of a point in the phenotypic space maximizing the difference $\lambda_2(O_1, O_2) - \lambda_3(O_1, O_2, O_3)$.

**Corollary 4.9** (Minimizing $\lambda_3(O_1, O_2, \cdot)$). Assume that $\delta > 0$, $\alpha > 0$, $\mu > 0$, $O_1 \in \mathbb{R}^n$ and $O_2 \in \mathbb{R}^n$ are fixed. Then the function $O_3 \mapsto \lambda_3(O_1, O_2, O_3)$ has a minimum in $\mathbb{R}^n$, equal to $O_1 = O_2$ if these two optima are identical, or lying in the line $(O_1 O_2)$ if $O_1 \neq O_2$.

*Proof.* Without loss of generality, one can assume that $O_1$ and $O_2$ are given as in (19). On the one hand, if $\beta = 0$, then $O_1 = O_2 = \mathcal{O}$ and $\lambda_3(O_1, O_2, \mathcal{O}) = \lambda_2(O_1, O_2) = \lambda_1$, by Proposition 3.1, while $\lambda_3(O_1, O_2, O_3) > \lambda_1$ for all $O_3 \neq \mathcal{O}$, by Proposition 3.1 again. The conclusion of Corollary 4.9 is then immediate in this case. On the other hand, if $\beta > 0$, the conclusion easily follows from the continuity property in Proposition 3.1, together with Propositions 4.5, 4.6 and 4.8. $\square$

**Remark 4.10.** From Proposition 4.7, it follows that the middle optimum $(O_1 + O_2)/2$ is not a minimum of the map $O_3 \mapsto \lambda_3(O_1, O_2, O_3)$ when $\|O_1 - O_2\|$ is large enough. Hence, for all $\|O_1 - O_2\|$ large enough, by symmetry, the map $\lambda_3(O_1, O_2, \cdot)$ has at least two different minima, of the type $(\pm a_\beta, 0, \dots, 0)$, for a certain $a_\beta > 0$, under the notation (19). For all $\|O_1 - O_2\|$ large enough, these optima maximise the difference $\lambda_2(O_1, O_2) - \lambda_3(O_1, O_2, \cdot)$ (which is positive, by Proposition 4.6).

4.3 Numerical simulations

We present in Figure 1 a map of the positions of the third optimum $O_3$ in the phenotypic plane $\mathbb{R}^2$ with a colorbar picturing the corresponding difference $\lambda_2(O_1, O_2) - \lambda_3(O_1, O_2, O_3)$, with $O_1$ and $O_2$ as in (19). We recall that the difference $\lambda_2(O_1, O_2) - \lambda_3(O_1, O_2, O_3)$ measures the gain (or loss, if negative) in the chances of persistence, when the third host is added to the system. All Matlab source code used to generate the analyses performed here is available at https://doi.org/10.17605/OSF.IO/QAV2M.

This figure serves as an illustration of the different results of this section. In each figure, the set $\Lambda_\delta$, defined in (22), corresponds to the green regions. The last column illustrates the behaviour $\delta \to +\infty$ of Proposition 4.1. In particular, on the first line, for $\beta = 1/2$, there is a perfect match between the ball given in the asymptotic regime $\delta \to +\infty$ and the numerical simulations. For the other values of parameter $\beta$ pictured here, $\delta = 10$ seems not large enough to be considered infinite, but one can still observe on the last two lines, going from the left to the right a “convexification” of the set $\Lambda_\delta$, for the lack of the better word, that points in the direction of the convergence towards the eventual ball, denoted $\Lambda_\infty$ in Proposition 4.1.
$\delta = 0.01$

$\delta = 1$

$\delta = 10$

<figure><img src="image_3.png" /><figcaption>Figure 1: **Effect of adding a third host on persistence.** Each point represent the position of a third optimum, while the two other optima are fixed, at the respective positions $(-\beta, 0)$ and $(\beta, 0)$, pictured in blue. The colormap corresponds to the value taken by the difference $\lambda_2 - \lambda_3$, depending on the position of the third optimum. Each row corresponds to a different value of $\delta$ and each line to a different value of $\beta$. We assumed here that $\mu^2/2 = \alpha = 1$.</figcaption></figure>
On the other hand, the first column represents the regime $\delta \to 0^+$, detailed in Proposition 4.3. In particular, it can be observed that, for small $\beta$, the set $\Lambda_\delta$ has the shape of an “athletic stadium”. Then, for larger $\beta$, the set $\Lambda_\delta$ becomes “oxygen molecule shaped”, and thus no longer convex. Interestingly, increasing $\beta$ leads to a loss of connectedness of the set $\Lambda_\delta$, see the picture for $\delta = 0.01$ and $\beta = 2$. In particular, in agreement with Proposition 4.7, this shows that the central point $(0,0)$ is not optimal in terms of persistence.

Actually, this (numerically observed) loss of connectedness suggests that the introduction of the third host “in the middle” may lead to a decrease in pathogen fitness, compared to the case of only two hosts. One should be careful that this rather surprising result is not implied by Proposition 4.7, which compares the introduction of a third host having an optimum in the middle of the first two hosts with that of a third host having the same optimum as one of the first two hosts. Notice that, as seen in the computations after Proposition 4.3, $O_3 = (0,0)$ always belongs to $\mathcal{P}$ defined in (25). However, from the numerical simulations, it is natural to conjecture that $(0,0) \notin \Lambda_\delta$ should occur for some “particular combinations of large $\beta$-small $\delta$” that should depend on $\mu$ (mutations) and $\alpha$ (selection). This makes the proof of this conjecture very challenging.

Finally, we can verify that, as proved in Proposition 4.6, if the third optimum is very close (or indeed at the same position) than one of the other two optima, then the fitness of the system with three optima is better than the one with two. This corresponds to the concept of “evolutionary springboard”, as adaptation to one optimum is favoured by the presence of another one nearby. Incidentally, in accordance with Proposition 4.8, the closer the third optimum to the axis $(O_1O_2)$, the larger the chances of persistence. This means that on the $x_2$-axis the third optimum is best positioned at the same phenotypic height than the two fixed optima. On the other axis (the $x_1$-axis), despite the symmetry of the two fixed optima, it is much harder to get qualitative information as exemplified above by the case where the third host is in the middle.

# 5 Discussion

We proposed a model to describe the adaptation dynamics of a phenotypically structured population in a $H$-patch environment, with each patch associated with a different phenotypic optimum, and performed a rigorous mathematical analysis of this model. This analysis sheds new light on the effect of increasing the number of hosts on the persistence of a (say, pathogenic) population. As is already known with two hosts, a pathogen population will always have more difficulty persisting on multiple hosts connected by migration than on a single host (Proposition 3.2). Furthermore, increasing the migration rate $\delta > 0$ further reduces the chances of population persistence (Proposition 3.1). However, in some configurations where pathogen persistence is not possible with two hosts, adding a third host achieves persistence. In such case, the presence of a third host causes a springboard effect, by modifying the geometric configuration of the phenotypic space. Thus, as we already know experimentally with the examples of zoonoses mentioned in the Introduction, the increase in the number of hosts does not necessarily mean a decrease in the chances of persistence as long as the initial number of hosts is at least equal to 2. However, the occurrence of this springboard effect, or on the contrary of a detrimental effect of increasing the number of hosts on the persistence of the pathogen, depends in a rather complex way on the respective positions of the optimal phenotypes associated with each host. Our results, which are based on a fixed point theorem, comparison principles, integral estimates, variational arguments and rearrangement techniques, provide a better understanding of these dependencies.
In the asymptotic regimes of small and large migration rates, we were able to charac-
terise very precisely the situations where the addition of a third host increases or decreases
the chances of persistence of the pathogen, compared to a situation with only two hosts.
Proposition 3.4 shows that when the migration rate is large, the system with $H$ host be-
haves as a single population, living in a single host and with “effective” fitness equal to the
mean of the $H$ original fitness functions. In this case, where the population can be quali-
fied as “generalist”, we were able to show that the addition of a third host is favourable to
persistence if and only if the optimum associated with this host is inside a certain ball con-
taining the two other optima (Proposition 4.1). When the migration rate is low, as shown
in Proposition 4.3 (and Figure 1), the shape of the region $\Lambda_\delta$ in the phenotypical space
where the presence of a third host leads to higher chances of persistence of the pathogen is
more complex and depends on the parameters. In particular, we see that this region is not
necessarily convex especially if the two initial optima are far from each other. Moreover,
while the addition of a third host always decreases the chances of persistence when the
optima associated with the three hosts form an equilateral triangle in the phenotypic space
(Proposition 4.4), we observe that when migration is low this equilateral configuration is
very close to the $\Lambda_\delta$ region. Of course, some of these observations are only possible if the
dimension of the phenotype space is at least 2. This illustrates the importance of taking
pleiotropy into account (in the sense that a single mutation can affect multiple pheno-
types), as already noted in [24] in the case $H = 1$ (but with other assumptions on the
mutational operator).

In the general case (i.e., with an arbitrary migration rate $\delta$), and when the positions $\mathbf{O}_1$ and $\mathbf{O}_2$ of the two other hosts are fixed, we proved the existence of a position in the phenotypic space of the optimum associated with the third host maximizing the chances of pathogen persistence (Corollary 4.9). This position is necessarily on the line ($\mathbf{O}_1 \mathbf{O}_2$) (Proposition 4.8). However, contrary to a natural intuition, it is not necessarily located in the exact mean position between $\mathbf{O}_1$ and $\mathbf{O}_2$ (Proposition 4.7). In particular, if the two initial optima $\mathbf{O}_1$ and $\mathbf{O}_2$ are far apart, it is more beneficial for the persistence of the pathogen to introduce a third host associated with an optimum $\mathbf{O}_3$ that is in a neighbourhood of $\mathbf{O}_1$ or $\mathbf{O}_2$. This result which is consistent with the simulations in Figure 1 shows the richness of the outcomes.

From a mathematical point of view, an important part of the difficulty in dealing with $H \ge 3$ hosts, compared to $H = 1$ or $H = 2$ hosts, comes from the lack of symmetry. In the $H = 2$ case, provided that the initial conditions are themselves symmetric about $(\mathbf{O}_1 + \mathbf{O}_2)/2$, the solutions $u_1$ and $u_2$ remain symmetric about this point. This type of symmetry argument was used in [12] to reduce the study of the nonlinear system (1) to that of a linear scalar equation. Such arguments cannot be applied in the general case $H \ge 3$ where we had to treat the complete nonlinear and nonlocal system. This has led to very different proofs, starting with the well-posedness of the Cauchy problem (Theorem 2.1). The stationary states of the system illustrate well the complexity induced by the absence of symmetry: when there is persistence and $H = 1, 2$, the stationary states are proportional to the principal eigenfunctions, but this is no longer true in general in the case $H \ge 3$, because of a possible symmetry breaking (Proposition 2.3). Again because of the lack of symmetry, the populations at equilibrium have no reason to be equal in each host. Even when $H = 2$, the eigenfunctions are not radially symmetric (they are biased towards the other optimum, see formula 46 in [12]). When $H \ge 3$, the qualitative description of their profile is even more involved, leading to new difficulties in the comparison between different geometric configurations of the positions of the optima. Although we were able to show that the “best position” of $\mathbf{O}_3$ when $\mathbf{O}_1$ and $\mathbf{O}_2$ are fixed should belong to the
line $(O_1 O_2)$, the rearrangement techniques that we used for this result did not allow to determine the best position exactly. This is because the rearrangement result is only true in the directions where $O_1$ and $O_2$ are stable by rearrangement (and thus not in the direction $x_1$ with our conventions; otherwise the position $O_3 = (O_1 + O_2)/2$ would always be optimal, which would contradict the result of the Proposition 4.7). Note that the rearrangement inequalities used in the proof of Proposition 4.8 are not standard, because of the unbounded domain and of the sign changes of $r_i$.

As mentioned in the Introduction, one of the fundamental principles in agroecology is to promote diversified agroecosystems rather than uniform cultures. Our results show, with a rigorous mathematical analysis, that what is true when moving from a monoculture to a varietal mixture is not always true when moving from a varietal mixture with $H \ge 2$ species to a mixture with $H + 1$ species. In such a situation, it is necessary to know the respective positions of the optima associated with the different hosts, from the point of view of the pathogen, to predict the positive or negative impact of diversification. From an experimental point of view, it would therefore be necessary to be able to estimate these positions, by comparing, for example, via cross-inoculation experiments, the effect of pre-adaptation of a pathogen on each host $i$ before its introduction on a host $j$, for each pair $i, j$. Of course, the vision here, where each host is described by a unique optimum (FGM model) and where the different parameters do not depend on the host $i$ is very schematic, and the comparison with quantitative data would require more complex approaches. However, the model (1) has the advantage of being mathematically tractable while capturing the antagonistic effects of host diversification.

# 6 Proofs of the main results

This section is devoted to the proofs of all main results of the paper which have not been shown along the way in the previous sections. The results of Section 2 are proved in Section 6.1, while Section 6.2 contains the proofs of the results of Section 3, and Section 6.3 those of Section 4.

## 6.1 Proofs of Theorem 2.1 and Propositions 2.2-2.3

*Proof of Theorem 2.1.* Our strategy involves firstly a fixed point argument to construct a unique local-in-time solution, and secondly the derivation of some estimates to conclude that the solution is defined for all time $t \ge 0$. One recalls that $u^0 = (u_1^0, \dots, u_H^0)$ satisfies (4)-(5). Having in mind the assumption (5) and the desired inequality (6), we can assume without loss of generality that

$$0 < \theta \le \frac{1}{\mu}$$

in (5). Set $r_{\max}^+ := \max(r_{\max}, 0)$, let $\omega_n$ be the $(n-1)$-dimensional Lebesgue measure of the unit Euclidean sphere of $\mathbb{R}^n$, let $B := \max_{1 \le i \le H} \|O_i\|$ and

$$(29) \quad C_1 := \frac{e}{3} \theta^{-n} \omega_n(n-1)! \times [|r_{\max}| + \alpha B^2 + \alpha \theta^{-2} n(n+1) + K e^{r_{\max}^+} \theta^{-n} \omega_n(n-1)!(e+3)+2\delta].$$

Denote

$$(30) \quad T_1 := \left(1 + \frac{\theta^2 \mu^2}{2} + r_{\max}^+ + 3 \theta^{-n} \omega_n (n-1)! K e^{r_{\max}^+ + 1} + 9 C_1^2\right)^{-1} \in (0, 1),$$
and

$$
\begin{equation} \tag{31}
\begin{aligned}
E &:= \left\{ \mathbf{u} = (u_1, \dots, u_H) \in C([0, T_1] \times \mathbb{R}^n, \mathbb{R}^H) \text{ such that } u_i(0, \cdot) = u_i^0 \text{ in } \mathbb{R}^n, \right. \\
  &\qquad \left. 0 \le u_i(t, x) \le K e^{(r_{\max}+1)t-\theta\|\mathbf{x}\|} \text{ for all } (t, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n, \right. \\
  &\qquad \left. \text{and } t \mapsto \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, \mathrm{d}\mathbf{x} \in C^{0,1/2}([0, T_1]), \text{ for all } 1 \le i \le H \right\}.
\end{aligned}
\end{equation}
$$

We equip E with the distance

$$
\begin{equation} \tag{32}
\begin{aligned}
d_E(\mathbf{u}, \mathbf{v}) &:= \max \left( \max_{1 \le i \le H} \|(t, \mathbf{x}) \mapsto (u_i(t, \mathbf{x}) - v_i(t, \mathbf{x})) e^{\theta \|\mathbf{x}\|}\| \right)_{L^\infty([0, T_1] \times \mathbb{R}^n)}, \\
&\quad \max_{1 \le i \le H} \left\| t \mapsto \int_{\mathbb{R}^n} (u_i(t, \mathbf{x}) - v_i(t, \mathbf{x})) \mathrm{d}\mathbf{x} \right\|_{C^{0,1/2}([0, T_1])},
\end{aligned}
\end{equation}
$$

where $\|g\|_{C^{0,1/2}([0,T_1])} = \max_{t \in [0,T_1]} |g(t)| + \sup_{t' \neq t \in [0,T_1]} |g(t) - g(t')|/|t - t'|^{1/2}$ for $g \in C^{0,1/2}([0,T_1])$. The space $(E, d_E)$ is a non-empty complete metric space.

Define now a mapping $\mathcal{F}$ as follows:

$$
\begin{array}{rcl}
\mathcal{F} : & E & \to E \\
& \mathbf{u} & \mapsto \mathcal{F}(\mathbf{u}) := \mathbf{U},
\end{array}
$$

where, for a given $\mathbf{u} \in E$, $\mathbf{U} = (U_1, \dots, U_H) \in C([0, T_1] \times \mathbb{R}^n, \mathbb{R}^H) \cap C_{t;\mathbf{x}}^{1/2}((0, T_1] \times \mathbb{R}^n, \mathbb{R}^H)$
is the classical bounded nonnegative solution of the linear system

$$
\partial_t U_i(t, \mathbf{x}) = \frac{\mu^2}{2} \Delta U_i(t, \mathbf{x}) + \left( r_i(\mathbf{x}) - \int_{\mathbb{R}^n} u_i(t, \mathbf{y}) \, \mathrm{d}\mathbf{y} \right) U_i(t, \mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} (U_k(t, \mathbf{x}) - U_i(t, \mathbf{x}))
$$

in $(0, T_1] \times \mathbb{R}^n$, for $1 \le i \le H$, with initial condition $U_i(0, \cdot) = u_i^0$ in $\mathbb{R}^n$. Using the Hölder
continuity of the maps $t \mapsto \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, \mathrm{d}\mathbf{x}$ in $[0, T_1]$, such a solution $\mathbf{U}$ exists and is unique,
by [3].

To show that $\mathcal{F}(E) \subset E$, let us first check that, for each $\mathbf{u} \in E$, the image $\mathbf{U} = \mathcal{F}(\mathbf{u})$
satisfies the pointwise bounds of $E$. Notice that

$$
(33) \quad \partial_t U_i(t, \mathbf{x}) \le \frac{\mu^2}{2} \Delta U_i(t, \mathbf{x}) + r_{\max} U_i(t, \mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} (U_k(t, \mathbf{x}) - U_i(t, \mathbf{x}))
$$

for every $1 \le i \le H$ and $(t, \mathbf{x}) \in (0, T_1] \times \mathbb{R}^n$. Take any vector $\xi \in \mathbb{R}^n$ such that $\|\xi\| = 1$
and define

$$
(34) \quad \bar{U}_i(t, \mathbf{x}) := K e^{(r_{\max}+1)t - \theta \xi \cdot \mathbf{x}}, \quad 1 \le i \le H, \ t \in \mathbb{R}, \ \mathbf{x} \in \mathbb{R}^n.
$$

Notice that the positive functions $\bar{U}_i$ actually do not depend on $i$, and satisfy

$$
(35) \quad \partial_t \bar{U}_i(t, \mathbf{x}) = \frac{\mu^2}{2} \Delta \bar{U}_i(t, \mathbf{x}) + (r_{\max} + 1) \bar{U}_i(t, \mathbf{x}) - \frac{\theta^2 \mu^2}{2} \bar{U}_i(t, \mathbf{x})
$$

in $\mathbb{R} \times \mathbb{R}^n$. Since we had assumed that $0 < \theta \le 1/\mu$ without loss of generality, we get that

$$
(36) \quad \partial_t \bar{U}_i(t, \mathbf{x}) \ge \frac{\mu^2}{2} \Delta \bar{U}_i(t, \mathbf{x}) + r_{\max} \bar{U}_i(t, \mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} (\bar{U}_k(t, \mathbf{x}) - \bar{U}_i(t, \mathbf{x}))
$$
for every $1 \le i \le H$ and $(t, \mathbf{x}) \in \mathbb{R} \times \mathbb{R}^n$ (recall that $\bar{U}_k - \bar{U}_i \equiv 0$). Moreover, we know from (5) that $U_i(0, \cdot) = u_i^0 \le \bar{U}_i(0, \cdot)$ in $\mathbb{R}^n$. Therefore, by (33) and (36), $\mathbf{U}$ and $\bar{\mathbf{U}} := (\bar{U}_1, \dots, \bar{U}_H)$ are sub- and super-solutions of a linear cooperative system, and the parabolic maximum principle implies that $U_i(t, \mathbf{x}) \le \bar{U}_i(t, \mathbf{x})$ for all $1 \le i \le H$ and $(t, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n$, that is, $U_i(t, \mathbf{x}) \le K e^{(r_{\max}+1)t-\theta\xi \cdot \mathbf{x}}$. Since $\xi$ was arbitrary in the Euclidean unit sphere of $\mathbb{R}^n$, and remembering the nonnegativity of each $U_i$, we get that

$$ (37) \qquad 0 \le U_i(t, \mathbf{x}) \le K e^{(r_{\max}+1)t-\theta\|\mathbf{x}\|} $$

for all $1 \le i \le H$ and $(t, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n$. The Lebesgue dominated convergence theorem then implies that each map $t \mapsto \int_{\mathbb{R}^n} U_i(t, \mathbf{x})\mathbf{d}\mathbf{x}$ is well defined and continuous in $[0, T_1]$. It remains to show that each such map is in $C^{0,1/2}([0, T_1])$. Since the functions $U_i$ and $r_i U_i$ decay exponentially to 0 as $\|\mathbf{x}\| \to +\infty$ uniformly in $t \in [0, T_1]$, it follows from standard parabolic estimates that, for each $0 < a \le b \le T_1$, the functions $\partial_{x_j} U_i$ decay exponentially to 0 as $\|\mathbf{x}\| \to +\infty$ uniformly in $t \in [a, b]$, for every $1 \le i \le H$ and $1 \le j \le n$. Therefore, for every $0 < t \le t' \le T_1$, by integrating the equation satisfied by $U_i$ over $[t, t'] \times B(\mathcal{O}, R)$ and passing to the limit as $R \to +\infty$, one infers that

$$ (38) \qquad \int_{\mathbb{R}^n} U_i(t', \mathbf{x})\mathbf{d}\mathbf{x} - \int_{\mathbb{R}^n} U_i(t, \mathbf{x})\mathbf{d}\mathbf{x} = \int_t^{t'} \left( \int_{\mathbb{R}^n} r_i(\mathbf{x})U_i(s, \mathbf{x})\mathbf{d}\mathbf{x} \right) \mathbf{d}s \\ - \int_t^{t'} \left( \int_{\mathbb{R}^n} u_i(s, \mathbf{x})\mathbf{d}\mathbf{x} \right) \left( \int_{\mathbb{R}^n} U_i(s, \mathbf{x})\mathbf{d}\mathbf{x} \right) \mathbf{d}s \\ + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} \int_t^{t'} \int_{\mathbb{R}^n} (U_k(s, \mathbf{x}) - U_i(s, \mathbf{x}))\mathbf{d}\mathbf{x}\mathbf{d}s. $$

By (37) and the definitions of $r_i$ and $E$, the integrals $\int_{\mathbb{R}^n} r_i(\mathbf{x})U_i(s, \mathbf{x})\mathbf{d}\mathbf{x}$, $\int_{\mathbb{R}^n} u_i(s, \mathbf{x})\mathbf{d}\mathbf{x}$ and $\int_{\mathbb{R}^n} U_i(s, \mathbf{x})\mathbf{d}\mathbf{x}$ are all bounded in absolute value, say by a constant $C$, uniformly in $s \in [0, T_1]$ and $1 \le i \le H$. Thus,

$$ (39) \qquad \left| \int_{\mathbb{R}^n} U_i(t', \mathbf{x})\mathbf{d}\mathbf{x} - \int_{\mathbb{R}^n} U_i(t, \mathbf{x})\mathbf{d}\mathbf{x} \right| \le (C + C^2 + 2\delta C) |t - t'| $$

for every $1 \le i \le H$ and $0 < t < t' \le T_1$ (and also with $t = 0$ by continuity). As a consequence, the maps $t \mapsto \int_{\mathbb{R}^n} U_i(t, \mathbf{x})\mathbf{d}\mathbf{x}$ are of class $C^{0,1}([0, T_1]) \subset C^{0,1/2}([0, T_1])$. Finally, $\mathbf{U} \in \mathbf{E}$.

Next, we show that $\mathcal{F}$ is a contraction mapping, that is, there exists $0 < \kappa < 1$ such that

$$ (40) \qquad \forall (\mathbf{u}, \mathbf{v}) \in E \times E, \quad d_E(\mathcal{F}(\mathbf{u}), \mathcal{F}(\mathbf{v})) \le \kappa d_E(\mathbf{u}, \mathbf{v}). $$

Take any $\mathbf{u}, \mathbf{v} \in E$, and denote $\mathbf{U} := \mathcal{F}(\mathbf{u})$, $\mathbf{V} := \mathcal{F}(\mathbf{v})$, $\mathbf{w} := \mathbf{u} - \mathbf{v}$ and $\mathbf{W} := \mathbf{U} - \mathbf{V}$. With these notations, for each $1 \le i \le H$, the function $W_i$ belongs to $C([0, T_1] \times \mathbb{R}^n, \mathbb{R}^H) \cap C_{t;\mathbf{x}}^{1;2}((0, T_1] \times \mathbb{R}^n, \mathbb{R}^H)$ and solves

$$ (41) \qquad \partial_t W_i(t, \mathbf{x}) = \frac{\mu^2}{2} \Delta W_i(t, \mathbf{x}) + \left( r_i(\mathbf{x}) - \int_{\mathbb{R}^n} u_i(t, \mathbf{y}) \mathbf{d}\mathbf{y} \right) W_i(t, \mathbf{x}) \\ - V_i(t, \mathbf{x}) \int_{\mathbb{R}^n} w_i(t, \mathbf{y}) \mathbf{d}\mathbf{y} + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} (W_k(t, \mathbf{x}) - W_i(t, \mathbf{x})) $$
in $(0, T_1] \times \mathbb{R}^n$, together with $W_i(0, \cdot) = 0$ in $\mathbb{R}^n$. Consider any $\xi \in \mathbb{R}^n$ with $\|\xi\| = 1$. Our supersolution candidate for this system will be $\overline{W} = (\overline{W}_1, \dots, \overline{W}_H)$ defined by

$$ \overline{W}_i(t, \mathbf{x}) := \frac{1}{3} d_E(\mathbf{u}, \mathbf{v}) e^{t/T_1 - \theta \xi \cdot \mathbf{x}}, \quad 1 \le i \le H, \ 0 < t \le T_1, \ \mathbf{x} \in \mathbb{R}^n. $$

Notice that the functions $\overline{W}_i$ actually do not depend on $i$. First of all, we recall that $\omega_n$ is the $(n-1)$-dimensional Lebesgue measure of the unit Euclidean sphere of $\mathbb{R}^n$. Now, since $\mathbf{V} = \mathcal{F}(\mathbf{v}) \in E$, one has, for every $1 \le i \le H$ and $(t, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n$,

$$ \begin{aligned} \left| V_i(t, \mathbf{x}) \int_{\mathbb{R}^n} w_i(t, \mathbf{y}) \mathrm{d}\mathbf{y} \right| & \le V_i(t, \mathbf{x}) \int_{\mathbb{R}^n} |u_i(t, \mathbf{y}) - v_i(t, \mathbf{y})| \mathrm{d}\mathbf{y} \\ & \le d_E(\mathbf{u}, \mathbf{v}) V_i(t, \mathbf{x}) \int_{\mathbb{R}^n} e^{-\theta \|\mathbf{y}\|} \mathrm{d}\mathbf{y} \\ & = \theta^{-n} \omega_n (n-1)! d_E(\mathbf{u}, \mathbf{v}) V_i(t, \mathbf{x}) \\ & \le \theta^{-n} \omega_n (n-1)! d_E(\mathbf{u}, \mathbf{v}) K e^{(r_{\max}+1)t - \theta \|\mathbf{x}\|} \\ & \le \theta^{-n} \omega_n (n-1)! d_E(\mathbf{u}, \mathbf{v}) K e^{(r_{\max}^+ + 1)T_1 - \theta \xi \cdot \mathbf{x}}. \end{aligned} \tag{42} $$

On the other hand, the definition (30) of $T_1$ yields

$$ \frac{1}{T_1} \ge \frac{\theta^2 \mu^2}{2} + r_{\max}^+ + 3 \theta^{-n} \omega_n (n-1)! K e^{(r_{\max}^+ + 1)T_1 - t/T_1} \tag{43} $$

for all $t \in [0, T_1]$. Since $r_i \le r_{\max} \le r_{\max}^+$ in $\mathbb{R}^n$ and $u_i \ge 0$ in $[0, T_1] \times \mathbb{R}^n$, it then follows from (42)-(43) that each nonnegative function $\overline{W}_i$ satisfies

$$ \begin{aligned} \partial_t \overline{W}_i(t, \mathbf{x}) & \ge \frac{\mu^2}{2} \Delta \overline{W}_i(t, \mathbf{x}) + r_{\max}^+ \overline{W}_i(t, \mathbf{x}) - V_i(t, \mathbf{x}) \int_{\mathbb{R}^n} w_i(t, \mathbf{y}) \mathrm{d}\mathbf{y} \\ & \ge \frac{\mu^2}{2} \Delta \overline{W}_i(t, \mathbf{x}) + \left( r_i(\mathbf{x}) - \int_{\mathbb{R}^n} u_i(t, \mathbf{y}) \mathrm{d}\mathbf{y} \right) \overline{W}_i(t, \mathbf{x}) \\ & \phantom{\ge} - V_i(t, \mathbf{x}) \int_{\mathbb{R}^n} w_i(t, \mathbf{y}) \mathrm{d}\mathbf{y} + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \ne i}}^H (\overline{W}_k(t, \mathbf{x}) - \overline{W}_i(t, \mathbf{x})) \end{aligned} $$

for all $(t, \mathbf{x}) \in (0, T_1] \times \mathbb{R}^n$ (recall that $\overline{W}_k - \overline{W}_i \equiv 0$). Therefore, the functions $Z_i := \overline{W}_i - W_i$ are continuous in $[0, T_1] \times \mathbb{R}^n$, nonnegative at $t = 0$, and satisfy

$$ \partial_t Z_i(t, \mathbf{x}) \ge \frac{\mu^2}{2} \Delta Z_i(t, \mathbf{x}) + \left( r_i(\mathbf{x}) - \int_{\mathbb{R}^n} u_i(t, \mathbf{y}) \mathrm{d}\mathbf{y} \right) Z_i(t, \mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \ne i}}^H (Z_k(t, \mathbf{x}) - Z_i(t, \mathbf{x})) $$

in $(0, T_1] \times \mathbb{R}^n$. The comparison principle for this linear cooperative system implies that $Z_i(t, \mathbf{x}) \ge 0$, that is, $W_i(t, \mathbf{x}) \le \overline{W}_i(t, \mathbf{x})$, for every $1 \le i \le H$ and $(t, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n$. By changing $W_i$ into $-W_i$ in (41) (and the sign of the first integral in the second line) and arguing similarly, one gets that $|W_i(t, \mathbf{x})| \le \overline{W}_i(t, \mathbf{x}) = (d_E(\mathbf{u}, \mathbf{v})/3) \times e^{t/T_1 - \theta \xi \cdot \mathbf{x}}$ for every $1 \le i \le H$ and $(t, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n$. Since $\xi$ was arbitrary in the unit Euclidean sphere of $\mathbb{R}^n$, one concludes that

$$ |U_i(t, \mathbf{x}) - V_i(t, \mathbf{x})| = |W_i(t, \mathbf{x})| \le \frac{1}{3} d_E(\mathbf{u}, \mathbf{v}) e^{t/T_1 - \theta \|\mathbf{x}\|} \le \frac{e}{3} d_E(\mathbf{u}, \mathbf{v}) e^{-\theta \|\mathbf{x}\|} \tag{44} $$

for every $1 \le i \le H$ and $(t, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n$.
In order to establish (40), it remains to estimate the $C^{0,1/2}([0, T_1])$ norm of the functions $t \mapsto \int_{\mathbb{R}^n}(U_i(t, \mathbf{x}) - V_i(t, \mathbf{x})) d\mathbf{x}$. For every $1 \le i \le H$ and $0 < t < t' \le T_1$, by using (38) for both $U_i$ and $V_i$, one infers that

$$ (45) \qquad \int_{\mathbb{R}^n} W_i(t', \mathbf{x}) d\mathbf{x} - \int_{\mathbb{R}^n} W_i(t, \mathbf{x}) d\mathbf{x} = \int_t^{t'} \left( \int_{\mathbb{R}^n} r_i(\mathbf{x}) W_i(s, \mathbf{x}) d\mathbf{x} \right) ds \\ - \int_t^{t'} \left( \int_{\mathbb{R}^n} u_i(s, \mathbf{x}) d\mathbf{x} \right) \left( \int_{\mathbb{R}^n} W_i(s, \mathbf{x}) d\mathbf{x} \right) ds \\ - \int_t^{t'} \left( \int_{\mathbb{R}^n} w_i(s, \mathbf{x}) d\mathbf{x} \right) \left( \int_{\mathbb{R}^n} V_i(s, \mathbf{x}) d\mathbf{x} \right) ds \\ + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} \int_t^{t'} \int_{\mathbb{R}^n} (W_k(s, \mathbf{x}) - W_i(s, \mathbf{x})) d\mathbf{x} ds. $$

Let us now bound in absolute value each term of the right-hand side of (45). Since $|r_i(\mathbf{x})| \le |r_{\max}| + \alpha \|\mathbf{x}\|^2 + \alpha B^2$ by (2) (we recall that $B = \max_{1 \le i \le H} \|\mathbf{O}_i\|$), it follows from (44) that

$$ \left| \int_t^{t'} \left( \int_{\mathbb{R}^n} r_i(\mathbf{x}) W_i(s, \mathbf{x}) d\mathbf{x} \right) ds \right| \le \frac{e}{3} d_E(\mathbf{u}, \mathbf{v}) \left[ \left( |r_{\max}| + \alpha B^2 \right) \theta^{-n} \omega_n(n-1)! + \alpha \theta^{-n-2} \omega_n(n+1)! \right] \times |t-t'| $$

for all $0 < t < t' \le T_1$. Similarly, using (44) and $\mathbf{u} \in E$, together with $T_1 \le 1$, the second term of the right-hand side of (45) can be bounded by

$$ \left| \int_t^{t'} \left( \int_{\mathbb{R}^n} u_i(s, \mathbf{x}) d\mathbf{x} \right) \left( \int_{\mathbb{R}^n} W_i(s, \mathbf{x}) d\mathbf{x} \right) ds \right| \le \frac{1}{3} d_E(\mathbf{u}, \mathbf{v}) K e^{r_{\max}^+ + 2} \theta^{-2n} \omega_n^2(n-1)!^2 \times |t-t'|, $$

for all $0 < t < t' \le T_1$. Furthermore, since $|w_i(s, \mathbf{x})| \le d_E(\mathbf{u}, \mathbf{v}) e^{-\theta \|\mathbf{x}\|}$ and $|V_i(s, \mathbf{x})| \le K e^{(r_{\max}+1)s - \theta \|\mathbf{x}\|} \le K e^{r_{\max}^+ + 1 - \theta \|\mathbf{x}\|}$ for every $1 \le i \le H$ and $(s, \mathbf{x}) \in [0, T_1] \times \mathbb{R}^n$ (because $\mathbf{V} = \mathcal{F}(\mathbf{v}) \in E$ and $T_1 \le 1$), the integral in the third line of the right-hand side of (45) can be estimated as

$$ \left| \int_t^{t'} \left( \int_{\mathbb{R}^n} w_i(s, \mathbf{x}) d\mathbf{x} \right) \left( \int_{\mathbb{R}^n} V_i(s, \mathbf{x}) d\mathbf{x} \right) ds \right| \le d_E(\mathbf{u}, \mathbf{v}) K e^{r_{\max}^+ + 1} \theta^{-2n} \omega_n^2(n-1)!^2 \times |t-t'|. $$

Lastly, (44) implies that

$$ \left| \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} \int_t^{t'} \int_{\mathbb{R}^n} (W_k(s, \mathbf{x}) - W_i(s, \mathbf{x})) d\mathbf{x} ds \right| \le \frac{2}{3} \delta e d_E(\mathbf{u}, \mathbf{v}) \theta^{-n} \omega_n(n-1)! \times |t-t'| $$

for every $1 \le i \le H$ and $0 < t < t' \le T_1$. Finally, putting all the above inequalities into (45) and remembering the definition (29) of $C_1$, it follows that

$$ \left| \int_{\mathbb{R}^n} W_i(t', \mathbf{x}) d\mathbf{x} - \int_{\mathbb{R}^n} W_i(t, \mathbf{x}) d\mathbf{x} \right| \le C_1 d_E(\mathbf{u}, \mathbf{v}) \times |t-t'| $$

for all $0 < t < t' \le T_1$, and then for all $0 \le t \le t' \le T_1$ by continuity of $s \mapsto \int_{\mathbb{R}^n} W_i(s, \mathbf{x}) d\mathbf{x}$ in $[0, T_1]$. Remember now that $C_1T_1 \le C_1\sqrt{T_1} \le 1/3$ by (30), and that $W_i(0, \cdot) = 0$ in $\mathbb{R}^N$ for every $1 \le i \le H$. As a consequence,

$$ \left\| t \mapsto \int_{\mathbb{R}^n} W_i(t, \mathbf{x}) d\mathbf{x} \right\|_{C^{0,1/2}([0,T_1])} \le \frac{2}{3} d_E(\mathbf{u}, \mathbf{v}) $$
for all $1 \le i \le H$. Together with (44), and using the definition (32) of $d_E$, one obtains
$d_E(\mathbf{U}, \mathbf{V}) \le \kappa d_E(\mathbf{u}, \mathbf{v})$, with $\kappa := \max(e/3, 2/3) = e/3 < 1$, that is, (40) is proved.

The Banach-Picard fixed point theorem then implies the existence and uniqueness of $\mathbf{u} \in E$ such that $\mathcal{F}(\mathbf{u}) = \mathbf{u}$. Therefore, $\mathbf{u}$ belongs to $C([0, T_1] \times \mathbb{R}^n, \mathbb{R}^H) \cap C_{t;\mathbf{x}}^{1;2}((0, T_1] \times \mathbb{R}^n, \mathbb{R}^H)$ and it solves (1) in $(0, T_1] \times \mathbb{R}^n$ with initial condition $\mathbf{u}^0$. Furthermore, owing to the definition of $E$ in (31) and the inequality $\theta \le 1/\mu$, $\mathbf{u}$ satisfies (6) in $[0, T_1] \times \mathbb{R}^n$ and, by (39), the functions $t \mapsto \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, \mathbf{d}\mathbf{x}$ are Lipschitz-continuous in $[0, T_1]$.

Let us now show that $\mathbf{u}$ can be extended to a global solution, that is, defined for all $t \ge 0$. First of all, observe that $T_1$ in (30) can be written as

$$T_1 = (A_1 + A_2 K + A_3 K^2)^{-1},$$

for some positive constants $A_1, A_2$ and $A_3$ that are independent of $K$. Consider then any $T'_1 \in (0, T_1)$. From (6) satisfied by $\mathbf{u}$ at time $T'_1$, one has $0 \le u_i(T'_1, \mathbf{x}) \le K e^{(r_{\max}+1)T'_1} e^{-\theta \|\mathbf{x}\|}$ for all $1 \le i \le H$ and $\mathbf{x} \in \mathbb{R}^n$. In other words, the continuous functions $u_i(T'_1, \cdot)$ satisfy properties similar to (4)-(5), with $K$ replaced by $K e^{(r_{\max}+1)T'_1}$. From the arguments of the previous paragraphs, there exists a unique solution $\tilde{\mathbf{u}} = (\tilde{u}_1, \dots, \tilde{u}_H) \in C([T'_1, T'_2] \times \mathbb{R}^n, \mathbb{R}^H) \cap C_{t;\mathbf{x}}^{1;2}((T'_1, T'_2] \times \mathbb{R}^n, \mathbb{R}^H)$ such that, for every $1 \le i \le H$, $\tilde{u}_i(T'_1, \cdot) = u_i(T'_1, \cdot)$ in $\mathbb{R}^n$,

$$0 \le \tilde{u}_i(t, x) \le K e^{(r_{\max}+1)T'_1} e^{(r_{\max}+1)(t-T'_1)-\theta \|\mathbf{x}\|} = K e^{(r_{\max}+1)t-\theta \|\mathbf{x}\|}$$

for all $(t, \mathbf{x}) \in [T'_1, T'_2] \times \mathbb{R}^n$ and the functions $t \mapsto \int_{\mathbb{R}^n} \tilde{u}_i(t, \mathbf{x}) \, \mathbf{d}\mathbf{x}$ are in $C^{0,1}([T'_1, T'_2])$, with

$$T'_2 - T'_1 = (A_1 + A_2 K e^{(r_{\max}+1)T'_1} + A_3 K^2 e^{2(r_{\max}+1)T'_1})^{-1}.$$

Since the uniqueness also holds in any interval $[T'_1, T''_2]$ with $T'_1 < T''_2 \le T'_2$ and since $\mathbf{u}$ restricted to $[T'_1, \min(T_1, T'_2)] \times \mathbb{R}^n$ satisfies the same properties as the restriction of $\tilde{\mathbf{u}}$ to this set, one gets that these two restrictions are equal. Therefore, the function $\mathbf{u} = (u_1, \dots, u_H)$ extended by $\tilde{\mathbf{u}}$ in $(T_1, T'_2] \times \mathbb{R}^n$ (provided this set is not empty), is then of class $C([0, T'_2] \times \mathbb{R}^n, \mathbb{R}^H) \cap C_{t;\mathbf{x}}^{1;2}((0, T'_2] \times \mathbb{R}^n, \mathbb{R}^H)$. Furthermore, $0 \le u_i(t, x) \le K e^{(r_{\max}+1)t-\theta \|\mathbf{x}\|}$ for all $(t, \mathbf{x}) \in [0, T'_2] \times \mathbb{R}^n$ and $1 \le i \le H$, and the functions $t \mapsto \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, \mathbf{d}\mathbf{x}$ are in $C^{0,1}([0, T'_2])$. Since this argument holds for every $T'_1 \in (0, T_1)$, it follows that $\mathbf{u}$ can then be extended in $[0, T_2) \times \mathbb{R}^n$, with

$$T_2 - T_1 = (A_1 + A_2 K e^{(r_{\max}+1)T_1} + A_3 K^2 e^{2(r_{\max}+1)T_1})^{-1},$$

and the above properties hold with $[0, T'_2]$ and $(0, T'_2]$ replaced by $[0, T_2)$ and $(0, T_2)$, respectively, with the functions $t \mapsto \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, \mathbf{d}\mathbf{x}$ belonging to $C_{\text{loc}}^{0,1}([0, T_2])$.

By an immediate induction, there exists a sequence $(T_m)_{m \ge 1}$ of positive real numbers such that, for every $m \ge 1$,

$$(46) \qquad T_{m+1} - T_m = (A_1 + A_2 K e^{(r_{\max}+1)T_m} + A_3 K^2 e^{2(r_{\max}+1)T_m})^{-1}$$

and $\mathbf{u}$ can be extended as a $C([0, T_m] \times \mathbb{R}^n, \mathbb{R}^H) \cap C_{t;\mathbf{x}}^{1;2}((0, T_m] \times \mathbb{R}^n, \mathbb{R}^H)$ function satisfying $0 \le u_i(t, x) \le K e^{(r_{\max}+1)t-\theta \|\mathbf{x}\|}$ for all $(t, \mathbf{x}) \in [0, T_m] \times \mathbb{R}^n$ and $1 \le i \le H$, and the functions $t \mapsto \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, \mathbf{d}\mathbf{x}$ are in $C_{\text{loc}}^{0,1}([0, T_m])$. The sequence $(T_m)_{m \ge 1}$ is increasing and is not bounded (otherwise, it would converge to a positive real number while $T_{m+1} - T_m$ would converge to 0, which would contradict (46)). Therefore, $T_m \to +\infty$ and the function $\mathbf{u}$ is then defined in $[0, +\infty) \times \mathbb{R}^n$, and it satisfies all the desired properties of the conclusion of Theorem 2.1.
Finally, if $\mathbf{v}$ is another solution satisfying all the properties of $\mathbf{u}$ listed in Theorem 2.1, it follows that the restrictions of $\mathbf{v}$ and $\mathbf{u}$ to $[0, T_1] \times \mathbb{R}^n$ belong to $E$ defined in (31), and that $\mathcal{F}(\mathbf{v}) = \mathbf{v}$, from the equations satisfied by the functions $v_i$. Therefore, by uniqueness, the restrictions of $\mathbf{v}$ and $\mathbf{u}$ to $[0, T_1] \times \mathbb{R}^n$ are equal. By an immediate induction, using the construction of the above sequence $(T_m)_{m \ge 1}$, one gets that $\mathbf{v}$ and $\mathbf{u}$ are equal in $[0, T_m) \times \mathbb{R}^n$ for every $m \ge 1$, hence $\mathbf{u} = \mathbf{v}$ in $[0, +\infty) \times \mathbb{R}^n$. The proof of Theorem 2.1 is thereby complete. $\square$

**Remark 6.1.** It follows from the above proof that, for any $\varepsilon \in (0,1)$, the quantity $r_{\max} + 1$ could be replaced by $r_{\max} + \varepsilon$ in (6), at the expense of replacing $\min(\theta, 1/\mu)$ by $\min(\theta, \sqrt{\varepsilon}/\mu)$. Indeed, in the above proof, by assuming that $0 < \theta \le \sqrt{\varepsilon}/\mu$ without loss of generality, the functions $\bar{U}_i$ in (34) would now be defined as $\bar{U}_i(t, \mathbf{x}) := K e^{(r_{\max}+\varepsilon)t-\theta\xi \cdot \mathbf{x}}$ and would satisfy (35), with $r_{\max} + \varepsilon$ instead of $r_{\max} + 1$, but they would still fulfill (36). The rest of the proof is unchanged, with just $(r_{\max} + 1, r_{\max}^+ + 1, r_{\max}^+ + 2)$ replaced by $(r_{\max} + \varepsilon, r_{\max}^+ + \varepsilon, r_{\max}^+ + \varepsilon + 1)$ in the definitions or calculations. As a consequence, if $r_{\max} < 0$, the solution $\mathbf{u}$ constructed in Theorem 2.1 is such that, for every $1 \le i \le H$, $\|u_i(t, \cdot)\|_{L^{\infty}(\mathbb{R}^n)} \to 0$ as $t \to +\infty$ exponentially fast.

*Proof of Proposition 2.2.* (i) We first assume $\lambda_H > 0$. Consider the following linear system, for $1 \le i \le H$:

$$ (47) \quad \partial_t v_i(t, \mathbf{x}) = \frac{\mu^2}{2} \Delta v_i(t, \mathbf{x}) + r_i(\mathbf{x}) v_i(t, \mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \ne i}}^{H} (v_k(t, \mathbf{x}) - v_i(t, \mathbf{x})). $$

It corresponds to a Malthusian growth of the population. Since $N_i(t) = \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, d\mathbf{x} \ge 0$ for all $t \ge 0$ by (6), the solution $\mathbf{u}$ of (1)-(5) is a subsolution of (47), while the function $(t, \mathbf{x}) \mapsto \Phi(\mathbf{x}) e^{-\lambda_H t}$ is a solution, by definition of the principal eigenpair $(\lambda_H, \Phi)$ of (8). Since $\mathbf{u}^0$ is assumed to be compactly supported and $\Phi$ is continuous and positive componentwise, there is a constant $A > 0$ such that $\mathbf{u}^0 \le A \Phi$ componentwise in $\mathbb{R}^n$. The comparison principle applied to the cooperative system (47) implies that $\mathbf{u}(t, \mathbf{x}) \le A \Phi(\mathbf{x}) e^{-\lambda_H t}$ componentwise for all $(t, \mathbf{x}) \in [0, +\infty) \times \mathbb{R}^n$, hence $N_i(t) \to 0$ as $t \to +\infty$ for every $1 \le i \le H$.

(ii) We now assume that $\lambda_H = 0$. Notice that the previous arguments imply that the continuous functions $N_i$ are bounded in $[0, +\infty)$. More precisely, as in the previous paragraph, since $\mathbf{u}^0$ is assumed to be compactly supported, there is a constant $A > 0$ such that $\mathbf{u}^0 \le A \Phi$ componentwise in $\mathbb{R}^n$, hence $\mathbf{u}(t, \mathbf{x}) \le A \Phi(\mathbf{x})$ componentwise for all $(t, \mathbf{x}) \in [0, +\infty) \times \mathbb{R}^n$. Assume now by way of contradiction that $\liminf_{t \to +\infty} \min_{1 \le i \le H} N_i(t) > 0$. Then there are $T > 0$ and $\varepsilon > 0$ such that $N_i(t) = \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, d\mathbf{x} \ge \varepsilon$ for every $1 \le i \le H$ and $t \ge T$. For $t \ge T$, the solution $\mathbf{u}$ is a subsolution of the linear system

$$ \partial_t v_i(t, \mathbf{x}) = \frac{\mu^2}{2} \Delta v_i(t, \mathbf{x}) + (r_i(\mathbf{x}) - \varepsilon) v_i(t, \mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \ne i}}^{H} (v_k(t, \mathbf{x}) - v_i(t, \mathbf{x})), $$

whereas $(t, \mathbf{x}) \mapsto A \Phi(\mathbf{x}) e^{-\varepsilon(t-T)}$ is a solution (since $\lambda_H = 0$). The maximum principle applied to this cooperative system implies that $\mathbf{u}(t, \mathbf{x}) \le A \Phi(\mathbf{x}) e^{-\varepsilon(t-T)}$ for all $(t, \mathbf{x}) \in [T, +\infty) \times \mathbb{R}^n$. Hence, $N_i(t) \to 0$ as $t \to +\infty$, for every $1 \le i \le H$, contradicting our assumption. As a consequence, $\liminf_{t \to +\infty} \min_{1 \le i \le H} N_i(t) = 0$.

(iii) Assume finally that $\lambda_H < 0$ and assume that the desired conclusion (12) is not satisfied. Then there are $T > 0$ and $\varepsilon \in (0, -\lambda_H)$ such that $N_i(t) = \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, d\mathbf{x} \le$
$-\lambda_H - \varepsilon$ for every $1 \le i \le H$ and $t \ge T$. Since $\lambda_H = \lim_{R \to +\infty} \lambda_H^R$, where $\lambda_H^R$ denotes the principal eigenvalue of the operator $\mathcal{A}$ acting on $(C_0^\infty(\overline{B(\mathcal{O}, R)))^H$, there is $R > 0$ such that $\lambda_H^R \le \lambda_H + \varepsilon/2$, hence $N_i(t) = \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, d\mathbf{x} \le -\lambda_H^R - \varepsilon/2$ for every $1 \le i \le H$ and $t \ge T$. Now, each function $u_i(T, \cdot)$ is continuous and positive in $\mathbb{R}^n$, hence there is $\eta > 0$ such that $\mathbf{u}(T, \mathbf{x}) \ge \eta \Phi^R(\mathbf{x})$ componentwise for all $\mathbf{x} \in \overline{B(\mathcal{O}, R)$, where $\Phi^R \in (C_0^\infty(\overline{B(\mathcal{O}, R)))^H$ denotes a principal eigenvector of $\mathcal{A}$, associated to the principal eigenvalue $\lambda_H^R$ (that is, $\Phi^R$ is positive componentwise in $B(\mathcal{O}, R)$, it vanishes on $\partial B(\mathcal{O}, R)$ and it satisfies $\mathcal{A}\Phi^R = \lambda_H^R\Phi^R$ in $\overline{B(\mathcal{O}, R)$). The function $(t, \mathbf{x}) \mapsto \eta \Phi^R(\mathbf{x}) e^{(\varepsilon/2)(t-T)}$ is a solution of the linear system

$$ \partial_t v_i(t, \mathbf{x}) = \frac{\mu^2}{2} \Delta v_i(t, \mathbf{x}) + (r_i(\mathbf{x}) + \lambda_H^R + \varepsilon/2) v_i(t, \mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} (v_k(t, \mathbf{x}) - v_i(t, \mathbf{x})) $$

in $\mathbb{R} \times \overline{B(\mathcal{O}, R)$, with Dirichlet boundary conditions on $\mathbb{R} \times \partial B(\mathcal{O}, R)$, whereas $\mathbf{u}$ is a supersolution in $[T, +\infty) \times \overline{B(\mathcal{O}, R)$, and the functions are compared at time $T$. The maximum principle applied to this cooperative system implies that $\mathbf{u}(t, \mathbf{x}) \ge \eta \Phi^R(\mathbf{x}) e^{(\varepsilon/2)(t-T)}$ for all $(t, \mathbf{x}) \in [T, +\infty) \times \overline{B(\mathcal{O}, R)$. Hence, $N_i(t) \ge \int_{B(\mathcal{O}, R)} u_i(t, \mathbf{x}) \, d\mathbf{x} \to +\infty$ as $t \to +\infty$, for every $1 \le i \le H$, contradicting our assumption. As a consequence, (12) is proved, and the proof of Proposition 2.2 is complete. $\square$

**Remark 6.2.** We here show that, when $\lambda_H \ge 0$, the system (1) does not admit any positive bounded stationary state. Assume by way of contradiction that such a stationary state $\mathbf{p} = (p_1, \dots, p_H)$ exists. It necessarily converges to 0 exponentially as $\|\mathbf{x}\| \to +\infty$, from the confining properties of the fitnesses $r_i$. Let $R > 0$ be such that $r_i(\mathbf{x}) < 0$ for all $1 \le i \le H$ and $\|\mathbf{x}\| \ge R$. Let also $A > 0$ be such that $\mathbf{p}(\mathbf{x}) \le A \Phi(\mathbf{x})$ componentwise for all $\|\mathbf{x}\| \le R$, where $\Phi = (\varphi_1, \dots, \varphi_H)$ denotes as usual the, positive, principal eigenvector of (8). Notice that, from its positivity, the stationary state $\mathbf{p}$ of (1) is a subsolution of the following linear system

$$ \frac{\mu^2}{2} \Delta p_i(\mathbf{x}) + r_i(\mathbf{x}) p_i(\mathbf{x}) + \frac{\delta}{H-1} \sum_{\substack{k=1 \\ k \neq i}}^{H} (p_k(\mathbf{x}) - p_i(\mathbf{x})) \ge 0, $$

while $\Phi$ is a supersolution of the same linear system, because $\lambda_H \ge 0$. The weak maximum principle for this system then yields $\mathbf{p}(\mathbf{x}) \le A \Phi(\mathbf{x})$ componentwise for all $\|\mathbf{x}\| \ge R$. Hence, $\mathbf{p}(\mathbf{x}) \le A \Phi(\mathbf{x})$ componentwise for all $\mathbf{x} \in \mathbb{R}^n$. Let finally $\mathbf{u}$ be the solution of (1) given by Theorem 2.1, with initial condition $\mathbf{u}^0 := \mathbf{p}$. On the one hand, by the uniqueness property of that theorem, there holds $\mathbf{u}(t, \mathbf{x}) = \mathbf{p}(\mathbf{x})$ for all $t \ge 0$ and $\mathbf{x} \in \mathbb{R}^n$. On the other hand, the same arguments as in the proof of Proposition 2.2 in the cases $\lambda_H > 0$ or $\lambda_H = 0$ then imply that, in both cases, one has at least $\liminf_{t \to +\infty} (\min_{1 \le i \le H} \int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, d\mathbf{x}) = 0$. This leads to a contradiction, since $\int_{\mathbb{R}^n} u_i(t, \mathbf{x}) \, d\mathbf{x} = \int_{\mathbb{R}^n} p_i(\mathbf{x}) \, d\mathbf{x} > 0$ is independent of $t \ge 0$, for every $1 \le i \le H$. As a conclusion, (1) does not admit any positive bounded stationary state when $\lambda_H \ge 0$.

**Proof of Proposition 2.3.** (i) Assume that $H = 1$ and $\lambda_1 < 0$. In this case, any positive bounded stationary state $p_1$ necessarily decays to 0 at least exponentially as $\|\mathbf{x}\| \to +\infty$, and then belongs to $H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$ from standard elliptic estimates. It is also a principal eigenfunction of the operator $\mathcal{A}$ defined in (7), with associated principal eigenvalue $-\int_{\mathbb{R}^n} p_1$. By uniqueness of the eigenpair $(\lambda_1, \varphi_1)$, up to a multiplication of $\varphi_1$ by positive constants, we deduce that there exists $\kappa > 0$ such that $p_1 = \kappa \varphi_1$ in $\mathbb{R}^n$, and that
$\int_{\mathbb{R}^n} p_1 = -\lambda_1$. We get that $\kappa = -\lambda_1 / (\int_{\mathbb{R}^n} \varphi_1)$. As it is immediate to check that $\kappa \varphi_1$ is a positive bounded stationary state of (1), this shows (i).

(ii) Assume that $H = 2$ and $\lambda_2 < 0$. Without loss of generality, up to rotation and translation of the phenotypic space, we can fix the optima $\mathbf{O}_1$ and $\mathbf{O}_2$ at $\mathbf{O}_1 := (-\beta, 0, \dots, 0)$, $\mathbf{O}_2 := (\beta, 0, \dots, 0)$, for some $\beta \ge 0$. First, we recall from [12] that, by uniqueness of the principal eigenfunctions (up to multiplication) and by symmetry of the problem (1) with respect to the origin, we have $\varphi_2(\mathbf{x}) = \varphi_1(\iota(\mathbf{x})) > 0$ for all $\mathbf{x} \in \mathbb{R}^n$, with

$$ \iota(\mathbf{x}) := \iota(x_1, x_2, \dots, x_n) = (-x_1, x_2, \dots, x_n). $$

It follows that $\int_{\mathbb{R}^n} \varphi_1 = \int_{\mathbb{R}^n} \varphi_2 > 0$. Furthermore, the principal eigenfunctions satisfy:

$$ (48) \qquad -\frac{\mu^2}{2} \Delta\varphi_i(\mathbf{x}) - r_i(\mathbf{x})\varphi_i(\mathbf{x}) + \delta(\varphi_i(\mathbf{x}) - \varphi_j(\mathbf{x})) = \lambda_2 \varphi_i, \quad \mathbf{x} \in \mathbb{R}^n, $$

for $(i, j) = (1, 2)$ and $(i, j) = (2, 1)$. Define $(p_1, p_2)$ as in (13) and $\bar{r}_i$ as in (11), with $u_i(t, \mathbf{y})$ replaced by $p_i(\mathbf{y})$. Owing to the definition (13), integrating (48) over $\mathbb{R}^n$ and using that the $\varphi_i$'s and its first order derivatives decay at least exponentially as $\|\mathbf{x}\| \to +\infty$, we infer that $\bar{r}_1 = \bar{r}_2 = -\lambda_2$. We also note that $\int_{\mathbb{R}^n} p_1 = \int_{\mathbb{R}^n} p_2 = -\lambda_2$. Dividing (48) by $\int_{\mathbb{R}^n} \varphi_i = \int_{\mathbb{R}^n} \varphi_j$, and multiplying by $-\lambda_2$ we then get:

$$ -\frac{\mu^2}{2} \Delta p_i(\mathbf{x}) - r_i(\mathbf{x})p_i(\mathbf{x}) + \delta(p_i(\mathbf{x}) - p_j(\mathbf{x})) = -\left(\int_{\mathbb{R}^n} p_i(\mathbf{y}) \, d\mathbf{y}\right) p_i(\mathbf{x}), \quad \mathbf{x} \in \mathbb{R}^n, $$

for $(i, j) = (1, 2)$ and $(i, j) = (2, 1)$. This proves (ii).

(iii) Assume that $H \ge 3$ and consider the optima $\mathbf{O}_j$'s located at $\mathbf{O}_1 := (-\beta, 0, \dots, 0)$ and $\mathbf{O}_i := (\beta, 0, \dots, 0)$ for all $2 \le i \le H$, with $\beta > 0$. Let $(\lambda_H, \Phi)$ be the principal eigenpair of (8), with $\Phi = (\varphi_1, \dots, \varphi_H)$ here normalised so that $\int_{\mathbb{R}^n} \varphi_1 = 1$. First of all, by uniqueness, we readily check that $\varphi_i = \varphi_2$ for all $2 \le i \le H$, hence the eigenfunctions satisfy:

$$ (49) \qquad \begin{cases} -\frac{\mu^2}{2} \Delta\varphi_1(\mathbf{x}) - r_1(\mathbf{x})\varphi_1(\mathbf{x}) + \delta(\varphi_1(\mathbf{x}) - \varphi_2(\mathbf{x})) = \lambda_H \varphi_1(\mathbf{x}), \\ -\frac{\mu^2}{2} \Delta\varphi_2(\mathbf{x}) - r_2(\mathbf{x})\varphi_2(\mathbf{x}) + \frac{\delta}{H-1}(\varphi_2(\mathbf{x}) - \varphi_1(\mathbf{x})) = \lambda_H \varphi_2(\mathbf{x}). \end{cases} $$

We now claim that, for all $\beta > 0$ large enough, no positive multiple of $\Phi$ is a stationary state of the system (1). Assume by way of contradiction that there exist a sequence $(\beta_m)_{m \in \mathbb{N}}$ diverging to $+\infty$ and a sequence $(\kappa_m)_{m \in \mathbb{N}}$ of positive real numbers such that $(\kappa_m \varphi_{1,m}, \kappa_m \varphi_{2,m}, \dots, \kappa_m \varphi_{2,m})$ is a positive bounded stationary state of (1), associated with the optima $\mathbf{O}_{1,m} := (-\beta_m, 0, \dots, 0)$ and $\mathbf{O}_{2,m} = \dots = \mathbf{O}_{H,m} := (\beta_m, 0, \dots, 0)$. Call $(\lambda_{H,m})_{m \in \mathbb{N}}$ the sequence of associated principal eigenvalues, and normalise the principal eigenvectors so that $\int_{\mathbb{R}^n} \varphi_{1,m}(\mathbf{x})d\mathbf{x} = 1$. For each $m \in \mathbb{N}$, plugging the solution $(\kappa_m \varphi_{1,m}, \kappa_m \varphi_{2,m}, \dots, \kappa_m \varphi_{2,m})$ into (1) yields $\kappa_m \int_{\mathbb{R}^n} \varphi_{i,m}(\mathbf{x})d\mathbf{x} = -\lambda_H$ for every $1 \le i \le H$, hence $\int_{\mathbb{R}^n} \varphi_{2,m} = \int_{\mathbb{R}^n} \varphi_{1,m} = 1$ in (49) (with $\varphi_1$ and $\varphi_2$ replaced by $\varphi_{1,m}$ and $\varphi_{2,m}$). From the conclusion of Proposition 3.2 (whose proof is independent of this one), each $\lambda_{H,m}$ satisfies

$$ \lambda_1 \le \lambda_{H,m} \le \lambda_1 + \delta = \frac{\mu n \sqrt{\alpha}}{2} - r_{\max} + \delta. $$
Therefore, up to extraction of a subsequence one can assume that

$$ \lambda_{H,m} \to \lambda_{\infty} \in \mathbb{R} \text{ as } m \to +\infty. $$

Call now $\phi_m := \varphi_{1,m}(\cdot + \mathbf{O}_{1,m})$ and $\psi_m := \varphi_{2,m}(\cdot + \mathbf{O}_{1,m}) = \varphi_{2,m}(\cdot - \mathbf{O}_{2,m})$. From (49), one has

$$ (50) \quad -\frac{\mu^2}{2} \Delta \phi_m(\mathbf{x}) - \left( r_{\max} - \alpha \frac{\|\mathbf{x}\|^2}{2} \right) \phi_m(\mathbf{x}) + \delta(\phi_m(\mathbf{x}) - \psi_m(\mathbf{x})) = \lambda_{H,m} \phi_m(\mathbf{x}) $$

for all $m \in \mathbb{N}$ and $\mathbf{x} \in \mathbb{R}^n$. Since $\int_{\mathbb{R}^n} \phi_m = \int_{\mathbb{R}^n} \psi_m = 1$ for all $m \in \mathbb{N}$, it follows that the sequence $(\Delta\phi_m)_{m \in \mathbb{N}}$ is bounded in $L_{\text{loc}}^1(\mathbb{R}^n)$ and then (see [22]<fnref n="3" />) that the sequence $(\phi_m)_{m \in \mathbb{N}}$ is bounded in $W_{\text{loc}}^{1,p}(\mathbb{R}^n)$ for every $1 \le p < n/(n-1)$ if $n \ge 2$ (respectively bounded in $W_{\text{loc}}^{1,\infty}(\mathbb{R})$ if $n=1$). Therefore, up to extraction of a subsequence, there is a nonnegative function $\phi_\infty$ in $W_{\text{loc}}^{1,p}(\mathbb{R}^n)$ for every $1 \le p < n/(n-1)$ if $n \ge 2$ (respectively for every $1 \le p \le +\infty$ if $n=1$) such that $\phi_m \to \phi_\infty$ almost everywhere in $\mathbb{R}^n$ and in $L_{\text{loc}}^q(\mathbb{R}^n)$ for every $1 \le q < n/(n-2)$ if $n \ge 3$ (respectively for every $1 \le q < +\infty$ if $n=2$, respectively in $L_{\text{loc}}^\infty(\mathbb{R})$ if $n=1$), and $\partial_{x_j}\phi_m \to \partial_{x_j}\phi_\infty$ weakly in $L_{\text{loc}}^p(\mathbb{R}^n)$ for every $1 \le p < n/(n-1)$ if $n \ge 2$ (respectively for every $1 \le p < +\infty$ if $n=1$). Since $\int_{B(\mathcal{O},R)} \phi_\infty = \lim_{m \to +\infty} \int_{B(\mathcal{O},R)} \phi_m \le 1$ for every $R > 0$, it follows from the monotone convergence theorem that $\int_{\mathbb{R}^n} \phi_\infty \le 1$. On the other hand, by integrating (50) over $\mathbb{R}^n$ and using that $\int_{\mathbb{R}^n} \phi_m = \int_{\mathbb{R}^n} \psi_m = 1$, one gets that

$$ (51) \quad \frac{\alpha R^2}{2} \int_{\mathbb{R}^n \setminus B(\mathcal{O},R)} \phi_m(\mathbf{x}) \, d\mathbf{x} \le \frac{\alpha}{2} \int_{\mathbb{R}^n} \|\mathbf{x}\|^2 \phi_m(\mathbf{x}) \, d\mathbf{x} = \lambda_{H,m} + r_{\max} \le \frac{\mu n \sqrt{\alpha}}{2} + \delta $$

for all $m \in \mathbb{N}$ and $R > 0$. Since the nonnegative functions $\phi_m$ converge to $\phi_\infty$ in $L_{\text{loc}}^1(\mathbb{R}^n)$ and have unit $L^1(\mathbb{R}^n)$ norm, it follows that $\int_{\mathbb{R}^n} \phi_\infty(\mathbf{x}) \, d\mathbf{x} = 1$. As in (51), by using now $\psi_m(\cdot + 2\mathbf{O}_{2,m}) = \varphi_{2,m}(\cdot + \mathbf{O}_{2,m})$ instead of $\phi_m = \varphi_{1,m}(\cdot + \mathbf{O}_{1,m})$, one also has

$$ \frac{\alpha R^2}{2} \int_{\mathbb{R}^n \setminus B(2\mathbf{O}_{2,m}, R)} \psi_m(\mathbf{x}) \, d\mathbf{x} \le \frac{\alpha}{2} \int_{\mathbb{R}^n} \|\mathbf{x}\|^2 \psi_m(\mathbf{x} + 2\mathbf{O}_{2,m}) \, d\mathbf{x} = \lambda_{H,m} + r_{\max} \le \frac{\mu n \sqrt{\alpha}}{2} + \delta $$

for all $m \in \mathbb{N}$ and $R > 0$. Since $\psi_m \ge 0$ in $\mathbb{R}^n$ and since $\mathbf{O}_{2,m} = (\beta_m, 0, \dots, 0)$ with $\lim_{m \to +\infty} \beta_m = +\infty$, it follows that $\psi_m \to 0$ as $m \to +\infty$ in $L_{\text{loc}}^1(\mathbb{R}^n)$. Passing to the limit as $m \to +\infty$ in (50) after multiplying against any $C^\infty(\mathbb{R}^n)$ function with compact support, it follows that $\phi_\infty$ is a distributional solution of

$$ (52) \quad -\frac{\mu^2}{2} \Delta \phi_\infty - \left( r_{\max} - \alpha \frac{\|\mathbf{x}\|^2}{2} \right) \phi_\infty(\mathbf{x}) + \delta \phi_\infty = \lambda_\infty \phi_\infty $$

in $\mathbb{R}^n$. Since $\phi_\infty \in W_{\text{loc}}^{1,p}(\mathbb{R}^n)$ for every $1 \le p < n/(n-1)$ if $n \ge 2$ (respectively for every $1 \le p \le +\infty$ if $n=1$), one infers from the weak formulation of (52) that $\partial_{x_j}\phi_\infty$ actually belongs to $W_{\text{loc}}^{1,p}(\mathbb{R}^n)$ for every $1 \le j \le n$ and then $\phi_\infty \in W_{\text{loc}}^{2,p}(\mathbb{R}^n)$, for every $1 \le p < n/(n-1)$ if $n \ge 2$ (respectively for every $1 \le p \le +\infty$ if $n=1$). Therefore, by bootstrapping and differentiating, it follows that $\phi_\infty$ is a $C^\infty(\mathbb{R}^n)$ solution of (52). Since it is nonnegative with unit $L^1(\mathbb{R}^n)$ mass, the strong elliptic maximum principle implies that $\phi_\infty > 0$ in $\mathbb{R}^n$. In particular, there is $R_0 > 0$ such that the function $\phi_\infty$ is subharmonic in $\mathbb{R}^n \setminus B(\mathcal{O}, R_0)$. Hence,

$$ 0 < \phi_{\infty}(\mathbf{x}) \le n \omega_n^{-1} \int_{B(\mathbf{x}, 1)} \phi_{\infty}(\mathbf{y}) \, d\mathbf{y} \le n \omega_n^{-1} \int_{\mathbb{R}^n \setminus B(\mathcal{O}, \|\mathbf{x}\| - 1)} \phi_{\infty}(\mathbf{y}) \, d\mathbf{y} $$

<footnote n="3">³The authors are grateful to L. Dupaigne and F. Murat for helpful discussions on this point.</footnote>
for all $||x|| \ge R_0 + 1$, where $\omega_n/n$ is the Lebesgue measure of the Euclidean unit ball of $\mathbb{R}^n$. Since $\phi_\infty \in L^1(\mathbb{R}^n)$, one then deduces that $\phi_\infty(\mathbf{x}) \to 0$ as $||\mathbf{x}|| \to +\infty$. On the other hand, for every $\sigma > 0$, there is $R_1 > 0$ such that $\alpha||\mathbf{x}||^2/2 - r_{\max} + \delta \ge \lambda_\infty$ for all $||\mathbf{x}|| \ge R_1$ and such that the function $\mathbf{x} \mapsto e^{-\sigma||\mathbf{x}||}$ is a supersolution of (52) in $\mathbb{R}^n \setminus B(\mathcal{O}, R_1)$. Therefore, by choosing $A > 0$ such that $\max_{\partial B(\mathcal{O}, R_1)} \phi_\infty \le Ae^{-\sigma R_1}$, the weak maximum principle implies that $0 < \phi_\infty(\mathbf{x}) \le Ae^{-\sigma||\mathbf{x}||}$ for all $||\mathbf{x}|| \ge R_1$. As a consequence, $\phi_\infty$ converges to 0 as $||\mathbf{x}|| \to +\infty$ faster than any exponentially decaying function. So does then the function $\mathbf{x} \mapsto ||\mathbf{x}||^2\phi_\infty(\mathbf{x})$, hence standard elliptic estimates imply that $||\nabla\phi_\infty||$ also converges to 0 as $||\mathbf{x}|| \to +\infty$ faster than any exponentially decaying function. Finally, $\phi_\infty \in H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$ and, since it is positive, it is therefore the principal eigenfunction of the operator $-(\mu^2/2)\Delta - (r_{\max} - \alpha||\mathbf{x}||^2/2) + \delta$, with principal eigenvalue $\lambda_\infty$. Hence,

$$ (53) \qquad \lambda_{\infty} = \lambda_{1} + \delta = \frac{\mu n \sqrt{\alpha}}{2} - r_{\max} + \delta. $$

Finally, by arguing similarly with the functions $\tilde{\phi}_m := \varphi_{1,m}(\cdot + \mathbf{O}_{2,m})$ and $\tilde{\psi}_m := \varphi_{2,m}(\cdot + \mathbf{O}_{2,m})$, and passing to the limit as $m \to +\infty$ in the equation

$$ -\frac{\mu^2}{2} \Delta \tilde{\psi}_m(\mathbf{x}) - \left( r_{\max} - \alpha \frac{||\mathbf{x}||^2}{2} \right) \tilde{\psi}_m(\mathbf{x}) + \frac{\delta}{H-1} (\tilde{\psi}_m(\mathbf{x}) - \tilde{\phi}_m(\mathbf{x})) = \lambda_{H,m} \tilde{\psi}_m(\mathbf{x}), $$

one similarly gets that

$$ (54) \qquad \lambda_{\infty} = \lambda_{1} + \frac{\delta}{H-1} = \frac{\mu n \sqrt{\alpha}}{2} - r_{\max} + \frac{\delta}{H-1}, $$

contradicting (53), since $\delta > 0$ and $H \ge 3$. The proof of Proposition 2.3 is thereby complete. $\square$

## 6.2 Proofs of Propositions 3.1-3.2 and 3.4

We actually start with the proof of Proposition 3.2, since its conclusion is used in the proof of Proposition 3.1.

*Proof of Proposition 3.2.* First of all, since

$$ (55) \qquad 2 \int_{\mathbb{R}^n} \phi(\mathbf{x}) \psi(\mathbf{x}) \, d\mathbf{x} \le \int_{\mathbb{R}^n} (\phi(\mathbf{x}))^2 \, d\mathbf{x} + \int_{\mathbb{R}^n} (\psi(\mathbf{x}))^2 \, d\mathbf{x} $$

for any two $L^2(\mathbb{R}^n)$ functions $\phi$ and $\psi$, and since the principal eigenvector $\Phi = (\varphi_1, \dots, \varphi_H)$ of $\mathcal{A}$ is normalised in $L^2(\mathbb{R}^n)^H$, the first inequality in (18) follows immediately. The second inequality follows from $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) = Q_H(\Phi)$ and the fact that $\lambda_1$ is the minimum of $Q_H$ when $\delta = 0$:

$$ \sum_{i=1}^{H} \left( \frac{\mu^2}{2} \int_{\mathbb{R}^n} ||\nabla \varphi_i(\mathbf{x})||^2 d\mathbf{x} - \int_{\mathbb{R}^n} r_i(\mathbf{x}) (\varphi_i(\mathbf{x}))^2 d\mathbf{x} \right) \ge \sum_{i=1}^{H} \lambda_1(\alpha, \mu) \int_{\mathbb{R}^n} (\varphi_i(\mathbf{x}))^2 d\mathbf{x} = \lambda_1(\alpha, \mu). $$

Now, recalling that $G_i = G(\cdot - \mathbf{O}_i)$ and testing the Rayleigh quotient $Q_H$ at $\Psi := (G_1/\sqrt{H}, \dots, G_H/\sqrt{H})$, we get that

$$ \begin{aligned} \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) &\le Q_H(\Psi) \\ &= \lambda_1(\alpha, \mu) + \delta \left( 1 - \sum_{1 \le i < j \le H} \frac{2}{H(H-1)} \int_{\mathbb{R}^n} G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x} \right) \\ &\le \lambda_1(\alpha, \mu) + \delta. \end{aligned} $$
Lastly, for any $\Gamma \in H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$ such that $\int_{\mathbb{R}^n} \Gamma^2 = 1/H$, we observe that the quantity $Q_H(\Gamma, \dots, \Gamma)$ is independent of $\delta$ and therefore the Rayleigh formula (9)-(10) (together with the lower bound $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \ge \mu n \sqrt{\alpha}/2 - r_{\max}$) implies that $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is bounded independently of $\delta$ (for every given $\alpha > 0, \mu > 0$ and $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$). $\square$

*Proof of Proposition 3.1.* Let us begin with the proof of the continuity of the map

$$ (\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \mapsto \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) $$

in $(0, +\infty)^3 \times (\mathbb{R}^n)^H$. Fix $(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \in (0, +\infty)^3 \times (\mathbb{R}^n)^H$ and consider any sequence $(\delta^m, \alpha^m, \mu^m, \mathbf{O}_1^m, \dots, \mathbf{O}_H^m)_{m \in \mathbb{N}}$ in $(0, +\infty)^3 \times (\mathbb{R}^n)^H$ and converging to $(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$. Since the sequence $(\lambda_H(\delta^m, \alpha^m, \mu^m, \mathbf{O}_1^m, \dots, \mathbf{O}_H^m))_{m \in \mathbb{N}}$ is bounded by Proposition 3.2, one can assume that it converges to some real number $\lambda$, up to extraction of a subsequence. For each $m \in \mathbb{N}$, call $\Phi^m := (\varphi_1^m, \dots, \varphi_H^m)$ the normalised principal eigenvector associated to the principal eigenvalue $\lambda_H(\delta^m, \alpha^m, \mu^m, \mathbf{O}_1^m, \dots, \mathbf{O}_H^m)$. As at the beginning of the proof of Proposition 3.2, the inequality (55) applied to $\int_{\mathbb{R}^n} \varphi_i^m \varphi_j^m$ implies that

$$ \sum_{i=1}^{H} \left( \frac{(\mu^m)^2}{2} \int_{\mathbb{R}^n} \|\nabla \varphi_i^m(\mathbf{x})\|^2 d\mathbf{x} + \frac{\alpha^m}{2} \int_{\mathbb{R}^n} \|\mathbf{x} - \mathbf{O}_i^m\|^2 (\varphi_i^m(\mathbf{x}))^2(\mathbf{x}) d\mathbf{x} \right) \\ \le Q_H(\Phi^m) + r_{\max} = \lambda_H(\delta^m, \alpha^m, \mu^m, \mathbf{O}_1^m, \dots, \mathbf{O}_H^m) + r_{\max} $$

for all $m \in \mathbb{N}$. Together with $\int_{\mathbb{R}^n} \|\Phi^m(\mathbf{x})\|^2 d\mathbf{x} = 1$, the sequence $(\Phi^m)_{m \in \mathbb{N}} = ((\varphi_1^m, \dots, \varphi_H^m))_{m \in \mathbb{N}}$ is bounded in $(H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$. Hence, up to extraction of a sub-sequence, it converges weakly in $H^1(\mathbb{R}^n)^H$, strongly in $L_{\text{loc}}^2(\mathbb{R}^n)^H$, and almost everywhere in $\mathbb{R}^n$, to a certain limit $\Phi^\infty := (\varphi_1^\infty, \dots, \varphi_H^\infty) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$, which is nonnegative componentwise. The boundedness of the sequences $(\int_{\mathbb{R}^n} \|\mathbf{x}\|^2 (\varphi_i^m(\mathbf{x}))^2 d\mathbf{x})_{m \in \mathbb{N}}$ also implies (as for the function $\phi_\infty$ after (51)) that the sequence $(\Phi^m)_{m \in \mathbb{N}}$ converges to $\Phi^\infty$ in $L^2(\mathbb{R}^n)^H$, and that

$$ \int_{\mathbb{R}^n} \|\Phi^\infty(\mathbf{x})\|^2 d\mathbf{x} = 1. $$

By passing to the limit weakly in (8) (with parameters $(\delta^m, \alpha^m, \mu^m, \mathbf{O}_1^m, \dots, \mathbf{O}_H^m)$), it follows that the pair $(\lambda, \Phi^\infty)$ is a weak solution of (8), with $\Phi^\infty$ being normalised in $L^2(\mathbb{R}^n)^H$ and nonnegative componentwise. From elliptic regularity, $\Phi^\infty$ is actually a $C^\infty(\mathbb{R}^n)^H$ solution of this problem and, from the strong elliptic maximum principle, it is positive componentwise in $\mathbb{R}^n$. By uniqueness of the principal eigenpair of (8), one infers that $\Phi^\infty = \Phi$ and $\lambda = \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$. By uniqueness of the limit, it also follows that $\Phi^m \to \Phi$ in $L^2(\mathbb{R}^n)^H$ as $m \to +\infty$. One has then shown the continuity of the map $(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \mapsto \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ in $(0, +\infty)^3 \times (\mathbb{R}^n)^H$. Together with the inequalities (18) and the definition (16), one concludes that the map $(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H) \mapsto \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is continuous in $[0, +\infty) \times (0, +\infty)^2 \times (\mathbb{R}^n)^H$.

Next, since the quantity $Q_H(\Psi)$ in (10) is affine with respect to $(\delta, \alpha, \mu) \in (0, +\infty)^3$ for each $\Psi \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$, one gets from (9) that the map $(\delta, \alpha, \mu) \mapsto \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is concave in $(0, +\infty)^3$, for each $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$. Together with the continuity proved in the previous paragraph, one infers that the map $(\delta, \alpha, \mu) \mapsto \lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ is concave in $[0, +\infty) \times (0, +\infty)^2$.

Consider now $\delta \ge 0, \alpha > 0, 0 < \mu_1 < \mu_2$, $(\mathbf{O}_i)_{1 \le i \le H} \in (\mathbb{R}^n)^H$, and let us show that $\lambda_H(\delta, \alpha, \mu_1, \mathbf{O}_1, \dots, \mathbf{O}_H) < \lambda_H(\delta, \alpha, \mu_2, \mathbf{O}_1, \dots, \mathbf{O}_H)$. The case $\delta = 0$ is obvious by (16),
so let us assume that $\delta > 0$. Let $\Phi_1$ and $\Phi_2$ be the principal eigenvectors associated to the
principal eigenvalues $\lambda_H(\delta, \alpha, \mu_1, \mathbf{O}_1, \dots, \mathbf{O}_H)$ and $\lambda_H(\delta, \alpha, \mu_2, \mathbf{O}_1, \dots, \mathbf{O}_H)$, respectively,
and let $Q_{H,1}$ and $Q_{H,2}$ be the functionals defined in (10), associated with $\mu_1$ and $\mu_2$,
respectively. Since no component of the vector $(\varphi_1, \dots, \varphi_H) := \Phi_2$ is constant, one has
$\int_{\mathbb{R}^n} \|\nabla \varphi_i(\mathbf{x})\|^2 d\mathbf{x} > 0$ for each $1 \le i \le H$, hence

$$
\lambda_H(\delta, \alpha, \mu_2, \mathbf{O}_1, \dots, \mathbf{O}_H) = Q_{H,2}(\Phi_2) > Q_{H,1}(\Phi_2) \ge \lambda_H(\delta, \alpha, \mu_1, \mathbf{O}_1, \dots, \mathbf{O}_H).
$$

The strict increasing monotonicity of $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ with respect to $\alpha > 0$
can be shown similarly. Lastly, the monotonicity of $\lambda_H(\delta, \alpha, \mu, \mathbf{O}_1, \dots, \mathbf{O}_H)$ with respect
to $\delta \ge 0$ will actually be shown independently in Proposition 4.3 in Section 6.3. $\square$

Proof of Proposition 3.4. The first part leading to (20) is actually a straightforward adaptation of the proof of [12, Proposition 4 (i)] performed in the case $H = 2$. Details are omitted. Next, since $R_H(\mathbf{x}) = r_{\max} - \alpha \|\mathbf{x} - \mathbf{M}\|^2/2 + R_H(\mathbf{M}) - r_{\max}$, we get (21). $\square$

6.3 Proofs of Propositions 4.3-4.8

*Proof of Proposition* 4.3. The parameters $\alpha > 0$ and $\mu > 0$ are fixed throughout the proof, as are the optima $(\mathbf{O}_i)_{1\le i\le H} \in (\mathbb{R}^n)^H$. Let us first show that the map $\delta \mapsto \lambda_H(\delta)$ is differentiable in $(0, +\infty)$. For $\delta > 0$ and $\Psi = (\psi_1, \dots, \psi_H) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$, we denote by $Q_H(\Psi, \delta)$ the Rayleigh quotient (10). We also denote $\Phi^\delta = (\varphi_1^\delta, \dots, \varphi_H^\delta)$ the principal (normalised) eigenvector associated with $\lambda_H(\delta)$, so that $\lambda_H(\delta) = Q_H(\Phi^\delta, \delta) = \min_{\Psi \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H} Q_H(\Psi, \delta)$.

Fix any $\delta > 0$. We notice that, for $\varepsilon > 0$,

$$
(56) \quad \frac{Q_H(\Phi^{\delta+\varepsilon}, \delta+\varepsilon) - Q_H(\Phi^{\delta+\varepsilon}, \delta)}{\varepsilon} \le \frac{\lambda_H(\delta+\varepsilon) - \lambda_H(\delta)}{\varepsilon} \le \frac{Q_H(\Phi^{\delta}, \delta+\varepsilon) - Q_H(\Phi^{\delta}, \delta)}{\varepsilon},
$$

together with

$$
(57) \quad \frac{Q_H(\Phi^{\delta+\epsilon}, \delta+\epsilon) - Q_H(\Phi^{\delta+\epsilon}, \delta)}{\epsilon} = 1 - \sum_{1\le i<j\le H} \frac{2}{H-1} \int_{\mathbb{R}^n} \varphi_i^{\delta+\epsilon}(\mathbf{x}) \varphi_j^{\delta+\epsilon}(\mathbf{x}) \, d\mathbf{x}
$$

and

$$
(58) \quad \frac{Q_H(\Phi^\delta, \delta + \epsilon) - Q_H(\Phi^\delta, \delta)}{\epsilon} = 1 - \sum_{1 \le i < j \le H} \frac{2}{H-1} \int_{\mathbb{R}^n} \varphi_i^\delta(\mathbf{x}) \varphi_j^\delta(\mathbf{x}) \, d\mathbf{x}.
$$

From the proof of the continuity property in Proposition 3.1, one knows that $\Phi^{\delta+\epsilon} \to \Phi^\delta$
in $L^2(\mathbb{R}^n)^H$ as $\epsilon \to 0$. Thus, inserting (57) and (58) into (56) and passing to the limit
as $\epsilon \to 0^+$, we get

$$
\frac{\lambda_H(\delta + \epsilon) - \lambda_H(\delta)}{\epsilon} \underset{\epsilon \to 0^+}{\longrightarrow} 1 - \sum_{1 \le i < j \le h} \frac{2}{H-1} \int_{\mathbb{R}^n} \varphi_i^\delta(\mathbf{x}) \varphi_j^\delta(\mathbf{x}) \, d\mathbf{x}.
$$

By reverting the inequalities in (56) when $\epsilon \in (-\delta, 0)$, and using the arbitrariness of $\delta > 0$,
one then concludes that the map $\delta \mapsto \lambda_H(\delta)$ is differentiable in $(0, +\infty)$, with

$$
(59) \quad \lambda'_H(\delta) = 1 - \sum_{1 \le i < j \le H} \frac{2}{H-1} \int_{\mathbb{R}^n} \varphi_i^\delta(\mathbf{x}) \varphi_j^\delta(\mathbf{x}) \, d\mathbf{x} < 1 \quad \text{for all } \delta > 0.
$$
Notice also that, if the optima $\mathbf{O}_i$ are all identical, then by uniqueness the functions $\varphi_i^\delta$ are also all identical, hence $\int_{\mathbb{R}^n} \varphi_i^\delta(\mathbf{x}) \varphi_j^\delta(\mathbf{x}) \mathrm{d}\mathbf{x} = 1/H$ for all $1 \le i \le H$ (remember that $\Phi^\delta = (\varphi_1^\delta, \dots, \varphi_H^\delta)$ has unit $L^2(\mathbb{R}^n)^H$ norm) and $\lambda'_H(\delta) = 0$ for all $\delta > 0$. On the other hand, if the optima $\mathbf{O}_i$ are not all identical, then, for any $\delta > 0$, the functions $\varphi_i^\delta$ are not all identical: otherwise, each $\varphi_i^\delta \in H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$ would be a classical positive solution of $-(\mu^2/2)\Delta\varphi_i^\delta(\mathbf{x}) - r_i(\mathbf{x})\varphi_i^\delta(\mathbf{x}) = \lambda_H(\delta)\varphi_i^\delta(\mathbf{x})$ in $\mathbb{R}^n$ and then $r_i(\mathbf{x}) \equiv r_j(\mathbf{x})$ in $\mathbb{R}^n$ and $\mathbf{O}_i = \mathbf{O}_j$ for all $1 \le i, j \le H$, a contradiction. Thus, for every $\delta > 0$, there are $1 \le i_0 < j_0 \le H$ such that $\varphi_{i_0}^\delta \ne \varphi_{j_0}^\delta$ in $\mathbb{R}^n$, hence

$$ \int_{\mathbb{R}^n} \varphi_{i_0}^\delta(\mathbf{x}) \varphi_{j_0}^\delta(\mathbf{x}) \mathrm{d}\mathbf{x} < \frac{1}{2} \int_{\mathbb{R}^n} (\varphi_{i_0}^\delta(\mathbf{x}))^2 \mathrm{d}\mathbf{x} + \frac{1}{2} \int_{\mathbb{R}^n} (\varphi_{j_0}^\delta(\mathbf{x}))^2 \mathrm{d}\mathbf{x}, $$

while the large inequality always holds for every $1 \le i, j \le H$. Therefore, (59) and the normalisation of $\Phi^\delta = (\varphi_1^\delta, \dots, \varphi_H^\delta)$ imply that $\lambda'_H(\delta) > 0$ for all $\delta > 0$.

Let us now show the differentiability of $\delta \mapsto \lambda_H(\delta)$ at $\delta = 0$. We recall that

$$ \lambda_H(0) = \lambda_1 = \frac{\mu n \sqrt{\alpha}}{2} - r_{\max} = \frac{\mu^2}{2} \int_{\mathbb{R}^n} \|\nabla G_i(\mathbf{x})\|^2 \mathrm{d}\mathbf{x} - \int_{\mathbb{R}^n} r_i(\mathbf{x}) (G_i(\mathbf{x}))^2 \mathrm{d}\mathbf{x} $$

for every $1 \le i \le H$, with $G_i(\mathbf{x}) = G(\mathbf{x} - \mathbf{O}_i)$ and $G$ defined in (15). We also observe that Proposition 3.2 implies that the function $\delta \mapsto \lambda_H(\delta)$ is continuous at 0, and that $0 \le \liminf_{\delta \to 0^+} (\lambda_H(\delta) - \lambda_H(0))/\delta \le \limsup_{\delta \to 0^+} (\lambda_H(\delta) - \lambda_H(0))/\delta < 1$ (the last strict inequality comes from the third inequality in (18)).

Consider now any sequence $(\varepsilon_m)_{m \in \mathbb{N}}$ of positive real numbers converging to 0. The same arguments as in the proof of the continuity property in Proposition 3.1 imply that, up to extraction of a subsequence, the sequence $(\Phi^{\varepsilon_m})_{m \in \mathbb{N}} = ((\varphi_1^{\varepsilon_m}, \dots, \varphi_H^{\varepsilon_m}))_{m \in \mathbb{N}}$ converges weakly in $H^1(\mathbb{R}^n)^H$, strongly in $L^2(\mathbb{R}^n)^H$, and almost everywhere in $\mathbb{R}^n$, to a certain limit $\Phi := (\varphi_1, \dots, \varphi_H) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$, which is nonnegative componentwise and satisfies $\int_{\mathbb{R}^n} \|\Phi(\mathbf{x})\|^2 \mathrm{d}\mathbf{x} = 1$. Furthermore, each function $\varphi_i$ is a $C^\infty(\mathbb{R}^n) \cap H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$ nonnegative solution of

$$ -\frac{\mu^2}{2} \Delta \varphi_i(\mathbf{x}) - r_i(\mathbf{x}) \varphi_i(\mathbf{x}) = \lambda_1 \varphi_i(\mathbf{x}) \quad \text{in } \mathbb{R}^n. $$

Therefore, by uniqueness of the principal eigenfunction for this equation, it follows that there is $p_i \in [0, +\infty)$ such that

$$ \varphi_i(\mathbf{x}) = p_i G_i(\mathbf{x}) = p_i G(\mathbf{x} - \mathbf{O}_i) \quad \text{for all } \mathbf{x} \in \mathbb{R}^n, $$

hence

$$ (60) \qquad (\varphi_1^{\varepsilon_m}, \dots, \varphi_H^{\varepsilon_m}) \to (p_1 G_1, \dots, p_H G_H) \quad \text{in } L^2(\mathbb{R}^n)^H \quad \text{as } m \to +\infty. $$

Since $\Phi$ is normalised in $L^2(\mathbb{R}^n)^H$, one also has

$$ p_1^2 + \dots + p_H^2 = 1. $$

Now, consider any $\mathbf{q} = (q_1, \dots, q_H) \in \mathbb{R}^H$ with $\|\mathbf{q}\| = 1$, and call $\Psi := (q_1 G_1, \dots, q_H G_H) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^H$. Notice that $\int_{\mathbb{R}^n} \|\Psi(\mathbf{x})\|^2 \mathrm{d}\mathbf{x} = 1$. Observe also that each $\varphi_{H,i}^{\varepsilon_m}$ belongs to $H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$, hence

$$ \lambda_1 \int_{\mathbb{R}^n} (\varphi_i^{\varepsilon_m}(\mathbf{x}))^2 \mathrm{d}\mathbf{x} \le \frac{\mu^2}{2} \int_{\mathbb{R}^n} \|\nabla \varphi_i^{\varepsilon_m}(\mathbf{x})\|^2 \mathrm{d}\mathbf{x} - r_{\max} \int_{\mathbb{R}^n} (\varphi_i^{\varepsilon_m}(\mathbf{x}))^2 \mathrm{d}\mathbf{x} + \frac{\alpha}{2} \int_{\mathbb{R}^n} \|\mathbf{x} - \mathbf{O}_i\|^2 (\varphi_i^{\varepsilon_m}(\mathbf{x}))^2(\mathbf{x}) \mathrm{d}\mathbf{x} \\ \underbrace{=}_{:=Q_1(\varphi_i^{\varepsilon_m}, 0)} Q_1(\varphi_i^{\varepsilon_m}, 0) $$
and $\lambda_H(0) = \lambda_1 \le Q_H(\Phi^{\epsilon_m}, 0)$ by summing the inequalities from $i=1$ to $H$. The inequality $\lambda_H(\epsilon_m) = Q_H(\Phi^{\epsilon_m}, \epsilon_m) \le Q_H(\Psi, \epsilon_m)$ and the equality $Q_H(\Psi, 0) = \lambda_1 = \lambda_H(0)$ then imply that

$$ \frac{Q_H(\Phi^{\epsilon_m}, \epsilon_m) - Q_H(\Phi^{\epsilon_m}, 0)}{\epsilon_m} \le \frac{\lambda_H(\epsilon_m) - \lambda_H(0)}{\epsilon_m} \le \frac{Q_H(\Psi, \epsilon_m) - Q_H(\Psi, 0)}{\epsilon_m} $$

for every $m \in \mathbb{N}$, hence, as in (57)-(58),

$$ (61) \qquad 1 - \sum_{1 \le i < j \le H} \frac{2}{H-1} \int_{\mathbb{R}^n} \varphi_i^{\epsilon_m}(\mathbf{x}) \varphi_j^{\epsilon_m}(\mathbf{x}) \, d\mathbf{x} \\ \le \frac{\lambda_H(\epsilon_m) - \lambda_H(0)}{\epsilon_m} \le 1 - \sum_{1 \le i < j \le H} \frac{2}{H-1} \int_{\mathbb{R}^n} q_i q_j G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x}. $$

Using (60) and passing to the limit as $m \to +\infty$ in the left-hand side of the previous inequality leads to

$$ \sum_{1 \le i < j \le H} \int_{\mathbb{R}^n} q_i q_j G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x} \le \sum_{1 \le i < j \le H} \int_{\mathbb{R}^n} p_i p_j G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x}, $$

hence

$$ \sum_{1 \le i, j \le H} \int_{\mathbb{R}^n} q_i q_j G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x} \le \sum_{1 \le i, j \le H} \int_{\mathbb{R}^n} p_i p_j G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x}, $$

since each $G_i = G(\cdot - \mathbf{O}_i)$ has unit $L^2(\mathbb{R}^n)$ norm and since both vectors $\mathbf{p} := (p_1, \dots, p_H)$ and $\mathbf{q} = (q_1, \dots, q_H)$ of $\mathbb{R}^H$ have unit Euclidean norm. Denote $A = (a_{ij})_{1 \le i, j \le H}$ the symmetric matrix whose entries are

$$ (62) \qquad a_{ij} := \int_{\mathbb{R}^n} G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x} = e^{-\sqrt{\alpha} \|\mathbf{O}_i - \mathbf{O}_j\|^2 / (4\mu)} > 0, $$

by (17). One has then obtained that $\mathbf{q}A\mathbf{q} \le \mathbf{p}A\mathbf{p}$ for all $\mathbf{q} \in \mathbb{R}^H$ with unit Euclidean norm (namely the quadratic form associated to $A$ is maximal on the unit Euclidean sphere of $\mathbb{R}^H$ at the normalised vector $\mathbf{p}$). From Perron-Frobenius theorem, since $a_{ij} > 0$ and $p_i \ge 0$ for all $1 \le i, j \le H$, it follows that $\mathbf{p}A\mathbf{p}$ is the largest eigenvalue $\mu_A$ of $A$, that $\mathbf{p}$ is the unique normalised principal eigenvector of the matrix $A$, and that $p_i > 0$ for all $1 \le i \le H$. By uniqueness of the limit, one concludes that $\Phi^{\epsilon} \to (p_1 G_1, \dots, p_H G_H)$ in $L^2(\mathbb{R}^n)^H$ as $\epsilon \to 0^+$. Furthermore, by using $\mathbf{q} = \mathbf{p}$ and $\epsilon$ instead of $\epsilon_m$ in (61), and by passing to the limit as $\epsilon \to 0^+$, one gets that $\delta \mapsto \lambda_H(\delta)$ is differentiable at 0, with

$$ (63) \qquad \lambda'_H(0) = 1 - \sum_{1 \le i < j \le H} \frac{2}{H-1} \int_{\mathbb{R}^n} p_i p_j G_i(\mathbf{x}) G_j(\mathbf{x}) \, d\mathbf{x} = 1 - \frac{1}{H-1} \sum_{1 \le i \ne j \le H} p_i p_j a_{ij}. $$

But $\sum_{1 \le i \ne j \le H} p_i p_j a_{ij} = \mathbf{p}A\mathbf{p} - \sum_{1 \le i \le H} p_i^2 = \mu_A - 1$. Finally,

$$ (64) \qquad \lambda'_H(0) = 1 - \frac{\mu_A - 1}{H - 1} = \frac{H - \mu_A}{H - 1}. $$

Observe that (63) confirms that $\lambda'_H(0) < 1$ (since $a_{ij} > 0$ and $p_i > 0$ for all $1 \le i, j \le H$). Furthermore, if the optima $\mathbf{O}_i$ are all identical, then we already know that $\lambda'_H(\delta) = 0$ for all $\delta > 0$, hence the function $\delta \mapsto \lambda_H(\delta)$, which is continuous in $[0, +\infty)$ is constant in $[0, +\infty)$, thus $\lambda'_H(0) = 0$ too and $\lambda_H(\delta) = \lambda_H(0) = \lambda_1 = \mu n \sqrt{\alpha}/2 - r_{\max}$ for all $\delta \ge 0$.
On the other hand, if the optima $O_i$ are not all identical, then the functions $G_i$ are not all identical and there are $1 \le i_0 < j_0 \le H$ such that $G_{i_0} \ne G_{j_0}$, hence $p_{i_0} G_{i_0} \ne p_{j_0} G_{j_0}$ (remember that $p_i > 0$ for all $1 \le i \le H$ and that the functions $G_i$ have all unit $L^2(\mathbb{R}^n)$ norm) and

$$ \int_{\mathbb{R}^n} p_{i_0} p_{j_0} G_{i_0}(\mathbf{x}) G_{j_0}(\mathbf{x}) \, d\mathbf{x} < \frac{1}{2} \int_{\mathbb{R}^n} (p_{i_0} G_{i_0}(\mathbf{x}))^2 \, d\mathbf{x} + \frac{1}{2} \int_{\mathbb{R}^n} (p_{j_0} G_{j_0}(\mathbf{x}))^2 \, d\mathbf{x} = \frac{p_{i_0}^2 + p_{j_0}^2}{2}, $$

while the large inequality holds for all $1 \le i, j \le H$. Therefore, (63) and the equality $p_1^2 + \cdots + p_H^2 = 1$ imply that $\lambda'_H(0) > 0$.

Lastly, with $\mathbf{O}_1, \mathbf{O}_2$ as in (19), formula (64) implies that

$$ \lambda'_2(0) = 1 - e^{-\sqrt{\alpha\beta^2}/\mu} $$

and

$$ \lambda'_3(0) = 1 - \sqrt{\frac{a_{12}^2 + a_{13}^2 + a_{23}^2}{3}} \times \cos \left[ \frac{1}{3} \arccos \left( a_{12}a_{13}a_{23} \sqrt{\frac{27}{(a_{12}^2 + a_{13}^2 + a_{23}^2)^3}} \right) \right], $$

from the del Ferro-Cardan formula for the roots of the characteristic polynomial of the matrix $A$ (notice that the argument inside the arccos ranges in $(0, 1]$ from the positivity of the $a_{ij}$'s and the arithmetic-geometric mean inequality). After some algebraic transformations based on (62), one infers that $\mathcal{P} := \{\mathbf{O}_3 = (a_1, \dots, a_n) \in \mathbb{R}^n : \lambda'_3(0) < \lambda'_2(0)\}$ is given by formula (25), and that $\{\mathbf{O}_3 = (a_1, \dots, a_n) \in \mathbb{R}^n : \lambda'_3(0) > \lambda'_2(0)\} = \mathbb{R}^n \setminus \bar{\mathcal{P}}$. The last conclusion of Proposition 4.3 immediately follows, and the proof is thereby complete. $\square$

**Proof of Proposition 4.4.** Denote $\lambda_3^e = \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$, with $\mathbf{O}_1 = (-\beta, 0, \dots, 0)$, $\mathbf{O}_2 = (\beta, 0, \dots, 0)$, $\mathbf{O}_3 = (0, \sqrt{3}\beta, 0, \dots, 0)$, and $\beta > 0$. We first show some symmetry properties. Let $(\varphi_1, \varphi_2, \varphi_3)$ be the normalised principal eigenvector associated with the principal eigenvalue $\lambda_3^e$. Define the symmetry

$$ (65) \qquad \iota(\mathbf{x}) = \iota(x_1, x_2, \dots, x_n) := (-x_1, x_2, \dots, x_n). $$

As $r_2 \circ \iota = r_1$ and $r_3 \circ \iota = r_3$, we get that $(\varphi_2 \circ \iota, \varphi_1 \circ \iota, \varphi_3 \circ \iota)$ is also a normalised positive eigenvector. By uniqueness, it follows that

$$ \varphi_2 = \varphi_1 \circ \iota \quad \text{and} \quad \varphi_3 = \varphi_3 \circ \iota. $$

Now, define

$$ \kappa(\mathbf{x}) = \kappa(x_1, x_2, \dots, x_n) := \left( \frac{1}{2}(x_1 - \beta) + \frac{\sqrt{3}}{2}x_2, \frac{\sqrt{3}}{2}(x_1 + \beta) - \frac{1}{2}x_2, x_3, \dots, x_n \right), $$

which is the orthogonal affine reflection with respect to the hyperplane parallel to the directions $x_3, \dots, x_n$ and containing the line joining $\mathbf{O}_1$ and $(\mathbf{O}_2 + \mathbf{O}_3)/2$. We have

$$ r_1 \circ \kappa = r_1, \quad r_2 \circ \kappa = r_3, \quad r_3 \circ \kappa = r_2. $$

Additionally, one has $\Delta(\varphi \circ \kappa)(\mathbf{x}) = (\Delta\varphi)(\kappa(\mathbf{x}))$ for all $\varphi \in C^2(\mathbb{R}^n)$ and $\mathbf{x} \in \mathbb{R}^n$. Thus, $(\varphi_1 \circ \kappa, \varphi_3 \circ \kappa, \varphi_2 \circ \kappa)$ is also a normalised positive eigenvector. Again, by uniqueness, one infers that

$$ \varphi_1 = \varphi_1 \circ \kappa, \quad \varphi_2 = \varphi_3 \circ \kappa, \quad \text{and} \quad \varphi_3 = \varphi_2 \circ \kappa. $$
Consider now the equation satisfied by $\varphi_1$, namely

$$ (66) \quad -\frac{\mu^2}{2} \Delta \varphi_1(\mathbf{x}) - r_1(\mathbf{x}) \varphi_1(\mathbf{x}) - \delta \left( \frac{\varphi_2(\mathbf{x}) + \varphi_3(\mathbf{x})}{2} - \varphi_1(\mathbf{x}) \right) = \lambda_3^e \varphi_1(\mathbf{x}). $$

Multiplying by $\varphi_1$ and integrating by parts, we get

$$ \int_{\mathbb{R}^n} \left( \frac{\mu^2}{2} \|\nabla \varphi_1(\mathbf{x})\|^2 - r_1(\mathbf{x}) (\varphi_1(\mathbf{x}))^2 \right) d\mathbf{x} + \delta \int_{\mathbb{R}^n} (\varphi_1(\mathbf{x}))^2 d\mathbf{x} \\ - \frac{\delta}{2} \int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_2(\mathbf{x}) d\mathbf{x} - \frac{\delta}{2} \int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_3(\mathbf{x}) d\mathbf{x} = \lambda_3^e \int_{\mathbb{R}^n} (\varphi_1(\mathbf{x}))^2 d\mathbf{x}. $$

Moreover,

$$ \int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_3(\mathbf{x}) d\mathbf{x} = \int_{\mathbb{R}^n} (\varphi_1 \circ \kappa)(\mathbf{x}) (\varphi_2 \circ \kappa)(\mathbf{x}) d\mathbf{x} = \int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_2(\mathbf{x}) d\mathbf{x}. $$

Finally, we have

$$ (67) \quad \int_{\mathbb{R}^n} \left( \frac{\mu^2}{2} \|\nabla \varphi_1(\mathbf{x})\|^2 - r_1(\mathbf{x}) (\varphi_1(\mathbf{x}))^2 \right) d\mathbf{x} + \delta \int_{\mathbb{R}^n} (\varphi_1(\mathbf{x}))^2 d\mathbf{x} \\ - \delta \int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_2(\mathbf{x}) d\mathbf{x} = \lambda_3^e \int_{\mathbb{R}^n} (\varphi_1(\mathbf{x}))^2 d\mathbf{x}, $$

and, by symmetry (using a composition with the function $\iota$),

$$ (68) \quad \int_{\mathbb{R}^n} \left( \frac{\mu^2}{2} \|\nabla \varphi_2(\mathbf{x})\|^2 - r_2(\mathbf{x}) (\varphi_2(\mathbf{x}))^2 \right) d\mathbf{x} + \delta \int_{\mathbb{R}^n} (\varphi_2(\mathbf{x}))^2 d\mathbf{x} \\ - \delta \int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_2(\mathbf{x}) d\mathbf{x} = \lambda_3^e \int_{\mathbb{R}^n} (\varphi_2(\mathbf{x}))^2 d\mathbf{x}. $$

Adding (67) and (68), we observe that

$$ (69) \quad Q_2 \left( \frac{\varphi_1}{\sqrt{\int_{\mathbb{R}^n} (\varphi_1^2 + \varphi_2^2)(\mathbf{x}) d\mathbf{x}}}, \frac{\varphi_2}{\sqrt{\int_{\mathbb{R}^n} (\varphi_1^2 + \varphi_2^2)(\mathbf{x}) d\mathbf{x}}} \right) = \lambda_3^e, $$

where $Q_2$ is the Rayleigh quotient defined by (10) with $H = 2$. Since $(\varphi_1, \varphi_2) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^2$, it follows from (9) that $\lambda_3^e \ge \lambda_2(\mathbf{O}_1, \mathbf{O}_2)$.

Finally, if $\lambda_3^e$ were equal to $\lambda_2(\mathbf{O}_1, \mathbf{O}_2)$, since the principal eigenvector associated with the principal eigenvalue $\lambda_2(\mathbf{O}_1, \mathbf{O}_2)$ with $H = 2$ is the unique positive normalised minimum of $Q_2$ in $(H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^2$, it would follow from (69) that $(\varphi_1, \varphi_2)$ would be (up to normalisation) a principal eigenvector of (8) with $H = 2$. Together with (66) and using the positivity of $\delta$, one would infer that $\varphi_3 \equiv \varphi_2$ in $\mathbb{R}^n$. Similarly, by using the equation similar to (66) satisfied by $\varphi_2$, one would get that $\varphi_3 \equiv \varphi_1$ in $\mathbb{R}^n$. Finally, $\varphi_1 \equiv \varphi_2 \equiv \varphi_3$ in $\mathbb{R}^n$ and (66) implies that $\varphi_1$ would be a positive multiple of the normalised principal eigenfunction $G_1 = G(\cdot - \mathbf{O}_1)$ of the case of only one host. But, similarly, $\varphi_2$ would be a positive multiple of $G_2 = G(\cdot - \mathbf{O}_2)$. This leads to a contradiction, since $\mathbf{O}_1 \ne \mathbf{O}_2$ (because $\beta > 0$). Therefore, $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3) = \lambda_3^e > \lambda_2(\mathbf{O}_1, \mathbf{O}_2)$ and the proof of Proposition 4.4 is thereby complete. $\square$

**Proof of Proposition 4.5.** With the notation (26), evaluating the Rayleigh quotient associated with $\lambda_3 = \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$ at $\tilde{\Phi} := (\tilde{\varphi}_1, \tilde{\varphi}_2, 0)$, we get

$$ (70) \quad \lambda_3 \le Q_3(\tilde{\Phi}) = \tilde{\lambda}_2 = \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_2). $$
On the other hand, as the minimum in the Rayleigh quotient $Q_3$ is reached at $\Phi = (\varphi_1, \varphi_2, \varphi_3)$, the normalised eigenvector associated with $\lambda_3$, we have

$$ \lambda_3 \geq \min_{(\psi_1, \psi_2) \in \mathcal{E}^2} \left\{ \frac{\mu^2}{2} \int_{\mathbb{R}^n} \|\nabla \psi_1\|^2 + \|\nabla \psi_2\|^2 - \int_{\mathbb{R}^n} (r_1 \psi_1^2 + r_2 \psi_2^2) - \delta \int_{\mathbb{R}^n} \psi_1 \psi_2 \right\} \\ + \min_{\psi \in \mathcal{E}} \left\{ \frac{\mu^2}{2} \int_{\mathbb{R}^n} \|\nabla \psi\|^2 - \int_{\mathbb{R}^n} r_3 \psi^2 \right\} + \delta - \delta \int_{\mathbb{R}^n} \varphi_3 (\varphi_1 + \varphi_2), $$

where $\mathcal{E} = H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$. As the first minimum in the above formula is precisely $(\tilde{\lambda}_2 - \delta)(1 - \int_{\mathbb{R}^n} \varphi_3^2)$, and the second minimum is $\lambda_1 \int_{\mathbb{R}^n} \varphi_3^2$, we get

$$ (71) \qquad \lambda_3 \geq (\tilde{\lambda}_2 - \delta) \left(1 - \int_{\mathbb{R}^n} \varphi_3^2\right) + \lambda_1 \int_{\mathbb{R}^n} \varphi_3^2 + \delta - \delta \int_{\mathbb{R}^n} \varphi_3 (\varphi_1 + \varphi_2). $$

Our goal is now to show that $\int_{\mathbb{R}^n} \varphi_3 (\varphi_1 + \varphi_2)$ becomes small as $\|\mathbf{O}_3\| \to +\infty$. Since $\lambda_3 = Q_3(\Phi)$, we have

$$ \lambda_3 > - \sum_{i=1}^{3} \int_{\mathbb{R}^n} r_i \varphi_i^2 + \delta - \delta \sum_{1 \le i < j \le 3} \int_{\mathbb{R}^n} \varphi_i \varphi_j \ge - \sum_{i=1}^{3} \int_{\mathbb{R}^n} r_i \varphi_i^2 $$

since $\int_{\mathbb{R}^n} \|\Phi\|^2 = 1$. Using (70) and the definition (2) of $r_i(\mathbf{x})$, we get

$$ \alpha \sum_{i=1}^{3} \int_{\mathbb{R}^n} \frac{\|\mathbf{x} - \mathbf{O}_i\|^2}{2} (\varphi_i(\mathbf{x}))^2 d\mathbf{x} < \tilde{\lambda}_2 + r_{\max}. $$

This implies that, for any radius $R > 0$, and $i = 1, 2, 3$,

$$ \|\varphi_i\|_{L^2(\mathbb{R}^n \setminus B(\mathbf{O}_i, R))}^2 < \frac{2(\tilde{\lambda}_2 + r_{\max})}{\alpha R^2}. $$

Next, we have, for $i = 1, 2$,

$$ \int_{\mathbb{R}^n} \varphi_3 \varphi_i = \int_{B(\mathbf{O}_i, R)} \varphi_3 \varphi_i + \int_{\mathbb{R}^n \setminus B(\mathbf{O}_i, R)} \varphi_3 \varphi_i \\ \le \|\varphi_3\|_{L^2(B(\mathbf{O}_i, R))} \|\varphi_i\|_{L^2(B(\mathbf{O}_i, R))} + \|\varphi_3\|_{L^2(\mathbb{R}^n \setminus B(\mathbf{O}_i, R))} \|\varphi_i\|_{L^2(\mathbb{R}^n \setminus B(\mathbf{O}_i, R))}. $$

Taking $R = \min(\|\mathbf{O}_3 - \mathbf{O}_1\|/2, \|\mathbf{O}_3 - \mathbf{O}_2\|/2)$, we have $B(\mathbf{O}_i, R) \subset \mathbb{R}^n \setminus B(\mathbf{O}_3, R)$ for $i = 1, 2$. Thus, we get

$$ \begin{align} \int_{\mathbb{R}^n} \varphi_3 (\varphi_1 + \varphi_2) &\le \|\varphi_3\|_{L^2(\mathbb{R}^n \setminus B(\mathbf{O}_3, R))} (\|\varphi_1\|_{L^2(\mathbb{R}^n)} + \|\varphi_2\|_{L^2(\mathbb{R}^n)}) \nonumber \\ &\quad + \|\varphi_3\|_{L^2(\mathbb{R}^n)} (\|\varphi_1\|_{L^2(\mathbb{R}^n \setminus B(\mathbf{O}_1, R))} + \|\varphi_2\|_{L^2(\mathbb{R}^n \setminus B(\mathbf{O}_2, R))}) \nonumber \\ \tag{72} &\le \sqrt{\frac{2(\tilde{\lambda}_2 + r_{\max})}{\alpha R^2}} (\|\varphi_1\|_{L^2(\mathbb{R}^n)} + \|\varphi_2\|_{L^2(\mathbb{R}^n)} + 2\|\varphi_3\|_{L^2(\mathbb{R}^n)}) \nonumber \\ &\le \sqrt{\frac{12(\tilde{\lambda}_2 + r_{\max})}{\alpha R^2}}, \end{align} $$

where the last inequality uses $\int_{\mathbb{R}^n} \|\Phi\|^2 = \int_{\mathbb{R}^n} \varphi_1^2 + \int_{\mathbb{R}^n} \varphi_2^2 + \int_{\mathbb{R}^n} \varphi_3^2 = 1$ together with the Cauchy-Schwarz inequality $a + b + 2c \le \sqrt{6(a^2 + b^2 + c^2)}$ for any $(a, b, c) \in \mathbb{R}^3$.
We come back to (71). Since $\int_{\mathbb{R}^n} \varphi_3^2 \in (0, 1)$, we have

$$ \lambda_3 \geq \inf_{a \in (0,1)} \{ (\tilde{\lambda}_2 - \delta)(1-a) + \lambda_1 a \} + \delta - \delta \int_{\mathbb{R}^n} \varphi_3 (\varphi_1 + \varphi_2). $$

Using (27), namely $\tilde{\lambda}_2 - \delta \leq \lambda_1$, we get

$$ \lambda_3 \geq \tilde{\lambda}_2 - \delta \int_{\mathbb{R}^n} \varphi_3 (\varphi_1 + \varphi_2). $$

Finally, we recall that $R = \min(\|\mathbf{O}_3 - \mathbf{O}_1\|/2, \|\mathbf{O}_3 - \mathbf{O}_2\|/2)$ in (72) and, together with (70), we get (28). The last conclusion of Proposition 4.5 then follows from (27). □

**Proof of Proposition 4.6.** We denote $(\varphi_1, \varphi_2)$ the normalised principal eigenvector associated with $\lambda_2 = \lambda_2(\mathbf{O}_1, \mathbf{O}_2)$. By symmetry we have $\varphi_1 = \varphi_2 \circ \iota$, with $\iota$ as in (65), hence $\int_{\mathbb{R}^n} \varphi_1^2 = \int_{\mathbb{R}^n} \varphi_2^2 = 1/2$. Next, we take $\mathbf{O}_3 \in B(\mathbf{O}_2, \rho)$ for some $\rho > 0$ (the case $\mathbf{O}_3 \in B(\mathbf{O}_1, \rho)$ can be handled similarly), we call $\lambda_3 = \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$, and we test the Rayleigh quotient $Q_3$ at $\Phi := (a\varphi_1, a\varphi_2, b\varphi_2)$ with

$$ a^2 + \frac{b^2}{2} = 1, $$

so that $\Phi$ is normalised. We obtain after some straightforward computations

$$ \begin{aligned} \lambda_3 \le Q_3(\Phi) &= a^2\lambda_2 + b^2 \left( \frac{\mu^2}{2} \int_{\mathbb{R}^n} \|\nabla\varphi_2\|^2 - r_2\varphi_2^2 \right) + b^2 \int_{\mathbb{R}^n} (r_2 - r_3)\varphi_2^2 \\ &\quad + \delta \left( \int_{\mathbb{R}^n} b^2\varphi_2^2 + a^2\varphi_1\varphi_2 - ab\varphi_1\varphi_2 - ab\varphi_2^2 \right) \\ &= a^2\lambda_2 + b^2 \left( \lambda_2 \int_{\mathbb{R}^n} \varphi_2^2 + \delta \int_{\mathbb{R}^n} (\varphi_1 - \varphi_2)\varphi_2 \right) + b^2 \int_{\mathbb{R}^n} (r_2 - r_3)\varphi_2^2 \\ &\quad + \delta \left( \int_{\mathbb{R}^n} b^2\varphi_2^2 + a^2\varphi_1\varphi_2 - ab\varphi_1\varphi_2 - ab\varphi_2^2 \right) \\ &= \lambda_2 + \delta \left( (b^2 + a^2 - ab) \int_{\mathbb{R}^n} \varphi_1\varphi_2 - ab \int_{\mathbb{R}^n} \varphi_2^2 \right) + b^2 \int_{\mathbb{R}^n} (r_2 - r_3)\varphi_2^2. \end{aligned} $$

We now select $a = b = \sqrt{2/3}$ and thus obtain

$$ \lambda_3 - \lambda_2 \le \frac{2\delta}{3} \left( \int_{\mathbb{R}^n} \varphi_1\varphi_2 - \int_{\mathbb{R}^n} \varphi_2^2 \right) + \frac{2}{3} \int_{\mathbb{R}^n} (r_2 - r_3)\varphi_2^2. $$

Since the Cauchy-Schwarz inequality yields $\int_{\mathbb{R}^n} \varphi_1\varphi_2 - \int_{\mathbb{R}^n} \varphi_2^2 < 0$ ($\varphi_1$ and $\varphi_2$ are not colinear because $\beta \ne 0$) and since the dominated convergence theorem implies that $\int_{\mathbb{R}^n} (r_2 - r_3)\varphi_2^2 \to 0$ as $\rho \to 0$, we have $\lambda_3 - \lambda_2 < 0$ if $\rho > 0$ is small enough. □

**Proof of Proposition 4.7.** First, from Proposition 3.2 and Corollary 3.3, we know that, whatever the position of $\mathbf{O}_3$, we have

$$ (73) \quad \lambda_1 \le \lambda_1 + \delta \left( 1 - \sum_{1 \le i < j \le 3} \int_{\mathbb{R}^n} \varphi_i(\mathbf{x})\varphi_j(\mathbf{x})d\mathbf{x} \right) \le \lambda_3 \le \lambda_1 + \delta \left( 1 - \frac{e^{-\sqrt{\alpha}\beta^2/\mu}}{3} \right). $$
Set $\psi_i = \varphi_i(\cdot + \mathbf{O}_i)$. We have

$$
\begin{align*}
\lambda_3 &= -r_{\max} + \sum_{i=1}^{3} \int_{\mathbb{R}^n} \left( \frac{\mu^2}{2} \|\nabla \psi_i(\mathbf{x})\|^2 + \frac{\alpha}{2} \|\mathbf{x}\|^2 (\psi_i(\mathbf{x}))^2 \right) \mathrm{d}\mathbf{x} + \delta - \delta \sum_{1 \le i < j \le 3} \int_{\mathbb{R}^n} \varphi_i(\mathbf{x}) \varphi_j(\mathbf{x}) \mathrm{d}\mathbf{x} \\
&\ge -r_{\max} + \sum_{i=1}^{3} \int_{\mathbb{R}^n} \left( \frac{\mu^2}{2} \|\nabla \psi_i(\mathbf{x})\|^2 + \frac{\alpha}{2} \|\mathbf{x}\|^2 (\psi_i(\mathbf{x}))^2 \right) \mathrm{d}\mathbf{x}.
\end{align*}
$$

Thus, from (73), we obtain that

$$
\sum_{i=1}^{3} \int_{\mathbb{R}^n} \frac{\alpha}{2} \|\mathbf{x}\|^2 (\psi_i(\mathbf{x}))^2 \mathrm{d}\mathbf{x} \leq r_{\max} + \lambda_1 + \delta \left( 1 - \frac{e^{-\sqrt{\alpha}\beta^2/\mu}}{3} \right).
$$

Considering any radius $R > 0$, the last inequality implies that

$$
\|\psi_i\|_{L^2(\mathbb{R}^n \setminus B(\mathcal{O}, R))}^2 \le \frac{2(r_{\max} + \lambda_1 + \delta)}{\alpha R^2}
$$

for every $\beta > 0$ and $1 \le i \le 3$.

Now, consider the particular case where $\mathbf{O}_3 = \mathcal{O} = (0, \dots, 0)$, so that the Euclidean
distance between two optima is at least $\beta$. In this case, $\psi_3 = \varphi_3$ and we have, for every
$R > 0$ and $\beta \ge 2R$,

$$
\begin{align*}
\int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_3(\mathbf{x}) \, \mathrm{d}\mathbf{x} &= \int_{\mathbb{R}^n} \psi_1(\mathbf{x} - \mathbf{O}_1) \psi_3(\mathbf{x}) \, \mathrm{d}\mathbf{x} \\
&= \int_{B(\mathcal{O}, R)} \psi_1(\mathbf{x} - \mathbf{O}_1) \psi_3(\mathbf{x}) \, \mathrm{d}\mathbf{x} + \int_{\mathbb{R}^n \setminus B(\mathcal{O}, R)} \psi_1(\mathbf{x} - \mathbf{O}_1) \psi_3(\mathbf{x}) \, \mathrm{d}\mathbf{x} \\
&\le \|\psi_1\|_{L^2(\mathcal{O}_1, R)} \|\psi_3\|_{L^2(\mathbb{R}^n)} + \|\psi_1\|_{L^2(\mathbb{R}^n)} \|\psi_3\|_{L^2(\mathbb{R}^n \setminus B(\mathcal{O}, R))} \\
&\le \|\psi_1\|_{L^2(\mathbb{R}^n \setminus B(\mathcal{O}, R))} \|\psi_3\|_{L^2(\mathbb{R}^n)} + \|\psi_1\|_{L^2(\mathbb{R}^n)} \|\psi_3\|_{L^2(\mathbb{R}^n \setminus B(\mathcal{O}, R))} \\
&\le \|\psi_1\|_{L^2(\mathbb{R}^n \setminus B(\mathcal{O}, R))} + \|\psi_3\|_{L^2(\mathbb{R}^n \setminus B(\mathcal{O}, R))}
\end{align*}
$$

since $B(\mathbf{O}_1, R) \subset \mathbb{R}^n \setminus B(\mathcal{O}, R)$ and $\|\psi_i\|_{L^2(\mathbb{R}^n)} = \|\varphi_i\|_{L^2(\mathbb{R}^n)} \le 1$ for each $1 \le i \le 3$. As a
result

$$
\int_{\mathbb{R}^n} \varphi_1(\mathbf{x}) \varphi_3(\mathbf{x}) d\mathbf{x} \le \frac{\sqrt{8(r_{\max} + \lambda_1 + \delta)}}{R\sqrt{\alpha}}.
$$

Applying the same arguments, we finally get

$$
\sum_{1 \le i < j \le 3} \int_{\mathbb{R}^n} \varphi_i(\mathbf{x}) \varphi_j(\mathbf{x}) \, \mathrm{d}\mathbf{x} \le \frac{\sqrt{72(r_{\max} + \lambda_1 + \delta)}}{R\sqrt{\alpha}}.
$$

Together with (73), this shows that

$$
\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathcal{O}) \to \lambda_1 + \delta \text{ as } \beta \to +\infty.
$$

Consider now the case $\mathbf{O}_3 = \mathbf{O}_1$. We note that $\lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_1) = \lambda_3(\mathbf{O}_1, \mathbf{O}_1, \mathbf{O}_2)$.
Then, using the result of Proposition 4.5, we obtain that $\lambda_3(\mathbf{O}_1, \mathbf{O}_1, \mathbf{O}_2) \to \lambda_2(\mathbf{O}_1, \mathbf{O}_1)$
as $\beta \to +\infty$, where $\tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_1)$ is defined by (26), in the particular case where the two
optima are at the same position. In such case, $\tilde{\varphi}_1 = \tilde{\varphi}_2$ by uniqueness and therefore, from
(26),

$$
-\frac{\mu^2}{2} \Delta \tilde{\varphi}_1 - r_1(\mathbf{x}) \tilde{\varphi}_1 + \frac{\delta}{2} \tilde{\varphi}_1 = \tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_1) \tilde{\varphi}_1,
$$

so that $\tilde{\lambda}_2(\mathbf{O}_1, \mathbf{O}_1) = \lambda_1 + \delta/2$, and we are done.

<figure><img src="image_4.png" /></figure>
**Remark 6.3.** The conclusion $\lim_{\beta \to +\infty} \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_1) = \lambda_1 + \delta/2$ of Proposition 4.7 looks at first glance similar to the conclusion (54) derived in the proof of Proposition 2.3 (iii). However, the arguments used in both proofs are not the same, since the proof of Proposition 4.7 mainly relies on $L^2$ estimates, while that of Proposition 2.3 (iii) is based on $L^1$ estimates, with a reasoning done by contradiction and actually leading to the contradicting conclusions (53) and (54). Notice also that the system (26) used in the proof of Proposition 4.7 is symmetric in $(\tilde{\varphi}_1, \tilde{\varphi}_2)$, whereas the system (49) derived in the argument by contradiction in the proof of Proposition 2.3 (iii) is not symmetric in $(\varphi_1, \varphi_2)$, complicating its analysis.

*Proof of Proposition 4.8.* We need to recall the notion of Steiner symmetrisation of a function with respect to the variable $x_k$. Consider first a measurable function $h : \mathbb{R} \to \mathbb{R}$, $x \mapsto h(x)$, which is either nonnegative and belongs to $L^p(\mathbb{R})$ for some $1 \le p < \infty$, or which is such that $h(x) \to \inf_{\mathbb{R}} h \in [-\infty, +\infty)$ as $|x| \to +\infty$. Then there exists a unique (in the class of functions which are identical almost everywhere) function $h^\sharp : \mathbb{R} \to \mathbb{R}$, $x \mapsto h^\sharp(x)$, such that: (i) $h^\sharp$ is symmetric with respect to $x = 0$ and nonincreasing in $[0, +\infty)$, i.e., for all $x, y \in \mathbb{R}$, $h^\sharp(x) \ge h^\sharp(y)$ if $|x| \le |y|$; (ii) $h^\sharp$ has the same distribution function as $h$, that is,

$$ \mathrm{meas}\{x \in \mathbb{R} : h^{\sharp}(x) > \alpha\} = \mathrm{meas}\{x \in \mathbb{R} : h(x) > \alpha\} $$

for all $\alpha \in \mathbb{R}$, where meas denotes the one-dimensional Lebesgue measure.

Consider now a measurable function $h : \mathbb{R}^n \to \mathbb{R}$, $\mathbf{x} \mapsto h(\mathbf{x})$, which is either non-negative and belongs to $L^p(\mathbb{R}^n)$ for some $1 \le p < \infty$, or which is such that $h(\mathbf{x}) \to \inf_{\mathbb{R}^n} h \in [-\infty, +\infty)$ as $\|\mathbf{x}\| \to +\infty$. For $1 \le k \le n$ and $(x_1, \dots, x_{k-1}, x_{k+1}, \dots, x_n) \in \mathbb{R}^{n-1}$, the measurable function $x_k \mapsto h(x_1, \dots, x_k, \dots, x_n)$ is either nonnegative and belongs to $L^p(\mathbb{R})$ (for almost every $(x_1, \dots, x_{k-1}, x_{k+1}, \dots, x_n) \in \mathbb{R}^{n-1}$), or is such that $h(x_1, \dots, x_k, \dots, x_n) \to \inf_{\mathbb{R}^n} h = \inf_{\mathbb{R}} h(x_1, \dots, x_{k-1}, \dots, x_{k+1}, \dots, x_n) \in [-\infty, +\infty)$ as $|x_k| \to +\infty$. We can then rearrange the function $x_k \mapsto h(x_1, \dots, x_k, \dots, x_n)$ as above. This corresponds to the Steiner rearrangement of $h$ with respect to the variable $x_k$, and we denote it by $h^{\sharp_k} : \mathbb{R}^n \to \mathbb{R}$, $\mathbf{x} \mapsto h^{\sharp_k}(\mathbf{x})$.

In the sequel, we use some known properties of this rearrangement, which we briefly recall. Firstly, for any nonnegative function $h \in L^p(\mathbb{R}^n)$ with $1 \le p < \infty$, the nonnegative function $h^{\sharp_k}$ belongs to $L^p(\mathbb{R}^n)$ too, and

$$ (74) \qquad \int_{\mathbb{R}^n} h^p = \int_{\mathbb{R}^n} (h^{\sharp_k})^p. $$

Secondly, for any nonnegative $L^2(\mathbb{R}^n)$ functions $h$ and $j$, the Hardy-Littlewood inequality asserts that

$$ (75) \qquad 0 \le \int_{\mathbb{R}^n} h \, j \le \int_{\mathbb{R}^n} h^{\sharp_k} \, j^{\sharp_k}. $$

Thirdly, the Pólya-Szegö inequality says that, for any nonnegative function $h \in W^{1,p}(\mathbb{R}^n)$ with $1 \le p < \infty$, the function $h^{\sharp_k}$ belongs to $W^{1,p}(\mathbb{R}^n)$ too, and

$$ (76) \qquad \int_{\mathbb{R}^n} \|\nabla h\|^p \ge \int_{\mathbb{R}^n} \|\nabla h^{\sharp_k}\|^p. $$

Equipped with this, we are in position to prove the result in Proposition 4.8. For the proof, without loss of generality; using the notation (19) and the fact that the map $\mathbf{O}_3 \mapsto \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$ is invariant by rotation around the axis $\mathbb{R} \times \{0\}^{n-1}$, we can assume
that $\mathbf{O}_3 \in \mathbb{R}^2 \times \{0\}^{n-2}$. Let $(\varphi_1, \varphi_2, \varphi_3) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^3$ be the normalised positive eigenvector associated with $\lambda_3 = \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3)$. Remember from (9) that

$$ (77) \qquad \lambda_3 = \sum_{i=1}^{3} \int_{\mathbb{R}^n} \left( \frac{\mu^2}{2} \|\nabla \varphi_i\|^2 - r_i \varphi_i^2 \right) + \delta - \delta \int_{\mathbb{R}^n} (\varphi_1 \varphi_2 + \varphi_2 \varphi_3 + \varphi_1 \varphi_3). $$

Consider $(\varphi_1^{\sharp 2}, \varphi_2^{\sharp 2}, \varphi_3^{\sharp 2})$ the Steiner symmetrisations of the eigenfunctions with respect to the variable $x_2$. From (74), we know that

$$ (78) \qquad \int_{\mathbb{R}^n} (\varphi_1^{\sharp 2})^2 + (\varphi_2^{\sharp 2})^2 + (\varphi_3^{\sharp 2})^2 = \int_{\mathbb{R}^n} \varphi_1^2 + \varphi_2^2 + \varphi_3^2 = 1. $$

Additionally, from (75), we have

$$ (79) \qquad \int_{\mathbb{R}^n} (\varphi_1 \varphi_2 + \varphi_2 \varphi_3 + \varphi_1 \varphi_3) \le \int_{\mathbb{R}^n} (\varphi_1^{\sharp 2} \varphi_2^{\sharp 2} + \varphi_2^{\sharp 2} \varphi_3^{\sharp 2} + \varphi_1^{\sharp 2} \varphi_3^{\sharp 2}). $$

Now, we have to prove that, for each $i = 1, 2, 3$, the function $\varphi_i^{\sharp 2}$ belongs to $H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$. First of all, it belongs to $H^1(\mathbb{R}^n)$, from (74) and (76). Let us now show that the integral $\int_{\mathbb{R}^n} \|\mathbf{x}\|^2 (\varphi_i^{\sharp 2}(\mathbf{x}))^2 \mathrm{d}\mathbf{x}$ converges. To do so, for any $R > 0$, call $b_R$ the nonnegative $L^2(\mathbb{R}^n)$ function defined by $b_R(\mathbf{x}) := \max(R^2 - \|\mathbf{x}\|^2, 0)$ and observe that

$$ \begin{aligned} \int_{\mathbb{R}^n} \min(\|\mathbf{x}\|^2, R^2) (\varphi_i(\mathbf{x}))^2 \mathrm{d}\mathbf{x} &= \int_{\mathbb{R}^n} R^2 (\varphi_i(\mathbf{x}))^2 \mathrm{d}\mathbf{x} - \int_{\mathbb{R}^n} b_R(\mathbf{x}) (\varphi_i(\mathbf{x}))^2 \mathrm{d}\mathbf{x} \\ &\ge \int_{\mathbb{R}^n} R^2 (\varphi_i^{\sharp 2}(\mathbf{x}))^2 \mathrm{d}\mathbf{x} - \int_{\mathbb{R}^n} b_R^{\sharp 2}(\mathbf{x}) (\varphi_i^{\sharp 2}(\mathbf{x}))^2 \mathrm{d}\mathbf{x} \end{aligned} $$

from (74)-(75) (we here use the fact that $\varphi_i^2$ belongs to $L^2(\mathbb{R}^n)$, since it is continuous and decays to 0 faster than exponentially as $\|\mathbf{x}\| \to +\infty$). Since $b_R^{\sharp 2} = b_R$, one then gets that

$$ \begin{aligned} \int_{\mathbb{R}^n} \|\mathbf{x}\|^2 (\varphi_i(\mathbf{x}))^2 \mathrm{d}\mathbf{x} &\ge \int_{\mathbb{R}^n} \min(\|\mathbf{x}\|^2, R^2) (\varphi_i(\mathbf{x}))^2 \mathrm{d}\mathbf{x} \\ &\ge \int_{\mathbb{R}^n} (R^2 - b_R(\mathbf{x})) (\varphi_i^{\sharp 2}(\mathbf{x}))^2 \mathrm{d}\mathbf{x} = \int_{\mathbb{R}^n} \min(\|\mathbf{x}\|^2, R^2) (\varphi_i^{\sharp 2}(\mathbf{x}))^2 \mathrm{d}\mathbf{x}. \end{aligned} $$

The monotone convergence theorem then yields the convergence of $\int_{\mathbb{R}^n} \|\mathbf{x}\|^2 (\varphi_i^{\sharp 2}(\mathbf{x}))^2 \mathrm{d}\mathbf{x}$. Thus, $\varphi_i^{\sharp 2} \in H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$, for each $i = 1, 2, 3$.

Next, we claim that, for each $i = 1, 2, 3$, the Hardy-Littlewood inequality can be applied to the couple $(r_i, \varphi_i^2)$, although $r_i$ does not have a constant sign in $\mathbb{R}^n$, namely we claim that

$$ (80) \qquad \int_{\mathbb{R}^n} r_i \varphi_i^2 \le \int_{\mathbb{R}^n} r_i^{\sharp 2} (\varphi_i^2)^{\sharp 2} = \int_{\mathbb{R}^n} r_i^{\sharp 2} (\varphi_i^{\sharp 2})^2. $$

Notice first that all integrals in (80) converge since the functions $\varphi_i$ and $\varphi_i^{\sharp 2}$ belong to $L^2(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n)$ and since $(r_1^{\sharp 2}, r_2^{\sharp 2}, r_3^{\sharp 2}) = (r_1, r_2, r_{\mathbf{O}_3^{\sharp}})$, where

$$ r_{\mathbf{O}_3^{\sharp}}(\mathbf{x}) = r_{\max} - \alpha \frac{\|\mathbf{x} - \mathbf{O}_3^{\sharp}\|^2}{2} $$

(the equality $r_3^{\sharp 2} = r_{\mathbf{O}_3^{\sharp}}$ holds because $\mathbf{O}_3 \in \mathbb{R}^2 \times \{0\}^{n-2}$). To show (80), set, for any arbitrary $K > 0$, $r_i^K(\mathbf{x}) := r_i(\mathbf{x}) + K$ if $r_i(\mathbf{x}) + K > 0$ and $r_i^K(\mathbf{x}) := 0$ otherwise, that
is, $r_i^K(\mathbf{x}) = \max(r_i(\mathbf{x}) + K, 0)$. The function $r_i^K$ is continuous and compactly supported, hence it is in $L^2(\mathbb{R}^n)$. From (75) applied to the nonnegative $L^2(\mathbb{R}^n)$ functions $r_i^K$ and $\varphi_i^2$, together with the property $(\varphi_i^2)^{\sharp 2} = (\varphi_i^{\sharp 2})^2$, we infer that

$$0 \le \int_{\mathbb{R}^n} r_i^K \varphi_i^2 \le \int_{\mathbb{R}^n} (r_i^K)^{\sharp 2} (\varphi_i^{\sharp 2})^2.$$

As $(r_i^K)^{\sharp 2} = \max(r_i^{\sharp 2} + K, 0)$ by definition of $r_i$ in (2), we get that

$$(81) \quad \int_{r_i+K>0} r_i \varphi_i^2 + K \int_{r_i+K>0} \varphi_i^2 - K \int_{r_i^{\sharp 2}+K>0} (\varphi_i^{\sharp 2})^2 \le \int_{r_i^{\sharp 2}+K>0} r_i^{\sharp 2} (\varphi_i^{\sharp 2})^2.$$

It also follows from (74) that the functions $\varphi_i$ and $\varphi_i^{\sharp 2}$ have the same $L^2(\mathbb{R}^n)$ norm, hence

$$(82) \quad K \int_{r_i+K>0} \varphi_i^2 - K \int_{r_i^{\sharp 2}+K>0} (\varphi_i^{\sharp 2})^2 = -K \int_{r_i+K \le 0} \varphi_i^2 + K \int_{r_i^{\sharp 2}+K \le 0} (\varphi_i^{\sharp 2})^2.$$

Next, we note that there are some positive constants $C$ and $K_0$ such that

$$\{\mathbf{x} \in \mathbb{R}^n : r_i(\mathbf{x}) + K \le 0\} \cup \{\mathbf{x} \in \mathbb{R}^n : r_i^{\sharp 2}(\mathbf{x}) + K \le 0\} \subset \mathbb{R}^n \setminus B(\mathcal{O}, C\sqrt{K})$$

for all $K \ge K_0$ and $1 \le i \le 3$. Thus, (82) implies that, for all $K \ge K_0$,

$$(83) \quad \left| K \int_{r_i+K>0} \varphi_i^2 - K \int_{r_i^{\sharp 2}+K>0} (\varphi_i^{\sharp 2})^2 \right| \le K \int_{\mathbb{R}^n \setminus B(\mathcal{O}, C\sqrt{K})} \varphi_i^2 + (\varphi_i^{\sharp 2})^2.$$

We then use the fact that the functions $x \mapsto \|\mathbf{x}\| \varphi_i(\mathbf{x})$ and $x \mapsto \|\mathbf{x}\| \varphi_i^{\sharp 2}(\mathbf{x})$ are in $L^2(\mathbb{R}^n)$ (the property for $\varphi_i^{\sharp 2}$ follows from the previous paragraph). Hence,

$$C^2 K \int_{\mathbb{R}^n \setminus B(\mathcal{O}, C\sqrt{K})} (\varphi_i(\mathbf{x}))^2 d\mathbf{x} \le \int_{\mathbb{R}^n \setminus B(\mathcal{O}, C\sqrt{K})} \|\mathbf{x}\|^2 (\varphi_i(\mathbf{x}))^2 d\mathbf{x} \to 0 \quad \text{as } K \to +\infty,$$

and a similar property holds for $\varphi_i^{\sharp 2}$. Finally, (83) shows that

$$K \int_{r_i+K>0} \varphi_i^2 - K \int_{r_i^{\sharp 2}+K>0} (\varphi_i^{\sharp 2})^2 \to 0 \quad \text{as } K \to +\infty.$$

Passing to the limit $K \to +\infty$ in (81) and using the dominated convergence theorem, we obtain (80).

Finally, using (76), we get that

$$(84) \quad \sum_{i=1}^3 \int_{\mathbb{R}^n} \|\nabla \varphi_i\|^2 \ge \sum_{i=1}^3 \int_{\mathbb{R}^n} \|\nabla \varphi_i^{\sharp 2}\|^2.$$

Coming back to (77), and using (79), (80) and (84), it follows that

$$(85) \quad \lambda_3 \ge \sum_{i=1}^3 \int_{\mathbb{R}^n} \left( \frac{\mu^2}{2} \|\nabla \varphi_i^{\sharp 2}\|^2 - r_i^{\sharp 2} (\varphi_i^{\sharp 2})^2 \right) + \delta - \delta \int_{\mathbb{R}^n} (\varphi_1^{\sharp 2} \varphi_2^{\sharp 2} + \varphi_2^{\sharp 2} \varphi_3^{\sharp 2} + \varphi_1^{\sharp 2} \varphi_3^{\sharp 2}).$$

Using (78) and $(\varphi_1^{\sharp 2}, \varphi_2^{\sharp 2}, \varphi_3^{\sharp 2}) \in (H^1(\mathbb{R}^n) \cap L_w^2(\mathbb{R}^n))^3$, together with $(r_1^{\sharp 2}, r_2^{\sharp 2}, r_3^{\sharp 2}) = (r_1, r_2, r_{\mathbf{O}_3^{\sharp}})$, it follows that $(\varphi_1^{\sharp 2}, \varphi_2^{\sharp 2}, \varphi_3^{\sharp 2})$ is an admissible triplet in the Rayleigh quotient associated with the optima $(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3^{\sharp})$. As a conclusion, (85) yields $\lambda_3 \ge \lambda_3(\mathbf{O}_1, \mathbf{O}_2, \mathbf{O}_3^{\sharp})$, which proves Proposition 4.8. $\square$
References

[1] S. Agmon. *Lectures on exponential decay of solutions of second-order elliptic equations: bounds on eigenfunctions of N-body Schrödinger operators*, vol. 29 of *Mathematical Notes*. Princeton University Press, Princeton, NJ, 1982.

[2] M. Alfaro and M. Veruete. Evolutionary branching via replicator-mutator equations, *J. Dynamics Diff. Equations*, 31:2029–2052, 2019.

[3] P. Besala. Fundamental solution and Cauchy problem for a parabolic system with unbounded coefficients. *J. Diff. Equations*, 33:26–38, 1979.

[4] T. Caquet, C. Gascuel, and M. Tixier-Boichard. *Agroécologie : des recherches pour la transition des filières et des territoires*. Quae, 2020.

[5] L. Cardoulis. Systems of Schrödinger equations in the whole space. *Rostock. Math. Kolloq.*, 65:29–49, 2009.

[6] L. Cardoulis. Principal eigenvalues for systems of Schrödinger equations defined in the whole space with indefinite weights. *Math. Slovaca*, 65:1079, 2015.

[7] L. Damascelli and F. Pacella. Symmetry results for cooperative elliptic systems via linearization. *SIAM J. Math. Anal.*, 45:1003–1026, 2013.

[8] F. Débarre, O. Ronce, and S. Gandon. Quantifying the effects of migration and mutation on adaptation and demography in spatially heterogeneous environments. *J. Evolutionary Biology*, 26:1185–1202, 2013.

[9] R. Djidjou-Demasse, A. Ducrot, and F. Fabre. Steady state concentration for a phenotypic structured problem modeling the evolutionary epidemiology of spore producing pathogens. *Math. Models Methods Appl. Sci.*, 27:385–426, 2017.

[10] Food and Agriculture Organization of the United Nations. *The 10 Elements of Agroecology. Guiding the Transition to Sustainable Food and Agricultural Systems*. FAO, 2018.

[11] G. García-Ramos and M. Kirkpatrick. Genetic models of adaptation and gene flow in peripheral populations. *Evolution*, 51:21–28, 1997.

[12] F. Hamel, F. Lavigne, and L. Roques. Adaptation in a heterogeneous environment. I: Persistence versus extinction. *J. Math. Biology*, 83:14, 2021.

[13] A. Latinne, B. Hu, K. J. Olival, G. Zhu, L. Zhang, H. Li, A. A. Chmura, H. E. Field, C. Zambrana-Torrelio, J. H. Epstein, B. Li, W. Zhang, L.-F. Wang, Z.-L. Shi, and P. Daszak. Origin and cross-species transmission of bat coronaviruses in China. *Nature Communications*, 11:1–15, 2020.

[14] S. K. P. Lau, P. C. Y. Woo, K. S. M. Li, Y. Huang, H.-W. Tsoi, B. H. L. Wong, S. S. Y. Wong, S.-Y. Leung, K.-H. Chan, and K.-Y. Yuen. Severe acute respiratory syndrome coronavirus-like virus in Chinese horseshoe bats. *Proc. Nat. Acad. Sciences*, 102:14040–14045, 2005.

[15] F. Lavigne, G. Martin, Y. Anciaux, J. Papaïx, and L. Roques. When sinks become sources: adaptive colonization in asexuals. *Evolution*, 74:29–42, 2020.

[16] G. Martin and T. Lenormand. The fitness effect of mutations across environments: Fisher's geometrical model with multiple optima. *Evolution*, 69:1433–1447, 2015.

[17] G. Meszéna, I. Czibula, and S. Geritz. Adaptive dynamics in a 2-patch environment: a toy model for allopatric and parapatric speciation. *J. Biol. Systems*, 5:265–284, 1997.

[18] S. Mirrahimi. A Hamilton–Jacobi approach to characterize the evolutionary equilibria in heterogeneous environments. *Math. Models Methods Appl. Sci.*, 27:2425–2460, 2017.
[19] S. Mirrahimi and S. Gandon. Evolution of specialization in heterogeneous environments: equilibrium between selection, mutation and migration. *Genetics*, 214:479-491, 2020.

[20] J. Papaïx, O. David, C. Lannou, and H. Monod. Dynamics of adaptation in spatially heterogeneous metapopulations. *PloS one*, 8:e54697, 2013.

[21] C. R. Parrish, P. R. Murcia, and E. C. Holmes. Influenza virus reservoirs and intermediate hosts: dogs, horses, and new possibilities for influenza virus exposure of humans. *J. Virology*, 89:2990-2994, 2015.

[22] A. C. Ponce. *Elliptic PDEs, measures and capacities, from the Poisson equation to nonlinear Thomas-Fermi problems*, vol. 23 of *Tracts in Mathematics*. European Math. Society, 2016.

[23] O. Tenaillon. The utility of Fisher's geometric model in evolutionary genetics. *Annual Rev. Ecology, Evolution, and Systematics*, 45:179-201, 2014.

[24] D. Waxman and J. R. Peck. Pleiotropy and the preservation of perfection. *Science*, 279:1210-1213, 1998.

[25] R. J. Webby and R. G. Webster. Emergence of influenza A viruses. *Phil. Trans. Royal Soc. London, Ser. B: Biol. Sciences*, 356:1817-1828, 2001.