Pulse solutions for an extended Klausmeier model with spatially varying coefficients
Robbin Bastiaansen, Martina Chirilus-Bruckner, Arjen Doelman
PPulse solutions for an extended Klausmeier model with spatiallyvarying coefficients ∗ Robbin Bastiaansen † Martina Chirilus-Bruckner † Arjen Doelman † December 20, 2018
Abstract
Motivated by its application in ecology, we consider an extended Klausmeier model, a singularly perturbedreaction-advection-diffusion equation with spatially varying coefficients. We rigorously establish existence ofstationary pulse solutions by blending techniques from geometric singular perturbation theory with boundsderived from the theory of exponential dichotomies. Moreover, the spectral stability of these solutions isdetermined, using similar methods. It is found that, due to the break-down of translation invariance, thepresence of spatially varying terms can stabilize or destabilize a pulse solution. In particular, this leads tothe discovery of a pitchfork bifurcation and existence of stationary multi-pulse solutions.
Since Alan Turing’s revolutionary insight that patterns can emerge spontaneously in systems with multiplespecies if these diffuse at different rates [43], systems of reaction-diffusion equations have served as prototypicalpattern forming models. Scientists have been using these reaction-diffusion models successfully to describe forinstance animal markings [28], embryo development [32] and the faceted eye of
Drosophila [31]. Special interesthas been given to localized solutions (e.g. pulses, fronts), that arise when the diffusivity of species involved isvery different. The prototypical (two-component) model (in one spatial dimensional) is a singularly perturbedequation of the (scaled) form (cid:26) ∂ t U = ∂ x U + H ( x, u, u x , v, v x ; ˜ ε ) ,∂ t V = ˜ ε ∂ x V + H ( x, u, u x , v, v x ; ˜ ε ) , (1)where 0 < ˜ ε (cid:28) H , H 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-Scottmodel [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 focusedon the introduction of localized spatial inhomogeneities [44, 34, 35, 48, 49, 21]; also (often formal) research hasbeen done on reaction-diffusion equations with (less restricted) spatially varying coefficients [9, 8, 2, 47, 46, 7]. Inthis article, we aim to expand the knowledge of such systems, by studying a reaction-diffusion system with fairlygeneric spatially varying coefficients rigorously; motivated by its use in ecology (see Remark 2), we consider thefollowing extended Klausmeier model with spatially varying coefficients [27, 4]: (cid:40) ∂ t U = ∂ x U + f ( x ) ∂ x U + g ( x ) U + a − U − U V ,∂ t V = D ∂ x V − mV + U V , (2)with x ∈ R , t ≥ , U = U ( x, t ) , V = V ( x, t ) ∈ R , parameters D, a, m > f, g ∈ C b ( R ). Certainconditions are imposed on the parameters and functions f and g – these will be explained in section 1.1. Remark 1.
The model (2) can be brought into the form of (1) by a series of scalings – see section 2 and [17]. ∗ Last edited: December 20, 2018. † Mathematical Institute, Leiden University, 2300 RA Leiden, The Netherlands ([email protected], [email protected], [email protected]). a r X i v : . [ m a t h . A P ] D ec . . . . x U ( x ) − − . . V ( x ) (a) h ( x ) = 0 − − . . . . x U ( x ) . . V ( x ) (b) h ( x ) = exp( − x / . . . . U ( x ) − π − π π π . . x V ( x ) (c) h ( x ) = 0 . x ) Figure 1 – Numerical simulation resulting in a stationary pulse solution for (2) with f ( x ) = h (cid:48) ( x ), g ( x ) = h (cid:48)(cid:48) ( x ),where h ( x ) = 0 (a), h ( x ) = exp( − x /
2) (b) and h ( x ) = 0 . x ) (c). U, V components are blue and red respectively,while the orange curve depicts the bounded solution u b to which the U -component converges for | x | → ∞ . Remark 2 (Application of the extended Klausmeier model) . This system of equations is used as a model inecology to describe the dynamics of vegetation ( U ) and water ( V ). The extended Klausmeier model (2) takes intoaccount the amount of rainfall ( a > ) and mortality rate of the vegetation ( m > ) and goes beyond its classicalversion by modeling a smooth, spatially varying terrain h = h ( x ) which then enters (2) as f ( x ) = h (cid:48) ( x ) , g ( x ) = h (cid:48)(cid:48) ( x ) (see [4]). Variants of the Klausmeier model have been studied in various articles ranging from ecologicalstudies [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 ofmore complicated patterns is a thorough understanding of their building blocks, namely, localized solutions. Thepresent paper is motivated by observations – both in numerical simulations and in real ecosystems [4, 5] – of theimpact 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 solutionsto (2). The presence of spatially varying coefficients, however, alters the approach that usually is taken in the caseof constant coefficients models. For one, with spatially constant coefficients, (2) possesses a uniform stationarystate, with V ≡
0, to which pulse solutions converge for x → ±∞ . In the case of spatially varying coefficients,however, typically such uniform stationary state does not exist; instead, a bounded solution ( u, v ) = ( u b , x → ±∞ – see Figure 1. Moreover, standardproofs using geometric singular perturbation theory typically rely on the availability of closed form expressionsfor orbits of subsystems of (2) – see below. These are no longer available in case of generic spatially varyingcoefficients, and only bounds can be found. Indeed, the core contribution of the present work is to overcome thesedifficulties, which we do by blending geometric singular perturbation theory [24] with the theory of exponentialdichotomies [11] in a new way.In this article, we initially follow the ‘standard’ approach of geometric singular perturbation theory. Thatis, we introduce a small parameter ε := am – see assumption (A1) in section 1.1– and construct a stationarypulse solution to (2) in the limit ε = 0, which present itself as a homoclinic orbit in the related stationaryfast-slow ODE system – in case of spatially varying coefficients it is homoclinic to the bounded solution. Forthis construction, the full system is split into a fast subsystem, and a (super)slow subsystem on a so-called slowmanifold M that consists of fixed points of the fast subsystem. We establish fast connections to and from M that take off from submanifold T o ⊂ M and touch down on submanifold T d ⊂ M . On M , we construct stableand unstable submanifolds W s/u ( u b ) ⊂ M that consists of points on M that converge to the bounded solutionfor x → ∞ respectively x → −∞ . Intersections between these unstable/stable manifolds and take-off/touch-down submanifolds (and a symmetry assumption) then establish the existence of pulse solutions to (2). Finally,persistence of these pulse solutions for ε > U ( x, t ) , V ( x, t )) = (˜ u ( x ) , ˜ v ( x )) of (2) fulfill the system of ODEs (cid:40) u xx + f ( x )˜ u x + g ( x )˜ u + a − ˜ u − ˜ u ˜ v , D m ˜ v xx − ˜ v + m ˜ u ˜ v . (3)After a sequence of (re)scalings, it can be seen that the associated fast subsystem is not affected by the spatiallyvarying terms and can be studied using standard methods. However, the slow subsystem, on the slow manifold M , is affected by the spatially varying terms. This subsystem is given (when rescaling ˆ u = a ˜ u ) by (cid:26) ∂ x ˆ u = ˆ p,∂ x ˆ p = − f ( x )ˆ p − g ( x )ˆ u − u. (4)2 igure 2 – Sketches of the bounded solution (blue) and its stable (green) respectively unstable (red) manifolds incase of constant coefficients (left) and varying coefficients (center and right). ˜ p ˜ u (a) Constant coefficient case ˜ p ˜ u (b) Strong enough bounds ˜ p ˜ u (c) Too weak bounds
Figure 3 – Sketches of a crosssection of M that illustrate the heart of the existence proof. In green the takeoffand touchdown curves are shown, the solid blue lines indicate (possible) l s/u (0), the dashed blue lines l s/u (0) for theconstant coefficient case f = 0 , g = 0. The shaded blue area indicates all possible locations of l s/u (0); the shadedred region the possible locations of the bounded solution. The existence proof works when bounds on u b and l s/u (0)are strong enough such that l u (0) necessarily intersects with T o (0) – this happens when all straight lines that startfrom the red region and stay within the blue region intersect the green curves. If bounds are strong enough this isthe case – as illustrated in (b) – but when bounds are too weak this is not the case and existence is not guaranteedby this method – as illustrated in (c). In (a) the situation for the constant coefficient case is shown. For f and g constant, (4) can be solved explicitly and the stable and unstable manifolds W s,u ( u b ) are knownexplicitly. In case of (spatially) varying f and g , typically no closed form solutions are available; however, whenthese varying coefficients are sufficiently small – specifically, when δ := sup x ∈ R (cid:112) f ( x ) + g ( x ) < (so δ canbe O (1) with respect to ε ); see section 2.3 – the dynamics of (4) can be related to the constant coefficient case f, g ≡ f, g ≡ D family of solutions that converge to the (unique) bounded solution to (4) for x → ∞ and a 1 D family of solutions that converge to the bounded solution for x → −∞ . These families of solutions essentiallyform the stable and unstable manifolds W s,u ( u b ). Due to the linear nature of (4), these (un)stable manifoldsare made up of straight lines, i.e. W s,u ( u b ) = ∪ x ∈ R ( x, l s,u ( x )) where l s,u ( x ) describes a straight line in R . Animportant difference now arises between the cases of constant and varying coefficients: when f, g ≡
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 caseof varying coefficients – see Figure 2. The theory of exponential dichotomies enables us to bound the variationof the lines l s,u ( x ); if δ is small enough (i.e. δ < δ c ( a, m, D ), where δ c ≤ / O (1) with respect to ε ), thesebounds are strict enough that a non-empty intersection (0 , l u (0)) ∩ T o is guaranteed – thus establishing existenceof a (symmetric) pulse solution to (2). See Figure 3 for a sketch.Next, the spectral stability of the thus created pulse solutions is studied. Using similar bounds as in theexistence problem, it is shown that eigenvalues are δ -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 boundedto the stable half-plane { λ ∈ C : Re λ < } . For the ‘small’ eigenvalue – located close to the origin – it is moresubtle. In case of f, g ≡ λ Im λ Figure 4 – Sketch of the spectral bounds obtained in this paper. The shaded areas indicate the possible locationsof spectra in the case of varying coefficients. The solid lines and crosses indicate the location of the essential andpoint spectra in the case of constant coefficients: the essential spectrum (orange), the ‘large’ eigenvalues (red) andthe ‘small’ eigenvalue (green). of (2). The introduction of spatially varying coefficients to the system breaks this invariance and as a result thesmall 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 ofthe functions f and g . In particular, when taking f = h (cid:48) , g = h (cid:48)(cid:48) , the location of the small eigenvalue is relatedto the curvature g = h (cid:48)(cid:48) of h : when the curvature is weak, the pulse solution is stable if g (0) = h (cid:48)(cid:48) (0) < g (0) = h (cid:48)(cid:48) (0) >
0; for strong curvature, this is flipped, due to a pitchfork bifurcation.Finally, the break-down of the translation invariance in (2) has another novel effect. In case of constantcoefficients, stationary multi-pulse solutions – solutions with multiple fast excursions – do not exist, due to thepresence 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 theirexistence 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 pulsesolutions to (2); here we first consider the case f, g ≡ f and g . Then, using the theory of exponential dichotomies, both cases are related to each other, resulting in boundsfor the generic case that allow us to prove existence. In section 3 we study the spectral stability of found pulsesolutions, again by relating the generic case to the constant coefficient case of f, g ≡
0. Then, in section 4 weconsider the small eigenvalues more in-depth using formal and numerical techniques, focusing on the possibleoccurrence of bifurcations; we also present stationary multi-pulse solutions. We conclude with a discussion ofthe results in section 5.
We will make several assumptions throughout the manuscript. Some are crucial, while some serve to simplifythe exposition. ( A1 ) : ε := am (cid:28)
1; (5)( A2 ) : f ( − x ) = − f ( x ) , g ( − x ) = g ( x ) , for all x ∈ R ; (6)( A3 ) : sup x ∈ R (cid:112) f ( x ) + g ( x ) <
14 ; (7)( A4 ) : lim x →±∞ f ( x ) , g ( x ) = 0 ; (8)( A5 ) : || f || C b = O (1) , || g || C b = O (1) (cid:16) w.r.t. am (cid:17) (9)Assumption (A1) ensures the presence of a small parameter, necessary to use geometric singular perturbationtheory [37, 4]. (A2) is a symmetry assumption, that ensures (2) possesses a (point) symmetry in x = 0;this technicality significantly simplifies our rigorous proof; pulse solutions can also be found formally and/or4umerically when (A2) does not hold (and we expect that their existence can be established rigorously byextending our methods). Then, assumption (A3) stems from the theory of exponential dichotomies: when thisholds, solutions to (4) for generic f and g can be linked to solutions of (4) with f, g ≡
0; when (A3) does nothold, this link is not provided by the theory of exponential dichotomies. Assumption (A4) is a technicality thatis only needed in the stability section (specifically for the elephant-trunk method to work); for the existencetheorems it is not necessary; in fact, it is suspected that even stability results continue to hold when (A4) isviolated – see also Remarks 20 and 21. Finally, assumption (A5) is needed to pass limits in the treatment of thefast-slow system.
A crucial step for making the stationary ODE (3) amenable to analytic considerations is to find a parameterregime convenient for rigorous perturbation techniques. While there are various choices, we pick a specific onefor clarity, since our focus is on novel phenomena due to the non-autonomous character of the system and notto classify all possible dynamics across parameter regimes.Following [14, 10, 4], we rescale the spatial coordinate (motivated by the diffusivity of the v -component) andthe amplitudes of the unknowns by ξ := √ mD x , ˜ u = m √ mDa u , ˜ v = a √ mD v , (10)to get (cid:40) u ξξ = a m (cid:104) D ma u − Dm √ ma f (cid:16) D √ m ξ (cid:17) u ξ − D ma g (cid:16) D √ m ξ (cid:17) u − D √ m + uv (cid:105) ,v ξξ = v − uv . (11)It is now convenient to introduce 0 < ε := am , < µ := m √ mDa , (12)and write the above ODEs as the first order system of ODEs ˙ u = εp , ˙ p = ε (cid:2) ε µ u − εµf (cid:0) ε µξ (cid:1) p − ε µ g (cid:0) ε µξ (cid:1) u − ε µ + uv (cid:3) , ˙ v = q , ˙ q = v − uv . (13)In order to use geometric singular perturbation theory, we make the customary assumption (A1), that is,0 < ε (cid:28) . (14)and stipulate assumption (A5) so we can pass to limits.In the autonomous case f ≡ g ≡
0, system (13) has a fixed point (1 /µ, , ,
0) and stationary pulsesolutions of (2) correspond to orbits that are homoclinic to (1 /µ, , , f (cid:54) = 0 , g (cid:54) = 0 there is no fixed point, but instead a unique bounded solution ( u b , p b , , u b , p b , ,
0) is established inthe following proposition proven later in section 2.3 (in the proof of Proposition 4).
Proposition 1 (Existence of a bounded solution for (13)) . Let assumptions (A3) and (A4) be fulfilled. Then (13) has a unique bounded solution ( u b , p b , , that satisfies lim ξ ←±∞ ( u b , p b , ,
0) = (1 /µ, , , . (15) Remark 3 (Orbits homoclinic to bounded solutions) . Note that the assumption lim x →±∞ 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 onlybounded without approaching a constant state when | x | → ∞ , the corresponding constructed pulse solution is alsoa homoclinic to the respective bounded solution. An illustration of such a case is given in Figure 1c, where, dueto the periodicity of the coefficients f, g , the bounded background solution is periodic and so is the pulse solutionin its tails. To highlight the novelty of the presented approach, we first briefly explain how the construction is carriedout in the constant coefficient case f = g = 0, to then proceed to the non-autonomous case.5 .1 Stationary pulse solutions for f = 0 and g = 0 The fast system reads ˙ u = εp , ˙ p = ε (cid:2) ε µ u − ε µ + uv (cid:3) , ˙ v = q , ˙ q = v − uv . (16)Note that this system possesses the symmetry ( ξ, u, p, v, q ) → ( − ξ, u, − p, v, − q ). The corresponding slow systemin the slow scaling η = εξ is given by u (cid:48) = p ,p (cid:48) = ε µ u − ε µ + uv ,εv (cid:48) = q ,εq (cid:48) = v − uv . (17)Restricted to the invariant manifold (cid:102) M := { ( u, p, , | u > } (18)it reads (cid:40) u (cid:48) = p ,p (cid:48) = ε µ u − ε µ , (19)which has a saddle structure around the fixed point (cid:16) µ , (cid:17) with stable and unstable eigenspaces given by˜ l u/s := (cid:26) ( u, p ) | p = εµ ( u − µ ) (cid:27) . (20) Remark 4.
Note that this step is much more intricate in the case of varying coefficients f, g where explicitsolutions are possible only for very specific choices of coefficients. Therefore, one must resort to estimationtechniques for the general case. Overcoming this difficulty using exponential dichotomies is the core contributionof the present work.
The reduced fast system has the form ˙ u = 0 , ˙ p = 0 , ˙ v = q , ˙ q = v − uv . (21)A sketch of its planar subsystem ˙ v = q, ˙ q = v − uv can be found in Figures 5a; this planar subsystem is aHamiltonian system with Hamiltonian H ( v, q ; u ) = 12 q − v + 13 uv . (22)Its fixed point ( v, q ) = (0 ,
0) features a saddle structure and a family of homoclinic orbits (cid:40) v (0) hom ( ξ ; u ) = u ω ( ξ ) , ω ( ξ ) := sech ( ξ/ ,q (0) hom ( ξ ; u ) = ˙ v hom ( ξ ; u ) , u ∈ R \{ } , (23)connecting its stable and unstable manifolds. Hence, (21) is a Hamiltonian system with Hamiltonian (cid:101) K ( u, p, v, q ) = H ( v, q ; u ) . (24)The invariant manifold (cid:102) M from (18) is the collection of saddle points ( u, p, , , u > , p ∈ R , for (21) and is,hence, normally hyperbolic. For its stable and unstable manifolds W s/u ( (cid:102) M ) it holds true that dim[ W s/u ( (cid:102) M )] = 36nd, in fact, W s ( (cid:102) M ) and W u ( (cid:102) M ) (partly) coincide, where the intersection is simply given by the family ofhomoclinic orbits. Moreover, we have that (cid:101) K ( u, p, v, q ) | ( u,p,v,q ) ∈ (cid:102) M = 0.For ε >
0, we note that (cid:102) M is still an invariant manifold of the full system (16). It is a standard result ingeometric singular perturbation theory (see, e.g. the classic articles [42, 24, 26] or, more recent, [30]) that, for ε sufficiently small, its stable and unstable manifolds persist as W s/uε ( (cid:102) M ) with dim[ W s/uε ( (cid:102) M )] = 3, but do notnecessarily coincide anymore. In fact, they generically meet in a 2D intersection in R .In order to analyze the persistence of homoclinic orbits we measure the distance of W sε ( (cid:102) M ) and W uε ( (cid:102) M ) inthe hyperplane (cid:101) 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 , p , v max , (cid:101) K and analyze itsdifference during the jump of the orbit through the fast field I f := (cid:18) − √ ε , √ ε (cid:19) , (25)by setting up ∆ I f (cid:101) K = (cid:101) K (1 / √ ε ) − (cid:101) K ( − / √ ε ) = (cid:90) I f ddξ (cid:101) K ( ξ ) dξ = 13 ε (cid:90) I f p ( ξ ) v hom ( ξ ) dξ + h.o.t. (26)where we used that ddξ (cid:101) K = ∂∂u H ( v, q ; u )( dudξ ) + ddξ H ( v, q ; u ) = v ( dudξ ) + 0 = εv p . We may set (using the factthat p is constant to leading order) p ( ξ ) = p (0) + εp (1) ( ξ ) + h.o.t. Therefore, in order to make this differencevanish 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 (cid:102) M is established through the intersection W sε ( (cid:102) M ) ∩ W uε ( (cid:102) M ) ∩ R , the remaining task is to determine possible take-off and touch-down points on (cid:102) M andinvestigate if these intersect the stable and unstable eigenspaces l s/u appropriately to form a homoclinic. To thisend we observe that∆ I f u = u (1 / √ ε ) − u ( − / √ ε ) = (cid:90) I f ddξ u ( ξ ) dξ = ε (cid:90) I f p (1) ( ξ ) dξ = O ( ε / ) , (27)∆ I f p = p (1 / √ ε ) − p ( − / √ ε ) = (cid:90) I f ddξ p ( ξ ) dξ = εu (cid:90) I f v (0) hom ( ξ ) dξ = 6 u ε + h.o.t. , (28)so, to leading order, only the p -variable changes during the fast jump, and therefore, the take-off and touch-downcurves on (cid:102) M are to leading order given by (cid:101) T o/d := (cid:26) ( u, p, , | p = ∓ εu , u > (cid:27) , (29)where we used that, by symmetry, to leading order p ( ± / √ ε ) = p (0) ±
12 ∆ I f p = ε (cid:18) p (1) (0) ± u (cid:19) . (30)Finally, a straightforward computation of the intersection points of these with the stable and unstable eigenspaces l s/u gives two possible homoclinics when µ ≤ , with u ± = 1 ± √ − µ µ (cid:18) for µ ≤ (cid:19) . (31) Remark 5.
When µ (cid:28) , the expression for u ± , (31) , can be expanded in terms of µ ; this yields for u ± thefollowing expansions u − = 3 + 9 µ + O ( µ ) u +0 = µ − − µ + O ( µ ) (32)A conceptual sketch of the dynamics on (cid:102) M , along with an excursion through the fast field, is given inFigure 5b. Moreover, in Figures 6a and 6b, the evolution of a homoclinic solution is projected onto manifold (cid:102) M .7 q u (cid:16) v (0) hom , q (0) hom (cid:17) (a) Sketch of fast reduced system M l s l u (b) Sketch of homoclinic solution
Figure 5 – Sketches of the fast reduced system (21) (a) and the dynamics on the slow manifold M along with, inred, the excursion through the fast field (b). f and g First, we convert the non-autonomous system into an autonomous one by setting s ( ξ ) := D √ m ξ = ε µξ , (33)which gives the extended (autonomous) fast system ˙ s = ε µ , ˙ u = εp , ˙ p = ε (cid:2) ε µ u − εµf ( s ) p − ε µ g ( s ) u − ε µ + uv (cid:3) , ˙ v = q , ˙ q = v − uv . (34)It is important to note that the symmetry assumptions (A2) on f and g translate directly into a symmetry for(34) which is crucial for the construction of a homoclinic. Lemma 1 (Symmetry of (34)) . Let the symmetry assumptions (A2) be fulfilled, that is, let f be an odd functionand g be an even function. Then we have for (34) the symmetry ( s, u, p, v, q ) → ( − s, u, − p, v, − q ) . The slow system corresponding to (34) in the slow variable η = εξ is given by s (cid:48) = εµ ,u (cid:48) = p ,p (cid:48) = ε µ u − εµf ( s ) p − ε µ g ( s ) u − ε µ + uv ,εv (cid:48) = q ,εq (cid:48) = v − uv . (35)It possesses a three-dimensional invariant manifold M := { ( s, u, p, , | u > , s, p ∈ R } ⊂ R , (36)on which it takes the form s (cid:48) = εµ ,u (cid:48) = p ,p (cid:48) = ε µ u − εµf ( s ) p − ε µ g ( s ) u − ε µ . (37)which is an extension of the non-autonomous system (cid:40) u (cid:48) = p ,p (cid:48) = ε µ u − εµf ( εµη ) p − ε µ g ( εµη ) u − ε µ . (38)8 − . l u (0) l s (0) To (0) Td (0) xU x U (a) ( x, U, U x )-plane for h ( x ) = 0 −
50 500 . l u (0) l s (0) x → ±∞ x = 0 x = − x = 1 T o T d U x U (b) ( U, U x )-plane for h ( x ) = 0 − − . l u (0) l s (0) To (0) Td (0) xU x U (c) ( x, U, U x )-plane for h ( x ) = exp( − x / −
50 500 . x → ±∞ x = 0 x = 0 x = − x = − x = 1 x = 1 T o T d l u (0) l s (0) U x U (d) ( U, U x )-plane for h ( x ) = exp( − x / − π − π π π − . . . T o T d l u (0) l s (0) xU x U (e) ( x, U, U x )-plane for h ( x ) = 0 . x ) −
10 100 . . . x = 0 x = 0 x = − x = − x = 1 x = 1 T o T d l u (0) l s (0) U x U (f ) ( U, U x )-plane for h ( x ) = 0 . x ) Figure 6 – Numerical simulations resulting in a stationary pulse solution for (2) with f ( x ) = h (cid:48) ( x ), g ( x ) = h (cid:48)(cid:48) ( x ),where h ( x ) = 0 (a,b), h ( x ) = exp( − x /
2) (c,d) and h ( x ) = 0 . x ) (e,f). Shown are projections to the ( x, U, U x )-plane (a,c,e) and the ( U, U x )-plane (b,d,f) of a stationary pulse solution (blue) and the bounded solution u b to whichthe U -component converges for | x | → ∞ . Parts of the take-off and touch-down curves ( T o/d ) along with stable andunstable manifolds at x = 0 are also sketched in green respectively red. Note that the plots in this figure correspondto the plots in Figure 1.
9t is now convenient to introduce (or, actually, return to) the super-slow variable x = εµη . We set u ( η ) = µ ˆ u ( εµη ) = µ ˆ u ( x ) and return to the second order non-autonomous setting (cid:40) ddx ˆ u = ˆ p , ddx ˆ p = ˆ u − f ( x ) ˆ p − g ( x ) ˆ u − . (39) Lemma 2 (Symmetry of (39)) . Let the symmetry assumptions (A2) be fulfilled, that is, let f be an odd functionand g be an even function. Then we have for (39) the symmetry ( x, ˆ u, ˆ p ) → ( − x, ˆ u, − ˆ p ) . Remark 6.
For conciseness, we note that we have three different scales: fast scale ξ , slow scale η = εξ , super-slow scale x = εµη = ε µξ The construction that we illustrate in this article therefore relies heavily on assumption (A1). The specificdefinition of the small parameter is convenient since the fast reduced system is an ODE which is known to havehomoclinic solutions and the slow system on the critical manifold M is a linear planar system. Remark 7.
Note the difference between p = dudη and ˆ p = d ˆ udx . Hence, p = ε ˆ p . Proposition 2 (Dynamics on M ) . Consider the slow system on M (37) with f, g fulfilling (A3). Then thereexists a unique bounded solution (ˆ u b , ˆ p b ) of (39) and corresponding connected set Γ ⊂ R ∪ {∞} such that thefollowing holds true: For each fixed x ∈ R there exists C s/u ( x ) ∈ Γ and lines l s/u ( x ) := { (ˆ u, ˆ p ) | ˆ p − ˆ u (cid:48) b ( x ) = C s/u ( x )(ˆ u − ˆ u b ( x )) } , (40) such that the solution to the initial value problem (39) with (ˆ u ( x ) , ˆ p ( x )) = (ˆ u , ˆ p ) ∈ l s ( x ) converges to (ˆ u b , ˆ p b ) for x → ∞ , while with (ˆ u ( x ) , ˆ p ( x )) = (ˆ u , ˆ p ) ∈ l u ( x ) it converges to (ˆ u b , ˆ p b ) for x → −∞ . 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 ∈ 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 ˆ u b ( x ) = µu b ( ξ ). Remark 8.
When lim x →±∞ f ( x ) , g ( x ) = 0 (i.e. assumption (A4)), the unique bounded solution (ˆ u b , ˆ p b ) limitsto the fixed point of the autonomous equation. That is, lim x →±∞ (ˆ u b ( x ) , ˆ p b ( x )) = (1 , . (41)This result implies that there are trajectories on M that lead to and away from the bounded solution (ˆ u b , ˆ p b ).Hence, the only remaining construction steps are the analysis of persistence of orbits biasymptotic to M andtheir touch-down/take-off locations. We therefore switch back to the fast system and examine the dynamicsduring the jump of an orbit through the fast field. In order to pass to the reduced fast system, we use theassumption (A5) so, in the limit ε →
0, we get the reduced fast system ˙ s = 0 , ˙ u = 0 , ˙ p = 0 , ˙ v = q , ˙ q = v − uv . (42)Note that in the reduced fast system the non-autonomous character of our problem is not visible. The onlydifference is the added trivial equation ˙ s = 0. As alluded to in the constant coefficient case in section 2.1, theplanar subsystem ˙ v = q, ˙ q = v − uv is known to be Hamiltonian and features a homoclinic to the saddle point( v, q ) = (0 ,
0) which can be specified explicitly (see (23)). As a result, also (42) is Hamiltonian with K ( s, u, p, v, q ) = H ( v, q ; u ) . (43)The invariant manifold M from (36) is the collection of saddle points ( s, u, p, , , u > , s, p ∈ R , for (42) and is,hence, normally hyperbolic. For its stable and unstable manifolds W s/u ( M ) it holds true that dim[ W s/u ( M )] = 4and, in fact, W s ( M ) and W u ( M ) (partly) coincide, where the intersection is simply given by the family ofhomoclinic orbits. Moreover, we have that K ( s, u, p, v, q ) | ( s,u,p,v,q ) ∈M = 0.The analogy with the constant coefficient case continues for ε > M isan invariant manifold of the full system (34) and that its stable and unstable manifolds persist as W s/uε ( M ) withdim[ W s/uε ( M )] = 4, but do not necessarily coincide anymore. In fact, they generically meet in a 3-D intersectionin R . 10 roposition 3 (Persistence of a homoclinic connection) . Let ε be sufficiently small.1. Define the hyperplane R = { ( s, u, p, v, q ) | q = 0 } . Then dim[ W sε ( M ) ∩ W uε ( M ) ∩ R ] = 2 and orbits in thisintersection fulfill p ( ξ ) = εp (1) ( ξ ) + h.o.t. , that is, the leading order constant term p (0) vanishes.2. The take-off and touch-down surfaces on M of orbits in the intersection W sε ( M ) ∩ W uε ( M ) ∩ R are toleading order given by T o/d ( s ) := (cid:26) ( s, u, p, , | p = ∓ εu , u > (cid:27) . (44)
3. For orbits in the intersection W sε ( M ) ∩ W uε ( M ) ∩ R the touch-down curve T d (0) and stable line l s (0) from (40) intersect in at most two points u ± = u b (0) ± (cid:112) u b (0) + 12 / ( µC s (0))2 , (45) where C s (0) is the slope of the stable line l s (0) from (40) and ˆ u b = µu b is the (rescaled) bounded backgroundsolution from Proposition 2. By symmetry, the analogous is true for the take-off curve T o (0) and unstableline l u (0) from (40) . In particular, the thus computed u ± -values coincide by the aforementioned symmetry C u (0) = − C s (0) – see Proposition 2.4. There are two even homoclinic orbits for (3) with u ± > in case u b (0) + 12 / ( µC s (0)) > and u b (0) − (cid:112) u b (0) + 12 / ( µC s (0)) > . Remark 9.
If we set u b (0) = µ and C s (0) = − in (45) , we recover (31) .Proof. Measuring the distance of W sε ( M ) and W uε ( M ) in the hyperplane R can again be accomplished usingthe difference of the Hamiltonian K during the fast the jump of the orbit through the fast field (25). We haveexactly as in the constant coefficient case (26) where (using that p is constant to leading order) we have set p ( ξ ) = p (0) + εp (1) ( ξ ) + h.o.t. , and used that ddξ K = ∂∂s K ( s, u, p, v, q )( dsdξ ) + ∂∂u H ( v, q ; u )( dudξ ) + ddξ H ( v, q ; u ) =0 + v ( dudξ ) + 0 = εv p . In order to make this difference vanish to leading order, we evidently need that p (0) = 0and 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 variablesduring the jump through the fast field:∆ I f s = s (1 / √ ε ) − s ( − (1 / √ ε )) = (cid:90) I f ddξ s ( ξ ) dξ = 2 √ ε ε µ = O ( ε / ) , (46)∆ I f u = u (1 / √ ε ) − u ( − (1 / √ ε )) = (cid:90) I f ddξ u ( ξ ) dξ = ε (cid:90) I f p (1) ( ξ ) dξ = O ( ε / ) , (47)∆ I f p = p (1 / √ ε ) − p ( − (1 / √ ε )) = (cid:90) I f ddξ p ( ξ ) dξ = εu (cid:90) I f v (0) hom ( ξ ) dξ = 6 u ε + h.o.t. , (48)Hence, to leading order, only the p -variable changes during the fast jump, and therefore, the take-off and touch-down curves on M are to leading order given by (44) where we used that, by symmetry, p ( ± / √ ε ) = p (0) ± ∆ I f p . This proves the second statement.Equating (44) and (40) (where we used that p = ε ˆ p – see Remark 7) gives the equality εµC s (0) ( u − u b (0)) = 3 εu ; (49)the solutions of which give the claimed expression (45) in the third statement. Finally, the fourth statementfollows from inspecting (45).Two examples of homoclinic solutions for varying f and g can be found in Figures 6c–6f. In these figures theevolution of a homoclinic solution is projected onto the manifold M , which shows the essence of Proposition 3.Proposition 3 thus establishes existence of homoclinic solutions for (3) under the conditions stated in Propo-sition 3(4). However, in the case of varying coefficients, there typically are no explicit expressions available forthe bounded solution u b (0) and the constant C s (0). To circumvent this, in the next section we derive bounds onthese using the theory of exponential dichotomy, which simultaneously forms the proof of Proposition 2.11 .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 M inexplicit expressions. Instead, our main tools for constructing a saddle-like structure on M are from the theoryof exponential dichotomies. To fix notation and keep the exposition self-contained, we state (following [11]) thedefinition of exponential dichotomies along with a selection of results that we use here. Definition 1 (Exponential Dichotomy) . Consider the planar ODE ddx Y = B ( x ) Y for the unknown Y : R → R and with B : R → R × a matrix-valued function which is continuous on R . Let Φ = Φ( x ) be the associatedcanonical solution operator. This ODE is said to have an exponential dichotomy if there is a projection matrix P and positive constants K and ρ such that (cid:107) Φ( x ) P Φ − (˜ x ) (cid:107) ≤ Ke − ρ ( x − ˜ x ) , x ≥ ˜ x , (cid:107) Φ( x ) ( I − P ) Φ − (˜ x ) (cid:107) ≤ Ke + ρ ( x − ˜ x ) , x ≤ ˜ x . In the next section we will be interested in first order ODEs of the form ddx Y = [ A + A ( x )] Y + F , (50)with x ∈ R , Y : R → R , A ∈ R × , A : R → R × , F ∈ R . In particular, we would like to corroborate knowledgeof the autonomous version (which is often available in terms of explicit solutions) to deduce qualitative resultsfor the full non-autonomous one. For the sake of clarity, we assemble first all auxiliary systems in one place:First, we have the homogeneous, autonomous system ddx Z h = A Z h . (51)Then, there is the homogeneous, non-autonomous system ddx Y h = [ A + A ( x )] Y h . (52)Finally, we have the inhomogeneous, autonomous system ddx Z = A Z + F. (53) Proposition 4 (Roughness and closeness of bounded solutions) . Let K aut , ρ aut > be the exponential dichotomyconstants of the homogeneous, autonomous ODE (51) and Φ aut , P aut the corresponding solution and projectionoperators. If δ := sup x ∈ R ||| A ( x ) ||| < ρ aut K aut , (54) the non-autonomous ODE (52) has an exponential dichotomy for which the following holds true.1. (Roughness) The exponential dichotomy constants of the homogeneous, non-autonomous ODE (52) are K = K aut and ρ = ρ aut − K aut δ , and concerning the solution and projection operators Φ , P of (52) wehave upon defining Q ( x ) := Φ( x ) P Φ − ( x ) , Q aut ( x ) := Φ aut ( x ) P aut Φ aut − ( x ) (55) the estimate ||| Q ( x ) − Q aut ( x ) ||| ≤ K aut δρ aut , x ∈ R . (56)
2. (Closeness of bounded solutions) There exist unique bounded solutions Z b,aut , Y b of the inhomogeneous,autonomous and non-autonomous ODEs (53) and (50) . In particular, they satisfy sup x ∈ R ||| Y b ( x ) − Z b,aut ( x ) ||| ≤ δK aut Kρ aut ρ (cid:107) F (cid:107) . (57)12 roof. The first statement is the persistence of exponential dichotomies, known as “roughness”, and is a stan-dard result (see [11, Ch.4, Prop.1]). Moreover, another standard result from the theory of exponential di-chotomies stipulates that inhomogeneous equations have unique bounded solutions, when the homogeneousequations 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 (cid:48) ( x ) = A W ( x ) + G ( x ) with G ( x ) = A ( x ) Y b ( x ). The unique bounded solution W b of this ODE satisfies theestimate sup x ∈ R (cid:107) W b ( x ) (cid:107) ≤ K aut ρ aut sup x ∈ R (cid:107) G ( x ) (cid:107) ≤ δK aut Kρ aut ρ (cid:107) F (cid:107) , where we used that sup x ∈ R (cid:107) Y b ( x ) (cid:107) ≤ Kρ (cid:107) F (cid:107) . M (Proof of Proposition 2) Let us introduce the more concise notation Y = (cid:0) ˆ u, ddx ˆ u (cid:1) T such that (39) has the form of (50) from the previoussection; that is, ddx Y = [ A + A ( x )] Y + F , (58)with A = (cid:18) (cid:19) , A ( x ) = (cid:18) − g ( x ) − f ( x ) (cid:19) , F = (cid:18) − (cid:19) . (59) Lemma 3 (Exponential Dichotomy Constants and Roughness) . With the notation of Proposition 4, let δ = sup x ∈ R (cid:112) f ( x ) + g ( x ) < . (60) Then we have ρ aut = K aut = 1 , ρ = 1 − δ, K = 5 / and ||| Q ( x ) − Q aut ( x ) ||| ≤ δ , x ∈ R . (61) Proof.
We have the canonical solution operator Φ( x ) = e A x . The eigenvalues of the matrix A are ± v = √ (1 , T , w = √ (1 , − T . Thus the fixed point Y = (0 , T is asaddle. From this it is clear that we can choose P = ww T = 12 (cid:18) − − (cid:19) . With the basis transformation matrix B = ( v | w ) and the diagonal matrix D = diag(1 , −
1) we then get (cid:107) Φ( x ) P Φ − ( s ) (cid:107) = (cid:107) Be Dx B − P Be − Ds B − (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) − − (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) e − ( x − s ) e − ( x − s ) . A similar reasoning – where one can use that I − P = vv T – gives (cid:107) Φ( x )( I − P )Φ − ( s ) (cid:107) = e ( x − s ) . Thus we have the estimate for exponential dichotomies from Definition 1 with ρ aut = 1 and K aut = 1. Theremaining statements can now be read off Proposition 4.The roughness of exponential dichotomies established in Lemma 3 provides a bound on the projection op-erator Q ( x ) of the non-autonomous system. However, this bound cannot be used directly to prove existence ofhomoclinic 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 (58).These can be defined as W s ( Y b ) := (cid:8) ( x, Y s ( x )) | Y s ( x ) = Y b ( x ) + Φ( x ) P Φ − ( x ) r , r ∈ R (cid:9) , (62) W u ( Y b ) := (cid:8) ( x, Y u ( x )) | Y u ( x ) = Y b ( x ) + Φ( x )( Id − P )Φ − ( x ) r , r ∈ R (cid:9) , (63)where Φ , P are the solution and projection operator for (58). For the construction that we have in mind, it isconvenient to notice that W s/u ( Y b ) = (cid:91) x ∈ R ( x, l s/u ( x )) , (64)13ith lines l s ( x ) = (cid:8) Y s ( x ) | Y s ( x ) = Y b ( x ) + Φ( x ) P Φ − ( x ) r , r ∈ R (cid:9) , (65) l u ( x ) = (cid:8) Y u ( x ) | Y u ( x ) = Y b ( x ) + Φ( x )( I − P )Φ − ( x ) r , r ∈ R (cid:9) . (66)While, in general, it is not possible to find explicit expressions for these objects, we can derive estimates for theirlocations. For this we first observe that the line l s can be written equivalently as l s ( x ) = { (ˆ u, ˆ p ) | ˆ p − ˆ u (cid:48) b ( x ) = C ( x )(ˆ u − ˆ u b ( x )) } , (67)where C ( x ) is the slope of the line. Starting from the bound on the projection operator Q ( x ) = Φ( x ) P Φ − ( x )derived in Lemma 3, a bound on the projection lines will be established in Lemma 4, which is then subsequentlyused to find a bound on the slope C ( x ) via the angle θ ( x ) of the line in Lemma 5.In particular, for the case of (39), we thus obtain l s ( x ) = (cid:110) (ˆ u, ˆ p ) | ˆ p − ˆ u (cid:48) b ( x ) = ( − C ( x ))(ˆ u − ˆ u b ( x )) (cid:111) , (68)with ˜ C ( x ) as in Lemma 5 taking into account that the projection operator depends on x , that is, Q = Q ( x ) andso does the angle θ = θ ( x ), which defines C = C ( x ) and, hence, also ˜ C = ˜ C ( x ).The rest of this section consists of the two technical lemmas that ultimately derive a bound for ˜ C . Lemma 4 (Closeness of projection lines) . Let Q and Q aut be the projection matrices with rank as defined inProposition 4(i), i.e. there are unit vectors q and q aut such that Q = qq T and Q aut = q aut q Taut , and (cid:107) Q − Q aut (cid:107) < δ holds true. Then either (cid:107) q − q aut (cid:107) < √ δ or (cid:107) q + q aut (cid:107) < √ δ .Proof. We prove the equivalent statement that from (cid:107) q − q aut (cid:107) ≥ √ δ and (cid:107) q + q aut (cid:107) ≥ √ δ it follows that (cid:107) Q − Q aut (cid:107) ≥ δ . First we observe that( q − q aut )( q T + q Taut )( q + q aut ) = ( qq T − q aut q Taut )( q + q aut ) + ( qq Taut − q aut q T )( q + q aut )= 2( qq T − q aut q Taut )( q + q aut ) = 2( Q − Q aut )( q + q aut ) . (69)Therefore, by assumption (cid:107) Q − Q aut (cid:107) (cid:107) q + q aut (cid:107) ≥ (cid:107) ( Q − Q aut )( q + q aut ) (cid:107) = 12 (cid:107) ( q − q aut )( q T + q Taut )( q + q aut ) (cid:107) = 12 (cid:107) q + q aut (cid:107) (cid:107) q − q aut (cid:107) ≥ δ (cid:107) q + q aut (cid:107) , from which it follows that (cid:107) Q − Q aut (cid:107) ≥ δ .The previous lemma establishes closeness of projection lines of the autonomous and the non-autonomouscase. The thus obtained bounds on norms can be transferred to bounds on the slope C by use of elementarygeometry. Note that transforming the norm bounds in this way leads to singularities when a projection linepasses the vertical axis (which also leads to a seemingly disjoint set of admittable slopes). A visualisation of theresults 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 as defined in Proposi-tion 4(i), i.e. there are unit vectors q and q aut such that Q = qq T and Q aut = q aut q Taut , and (cid:107) Q − Q aut (cid:107) < δ holds true. Furthermore, let θ, θ aut ∈ [ − π, π ) be defined by q =: (cos( θ ) , sin( θ )) , q aut = (cos( θ aut ) , sin( θ aut )) suchthat the slopes of the lines spans by q and q aut are given by C := tan( θ ) , C aut := tan( θ aut ) . (70) Then there exist constants C min/max ( δ, C aut ) defined by C min ( δ, C aut ) := − (1 + C aut ) √ √ δ √ − δ (1 − δ )+2 C aut √ √ δ √ − δ , if δ (cid:54) = (cid:18) C aut √ C aut (cid:19) −∞ , if δ = (cid:18) C aut √ C aut (cid:19) (71) C max ( δ, C aut ) := +(1 + C aut ) √ √ δ √ − δ (1 − δ ) − C aut √ √ δ √ − δ , if δ (cid:54) = (cid:18) − C aut √ C aut (cid:19) ;+ ∞ , if δ = (cid:18) − C aut √ C aut (cid:19) , (72)14 −√
28 18 14 − − δ ˜ C (a) Plots of C min (blue) in (71) and C max (red) in (72) as functions of δ for C aut = − ˆ p ˆ u (b) Bounds on slope for C aut = − δ < −√ . ˆ p ˆ u (c) Bounds on slope for C aut = − δ > −√ . Figure 7 – Visualisation of the results of Lemma 5. In (a) plots of C min (blue) and C max (red) are shown as functionof δ for C aut = −
1, i.e. the set Γ( δ, − C − C aut . In (b) and (c) plots of the possible slopes C are shown for some δ < −√ (b) and δ > −√ (c).The green line indicates the slope value C aut = − such that C − C aut ∈ Γ ( δ, C aut ) , where Γ ( δ, C aut ) := (cid:16) C min ( δ, C aut ) , C max ( δ, C aut ) (cid:17) , if C min ( δ, C aut ) < C max ( δ, C aut ) ; (cid:16) − ∞ , C max ( δ, C aut ) (cid:17) ∪ (cid:16) C min ( δ, C aut ) , + ∞ (cid:17) , if C max ( δ, C aut ) < C min ( δ, C aut ) . (73) In particular, for q aut = √ (1 , − T we have C aut = − and, hence, C = − C , ˜ C ∈ Γ( δ, − . (74) Proof.
For technical reasons we assume that (cid:107) q − q aut (cid:107) ≤ (cid:107) q + q aut (cid:107) ; if this inequality does not hold, we canscale q → − q without changing the projection matrix Q . Then, with∆ θ := θ − θ aut , (75)we have C − C aut = tan( θ ) − tan( θ aut ) = tan(∆ θ + θ aut ) − tan( θ aut ) = (1 + C aut ) (cid:18) tan(∆ θ )1 − C aut tan(∆ θ ) (cid:19) . (76)From (cid:107) Q − Q aut (cid:107) < δ we know by the previous lemma that (cid:107) q − q aut (cid:107) < √ δ and, hence, since q and q aut areunit vectors, we have 0 ≤ − q T q aut ) = (cid:107) q − q aut (cid:107) < δ = ⇒ − δ < q T q aut . (77)Since arccos( z ) is monotonically decreasing, we hence get from | ∆ θ | = arccos( q T q aut ) that − arccos(1 − δ ) < ∆ θ < arccos(1 − δ ) . (78)Furthermore, since tan( z )1 − C aut tan( z ) is monotonically increasing in z , we have the claimed result by usingtan( ± arccos( z )) = ± √ − z z and some simplifications in (76). Here, we first state our main existence results in detail. Their proofs are given in section 2.6.15 heorem 1 (Existence for general f, g ) . Let assumptions (A1), (A2), (A3) and (A5) be satisfied. Then thereis a µ ∗ with < µ ∗ < and corresponding ε ∗ = ε ∗ ( µ ) > , < δ ∗ = δ ∗ ( µ ) < −√ such that the following holdstrue: For any ε, µ, δ with < µ < µ ∗ , < ε < ε ∗ = ε ∗ ( µ ) , δ = sup x ∈ R (cid:112) f ( x ) + g ( x ) < δ ∗ = δ ∗ ( µ ) , (79) the stationary wave ODE (34) has (two) orbits ( s p ( ξ ) , u p ( ξ ) , p p ( ξ ) , v p ( ξ ) , q p ( ξ )) , that are homoclinic to thebounded solution (cid:16) ξ, ˆ u b ( ε µξ ) µ , ε ˆ u (cid:48) b ( ε µξ ) , , (cid:17) , with ( u p ( ξ ) , v p ( ξ )) to leading order given by (cid:34) (ˆ u b ( ε µξ ) − (ˆ u b (0) − µu ) ˆ u − ( ε µξ )) µ (cid:35) χ − s ( ξ ) + (cid:34) u u sech (cid:16) ξ (cid:17) (cid:35) χ f ( ξ ) + (cid:34) (ˆ u b ( ε µξ ) − (ˆ u b (0) − µu ) ˆ u + ( ε µξ )) µ (cid:35) χ + s ( ξ ) (80) with u = u − or u = u +0 from (45) , i.e. u = ˆ u b (0) ± (cid:112) ˆ u b (0) + 12 µ/C s (0)2 µ ; (81)ˆ u b the bounded solution from Proposition 2 and where the indicator functions χ − s ( ξ ) = χ ( −∞ , − / √ ε ) , χ f ( ξ ) = χ ( − / √ ε, / √ ε ) , χ + s ( ξ ) = χ ( / √ ε, ∞ ) (82) distinguishes the behavior of the solution in the fast and super-slow fields. Furthermore, for ˆ u ± we have theestimates | ˆ u ± ( x ) | ≤ Ce − (1 − δ ) | x | , x ≷ , for some C > , the bounded solution u b obeys sup x ∈ R (cid:113) (ˆ u b ( x ) − + ˆ u (cid:48) b ( x ) ≤ δ − δ . Finally, this homoclinic orbit gives rise to a stationary pulse solution (cid:34) U p ( x, t ) V p ( x, t ) (cid:35) = m √ mDa u (cid:16) √ mD x (cid:17) aD √ m v (cid:16) √ mD x (cid:17) (83) for the Klausmeier model (2) that is biasymptotic to the bounded state (cid:16) a ˆ u b (cid:16) √ mD x (cid:17) , (cid:17) . Corollary 1 (Existence for f, g = 0) . Let f, g = 0 , and the conditions from Theorem 1 be fulfilled. Then ˆ u ± ( x ) = e ∓ x , ˆ u b ≡ . Corollary 2 (Existence for small f, g ) . Let the conditions from Theorem 1 be fulfilled and f = δ ˜ f , g = δ ˜ g where ˜ f , ˜ g = O (1) , < δ (cid:28) (i.e. sup x ∈ R (cid:113) ˜ f ( x ) + ˜ g ( x ) = 1 . Then ˆ u + ( x ) = e − x + δ (cid:20) − e x (cid:90) ∞ x ( ˜ f ( z ) − ˜ g ( z )) e − z dz + e − x (cid:18)(cid:90) ∞ ( ˜ f ( z ) − ˜ g ( z )) e − z dz + (cid:90) x ( ˜ f ( z ) − ˜ g ( z )) ds (cid:19)(cid:21) + h.o.t. , ˆ u − ( x ) = e x + δ (cid:20) e − x (cid:90) x −∞ ( ˜ f ( z ) + ˜ g ( z )) e − z ds − e − x (cid:18)(cid:90) −∞ ( ˜ f ( z ) + ˜ g ( z )) e − z dz + (cid:90) x ( ˜ f ( z ) + ˜ g ( z )) dz (cid:19)(cid:21) + h.o.t. , ˆ u b ( x ) = 1 + δ (cid:20) e x (cid:90) ∞ x ˜ g ( z ) e − z dz + e − x (cid:90) x −∞ ˜ g ( z ) e z dz (cid:21) + h.o.t. . Moreover, u as in (81) can be expressed in terms of δ as u = u + δu + h.o.t., (84) where u corresponds to the u -value for the autonomous case, i.e. u is given by (31) . orollary 3 (Existence for h ( x ) = − βx )) . Let h ( x ) = − βx ) , β > , f = h (cid:48) , g = h (cid:48)(cid:48) , and theconditions from Theorem 1 be fulfilled. Then ˆ u ± ( x ) = e ∓ √ β x cosh( βx ) , ˆ u b ( x ) = u − ( x )2 (cid:112) β (cid:90) ∞ x e − √ β z sech( βz ) dz + u + ( x )2 (cid:112) β (cid:90) x −∞ e √ β z sech( βz ) dz . Remark 10.
Pulses solutions as in Corollary 3 exist for any β > without the need of the general assumptionon δ as in Theorem 1; since the flow on M can be solved explicitly for these functions f and g , no condition on δ is needed. Remark 11.
Since the flow on M can be solved explicitly for the functions f and g as in Corollary 3, it is alsopossible to prove existence of symmetric, stationary -pulse solutions (and, in fact, any symmetric, stationary N -pulse solution). Note that normally, for f, g ≡ , these do not exist, since pulses in (2) repel each other [13, 4];this repulsive force can only be overcome by driving forces due to the spatially varying functions f and g . Wecome back to these multi-pulse solutions in section 4.5. The proofs of the existence results in section 2.5 follow from the theory developed in the preceding sections. Theheart 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 δ 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 therigorous 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 Propo-sition 3(4) are satisfied. Since u b (0) = ˆ u b (0) /µ , these hold if and only if the following three bounds hold true:(i) ˆ u b (0) > C s (0) < u b (0) + 12 µ/C s (0) > u b (0) > − δ − δ , (85)and by Lemma 5 we have C s (0) = − C, ˜ C ∈ Γ( δ, − , (86)where Γ is as in (73). Using these, bound (i) is satisfied when δ < and bound (ii) when δ < −√ . Since thebound (iii) holds true when δ = 0 and µ < , continuity of mentioned bounds on ˆ u b (0) and C s (0) guaranteesthe existence of the critical value 0 < δ ∗ ( µ ) < −√ . Proof of Corollary 1.
This follows immediately from solving (39) with f, g ≡
0, and is also carried out in moredetail in section 2.1.
Proof of Corollary 2.
The super-slow system on M in (39) can be solved using a regular expansion in 0 < δ (cid:28) x →∞ ˆ u + ( x ) and lim x →−∞ ˆ u − ( x ) exist, the results follow by a straightforward calculation. Proof of Corollary 3.
One can easily verify that ˆ u ± solve (39), and that lim x →±∞ ˆ u ± ( x ) = 0. The boundedsolution ˆ u b follows from a standard variation of constants method.17 Linear stability analysis
In the previous section, we proved the existence of stationary 1-pulse solutions to (2). In this section we studythe linear stability of these solutions. For ( U p , V p ) a pulse solution from Theorem 1 we define the linear operator L (cid:18) ¯ U ¯ V (cid:19) = (cid:18) ∂ x ¯ U + f ( x ) ∂ x ¯ U + g ( x ) ¯ U − ¯ U − V p ¯ U − U p V p ¯ VD ∂ x ¯ V − m ¯ V + V p ¯ U + 2 U p V p ¯ V . (cid:19) , (87)with L : H ( R ) × H ( R ) ⊂ L ( R ) × L ( R ) → L ( R ) × L ( R ) and its spectrum by Σ( L ), where we distinguishbetween the point spectrum Σ pt ( L ) and the essential spectrum Σ ess ( L ) = Σ( L ) \ Σ pt ( L ) – we denote the elementsof Σ ess ( L ) by λ . As customary, we say that ( U p , V p ) is linearly stable if there is no spectrum in the right halfplane. In order to keep the exposition at reasonable length, we will concentrate here on characterizing parameterregimes where the only instability that can occur is through the (translational) zero eigenvalue which startsmoving 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 (2) as in Theorem 1. Then the essential spectrum of L from (87) is Σ ess ( L ) = ( −∞ , max {− m, − } ] , (88) and, hence, lies in the left half-plane.Proof. The limiting operator of L at x → ±∞ is L ∞ := diag[ ∂ x − , D ∂ x − m ] (note that we thus explicitlyuse assumption (A4)). Therefore, we have that the boundaries of the essential spectrum are λ ( k ) = − ( k + 1), λ ( k ) = − ( D k + m ), k ∈ R , which immediately gives the claimed result.The assumptions on f, g allow (again through the use of exponential dichotomies) the derivation of boundson the location of the point spectrum, which, under the assumption that f, g are chosen ‘small’, can be furtherrefined to track the one small eigenvalue that can possibly lead to bifurcations. The proof of the followingstatements 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 (2) with u = u − as in (81) . Then there exist constants m c , µ ∗ , ν ∗ > such that if either(i) m < m c and µ < µ ∗ or (ii) m > m c and µ √ m < ν ∗ , then there exists a δ c > such that if ≤ δ < δ c precisely one eigenvalue λ is O ( ε ) -close to and all other eigenvalues of 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 theoryof exponential dichotomies.
Remark 12.
Note that Theorem 2 only holds for pulse solutions with u = u − ; pulse solutions with u = u +0 are always unstable. See also Remark 18. Remark 13.
The constants m c , µ ∗ and ν ∗ in Theorem 2 can be computed explicitly (see Lemma 10). Theorem 3 (Small eigenvalue close to λ = 0 for small f , g ) . Assuming that f = δ ˜ f , g = δ ˜ g with < δ (cid:28) , ˜ f , ˜ g = O (1) (i.e. sup x ∈ R (cid:113) ˜ f ( x ) + ˜ g ( x ) = 1 ), there exists a constant τ ∗ > such that if τ := ε µm < τ ∗ thesmall eigenvalue λ close to λ = 0 is located, to leading order, at λ = 2 τ δu − τ (1 − µu ) (cid:90) + ∞ e − x (cid:16) ˜ f (cid:48) ( x )(1 − µu ) + ˜ g (cid:48) ( x )[ e x + µu − (cid:17) dx, (89) where u is as in (81) and Corollary 2.Proof. This statement is derived in section 3.2 by employing a regular expansion in δ . Corollary 4.
Let the conditions of Theorem 3 be fulfilled. Then, in the double asymptotic limit µ (cid:28) and τ := ε µm (cid:28) the leading order expression for λ becomes λ = 23 τ (cid:90) ∞ e − x (cid:16) ˜ f (cid:48) ( x ) + ˜ g (cid:48) ( x )[ e x − (cid:17) dx. (90) Remark 14.
When the term τ = ε µm = a Dm √ m in (89) becomes too large (larger than τ ∗ ), the pulse becomesunstable due to a traveling wave bifurcation/drift instability [10, 14]. .1 Qualitative description of the point spectrum location (Proof of Theorem 2) This section is devoted to finding the point spectrum of the operator L . For that, we use a decompositionmethod for the Evans function, first developed in [1, 16], which is supplemented by the theory of exponentialdichotomies to treat the varying coefficients in (2). As before, the following computations will again heavily relyon the singularly perturbed structure. Therefore, we introduce for the eigenvalue problem ( L − λI )( ¯ U , ¯ V ) T = 0,that is, λ ¯ U = d dx ¯ U + f ( x ) ddx ¯ U + g ( x ) ¯ U − ¯ U − V p ¯ U − U p V P ¯ V , m λ ¯ V = D m d dx ¯ V − ¯ V + m V p ¯ U + m U p V p ¯ V , (91)and the scalings (analogous to (10) and (12)) ξ = D √ m = ε µx , ¯ U = mεµ ¯ u , U p = mεµu p , ¯ V = 1 εµ ¯ v , V p = 1 εµ v p , (92)to get the fast eigenvalue problem (cid:26) ε µ λ ¯ u = ¨¯ u − ε [2 u p v p ¯ v + v p ¯ u ] − ε µ ¯ u + ε µf ( ε µξ ) ˙¯ u + ε µ g ( ε µξ )¯ u , m λ ¯ v = ¨¯ v − ¯ v + [2 u p v p ¯ v + v p ¯ u ] , (93)which suggests (just as in [4, 10, 14]) the introduction of the scaled eigenvalue parameter λ = mλ , (94)so, finally, (cid:26) ε µ mλ ¯ u = ¨¯ u − ε [2 u p v p ¯ v + v p ¯ u ] − ε µ ¯ u + ε µf ( ε µξ ) ˙¯ u + ε µ g ( ε µξ )¯ u ,λ ¯ v = ¨¯ v − ¯ v + [2 u p v p ¯ v + v p ¯ u ] . (95)It is convenient to introduce φ := (cid:0) ¯ u, ˙¯ u/ ( ε µ ) , ¯ v, ˙¯ v (cid:1) and to write the above ODEs as the system of first orderODEs ˙ φ = A ( ξ ; λ, ε, µ, m ) φ, (96)where A ( ξ ; λ, ε, µ, m ) = ε µ v p /µ + ε µ (cid:2) mλ − g ( ε µξ ) (cid:3) − ε µf ( ε µξ ) 2 u p v p /µ
00 0 0 1 − v p λ − u p v p . (97)From the existence analysis in section 2, we have seen that the real line R can be split in one fast region, I f ,near the pulse location and two super slow fields I ± s to both sides of the fast field: I − s := (cid:18) −∞ , − √ ε (cid:19) , I f := (cid:20) − √ ε , √ ε (cid:21) , I + s := (cid:18) √ ε , ∞ (cid:19) . Since we know that v p vanished to leading order in the slow fields, we have in those regions the system matrix A s ( ξ ; λ, ε, µ, m ) := ε µ ε µ (cid:2) mλ − g ( ε µξ ) (cid:3) − ε µf ( ε µξ ) 0 00 0 0 10 0 1 + λ , (98)that is, the dynamics for slow and fast variables are decoupled. Any value λ ∈ C for which this system of ODEshas a non-trivial solution in L ( R ) × L ( R ) corresponds to an eigenvalue λ = mλ of L . A mechanism (that isby now standard) for detecting eigenvalues is the construction of an Evans function, whose roots coincide withthe eigenvalues of L . Although the Evans function can also be extended into the essential spectrum, we do notneed this in the present work and rather restrict λ to C e := C \ { λ ∈ R : λ ≤ max {− , − /m }} = (cid:26) λ = λm : λ / ∈ Σ ess ( L ) (cid:27) , (99)on which the Evans function is analytic. 19 .1.1 Evans function construction By (conditions and results of) Theorem 1 and assumption (A4), we know that the limiting matrix for | ξ | → ∞ is given by A ∞ ( λ, ε, µ, m ) := ε µ ε µ [1 + mλ ] 0 0 00 0 0 10 0 1 + λ . (100)Its eigenvalues Λ , , , and eigenvectors E , , , areΛ , ( λ ) = ±√ λ, Λ , ( λ ) = ± ε µ √ mλE , ( λ ) = (0 , , , Λ , ) T , E , ( λ ) = (cid:0) , ±√ mλ, , (cid:1) T . (101)where Re (Λ ( λ )) < Re (Λ ( λ )) < < Re (Λ ( λ )) < Re (Λ ( λ )) for λ ∈ C e .The system ˙ φ ∞ = A ∞ ( λ, ε, µ, m ) φ ∞ admits exponential dichotomies on C e . Since A ∞ is exponentially closeto A for large | ξ | , the stable and unstable subspaces of ˙ φ = A ( ξ ; λ, ε, µ, m ) φ and ˙ φ ∞ = A ∞ ( λ, ε, µ, m ) φ ∞ aresimilar when | ξ | → ∞ . In particular, for all λ ∈ C e there is a two-dimensional family of solutions, Φ −∞ ( λ ), to˙ φ ∞ = A ∞ ( λ, ε, µ, m ) φ ∞ such that lim ξ →−∞ φ −∞ ( ξ ) = 0 for all φ −∞ ∈ Φ −∞ ( λ ), and a two-dimensional family ofsolutions, Φ + ∞ ( λ ), to ˙ φ ∞ = A ∞ ( λ, ε, µ, m ) φ ∞ such that lim ξ →∞ φ + ∞ ( ξ ) = 0 for all φ + ∞ ∈ Φ + ∞ ( λ ), which impliesthat the system ˙ φ = A ( ξ ; λ, ε, µ, m ) φ also possesses two two-dimensional families of solutions, Φ − ( λ ) and Φ + ( λ )with the same properties.For the system ˙ φ = A ( ξ ; λ, ε, µ, m ) φ , however, it is possible that the intersection Φ + ( λ ) ∩ Φ − ( λ ) is nonempty.The values λ ∈ C e for which this happens correspond to λ = mλ in the point spectrum Σ pt . To find these, weuse a Evans function [1, 16], which is defined as D ( λ ) = det [ φ (0; λ ) , φ (0; λ ) , φ (0; λ ) , φ (0; λ )] , (102)where { φ ( · ; λ ) , φ ( · ; λ ) } spans the space Φ − ( λ ) and { φ ( · ; λ ) , φ ( · ; λ ) } spans the space Φ + ( λ ). For notationalclarity we have suppressed the dependence on the other parameters. Essentially, the Evans function D ( λ )measures the linear independence of the solution functions φ ,..., . Therefore, zeros of D ( λ ) correspond to valuesof λ for which Φ + ( λ ) ∩ Φ − ( λ ) (cid:54) = ∅ , and thus to eigenvalues in the point spectrum [1].In (102) the solutions φ ,..., 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 inthese 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 (2) as in Theorem 1.Then all eigenvalues λ ∈ Σ pt associated to (95) are roots of the Evans function D ( λ ) = t ( λ ) t ( λ )(1 + mλ )(1 + λ ) exp (cid:18)(cid:90) ∞ f ( x ) dx (cid:19) , (103) where t and t are analytic (transmission) functions of λ , defined by lim ξ →∞ φ ( ξ ; λ ) e − Λ ( λ ) ξ = t E ; (104)lim ξ →∞ φ ( ξ ; λ ) e − Λ ( λ ) ξ = t E , (105) where φ is the (unique) solution to (96) for which lim t →−∞ φ ( ξ ; λ ) e − Λ ( λ ) ξ = E ; (106) and φ is the (unique) solution to (96) (if t ( λ ) (cid:54) = 0 ) for which lim t →−∞ φ ( ξ ; λ ) e − Λ ( λ ) ξ = E ; (107)lim t →∞ φ ( ξ ; λ ) e Λ ( λ ) ξ = 0 (108)20 roof. The proof is heavily based on [16, Section 3.2]. Therefore, we present here only an outline of the proofand refer the interested reader to [16] for more details.The heart of the proof is based on choosing φ ,..., in such way that the scale separation of (2) can be exploited.Because A and A ∞ are exponentially close when ξ → −∞ , there is a unique solution φ such that φ closely follows E ( λ ) e Λ ( λ ) ξ as ξ → −∞ . More precisely, we define φ uniquely such that lim ξ →−∞ φ ( ξ ; λ ) e − Λ ( λ ) ξ = E ( λ ).For ξ → ∞ , we do not know the precise form of φ , but we do know that, asymptotically, it is a combinationof the eigenfunctions of the system ˙ φ ∞ = A ∞ φ ∞ . That is, φ ( ξ ; λ ) → t ( λ ) E e Λ ( λ ) ξ + t ( λ ) E e Λ ( λ ) ξ + t ( λ ) E e Λ ( λ ) ξ + t ( λ ) e Λ ( λ ) ξ as ξ → ∞ , where t , . . . , t are analytic transmission functions.Next, φ must be chosen such that { φ ( · , λ ) , φ ( · , λ ) } spans Φ − ( λ ). As this does not determine φ uniquely,we may, additionally, require that φ grows, at most, as E ( λ ) e Λ ( λ ) ξ for ξ → ∞ . More precisely, we define φ uniquely such that lim ξ →−∞ φ ( ξ ; λ ) e − Λ ( λ ) ξ = E and lim ξ → + ∞ φ ( ξ ; λ ) e − Λ ( λ ) ξ = 0 (note that this constructionis based on insight in t – that may not be 0 – that is obtained by the ‘elephant trunk procedure’, see [16, 25] andRemark 20). For ξ → ∞ , φ is then asymptotically given by φ ( ξ ; λ ) → t ( λ ) E ( λ ) e Λ ( λ ) ξ + t ( λ ) E ( λ ) e Λ ( λ ) ξ + t ( λ ) e Λ ( λ ) ξ as ξ → ∞ , where t , t , t are analytical transmission functions.In a similar vein the solutions φ and φ can be defined such that lim ξ →∞ φ ( ξ ; λ ) e − Λ ( λ ) = E ( λ ) andlim ξ →∞ φ ( ξ ; λ ) e − Λ ( λ ) = E ( λ ).Then, using that (cid:80) j =1 Λ j ( λ ) = 0 and by Liouville’s formula, the Evans function (102) can be rewritten: D ( λ ) = lim ξ →∞ det [ φ ( ξ ; λ ) , φ ( ξ ; λ ) , φ ( ξ ; λ ) , φ ( ξ ; λ )] exp (cid:32) − (cid:90) ξ Tr A ( z ) dz (cid:33) = lim ξ →∞ det (cid:104) φ ( ξ ; λ ) e − Λ ( λ ) ξ , φ ( ξ ; λ ) e − Λ ( λ ) ξ , φ ( ξ ; λ ) e − Λ ( λ ) ξ , φ ( ξ ; λ ) e − Λ ( λ ) ξ (cid:105) exp (cid:32) − (cid:90) ξ Tr A ( z ) dz (cid:33) = det [ t ( λ ) E ( λ ) , t ( λ ) E ( λ ) , E ( λ ) , E ( λ )] exp (cid:18)(cid:90) ∞ f ( x ) dx (cid:19) = t ( λ ) t ( λ )(1 + mλ )(1 + λ ) exp (cid:18)(cid:90) ∞ f ( x ) dx (cid:19) . The roots λ ∈ C e of D ( λ ) thus correspond to the roots of t ( λ ) t ( λ ). The next goal, therefore, is to determinethe roots of these transmission functions. t The transmission function t is closely related to the linearization around the pulse in the fast field,( L r − λ ) v = 0 , L r v := ∂ ξ v − [1 − ξ/ ] v. (109)The eigenvalues of L r are well-known to be λ r0 = 5 / λ r1 = 0 and λ r2 = − /
4. By a standard winding numberargument, it follows that roots of t lie O ( ε )-close to these eigenvalues λ r0 , λ r1 and λ r2 . Lemma 8 (Properties of t ) . Let the conditions of Proposition 7 be fulfilled. The roots of t lie O ( ε ) close tothe eigenvalues (counting multiplicity) of L r , i.e. close to λ r0 = 5 / , λ r1 = 0 and λ r2 = − / .Proof. See [16, Lemma 4.1].Although t has a root (with multiplicity 1) close to λ r = 5 /
4, this does not mean that D ( λ ) has a root forthe same value of λ , since – as will be discussed in the next section – the transmission function t has a pole oforder 1 for the same λ , thus preventing it from being an eigenvalue of L – in the literature, this is known as the‘NLEP paradox’.In studies of autonomous systems, the root of t close to λ = 0 is actually located precisely at λ = 0 becauseof the translation invariance of those autonomous systems. However, (2) is non-autonomous and therefore thisreasoning no longer holds and the eigenvalue close to λ r = 0 can have negative or positive real part. As t does not have a pole for this λ – as will be discussed in the next section – the Evans function D ( λ ) has a root forthis value; it thus corresponds to an eigenvalue of L . To our best knowledge, it is, in general, not possible todetermine the precise location of this eigenvalue; in section 3.2 we compute its location using standard regularperturbation techniques when the non-autonomous terms are small.21 .1.3 Slow transmission function t To determine the transmission function t , we focus on the function φ , as defined in Proposition 7. Perconstruction, we know that φ ( ξ ; λ ) → t ( λ ) E ( λ ) e Λ ( λ ) ξ + t ( λ ) E ( λ ) e Λ ( λ ) ξ + t ( λ ) e Λ ( λ ) ξ as ξ → ∞ . As | Λ ( λ ) | (cid:29) | Λ , ( λ ) | for λ ∈ C e , the term e Λ ( λ ) ξ is exponentially small in the slow fields I ± s . Therefore, we have φ ( ξ ; λ ) ≈ t ( λ ) E ( λ ) e Λ ( λ ) ξ + t ( λ ) E ( λ ) e Λ ( λ ) ξ for ξ ∈ I + s sufficiently large. In this way, φ in the slow fieldsis related to the properties of the exponentially asymptotic constant-coefficient system ˙ φ ∞ = A ∞ ( λ, ε, µ, m ) φ ∞ .However, we need to relate φ in the slow fields to the exponentially asymptotic non-autonomous system ˙ φ s = A s ( ξ ; λ, ε, µ, m ) φ s to determine t .In the slow fields the system ˙ φ s = A s ( ξ ; λ, ε, µ, m ) φ s has the dynamics for the (¯ u, ¯ p ) part completely separatedfrom the dynamics of the (¯ v, ¯ q ) part. The (¯ u, ¯ p ) part is governed by the non-autonomous ODE (cid:18) ˙¯ u ˙¯ p (cid:19) = ε µ [ B ( λ ) + B ( ξ )] (cid:18) ¯ u ¯ p (cid:19) , (110)where B ( λ ) = (cid:18) mλ (cid:19) ; B ( ξ ) = (cid:18) − g ( ε µξ ) − f ( ε µξ ) (cid:19) . Here, only the matrix B carries the non-autonomous part of the differential equation and the system with-out B corresponds to the (¯ u, ¯ p ) part of the system ˙ φ ∞ = A ∞ ( λ, ε, µ, m ) φ ∞ , which has spatial eigenvaluesΛ , = ± ε µ √ mλ . When λ ∈ C e this autonomous system admits an exponential dichotomy on R and, there-fore, by roughness the non-autonomous system (110) does so as well, provided that δ = sup x ∈ R (cid:112) f ( x ) + g ( x ) =sup x ∈ R (cid:107) B ( x ) (cid:107) is sufficiently small. Under these conditions, there exist ˜ ψ ( ξ ; λ ) = ( u ( ξ ; λ ) , p ( ξ ; λ ) , , T and ˜ ψ ( ξ ; λ ) = ( u ( ξ ; λ ) , p ( ξ ; λ ) , , T such that ˜ ψ ( ξ ; λ ) → E ( λ ) e Λ ( λ ) ξ and ˜ ψ ( ξ ; λ ) → E ( λ ) e Λ ( λ ) ξ as | ξ | → ∞ . The same reasoning as before can now be used to deduce that φ ( ξ ; λ ) ≈ ˜ ψ ( ξ ; λ ) for ξ ∈ I − s and φ ( ξ ; λ ) ≈ t ( λ ) ˜ ψ ( ξ ; λ ) + t ( λ ) ˜ ψ ( ξ ; λ ) for ξ ∈ I + s .To compute t we need to track the changes of ¯ u and ¯ p during the fast transition when ξ ∈ I f . From (95),it follows that ¯ u stays constant to leading order. Hence, matching φ at the ends of both super-slow fields I ± s gives the leading order matching condition u (0; λ ) = t ( λ ) u (0; λ ) + t ( λ ) u (0; λ ) . (111)The ¯ p component changes in the fast field. On the one hand, this change is given by the difference of ¯ p valuesat both ends of the slow fields I ± s , i.e.∆ s ¯ p = t ( λ ) p (0; λ ) + t ( λ ) p (0; λ ) − p (0; λ ) . (112)On the other hand, the accumulated jump over the fast field is∆ f ¯ p = 1 µ (cid:90) I f (cid:0) v p ( ξ ) u (0; λ ) + 2 u p ( ξ ) v p ( ξ )¯ v ( ξ ; λ ) (cid:1) dξ, (113)where ¯ v satisfies ( L r − λ ) ¯ v = − u (0; λ ) v p ( ξ ) . We recall that, in the fast field, to leading order, u p = u and v p = ωu , where ω ( ξ ) = sech( ξ/ . We rescale ¯ v ( ξ ; λ ) = − u (0; λ ) u V in ( ξ ; λ ). Then (113) becomes∆ f ¯ p = 1 µ u (0; λ ) u (cid:90) I f (cid:0) ω ( ξ ) − ω ( ξ ) V in ( ξ ; λ ) (cid:1) dξ = 1 µ u (0; λ ) u (6 − R ( λ )) + h.o.t. (114)where R ( λ ) := (cid:90) ∞−∞ ω ( ξ ) V in ( ξ ; λ ) dξ (115)and V in satisfies ( L r − λ ) V in ( ξ ; λ ) = ω ( ξ ) . (116)Equating ∆ s ¯ p = ∆ f ¯ p and by (111) one readily derives (at leading order in ε ) t ( λ ) = 1 + 1 µ u − R ( λ ) p (0; λ ) u (0; λ ) − p (0; λ ) u (0; λ ) . (117)22ecause of the symmetry f ( x ) = f ( − x ), g ( x ) = − g ( − x ), it follows that u (0; λ ) = u (0; λ ) and p (0; λ ) = − p (0; λ ). Hence t ( λ ) = 1 + 1 µ u − R ( λ ) p (0; λ ) u (0; λ ) . (118)The inhomogeneous ODE ( L r − λ ) V in = ω admits bounded solutions for all λ that are not eigenvalues of L r .When λ is an eigenvalue, though, a bounded solution only exists if the following Fredholm condition is satisfied: (cid:90) ∞−∞ ω v ∗ dξ = 0 , (119)where v ∗ is the corresponding eigenfunction. Therefore, by Sturm-Liouville theory, it is clear that there is abounded solution for λ r = 0, but not for λ r = 5 / λ r = − /
4. That is, R ( λ ), and therefore t , has poles oforder 1 at λ r and λ 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 (2) as described in Theorem 1. It then holds true that the eigenvalues of the operator L in (87) arising from linearization around the pulse solution ( U p , V p ) coincide on C e with the roots of the Evansfunction D ( λ ) = t ( λ ) t ( λ ) (cid:101) D ( λ ) , (120) with (cid:101) D ( λ ) (cid:54) = 0 , λ ∈ C e and where the so-called fast transmission function is given by t ( λ ) = C (cid:16) λ − λ f (cid:17) (cid:16) λ − λ f (cid:17) (cid:16) λ − λ f (cid:17) , (121) with λ f = O ( ε ) , while the so-called slow transmission function is given by t ( λ ) = C (cid:101) t ( λ ) (cid:16) λ − λ f (cid:17) (cid:16) λ − λ f (cid:17) , (122) with some C , C , λ f , λ f ∈ R \ { } and (cid:101) t an analytic function on C e . In particular, t ( λ ) = 1 + 1 u µ (cid:18) − R ( λ ) p (0; λ ) /u (0; λ ) (cid:19) , (123) where p (0; λ ) /u (0; λ ) is the slope of the unstable manifold of the trivial solution to (110) at x = 0 , and R isgiven (at leading order in ε ) by R ( λ ) = (cid:90) ∞−∞
32 sech( ξ/ V in ( ξ ; λ ) dξ , (124) where V in satisfies ( L r − λ ) V in = sech( ξ/ . Remark 15.
The function R has been extensively studied in [4, Section 3.1.1], [20, Section 4.1] and [19, Section5]. We would like to stress, however, that R in this article has a different factor in front of it and is definedin terms of λ , whereas in [20, 19] it is defined as function of P := 2 √ λ . A plot of R has been included inFigure 8. Remark 16.
The eigenvalue problem is often written as a nonlocal eigenvalue problem (NLEP). This can beachieved via the transformation V in ( ξ ; λ ) = 3 − µu p (0; λ ) u (0; λ ) (cid:82) ∞−∞ ω ( ξ ) f ( ξ ; λ ) dξ z ( ξ ; λ ) , which results in the NLEP ( L r − λ ) z = ω (cid:82) ∞−∞ ωz dξ − µu p (0; λ ) u (0; λ ) . ˆ λ R (ˆ λ ) Figure 8 – A plot of the function R ( λ ). The red lines show the form of R ( λ ) for real-valued λ , whereas the bluelines also show the complex λ for which R ( λ ) is real-valued; the green, dashed lines indicate the poles of the R ( λ ). -1 1 2-2525 Re λ R ( λ ) − √ mλ (a) m = 0 . -1 1 2-2525 Re λ R ( λ ) − √ mλ (b) m = 1 . -1 1 2-2525 Re λ R ( λ ) − √ mλ (c) m = 10 Figure 9 – Plots of the right-hand side of (125) for various m . The red lines indicate the values for real-valued λ ,whereas the blue lines indicate complex λ for which the right-hand side of (125) is real-valued; in green the poles areshown; see [4] for more details. t In the constant coefficient case f, g ≡
0, we have that p (0; λ ) /u (0; λ ) = √ mλ and so t ( λ ) = 0 reduces to µu = R ( λ ) − √ mλ , (125)with u as in (81), and eigenvalues can be readily extracted from this condition – see [4]; in Figure 9, we showplots of the right-hand side for various m . With additional asymptotic approximations, m (cid:28) m (cid:29)
1, thiscan be reduced even further, to leading order to, µu = R ( λ ) − , when m (cid:28) νu = R ( λ ) − √ λ , when m (cid:29)
1; (126)where ν = m Da = µ √ m . (127)Now, when µ (cid:28)
1, respectively ν (cid:28)
1, the left-hand side of these expressions becomes asymptotically small(since u = u − = O (1), see (81) and Remark 5), but stays positive. Hence solutions λ accumulate at points forwhich R ( λ ) − ≈
0, which happens to be at the tip of the essential spectrum, i.e. λ = λ/m ≈ −
1, see Figure 9and [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, werelate the non-autonomous equation to the autonomous equation. Here, it is useful to rescale (110) such that it24as the form of (58). Specifically, we set ˜ x = ε µ |√ mλ | ξ and ¯ p = |√ mλ | ˜ p , under which (110) turns intothe system (cid:18) ¯ u (cid:48) ˜ p (cid:48) (cid:19) = (cid:34)(cid:18) (cid:19) + (cid:32) − g (˜ x/ |√ mλ | ) | mλ | − f (˜ x/ |√ mλ | ) |√ mλ | (cid:33)(cid:35) (cid:18) ¯ u ˜ p (cid:19) . (128)The autonomous part of this equation corresponds to the autonomous part for the existence problem – seesection 2.4 – and thus possesses an exponential dichotomy with constants K = 1 and ρ = 1. Therefore, for agiven λ ∈ C e , by roughness (Proposition 4) it follows that the full non-autonomous equation has an exponentialdichotomy as well when sup x ∈ R |√ mλ | (cid:115) g ( x ) | mλ | + f ( x ) < . (129)It is easily verified that this condition is satisfied when δ = sup x ∈ R (cid:112) f ( x ) + g ( x ) < δ c ( λ ) := 14 |√ mλ | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:114) mλ mλ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (130)Thus, for all λ ∈ C e , we obtain a (different) bound δ c ( λ ). Since δ c ( λ ) ↓ |√ mλ | ↓ λ approaches − /m – we cannot take the infimum over the region C e . Instead, we further restrict λ to λ ∈ ˜ C e := C e ∩ (cid:8) λ ∈ C : | λ + m | > m (cid:9) . Note that C + ⊂ ˜ C e . Then the infimum of δ c ( λ ) over this region exists, and wedefine it as δ c := inf λ ∈ ˜ C e δ c ( λ ) = √ ≈ . δ < δ c , (128) possesses an exponential dichotomy for all λ ∈ ˜ C e .Moreover, for all λ ∈ ˜ C e and δ < δ c , the slope p (0; λ ) /u (0; λ ) of the non-autonomous case can be relatedto that of the autonomous case, along the same lines as in the existence proof in section 2.4 (specifically, asin Lemma 5). That is, there are O (1) constants 0 < C − ( δ ) ≤ ≤ C + ( δ ) such that ˜ p (0; λ ) = C ¯ u (0; λ ) forsome C ∈ ( C − ( δ ) , C + ( δ )). Rescaling back to the original variables then yields p (0; λ ) /u (0; λ ) = C √ mλ .Therefore t ( λ ) = 0 reduces to Cµu = R ( λ ) − √ mλ . (131)The asymptotic arguments for the autonomous case can now be repeated and it readily follows that no solutionsare found with λ ∈ ˜ C e . In particular t ( λ ) = 0 does not have solutions with Re λ >
0. We, hence, have thefollowing result.
Proposition 5 (Roots of the slow transmission function) . Let t be the slow transmission function from Lemma 9.Then, for λ ∈ (cid:8) λ ∈ C e : (cid:107) λ + m (cid:107) > m (cid:9) , t ( λ ) = 1 + 1 u µ (cid:18) − R ( λ ) C √ mλ (cid:19) , (132) with u = u − as in (81) and for some C ∈ R with < C min ( δ ) < C < C max ( δ ) < ∞ (133) and C min / max ( δ ) defined as in Lemma 5.Moreover, if either of the following two asymptotic approximations hold true,(i) m (cid:28) and µ (cid:28) ;(ii) m (cid:29) and ν (cid:28) ,then t ( λ ) = 0 does not have any solution λ ∈ C e with Re λ > . Combining Lemma 9 with Proposition 5 readily demonstrates Theorem 2.
If the asymptotic conditions on m , µ and ν from Proposition 5 do not hold, equation (131) still holds. Byrestricting δ further (i.e. taking a lower bound δ c ) stronger bounds on the constant C + can be enforced thatguarantee all roots of t lie to the left of the imaginary axis. The proof of this heavily relies on the proof forthe autonomous case (see e.g. [4]) and a careful estimation of the constant C + . Specifically, the following lemmacan be established: 25 λ Im λ (a) m = 0 . -1 1 2-1-0.50.51 Re λ Im λ (b) m = 3 -1 1 2-1-0.50.51 Re λ Im λ (c) m = 10 Figure 10 – Plots of skeletons on which λ that satisfy (131) necessarily need to lie. Lemma 10.
Let the conditions of Proposition 7 be fulfilled. Then there exists critical values m c = 3 , <µ ∗ ( m ) < (see Theorem 1) and ν ∗ ( m ) > such that if either of the following holds(i) m < m c and µ < µ ∗ ( m ) ;(ii) m > m c and ν < ν ∗ ( m ) ;(iii) m = m c and µ < µ ∗ ( m ) and ν < ν ∗ ( m ) ,then there exists a δ c > such that if δ < δ c the condition (131) has no solutions with Re λ > ; that is, t hasnot roots with positive real part. Remark 17. In (131) , the left-hand side is always real-valued. Hence, only λ ∈ C for which the right-hand sideis real-valued can satisfy (131) . Due to this, eigenvalues can only appear on a skeleton in C , of which the formonly depends on m . In Figure 10 we show several skeletons for different m . Note that this is the reason for (theshape 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 = u − (see (81) ; u − asin (45) and (31) ). There also exist pulse solutions with u = u +0 (with u +0 as in (45) and (31) ) and the reasoningalso holds for these, up to equation (131) . However, u +0 = O (cid:16) µ (cid:17) for these solutions (see Remark 5) and asan effect the left-hand side of (131) thus is asymptotically large (for µ (cid:28) ). As result, eigenvalues accumulatearound the poles of the right-hand side. In particular, because of this, these alternative pulse solution necessarilyhave an eigenvalue close to λ = 5 / > , making these pulse solutions unstable. Remark 19. If δ (cid:28) , a direct application of roughness of exponential dichotomies can be used to directly provethat eigenvalues of (96) necessarily lie O ( δ ) close to eigenvalues of the problem with f ≡ , g ≡ . Remark 20. If lim x →±∞ f ( x ) , g ( x ) exist but are not (all) equal to zero, a similar result can be found with minorchanges to the proof – provided that the essential spectrum lies to the left of the imaginary axis. Remark 21. If lim x →±∞ f ( x ) , g ( x ) do not exists, the outlined proof fails because the ‘elephant trunk’ procedureused in the proof of Lemma 7 does no longer work. If f and g approach (possibly different) period functions for x → ∞ a variant of this proof using a Ricatti transformation such as in [12] seems possible. λ = 0 (Proof of Theorem 3) In this section we assume that f ( x ) = δ (cid:101) f ( x ) , g ( x ) = δ (cid:101) g ( x ) , < δ (cid:28) , (cid:101) f , (cid:101) g = O (1) , sup x ∈ R (cid:113) (cid:101) f ( x ) + (cid:101) g ( x ) = 1 , (134)26hich will ease the derivation of a more detailed estimate (as given in Theorem 3) of the location of the smalleigenvalue around λ = 0 (in terms of δ ), so we set λ = δ (cid:101) λ . (135)The strategy to derive such an estimate is to relate the eigenvalue and existence problems in an appropriate wayand then use the Fredholm alternative. To this end, let us write the eigenvalue problem in the fast field (95) inthe more concise form δ ˜ λ (cid:18) ε µ m
00 1 (cid:19) (cid:18) ¯ u ¯ v (cid:19) = L u p ,v p (cid:18) ¯ u ¯ v (cid:19) , (136)and the existence problem in the fast field (11) as0 = L h (cid:18) u p v p (cid:19) + δL in ( ξ ) (cid:18) u p v p (cid:19) + N (cid:18) u p v p (cid:19) + (cid:18) a (cid:19) , (137)with (the linear part with constant coefficients) L h = (cid:18) ∂ ξ − ε µ ∂ ξ − (cid:19) , (138)and L in ( ξ ) = (cid:18) ε µ (cid:101) f ( ε µξ ) ∂ ξ + ε µ (cid:101) g ( ε µξ ) 00 0 (cid:19) , (139)and N the nonlinear terms. Recall that in the autonomous case the derivative of the pulse solution is aneigenfunction for the zero eigenvalue. Motivated by this, we take a derivative w.r.t. ξ of the non-autonomousexistence problem which gives0 = [ L h + δL in ( ξ ) + DN ( u p , v p )] (cid:124) (cid:123)(cid:122) (cid:125) = L up,vp (cid:18) ˙ u p ˙ v p (cid:19) + δ (cid:18) ddξ L in ( ξ ) (cid:19) (cid:18) u p v p (cid:19) , (140)and plug into the above eigenvalue problem (136) the ansatz (cid:18) ¯ u ¯ v (cid:19) = (cid:18) ˙ u p ˙ u p (cid:19) + δ (cid:18) (cid:101) u (cid:101) v (cid:19) , (141)which results in δ ˜ λ (cid:18) ε µ m
00 1 (cid:19) (cid:18) ˙ u p ˙ u p (cid:19) + δ ˜ λ (cid:18) ε µ m
00 1 (cid:19) (cid:18) (cid:101) u (cid:101) v (cid:19) = L u p ,v p (cid:18) ˙ u p ˙ u p (cid:19) + δ L u p ,v p (cid:18) (cid:101) u (cid:101) v (cid:19) . (142)Upon using (140) to replace the term featuring L u p ,v p ( ˙ u p , ˙ v p ) T , we get δ ˜ λ (cid:18) ε µ m
00 1 (cid:19) (cid:18) ˙ u p ˙ v p (cid:19) + δ ˜ λ (cid:18) ε µ m
00 1 (cid:19) (cid:18) (cid:101) u (cid:101) v (cid:19) = − δ (cid:18) ddξ L in ( ξ ) (cid:19) (cid:18) u p v p (cid:19) + δ L u p ,v p (cid:18) (cid:101) u (cid:101) v (cid:19) (143)For the perturbation analysis to follow we will use the notation u p, , v p, , ¯ u , ¯ v to indicate the leading order in δ of the corresponding terms. In particular, u p, , v p, are the pulse solutions for the homogeneous case f = g = 0as described in Corollary 1. We, hence, arrive at the leading order in δ of the previous equation L (cid:18) (cid:101) u (cid:101) v (cid:19) = (cid:18) αβ (cid:19) (144)with L := L u p, ,v p, = (cid:18) ∂ ξ − ε µ − ε v p, − ε u p, v p, v p, ∂ ξ − u p, v p, (cid:19) , (145)27nd (cid:18) αβ (cid:19) := ˜ λ (cid:18) ε µ m
00 1 (cid:19) (cid:18) ˙ u p, ˙ v p, (cid:19) + (cid:18) ddξ L in ( ξ ) (cid:19) (cid:18) u p, v p, (cid:19) = (cid:18) ε µ m ˜ λ ˙ u p, + ε µ ˜ f (cid:48) ( ε µξ ) ˙ u p, + ε µ ˜ g (cid:48) ( ε µξ ) u p, ˜ λ ˙ v p, (cid:19) . (146) In order to find an expression for the eigenvalue correction (cid:101) λ , we will make use of the Fredholm alternative for(144). Hence, we first need to study the kernel of the adjoint operator L ∗ = (cid:18) ∂ ξ − ε µ − ε v p, v p, − ε u p, v p, ∂ ξ − u p, v p, (cid:19) , that is, to find ( u ∗ , v ∗ ) T with L ∗ (cid:18) u ∗ v ∗ (cid:19) = 0 , (147)and rearrange the solvability condition (cid:28)(cid:18) u ∗ v ∗ (cid:19) , (cid:18) αβ (cid:19)(cid:29) L × L = 0 , (148)to get an expression for ˜ λ . Since (147) is again a singularly perturbed problem (in ε ), we split this problem intothree regions: two slow regions, I ± s , and one fast region, I f . As described in Theorem 1 and Corollary 1, we have u p, , ( ξ ) = µ (cid:104) − (1 − µu ) e + ε µξ (cid:105) , ξ ∈ I − s ; u , ξ ∈ I f ; µ (cid:104) − (1 − µu ) e − ε µξ (cid:105) , ξ ∈ I + s , v p, , ( ξ ) = , ξ ∈ I − s ; u ω ( ξ ) , ξ ∈ I f ;0 , ξ ∈ I + s , , (149)where ω ( ξ ) = sech( ξ/ and the notation “ p, ,
0” indicates that this the leading order in both, δ and ε . Inthe slow regions we have v p, , = 0 to leading order and therefore (again to leading order) u ∗ ( ξ ) = (cid:40) C − e ε µξ , ξ ∈ I − s ; C + e − ε µξ , ξ ∈ I + s ; v ∗ ( ξ ) = (cid:40) D − e ξ , ξ ∈ I − s ; D + e − ξ , ξ ∈ I + s , (150)where C ± and D ± are constants that need to be found via matching with the fast field at ξ = ± / √ ε . In thefast region, the adjoint problem is to leading order given by (cid:26) u ∗ + u ω v ∗ , v ∗ − v ∗ + 2 ωv ∗ . Up to a multiplicative constant, the only bounded solution to the v ∗ -equation is v ∗ = u ω (cid:48) . Matching with theslow fields indicates D ± = 0. The expression for u ∗ in I f can be found by integrating twice, which reveals u ∗ ( ξ ) = − u (cid:90) ξ ω ( z ) dz + C = − u
920 [6 cosh( ξ ) + cosh(2 ξ ) + 8] tanh( ξ/
2) sech( ξ/ + C =: σ ( ξ ) . The value of C turns out to be irrelevant and therefore we choose C = 0 for simplicity of presentation. Matchingwith the slow fields then gives C − = u and C + = − u . In summary, we have to leading order in εu ∗ ( ξ ) = + u e + ε µξ , ξ ∈ I − s ; σ ( ξ ) , ξ ∈ I f ; − u e − ε µξ , ξ ∈ I + s , v ∗ ( ξ ) = , ξ ∈ I − s ; u ω (cid:48) ( ξ ) , ξ ∈ I f ;0 , ξ ∈ I + s , , (151)and α ( ξ ) = ε µ e + ε µξ (cid:104) − m (cid:101) λ (1 − µu ) − (cid:101) f (cid:48) ( ε µξ )(1 − µu ) + (cid:101) g (cid:48) ( ε µξ ) (cid:16) e − ε µξ + µu − (cid:17)(cid:105) , ξ ∈ I − s ; ε µ ˜ g (cid:48) ( ε µξ ) u , ξ ∈ I f ; ε µ e − ε µξ (cid:104) m (cid:101) λ (1 − µu ) + (cid:101) f (cid:48) ( ε µξ )(1 − µu ) + (cid:101) g (cid:48) ( ε µξ ) (cid:16) e + ε µξ + µu − (cid:17)(cid:105) , ξ ∈ I + s , (152)28 ( ξ ) = , ξ ∈ I − s ; (cid:101) λu ω (cid:48) ( ξ ) , ξ ∈ I f ;0 , ξ ∈ I + s ; . , (153)We can now assemble the different terms for the solvability condition (cid:28)(cid:18) u ∗ v ∗ (cid:19) , (cid:18) αβ (cid:19)(cid:29) L × L = (cid:90) I − s ∪ I f ∪ I + s u ∗ ( ξ ) α ( ξ ) dξ + (cid:90) I − s ∪ I f ∪ I + s v ∗ ( ξ ) β ( ξ ) dξ (154)Using that f is odd and g is even, which makes f (cid:48) even and g (cid:48) odd, we get to leading order (cid:90) I − s u ∗ ( ξ ) α ( ξ ) dξ = + ε µ (cid:18) u (cid:19) (cid:90) I − s e +2 ε µξ (cid:16) − m ˜ λ (1 − µu ) − − ˜ f (cid:48) ( ε µξ )(1 − µu ) − +˜ g (cid:48) ( ε µξ )[ e − ε µξ + µu − (cid:17) dξ = + ε µ (cid:18) u (cid:19) (cid:90) + ∞ e − x (cid:16) − m ˜ λ (1 − µu ) − ˜ f (cid:48) ( x )(1 − µu ) − ˜ g (cid:48) ( x )[ e x + µu − (cid:17) dx + h.o.t. = − ε µ (cid:18) u (cid:19) (cid:90) + ∞ e − x (cid:16) m ˜ λ (1 − µu ) + ˜ f (cid:48) ( x )(1 − µu ) + ˜ g (cid:48) ( x )[ e x + µu − (cid:17) dx + h.o.t. = − ε µ (cid:18) u (cid:19) (cid:18) m (1 − µu )˜ λ + (cid:90) + ∞ e − x (cid:16) ˜ f (cid:48) ( x )(1 − µu ) + ˜ g (cid:48) ( x )[ e x + µu − (cid:17) dx (cid:19) + h.o.t. (cid:90) I + s u ∗ ( ξ ) α ( ξ ) dξ = − ε µ (cid:18) u (cid:19) (cid:90) I + s e − ε µξ (cid:16) m ˜ λ (1 − µu ) + ˜ f (cid:48) ( ε µξ )(1 − µu ) + ˜ g (cid:48) ( ε µξ )[ e + ε µξ + µu − (cid:17) dξ = − ε µ (cid:18) u (cid:19) (cid:18) m (1 − µu )˜ λ + (cid:90) + ∞ e − x (cid:16) ˜ f (cid:48) ( x )(1 − µu ) + ˜ g (cid:48) ( x )[ e x + µu − (cid:17) dx (cid:19) + h.o.t. (cid:90) I f u ∗ ( ξ ) α ( ξ ) dξ = (cid:90) I f ε µ ˜ g (cid:48) ( ε µξ ) u dξ = O ( ε − / µ ) (cid:90) I ± s v ∗ ( ξ ) β ( ξ ) dξ = h.o.t (cid:90) I f v ∗ ( ξ ) β ( ξ ) dξ = (cid:90) I f ˜ λ u ω (cid:48) ( ξ ) dξ = ˜ λu (cid:18) u (cid:19) + h.o.t. . Putting all pieces together, the solvability condition reads (cid:28)(cid:18) u ∗ v ∗ (cid:19) , (cid:18) αβ (cid:19)(cid:29) L × L = (cid:18) u (cid:19) (cid:104) ˜ λu − ε µ (cid:16) m ˜ λ (1 − µu )+2 (cid:90) + ∞ e − x (cid:16) ˜ f (cid:48) ( x )(1 − µu ) + ˜ g (cid:48) ( x )[ e x + µu − (cid:17) dx (cid:19)(cid:21) + h.o.t. = 0 , which can be rearranged to˜ λ = 2 ε µu − ε µm (1 − µu ) (cid:90) + ∞ e − x (cid:16) ˜ f (cid:48) ( x )(1 − µu ) + ˜ g (cid:48) ( x )[ e x + µu − (cid:17) dx + h.o.t. (155)Since the problem is solved by a regular perturbation approach, the asymptotic analysis may be validatedrigorously by classical methods (i.e. by rigorously controlling the higher order terms); alternatively a geometricalapproach based on Lin’s method may be employed (see e.g. [3]).To show Corollary 4, we observe that in the double asymptotic limit τ := ε µm (cid:28) µ (cid:28)
1, the leadingorder expression for ˜ λ becomes˜ λ = 2 ε µ (cid:90) ∞ e − x (cid:16) ˜ f (cid:48) ( x ) + ˜ g (cid:48) ( x )[ e x − (cid:17) dx + h.o.t. (156)where we used that u = u − ( µ ) → µ → Going back to the ecological application, we set f ( x ) = h (cid:48) ( x ) and g ( x ) = h (cid:48)(cid:48) ( x ). Depending on the rate oftopographical variation, several different simplifications can be made to Theorem 3, that allow us to makegeneric statements about stability of pulse solutions on these terrains.First, if the topographical changes are small, i.e. when h = O ( δ ), we can write h ( x ) = δ ˜ h ( x ) and then (89)can be simplified (via integration by parts): 29 orollary 5 (small eigenvalue for height function h ) . Let the conditions of Theorem 3 be fulfilled. If ˜ f ( x ) =˜ h (cid:48) ( x ) and ˜ g ( x ) = ˜ h (cid:48)(cid:48) ( x ) , then (89) becomes λ = 2 δτu − τ (1 − µu ) (cid:20) − µu ˜ h (cid:48)(cid:48) (0) + ˜ h (0)(1 − µu ) + (cid:90) ∞ ˜ h ( x ) (cid:0) e − x − − µu ) e − x (cid:1) dx (cid:21) ; (157) additionally, in the double asymptotic limit τ := ε µm (cid:28) , µ (cid:28) this further reduces to λ = 23 δτ (cid:20) ˜ h (0) + (cid:90) ∞ ˜ h ( x ) (cid:0) e − x − e − x (cid:1) dx. (cid:21) + h.o.t. (158) Remark 22.
Note that ˜ h appears in (157) , while it does not appear in the original PDE (2) , where only itsderivatives appear. Thus, increasing ˜ h by an additive constant does not affect the system, and in particular shouldnot affect (157) . Since (cid:82) ∞ (cid:0) e − x − − µu ) e − x (cid:1) dx = − (1 − µu ) the result in (157) is indeed not changedwhen adding a constant to the height function ˜ h . Second, if topographical variation happens only over long spatial scales (i.e. for terrains with weak curvature),we can write ˜ h ( x ) = ˆ h ( σx ), where 0 < σ (cid:28) f ( x ) = σ ˆ h (cid:48) ( σx ) = O ( σ ) and ˜ g ( x ) = σ ˆ h (cid:48)(cid:48) ( σx ) = O ( σ ). Because of the difference in size of ˜ f and ˜ g , the sign of λ canbe related to the sign of ˆ h (cid:48)(cid:48) (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 ˜ f ( x ) = σ ˆ h (cid:48) ( σx ) and ˜ g ( x ) = σ ˆ h (cid:48)(cid:48) ( σx ) with < σ (cid:28) , the leading order expansion of (89) becomes λ = τ δσ (1 − µu ) u − τ (1 − µu ) ˆ h (cid:48)(cid:48) (0); (159) additionally, in the double asymptotic limit τ := ε µm (cid:28) , µ (cid:28) , this further reduces to λ = 13 τ δσ ˆ h (cid:48)(cid:48) (0) + h.o.t. (160) Furthermore, it follows that sgn λ = sgn ˆ h (cid:48)(cid:48) (0) , i.e. (vegetation) pulses on hilltops are stable and in valleys areunstable.Proof. Since | ˜ f (cid:48) ( x ) | (cid:29) | ˜ g (cid:48) ( x ) | we can neglect the terms with ˜ g (cid:48) ( x ) in (89), thus obtaining λ = 2 τ δ (1 − µu ) u − τ (1 − µu ) (cid:90) ∞ ˜ f (cid:48) ( x ) e − x dx. (161)Substitution of ˜ f (cid:48) ( x ) = σ ˆ h (cid:48)(cid:48) ( σx ) and Taylor expanding ˆ h (cid:48)(cid:48) as ˆ h (cid:48)(cid:48) ( x ) = ˆ h (cid:48)(cid:48) (0) + O ( σ ) immediately yields (159);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 ˜ h ( x ) = ˘ h ( x/σ ), where 0 < σ (cid:28) f ( x ) = ˘ h (cid:48) ( x/σ ) /σ = O (1 /σ ) and ˜ g ( x ) = ˘ h (cid:48)(cid:48) ( x/σ ) /σ = O (1 /σ ). Again, the sign of λ can be related to the sign of ˘ h (cid:48)(cid:48) (0), thoughthe results are now flipped: Corollary 7 (small eigenvalue for terrains with strong curvature) . Let the conditions of Theorem 3 be fulfilled.If ˜ f ( x ) = ˘ h (cid:48) ( x/σ ) /σ and ˜ g ( x ) = ˘ h (cid:48)(cid:48) ( x/σ ) /σ with < σ (cid:28) and ˘ h ( y ) , ˘ h (cid:48) ( y ) , ˘ h (cid:48)(cid:48) ( y ) → exponentially fast for | y | → ∞ , the leading (and next-leading) order expansion of (89) becomes λ = 2 τ δu − τ (1 − µu ) (cid:20) − µu σ ˘ h (cid:48)(cid:48) (0) + (1 − µu ) ˘ h (0) (cid:21) ; (162) additionally, in the double asymptotic limit τ := ε µm (cid:28) , µ (cid:28) , this further reduces to λ = 23 τ δ ˘ h (0) . (163) Furthermore, it follows that sgn λ = − sgn ˘ h (cid:48)(cid:48) (0) when µ (cid:54) = 0 , i.e. (vegetation) pulses on hilltops are unstableand in valleys are stable; and sgn λ = sgn ˘ h (0) when µ = 0 . roof. Substitution of ˜ h ( x ) = ˘ h ( x/σ ) and the use of the transformation y = x/σ in (157) yields λ = 2 δτu − τ (1 − µu ) (cid:20) − µu σ ˘ h (cid:48)(cid:48) (0) + (1 − µu ) ˘ h (0) + σ (cid:90) ∞ ˘ h ( y ) (cid:0) e − σy − − µu ) e − σy (cid:1) dy. (cid:21) (164)Taylor expanding the exponential functions then indicates the integral contributes only at order O ( δτ σ ). Hencethe claimed results follow.Thus, the corollaries in this section indicate that – under certain assumptions on the limiting behavior of thetopography function h – vegetation patterns concentrated on hilltops are stable if the terrain has weak curvatureand unstable if the terrain has strong curvature; similarly, patterns concentrated in valleys are unstable forterrains with weak curvature, but they become stable if the terrain has strong curvature. A more in-depthinspection of this phenomena can be found in section 4.4, where a few explicit terrain functions h are studiednumerically. In the previous section we found that, under certain ‘standard’ assumptions on the system’s parameters, all largeeigenvalues of a homoclinic pulse solution reside to the left of the imaginary axis. Only one small eigenvaluecan lead to destabilization of the pulse solution. Since this small eigenvalue is closely related to the translationinvariance of the system without spatially varying coefficients, it is possible to study its effects by projecting thewhole system unto the corresponding eigenspace.This derivation enables us to reduce the full PDE dynamics of (2) to a simpler ODE that describes themovement of the pulse’s location. Concretely, let P denote the location of the center of the pulse. Then thetime-evolution of P is given by dPdt = τ (cid:2) ˜ u x ( P + ) − ˜ u x ( P − ) (cid:3) , (165)where the superscripts ± denote taking the upper respectively lower limit, τ := ε µm = Da m √ m and ˜ u solves thedifferential-algebraic equation ˜ u xx + f ( x )˜ u x + g ( x )˜ u + 1 − ˜ u = 0˜ u ( P ) = µu ˜ u x ( P + ) − ˜ u x ( P − ) = u (166)We follow [4] and only give a short formal derivation of this PDE-to-ODE reduction, in section 4.1. Werefrain from going into the details of (proving) the validity of this reduction. Although the renormalization groupapproach of [6, 18] for semi-strong pulse interactions has not yet been applied to systems with inhomogeneousterms, it can naturally be extended to include these effects. However, it should be noted that, so far, the resultsand techniques of [6, 18] only cover strongly restricted region in parameter space: the general issue of validity ofthe reduction of semi-strong pulse interactions to finite dimensional settings still largely remains an open questionin the field – see also [4]. As a consequence, we formulate the main results of this section as Propositions andonly provide their formal derivations.Using the pulse location ODE (165) we use formal analysis in section 4.2 to present a scheme by which we candetermine the stability of the homoclinic pulse patterns of Theorem 2.5 for any functions f and g , i.e. withoutthe restriction on their size by which we obtained Theorem 3; in section 4.3 we (formally) validate this schemeby reducing it to the setting of Theorem 3, i.e. by assuming that f, g = O ( δ ) (with δ (cid:28) In this section we formally derive the pulse location ODE (165). Mathematically, this amounts to tracking per-turbations along translational eigenvalues; this approach is sometimes called the ‘collective coordinate method’.Specifically, in this section, we show
Proposition 6.
Let ε = am (cid:28) , τ = Da m √ m (cid:28) and µ = Dm √ ma ≤ O (1) (w.r.t. ε ). Let P denote the locationof the homoclinic pulse’s center. Then the evolution of P is described by the pulse location ODE (165) . ormal derivation, cf. [4]. We introduce the stretched travelling-wave coordinate ξ = √ mD ( x − P ( t )) = √ mD (cid:18) x − P (0) − (cid:90) t dPdt ( s ) ds (cid:19) , scale dPdt = Da m √ m c ( t ) and use scalings (10) to transform (2) to get (cid:40) − a m Dm √ ma Da m √ m c ( t ) u ξ = u ξξ − a m (cid:104) D ma u − Dm √ ma f (cid:16) D √ m ξ (cid:17) u ξ − D ma g (cid:16) D √ m ξ (cid:17) u − D √ m + uv (cid:105) − a m c ( t ) v ξ = v ξξ − v + uv (167)To find the solution in the fast region I f = [ − / √ ε, / √ ε ], close to the pulse location, we expand u and v interms of ε and look for solution of the form (cid:40) u = u + ε u + . . .v = v + ε v + . . . (168)To leading order (167) is given by (cid:26) u (cid:48)(cid:48) , v (cid:48)(cid:48) − v + u v . (169)Hence we find u to be constant and v ( ξ ) = 32 1 u sech( ξ/ . (170)The next order of (167) is (cid:26) u (cid:48)(cid:48) = u v ,v (cid:48)(cid:48) − v + 2 u v v = − c ( t ) v (cid:48) − v u . (171)It is not a priori clear whether the v -equation is solvable; the self-adjoint operator L := ∂ ξ − u v has anon-empty kernel, since L v (cid:48) = 0, and therefore the inhomogeneous v -equation is only solvable when the followingFredholm condition holds (cid:90) I f c ( t ) v (cid:48) ( η ) dη = − (cid:90) I f v ( η ) u ( η ) v (cid:48) ( η ) dη. (172)Upon integrating by parts twice on the right-hand side we obtain (cid:90) I f c ( t ) v (cid:48) ( η ) dη = − (cid:20) u (cid:48) ( η ) (cid:90) η v ( y ) dy (cid:21) η =+1 / √ εη = − / √ ε + 13 (cid:90) I f u (cid:48)(cid:48) ( η ) (cid:90) η v ( y ) dydη + h.o.t. (173)Since v is an even function, u (cid:48)(cid:48) is an even function and η (cid:55)→ (cid:82) η v ( y ) dy is an odd function. Therefore the lastintegral vanishes and we obtain c ( t ) (cid:90) I f v (cid:48) ( η ) dη = 16 (cid:20) u (cid:48) (cid:18) √ ε (cid:19) + u (cid:48) (cid:18) − √ ε (cid:19)(cid:21) (cid:90) I f v ( η ) dη. (174)The integrals over the fast field I f can be approximated by integrals over R , since v decays exponentially withinfast field. Hence we find c ( t ) = 1 u (cid:20) u (cid:48) (cid:18) √ ε (cid:19) + u (cid:48) (cid:18) − √ ε (cid:19)(cid:21) . (175)Finally, it follows from the u -equation in (171) that u (cid:48) (cid:18) √ ε (cid:19) − u (cid:48) (cid:18) − √ ε (cid:19) = (cid:90) I f u (cid:48)(cid:48) ( η ) dη = (cid:90) I f u v ( η ) dη = 6 u + h.o.t. (176)Combining this with (175) we obtain c ( t ) = 16 (cid:34) u (cid:48) (cid:18) √ ε (cid:19) − u (cid:48) (cid:18) − √ ε (cid:19) (cid:35) (177)The values of u (cid:48) ( ± / √ ε ) can be matched to the solutions ˆ u in the slow fields. Careful inspection of thescalings involved reveals u (cid:48) ( ± / √ ε ) = ˆ u x ( P ± ), where ˆ u satisfies the differential-algebraic equation (166). Since dPdt = τ c ( t ) this concludes the proof. Remark 23.
Note the link with the notation in section 2: u (cid:48) = ˆ p . See also Remark 7. .2 Stability of fixed points of pulse location ODE (165) The pulse location ODE (165) describes the movement of a pulse over time. In general, for generic functions f and g , it is not possible to solve (166) in closed form, and therefore the pulse location ODE (165) cannot beexpressed more explicitly for generic functions f and g . Thus, in general, (165) can only be solved numerically –for instance using the numerical scheme developed in [4]. Moreover, for generic f and g fixed points of (165) canonly be obtained numerically. However, when f and g obey the symmetry assumptions (A2), one can readilyobtain that P ∗ = 0 is a fixed point. It is possible to determine the stability of fixed points using (165) via directnumerics, but this can be rather time-intensive and is prone to errors close to bifurcation points. Instead, it isbetter to first use asymptotic expansions to derive a stability condition that can be checked (numerically) moreeasily. Proposition 7.
Let the conditions of Proposition 6 be satisfied, let µ (cid:28) and let P ∗ be a fixed point of (165) .Then, the eigenvalue λ – where λ = mλ , see (94) – corresponding to the pulse solution with a pulse located atthe fixed point P ∗ is given by λ = τ (cid:8) u (cid:48) ( P + ∗ ) (cid:2) ˜ u (cid:48)(cid:48) ( P + ∗ ) + ˜ w (cid:48) ( P + ∗ ) (cid:3) − u (cid:48) ( P −∗ ) (cid:2) ˜ u (cid:48)(cid:48) ( P −∗ ) + ˜ w (cid:48) ( P −∗ ) (cid:3)(cid:9) . (178) Here ˜ u and ˜ w solve the coupled ODE system u (cid:48)(cid:48) + f ˜ u (cid:48) + g ˜ u − ˜ u + 1 , w (cid:48)(cid:48) + f ˜ w (cid:48) + g ˜ w − ˜ w, ˜ u ( P ∗ ) = 0 , ˜ w ( P ±∗ ) = − ˜ u (cid:48) ( P ±∗ ) . (179) 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 ˜ u (cid:48) ( P + ∗ ) = − ˜ u (cid:48) ( P −∗ ) , ˜ u (cid:48)(cid:48) ( P + ∗ ) = ˜ u (cid:48)(cid:48) ( P −∗ ) and ˜ w (cid:48) ( P + ∗ ) = ˜ w ( P −∗ ) . Therefore, (178) reduces to λ = 2 τ u (cid:48) ( P + ∗ ) (cid:2) ˜ u (cid:48)(cid:48) ( P + ∗ ) + ˜ w (cid:48) ( P + ∗ ) (cid:3) . (180) Remark 25.
The condition µ (cid:28) in Theorem (7) is not strictly necessary. When this condition holds, thedifferential-algebraic system (166) simplifies to a normal boundary value problem, since ˜ u ( P ) = 0 to leadingorder. However, when µ = O (1) (w.r.t. ε ) the procedure explained below is still applicable and one can derivea similar result; only this time, u in (166) needs to be expanded as well and ˜ u and ˜ w satisfy the coupleddifferential-algebraic system u (cid:48)(cid:48) + f ˜ u (cid:48) + g ˜ u − ˜ u + 1 , w (cid:48)(cid:48) + f ˜ w (cid:48) + g ˜ w − ˜ w, ˜ u ( P ∗ ) = µu , ˜ w ( P ±∗ ) = − ˜ u (cid:48) ( P ±∗ ) + µw , ˜ u (cid:48) ( P + ∗ ) − ˜ u (cid:48) ( P −∗ ) = u , ˜ w (cid:48) ( P + ∗ ) − ˜ w (cid:48) ( P −∗ ) = w u + ˜ u (cid:48)(cid:48) ( P −∗ ) − ˜ u (cid:48)(cid:48) ( P + ∗ ) . (181) Formal derivation.
To find the eigenvalue λ we need to evaluate the derivative of the right-hand side of (165)at the fixed point P ∗ . That is, λ = ddP (cid:104) τ (cid:0) ˜ u (cid:48) ( P + ) − ˜ u (cid:48) ( P − ) (cid:1)(cid:105) P = P ∗ = τ (cid:34) u (cid:48) ( P + ∗ ) (cid:18) ddP ˜ u (cid:48) ( P + ) (cid:19) P = P ∗ − u (cid:48) ( P −∗ ) (cid:18) ddP ˜ u (cid:48) ( P − ) (cid:19) P = P ∗ (cid:35) . (182)By definition of the derivative ddP (cid:2) ˜ u (cid:48) ( P ± ) (cid:3) = lim φ → ˜ u (cid:48) φ (( P + φ ) ± ) − ˜ u (cid:48) ( P ± ) φ , (183)where ˜ u φ solves (166) with every P replaced by P + φ . For small φ , ˜ u φ can be related to ˜ u via a regular expansion.Specifically, let | φ | (cid:28)
1, and expand ˜ u φ = ˜ u + φ ˜ w . Substitution in (166) and careful bookkeeping readiliy shows33hat ˜ u and ˜ w satisfy (179). Finally, upon substituting the expansion for ˜ u φ into (183) and the use of a Taylorexpansion we obtain ddP (cid:2) ˜ u (cid:48) ( P ± ) (cid:3) = lim φ → ˜ u (cid:48) (( P + φ ) ± ) + φ ˜ w (cid:48) (( P + φ ) ± ) − ˜ u ( P ± ) φ = lim φ → ˜ u (cid:48) ( P ± ) + φ ˜ u (cid:48)(cid:48) ( P ± ) + φ ˜ w (cid:48) ( P ± ) − ˜ u (cid:48) ( P ± ) φ = ˜ u (cid:48)(cid:48) ( P ± ) + ˜ w (cid:48) ( P ± ) . Finally, substitution into (182) gives (178).
As an example of the use of Proposition 7, in this section we use Proposition 7 to give another proof for Theorem 3in the limit µ (cid:28)
1. This not only shows the applicability of Proposition 7 but especially the relevance of thepulse location ODE (165). Moreover, it also provides a confirmation of the validity of the formal results in thissection.
Alternative formal derivation of Theorem 3 for µ (cid:28) . Since f and g satisfy the symmetry assumption (A2),the eigenvalue λ is given by (180). Therefore, it suffices to only look at the solutions ˜ u and ˜ w to (179) for x > f, g = O ( δ ) with δ (cid:28)
1, we use regular expansions for ˜ u and ˜ w ; that is, we set˜ u = ˜ u + δ ˜ u + . . . , ˜ w = ˜ w + δ ˜ w + . . . . Substitution in (179) gives at leading order u (cid:48)(cid:48) − ˜ u + 1 , w (cid:48)(cid:48) − ˜ u , ˜ u (0) = 0 , ˜ w (0 + ) = − ˜ w (cid:48) (0 + ); (184)and at the next order, O ( δ ), we find ˜ u (cid:48)(cid:48) − ˜ u = − ˜ f ˜ u (cid:48) − ˜ g ˜ u , ˜ w (cid:48)(cid:48) − ˜ w = − ˜ f ˜ w (cid:48) − ˜ g ˜ w , ˜ u (0) = 0 , ˜ w (0 + ) = − ˜ u (cid:48) (0 + ) . (185)Using the usual techniques to solve these ODEs, one can verify that˜ u ( x ) = 1 − e − x (186)˜ u ( x ) = 12 e x (cid:90) ∞ x F ( z ) e − z dz − (cid:90) ∞ F ( z ) e − z dz + 12 e − x (cid:90) x F ( z ) e z dz (187)˜ w ( x ) = − e − x (188)˜ w ( x ) = 12 e x (cid:90) ∞ x G ( z ) e − z dz − e − x (cid:90) ∞ G ( z ) e − z dz + 12 e − x (cid:90) x G ( z ) e z dz − e − x (cid:90) ∞ F ( z ) e − z dz (189)where F ( z ) := ˜ f ( z ) e − z + ˜ g ( z )(1 − e − z ) , (190) G ( z ) := ˜ f ( z ) e − z − ˜ g ( z ) e − z . (191)Substitution of these expansions in (180) then yields λ = 23 τ [˜ u (cid:48) (0) + δ ˜ u (cid:48) (0)] [˜ u (cid:48)(cid:48) (0) + δ ˜ u (cid:48)(cid:48) (0) + ˜ w (cid:48) (0) + δ ˜ w (cid:48) (0)] + O ( δ )= 23 τ (cid:20) δ (cid:90) ∞ F ( z ) e − z e − z (cid:21) (cid:20) − δ (cid:90) ∞ ( F ( z ) + G ( z )) e − z dz (cid:21) + O ( δ )= 23 δτ (cid:90) ∞ ( F ( z ) + G ( z )) e − z dz + O ( δ )= 23 δτ (cid:90) ∞ (cid:16) f ( z ) e − z + ˜ g ( z )[1 − e − z ] e − z (cid:17) dz + O ( δ )= 23 δτ (cid:90) ∞ (cid:16) ˜ f (cid:48) ( z ) e − z + ˜ g (cid:48) ( z )(1 − e − z ) e − z (cid:17) dz + O ( δ ) . − . . . . P ∗ B (a) Bifurcation diagram for A =1 − − . . . . P ∗ B (b) Bifurcation diagram for A = − − − − . . A B c (c) Bifurcation value B c ( A ) , , , , t v − − . . . x h ( x ) (d) A = 1, B = 0 . , , , , , t v − − . . . x h ( x ) (e) A = 1, B = 1 . , , , , , t v − − . . − − . x h ( x ) (f ) A = − B = 0 . , , , , t v − − . . − − . x h ( x ) (g) A = − B = 1 . Figure 11 – Numerical results for h ( x ) = Ae − Bx . Shown are bifurcation diagrams for A = 1 (a) and A = − B c ( A ) of the pitchfork bifurcation (c), and (parts of) various simulations of the full PDEillustrating the change of stability along with a plot of the function h ( x ) (d-g). The green areas in (c) indicate theparameter region in which the fixed point P ∗ = 0 is stable. In the PDE simulations we have used parameters a = 0 . m = 0 . D = 0 .
01 and taken x ∈ [ − , Finally, we note that the eigenvalue has been rescaled as λ = mλ in Theorem 2. Since τ /m = ε µ and u = 3in the limit µ (cid:28)
1, we have indeed recovered (156), i.e. Theorem 3, in the case µ (cid:28) In this section, we study a few explicit functions f and g ; in all examples we specify a function h and take f = h (cid:48) , g = h (cid:48)(cid:48) . Not all functions we consider here limit to 0 as | x | → ∞ ; that is, some violate assumption (A4).Therefore, these examples also form an outlook, illustrating how the results in this paper are expected to extendbeyond the imposed assumptions on functions f and g . Specifically, we consider the following four examples:(i) h ( x ) = Ae − Bx , ( A ∈ R , B > h ( x ) = A sech( Bx ), ( A ∈ R , B > h ( x ) = A cos( Bx ), ( A ∈ R , B > h ( x ) = − βx )), ( β > | x |→∞ f ( x ) , g ( x ) = 0 in cases (i)–(ii), which therefore satisfy assumption (A4). In case (iii) f and g are periodic when | x | (cid:29)
1; in case (iv) f and g do have well-defined (though non-zero) limits for | x | → ∞ . Remark 26.
Note that
A > in (i)–(ii) corresponds to ‘hill-like’ topographies and A < 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 (165) and Proposition 7, we have tracked the fixed points and their stabilityfor these examples in the limit µ (cid:28)
1, using numerical continuation methods. The resulting bifurcation diagramsfor (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.35 − . . . . . P ∗ B (a) Bifurcation diagram for A =1 − − . . . . P ∗ B (b) Bifurcation diagram for A = − − . . A B c (c) Bifurcation value B c ( A ) , , , , x t v − − . . . x h ( x ) (d) A = 1, B = 0 . , , , , x t v − − . . . x h ( x ) (e) A = 1, B = 2 . , , , , t v − − . . − . − x h ( x ) (f ) A = − B = 0 . , , , , t v − − . . x h ( x ) (g) A = − B = 2 . Figure 12 – Numerical results for h ( x ) = A sech( Bx ). Shown are bifurcation diagrams (solid for stable; dashedfor unstable fixed points) for A = 1 (a) and A = − B c ( A ) of the pitchfork bifurcation(c), and (parts of) various simulations of the full PDE illustrating the change of stability along with a plot of thefunction h ( x ) (d-g). The green areas in (c) indicate the parameter region in which the fixed point P ∗ = 0 is stable.In the PDE simulations we have used parameters a = 0 . m = 0 . D = 0 .
01 and taken x ∈ [ − , for weak curvature topographies – the stability of these fixed points is determined by the sign of A : A >
A < B –i.e. topographies with strong curvature– the stability of those fixed points changes through a pitchfork bifurcation and new behavior is observed. Incase (iii) this even leads to the possibility that both the tops ( BP = 0) as well as the valleys ( BP = ± π ) formstable fixed points of (165). 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) and13(b) (for topographies (i), (ii) and (iii)). Remark 27.
Theorem 3, and in particular (156) and (158) , provide a leading order analytic expression for B c (0) . Evaluating these yields B c (0) ≈ . (i), B c (0) ≈ . (ii) and B c (0) = √ (iii), which is confirmed bythe numerical continuation that indicate B c (0) ≈ . (i), B c (0) ≈ . (ii) and B c (0) = 1 . (iii). Note that A = 0 is, indeed, just the flat terrain h ( x ) ≡ ; 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 pointsand, for well-chosen parameter values, these simulations show convergence to fixed points not located at thepoint of symmetry. Note also that in the case of periodic topography (i.e. case (iii)), there indeed is a regionof 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 thesame 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 µ (cid:28) dPdt = τ (cid:104) (cosh( βP ) I ( P )) − (cosh( βP ) I ( P )) (cid:105) , (192)36 π − π π π BP ∗ B (a) Bifurcation diagram for A = 1 − − − A B c (b) Bifurcation value B c ( A ) , , , , x t v − π − π/ π/ π − Bx h ( x ) (c) A = 1, B = 1 , , , t v − π − π/ π/ π − Bx h ( x ) (d) A = 1, B = 1 . , , , , t v − π − π/ π/ π − Bx h ( x ) (e) A = 1, B = 1 . , , , , t v − π − π/ π/ π − Bx h ( x ) (f ) A = 1, B = 2 Figure 13 – Numerical results for h ( x ) = A cos( Bx ). Shown are the bifurcation diagram (solid for stable; dashedfor unstable fixed points) for A = 1 (a), the bifurcation value B c ( A ) of the pitchfork bifurcation at x = 0 (b), and(parts of) various simulations of the full PDE illustrating the change of stability along with a plot of the function h ( x ) (c-f). The green areas in (b) indicate the parameter region in which the fixed point P ∗ = 0 is stable. In thePDE simulations we have used parameters a = 0 . m = 0 . D = 0 .
002 and taken x ∈ [ − , where I ( P ) := (cid:90) ∞ P e r ( P − z ) sech( βz ) dz ; I ( P ) := (cid:90) P −∞ e − r ( P − z ) sech( βz ) dz. (193)Thus, a point P ∗ is a fixed point if and only if I ( P ∗ ) = I ( P ∗ ). Straightforward inspection reveals that P ∗ = 0therefore is the unique fixed point in case (iv) for all values of β >
0. By Proposition 7 and equation (180) thecorresponding (small) eigenvalue λ can be approximated by λ = 2 τ I (0) ( r I (0) − . (194)Upon noting that r I (0) − − β (cid:90) ∞ sech( βz ) tanh( βz ) e − rz dz < , (195)it is clear that λ <
0. Hence, P ∗ = 0 is the only fixed point of (192) in case (iv), which is (globally) stable –for all β >
0. Direct PDE simulations verify this – even when the asymptotic limit µ (cid:28) The focus in this article has been on single pulse solutions to (2). As a short encore we briefly discuss thepossibility of stationary multi-pulse solutions – i.e. solutions with multiple fast excursions. The movement ofthese solutions can be captured in an ODE much akin to 165. Specifically, let P , . . . , P N denote the location of N pulses. Then their movement is described by the ODE dP j dt = τ (cid:2) ˜ u x ( P + j ) − ˜ u x ( P − j ) (cid:3) , ( j = 1 , . . . , N ) (196)37 , , , , t v − − . . − − . x h ( x ) Figure 14 – Direct numerical PDE simulation for h ( x ) = − βx ) for β = 1 along with a plot of the function h ( x ). In the PDE simulation we have used the parameters a = 0 . , m = 0 . , D = 0 .
01 and taken x ∈ [ − , , t v − − x h ( x ) (a) h ( x ) = 0 , , , , , x t v − − . . − − . x h ( x ) (b) h ( x ) = − x ) , , , , , x t v − − − . x h ( x ) (c) h ( x ) = e − x / , , , , , x t v − − . . . x h ( x ) (d) h ( x ) = sech( x/ Figure 15 – Numerical simulation of several multi-pulse solutions to (2) for various h , with f = h (cid:48) and g = h (cid:48)(cid:48) . (a) h ( x ) = 0: no stable stationary multi-pulse solution is found; (b) h ( x ) = − x ): the existence of a stable two-pulse solution; (c) h ( x ) = e − x / : a stable three-pulse solution; (d) h ( x ) = sech( x/ x -domain is shown for clarity. Also note that,using (199), it is found that P ∗ ≈ .
51 in (b). where ˜ u satisfies the differential-algebraic system ˜ u xx + f ( x )˜ u x + g ( x )˜ u + 1 − ˜ u = 0˜ u ( P j ) = µu j ( j = 1 , . . . , N )˜ u x ( P + j ) − ˜ u x ( P − j ) = u j ( j = 1 , . . . , N ) (197)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 ≡
0, it is well-known that stationary multi-pulse solutions do not exist [13,4]. In fact, from (196) one can verify that in 2-pulse solutions the pulses typically move away from each otherwith a speed proportional to e − ∆ P , where ∆ P := P − P 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 move-ment. Therefore stationary pulse solutions do exist in (2) for well-chosen f and g . In Figure 15 we show severalnumerical 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 O ( f, g ) effect on the movement speed of the pulses. Find-ing fixed points of (196) – i.e. finding stationary multi-pulse solutions to (2) – thus boils down to balancing twoeffects of different size. In particular, if f, g = O ( δ ) , only multi-pulse solutions exist with ∆ P = O ( − ln( δ )) (cid:29) .In this case, existence of stationary multi-pulse solutions can be established rigorously by asymptotic analysis andthe 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 canbe stable and unstable depending on the parameter values – similar to the one-pulse variants. A description ofhow 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-pulsesolutions (however, see Remark 28 for the case of small f , g ). We do remark however that stationary multi-pulsesolutions can be constructed for f and g such that (197) can be solved explicitly, as illustrated by the followingproposition. 38 roposition 8. Let h ( x ) = − βx ) , β > , f = h (cid:48) , g = h (cid:48)(cid:48) and let µ (cid:28) . Then there exists a P ∗ > such that (2) admits a stationary symmetric two-pulse solutions with pulses at P = − P ∗ and P = P ∗ .Formal derivation. By symmetry of the desired two-pulse solution, we may set P = P , P = − P . Moreover,necessarily ˜ u (cid:48) (0) = 0. Since µ (cid:28)
1, to leading order we have ˜ u ( P ) = ˜ u ( − P ) = 0. Therefore ˜ u is given to leadingorder by ˜ u ( x ) = ˆ u b ( x ) − ˆ u b ( − P )ˆ u − ( − P ) ˆ u − ( x ) , x < − P, ˆ u b ( x ) − ˆ u b ( P )ˆ u + ( P )+ˆ u − ( P ) (ˆ u + ( x ) + ˆ u − ( x )) , − P < x < P ;ˆ u b ( x ) − ˆ u b ( P )ˆ u + ( P ) ˆ u + ( x ) , x > P ; (198)where ˜ u ± and ˜ u b are as in Corollary 3. To have stationary pulse solutions, by (196) we need to have T ( P ) := ˆ u (cid:48) b ( P ) − ˆ u b ( P ) (cid:34) (cid:112) β (cid:16) tanh( (cid:112) β P ) − (cid:17) + β tanh( βP ) (cid:35) = 0 , (199)Upon noting that T (0) = 12 (cid:90) ∞ e − √ β z sech( βz ) dz > , (200)and, since lim P →∞ ˆ u b ( P ) = 1 and lim P →∞ ˆ u (cid:48) b ( P ) = 0,lim P →∞ T ( P ) = − β < , (201)continuity of T guarantees the existence of P ∗ > Remark 30.
This result can be established rigorously by geometric singular perturbation theory, using the meth-ods detailed in section 2. We refrain from giving the details of this procedure.
In this paper, we studied pulse solutions in a reaction-advection-diffusion system with spatially varying coef-ficients. The existence of stationary pulse solutions at a point of symmetry was established by combining theusual techniques from geometric singular perturbation theory with the tools from the theory of exponentialdichotomies. The latter has been used to generate a saddle-like structure in the slow subsystem, and to obtainbounds on the stable/unstable manifolds of this subsystem. These techniques have also been used to determinethe spectral stability of these pulse solutions. None of these concepts or ideas are model-dependent and thereforecould 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 boundedto the stable half-plane, under conditions similar to the usual, constant coefficient case. Although we did notfocus on the dynamics of solutions when a large eigenvalue crosses the imaginary axis, simulations show the usualpulse annihilation and pulse splitting phenomena. However, the introduction of spatially varying coefficients doeshave a significant effect on the so-called ‘small’ eigenvalues (close to λ = 0) because of the break-down of thetranslation 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]. Weexpect that a careful tuning of f and g can either prevent or force these subsequent bifurcations, which mayhave 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 (cid:48) , g = h (cid:48)(cid:48) (where h is used to modelthe topography of a dryland ecosystem). Here, we were able to link the stability of (stationary) pulse solutionto the curvature of h . If the curvature is weak, the pulse is stable if h (cid:48)(cid:48) (0) < h (cid:48)(cid:48) (0) >
0; forstrong curvature the opposite is true: the pulse is stable if h (cid:48)(cid:48) (0) > h (cid:48)(cid:48) (0) <
0. We found thatthis change in stability typically happens via a pitchfork bifurcation, and showed that the associated parametercombinations 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 (2) does not posses the symmetry ( x, u ) → ( − 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 ofsubsequent work. 39oreover, in case of spatially varying coefficients, the system (2) can also posses stationary multi-pulsesolutions – i.e. solutions that have multiple fast excursions. When f, g ≡
0, these solutions do not exist. Becausethe spatially varying coefficients break the translation invariance of the system, these multi-pulse solutions canexist – for well-chosen functions f and g . In this article we gave numerical evidence for this and showed theirexistence for a specific choice of functions. We do not think their existence can be proven in as much generalityas the existence of stationary one pulse solutions – certainly, the bounds used in this paper, provided by thetheory 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 subsequentpulses then becomes asymptotically large and asymptotic analysis needs to be done with great care to keep trackof 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 usedto 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 thisapplication. Specifically, the treatment of a spatially varying height function h is new and is inherently morerealistic than taking a constant topography (or a constantly sloped topography) as has been done in the past (seee.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 downhillmovement of pulses. This aligns better with measurements, where also both uphill and downhill movement canbe observed – even within the same general region [23, 5]. In this regard, the study in this paper can be seen asa 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 onperiodic topographies. This work was funded by NWO’s Mathematics of Planet Earth program.
References [1] J Alexander, Robert Gardner, and CKRT Jones. A topological invariant arising in the stability analysis oftravelling waves.
J. reine angew. Math , 410(167-212):143, 1990.[2] Daniele Avitabile, Victor F Bren˜a, and Michael J Ward. Spot dynamics in a reaction-diffusion model ofplant root hair initiation.
SIAM Journal on Applied Mathematics , 78(1):291–319, 2018.[3] Robbin Bastiaansen, Paul Carter, and Arjen Doelman. Stable planar vegetation stripe patterns on slopedterrain in dryland ecosystems. submitted , 2018.[4] Robbin Bastiaansen and Arjen Doelman. The dynamics of disappearing pulses in a singularly perturbedreaction-diffusion system with parameters that vary in time and space.
Physica D , ?:?, ?[5] Robbin Bastiaansen, Olfa Ja¨ıbi, Vincent Debaluwe, Maarten Eppinga, Koen Siteur, Eric Siero, St´ephaneMermozh, Alexandre Bouvet, Arjen Doelman, and Max Rietkerk. Multi-stability of model and realdryland ecosystems through spatial self-organization.
Proceedings of the National Academy of Sciences ,115(44):11256–11261, 2018.[6] Thomas Bellsky, Arjen Doelman, Tasso J Kaper, and Keith Promislow. Adiabatic stability under semi-strong interactions: the weakly damped regime.
Indiana University Mathematics Journal , pages 1809–1859,2013.[7] Henri Berestycki, Juncheng Wei, and Matthias Winter. Existence of symmetric and asymmetric spikes fora crime hotspot model.
SIAM Journal on Mathematical Analysis , 46(1):691–719, 2014.[8] Victor Bren˜a Medina, Alan R Champneys, C Grierson, and Michael J Ward. Mathematical modeling ofplant root hair initiation: Dynamics of localized patches.
SIAM Journal on Applied Dynamical Systems ,13(1):210–248, 2014.[9] Victor F Bren˜a Medina, Daniele Avitabile, Alan R Champneys, and Michael J Ward. Stripe to spottransition in a plant root hair initiation model.
SIAM Journal on Applied Mathematics , 75(3):1090–1119,2015. 4010] Wan Chen and Michael J Ward. Oscillatory instabilities and dynamics of multi-spike patterns for theone-dimensional Gray-Scott model.
European Journal of Applied Mathematics , 20(2):187–214, 2009.[11] W.A. Coppel.
Dichotomies in Stability Theory , volume 629 of
Lecture Notes in Mathematics . Springer-Verlag, Berlin-New York, 1978.[12] Bj¨orn de Rijk, Arjen Doelman, and Jens Rademacher. Spectra and stability of spatially periodic pulsepatterns: Evans function factorization via Riccati transformation.
SIAM J. Math. Anal. , 48(1):61–121,2016.[13] Arjen Doelman, Wiktor Eckhaus, and Tasso J. Kaper. Slowly modulated two-pulse solutions in the Gray-Scott model i: Asymptotic construction and stability.
SIAM Journal on Applied Mathematics , 61(3):1080–1102, 2000.[14] Arjen Doelman, Wiktor Eckhaus, and Tasso J Kaper. Slowly modulated two-pulse solutions in the Gray-Scott model ii: Geometric theory, bifurcations, and splitting dynamics.
SIAM Journal on Applied Mathe-matics , 61(6):2036–2062, 2001.[15] Arjen Doelman, Robert A Gardner, and Tasso J Kaper. Stability analysis of singular patterns in the 1DGray-Scott model: a matched asymptotics approach.
Physica D: Nonlinear Phenomena , 122(1-4):1–36,1998.[16] Arjen Doelman, Robert A Gardner, and Tasso J Kaper. Large stable pulse solutions in reaction-diffusionequations.
Indiana University Mathematics Journal , 50(1):443–507, 2001.[17] Arjen Doelman and Tasso J Kaper. Semistrong pulse interactions in a class of coupled reaction-diffusionequations.
SIAM Journal on Applied Dynamical Systems , 2(1):53–96, 2003.[18] Arjen Doelman, Tasso J Kaper, and Keith Promislow. Nonlinear asymptotic stability of the semistrong pulsedynamics in a regularized Gierer–Meinhardt model.
SIAM Journal on Mathematical Analysis , 38(6):1760–1787, 2007.[19] Arjen Doelman, Jens D. M. Rademacher, and Sjors van der Stelt. Hopf dances near the tips of Busseballoons.
Discrete and Continuous Dynamical Systems-Series S , 5(1):61–92, 2012.[20] Arjen Doelman and Harmen van der Ploeg. Homoclinic stripe patterns.
SIAM Journal on Applied DynamicalSystems , 1(1):65–104, 2002.[21] Arjen Doelman, Peter van Heijster, and F Xie. A geometric approach to stationary defect solutions in onespace dimension.
SIAM Journal on Applied Dynamical Systems , 15(2):655–712, 2016.[22] Arjen Doelman and Frits Veerman. An explicit theory for pulses in two component, singularly perturbed,reaction–diffusion equations.
Journal of Dynamics and Differential Equations , 27(3-4):555–595, 2015.[23] David L Dunkerley. Vegetation mosaics of arid western New South Wales, Australia: Considerations oftheir origin and persistence. In
Patterns of Land Degradation in Drylands , pages 315–345. Springer, 2014.[24] Neil Fenichel. Geometric singular perturbation theory for ordinary differential equations.
Journal of differ-ential equations , 31(1):53–98, 1979.[25] Robert Gardner and Christopher KRT Jones. Stability of travelling wave solutions of diffusive predator-preysystems.
Transactions of the American Mathematical Society , 327(2):465–524, 1991.[26] Christopher KRT Jones. Geometric singular perturbation theory. In
Dynamical systems , pages 44–118.Springer, 1995.[27] Christopher A Klausmeier. Regular and irregular patterns in semiarid vegetation.
Science , 284(5421):1826–1828, 1999.[28] AJ Koch and Hans Meinhardt. Biological pattern formation: from basic mechanisms to complex structures.
Reviews of modern physics , 66(4):1481, 1994. 4129] Theodore Kolokolnikov, Michael J. Ward, and Juncheng Wei. The existence and stability of spike equilibriain the one-dimensional GrayScott model on a finite domain.
Applied Mathematics Letters , 18(8):951 – 956,2005.[30] Christian Kuehn.
Multiple time scale dynamics , volume 191 of
Applied Mathematical Sciences . Springer,Cham, 2015.[31] Philip K Maini. Applications of mathematical modelling to biological pattern formation. In
Coherentstructures in complex systems , pages 205–217. Springer, 2001.[32] Hans Meinhardt. Models of biological pattern formation: from elementary steps to the organization ofembryonic axes.
Current topics in developmental biology , 81:1–63, 2008.[33] I Moyles, WH Tse, and MJ Ward. Explicitly solvable nonlocal eigenvalue problems and the stability oflocalized stripes in reaction-diffusion systems.
Studies in Applied Mathematics , 136(1):89–136, 2016.[34] Yasumasa Nishiura, Yoshihito Oyama, Kei-ichi Ueda, et al. Dynamics of traveling pulses in heterogeneousmedia of jump type.
Hokkaido mathematical journal , 36(1):207–242, 2007.[35] Yasumasa Nishiura, Takashi Teramoto, Xiaohui Yuan, and Kei-Ichi Ueda. Dynamics of traveling pulses inheterogeneous media.
Chaos: An interdisciplinary journal of nonlinear science , 17(3):037104, 2007.[36] V Rottsch¨afer, JC Tzou, and MJ Ward. Transition to blow-up in a reaction–diffusion model with localizedspike solutions.
European Journal of Applied Mathematics , 28(6):1015–1055, 2017.[37] Lotte Sewalt and Arjen Doelman. Spatially periodic multipulse patterns in a generalized Klausmeier-Gray-Scott model.
SIAM J. Appl. Dyn. Syst. , 16(2):1113–1163, 2017.[38] Jonathan A Sherratt. History-dependent patterns of whole ecosystems.
Ecological Complexity , 14:8–20,2013.[39] Jonathan A Sherratt. Using wavelength and slope to infer the historical origin of semiarid vegetation bands.
Proceedings of the National Academy of Sciences , page 201420171, 2015.[40] Koen Siteur, Eric Siero, Maarten B Eppinga, Jens DM Rademacher, Arjen Doelman, and Max Rietkerk.Beyond Turing: The response of patterned ecosystems to environmental change.
Ecological Complexity ,20:81–96, 2014.[41] Wentao Sun, Michael J. Ward, and Robert Russell. The slow dynamics of two-spike solutions for the Gray-Scott and Gierer–Meinhardt systems: Competition and oscillatory instabilities.
SIAM Journal on AppliedDynamical Systems , 4(4):904–953, 2005.[42] Andrei Tikhonov. On the dependence of the solutions of differential equations on a small parameter.
Matematicheskii sbornik , 64(2):193–204, 1948.[43] Alan Mathison Turing. The chemical basis of morphogenesis.
Phil. Trans. R. Soc. Lond. B , 237(641):37–72,1952.[44] Peter Van Heijster, Arjen Doelman, Tasso J Kaper, Yasumasa Nishiura, and Kei-Ichi Ueda. Pinned frontsin heterogeneous media of jump type.
Nonlinearity , 24(1):127, 2010.[45] Frits Veerman and Arjen Doelman. Pulses in a Gierer–Meinhardt equation with a slow nonlinearity.
SIAMJournal on Applied Dynamical Systems , 12(1):28–60, 2013.[46] Juncheng Wei and Matthias Winter. Stable spike clusters for the one-dimensional Gierer–Meinhardt system.
European Journal of Applied Mathematics , 28(4):576–635, 2017.[47] Juncheng Wei, Matthias Winter, and Wen Yang. Stable spike clusters for the precursor Gierer–Meinhardtsystem in R . Calculus of Variations and Partial Differential Equations , 56(5):142, 2017.[48] Jack Xin. Front propagation in heterogeneous media.
SIAM review , 42(2):161–230, 2000.[49] Xiaohui Yuan, Takashi Teramoto, and Yasumasa Nishiura. Heterogeneity-induced defect bifurcation andpulse dynamics for a three-component reaction-diffusion system.