# Travelling waves in a PDE–ODE coupled system with nonlinear diffusion

K. Mitra111email<EMAIL_ADDRESS>J.M. Hughes  University of British
Columbia, Canada  S. Sonner  H.J. Eberl  University of Guelph, Canada  J.D.
Dockery  Montana State University, USA

###### Abstract

We analyze travelling wave (TW) solutions for nonlinear systems consisting of
an ODE coupled to a degenerate PDE with a diffusion coefficient that vanishes
as the solution tends to zero and blows up as it approaches its maximum value.
Stable TW solutions for such systems have previously been observed numerically
as well as in biological experiments on the growth of cellulolytic biofilms.
In this work, we provide an analytical justification for these observations
and prove existence and stability results for TW solutions of such models.

Using the TW ansatz and a first integral, the system is reduced to an
autonomous dynamical system with two unknowns. Analysing the system in the
corresponding phase–plane, the existence of a unique TW is shown, which
possesses a sharp front and a diffusive tail, and is moving with a constant
speed. The linear stability of the TW in two space dimensions is proven under
suitable assumptions on the initial data. Finally, numerical simulations are
presented that affirm the theoretical predictions on the existence, stability
and parametric dependence of the TWs.

Keywords: travelling waves, stability, degenerate diffusion, PDE–ODE system,
biofilm

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries and main result
    1. 2.1 Auxiliary quantities and reductions
    2. 2.2 Main theorem
  3. 3 The existence of travelling waves
    1. 3.1 The phase–plane
    2. 3.2 The orbit connecting with $(0,1)$
    3. 3.3 Uniqueness of the orbit connecting with $(0,1)$ at $\xi=0$
    4. 3.4 Ordering of the orbits with respect to $v$
    5. 3.5 Existence/non-existence of travelling waves
  4. 4 Linear stability of the travelling waves in two space dimensions
    1. 4.1 Asymtotic expansion and linearization
    2. 4.2 Stability in $L^{1}$-norm
  5. 5 Numerical results
    1. 5.1 PDE simulations
      1. 5.1.1 One dimensional results: transience and stability
      2. 5.1.2 Two dimensional results: transience and stability
    2. 5.2 ODE simulations
      1. 5.2.1 TW algorithm
      2. 5.2.2 Validation of the TW algorithm using the PDE scheme
      3. 5.2.3 Parametric study in the limiting cases
    3. 5.3 Numerical continuation method
  6. 6 Interpretation of the results
  7. A Grid independence study for the PDE simulations

##  1 Introduction

Motivated by mathematical models for biofilm growth, we investigate travelling
wave (TW) solutions for coupled PDE–ODE systems of the form

$\displaystyle\partial_{t}M=\partial_{x}[D(M)\,\partial_{x}M]+\left(f(S)-\lambda\right)M,$
(1.1a) $\displaystyle\partial_{t}S=-\gamma\,f(S)\,M,$ (1.1b)

where $x\in{\mathbb{R}}$ and $t>0$ represent the space and time coordinates
respectively. The functions $M$ and $S$ are normalized, $S$ takes values in
the interval $[0,\infty)$ and $M$ in $[0,1)$. The growth and decay
characteristics of the system are represented by the constants
$\lambda,\,\gamma\in(0,1)$. The diffusion coefficient $D:[0,1)\to[0,\infty)$
has a singularity as well as a degeneracy. More specifically, it satisfies
$D(m)\nearrow\infty$ as $m\nearrow 1$, and $D(0)=0$. The source function
$f:[0,\infty)\to[0,1]$ is Lipschitz continuous and increasing. The exact
assumptions on $D$ and $f$ are stated in Section 2.

For $M$ we consider the following boundary conditions

$\displaystyle M(\pm\infty,t)=0,\quad(D(M)\,\partial_{x}M)(\pm\infty,t)=0,$
(1.2a) and initial conditions that are consistent with the above. In
particular, we assume that $\displaystyle\begin{cases}M(x,0)\in[0,1),\quad
S(x,0)\in[0,\infty)&\text{ for all }x<0,\\\ M(x,0)=0,\quad S(x,0)=1&\text{ for
all }x\geq 0.\end{cases}$ (1.2b)

Standard approaches to prove the existence and stability of TW solutions for
equations such as the Fischer-KPP equation do not generalize to the system
(1.1). Difficulties arise through the degeneracy and singularity of the
diffusion coefficient, and the nonlinear coupling between the PDE and the ODE,
which leads to a non-monotone profile for $M$. In this paper, we use ordering
of orbits in the $(M,S)$ phase–plane to prove the existence and uniqueness of
a TW solution for the system (1.1) satisfying suitable boundary conditions. We
further derive stability results using asymptotic expansions. Moreover, we
present numerical simulations that affirm the theoretical predictions on the
existence, stability and parametric dependence of the TWs.

Several examples of semilinear evolution equations coupled to an ODE through
the source–term can be found in [18] as models of reactive transport through
the subsurface. In mathematical biology, examples of similar systems are found
in [23, Chapter 13], as well as a discussion on their TW solutions. PDE–ODE
coupled systems of various type with nonlinear diffusion coefficients are used
to model variety of other physical or biological processes, ranging from
hysteretic flow through porous media [19], to tumor growth [15], and wound
healing [16]. The particular motivation for our analysis is the model for
cellulolytic biofilm growth in [11] for which TW solutions have been observed
numerically. Cellulolytic biofilms play an important role in the production of
cellulosic ethanol, a renewable biofuel that can be implemented in the
existing transportation infrastructure. In contrast to more traditional
biofilms that form on mostly abiotic surfaces and develop colonies that grow
into the surrounding aqueous phase, many cellulolytic biofilms degrade and
consume the biological material that supports them and form crater like
depressions, a phenomenon known as inverted colony formation. Since the
nutrients are immobile, whereas, the biofilm expands spatially, a PDE–ODE
model was proposed for cellulolytic biofilms in [11] which is a special case
of (1.1). In the model, $M$ represented the biomass density and $S$ the
nutrient concentration. The bacteria consume nutrients and degrade the
biological material which results in the production of biomass. This was
described in [11] using the Monod reaction function $f$, whereas, the spatial
spreading of biomass was modelled by a density-dependent diffusion coefficient
$D$, with the corresponding expressions being

$\displaystyle D(m)$ $\displaystyle=\frac{\delta\,m^{a}}{(1-m)^{b}},$
$\displaystyle\delta>0,\ a,\,b>1,$ (1.3a) $\displaystyle f(s)$
$\displaystyle=\frac{s}{\kappa+s},$ $\displaystyle 0<\kappa\ll 1.$ (1.3b)

Cellulolytic biofilm formation has been studied both numerically and
experimentally. In [31], an agent based stochastic discrete cellular automaton
model was used to study the system. The model in [11], on the other hand, is
the deterministic continuum model (1.1) with the specific functions $f$ and
$D$ in (1.3). Extensions of this model that account for attachment of cells
from the aqueous phase to the biofilm were presented in [24, 17] using either
Itô stochastic differential equations or random differential equations.
Numerical simulations of the original deterministic model in [11] suggest the
existence of TW solutions, which describe a constant rate of degradation of
the cellulosic material that is utilised by the bacteria. Additional numerical
evidence for this is given in [17] where a different time integration method
is used, along with an independent implementation. Celluloytic biofilm systems
are very difficult to observe experimentally with time-lapse microscopy
techniques. Nevertheless, in [31] experiments were reported that suggest
degradation of paper chads by cellulosic biofilms at a constant speed, which
gives indirect evidence for TW like degradation in the biological system.
Furthermore, experimental observations in [31] and in [9] indicate that the
width of the microbially active band remains constant as the wave of microbial
crater formation spreads. A rigorous proof of the existence and features of TW
solutions of the model in [11], or an answer to the question ‘under which
conditions on parameters such TWs can be found’, have so far been open
problems. In this study we provide the answers. For this purpose, the initial
and boundary conditions (1.2) are chosen to be consistent with the physical
setting of the numerical experiments in [11]. We emphasize that, although our
study is motivated and prompted by the cellulosic biofilm system, our results
are valid for a significantly wider class of problems.

The existence of TWs for the (scalar) porous medium equation with a nonlinear
source term was investigated in [5]. It was shown that under certain
conditions on the coefficients, there exists a minimum speed for which TW
solutions exist. The results were extended to include nonlinear advection
terms in [3]. Furthermore, in [1] the stability of such TW solutions was shown
in one space dimension. For the porous medium equation with Fischer type
reaction term, the existence of TWs was shown in [25], again for wave-speeds
larger than a minimum value. In [4], further qualitative properties of the TW
solutions were analysed. These results were generalized in [12] for equations
with the biofilm diffusion coefficient $D$ in (1.3a). However, the
aforementioned results are limited to scalar equations with Fisher-type
nonlinear reaction term, and thus, exclude the complex interplay between the
ODE and the PDE solutions. While the TW profile is a monotone function with
respect to $x$ in all the mentioned results, this is not the case for the TW
profile of $M$ in our system.

TWs for PDE–ODE coupled systems have been studied for multiphase flow through
porous media in [28, 22, 20], where non-monotone profiles of $M$ have been
observed. The ordering of orbits in the phase–plane is also used in these
papers to predict the behaviour of the TWs. The existence of TWs in two-
dimensions for a PDE–ODE model was investigated in [21] in the context of
hysteretic flow through porous media, and non-planar TWs were shown to exist.
However, TWs in these cases originate from the advection term, rather than the
source term. Non-monotone profiles have also been observed for the TWs of
PDE–ODE systems arising in biology, see [15, 16] for examples. TWs for a
PDE–PDE coupled model of bacteria spreading in an aqueous phase were analysed
in [26]. Nevertheless, these systems differ fundamentally in their structure
from (1.1). In our setting, as will be observed later, the TWs have
distinctive features which distinguish them from the examples above. They
inherit a sharp front and a minimum speed of propagation like TWs of the
porous medium equation. However, due to the coupling with the ODE, they
exhibit a non-monotone profile with a diffusive tail.

The outline of our paper is as follows: In Section 2, we state the assumptions
on the associated functions, and using the TW ansatz, the system (1.1) is
reduced to a dynamical system with two unknowns. The existence result for the
TWs is also stated, see Theorem 2.1. In Section 3, using phase–plane analysis
we develop the auxiliary results which are then used to prove Theorem 2.1.
Section 4 is dedicated to proving a linear stability result for the TWs in two
space dimensions using asymptotic expansion. In Section 5, numerical results
are presented for a discretization of the full PDE–ODE system, an ODE approach
inspired by the TW analysis, and a numerical continuation approach. All three
different approaches concur about the existence/non-existence of TWs in a
parametric regime indicated by our theory. Furthermore, numerical results
showing the influence of the parameters, and the stability of the TWs are
presented. In Section 6, we interpret the analytical and numerical results in
the context of cellulolytic biofilms, and discuss possible generalizations and
future applications.

##  2 Preliminaries and main result

We investigate the existence of solutions for the system (1.1) with the
boundary and initial conditions (1.2). To this end, we assume that $D$ and $f$
have the following properties:

  1. (P1)

The diffusion coefficient $D:[0,1)\to[0,\infty)$ is an increasing function in
$C^{1}([0,1))$ which satisfies for constants $a,\,b>1$

$\lim\limits_{m\searrow 0}\frac{D(m)}{m^{a}}\in{\mathbb{R}}^{+},\;\;\text{ and
}\;\;\lim\limits_{m\nearrow 1}(1-m)^{b}\,D(m)\in{\mathbb{R}}^{+}.$

  2. (P2)

The source function $f:[0,1]\to[0,1]$ is an increasing function in
$C^{1}([0,1])$ which satisfies for a constant $\kappa\in(0,1]$,

$f(0)=0,\quad f^{\prime}(0)=\frac{1}{\kappa},\;\;\text{ and
}\;\;f(1)\in(\lambda,1].$

The TW ansatz for the subsequent theory is:

######  Assumption 2.1 (TW ansatz).

For a wave-speed $v>0$, and the travelling wave coordinate $\xi=x-vt$, there
exist $M,S:{\mathbb{R}}\times{\mathbb{R}}^{+}\to[0,1]$ which satisfy
(1.1)–(1.2) and

$M(x,t)=M(\xi),\quad S(x,t)=S(\xi).$

With this ansatz, the system (1.1) is written as

$\displaystyle-$ $\displaystyle
v\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=\tfrac{\mathrm{d}}{\mathrm{d}\xi}\left[D(M)\,\tfrac{\mathrm{d}}{\mathrm{d}\xi}M\right]+\left(f(S)-\lambda\right)M,$
(2.1a) $\displaystyle v\tfrac{\mathrm{d}}{\mathrm{d}\xi}S=\gamma\,f(S)\,M.$
(2.1b)

The initial and boundary conditions (1.2) are transformed into

$\displaystyle M(-\infty)$
$\displaystyle=(D(M)\,\tfrac{\mathrm{d}}{\mathrm{d}\xi}M)(-\infty)=0,$ (2.2a)
$\displaystyle M(\xi)$
$\displaystyle=(D(M)\,\tfrac{\mathrm{d}}{\mathrm{d}\xi}M)(\xi)=0,\quad
S(\xi)=1\text{ for all }\xi\geq 0.$ (2.2b)

###  2.1 Auxiliary quantities and reductions

For a given solution $(M,S)$ of (2.1)–(2.2), the accumulated biomass
$\omega:{\mathbb{R}}\to{\mathbb{R}}^{+}$ until $\xi\in{\mathbb{R}}$, is
defined as

$\displaystyle\omega(\xi):=\int_{-\infty}^{\xi}M,\quad\text{ implying
}\quad\tfrac{\mathrm{d}}{\mathrm{d}\xi}\omega=M\text{ and }\omega(-\infty)=0.$
(2.3)

Assuming $M$, $S$ and $\omega$ are smooth enough, (2.1a) is then rewritten
using (2.1b) and (2.3) as

$-v\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=\tfrac{\mathrm{d}}{\mathrm{d}\xi}\left[D(M)\,\tfrac{\mathrm{d}}{\mathrm{d}\xi}M\right]+\tfrac{\mathrm{d}}{\mathrm{d}\xi}\left(\tfrac{v}{\gamma}S-\lambda\omega\right).$

Integrating the above equation from $-\infty$ to $\xi$, we have using (2.2)
that

$-vM=D(M)\,\tfrac{\mathrm{d}}{\mathrm{d}\xi}M+\tfrac{v}{\gamma}(S-S(-\infty))-\lambda\omega.$

Observe that $S$ is a non-decreasing function by (2.1b) with $S(0)=1$, and
$S(\xi)\geq 0$ for all $\xi\in{\mathbb{R}}$. Therefore, $S(-\infty)\geq 0$ is
well-defined. Upon rearranging the above equation one has

$\displaystyle
D(M)\,\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=\lambda\omega-v\,(M+\tfrac{1}{\gamma}(S-S(-\infty))).$
(2.4)

Passing $\xi\to\infty$ and using (2.2) one further has

$\displaystyle\omega(+\infty)=\tfrac{v}{\lambda\,\gamma}(1-S(-\infty)),$ (2.5)

which serves as a kind of Rankine-Hugoniot condition for the wave-speed $v$.
Finally, using the relations above, (2.1) is rewritten as an autonomous
dynamical system for $M$, $S$ and $\omega$,

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=\tfrac{1}{D(M)}[\lambda\omega-v\,(M+\tfrac{1}{\gamma}(S-S(-\infty)))],$
(2.6a)
$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\xi}S=\tfrac{\gamma}{v}\,f(S)\,M,$
(2.6b) $\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\xi}\omega=M.$ (2.6c)

Observing that $S(\xi)\geq S(-\infty)\geq 0$ for all $\xi\in{\mathbb{R}}$,
equation (2.6b) is rewritten using (2.6c) as

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\xi}S=\tfrac{\gamma}{v}f(S)\tfrac{\mathrm{d}}{\mathrm{d}\xi}\omega\quad\text{
or }\quad\tfrac{v}{\gamma
f(S)}\tfrac{\mathrm{d}}{\mathrm{d}\xi}S=\tfrac{\mathrm{d}}{\mathrm{d}\xi}\omega.$
(2.7)

We introduce the function $F:(0,1]\to{\mathbb{R}}^{+}$ as

$\displaystyle
F(s):=\int^{1}_{s}\frac{\mathrm{d}\varrho}{f(\varrho)}.\quad\text{ It follows
from \ref{prop:f} that }\;F^{\prime}(s)<0,\;\lim\limits_{s\searrow
0}F(s)=\infty,\;F(1)=0.$ (2.8)

The limit $F(s)\to\infty$ for $s\searrow 0$ follows from (P2) since $f(s)\sim
s/\kappa$ in a right neighbourhood of $s=0$, and consequently
$F(s)=\int_{s}\frac{1}{f}\sim-\kappa\log(s)\to\infty$ as $s\searrow 0$.
Integrating (2.7) from $\xi\leq 0$ to $+\infty$ and using (2.2) one has
$\tfrac{v}{\gamma}F(S)=\omega(+\infty)-\omega$, which upon rearranging and
using (2.5) gives

$\displaystyle\omega=\tfrac{v}{\gamma}\left[\tfrac{1}{\lambda}(1-S(-\infty))-F(S)\right].$
(2.9)

Passing $\xi\to-\infty$ in the above equation, using $\omega(-\infty)=0$ from
(2.3) and cancelling equal terms, we get

$S(-\infty)+\lambda F(S(-\infty))=1.$

Since $F$ is strictly decreasing and convex as evident from (P2), there can at
most be two solutions of the equation $g(s)=s+\lambda F(s)=1$ in $(0,1]$. One
trivial solution is $s=1$. Since $g(0)=\infty$ from (2.8), the existence of
the second solution is guaranteed if $g^{\prime}(1)>0$, or
$F^{\prime}(1)=-1/f(1)>-1/\lambda$ which holds due to (P2). Hence, we define
$s_{-\infty}\in(0,1)$ as the solution of

$\displaystyle g(s_{-\infty})=s_{-\infty}+\lambda F(s_{-\infty})=1.$ (2.10)

######  Remark 2.1 (The value of $s_{-\infty}\in(0,1)$).

Let $f$ be given by the expression in (1.3). Then, $F(s)=1-s-\kappa\log(s)$.
Thus, for $\kappa\ll 1$ and $\lambda\sim 1$, one has

$s_{-\infty}\approx\exp\left(-\tfrac{1-\lambda}{\kappa\lambda}\right).$

For the parameters $\lambda=0.42$, $\gamma=0.4$ and $\kappa=0.01$ used in
[11], we estimate that

$s_{-\infty}\approx 10^{-60},$

which is negligible for all practical purposes. This suggests that in this
parameter regime, the substrate is fully depleted after the TW has passed. In
Section 5.3, we provide an example of parameters for which $s_{-\infty}\approx
0.11$, and hence, there remains a significant level of residual substrates.

Finally, substituting (2.9) into (2.6) we get the reduced autonomous dynamical
system with only two unknowns $M$ and $S$, i.e.,

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=\tfrac{v}{\gamma\,D(M)}[(1-S-\lambda
F(S))-\gamma M],$ (2.11a)
$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\xi}S=\tfrac{\gamma}{v}\,f(S)\,M.$
(2.11b)

The revised boundary conditions for this system are

$\displaystyle M(-\infty)=0,\;S(-\infty)=s_{-\infty}\text{ and
}M(\xi)=0,\;S(\xi)=1\text{ for all }\xi\geq 0.$ (2.12)

This will be the main system analysed in this paper.

######  Remark 2.2 (The flux conditions at $\xi=0$ and $\xi=-\infty$).

Observe that any solution $(M,S)$ of the system (2.11)–(2.12), automatically
satisfies the boundary condition for the flux, i.e.,
$D(M)\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=0$ at $\xi=0$ and $\xi=-\infty$.

###  2.2 Main theorem

For the rest of this study, we focus on the following parametric regime: for
$F$ defined in (2.8), and $g=\mathbb{I}+\lambda F$, let
$\lambda,\gamma\in(0,1)$ be such that

$\displaystyle g(f^{-1}(\lambda))=f^{-1}(\lambda)+\lambda
F(f^{-1}(\lambda))<1-\gamma.$ (2.13)

Observe that, the function $g(y)$ takes its minimum value $g_{\min}$ at
$y=f^{-1}(\lambda)$. The condition (P2) then guarantees that
$f^{-1}(\lambda)\in(0,1)$. The existence of $s_{-\infty}\in(0,1)$ in (2.10)
proves that $g_{\min}<1$. Assumption (2.13) enforces a stronger condition,
i.e., that $g_{\min}<1-\gamma$.

For the function $f(s)=s/(\kappa+s)$ the condition can be stated in a more
compact form as

$\displaystyle 0<\gamma+\lambda+\kappa\lambda\,(1-\log(\kappa\lambda))\leq 1.$
(2.14)

We introduce the following important integral:

$\displaystyle\mathcal{G}(s)$ $\displaystyle:=\int_{s}^{1}(\varrho+\lambda
F(\varrho)-(1-\gamma))\,\frac{\mathrm{d}\varrho}{f(\varrho)}$ (2.15a)
$\displaystyle\overset{\eqref{eq:DefF}}{=}\int_{s}^{1}\frac{\varrho}{f(\varrho)}\,\mathrm{d}\varrho+\frac{\lambda}{2}\,F^{2}(s)-(1-\gamma)\,F(s).$
(2.15b)

This representation follows using $F^{\prime}=1/f$ and $F(1)=0$ from (2.8).
The shape of the $\mathcal{G}$-integral is shown in Figure 1. Observe that,
$\mathcal{G}(1)=0$ and $\mathcal{G}(s)$ is a decreasing function for
$s>s^{*}$, where $s=s^{*}$ solves $s+\lambda F(s)=1-\gamma$. Hence,
$\mathcal{G}>0$ in an interval $(s_{g},1)$ where $s_{g}\in(0,s^{*})$. Also
from (2.15b), $\mathcal{G}(s)\to+\infty$ as $s\searrow 0$ since
$F(s)\to\infty$ in this case, see (2.8). Hence, depending on the parameter
values, $\mathcal{G}$ might or might not have a negative part. This has a
profound effect on the existence of TWs as stated in our main theorem below.

$s$$\mathcal{G}$$0$$\lambda_{1}=\gamma_{1}=.4$$\lambda_{2}=\gamma_{2}=.3$$s^{*}_{2}$$s^{*}_{1}$
Figure 1: The plot of $\mathcal{G}$ as defined in (2.15) for
$f(s)=s/(\kappa+s)$, $\lambda_{1}=\gamma_{1}=0.4$ and
$\lambda_{2}=\gamma_{2}=0.3$. The points $s=s^{*}\in(0,1)$ marked are the
solutions of $s+\lambda F(s)=1-\gamma$.

######  Theorem 2.1 (Existence of the TW solution).

Assume (P1)–(P2). Let (2.13) be satisfied and let $\mathcal{G}(s)>0$ for all
$s\in(s_{-\infty},1)$. Then there exists a unique $v>0$ such that a travelling
wave solution $(M,S):{\mathbb{R}}\to[0,1]^{2}$ with
$D(M)\tfrac{\mathrm{d}}{\mathrm{d}\xi}M\in C({\mathbb{R}})$ and $S\in
C^{1}({\mathbb{R}})$ exists satisfying (2.11)–(2.12).

######  Remark 2.3 (Conditions on existence).

The condition (2.13), used in Theorem 2.1, provides upper bounds for $\lambda$
and $\gamma$, whereas, the condition $\mathcal{G}(s)>0$ for
$s\in(s_{-\infty},1)$ provides lower bounds for $\lambda$ and $\gamma$ for the
existence of the TW solutions. We show in Proposition 3.6 that the latter
condition is also a necessary condition.

##  3 The existence of travelling waves

In this section we prove Theorem 2.1 by analysing the dynamical system (2.11).

###  3.1 The phase–plane

For a given orbit $\xi\mapsto(M,S)$ satisfying (2.11), the scaled TW
coordinate $\tau$ is defined by the coordinate transform

$\displaystyle\tau(\xi):=\int_{0}^{\xi}\tfrac{\mathrm{d}\varrho}{D(M(\varrho))},\quad\text{
implying }\tfrac{\mathrm{d}}{\mathrm{d}\xi}\tau=\tfrac{1}{D(M(\xi))}.$ (3.1)

Moreover, to shorten notation, we introduce

$\displaystyle\ell(s;\lambda,\gamma):=\frac{1}{\gamma}[(1-s)-\lambda F(s)],$
(3.2)

which from (P2) and (2.8) has the properties

$\displaystyle\ell\in C^{1}((0,1]),\;\ell^{\prime}(s)>0\text{ for
}s<s_{\mathrm{M}}:=f^{-1}(\lambda),\text{ and }\ell^{\prime}(s)<0\text{ for
}s>s_{\mathrm{M}};$ (3.3a) $\displaystyle\ell(1)=0,\text{ and
}\lim\limits_{s\searrow 0}\ell(s)=-\infty.$ (3.3b)

Observe that, in terms of the function $g$ introduced in Section 2,
$\ell(s)=\gamma^{-1}[1-g(s)]$ (recall that $g=\mathbb{I}+\lambda\,F$). In the
following sections, we will only use the properties (3.3) of $\ell$, and will
not further use $F$ or $g$, to keep the notation as clear as possible. From
(2.15a), we additionally have that

$\displaystyle\mathcal{G}(s)=\gamma\int_{s}^{1}\left(\frac{1-\ell(\varrho)}{f(\varrho)}\right)\mathrm{d}\varrho.$
(3.4)

With the coordinate transform (3.1) and the definition (3.2), system (2.11) is
re-written as

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\tau}M=v\,[\ell(S)-M],$ (3.5a)
$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\tau}S=\tfrac{\gamma}{v}\,f(S)\,MD(M).$
(3.5b)

$m$$s$$0$$1$$1$$s^{*}$$s_{*}$$m=\ell(s)$$(M,S)$

$m$$s$$0$$1$$(M,S)$$s_{-\infty}$$m=\ell(s)$

Figure 2: (left) The direction of orbits of the dynamical system (2.11) in the
phase–plane $[0,1)\times(0,1]$. (right) A zoomed view into the phase–plane
near the equilibrium point $(0,s_{-\infty})$. An orbit $(M,S)$ originating
from $(0,s_{-\infty})$ is also shown.

Since this is an autonomous system, we look into the phase–plane
$[0,1)\times(0,1]$. Henceforth, $(m,s)$ will represent a point in this
phase–plane, and $(M,S)$ will denote the orbits in this plane. The directions
of the orbits are shown in Figure 2. The line $m=\ell(s)$ is highlighted and
represents the nullcline of $M$, i.e., the points where
$\tfrac{\mathrm{d}}{\mathrm{d}\tau}M=0$. Due to the restriction (2.13)
imposed, this nullcline intersects the line $m=1$ at precisely two points
$s_{*}$ and $s^{*}$. This follows from the properties of $\ell$ in (3.3) along
with the observation that

$\ell(s_{\mathrm{M}})=\gamma^{-1}[1-s_{\mathrm{M}}-\lambda
F(s_{\mathrm{M}})]\overset{\eqref{eq:parameter}}{>}1,\quad\text{ where
}s_{\mathrm{M}}=f^{-1}(\lambda).$

Hence, the $s$-coordinates $s_{*},\,s^{*}\in(0,1)$ satisfying

$\displaystyle s_{*}:=\min\\{s\in(0,1):\ell(s)=1\\}<s_{\mathrm{M}},\quad
s^{*}:=\max\\{s\in(0,1):\ell(s)=1\\}>s_{\mathrm{M}},$ (3.6)

exist (consistent with the definition of $s^{*}\in(0,1)$ below (2.15)).
Consequently, we have the following:

######  Lemma 3.1 (The existence of orbits).

Let $v>0$ be given, (2.13) be satisfied and let
$(M_{0},S_{0})\in\mathfrak{R}:=[0,1)\times[s_{-\infty},1].$ Then, there exists
a unique orbit $\tau\mapsto(M,S)\in(C^{1}({\mathbb{R}}))^{2}$, satisfying
(3.5) with $(M,S)(0)=(M_{0},S_{0})$. The equilibrium points of the system
(3.5) are $(0,s_{-\infty})$ and $(0,1)$ with $s_{-\infty}\in(0,1)$ satisfying
(2.10). If $(M_{0},S_{0})$ is not an equilibrium point, then

  1. (i)

for $\tau>0$, the orbit either exits $\mathfrak{R}$ through the line
$\\{s=1\\}$ or ends at $(0,1)$.

  2. (ii)

for $\tau<0$, the orbit enters $\mathfrak{R}$ either though the line segment
$\\{m=0,s\geq s_{-\infty}\\}$ or through $\\{s=s_{-\infty}\\}$.

###### Proof.

Observe that, orbits satisfying (3.5) are locally well-posed at any point in
$[0,1)\times[s_{-\infty},1]$. This follows from the Picard–Lindelöf theorem
since the right hand sides of (3.5) are locally–Lipschitz with respect to $M$
and $S$ for $M<1$ and $S\in[s_{-\infty},1]$.

_(i)_ The direction of the orbits implies that $(M,S)$ can exit $\mathfrak{R}$
through the boundaries $\\{m=1\\}$ or $\\{s=1\\}$. To rule out the line
$\\{m=1\\}$, we assume the contrary, i.e., we suppose there exists
$\tau_{1}>0$ such that $M(\tau_{1})=1$ and $M(\tau)<1$ for all
$0<\tau<\tau_{1}$. Then, $(M,S)$ satisfies (3.5) for all $\tau<\tau_{1}$. It
is straightforward to see that $S(\tau_{1})\in[s_{*},s^{*}]$ since
$\lim_{\tau\nearrow\tau_{1}}\tfrac{\mathrm{d}}{\mathrm{d}\tau}M=v\,[\ell(S(\tau_{1}))-1]$
and the definition of $\tau_{1}$ demands that
$\lim_{\tau\nearrow\tau_{1}}\tfrac{\mathrm{d}}{\mathrm{d}\tau}M\geq 0$. Using
the intermediate value theorem, for a given $\varepsilon<1-M_{0}$, there
exists $\tau_{\varepsilon}\in(0,\tau_{1})$ such that
$M(\tau_{\varepsilon})=1-\varepsilon$. Moreover, observe that there exists a
constant $C>0$ such that

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\tau}M\leq C\text{ in
}(\tau_{\varepsilon},\tau_{1}),\;\text{ or, integrating in
}(\tau_{\varepsilon},\tau_{1}),\quad\tau_{1}-\tau_{\varepsilon}\geq\frac{\varepsilon}{C}.$
(3.7)

From (3.5b) and using (3.7), for some constants $C_{1/2}>0$ independent of
$\varepsilon$, one has

$S(\tau_{1})-S(\tau_{\varepsilon})=\int_{\tau_{\varepsilon}}^{\tau_{1}}\tfrac{\gamma}{v}f(S)MD(M)\mathrm{d}\tau\geq
C_{1}\int_{\tau_{\varepsilon}}^{\tau_{1}}D(1-\varepsilon)\,\mathrm{d}\tau\overset{\ref{prop:D}}{\geq}C_{2}\frac{\tau_{1}-\tau_{\varepsilon}}{\varepsilon^{b}}\overset{\eqref{eq:Mleq1contracdiction}}{\geq}\frac{C_{2}}{C\varepsilon^{b-1}}\to\infty,$

as $\varepsilon\to 0$ (as $b>1$). Since $S(\tau_{\varepsilon})>0$, this
contradicts $S(\tau_{1})\in[s_{*},s^{*}]$, thus proving that $(M,S)$ cannot
exit through $\\{m=1\\}$.

_(ii)_ For $\tau<0$, the fact that the orbits can enter through the mentioned
segments is clear. The orbit cannot enter through the boundary $\\{s=1\\}$
since $S$ is strictly increasing in $\mathfrak{R}$ for $m>0$, and $(0,1)$ is
an equilibrium point. The fact that the orbit cannot enter through $\\{m=1\\}$
follows similarly as the proof of point (i). ∎

######  Remark 3.1 (The travelling wave solutions avoid the degeneracy at
$m=1$).

The proof above shows that, for any TW solution $(M,S)$ satisfying
(2.11)–(2.12), there exists a constant $\varepsilon>0$ such that

$0\leq M(\xi)\leq 1-\varepsilon,\text{ for all }\xi\in{\mathbb{R}}.$

Hence, the degeneracy of the diffusion coefficient $D$, due to the possibility
of $D(M)\to\infty$ as $M\nearrow 1$, is avoided.

###  3.2 The orbit connecting with $(0,1)$

For any orbit $(M,S)$ described in Lemma 3.1, $S$ is strictly increasing for
all $\tau\in{\mathbb{R}}$ provided $(M,S)\in\mathfrak{R}$ and $M>0$. Hence,
for a given $s\in(s_{-\infty},1)$ there can exist at most one
$\tau\in{\mathbb{R}}$ such that $S(\tau)=s$. This allows us to introduce the
unique mapping $S\mapsto M$ through the following function.

######  Definition 1 (The $\mathfrak{M}$–map).

For a given $(M_{0},S_{0})\in\mathfrak{R}$, let $(M,S)$ be the unique orbit
$\tau\mapsto(M,S)\in(C^{1}({\mathbb{R}}))^{2}$, satisfying (3.5) and
$(M,S)(0)=(M_{0},S_{0})$. Then the continuous function
$\mathfrak{M}:[s_{-\infty},1]\to[0,1)$ is defined as

$\displaystyle\mathfrak{M}(s):=\begin{cases}M(\tau)&\text{ if there exists
$\tau\in{\mathbb{R}}$ such that }S(\tau)=s,\\\ 0&\text{
otherwise}.\end{cases}$ (3.8)

Let us introduce the function

$\displaystyle\Phi(m):=\int_{0}^{m}\varrho\,D(\varrho)\,\mathrm{d}\varrho,\text{
such that }\Phi^{\prime}(m)=mD(m)\geq 0.$ (3.9)

Observe that $\mathfrak{M}$ satisfies $\mathfrak{M}(S_{0})=M_{0}$, and for all
$\mathfrak{M}>0$,

$\displaystyle\dfrac{\mathrm{d}\mathfrak{M}}{\mathrm{d}s}=\dfrac{v^{2}}{\gamma}\dfrac{\ell(s)-\mathfrak{M}(s)}{f(s)\,\mathfrak{M}\,D(\mathfrak{M})}.$
(3.10a) Using (3.9), we alternatively rewrite the equation above as
$\displaystyle\dfrac{\mathrm{d}\Phi(\mathfrak{M})}{\mathrm{d}s}=\dfrac{v^{2}}{\gamma
f(s)}[\ell(s)-\mathfrak{M}(s)].$ (3.10b)

Our focus will be on a specific group of maps $\mathfrak{M}$ which originate
from $(\varepsilon,1)$.

######  Lemma 3.2.

For fixed $v,\,\varepsilon>0$ and $(M_{0},S_{0})=(\varepsilon,1)$, let
$\mathfrak{M}^{\varepsilon}$ denote the $\mathfrak{M}$-mapping in the sense of
Definition 1. For $a>1$ introduced in (P1), let
$\underline{\mathfrak{M}}:(0,1]\to[0,1]$ solve

$\int_{0}^{\underline{\mathfrak{M}}(s)}\frac{D(\varrho)}{\varrho^{a-1}}\,\mathrm{d}\varrho=\frac{v^{2}}{\gamma}F(s).$

Then there exists $\underline{s}\in(0,1)$ independent of $\varepsilon$, such
that

$\ell(s)<\underline{\mathfrak{M}}(s)<\mathfrak{M}^{\varepsilon}(s)\text{ for
all }\underline{s}\leq s<1.$

###### Proof.

Observe from (P1) that $\underline{\mathfrak{M}}$ is well-defined and
satisfies the ODE,

$\displaystyle\dfrac{\mathrm{d}\Phi(\underline{\mathfrak{M}})}{\mathrm{d}s}=-\frac{v^{2}}{\gamma}\dfrac{\underline{\mathfrak{M}}^{a}}{f(s)}$
(3.11)

and $\underline{\mathfrak{M}}(1)=0$. Subtracting (3.10b) and integrating in
$(s,1)$ one gets

$\displaystyle\Phi(\mathfrak{M}^{\varepsilon}(s))-\Phi(\underline{\mathfrak{M}}(s))=\Phi(\varepsilon)+\frac{v^{2}}{\gamma}\int_{s}^{1}\left[(\mathfrak{M}^{\varepsilon}(\varrho)-\underline{\mathfrak{M}}(\varrho))+(\underline{\mathfrak{M}}(\varrho)-\underline{\mathfrak{M}}^{a}(\varrho)-\ell(\varrho))\right]\frac{\mathrm{d}\varrho}{f(\varrho)}$
$\displaystyle=\Phi(\varepsilon)+\frac{v^{2}}{\gamma}\int_{s}^{1}\frac{1}{f(\varrho)}(\mathfrak{M}^{\varepsilon}-\underline{\mathfrak{M}})(\varrho)+\frac{v^{2}}{\gamma}\int_{s}^{1}\frac{\underline{\mathfrak{M}}(\varrho)}{f(\varrho)}\left(1-\underline{\mathfrak{M}}^{a-1}(\varrho)-\tfrac{\ell(\varrho)}{\underline{\mathfrak{M}}(\varrho)}\right).$
(3.12)

Note that,
$\ell^{\prime}(1)=-\gamma^{-1}(1-\lambda/f(1))\overset{\ref{prop:f}}{<}0$,
whereas $\underline{\mathfrak{M}}(s)\to 0$ as $s\to 1$, and therefore,

$\frac{d\underline{\mathfrak{M}}}{ds}=-\frac{v^{2}}{\gamma}\frac{\underline{\mathfrak{M}}^{a-1}}{f(s)\,D(\underline{\mathfrak{M}})}\to-\infty\text{
as }s\to 1.$

Hence, using L’Hôpital’s rule, there exists $\underline{s}\in(0,1)$
independent of $\varepsilon$, such that

$\displaystyle\left(1-\underline{\mathfrak{M}}^{a-1}(s)-\frac{\ell(s)}{\underline{\mathfrak{M}}(s)}\right)\geq
0\text{ for all }\underline{s}\leq s<1.$ (3.13)

Observe that $\mathfrak{M}^{\varepsilon}(s)>\underline{\mathfrak{M}}(s)$ in a
left neighbourhood of $s=1$, simply because
$\mathfrak{M}^{\varepsilon}(1)=\varepsilon>0=\underline{\mathfrak{M}}(1)$.
Then, (3.12)–(3.13) imply that
$\mathfrak{M}^{\varepsilon}(s)>\underline{\mathfrak{M}}(s)$ for all
$\underline{s}\leq s<1$. To see this, assume the contrary, i.e.,
$\mathfrak{M}^{\varepsilon}(s_{1})=\underline{\mathfrak{M}}(s_{1})$ for some
$s_{1}\in(\underline{s},1)$ and
$\mathfrak{M}^{\varepsilon}(s)>\underline{\mathfrak{M}}(s)$ for $s_{1}<s<1$.
Then, from (3.12)–(3.13) we have
$\Phi(\mathfrak{M}^{\varepsilon}(s_{1}))>\Phi(\underline{\mathfrak{M}}(s_{1}))$,
thus contradicting our assumption. This concludes the proof. ∎

######  Theorem 3.3 (Existence of an orbit connecting with $(0,1)$).

Let $v>0$ be fixed and (2.13) be satisfied. Let $\mathfrak{M}^{\varepsilon}$
denote the $\mathfrak{M}$-mapping introduced in Lemma 3.2 with
$\mathfrak{M}^{\varepsilon}(1)=\varepsilon$. Then there exists a function
$\mathfrak{M}:[s_{-\infty},1]\to[0,1)$ which satisfies (3.10) with
$\mathfrak{M}(1)=0$ and for all $s\in[s_{-\infty},1]$,
$\mathfrak{M}^{\varepsilon}(s)\to\mathfrak{M}(s)$ as $\varepsilon\to 0$ (see
Figure 3). Moreover, define the function $\zeta:(s_{-\infty},1]\to(-\infty,0]$
$($the $\zeta$-$\mathrm{map})$ as

$\displaystyle\zeta(s):=-\frac{v}{\gamma}\int_{s}^{1}\frac{\mathrm{d}\varrho}{f(\varrho)\,\mathfrak{M}(\varrho)}.$
(3.14)

Then $\zeta$ is differentiable and increasing with $s$ whenever
$\mathfrak{M}(s)>0$, and $\lim_{s\nearrow 1}\zeta(s)=0$. For any
$s\in(s_{-\infty},1]$ and $\zeta(s)\in{\mathbb{R}}^{-}$, defining
$\xi=\zeta(s)$, $M=\mathfrak{M}(s)$ and $S=s$, the mapping $\xi\mapsto(M,S)$
solves (2.11).

$m$$s$$0$$1$$1$$m=\ell(s)$$0<\varepsilon_{1}<\varepsilon_{2}$$(\varepsilon_{1},1)$$(\varepsilon_{2},1)$$\underline{\mathfrak{M}}$$\underline{s}$$\mathfrak{M}$$\mathfrak{M}^{\varepsilon_{1}}$$\mathfrak{M}^{\varepsilon_{2}}$
Figure 3: The $\mathfrak{M}$-mapping, introduced in Theorem 3.3 with
$\mathfrak{M}(1)=0$, and two $\mathfrak{M}^{\varepsilon}$-mappings, introduced
in Lemma 3.2 with $\mathfrak{M}^{\varepsilon}(1)=\varepsilon$ for
$\varepsilon=\varepsilon_{1}>0$ and
$\varepsilon=\varepsilon_{2}>\varepsilon_{1}$. The function
$\underline{\mathfrak{M}}$ defined in Lemma 3.2, providing a lower bound, is
also shown.

###### Proof.

For a fixed $v>0$ and $s\in(s_{-\infty},1]$,

$\displaystyle\mathfrak{M}_{\varepsilon_{1}}(s)\leq\mathfrak{M}_{\varepsilon_{2}}(s)<1$
$\displaystyle\text{ if }0<\varepsilon_{1}<\varepsilon_{2}<1,$ and the
equality holds only if
$\displaystyle\mathfrak{M}_{\varepsilon_{1}}(s)=\mathfrak{M}_{\varepsilon_{2}}(s)=0.$
(3.15)

This is evident since orbits, corresponding to initial values
$(\varepsilon_{1},1)$ and $(\varepsilon_{2},1)$, do not intersect in the
interior of $\mathfrak{R}$ due to uniqueness of solutions, see Lemma 3.1. The
lemma further yields that $\mathfrak{M}^{\varepsilon}<1$ if $\varepsilon<1$.
Since, $\mathfrak{M}^{\varepsilon}(s)$ is bounded below by
$\underline{\mathfrak{M}}(s)$ for $s\in(\underline{s},1)$, see Lemma 3.2,
there exists $\mathfrak{M}(s)$ such that

$\displaystyle\mathfrak{M}(s):=\lim\limits_{\varepsilon\searrow
0}\mathfrak{M}^{\varepsilon}(s)\;\;\text{ for }s\in(s_{-\infty},1),\;\;\text{
and }\;\;\mathfrak{M}(s)\geq\underline{\mathfrak{M}}(s)>0\text{ for
}s\in(\underline{s},1),$ (3.16)

see Figure 3. Let us take $s\in(\underline{s},1)$. Observe that by (3.10b),
$\mathfrak{M}^{\varepsilon}$ satisfies in this interval

$\displaystyle\Phi(\mathfrak{M}^{\varepsilon}(s))=\Phi(\varepsilon)+\frac{v^{2}}{\gamma}\int_{s}^{1}(\mathfrak{M}^{\varepsilon}(\varrho)-\ell(\varrho))\,\frac{\mathrm{d}\varrho}{f(\varrho)}.$
(3.17)

Due to (3.15), if $\varepsilon_{0}\in(0,1)$, then $\mathfrak{M}^{\varepsilon}$
is uniformly bounded away from 1 in $[s_{-\infty},1]$ for all
$\varepsilon\in(0,\varepsilon_{0}]$. Thus, $\Phi$ can be assumed to be locally
Lipschitz in the above equation. Hence, passing the limit $\varepsilon\to 0$,
we get for all $s\in(\underline{s},1)$,

$\displaystyle\Phi(\mathfrak{M}(s))=\frac{v^{2}}{\gamma}\int_{s}^{1}(\mathfrak{M}(\varrho)-\ell(\varrho))\,\frac{\mathrm{d}\varrho}{f(\varrho)},$
(3.18)

which upon differentiation proves that $\mathfrak{M}$ satisfies (3.10) in
$(\underline{s},1)$ with $(M_{0},S_{0})=(0,1)$. For $s<\underline{s}$, the
mapping $\mathfrak{M}$ is simply extended by solving the equation (3.10) with
$(M_{0},S_{0})=(\mathfrak{M}(\underline{s}),\underline{s})\in\mathfrak{R}$.
The existence of $\mathfrak{M}$ in this case also follows from the existence
of orbits, i.e., Lemma 3.1.

The differentiability and monotonicity of the function $\zeta$ is obvious from
(3.14). To prove that $\lim_{s\nearrow 1}\zeta(s)=0$, we estimate for any
$s\in(\underline{s},1)$,

$\displaystyle
0<\int_{s}^{1}\frac{\mathrm{d}\varrho}{f(\varrho)\,\mathfrak{M}(\varrho)}\overset{\eqref{eq:MepsOrder}}{\leq}\int_{s}^{1}\frac{\mathrm{d}\varrho}{f(\varrho)\,\underline{\mathfrak{M}}(\varrho)}\overset{\eqref{eq:McLB}}{=}\frac{\gamma}{v^{2}}\int^{\underline{\mathfrak{M}}(s)}_{0}\frac{D(m)}{m^{a}}\mathrm{d}m<\infty.$
(3.19)

Hence, $\zeta(s)>-\infty$ for all $s\in(\underline{s},1]$, and passing to the
limit $s\to 1$ one obtains that $\zeta(1^{-})=0$.

Differentiating (3.14) and using (3.10a), it immediately follows that the
mapping $\xi\mapsto(M,S)$ solves (2.11). ∎

######  Corollary 3.3.1 (Behaviour of the orbit connecting to $(0,1)$).

Let $(M,S)$ be the orbit defined in Theorem 3.3. Then there exists
$\underline{\xi}\in{\mathbb{R}}^{-}\cup\\{-\infty\\}$ such that
$(M,S)\in\mathfrak{R}=[0,1)\times[s_{-\infty},1]$ for all
$\xi\geq\underline{\xi}$. In $\mathfrak{R}\cap\\{s\geq s_{*}\\}$, both $M$ and
$S$ increase with $\xi$ until $M=\ell(S)$ is satisfied for some $S>s^{*}$
$(s_{*},\;s^{*}\in(s_{-\infty},1)$ defined in (3.6)$)$, after which $M$
decreases and $S$ remains increasing. Finally, $(M,S)=(0,1)$ for all $\xi\geq
0$.

The statement is evident from the direction of orbits in the phase–plane,
Figure 2, Theorem 3.3 and Lemma 3.1. Below, we prove that $(M,S)$ is the
unique orbit which connects with $(0,1)$ at a finite $\xi$-coordinate.

###  3.3 Uniqueness of the orbit connecting with $(0,1)$ at $\xi=0$

For a given $v>0$, there are in fact infinitely many orbits $(M,S)$ that
connect to $(0,1)$ as $\xi\to+\infty$. However, only one orbit, the one
constructed in Section 3.2, connects to $(0,1)$ at $\xi=0$. This statement
will be proved below. This is a common phenomenon for TW solutions of
degenerate diffusion equations [4, 5] and the unique orbit corresponds to the
TW with minimum speed in these cases.

######  Proposition 3.4 (Uniqueness of the orbit connecting to $(0,1)$ for
some $\xi\in{\mathbb{R}}$).

For a fixed $v>0$, let
$(\widetilde{M},\widetilde{S})\in(C^{1}({\mathbb{R}}))^{2}$ be an orbit
satisfying (2.11) and connecting with $(0,1)$ from $\\{s<1\\}$. Let
$\widetilde{\mathfrak{M}}$ and $\widetilde{\zeta}$ denote the corresponding
$\mathfrak{M}$-mapping (Definition 1) and $\zeta$-mapping (Theorem 3.3) of
$(\widetilde{M},\widetilde{S})$ with $\widetilde{\xi}=\widetilde{\zeta}(1)$.
Then, $\widetilde{\xi}<\infty$ if and only if $(\widetilde{M},\widetilde{S})$
is the orbit defined in Theorem 3.3.

###### Proof.

We first show that

$\displaystyle\widetilde{\xi}=\widetilde{\zeta}(1)<\infty\quad\text{ implies
}\quad\frac{\mathrm{d}}{\mathrm{d}s}\widetilde{\mathfrak{M}}(1)=-\infty.$
(3.20)

Assume the contrary, i.e,
$\frac{\mathrm{d}}{\mathrm{d}s}\widetilde{\mathfrak{M}}(1)>-\infty$. Then,
there exists $s\in(s_{-\infty},1)$ and $L>0$ such that

$\widetilde{\mathfrak{M}}(\varrho)\leq L\,(1-\varrho)\text{ for all
}\varrho\in[s,1].$

This implies $\widetilde{\xi}=\widetilde{\zeta}(1)=\infty$ since, writing
formally,

$\widetilde{\zeta}(1)-\widetilde{\zeta}(s)=\frac{v}{\gamma}\int^{1}_{s}\frac{\mathrm{d}\varrho}{f(\varrho)\,\widetilde{\mathfrak{M}}(\varrho)}\geq\frac{v}{\gamma}\int^{1}_{s}\frac{\mathrm{d}\varrho}{L\,f(1)\,(1-\varrho)}=\infty.$

Now, let us assume
$\frac{\mathrm{d}}{\mathrm{d}s}\widetilde{\mathfrak{M}}(1)=-\infty$. Since
$\widetilde{\mathfrak{M}}$ is continuous and differentiable in a left
neighbourhood of $s=1$, and $\ell^{\prime}(1)<0$ is bounded, for any given
$\nu\in(0,1)$, by the L’Hôpital’s rule, there exists
$\widetilde{s}_{\nu}\in(0,1)$ such that

$\displaystyle\frac{\ell(s)}{\widetilde{\mathfrak{M}}(s)}<1-\nu\text{ for all
}s\in[\widetilde{s}_{\nu},1].$ (3.21)

Let $(M,S)$ be the orbit defined in Theorem 3.3 with $\mathfrak{M}$ and
$\zeta$ being the corresponding mappings. To shorten notations we introduce

$\displaystyle\phi=\Phi(\mathfrak{M}),\quad\widetilde{\phi}=\Phi(\widetilde{\mathfrak{M}}),\;\text{
and }\;\eth\phi=\phi-\widetilde{\phi}.$ (3.22)

Since $\mathfrak{M}(s)=\lim_{\varepsilon\searrow
0}\mathfrak{M}^{\varepsilon}(s)$, see Theorem 3.3, the ordering of orbits for
a fixed $v>0$ implies that

$\displaystyle\mathfrak{M}(s)\geq\widetilde{\mathfrak{M}}(s),\quad\phi\geq\widetilde{\phi},\quad\eth\phi\geq
0,\;\text{ for all }s\in(s_{-\infty},1).$ (3.23)

From (3.10b) applied to $\mathfrak{M}$ and $\widetilde{\mathfrak{M}}$, we then
obtain

$\displaystyle\frac{\mathrm{d}(\eth\phi)}{\mathrm{d}s}=-\frac{v^{2}}{\gamma
f(s)}[\Phi^{-1}(\phi)-\Phi^{-1}(\widetilde{\phi})].$ (3.24)

Since $\Phi$ is convex and strictly increasing, $\Phi^{-1}$ exists and is
concave, and
$\\{\Phi^{-1}\\}^{\prime}(\Phi(m))=1/\Phi^{\prime}(m)=1/(m\,D(m))$. Hence,
(3.23) implies
$\Phi^{-1}(\phi)-\Phi^{-1}(\widetilde{\phi})\leq\\{\Phi^{-1}\\}^{\prime}(\widetilde{\phi})\,\eth\phi={\eth\phi}/{(\widetilde{\mathfrak{M}}\,D(\widetilde{\mathfrak{M}}))}$.
Thus, (3.24) yields

$\displaystyle\frac{\mathrm{d}(\eth\phi)}{\mathrm{d}s}+\frac{v^{2}}{\gamma
f(s)}\frac{\eth\phi}{\widetilde{\mathfrak{M}}\,D(\widetilde{\mathfrak{M}})}\geq
0.$ (3.25)

Using the integrating factor
$\exp\left(\frac{v^{2}}{\gamma}\int_{\widetilde{s}_{\nu}}^{s}\frac{\mathrm{d}\varrho}{f\,\widetilde{\mathfrak{M}}\,D(\widetilde{\mathfrak{M}})}\right)$
and integrating in $(\widetilde{s}_{\nu},s)$ one has

$\displaystyle\eth\phi(s)\geq\eth\phi(\widetilde{s}_{\nu})\,\exp\left(-\frac{v^{2}}{\gamma}\int_{\widetilde{s}_{\nu}}^{s}\frac{\mathrm{d}\varrho}{f(\varrho)\,\widetilde{\mathfrak{M}}(\varrho)\,D(\widetilde{\mathfrak{M}}(\varrho))}\right).$
(3.26)

Observe from (3.21) and (3.10a) that
$\frac{\mathrm{d}\widetilde{\mathfrak{M}}}{\mathrm{d}s}<0$ for
$s\in(\widetilde{s}_{\nu},1)$, which gives

$\displaystyle-\frac{v^{2}}{\gamma}\frac{1}{f\,\widetilde{\mathfrak{M}}\,D(\widetilde{\mathfrak{M}})}\overset{\eqref{eq:McEqA}}{=}\frac{1}{\widetilde{\mathfrak{M}}(s)-\ell(s)}\frac{\mathrm{d}\widetilde{\mathfrak{M}}}{\mathrm{d}s}=\tfrac{1}{\left(1-\frac{\ell(s)}{\widetilde{\mathfrak{M}}(s)}\right)}\frac{1}{\widetilde{\mathfrak{M}}}\frac{\mathrm{d}\widetilde{\mathfrak{M}}}{\mathrm{d}s}\overset{\eqref{eq:elltildeM}}{\geq}\frac{1}{\nu}\frac{1}{\widetilde{\mathfrak{M}}}\frac{\mathrm{d}\widetilde{\mathfrak{M}}}{\mathrm{d}s}.$
(3.27)

Putting this in (3.26), we have

$\eth\phi(s)\geq\eth\phi(\widetilde{s}_{\nu})\,\exp\left(\frac{1}{\nu}\int_{\widetilde{s}_{\nu}}^{s}\frac{\mathrm{d}\widetilde{\mathfrak{M}}}{\widetilde{\mathfrak{M}}}\right)=\eth\phi(\widetilde{s}_{\nu})\left(\frac{\widetilde{\mathfrak{M}}(s)}{\widetilde{\mathfrak{M}}(\widetilde{s}_{\nu})}\right)^{\frac{1}{\nu}}.$

Rearranging the above relation using (3.22) one obtains for a constant
$\widetilde{C}^{1}_{\nu}>0$ only dependent on $\widetilde{s}_{\nu}<1$ that

$\displaystyle\Phi(\widetilde{\mathfrak{M}})+(\widetilde{C}^{1}_{\nu}\widetilde{\mathfrak{M}})^{\frac{1}{\nu}}\leq\Phi(\mathfrak{M}),\quad\text{
or
}\quad\widetilde{C}^{1}_{\nu}\widetilde{\mathfrak{M}}\leq(\Phi(\mathfrak{M}))^{\nu}.$
(3.28)

Recalling (P1) and the definition of $\Phi$ in (3.9), we choose $\nu\in(0,1)$
such that there exists a constant $\widetilde{C}^{2}_{\nu}>0$ for which

$\displaystyle\Phi(m)^{\nu}\leq\widetilde{C}^{2}_{\nu}\Phi^{\prime}(m)\overset{\eqref{eq:DefPhi}}{=}\widetilde{C}^{2}_{\nu}\,m\,D(m)\text{
for all }m\in(0,1).$ (3.29)

In a right neighbourhood $[0,\varepsilon]$ of $m=0$ ($\epsilon>0$), where
$D(m)\sim Cm^{a}$ (see (P1)) for some constant $C>0$, we have
$\Phi^{\prime}(m)\sim Cm^{a+1}$ and $\Phi(m)\sim\frac{C}{a+2}m^{a+2}$. Hence,
for any constant $\nu\in[(1+a)/(2+a),1)$ there exists
$\widetilde{C}^{2}_{\nu}>0$ satisfying (3.29) in $[0,\varepsilon]$. For
$m\in(0,1)$, one can take $\widetilde{C}^{2}_{\nu}$ to be the minimum of the
value of $\widetilde{C}^{2}_{\nu}$ for $m\in[0,\varepsilon]$, and
$\min_{m\in(\varepsilon,1)}\Phi(m)^{\nu}(\Phi^{\prime}(m))^{-1}>0$. Then, for
$s\in(\widetilde{s}_{\nu},1)$ one has

$\displaystyle\widetilde{\zeta}(s)-\widetilde{\zeta}(\widetilde{s}_{\nu})$
$\displaystyle=\frac{v}{\gamma}\int_{\widetilde{s}_{\nu}}^{s}\frac{\mathrm{d}\varrho}{f(\varrho)\,\widetilde{\mathfrak{M}}(\varrho)}\overset{\eqref{eq:TildePhiIneq}}{\geq}\frac{v}{\gamma\widetilde{C}^{1}_{\nu}}\int_{\widetilde{s}_{\nu}}^{s}\frac{\mathrm{d}\varrho}{f(\varrho)\,\Phi(\mathfrak{M}(\varrho))^{\nu}}$
$\displaystyle\overset{\eqref{eq:PhiGrowthRate}}{\geq}\frac{v}{\gamma\widetilde{C}^{1}_{\nu}\widetilde{C}^{2}_{\nu}}\int_{\widetilde{s}_{\nu}}^{s}\frac{\mathrm{d}\varrho}{f(\varrho)\,\mathfrak{M}(\varrho)D(\mathfrak{M}(\varrho))}\overset{\eqref{eq:GEs}}{=}\frac{1}{\widetilde{C}^{1}_{\nu}\widetilde{C}^{2}_{\nu}}\int_{\widetilde{s}_{\nu}}^{s}\left(\frac{\mathrm{d}\tau}{\mathrm{d}S}\right)\mathrm{d}S.$
(3.30)

The orbit $(M,S)\to(0,1)$ only as $\tau\to\infty$ since $(0,1)$ is an
equilibrium point of (3.5) and the right hand side of (3.5) is locally
Lipschitz with respect to $M$ and $S$. Hence, the right hand side of (3.3)
tends to $+\infty$ as $s\nearrow 1$. This contradicts
$\widetilde{\xi}=\widetilde{\zeta}(1)$ being bounded, which concludes the
proof of the proposition.

∎

By this point we have found a unique orbit which satisfies half of the
boundary conditions in (2.12), i.e. $(M,S)(0)=(0,1)$. For the boundary
condition at the other end, i.e. $\xi\to-\infty$, we investigate how the orbit
varies with $v$.

###  3.4 Ordering of the orbits with respect to $v$

$m$$s$$0$$1$$1$$s^{*}$$s_{*}$$\;m=\ell(s)$$v_{1}<v_{2}<v_{3}<v_{4}$$\mathfrak{M}_{v_{1}}$$\mathfrak{M}_{v_{2}}$$\mathfrak{M}_{v_{3}}$$\mathfrak{M}_{v_{4}}$
Figure 4: The ordering of orbits as stated in Theorem 3.5 computed for the
parameter set provided in [11].

######  Theorem 3.5 (The ordering of the orbits with respect to $v$).

For $v>0$, let the function $\mathfrak{M}_{v}:[s_{-\infty},1)\to[0,1]$ denote
the $\mathfrak{M}$-mapping introduced in Theorem 3.3 that satisfies
$\mathfrak{M}_{v}(1)=0$. Then, for a fixed $s\in[s_{-\infty},1)$,
$\mathfrak{M}_{v}(s)$ varies continuously with $v$. Moreover,

  1. (i)

$\mathfrak{M}_{v}(s)$ is strictly increasing with $v$.

  2. (ii)

$\mathfrak{M}_{v}(s)\to 0$ as $v\to 0$.

  3. (iii)

Let $\mathcal{G}(s)>0$ for all $s\in(s_{-\infty},1)$. Then,
$\mathfrak{M}_{v}(s)\to 1$ as $v\to\infty$.

Schematically the ordering of the orbits is shown in Figure 4.

###### Proof.

_(i)_ For a fixed $v>0$, let $\mathfrak{M}^{\varepsilon}_{v}$ denote the
$\mathfrak{M}$-mapping discussed in Lemma 3.2 with
$\mathfrak{M}^{\varepsilon}_{v}(1)=\varepsilon>0$. Let $0<v_{1}<v_{2}$.
Observe that

$\frac{\mathrm{d}}{\mathrm{d}s}\Phi(\mathfrak{M}^{\varepsilon}_{v_{1}}(1))=-\frac{v^{2}_{1}}{\gamma\,f(1)}\varepsilon>-\frac{v^{2}_{2}}{\gamma\,f(1)}\,\varepsilon=\frac{\mathrm{d}}{\mathrm{d}s}\Phi(\mathfrak{M}^{\varepsilon}_{v_{2}}(1)),$

implying
$\Phi(\mathfrak{M}^{\varepsilon}_{v_{1}}(s))<\Phi(\mathfrak{M}^{\varepsilon}_{v_{2}}(s))$
in a left neighbourhood of $s=1$. We show below that
$\Phi(\mathfrak{M}_{v_{1}}^{\varepsilon}(s))<\Phi(\mathfrak{M}_{v_{2}}^{\varepsilon}(s))$
for all $s\in[s_{-\infty},1)$ provided
$\mathfrak{M}_{v_{2}}^{\varepsilon}(s)>0$. To show this, let us assume the
contrary, i.e.
$\mathfrak{M}_{v_{1}}^{\varepsilon}(s_{1})=\mathfrak{M}_{v_{2}}^{\varepsilon}(s_{1})>0$
for some $s_{1}\in[s_{-\infty},1)$ with
$\mathfrak{M}^{\varepsilon}_{v_{1}}(s)<\mathfrak{M}^{\varepsilon}_{v_{2}}(s)$
for all $s\in(s_{1},1)$. Subtracting two versions of (3.10b), we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}(\Phi(\mathfrak{M}_{v_{1}}^{\varepsilon}(s))-\Phi(\mathfrak{M}_{v_{2}}^{\varepsilon}(s)))$
$\displaystyle=\frac{1}{f(s)}\left[-\frac{v^{2}_{1}}{\gamma}(\mathfrak{M}_{v_{1}}^{\varepsilon}-\mathfrak{M}_{v_{2}}^{\varepsilon})-\dfrac{v_{2}^{2}-v_{1}^{2}}{\gamma}[\ell(s)-\mathfrak{M}_{v_{2}}^{\varepsilon}]\right]$
$\displaystyle=$
$\displaystyle-\frac{v^{2}_{1}}{\gamma\,f(s)}(\mathfrak{M}_{v_{1}}^{\varepsilon}-\mathfrak{M}_{v_{2}}^{\varepsilon})-\left(1-\frac{v_{1}^{2}}{v_{2}^{2}}\right)\frac{\mathrm{d}}{\mathrm{d}s}\Phi(\mathfrak{M}_{v_{2}}^{\varepsilon}(s)).$
(3.31)

Integrating in $(s_{1},1)$ we get

$\displaystyle\Phi(\mathfrak{M}_{v_{2}}^{\varepsilon}(s_{1}))-\Phi(\mathfrak{M}_{v_{1}}^{\varepsilon}(s_{1}))=\frac{v^{2}_{1}}{\gamma}\int_{s_{1}}^{1}(\mathfrak{M}_{v_{2}}^{\varepsilon}-\mathfrak{M}_{v_{1}}^{\varepsilon})\frac{\mathrm{d}s}{f(s)}+\left(1-\frac{v_{1}^{2}}{v_{2}^{2}}\right)\Phi(\mathfrak{M}_{v_{2}}^{\varepsilon}(s_{1}))>0,$

thus contradicting our assumption. Hence,
$\Phi(\mathfrak{M}_{v_{1}}^{\varepsilon}(s))<\Phi(\mathfrak{M}_{v_{2}}^{\varepsilon}(s))$
for all $s\in[s_{-\infty},1)$ and passing to the limit $\varepsilon\to 0$ we
obtain

$\displaystyle\Phi(\mathfrak{M}_{v_{1}}(s))\leq\Phi(\mathfrak{M}_{v_{2}}(s))\text{
for all }s\in[s_{-\infty},1).$ (3.32)

To prove that $\mathfrak{M}_{v}$ strictly increases with $v$, we observe
integrating (3.10b) in $(s,1)$ that

$\frac{\Phi(\mathfrak{M}_{v_{2}}(s))}{v_{2}^{2}}=\frac{1}{\gamma}\int_{s}^{1}[\mathfrak{M}_{v_{2}}-\ell](\varrho)\frac{\mathrm{d}\varrho}{f(\varrho)}\overset{\eqref{eq:McMonotone}}{\geq}\frac{1}{\gamma}\int_{s}^{1}[\mathfrak{M}_{v_{1}}-\ell](\varrho)\frac{\mathrm{d}\varrho}{f(\varrho)}=\frac{\Phi(\mathfrak{M}_{v_{1}}(s))}{v_{1}^{2}},$

which yields the important inequality

$\displaystyle\Phi(\mathfrak{M}_{v_{2}}(s))\geq\frac{v_{2}^{2}}{v_{1}^{2}}\Phi(\mathfrak{M}_{v_{1}}(s))>\Phi(\mathfrak{M}_{v_{1}}(s)).$
(3.33)

This proves that $\mathfrak{M}_{v}(s)$ is strictly increasing with respect to
$v$.

_(ii)_ Integrating (3.10b) in $(s,1)$ we have

$\Phi(\mathfrak{M}_{v}(s))=\frac{v^{2}}{\gamma}\int_{s}^{1}(\mathfrak{M}_{v}(\varrho)-\ell(\varrho))\frac{\mathrm{d}\varrho}{f(\varrho)}<\frac{v^{2}}{\gamma}\int_{s}^{1}(1-\ell(\varrho))\frac{\mathrm{d}\varrho}{f(\varrho)}\overset{\eqref{eq:relationGell}}{=}\left(\frac{v}{\gamma}\right)^{2}\mathcal{G}(s).$

This proves point (ii).

_(iii)_ We prove the statement first in $(s^{*},1)$. Recall that
$s^{*}\in(0,1)$ satisfies $\ell(s^{*})=1$. Let $S^{\dagger}_{v}\in(s^{*},1)$
be such that $\mathfrak{M}_{v}(S^{\dagger}_{v})=\ell(S^{\dagger}_{v})$. The
existence of $S^{\dagger}_{v}$ is guaranteed by Corollary 3.3.1. Also, point
(i) implies that $S^{\dagger}_{v}$ decreases as $v$ increases. Let
$S^{\dagger}=\lim\limits_{v\to\infty}S^{\dagger}_{v}$. Then Corollary 3.3.1
yields again

$\displaystyle\max_{s\in(s_{*},1)}\mathfrak{M}_{v}(s)=\ell(S^{\dagger}_{v}),\text{
and hence,
}\sup_{v>0}\max_{s\in(s_{*},1)}\mathfrak{M}_{v}(s)=\ell(S^{\dagger}).$ (3.34)

Let $S^{\dagger}>s^{*}$, yielding $\ell(S^{\dagger})<1$. Then using (3.33) one
has for some $v_{1}=v>0$ and $v_{2}=v^{2}$ that

$\Phi(\mathfrak{M}_{v_{2}}(S^{\dagger}_{v_{1}}))\geq\frac{v_{2}^{2}}{v_{1}^{2}}\Phi(\mathfrak{M}_{v_{1}}(S^{\dagger}_{v_{1}}))=v^{2}\Phi(\ell(S^{\dagger}_{v}))\to\infty,\text{
as }v\to\infty,$

which contradicts the upper bound in (3.34) since $\ell(S^{\dagger})<1$.
Hence, $\lim\limits_{v\to\infty}S^{\dagger}_{v}=s^{*}$.

Now let $s\in(s^{*},1)$. Then, there exists $v_{1}>0$ such that
$S^{\dagger}_{v_{1}}\geq s$. Using (3.33) for $v>v_{1}$,

$\displaystyle\Phi(\mathfrak{M}_{v}(s))\geq\frac{v^{2}}{v_{1}^{2}}\Phi(\mathfrak{M}_{v_{1}}(s))\geq\frac{v^{2}}{v_{1}^{2}}\Phi(\ell(s))\to\infty\text{
as }v\to\infty.$ (3.35)

Hence, $\mathfrak{M}_{v}(s)\to 1$ for all $s\in(s^{*},1)$.

Now, we extend the result to $(s_{-\infty},s^{*})$. Set $\check{s}_{0}=s^{*}$
and let $\check{s}_{k}\in[s_{-\infty},\check{s}_{k-1}]$, $k\in{\mathbb{N}}$,
be recursively defined by the formula

$\int_{\check{s}_{k}}^{\check{s}_{k-1}}\frac{\ell(\varrho)}{f(\varrho)}\,\mathrm{d}\varrho=\int^{1}_{\check{s}_{k-1}}\frac{1-\ell(\varrho)}{f(\varrho)}\,\mathrm{d}\varrho\overset{\eqref{eq:DefG}}{=}\gamma^{-1}\mathcal{G}(\check{s}_{k-1})>0.$

If such $\check{s}_{k}$ does not exist in $[s_{-\infty},\check{s}_{k-1}]$ then
we set $\check{s}_{k}=s_{-\infty}$. Assume that
$\lim_{v\to\infty}\mathfrak{M}_{v}(s)\to 1$ holds for all
$s\in(\check{s}_{k-1},1)$. This is certainly true for $k=1$. Then we show that
it holds for all $s\in(\check{s}_{k},1)$. For all
$s\in(\check{s}_{k},\check{s}_{k-1})$ one has

$\displaystyle\Phi(\mathfrak{M}_{v}(s))$
$\displaystyle=\frac{v^{2}}{\gamma}\int_{s}^{1}\frac{1}{f}(\mathfrak{M}_{v}-\ell)=\frac{v^{2}}{\gamma}\left[\int_{\check{s}_{k-1}}^{1}\frac{1}{f}(\mathfrak{M}_{v}-\ell)+\int^{\check{s}_{k-1}}_{s}\frac{1}{f}(\mathfrak{M}_{v}-\ell)\right]$
$\displaystyle\geq\frac{v^{2}}{\gamma}\left[\int_{\check{s}_{k-1}}^{1}\frac{1}{f}(\mathfrak{M}_{v}-\ell)-\int^{\check{s}_{k-1}}_{s}\frac{\ell}{f}\right].$
(3.36)

The integral
$\int_{\check{s}_{k-1}}^{1}\frac{1}{f}(\mathfrak{M}_{v}-\ell)\to\int_{\check{s}_{k-1}}^{1}\frac{1}{f}(1-\ell)\overset{\eqref{eq:relationGell}}{=}\gamma^{-1}\mathcal{G}(\check{s}_{k-1})$
as $v\to\infty$. Thus
$\int_{\check{s}_{k-1}}^{1}\frac{1}{f}(\mathfrak{M}_{v}-\ell)-\int^{\check{s}_{k-1}}_{s}\frac{\ell}{f}>0$
for large $v$ and $s>\check{s}_{k}$. Hence, passing $v\to\infty$ in (3.4) one
gets, $\mathfrak{M}_{v}(s)\to 1$ for all $s>\check{s}_{k}$, thus proving the
statement.

_Continuity:_ Finally, we prove that $\lim\limits_{v\to
v_{0}}\mathfrak{M}_{v}(s)=\mathfrak{M}_{v_{0}}(s)$ for any $v_{0}>0$. Let us
handle the case $v\searrow v_{0}$ first. From point (i),
$\lim\limits_{v\searrow v_{0}}\mathfrak{M}_{v}(s)\geq\mathfrak{M}_{v_{0}}(s)$.
Let us assume

$\lim\limits_{v\searrow
v_{0}}\mathfrak{M}_{v}(s)=:\bar{\mathfrak{M}}_{v_{0}}(s)>\mathfrak{M}_{v_{0}}(s)\text{
for some }s\in(s_{-\infty},1).$

Just as in the proof of Theorem 3.3, $\bar{\mathfrak{M}}_{v_{0}}$ satisfies
(3.10) with $v=v_{0}$. Since
$\mathfrak{M}^{\varepsilon}_{v_{0}}(s)\searrow\mathfrak{M}_{v_{0}}(s)$ as
$\varepsilon\searrow 0$ ($\mathfrak{M}^{\varepsilon}_{v_{0}}$ as defined in
Lemma 3.2 for $v=v_{0}$), one can choose $\varepsilon>0$ small enough such
that
$\mathfrak{M}_{v_{0}}(s)<\mathfrak{M}_{v_{0}}^{\varepsilon}(s)<\bar{\mathfrak{M}}_{v_{0}}(s)$.
Then $\mathfrak{M}^{\varepsilon}_{v_{0}}$ and $\bar{\mathfrak{M}}_{v_{0}}$
both satisfy (3.10) and

$(\mathfrak{M}^{\varepsilon}_{v_{0}}-\bar{\mathfrak{M}}_{v_{0}})(1)=\varepsilon>0,\text{
whereas
}(\mathfrak{M}^{\varepsilon}_{v_{0}}-\bar{\mathfrak{M}}_{v_{0}})(s)<0,$

implying that they intersect at some intermediate point. Since both
$\mathfrak{M}^{\varepsilon}_{v_{0}}$ and $\bar{\mathfrak{M}}_{v_{0}}$
correspond to orbits in the phase–plane with the same $v=v_{0}$, and the
orbits cannot intersect in $\\{m>0\\}$ we have our contradiction.

The proof of the case $v\nearrow v_{0}$ follows from Proposition 3.4 since
$\bar{\mathfrak{M}}_{v_{0}}=\lim_{v\nearrow
v_{0}}\mathfrak{M}_{v}<\mathfrak{M}_{v_{0}}$ would imply that there are two
mappings, $\mathfrak{M}_{v_{0}}$ and $\bar{\mathfrak{M}}_{v_{0}}$, which
satisfy $\mathfrak{M}(1)=0$ with a corresponding $\zeta(1)=0$, thus
contradicting Proposition 3.4. ∎

From the proof of Theorem 3.5 (ii), we obtain the following property:

######  Corollary 3.5.1 (Bounds on $\mathfrak{M}_{v}$).

For $v>0$, let the function $\mathfrak{M}_{v}:[s_{-\infty},1]\to[0,1]$ denote
the $\mathfrak{M}$-mapping introduced in Theorem 3.3 that satisfies
$\mathfrak{M}_{v}(1)=0$. Then, for all $v_{0}<v$ one has

$\frac{v^{2}}{v_{0}^{2}}\Phi(\mathfrak{M}_{v_{0}}(s))\leq\Phi(\mathfrak{M}_{v}(s))\leq\left(\frac{v}{\gamma}\right)^{2}\mathcal{G}(s).$

###  3.5 Existence/non-existence of travelling waves

###### Proof of Theorem 2.1.

Theorem 3.3 and Proposition 3.4 prove the existence of a unique orbit
$(M_{v},S_{v})$ satisfying (2.11) and connecting with $(0,1)$ at $\xi=0$.
Corollary 3.3.1 shows that the orbit crosses the line $m=\ell(s)$ for some
$s\in(s^{*},1)$ and Lemma 3.1 shows that it enters
$\mathfrak{R}=[0,1)\times(s_{-\infty},1]$ through either
$\\{m=0,s\in[s_{-\infty},1]\\}$ or through $\\{s=s_{-\infty}\\}$. Theorem 3.5
proves that the orbit varies continuously with $v$, and for large $v$, the
orbit enters through $s=s_{-\infty}$. Hence, it remains to be shown that for
small $v$, the orbit enters through $\\{m=0,s\in[s_{-\infty},1]\\}$. This will
prove, by continuity, the existence of $v=\bar{v}>0$ such that the
corresponding orbit connects with $(0,s_{-\infty})$ which is the intersection
point of the two segments.

We show that for any $\hat{s}\in[s_{-\infty},1)$ there exists a corresponding
$\hat{v}>0$ such that

$\mathfrak{M}_{\hat{v}}(\hat{s})=0\text{ for all }v\leq\hat{v}.$

Assume that no such $\hat{v}>0$ exists. Then $\mathfrak{M}_{v}(\hat{s})>0$ for
all $v>0$. Integrating (3.10b) in $(\hat{s},1)$ one then has

$0<\Phi(\mathfrak{M}_{v}(\hat{s}))=\frac{v^{2}}{\gamma}\int_{\hat{s}}^{1}(\mathfrak{M}_{v}(\varrho)-\ell(\varrho))\,\frac{\mathrm{d}\varrho}{f(\varrho)}.$

However, Corollary 3.5.1 implies that
$\int_{\hat{s}}^{1}(\mathfrak{M}_{v}(\varrho)-\ell(\varrho))\,\frac{\mathrm{d}\varrho}{f(\varrho)}<0$
for small enough $v>0$, since $\mathfrak{M}_{v}(\varrho)\searrow 0$ uniformly
as $v\searrow 0$. This is a contradiction to
$\Phi(\mathfrak{M}_{v}(\hat{s}))>0$. Hence, the hypothesized $\hat{v}$ exists.
Setting $\hat{s}=s_{-\infty}$ proves Theorem 2.1. ∎

The profiles of $(M,S)$ as functions of $\xi\leq 0$ are shown in Figure 5. An
interesting feature of this TW is that the profile of $M$ has a sharp front at
$\xi=0$, whereas, it has a diffused tail at the rear. The TW is in fact a
travelling pulse since it connects an equilibrium state with $M=0$ to another
equilibrium state with $M=0$, in contrast to TWs for nonlinear diffusion
problems with Fischer type source terms, see [12, 5, 4].

Figure 5: The TW profile for the default parameter set given in Table 1
obtained using Algorithm 5.1. In the (left) plot the graphs of $M(\xi)$,
$S(\xi)$ and $\omega(\xi)=\int_{-\infty}^{\xi}M$ are shown. The (right) plot
shows the orbit $(M,S)$ in the phase–plane (note that $s_{-\infty}\approx
10^{-60}$ in this case). The corresponding wave-speed is $v=0.01454$.

######  Remark 3.2 (Consistency of the TW solution).

Let $(M,S)$ be the TW solution described in Theorem 2.1 and let $\mathfrak{M}$
denote its $\mathfrak{M}$-mapping. Then $\ell(s)/\mathfrak{M}\to 0$ as
$s\nearrow 1$ from (3.20), and consequently recalling the TW Assumption 2.1
one has

$\displaystyle\begin{cases}-D(M)\partial_{x}M=-D(M)\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=v(M-\ell(S))\to
0,\\\
-\tfrac{D(M)}{M}\partial_{x}M=-\tfrac{D(M)}{M}\tfrac{\mathrm{d}}{\mathrm{d}\xi}M=v\left(1-\tfrac{\ell(S)}{M}\right)\to
v,\end{cases}\quad\text{ as }x-vt=\xi\nearrow 0.$ (3.37)

Here, the term $-D(M)\partial_{x}M$ represents the flux and the term
$-\tfrac{D(M)}{M}\partial_{x}M$ is commonly referred to as the particle speed
in the literature of the porous medium equation, see [30, Chapter 19]. The
fact that the flux vanishes at $\xi=0$ and the particle velocity is $v$
implies that the TW is physically consistent with the model (1.1).

######  Remark 3.3 (Regularity of the sharp front).

From the limits in (3.37), one can also deduce the regularity of $M(\xi)$ as
$\xi\nearrow 0$. Taking the expression (1.3a) for $D$, one has in the limit
$\xi\nearrow 0$, observing that $M(\xi)\searrow 0$, that
$\delta\,M^{a-1}\,\tfrac{\mathrm{d}}{\mathrm{d}\xi}M\sim-v$, or

$\displaystyle M(\xi)\sim\tfrac{a\,v}{\delta}|\xi|^{\frac{1}{a}},\text{ for
}\xi\text{ in a left neighbourhood of }0.$ (3.38)

Consequently, a sharp front exists at $\xi=0$, however, the TW profile still
remains Hölder continuous.

Finally, the techniques developed in this section also allow us to give
necessary conditions for existence of the TWs.

######  Proposition 3.6 (Non-existence of travelling wave solutions).

Let (2.13) be satisfied and $\mathcal{G}(\bar{s})<0$ for some
$\bar{s}\in(s_{-\infty},1)$. Then there exists no TW orbit $(M,S)$ satisfying
both (2.11) and (2.12). Similarly, if $\lambda\geq 1$ then no TW exists.

###### Proof.

For $\mathcal{G}(\bar{s})<0$, assume that such a solution exists for some
$v>0$. Then, integrating (3.10) in $(\bar{s},1)$ one has

$0<\frac{\gamma}{v^{2}}\Phi(\mathfrak{M}_{v}(\bar{s}))=\int_{\bar{s}}^{1}(\mathfrak{M}_{v}(\varrho)-\ell(\varrho))\,\frac{\mathrm{d}\varrho}{f(\varrho)}<\int_{\bar{s}}^{1}(1-\ell)\,\frac{\mathrm{d}\varrho}{f(\varrho)}=\gamma^{-1}\mathcal{G}(\bar{s})<0,$

which is a contradiction. For the case $\lambda\geq 1$, we have from (3.2)
that $\ell^{\prime}(s)\overset{\eqref{eq:DefF}}{=}\frac{1}{\gamma
f(s)}[f(s)-\lambda]\overset{\ref{prop:f}}{\leq}0$ and hence no
$s_{-\infty}\in(0,1)$ satisfying (2.10) exists, thus giving the result. ∎

##  4 Linear stability of the travelling waves in two space dimensions

In this section, we analyse the linear stability of the TW solutions in two
space dimensions. It was observed numerically in [11, 17] that the TWs are
stable with respect to random perturbations, and that solutions resulting from
arbitrary initial data converge to a TW solution in a long rectangular domain.
The stability of the TWs is further investigated numerically in Section 5. In
this section, we use asymtotic expansions to prove linear $L^{1}$-stability of
the TWs under certain (suitable) assumptions.

###  4.1 Asymtotic expansion and linearization

Let $(M^{\langle{0}\rangle},S^{\langle{0}\rangle})$ denote the TWs introduced
in Theorem 2.1. For a given $L>0$, we consider the infinite-cylinder
$\Omega:={\mathbb{R}}\times(0,L)$ as the domain. Let $(x,y)\in\Omega$ denote
the spatial coordinates. Then, we consider the two dimensional version of
(1.1) in $\Omega$:

$\displaystyle\partial_{t}M=\partial_{x}[D(M)\,\partial_{x}M]+\partial_{y}[D(M)\,\partial_{y}M]+\left(f(S)-\lambda\right)M,$
(4.1a) $\displaystyle\partial_{t}S=-\gamma\,f(S)\,M.$ (4.1b)

We consider the following boundary conditions for (4.1):

$\displaystyle M(\pm\infty,y,t)=0,\;\;(D(M)\,\partial_{x}M)(\pm\infty,y,t)=0$
$\displaystyle\text{ for all }y\in(0,L)\text{ and }t>0,$ (4.2a)
$\displaystyle(D(M)\,\partial_{y}M)(x,0,t)=(D(M)\,\partial_{y}M)(x,L,t)=0$
$\displaystyle\text{ for all }x\in{\mathbb{R}}\text{ and }t>0.$ (4.2b)

The boundary condition (4.2a) generalizes (1.2) to two dimensions, whereas,
(4.2b) introduces homogeneous Neumann conditions for the lateral boundaries.
As initial condition we choose for an arbitrarily small $\varepsilon>0$,

$\displaystyle\begin{cases}M(x,y,0)&=M^{\langle{0}\rangle}(x)+\varepsilon\sum_{n=0}^{\infty}M^{\langle{1}\rangle}_{n,0}(x)\cos\left(\Lambda_{n}\,y\right),\\\\[5.0pt]
S(x,y,0)&=S^{\langle{0}\rangle}(x)+\varepsilon\sum_{n=0}^{\infty}S^{\langle{1}\rangle}_{n,0}(x)\cos\left(\Lambda_{n}\,y\right),\end{cases}\text{
where }\Lambda_{n}:=\tfrac{2\pi n}{L},$ (4.3a) and the functions
$M^{\langle{1}\rangle}_{n,0},\;S^{\langle{1}\rangle}_{n,0}\in
L^{1}({\mathbb{R}}^{-})$ are bounded, smooth, and satisfy $\displaystyle
M^{\langle{1}\rangle}_{n,0}(0)=S^{\langle{1}\rangle}_{n,0}(0)=0,\text{ and
}\lim\limits_{x\searrow-\infty}M^{\langle{1}\rangle}_{n,0}(x)=\lim\limits_{x\searrow-\infty}S^{\langle{1}\rangle}_{n,0}(x)=0\text{
for all }n\in{\mathbb{N}}_{0}.$ (4.3b)

######  Remark 4.1 (Generality of the initial condition).

The initial condition (4.3) can be generalized to include also
$\sin(\Lambda_{n}y)$ components which then covers all smooth initial
conditions having periodic boundaries at $y\in\\{0,\,L\\}$ (by Fourier series
expansion). The main result of this section (Theorem 4.1) remains unchanged
provided that $\sin(\Lambda_{n}y)$ contributions are added in (4.8). For
simplicity, we have used only the $\cos(\Lambda_{n}y)$ components here.
Observe that the initial condition (4.3) is consistent with the zero Neumann
conditions in (4.2b).

We assume that the solution of (1.1) can be described in this case by the
asymptotic expansion

$\displaystyle\begin{cases}M(x,y,t)=M^{\langle{0}\rangle}(x-vt)+\varepsilon
M^{\langle{1}\rangle}(x,y,t)+\dots,\\\
S(x,y,t)=S^{\langle{0}\rangle}(x-vt)+\varepsilon
S^{\langle{1}\rangle}(x,y,t)+\dots,\end{cases}$ (4.4)

where $v>0$ represents the wave speed, and $M^{\langle{1}\rangle}_{n}$,
$S^{\langle{1}\rangle}_{n}$ are differentiable functions such that the
boundary conditions (4.2) and initial conditions (4.3) are satisfied. Then by
inserting the expansion (4.4) in (4.1) and equating the $\varepsilon$-order
terms one obtains the system

$\displaystyle\partial_{t}M^{\langle{1}\rangle}$
$\displaystyle=\partial_{xx}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle})+\partial_{yy}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle})$
$\displaystyle\quad+(f(S^{\langle{0}\rangle})-\lambda)M^{\langle{1}\rangle}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle},$
(4.5a) $\displaystyle\partial_{t}S^{\langle{1}\rangle}$
$\displaystyle=-\gamma[f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}].$
(4.5b)

Since the TW coordinate, $\xi=x-vt$, is a more natural space-coordinate
compared to $x$ to analyse this problem, we use the coordinate transform
$(x,y,t)\mapsto(\xi,y^{\prime},t^{\prime})$, where

$\displaystyle\xi=x-vt,\;y^{\prime}=y\text{ and }t^{\prime}=t,\text{ giving
}\partial_{t}=\partial_{t^{\prime}}-v\partial_{\xi},\;\partial_{x}=\partial_{\xi}\text{
and }\partial_{y^{\prime}}=\partial_{y}.$ (4.6)

System (4.5) is modified accordingly. For simplicity, referring to
$(y^{\prime},t^{\prime})$ as $(y,t)$, we then have the following problem

$\displaystyle\partial_{t}M^{\langle{1}\rangle}-v\,\partial_{\xi}M^{\langle{1}\rangle}$
$\displaystyle=\partial_{\xi\xi}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle})+\partial_{yy}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle})$
$\displaystyle\quad+(f(S^{\langle{0}\rangle})-\lambda)M^{\langle{1}\rangle}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle},$
(4.7a)
$\displaystyle\partial_{t}S^{\langle{1}\rangle}-v\,\partial_{\xi}S^{\langle{1}\rangle}$
$\displaystyle=-\gamma[f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}].$
(4.7b)

Due to the form of the initial condition prescribed, we look for
$M^{\langle{1}\rangle}$ and $S^{\langle{1}\rangle}$ of the following form

$\displaystyle\begin{cases}M^{\langle{1}\rangle}(x,y,t)&=\sum\limits_{n=0}^{\infty}M^{\langle{1}\rangle}_{n}(\xi,t)\cos\left(\Lambda_{n}\,y\right),\\\\[5.0pt]
S^{\langle{1}\rangle}(x,y,t)&=\sum\limits_{n=0}^{\infty}S^{\langle{1}\rangle}_{n}(\xi,t)\cos\left(\Lambda_{n}\,y\right),\end{cases}\text{
with }\Lambda_{n}=\tfrac{2\pi n}{L}.$ (4.8)

Substituting (4.8) in (4.7), observing that
$(M^{\langle{0}\rangle},S^{\langle{0}\rangle})$ only depends on $\xi$, and
equating the $\cos\left(\Lambda_{n}y\right)$ terms, we have for each
$n\in{\mathbb{N}}_{0}$ the system

$\displaystyle\partial_{t}M^{\langle{1}\rangle}_{n}-v\,\partial_{\xi}M^{\langle{1}\rangle}_{n}$
$\displaystyle=\partial_{\xi\xi}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n})-\Lambda_{n}^{2}M^{\langle{1}\rangle}_{n}$
$\displaystyle\quad+(f(S^{\langle{0}\rangle})-\lambda)M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n},$
(4.9a)
$\displaystyle\partial_{t}S^{\langle{1}\rangle}_{n}-v\,\partial_{\xi}S^{\langle{1}\rangle}_{n}$
$\displaystyle=-\gamma[f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n}].$
(4.9b)

In order to satisfy the boundary and initial conditions (4.2)–(4.3) for each
$n\in{\mathbb{N}}_{0}$ one must have

$\displaystyle M^{\langle{1}\rangle}_{n}(0,t)$
$\displaystyle=S^{\langle{1}\rangle}_{n}(0,t)=0$ $\displaystyle\text{ for all
}t>0,$ (4.10a)
$\displaystyle\lim\limits_{\xi\searrow-\infty}M^{\langle{1}\rangle}_{n}(\xi)$
$\displaystyle=\lim\limits_{\xi\searrow-\infty}S^{\langle{1}\rangle}_{n}(\xi)=0,$
(4.10b) $\displaystyle M^{\langle{1}\rangle}_{n}(\xi,0)$
$\displaystyle=M^{\langle{1}\rangle}_{n,0}(\xi),\quad
S^{\langle{1}\rangle}_{n}(\xi,0)=S^{\langle{1}\rangle}_{n,0}(\xi)$
$\displaystyle\text{ for all }\xi<0.$ (4.10c)

######  Remark 4.2 (Choice of boundary conditions (4.10)).

The boundary conditions at $\xi=0$ prescribed in (4.10) imply that
$\partial_{\xi}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n})=D(M^{\langle{0}\rangle})\partial_{\xi}M^{\langle{1}\rangle}_{n}+\partial_{\xi}D(M^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n}=0$
at $\xi=0$ since $\partial_{\xi}D(M^{\langle{0}\rangle})$ is bounded. Thus,
the flux is zero at $\xi=0$ which ensures that
$(M^{\langle{1}\rangle}_{n}(t),S^{\langle{1}\rangle}_{n}(t))$, as a solution
to (4.9), can be extended to ${\mathbb{R}}^{+}$ by setting

$M^{\langle{1}\rangle}_{n}=S^{\langle{1}\rangle}_{n}=0,\text{ for all
}\xi>0\text{ and }t>0.$

Hence, (4.9) is satisfied for all $\xi\in{\mathbb{R}}$ and $t>0$. Finally, the
boundary conditions $M^{\langle{1}\rangle}_{n}=S^{\langle{1}\rangle}_{n}=0$ at
$\xi=-\infty$ is consistent with the initial conditions in (4.3) and make it
possible to have absolutely integrable solutions
$(M^{\langle{1}\rangle}_{n},S^{\langle{1}\rangle}_{n})$.

###  4.2 Stability in $L^{1}$-norm

######  Theorem 4.1 (Stability of the travelling waves).

Let
$(M^{\langle{0}\rangle},S^{\langle{0}\rangle}):{\mathbb{R}}^{-}\to[0,1]^{2}$
be a travelling wave solution satisfying (2.11)–(2.12) with a given wave-speed
$v>0$. Assume that for all $n\in{\mathbb{N}}_{0}$, a continuously
differentiable and absolutely integrable solution
$(M^{\langle{1}\rangle}_{n},S^{\langle{1}\rangle}_{n}):{\mathbb{R}}^{-}\times[0,\infty)\to{\mathbb{R}}^{2}$
exists of the problem (4.9) which satisfies the initial and boundary
conditions (4.10). Then, for any given $t>0$, one has

$\displaystyle\int_{{\mathbb{R}}^{-}}\left[|M^{\langle{1}\rangle}_{n}(t)|+\tfrac{1}{\gamma}|S^{\langle{1}\rangle}_{n}(t)|\right]+\left(\lambda+\Lambda_{n}^{2}\right)\int_{0}^{t}\int_{{\mathbb{R}}^{-}}|M^{\langle{1}\rangle}_{n}|$
$\displaystyle\leq\int_{{\mathbb{R}}^{-}}\left[|M^{\langle{1}\rangle}_{n,0}|+\tfrac{1}{\gamma}|S^{\langle{1}\rangle}_{n,0}|\right]+\int_{0}^{t}\int_{\\{M^{\langle{1}\rangle}_{n}\cdot
S^{\langle{1}\rangle}_{n}\leq
0\\}}[2f(S^{\langle{0}\rangle})|M^{\langle{1}\rangle}_{n}|-f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}|S^{\langle{1}\rangle}_{n}|].$
(4.11)

Consequently, if $\Lambda_{n}^{2}>2f(1)-\lambda>0$, then
$\int_{0}^{\infty}\int_{{\mathbb{R}}^{-}}|M^{\langle{1}\rangle}_{n}|<\infty$,
and

$\displaystyle\int_{{\mathbb{R}}^{-}}\left[|M^{\langle{1}\rangle}_{n}(t)|+\tfrac{1}{\gamma}|S^{\langle{1}\rangle}_{n}(t)|\right]\text{
is strictly decreasing with respect to }t>0.$ (4.12)

######  Remark 4.3 (Stability of the travelling waves).

$\text{If }\;L<\frac{2\pi}{\sqrt{2f(1)-\lambda}},\text{ then \eqref{eq:Decay}
holds for all }n\in{\mathbb{N}},$

thus proving stability of the TWs in two dimensions. However, for $n=0$,
$\Lambda_{n}=0$ and the stability of the TW is not guaranteed. In practice,
this means that perturbations to the TW that have fast transverse variations
decay. The case $n=0$ represents no perturbation in the transverse direction
but only in the direction of the TW. Unfortunately, only conditional stability
can be proven for this case using our analysis. Increased stability due to
transverse variations, similar to Theorem 4.1, has been studied earlier, for
example in [29]. Using numerical simulations we show in Section 5.1.1 that the
TWs are also stable in terms of large longitudinal perturbations.

###### Proof.

The proof uses the well-known $L^{1}$-contraction principle. We reproduce a
formal version here for the sake of brevity.

Let ${\rm sign}_{\varepsilon}:{\mathbb{R}}\to[-1,1]$ denote a regularised
version of the signum function for $\varepsilon>0$ with
$\mathcal{U}_{\varepsilon}$ as its primitive. More precisely,

$\displaystyle{\rm sign}_{\varepsilon}(u):=\begin{cases}1&\text{ if
}u>\varepsilon,\\\ u/\varepsilon&\text{ if }|u|\leq\varepsilon,\\\ -1&\text{
if }u<-\varepsilon,\\\ \end{cases}\text{ and
}\;\mathcal{U}_{\varepsilon}(u)=\int_{0}^{u}{\rm sign}_{\varepsilon}.$ (4.13)

Note the following properties of these functions for future use

$\displaystyle|u\,{{\rm sign}_{\varepsilon}}^{\prime}(u)|\begin{cases}<1\text{
for }|u|\leq\varepsilon,\\\ =0\text{ for }|u|>\varepsilon,\\\
\end{cases}\text{ and }\;\mathcal{U}_{\varepsilon}(u)\geq 0\text{ with
equality only for }u=0,$ (4.14a) $\displaystyle{\rm
sign}_{\varepsilon}(u)\to{\rm sign}(u),\;u\,{\rm
sign}_{\varepsilon}(u)\to|u|,\text{ and
}\mathcal{U}_{\varepsilon}(u)\to|u|\text{ pointwise as }\varepsilon\searrow
0.$ (4.14b)

We use ${\rm sign}_{\varepsilon}(M^{\langle{1}\rangle}_{n})$ as a test
function in (4.9a). Multiplying (4.9a) by ${\rm
sign}_{\varepsilon}(M^{\langle{1}\rangle}_{n})$ and integrating in
${\mathbb{R}}^{-}$ one has term by term

$\displaystyle\int_{{\mathbb{R}}^{-}}{\rm
sign}_{\varepsilon}(M^{\langle{1}\rangle}_{n})\,\partial_{t}M^{\langle{1}\rangle}_{n}=\partial_{t}\left(\int_{{\mathbb{R}}^{-}}\mathcal{U}_{\varepsilon}(M^{\langle{1}\rangle}_{n})\right)\overset{\eqref{eq:propSignB}}{\to}\partial_{t}\left(\int_{{\mathbb{R}}^{-}}|M^{\langle{1}\rangle}_{n}|\right)\text{
as }\varepsilon\searrow 0,$ (4.15a) $\displaystyle-v$
$\displaystyle\int_{{\mathbb{R}}^{-}}{\rm
sign}_{\varepsilon}(M^{\langle{1}\rangle}_{n})\,\partial_{\xi}M^{\langle{1}\rangle}_{n}=-v\int_{{\mathbb{R}}^{-}}\partial_{\xi}\,\mathcal{U}_{\varepsilon}(M^{\langle{1}\rangle}_{n})=v[\mathcal{U}_{\varepsilon}(M^{\langle{1}\rangle}_{n}(-\infty,t))-\mathcal{U}_{\varepsilon}(0)]\overset{\eqref{eq:stabBCpde}}{=}0.$
(4.15b) From the second order term, one has using integration by parts, and
the boundary conditions in (4.10) that
$\displaystyle\int_{{\mathbb{R}}^{-}}{\rm
sign}_{\varepsilon}(M^{\langle{1}\rangle}_{n})\,\partial_{\xi\xi}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n})=-\int_{{\mathbb{R}}^{-}}\partial_{\xi}\left({\rm
sign}_{\varepsilon}(M^{\langle{1}\rangle}_{n})\right)\partial_{\xi}(D(M^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n})$
$\displaystyle=-\int_{{\mathbb{R}}^{-}}{{\rm
sign}_{\varepsilon}}^{\prime}(M^{\langle{1}\rangle}_{n})D(M^{\langle{0}\rangle})|\partial_{\xi}M^{\langle{1}\rangle}_{n}|^{2}-\int_{{\mathbb{R}}^{-}}{{\rm
sign}_{\varepsilon}}^{\prime}(M^{\langle{1}\rangle}_{n})\,M^{\langle{1}\rangle}_{n}\,\partial_{\xi}\left(D(M^{\langle{0}\rangle})\right)\partial_{\xi}M^{\langle{1}\rangle}_{n}$
$\displaystyle\overset{\eqref{eq:propSignA}}{\leq}\int_{\\{|M^{\langle{1}\rangle}_{n}|<\varepsilon\\}}|\partial_{\xi}D(M^{\langle{0}\rangle})\,\partial_{\xi}M^{\langle{1}\rangle}_{n}|\to
0\qquad\text{ as }\varepsilon\searrow 0.$ (4.15c) Finally, for the source
terms, one has as $\varepsilon\to 0$,
$\displaystyle\int_{{\mathbb{R}}^{-}}{\rm
sign}_{\varepsilon}(M^{\langle{1}\rangle}_{n})\,[(f(S^{\langle{0}\rangle})-\lambda-\Lambda_{n}^{2})M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n}]$
$\displaystyle\overset{\eqref{eq:propSignB}}{\to}-(\lambda+\Lambda_{n}^{2})\int_{{\mathbb{R}}^{-}}|M^{\langle{1}\rangle}_{n}|+\int_{{\mathbb{R}}^{-}}{\rm
sign}(M^{\langle{1}\rangle}_{n})\,[f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n}].$
(4.15d)

Similarly, multiplying (4.9b) by ${\rm
sign}_{\varepsilon}(S^{\langle{1}\rangle}_{n})$, integrating in
${\mathbb{R}}^{-}$ and following the steps of (4.15) one has

$\displaystyle\partial_{t}\left(\int_{{\mathbb{R}}^{-}}|S^{\langle{1}\rangle}_{n}|\right)\leq-\gamma\int_{{\mathbb{R}}^{-}}{\rm
sign}(S^{\langle{1}\rangle}_{n})\,[f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n}].$
(4.16)

Adding (4.15)–(4.16) one thus obtains

$\displaystyle\partial_{t}\left(\int_{{\mathbb{R}}^{-}}\left[|M^{\langle{1}\rangle}_{n}|+\tfrac{1}{\gamma}|S^{\langle{1}\rangle}_{n}|\right]\right)+\left(\lambda+\Lambda_{n}^{2}\right)\int_{{\mathbb{R}}^{-}}|M^{\langle{1}\rangle}_{n}|$
$\displaystyle\leq\int_{{\mathbb{R}}^{-}}({\rm
sign}(M^{\langle{1}\rangle}_{n})-{\rm
sign}(S^{\langle{1}\rangle}_{n}))\,[f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n}].$
(4.17)

Note that if both $M^{\langle{1}\rangle}_{n},\,S^{\langle{1}\rangle}_{n}>0$ or
$M^{\langle{1}\rangle}_{n},\,S^{\langle{1}\rangle}_{n}<0$ then the right hand
side vanishes. On the other hand, if $M^{\langle{1}\rangle}_{n}>0$ and
$S^{\langle{1}\rangle}_{n}<0$ then ${\rm sign}(M^{\langle{1}\rangle}_{n})-{\rm
sign}(S^{\langle{1}\rangle}_{n})=2$ and
$f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n}=f(S^{\langle{0}\rangle})|M^{\langle{1}\rangle}_{n}|-f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}|S^{\langle{1}\rangle}_{n}|$.
Hence, the integrand on the right hand side becomes

$2(f(S^{\langle{0}\rangle})|M^{\langle{1}\rangle}_{n}|-f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}|S^{\langle{1}\rangle}_{n}|)<2f(S^{\langle{0}\rangle})|M^{\langle{1}\rangle}_{n}|-f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}|S^{\langle{1}\rangle}_{n}|.$

By symmetry, we have the same inequality when $M^{\langle{1}\rangle}_{n}<0$
and $S^{\langle{1}\rangle}_{n}>0$. Including the trivial cases of
$M^{\langle{1}\rangle}_{n}=0$ and/or $S^{\langle{1}\rangle}_{n}=0$ which in
both cases yield

$({\rm sign}(M^{\langle{1}\rangle}_{n})-{\rm
sign}(S^{\langle{1}\rangle}_{n}))\,[f(S^{\langle{0}\rangle})M^{\langle{1}\rangle}_{n}+f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}S^{\langle{1}\rangle}_{n}]\leq
2f(S^{\langle{0}\rangle})|M^{\langle{1}\rangle}_{n}|-f^{\prime}(S^{\langle{0}\rangle})M^{\langle{0}\rangle}|S^{\langle{1}\rangle}_{n}|,$

we have (4.11) by integrating (4.2) in time.

Since $f(S^{\langle{0}\rangle})<f(1)$ in ${\mathbb{R}}^{-}$, we have
$\lambda+\Lambda_{n}^{2}>2f(S^{\langle{0}\rangle})$ in (4.12). This implies
that there exists a constant $c_{1}>0$ such that for any $0<t_{1}<t_{2}$,

$\int_{{\mathbb{R}}^{-}}\left[|M^{\langle{1}\rangle}_{n}(t_{2})|+\tfrac{1}{\gamma}|S^{\langle{1}\rangle}_{n}(t_{2})|\right]+c_{1}\int_{t_{1}}^{t_{2}}\int_{{\mathbb{R}}^{-}}|M^{\langle{1}\rangle}_{n}|\leq\int_{{\mathbb{R}}^{-}}\left[|M^{\langle{1}\rangle}_{n}(t_{1})|+\tfrac{1}{\gamma}|S^{\langle{1}\rangle}_{n}(t_{1})|\right].$

Hence,
$\int_{{\mathbb{R}}^{-}}\left[|M^{\langle{1}\rangle}_{n}(t)|+\tfrac{1}{\gamma}|S^{\langle{1}\rangle}_{n}(t)|\right]$
is a decreasing function having a limit, and
$\int_{0}^{t}\int_{{\mathbb{R}}^{-}}|M^{\langle{1}\rangle}_{n}|$ is bounded
uniformly for all $t>0$. ∎

##  5 Numerical results

In this section, we verify the analytical predictions of Sections 3 and 4
numerically by computing solutions of the PDE systems (1.1) and (4.1). It is
shown that for an arbitrary initial condition, the PDE solutions indeed
converge to a profile moving with constant speed both in one and two space
dimensions, which further shows numerically the stability of the TWs. The TWs
are also obtained directly by solving (2.11) and finding the correct wave-
speed by a bisection algorithm. The TWs produced by the PDE simulations and
the bisection algorithm are shown to coincide. The algorithm is then used to
further investigate the parametric dependence of the TW profiles nd wave-speed
for a greater range of parameters. Finally, a numerical continuation approach
is considered which enables us to study the limiting cases for which the
assumptions in Theorem 2.1 are satisfied.

###  5.1 PDE simulations

Here, we solve the PDE systems (1.1) (one space dimension) and (4.1) (two
space dimensions) with expressions (1.3) for $D,\,f$, on finite rectangular
domains using homogeneous Neumann boundary conditions. The backward Euler
method with uniform time-step $\Delta t>0$ is used for the temporal
discretization. A uniform square grid with mesh-size $\Delta x>0$ is used to
discretize the domain. The PDEs are solved using the standard two-point flux
approximation finite volume method. The solution is approximated at the
centres of each grid cell. The diffusion across the cell interfaces are
approximated via arithmetic averaging as described in [10]. For time
integration we use the trapezoidal method. The large system of arithmetic
equations generated by this approach is solved via a fixed point iteration
scheme described in [17]. The default set of simulation parameters is given in
Table 1. These parameters will be used throughout this section unless stated
otherwise.

Model parameter | Symbol | Value | Reference
---|---|---|---
Motility coefficient in (1.3) | $\delta$ | $10^{-6}$ | [11]
Diffusion exponent 1 in (1.3) | $a$ | 4.0 | [11]
Diffusion exponent 2 in (1.3) | $b$ | 4.0 | [11]
Half saturation concentration in (1.3) | $\kappa$ | 0.01 | [11]
Maximum consumption rate | $\gamma$ | 0.4 | [11]
Cell-loss rate | $\lambda$ | 0.42 | [11]
Table 1: Default model parameters (their names refer to the context of
cellulolytic biofilm models) used for the simulations of (1.1), (4.1) and
(2.11). The functions $D$ and $f$ are as in (1.3). All values in the table are
dimensionless.

####  5.1.1 One dimensional results: transience and stability

In this case, our spatial domain is $(0,H)$ for $H>0$. The default initial
conditions for the system are for $h\in(0,1)$ and $d\in(0,H)$,

$\displaystyle
M(x,0)=\begin{cases}h\left(1-\tfrac{x^{4}}{d^{4}}\right)&\text{if }x\leq d,\\\
0&\text{otherwise},\end{cases}\quad S(x,0)\equiv 1.$ (5.1)

The numerical parameters used for the simulation are

$\displaystyle H=1,\;\;\Delta x=2^{-N},\;\;\Delta t=10^{-2}\,(2^{9}\Delta
x)^{2},\;\;d=\tfrac{5}{127},\;\;h=0.1.$ (5.2)

In the above, $N$ is an integer between 9 and 16. A grid independence study is
done in Appendix A and based on the result the default value of $N=14$ is
chosen.

Figure 6: Transient behaviour of the $M$-profiles obtained from the PDE
simulations. The parameters are taken from Table 1 except in (left)
$\lambda=0.42$, (center) $\lambda=0.60$, and (right) $a=b=2$.

Observe that the initial condition (5.1) is arbitrary and has no relation to
the TWs. Nevertheless, the numerical solutions develop into profiles that move
with a constant speed. This is shown in Figure 6 using three sets of
parameters. To estimate the wave-speed of the developed profile, we calculate
the interface of the biomass wave by taking the largest $x$-coordinate such
that $M(x,t)>10^{-2}$. We use $10^{-2}$ as an approximation for $0$ to avoid
any numerical noise generated due to the degeneracy. Once we have the wave
interface, we can estimate the wave speed by fitting a linear function through
the data points corresponding to the wave interface. This is done using the
built-in function fit from GNUPLOT which gives the wave-speed. Figure 7 (left)
illustrates the process for the parameter set in Table 1. To calculate the
wave profile, we plot the solution $M(x,t)$ after the profiles have developed.
To verify that a TW solution exists, we plot the wave profiles at various time
points $t_{i}$ and horizontally translate the waves by $v(t_{0}-t_{i})$. If
the wave profiles coincide, then we conclude that a TW solution exists. Unless
otherwise stated, all the PDE-simulations are verified to permit a TW
solution.

Figure 7: (left) Wave-speed calculation for the default parameters in Table 1
using linear fitting. The resulting wave-speed estimated is $v=0.01494$.
(right) The $M$-profiles horizontally shifted by $v(t_{0}-t_{i})$ for various
time points $t_{i}\in\\{56,57,58,59,60,64\\}$.

####  5.1.2 Two dimensional results: transience and stability

For the two dimensional results our spatial domain is the rectangle
$[0,H]\times[0,L]$. The initial condition prescribed is

$\displaystyle
M(x,y,0)=\begin{cases}h\left(1-\tfrac{x^{4}}{d^{4}}\right)\left(1+\tfrac{1}{5}\cos\left(\frac{2\pi
y}{L}\right)\right)&\text{if }x\leq d,\\\ 0&\text{otherwise},\end{cases}\quad
S(x,y,0)\equiv 1.$ (5.3)

The numerical parameters for this case are

$\displaystyle H=2,\;\;L=1,\;\;\Delta x=2^{-9},\;\;\Delta
t=0.01,\;\;d=\tfrac{5}{127},\;\;h=0.1.$ (5.4)

As before, the initial condition in this case is arbitrary and deviates
largely from any TW profile both in terms of its $x$ and $y$ variations. As
such, it does not satisfy the restrictions imposed in Section 4 and Theorem
4.1 for proving linear stability. Nevertheless, Figure 8 shows that the
numerical solution still develops slowly into a planar front that moves with a
constant speed. This speed, estimated in Figure 9, is very close to the wave-
speed computed at the same level of discretization, i.e., $\Delta x=2^{-9}$,
for the one dimensional case, see Table 2 of Appendix A.

Figure 8: The transient behaviour of the numerical solution of (4.1) subjected
to the initial condition (5.3). The parameter values are taken from Table 1.
Figure 9: Wave-speed estimation for the two dimensional simulation presented
in Figure 8.

###  5.2 ODE simulations

Upon verifying that the solutions indeed converge to TW solutions, we device a
more direct and faster method to compute the TW profiles.

####  5.2.1 TW algorithm

To directly find the TW profile and the corresponding wave-speed $v>0$ for a
given set of parameters, we use a semi-analytic approach that utilizes the
monotonicity property described in Theorem 3.5 with respect to $v>0$. We solve
the dynamical system (3.5) for a given $v>0$, except we solve the negation of
the vector field, i.e.

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\tau}M^{v}=-v\,[\ell(S^{v})-M^{v}],$
(5.5a)
$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\tau}S^{v}=-\tfrac{\gamma}{v}{f(S^{v})M^{v}D(M^{v})},$
(5.5b)

with $D$ and $f$ as in (1.3). The above system is solved numerically for
$\tau>0$ using the 4th-order Runge-Kutta method [27] and the following initial
condition is used for the computation:

$\displaystyle(M^{v},S^{v})(0)=(\varepsilon,1),$ $\displaystyle\text{ for
}\varepsilon=10^{-3}.$ (5.6)

This avoids the problem of starting from the degenerate equilibrium point
$(0,1)$ and gives a close approximation to the TW profile as indicated by
Theorem 3.3.

Then, the bisection method is used to determine the wave-speed for which
$(M^{v},S^{v})$ connects to the other equilibrium point $(0,s_{-\infty})$. Let
$\bar{v}>0$ be large enough so that $(M^{\bar{v}},S^{\bar{v}})$ exits the
region $\mathfrak{R}:=[0,1)\times[s_{-\infty},1]$ through the line
$\\{s=s_{-\infty}\\}$. The existence of such $\bar{v}>0$ is ensured by Theorem
3.5. Similarly, let $\underline{v}\in(0,\bar{v})$ be such that
$(M^{\underline{v}},S^{\underline{v}})$ exits $\mathfrak{R}$ through the line
$\\{m=0\\}$. Then we use the following algorithm to determine $v>0$:

######  Algorithm 5.1 (Bisection iteration to determine the travelling wave).

  1. 1.

Set $v=\frac{1}{2}|\bar{v}+\underline{v}|$. Solve $(M^{v},S^{v})$ satisfying
(5.5) and (5.6) numerically.

  2. 2.

If $(M^{v},S^{v})$ exits $\mathfrak{R}$ through $s=s_{-\infty}$, then set
$\bar{v}=v$. Else, set $\underline{v}=v$.

  3. 3.

If $|\bar{v}-\underline{v}|<10^{-4}\,v$ then stop; otherwise go to Step 1.

For the default parameter set in Table 1, Algorithm 5.1 is over 1000 times
faster than the PDE computation under the same computational set-up.

####  5.2.2 Validation of the TW algorithm using the PDE scheme

Figure 10: The comparison of the $M$-profiles obtained by using Algorithm 5.1
with the PDE simulations. The parameters are taken from Table 1 except in
(left) $\lambda=0.42$, (center) $\lambda=0.60$, and (right) $a=b=2$.

We compare the profiles obtained from Algorithm 5.1 with the profiles of the
full PDE simulations. The same three test cases are chosen as in Figure 6. The
results are shown in Figure 10. We consider the agreement to be excellent with
the wave-speed varying only by $2.72\%$ in the default case of Table 1 for
$N=14$ in (5.2), and the profiles are barely distinguishable from each other
for the two methods.

####  5.2.3 Parametric study in the limiting cases

Figure 10 shows us that when $\lambda$ changes from $0.42$ to $0.6$, the wave-
speed and the perceived width of the TW become about 20 times smaller. Hence,
PDE simulations become impractical both for large $\lambda$ values, due to
very small mesh-sizes required, and for small $\lambda$ values, due to the
increase in the required domain size. Similar problems occur for the
coefficients $a,\,b$ (see Figure 10) and the parameters $\gamma$ and $\kappa$.
Since Algorithm 5.1 is much faster compared to the PDE computations, these
cases can be better explored thoroughly using the ODE simulations.

An interesting case-study is the investigation of the conditions provided in
Theorem 2.1 for the existence of TWs. They are, Condition 1:
$\mathcal{G}(s)>0$ for all $s\in(s_{-\infty},1)$, which is also a necessary
condition due to Proposition 3.6; and Condition 2: (2.13) is satisfied. For
$\gamma$ and $\kappa$ given in Table 1, Condition 1 is satisfied if
$\lambda\geq 0.26$ whereas Condition 2 is satisfied if $\lambda\leq 0.56$.
Here we are interested in exploring the limits of $\lambda$ for which the TWs
exist. We already saw from Figures 6 and 10 that a TW solution exists for
$\lambda=0.6$, thus indicating that Condition 2 is not a necessary condition.
Using the TW algorithm, we also find TW solutions up to $\lambda=0.8$. The
profile is much narrower in this case and has a minuscule $v=0.000033$, see
Figure 11 (left). However, we were unable to find any TW solutions for
$\lambda\geq 0.9$ which is still smaller than the absolute limit of
$\lambda=1$ for which it is guaranteed that TW solutions do not exist, see
Proposition 3.6. Theorem 2.1 does not guarantee the existence of TW solutions
in these cases since the nullcline $m=\ell(s)$ does not intersect the line
$\\{m=1\\}$, see Figure 11 (right).

Figure 11: (left) The plots of $M^{v}$ and $S^{v}$ vs. $\xi=x-vt$ computed
using Algorithm 5.1 for $\lambda=0.8$. (right) The orbit $(M^{v},S^{v})$ in
the $m$–$s$ phase–plane.

On the other hand, if $\lambda$ is chosen smaller than $\lambda=0.42$, then we
see a rapid increase in wave-speed and a widening of the profile, see Figure
12. In fact, varying between $\lambda=0.38$ and $\lambda=0.36$ the wave-speed
and the profile width increase ten fold. As such, the computational time
required for the algorithm to converge increases exponentially, and we were
unable to obtain TW solutions for values of $\lambda$ smaller than $0.36$.
However, the trend in behaviour when varying $\lambda$ is evident from the
simulations.

Figure 12: The plots of $M^{v}$ and $S^{v}$ vs. $\xi=x-vt$ computed using
Algorithm 5.1. For the (left) plot $\lambda=0.38$, and the (right) plot
$\lambda=0.36$.

###  5.3 Numerical continuation method

Another approach to study parameter regimes for which the system (3.5) has a
connecting orbit is to use numerical continuation techniques [6]. The approach
is to find a connecting orbit on the center unstable manifold of the rest
state $(0,s_{-\infty})$ to the strong stable manifold of the rest state
$(0,1)$.

Following [6], we set this up as a two point boundary value problem:

$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\tau}M=v\,T\,[\ell(S)-M],$
$\displaystyle\mbox{~{}~{}~{}for $\tau\in(0,1)$},$ (5.7a)
$\displaystyle\tfrac{\mathrm{d}}{\mathrm{d}\tau}S=\tfrac{\gamma}{v}\,T\,\,f(S)\,MD(M),$
$\displaystyle\mbox{~{}~{}~{}for $\tau\in(0,1)$},$ (5.7b)
$\displaystyle(M(0),S(0))=(\varepsilon_{0},s_{\varepsilon_{0}}),$ (5.7c)
$\displaystyle(M(1),S(1)=(\varepsilon_{1},1),$ (5.7d)

where $s_{\varepsilon_{0}}$ is the smallest root of the equation

$\ell(s)=\varepsilon_{0},$ (5.8)

i.e., $\ell(s_{\varepsilon_{0}})=\varepsilon_{0}$. Note that (5.7d) is a point
on the linearized strong stable manifold for the rest state $(0,1)$. The
center unstable manifold of the rest state $(0,s_{-\infty})$ is well
approximated by $(\ell(s),s)$ for $s$ near $s_{-\infty}$ since the homological
equation for the center manifold is given by

$M\,D(M)\,\dfrac{\mathrm{d}M}{\mathrm{d}S}=\frac{v^{2}}{\gamma
f(S)}\,[\ell(S)-M],$

and $MD(M)={\cal O}(M^{1+a})$ for $M$ near zero. This is why we use the
boundary condition (5.7c). The parameter $T$ is the time of travel and is
taken to be very large in order to approximate the heteroclinic orbit. An
integral condition (see, e.g. [6, 7]) is added to (5.7) to facilitate adaptive
mesh selection when computing a branch of solutions. We use AUTO-07P [8] for
these computations.

To further illustrate the conclusions of Theorem 2.1, and to investigate the
limiting cases where the condition $\mathcal{G}(s)>0$ for all
$s\in(s_{-\infty},1)$ is violated, we compute the heteroclinic connections
using the parameter values $a=b=2$, $\kappa=1$, and $\lambda=0.3$. With these
values of $\kappa$ and $\lambda$, we find numerically that $s_{-\infty}\approx
0.1319$ and $\mathcal{G}(s)>0$ for all $s\in(s_{-\infty},1)$ provided
$\gamma>\approx 0.1093$. Hence, TWs cannot exist for $\gamma$ below this
threshold. This is exactly what is observed numerically from Figure 13 (left)
where the wave-speed $v$ is plotted against corresponding $\gamma$ values. The
horizontal asymptote shows that no solution is possible below a certain
threshold of $\gamma$ close to the predicted value. In the (right) figure we
see how the orbits vary in the phase–plane when $\gamma$ is decreased, tending
towards the line $\\{m=1\\}$ uniformly. This is similar to what was observed
in the case when $\lambda$ was lowered, see Figure 12.

Figure 13: The results for the numerical continuation method. The (left) plot
is the bifurcation diagram of $\gamma$ versus the wave speed $v$. In the
(right) plot are the corresponding solutions in the phase–plane.

In Figure 14 (left) we see a representative solution in the phase–plane along
with the nullcline, $m=\ell(s)$. The effect of using the phase-condition as
described in [7] is illustrated in Figure 14 (right) where solutions at two
different parameter values are depicted.

Figure 14: (left) Plot of one solution (blue) in the plase–plane along with
the nullcline $m=\ell(s)$ (red). Here $\gamma=0.1099$ and $v=16.227$. (right)
Plot of two different solutions as functions of $\tau$, one for
$\gamma=0.3893$ (red) and one for $\gamma=0.1099$ (blue).

##  6 Interpretation of the results

The analysis in Section 3 shows that TWs for the system (1.1) exist as
expected from the numerical experiments in [11] and observations on
cellulolytic biofilms made in [31]. Moreover, we predicted theoretically and
verified numerically that they exist for a wide range of parameters. The
resulting TWs are found to be stable even for large perturbations (see Figures
6 and 8) which agrees with the numerical observations in [11, 17]. In the
context of biofilm growth, the TWs describe the formation of crater like
structures (inverted colonies), i.e. the invasion and degradation of the
undisturbed cellulosic environment by microbes. At the leading edge of this
invading front is a microbially active layer that consumes the substrate. In
the wake of this layer, decay terms dominate over growth terms, and thus, this
region is dominated by the dynamics of how fast substrates degrade and how
fast bacteria decay if growth cannot be sustained. This complex interplay
gives the TW its distinct shape.

Our analysis also suggests that the TWs do not exist if for example, either
the cell-loss rate $\lambda\in(0,1)$ is too small or too large (Proposition
3.6). From the numerical experiments in Section 5.2, we can postulate the
reason behind the non-existence of TWs. We see from Figure 12 that as
$\lambda$ decreases, the TW profile becomes wider very quickly, and the
biofilm concentration $M$ approaches a value close to 1 in a large interval.
This is expected since small $\lambda$ implies less decay of $M$, and for very
small $\lambda$ we expect that the biofilm concentration would grow
monotonically with time and reach 1 in every point of the domain, implying
that a TW solution cannot exist. On the other hand, for $\lambda$ large, the
profiles become narrower and their amplitudes decrease as seen from Figure 11.
Hence, one expects that for $\lambda$ large enough, the initial biofilm
profile would decay to 0 monotonically with time.

The variation of the wave-speed can also be explained through these
observations. If $\lambda$ is small, then the profile is wider. Hence, the
bacteria consume the substrate faster, which results in a higher wave-speed.
For $\lambda$ large the effect is reversed.

Similarly, the effects of the consumption rate $\gamma>0$ on the existence,
profile-width and wave-speed of the TW can be explained, see Figure 13. Higher
values of $\gamma$ result in a faster consumption of the substrate, which
leads to a decrease in the production of biomass. Hence, the effects of
increasing $\gamma$ are analogous to the effects of increasing $\lambda$.

Lastly, we remark that our approach can also be applied to PDE–ODE systems
with a porous media type diffusion coefficient (i.e. $D(M)=M^{a}$), although
boundedness of $M\in[0,1)$ cannot be expected in this case. Furthermore, we
expect that the results in this paper can be extended to study TW solutions of
PDE–ODE systems with multiple substrates, see the system in [13] for instance.
Due to the structure of the TW, as a pulse with sharp front and diffusive tail
travelling at a constant speed, systems of the form (1.1) can likely be used
to model several other biological and physical processes with immobile
substrates, such as tumor growth, fungal growth, and the spreading of
wildfire.

### Acknowledgements

K. Mitra and S. Sonner would like to thank the Nederlandse Organisatie voor
Wetenschappelijk Onderzoek (NWO) for their support through the grant
OCENW.KLEIN.358. J. Hughes is supported by the Natural Sciences and
Engineering Research Council (NSERC) of Canada. The data and models used in
this paper are publicly available and have been properly referenced.

## References

  * [1] Z. Biró.  Stability of travelling waves for degenerate reaction-diffusion equations of KPP-type.  Advanced Nonlinear Studies, 2(4): 357–371, 2002.
  * [2] R.L. Burden, and J.D. Faires.  Numerical Analysis.  Brooks Cole, Boston, 2010.
  * [3] A. de Pablo, and A. Sánchez.  Global travelling waves in reaction-convection-diffusion equations.  Journal of Differential Equations 165(2): 377–413, 2000.
  * [4] A. de Pablo, and A. Sánchez.  Travelling wave behaviour for a porous-Fisher equation.  European Journal of Applied Mathematics, 9(3): 285–304, 1998.
  * [5] A. de Pablo, and J.L. Vázquez.  Travelling waves and finite propagation in a reaction-diffusion equation.  Journal of Differential Equations, 93(1): 19–61, 1991.
  * [6] E. Doedel and M.J. Friedman.  Numerical Computation of heteroclinic orbits.  Journal of Computational and Applied Mathematics, 26: 155–170, 1989\.
  * [7] E. Doedel, H. Keller, J. Kernevez,  Numerical Analysis and Control of Bifurction Problems (II): Bifurction in Infinite Dimensions.  International Journal of Bifurcation and Chaos, 1(4):745-772, 1991.
  * [8] E. Doedel, et al.  AUTO-07P: Continuation and bifurcation software for ordinary differential equations. Numerical Computation of heteroclinic orbits.  GitHub repository, 2007,  https://github.com/auto-07p/auto-07p.
  * [9] A. Dumitrache, H.J. Eberl, D.G. Allen, and G.M. Wolfaardt.  Mathematical modeling to validate on–line $\text{CO}_{2}$ measurements as a metric for cellulolytic biofilm activity in continuous–flow bioreactors.  Biochemical Engineering Journal 101: 55–67, 2015.
  * [10] H.J. Eberl, and L. Demaret.  A finite difference scheme for a degenerated diffusion equation arising in microbial ecology.  Electronic Journal of Differential Equations, 15: 77–95, 2007.
  * [11] H.J. Eberl, E.M. Jalbert, A. Dumitrache, and G.M. Wolfaardt.  A spatially explicit model of inverse colony formation of cellulolytic biofilms.  Biochemical Engineering Journal, 122: 141–151, 2017.
  * [12] M. Efendiev, and J. Muller.  Classification of existence and non-existence of running fronts in case of fast diffusion.  Advances in Mathematical Sciences and Applications, 19(1): 285, 2009\.
  * [13] M. A. Efendiev, M. Otani, and H.J. Eberl.  Mathematical Analysis of a PDE-ODE Coupled Model of Mitochondrial Swelling with Degenerate Calcium Ion Diffusion.  SIAM Journal on Mathematical Analysis, 52(1): 543–569, 2020.
  * [14] M. Ghasemi, and H.J. Eberl.  Time adaptive numerical solution of a highly degenerate diffusion-reaction biofilm model based on regularisation.  Journal of Scientific Computing, 74(2): 1060–1090, 2018.
  * [15] K. Harley, P. van Heijster, R. Marangell, G.J Pettet, and M. Wechselberger.  Existence of traveling wave solutions for a model of tumor invasion.  SIAM Journal on Applied Dynamical Systems, 13(1): 366–396, 2014.
  * [16] K. Harley, P. van Heijster, R. Marangell, G.J. Pettet, and M. Wechselberger.  Novel solutions for a model of wound healing angiogenesis.  Nonlinearity, 27(12): 2975, 2014.
  * [17] J.M. Hughes, H.J. Eberl, and S. Sonner.  A mathematical model of discrete attachment to a cellulolytic biofilm using random DEs  submitted
  * [18] J.D. Logan.  Transport modeling in hydrogeochemical systems (Vol. 15).  Springer Science and Business Media, 2001.
  * [19] K. Mitra.  Existence and properties of solutions of the extended play-type hysteresis model.  Journal of Differential Equations, 288: 118–140, 2021.
  * [20] K. Mitra, T. Köppl, I.S. Pop, C.J. van Duijn, and R. Helmig.  Fronts in two–phase porous media flow problems: The effects of hysteresis and dynamic capillarity.  Studies in Applied Mathematics, 144(4): 449–492, 2020.
  * [21] K. Mitra, A. Rätz, and B. Schweizer.  Travelling wave solutions for gravity fingering in porous media flows.  arXiv preprint, arXiv:2011.10792, 2020.
  * [22] K. Mitra, and C. J. van Duijn.  Wetting fronts in unsaturated porous media: The combined case of hysteresis and dynamic capillary pressure.  Nonlinear Analysis: Real World Applications, 50: 316–341, 2019.
  * [23] J.D. Murray.  Mathematical biology I. An introduction.  Springer, 2002.
  * [24] Y. Rohanizadegan, S. Sonner, and H.J. Eberl.  Discrete attachment to a cellulolytic biofilm modeled by an Itô stochastic differential equation.  Mathematical Biosciences and Engineering, 17: 2236–2271, 2020.
  * [25] F.A. Sánchez-Gardu$\tilde{\text{n}}$o, P.K. Maini, and M.E. Kappos.  A shooting argument approach to a sharp-type solution for nonlinear degenerate Fisher-KPP equations.  IMA journal of applied mathematics, 57(3): 211–221, 1996.
  * [26] R.A. Satnoianu, P.K. Maini, F.A. Sánchez-Gardu$\tilde{\text{n}}$o, and J.P. Armitage.  Travelling waves in a nonlinear degenerate diffusion model for bacterial pattern formation.  Discrete and Continuous Dynamical Systems–B, 1(3): 339, 2001.
  * [27] E. Süli and D.F. Mayers.  An introduction to numerical analysis.  Cambridge University Press, New York, 2003.
  * [28] C.J. van Duijn, K. Mitra, and I.S. Pop.  Travelling wave solutions for the Richards equation incorporating non-equilibrium effects in the capillarity pressure  Nonlinear Analysis: Real World Applications, 41 (C): 232–268, 2018.
  * [29] C.J. van Duijn, G.J.M. Pieters, and P.A.C. Raats.  On the stability of density stratified flow below a ponded surface.  Transport in Porous Media 127(3): 507–548, 2019.
  * [30] J.L. Vázquez. The porous medium equation: mathematical theory. Oxford University Press, 2007.
  * [31] Z.W. Wang, S.H. Lee, J.G. Elkins, and J.L. Morrell-Falvey.  Spatial and temporal dynamics of cellulose degradation and biofilm formation by Caldicellulosiruptor obsidiansis and Clostridium thermocellum.  AMB Express 1: 1–10, 2011.

##  Appendix A Grid independence study for the PDE simulations

In this section, we study how the mesh-size $\Delta x$ affects the TW solution
obtained numerically from solving the PDE. The focus is to investigate the
convergence of the wave-speed and profile as $\Delta x\to 0$. Let $v_{N}$
denote the wave-speed, numerically measured using the scheme described in
Section 5.1.1, for the $2^{N}$-grid simulation. We quantify the relative
differences between the wave-speeds across grid resolutions via the ratios

$\displaystyle{\rm ReD}_{N}=\frac{|v_{N}-v_{16}|}{|v_{16}|},\text{ and }{\rm
GRR}_{N}=\frac{{\rm ReD}_{N-1}}{{\rm ReD}_{N}}.$ (A.1)

Here, the quantity ${\rm GRR}_{N}$ stands for the grid refinement ratio which
gives a quantitative comparison between the successive relative differences.
The wave-speeds for each simulation along with the relative difference
calculations are given in Table 2.

N | Grid Size | Wave-speed $v_{N}$ | ReDN | GRRN
---|---|---|---|---
9 | $512$ | 0.019450 | 0.308242 | –
10 | $1024$ | 0.017100 | 0.150125 | 2.05
11 | $2048$ | 0.015909 | 0.070031 | 2.14
12 | $4096$ | 0.015335 | 0.031434 | 2.23
13 | $8192$ | 0.015064 | 0.013201 | 2.38
14 | $16384$ | 0.014944 | 0.005165 | 2.56
Table 2: Wave-speeds calculated from the various grid resolutions along with
the relative difference calculations ReDN and GRRN from (A.1). The wave speed
for the $2^{16}$ grid simulation is $v_{16}=0.014868$.

Figure 15: Wave profiles of the various grid refinement simulations. The
simulation parameters are given in Table 1.

It is seen from Table 2 that the wave-speeds decrease as the grid is refined
and relative differences tend to zero. These observations indicate that we
have convergence with respect to the wave-speed. Figure 15 shows the
convergence of the $M$ and $S$ profiles as the grid is made finer. A narrower
overall wave profile and a steeper wave front is observed with finer grid
resolutions. This is due to less numerical diffusion at smaller discretization
levels.

The convergence of both the wave-speed and the wave profile guarantees that
the TWs can indeed be reproduced using the numerical scheme described in
Section 5.1. Given the results in Table 2 and considerations on simulation
runtime, we use a $2^{14}$-grid in our simulations and a time step of $\Delta
t=10^{-2}(2^{9}(1/2^{14}))^{2}$ as presented in (5.2).