Pulse solutions for an extended Klausmeier model with spatially varying coefficients

TL;DR

This paper analyzes an extended Klausmeier ecological model with spatially varying coefficients, establishing the existence and stability of pulse solutions, and revealing how spatial heterogeneity influences their stability and bifurcations.

Contribution

It provides a rigorous mathematical framework for existence and spectral stability of pulse solutions in a spatially heterogeneous Klausmeier model, including bifurcation analysis.

Findings

Spatial heterogeneity can stabilize or destabilize pulse solutions.

Existence of stationary multi-pulse solutions due to bifurcations.

Identification of a pitchfork bifurcation affecting pulse stability.

Abstract

Motivated by its application in ecology, we consider an extended Klausmeier model, a singularly perturbed reaction-advection-diffusion equation with spatially varying coefficients. We rigorously establish existence of stationary pulse solutions by blending techniques from geometric singular perturbation theory with bounds derived from the theory of exponential dichotomies. Moreover, the spectral stability of these solutions is determined, using similar methods. It is found that, due to the break-down of translation invariance, the presence of spatially varying terms can stabilize or destabilize a pulse solution. In particular, this leads to the discovery of a pitchfork bifurcation and existence of stationary multi-pulse solutions.

Equations

$$\left\{\begin{array}[]{rcrcl}\partial_{t}U&=&\partial_{x}^{2}U&+&\mathcal{H}_{1}\left(x,u,u_{x},v,v_{x};\tilde{\varepsilon}\right),\\ \partial_{t}V&=&\tilde{\varepsilon}^{2}\partial_{x}^{2}V&+&\mathcal{H}_{2}\left(x,u,u_{x},v,v_{x};\tilde{\varepsilon}\right),\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{rcrl}\partial_{t}U&=&\partial_{x}^{2}U&+f(x)\partial_{x}U+g(x)U+a-U-UV^{2}\,,\\[5.69046pt] \partial_{t}V&=&D^{2}\partial_{x}^{2}V&-\ mV+UV^{2}\,,\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{rcrl}0&=&\tilde{u}_{xx}&+f(x)\tilde{u}_{x}+g(x)\tilde{u}+a-\tilde{u}-\tilde{u}\tilde{v}^{2}\,,\\[2.84544pt] 0&=&\frac{D^{2}}{m}\tilde{v}_{xx}&-\tilde{v}+\frac{1}{m}\tilde{u}\tilde{v}^{2}\,.\end{array}\right.$$

$$\left\{\begin{array}[]{rcl}\partial_{x}\hat{u}&=&\hat{p},\\ \partial_{x}\hat{p}&=&-f(x)\hat{p}-g(x)\hat{u}-1+\hat{u}.\end{array}\right.$$

$$(A1):$$

$$(A2):$$

$$(A3):$$

$$(A4):$$

$$(A5):$$

$$\xi :=\frac{\sqrt{m}}{D}x,\tilde{u}=\frac{m\sqrt{m}D}{a}u,\tilde{v}=\frac{a}{\sqrt{m}D}v,$$

$$\displaystyle\left\{\begin{array}[]{rcl}u_{\xi\xi}&=&\frac{a^{2}}{m^{2}}\left[\frac{D^{2}m}{a^{2}}u-\frac{Dm\sqrt{m}}{a^{2}}f\left(\frac{D}{\sqrt{m}}\xi\right)u_{\xi}-\frac{D^{2}m}{a^{2}}g\left(\frac{D}{\sqrt{m}}\xi\right)u-\frac{D}{\sqrt{m}}+uv^{2}\right]\,,\\[2.84544pt] v_{\xi\xi}&=&v-uv^{2}\,.\end{array}\right.$$

$$0<\epsilon :=\frac{a}{m},0<\mu :=\frac{m\sqrt{m}D}{a^2},$$

$$\displaystyle\left\{\begin{array}[]{rcl}\dot{u}&=&\varepsilon p\,,\\[2.84544pt] \dot{p}&=&\varepsilon\left[\varepsilon^{2}\mu^{2}u-\varepsilon\mu f\left(\varepsilon^{2}\mu\xi\right)p-\varepsilon^{2}\mu^{2}g\left(\varepsilon^{2}\mu\xi\right)u-\varepsilon^{2}\mu+uv^{2}\right]\,,\\[2.84544pt] \dot{v}&=&q\,,\\[2.84544pt] \dot{q}&=&v-uv^{2}\,.\end{array}\right.$$

$$0<\epsilon \ll 1.$$

$$\underset{\xi \leftarrow \pm \infty }{\lim }(u_b,p_b,0,0)=\left(1/\mu ,0,0,0\right).$$

$$\displaystyle\left\{\begin{array}[]{rcl}\dot{u}&=&\varepsilon p\,,\\[2.84544pt] \dot{p}&=&\varepsilon\left[\varepsilon^{2}\mu^{2}u-\varepsilon^{2}\mu+uv^{2}\right]\,,\\[2.84544pt] \dot{v}&=&q\,,\\[2.84544pt] \dot{q}&=&v-uv^{2}\,.\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{rcl}u^{\prime}&=&p\,,\\[2.84544pt] p^{\prime}&=&\varepsilon^{2}\mu^{2}u-\varepsilon^{2}\mu+uv^{2}\,,\\[2.84544pt] \varepsilon v^{\prime}&=&q\,,\\[2.84544pt] \varepsilon q^{\prime}&=&v-uv^{2}\,.\end{array}\right.$$

$$\tilde{\mathcal{M}}:=\{(u,p,0,0)|u>0\}$$

$$\displaystyle\left\{\begin{array}[]{rcl}u^{\prime}&=&p\,,\\[2.84544pt] p^{\prime}&=&\varepsilon^{2}\mu^{2}u-\varepsilon^{2}\mu\,,\end{array}\right.$$

$${\tilde{l}}^{u/s}:=\left\{(u,p)|p=\epsilon \mu (u-\frac{1}{\mu })\right\}.$$

$$\displaystyle\left\{\begin{array}[]{rcl}\dot{u}&=&\ 0\,,\quad\dot{p}\ =\ 0\,,\\[2.84544pt] \dot{v}&=&q\,,\\[2.84544pt] \dot{q}&=&v-uv^{2}\,.\end{array}\right.$$

$$H(v,q;u)=\frac{1}{2}{q}^2-\frac{1}{2}{v}^2+\frac{1}{3}u{v}^3.$$

$$\displaystyle\left\{\begin{array}[]{rcl}v_{hom}^{(0)}(\xi;u_{0})&=&\frac{1}{u_{0}}\,\omega(\xi)\,,\quad\omega(\xi):=\frac{3}{2}\,\mathrm{sech}^{2}(\xi/2)\,,\\ q_{hom}^{(0)}(\xi;u_{0})&=&\dot{v}_{hom}(\xi;u_{0})\,,\quad u_{0}\in\mathbb{R}\backslash\{0\}\,,\end{array}\right.$$

$$\tilde{K}(u,p,v,q)=H(v,q;u).$$

$$I_f:=\left(-\frac{1}{\sqrt{\epsilon }},\frac{1}{\sqrt{\epsilon }}\right),$$

$${\Delta }_{I_f}\tilde{K}=\tilde{K}(1/\sqrt{\epsilon })-\tilde{K}(-1/\sqrt{\epsilon })={\int }_{I_f}\frac{d}{d\xi }\tilde{K}(\xi )d\xi =\frac{1}{3}\epsilon {\int }_{I_f}p(\xi )v_{hom}(\xi {)}^{3}d\xi +h.o.t.$$

$${\Delta }_{I_f}u$$

$${\Delta }_{I_f}p$$

$${\tilde{T}}_{o/d}:=\left\{\left(u,p,0,0\right)|p=\mp \frac{3\epsilon }{u},u>0\right\},$$

$$p(\pm 1/\sqrt{\epsilon })=p(0)\pm \frac{1}{2}{\Delta }_{I_f}p=\epsilon \left(p^{(1)}(0)\pm \frac{3}{u_0}\right).$$

$$u_0^{\pm }=\frac{1\pm \sqrt{1-12\mu }}{2\mu }\left(\text{mbox}for\mu \le \frac{1}{12}\right).$$

$$\begin{array}[]{rcrcrcrcrc}u_{0}^{-}&=&&&3&+&9\mu&+&\mathcal{O}(\mu^{2})\\ u_{0}^{+}&=&\frac{1}{\mu}&-&3&-&9\mu&+&\mathcal{O}(\mu^{2})\end{array}$$

$$s(\xi ):=\frac{D}{\sqrt{m}}\xi ={\epsilon }^2\mu \xi ,$$

Full text

1 Introduction

Since Alan Turing’s revolutionary insight that patterns can emerge spontaneously in systems with multiple species if these diffuse at different rates [43], systems of reaction-diffusion equations have served as prototypical pattern forming models. Scientists have been using these reaction-diffusion models successfully to describe for instance animal markings [28], embryo development [32] and the faceted eye of Drosophila [31]. Special interest has been given to localized solutions (e.g. pulses, fronts), that arise when the diffusivity of species involved is very different. The prototypical (two-component) model (in one spatial dimensional) is a singularly perturbed equation of the (scaled) form

[TABLE]

where $0<\tilde{\varepsilon}\ll 1$ is a measure for the ratio of diffusion constants, and $\mathcal{H}_{1}$, $\mathcal{H}_{2}$ are sufficiently smooth functions. Because of the singular perturbed nature of (1), it is possible to establish existence and determine (linear) stability of localized patterns in these models. In the past, this has been done successfully for the Gray-Scott model [14, 15, 17, 10, 29, 41], the Gierer-Meinhardt model [16, 45, 17, 41], and in several other settings [12, 22, 36, 33]. However, these studies are usually limited to models with constant coefficients. Some research has focused on the introduction of localized spatial inhomogeneities [44, 34, 35, 48, 49, 21]; also (often formal) research has been done on reaction-diffusion equations with (less restricted) spatially varying coefficients [9, 8, 2, 47, 46, 7]. In this article, we aim to expand the knowledge of such systems, by studying a reaction-diffusion system with fairly generic spatially varying coefficients rigorously; motivated by its use in ecology (see Remark 2), we consider the following extended Klausmeier model with spatially varying coefficients [27, 4]:

[TABLE]

with $x\in\mathbb{R},t\geq 0,U=U(x,t),V=V(x,t)\in\mathbb{R}$, parameters $D,a,m>0$ and functions $f,g\in C^{1}_{b}(\mathbb{R})$. Certain conditions are imposed on the parameters and functions $f$ and $g$ – these will be explained in section 1.1.

Remark 1**.**

The model (4) can be brought into the form of (1) by a series of scalings – see section 2 and [17].

Remark 2** (Application of the extended Klausmeier model).**

This system of equations is used as a model in ecology to describe the dynamics of vegetation ($U$) and water ($V$). The extended Klausmeier model (4) takes into account the amount of rainfall ($a>0$) and mortality rate of the vegetation ($m>0$) and goes beyond its classical version by modeling a smooth, spatially varying terrain $h=h(x)$ which then enters (4) as $f(x)=h^{\prime}(x),g(x)=h^{\prime\prime}(x)$ (see [4]). Variants of the Klausmeier model have been studied in various articles ranging from ecological studies [27, 5] to mathematical analysis [4, 40, 38, 39]. The focus of all these studies are vegetation patterns, which have been found to play a crucial role in the process of desertification. A starting point for the analysis of more complicated patterns is a thorough understanding of their building blocks, namely, localized solutions. The present paper is motivated by observations – both in numerical simulations and in real ecosystems [4, 5] – of the impact of nontrivial topographies on the dynamics of localized vegetation patterns.

The focus of this article is to analyze existence, stability and (some) bifurcations of stationary pulse solutions to (4). The presence of spatially varying coefficients, however, alters the approach that usually is taken in the case of constant coefficients models. For one, with spatially constant coefficients, (4) possesses a uniform stationary state, with $V\equiv 0$, to which pulse solutions converge for $x\rightarrow\pm\infty$. In the case of spatially varying coefficients, however, typically such uniform stationary state does not exist; instead, a bounded solution $(u,v)=(u_{b},0)$ exists and pulse solutions converge to this bounded solution for $x\rightarrow\pm\infty$ – see Figure 1. Moreover, standard proofs using geometric singular perturbation theory typically rely on the availability of closed form expressions for orbits of subsystems of (4) – see below. These are no longer available in case of generic spatially varying coefficients, and only bounds can be found. Indeed, the core contribution of the present work is to overcome these difficulties, which we do by blending geometric singular perturbation theory [24] with the theory of exponential dichotomies [11] in a new way.

In this article, we initially follow the ‘standard’ approach of geometric singular perturbation theory. That is, we introduce a small parameter $\varepsilon:=\frac{a}{m}$ – see assumption (A1) in section 1.1– and construct a stationary pulse solution to (4) in the limit $\varepsilon=0$, which present itself as a homoclinic orbit in the related stationary fast-slow ODE system – in case of spatially varying coefficients it is homoclinic to the bounded solution. For this construction, the full system is split into a fast subsystem, and a (super)slow subsystem on a so-called slow manifold $\mathcal{M}$ that consists of fixed points of the fast subsystem. We establish fast connections to and from $\mathcal{M}$ that take off from submanifold $T_{o}\subset\mathcal{M}$ and touch down on submanifold $T_{d}\subset\mathcal{M}$. On $\mathcal{M}$, we construct stable and unstable submanifolds $W^{s/u}(u_{b})\subset\mathcal{M}$ that consists of points on $\mathcal{M}$ that converge to the bounded solution for $x\rightarrow\infty$ respectively $x\rightarrow-\infty$. Intersections between these unstable/stable manifolds and take-off/touch-down submanifolds (and a symmetry assumption) then establish the existence of pulse solutions to (4). Finally, persistence of these pulse solutions for $\varepsilon>0$ is guaranteed by geometric singular perturbation theory [24].

Specifically, stationary solutions $(U(x,t),V(x,t))=(\tilde{u}(x),\tilde{v}(x))$ of (4) fulfill the system of ODEs

[TABLE]

After a sequence of (re)scalings, it can be seen that the associated fast subsystem is not affected by the spatially varying terms and can be studied using standard methods. However, the slow subsystem, on the slow manifold $\mathcal{M}$, is affected by the spatially varying terms. This subsystem is given (when rescaling $\hat{u}=a\tilde{u}$) by

[TABLE]

For $f$ and $g$ constant, (8) can be solved explicitly and the stable and unstable manifolds $W^{s,u}(u_{b})$ are known explicitly. In case of (spatially) varying $f$ and $g$, typically no closed form solutions are available; however, when these varying coefficients are sufficiently small – specifically, when $\delta:=\sup_{x\in\mathbb{R}}\sqrt{f(x)^{2}+g(x)^{2}}<\frac{1}{4}$ (so $\delta$ can be $\mathcal{O}(1)$ with respect to $\varepsilon$); see section 2.3 – the dynamics of (8) can be related to the constant coefficient case $f,g\equiv 0$ using the theory of exponential dichotomies.

In particular, the saddle structure – present for $f,g\equiv 0$ – persists as exponential dichotomy. Therefore, (8) possesses a $1D$ family of solutions that converge to the (unique) bounded solution to (8) for $x\rightarrow\infty$ and a $1D$ family of solutions that converge to the bounded solution for $x\rightarrow-\infty$. These families of solutions essentially form the stable and unstable manifolds $W^{s,u}(u_{b})$. Due to the linear nature of (8), these (un)stable manifolds are made up of straight lines, i.e. $W^{s,u}(u_{b})=\cup_{x\in\mathbb{R}}(x,l^{s,u}(x))$ where $l^{s,u}(x)$ describes a straight line in $\mathbb{R}^{2}$. An important difference now arises between the cases of constant and varying coefficients: when $f,g\equiv 0$, the lines $l^{s,u}(x)$ do not depend on $x$; when $f$ and $g$ are spatially varying, they do. Hence, $W^{s,u}(u_{b})$ appears wiggly in case of varying coefficients – see Figure 2. The theory of exponential dichotomies enables us to bound the variation of the lines $l^{s,u}(x)$; if $\delta$ is small enough (i.e. $\delta<\delta_{c}(a,m,D)$, where $\delta_{c}\leq 1/4$ is $\mathcal{O}(1)$ with respect to $\varepsilon$), these bounds are strict enough that a non-empty intersection $(0,l^{u}(0))\cap T_{o}$ is guaranteed – thus establishing existence of a (symmetric) pulse solution to (4). See Figure 3 for a sketch.

Next, the spectral stability of the thus created pulse solutions is studied. Using similar bounds as in the existence problem, it is shown that eigenvalues are $\delta$-close to their counterparts in case of constant coefficients – see Figure 4. That is, under several conditions, typical for these systems, the ‘large’ eigenvalues can be bounded to the stable half-plane $\{\lambda\in\mathbb{C}:\mbox{Re}\lambda<0\}$. For the ‘small’ eigenvalue – located close to the origin – it is more subtle. In case of $f,g\equiv 0$ this small eigenvalue is located precisely at the origin due to the translation invariance of (4). The introduction of spatially varying coefficients to the system breaks this invariance and as a result the small eigenvalue moves to the stable or the unstable half-plane.

Tracking of this eigenvalue indicates that it can, indeed, move to either half-plane, depending on the form of the functions $f$ and $g$. In particular, when taking $f=h^{\prime}$, $g=h^{\prime\prime}$, the location of the small eigenvalue is related to the curvature $g=h^{\prime\prime}$ of $h$: when the curvature is weak, the pulse solution is stable if $g(0)=h^{\prime\prime}(0)<0$ and unstable if $g(0)=h^{\prime\prime}(0)>0$; for strong curvature, this is flipped, due to a pitchfork bifurcation.

Finally, the break-down of the translation invariance in (4) has another novel effect. In case of constant coefficients, stationary multi-pulse solutions – solutions with multiple fast excursions – do not exist, due to the presence of the translation invariance. If this invariance is broken, they can exist; the introduction of functions $f$ and $g$ now allows for these stationary multi-pulse solutions (under some conditions on $f$ and $g$) and their existence can be established (although we refrain from going in the details).

The set-up for the rest of this paper is as follows. In section 2, we establish existence of stationary pulse solutions to (4); here we first consider the case $f,g\equiv 0$ and subsequently the case of generic (bounded) $f$ and $g$. Then, using the theory of exponential dichotomies, both cases are related to each other, resulting in bounds for the generic case that allow us to prove existence. In section 3 we study the spectral stability of found pulse solutions, again by relating the generic case to the constant coefficient case of $f,g\equiv 0$. Then, in section 4 we consider the small eigenvalues more in-depth using formal and numerical techniques, focusing on the possible occurrence of bifurcations; we also present stationary multi-pulse solutions. We conclude with a discussion of the results in section 5.

1.1 Assumption

We will make several assumptions throughout the manuscript. Some are crucial, while some serve to simplify the exposition.

[TABLE]

Assumption (A1) ensures the presence of a small parameter, necessary to use geometric singular perturbation theory [37, 4]. (A2) is a symmetry assumption, that ensures (4) possesses a (point) symmetry in $x=0$; this technicality significantly simplifies our rigorous proof; pulse solutions can also be found formally and/or numerically when (A2) does not hold (and we expect that their existence can be established rigorously by extending our methods). Then, assumption (A3) stems from the theory of exponential dichotomies: when this holds, solutions to (8) for generic $f$ and $g$ can be linked to solutions of (8) with $f,g\equiv 0$; when (A3) does not hold, this link is not provided by the theory of exponential dichotomies. Assumption (A4) is a technicality that is only needed in the stability section (specifically for the elephant-trunk method to work); for the existence theorems it is not necessary; in fact, it is suspected that even stability results continue to hold when (A4) is violated – see also Remarks 20 and 21. Finally, assumption (A5) is needed to pass limits in the treatment of the fast-slow system.

2 Analysis of stationary pulse solutions

A crucial step for making the stationary ODE (7) amenable to analytic considerations is to find a parameter regime convenient for rigorous perturbation techniques. While there are various choices, we pick a specific one for clarity, since our focus is on novel phenomena due to the non-autonomous character of the system and not to classify all possible dynamics across parameter regimes.

Following [14, 10, 4], we rescale the spatial coordinate (motivated by the diffusivity of the $v$-component) and the amplitudes of the unknowns by

[TABLE]

to get

[TABLE]

It is now convenient to introduce

[TABLE]

and write the above ODEs as the first order system of ODEs

[TABLE]

In order to use geometric singular perturbation theory, we make the customary assumption (A1), that is,

[TABLE]

and stipulate assumption (A5) so we can pass to limits.

In the autonomous case $f\equiv 0$ and $g\equiv 0$, system (23) has a fixed point $\left(1/\mu,0,0,0\right)$ and stationary pulse solutions of (4) correspond to orbits that are homoclinic to $\left(1/\mu,0,0,0\right)$; see Figure 1(a) for an example. In the non-autonomous case $f\neq 0,g\neq 0$ there is no fixed point, but instead a unique bounded solution $(u_{b},p_{b},0,0)$. In this case, stationary pulse solutions of (4) correspond to orbits that are homoclinic to this bounded solutions; see Figures 1(b) and 1(c) for examples. The existence of said unique bounded solution $(u_{b},p_{b},0,0)$ is established in the following proposition proven later in section 2.3 (in the proof of Proposition 4).

Proposition 1** (Existence of a bounded solution for (23)).**

Let assumptions (A3) and (A4) be fulfilled. Then (23) has a unique bounded solution $(u_{b},p_{b},0,0)$ that satisfies

[TABLE]

Remark 3** (Orbits homoclinic to bounded solutions).**

Note that the assumption $\lim_{x\rightarrow\pm\infty}f(x),g(x)=0$ in (A4) is not necessary for the existence proof, but will be used in the stability analysis. In case $f,g$ are only bounded without approaching a constant state when $|x|\rightarrow\infty$, the corresponding constructed pulse solution is also a homoclinic to the respective bounded solution. An illustration of such a case is given in Figure 1(c), where, due to the periodicity of the coefficients $f,g$, the bounded background solution is periodic and so is the pulse solution in its tails.

To highlight the novelty of the presented approach, we first briefly explain how the construction is carried out in the constant coefficient case $f=g=0$, to then proceed to the non-autonomous case.

2.1 Stationary pulse solutions for $f=0$ and $g=0$

The fast system reads

[TABLE]

Note that this system possesses the symmetry $(\xi,u,p,v,q)\rightarrow(-\xi,u,-p,v,-q)$. The corresponding slow system in the slow scaling $\eta=\varepsilon\xi$ is given by

[TABLE]

Restricted to the invariant manifold

[TABLE]

it reads

[TABLE]

which has a saddle structure around the fixed point $\left(\frac{1}{\mu},0\right)$ with stable and unstable eigenspaces given by

[TABLE]

Remark 4**.**

Note that this step is much more intricate in the case of varying coefficients $f,g$ where explicit solutions are possible only for very specific choices of coefficients. Therefore, one must resort to estimation techniques for the general case. Overcoming this difficulty using exponential dichotomies is the core contribution of the present work.

The reduced fast system has the form

[TABLE]

A sketch of its planar subsystem $\dot{v}=q,\dot{q}=v-uv^{2}$ can be found in Figures 5(a); this planar subsystem is a Hamiltonian system with Hamiltonian

[TABLE]

Its fixed point $(v,q)=(0,0)$ features a saddle structure and a family of homoclinic orbits

[TABLE]

connecting its stable and unstable manifolds. Hence, (44) is a Hamiltonian system with Hamiltonian

[TABLE]

The invariant manifold $\widetilde{\mathcal{M}}$ from (36) is the collection of saddle points $(u,p,0,0),u>0,p\in\mathbb{R},$ for (44) and is, hence, normally hyperbolic. For its stable and unstable manifolds $W_{0}^{s/u}(\widetilde{\mathcal{M}})$ it holds true that $\dim[W_{0}^{s/u}(\widetilde{\mathcal{M}})]=3$ and, in fact, $W_{0}^{s}(\widetilde{\mathcal{M}})$ and $W_{0}^{u}(\widetilde{M})$ (partly) coincide, where the intersection is simply given by the family of homoclinic orbits. Moreover, we have that $\widetilde{K}(u,p,v,q)|_{(u,p,v,q)\in\widetilde{\mathcal{M}}}=0$.

For $\varepsilon>0$, we note that $\widetilde{\mathcal{M}}$ is still an invariant manifold of the full system (30). It is a standard result in geometric singular perturbation theory (see, e.g. the classic articles [42, 24, 26] or, more recent, [30]) that, for $\varepsilon$ sufficiently small, its stable and unstable manifolds persist as $W_{\varepsilon}^{s/u}(\widetilde{\mathcal{M}})$ with $\dim[W_{\varepsilon}^{s/u}(\widetilde{\mathcal{M}})]=3$, but do not necessarily coincide anymore. In fact, they generically meet in a 2D intersection in $\mathbb{R}^{4}$.

In order to analyze the persistence of homoclinic orbits we measure the distance of $W_{\varepsilon}^{s}(\widetilde{\mathcal{M}})$ and $W_{\varepsilon}^{u}(\widetilde{\mathcal{M}})$ in the hyperplane $\widetilde{R}=\{(u,p,v,q)~{}|~{}q=0\}$, that is, we fix an even homoclinic orbit $(u_{hom},p_{hom},v_{hom},q_{hom})$ with $(u_{hom}(0),p_{hom}(0),v_{hom}(0),q_{hom}(0))=(u_{0},p_{0},v_{max},0)$. To this end we use the Hamiltonian $\widetilde{K}$ and analyze its difference during the jump of the orbit through the fast field

[TABLE]

by setting up

[TABLE]

where we used that $\frac{d}{d\xi}\widetilde{K}=\frac{\partial}{\partial u}H(v,q;u)(\frac{du}{d\xi})+\frac{d}{d\xi}H(v,q;u)=\frac{1}{3}v^{3}(\frac{du}{d\xi})+0=\frac{1}{3}\varepsilon v^{3}p$. We may set (using the fact that $p$ is constant to leading order) $p(\xi)=p^{(0)}+\varepsilon p^{(1)}(\xi)+h.o.t.$ Therefore, in order to make this difference vanish to leading order, we evidently need that $p^{(0)}=0$ and $p^{(1)}(0)=0$.

Now that a departure and return mechanism from and back to $\widetilde{\mathcal{M}}$ is established through the intersection $W_{\varepsilon}^{s}(\widetilde{\mathcal{M}})\cap W_{\varepsilon}^{u}(\widetilde{\mathcal{M}})\cap R$, the remaining task is to determine possible take-off and touch-down points on $\widetilde{\mathcal{M}}$ and investigate if these intersect the stable and unstable eigenspaces $l^{s/u}$ appropriately to form a homoclinic. To this end we observe that

[TABLE]

so, to leading order, only the $p$-variable changes during the fast jump, and therefore, the take-off and touch-down curves on $\widetilde{\mathcal{M}}$ are to leading order given by

[TABLE]

where we used that, by symmetry, to leading order

[TABLE]

Finally, a straightforward computation of the intersection points of these with the stable and unstable eigenspaces $l^{s/u}$ gives two possible homoclinics when $\mu\leq\frac{1}{12}$, with

[TABLE]

Remark 5**.**

When $\mu\ll 1$, the expression for $u_{0}^{\pm}$, (56), can be expanded in terms of $\mu$; this yields for $u_{0}^{\pm}$ the following expansions

[TABLE]

A conceptual sketch of the dynamics on $\widetilde{\mathcal{M}}$, along with an excursion through the fast field, is given in Figure 5(b). Moreover, in Figures 6(a) and 6(b), the evolution of a homoclinic solution is projected onto manifold $\widetilde{\mathcal{M}}$.

2.2 Stationary pulse solutions for varying $f$ and $g$

First, we convert the non-autonomous system into an autonomous one by setting

[TABLE]

which gives the extended (autonomous) fast system

[TABLE]

It is important to note that the symmetry assumptions (A2) on $f$ and $g$ translate directly into a symmetry for (64) which is crucial for the construction of a homoclinic.

Lemma 1** (Symmetry of (64)).**

Let the symmetry assumptions (A2) be fulfilled, that is, let $f$ be an odd function and $g$ be an even function. Then we have for (64) the symmetry $(s,u,p,v,q)\rightarrow(-s,u,-p,v,-q)$.

The slow system corresponding to (64) in the slow variable $\eta=\varepsilon\xi$ is given by

[TABLE]

It possesses a three-dimensional invariant manifold

[TABLE]

on which it takes the form

[TABLE]

which is an extension of the non-autonomous system

[TABLE]

It is now convenient to introduce (or, actually, return to) the super-slow variable $x=\varepsilon\mu\eta$. We set $u(\eta)=\frac{1}{\mu}\hat{u}(\varepsilon\mu\eta)=\frac{1}{\mu}\hat{u}(x)$ and return to the second order non-autonomous setting

[TABLE]

Lemma 2** (Symmetry of (81)).**

Let the symmetry assumptions (A2) be fulfilled, that is, let $f$ be an odd function and $g$ be an even function. Then we have for (81) the symmetry $(x,\hat{u},\hat{p})\rightarrow(-x,\hat{u},-\hat{p})$.

Remark 6**.**

*For conciseness, we note that we have three different scales:

**fast scale $\xi\,,\quad$ slow scale $\eta=\varepsilon\xi\,,\quad$ super-slow scale $x=\varepsilon\mu\eta=\varepsilon^{2}\mu\xi$

**

The construction that we illustrate in this article therefore relies heavily on assumption (A1). The specific definition of the small parameter is convenient since the fast reduced system is an ODE which is known to have homoclinic solutions and the slow system on the critical manifold $\mathcal{M}$ is a linear planar system.

Remark 7**.**

Note the difference between $p=\frac{du}{d\eta}$ and $\hat{p}=\frac{d\hat{u}}{dx}$. Hence, $p=\varepsilon\hat{p}$.

Proposition 2** (Dynamics on $\mathcal{M}$).**

Consider the slow system on $\mathcal{M}$ (75) with $f,g$ fulfilling (A3). Then there exists a unique bounded solution $(\hat{u}_{b},\hat{p}_{b})$ of (81) and corresponding connected set $\Gamma\subset\mathbb{R}\cup\{\infty\}$ such that the following holds true: For each fixed $x\in\mathbb{R}$ there exists $C^{s/u}(x)\in\Gamma$ and lines

[TABLE]

such that the solution to the initial value problem (81) with $(\hat{u}(x),\hat{p}(x))=(\hat{u}_{0},\hat{p}_{0})\in l^{s}(x)$ converges to $(\hat{u}_{b},\hat{p}_{b})$ for $x\rightarrow\infty$, while with $(\hat{u}(x),\hat{p}(x))=(\hat{u}_{0},\hat{p}_{0})\in l^{u}(x)$ it converges to $(\hat{u}_{b},\hat{p}_{b})$ for $x\rightarrow-\infty$. Moreover, if $f$ and $g$ fulfill the symmetry assumption (A2), $C^{s/u}$ posses the symmetry $C^{s}(x)=-C^{u}(-x)$ for all $x\in\mathbb{R}$. In particular, $C^{s}(0)=-C^{u}(0)$.

The proof of Proposition 2 constitutes the contents of section 2.4. Also note the similarities with Proposition 1, since the bounded solutions mentioned in both Propositions are identical up ot the scaling $\hat{u}_{b}(x)=\mu u_{b}(\xi)$.

Remark 8**.**

When $\lim_{x\rightarrow\pm\infty}f(x),g(x)=0$ (i.e. assumption (A4)), the unique bounded solution $(\hat{u}_{b},\hat{p}_{b})$ limits to the fixed point of the autonomous equation. That is,

[TABLE]

This result implies that there are trajectories on $\mathcal{M}$ that lead to and away from the bounded solution $(\hat{u}_{b},\hat{p}_{b})$. Hence, the only remaining construction steps are the analysis of persistence of orbits biasymptotic to $\mathcal{M}$ and their touch-down/take-off locations. We therefore switch back to the fast system and examine the dynamics during the jump of an orbit through the fast field. In order to pass to the reduced fast system, we use the assumption (A5) so, in the limit $\varepsilon\rightarrow 0$, we get the reduced fast system

[TABLE]

Note that in the reduced fast system the non-autonomous character of our problem is not visible. The only difference is the added trivial equation $\dot{s}=0$. As alluded to in the constant coefficient case in section 2.1, the planar subsystem $\dot{v}=q,\dot{q}=v-uv^{2}$ is known to be Hamiltonian and features a homoclinic to the saddle point $(v,q)=(0,0)$ which can be specified explicitly (see (48)). As a result, also (87) is Hamiltonian with

[TABLE]

The invariant manifold $\mathcal{M}$ from (71) is the collection of saddle points $(s,u,p,0,0),u>0,s,p\in\mathbb{R},$ for (87) and is, hence, normally hyperbolic. For its stable and unstable manifolds $W_{0}^{s/u}(\mathcal{M})$ it holds true that $\dim[W_{0}^{s/u}(\mathcal{M})]=4$ and, in fact, $W_{0}^{s}(\mathcal{M})$ and $W_{0}^{u}(\mathcal{M})$ (partly) coincide, where the intersection is simply given by the family of homoclinic orbits. Moreover, we have that $K(s,u,p,v,q)|_{(s,u,p,v,q)\in{\mathcal{M}}}=0$.

The analogy with the constant coefficient case continues for $\varepsilon>0$ sufficiently small; we still have that ${\mathcal{M}}$ is an invariant manifold of the full system (64) and that its stable and unstable manifolds persist as $W_{\varepsilon}^{s/u}(\mathcal{M})$ with $\dim[W_{\varepsilon}^{s/u}(\mathcal{M})]=4$, but do not necessarily coincide anymore. In fact, they generically meet in a 3-D intersection in $\mathbb{R}^{5}$.

Proposition 3** (Persistence of a homoclinic connection).**

Let $\varepsilon$ be sufficiently small.

Define the hyperplane $R=\{(s,u,p,v,q)~{}|~{}q=0\}$. Then $\dim[W_{\varepsilon}^{s}(\mathcal{M})\cap W_{\varepsilon}^{u}(\mathcal{M})\cap R]=2$ and orbits in this intersection fulfill $p(\xi)=\varepsilon p^{(1)}(\xi)+h.o.t.$, that is, the leading order constant term $p^{(0)}$ vanishes. 2. 2.

The take-off and touch-down surfaces on $\mathcal{M}$ of orbits in the intersection $W_{\varepsilon}^{s}(\mathcal{M})\cap W_{\varepsilon}^{u}(\mathcal{M})\cap R$ are to leading order given by

[TABLE] 3. 3.

For orbits in the intersection $W_{\varepsilon}^{s}(\mathcal{M})\cap W_{\varepsilon}^{u}(\mathcal{M})\cap R$ the touch-down curve ${T}_{d}(0)$ and stable line $l^{s}(0)$ from (82) intersect in at most two points

[TABLE]

where $C^{s}(0)$ is the slope of the stable line $l^{s}(0)$ from (82) and $\hat{u}_{b}=\mu u_{b}$ is the (rescaled) bounded background solution from Proposition 2. By symmetry, the analogous is true for the take-off curve ${T}_{o}(0)$ and unstable line $l^{u}(0)$ from (82). In particular, the thus computed $u_{0}^{\pm}$-values coincide by the aforementioned symmetry $C^{u}(0)=-C^{s}(0)$ – see Proposition 2. 4. 4.

There are two even homoclinic orbits for (7) with $u_{0}^{\pm}>0$ in case $u_{b}(0)^{2}+12/(\mu C^{s}(0))>0$ and $u_{b}(0)-\sqrt{u_{b}(0)^{2}+12/(\mu C^{s}(0))}>0$.

Remark 9**.**

If we set $u_{b}(0)=\frac{1}{\mu}$ and $C^{s}(0)=-1$ in (90), we recover (56).

Proof.

Measuring the distance of $W_{\varepsilon}^{s}({\mathcal{M}})$ and $W_{\varepsilon}^{u}({\mathcal{M}})$ in the hyperplane $R$ can again be accomplished using the difference of the Hamiltonian $K$ during the fast the jump of the orbit through the fast field (50). We have exactly as in the constant coefficient case (51) where (using that $p$ is constant to leading order) we have set $p(\xi)=p^{(0)}+\varepsilon p^{(1)}(\xi)+h.o.t.\,,$ and used that $\frac{d}{d\xi}{K}=\frac{\partial}{\partial s}K(s,u,p,v,q)(\frac{ds}{d\xi})+\frac{\partial}{\partial u}H(v,q;u)(\frac{du}{d\xi})+\frac{d}{d\xi}H(v,q;u)=0+\frac{1}{3}v^{3}(\frac{du}{d\xi})+0=\frac{1}{3}\varepsilon v^{3}p$. In order to make this difference vanish to leading order, we evidently need that $p^{(0)}=0$ and $p^{(1)}(0)=0$. This proves the first statement.

In order to construct the take-off and touch-down curves, we again investigate the change of the fast variables during the jump through the fast field:

[TABLE]

Hence, to leading order, only the $p$-variable changes during the fast jump, and therefore, the take-off and touch-down curves on ${\mathcal{M}}$ are to leading order given by (89) where we used that, by symmetry, $p(\pm 1/\sqrt{\varepsilon})=p(0)\pm\frac{1}{2}\Delta_{I_{f}}p\,.$ This proves the second statement.

Equating (89) and (82) (where we used that $p=\varepsilon\hat{p}$ – see Remark 7) gives the equality

[TABLE]

the solutions of which give the claimed expression (90) in the third statement. Finally, the fourth statement follows from inspecting (90). ∎

Two examples of homoclinic solutions for varying $f$ and $g$ can be found in Figures 6(c)–6(f). In these figures the evolution of a homoclinic solution is projected onto the manifold $\mathcal{M}$, which shows the essence of Proposition 3.

Proposition 3 thus establishes existence of homoclinic solutions for (7) under the conditions stated in Proposition 3(4). However, in the case of varying coefficients, there typically are no explicit expressions available for the bounded solution $u_{b}(0)$ and the constant $C^{s}(0)$. To circumvent this, in the next section we derive bounds on these using the theory of exponential dichotomy, which simultaneously forms the proof of Proposition 2.

2.3 Some basic results from the theory of exponential dichotomies

When $f$ and/or $g$ are non-constant, generically it is not possible to capture the dynamics on manifold $\mathcal{M}$ in explicit expressions. Instead, our main tools for constructing a saddle-like structure on $\mathcal{M}$ are from the theory of exponential dichotomies. To fix notation and keep the exposition self-contained, we state (following [11]) the definition of exponential dichotomies along with a selection of results that we use here.

Definition 1** (Exponential Dichotomy).**

Consider the planar ODE $\frac{d}{dx}Y=B(x)Y$ for the unknown $Y:\mathbb{R}\rightarrow\mathbb{R}^{2}$ and with $B:\mathbb{R}\rightarrow\mathbb{R}^{2\times 2}$ a matrix-valued function which is continuous on $\mathbb{R}$. Let $\Phi=\Phi(x)$ be the associated canonical solution operator. This ODE is said to have an exponential dichotomy if there is a projection matrix $P$ and positive constants $K$ and $\rho$ such that

[TABLE]

In the next section we will be interested in first order ODEs of the form

[TABLE]

with $x\in\mathbb{R},Y:\mathbb{R}\rightarrow\mathbb{R}^{2},A_{0}\in\mathbb{R}^{2\times 2},A:\mathbb{R}\rightarrow\mathbb{R}^{2\times 2},F\in\mathbb{R}^{2}$. In particular, we would like to corroborate knowledge of the autonomous version (which is often available in terms of explicit solutions) to deduce qualitative results for the full non-autonomous one. For the sake of clarity, we assemble first all auxiliary systems in one place:

[TABLE]

Proposition 4** (Roughness and closeness of bounded solutions).**

Let $K_{aut},\rho_{aut}>0$ be the exponential dichotomy constants of the homogeneous, autonomous ODE (97) and $\Phi_{aut},P_{aut}$ the corresponding solution and projection operators. If

[TABLE]

the non-autonomous ODE (98) has an exponential dichotomy for which the following holds true.

(Roughness)* The exponential dichotomy constants of the homogeneous, non-autonomous ODE (98) are $K=\frac{5}{2}K_{aut}^{2}$ and ${\rho}=\rho_{aut}-2K_{aut}\delta$, and concerning the solution and projection operators ${\Phi},{P}$ of (98) we have upon defining*

[TABLE]

the estimate

[TABLE] 2. 2.

(Closeness of bounded solutions)* There exist unique bounded solutions $Z_{b,aut},Y_{b}$ of the inhomogeneous, autonomous and non-autonomous ODEs (99) and (95). In particular, they satisfy*

[TABLE]

Proof.

The first statement is the persistence of exponential dichotomies, known as “roughness”, and is a standard result (see [11, Ch.4, Prop.1]). Moreover, another standard result from the theory of exponential dichotomies stipulates that inhomogeneous equations have unique bounded solutions, when the homogeneous equations have an exponential dichotomy and the inhomogeneous terms are bounded (see [11, Ch.8, Prop.2]). Then, to demonstrate the rest of the second statement, we define $W(x)={Y}_{b}(x)-Z_{b,aut}(x)$ which gives $W^{\prime}(x)=A_{0}W(x)+G(x)$ with $G(x)=A(x){Y}_{b}(x)$. The unique bounded solution $W_{b}$ of this ODE satisfies the estimate $\sup_{x\in\mathbb{R}}\|W_{b}(x)\|\leq\frac{2{K}_{aut}}{{\rho_{aut}}}\,\sup_{x\in\mathbb{R}}\|G(x)\|\leq\frac{4\delta K_{aut}{K}}{\rho_{aut}{\rho}}\,\|F\|\,,$ where we used that $\sup_{x\in\mathbb{R}}\|{Y}_{b}(x)\|\leq\frac{2{K}}{{\rho}}\|F\|\,.$

∎

2.4 Dynamics on $\mathcal{M}$ (Proof of Proposition 2)

Let us introduce the more concise notation $Y=\left(\hat{u},\frac{d}{dx}\hat{u}\right)^{T}$ such that (81) has the form of (95) from the previous section; that is,

[TABLE]

with

[TABLE]

Lemma 3** (Exponential Dichotomy Constants and Roughness).**

With the notation of Proposition 4, let

[TABLE]

Then we have $\rho_{aut}=K_{aut}=1,\rho=1-2\delta,K=5/2$ and

[TABLE]

Proof.

We have the canonical solution operator $\Phi(x)=e^{A_{0}x}$. The eigenvalues of the matrix $A_{0}$ are $\pm 1$ and the corresponding normed eigenvectors are $v=\frac{1}{\sqrt{2}}(1,1)^{T},w=\frac{1}{\sqrt{2}}(1,-1)^{T}$. Thus the fixed point $Y=(0,0)^{T}$ is a saddle. From this it is clear that we can choose

[TABLE]

With the basis transformation matrix $B=\left(v~{}|~{}w\right)$ and the diagonal matrix $D=\mbox{diag}(1,-1)$ we then get

[TABLE]

A similar reasoning – where one can use that $I-P=vv^{T}$ – gives

[TABLE]

Thus we have the estimate for exponential dichotomies from Definition 1 with $\rho_{aut}=1$ and $K_{aut}=1$. The remaining statements can now be read off Proposition 4. ∎

The roughness of exponential dichotomies established in Lemma 3 provides a bound on the projection operator $Q(x)$ of the non-autonomous system. However, this bound cannot be used directly to prove existence of homoclinic solutions using geometric singular perturbation theory, as geometric properties need to be derived. In particular, we need to find the stable and unstable manifolds for the unique bounded solution $Y_{b}$ of (104). These can be defined as

[TABLE]

where $\Phi,P$ are the solution and projection operator for (104). For the construction that we have in mind, it is convenient to notice that

[TABLE]

with lines

[TABLE]

While, in general, it is not possible to find explicit expressions for these objects, we can derive estimates for their locations. For this we first observe that the line $l^{s}$ can be written equivalently as

[TABLE]

where $C(x)$ is the slope of the line. Starting from the bound on the projection operator $Q(x)=\Phi(x)P\Phi^{-1}(x)$ derived in Lemma 3, a bound on the projection lines will be established in Lemma 4, which is then subsequently used to find a bound on the slope $C(x)$ via the angle $\theta(x)$ of the line in Lemma 5.

In particular, for the case of (81), we thus obtain

[TABLE]

with $\tilde{C}(x)$ as in Lemma 5 taking into account that the projection operator depends on $x$, that is, $Q=Q(x)$ and so does the angle $\theta=\theta(x)$, which defines $C=C(x)$ and, hence, also $\tilde{C}=\tilde{C}(x)$.

The rest of this section consists of the two technical lemmas that ultimately derive a bound for $\tilde{C}$.

Lemma 4** (Closeness of projection lines).**

Let $Q$ and ${Q}_{aut}$ be the projection matrices with rank $1$ as defined in Proposition 4(i), i.e. there are unit vectors $q$ and $q_{aut}$ such that $Q=qq^{T}$ and ${Q}_{aut}={q}_{aut}{q}_{aut}^{T}$, and $\|Q-{Q}_{aut}\|<4\delta$ holds true. Then either $\|q-{q}_{aut}\|<\sqrt{8\delta}$ or $\|q+{q}_{aut}\|<\sqrt{8\delta}$.

Proof.

We prove the equivalent statement that from $\|q-q_{aut}\|\geq\sqrt{8\delta}$ and $\|q+q_{aut}\|\geq\sqrt{8\delta}$ it follows that $\|Q-Q_{aut}\|\geq 4\delta$. First we observe that

[TABLE]

Therefore, by assumption

[TABLE]

from which it follows that $\|Q-Q_{aut}\|\geq 4\delta$. ∎

The previous lemma establishes closeness of projection lines of the autonomous and the non-autonomous case. The thus obtained bounds on norms can be transferred to bounds on the slope $C$ by use of elementary geometry. Note that transforming the norm bounds in this way leads to singularities when a projection line passes the vertical axis (which also leads to a seemingly disjoint set of admittable slopes). A visualisation of the results of Lemma 5 are given in Figure 7. In particular, the resulting bounds for the slope are shown.

Lemma 5** (Closeness of slopes).**

Let $Q$ and ${Q}_{aut}$ be the projection matrices with rank $1$ as defined in Proposition 4(i), i.e. there are unit vectors $q$ and $q_{aut}$ such that $Q=qq^{T}$ and ${Q}_{aut}={q}_{aut}{q}_{aut}^{T}$, and $\|Q-Q_{aut}\|<4\delta$ holds true. Furthermore, let $\theta,\theta_{aut}\in[-\pi,\pi)$ be defined by $q=:(\cos(\theta),\sin(\theta)),q_{aut}=(\cos(\theta_{aut}),\sin(\theta_{aut}))$ such that the slopes of the lines spans by $q$ and $q_{aut}$ are given by

[TABLE]

Then there exist constants $C_{min/max}(\delta,C_{aut})$ defined by

[TABLE]

such that $C-C_{aut}\in\Gamma\left(\delta,C_{aut}\right)$, where

[TABLE]

In particular, for $q_{aut}=\frac{1}{\sqrt{2}}(1,-1)^{T}$ we have $C_{aut}=-1$ and, hence,

[TABLE]

Proof.

For technical reasons we assume that $\|q-q_{aut}\|\leq\|q+q_{aut}\|$; if this inequality does not hold, we can scale $q\rightarrow-q$ without changing the projection matrix $Q$. Then, with

[TABLE]

we have

[TABLE]

From $\|Q-{Q}_{aut}\|<4\delta$ we know by the previous lemma that $\|q-{q}_{aut}\|<\sqrt{8\delta}$ and, hence, since $q$ and $q_{aut}$ are unit vectors, we have

[TABLE]

Since $\mathrm{arccos}(z)$ is monotonically decreasing, we hence get from $|\Delta\theta|=\mathrm{arccos}(q^{T}q_{aut})$ that

[TABLE]

Furthermore, since $\frac{\tan(z)}{1-C_{aut}\tan(z)}$ is monotonically increasing in $z$, we have the claimed result by using

[TABLE]

and some simplifications in (122). ∎

2.5 Existence results

Here, we first state our main existence results in detail. Their proofs are given in section 2.6.

Theorem 1** (Existence for general $f,g$).**

Let assumptions (A1), (A2), (A3) and (A5) be satisfied. Then there is a $\mu^{*}$ with $0<\mu^{*}<\frac{1}{12}$ and corresponding $\varepsilon^{*}=\varepsilon^{*}(\mu)>0,0<\delta^{*}=\delta^{*}(\mu)<\frac{2-\sqrt{2}}{8}$ such that the following holds true: For any $\varepsilon,\mu,\delta$ with

[TABLE]

the stationary wave ODE (64) has (two) orbits $\left(s_{p}(\xi),u_{p}(\xi),p_{p}(\xi),v_{p}(\xi),q_{p}(\xi)\right)$, that are homoclinic to the bounded solution $\left(\xi,\frac{\hat{u}_{b}(\varepsilon^{2}\mu\xi)}{\mu},\varepsilon\hat{u}_{b}^{\prime}(\varepsilon^{2}\mu\xi),0,0\right)$, with $\left(u_{p}(\xi),v_{p}(\xi)\right)$ to leading order given by

[TABLE]

with $u_{0}=u_{0}^{-}$ or $u_{0}=u_{0}^{+}$ from (90), i.e.

[TABLE]

$\hat{u}_{b}$* the bounded solution from Proposition 2 and where the indicator functions*

[TABLE]

distinguishes the behavior of the solution in the fast and super-slow fields. Furthermore, for $\hat{u}_{\pm}$ we have the estimates

[TABLE]

for some $C>0$, the bounded solution $u_{b}$ obeys

[TABLE]

Finally, this homoclinic orbit gives rise to a stationary pulse solution

[TABLE]

for the Klausmeier model (4) that is biasymptotic to the bounded state $\left(a\hat{u}_{b}\left(\frac{\sqrt{m}}{D}x\right),0\right)$.

Corollary 1** (Existence for $f,g=0$).**

Let $f,g=0$, and the conditions from Theorem 1 be fulfilled. Then

[TABLE]

Corollary 2** (Existence for small $f,g$).**

Let the conditions from Theorem 1 be fulfilled and $f=\delta\tilde{f}$, $g=\delta\tilde{g}$ where $\tilde{f},\tilde{g}=\mathcal{O}(1),0<\delta\ll 1$ (i.e. $\sup_{x\in\mathbb{R}}\sqrt{\tilde{f}(x)^{2}+\tilde{g}(x)^{2}}=1$. Then

[TABLE]

Moreover, $u_{0}$ as in (133) can be expressed in terms of $\delta$ as

[TABLE]

where $u_{00}$ corresponds to the $u_{0}$-value for the autonomous case, i.e. $u_{00}$ is given by (56).

Corollary 3** (Existence for $h(x)=-2\ln\cosh(\beta x)$).**

Let $h(x)=-2\ln\cosh(\beta x),\beta>0,f=h^{\prime},g=h^{\prime\prime}$, and the conditions from Theorem 1 be fulfilled. Then

[TABLE]

Remark 10**.**

Pulses solutions as in Corollary 3 exist for any $\beta>0$ without the need of the general assumption on $\delta$ as in Theorem 1; since the flow on $\mathcal{M}$ can be solved explicitly for these functions $f$ and $g$, no condition on $\delta$ is needed.

Remark 11**.**

Since the flow on $\mathcal{M}$ can be solved explicitly for the functions $f$ and $g$ as in Corollary 3, it is also possible to prove existence of symmetric, stationary $2$-pulse solutions (and, in fact, any symmetric, stationary $N$-pulse solution). Note that normally, for $f,g\equiv 0$, these do not exist, since pulses in (4) repel each other [13, 4]; this repulsive force can only be overcome by driving forces due to the spatially varying functions $f$ and $g$. We come back to these multi-pulse solutions in section 4.5.

2.6 Proof of existence results

The proofs of the existence results in section 2.5 follow from the theory developed in the preceding sections. The heart of these proofs is formed by Proposition 3 and the bounds on the bounded solution $u_{b}$ and the slopes $C^{s/u}$ as found in Proposition 2. Ultimately, it boils down to taking $\delta$ small enough such that an intersection between $l^{s}(0)$ and $T_{o}(0)$ is guaranteed. A sketch of this idea is given in Figure 3; the rest of this section is devoted to the rigorous proof of the existence theorem and the corollories in section 2.5.

Proof of Theorem 1.

Existence of the homoclinic orbits is established by Proposition 3 if the conditions in Proposition 3(4) are satisfied. Since $u_{b}(0)=\hat{u}_{b}(0)/\mu$, these hold if and only if the following three bounds hold true:

• (i)

$\hat{u}_{b}(0)>0$;

• (ii)

$C^{s}(0)<0$;

• (iii)

$\hat{u}_{b}(0)^{2}+12\mu/C^{s}(0)>0$.

By Proposition 4 and Lemma 3, we have

[TABLE]

and by Lemma 5 we have

[TABLE]

where $\Gamma$ is as in (119). Using these, bound (i) is satisfied when $\delta<\frac{1}{12}$ and bound (ii) when $\delta<\frac{2-\sqrt{2}}{8}$. Since the bound (iii) holds true when $\delta=0$ and $\mu<\frac{1}{12}$, continuity of mentioned bounds on $\hat{u}_{b}(0)$ and $C^{s}(0)$ guarantees the existence of the critical value $0<\delta^{*}(\mu)<\frac{2-\sqrt{2}}{8}$. ∎

Proof of Corollary 1.

This follows immediately from solving (81) with $f,g\equiv 0$, and is also carried out in more detail in section 2.1. ∎

Proof of Corollary 2.

The super-slow system on $\mathcal{M}$ in (81) can be solved using a regular expansion in $0<\delta\ll 1$. By requiring that $\lim_{x\rightarrow\infty}\hat{u}_{+}(x)$ and $\lim_{x\rightarrow-\infty}\hat{u}_{-}(x)$ exist, the results follow by a straightforward calculation. ∎

Proof of Corollary 3.

One can easily verify that $\hat{u}_{\pm}$ solve (81), and that $\lim_{x\rightarrow\pm\infty}\hat{u}_{\pm}(x)=0$. The bounded solution $\hat{u}_{b}$ follows from a standard variation of constants method. ∎

3 Linear stability analysis

In the previous section, we proved the existence of stationary $1$-pulse solutions to (4). In this section we study the linear stability of these solutions. For $(U_{p},V_{p})$ a pulse solution from Theorem 1 we define the linear operator

[TABLE]

with $\mathcal{L}:H^{2}(\mathbb{R})\times H^{2}(\mathbb{R})\subset L^{2}(\mathbb{R})\times L^{2}(\mathbb{R})\rightarrow L^{2}(\mathbb{R})\times L^{2}(\mathbb{R})$ and its spectrum by $\Sigma(\mathcal{L})$, where we distinguish between the point spectrum $\Sigma_{\mathrm{pt}}(\mathcal{L})$ and the essential spectrum $\Sigma_{\mathrm{ess}}(\mathcal{L})=\Sigma(\mathcal{L})\setminus\Sigma_{\mathrm{pt}}(\mathcal{L})$ – we denote the elements of $\Sigma_{\mathrm{ess}}(\mathcal{L})$ by $\underline{\lambda}$. As customary, we say that $(U_{p},V_{p})$ is linearly stable if there is no spectrum in the right half plane. In order to keep the exposition at reasonable length, we will concentrate here on characterizing parameter regimes where the only instability that can occur is through the (translational) zero eigenvalue which starts moving due to the introduction of spatially varying $f$ and/or $g$. In particular, there are no essential instabilities:

Lemma 6** (Essential spectrum).**

Let the conditions of Theorem 1 and assumption (A4) be fulfilled, and let $(U_{p},V_{p})$ be a pulse solution to (4) as in Theorem 1. Then the essential spectrum of $\mathcal{L}$ from (143) is

[TABLE]

and, hence, lies in the left half-plane.

Proof.

The limiting operator of $\mathcal{L}$ at $x\rightarrow\pm\infty$ is $\mathcal{L}_{\infty}:=\mathrm{diag}[\partial_{x}^{2}-1,D^{2}\partial_{x}^{2}-m]$ (note that we thus explicitly use assumption (A4)). Therefore, we have that the boundaries of the essential spectrum are $\underline{\lambda}_{1}(k)=-(k^{2}+1)$, $\underline{\lambda}_{2}(k)=-(D^{2}k^{2}+m)$, $k\in\mathbb{R}$, which immediately gives the claimed result. ∎

The assumptions on $f,g$ allow (again through the use of exponential dichotomies) the derivation of bounds on the location of the point spectrum, which, under the assumption that $f,g$ are chosen ‘small’, can be further refined to track the one small eigenvalue that can possibly lead to bifurcations. The proof of the following statements will be the subject of the next sections.

Theorem 2** (Point spectrum).**

Let the conditions of Theorem 1 and assumption (A4) be fulfilled, and let $(U_{p},V_{p})$ be a pulse solution to (4) with $u_{0}=u_{0}^{-}$ as in (133). Then there exist constants $m_{c},\mu^{*},\nu^{*}>0$ such that if either (i) $m<m_{c}$ and $\mu<\mu^{*}$ or (ii) $m>m_{c}$ and $\mu\sqrt{m}<\nu^{*}$, then there exists a $\delta_{c}>0$ such that if $0\leq\delta<\delta_{c}$ precisely one eigenvalue $\underline{\lambda}_{0}$ is $\mathcal{O}(\varepsilon)$-close to [math] and all other eigenvalues of $\mathcal{L}$ lie in the left-half plane.

Proof.

The statement is demonstrated in section 3.1 by combining the setup of an Evans function and the theory of exponential dichotomies. ∎

Remark 12**.**

Note that Theorem 2 only holds for pulse solutions with $u_{0}=u_{0}^{-}$; pulse solutions with $u_{0}=u_{0}^{+}$ are always unstable. See also Remark 18.

Remark 13**.**

The constants $m_{c}$, $\mu^{*}$ and $\nu^{*}$ in Theorem 2 can be computed explicitly (see Lemma 10).

Theorem 3** (Small eigenvalue close to $\underline{\lambda}=0$ for small $f$, $g$).**

Assuming that $f=\delta\tilde{f}$, $g=\delta\tilde{g}$ with $0<\delta\ll 1$, $\tilde{f},\tilde{g}=\mathcal{O}(1)$ (i.e. $\sup_{x\in\mathbb{R}}\sqrt{\tilde{f}(x)^{2}+\tilde{g}(x)^{2}}=1$), there exists a constant $\tau^{*}>0$ such that if $\tau:=\varepsilon^{4}\mu m<\tau^{*}$ the small eigenvalue $\underline{\lambda}_{0}$ close to $\underline{\lambda}=0$ is located, to leading order, at

[TABLE]

where $u_{0}$ is as in (133) and Corollary 2.

Proof.

This statement is derived in section 3.2 by employing a regular expansion in $\delta$. ∎

Corollary 4**.**

Let the conditions of Theorem 3 be fulfilled. Then, in the double asymptotic limit $\mu\ll 1$ and $\tau:=\varepsilon^{4}\mu m\ll 1$ the leading order expression for $\underline{\lambda}_{0}$ becomes

[TABLE]

Remark 14**.**

When the term $\tau=\varepsilon^{4}\mu m=\frac{a^{2}D}{m\sqrt{m}}$ in (145) becomes too large (larger than $\tau^{*}$), the pulse becomes unstable due to a traveling wave bifurcation/drift instability [10, 14].

3.1 Qualitative description of the point spectrum location (Proof of Theorem 2)

This section is devoted to finding the point spectrum of the operator $\mathcal{L}$. For that, we use a decomposition method for the Evans function, first developed in [1, 16], which is supplemented by the theory of exponential dichotomies to treat the varying coefficients in (4). As before, the following computations will again heavily rely on the singularly perturbed structure. Therefore, we introduce for the eigenvalue problem $(\mathcal{L}-\underline{\lambda}I)(\bar{U},\bar{V})^{T}=0$, that is,

[TABLE]

and the scalings (analogous to (14) and (18))

[TABLE]

to get the fast eigenvalue problem

[TABLE]

which suggests (just as in [4, 10, 14]) the introduction of the scaled eigenvalue parameter

[TABLE]

so, finally,

[TABLE]

It is convenient to introduce $\phi:=\left(\bar{u},\dot{\bar{u}}/(\varepsilon^{2}\mu),\bar{v},\dot{\bar{v}}\right)$ and to write the above ODEs as the system of first order ODEs

[TABLE]

where

[TABLE]

From the existence analysis in section 2, we have seen that the real line $\mathbb{R}$ can be split in one fast region, $I_{f}$, near the pulse location and two super slow fields $I_{s}^{\pm}$ to both sides of the fast field:

[TABLE]

Since we know that $v_{p}$ vanished to leading order in the slow fields, we have in those regions the system matrix

[TABLE]

that is, the dynamics for slow and fast variables are decoupled. Any value $\lambda\in\mathbb{C}$ for which this system of ODEs has a non-trivial solution in $L^{2}(\mathbb{R})\times L^{2}(\mathbb{R})$ corresponds to an eigenvalue $\underline{\lambda}=m\lambda$ of $\mathcal{L}$. A mechanism (that is by now standard) for detecting eigenvalues is the construction of an Evans function, whose roots coincide with the eigenvalues of $\mathcal{L}$. Although the Evans function can also be extended into the essential spectrum, we do not need this in the present work and rather restrict $\lambda$ to

[TABLE]

on which the Evans function is analytic.

3.1.1 Evans function construction

By (conditions and results of) Theorem 1 and assumption (A4), we know that the limiting matrix for $|\xi|\rightarrow\infty$ is given by

[TABLE]

Its eigenvalues $\Lambda_{1,2,3,4}$ and eigenvectors $E_{1,2,3,4}$ are

[TABLE]

where $\mbox{Re}\left(\Lambda_{1}(\lambda)\right)<\mbox{Re}\left(\Lambda_{2}(\lambda)\right)<0<\mbox{Re}\left(\Lambda_{3}(\lambda)\right)<\mbox{Re}\left(\Lambda_{4}(\lambda)\right)$ for $\lambda\in\mathcal{C}_{e}$.

The system $\dot{\phi}_{\infty}=A_{\infty}(\lambda,\varepsilon,\mu,m)\phi_{\infty}$ admits exponential dichotomies on $\mathcal{C}_{e}$. Since $A_{\infty}$ is exponentially close to $A$ for large $|\xi|$, the stable and unstable subspaces of $\dot{\phi}=A(\xi;\lambda,\varepsilon,\mu,m)\phi$ and $\dot{\phi}_{\infty}=A_{\infty}(\lambda,\varepsilon,\mu,m)\phi_{\infty}$ are similar when $|\xi|\rightarrow\infty$. In particular, for all $\lambda\in\mathcal{C}_{e}$ there is a two-dimensional family of solutions, $\Phi_{\infty}^{-}(\lambda)$, to $\dot{\phi}_{\infty}=A_{\infty}(\lambda,\varepsilon,\mu,m)\phi_{\infty}$ such that $\lim_{\xi\rightarrow-\infty}\phi_{\infty}^{-}(\xi)=0$ for all $\phi_{\infty}^{-}\in\Phi_{\infty}^{-}(\lambda)$, and a two-dimensional family of solutions, $\Phi_{\infty}^{+}(\lambda)$, to $\dot{\phi}_{\infty}=A_{\infty}(\lambda,\varepsilon,\mu,m)\phi_{\infty}$ such that $\lim_{\xi\rightarrow\infty}\phi_{\infty}^{+}(\xi)=0$ for all $\phi_{\infty}^{+}\in\Phi_{\infty}^{+}(\lambda)$, which implies that the system $\dot{\phi}=A(\xi;\lambda,\varepsilon,\mu,m)\phi$ also possesses two two-dimensional families of solutions, $\Phi^{-}(\lambda)$ and $\Phi^{+}(\lambda)$ with the same properties.

For the system $\dot{\phi}=A(\xi;\lambda,\varepsilon,\mu,m)\phi$, however, it is possible that the intersection $\Phi^{+}(\lambda)\cap\Phi^{-}(\lambda)$ is nonempty. The values $\lambda\in\mathcal{C}_{e}$ for which this happens correspond to $\underline{\lambda}=m\lambda$ in the point spectrum $\Sigma_{\mathrm{pt}}$. To find these, we use a Evans function [1, 16], which is defined as

[TABLE]

where $\{\phi_{1}(\cdot;\lambda),\phi_{2}(\cdot;\lambda)\}$ spans the space $\Phi^{-}(\lambda)$ and $\{\phi_{3}(\cdot;\lambda),\phi_{4}(\cdot;\lambda)\}$ spans the space $\Phi^{+}(\lambda)$. For notational clarity we have suppressed the dependence on the other parameters. Essentially, the Evans function $\mathcal{D}(\lambda)$ measures the linear independence of the solution functions $\phi_{1,\ldots,4}$. Therefore, zeros of $\mathcal{D}(\lambda)$ correspond to values of $\lambda$ for which $\Phi^{+}(\lambda)\cap\Phi^{-}(\lambda)\neq\varnothing$, and thus to eigenvalues in the point spectrum [1].

In (160) the solutions $\phi_{1,\ldots,4}$ are not uniquely defined, and any choice leads to the same eigenvalues. However, for singularly perturbed partial differential equations a specific choice enables the use of the scale separation in these equations, which in turn makes it possible to determine the eigenvalues.

Lemma 7**.**

Let the conditions of Theorem 1 be fulfilled and let $(U_{p},V_{p})$ be a pulse solution to (4) as in Theorem 1. Then all eigenvalues $\lambda\in\Sigma_{\mathrm{pt}}$ associated to (153) are roots of the Evans function

[TABLE]

where $t_{11}$ and $t_{22}$ are analytic (transmission) functions of $\lambda$, defined by

[TABLE]

where $\phi_{1}$ is the (unique) solution to (154) for which

[TABLE]

Proof.

The proof is heavily based on [16, Section 3.2]. Therefore, we present here only an outline of the proof and refer the interested reader to [16] for more details.

The heart of the proof is based on choosing $\phi_{1,\ldots,4}$ in such way that the scale separation of (4) can be exploited. Because $A$ and $A_{\infty}$ are exponentially close when $\xi\rightarrow-\infty$, there is a unique solution $\phi_{1}$ such that $\phi_{1}$ closely follows $E_{1}(\lambda)e^{\Lambda_{1}(\lambda)\xi}$ as $\xi\rightarrow-\infty$. More precisely, we define $\phi_{1}$ uniquely such that $\lim_{\xi\rightarrow-\infty}\phi_{1}(\xi;\lambda)e^{-\Lambda_{1}(\lambda)\xi}=E_{1}(\lambda)$. For $\xi\rightarrow\infty$, we do not know the precise form of $\phi_{1}$, but we do know that, asymptotically, it is a combination of the eigenfunctions of the system $\dot{\phi}_{\infty}=A_{\infty}\phi_{\infty}$. That is, $\phi_{1}(\xi;\lambda)\rightarrow t_{11}(\lambda)E_{1}e^{\Lambda_{1}(\lambda)\xi}+t_{12}(\lambda)E_{2}e^{\Lambda_{2}(\lambda)\xi}+t_{13}(\lambda)E_{3}e^{\Lambda_{3}(\lambda)\xi}+t_{14}(\lambda)e^{\Lambda_{4}(\lambda)\xi}$ as $\xi\rightarrow\infty$, where $t_{11},\ldots,t_{14}$ are analytic transmission functions.

Next, $\phi_{2}$ must be chosen such that $\{\phi_{1}(\cdot,\lambda),\phi_{2}(\cdot,\lambda)\}$ spans $\Phi^{-}(\lambda)$. As this does not determine $\phi_{2}$ uniquely, we may, additionally, require that $\phi_{2}$ grows, at most, as $E_{2}(\lambda)e^{\Lambda_{2}(\lambda)\xi}$ for $\xi\rightarrow\infty$. More precisely, we define $\phi_{2}$ uniquely such that $\lim_{\xi\rightarrow-\infty}\phi_{2}(\xi;\lambda)e^{-\Lambda_{2}(\lambda)\xi}=E_{2}$ and $\lim_{\xi\rightarrow+\infty}\phi_{2}(\xi;\lambda)e^{-\Lambda_{1}(\lambda)\xi}=0$ (note that this construction is based on insight in $t_{11}$ – that may not be [math] – that is obtained by the ‘elephant trunk procedure’, see [16, 25] and Remark 20). For $\xi\rightarrow\infty$, $\phi_{2}$ is then asymptotically given by $\phi_{2}(\xi;\lambda)\rightarrow t_{22}(\lambda)E_{2}(\lambda)e^{\Lambda_{2}(\lambda)\xi}+t_{23}(\lambda)E_{3}(\lambda)e^{\Lambda_{3}(\lambda)\xi}+t_{24}(\lambda)e^{\Lambda_{4}(\lambda)\xi}$ as $\xi\rightarrow\infty$, where $t_{21}$, $t_{23}$, $t_{24}$ are analytical transmission functions.

In a similar vein the solutions $\phi_{3}$ and $\phi_{4}$ can be defined such that $\lim_{\xi\rightarrow\infty}\phi_{4}(\xi;\lambda)e^{-\Lambda_{4}(\lambda)}=E_{4}(\lambda)$ and $\lim_{\xi\rightarrow\infty}\phi_{3}(\xi;\lambda)e^{-\Lambda_{3}(\lambda)}=E_{3}(\lambda)$.

Then, using that $\sum_{j=1}^{4}\Lambda_{j}(\lambda)=0$ and by Liouville’s formula, the Evans function (160) can be rewritten:

[TABLE]

∎

The roots $\lambda\in\mathcal{C}_{e}$ of $\mathcal{D}(\lambda)$ thus correspond to the roots of $t_{11}(\lambda)t_{22}(\lambda)$. The next goal, therefore, is to determine the roots of these transmission functions.

3.1.2 Fast transmission function $t_{11}$

The transmission function $t_{11}$ is closely related to the linearization around the pulse in the fast field,

[TABLE]

The eigenvalues of $\mathcal{L}^{\mathrm{r}}$ are well-known to be $\lambda_{0}^{\mathrm{r}}=5/4$, $\lambda_{1}^{\mathrm{r}}=0$ and $\lambda_{2}^{\mathrm{r}}=-3/4$. By a standard winding number argument, it follows that roots of $t_{11}$ lie $\mathcal{O}(\varepsilon)$-close to these eigenvalues $\lambda_{0}^{\mathrm{r}}$, $\lambda_{1}^{\mathrm{r}}$ and $\lambda_{2}^{\mathrm{r}}$.

Lemma 8** (Properties of $t_{11}$).**

Let the conditions of Proposition 7 be fulfilled. The roots of $t_{11}$ lie $\mathcal{O}(\varepsilon)$ close to the eigenvalues (counting multiplicity) of $\mathcal{L}^{\mathrm{r}}$, i.e. close to $\lambda_{0}^{\mathrm{r}}=5/4$, $\lambda_{1}^{\mathrm{r}}=0$ and $\lambda_{2}^{\mathrm{r}}=-3/4$.

Proof.

See [16, Lemma 4.1]. ∎

Although $t_{11}$ has a root (with multiplicity $1$) close to $\lambda^{r}_{0}=5/4$, this does not mean that $\mathcal{D}(\lambda)$ has a root for the same value of $\lambda$, since – as will be discussed in the next section – the transmission function $t_{22}$ has a pole of order $1$ for the same $\lambda$, thus preventing it from being an eigenvalue of $\mathcal{L}$ – in the literature, this is known as the ‘NLEP paradox’.

In studies of autonomous systems, the root of $t_{11}$ close to $\lambda=0$ is actually located precisely at $\lambda=0$ because of the translation invariance of those autonomous systems. However, (4) is non-autonomous and therefore this reasoning no longer holds and the eigenvalue close to $\lambda^{r}_{1}=0$ can have negative or positive real part. As $t_{22}$ does not have a pole for this $\lambda$ – as will be discussed in the next section – the Evans function $\mathcal{D}(\lambda)$ has a root for this value; it thus corresponds to an eigenvalue of $\mathcal{L}$. To our best knowledge, it is, in general, not possible to determine the precise location of this eigenvalue; in section 3.2 we compute its location using standard regular perturbation techniques when the non-autonomous terms are small.

3.1.3 Slow transmission function $t_{22}$

To determine the transmission function $t_{22}$, we focus on the function $\phi_{2}$, as defined in Proposition 7. Per construction, we know that $\phi_{2}(\xi;\lambda)\rightarrow t_{22}(\lambda)E_{2}(\lambda)e^{\Lambda_{2}(\lambda)\xi}+t_{23}(\lambda)E_{3}(\lambda)e^{\Lambda_{3}(\lambda)\xi}+t_{24}(\lambda)e^{\Lambda_{4}(\lambda)\xi}$ as $\xi\rightarrow\infty$. As $|\Lambda_{4}(\lambda)|\gg|\Lambda_{2,3}(\lambda)|$ for $\lambda\in\mathcal{C}_{e}$, the term $e^{\Lambda_{4}(\lambda)\xi}$ is exponentially small in the slow fields $I_{s}^{\pm}$. Therefore, we have $\phi_{2}(\xi;\lambda)\approx t_{22}(\lambda)E_{2}(\lambda)e^{\Lambda_{2}(\lambda)\xi}+t_{23}(\lambda)E_{3}(\lambda)e^{\Lambda_{3}(\lambda)\xi}$ for $\xi\in I_{s}^{+}$ sufficiently large. In this way, $\phi_{2}$ in the slow fields is related to the properties of the exponentially asymptotic constant-coefficient system $\dot{\phi}_{\infty}=A_{\infty}(\lambda,\varepsilon,\mu,m)\phi_{\infty}$. However, we need to relate $\phi_{2}$ in the slow fields to the exponentially asymptotic non-autonomous system $\dot{\phi}_{s}=A_{s}(\xi;\lambda,\varepsilon,\mu,m)\phi_{s}$ to determine $t_{22}$.

In the slow fields the system $\dot{\phi}_{s}=A_{s}(\xi;\lambda,\varepsilon,\mu,m)\phi_{s}$ has the dynamics for the $(\bar{u},\bar{p})$ part completely separated from the dynamics of the $(\bar{v},\bar{q})$ part. The $(\bar{u},\bar{p})$ part is governed by the non-autonomous ODE

[TABLE]

where

[TABLE]

Here, only the matrix $B_{1}$ carries the non-autonomous part of the differential equation and the system without $B_{1}$ corresponds to the $(\bar{u},\bar{p})$ part of the system $\dot{\phi}_{\infty}=A_{\infty}(\lambda,\varepsilon,\mu,m)\phi_{\infty}$, which has spatial eigenvalues $\Lambda_{2,3}=\pm\varepsilon^{2}\mu\ \sqrt{1+m\lambda}$. When $\lambda\in\mathcal{C}_{e}$ this autonomous system admits an exponential dichotomy on $\mathbb{R}$ and, therefore, by roughness the non-autonomous system (168) does so as well, provided that $\delta=\sup_{x\in\mathbb{R}}\sqrt{f(x)^{2}+g(x)^{2}}=\sup_{x\in\mathbb{R}}\|B_{1}(x)\|$ is sufficiently small. Under these conditions, there exist $\tilde{\psi}_{2}(\xi;\lambda)=(u_{2}(\xi;\lambda),p_{2}(\xi;\lambda),0,0)^{T}$ and $\tilde{\psi}_{3}(\xi;\lambda)=(u_{3}(\xi;\lambda),p_{3}(\xi;\lambda),0,0)^{T}$ such that $\tilde{\psi}_{2}(\xi;\lambda)\rightarrow E_{2}(\lambda)e^{\Lambda_{2}(\lambda)\xi}$ and $\tilde{\psi}_{3}(\xi;\lambda)\rightarrow E_{3}(\lambda)e^{\Lambda_{3}(\lambda)\xi}$ as $|\xi|\rightarrow\infty$. The same reasoning as before can now be used to deduce that $\phi_{2}(\xi;\lambda)\approx\tilde{\psi}_{2}(\xi;\lambda)$ for $\xi\in I_{s}^{-}$ and $\phi_{2}(\xi;\lambda)\approx t_{22}(\lambda)\tilde{\psi}_{2}(\xi;\lambda)+t_{23}(\lambda)\tilde{\psi}_{3}(\xi;\lambda)$ for $\xi\in I_{s}^{+}$.

To compute $t_{22}$ we need to track the changes of $\bar{u}$ and $\bar{p}$ during the fast transition when $\xi\in I_{f}$. From (153), it follows that $\bar{u}$ stays constant to leading order. Hence, matching $\phi_{2}$ at the ends of both super-slow fields $I_{s}^{\pm}$ gives the leading order matching condition

[TABLE]

The $\bar{p}$ component changes in the fast field. On the one hand, this change is given by the difference of $\bar{p}$ values at both ends of the slow fields $I_{s}^{\pm}$, i.e.

[TABLE]

On the other hand, the accumulated jump over the fast field is

[TABLE]

where $\bar{v}$ satisfies $\left(\mathcal{L}^{r}-\lambda\right)\bar{v}=-u_{2}(0;\lambda)v_{p}(\xi)^{2}$. We recall that, in the fast field, to leading order, $u_{p}=u_{0}$ and $v_{p}=\frac{\omega}{u_{0}}$, where $\omega(\xi)=\frac{3}{2}\operatorname{sech}(\xi/2)^{2}$. We rescale $\bar{v}(\xi;\lambda)=-\frac{u_{2}(0;\lambda)}{u_{0}^{2}}V_{\mathrm{in}}(\xi;\lambda)$. Then (171) becomes

[TABLE]

where

[TABLE]

and $V_{\mathrm{in}}$ satisfies

[TABLE]

Equating $\Delta_{\mathrm{s}}\ \bar{p}=\Delta_{\mathrm{f}}\ \bar{p}$ and by (169) one readily derives (at leading order in $\varepsilon$)

[TABLE]

Because of the symmetry $f(x)=f(-x)$, $g(x)=-g(-x)$, it follows that $u_{2}(0;\lambda)=u_{3}(0;\lambda)$ and $p_{2}(0;\lambda)=-p_{3}(0;\lambda)$. Hence

[TABLE]

The inhomogeneous ODE $\left(\mathcal{L}^{r}-\lambda\right)V_{in}=\omega^{2}$ admits bounded solutions for all $\lambda$ that are not eigenvalues of $\mathcal{L}^{r}$. When $\lambda$ is an eigenvalue, though, a bounded solution only exists if the following Fredholm condition is satisfied:

[TABLE]

where $v^{*}$ is the corresponding eigenfunction. Therefore, by Sturm-Liouville theory, it is clear that there is a bounded solution for $\lambda^{r}_{1}=0$, but not for $\lambda_{0}^{r}=5/4$ or $\lambda_{2}^{r}=-3/4$. That is, $\mathcal{R}(\lambda)$, and therefore $t_{22}$, has poles of order $1$ at $\lambda_{0}^{r}$ and $\lambda_{2}^{r}$.

We have, hence, demonstrated the following:

Lemma 9** (Evans function).**

Let the conditions of Theorem 1 and assumption (A4) be fulfilled, and let $(U_{p},V_{p})$ be a pulse solution to (4) as described in Theorem 1. It then holds true that the eigenvalues of the operator $\mathcal{L}$ in (143) arising from linearization around the pulse solution $(U_{p},V_{p})$ coincide on $\mathcal{C}_{e}$ with the roots of the Evans function

[TABLE]

with $\widetilde{\mathcal{D}}(\lambda)\neq 0,\lambda\in\mathcal{C}_{e}$ and where the so-called fast transmission function is given by

[TABLE]

with $\lambda_{1}^{f}=\mathcal{O}(\varepsilon)$, while the so-called slow transmission function is given by

[TABLE]

with some $C_{1},C_{2},\lambda_{0}^{f},\lambda_{2}^{f}\in\mathbb{R}\setminus\{0\}$ and $\widetilde{t}_{22}$ an analytic function on $\mathcal{C}_{e}$. In particular,

[TABLE]

where $p_{2}(0;\lambda)/u_{2}(0;\lambda)$ is the slope of the unstable manifold of the trivial solution to (168) at $x=0$, and $\mathcal{R}$ is given (at leading order in $\varepsilon$) by

[TABLE]

where $V_{\mathrm{in}}$ satisfies $\left(\mathcal{L}^{r}-\lambda\right)V_{\mathrm{in}}=\frac{9}{4}\operatorname{sech}(\xi/2)^{4}$.

Remark 15**.**

The function $\mathcal{R}$ has been extensively studied in [4, Section 3.1.1], [20, Section 4.1] and [19, Section 5]. We would like to stress, however, that $\mathcal{R}$ in this article has a different factor in front of it and is defined in terms of $\lambda$, whereas in [20, 19] it is defined as function of $P:=2\sqrt{1+\lambda}$. A plot of $\mathcal{R}$ has been included in Figure 8.

Remark 16**.**

The eigenvalue problem is often written as a nonlocal eigenvalue problem (NLEP). This can be achieved via the transformation

[TABLE]

which results in the NLEP

[TABLE]

3.1.4 Roots of transmission function $t_{22}$

In the constant coefficient case $f,g\equiv 0$, we have that $p_{2}(0;\lambda)/u_{2}(0;\lambda)=\sqrt{1+m\lambda}$ and so $t_{22}(\lambda)=0$ reduces to

[TABLE]

with $u_{0}$ as in (133), and eigenvalues can be readily extracted from this condition – see [4]; in Figure 9, we show plots of the right-hand side for various $m$. With additional asymptotic approximations, $m\ll 1$ and $m\gg 1$, this can be reduced even further, to leading order to,

[TABLE]

where

[TABLE]

Now, when $\mu\ll 1$, respectively $\nu\ll 1$, the left-hand side of these expressions becomes asymptotically small (since $u_{0}=u_{0}^{-}=\mathcal{O}(1)$, see (133) and Remark 5), but stays positive. Hence solutions $\lambda$ accumulate at points for which $\mathcal{R}(\lambda)-3\approx 0$, which happens to be at the tip of the essential spectrum, i.e. $\lambda=\underline{\lambda}/m\approx-1$, see Figure 9 and [4]. Certainly, no eigenvalues with positive real parts are found.

This idea can be expanded to include the non-autonomous cases. For this, as in the existence problem, we relate the non-autonomous equation to the autonomous equation. Here, it is useful to rescale (168) such that it has the form of (104). Specifically, we set $\tilde{x}=\varepsilon^{2}\mu|\sqrt{1+m\lambda}|\xi$ and $\bar{p}=|\sqrt{1+m\lambda}|\tilde{p}$, under which (168) turns into the system

[TABLE]

The autonomous part of this equation corresponds to the autonomous part for the existence problem – see section 2.4 – and thus possesses an exponential dichotomy with constants $K=1$ and $\rho=1$. Therefore, for a given $\lambda\in\mathcal{C}_{e}$, by roughness (Proposition 4) it follows that the full non-autonomous equation has an exponential dichotomy as well when

[TABLE]

It is easily verified that this condition is satisfied when

[TABLE]

Thus, for all $\lambda\in\mathcal{C}_{e}$, we obtain a (different) bound $\delta_{c}(\lambda)$. Since $\delta_{c}(\lambda)\downarrow 0$ as $|\sqrt{1+m\lambda}|\downarrow 0$ – i.e. when $\lambda$ approaches $-1/m$ – we cannot take the infimum over the region $\mathcal{C}_{e}$. Instead, we further restrict $\lambda$ to $\lambda\in\tilde{C}_{e}:=\mathcal{C}_{e}\cap\left\{\lambda\in\mathbb{C}:|\lambda+\frac{1}{m}|>\frac{1}{2m}\right\}$. Note that $\mathbb{C}^{+}\subset\tilde{C}_{e}$. Then the infimum of $\delta_{c}(\lambda)$ over this region exists, and we define it as $\delta_{c}:=\inf_{\lambda\in\tilde{C}_{e}}\delta_{c}(\lambda)=\frac{\sqrt{6}}{24}\approx 0.102$. Thus, if $\delta<\delta_{c}$, (186) possesses an exponential dichotomy for all $\lambda\in\tilde{C}_{e}$.

Moreover, for all $\lambda\in\tilde{C}_{e}$ and $\delta<\delta_{c}$, the slope $p_{2}(0;\lambda)/u_{2}(0;\lambda)$ of the non-autonomous case can be related to that of the autonomous case, along the same lines as in the existence proof in section 2.4 (specifically, as in Lemma 5). That is, there are $\mathcal{O}(1)$ constants $0<C_{-}(\delta)\leq 1\leq C_{+}(\delta)$ such that $\tilde{p}(0;\lambda)=C\bar{u}(0;\lambda)$ for some $C\in\left(C_{-}(\delta),C_{+}(\delta)\right)$. Rescaling back to the original variables then yields $p_{2}(0;\lambda)/u_{2}(0;\lambda)=C\sqrt{1+m\lambda}$. Therefore $t_{22}(\lambda)=0$ reduces to

[TABLE]

The asymptotic arguments for the autonomous case can now be repeated and it readily follows that no solutions are found with $\lambda\in\tilde{C}_{e}$. In particular $t_{22}(\lambda)=0$ does not have solutions with $\mbox{Re}\lambda>0$. We, hence, have the following result.

Proposition 5** (Roots of the slow transmission function).**

Let $t_{22}$ be the slow transmission function from Lemma 9. Then, for $\lambda\in\left\{\lambda\in\mathcal{C}_{e}:\|\lambda+\frac{1}{m}\|>\frac{1}{2m}\right\}$,

[TABLE]

with $u_{0}=u_{0}^{-}$ as in (133) and for some $C\in\mathbb{R}$ with

[TABLE]

*and $C_{\mathrm{min}/\mathrm{max}}(\delta)$ defined as in Lemma 5.

Moreover, if either of the following two asymptotic approximations hold true,*

• (i)

$m\ll 1$* and $\mu\ll 1$;*

• (ii)

$m\gg 1$* and $\nu\ll 1$,*

then $t_{22}(\lambda)=0$ does not have any solution $\lambda\in\mathcal{C}_{e}$ with $\mbox{Re}\lambda>0$.

Combining Lemma 9 with Proposition 5 readily demonstrates Theorem 2.

3.1.5 Further remarks

If the asymptotic conditions on $m$, $\mu$ and $\nu$ from Proposition 5 do not hold, equation (189) still holds. By restricting $\delta$ further (i.e. taking a lower bound $\delta_{c}$) stronger bounds on the constant $C_{+}$ can be enforced that guarantee all roots of $t_{22}$ lie to the left of the imaginary axis. The proof of this heavily relies on the proof for the autonomous case (see e.g. [4]) and a careful estimation of the constant $C_{+}$. Specifically, the following lemma can be established:

Lemma 10**.**

Let the conditions of Proposition 7 be fulfilled. Then there exists critical values $m_{c}=3$, $0<\mu^{*}(m)<\frac{1}{12}$ (see Theorem 1) and $\nu^{*}(m)>0$ such that if either of the following holds

• (i)

$m<m_{c}$* and $\mu<\mu^{*}(m)$;*

• (ii)

$m>m_{c}$* and $\nu<\nu^{*}(m)$;*

• (iii)

$m=m_{c}$* and $\mu<\mu^{*}(m)$ and $\nu<\nu^{*}(m)$,*

then there exists a $\delta_{c}>0$ such that if $\delta<\delta_{c}$ the condition (189) has no solutions with $\mbox{Re}\lambda>0$; that is, $t_{22}$ has not roots with positive real part.

Remark 17**.**

In (189), the left-hand side is always real-valued. Hence, only $\lambda\in\mathbb{C}$ for which the right-hand side is real-valued can satisfy (189). Due to this, eigenvalues can only appear on a skeleton in $\mathbb{C}$, of which the form only depends on $m$. In Figure 10 we show several skeletons for different $m$. Note that this is the reason for (the shape of) the bounds on the ‘large’ eigenvalues shown in Figure 4 (in red).

Remark 18**.**

The arguments in this section have been applied to pulse solutions with $u_{0}=u_{0}^{-}$ (see (133); $u_{0}^{-}$ as in (90) and (56)). There also exist pulse solutions with $u_{0}=u_{0}^{+}$ (with $u_{0}^{+}$ as in (90) and (56)) and the reasoning also holds for these, up to equation (189). However, $u_{0}^{+}=\mathcal{O}\left(\frac{1}{\mu}\right)$ for these solutions (see Remark 5) and as an effect the left-hand side of (189) thus is asymptotically large (for $\mu\ll 1$). As result, eigenvalues accumulate around the poles of the right-hand side. In particular, because of this, these alternative pulse solution necessarily have an eigenvalue close to $\lambda=5/4>0$, making these pulse solutions unstable.

Remark 19**.**

If $\delta\ll 1$, a direct application of roughness of exponential dichotomies can be used to directly prove that eigenvalues of (154) necessarily lie $\mathcal{O}(\delta)$ close to eigenvalues of the problem with $f\equiv 0$, $g\equiv 0$.

Remark 20**.**

If $\lim_{x\rightarrow\pm\infty}f(x),g(x)$ exist but are not (all) equal to zero, a similar result can be found with minor changes to the proof – provided that the essential spectrum lies to the left of the imaginary axis.

Remark 21**.**

If $\lim_{x\rightarrow\pm\infty}f(x),g(x)$ do not exists, the outlined proof fails because the ‘elephant trunk’ procedure used in the proof of Lemma 7 does no longer work. If $f$ and $g$ approach (possibly different) period functions for $x\rightarrow\infty$ a variant of this proof using a Ricatti transformation such as in [12] seems possible.

3.2 Small eigenvalue close to $\lambda=0$ (Proof of Theorem 3)

In this section we assume that

[TABLE]

which will ease the derivation of a more detailed estimate (as given in Theorem 3) of the location of the small eigenvalue around $\lambda=0$ (in terms of $\delta$), so we set

[TABLE]

The strategy to derive such an estimate is to relate the eigenvalue and existence problems in an appropriate way and then use the Fredholm alternative. To this end, let us write the eigenvalue problem in the fast field (153) in the more concise form

[TABLE]

and the existence problem in the fast field (17) as

[TABLE]

with (the linear part with constant coefficients)

[TABLE]

and

[TABLE]

and $N$ the nonlinear terms. Recall that in the autonomous case the derivative of the pulse solution is an eigenfunction for the zero eigenvalue. Motivated by this, we take a derivative w.r.t. $\xi$ of the non-autonomous existence problem which gives

[TABLE]

and plug into the above eigenvalue problem (200) the ansatz

[TABLE]

which results in

[TABLE]

Upon using (220) to replace the term featuring $\mathbb{L}_{u_{p},v_{p}}(\dot{u}_{p},\dot{v}_{p})^{T}$, we get

[TABLE]

For the perturbation analysis to follow we will use the notation $u_{p,0},v_{p,0},\bar{u}_{0},\bar{v}_{0}$ to indicate the leading order in $\delta$ of the corresponding terms. In particular, $u_{p,0},v_{p,0}$ are the pulse solutions for the homogeneous case $f=g=0$ as described in Corollary 1. We, hence, arrive at the leading order in $\delta$ of the previous equation

[TABLE]

with

[TABLE]

and

[TABLE]

In order to find an expression for the eigenvalue correction $\widetilde{\lambda}$, we will make use of the Fredholm alternative for (258). Hence, we first need to study the kernel of the adjoint operator

[TABLE]

that is, to find $(u^{*},v^{*})^{T}$ with

[TABLE]

and rearrange the solvability condition

[TABLE]

to get an expression for $\tilde{\lambda}$. Since (275) is again a singularly perturbed problem (in $\varepsilon$), we split this problem into three regions: two slow regions, $I_{s}^{\pm}$, and one fast region, $I_{f}$. As described in Theorem 1 and Corollary 1, we have

[TABLE]

where $\omega(\xi)=\frac{3}{2}\operatorname{sech}(\xi/2)^{2}$ and the notation “$p,0,0$” indicates that this the leading order in both, $\delta$ and $\varepsilon$. In the slow regions we have $v_{p,0,0}=0$ to leading order and therefore (again to leading order)

[TABLE]

where $C^{\pm}$ and $D^{\pm}$ are constants that need to be found via matching with the fast field at $\xi=\pm 1/\sqrt{\varepsilon}$. In the fast region, the adjoint problem is to leading order given by

[TABLE]

Up to a multiplicative constant, the only bounded solution to the $v^{*}$-equation is $v^{*}=\frac{1}{u_{0}}\omega^{\prime}$. Matching with the slow fields indicates $D^{\pm}=0$. The expression for $u^{*}$ in $I^{f}$ can be found by integrating twice, which reveals

[TABLE]

The value of $C_{2}$ turns out to be irrelevant and therefore we choose $C_{2}=0$ for simplicity of presentation. Matching with the slow fields then gives $C^{-}=\frac{6}{5u_{0}^{3}}$ and $C^{+}=-\frac{6}{5u_{0}^{3}}$. In summary, we have to leading order in $\varepsilon$

[TABLE]

and

[TABLE]

[TABLE]

We can now assemble the different terms for the solvability condition

[TABLE]

Using that $f$ is odd and $g$ is even, which makes $f^{\prime}$ even and $g^{\prime}$ odd, we get to leading order

[TABLE]

Putting all pieces together, the solvability condition reads

[TABLE]

which can be rearranged to

[TABLE]

Since the problem is solved by a regular perturbation approach, the asymptotic analysis may be validated rigorously by classical methods (i.e. by rigorously controlling the higher order terms); alternatively a geometrical approach based on Lin’s method may be employed (see e.g. [3]).

To show Corollary 4, we observe that in the double asymptotic limit $\tau:=\varepsilon^{4}\mu m\ll 1$ and $\mu\ll 1$, the leading order expression for $\tilde{\lambda}$ becomes

[TABLE]

where we used that $u_{0}=u_{0}^{-}(\mu)\rightarrow 3$ for $\mu\rightarrow 0$ (see Corollary 2 and (57)).

3.2.1 Interpretation of results for ecological applications

Going back to the ecological application, we set $f(x)=h^{\prime}(x)$ and $g(x)=h^{\prime\prime}(x)$. Depending on the rate of topographical variation, several different simplifications can be made to Theorem 3, that allow us to make generic statements about stability of pulse solutions on these terrains.

First, if the topographical changes are small, i.e. when $h=\mathcal{O}(\delta)$, we can write $h(x)=\delta\tilde{h}(x)$ and then (145) can be simplified (via integration by parts):

Corollary 5** (small eigenvalue for height function $h$).**

Let the conditions of Theorem 3 be fulfilled. If $\tilde{f}(x)=\tilde{h}^{\prime}(x)$ and $\tilde{g}(x)=\tilde{h}^{\prime\prime}(x)$, then (145) becomes

[TABLE]

additionally, in the double asymptotic limit $\tau:=\varepsilon^{4}\mu m\ll 1$, $\mu\ll 1$ this further reduces to

[TABLE]

Remark 22**.**

Note that $\tilde{h}$ appears in (289), while it does not appear in the original PDE (4), where only its derivatives appear. Thus, increasing $\tilde{h}$ by an additive constant does not affect the system, and in particular should not affect (289). Since $\int_{0}^{\infty}\left(e^{-x}-4(1-\mu u_{0})e^{-2x}\right)\ dx=-(1-2\mu u_{0})$ the result in (289) is indeed not changed when adding a constant to the height function $\tilde{h}$.

Second, if topographical variation happens only over long spatial scales (i.e. for terrains with weak curvature), we can write $\tilde{h}(x)=\hat{h}(\sigma x)$, where $0<\sigma\ll 1$ to indicate the large-scale spatial variability. Hence, $\tilde{f}(x)=\sigma\hat{h}^{\prime}(\sigma x)=\mathcal{O}(\sigma)$ and $\tilde{g}(x)=\sigma^{2}\hat{h}^{\prime\prime}(\sigma x)=\mathcal{O}(\sigma^{2})$. Because of the difference in size of $\tilde{f}$ and $\tilde{g}$, the sign of $\underline{\lambda}_{0}$ can be related to the sign of $\hat{h}^{\prime\prime}(0)$, i.e. to the local curvature at the location of the pulse.

Corollary 6** (small eigenvalue for terrains with weak curvature).**

Let the conditions of Theorem 3 be fulfilled. If $\tilde{f}(x)=\sigma\hat{h}^{\prime}(\sigma x)$ and $\tilde{g}(x)=\sigma^{2}\hat{h}^{\prime\prime}(\sigma x)$ with $0<\sigma\ll 1$, the leading order expansion of (145) becomes

[TABLE]

additionally, in the double asymptotic limit $\tau:=\varepsilon^{4}\mu m\ll 1$, $\mu\ll 1$, this further reduces to

[TABLE]

Furthermore, it follows that $\mbox{sgn}\ \underline{\lambda}_{0}=\mbox{sgn}\ \hat{h}^{\prime\prime}(0)$, i.e. (vegetation) pulses on hilltops are stable and in valleys are unstable.

Proof.

Since $|\tilde{f}^{\prime}(x)|\gg|\tilde{g}^{\prime}(x)|$ we can neglect the terms with $\tilde{g}^{\prime}(x)$ in (145), thus obtaining

[TABLE]

Substitution of $\tilde{f}^{\prime}(x)=\sigma^{2}\hat{h}^{\prime\prime}(\sigma x)$ and Taylor expanding $\hat{h}^{\prime\prime}$ as $\hat{h}^{\prime\prime}(x)=\hat{h}^{\prime\prime}(0)+\mathcal{O}(\sigma^{3})$ immediately yields (291); the rest of the statement follows straightforwardly. ∎

Third, if topographical variation happens over short spatial scales (i.e. for terrains with strong curvature), we can write $\tilde{h}(x)=\breve{h}\left(x/\sigma\right)$, where $0<\sigma\ll 1$ to indicate the short spatial scales. Hence, $\tilde{f}(x)=\breve{h}^{\prime}\left(x/\sigma\right)/\sigma=\mathcal{O}(1/\sigma)$ and $\tilde{g}(x)=\breve{h}^{\prime\prime}\left(x/\sigma\right)/\sigma^{2}=\mathcal{O}(1/\sigma^{2})$. Again, the sign of $\underline{\lambda}_{0}$ can be related to the sign of $\breve{h}^{\prime\prime}(0)$, though the results are now flipped:

Corollary 7** (small eigenvalue for terrains with strong curvature).**

Let the conditions of Theorem 3 be fulfilled. If $\tilde{f}(x)=\breve{h}^{\prime}\left(x/\sigma\right)/\sigma$ and $\tilde{g}(x)=\breve{h}^{\prime\prime}\left(x/\sigma\right)/\sigma^{2}$ with $0<\sigma\ll 1$ and $\breve{h}(y),\breve{h}^{\prime}(y),\breve{h}^{\prime\prime}(y)\rightarrow 0$ exponentially fast for $|y|\rightarrow\infty$, the leading (and next-leading) order expansion of (145) becomes

[TABLE]

additionally, in the double asymptotic limit $\tau:=\varepsilon^{4}\mu m\ll 1$, $\mu\ll 1$, this further reduces to

[TABLE]

Furthermore, it follows that $\mbox{sgn}\ \underline{\lambda}_{0}=-\mbox{sgn}\ \breve{h}^{\prime\prime}(0)$ when $\mu\neq 0$, i.e. (vegetation) pulses on hilltops are unstable and in valleys are stable; and $\mbox{sgn}\ \underline{\lambda}_{0}=\mbox{sgn}\ \breve{h}(0)$ when $\mu=0$.

Proof.

Substitution of $\tilde{h}(x)=\breve{h}(x/\sigma)$ and the use of the transformation $y=x/\sigma$ in (289) yields

[TABLE]

Taylor expanding the exponential functions then indicates the integral contributes only at order $\mathcal{O}(\delta\tau\sigma)$. Hence the claimed results follow. ∎

Thus, the corollaries in this section indicate that – under certain assumptions on the limiting behavior of the topography function $h$ – vegetation patterns concentrated on hilltops are stable if the terrain has weak curvature and unstable if the terrain has strong curvature; similarly, patterns concentrated in valleys are unstable for terrains with weak curvature, but they become stable if the terrain has strong curvature. A more in-depth inspection of this phenomena can be found in section 4.4, where a few explicit terrain functions $h$ are studied numerically.

4 The effect of the small eigenvalue: movement of pulses

In the previous section we found that, under certain ‘standard’ assumptions on the system’s parameters, all large eigenvalues of a homoclinic pulse solution reside to the left of the imaginary axis. Only one small eigenvalue can lead to destabilization of the pulse solution. Since this small eigenvalue is closely related to the translation invariance of the system without spatially varying coefficients, it is possible to study its effects by projecting the whole system unto the corresponding eigenspace.

This derivation enables us to reduce the full PDE dynamics of (4) to a simpler ODE that describes the movement of the pulse’s location. Concretely, let $P$ denote the location of the center of the pulse. Then the time-evolution of $P$ is given by

[TABLE]

where the superscripts $\pm$ denote taking the upper respectively lower limit, $\tau:=\varepsilon^{4}\mu m=\frac{Da^{2}}{m\sqrt{m}}$ and $\tilde{u}$ solves the differential-algebraic equation

[TABLE]

We follow [4] and only give a short formal derivation of this PDE-to-ODE reduction, in section 4.1. We refrain from going into the details of (proving) the validity of this reduction. Although the renormalization group approach of [6, 18] for semi-strong pulse interactions has not yet been applied to systems with inhomogeneous terms, it can naturally be extended to include these effects. However, it should be noted that, so far, the results and techniques of [6, 18] only cover strongly restricted region in parameter space: the general issue of validity of the reduction of semi-strong pulse interactions to finite dimensional settings still largely remains an open question in the field – see also [4]. As a consequence, we formulate the main results of this section as Propositions and only provide their formal derivations.

Using the pulse location ODE (297) we use formal analysis in section 4.2 to present a scheme by which we can determine the stability of the homoclinic pulse patterns of Theorem 2.5 for any functions $f$ and $g$, i.e. without the restriction on their size by which we obtained Theorem 3; in section 4.3 we (formally) validate this scheme by reducing it to the setting of Theorem 3, i.e. by assuming that $f,g=\mathcal{O}(\delta)$ (with $\delta\ll 1$), and showing that this indeed confirms the results of Theorem 3. Next, we study a few explicit functions in section 4.4 – focusing on what happens when the pulse solution changes stability type. Finally, we briefly consider multi-pulse dynamics in section 4.5.

4.1 Formal derivation of pulse location ODE

In this section we formally derive the pulse location ODE (297). Mathematically, this amounts to tracking perturbations along translational eigenvalues; this approach is sometimes called the ‘collective coordinate method’. Specifically, in this section, we show

Proposition 6**.**

Let $\varepsilon=\frac{a}{m}\ll 1$, $\tau=\frac{Da^{2}}{m\sqrt{m}}\ll 1$ and $\mu=\frac{Dm\sqrt{m}}{a^{2}}\leq\mathcal{O}(1)$ (w.r.t. $\varepsilon$). Let $P$ denote the location of the homoclinic pulse’s center. Then the evolution of $P$ is described by the pulse location ODE (297).

Formal derivation, cf. [4]. We introduce the stretched travelling-wave coordinate

[TABLE]

scale $\frac{dP}{dt}=\frac{Da^{2}}{m\sqrt{m}}c(t)$ and use scalings (14) to transform (4) to get

[TABLE]

To find the solution in the fast region $I_{f}=\left[-1/\sqrt{\varepsilon},1/\sqrt{\varepsilon}\right]$, close to the pulse location, we expand $u$ and $v$ in terms of $\varepsilon$ and look for solution of the form

[TABLE]

To leading order (299) is given by

[TABLE]

Hence we find $u_{0}$ to be constant and

[TABLE]

The next order of (299) is

[TABLE]

It is not a priori clear whether the $v$-equation is solvable; the self-adjoint operator $\mathcal{L}:=\partial_{\xi}^{2}-1+2u_{0}v_{0}$ has a non-empty kernel, since $\mathcal{L}v_{0}^{\prime}=0$, and therefore the inhomogeneous $v$-equation is only solvable when the following Fredholm condition holds

[TABLE]

Upon integrating by parts twice on the right-hand side we obtain

[TABLE]

Since $v_{0}$ is an even function, $u_{1}^{\prime\prime}$ is an even function and $\eta\mapsto\int_{0}^{\eta}v_{0}(y)^{3}dy$ is an odd function. Therefore the last integral vanishes and we obtain

[TABLE]

The integrals over the fast field $I_{f}$ can be approximated by integrals over $\mathbb{R}$, since $v_{0}$ decays exponentially within fast field. Hence we find

[TABLE]

Finally, it follows from the $u$-equation in (303) that

[TABLE]

Combining this with (307) we obtain

[TABLE]

The values of $u_{1}^{\prime}(\pm 1/\sqrt{\varepsilon})$ can be matched to the solutions $\hat{u}$ in the slow fields. Careful inspection of the scalings involved reveals $u_{1}^{\prime}(\pm 1/\sqrt{\varepsilon})=\hat{u}_{x}(P^{\pm})$, where $\hat{u}$ satisfies the differential-algebraic equation (298). Since $\frac{dP}{dt}=\tau c(t)$ this concludes the proof.

Remark 23**.**

Note the link with the notation in section 2: $u_{1}^{\prime}=\hat{p}$. See also Remark 7.

4.2 Stability of fixed points of pulse location ODE (297)

The pulse location ODE (297) describes the movement of a pulse over time. In general, for generic functions $f$ and $g$, it is not possible to solve (298) in closed form, and therefore the pulse location ODE (297) cannot be expressed more explicitly for generic functions $f$ and $g$. Thus, in general, (297) can only be solved numerically – for instance using the numerical scheme developed in [4]. Moreover, for generic $f$ and $g$ fixed points of (297) can only be obtained numerically. However, when $f$ and $g$ obey the symmetry assumptions (A2), one can readily obtain that $P_{*}=0$ is a fixed point. It is possible to determine the stability of fixed points using (297) via direct numerics, but this can be rather time-intensive and is prone to errors close to bifurcation points. Instead, it is better to first use asymptotic expansions to derive a stability condition that can be checked (numerically) more easily.

Proposition 7**.**

Let the conditions of Proposition 6 be satisfied, let $\mu\ll 1$ and let $P_{*}$ be a fixed point of (297). Then, the eigenvalue $\underline{\lambda}$ – where $\underline{\lambda}=m\lambda$, see (152) – corresponding to the pulse solution with a pulse located at the fixed point $P_{*}$ is given by

[TABLE]

Here $\tilde{u}$ and $\tilde{w}$ solve the coupled ODE system

[TABLE]

Remark 24**.**

If $f$ and $g$ satisfy the symmetry assumption (A2) and $P_{*}$ is located at the point of symmetry, i.e. $P_{*}=0$, then symmetry forces $\tilde{u}^{\prime}(P_{*}^{+})=-\tilde{u}^{\prime}(P_{*}^{-})$, $\tilde{u}^{\prime\prime}(P_{*}^{+})=\tilde{u}^{\prime\prime}(P_{*}^{-})$ and $\tilde{w}^{\prime}(P_{*}^{+})=\tilde{w}(P_{*}^{-})$. Therefore, (310) reduces to

[TABLE]

Remark 25**.**

The condition $\mu\ll 1$ in Theorem (7) is not strictly necessary. When this condition holds, the differential-algebraic system (298) simplifies to a normal boundary value problem, since $\tilde{u}(P)=0$ to leading order. However, when $\mu=\mathcal{O}(1)$ (w.r.t. $\varepsilon$) the procedure explained below is still applicable and one can derive a similar result; only this time, $u_{0}$ in (298) needs to be expanded as well and $\tilde{u}$ and $\tilde{w}$ satisfy the coupled differential-algebraic system

[TABLE]

Formal derivation. To find the eigenvalue $\underline{\lambda}$ we need to evaluate the derivative of the right-hand side of (297) at the fixed point $P_{*}$. That is,

[TABLE]

By definition of the derivative

[TABLE]

where $\tilde{u}_{\phi}$ solves (298) with every $P$ replaced by $P+\phi$. For small $\phi$, $\tilde{u}_{\phi}$ can be related to $\tilde{u}$ via a regular expansion. Specifically, let $|\phi|\ll 1$, and expand $\tilde{u}_{\phi}=\tilde{u}+\phi\tilde{w}$. Substitution in (298) and careful bookkeeping readiliy shows that $\tilde{u}$ and $\tilde{w}$ satisfy (311). Finally, upon substituting the expansion for $\tilde{u}_{\phi}$ into (315) and the use of a Taylor expansion we obtain

[TABLE]

Finally, substitution into (314) gives (310).

4.3 Small eigenvalue in case of small spatially varying coefficients

As an example of the use of Proposition 7, in this section we use Proposition 7 to give another proof for Theorem 3 in the limit $\mu\ll 1$. This not only shows the applicability of Proposition 7 but especially the relevance of the pulse location ODE (297). Moreover, it also provides a confirmation of the validity of the formal results in this section.

Alternative formal derivation of Theorem 3 for $\mu\ll 1$. Since $f$ and $g$ satisfy the symmetry assumption (A2), the eigenvalue $\underline{\lambda}$ is given by (312). Therefore, it suffices to only look at the solutions $\tilde{u}$ and $\tilde{w}$ to (311) for $x>0$. Since $f,g=\mathcal{O}(\delta)$ with $\delta\ll 1$, we use regular expansions for $\tilde{u}$ and $\tilde{w}$; that is, we set

[TABLE]

Substitution in (311) gives at leading order

[TABLE]

and at the next order, $\mathcal{O}(\delta)$, we find

[TABLE]

Using the usual techniques to solve these ODEs, one can verify that

[TABLE]

where

[TABLE]

Substitution of these expansions in (312) then yields

[TABLE]

Finally, we note that the eigenvalue has been rescaled as $\underline{\lambda}=m\lambda$ in Theorem 2. Since $\tau/m=\varepsilon^{4}\mu$ and $u_{0}=3$ in the limit $\mu\ll 1$, we have indeed recovered (288), i.e. Theorem 3, in the case $\mu\ll 1$.

4.4 Examples of stationary single-pulse solutions

In this section, we study a few explicit functions $f$ and $g$; in all examples we specify a function $h$ and take $f=h^{\prime}$, $g=h^{\prime\prime}$. Not all functions we consider here limit to [math] as $|x|\rightarrow\infty$; that is, some violate assumption (A4). Therefore, these examples also form an outlook, illustrating how the results in this paper are expected to extend beyond the imposed assumptions on functions $f$ and $g$. Specifically, we consider the following four examples:

• (i)

$h(x)=Ae^{-Bx^{2}}$, ($A\in\mathbb{R}$, $B>0$);

• (ii)

$h(x)=A\operatorname{sech}(Bx)$, ($A\in\mathbb{R}$, $B>0$);

• (iii)

$h(x)=A\cos(Bx)$, ($A\in\mathbb{R}$, $B>0$);

• (iv)

$h(x)=-2\ln(\cosh(\beta x))$, ($\beta>0$).

Note that $\lim_{|x|\rightarrow\infty}f(x),g(x)=0$ in cases (i)–(ii), which therefore satisfy assumption (A4). In case (iii) $f$ and $g$ are periodic when $|x|\gg 1$; in case (iv) $f$ and $g$ do have well-defined (though non-zero) limits for $|x|\rightarrow\infty$.

Remark 26**.**

Note that $A>0$ in (i)–(ii) corresponds to ‘hill-like’ topographies and $A<0$ to ‘valley-like’ topographies. The value of $B$ in (i)–(iii) is a measure of the curvature of the terrain; the higher the value of $B$, the stronger the curvature of the terrain modeled by the function $h$.

Using the pulse location ODE (297) and Proposition 7, we have tracked the fixed points and their stability for these examples in the limit $\mu\ll 1$, using numerical continuation methods. The resulting bifurcation diagrams for (i) are shown in Figure 11(a-b), for (ii) in Figure 12(a-b) and for (iii) in Figure 13(a). In all of these cases, we find fixed points at the point of symmetry, corroborating the results in section 2. For small $B$ values – i.e. for weak curvature topographies – the stability of these fixed points is determined by the sign of $A$: $A>0$ leads to stable and $A<0$ to unstable fixed points – corroborating previous intuition indicating that pulses migrate in uphill direction [40, 37, 4]. However, for sufficiently large values of $B$ –i.e. topographies with strong curvature – the stability of those fixed points changes through a pitchfork bifurcation and new behavior is observed. In case (iii) this even leads to the possibility that both the tops ($BP=0$) as well as the valleys ($BP=\pm\pi$) form stable fixed points of (297). The bifurcation value of the pitchfork bifurcation, $B_{c}(A)$, depends on the value of $A$. Using numerical continuation methods we also tracked this value; the results are in Figures 11(c), 12(c) and 13(b) (for topographies (i), (ii) and (iii)).

Remark 27**.**

Theorem 3, and in particular (288) and (290), provide a leading order analytic expression for $B_{c}(0)$. Evaluating these yields $B_{c}(0)\approx 0.75$ (i), $B_{c}(0)\approx 1.23$ (ii) and $B_{c}(0)=\sqrt{2}$ (iii), which is confirmed by the numerical continuation that indicate $B_{c}(0)\approx 0.75$ (i), $B_{c}(0)\approx 1.24$ (ii) and $B_{c}(0)=1.43$ (iii). Note that $A=0$ is, indeed, just the flat terrain $h(x)\equiv 0$; however, these results for $A=0$ should be interpreted to apply to ‘small’ topographical functions only, where $A$ is asymptotically small.

Moreover, these observations are validated by numerical simulation of the full PDE – see Figure 11(d-g) for (i), Figure 12(d-g) for (ii) and Figure 13(c-f) for (iii). Here, we observe the change in stability of the fixed points and, for well-chosen parameter values, these simulations show convergence to fixed points not located at the point of symmetry. Note also that in the case of periodic topography (i.e. case (iii)), there indeed is a region of $B$-values for which both a pulse at the top of a hill and one at the bottom of a valley can be stable (for the same $B$ value). Thus, we are led to conclude that a pitchfork bifurcation occurs at the critical values $B_{c}(0)$. Simulations indicate that these exist also when the asymptotic limit $\mu\ll 1$ does not hold.

For the last function, (iv), it is possible to derive the pulse location ODE (297) explicitly, since (298) can be solved explicitly – see Corollary 3. Using the expressions given in Corollary 3, a straightforward computation reduces (297) to

[TABLE]

where

[TABLE]

Thus, a point $P_{*}$ is a fixed point if and only if $\mathcal{I}_{1}(P_{*})=\mathcal{I}_{2}(P_{*})$. Straightforward inspection reveals that $P_{*}=0$ therefore is the unique fixed point in case (iv) for all values of $\beta>0$. By Proposition 7 and equation (312) the corresponding (small) eigenvalue $\underline{\lambda}$ can be approximated by

[TABLE]

Upon noting that

[TABLE]

it is clear that $\underline{\lambda}<0$. Hence, $P_{*}=0$ is the only fixed point of (324) in case (iv), which is (globally) stable – for all $\beta>0$. Direct PDE simulations verify this – even when the asymptotic limit $\mu\ll 1$ does not hold – see Figure 14.

4.5 Stationary multi-pulse solutions

The focus in this article has been on single pulse solutions to (4). As a short encore we briefly discuss the possibility of stationary multi-pulse solutions – i.e. solutions with multiple fast excursions. The movement of these solutions can be captured in an ODE much akin to 297. Specifically, let $P_{1},\ldots,P_{N}$ denote the location of $N$ pulses. Then their movement is described by the ODE

[TABLE]

where $\tilde{u}$ satisfies the differential-algebraic system

[TABLE]

The derivation is similar to that of Proposition 7; we omit the details here and refer the interested reader to [4] for a full coverage.

In case of constant coefficients $f,g\equiv 0$, it is well-known that stationary multi-pulse solutions do not exist [13, 4]. In fact, from (328) one can verify that in $2$-pulse solutions the pulses typically move away from each other with a speed proportional to $e^{-\Delta P}$, where $\Delta P:=P_{2}-P_{1}$ is the distance between the pulses – see [13, 4].

However, the non-autonomous terms $f$ and $g$ affect the movement speed and can cancel this repulsive movement. Therefore stationary pulse solutions do exist in (4) for well-chosen $f$ and $g$. In Figure 15 we show several numerical examples of (stable) stationary multi-pulse solutions for various choices of $f$ and $g$.

Remark 28**.**

The spatially varying $f$ and $g$ have a order $\mathcal{O}(f,g)$ effect on the movement speed of the pulses. Finding fixed points of (328) – i.e. finding stationary multi-pulse solutions to (4) – thus boils down to balancing two effects of different size. In particular, if $f,g=\mathcal{O}(\delta)$, only multi-pulse solutions exist with $\Delta P=\mathcal{O}\left(-\ln(\delta)\right)\gg 1$. In this case, existence of stationary multi-pulse solutions can be established rigorously by asymptotic analysis and the methods of geometric singular perturbation theory.

Remark 29**.**

We do not present a full analysis of the spectrum of (evolving) multi-pulse solutions here; they can be stable and unstable depending on the parameter values – similar to the one-pulse variants. A description of how to find the spectrum of multi-pulse solutions can be found in [4].

For generic functions $f$ and $g$ it is, at the moment, not possible to prove existence of stationary multi-pulse solutions (however, see Remark 28 for the case of small $f$, $g$). We do remark however that stationary multi-pulse solutions can be constructed for $f$ and $g$ such that (329) can be solved explicitly, as illustrated by the following proposition.

Proposition 8**.**

Let $h(x)=-2\ln\cosh(\beta x)$, $\beta>0$, $f=h^{\prime}$, $g=h^{\prime\prime}$ and let $\mu\ll 1$. Then there exists a $P_{*}>0$ such that (4) admits a stationary symmetric two-pulse solutions with pulses at $P_{1}=-P_{*}$ and $P_{2}=P_{*}$.

Formal derivation. By symmetry of the desired two-pulse solution, we may set $P_{2}=P$, $P_{1}=-P$. Moreover, necessarily $\tilde{u}^{\prime}(0)=0$. Since $\mu\ll 1$, to leading order we have $\tilde{u}(P)=\tilde{u}(-P)=0$. Therefore $\tilde{u}$ is given to leading order by

[TABLE]

where $\tilde{u}_{\pm}$ and $\tilde{u}_{b}$ are as in Corollary 3. To have stationary pulse solutions, by (328) we need to have

[TABLE]

Upon noting that

[TABLE]

and, since $\lim_{P\rightarrow\infty}\hat{u}_{b}(P)=1$ and $\lim_{P\rightarrow\infty}\hat{u}_{b}^{\prime}(P)=0$,

[TABLE]

continuity of $\mathcal{T}$ guarantees the existence of $P_{*}>0$ as claimed.

Remark 30**.**

This result can be established rigorously by geometric singular perturbation theory, using the methods detailed in section 2. We refrain from giving the details of this procedure.

5 Discussion

In this paper, we studied pulse solutions in a reaction-advection-diffusion system with spatially varying coefficients. The existence of stationary pulse solutions at a point of symmetry was established by combining the usual techniques from geometric singular perturbation theory with the tools from the theory of exponential dichotomies. The latter has been used to generate a saddle-like structure in the slow subsystem, and to obtain bounds on the stable/unstable manifolds of this subsystem. These techniques have also been used to determine the spectral stability of these pulse solutions. None of these concepts or ideas are model-dependent and therefore could be used in a wider variety of models, including Gierer-Meinhardt type models.

Analysis of the spectrum associated to these pulse solutions showed that ‘large’ eigenvalues can be bounded to the stable half-plane, under conditions similar to the usual, constant coefficient case. Although we did not focus on the dynamics of solutions when a large eigenvalue crosses the imaginary axis, simulations show the usual pulse annihilation and pulse splitting phenomena. However, the introduction of spatially varying coefficients does have a significant effect on the so-called ‘small’ eigenvalues (close to $\lambda=0$) because of the break-down of the translation invariance in the system. Therefore, well-chosen $f$ and $g$ can either stabilize or destabilize solutions. When the small eigenvalue is in the unstable half-plane, the pulse solution is unstable and as an effect its position changes. In some cases, this in turn can subsequently lead to a pulse annihilation or a pulse splitting [4]. We expect that a careful tuning of $f$ and $g$ can either prevent or force these subsequent bifurcations, which may have a relevance in the maintenance of vegetation patterns in semi-arid climates.

The small eigenvalues were studied more in-depth in the case of $f=h^{\prime}$, $g=h^{\prime\prime}$ (where $h$ is used to model the topography of a dryland ecosystem). Here, we were able to link the stability of (stationary) pulse solution to the curvature of $h$. If the curvature is weak, the pulse is stable if $h^{\prime\prime}(0)<0$ and unstable if $h^{\prime\prime}(0)>0$; for strong curvature the opposite is true: the pulse is stable if $h^{\prime\prime}(0)>0$ and unstable if $h^{\prime\prime}(0)<0$. We found that this change in stability typically happens via a pitchfork bifurcation, and showed that the associated parameter combinations can be obtained numerically. However, we did not consider a fully general class of functions $f$ and $g$, and we do not know in which way these results generalize to other functions $f$ and $g$ – although for choices $f$ and $g$ for which (4) does not posses the symmetry $(x,u)\rightarrow(-x,u)$ (i.e. when assumption (A2) does not hold), the pitchfork bifurcation will break down. A precise treatment of such generic functions could be the topic of subsequent work.

Moreover, in case of spatially varying coefficients, the system (4) can also posses stationary multi-pulse solutions – i.e. solutions that have multiple fast excursions. When $f,g\equiv 0$, these solutions do not exist. Because the spatially varying coefficients break the translation invariance of the system, these multi-pulse solutions can exist – for well-chosen functions $f$ and $g$. In this article we gave numerical evidence for this and showed their existence for a specific choice of functions. We do not think their existence can be proven in as much generality as the existence of stationary one pulse solutions – certainly, the bounds used in this paper, provided by the theory of exponential dichotomies, are not sufficient in the regions between pulses. For sufficiently small $f$ and $g$, an asymptotic analysis can be developed to overcome this issue, although the distance between subsequent pulses then becomes asymptotically large and asymptotic analysis needs to be done with great care to keep track of the right scalings; this is topic of ongoing research.

Finally, the extended Klausmeier model studied in this paper has its application in ecology, where it is used to model dryland ecosystems. The studied pulse solutions in this model correspond to vegetation ‘patches’ that are typically found in those ecosystems. Naturally, the results in this paper can therefore be used for this application. Specifically, the treatment of a spatially varying height function $h$ is new and is inherently more realistic than taking a constant topography (or a constantly sloped topography) as has been done in the past (see e.g. [40, 5, 27, 3, 3]). Typically, the constant coefficient models exhibit pulses that only move uphill. However, as illustrated with numerics, we have shown that a varying topography can lead to both uphill and downhill movement of pulses. This aligns better with measurements, where also both uphill and downhill movement can be observed – even within the same general region [23, 5]. In this regard, the study in this paper can be seen as a first step to better understand the role of topographic variability in pattern formation.

Acknowledgements

We like to thank Marco Wolters for his exploratory (bachelor) research on the migration of vegetation pulses on periodic topographies. This work was funded by NWO’s Mathematics of Planet Earth program.