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

HAL
open science

# Steady-state solutions for a reaction-diffusion equation with Robin boundary conditions: Application to the control of dengue vectors

Luís Almeida, Pierre-Alexandre Bliman, Nga Nguyen, Nicolas Vauchelet

► To cite this version:

Luís Almeida, Pierre-Alexandre Bliman, Nga Nguyen, Nicolas Vauchelet. Steady-state solutions for a reaction-diffusion equation with Robin boundary conditions: Application to the control of dengue vectors. 2023. [hal-03775887v2](https://hal.science/hal-03775887v2)

HAL Id: hal-03775887

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

Preprint submitted on 25 Sep 2023

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Steady-state solutions for a reaction-diffusion equation with Robin boundary conditions: Application to the control of dengue vectors.

L. Almeida<fnref n="1" />, P.A. Bliman<fnref n="2" />,<fnref n="1" />, N. Nguyen<fnref n="3" />,<fnref n="1" />, N. Vauchelet<fnref n="3" />

## Abstract

In this paper, we investigate an initial-boundary-value problem of a reaction-diffusion equation in a bounded domain with a Robin boundary condition and introduce some particular parameters to consider the non-zero flux on the boundary. This problem arises in the study of mosquito populations under the intervention of the population replacement method, where the boundary condition takes into account the inflow and outflow of individuals through the boundary. Using phase-plane analysis, the present paper studies the existence and properties of non-constant steady-state solutions depending on several parameters. Then, we prove some sufficient conditions for their stability. We show that the long-time efficiency of this control method depends strongly on the size of the treated zone and the migration rate. To illustrate these theoretical results, we provide some numerical simulations in the framework of mosquito population control.

## 1 Introduction

The study of scalar reaction-diffusion equations $\partial_t p - \Delta p = f(p)$ with a given nonlinearity $f$ has a long history. For suitable choices of $f$, this equation can be used to model some phenomena in biology such as population dynamics (see e.g. [4], [16], [25]). To investigate the structure of the steady-state solutions, the semilinear elliptic equation $\Delta p + f(p) = 0$ has been studied extensively.

Many results about the multiplicity of positive solutions for the parametrized version $\Delta p + \lambda f(p) = 0$ in a bounded domain are known. Here, $\lambda$ is a positive parameter. Various works investigated the number of solutions and the global bifurcation diagrams of this equation according to different classes of the nonlinearity $f$ and boundary conditions. For Dirichlet problems, in [15], Lions used many “bifurcation diagrams” to describe the solution set of this equation with several kinds of nonlinearities $f$, and gave nearly optimal multiplicity results in each case. The exact number of solutions and the precise bifurcation diagrams with cubic-like nonlinearities $f$ were given in the works of Korman et. al. [13], [14], Ouyang and Shi [18] and references therein. In these works, the authors developed a global bifurcation approach to obtain the exact multiplicity of positive solutions. In the case of one-dimensional space with a two-point boundary, Korman gave a survey of this approach in [12]. Another approach was given by Smoller and Wasserman in [24] using phase-plane analysis and the time mapping method. This method was completed and applied in the works of Wang [28], [29]. While the bifurcation approach is convenient to solve the problem with more general cubic nonlinearities $f$, the phase-plane method is more intuitive and easier to compute.

Although many results were obtained concerning the number of solutions for Dirichlet problems, relatively little seems to be known concerning the results for other kinds of boundary conditions. For the Neumann problem, the works of Smoller and Wasserman [24], Schaaf [21], and Korman [11] dealt with cubic-like nonlinearities $f$ in one dimension. Recently, more works have been done for Robin boundary conditions (see e.g. [3], [22], [33]), Neumann-Robin boundary conditions (see e.g. [27]), or even nonlinear boundary conditions (see e.g. [6], [7] and references therein). However, those works only focused on

<footnote n="1">¹LJLL, Sorbonne University, 5 Place Jussieu, 75005 Paris, CNRS University Paris Cite, France</footnote>

<footnote n="2">²MAMBA, Inria Paris; LJLL, Sorbonne University, 5 Place Jussieu, 75005 Paris, France</footnote>

<footnote n="3">³LAGA, CNRS UMR 7539, Institut Galilée, Université Sorbonne Paris Nord, 99 avenue Jean-Baptiste Clément, 93430 Villetaneuse, France</footnote>
other types of nonlinearities such as positive or monotone $f$. An analogous problem with advection term was studied in [31], [30] for cubic-like nonlinearities, but in these works, they used a homogeneous non-symmetric Robin boundary condition to characterize the open or closed environment boundary. To the best of our knowledge, the study of inhomogeneous symmetric Robin problems with cubic-like nonlinearities remains quite open.

In this paper, we study the steady-state solutions with values in $[0, 1]$ of a reaction-diffusion equation in one dimension with inhomogeneous Robin boundary conditions

$$ \partial_t p^0 - \partial_{xx} p^0 = f(p^0), \quad (t, x) \in (0, T) \times \Omega, \quad (1a) $$

$$ \frac{\partial p^0}{\partial \nu} = -D(p^0 - p^{\text{ext}}), \quad (t, x) \in (0, T) \times \partial\Omega, \quad (1b) $$

$$ p^0(0, x) = p^{\text{init}}(x), \quad x \in \Omega, \quad (1c) $$

where $\Omega = (-L, L)$ is a bounded domain in $\mathbb{R}$, time $T > 0$. The steady-state solutions satisfy the following elliptic boundary-value problem,

$$ -p''(x) = f(p(x)), \quad x \in (-L, L), \quad (2a) $$

$$ p'(L) = -D(p(L) - p^{\text{ext}}), \quad (2b) $$

$$ -p'(-L) = -D(p(-L) - p^{\text{ext}}), \quad (2c) $$

where $L > 0$, $D > 0$, $p^{\text{ext}} \in (0, 1)$ are constants. The reaction term $f: [0, 1] \to \mathbb{R}$ is of class $C^1$, with three roots $\{0, \theta, 1\}$ where $0 < \theta < 1$ (see Figure 1a). The dynamics of (1) can be determined by the structure of steady-state solutions which satisfy (2). Note that, by changing variable from $x$ to $y = x/L$, then (2) becomes $p''(y) + L^2 f(p(y)) = 0$ on $(-1, 1)$ with parameter $L^2$. Thus, we study problem (2) with three parameters $L > 0$, $D > 0$, and $p^{\text{ext}} \in (0, 1)$.

The Robin boundary condition considered in (1) and (2) means that the flow across the boundary points is proportional to the difference between the surrounding density and the density just inside the interval. Here we assume that $p^{\text{ext}}$ does not depend on space variable $x$ nor time variable $t$.

The existence of classical solutions for such problems was studied widely in the theory of elliptic and parabolic differential equations (see, for example, [19]). In our problem, due to difficulties caused by the inhomogeneous Robin boundary condition and the variety of parameters, we cannot obtain the exact multiplicity of solutions. However, our main results in Theorems 2.1 and 2.2 show how the existence of solutions and their “shapes” depend on parameters $D$, $p^{\text{ext}}$ and $L$. The idea of phase plane analysis and time-map method as in [24] are extended to prove these results.

Since the solutions of (2) are equilibria of (1), their stability and instability are the next problems that we want to investigate. The stability analysis of the non-constant steady-state solutions is a delicate problem especially when the system under consideration has multiple steady-state solutions. In Theorem 2.3, we use the principle of linearized stability to give some sufficient conditions for stability. Finally, as a consequence of these theorems, we obtain Corollary 2.1 which provides a comprehensive result about existence and stability of the steady-state solutions when the size $L$ is small.

The main biological application of our results is the control of dengue vectors. *Aedes* mosquitoes are vectors of many vector-borne diseases, including dengue. Recently, a biological control method using an endosymbiotic bacterium called *Wolbachia* has gathered a lot of attention. *Wolbachia* helps reduce the vectorial capacity of mosquitoes and can be passed to the next generation. Massive release of mosquitoes carrying this bacterium in the field is thus considered as a possible method to replace wild mosquitoes and prevent dengue epidemics. Reaction-diffusion equations have been used in previous works to model this replacement strategy (see [1, 2, 26]). In this work, we introduce the Robin boundary condition to describe the migration of mosquitoes through the boundary. Since inflows of wild mosquitoes and outflows of mosquitoes carrying *Wolbachia* may affect the efficiency of the method, the study of existence and stability of steady-state solutions depending on parameters $D$, $p^{\text{ext}}$ and $L$ as in (2), (1) will provide necessary information to maintain the success of the control method using *Wolbachia* under the effects of migration.

Problem (1) arises often in the study of population dynamics. $p^0$ is usually considered as the relative proportion of one population when there are two populations in competition. This is why, we only focus on solutions with values that belong to the interval $[0, 1]$. Problem (1) is derived from the idea in paper [26], where the authors reduce a reaction-diffusion system modeling the competition between two
populations $n_1$ and $n_2$ to a scalar equation on the proportion $p = \frac{n_1}{n_1+n_2}$. More precisely, they consider two populations with a very high fecundity rate scaled by a parameter $\epsilon > 0$ and propose the following system depending on $\epsilon$ for $t > 0, x \in \mathbb{R}^d$,

$$ \partial_t n_1^\epsilon - \Delta n_1^\epsilon = n_1^\epsilon f_1(n_1^\epsilon, n_2^\epsilon), \qquad (3a) $$

$$ \partial_t n_2^\epsilon - \Delta n_2^\epsilon = n_2^\epsilon f_2(n_1^\epsilon, n_2^\epsilon). \qquad (3b) $$

The authors obtained that under some appropriate conditions, the proportion $p^\epsilon = \frac{n_1^\epsilon}{n_1^\epsilon+n_2^\epsilon}$ converges strongly in $L^2(0,T; L^2(\mathbb{R}^d))$, and weakly in $L^2(0,T; H^1(\mathbb{R}^d))$ to the solution $p^0$ of the scalar reaction-diffusion equation $\partial_t p^0 - \Delta p^0 = f(p^0)$ when $\epsilon \to 0$, where $f$ can be given explicitly from $f_1, f_2$. Now, in order to describe and study the migration phenomenon, we aim here at considering system (3) in a bounded domain $\Omega$ and introduce the boundary conditions to characterize the inflow and outflow of individuals as follows

$$ \frac{\partial n_1^\epsilon}{\partial \nu} = -D(n_1^\epsilon - n_1^{\text{ext}, \epsilon}), \quad (t, x) \in (0, T) \times \partial \Omega, \qquad (4a) $$

$$ \frac{\partial n_2^\epsilon}{\partial \nu} = -D(n_2^\epsilon - n_2^{\text{ext}, \epsilon}), \quad (t, x) \in (0, T) \times \partial \Omega, \qquad (4b) $$

where $n_1^{\text{ext}, \epsilon}$, $n_2^{\text{ext}, \epsilon}$ depend on $\epsilon$ but do not depend on time $t$ and position $x$. (4) models the tendency of the population to cross the boundary, with rates proportional to the difference between the surrounding density and the density just inside $\Omega$. Reusing the idea in [26], we prove in Section A that the proportion $p^\epsilon = \frac{n_1^\epsilon}{n_1^\epsilon+n_2^\epsilon}$ converges on any bounded time-domain to the solution of (1) when $\epsilon$ goes to zero. Hence, we can reduce the system (3), (4) to a simpler setting as in (1). The proof is based on a relative compactness argument that was also used in previous works about singular limits (e.g. [8,9,26]), but here, the use of the trace theorem is necessary to prove the limit on the boundary.

The outline of this work is the following. In the next section, we present the setting of the problem and the main results. In Section 3, we provide detailed proof of these results. Section 4 is devoted to an application to the biological control of mosquitoes. We also present numerical simulations to illustrate the theoretical results we obtained. Section A is devoted to proving the asymptotic limit of a 2-by-2 reaction-diffusion system when the reaction rate goes to infinity. Finally, we end this article with a conclusion and perspectives section.

# 2 Results on the steady-state solutions

## 2.1 Setting of the problem

In one-dimensional space, consider the system (1) in a bounded domain $\Omega = (-L, L) \subset \mathbb{R}$. Let $D > 0$, $p^{\text{ext}} \in (0, 1)$ be some constant and $p^{\text{init}}(x) \in [0, 1]$ for all $x \in (-L, L)$. The reaction term $f$ satisfies the following assumptions

**Assumption 2.1 (bistability).** Function $f: [0, 1] \to \mathbb{R}$ is of class $C^1([0, 1])$ and $f(0) = f(\theta) = f(1) = 0$ with $\theta \in (0, 1)$, $f(q) < 0$ for all $q \in (0, \theta)$, and $f(q) > 0$ for all $q \in (\theta, 1)$. Moreover, $\int_0^1 f(s)ds > 0$.

**Assumption 2.2 (convexity).** There exist $\alpha_1 \in (0, \theta)$ and $\alpha_2 \in (\theta, 1)$ such that $f'(\alpha_1) = f'(\alpha_2) = 0$, $f'(q) < 0$ for any $q \in [0, \alpha_1) \cup (\alpha_2, 1]$, and $f'(q) > 0$ for $q \in (\alpha_1, \alpha_2)$. Moreover, $f$ is convex on $(0, \alpha_1)$ and concave on $(\alpha_2, 1)$.

A function $f$ satisfying Assumptions 2.1 and 2.2 is illustrated in Figure 1a.

**Remark 2.1.**

1. Due to Assumption 2.1 and the fact that $p^{\text{ext}} \in (0, 1)$, $p^{\text{init}}(x) \in [0, 1]$ for any $x$, one has that 0 and 1 are respectively sub- and super-solution of problem (1). Since $f$ is Lipschitz continuous on $(0, 1)$ then by Theorem 4.1, Section 2.4 in [19], we obtain that problem (1) has a unique solution $p^0$ that is in $C^{1,2}((0, T] \times \Omega)$ with $0 \le p^0(t, x) \le 1$ for all $x \in (-L, L)$, $t > 0$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the functions f and F</figcaption></figure>

2. Again by Assumption 2.1, 0 and 1 are respectively sub- and super-solutions of (2). For fixed values of D, $p^{\text{ext}}$ and L, we use the same method as in [19] to obtain that there exists a $C^2$ solution of (2) with values in $[0, 1]$. However, Assumption 2.1 and 2.2 on f are not enough to conclude the uniqueness of the solution. In the following section, we prove that the stationary problem (2) may have multiple solutions and their existence depends on the values of the parameters.

The following proposition shows that solutions of system (2) always have at least one extreme value in $(-L, L)$.

**Proposition 2.1.** For any $p^{\text{ext}} \in (0, 1)$ and $p^{\text{ext}} \neq \theta$, system (2) does not have any non-constant monotone solution on the whole interval $(-L, L)$.

*Proof.* Assume that (2) admits an increasing solution $p$ on $(-L, L)$ (the case when $p$ is decreasing on $(-L, L)$ is analogous). Thus, we have $p'(x) \ge 0$ for all $x \in [-L, L]$ and $p(L) > p(-L)$. So thanks to the boundary condition of (2), one has

$$Dp^{\text{ext}} = p'(L) + Dp(L) \ge Dp(L) > Dp(-L) \ge -p'(-L) + Dp(-L) = Dp^{\text{ext}},$$

which is impossible. Therefore, we can deduce that the solutions of system (2) always admit at least one local extremum on the open interval $(-L, L)$. $\square$

To study system (2), we define function $F$ (see Figure 1b) as below

$$F(q) = \int_{0}^{q} f(s)ds, \qquad (5)$$

then $F'(q) = f(q)$ and $F(0) = 0$. From Assumption 2.1, $F$ reaches the minimal value at $q = \theta$ and the (locally) maximal values at $q = 0$ and $q = 1$. Since $\int_{0}^{1} f(s)ds > 0$, then $F(1) > F(0)$, it implies that $F(1) = \max_{[0,1]} F; F(\theta) = \min_{[0,1]} F$. Moreover, since $F(\theta) < F(0)$ and function $F$ is monotone in $(\theta, 1)$ ($F'(q) = f(q) > 0$ for any $q \in (\theta, 1)$). Thus, there exists a unique value $\beta \in (\theta, 1)$ such that

$$F(\beta) = F(0) = 0. \qquad (6)$$

The main results of the present work concern the existence and stability of steady-state solutions of (1), i.e. solutions of (2).

## 2.2 Existence of steady-state solutions

In our result, we first focus on two types of steady-state solutions defined as follows
<figure><img src="image_3.png" /><figcaption>Figure 2: Sketch of the symmetric steady-state solutions p</figcaption></figure>

**Definition 2.1.** Consider a steady-state solution $p(x)$,

p is called a symmetric-decreasing (SD) solution when p is symmetric on $(-L, L)$ with values in $[0, 1]$, decreasing on $(0, L)$ and $p'(0) = 0$ (see Figure 2a).

Similarly, p is called a symmetric-increasing (SI) solution when p is symmetric on $(-L, L)$ with values in $[0, 1]$, increasing on $(0, L)$ and $p'(0) = 0$ (see Figure 2b).

Any solution which is either (SD) or (SI) is called a symmetric-monotone (SM) solution.

The following theorems present the main result of the existence of (SM) solutions depending upon the parameters. For each value of $p^{\text{ext}} \in (0, 1)$ and $D > 0$, we find the critical values of L such that (2) admits solutions.

**Theorem 2.1.** In a bounded domain $\Omega = (-L, L) \subset \mathbb{R}$, consider the stationary problem (2). Assume that the reaction term f satisfies Assumptions 2.1 and 2.2. Then, there exist two functions

$$ M_d, M_i: (0,1) \times (0, +\infty) \rightarrow [0, +\infty], \\ (p^{\text{ext}}, D) \mapsto M_d(p^{\text{ext}}, D), M_i(p^{\text{ext}}, D), \qquad (7) $$

such that for any $p^{\text{ext}} \in (0, 1)$, $D > 0$, problem (2) admits at least one (SD) solution (resp., (SI) solution) if and only if $L > M_d(p^{\text{ext}}, D)$ (resp., $L > M_i(p^{\text{ext}}, D)$) and the values of these solutions are in $[p^{\text{ext}}, 1]$ (resp., $[0, p^{\text{ext}}]$). More precisely,

1. If $0 < p^{\text{ext}} < \theta$, then for any $D > 0$, $M_i(p^{\text{ext}}, D) = 0$ and $M_d(p^{\text{ext}}, D) \in (0, +\infty)$.

Moreover, if $p^{\text{ext}} \le \alpha_1$, the (SI) solution is unique.

2. If $\theta < p^{\text{ext}} < 1$, then for any $D > 0$, $M_d(p^{\text{ext}}, D) = 0$. If $\alpha_2 \le p^{\text{ext}}$, the (SD) solution is unique.

Moreover, consider $\beta$ as in (6),

* if $p^{\text{ext}} \le \beta$, then $M_i(p^{\text{ext}}, D) \in (0, +\infty)$ for any $D > 0$;

* if $p^{\text{ext}} > \beta$, then there exists a constant $D_* > 0$ such that $M_i(p^{\text{ext}}, D) \in (0, +\infty)$ for any $D < D_*$, and $M_i(p^{\text{ext}}, D) = +\infty$ for $D \ge D_*$.

3. If $p^{\text{ext}} = \theta$, then $M_d(\theta, D) = M_i(\theta, D) = 0$. Moreover, there exists a constant solution $p \equiv \theta$.

In the statement of the above result, $M_i = 0$ means that for any $L > 0$, (2) always admits (SI) solutions. $M_i = +\infty$ means that there is no (SI) solution even when L is large. The same interpretation applies for $M_d$.

Besides, problem (2) can also admit solutions that are neither (SD) nor (SI). The following theorem provides an existence result for those solutions.

**Theorem 2.2.** In a bounded domain $\Omega = (-L, L) \subset \mathbb{R}$, consider the stationary problem (2). Assume that the reaction term f satisfies Assumption 2.1 and 2.2. Then, there exists a function

$$ M_*: (0,1) \times (0, +\infty) \rightarrow [0, +\infty], \\ (p^{\text{ext}}, D) \mapsto M_*(p^{\text{ext}}, D), \qquad (8) $$
<table><caption>Table 1: The existence of steady-state solutions corresponding to values of parameters.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="3">0 &lt; p<sup>ext</sup> &lt; &theta;, D &gt; 0</th></tr><tr><th>0 = M<sub>i</sub> &lt; L &lt; M<sub>d</sub></th><th>M<sub>d</sub> &le; L &lt; M<sub>*</sub></th><th>M<sub>*</sub> &le; L &lt; +&infin;</th></tr></thead><tbody><tr><th>Types of solutions</th><td>(SI)</td><td>(SI), (SD)</td><td>(SI), (SD), non-(SM)</td></tr><tr><th rowspan="2">Parameters</th><td colspan="3">p<sup>ext</sup> = &theta;, D &gt; 0</td></tr><tr><td colspan="3">0 = M<sub>d</sub> = M<sub>i</sub> = M<sub>*</sub> &lt; L</td></tr><tr><th>Types of solutions</th><td colspan="3">(SD), (SI), non-(SM)</td></tr><tr><th rowspan="2">Parameters</th><td colspan="3">&theta; &lt; p<sup>ext</sup> &le; &beta;, D &gt; 0 or p<sup>ext</sup> &gt; &beta;, 0 &lt; D &lt; D<sub>*</sub></td></tr><tr><th>0 = M<sub>d</sub> &lt; L &lt; M<sub>i</sub></th><th>M<sub>i</sub> &le; L &lt; M<sub>*</sub></th><th>M<sub>*</sub> &le; L &lt; +&infin;</th></tr><tr><th>Types of solutions</th><td>(SD)</td><td>(SI), (SD)</td><td>(SI), (SD), non-(SM)</td></tr><tr><th rowspan="2">Parameters</th><td colspan="3">&beta; &lt; p<sup>ext</sup> &lt; 1, D &gt; D<sub>*</sub></td></tr><tr><td colspan="3">0 = M<sub>d</sub> &lt; L &lt; M<sub>i</sub> = M<sub>*</sub> = +&infin;</td></tr><tr><th>Types of solutions</th><td colspan="3">(SD)</td></tr></tbody></table>

such that for any $p^{\text{ext}} \in (0,1)$, $D > 0$, problem (2) admits at least one solution which is not (SM) if and only if $L \ge M_*(p^{\text{ext}}, D)$. Moreover,

• If $p^{\text{ext}} \le \beta$, then for any $D > 0$, one has

$$0 < M_i(p^{\text{ext}}, D) + M_d(p^{\text{ext}}, D) < M_*(p^{\text{ext}}, D) < +\infty. \quad (9)$$

• If $p^{\text{ext}} > \beta$, then for any $D < D_*$, one has $0 < M_i(p^{\text{ext}}, D) < M_*(p^{\text{ext}}, D) < +\infty$. Otherwise, for $D \ge D_*$, $M_*(p^{\text{ext}}, D) = +\infty$. Here, $D_*$ was defined in Theorem 2.1.

The construction of $M_i$, $M_d$, $M_*$ will be done in the proof in Section 3. The idea of the proof is based on a careful study of the phase portrait of (2). To make the results more reader-friendly, we present the types of steady-state solutions corresponding to different parameters in Table 1.

In the next section, we present a result about the stability and instability of steady-state solutions of (2).

## 2.3 Stability of steady-state solutions

The definition of stability and instability used in the present work comes from Lyapunov stability

**Definition 2.2.** A steady-state solution $p(x)$ of (1) is called stable if for any constant $\epsilon > 0$, there exists a constant $\delta > 0$ such that when $||p^{\text{init}} - p||_{\infty} < \delta$, one has

$$||p^0(t, \cdot) - p||_{\infty} < \epsilon, \quad \text{for all } t > 0 \qquad (10)$$

where $p^0(t, x)$ is the unique solution of (1). If, in addition,

$$\lim_{t \to \infty} ||p^0(t, \cdot) - p||_{\infty} = 0, \qquad (11)$$

then $p$ is called asymptotically stable. The steady-state solution $p$ is called unstable if it is not stable.

The following theorem provides sufficient conditions for the stability of steady-state solutions given in Section 2.2.

**Theorem 2.3.** In the bounded domain $\Omega = (-L, L) \subset \mathbb{R}$, consider the problem (1) with the reaction term satisfying Assumptions 2.1 and 2.2. There exists a constant $\lambda_1 \in (0, \frac{\pi^2}{4L^2})$ such that for any steady-state solution $p$ of (1),

• If $f'(p(x)) > \lambda_1$ for all $x \in (-L, L)$, then $p$ is unstable.

• If $f'(p(x)) < \lambda_1$ for all $x \in (-L, L)$, then $p$ is asymptotically stable.

More precisely, $\lambda_1$ is the principal eigenvalue of the linear problem (12)

$$- \phi''(x) = \lambda \phi(x) \quad x \in (-L, L), \qquad (12a)$$

$$\phi'(L) = -D\phi(L), \qquad (12b)$$

$$\phi'(-L) = D\phi(-L), \qquad (12c)$$
where $\lambda$ is an eigenvalue with associated eigenfunction $\phi$. It may be proved that its value is the smallest positive solution of equation $\sqrt{\lambda} \tan(L\sqrt{\lambda}) = D$ (see more details in Section 3).

Note that we cannot apply the first statement if $\sup_{q \in (0,1)} f'(q) \le \lambda_1$. However, due to the fact that $\lambda_1 \in (0, \frac{\pi^2}{4L^2})$, when $L$ gets larger, the value of $\lambda_1$ gets closer to zero and the inequality in the first statement becomes valid.

**Remark 2.2.** By Assumption 2.2, $f'(q) \le 0 < \lambda_1$ for all $q \in [0, \alpha_1] \cup [\alpha_2, 1]$, we can deduce that the steady-state solutions with values smaller than $\alpha_1$ or larger than $\alpha_2$ are asymptotically stable.

As a consequence of Theorems 2.1, 2.2, and 2.3, the following important result provides complete information about the existence and stability of steady-state solutions in some special cases.

**Corollary 2.1.** In the bounded domain $\Omega = (-L, L) \subset \mathbb{R}$, consider the problem (1) with the reaction term satisfying Assumption 2.1 and 2.2. Then for any $D > 0$, we have

* If $p^{\text{ext}} \le \alpha_1$, for any $L > 0$, there exists exactly one (SI) steady-state solution $p$ and it is asymptotically stable. Moreover, if $L < M_d(p^{\text{ext}}, D)$, then $p$ is the unique steady-state solution of (1).

* If $p^{\text{ext}} \ge \alpha_2$, for any $L > 0$, there exists exactly one (SD) steady-state solution $p$ and it is asymptotically stable. Moreover, if $L < M_i(p^{\text{ext}}, D)$, then $p$ is the unique steady-state solution of (1).

**Remark 2.3.** This corollary gives us a comprehensive view of the long-time behavior of solutions of (1) when the size $L$ of the domain is small. In this case, the unique steady-state solution $p$ is symmetric, monotone on each half of $\Omega$, and asymptotically stable. Its values will be close to 0 if $p^{\text{ext}}$ is small and close to 1 if $p^{\text{ext}}$ is large. We discuss an essential application of this result in Section 4.

# 3 Proof of the theorems

## 3.1 Proof of existence

In this section, we use phase-plane analysis to prove the existence of both (SM) and non-(SM) steady-state solutions depending on the parameters. The studies of (SD) and (SI) solutions will be presented respectively in Sections 3.1.1 and 3.1.2. Then, using these results, we prove Theorem 2.1. The proof of Theorem 2.2 will be presented after that using the same technique.

First, we introduce the following function

$$E(p, p') = \frac{(p')^2}{2} + F(p). \qquad (13)$$

Since $\frac{d}{dx}E(p, p') = p'(p'' + f(p)) = 0$, then $E(p, p')$ is constant along the orbit of (2). From Proposition 2.1, we can deduce that there exists an $x_0 \in (-L, L)$ such that $p'(x_0) = 0$, thus one has

$$E(p(x_0), 0) = E(p(x), p'(x)), \qquad (14)$$

for all $x \in (-L, L)$. Therefore, the relation between $p'$ and $p$ is as below

$$p' = \pm\sqrt{2F(p(x_0)) - 2F(p)}. \qquad (15)$$

According to this relation, one has a phase plane as in Figure 3a, in which the curves illustrate the relation between $p'(x)$ and $p(x)$ in (15) with respect to different values of $p(x_0)$. We can see that some curves do not end on the axis $p = 0$ but wrap around the point $(\theta, 0)$. This is due to the fact that for any $p_1 \in [\theta, \beta]$, there exists a value $p_2 \in [0, \theta]$ such that $F(p_1) = F(p_2)$. Thus, if the curve passes through the point $(p_1, 0)$, it will also pass through the point $(p_2, 0)$ on the axis $p' = 0$. Moreover, those curves only exist if their intersection with the axis $p' = 0$ has p-coordinate less than or equal to $\beta$. Besides, the two straight lines show the relation between $p'$ and $p$ at the boundary points. Solutions of (2) correspond to those orbits that connect the intersection of the curves with the line $p' = D(p - p^{\text{ext}})$ to the intersection of the curves with the line $p' = -D(p - p^{\text{ext}})$.

In the phase plane in Figure 3a, orbit $T_1$ describes an (SD) solution, while orbit $T_2$ corresponds to a (SI) solution. On the other hand, the solid curve $T_3$ shows the orbit of a steady-state solution that is not symmetric-monotone.
<figure><img src="image_4.png" /><figcaption>Figure 3: Phase portraits of (2): straight lines illustrate the boundary conditions, solid curves show relations between $p'$ and $p$. Figure (a): curves $T_1$, $T_2$, and $T_3$ correspond to orbits of (SD), (SI), and non-(SM) solutions, respectively. Figure (b): curve $T_4$ corresponds to an orbit of a non-(SM) solution.</figcaption></figure>

**Remark 3.1.** (Graphical interpretation of $D_*$) The (SI) solutions (see Figure 2b) have orbit as $T_2$ in Figure 3a. This type of orbits only exists when the lines $p = \pm D(p - p^{\text{ext}})$ intersect the curves wrapping around the point $(\theta, 0)$. In the case when $p^{\text{ext}} > \beta$, the constant $D_* > 0$ in Theorem 2.1 is the slope of the tangent line to the curve passing through $(\beta, 0)$ as in Figure 3b. Hence, if $D > D_*$, there exists no (SI) solution. We construct explicitly the value of $D_*$ in Proposition 3.2 below.

Next, we establish some relations between the solution $p$ and the parameters based on the phase portrait above. For any $x > x_0$, if $p$ is monotone on $(x_0, x)$, we can invert $x \mapsto p(x)$ into function $p \mapsto X(p)$. We obtain $X'(p) = \frac{\pm 1}{\sqrt{2F(p(x_0))-2F(p)}}$. By integrating this equation, we obtain that

$$x - x_0 = \int_{p(x_0)}^{p(x)} \frac{(-1)^k ds}{\sqrt{2F(p(x_0)) - 2F(s)}}, \quad (16)$$

where $k = 1$ if $p$ is decreasing and $k = 2$ if $p$ is increasing on $(x_0, x)$. We can obtain the analogous formula for $x < x_0$.

First, we focus on symmetric-monotone (SM) solutions for which $p'(0) = 0$, then we analyze the integral in (16) with $x = L, x_0 = 0$. For any $p^{\text{ext}} \in (0, 1)$, using (15), we have

$$F(p(0)) = F(p(L)) + \frac{1}{2} D^2 (p(L) - p^{\text{ext}})^2 = G(p(L)), \quad (17)$$

for $F$ defined in (5) and

$$G(q) := F(q) + \frac{1}{2} D^2 (q - p^{\text{ext}})^2, \quad (18)$$

and from (16) with $x = L, x_0 = 0$, we have

$$L = \int_{p(0)}^{p(L)} \frac{(-1)^k ds}{\sqrt{2F(p(0)) - 2F(s)}}, \quad (19)$$

where $k = 1$ if $p$ is decreasing on $(0, L)$, $k = 2$ if $p$ is increasing on $(0, L)$.

Thus, the (SM) solution of (2) exists if there exist values $p(L)$ and $p(0)$ that satisfy (17) and (19). When such values exist, we can assess the value of $p(x)$ for any $x$ in $(-L, L)$ using (16).

Before proving the existence of such values of $p(0)$ and $p(L)$, we establish some useful properties of the function $G$ defined in (18). It is continuous in $[0, 1]$ and $G(q) \ge F(q)$ for all $q \in [0, 1]$. Moreover, the following lemma shows that $G$ has a unique minimum point.

**Lemma 3.1.** For any $p^{\text{ext}} \in (0, 1)$, there exists a unique value $\bar{q} \in (0, 1)$ such that $G'(\bar{q}) = 0$, $G'(q) < 0$ for all $q \in [0, \bar{q})$ and $G'(q) > 0$ for all $q \in (\bar{q}, 1]$. Particularly, $G(\bar{q}) = \min_{[0,1]} G$.
Proof. We have $G'(q) = f(q) + D^2(q - p^{\text{ext}})$. We consider the following cases.

**Case 1:** When $p^{\text{ext}} = \theta$, we have $G'(p^{\text{ext}}) = G'(\theta) = f(\theta) = 0$, $G'(q) < 0$ for all $q \in (0, \theta)$ and $G'(q) > 0$ for all $q \in (\theta, 1)$. Thus $\bar{q} = \theta = p^{\text{ext}}$.

**Case 2:** When $p^{\text{ext}} < \theta$, we have $G'(q) < 0$ for all $q \in [0, p^{\text{ext}}]$ and $G'(q) > 0$ for all $q \in [\theta, 1]$. So there exists at least one value $\bar{q} \in (p^{\text{ext}}, \theta)$ such that $G'(\bar{q}) = 0$.

For any $\bar{q} \in (p^{\text{ext}}, \theta)$ such that $G'(\bar{q}) = 0$, we have $f(\bar{q}) + D^2(\bar{q} - p^{\text{ext}}) = 0$ so that $D^2 = -\frac{f(\bar{q})}{\bar{q}-p^{\text{ext}}}$. We can prove that $G''(\bar{q})$ is strictly positive. Indeed, from Assumption 2.2 we have that $\alpha_1$ is the unique value in $(0, \theta)$ such that $f'(\alpha_1) = 0$, thus $f(\alpha_1) = \min_{[0,\theta]} f < 0$.

If $\alpha_1 \le \bar{q} < \theta$ then $f'(\bar{q}) \ge 0$. One has $G''(\bar{q}) = f'(\bar{q}) + D^2 > 0$.

If $p^{\text{ext}} < \bar{q} < \alpha_1$, due to the fact that $f$ is convex in $(0, \alpha_1)$ one has $f'(\bar{q}) \ge \frac{f(\bar{q})-f(p^{\text{ext}})}{\bar{q}-p^{\text{ext}}}$. Since $f(p^{\text{ext}}) < 0$, one has $G''(\bar{q}) = f'(\bar{q}) + D^2 = f'(\bar{q}) - \frac{f(\bar{q})}{\bar{q}-p^{\text{ext}}} > f'(\bar{q}) + \frac{f(p^{\text{ext}})-f(\bar{q})}{\bar{q}-p^{\text{ext}}} \ge 0$. One can deduce that $\bar{q}$ is the unique value in $(0, 1)$ such that $G'(\bar{q}) = 0$ and $G(\bar{q}) = \min_{[0,1]} G$, so it satisfies Lemma 3.1.

**Case 3:** When $p^{\text{ext}} > \theta$, the proof is analogous to case 2 but using the concavity of $f$ in $(\alpha_2, 1)$. We obtain that there exists a unique value $\bar{q}$ in $(\theta, p^{\text{ext}})$ that satisfies Lemma 3.1. $\square$

When $p^{\text{ext}} = \theta$, it is easy to check that $p \equiv \theta$ is a solution of (2). We now analyze two types of (SM) solutions (see Figure 2) in the following parts.

### 3.1.1 Existence of (SD) solutions

In this part, the solution $p$ we study is symmetric on $(-L, L)$ and decreasing on $(0, L)$ (see Figure 2a). So $p(L) < p(x) < p(0)$ for any $x \in (0, L)$. But from (15), we have that $F(p(x)) \le F(p(0))$, so $F'(p(0)) \ge 0$. It implies that $p(0) \in [\theta, 1]$. Next, we use two steps to study the existence of (SD) solutions:

**Step 1: Rewriting as a non-linear equation on $p(L)$**

For any $q \in (\theta, 1)$, we have $F'(q) = f(q) > 0$ so $F|_{(\theta,1)} : (\theta, 1) \to (F(\theta), F(1))$ is invertible. Define $F_1^{-1} := (F|_{(\theta,1)})^{-1} : (F(\theta), F(1)) \to (\theta, 1)$, and $F_1^{-1}(F(\theta)) = \theta, F_1^{-1}(F(1)) = 1$. Then, $F_1^{-1}$ is continuous in $[F(\theta), F(1)]$. For any $y \in (F(\theta), F(1))$, one has $(F_1^{-1})'(y) = \frac{1}{F'(F_1^{-1}(y))} = \frac{1}{f(F_1^{-1}(y))} > 0$, so $F_1^{-1}$ is an increasing function in $(F(\theta), F(1))$. From (17) and (19), since $p$ is decreasing in $(0, L)$, we have $L = \int_{p(L)}^{p(0)} \frac{ds}{\sqrt{2G(p(L)) - 2F(s)}}$. Denote

$$ \mathcal{F}_1(q) := \int_q^{F_1^{-1}(G(q))} \frac{ds}{\sqrt{2G(q) - 2F(s)}}. \qquad (20) $$

Hence, a (SD) solution $p$ of system (2) has $p(0) = F_1^{-1}(G(p(L)))$, and $p(L)$ satisfies

$$ L = \mathcal{F}_1(p(L)). \qquad (21) $$

Moreover, one has $p'(x) \le 0$ for all $x \in (0, L)$ thus $-D(p(L) - p^{\text{ext}}) = p'(L) \le 0$. One can deduce that

$$ p(L) \ge p^{\text{ext}}. \qquad (22) $$

**Step 2: Solving (21) in $[p^{\text{ext}}, 1]$**

The following proposition states the existence of a solution of (21).

**Proposition 3.1.** For any $D > 0$, $p^{\text{ext}} \in (0, 1)$, we have

1. If $0 < p^{\text{ext}} < \theta$, then there exists a constant $M_1 > 0$ such that equation (21) has at least one solution $p(L) \ge p^{\text{ext}}$ if and only if $L \ge M_1$.

2. If $\theta \le p^{\text{ext}} < 1$, then equation (21) admits at least one solution $p(L) \ge p^{\text{ext}}$ for all $L > 0$. If $p^{\text{ext}} \ge \alpha_2$, then this solution is unique.
Proof. Since $F_1^{-1}$ is only defined in $[F(\theta), F(1)]$, we need to find $p(L) \in [p^{\text{ext}}, 1]$ such that $G(p(L)) \in [F(\theta), F(1)]$.

For all $q \in (0, 1)$, we have $G(q) \ge F(q) \ge F(\theta)$ and from Lemma 3.1, there exists a value $\bar{q} \in (0, 1)$ such that $\min_{[0,1]} G = G(\bar{q}) \le G(p^{\text{ext}}) = F(p^{\text{ext}}) < \max_{[0,1]} F = F(1)$. Moreover, one has $G(1) > F(1)$, thus there exists a value $p^* \in (p^{\text{ext}}, 1)$ such that $G(p^*) = F(1)$. Then, for all $q \in [p^{\text{ext}}, p^*]$, $G(q) \in [F(\theta), F(1)]$ and we will find $p(L)$ in $[p^{\text{ext}}, p^*]$. Since $F_1^{-1}$ increases in $(F(\theta), F(1))$, then $p(0) = F_1^{-1}(G(p(L))) \ge F_1^{-1}(F(p(L))) \ge p(L)$.

Function $\mathcal{F}_1$ in (20) is well-defined and continuous in $[p^{\text{ext}}, p^*)$, $\mathcal{F} \ge 0$ in $[p^{\text{ext}}, p^*)$. Moreover, since $F'(1) = 0$, one has $\lim_{p \to p^*} \mathcal{F}_1(p) = \int_{p^*}^{1} \frac{ds}{\sqrt{2F(1) - 2F(s)}} = +\infty$.

**Case 1:** If $0 < p^{\text{ext}} < \theta$, we will prove that $\mathcal{F}_1$ is strictly positive in $[p^{\text{ext}}, p^*)$. Indeed, for any $y \in [0, 1]$, if $y < \theta$, by the definition of $F_1^{-1}$, we have $F_1^{-1}(G(y)) \in [\theta, 1]$ so $F_1^{-1}(G(y)) > y$. If $y \ge \theta > p^{\text{ext}}$ then $G(y) = F(y) + \frac{1}{2}D^2(y - p^{\text{ext}})^2 > F(y)$ so again $F_1^{-1}(G(y)) > y$. Hence $\mathcal{F}_1(y) > 0$ for all $y \in [p^{\text{ext}}, p^*)$. We have $\mathcal{F}_1(p) \to +\infty$ when $p \to p^*$, so there exists $p \in [p^{\text{ext}}, p^*)$ such that $M_1 := \mathcal{F}_1(p) = \min_{[p^{\text{ext}}, p^*]} \mathcal{F}_1 > 0$, and system (21) admits at least one solution if and only if $L \ge M_1$.

**Case 2:** If $\theta \le p^{\text{ext}} < 1$, one has $G(p^{\text{ext}}) = F(p^{\text{ext}})$, then $F_1^{-1}(G(p^{\text{ext}})) = p^{\text{ext}}$ so $\mathcal{F}_1(p^{\text{ext}}) = 0$. On the other hand, $\mathcal{F}_1(p) \to +\infty$ when $p \to p^*$. Thus, for any $L > 0$, there always exists at least one value $p(L) \in (p^{\text{ext}}, p^*)$ such that $\mathcal{F}_1(p(L)) = L$.

*Proof of uniqueness:* When $p^{\text{ext}} \ge \alpha_2$, we can prove that $\mathcal{F}'_1 > 0$ on $(p^{\text{ext}}, p^*)$. Indeed, denoting $\gamma(q) = F_1^{-1}(G(q))$, and changing the variable from $s$ to $t$ such that $s = t\gamma(q) + (1-t)q$, one has

$$ \mathcal{F}_1(q) = \int_0^1 \frac{[\gamma(q) - q]dt}{\sqrt{2F(\gamma(q)) - 2F(t\gamma(q) + (1-t)q)}} $$

To simplify, denote $s(q) = t\gamma(q) + (1-t)q$. For any $t \in (0, 1)$, one has $q < s(q) < \gamma(q)$. Let us define $\Delta F = F(\gamma(q)) - F(s(q))$, then one has

$$ \begin{aligned} \sqrt{2}\mathcal{F}'_1(q) &= \int_0^1 (\gamma'(q) - 1)(\Delta F)^{-1/2} dt - \frac{1}{2}\int_0^1 (\Delta F)^{-3/2}(\gamma(q) - q) \frac{d\Delta F}{dq} dt \\ &= \int_0^1 (\Delta F)^{-3/2} \left[ (\gamma'(q) - 1)\Delta F - \frac{1}{2}(\gamma(q) - q)(f(\gamma(q))\gamma'(q) - f(s(q))s'(q)) \right]. \end{aligned} $$

Let $P$ be the formula in the brackets, then

$$ \begin{aligned} P &= (\gamma' - 1)\Delta F - \frac{1}{2}(\gamma - q)[f(\gamma)\gamma' - f(s)(t\gamma' + 1 - t)] \\ &= (\gamma' - 1)[\Delta F - \frac{1}{2}(\gamma - q)f(\gamma) + \frac{1}{2}(s - q)f(s)] - \frac{1}{2}(\gamma - q)(f(\gamma) - f(s)), \end{aligned} $$

Define $\psi(y) := F(y) - \frac{1}{2}f(y)(y - q)$ for any $y \in [q, \gamma(q)]$, then one has $\psi'(y) = \frac{1}{2}[f(y) - f'(y)(y - q)] \ge \frac{f(q)}{2} > 0$ since $y \ge q > p^{\text{ext}} \ge \alpha_2$ and $f$ is concave in $(\alpha_2, 1)$, $f(q) > 0$. Moreover, $f$ is decreasing on $(\alpha_2, 1)$ so $0 < f(\gamma(q)) < f(s(q)) < f(q)$, and $\gamma'(q) = \frac{G'(q)}{f(F_1^{-1}(G(q)))} = \frac{f(q) + D^2(q - p^{\text{ext}})}{f(\gamma(q))} > 1$. Hence, we can deduce that $P = (\gamma' - 1)(\psi(\gamma) - \psi(s)) - \frac{1}{2}(\gamma - q)(f(\gamma) - f(s)) > 0$ for any $t \in (0, 1)$. This proves that function $\mathcal{F}_1$ is increasing on $(p^{\text{ext}}, p^*)$, so the solution of equation (21) is unique. $\square$

### 3.1.2 Existence of (SI) solutions

In this case, the technique we use to prove the existence of (SI) solutions is analogous to (SD) solutions except in the case when $p^{\text{ext}} > \beta$ (case 3 below). Since the proof is not straightforward, it is worth to re-establish this technique for (SI) solutions in two following steps:

#### Step 1: Rewriting as a non-linear equation on $p(L)$

Since now $p$ is symmetric on $(-L, L)$ and increasing in $(0, L)$ (see Figure 2b), then $p(0) < p(x) < p(L)$ for any $x \in (0, L)$. But from (15), we have that $F(p(x)) \le F(p(0))$, so $F'(p(0)) \le 0$. This implies that $p(0) \in [0, \theta]$.

For any $q \in (0, \theta)$, we have $F'(q) = f(q) < 0$ so $F|_{(0,\theta)}: (0, \theta) \to (F(\theta), F(0))$ is invertible. Define $F_2^{-1} := (F|_{(0,\theta)})^{-1}: (F(\theta), F(0)) \to (0, \theta)$, $F_2^{-1}(F(\theta)) = \theta$, $F_2^{-1}(F(0)) = 0$, and $F_2^{-1}$ is continuous in $[F(\theta), F(0)]$. For any $y \in (F(\theta), F(0))$, $(F_2^{-1})'(y) = \frac{1}{F'(F_2^{-1}(y))} = \frac{1}{f(F_2^{-1}(y))} < 0$, so $F_2^{-1}$ is a decreasing
function in $(F(\theta), F(0))$. From (17) and (19), we have $L = \int_{p(0)}^{p(L)} \frac{ds}{\sqrt{2G(p(L))-2F(s)}}$. Denote

$$ \mathcal{F}_2(q) := \int_{F_2^{-1}(G(q))}^{q} \frac{ds}{\sqrt{2G(q) - 2F(s)}}. \qquad (23) $$

Hence, a (SI) solution of system (2) has $p(0) = F_2^{-1}(G(p(L)))$, and $p(L)$ satisfies

$$ L = \mathcal{F}_2(p(L)), \qquad (24) $$

and in this case, one needs to find $p(L)$ in $[0, p^{\text{ext}}]$.

### Step 2: Solving of (24) in $[0, p^{\text{ext}}]$

The following proposition states the existence of a solution of (24).

**Proposition 3.2.** For any $p^{\text{ext}} \in (0, 1)$, considering the value $\beta$ as in (6), we have:

1. If $0 < p^{\text{ext}} \le \theta$, then equation (24) admits at least one solution $p$ with $p(L) \le p^{\text{ext}}$ for all $L > 0, D > 0$. If $p^{\text{ext}} \le \alpha_1$, this solution is unique.

2. If $\theta < p^{\text{ext}} \le \beta$, then for all $D > 0$, there exists a constant $M_2 > 0$ such that equation (24) has at least one solution $p$ with $p(L) \le p^{\text{ext}}$ if and only if $L \ge M_2$.

3. If $\beta < p^{\text{ext}} < 1$, then there exists a constant $D_* > 0$ such that when $D \ge D_*$, equation (24) has no solution. Otherwise, there exists a constant $M_3 > 0$ such that equation (24) has at least one solution $p$ with $p(L) \le p^{\text{ext}}$ if and only if $L \ge M_3$.

*Proof.* As we assume that $F(0) < F(1)$ and $F(\theta) < F(0)$ then, due to the continuity of $F$, one can deduce that there exists a value $\beta \in (\theta, 1)$ such that $F(\beta) = F(0) = 0$.

Since $F_2^{-1}$ is only defined in $[F(\theta), F(0)]$, we need to find $p(L) \in [0, p^{\text{ext}}]$ such that $G(p(L)) \in [F(\theta), F(0)]$. For all $q \in (0, 1)$, we have $G(q) \ge F(q) \ge F(\theta)$, thus equation (24) has solutions if and only if $\min_{[0,1]} G < F(0)$. Even when $\min_{[0,1]} G = G(\bar{q}) = F(0)$, $\mathcal{F}_2$ is still not defined in $[0, 1]$ since $\mathcal{F}_2(\bar{q}) = +\infty$.

One has the following cases:

**Case 1:** $0 < p^{\text{ext}} \le \theta$:

We have $\min_{[0,1]} G = G(\bar{q}) \le G(p^{\text{ext}}) = F(p^{\text{ext}}) < \max_{[0,\theta]} F = F(0)$, and $G(0) > F(0)$ so there is a value $p_* \in (0, p^{\text{ext}})$ such that $G(p_*) = F(0)$. Moreover $F'(0) = 0$, then $\lim_{p \to p^*} \mathcal{F}_2(p) = +\infty$. Thus, function $\mathcal{F}_2$ is only well-defined and continuous in $(p_*, p^{\text{ext}}]$.

When $0 < p^{\text{ext}} \le \theta$, $F_2^{-1}(G(p^{\text{ext}})) = F_2^{-1}(F(p^{\text{ext}})) = p^{\text{ext}}$ so $\mathcal{F}_2(p^{\text{ext}}) = 0$. We can deduce that for any $L > 0$, there always exists at least one value $p(L) \in (p_*, p^{\text{ext}})$ such that $\mathcal{F}_2(p(L)) = L$. When $p^{\text{ext}} \le \alpha_1$, arguing analogously to the second case of Proposition 3.1, one has $\mathcal{F}'_2 < 0$ on $(p_*, p^{\text{ext}})$, thus the solution is unique.

**Case 2:** $\theta < p^{\text{ext}} \le \beta$:

Since $F$ increases on $(\theta, 1)$, then $\min_{[0,1]} G = G(\bar{q}) < G(p^{\text{ext}}) = F(p^{\text{ext}}) \le F(\beta) = F(0)$. Analogously to the previous case, $\mathcal{F}_2$ is well-defined and continuous in $(p_*, p^{\text{ext}}]$, $\lim_{p \to p^*} \mathcal{F}_2(p) = +\infty$, and $\mathcal{F}_2$ is strictly positive in $(p_*, p^{\text{ext}}]$. Therefore, there exists $p \in (p_*, p^{\text{ext}}]$ such that

$$ M_2 := \mathcal{F}_2(p) = \min_{[p_*, p^{\text{ext}}]} \mathcal{F}_2 > 0, \qquad (25) $$

and system (24) admits as least one solution if and only if $L \ge M_2$.

**Case 3:** $\beta < p^{\text{ext}} < 1$:

Consider the function $H(q) = F(q) + \frac{1}{2}f(q)(p^{\text{ext}} - q)$ defined in an interval $[\theta, p^{\text{ext}}]$. For any $\theta < q < p^{\text{ext}}$, one can prove that $H'(q) \ge 0$.

Indeed, if $q \le \alpha_2$, then $f'(q) \ge 0$, and $f(q) > 0$. One has $H'(q) = \frac{1}{2}f(q) + \frac{1}{2}f'(q)(p^{\text{ext}} - q) > 0$. If $q > \alpha_2$, from Assumption 2.2, the function $f$ is concave in $(\alpha_2, 1)$, and hence $f'(q)(p^{\text{ext}} - q) \ge f(p^{\text{ext}}) - f(q)$. Thus,

$$ H'(q) = \frac{1}{2}(p^{\text{ext}} - q) \left( f'(q) + \frac{f(q)}{p^{\text{ext}} - q} \right) > \frac{1}{2}(p^{\text{ext}} - q) \left( f'(q) + \frac{f(q) - f(p^{\text{ext}})}{p^{\text{ext}} - q} \right) \ge 0. $$
Therefore, function H increases in $(\theta, p^{\text{ext}})$. Moreover $H(\theta) = F(\theta) < F(0)$ and $H(p^{\text{ext}}) = F(p^{\text{ext}}) > F(\beta) = F(0)$, and so there exists a unique value $\bar{p}_* \in (\theta, p^{\text{ext}})$ such that $H(\bar{p}_*) = F(0)$. Take $D_* > 0$ such that $D_*^2 = \frac{f(\bar{p}_*)}{p^{\text{ext}}-\bar{p}_*}$. Then, for any $D > 0$, from Lemma 3.1, there is a unique value $\bar{q} \in (\theta, p^{\text{ext}})$ such that $G'(\bar{q}) = 0$, $G(\bar{q}) = \min_{[0,1]} G$, and $D^2 = \frac{f(\bar{q})}{p^{\text{ext}}-\bar{q}}$. If $D < D_*$, then $\frac{f(\bar{q})}{p^{\text{ext}}-\bar{q}} < \frac{f(\bar{p}_*)}{p^{\text{ext}}-\bar{p}_*}$.

Let $h(q) = \frac{f(q)}{p^{\text{ext}}-q}$, then $h'(q) = \frac{1}{p^{\text{ext}}-q} \left( f'(q) + \frac{f(q)}{p^{\text{ext}}-q} \right) > 0$ for $q \in (\theta, p^{\text{ext}})$. So function h is increasing in $(\theta, p^{\text{ext}})$, and we can deduce that $\bar{q} < \bar{p}_*$. Hence, $\min_{[0,1]} G = G(\bar{q}) = F(\bar{q}) + \frac{1}{2} D^2 (p^{\text{ext}} - \bar{q})^2 = F(\bar{q}) + \frac{1}{2} f(\bar{q}) (p^{\text{ext}} - \bar{q}) = H(\bar{q}) < H(\bar{p}_*) = F(0)$.

Moreover, $G(p^{\text{ext}}) = F(p^{\text{ext}}) > F(\beta) = F(0)$, $G(0) > F(0)$. Thus, there exists a maximal interval $(q_*, q^*) \subset [0, p^{\text{ext}}]$ such that $G(q) \in (F(\theta), F(0))$ for all $q \in (q_*, q^*)$. We have $0 < q_* < \bar{q} < q^* < p^{\text{ext}}$ and $G(q_*) = G(q^*) = F(0)$. Therefore, $\mathcal{F}_2$ is well-defined and continuous in $(q_*, q^*)$, and $\lim_{p \to q^*} \mathcal{F}_2(p) = \lim_{p \to q_*} \mathcal{F}_2(p) = +\infty$. Reasoning like in the previous case, (24) admits solution if and only if $L \ge M_3$, where

$$M_3 := \min_{[q_*, q^*]} \mathcal{F}_2 > 0, \qquad (26)$$

On the other hand, if $D \ge D_*$, $\min_{[0,1]} G \ge F(0)$, and equation (24) has no solution. $\square$

**Proof of Theorem 2.1.** As we showed in Section 3.1.1, the (SD) steady-state solution $p$ of (2) has $p(L)$ satisfying equation (21). From Proposition 3.1, we can deduce that for fixed $p^{\text{ext}} \in (0, 1)$, $D > 0$, $M_d(p^{\text{ext}}, D) = \min_q \mathcal{F}_1(q)$. Thus, we obtain the results for (SD) steady-state solutions of (2) in Theorem 2.1.

Similarly, Proposition 3.2 provides that for fixed $p^{\text{ext}} \in (0, 1)$, $D > 0$, we have $M_i(p^{\text{ext}}, D) = \min_q \mathcal{F}_2(q)$ when $p^{\text{ext}} \le \beta$ or $D < D_*$. Otherwise, $M_i(p^{\text{ext}}, D) = +\infty$. $\square$

### 3.1.3 Existence of non-(SM) solutions

As we can see in the phase portrait in Figure 3, there exist some solutions of (2) which are neither (SD) nor (SI). These solutions can be non-symmetric or can have more than one (local) extremum. By studying these cases, we prove Theorem 2.2 as follows

**Proof of Theorem 2.2.** We can see from Figure 3a that for fixed $p^{\text{ext}} \le \beta$, $D > 0$, the non-(SM) solutions $p$ of (2) have more than one (local) extreme value because their orbits have at least two intersections with the axis $p' = 0$ (see e.g. $T_3$). Those solutions have the same local minimum values, denoted $p_{\min}$, and the same maximum values, denoted $p_{\max}$. Moreover, we have $p_{\min} < \theta < p_{\max}$, and $F(p_{\min}) = F(p_{\max})$.

Since the orbits make a round trip of distance $2L$, then the more extreme values a solution has, the larger $L$ is. Hence, to find the minimal value $M_*$, we study the case when $p$ has one local minimum and one local maximum with orbit as $T_3$ in Figure 3a. Then we have

$$G(p(-L)) = G(p(L)) = F(p_{\min}) = F(p_{\max}), \qquad (27)$$

and by using (16), we obtain

$$2L = \mathcal{F}_1((p(-L)) + \int_{p_{\min}}^{p_{\max}} \frac{ds}{\sqrt{2F(p_{\min}) - 2F(s)}} + \mathcal{F}_2(p(L)) = 2[\mathcal{F}_1(p(-L)) + \mathcal{F}_2(p(L))] + \int_{p(L)}^{p(-L)} \frac{ds}{\sqrt{2G(p(L)) - 2F(s)}}.$$

Using the same idea as above, we can show that $L$ depends continuously on $p(L)$. Moreover, we know that $M_d = \min \mathcal{F}_1$, $M_i = \min \mathcal{F}_2$, therefore there exists a constant $M_*$ such that (2) admits at least one non-(SM) solution $p$ if and only if $L \ge M_* > M_d + M_i$.
On the other hand, for fixed $p^{\text{ext}} > \beta, D < D_*$, it is possible that (2) admits a non-symmetric solution with only one minimum. The orbit of this solution is as $T_4$ in Figure 3b. In this case, we have $G(p(L)) = G(p(-L)) = F(p_{\min})$ with $p(-L) < p(L)$ and

$$2L = \mathcal{F}_2(p(-L)) + \mathcal{F}_2(p(L)) > 2M_i.$$

Hence, we only need $M_* > M_i$. $\square$

## 3.2 Stability analysis

We first study the principal eigenvalue and eigenfunction for the linear problem. Then by using these eigenelements, we construct the super- and sub-solution of (1) and prove the stability and instability corresponding to each case in Theorem 2.3.

**Proof of Theorem 2.3.** Consider the corresponding linear eigenvalue problem (12). We can see that $\phi = \cos(\sqrt{\lambda}x)$ is an eigenfunction iff $\sqrt{\lambda} \tan(L\sqrt{\lambda}) = D$. Denote $\lambda_1$ the smallest positive value of $\lambda$ which satisfies this equality, thus $L\sqrt{\lambda_1} \in (0, \frac{\pi}{2})$. Hence, $\lambda_1 \in (0, \frac{\pi^2}{4L^2})$. Moreover, for any $x \in (-L, L)$, the corresponding eigenfunction $\phi_1(x) = \cos(\sqrt{\lambda_1}x)$ takes values in $(0, 1)$.

**Proof of stability:** Now let $p$ be a steady-state solution of (1) governed by (2). First, we prove that if $f'(p(x)) < \lambda_1$ for any $x \in (-L, L)$ then $p$ is asymptotically stable. Indeed, since $f'(p(x)) < \lambda_1$, there exist positive constants $\delta, \gamma$ with $\gamma < \lambda_1$ such that for any $\eta \in [0, \delta]$,

$$f(p + \eta) - f(p) \leq (\lambda_1 - \gamma)\eta, \quad f(p) - f(p - \eta) \leq (\lambda_1 - \gamma)\eta, \qquad (28)$$

on $(-L, L)$. Now consider

$$\bar{p}(t, x) = p(x) + \delta e^{-\gamma t} \phi_1(x), \quad \underline{p}(t, x) = p(x) - \delta e^{-\gamma t} \phi_1(x).$$

Assume that $p^{\text{init}}(x) \leq p(x) + \delta\phi_1(x)$. Then by (28), we have that $\bar{p}$ is a super-solution of (1) because

$$\partial_t \bar{p} - \partial_{xx} \bar{p} = (\lambda_1 - \gamma) \delta e^{-\gamma t} \phi_1(x) + f(p) \geq f(p + \delta e^{-\gamma t} \phi_1(x)) = f(\bar{p}),$$

due to the fact that $0 < \delta e^{-\gamma t} \phi_1(x) < \delta$ for any $t > 0$, $x \in (-L, L)$. Moreover, at the boundary points one has $\frac{\partial \bar{p}}{\partial \nu} + D(\bar{p} - p^{\text{ext}}) = \frac{\partial p}{\partial \nu} + D(p - p^{\text{ext}}) = 0$.

Similarly, if we have $p^{\text{init}}(x) \geq p(x) - \delta\phi_1(x)$, and so $\underline{p}$ is a sub-solution of (1). Then, by the method of super- and sub-solution (see e.g. [19]), the solution $p^0$ of (1) satisfies $\underline{p} \leq p^0 \leq \bar{p}$. Hence, $|p^0(t, x) - p(x)| \leq \delta e^{-\gamma t} \phi_1(x)$. Therefore, we can conclude that, whenever $|p^{\text{init}}(x) - p(x)| \leq \delta\phi_1(x)$ for any $x \in (-L, L)$, the solution $p^0$ of (1) converges to the steady-state $p$ when $t \to +\infty$. This shows the stability of $p$.

**Proof of instability:** In the case when $f'(p(x)) > \lambda_1$, there exist positive constants $\delta, \gamma$, with $\gamma < \lambda_1$, such that for any $\eta \in [0, \delta]$,

$$f(p + \eta) - f(p) \geq (\lambda_1 + \gamma)\eta, \qquad (29)$$

on $(-L, L)$.

For any $p^{\text{init}} > p$, there exists a positive constant $\sigma < 1$ such that $p^{\text{init}} \geq p + \delta(1 - \sigma)$. Then $\tilde{p}(t, x) = p(x) + \delta(1 - \sigma e^{-\gamma't})\phi_1(x)$, with $\gamma' < \gamma$ small enough, is a sub-solution of (1). Indeed, by applying (29) with $\eta = \delta(1 - \sigma e^{-\gamma't})\phi_1(x) \in [0, \delta]$ for any $x \in (-L, L)$, we have

$$\partial_t \tilde{p} - \partial_{xx} \tilde{p} = \gamma' \delta \sigma e^{-\gamma' t} \phi_1(x) + \lambda_1 \delta (1 - \sigma e^{-\gamma' t}) \phi_1(x) + f(p) \leq f(p + \delta(1 - \sigma e^{-\gamma' t}) \phi_1(x))$$

if $\gamma \geq \frac{\gamma' \sigma e^{-\gamma' t}}{1 - \sigma e^{-\gamma' t}} = \frac{\gamma' \sigma}{e^{\gamma' t} - \sigma}$ for any $t \geq 0$. This inequality holds when we choose $\gamma' \leq \frac{\gamma(1-\sigma)}{\sigma}$. Now, we have that $\tilde{p}$ is a sub-solution of (1), thus for any $t \geq 0$, $x \in (-L, L)$, the corresponding solution $p^0$ satisfies

$$p^0(t, x) - p(x) \geq \tilde{p}(t, x) - p(x) \geq \delta(1 - \sigma e^{-\gamma't})\phi_1(x).$$

Hence, for a given positive $\epsilon < \delta \min_x \phi_1(x)$, when $t \to +\infty$, solution $p^0$ cannot remain in the $\epsilon$-neighborhood of $p$ even if $p^{\text{init}} - p$ is small. This implies the instability of $p$. $\square$
*Proof of Corollary 2.1.* For $p^{\text{ext}} \le \alpha_1 < \theta, D > 0$, from Theorem 2.1, the (SI) steady-state solution $p$ exists for any $L > 0$ and is unique, $p(x) \le p^{\text{ext}} \le \alpha_1$ for all $x \in (-L, L)$. Moreover, from Assumption 2.2, the reaction term $f$ has $f'(q) < 0$, for any $q \in (0, \alpha_1)$. Then, for any $x \in (-L, L)$, $f'(p(x)) \le 0 < \lambda_1$. Hence, $p$ is asymptotically stable.

Besides, from Theorems 2.1 and 2.2, for any $L > 0$ such that $L < M_d(p^{\text{ext}}, D) < M_*(p^{\text{ext}}, D)$, (1) has neither (SD) nor non-(SM) steady-state solutions. So the (SI) steady-state solution is the unique steady-state solution.

Using a similar argument for the case $p^{\text{ext}} \ge \alpha_2$, we obtain Corollary 2.1. □

# 4 Application to the control of dengue vectors by the introduction of the bacterium Wolbachia

## 4.1 Model

In this section, we show an application of our model to the control of mosquitoes using Wolbachia. Mosquitoes of genus *Aedes* are the vector of many dangerous arboviruses, such as dengue, zika, chikungunya and others. There exists neither effective treatment nor vaccine for these vector-borne diseases, and in such conditions, the main method to control them is to control the vector population. A biological control method using a bacterium called Wolbachia (see [10]) was discovered and developed with this purpose. Besides reducing the ability of mosquitoes to transmit viruses, Wolbachia also causes an important phenomenon called cytoplasmic incompatibility (CI) on mosquitoes. More precisely, if a wild female mosquito is fertilized by a male carrying Wolbachia, its eggs almost cannot hatch. For more details about CI, we refer to [32]. In the case of *Aedes* mosquitoes, Wolbachia reduces lifespan, changes fecundity, and blocks the development of the virus. However, it does not influence the way mosquitoes move.

In [26], model (3), (4) was considered with $n_1 = n_i$ the density of the mosquitoes which are infected by Wolbachia and $n_2 = n_u$ the density of wild uninfected mosquitoes. Consider the following positive parameters:

* $d_u, \delta d_u$: death rate of, respectively uninfected mosquitoes and infected mosquitoes, $\delta > 1$ since Wolbachia reduces the lifespan of the mosquitoes;

* $b_u, (1 - s_f)b_u$: birth rate of, respectively uninfected mosquitoes and infected ones. Here $s_f \in [0, 1)$ characterizes the fecundity decrease;

* $s_h \in (0, 1]$: the fraction of uninfected females' eggs fertilized by infected males that do not hatch, due to the cytoplasmic incompatibility (CI);

* $K$: carrying capacity, $A$: diffusion coefficient.

Parameters $\delta, s_f, s_h$ have been estimated in several cases and can be found in the literature (see [1] and references therein). We always assume that $s_f < s_h$ (in practice, $s_f$ is close to 0 while $s_h$ is close to 1).

Several models have been proposed using these parameters. In the present study, a system of Lotka-Volterra type is proposed, where the parameter $\epsilon > 0$ is used to characterize the high fertility as follows

$$ \partial_t n_i^{\epsilon} - A \partial_{xx} n_i^{\epsilon} = (1 - s_f) \frac{b_u}{\epsilon} n_i^{\epsilon} \left( 1 - \frac{n_i^{\epsilon} + n_u^{\epsilon}}{K} \right) - \delta d_u n_i^{\epsilon}, \quad (30a) $$

$$ \partial_t n_u^{\epsilon} - A \partial_{xx} n_u^{\epsilon} = \frac{b_u}{\epsilon} n_u^{\epsilon} \left( 1 - s_h \frac{n_i}{n_i^{\epsilon} + n_u^{\epsilon}} \right) \left( 1 - \frac{n_i^{\epsilon} + n_u^{\epsilon}}{K} \right) - d_u n_u^{\epsilon}, \quad (30b) $$

where the reaction term describes birth and death. The factor $(1 - s_h \frac{n_i^{\epsilon}}{n_i^{\epsilon} + n_u^{\epsilon}})$ characterizes the cytoplasmic incompatibility. Indeed, when $s_h = 1$, no egg of uninfected females fertilized by infected males can hatch, that is, there is complete cytoplasmic incompatibility. The factor becomes $\frac{n_i^{\epsilon}}{n_i^{\epsilon} + n_u^{\epsilon}}$ which means the birth rate of uninfected mosquitoes depends on the proportion of uninfected parents because only an uninfected couple can lay uninfected eggs. Whereas, $s_h = 0$ means that all the eggs of uninfected females hatch. In this case, the factor $(1 - s_h \frac{n_i^{\epsilon}}{n_i^{\epsilon} + n_u^{\epsilon}})$ becomes 1, so the growth rate of uninfected population is not altered by the pressure of the infected one.
In paper [26], the same model was studied in the entire space $\mathbb{R}$. In that case, the system (30) has exactly two stable equilibria, namely the *Wolbachia* invasion steady state and the *Wolbachia* extinction steady state. In this paper, the authors show that when $\epsilon \to 0$ and the reaction terms satisfy some appropriate conditions, the proportion $p^\epsilon = \frac{n_i^\epsilon}{n_i^\epsilon+n_u^\epsilon}$ converges to the solution $p^0$ of the scalar equation $\partial_t p^0 - A\partial_{xx}p^0 = f(p^0)$, with the reaction term

$$f(p) = \delta d_u s_h \frac{p(1-p)(p-\theta)}{s_h p^2 - (s_f + s_h)p + 1}, \quad (31)$$

with $\theta = \frac{s_f+\delta-1}{\delta s_h}$. We will always assume that $s_f + \delta(1 - s_h) < 1$, so $\theta \in (0, 1)$, and $f$ is a bistable function on $(0, 1)$. The two stable steady states 1 and 0 of (1) correspond to the success or failure of the biological control using *Wolbachia*.

## 4.2 Mosquito population in presence of migration

In this study, the migration of mosquitoes is taken into account. Typically, the inflow of wild unin- fected mosquitoes and the outflow of the infected ones may influence the efficiency of the method using *Wolbachia*. Here, to model this effect, system (30) is considered in a bounded domain with appropriate boundary conditions to characterize the migration of mosquitoes. In one-dimensional space, we consider $\Omega = (-L, L)$ and Robin boundary conditions as in (4) at points $x = -L$, and $x = L$

$$\frac{\partial n_i^\epsilon}{\partial \nu} = -D(n_i^\epsilon - n_i^{\text{ext},\epsilon}), \quad (32a)$$

$$\frac{\partial n_u^\epsilon}{\partial \nu} = -D(n_u^\epsilon - n_u^{\text{ext},\epsilon}), \quad (32b)$$

where $n_i^{\text{ext},\epsilon}$, $n_u^{\text{ext},\epsilon}$ do not depend on $t$ and $x$ but depend on parameter $\epsilon > 0$. Denote $p^\epsilon = \frac{n_i^\epsilon}{n_i^\epsilon+n_u^\epsilon}$, $n^\epsilon = \frac{1}{\epsilon}(1 - \frac{n_i^\epsilon+n_u^\epsilon}{K})$. In Section A, we prove that when $\epsilon \to 0$, up to extraction of sub-sequences, $n^\epsilon$ converges weakly to $n^0 = h(p^0)$ for some explicit function $h$, and $p^\epsilon$ converges strongly towards solution $p^0$ of (1) where $p^{\text{ext}}$ is the limit of $\frac{n_i^{\text{ext},\epsilon}}{n_i^{\text{ext},\epsilon}+n_u^{\text{ext},\epsilon}}$ when $\epsilon \to 0$, and the reaction term $f$ as in (31). Function $f$ satisfies Assumptions 2.1 and 2.2, so the results in Theorem 2.1 and 2.3 can be applied to this problem. By changing the spatial scale, we can normalize the diffusion coefficient into $A = 1$.

In this application, the parameters $L, D, p^{\text{ext}}$ correspond to the size of $\Omega$, the migration rate of mosquitoes, and the proportion of infected mosquitoes surrounding the boundary. The main results in the present paper give information about the existence and stability of equilibria depending upon different conditions for these parameters. Especially, from Corollary 2.1, we obtain that when the size $L$ of the domain is small, there exists a unique equilibrium for this problem and its values depend on the proportion of mosquitoes carrying *Wolbachia* outside the domain ($p^{\text{ext}}$). More precisely, when $p^{\text{ext}}$ is small (i.e., $p^{\text{ext}} \le \alpha_1$), the solution of (1) converges to the steady-state solution close to 0, which corresponds to the extinction of mosquitoes carrying *Wolbachia*. Therefore, in this situation, the replacement strategy fails because of the migration through the boundary. Otherwise, when the proportion outside the domain is high (i.e., $p^{\text{ext}} \ge \alpha_2$), then the long-time behavior of solutions of (1) has values close to 1, which means that the mosquitoes carrying *Wolbachia* can invade the whole population.

## 4.3 Numerical illustration

In this section, we present the numerical illustration for the above results. Parameters are fixed according to biologically relevant data (adapted from [5]). Time unit is the day, and parameters per day are in Table 2. Then, the reaction term $f$ in (31) has $\theta = 0.2375$, $\beta \approx 0.3633$, $\alpha_1 \approx 0.12$, $\alpha_2 \approx 0.7$. As proposed in section 3 of the modeling article [17], we may pick the value $830m^2$ per day for the diffusivity of *Aedes* mosquitoes. Choose $A = 1$, so the $x$-axis unit in the simulation corresponds to $\sqrt{830}/1 \approx 29$ m.

In the following parts, we check the convergence of $p^\epsilon$ when $\epsilon \to 0$ in 4.3.1. In 4.3.2, corresponding to different parameters, we compute numerically the solutions of (1) and (2) to check their existence and stability.
<table><caption>Table 2: Parameters for the numerical illustration</caption><thead><tr><th>Parameters</th><th>b<sub>u</sub></th><th>d<sub>u</sub></th><th>&delta;</th><th>&sigma;</th><th>s<sub>f</sub></th><th>s<sub>h</sub></th></tr></thead><tbody><tr><td>Values</td><td>1.12</td><td>0.27</td><td><sup>10</sup>&frasl;<sub>9</sub></td><td>1</td><td>0.1</td><td>0.8</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Convergence of $p^\epsilon$ to $p^0$ as $\epsilon$ goes to zero. The solid lines represent the solution $p^0(t, x)$ of (1) at $t = 50$ days. The dashed lines represent the proportion $p^\epsilon = \frac{n_i^\epsilon}{n_i^\epsilon+n_u^\epsilon}$ of solution $n_i^\epsilon$, $n_u^\epsilon$ of system (30), (32) at $t = 50$.</figcaption></figure>

### 4.3.1 Convergence to the scalar equation

Consider a mosquito population with a large fecundity rate that is, $\epsilon \ll 1$. Model (30) with boundary condition in (32) takes into account the migration of mosquitoes.

Fix $D = 0.05$, $p^{\text{ext}} = 0.1$ and $L = 2$, the system (30), (32) is solved numerically thanks to a semi-implicit finite difference scheme with 3 different values of the parameters $\epsilon$. The initial data are chosen such that $n_i^\epsilon(t=0) = n_u^\epsilon(t=0)$, that is, $p^{\text{init}} = 0.5$. In Figure 4, at time $t = 50$ days, the numerical solutions of (1) are plotted with blue solid lines, the proportions $p^\epsilon = \frac{n_i^\epsilon}{n_i^\epsilon+n_u^\epsilon}$ are plotted with dashed lines. We observe that when $\epsilon$ goes to 0, the proportion $p^\epsilon$ converges to the solution $p^0$ of system (1).

### 4.3.2 Steady-state solutions

For the different values of $p^{\text{ext}}$, the values of the integrals $F_1$ and $F_2$ as functions of $p(L)$ in (20) and (23) are plotted in Figure 5. For fixed values of $D$ and $p^{\text{ext}}$, Figure 5 can play the role of bifurcation diagrams that show the relation between the value $p(L)$ of symmetric solutions $p$ and parameter $L$. Then, we can obtain the critical values of parameter $L$. Next, we compute numerically the (SM) steady-state solutions of (1) with different values of $L > 0, D > 0, p^{\text{ext}} \in (0, 1)$.

**Numerical method:** To approximate the (SM) steady-state solution, we use the Newton method to solve nonlinear equations and follow these steps:

* Step 1: Solve $L = F_i(p(L))$ for $i = 1$ or 2, and obtain the values of $p(L)$.

* Step 2: Find $p(0)$ by solving (17).

* Step 3: For each $x$ in $(0, L)$, interpolate $p(x)$ by solving $x = \int_{p(0)}^{p(x)} \frac{(-1)^k ds}{\sqrt{2F(p(0))-2F(s)}}$ due to (16) with $k = 1$ if $p$ is decreasing and $k = 2$ if $p$ is increasing on $(0, x)$.

The construction of a non-(SM) steady-state solution is more sophisticated since it is hard to find $p(L)$ for a fixed $L$ like in step 1. We presented a numerical non-(SM) equilibrium in Figure 6c where we first fixed a value $p(L)$. Then similarly to step 2, we solved (27) to find all the extreme values of $p$. Finally, we applied step 3 with $p(0)$ replaced by $p_{\min}$ or $p_{\max}$.

We also plot the time dynamics of the solution $p^0(t, x)$ of (1) at $t = 10, 20, 40, 60, 100$ to verify the asymptotic stability of steady-state solutions. Next, we consider different values of $p^{\text{ext}}$ and present our
<figure><img src="image_6.png" /><figcaption>Figure 5: The blue and red solid lines represent respectively functions $F_1$ and $F_2$ of $p(L)$</figcaption></figure>

observation in each case.

• Case 1: $p^{\text{ext}} = 0.1 < \alpha_1$.

For $D = 0.05$ fixed, we observe in Figure 5a that for any $L > 0$, equation $F_2(p(L)) = L$ always admits exactly one solution. Thus, there always exists one (SI) steady-state solution with small values. We approximate that

$$M_d(0.1, 0.05) = M_1 \approx 0.8819, \quad M_*(0.1, 0.05) \approx 8.625.$$

Also from Figure 5a, we observe that when $L = M_1$, a bifurcation occurs and (1) admits an (SD) steady-state solution, and when $L > M_1$ one can obtain two (SD) solutions. Moreover, when $L \ge M_*$, there exist non-symmetric steady-state solutions. We do numerical simulations for two values of $L$ as follows.

For $L = 0.5 < M_1$, the unique equilibrium $\bar{p}_{21}$ is (SI) and has values close to 0 (see Figure 6a). Solution $p^0$ of (1) with any initial data converges to $\bar{p}_{21}$. This simulation is coherent with the asymptotic stability that we proved in Corollary 2.1.

For $L = 8.96 > M_* > M_1$, together with $\bar{p}_{21}$, there exist two more (SD) steady-state solutions, namely $\bar{p}_{11}$, $\bar{p}_{12}$, (see Figure 6b). This plot shows that these steady-state solutions are ordered, and the time-dependent solutions converge to either the largest one $\bar{p}_{11}$ or the smallest one $\bar{p}_{21}$, while $\bar{p}_{12}$ with intermediate values is not an attractor. In Figure 6c, we find numerically a non-symmetric solution $\bar{p}$ of (2) corresponding to orbit $T_3$ as in Figure 3a. Let the initial value $p^{\text{init}} = \bar{p}$, then we observe from Figure 6c that $p^0$ still converges to the symmetric equilibrium $\bar{p}_{21}$.

Moreover, the value $\lambda_1$ of Theorem 2.3 in this case is approximately equal to 0.0063. We also obtain that for any $x \in (-L, L)$,

$$f'(\bar{p}_{11}(x)) < 0, \quad f'(\bar{p}_{21}(x)) < 0, \quad f'(\bar{p}_{12}(x)) > 0.0462, \quad f'(\bar{p}(x)) > 0.022.$$

Therefore, by applying Theorem 2.3, we deduce that the steady-state solutions $\bar{p}_{11}$, $\bar{p}_{21}$ are asymptotically stable, $\bar{p}_{12}$ and the non-symmetric equilibrium $\bar{p}$ are unstable. Thus, the numerical simulations in Figure 6 are coherent to the theoretical results that we proved.

• Case 2: $p^{\text{ext}} = 0.8 > \alpha_2 > \beta$.

In this case, we obtain $D_* \approx 0.16$. We present numerical illustrations for two cases: $D = 0.05 < D_*$ and $D = 0.5 > D_*$.

◦ For $D = 0.05 < D_*$, we have $M_i(0.8, 0.05) = M_2 \approx 10.3646$ (see Figure 5b).

For $L = 2 < M_2$, the unique equilibrium $\bar{p}_{11}$ is (SD) and has values close to 1 (see Figure 7a). The time-dependent solution $p^0$ of (1) with any initial data converges to $\bar{p}_{11}$. This simulation is coherent to the asymptotic stability we obtained in Corollary 2.1.

For $L = 12 > M_2$, together with $\bar{p}_{11}$, there exist two more (SI) steady-state solutions, namely $\bar{p}_{21}$, $\bar{p}_{22}$, and they are ordered (see Figure 7b). In this case, we obtain approximately that $\lambda_1 \approx 0.0063$ and for any $x \in (-L, L)$, one has

$$f'(\bar{p}_{11}(x)) < 0, \quad f'(\bar{p}_{21}(x)) \in (-0.0398, 0.0368), \quad f'(\bar{p}_{22}(x)) \in (-0.0195, 0.0673).$$
<figure><img src="image_7.png" /><figcaption>Figure 6: Case $p^{\text{ext}} = 0.1, D = 0.05$: The solid lines illustrate the steady-state solutions. The dotted lines show the initial data of problem (1). The dashed lines represent the solution $p^0(t, x)$ with $t \in \{10, 20, 40, 60, 100\}$. The color of the dashed lines corresponds to the color of the equilibrium that they converge to.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Case $p^{\text{ext}} = 0.8$: The solid lines illustrate the steady-state solutions. The dotted lines show the initial data of problem (1). The dashed lines represent the solution $p^0(t, x)$ with $t \in \{10, 20, 40, 60, 100\}$. The color of the dashed lines corresponds to the color of the equilibrium that they converge to.</figcaption></figure>

By sufficient conditions in Theorem 2.3, we obtain that $\bar{p}_{11}$ is asymptotically stable but we can not conclude the stability for $\bar{p}_{21}$ and $\bar{p}_{22}$. The time dynamics of $p^0$ in Figure 7b suggests that the smallest steady-state solution $\bar{p}_{21}$ is asymptotically stable and $\bar{p}_{22}$ seems to be unstable.

◦ For $D = 0.5 > D_*$, function $F_2$ is not defined (see Figure 5c), so problem (2) admits only one (SD) steady-state, and we obtain that it is unique and asymptotically stable (see Figure 7c).

5 Conclusion and perspectives

We have studied the existence and stability of steady-state solutions with values in $[0, 1]$ of a reaction-diffusion equation

$$ \partial_t p - \partial_{xx} p = f(p) $$

on an interval $(-L, L)$ with cubic nonlinearity $f$ and inhomogeneous Robin boundary conditions

$$ \frac{\partial p}{\partial \nu} = D(p - p^{\text{ext}}), $$

where constant $p^{\text{ext}} \in (0, 1)$ is an analogue of $p$, and constant $D > 0$. We have shown how the analysis of this problem depends on the parameters $p^{\text{ext}}$, $D$, and $L$. More precisely, the main results say that there always exists a symmetric steady-state solution that is monotone on each half of the domain. For $p^{\text{ext}}$ large, the value of this steady-state solution is close to 1, otherwise, it is close to 0. Besides, the larger value of $L$, the more steady-state solutions this problem admits. We have found the critical values of $L$ so that when the parameters surpass these critical values, the number of steady-state solutions increases. We also provided some sufficient conditions for the stability and instability of the steady-state solutions.

We presented an application of our results on the control of dengue vector using Wolbachia bacterium that can be transmitted maternally. Since Wolbachia can help reduce vectorial capacity of the mosquitoes, the main goal of this method is to replace wild mosquitoes by mosquitoes carrying Wolbachia. In this
application, we considered $p$ as the proportion of mosquitoes carrying Wolbachia and used the equation above to model the dynamic of the mosquito population. The boundary condition describes the migration through the border of the domain. This replacement method only works when $p$ can reach an equilibrium close to 1. Therefore, the study of the existence and stability of the steady-state solution close to 1 is meaningful and depends strongly on the parameters $p^{\text{ext}}$, $D$, and $L$. In realistic situations, the proportion $p^{\text{ext}}$ of mosquitoes carrying Wolbachia outside the domain is usually low. Using the theoretical results proved in this article, one sees that, to have major chances of success, one should try to treat large regions ($L$ large), well isolated ($D$ small) and possibly applying a population replacement method in a zone outside $\Omega$ (to increase $p^{\text{ext}}$ by reducing its denominator).

As a natural continuation of the present work, higher dimension problems and more general boundary conditions can be studied. In more realistic cases, $p^{\text{ext}}$ can be considered to depend on space and the periodic solutions can be the next problem for our study. Besides, when an equilibrium close to 1 exists and is stable, one may consider multiple strategies using multiple releases of mosquitoes carrying Wolbachia. To optimize the number of mosquitoes released to guarantee the success of this method under the difficulties enlightened by this paper is an interesting problem for future works.

## A Asymptotic limit of reaction-diffusion systems

In [26], the authors reduced a 2-by-2 reaction-diffusion system of Lotka-Volterra type modeling two biological populations to a scalar equation as in (1) when the fecundity rate is very large. This limit problem was first proved in the whole domain. In the present study, we prove the limit for a system in a bounded domain with inhomogeneous Robin boundary conditions. In the following part, we recall the necessary assumptions and present results about this problem.

Although the main result of the paper is in one-dimensional space, the following result holds in any dimension $d$. Let $\Omega \subset \mathbb{R}^d$ be a bounded domain and consider the initial-boundary-value problem (33) depending on parameter $\epsilon > 0$,

$$ \partial_t n_1^\epsilon - \Delta n_1^\epsilon = n_1^\epsilon f_1^\epsilon(n_1^\epsilon, n_2^\epsilon), \quad (t, x) \in (0, T) \times \Omega, \quad (33a) $$

$$ \partial_t n_2^\epsilon - \Delta n_2^\epsilon = n_2^\epsilon f_2^\epsilon(n_1^\epsilon, n_2^\epsilon), \quad (t, x) \in (0, T) \times \Omega, \quad (33b) $$

$$ n_1^\epsilon(0, x) = n_1^{\text{init}, \epsilon}(x), \quad n_2^\epsilon(0, x) = n_2^{\text{init}, \epsilon}(x), \quad x \in \Omega, \quad (33c) $$

$$ \frac{\partial n_1^\epsilon}{\partial \nu} = -D(n_1^\epsilon - n_1^{\text{ext}, \epsilon}), \quad \frac{\partial n_2^\epsilon}{\partial \nu} = -D(n_2^\epsilon - n_2^{\text{ext}, \epsilon}), \quad (t, x) \in (0, T) \times \partial\Omega, \quad (33d) $$

where we assume that $f_1^\epsilon, f_2^\epsilon$ are smooth enough to guarantee existence and uniqueness of a classical solution for fixed $\epsilon$. More precisely, the following assumptions are made:

**Assumption A.1.** (Initial and boundary conditions). Assume that $n_1^{\text{init}, \epsilon}, n_2^{\text{init}, \epsilon} \in L^\infty(\Omega)$ with $n_1^{\text{init}, \epsilon}, n_2^{\text{init}, \epsilon} \ge 0$ and $n_2^{\text{init}, \epsilon}$ is not identical to 0.

$D > 0$ is constant, $n_1^{\text{ext}, \epsilon} \ge 0, n_2^{\text{ext}, \epsilon} > 0$ do not depend on time $t$ and position $x$.

To study the limit problem, we define the "rescaled total population" $n^\epsilon$ and proportion $p^\epsilon$, by

$$ n^\epsilon := \frac{1}{\epsilon} - n_1^\epsilon - n_2^\epsilon, \quad p^\epsilon := \frac{n_1^\epsilon}{n_1^\epsilon + n_2^\epsilon}. $$

Next, we recall some assumptions that were proposed in [26] on the families of functions $(f_1^\epsilon, f_2^\epsilon)_{\epsilon>0}$ to study the convergence of $p^\epsilon$ when $\epsilon \to 0$

**Assumption A.2.** Function $f_1^\epsilon, f_2^\epsilon$ are of class $C^2(\mathbb{R}_+^2\{0\})$, and for $i \in \{1, 2\}$ there exists $F_i \in C^2(\mathbb{R}^2)$ (independent of $\epsilon$) such that

$$ f_i^\epsilon(n_1^\epsilon, n_2^\epsilon) = F_i(n^\epsilon, p^\epsilon). \quad (34) $$

That is, we may write $f_i^\epsilon(n_1^\epsilon, n_2^\epsilon) = F_i(\frac{1}{\epsilon} - n_1^\epsilon - n_2^\epsilon, \frac{n_1^\epsilon}{n_1^\epsilon + n_2^\epsilon})$ for $i \in \{1, 2\}$.

Then, we can deduce that $p^\epsilon$ and $n^\epsilon$ satisfy system (35) as follows:

In $(0, T) \times \Omega$, we have

$$ \partial_t n^\epsilon - \Delta n^\epsilon = -(\frac{1}{\epsilon} - n^\epsilon) [p^\epsilon F_1(n^\epsilon, p^\epsilon) + (1-p^\epsilon) F_2(n^\epsilon, p^\epsilon)], \quad (35a) $$
$$ \partial_t p^\epsilon - \Delta p^\epsilon + \frac{2\epsilon A}{1 - \epsilon n^\epsilon} \nabla p^\epsilon \cdot \nabla n^\epsilon = p^\epsilon(1 - p^\epsilon)(F_1 - F_2)(n^\epsilon, p^\epsilon), \quad (35b) $$

on the boundary $(0, T) \times \partial\Omega$, we have

$$ \frac{\partial n^\epsilon}{\partial \nu} = -D(n^\epsilon - n^{\text{ext},\epsilon}), \quad \frac{\partial p^\epsilon}{\partial \nu} = -D(p^\epsilon - p^{\text{ext},\epsilon}) \frac{1 - \epsilon n^{\text{ext},\epsilon}}{1 - \epsilon n^\epsilon}, \quad (35c) $$

at time $t=0$, for any $x \in \Omega$, the initial data read

$$ n^{\epsilon}(0, x) = n^{\text{init},\epsilon}(x), \quad p^{\epsilon}(0, x) = p^{\text{init},\epsilon}(x), \quad (35d) $$

where $(F_1 - F_2)(n^\epsilon, p^\epsilon) = F_1(n^\epsilon, p^\epsilon) - F_2(n^\epsilon, p^\epsilon)$, and

$$ n^{\text{init},\epsilon} := \frac{1}{\epsilon} - n_1^{\text{init},\epsilon} - n_2^{\text{init},\epsilon}, \quad p^{\text{init},\epsilon} := \frac{n_1^{\text{init},\epsilon}}{n_1^{\text{init},\epsilon} + n_2^{\text{init},\epsilon}}, $$

$$ n^{\text{ext},\epsilon} := \frac{1}{\epsilon} - n_1^{\text{ext},\epsilon} - n_2^{\text{ext},\epsilon}, \quad p^{\text{ext},\epsilon} := \frac{n_1^{\text{ext},\epsilon}}{n_1^{\text{ext},\epsilon} + n_2^{\text{ext},\epsilon}}. $$

Let us denote $H(n, p) = -pF_1(n, p) - (1-p)F_2(n, p)$. The following assumption guarantees existence of zeros of $H$ given by $(n, p) = (h(p), p)$ for each $p \in [0, 1]$.

**Assumption A.3.** In addition to Assumption A.2,

(i) There exists $B > 0$ such that for all $n \ge 0$, $p \in [0, 1]$, $\partial_n H(n, p) \le -B$,

(ii) For all $p > 0$, $H(0, p) > 0$.

Conditions (i) and (ii) imply that for all $p \in [0, 1]$, there exists a unique $n =: h(p) \in \mathbb{R}_+^*$ such that $H(n, p) = 0$. We have $H \in C^2(\mathbb{R}_+^2)$ (from Assumption A.2) thus $h \in C^2(0, 1)$, with $H(h(p), p) = 0$ for all $p \in [0, 1]$.

The following assumptions are made for the initial data and boundary conditions

**Assumption A.4.** There exists a function $p^{\text{init}} \in L^2(\Omega)$ such that $p^{\text{init},\epsilon} \underset{\epsilon \to 0}{\rightharpoonup} p^{\text{init}}$ weakly in $L^2(\Omega)$.
Function $n^{\text{init},\epsilon} - h(0) \in L^2 \cap L^\infty(\Omega)$ is uniformly bounded in $\epsilon > 0$.

**Assumption A.5.** There exists positive constants $\tilde{\epsilon} > 0$, $\tilde{K} > 0$ such that for any $\epsilon \in (0, \tilde{\epsilon})$, we have $|n^{\text{ext},\epsilon}| < \tilde{K}$.

There exists a constant $p^{\text{ext}} \in (0, 1)$ not depending on $\epsilon$ such that $p^{\text{ext},\epsilon} \underset{\epsilon \to 0}{\rightharpoonup} p^{\text{ext}}$

**Convergence result.** For fixed $\epsilon > 0$, existence of solutions of (35) is classical (see, e.g. [20]). Following the idea in [26], we present the asymptotic limit of the proportion $p^\epsilon$ and $n^\epsilon$ in the following theorem.

**Theorem A.1.** Assume that Assumptions A.1-A.5 are satisfied and consider the solution $(n^\epsilon, p^\epsilon)$ of (35). Then, for all $T > 0$, we have the convergence

$$ \left\{ \begin{array}{ll} p^{\epsilon} \underset{\epsilon \to 0}{\longrightarrow} p^0 & \text{strongly in } L^2(0, T; L^2(\Omega)), \text{ weakly in } L^2(0, T; H^1(\Omega)), \\ n^{\epsilon} - h(p^{\epsilon}) \underset{\epsilon \to 0}{\longrightarrow} 0 & \text{strongly in } L^2(0, T; L^2(\Omega)), \text{ weakly in } L^2(0, T; H^1(\Omega)), \end{array} \right. $$

where $p^0$ is the unique solution of

$$ \left\{ \begin{array}{ll} \partial_t p^0 - \Delta p^0 = p^0(1-p^0)(F_1 - F_2)(h(p^0), p^0), & \text{in } (0, T) \times \Omega, \\ p^0(0, \cdot) = p^{\text{init}} & \text{in } \Omega, \\ \frac{\partial p^0}{\partial \nu} = -D(p^0 - p^{\text{ext}}) & \text{on } (0, T) \times \partial \Omega. \end{array} \right. $$

We recall the apriori estimates of [26] without proof and present some bounds on the boundary in Appendix A.1. Then we use the Aubin-Lions lemma and trace theorem to prove the limit in Appendix A.2.
## A.1 Uniform a priori estimates

First, we establish the uniform bound with respect to $\epsilon$ in $L^\infty$ in the following lemma

**Lemma A.1.** Under Assumptions A.1-A.5, for a given value $\epsilon > 0$, let $(n^\epsilon, p^\epsilon)$ be the unique solution of (35). Then, for any $T > 0$, $0 \le p^\epsilon \le 1$ in $[0, T] \times \bar{\Omega}$ for all $\epsilon > 0$. Also, there exists $\epsilon_0 > 0$, $K_0 > 0$ such that for any $\epsilon \in (0, \epsilon_0)$, $||n^\epsilon||_{L^\infty([0,T] \times \Omega)} \le K_0$.

Moreover, $n^\epsilon$ is uniformly bounded on $[0, T] \times \partial\Omega$.

*Proof.* Using the same method as in Lemma 5 of [26], we obtain the uniform bounds for $p^\epsilon$ in $[0, T] \times \bar{\Omega}$, and for $n^\epsilon$ in $L^\infty([0, T] \times \Omega)$.

Moreover, for any $x \in \partial\Omega$, let $\nu$ be the normal outward vector through $x$. Then, for $\delta > 0$ small enough, $x - \delta\nu \in \Omega$. From the boundary condition for $n^\epsilon$ in (35), one has for $t \in [0, T]$,

$$ \lim_{\delta \to 0^+} \frac{n^\epsilon(t, x) - n^\epsilon(t, x - \delta\nu)}{\delta} = -D(n^\epsilon(t, x) - n^{\text{ext}, \epsilon}). $$

So for any $\eta > 0$, there exists $\delta > 0$ small such that

$$ \left| \frac{n^\epsilon(t,x)-n^\epsilon(t,x-\delta\nu)}{\delta} + D(n^\epsilon(t,x)-n^{\text{ext},\epsilon}) \right| \le \eta. $$

Thus, $n^\epsilon(t, x)(1 + \delta D) \le n^\epsilon(t, x - \delta\nu) + \delta D n^{\text{ext},\epsilon} + \delta\eta$, then for $\eta$ and $\delta$ small enough, for any $\epsilon < \epsilon_0$, $t \in [0, T]$, $x \in \partial\Omega$, since $x - \delta\nu \in \Omega$, one has $|n^\epsilon(t, x)| \le K_0 + \delta D \tilde{K} + \delta\eta < K_1$. Then $n^\epsilon$ is uniformly bounded on $[0, T] \times \partial\Omega$ and $||n^\epsilon||_{L^\infty([0,T] \times \partial\Omega)} \le K_1$. $\square$

The following lemmas can be proved analogously to the proof in [26].

**Lemma A.2.** Under Assumptions A.1-A.5, for $\epsilon > 0$ small enough, let $(n^\epsilon, p^\epsilon)$ be the unique solution of (35). We have the following uniform estimates

$$ \epsilon \int_0^T \int_{\Omega} |\nabla n^\epsilon|^2 dxdt \le C_0, \quad \int_0^T \int_{\Omega} |\nabla p^\epsilon|^2 dxdt \le \bar{C}, \qquad (36) $$

for some positive constants $C_0$ and $\bar{C}$.

Denote $M^\epsilon := n^\epsilon - h(p^\epsilon)$ where $h$ is defined in Assumption A.3. The following provide the convergence of $M^\epsilon$.

**Lemma A.3.** Let $T > 0$, under Assumptions A.1-A.5, one has $M^\epsilon \to 0$ in $L^2(0, T; L^2(\Omega))$ when $\epsilon \to 0$.

Now, we provide a uniform estimate for $\partial_t p^\epsilon$ with respect to $\epsilon$ in the following lemma.

**Lemma A.4.** Under Assumptions A.1-A.5, for $\epsilon > 0$ small enough, $\partial_t p^\epsilon$ is uniformly bounded in $L^2(0, T; X')$ with respect to $\epsilon$, where $X = H^1(\Omega) \cap L^\infty(\Omega)$.

## A.2 Proof of convergence

The idea to prove Theorem A.1 is relied on the relative compactness obtained from the Aubin-Lions lemma below (see [23])

**Lemma A.5 (Aubin-Lions).** Let $T > 0$, $q \in (1, \infty)$, and $(\psi_n)_n$ a bounded sequence in $L^q(0, T; B)$, where $B$ is a Banach space. If $(\psi_n)$ is bounded in $L^q(0, T; X)$ and $X$ embeds compactly in $B$, and if $(\partial_t \psi_n)_n$ is bounded in $L^q(0, T; X')$ uniformly with respect to $n$, then $(\psi_n)_n$ is relatively compact in $L^q(0, T; B)$.

**Proof of Theorem A.1.** We use 3 steps to proof Theorem A.1. First, we obtain the relative compactness of $(p^\epsilon)$ by applying Aubin-Lions lemma, and prove that there exists (up to extracting subsequences) a limit function. Then, we study its behavior on the boundary using the trace theorem. Finally, thanks to our uniform bounds, we show that the limit function satisfies a problem whose solution is unique.

**Step 1:** In our problem, we need to apply the Lions-Aubin lemma with $q = 2, B = L^2(\Omega)$ and $X = H^1(\Omega) \cap L^\infty(\Omega)$ to $(\psi_\epsilon) = (p^\epsilon)_\epsilon$. The compact embedding from $X$ to $B$ is valid by the Rellich-Kondrachov theorem. In the previous section, we have already obtained uniform estimates that are sufficient to apply the Aubin-Lions lemma. The sequence $(p^\epsilon)_\epsilon$ is bounded in $L^2(0, T; L^2(\Omega))$ due to Lemma A.1
$$ ||p^{\epsilon}||_{L^2(0,T;L^2(\Omega))}^2 = \int_{0}^{T} \int_{\Omega} |p^{\epsilon}|^2 dx dt \le ||p^{\epsilon}||_{L^{\infty}(0,T;L^2(\Omega))}^2 \text{meas}(\Omega)T < \infty, $$

for $\epsilon < \epsilon_0$ small enough. Then, due to Lemma A.2, this sequence is bounded in $L^2(0, T; X)$. The sequence $(\partial_t p^\epsilon)_\epsilon$ is bounded in $L^2(0, T; X')$ by Lemma A.4. Thus, we can apply Aubin-Lions lemma and deduce that $(p^\epsilon)_\epsilon$ is strongly relatively compact in $L^2(0, T; L^2(\Omega))$. Therefore, there exists $p^0 \in L^2(0, T; H^1(\Omega))$ such that, up to extraction of subsequences, we have $p^\epsilon \to p^0$ strongly in $L^2((0, T) \times \Omega)$ and a.e., $\nabla p^\epsilon \to \nabla p^0$ weakly in $L^2((0, T) \times \Omega)$.

Moreover, by the triangle inequality we have $|n^\epsilon - h(p^0)| \le |n^\epsilon - h(p^\epsilon)| + |h(p^\epsilon) - h(p^0)| \le |n^\epsilon - h(p^\epsilon)| + ||h'||_{L^\infty([0,1])} |p^\epsilon - p^0|$. From the strong convergence of $p^\epsilon$ and $M^\epsilon$ in Lemma A.3 when $\epsilon \to 0$, we can deduce the following strong convergence in $L^2(0, T; L^2(\Omega))$

$$ n^{\epsilon} \to n^{0} := h(p^{0}) \qquad (37) $$

**Step 2:** Now, let us focus on the behavior on the boundary of the domain. Let the linear operator $\gamma$ be the trace operator on the boundary $(0, T) \times \partial\Omega$. For any $\epsilon \in (0, \epsilon_0)$ small enough, we have $\gamma(p^\epsilon) = p^\epsilon|_{(0,T)\times\partial\Omega}$, then by the trace theorem, one has

$$ ||\gamma(p^{\epsilon})||_{L^2(0,T;L^2(\partial\Omega))} \le C ||p^{\epsilon}||_{L^2(0,T;H^1(\Omega))} $$

where the constant $C$ only depends on $\Omega$. Then

$$ ||\gamma(p^{\epsilon})||_{L^2(0,T;L^2(\partial\Omega))}^2 \le C^2 \int_0^T \int_{\Omega} |p^{\epsilon}|^2 dx dt + C^2 \int_0^T \int_{\Omega} |\nabla p^{\epsilon}(t, \cdot)|^2 dx dt < \infty, $$

due to Lemma A.1 and A.2. Hence, we can deduce that $\gamma(p^\epsilon)$ is weakly convergent in $L^2((0, T) \times \partial\Omega)$. Let $\gamma^0 := \lim_{\epsilon \to 0} \gamma(p^\epsilon)$. For any function $\psi \in C^1(\Omega)$, and for $i = 1, \dots, d$, by Green's formula one has

$$ \int_{\Omega} \partial_i p^{\epsilon} \psi dx = - \int_{\Omega} p^{\epsilon} \partial_i \psi + \int_{\partial \Omega} \psi \gamma(p^{\epsilon}) \nu_i dS. $$

Since $p^\epsilon$ converges weakly to $p^0$ in $H^1(\Omega)$, when $\epsilon \to 0$, one has

$$ \int_{\Omega} \partial_i p^0 \psi dx = - \int_{\Omega} p^0 \partial_i \psi + \int_{\partial \Omega} \psi \gamma^0 \nu_i dS. $$

We can deduce that $\gamma^0 = \gamma(p^0)$.

**Step 3:** We pass to the limit in the weak formulation of (35), for any test function $\psi$ such that $\psi \in C^2([0, T] \times \bar{\Omega})$, $\psi(T, \cdot) = 0$ in $\Omega$, one has

$$ \begin{aligned} & -\underbrace{\int_0^T \int_\Omega p^\epsilon \partial_t \psi dx dt}_{\text{strong convergence}} + A \underbrace{\int_0^T \int_\Omega \nabla p^\epsilon \cdot \nabla \psi dx dt}_{\text{weak convergence}} = \underbrace{\int_\Omega p^{\text{init},\epsilon} \psi(0, \cdot) dx}_{\text{weak convergence}} \\ & -2\epsilon A \underbrace{\int_0^T \int_\Omega \frac{\psi}{1 - \epsilon n^\epsilon} \nabla p^\epsilon \nabla n^\epsilon dx dt}_{\text{bounded as } \epsilon \to 0} + \underbrace{\int_0^T \int_\Omega \psi p^\epsilon (1 - p^\epsilon) (F_1 - F_2)(n^\epsilon, p^\epsilon) dx dt}_{\text{strong convergence}} \\ & -DA \underbrace{\int_0^T \int_{\partial\Omega} (p^\epsilon - p^{\text{ext},\epsilon}) \frac{1 - \epsilon n^{\text{ext},\epsilon}}{1 - \epsilon n^\epsilon} dS}_{\text{weak convergence}} \end{aligned} $$

The weak convergence of the last term on the boundary is obtained from Lemma A.1 and Assumption A.5. When $\epsilon < \epsilon_0$, we have $n^{\text{ext},\epsilon}$, $n^\epsilon$ are uniformly bounded on $(0, T) \times \Omega$ with respect to $\epsilon$, then $\frac{1-\epsilon n^{\text{ext},\epsilon}}{1-\epsilon n^\epsilon}$ converges strongly to 1 when $\epsilon \to 0$. From the previous step, one has $p^\epsilon|_{\partial\Omega} = \gamma(p^\epsilon) \to \gamma(p^0)$ weakly in $L^2((0, T) \times \partial\Omega)$. Passing to the limit, we obtain that $p^0 \in L^2(0, T; H^1(\Omega))$ is a weak solution of the following problem

$$ \left\{ \begin{array}{ll} \partial_t p^0 - A \Delta p^0 = p^0(1-p^0)(F_1 - F_2)(n^0, p^0) & \text{in } (0,T) \times \Omega, \\ p^0(0, \cdot) = p^{\text{init}} & \text{in } \Omega, \\ \frac{\partial p^0}{\partial \nu} = -D(p^0 - p^{\text{ext}}) & \text{on } (0,T) \times \partial \Omega. \end{array} \right. $$
Using (37), we can deduce that this problem is a self-contained initial-boundary-value problem. Moreover, since 0 and 1 are respectively sub- and super-solutions of this problem, it admits a unique classical solution with values in $[0, 1]$. Hence, all the extracted sub-sequences converge to the same limit $p^0$ and $p^0|_{∂Ω} = γ(p^0)$. □

## Acknowledgments

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

This work has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 945322.

## Conflict of interests

There is no conflict of interest.

## References

[1] BARTON, N. H., AND TURELLI, M. Spatial waves of advance with bistable dynamics: cytoplasmic and genetic analogues of allee effects. *The American Naturalist* **178**, 3 (Sept. 2011), 48–75.

[2] CHAN, M. H. T., AND KIM, P. S. Modelling a wolbachia invasion using a slow–fast dispersal reaction–diffusion approach. *Bull Math Biol* **75**, 9 (June 2013).

[3] DANERS, D. Robin boundary value problems on arbitrary domains. *Trans. Amer. Math. Soc.* **352**, 9 (Mar. 2000), 4207–4236.

[4] FIFE, P. C. *Mathematical Aspects of Reacting and Diffusing Systems*, 1st ed. Springer, Berlin, Heidelberg, 1979.

[5] FOCKS, D. A., HAILE, D. G., DANIELS, E., AND MOUNT, G. A. Dynamic life table model for aedes aegypti (diptera: Culicidae): analysis of the literature and model development. *Journal of medical entomology* **30**, 6 (1993), 1003–1017.

[6] GODDARD, J., AND SHIVAJI, R. Stability analysis for positive solutions for classes of semilinear elliptic boundary-value problems with nonlinear boundary conditions. *Proceedings of the Royal Society of Edinburgh: Section A Mathematics* **147**, 5 (2017), 1019–1040.

[7] GORDON, P. V., KO, E., AND SHIVAJI, R. Multiplicity and uniqueness of positive solutions for elliptic equations with nonlinear boundary conditions arising in a theory of thermal explosion. *Nonlinear Analysis: Real World Applications* **15** (Jan. 2014), 51–57.

[8] HILHORST, D., IIDA, M., MIMURA, M., AND NINOMIYA, H. Relative compactness in lp of solutions of some 2m components competition-diffusion systems. *Discrete and Continuous Dynamical Systems* **21**, 1 (May 2008), 233–244.

[9] HILHORST, D., MARTIN, S., AND MIMURA, M. Singular limit of a competition-diffusion system with large interspecific interaction. *Journal of Mathematical Analysis and Applications* **390** (06 2012).

[10] HOFFMANN, A. A., MONTGOMERY, B. L., POPOVICI, J., ITURBE-ORMAETXE, I., JOHNSON, P. H., MUZZI, F., GREENFIELD, M., DURKAN, M., LEONG, Y. S., DONG, Y., COOK, H., AXFORD, J., CALLAHAN, A. G., KENNY, N., OMODEI, C., McGRAW, E. A., RYAN, P. A., RITCHIE, S. A., TURELLI, M., AND O'NEILL, S. L. Successful establishment of wolbachia in aedes populations to suppress dengue transmission. *Nature* **476**, 7361 (06 2011), 454–457.

[11] KORMAN, P. Exact multiplicity of solutions for a class of semilinear Neumann problems. *Communications on Applied Nonlinear Analysis* **9** (Jan. 2002).
[12] KORMAN, P. *Chapter 6 Global Solution Branches and Exact Multiplicity of Solutions for Two Point Boundary Value Problems*, vol. 3. North-Holland, Jan. 2006.

[13] KORMAN, P., LI, Y., AND OUYANG, T. Exact multiplicity results for boundary value problems with nonlinearities generalising cubic. *Proceedings of the Royal Society of Edinburgh Section A: Mathematics* 126, 3 (1996), 599–616. Publisher: Royal Society of Edinburgh Scotland Foundation.

[14] KORMAN, P., LI, Y., AND OUYANG, T. An Exact Multiplicity Result for a class of Semilinear Equations. *Communications in Partial Differential Equations* 22, 3-4 (Jan. 1997), 661–684.

[15] LIONS, P. On the Existence of Positive Solutions of Semilinear Elliptic Equations. *SIAM Review* 24 (1982), 441–467.

[16] MURRAY, J. D. *Mathematical Biology II: Spatial Models and Biomedical Applications*. Springer Science & Business Media, Feb. 2011.

[17] OTERO, M., SCHWEIGMANN, N., AND SOLARI, H. G. A stochastic spatial dynamical model for aedes aegypti. *Bulletin of mathematical biology* 70, 5 (2008), 1297–1325.

[18] OUYANG, T., AND SHI, J. Exact Multiplicity of Positive Solutions for a Class of Semilinear Problems. *Journal of Differential Equations* 146, 1 (June 1998), 121–156.

[19] PAO, C. V. *Nonlinear Parabolic and Elliptic Equations*, 1st ed. Springer, Boston, MA, 1992.

[20] PERTHAME, B. *Parabolic equations in biology*. Lecture Notes on Mathematical Modelling in the Life Sciences. Springer, 2015.

[21] SCHAAF, R. Global behaviour of solution branches for some Neumann problems depending on one or several parameters. *Journal für die reine und angewandte Mathematik* 1984, 346 (Jan. 1984), 1–31.

[22] SHI, S., AND LI, S. Existence of solutions for a class of semilinear elliptic equations with the Robin boundary value condition. *Nonlinear Analysis: Theory, Methods & Applications* 71, 7 (Oct. 2009), 3292–3298.

[23] SIMON, J. Compact sets in the space $l^p(o, t; b)$. *Annali di Matematica pura ed applicata* (Sept. 1986), 65–96.

[24] SMOLLER, J. *Shock waves and reaction–diffusion equations*, vol. 258. Springer Science & Business Media, 2012.

[25] SMOLLER, J., AND WASSERMAN, A. Global bifurcation of steady-state solutions. *Journal of Differential Equations* 39, 2 (1981), 269–290.

[26] STRUGAREK, M., AND VAUCHELET, N. Reduction to a single closed equation for 2 by 2 reaction-diffusion systems of Lotka-Volterra type. *SIAM Journal on Applied Mathematics* 76, 5 (2016), 2060–2080.

[27] TSAI, C.-C., WANG, S.-H., AND HUANG, S.-Y. Classification and evolution of bifurcation curves for a one-dimensional Neumann–Robin problem and its applications. *Electron. J. Qual. Theory Differ. Equ.*, 85 (2018), 1–30.

[28] WANG, S.-H. A correction for a paper by j. smoller and a. wasserman. *Journal of Differential Equations* 77, 1 (1989), 199–202.

[29] WANG, S.-H., AND KAZARINOFF, N. D. Bifurcation of steady-state solutions of a scalar reaction-diffusion equation in one space variable. *Journal of the Australian Mathematical Society* 52, 3 (June 1992), 343–355.

[30] WANG, Y., AND SHI, J. Persistence and Extinction of Population in Reaction-Diffusion-Advection Model with Weak Allee Effect Growth. *SIAM J. Appl. Math.* 79, 4 (Jan. 2019), 1293–1313. Publisher: Society for Industrial and Applied Mathematics.
[31] WANG, Y., SHI, J., AND WANG, J. Persistence and extinction of population in reaction-diffusion-advection model with strong Allee effect growth. *J. Math. Biol.* **78**, 7 (June 2019), 2093–2140.

[32] WERREN, J., BALDO, L., AND CLARK, M. *Wolbachia*: master manipulators of invertebrate biology. *Nat Rev Microbiol* **6** (Oct. 2008), 741–751.

[33] ZHANG, J., LI, S., AND XUE, X. Multiple solutions for a class of semilinear elliptic problems with Robin boundary condition. *Journal of Mathematical Analysis and Applications* **388**, 1 (Apr. 2012), 435–442.