Spectral stability of shock-fronted travelling waves under viscous relaxation

Regularised PDEs exhibiting travelling wave solutions can frequently be regarded as singularly perturbed dynamical systems, since the regularisation is modelled by the inclusion of strictly higher-order partial derivatives multiplied by a small constant. The singular perturbation in turn has a physical/dynamical interpretation reflecting an inherent hierarchy of spatiotemporal scales, which we refer to more succinctly as a fast-slow structure. We consider the following reaction-nonlinear diffusion PDE with a viscous relaxation term, exhibiting (regularised) shock-fronted travelling waves:¹¹1Throughout the paper we use bar notation to refer to phase space variables. Later on we will use unbarred variables to denote variables defined on the linearized subspaces along the curve, and carets to denote projectivisations of these variables.

for $(x,t)\in\mathbb{R}\times\mathbb{R}$ and $\varepsilon\geq 0$ the singular perturbation parameter characterising the strength of the the regularisation. The nonlinear diffusion and reaction terms given by the quadratic resp. cubic polynomials $D(\bar{U})$ and $R(\bar{U})$. We suppose that $D(\bar{U})<0$ within an interval $(a,b)\subset(0,1)$; the potential function $F(\bar{U}):=\int D(\bar{U})d\bar{U}$ is hence nonmonotone. Nonlinear diffusion processes with nonmonotone potential functions serve as prototypical models for the formation of coherent solutions with sharp fronts; see e.g. [37, 38]. Arguably the canonical system from a physical context is the Cahn-Hilliard model [27]. The system (1) for $\varepsilon=0$ is derived as a particular continuum limit of stochastic agent-based models of invasion processes [31, 32]. The choice of viscous regularisation in (1) appears in the shock regularisation literature [25, 26, 38], but there are other physically-motivated high-order regularising terms that give rise to (smoothed families of) shock-fronted travelling waves. We refer the reader to [22] for a comparison of the existence problems under viscous resp. nonlocal regularising terms relative to system (1).

In the case of existence, the effectiveness of geometric singular perturbation theory (GSPT) is now well established. When (1) is written in a coordinate frame that follows the travelling wave, it takes the structure of a closed fast-slow system of ODEs. A one-parameter family of travelling waves for $0<\varepsilon\ll 1$ is then constructed rigorously as a perturbation of a singular heteroclinic orbit, composed of segments defined by the slow flow along a so-called critical manifold, concatenated with fast jumps along the fast fibres, according to the dynamics of the layer problem. Li et. al. follow this approach in [22], numerically demonstrating the existence of smooth travelling waves of (1) for small values of $\varepsilon>0$. These one-parameter families of waves limit to ‘genuine’ (piecewise continuous) shock-fronted waves as $\varepsilon\to 0$; furthermore, such singular limits are nonunique, strongly depending on the regularisation chosen. Li and coauthors have also studied smooth travelling-wave solutions of the unperturbed problem arising as so-called ‘hole in the wall’ solutions; see [23] and the footnote on page 9. Generally speaking, one can attempt to construct shock-fronted limiting solutions directly as weak solutions, but there may be infinitely many; see eg. [13].

It is worth re-emphasizing the effect of the singular perturbation in (1), in terms of coordinate representations of the travelling wave. With respect to a natural Liénard-type representation, the unperturbed problem can be thought of as modelling the slow dynamics on the critical manifold itself. The introduction of a regularisation term is then tantamount to embedding this slow flow within a higher-dimensional space and introducing a rule (i.e. the fast layer flow) for connecting the wavefront smoothly. In [22], it is explicitly shown that viscous relaxation adds a one-dimensional layer flow, while a fourth-order nonlocal regularization term adds a two-dimensional layer flow.

The objective of this paper is to demonstrate the spectral stability of travelling waves arising in (1). To determine the spectral stability along a wave, we must find the spectrum $\sigma(L)$ of a corresponding linearized operator $L$. The total spectrum is decomposed into its continuous and point components as $\sigma(L)=\sigma_{e}(L)\cup\sigma_{p}(L)$. Typically, most of the work involves deducing the point spectrum $\sigma_{p}(L)$, by recasting the eigenvalue problem $(L-\lambda I)v=0$ as a bifurcation problem posed on the underlying linearized subspaces along the wave. The major point that we want to highlight is that although the existence problem remains amenable to the usual GSPT techniques, the analysis of the stability problem (in particular, the computation of the point spectrum) is quite distinct.

The literature on geometric and analytic techniques for linearized stability problems is substantial (see e.g. [18] and [28] for comprehensive surveys). In this paper we follow the geometric framework developed in a series of seminal papers by Alexander, Gardner, and Jones [1, 8, 14]. The key point is that $L$ also inherits a slow-fast structure, enabling a reduction via linearized subsystems defined along either the slow manifold or along the fast fibres. The typical approach is then to hope that the spectrum of the full system, for $\varepsilon>0$ sufficiently small, is close to the ‘fast’ and/or ‘slow’ spectra of these linearized subsystems, which are defined for $\varepsilon=0$.

Our principal result is Theorem 2.11. As a consequence of the abovementioned high order of the singular perturbation term representing the regularisation, the stability problem for small $\varepsilon>0$ is now close to a slow eigenvalue problem defined on the critical manifold. We give an outline of the proof in Section 6.2. To summarize, we show in Sections 6 and 7 that the eigenvalues for the perturbed problem are the same as the eigenvalues for the reduced problem, which are calculated in Section 8. Together with the fact that the essential spectrum lies in the left half complex plane as shown in Section 5, we can conclude that there is no spectrum of solutions to (1) in any fixed contour $K$ in the right half complex plane.

The principal estimate for constructing these reduced eigenvalue problems is the elephant trunk lemma, which characterises an attracting set over a fast unstable subbundle in the linearized space when $\lambda$ is not a fast eigenvalue; see the construction in the paper of Gardner & Jones [8]. We highlight that a nondegenerate fast eigenvalue problem is a requirement to construct linked elephant trunks over the entire wave. Furthermore, the construction of a slow eigenvalue problem, corresponding to the singular limit of some slow subbundle, is essentially auxiliary to that of the fast problem, and many of the key properties of the slow subbundle are indirectly enforced by the elephant trunk over the fast subbundle.

In order for us to define the promised slow eigenvalue problem, it is essential that we determine how exactly slow linear data on one branch of the slow manifold is transported across fast fibres to another branch of the slow manifold. When an elephant trunk lemma is available, it enforces certain nice properties, such as continuity of the slow subbundle across the fast layer in the singular limit. In other words, slow linear data from one branch of the slow manifolds is transported across fast fibres to the other branch identically, allowing a straightforward definition of slow eigenvalue problems across disjoint subsets of the slow manifolds.

We will show that our problem does not allow for elephant trunk-type estimates over the fast layer. The main obstruction is that, relative to a natural set of coordinates for the linearized system so that setting $\lambda=0$ recovers the usual variational equations along the flow, the parameter $\lambda$ enters the equations ‘weakly,’ i.e. only through $\mathcal{O}(\varepsilon\lambda)$ terms. The singular limit of the fast linearized system therefore degenerates such that there is ‘always a fast eigenvalue,’ i.e. there is a connection made between unstable and stable directions across the wavefront for all $\lambda\in\mathbb{C}$ as $\varepsilon\to 0$. In the present case the eigenvalue problem is low-dimensional enough that we can partly avoid this problem by concentrating on the slow problem only, but we return to this issue in the conclusion when we consider fourth-order regularizations.

Consequently, it must be the case that the slow linearized dynamics is responsible for the generation of eigenvalues near the tails of the wave. In this context, we think of our work as complementing Jones’ early result on the stability of a travelling pulse in the FitzHugh-Nagumo system, in which the eigenvalue problem of the full system is strongly controlled by the fast eigenvalue problems along the wave front and back [14]; as well as the result of Gardner and Jones on the stability of travelling waves in predator-prey systems, where both fast and slow eigenvalue problems are constructed [8].

While we do not have access to a fast elephant trunk over the entire wave, we nevertheless retain control of the linearized dynamics for small values of $\varepsilon$ since $\lambda$ enters the equations only weakly. The key observation here is that the eigenvalue problem can be thought of as a ‘weak’ ($\mathcal{O}(\varepsilon)$)-perturbation of the standard variational problem. We will show that the corresponding transport of slow data is not the identity map, but instead is replaced by a nontrivial jump map that we can write down explicitly in the limit as $\varepsilon\to 0$. The key technical tools that we use are desingularized linearized slow flows; $\varepsilon$-dependent rescalings along the fast fibres; and a continuous differentiability criterion that holds across the fold. We remark that this analysis holds also for the $\lambda=0$ case (corresponding to the regular variational dynamics carried by the wave), but to our knowledge such a statement about jump compatibility is new in the GSPT literature.

The exchange lemma, another fundamental GSPT technique (see e.g. [16, 17]),²²2Strictly speaking, in this paper we consider the case of an inclination lemma (see e.g. [5, 30]). We also point out that the adapation of the standard $(k+\sigma)$-exchange lemma to the case of zero unstable fast directions is straightforward. can also be adapted to our eigenvalue problem, since it is $\mathcal{O}(\varepsilon)$-close to the standard variational problem. Exchange lemma-type estimates allow us to track solutions of the eigenvalue problem as they leave the fast layer and enter small neighborhoods of the slow manifolds, where they enter partial elephant trunks and remain well-controlled along the tails of the wave. The key result here is that solutions entering neighborhoods of the slow manifolds in a generic way will be aligned closely to a slow subbundle defined on the slow manifold after sufficiently long times. This subbundle is defined from the eigenvalue problem, and is therefore not generally equivalent to the tangent bundle of the slow manifold (but it is nearby).

These new techniques cause no extra trouble for the topological arguments that allow us to calculate the point spectrum of the perturbed problem: the corresponding evaluation of first Chern numbers of certain augmented unstable bundles (defined in [1]) follows, as usual. Following the approach of Gardner & Jones in [8], we construct a homotopy between the augmented unstable bundle $\mathcal{E}_{\varepsilon}(K)$ of the ‘full’ problem with $\varepsilon>0$ and that of a reduced problem $\mathcal{E}_{0}(K)$, when $\varepsilon$ is sufficiently small. The reduced vector bundle is defined over separate hemispheres using the slow eigenvalue problem. The jump map now plays an essential role in defining the appropriate clutching function, which glues these vector bundles together along the edges of the hemispheres. We also highlight that the translational eigenvalue at $\lambda=0$ is now counted by the slow problem.

In the course of proving stability in our problem, we also demonstrate explicitly the geometric consequence of eigenvalue crossings in the full system for $0<\varepsilon\ll 1$: the projectivized solutions wind around according to the number of eigenvalues crossed, and these winds occur entirely within the slow dynamics of the critical manifold. The propagation of winding in the reduced problem, as well as the generation of new winds (eigenvalues) can be identified by locating the poles and zeroes, respectively, of a meromorphic Riccati-Evans function [10], defined with respect to a judiciously chosen section of the projective dynamics. This also gives us another opportunity to demonstrate the computational utility of the Riccati-Evans function to count eigenvalues via the argument principle.

We wish to highlight some results in the literature in the context of our present problem. Slow eigenvalues appear in a doubly-diffusive Fitzhugh-Nagumo system which is considered as an early application of the augmented unstable bundle theory in [1], but here the extra fast dynamics is relatively trivial; furthermore, the ‘slow’ problem is in fact fast-slow relative to another small parameter, with its eigenvalues arising from the corresponding fast front and back. The occurrence of such a nontrivial slow eigenvalue problem as described in our current analysis appears to be quite atypical; we highlight the work of Bose [3] in this vein as one of the only comparable examples we have found in the literature. Bose provides similar geometrical results for pulse (meta)stability in singularly-perturbed nonlocal RDEs. The singular perturbation occurs at the second order, so the resulting nonlocal eigenvalue problem can then be studied with traditional Sturm-Liouville techniques, and an oscillation theorem is proven. Relative to Bose’s result, we have two additional technical issues to settle: (i) we must first ‘reduce’ our singularly-perturbed third-order eigenvalue problem (see (24)) by showing closeness to a second-order eigenvalue problem defined on the critical manifold, and (ii) we must also find the corresponding compatibility condition to connect the eigenvalue problem across disjoint branches of the critical manifold.

We also wish to highlight de Rijk et. al.’s recent work on analytical slow-fast factorizations of Evans functions [6]; the relationship between the geometric obstructions discussed above and analytic conditions to produce reduced Evans functions in de Rijk et. al.’s paper bears further exploration, but we do not attempt this in our paper.

We find it sensible to closely follow the general structure in [8] in organizing our paper. In Sec. 2, we describe the construction of a one-parameter family of travelling waves limiting onto a singular shock-fronted travelling wave using geometric singular perturbation theory. In Sec. 3, we define the relevant geometric spaces in which our objects of interest lie and write down some key facts about them. In Sec. 4 we write down the spatial eigenvalue problem for the associated linearized operator of (1), and we write down some leading-order estimates for the eigenvalues and eigenvectors of the associated asymptotic systems. We also define the linearized slow and fast subsystems that we will use often to estimate the linear dynamics of the ‘full’ system. In Sec. 5, we describe some general stability theory for travelling waves and show that the essential spectrum is bounded away from the imaginary axis in the left-half complex plane. In Sec. 6, we construct the augmented unstable vector bundles over a complex contour enclosing the remaining point spectrum, and we also take the opportunity to outline the proof of the main stability theorem in the paper. In Sec. 7, we prove the main uniform estimates for the linearized dynamics in terms of the singular subsystems. We use these estimates to construct a homotopy between some augmented unstable bundles. In Sec. 8, we characterize the point spectrum of the reduced subsystem, finally allowing us to deduce the stability of the wave of the full system for $\varepsilon>0$ sufficiently small. We conclude in Sec. 9, and the most technical lemmas are relegated to the Appendices.

Our first task is to construct constant-speed travelling waves for the PDE (1). In this section we summarise the setup and analysis in [22]. We highlight that the techniques and terminology used here are entirely standard in the context of GSPT; we point to the usual references [15, 21] for general definitions.

Ahead of our analysis in the remaining sections, we explicitly describe the geometric spaces in which our objects of interest live. We will track the dynamics of fibers of complex vector bundles, corresponding to points in $\mathbf{Gr}(k,n)$, the Grassmannian of complex $k$-planes in $\mathbb{C}^{n}$. The general set-up is standard: we will study dynamics on the Grassmannian indirectly via embedded systems on complex projective space.

We first introduce some essential notation for complex projective space. The space $\mathbb{CP}^{n}$ is the quotient space $(\mathbb{C}^{n+1}-\{0\})/\sim{}$ subject to the equivalence relation that identifies complex rays, i.e. $y_{1},y_{2}\in\mathbb{C}^{n+1}-\{0\}$ are identified if $y_{1}=\alpha y_{2}$ for some $\alpha\in\mathbb{C}$. We denote the projection map by $\pi:\mathbb{C}^{n+1}-\{0\}\to\mathbb{CP}^{n}$, and for each $y\in\mathbb{C}^{n}$ we use the notation $[y]$ or $\hat{y}$ for the image $\pi(y)\in\mathbb{CP}^{n}$ (this notation also holds for sets $S\subset\mathbb{C}^{n}-\{0\}$).

Later, in Section 7.2 we will use the following elementary lemma to interpret metric neighborhoods on $\mathbb{CP}^{2}$ in terms of those on $\mathbb{C}^{3}$.

Proof of Lemma 3.1: Suppose that there is no such uniform constant for some fixed choice of $\delta>0$. Then we can find a sequence $\{y_{k}\}$ lying on the unit sphere in $\mathbb{C}^{n+1}$ such that for each $k$, we have $\hat{y}_{k}\notin N_{\delta}(\hat{P})$ and $|w_{k}|<\delta/k$, where $w_{k}$ denotes the orthogonal projection of $y_{k}$ with respect to $P$. Choose a convergent subsequence $y_{k}\to y_{*}$, relabeling $k$ if necessary. Then $y_{*}$ clearly has zero component in the direction orthogonal to $P$, implying that $y_{*}\in P$ and hence that $\hat{y}_{*}\in N_{\delta}(\hat{P})$. On the other hand, we have $\hat{y}_{k}\to\hat{y}_{*}\notin N_{\delta}(\hat{P})$ since the quotient map is continuous and the complement of $N_{\delta}(\hat{P})$ is closed, giving a contradiction. $\Box$

Let $v_{1},\cdots,v_{k}\in\mathbb{C}^{n}$ denote a spanning set of vectors for the $k$-plane $W$. The well-known Plücker embedding

provides homogeneous coordinates for the Grassmannian, called Plücker coordinates. The image $\psi(W)$ can be coordinatized by writing down the maximal minors/cofactors of the $(n\times k)$ matrix $M_{W}=[v_{1}~{}\cdots~{}v_{k}]$.

Example. Consider a complex 2-plane $W\in\mathbf{Gr}(2,3)$ with spanning complex vectors $v_{1}=(a_{1},b_{1},c_{1})^{T}$ and $v_{2}=(a_{2},b_{2},c_{2})^{T}$. A coordinate representation of $[v_{1}\wedge v_{2}]\in\mathbb{C}\mathbb{P}^{2}$ is given by the list of $2\times 2$ cofactors of the $3\times 2$ matrix $[v_{1}~{}~{}v_{2}]$, namely

Observe that the Plücker coordinate representation $N_{W}$ can be interpreted as a parametrisation of the plane $W$ by a (nontrivial) normal vector.

The Grassmannian spaces are metrizable. Here we construct an explicit metric, because we will need a way to measure ‘angles’ between lines and planes. Using a complex version of the construction given in [17], we take advantage of the Plücker embedding to equip $\mathbf{Gr}(k,n)$ with the Fubini-Study metric. Complex projective space can be embedded into the space of Hermitian $(n+1)\times(n+1)$ matrices $\mathcal{H}$ via the following isometric immersion:

where $x\in\mathbb{C}^{n+1}-\{0\}$ represents the ray $\mathbb{C}x$. The space of Hermitian matrices is in turn endowed with an inner product given by

which allows us to define a norm $\lVert P\rVert=d(P,P)^{1/2}$, and hence a metric induced by the norm given by $\text{dist}(P,Q):=\lVert P-Q\rVert$. We can pull back this metric on $\mathcal{H}$ to one on $\mathbb{CP}^{n}$, and hence on the Grassmannian:

We observe that $\text{dist}([x],[y])$ recovers the sine of the so-called Hermitian angle [29] between a corresponding pair of nontrivial complex vectors representing $[x],\,[y]$. For example, this metric on $\mathbf{Gr}(2,3)$ measures the Hermitian angle between a pair of normal vectors parametrising a corresponding pair of complex $2$-planes.

We will use the Plücker embedding to study (linear) dynamics which describe the evolution of $k$-planes in the Grassmannian. Concretely, let $y^{\prime}=b(\xi)y$ denote a nonautonomous linear system on $\mathbb{C}^{n+1}$. This linear system induces a nonlinear flow on $\mathbb{CP}^{n}$, and hence on the corresponding Grassmannian, which we denote by

We mention two elementary facts about about these induced flows. The first is that $y$ is an eigenvector of $b(\xi)$ if and only if $\hat{y}$ is a critical point of $\hat{b}(\hat{y},\xi)$. The second is that if $\{\lambda_{i}\}$ denotes the eigenvalues of the linearization of a (constant) matrix $b$ and if $y_{i}$ is an eigenvector corresponding to the eigenvalue $\lambda_{i}$, then the eigenvalues of the linearization of $\hat{b}$ near the corresponding critical point $\hat{y}_{i}$ are $\lambda_{j}-\lambda_{i}$ for $j\neq i$. As a consequence of these facts, if the largest eigenvalue of a linear system is simple, then the corresponding eigendirection in the projectivized system is a stable fixed point.

Proof: Choose an ordered basis for $V$ consisting of eigenvectors of $A$, and use it to construct an ordered basis for $\bigwedge^{2}V$. The conclusion follows from a straightforward calculation using the definition (22). $\Box$

If $x^{\prime}=Ax$ is a linear vector field, then $P^{\prime}=A^{(\sigma)}P$ defines a linear vector field on the $\sigma$th exterior power space. By using the Plücker embedding, the latter induced vector field can be interpreted as evolving points in $\mathbf{Gr}(\sigma,n)$ (see Sec. 3.2 in [17] for a complete description of this construction in the real setting; the extension to complex $k$-planes is direct).

We now assume that Hypotheses 2.7, 2.8, and 2.9 hold, and we proceed to define the spatial eigenvalue problem associated with the one-parameter family of waves $(\Gamma_{\varepsilon},c({\varepsilon}))$. Let us fix $\varepsilon>0$ sufficiently small and represent $\Gamma_{\varepsilon}$ by the solution $\bar{X}(\zeta,\varepsilon)=(\bar{U}(\zeta,\varepsilon),\bar{P}(\zeta,\varepsilon),\bar{V}(\zeta,\varepsilon))$. We consider perturbations of the form