# Nonlinear stability of shock-fronted travelling waves under nonlocal
regularization

Ian Lizarraga and Robert Marangell School of Mathematics and Statistics, The
University of Sydney<EMAIL_ADDRESS>(corresponding author)School
of Mathematics and Statistics, The University of Sydney,
<EMAIL_ADDRESS>
###### Abstract

We determine the nonlinear stability of shock-fronted travelling waves arising
in a reaction-nonlinear diffusion PDE, subject to a fourth-order spatial
derivative term multiplied by a small parameter $\varepsilon$ that models
nonlocal regularization. Motivated by the authors’ recent stability analysis
of shock-fronted travelling waves under viscous relaxation, our numerical
analysis is guided by the observation that there is a fast-slow decomposition
of the associated eigenvalue problem for the linearised operator. In
particular, we observe an astonishing reduction of the complex four-
dimensional eigenvalue problem into a real one-dimensional problem defined
along the slow manifolds; i.e. slow eigenvalues defined near the tails of the
shock-fronted wave for $\varepsilon=0$ govern the point spectrum of the
linearised operator when $0<\varepsilon\ll 1$.

## 1 Introduction

Reaction-nonlinear diffusion partial differential equations (hereafter RND
PDEs) have recently been attracting interest as a natural mathematical model
for coherent propagation with sharp fronts. On the one hand, they can be
derived as continuum limits of discrete models of biological competition and
exclusion [26, 27]; on the other, they are now known to admit shock-fronted
travelling wave solutions [19, 20].

At the same time, the technique of smoothing shock-type solutions in nonlinear
diffusion models via high-order regularization is by now well-motivated and
studied; a prototypical example in this context is the Cahn-Hilliard model
(see e.g. Pego and Penrose’s work in [24]). Regularized nonlinear diffusion
models with advection have also been considered by Witelski in [29, 30]. An
interesting feature of high-order regularizations is that they give rise to
distinct rules for shock selection. In [30], Witelski used singular
perturbation theory to derive generalised area rules for shock selection in
nonlinear diffusion models, corresponding to an interpolated family of viscous
and nonlocal regularizations. This result has recently been extended to the
RND PDE setting; see [4].

We consider the following RND PDE with a fourth-order spatial derivative term
modelling nonlocal regularization:111We follow the convention from our earlier
work [21]: barred functions of $(x,t)$ are used for the travelling wave
dynamics in phase space, and unbarred functions of $(x,t)$ are used for
linearized quantities.

$\displaystyle\frac{\partial\bar{U}}{\partial t}$
$\displaystyle=\frac{\partial}{\partial
x}\left(D(\bar{U})\frac{\partial\bar{U}}{\partial
x}\right)+R(\bar{U})-\varepsilon^{2}\frac{\partial^{4}\bar{U}}{\partial
x^{4}}$ (1)

for $(x,t)\in\mathbb{R}\times\mathbb{R}$. The small parameter $\varepsilon\geq
0$ characterizes the strength of the regularizing effect. We suppose that the
potential function $F(\bar{U}):=\int D(\bar{U})d\bar{U}$ is nonmonotone, i.e.
$D(\bar{U})<0$ within an interval $(a,b)\subset(0,1)$, and we also suppose
that the reaction term is bistable, being pinned at $U=0$ and $U=1$.

Using geometric singular perturbation theory (GSPT), Li and his coauthors [19]
showed the existence of a one-parameter family of (smoothed) shock-fronted
travelling waves for system (1). Using a coordinate frame that follows the
travelling wave, the PDE (1) can be written as a closed four-dimensional
system of ODEs in two fast and two slow variables. The shock-fronted
travelling waves are then found as heteroclinic orbits connecting a saddle-
type fixed point at $\bar{U}=1$ to one at $\bar{U}=0$. These fixed points lie
on disjoint branches of two-dimensional saddle-type slow manifolds, and the
shock segment corresponds to a fast jump between these manifolds. Li et al.
also derive the well-known equal area rule for these waves by utilizing the
Hamiltonian structure of the (fast) layer subsystem.

The objective of this paper is to determine the nonlinear stability of this
family of shock-fronted travelling waves. Following the standard strategy, we
calculate the spectrum $\sigma(\mathcal{L})$ of a corresponding sectorial
linearized operator $\mathcal{L}$ associated with the PDE (1). Our task is to
show that except for a simple translational eigenvalue that necessarily exists
at the origin, the entire spectrum is bounded within the left-half complex
plane in a ‘nice’ way, as we now describe. The total spectrum is decomposed
into its continuous and point components as
$\sigma(\mathcal{L})=\sigma_{e}(\mathcal{L})\cup\sigma_{p}(\mathcal{L})$. Our
first result concerns the location of the continuous spectrum and its
implications for nonlinear stability: completing the analysis that we
initiated in the Concluding Remarks of [21], we show in Sec. 3 that the
linearised operator is sectorial. We can already highlight a major difference
with the viscous relaxation limit; there, the continuous spectrum was
similarly shown to be bounded inside the left-half complex plane, but it was
asymptotically vertical, obstructing sectoriality of the operator.

It then remains to show that the translational eigenvalue is simple and that
the rest of the point spectrum is bounded within the left-half complex plane,
and this is the primary focus of this paper. Following standard geometric
stability theory for travelling waves (see e.g. [16] and [25]), we
characterize values $\lambda\in\sigma_{p}(\mathcal{L})$ in the point spectrum
via bifurcations of the eigenvalue problem $(\mathcal{L}-\lambda)v=0$ posed on
the underlying linearized subspaces along the travelling wave. Specifically,
we are interested in tracking the complex two-plane subbundle of unstable
subspaces along travelling wave solutions of (1), beginning at the asymptotic
steady-state $\bar{U}=1$ of (1). An eigenvalue
$\lambda\in\sigma_{p}(\mathcal{L})$ is found whenever this bundle has a
nontrivial intersection with the corresponding stable subbundle continued from
the stable subspace at $\bar{U}=0$.

Numerical methods to locate these intersections focus on robustly
approximating a suitable Evans function, whose roots coincide with values in
the point spectrum; we refer the reader to some standard references [16, 25].
A recent technique that has found much success in numerical stability
calculations is the computation of the winding number of a Riccati-Evans
function along a suitably chosen contour in the complex plane; see [9, 18] and
references therein. Unlike a traditional Evans function, the Riccati-Evans
function is typically not analytic but at best meromorphic, and is defined by
considering the induced flow of the linear eigenvalue problem on a suitable
(chart of the) Grassmannian (in the present case, we use $\text{Gr}(2,4)$, the
Grassmannian of complex 2-planes in $\mathbb{C}^{4}$). As demonstrated by
Harley et al. in [9], blowups of the induced matrix Riccati equation (i.e.
coordinate pole singularities) can be moved by selecting charts judiciously in
the domain of interest.

Our main result (see Sec. 4) is the use of a Riccati-Evans function to
calculate the point spectrum. We provide strong evidence that the point
spectrum consists of only a simple translational eigenvalue at the origin, and
another simple real eigenvalue of negative real part. We conclude that the
family of shock-fronted travelling waves is nonlinearly stable.

An equally significant objective of this paper is to explain the output of the
Riccati-Evans calculation by demonstrating an unambiguous fast-slow splitting
of the eigenvalue problem. We couch our result within the geometric framework
developed by Alexander, Gardner, and Jones in the series [1, 7, 12]. In recent
work by the authors [21], we showed rigorously that for $\varepsilon>0$
sufficiently small, the eigenvalue problem in the (three-dimensional) viscous
relaxation case was controlled by a slow eigenvalue problem defined on the
critical manifolds for $\varepsilon=0$, together with a newly-defined jump
map. Note that the unstable subbundle in the viscous relaxation case is a line
bundle (after a suitable time reversal). In the present case, we demonstrate
numerically that the unstable 2-plane bundle splits in such a way that a
‘fast’ line subbundle always makes a trivial unstable-to-unstable connection
for $\varepsilon>0$ sufficiently small, hence contributing no eigenvalues. On
the other hand, a ‘slow’ line subbundle continues to be controlled by a slow
eigenvalue problem.

One of the advantages of producing this splitting is that the analysis of the
slow problem in the nonlocal case carries over identically from that of the
viscous case: indeed, the reduced two-dimensional linearized system defined on
the critical manifold is exactly the same. For example, we can state
immediately that there are no slow eigenvalues with nonzero imaginary part
(see Sec. 8 in [21]). By studying the linear flow on projective space, we
therefore obtain the promised reduction of the original eigenvalue problem to
a real one-dimensional problem.

We wish to contrast our results with the travelling wave stability analysis of
another four-dimensional singularly-perturbed problem with two fast and two
slow variables considered by Gardner and Jones in [7]. In particular, we
report two crucial differences. The first is that the ‘fast unstable-to-
unstable’ connections that we find do not persist in the singular limit,
instead requiring a very delicate break in the (linearized) dynamics near a
heteroclinic orbit in the intermediate fast layer as soon as $\varepsilon>0$
(see Sec. 5). The second difference is the requirement of a nontrivial jump
map to define the slow problem, which can be clearly observed in our results
(see Sec. 6). These two aspects are related, as shown in the earlier setting
considered by Gardner and Jones: specifically, an elephant trunk lemma
enforces a trivial jump map for the slow dynamics (see Theorem 5.3(b) in [7]).
However, our fast eigenvalue problem degenerates in the singular limit
$\varepsilon\to 0$. As a consequence, we cannot construct a suitable ‘fast
elephant trunk’ that Gardner and Jones use in their analysis to track the fast
unstable solution of the eigenvalue problem through the shock layer.

The geometric issue of splitting the unstable 2-plane subbundle is also
intimately related to the analytic problem of factorizing an Evans function
into fast and slow reduced components over the entire travelling wave (see
[5]). Revisiting the present problem using the analytic factorization
procedure proposed by de Rijk et al. is likely to be fruitful.

The rest of this paper is organized as follows. In Sec. 2, we review the
construction of the one-parameter family of shock-fronted travelling waves in
the PDE (1). In Sec. 3 we write down the spatial eigenvalue problems for the
linearized operator of (1) and remind the reader of previous results
concerning the essential spectrum. In Sec. 4 we define the Riccati-Evans
function and compute its roots, thereby locating the point spectrum for the
singular perturbation parameter $\varepsilon=10^{-4}$. In Sec. 5 we describe
numerical observations of the fast (un)stable-to-(un)stable line bundle
connections in the eigenvalue problem, and in Sec. 6 we demonstrate that the
slow line subbundle is well-approximated by a reduced eigenvalue problem
defined near the slow manifolds of the system, together with a jump map. In
this section we also describe how the slow eigenvalue problem defined in the
singular limit $\varepsilon\to 0$ controls the generation of eigenvalues in
the point spectrum in the ‘full’ problem when $\varepsilon>0$. We conclude
with some open questions in Sec. 7.

## 2 Existence of shock-fronted travelling waves

We begin by reviewing how Li and his coauthors [19] used GSPT to construct a
family of shock-fronted travelling waves of (1) near the singular limit
$\varepsilon=0$ (i.e. ‘small’ regularisation). We freely use well-known
definitions from Fenichel theory (see [6] and the comprehensive introductions
in [13],[17]) when referring to concepts like normal hyperbolicity and locally
invariant slow manifolds.

We define the diffusion and reaction terms in (1) concretely as

$\displaystyle D(\bar{U})$ $\displaystyle=6(\bar{U}-7/12)(\bar{U}-5/6)$ (2)
$\displaystyle R(\bar{U})$ $\displaystyle=5\bar{U}(1-\bar{U})(\bar{U}-1/5).$
(3)

These definitions are selected for consistency with [19] and [21]. We use the
standard technique of expressing the PDE (1) in terms of the frame
$(\zeta,t^{\prime})=(x-ct,t)$, where $c$ is an unknown constant parameterizing
the wavespeed, and then seeking steady-states of the resulting transformed
PDE. By using a Liénard representation for the resulting closed four-
dimensional system of ODEs (see Sec. 2.1 in [21]), we obtain the following
system in two fast variables $(\bar{U},\bar{W})$ and two slow variables
$(\bar{P},\bar{V})$:

$\displaystyle\varepsilon\dot{\bar{U}}$ $\displaystyle=\bar{W}$ (4)
$\displaystyle\varepsilon\dot{\bar{W}}$ $\displaystyle=F(\bar{U})-\bar{V}$
$\displaystyle\dot{\bar{P}}$ $\displaystyle=-R(\bar{U})$
$\displaystyle\dot{\bar{V}}$ $\displaystyle=\bar{P}-c\bar{U}.$

Here $F(\bar{U}):=\int D(\bar{U})$ denotes an integral of $D(\bar{U})$, which
we fix as

$\displaystyle F(\bar{U})$
$\displaystyle=2\bar{U}^{3}-4\bar{U}^{2}+\frac{21}{8}\bar{U}.$ (5)

In terms of the stretched travelling wave variable $\xi:=\zeta/\varepsilon$,
we obtain the following system which is equivalent to (4) when
$\varepsilon>0$:222In keeping notational consistency with [7, 21] we use the
lowercase alphabet for quantities defined w.r.t. the stretched variable $\xi$
and uppercase otherwise, i.e. $a(\xi):=A(\zeta)=A(\varepsilon\xi)$ for an
arbitrary function $A(\zeta)$.

$\displaystyle\bar{u}^{\prime}$ $\displaystyle=\bar{w}$ (6)
$\displaystyle\bar{w}^{\prime}$ $\displaystyle=F(\bar{u})-\bar{v}$
$\displaystyle\bar{p}^{\prime}$ $\displaystyle=-\varepsilon R(\bar{u})$
$\displaystyle\bar{v}^{\prime}$ $\displaystyle=\varepsilon(\bar{p}-c\bar{u}).$

Note that the dot notation denotes the derivative with respect to $\zeta$ and
the prime notation denotes the derivative with respect to $\xi$. The system
(4) (or (6)) admits three fixed points corresponding to zeroes of the bistable
reaction term $R(\bar{u})$. The fixed points $q_{-}=(0,0,0,0)$ and
$q_{+}=(1,0,c,F(1))$ are saddle-type for each $\varepsilon>0$ and $c>0$, both
of which have two-dimensional local stable and unstable manifolds. The
travelling waves that we seek are now realized as heteroclinic orbits
connecting the saddle point at $\bar{u}=1$ (via a weak unstable direction) to
the one at $\bar{u}=0$ (via a weak stable direction).

The singular limit $\varepsilon\to 0$ of systems (4) and (6) are distinct. The
system (6) limits to the fast layer problem, defined as

$\displaystyle\bar{u}^{\prime}$ $\displaystyle=\bar{w}$ (7)
$\displaystyle\bar{w}^{\prime}$ $\displaystyle=F(\bar{u})-\bar{v}$
$\displaystyle\bar{p}^{\prime}$ $\displaystyle=0$
$\displaystyle\bar{v}^{\prime}$ $\displaystyle=0.$

Figure 1: Segments of $W^{s}(\bar{Q}_{-})$ and $W^{u}(\bar{Q}_{+})$ arising
from the slow flow (9), as $c$ is varied within the range $[0.17,0.22]$. The
singular homoclinic orbit $\Gamma_{0}$ is depicted by the red curve for
$c_{0}\approx 0.1968109995$, and the corresponding fast shock is also
overlaid. The fold lines at $\bar{u}=\bar{u}_{f,\pm}$ are depicted as gray
dashed lines, and the $\bar{u}$-values $\bar{u}=\bar{u}_{\pm}$ satisfying the
equal area rule are depicted as black dashed lines.

The slow variables $(\bar{p},\bar{v})$ can be interpreted as parameters in the
singular limit. On the stretched ‘timescale,’ we obtain a two-dimensional
problem containing a two-dimensional cubic manifold of equilibria called the
critical manifold:

$\displaystyle S$
$\displaystyle:=\\{(\bar{u},\bar{w},\bar{p},\bar{v})\in\mathbb{R}^{4}:\bar{v}=F(\bar{u}),~{}\bar{w}=0\\}.$

The layer problem (7) defines the dynamics along the fast fibers away from
this critical set. We can classify the normal hyperbolicity along points of
$S$ by studying the eigenvalues of the linearization of (7). There are always
two zero eigenvalues, called the trivial eigenvalues, corresponding to the
two-dimensional tangent spaces at every point in $S$. The signs of the
remaining two nontrivial eigenvalues at a point $x\in S$ determines the normal
hyperbolicity at that point. A calculation of the Jacobian shows us that $S$
is a saddle-type critical manifold on the sets
$S_{-}:=S\cap\\{\bar{u}<u_{f,-}:=3/4\\}$ and
$S_{+}:=S\cap\\{\bar{u}>\bar{u}_{f,+}:=5/6\\}$, and is repelling when
$\bar{u}_{f,-}<\bar{u}<\bar{u}_{f,+}$. The critical manifold contains two fold
lines at $\bar{u}=\bar{u}_{f,-}$ and $\bar{u}=\bar{u}_{f,+}$ along which it
loses normal hyperbolicity via a simple zero eigenvalue crossing. Note that
the fold lines occur at roots of the diffusion $D(\bar{u})$. The travelling
waves that we consider in this paper stay far away from the fold lines, in
contrast to the viscous case [21].

As highlighted in [19], the layer problem (7) has the structure of a
Hamiltonian system with

$\displaystyle H(u,w)$ $\displaystyle=\frac{1}{2}w^{2}-G(u)+vu,$

where $G(u)$ is any integral of $F(u)$. Fast connections adjoining the outer
saddle-type critical branches exist when there are heteroclinic orbits of (7)
connecting saddle points at values $u=u_{-}<u_{f,-}$ to those at
$u=u_{+}>u_{f,+}$. Using the fact that these are level sets of the
Hamiltonian, we obtain the equal area rule

$\displaystyle\int_{u_{-}}^{u_{+}}(v-F(u))\,du$ $\displaystyle=0$

which constrains such pairs of values $u_{\pm}$ . In view of the symmetry of
the graph of the cubic potential function (5) about its inflection point at
$\bar{u}=2/3$, we find the exact expressions

$\bar{u}_{\pm}=\frac{1}{12}(8\pm\sqrt{3}).$

(a)
(b)

Figure 2: (a) Heteroclinic orbit of (4) for the parameter set
$(\varepsilon,c)=(10^{-4},0.19686)$ on a
$(\bar{u},\bar{p},\bar{v})$-projection of the four-dimensional phase space
(red curve) together with segments of the singular homoclinic orbit on the
critical manifold $S$ (black dashed curves). The saddle-type branches
$S_{\pm}$ are also identified. (b) Magnified view near the folds of $S$. The
green dashed line depicts the equal area rule and the black dashed line
depicts the critical manifold. The orbit was computed using the bvp4c
boundary-value solver in MATLAB R2021a with a relative error tolerance of
$10^{-8}$.

We now write down the singular limit of (4). As $\varepsilon\to 0$, we obtain
the two-dimensional reduced problem, constrained to the critical manifold $S$:

$\displaystyle\bar{V}$ $\displaystyle=F(\bar{U})$ (8) $\displaystyle\bar{W}$
$\displaystyle=0$ $\displaystyle\dot{\bar{P}}$ $\displaystyle=-R(\bar{U})$
$\displaystyle\dot{\bar{V}}$ $\displaystyle=\bar{P}-c\bar{U}.$

The reduced problem defines the slow dynamics on the critical manifold in the
singular limit. The constraint $\dot{\bar{V}}=D(\bar{U})\dot{\bar{U}}$, which
specifies that the reduced flow should be everywhere tangent to $S$, provides
a convenient closed representation of the reduced problem as

$\displaystyle\dot{\bar{U}}$
$\displaystyle=\frac{1}{D(\bar{U})}(\bar{P}-c\bar{U})$ (9)
$\displaystyle\dot{\bar{P}}$ $\displaystyle=-R(\bar{U}).$

Note that the saddle singularities $\bar{q}_{-}$ and $\bar{q}_{+}$ for
$\varepsilon>0$ limit to saddle-type fixed points of (9) at
$\bar{Q}_{-}=(0,0)$ resp. $\bar{Q}_{+}=(1,c)$. Solution segments of the fast
and slow reduced problems (7) and (9) can then be formally concatenated with
the prospect of finding singular heteroclinic orbits connecting $Q_{-}$ and
$Q_{+}$.

In Fig. 1, we report the existence of such an orbit $\Gamma_{0}$, which is
also transverse with respect to variation in the wavespeed parameter $c$. One
of the important achievements of GSPT is the use of Fenichel theory together
with other estimates (such as the exchange lemma and blow-up theory, if
necessary) to rigorously construct heteroclinic orbits as robust perturbations
of these singular objects. We refer the reader to [19] for the detailed
construction of a one-parameter family of heteroclinic orbits
$\\{\Gamma_{\varepsilon}\\}_{\varepsilon\in(0,\bar{\varepsilon}]}$ for the
full system (4). Here $\bar{\varepsilon}$ denotes a sufficiently small upper
bound on $\varepsilon$. A member of this family is plotted in Fig. 2.

## 3 Point and essential spectra

We now turn to the stability analysis of the family of heteroclinic orbits
$\\{\Gamma_{\varepsilon}\\}_{\varepsilon\in(0,\bar{\varepsilon}]}$. For given
$\varepsilon>0$, we represent the corresponding member of this family by the
solution
$\bar{X}(\zeta,\varepsilon)=(\bar{U}(\zeta,\varepsilon),\bar{W}(\zeta,\varepsilon),\bar{P}(\zeta,\varepsilon),\bar{V}(\zeta,\varepsilon))$
and consider perturbed solutions of (1) of the form
$\tilde{U}(\zeta,t)=\bar{U}(\zeta,\varepsilon)+\delta e^{\lambda
t}U(\zeta)+\mathcal{O}(\delta^{2})$. At linear order we find the temporal
eigenvalue problem in the linearised operator $\mathcal{L}$, defined as:

$\displaystyle\lambda U$
$\displaystyle=cU_{\zeta}+(D(\bar{U})U)_{\zeta\zeta}+R^{\prime}(\bar{U})U_{\zeta}-\varepsilon^{2}U_{\zeta\zeta\zeta\zeta}=:\mathcal{L}(U)$
(10)

After collecting terms of linear order in $\delta$ and choosing a Liénard
coordinate representation for $\mathcal{L}$ that is compatible with that of
the travelling wave phase space dynamics, we arrive at the nonautonomous
linear dynamical system

$\displaystyle\varepsilon\dot{U}$ $\displaystyle=W$ (11)
$\displaystyle\varepsilon\dot{W}$ $\displaystyle=D(\bar{U})U-V$
$\displaystyle\dot{P}$ $\displaystyle=(\lambda-R^{\prime}(\bar{U}))U$
$\displaystyle\dot{V}$ $\displaystyle=P-cU$

on the (slow) $\zeta$ timescale and

$\displaystyle u^{\prime}$ $\displaystyle=w$ (12) $\displaystyle w^{\prime}$
$\displaystyle=D(\bar{u})u-v$ $\displaystyle p^{\prime}$
$\displaystyle=\varepsilon(\lambda-R^{\prime}(\bar{u}))u$ $\displaystyle
v^{\prime}$ $\displaystyle=\varepsilon(p-cu)$

on the stretched (fast) $\xi$ timescale.

Let us denote the linear system (11) as
$Z^{\prime}=A(\xi,\lambda,\varepsilon)Z$. We define the asymptotically
constant matrices associated with $A(\xi,\lambda,\varepsilon)$, given by

$A_{\pm}(\lambda,\varepsilon):=\lim_{\xi\to\pm\infty}A(\xi,\lambda,\varepsilon)$

(note: we analogously define the matrices $a_{\pm}$ in the natural way).

Now we want to compute the spectrum of the linearised operator $\mathcal{L}$.
We begin with the essential spectrum, $\sigma_{e}(\mathcal{L})$. Via Weyl’s
essential spectrum theorem (see, for example [16]), we can characterize the
essential spectrum of $\mathcal{L}$ by determining the values of $\lambda$ for
which the signatures of $A_{\pm}(\lambda,\varepsilon)$ differ; equivalently,
the Fredholm borders which bound the essential spectrum are characterized by
the condition that one of $A_{-}(\lambda,\varepsilon)$ and
$A_{+}(\lambda,\varepsilon)$ has an imaginary (spatial) eigenvalue $ik$. The
calculations are given in the Concluding Remarks of [21] (the nonlocal case we
consider is the case $a=0$ in that computation). Summarising, we have that the
dispersion relations forming the Fredholm borders are given by:

$\lambda_{\pm}=-\varepsilon^{2}k^{4}-D(\bar{U}^{\pm})k^{2}+R^{\prime}(\bar{U}^{\pm})+ick$
(13)

where

$\bar{U}^{\pm}:=\lim_{\zeta\to\pm\infty}\bar{U}$ (14)

(and analogously for the lower case letters).

These are seen to partition $\mathbb{C}$ into five disjoint sectors,
distinguished by the signatures of $A_{\pm}(\lambda,\varepsilon)$. We call
region to the right of the essential spectrum, where both
$A_{\pm}(\lambda,\varepsilon)$ have signature $(2,2)$, $\Omega$. That is

$\Omega:=\left\\{\lambda\in\mathbb{C}\bigg{|}\textrm{sgn}(A_{-})=\textrm{sgn}(A_{+})=(-,-,+,+)\right\\}.$

We use $\mathcal{A}_{i}$ with $i=1,2,3,4$ to denote the remaining regions. We
note that the essential spectrum is contained within the half plane
$\\{x+iy\in\mathbb{C}:x\leq-1=R^{\prime}(\bar{U}^{-})\\}$, and the behavior
for large $|\lambda|$ is dominated by the quartic term,
$-\varepsilon^{2}k^{4}$ (see Figure 3 and Table 1).

Figure 3: A plot of the essential spectrum (shaded region) of the operator $\mathcal{L}$ in (10). The five regions, $\Omega$ and $\mathcal{A}_{i}$ with $i=1,2,3,4$ are shown. The bounding lines are the Fredholm borders given by the dispersion relations in (13). For the figure, $\varepsilon=10^{-5}$, and $c=.19681$, with $D$ and $R$ as in (2) and (3).  Region | $\textrm{sgn}(A_{-}(\lambda,\varepsilon))$ | $\textrm{sgn}(A_{+}(\lambda,\varepsilon)f)$
---|---|---
$\Omega$ | $(-,-,+,+)$ | $(-,-,+,+)$
$\mathcal{A}_{1}$ | $(-,-,-,+)$ | $(-,-,+,+)$
$\mathcal{A}_{2}$ | $(-,-,++)$ | $(-,-,-,+)$
$\mathcal{A}_{3}$ | $(-,-,+,+)$ | $(-,-,-,+)$
$\mathcal{A}_{4}$ | $(-,-,-,+)$ | $(-,-,-,+)$
Table 1: The signatures of the asymptotically constant matrices which define
the regions $\Omega$ and $\mathcal{A}_{i}$.

For $\lambda\in\Omega$, the asymptotic operators
$A_{\pm}(\lambda,\varepsilon)$ are hyperbolic, having eigenvalues
$\mu^{\pm}_{i}$ (with $i=1,2,3,4$) that fall within the hierarchy

$\text{Re}(\mu^{\pm}_{4})\ll\text{Re}(\mu^{\pm}_{3})<0<\text{Re}(\mu^{\pm}_{2})\ll\text{Re}(\mu^{\pm}_{1}).$

This estimate can be shown to hold uniformly along the segments of the wave
near to the slow manifolds $S_{\varepsilon}$, i.e. by studying the
linearisations of $A(\xi,\lambda,\varepsilon)$ defined on the left and right
saddle branches of the slow manifold for $|\xi|$ sufficiently large, the
corresponding eigenvalues $\mu^{\pm}_{j}(\xi)$ inherit this hierarchy. See
e.g. [7, 21] for examples of such estimates in three- and four-dimensional
problems.

Note that this hierarchy implies the existence of fast unstable and fast
stable directions for the eigenvalue problem corresponding to the eigenspaces
of $\mu^{\pm}_{1}$ and $\mu^{\pm}_{4}$, respectively, together with a
$2$-plane slow subbundle given by the span of eigenspaces corresponding to
$\mu^{\pm}_{2,3}$.

We seek to determine whether the family of traveling waves
$\Gamma(c(\varepsilon),\varepsilon)$ is nonlinearly stable, i.e. whether for
fixed $\varepsilon>0$, every solution $U(t)$ of (1) with initial data
sufficiently near to the wave will tend to a translate of the wave at an
exponential rate as $t\to\infty$. General stability theory (see e.g. [11], as
well as Prop. 2.2 and Cor. 2.1 in [1]) provides broad hypotheses for which it
is sufficient to check the linearised stability of the wave. In fact, the
operator $\mathcal{L}$ is sectorial, as we shall see, so it is sufficient to
check the spectral stability of the wave. We then have linearised and
nonlinear stability of the wave due to standard results in the theory (see for
example [1, 11, 16]).

The dispersion relations in (13) show that the essential spectrum is contained
in the left half plane. We now give a geometric characterization of the point
spectrum. We restrict our attention to the set of eigenvalue parameters
$\lambda\in\Omega$, the subset of the complex plane to the right of the
Fredholm borders of the essential spectrum. For each such $\lambda\in\Omega$,
it is a direct calculation to check that for $\varepsilon>0$ sufficiently
small, the asymptotically constant matrices $A_{\pm}(\lambda,\varepsilon)$
each have two eigenvalues of positive real part and two eigenvalues of
negative real part (in particular, note that this is true for $\lambda=0$, as
indicated in the discussion of the hyperbolic fixed points in 2). By general
theory (see [1]), we may extend an unstable 2-plane bundle
$\varphi^{+}(\zeta,\lambda,\varepsilon)$ from the saddle point point $q_{+}$
at $\bar{u}=1$, using the linearised flow (11): it is the unique bundle that
converges to $\text{span}\\{v_{1}^{+},v_{2}^{+}\\}$ as $\zeta\to-\infty$,
where $v^{+}_{1,2}$ span the eigenspaces corresponding to the unstable
eigenvalues $\mu^{+}_{1,2}$ of (11) at the steady-state $\bar{U}=1$. There is
also an analogously defined stable 2-plane bundle
$\varphi^{-}(\zeta,\lambda,\varepsilon)$ extending outward from the stable
eigendirections of the fixed point $q_{-}$ perched at $\bar{u}=0$.

The general theory (see [1]) tells us that $\lambda\in\sigma_{p}(\mathcal{L})$
if and only if $\varphi^{-}$ and $\varphi^{+}$ have an intersection at some
common value $\zeta=\zeta_{0}$ (and hence for each $\zeta\in\mathbb{R}$) .
Defining the section

$\displaystyle\Sigma$ $\displaystyle=\\{\bar{u}=2/3\\}$ (15)

in the middle of the shock layer, we can without loss of generality translate
the coordinate $\zeta$ so that $\zeta=0$ when the travelling wave $\Gamma$
intersects $\Sigma$. It is well-known that an analytic Evans function can be
defined away from the essential spectrum, whose roots (and their multiplicity)
coincide with the values in the point spectrum (and their order,
respectively); see e.g. [16]. This defines a shooting problem posed with
respect to the section $\Sigma$.

For sectorial operators it can be shown (see [1, 7, 11]) that a sufficiently
large contour $K\subset\Omega$ can be selected so that the entire point
spectrum lies inside $K$. Indeed, we now mimic the argument given in the
proofs of Proposition 2.2 and Corollary 2.1 in [1] to show that the operator
$\mathcal{L}$ is sectorial. Setting

$Y:=-\varepsilon^{2}U_{\zeta\zeta}+D(\bar{U})U,$

we can re-write (10) as

$\displaystyle\varepsilon^{2}U_{\zeta\zeta}$ $\displaystyle=D(\bar{U})U-Y,$
(16) $\displaystyle Y_{\zeta\zeta}$ $\displaystyle=\lambda
U-(R^{\prime}(\bar{U})+c)U_{\zeta}.$

Writing $y:=\zeta|\lambda|^{\frac{1}{4}}$ and
$\tilde{Y}:=Y/|\lambda|^{\frac{1}{2}}$, we can rewrite (16) as

$\displaystyle\varepsilon^{2}U_{yy}$
$\displaystyle=\frac{D(\bar{U})U}{|\lambda|^{\frac{1}{2}}}-\tilde{Y}$ (17)
$\displaystyle\tilde{Y}_{yy}$
$\displaystyle=e^{i\arg\lambda}U-\frac{(R^{\prime}(\bar{U})+c)|\lambda|^{\frac{1}{4}}U_{y}}{|\lambda|}.$

As $|\lambda|\to+\infty$, the equation (18) limits to the constant coefficient
system

$\displaystyle\varepsilon^{2}U_{yy}$ $\displaystyle=-\tilde{Y}$ (18)
$\displaystyle\tilde{Y}_{yy}$ $\displaystyle=e^{i\arg\lambda}U.$

Writing (18) as a $4\times 4$ system

$\begin{pmatrix}U\\\ U_{y}\\\ \tilde{Y}\\\
\tilde{Y}_{y}\end{pmatrix}_{y}=\begin{pmatrix}0&1/\varepsilon&0&0\\\
0&0&-1/\varepsilon&0\\\ 0&0&0&1\\\
e^{i\arg\lambda}&0&0&0\end{pmatrix}\begin{pmatrix}U\\\ U_{y}\\\ \tilde{Y}\\\
\tilde{Y}_{y}\end{pmatrix}=:B\begin{pmatrix}U\\\ U_{y}\\\ \tilde{Y}\\\
\tilde{Y}_{y}\end{pmatrix},$ (19)

we find the following (spatial) eigenvalues of $B$:

$\mu=\pm\frac{e^{i\frac{\arg\lambda}{4}}}{\sqrt{2\varepsilon}}\left(1\pm
i\right).$ (20)

The key observations here are that if $0<\gamma<\pi/2$ and
$|\arg\lambda|<\pi-\gamma$, then there are two (spatial) eigenvalues of $B$
with positive real part and two with negative real part. The corresponding
unstable subspaces form an attracting set relative to the flow induced on the
Grassmanian $Gr(2,4)$ by equation (18). Therefore, if $|\lambda|\gg 1$ is
sufficiently large, the unstable 2-plane bundle
$\varphi^{+}(\zeta,\lambda,\varepsilon)$ stays close to this attractor, and
hence can not connect to $\varphi^{-}(\zeta,\lambda,\varepsilon)$ to form a
(temporal) eigenvalue. We also observe that the spatial eigenvalues (20)
remain uniformly separated as $\varepsilon\to 0$ (in contrast to the example
in [1]), and thus we can conclude the following (c.f. Proposition 2.2 and
Corollary 2.1 in [1]):

######  Proposition 3.1.

There are positive constants $M$ and $\delta$ ($0<\gamma<\pi/2$) so that if
$|\lambda|>M$ and $|arg\lambda|<\pi-\gamma$ (with $-\pi<\arg\lambda\leq\pi$),
then $\lambda\not\in\sigma(\mathcal{L})$.

We also have the following result asserting the existence of a maximal contour
for the point spectrum.

######  Corollary 3.2.

There exists a $\beta<0$ and a simple closed curve $K$, so that
$\sigma(\mathcal{L})\cap\left\\{\lambda:\textrm{Re}\lambda>\beta\right\\}$,
the spectrum of the operator $\mathcal{L}$ with real part greater than
$\beta$, is contained in the interior of $K$.

## 4 Computing the roots of a Riccati-Evans function

We now compute the point spectrum of the operator $\mathcal{L}$. The first
step is to formulate a Riccati-Evans function to solve the corresponding
shooting problem. Following Ledoux et al. [18] and Harley et al. [9], we think
of the flows of $\varphi^{\pm}(\zeta)$ as evolving elements in the
Grassmannian $Gr(2,4)$ of complex 2-planes in $\mathbb{C}^{4}$. The
Grassmannian has the structure of a complex manifold (which can be seen via
e.g. the well-known Plücker embedding). A $2$-plane $P$ in $\mathbb{C}^{4}$
can be represented by any $4\times 2$ matrix, called a frame for $P$, whose
columns are vectors in $\mathbb{C}^{4}$ that form a basis of the 2-plane. The
Grassmannian has a covering by local coordinate patches

$U_{\mathscr{i}}=\\{Y\in Gr(k,n):Y\cap Y_{i^{o}}=\emptyset\\},$

where $\mathscr{i}=\\{(i_{1},\cdots,i_{k}):1\leq i_{1}<i_{2}<\cdots<i_{k}\leq
n$ is a multi-index and $Y_{i^{o}}$ denotes the $(n-k)$-plane in
$\mathbb{C}^{n}$ spanned by the vectors $\\{e_{j}\\}_{j\notin\mathscr{i}}$.
Each element $A\in U_{\mathscr{i}}$ has a unique matrix representation
$A_{\mathscr{i}}$ with respect to the multi-index $\mathscr{i}$ of $A$, with
the property that the corresponding $k\times k$ submatrix of $A_{\mathscr{i}}$
is the identity matrix $I_{k}$. See Sec. 2 of [18] for details.

For concreteness, let us first assume that both subbundles
$\varphi^{\pm}(\zeta)$ remain on the coordinate patch $U_{\mathscr{i}}$ with
$\mathscr{i}=\\{1,\cdots,k\\}$ for all $\zeta\in\mathbb{R}$. Then for any
$2$-plane $P$ with arbitrarily chosen frame $V$, we arrive at the following
(unique) representation on $U_{\mathscr{i}}$ via right multiplication:

$\displaystyle V=\begin{pmatrix}X(\zeta)\\\ Y(\zeta)\end{pmatrix}$
$\displaystyle\sim\begin{pmatrix}I_{k}\\\ Y(\zeta)X(\zeta)^{-1}\end{pmatrix},$

where we have decomposed the original frame for $P$ into two complex-valued
$2\times 2$ matrices $X,\,Y$. Let us now specify the representation of the
$2$-plane dynamics on $\mathcal{U}_{\mathscr{i}}$ that is induced by the
eigenvalue problem (11). First perform a block decomposition of the
nonautonomous linear operator right-hand side into $2\times 2$ matrices:

$\displaystyle\left(\begin{array}[]{cc|cc}0&1/\varepsilon&0&0\\\
D(\bar{u})/\varepsilon&0&0&-1/\varepsilon\\\
\hline\cr\lambda-R^{\prime}(\bar{u})&0&0&0\\\ -c&0&1&0\end{array}\right)$
$\displaystyle=\left(\begin{array}[]{c|c}A&B\\\ \hline\cr
C&D\end{array}\right).$ (21)

Let us write $W(\zeta)=Y(\zeta)X(\zeta)^{-1}$. After some matrix algebra (see
[9]) we arrive at the $2\times 2$ complex-valued matrix Riccati equation

$\displaystyle W^{\prime}$ $\displaystyle=C+DW-WBW.$ (22)

Suppose $T$ is an invertible $4\times 4$ matrix representing a change in
chart, with $\det(T)=1$. Letting

$\displaystyle\begin{pmatrix}X_{T}\\\ Y_{T}\end{pmatrix}$
$\displaystyle:=T\begin{pmatrix}X\\\ Y\end{pmatrix}\text{ and}$ (23)
$\displaystyle\begin{pmatrix}A_{T}&B_{T}\\\ C_{T}&D_{T}\end{pmatrix}$
$\displaystyle:=T\begin{pmatrix}A&B\\\ C&D\end{pmatrix}T^{-1}$

and defining $W_{T}=X_{T}Y_{T}^{-1}$, we arrive at the corresponding matrix
Riccati equation with respect to the new chart:

$\displaystyle W_{T}^{\prime}$
$\displaystyle=C_{T}+D_{T}W_{T}-W_{T}B_{T}W_{T}.$ (24)

In this paper we fix the transformation

$\displaystyle T$ $\displaystyle=\begin{pmatrix}-i&0&1&0\\\ 0&i&0&1\\\
0&0&1&0\\\ 0&0&0&1\end{pmatrix},$

which we selected in an attempt to minimize singularities of the matrix
Riccati equation arising from leaving the ‘natural’ chart $U_{\mathscr{i}}$
(note the complex shear transformation $u\mapsto-iu+p$ in the first row). With
this choice of chart, we find a simple pole $\lambda_{p}$ which lies near the
real axis. We use this chart for our Riccati-Evans calculations throughout,
because it appears to give robust numerical results.

Let $z_{0}\in\Gamma_{\varepsilon}\cap\Sigma$ denote the (unique) intersection
of the travelling wave with the section $\Sigma$ defined in (15). The Riccati-
Evans function on this chart is then given by

$\displaystyle
E_{T}(z_{0},\lambda)=\det(W_{T}^{s}(z_{0};\lambda)-W^{u}_{T}(z_{0};\lambda)),$
(25)

where $W_{T}^{u,s}$ denote the corresponding values of $\varphi^{+,-}$,
respectively at the point of intersection (and written in the same chart
coordinates). Using the fact that $\det(T)=1$ and the argument principle, we
can assess the roots and poles of the Riccati-Evans function (and hence deduce
the zeroes of the Evans function) by computing its winding number along
judiciously chosen contours to the right of the essential spectrum.

(a)
(b)

Figure 4: Computations of the winding of the Riccati-Evans function (25) for
the parameter set $(\varepsilon,c)=(10^{-4},0.19686)$. (a) Semicircular
contour in the right-half complex plane of radius $10^{5}$ centered at the
origin, with its diameter aligning with the imaginary axis and also containing
a small semicircular detour of radius 1 to avoid the translational eigenvalue
at $\lambda=0$ (see inset). (b) Evaluation of the Riccati-Evans function along
this contour. The image is a closed curve that does not wind around the origin
(see inset)—thus, the Riccati-Evans function has a winding number of 0 along
the chosen contour.

(a)
(b)

Figure 5: Computations of the winding of the Riccati-Evans function along the
circle (a) $B(0.1,0)$ as well as along (b) $B(0.01,-0.29)$ and $B(0.02,-0.8)$.
The corresponding windings are $1$, $-1$, and $1$, respectively, which is
consistent with the existence of roots at $\lambda_{0}=0$ and
$\lambda_{1}\approx-0.29$ and a pole at $\lambda_{p}=-0.29$.

######  Remark 4.1.

Numerical integration of the matrix Riccati equation (24) appears to be very
stable relative to working directly with the linearised system (11) or its
induced two-plane dynamics via derivations (see e.g. [15]). The simpler case
of tracking an unstable line bundle is a useful starting point to understand
this numerical robustness. The obvious strategy there is to consider the
projectivization of the eigenvalue problem and to study the dynamics in
projective space over a suitable chart. Then the ‘most unstable direction’
becomes a curve of attracting fixed points over the slow manifolds (with
respect to the corresponding family of ‘frozen’ systems associated with the
projectivization of (11), where the time variable $\zeta$ in the nonautonomous
component is frozen). See e.g. [7, 21], where this observation is crucial in
the construction of so-called relatively invariant attracting sets over the
strong unstable direction. As a heuristic, attracting sets tend to be easier
to numerically approximate.

Turning to the present case, we observe numerically that the unstable
$2$-plane bundle similarly becomes an attracting set with respect to the
matrix Riccati equation (24). Let us recall that the induced derivation

$w^{\prime}=A^{(k)}w,$

specifying the $k$-plane dynamics of a linear system $z^{\prime}=Az$, has
eigenvectors corresponding to $k$-fold wedge products of the eigenvectors of
the original system, with the corresponding eigenvalues given by the $k$-fold
sums of the appropriate eigenvalues of the original system (see e.g. [3])—in
particular, the unstable subspace formed by taking the wedge product of two
unstable eigenvectors corresponds to an eigenvector of the new induced system,
corresponding to the new largest positive eigenvalue.

Then our present numerical observations are consistent with the fact that at
the projective level, these ‘unstable eigenplanes’ become attracting fixed
points, in analogy with the line bundle case. The numerical integration of
(24) is also stable when integrating the stable subbundle backwards in time,
for the same reason. We expect that the matrix-Riccati equation can be a
useful starting point to write down the corresponding estimates for relatively
invariant attracting sets over unstable $k$-plane bundles in $\mathbb{C}^{n}$.
$\triangle$

As shown in Fig. 4, an evaluation of the winding number of the Riccati-Evans
function (25) on a large contour in the right-half plane shows that there
appear to be no eigenvalues of positive real part. Our function appears to map
increasingly large semicircular contours onto a compact set bounded in the
left-half complex plane, and thus it does not appear possible to generate a
winding around the origin by choosing larger contours with the same shape.

We also compute a counterclockwise winding of $+1$ along large circular
contours containing the origin and to the right of the essential spectrum. In
view of the calculation in Fig. 4, we decide to investigate the winding along
much smaller circles centered on the real line near the origin. As shown in
Fig. 5, we evaluate the winding of the Riccati-Evans function about two small
circles $\lambda_{0}=0$ and $\lambda_{1}\approx-0.8$ to find an index of $+1$,
which is consistent with the existence of simple roots near these values. We
also find a pole $\lambda_{p}$ of order 1 (i.e. having a winding of $-1$),
with $\text{Re}(\lambda_{1})<\text{Re}(\lambda_{p})<\text{Re}(\lambda_{0})$.
The sum of the winding numbers of $\lambda_{0}$, $\lambda_{p}$, and
$\lambda_{1}$ gives a total index of $+1$, which is consistent with the
winding number evaluated about the large circles.

Further investigations along different choices of contour did not reveal the
existence of any other roots or poles. We conclude that the point spectrum
consists only of the simple eigenvalues $\lambda_{0},\,\lambda_{1}$ with

$\displaystyle\lambda_{0}$ $\displaystyle=0\text{ and}$ (26)
$\displaystyle\text{Re}(\lambda_{1})$ $\displaystyle\approx-0.8.$

Thus, the corresponding shock-fronted travelling wave is nonlinearly stable
for the singular perturbation value $\varepsilon=10^{-4}$. Investigations of a
few different small values of $\varepsilon$ ranging from $10^{-4}$ to
$10^{-2}$ gave very similar numerical results.

## 5 Fast (un)stable-to-(un)stable connections

We now show how the results in the previous section can be explained via a
fast-slow splitting of the (un)stable $2$-plane bundles into fast and slow
line subbundles. Such splittings were analyzed in great depth in the papers
[1, 7]; we relate the geometric construction here. The fundamental observation
is the topological characterization of the point spectrum within a contour via
the first Chern class of a complex vector bundle over a sphere called the
augmented unstable bundle. The sphere is formed by extending the contour along
the time direction to define an infinite cylinder, and then applying a
compactification; this sphere is endowed with the $2$-plane unstable bundle
defined using the flow of the eigenvalue problem.

The first Chern number is robust to small deformations of the complex vector
bundle and also distributes over Whitney sums of line bundles [2]. The latter
property potentially allows one to compute the windings of a $2$-plane bundle
along the contour $K$ by summing the windings of fast and slow line bundles
along $K$, and the former potentially allows a straightforward computation of
these ‘sub-windings’ for sufficiently small values of $\varepsilon>0$, by
means of reduced eigenvalue problems.

Our first task is demonstrate the existence of fast unstable-to-unstable
connections of line bundles for open sets of $\lambda\in\Omega$ and for
$\varepsilon>0$ sufficiently small. Recall the hyperbolic fast-slow eigenvalue
splittings $\mu_{i}^{\pm}$ of the eigenvalue problem near the tails of the
wave. We say that there is a fast unstable-to-unstable connection for
$\lambda$ if the fast unstable line bundle $\hat{e}_{1}$, defined as the
unique bundle that converges to the eigenspace of the eigenvalue of largest
real part $\text{Re}(\mu^{+}_{1})>0$ of $a_{+}(\lambda,\varepsilon)$ as
$\zeta\to-\infty$ under the flow of the eigenvalue problem, converges to the
eigenspace of the eigenvalue of largest real part $\text{Re}(\mu^{-}_{1})>0$
of $a_{-}(\lambda,\varepsilon)$ as $\zeta\to+\infty$. We can analogously
define a fast stable-to-stable connection $\hat{e}_{4}$ by reversing time and
asking for connections from the eigenspace of the smallest eigenvalue
$\mu^{-}_{4}<0$ for $A_{-}(\lambda,\varepsilon)$ to that corresponding to
$\mu^{+}_{4}<0$ for $A_{+}(\lambda,\varepsilon)$, respectively.

An important reduced system in this context is the fast reduced eigenvalue
problem defined from a singular limit of (12), i.e.

$\displaystyle u^{\prime}$ $\displaystyle=w$ (27) $\displaystyle w^{\prime}$
$\displaystyle=D(\bar{u})u.$

This singular limit is defined along the (singular) shock-front of the wave.
When the fast reduced eigenvalue problem has a fast unstable-to-unstable
connection for $\lambda\in\Omega$, a key estimate known as the elephant trunk
lemma (see [7]) allows this connection to be extended along the entire wave
for small values of $\varepsilon>0$. Unfortunately, since $\lambda$ does not
appear in our fast reduced eigenvalue problem, it degenerates so that there is
always a fast unstable-to-stable connection in the singular limit!333This is
true because (27) is the variational equation of the layer vector field (7)
along the singular shock front connecting two saddle points, and therefore
automatically has a nontrivial solution that vanishes exponentially quickly as
$\xi\to\pm\infty$. This motivates a numerical analysis to test whether this
connection indeed breaks into the desired fast unstable-to-unstable connection
when $\varepsilon>0$.

(a)

Figure 6: A depiction of the fast unstable-to-unstable connection (green) and
a fast stable-to-stable connection (magenta) when
$(\varepsilon,c,\lambda)=(10^{-4},0.19686,0)$, using the projectivization
$\beta_{1}=w/u$. The gray dashed lines denote the jump values
$\bar{u}=\bar{u}_{\pm}$ and the black dashed line segment depicts the singular
fast unstable-to-stable connection. Different choices of $\lambda$ give very
similar graphs.

In Fig. 6 we depict an example of such fast (un)stable-to-(un)stable
connections when $\lambda=0$. We choose a chart of projective space
$\mathbb{CP}^{3}$ with $(\beta_{1},\beta_{2},\beta_{3})=(w/u,p/u,v/u)$ to
depict the dynamics, because in the singular limit of (11) we can write down
the exact expressions

$\displaystyle v^{s}(\bar{u})=\begin{pmatrix}\sqrt{D(\bar{u})}\\\
-D(\bar{u})\\\ 0\\\ 0\end{pmatrix}$ $\displaystyle\text{ and
}v^{u}(\bar{u})=\begin{pmatrix}\sqrt{D(\bar{u})}\\\ D(\bar{u})\\\ 0\\\
0\end{pmatrix}$

characterizing the fast stable and unstable eigendirections for the
eigenvalues $\mu_{4}(\bar{u})=-\sqrt{D(\bar{u})}$ and
$\mu_{1}(\bar{u})=+\sqrt{D(\bar{u})}$, respectively. With respect to the chart
representation, we have $v^{s}(\bar{u})=(-\sqrt{D(\bar{u})},0,0)$ and
$v^{u}(\bar{u})=(\sqrt{D(\bar{u})},0,0)$. This accounts for the nearly linear
dynamics away from the jump in Fig. 6. We can also read off the kind of
connection made by checking whether the trajectory stays above $\beta_{1}>0$
(unstable-to-unstable) or below $\beta_{1}<0$ (stable-to-stable) away from the
shock.

######  Remark 5.1.

Recall that the case $\lambda=0$ corresponds to the variational problem, where
we already know there is an eigenvalue! The relevant connection in this case
occurs within the slow subbundle instead (see the discussion below (11)). We
discuss the relevant slow eigenvalue problem in the next section. $\triangle$

We test the possibility of a fast unstable-to-stable connection (i.e. a fast
eigenvalue) by defining a proxy fast reduced Riccati-Evans function for
$\varepsilon>0$ on the section $\Sigma$ as follows. Because we consider line
bundles in $\mathbb{C}^{3}$, our strategy is to consider a projectivisation of
the system (11) on a chart with $u\neq 0$, with
$(\beta_{1},\beta_{2},\beta_{3})=(w/u,p/u,v/u)$, and to define

$\displaystyle E_{f}^{\Sigma}(\lambda)$
$\displaystyle:=\beta_{1}^{u}(z_{0},\lambda)-\beta_{1}^{s}(z_{0},\lambda),$
(28)

where $z_{0}$ is defined on the section $\Sigma$ as before. This function is
‘Evans-like’ in the sense that if there are no roots of this function within a
contour, then it is hopeless for a fast unstable-to-stable connection to be
made for values of $\lambda$ within the contour: such a connection is
necessarily unique. We highlight, however, that a zero of $E_{f}^{\Sigma}$
does not necessarily imply the existence of a fast eigenvalue: the values of
$\beta_{2},\,\beta_{3}$ must also coincide. Regardless, it is meaningful to
test the winding number of $E_{f}^{\Sigma}(\lambda)$ evaluated along contours
in $\mathbb{C}$ in order to efficiently test candidates $\lambda$ for fast
connections.

We found only one simple zero of $E_{f}^{\Sigma}$ at the value $\lambda\approx
3718.025$. However, at this value it can be computed directly that
$\beta_{2}^{u}(z_{0},\lambda)-\beta_{2}^{s}(z_{0},\lambda)\approx-5.08$, and
thus there is still a fast unstable-to-unstable connection made at this value.
Using (25), we can verify that there is indeed no point spectrum in this
parameter region.

We found that the additional fast direction endowed by nonlocal regularization
(compared to viscous relaxation) appears to give us only trivial fast unstable
(resp. fast stable) connections, and thus the winding of the Evans function
for this two-plane bundle with this extra fast direction is not the same as
that without it, and this extra fast direction does not contribute to the
point spectrum.

## 6 The slow eigenvalue problem

Having shown that the fast connections do not contribute to the point
spectrum, we now demonstrate that a slow eigenvalue problem, defined near the
tails of the shock-fronted travelling wave, controls the point spectrum of the
full eigenvalue problem when $\varepsilon>0$ is sufficiently small. In analogy
to the definition of (27), we take the singular limit $\varepsilon\to 0$ of
(12) to obtain the reduced slow eigenvalue problem

$\displaystyle\dot{P}$
$\displaystyle=(\lambda-R^{\prime}(\bar{U}))(V/D(\bar{U}))$ (29)
$\displaystyle\dot{V}$ $\displaystyle=P-c_{0}V/D(\bar{U})$

defined with respect to the linearised constraint $V=D(\bar{U})U$. The reduced
slow eigenvalue problem is defined on the outer saddle-type branches of the
critical manifold $S_{0}$, in particular along the singular limit of the
travelling wave lying on these branches (i.e., the left and right slow
segments of the red trajectory depicted in Fig. 1). Reduced eigenvalue
problems of this sort appear already in [7]. We also point the reader to the
derivation in the viscous case [21]: we note that the reduced problem in the
present nonlocal case is topologically equivalent (modulo desingularisation
and a trivial orientation reversal in the ‘time’ variable $\zeta$).

For $\lambda\in\Omega$, the linear operator for the reduced slow eigenvalue
problem remains hyperbolic near the end states $\bar{U}=0$ and $\bar{U}=1$,
and so we retain a well-defined geometrical characterization of the eigenvalue
problem in terms of finding an unstable-to-stable connection. Specifically, we
can track the reduced slow unstable line bundle $\varphi^{+}_{s}$ from
$\bar{U}=1$ and check whether it connects with the reduced slow stable line
bundle $\varphi^{-}_{s}$ emanating from $\bar{U}=0$.

There is one technical issue to settle: since the reduced slow eigenvalue
problem is defined on disjoint segments of the critical manifold, we must also
define a jump map which specifies how data at the beginning of the shock on
$u=u_{+}$ is transported to $u=u_{-}$ at the end of the shock. This map is
derived as a singular limit in Sec. 7.3 of [21]: via an identical procedure,
we obtain the jump map $J_{\lambda}:\mathbb{C}^{2}\to\mathbb{C}^{2}$ given by

$\displaystyle J_{\lambda}\begin{pmatrix}P\\\ V\end{pmatrix}$
$\displaystyle=\begin{pmatrix}1&\frac{R(\bar{U})_{+}-R(\bar{U}_{-})-\lambda(\bar{U}_{+}-\bar{U}_{-})}{cU_{+}-\bar{P}_{+}}\\\
0&\frac{c\bar{U}_{+}-\bar{P}_{+}}{cU_{-}-\bar{P}_{-}}\end{pmatrix}\begin{pmatrix}P\\\
V\end{pmatrix}.$ (30)

(a) (b)

Figure 7: Slow eigenvalues (fast unstable-to-stable connections) of the
reduced slow eigenvalue problem (29)–(30) when (a) $\lambda_{0}=0$ and (b)
$\lambda_{1}\approx-0.80031$. Red point: attracting fixed point corresponding
to projectivized unstable direction at $\bar{U}=1$; black point: saddle fixed
point corresponding to projectivized stable direction at $\bar{U}=1$. The red
segments depict the flow of (29) and the intermediate magenta segments depict
the output of jump map (30).

Note that the jump map is linear, which implies that we can write down a well-
defined projectivization of the map. In the framework of [1, 7], the jump map
serves as a clutching function over the equator of the sphere when defining
the reduced slow line bundle. We also highlight that the jump map is defined
for $\lambda=0$, i.e. it is compatible with the evaluation of the reduced
vector field (9) at either end of the jump point.

Figure 8: Demonstration of a slow unstable-to-unstable connection for
$\lambda=100$. Red curve segments: solutions of the boundary-value problems
beginning on the weak unstable eigendirection (i.e. corresponding to the
eigenspace of $\mu^{+}_{2}$) for the saddle-point at $\bar{U}=1$ with the
section $\Sigma=\\{\bar{U}=2/3\\}$, and beginning on the weak stable
eigendirection (i.e. corresponding to the eigenspace of $\mu^{-}_{3}$ for the
saddle-point at $\bar{U}=0$) with that section. Red point: attracting fixed
point corresponding to projectivized unstable direction at $\bar{U}=1$; black
point: saddle fixed point corresponding to projectivized stable direction at
$\bar{U}=1$. Black dashed segments: solutions of the projectivization of (29)
on the chart $S=P/V$ (with $V\neq 0$) together with output from the jump map
(30). The red segments were computed using the bvp4c boundary-value solver
package in MATLAB, with a relative error tolerance of $10^{-8}$.

Via an identical analysis to that in Sec. 8 of [21], we may conclude that the
reduced slow eigenvalue problem (29) together with the jump map (30): (1) does
not have eigenvalues with nonzero imaginary part, and (2) does not have any
eigenvalues of positive real part. Because of these facts, we may restrict our
attention to the one-dimensional problem defined by projectivization of (29)
(by e.g. taking the chart $S=P/V$ with $V\neq 0$) and restricting the
parameter set to only real $\lambda$.

By direct computation (see Fig. 7) we find only two eigenvalues:

$\displaystyle\lambda_{0}$ $\displaystyle=0\text{ and}$ (31)
$\displaystyle\lambda_{1}$ $\displaystyle\approx-0.80031.$

By the analysis in Sec. 8 of [21], we may immediately conclude that both
eigenvalues are simple. Alternatively, we can also define a section, say
$\Sigma_{R}=\\{\bar{U}=0.7\\}$, and define a slow reduced Evans function

$\displaystyle E^{\Sigma_{R}}_{s}(\lambda)$
$\displaystyle:=S^{+}(z_{0},\lambda)-S^{-}(z_{0},\lambda),$

where $S^{+}(z_{0},\lambda)$ denotes the location of the unstable slow
subbundle $\varphi^{+}$ at the (first) intersection $z_{0}$ of the singular
limit of the travelling wave with $\Sigma_{R}$ (and analogously for
$S^{-}(z_{0},\lambda)$). These computations were carried out in Sec. 8 of [21]
for slightly different jump conditions and wavespeed, so we exclude them here
while pointing out that the numerical output is equivalent.

We emphasize that the slow eigenvalue count (31) is in excellent agreement
with the eigenvalue count of the ‘full’ problem (26) which was obtained with
Riccati-Evans computations. In fact, we can show even more: by using a
boundary value problem solver, we can set up a matching problem for the slow
line subbundles defined by the system (11) on the section $\Sigma$ and show
that they are in excellent agreement with the reduced dynamics (29)–(30). We
depict one such computation in Fig. 8, where we choose a real value
$\lambda=100$ for convenience.

We conclude our analysis by outlining the argument which relates the point
spectrum of the full problem with the point spectra of the fast and slow
reduced problems, according to the topological characterization given in [1,
7]. Fixing a large enough contour $K$ so that the entire point spectrum lies
inside it, the associated augmented unstable bundle of the full problem, which
we denote $\mathcal{E}_{\varepsilon}(K)$, can be split into a Whitney sum of
fast and slow line subbundles as

$\displaystyle\mathcal{E}_{\varepsilon}(K)$
$\displaystyle=\mathcal{E}^{f}_{\varepsilon}(K)\oplus\mathcal{E}^{s}_{\varepsilon}(K)$
(32)

when $\varepsilon>0$ is sufficiently small. It is well-known (see [2]) that
the first Chern number evaluation distributes over line bundle decompositions
as

$\displaystyle c_{1}(\mathcal{E}_{\varepsilon}(K))$
$\displaystyle=c_{1}(\mathcal{E}^{f}_{\varepsilon}(K))+c_{1}(\mathcal{E}^{s}_{\varepsilon}(K)),$

and by a theorem of Alexander, Gardner, and Jones [1], the first Chern number
of the augmented unstable bundle is equal to the eigenvalue count (with
multiplicity) of the linearised operator inside $K$. We can now focus on
computing separately the first Chern number of the fast and slow line bundles.

Our numerical observations (e.g. Fig. 8) suggest the existence of a
deformation

$\displaystyle\mathcal{E}^{s}_{\varepsilon}(K)$
$\displaystyle\cong\mathcal{E}^{s}_{0}(K),$

for $\varepsilon>0$ sufficiently small, where $\mathcal{E}^{s}_{0}(K)$ is a
reduced augmented line bundle defined by using the reduced eigenvalue problem
(29) together with a clutching operation defined from the jump map (30); see
rigorous examples of such constructions in [7, 21]. The resulting slow
eigenvalues can be counted as described in this section, i.e.
$c_{1}(\mathcal{E}^{s}_{\varepsilon}(D))\leq 2$ depending on the contour $D$
chosen.

It remains to evaluate $c_{1}(\mathcal{E}^{f}_{\varepsilon}(K))$ for
sufficiently small values of $\varepsilon>0$, and here we observe a surprising
departure from the techniques developed in [1, 7, 12]. We can ask whether
there exists a similar deformation onto some singular fast line bundle
$\mathcal{E}^{f}_{0}(K)$, i.e.

$\displaystyle\mathcal{E}^{f}_{\varepsilon}(K)$
$\displaystyle\cong\mathcal{E}^{f}_{0}(K).$

In [7], the singular fast line bundle is defined by finding singular fast
unstable-to-unstable connections in the fast reduced eigenvalue problem, and
then an elephant trunk lemma is proven to show that these connections persist
along the entire wave when $\varepsilon>0$. But our fast reduced eigenvalue
problem (27) degenerates! Nonetheless, it appears to be the case that
$c_{1}(\mathcal{E}^{f}_{\varepsilon}(K))=0$ for all sufficiently small
$\varepsilon>0$, in view of our analysis in Sec. 5.

This completes our analysis of the fourth-order nonlocal problem: we have

$(c_{1}(\mathcal{E}_{\varepsilon}(K))=2$

for all $\varepsilon>0$ sufficiently small, where the first Chern number is
exactly accounted for by the simple translational eigenvalue at the origin,
together with a nontrivial eigenvalue $\lambda_{1}$ with
$-1<\text{Re}(\lambda_{1})<0$. This nontrivial eigenvalue can be furthermore
be well-approximated by the corresponding slow eigenvalue $\lambda_{1}$ shown
in Fig. 7. Thus, the family of travelling waves
$\\{\Gamma(c(\varepsilon),\varepsilon)\\}_{\varepsilon\in(0,\bar{\varepsilon}]}$
is nonlinearly stable for each sufficiently small $\varepsilon>0$.

## 7 Concluding remarks

We have given a comprehensive analysis of the nonlinear (asymptotic) stability
of a family of shock-fronted travelling waves arising in the regularized RND
PDE (1). This result is shown by verifying that the linearised operator is
sectorial Sec. 3, and then calculating the point spectrum using a Riccati-
Evans function (Sec. 4). A fast-slow splitting of the corresponding eigenvalue
problem was also furnished (Secs. 5–6), allowing us to explain the emergence
of each eigenvalue and its order.

It would be interesting to analyze such fast-slow splittings in more general
models that incorporate both nonlocal and viscous regularizations as described
in [4], and in particular whether to assess whether the nonlinear stability of
the wave persists for all combinations of these two regularizations. The
corresponding stability result for the viscous relaxation limit in [21] is
slightly weaker because of the asymptotically vertical nature of the essential
spectrum; on the other hand, there are parameter regimes in the generalised
model for which viscous shocks can potentially form, and yet the linearised
operator is still sectorial. It is also of interest to determine to what
extent slow eigenvalue problems are responsible for the generation of
eigenvalues in the point spectrum for other more exotic choices of
regularization.

In view of the eigenvalue computation using first Chern numbers in the
previous section, further analysis (e.g. in the style of [7, 21]) that
explains the behavior of $\mathcal{E}_{\varepsilon}^{f}(K)$ for small values
of $\varepsilon$ is desirable. We suggest that a trivial ‘singular fast
bundle’ $\mathcal{E}^{f}_{0}(K)$ should exist, but it should be constructed on
the $\zeta$-timescale (and not on the $\xi$-timescale as is done in [7]), i.e.
it should be constructed using a clutching procedure, similarly to what is
done for the singular slow bundle. A rigorous approximation theorem for the
fast bundle (i.e. a proof that
$\mathcal{E}^{f}_{\varepsilon}(K)\cong\mathcal{E}^{f}_{0}(K)$ for
$\varepsilon>0$ sufficiently small) will require new estimates that do not use
elephant trunk lemma-style theory. We suggest that new estimates of exchange
lemma-type (see [21]) may prove useful; the standard $(k+\sigma)$–exchange
lemma is ideally suited to tracking plane bundles near the center-unstable
manifolds of saddle-type slow manifolds (see e.g. [14, 15]). Melnikov-type
analysis might also be useful to describe transversal breaking of degenerate
unstable-to-stable connections; see Lemma 1 and Theorem 1, part (ii) in [28]
for such a result in the case of planar viscous shockwaves.

## Acknowledgements

Both authors acknowledge support by the Australian Research Council Discovery
Project grant DP200102130. We would like to thank Björn de Rijk for a useful
discussion about analytical factorizations of Evans functions.

## References

  * [1] J. Alexander, R. Gardner, and C. Jones, A topological invariant arising in the stability analysis of travelling waves, J. reine angew. Math. 410 (1990), 167–212.
  * [2] M.F. Atiyah, K-Theory, W.A. Benjamin, Inc. (1967).
  * [3] N. Bourbaki, Algebra I, Springer-Verlag, Berlin (1989).
  * [4] Bronwyn Bradshaw-Hajek, Ian Lizarraga, Robert Marangell, and Martin Wechselberger. A geometric singular perturbation analysis of regularised reaction-nonlinear diffusion models including shocks, Proceedings of 47th Sapporo Symposium on Partial Differential Equations, (2022), pp. 53–64.
  * [5] B. de Rijk, A. Doelman, J.D.M. Rademacher, Spectra and stability of spatially periodic pulse patterns: Evans function factorization via Riccati transformation, SIAM J. Math. Anal. 48 (2016), 61–121.
  * [6] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations, J. Differential Equations 31 (1979), pp. 53–98.
  * [7] JR. Gardner and C. Jones, Stability of Travelling Wave Solutions of Diffusive Predator-Prey Systems, Trans. AMS 327 (1991), pp. 465–524.
  * [8] K.E. Harley, P. van Heijster, R. Marangell, G.J. Pettet, and M. Wechselberger, Novel Solutions for a model of wound healing angiogenesis, Nonlinearity 27 (2014), pp. 2975-3003.
  * [9] K.E. Harley, P. van Heijster, R. Marangell, G.J. Pettet, T.V. Roberts, and M. Wechselberger, (In)stability of Travelling Waves in a Model of Haptotaxis, SIAM J. Appl. Math. 80 (2020), pp. 1629–1653.
  * [10] A. Hatcher, Vector Bundles and K-Theory, unpublished manuscript version 2.2 available online at https://pi.math.cornell.edu/~hatcher/VBKT/VB.pdf (2017).
  * [11] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Mathematics, Springer-Verlag, Berlin (1981).
  * [12] C.K.R.T. Jones, Stability of the Travelling Wave Solution of the Fitzhugh-Nagumo System, Trans. AMS 286 (1984), pp. 431–469.
  * [13] C.K.R.T Jones, Geometric Singular Perturbation Theory, Lecture Notes in Mathematics, Dynamical Systems (Montecatini Terme), Springer, Berlin (2015).
  * [14] C.K.R.T. Jones and N. Kopell, Tracking Invariant Manifolds With Differential Forms in Singularly Perturbed Systems, J. Differential Equations (1994), pp. 64–88.
  * [15] C.K.R.T. Jones and S.K. Tin, Generalized Exchange Lemmas and Orbits Heteroclinic to Invariant Manifolds, DCDS-S (2009), pp. 967–1023.
  * [16] T. Kapitula and K. Promislow, Spectral and Dynamical Stability of Nonlinear Waves, Springer. Applied Mathematical Sciences Vol 185 (2013).
  * [17] C. Kuehn, Multiple Time Scale Dynamics, Springer International Publishing Switzerland (2015).
  * [18] V. Ledoux, S.J.A. Malham and V. Thümmler, Grassmannian spectral shooting, Math. Comp. 79 (2010), pp. 1585–1619.
  * [19] Y. Li, P. van Heijster, M.J. Simpson, and M. Wechselberger, Shock-fronted travelling waves in a reaction-diffusion model with nonlinear forward-backward-forward diffusion, Physica D 423 (2021), 132916.
  * [20] Y. Li, P. van Heijster, R. Marangell, and M.J. Simpson, Travelling wave solutions in a negative nonlinear diffusion–reaction model, J. Math. Biol. 81 (2020), pp. 1495–1522.
  * [21] I. Lizarraga and R. Marangell, Spectral stability of shock-fronted travelling waves under viscous relaxation, in review, arXiv:2208.10064 (2022), 70 pages.
  * [22] A. Novick-Cohen, R.L. Pego, Stable patterns in a viscous diffusion equation, Trans. AMS 324 (1991), pp. 331–351.
  * [23] V. Padrón, Effect of aggregation on population recovery modeled by a forward-backward pseudoparabolic equation, Trans. AMS 356 (2004), pp. 2739–2756.
  * [24] R.L. Pego, O. Penrose, Front migration in the nonlinear Cahn-Hilliard equation, P. Roy. Soc. A-Math. Phys. 422 (1989), pp. 261–278.
  * [25] B. Sandstede, Stability of travelling waves, Handbook of Dynamical Systems II, North-Holland (2002), 983–1055.
  * [26] M.J. Simpson, K.A. Landman, and B.D. Hughes, Cell invasion with proliferation mechanisms motivated by time-lapse data, Physica A 389 (2010), pp. 3779–3790.
  * [27] M.J. Simpson, K.A. Landman, B.D. Hughes, and A.E. Fernando, A model for mesoscale patterns in motile populations, Physica A 389 (2010), pp. 1412–1424.
  * [28] P. Szmolyan and H. Freistühler, Spectral stability of small-amplitude viscous shock waves in several space dimensions, Arch. Ration. Mech. Anal. 195-2 (2010), pp. 353–373.
  * [29] T.P. Witelski, Shocks in nonlinear diffusion, Appl. Math. Lett. 8 (1995), pp. 27–32.
  * [30] T.P. Witelski, The structure of internal layers for unstable nonlinear diffusion equations, Stud. Appl. Math. 97 (1996), pp. 277–300.