The existence of localized vegetation patterns in a systematically reduced model for dryland vegetation

In this paper we consider the 2-component reaction-diffusion model that was recently obtained by a systematic reduction of the 3-component Gilad et al. model for dryland ecosystem dynamics. The nonlinear structure of this model is more involved than other more conceptual models, such as the extended Klausmeier model, and the analysis a priori is more complicated. However, the present model has a strong advantage over these more conceptual models in that it can be more directly linked to ecological mechanisms and observations. Moreover, we find that the model exhibits a richness of analytically tractable patterns that exceeds that of Klausmeier-type models. Our study focuses on the 4-dimensional dynamical system associated with the reaction-diffusion model by considering traveling waves in 1 spatial dimension. We use the methods of geometric singular perturbation theory to establish the existence of a multitude of heteroclinic/homoclinic/periodic orbits that `jump' between (normally hyperbolic) slow manifolds, representing various kinds of localized vegetation patterns. The basic 1-front invasion patterns and 2-front spot/gap patterns that form the starting point of our analysis have a direct ecological interpretation and appear naturally in simulations of the model. By exploiting the rich nonlinear structure of the model, we construct many multi-front patterns that are novel, both from the ecological and the mathematical point of view. In fact, we argue that these orbits/patterns are not specific for the model considered here, but will also occur in a much more general (singularly perturbed reaction-diffusion) setting. We conclude with a discussion of the ecological and mathematical implications of our findings.


Introduction
Ecosystems consist of organisms that interact among themselves and with their environment.These interactions involve various kinds of feedback processes that may combine to form positive feedback loops and instabilities when environmental conditions change [48,49].In many ecosystems -drylands, peatlands, savannas, mussel beds, coral reefs, and ribbon forests -the leading feedback processes have different spatial scales: a short-range facilitation by local modification of the environment versus a long-range competition for resources [57].Like the well-established activator-inhibitor principle in bio-chemical systems [52], the combination of these scale-dependent feedback mechanisms can induce instabilities that result in large-scale spatial patterns, which are similar to a wide variety of vegetation patterns observed in drylands, peatlands, savannas and undersea [15,61,5,32,58,56,27].Varying climatic conditions and human disturbances may continue to propel ecosystem dynamics.Ecosystem response to decreasing rainfall, for example, may take the form of abrupt collapse to a nonproductive 'desert state' [59,70,56], or involve gradual desertification, consisting of a cascade of state transitions to sparser vegetation [65,4], or gradual vegetation retreat by front propagation [6,76].Understanding the dynamics of spatially extended ecosystems has become an active field of research in the last two decades -within communities of ecologists, environmental scientists, mathematicians and physicists.Apart from its obvious environmental and societal relevance, the phenomena exhibited pose fundamental challenges to the research field of pattern formation.
Several models of increasing complexity have been proposed in the past two decades.Of these, the models that have received most attention are the one-component model by Lefever and Lejeune [45], the two-component models by Klausmeier [42] and von Hardeberg et al. [70], and the three-component models by Rietkerk et al. [55] and Gilad et al. [28].A basic difference between these models is the manner by which they describe water dynamics.The Lefever-Lejeune model does not describe water dynamics at all, the Klausmeier model does describe water dynamics but does not make a clear distinction between soil water and surface water [69], while the von Hardenberg et al. model only takes soil water into account.The Rietkerk model and the Gilad et al. model describe both soil water and surface water dynamics and, therefore, capture more aspects of real dryland ecosystems.A major difference between these two models is the inclusion of water conduction by laterally spread roots, as an additional water-transport mechanism, in the Gilad et al. model.Despite these differences, all models appear to share a similar bifurcation structure, as analytical and numericalcontinuation studies reveal [46,33,16,75], except the Klausmeier model.This structure includes, in particular, a stationary uniform instability (i.e.involving the monotonic growth of spatially uniform perturbations) of the bare soil (zero biomass) state as the precipitation rate exceeds a threshold value.The Klasumeier model fails to capture that instability, leaving the bare soil state stable at all precipitation values.This behavior limits the applicability of the Klausmeier model to ecological contexts where the bare soil state is stabilized at relatively high precipitation rates, e.g. by high evaporation rates.Nevertheless, of all models, the Klausmeier model and its extension to include water diffusion have been studied to a greater extent [8,62,63,69], partly because the extended form coincides with the much studied Gray-Scott model for autocatalytic chemical reactions -see [3,10,60] and the references therein.
All models have been analyzed mathematically to various extents.Two main analytical approaches can be distinguished in these studies (see however Goto et al. [30]); linear stability and weakly nonlinear analysis near instability points [46,13,33,31,26,69], and singular perturbation analysis, based on the disparate length scales associated with biomass (short) and water (long) [8,3,10,60].Studies of the first category are strictly valid only near instability points, although they do capture essential parts of the bifurcation structure even far from these points and are quite insightful in this respect.By contrast, studies of the second category apply to the strongly nonlinear 'far-from-equilibrium' regime, where desertification transitions take place, and are, potentially, of higher ecological interest.So far, however, these studies have been limited to the simpler and less realistic Klausmeier model.
In this paper we apply a geometric singular perturbation analysis to a reduced version of the Gilad et al. model in order to study the existence of various forms of localized patterns.Singular perturbation theory has already been applied to three-component models -see for instance [23,68] -and could be applied, in principle, to the non-local three-component Gilad et al. model.Here we choose to consider ecological contexts that allow to reduce that model to a local two-component model for the vegetation biomass and the soil water content.Specifically, we assume soil types characterized by high infiltration rates of surface water into the soil, such as sandy soil, and plant species with laterally confined root zones (see A for more details).These conditions are met, for example, by Namibian grasslands showing localized and extended gap patterns ('fairy circles') [77].We further simplify the problem by assuming one space dimension.The reduced model reads: where B( X, T ) ≥ 0 and W ( X, T ) ≥ 0 represent areal densities of biomass and soil water, respectively, and X ∈ R, T ∈ R + are the space and time coordinates.In the biomass ( B) equation, Λ represents the biomass growth rate coefficient, K the maximal standing biomass, E is a measure for the root-to-shoot ratio, M the plant mortality rate and D B the seed-dispersal or clonal growth rate, while in the water ( W ) equation, P represents the precipitation rate, N the evaporation rate, R the reduction of the evaporation rate due to shading, Γ the water-uptake rate coefficient and D W the effective soil water diffusivity.Notice that the power of the factor (1 + E B) in both equations is unity, whereas in the reduced model in [77] the power is two.This difference stems from the consideration in this study of one space dimension rather than two (see A).
From the ecological point of view, the advantage in studying model (1.1) over the much analyzed Klausmeier model lies in the fact that it has been systematically derived from a more extended model that better captures relevant ecological processes, such as water uptake by plant roots (controlled by E), reduced evaporation by shading (controlled by R), and late-growth constraints, such as self-shading (controlled by K) -see [28,29,47,61].As a consequence, (mathematical) insights in (1.1) can be linked to ecological observations and mechanisms in a direct fashion.Naturally, there also is a disadvantage to analyzing a model that incorporates concrete ecological mechanisms: the more involved -algebraically more complex -nonlinear structure of (1.1) a priori makes it less suitable for an analytical study than the Klausmeier model (or other more conceptual models).However, that apparent disadvantage turned around into an advantage: we will find that the reduced model transcends by far the Klausmeier model in terms of richness of analytically tractable pattern solutions.
The model equations (1.1) represent a singularly perturbed system, because of the low seed-dispersal rate as compared with soil water diffusion, that is, D B D W [29,69,77].To make this explicit and to simplify (1.1) as much as possible, we introduce the following scalings, and set, By also introducing our main parameters, we arrive at, in which, A more detailed derivation of the scaled equations (1.5) from (1.1) is given in B. Since the signs of the parameters in (1.5) will turn out to be crucial in the upcoming analysis, we note explicitly that a, Ψ, Φ, Θ ≥ 0 while Ω ∈ R, i.e.Ω may be negative.
We study pattern formation in (1.1) by analyzing (1.5) using the methods of (geometric) singular perturbation theory [39,41] and thus 'exploit' the fact that ε 1 (1.4).In fact -apart from some observations in section 2.3 and the discussion section 4.2 -we focus completely on the 'spatial' 4-dimensional dynamical system that is obtained from (1.5) by considering 'simple' solutions that are stationary in a co-moving frame traveling with constant speed c.More specifically, in this paper we study the existence of traveling (and stationary) solutions -in particular localized (multi-)front solutions connecting a (uniform) bare soil state to a uniform vegetation state, or a bare soil state to itself (with a 'passage' along a vegetated state), etc. -by taking the classical approach of introducing a (uniformly) traveling coordinate ξ = x − ct, with speed c ∈ R an a priori free O(1) parameter (w.r.t. the asymptotically small parameter ε).By setting (B(x, t), W (x, t)) = (b(ξ), w(ξ)) and introducing p = b ξ and q = 1 ε w ξ , PDE (1.5) (1.7) Figure 2: Four sketches of 'higher order' localized patterns constructed in this paper.(a) A secondary traveling 1-front, the second one in a countable family of traveling 1-fronts between the bare soil state and a uniform vegetated state -all traveling with different speeds -that starts with the primary 1-front of Fig. 1(a) (Theorem 3.6).(b) The limiting orbit of the family sketched in (a): a 1-front connection between the bare soil state and a spatially periodic vegetation state (Theorem 3.5).(c,d) The first 2 representations of a (countable) 'higher order' family of localized (stationary, homoclinic 2-front) spot patterns with an increasing number of 'spatial oscillations' (Theorem 3.12).Fig. 1 shows four basic patterns that naturally appear in simulations of (1.5) and have identifiable ecological counterparts: vegetation fronts (ecotones), isolated vegetation spots and gaps, and periodic patterns [50,25,24,27].These patterns are rigorously constructed by the methods of singular perturbation theory in section 3. From the geometrical point of view, these constructions are natural and thus relatively straightforward: all patterns in Fig. 1 'jump' between two well-defined slow manifolds (of (1.7)) -see Theorems 3.4,3.11,3.13,and 3.15.Therefore, the main work in establishing these results lies in resolving technical issues (which can be done by the preparations of section 2).However, the preparations of section 2 also form the origin of the construction of a surprisingly rich 'space' of traveling and/or stationary patterns that goes way beyond those exhibited in Fig. 1 -see for instance the sketches of Fig. 2.These are novel patterns, at least from the point of view of explicit rigorous mathematical constructions in multi-component reaction-diffusion equations.However, similar patterns have been analyzed as (perturbations of) heteroclinic networks in a more abstract framework -see [53,54] and the references therein.Moreover, patterns similar to those of Fig. 2 have been observed in simulations of the Klausmeier-Gray-Scott model [74], although with parameter settings beyond that for which the mathematical singular-perturbation approach can be applied.
Here, our motivation to study these patterns is primarily ecological; however, we claim that patterns like these must also occur generically in the setting of a completely general class of singularly perturbed 2-component reactiondiffusion systems -as we will motivate in more detail in section 4.2.Thus, our explicit analysis of model (1.5) provides novel mathematical insights beyond that of the present ecological setting.The driving mechanism behind these patterns originates from the perturbed integrable flow on the slow manifolds associated with (1.7) -see sections 2.2 and 2.4.The perturbation terms are generically introduced by the O(ε) differences between the slow manifold and its ε → 0 limit, and they transform the (Hamiltonian) integrable reduced slow flow to a (planar) 'nonlinear oscillator with nonlinear friction' that can be studied by explicit Melnikov methods.Typically, one for instance expects (and finds: Theorem 2.4) persistent periodic solutions on the slow manifold.Associated with these persisting periodic solutions, one can subsequently construct heteroclinic 1-front connections between a critical point -representing the uniform bare soil state in the ecological setting -and such a periodic pattern (Theorems 3.5 and 3.9 and Fig. 2b) and a countable family of 'higher order' heteroclinic 1-fronts between critical points that limits on these orbits (Theorem 3.6 and Fig. 2a -where we note that Fig. 1a represents the very first -primary -member of this family).In the case of (stationary) localized spot patterns, one can construct a countable family of connections that follow the periodic orbit for arbitrarily many 'spatial oscillations' (Theorem 3.12 and Fig. 2c, 2d).Combining these insights with the ideas of [21], one may even construct many increasingly complex families of spatially periodic and aperiodic multi-spot/gap patterns (Corollary 3.16 and section 3.6).Moreover, we can explicitly study the associated bifurcation scenarios: in section 3.3 we present a scenario of cascading saddle-node bifurcations of heteroclinic 1-front connections starting from no such orbits to countably many -all traveling with different speed (Theorem 3.6 and Figs.1a, 2a and 2b) -back to 1 unique 1-front pattern (of the type presented in Fig. 1a) -see Fig. 9 in section 3.3.
Finally, we illustrate our analytical findings by several numerical simulations of PDE model (1.1)/(1.7)-see also Fig. 1.We did not systematically investigate the question whether all heteroclinic/homoclinic/periodic (multi-)front orbits of (1.7) constructed here indeed may be (numerically) observed as stable patterns in (1.7), either for general parameter combinations in (1.5) or for the more restricted class of ecologically relevant parameter combinations.This will be the subject of future work, as will be the analytical question about the spectral stability of the constructed patterns.These issues will be discussed more extensively in section 4.2, where we will also discuss further implications of our findings -both from the mathematical as well as from the ecological point of view.
The set-up of this paper is as follows.Section 2 is a preparatory section: in section 2.1 and 2.2 we consider the fast and slow reduced problems associated with (1.7), followed by a brief section -section 2.3 -in which we discuss the nature (and stability) of the critical points of (1.7) as uniform vegetated states in (1.5); in section 2. 4 we study the full, perturbed, slow flow on the slow manifolds (leading to Theorem 2.4).All localized patterns are constructed in section 3, which begins with (another) preparatory section -section 3.1 -in which we set up the geometry of orbits 'jumping' between slow manifolds.The primary traveling 1-front patterns of Fig. 1a are constructed in section 3.2, the associated higher order 1-fronts of Figs.2a and 2b in section 3.3.Stationary patterns are considered in 3.4 -on 1-fronts -and 3.4 -on 2-fronts of spot and gap type as shown in Figs.1b, 1c and Fig. 2(c,d); various families of spatially periodic multi-front patterns -including the basic ones of Fig. 1d -are constructed in section 3.6.Section 4 begins with section 4.1 in which we show various numerically obtained patterns -some of them beyond the analysis of the present paper -and ends with discussion section 4.2.
Remark 1.1.While the original model (1.1) has 8 parameters -(Λ, Γ, R, K, E, M, N, P ) -(neglecting D B , D W which are represented by ε), rescaled model (1.5) has only 5 parameters -(a, Ψ, Φ, Ω, Θ).We will formulate our results by stipulating conditions on (a, Ψ, Φ, Ω, Θ) and refrain from giving a corresponding range for the original parameters.Moreover, we notice that we have implicitly assumed that α > 0, i.e. that EK > 1 (1.3).This is a technical assumption (and not unrealistic from ecological point of view), the case 0 < EK < 1 can be treated in a completely analogous way -see B.

The fast reduced problem
The fast reduced limit problem associated to (1.7) is a two-parameter family of planar systems that is obtained by taking the limit ε → 0 in (1.7), These planar systems can have up to 3 (families of) critical points (parameterized by (w 0 , q 0 )) given by, Clearly, (b 0 , w 0 ) represents the (homogeneous) bare soil state B(x, t) ≡ 0, the other two solutions correspond to uniform vegetation states and only exist for w 0 > 4/(1 + 4a).The critical points also determine 3 two-dimensional invariant (slow) manifolds, M 0 0 and M ± 0 , (2.5) A straightforward analysis yields that the critical points (b + , 0) are saddles for all c ∈ R and that the points (b 0 , p 0 ) = (0, 0) are saddles for all c as long as w 0 < 1/a.Therefore, we consider in this paper w 0 such that, so that (parts of) the manifolds M 0 0 and M + 0 are normally hyperbolic for all w 0 that satisfy (2.6) (and thus persist as ε becomes nonzero [39,41]); moreover, all stable and unstable manifolds W s,u (M 0 0 ) and W s,u (M + 0 ) are 3-dimensional.(In this paper, we do not consider the manifold M − 0 for several reasons: (i) it is not normally hyperbolic in the crucial case of stationary patterns (i.e. for c = 0, under the -natural -assumption that the water concentration w 0 does not become negative), (ii) critical points for the full system (1.7) that limit on M − 0 as ε → 0 cannot correspond to stable homogeneous states of PDE (1.5) -see section 2.3.)The manifolds W s,u (M 0 0 ) and W s,u (M + 0 ) are determined by the stable and unstable manifolds of (0, 0) and (b + , 0).By the (relatively) simple cubic nature of (2.3) we do have explicit control over these manifolds in the relevant cases that they collide, i.e. that there is a heteroclinic connection between (0, 0) and (b + , 0).Although this is a classical procedure -see [52] -we provide a brief sketch here.
We assume that a heteroclinic solution of (2.3) between (0, 0) and (b + , 0) can also be written as a solution of the first order equation where n is a free pre-factor (and we know that this assumption provides all possible heteroclinic connections).Taking the derivative (w.r.t.ξ) yields an equation for b ξξ that must equal (2.3) -that we write as Working out the details yield explicit expressions for n and c, Thus, for a given c, there is a heteroclinic connection between M 0 0 and M + 0 at the 'level' w 0 = w ± h (c) if w 0 solves (2.8).A direct calculation yields that c ± (w 0 ) are strictly monotonic function of w 0 with inverse (2.9) We conclude that for a given c, there may be 'parabolic' -by the relation between b and p (2.7 (where we recall that q = q 0 ∈ R is still a free parameter).See Lemma 3.2 for a further discussion and analysis (for instance on the allowed c-intervals for which the heteroclinic connections exist: w ± h (c) must satisfy (2.6)).

The slow reduced limit problems
The slow reduced limit problem is obtained by taking the limit ε → 0 in (2.1).It is a planar problem in (w, q), w XX = −Ψ + Φ + Ωb + Θb 2 w. (2.12) i.e. (2.12) describes the (slow) flow on the (slow) manifolds M 0 0 and M ± 0 (2.5).The flow on M 0 0 is linear, with critical point P 0 0 = (0, 0, Ψ/Φ, 0) ∈ M 0 0 of saddle type -that corresponds to the uniform bare soil state (B(x, t), W (x, t)) ≡ (0, Ψ/Φ) of (1.5) -that has the stable and unstable manifolds (on M 0 0 ) given by Figure 3: Numerical simulations of dynamics of the fast reduced system (1.7) for a = 1 4 and two choices of w 0 ∈ U a (2.6), both featuring a heteroclinic orbit between the saddle points (0, 0) and (b + (w 0 ), 0): Since we focus on orbits -patterns -that 'jump' between M 0 0 and M + 0 (in the limit ε → 0), we do not consider the flow on M − 0 but focus on (the flow on) M + 0 , where and we notice explicitly that B and C may be negative (since Ω may be negative).For w satisfying (2.6), we define, and conclude that the critical points P +,j 0 = (b + (w +,j 0 ), 0, w +,j 0 , 0) ∈ M + 0 are determined as solutions of the quadratic equation, (2.18) Thus, the points P +,j 0 exist for parameter combinations such that C 2 − 4AD > 0. There are 2 critical points if additionally C < 0 and D > 0 and only 1 if D < 0.
Clearly, the flow (2.15) is integrable, with Hamiltonian given by with, for ã = a + 1 4 , Hence, if non-degenerate, the critical points P +,j 0 are either centers -P +,c 0 -or saddles -P +,s 0 .Notice that, except the uniform bare soil state (0, Ψ/Φ), all critical points correspond to uniform vegetation states (B(x, t), W (x, t)) ≡ ( B, W ) in (1.5) -see section 2.3.In the case that there is only 1 critical point P + 0 ∈ M + 0 , it can either be of saddle or center type: where W > 0 is the solution of (2.18).We notice that the stable and unstable manifolds (restricted to M + 0 ) of the saddle point P +,s 0 ∈ M + 0 are represented by, Figure 4: Phase portrait of the unperturbed flow (2.15) on M + 0 for parameters (a, Ψ, Φ, Ω, Θ) such that (2.32) holds.
In the upcoming analysis, we will be especially interested in the case of 2 critical points P +,s 0 , P +,c 0 ∈ M + 0 , therefore we investigate this situation on some more detail.First, we introduce D SN and σ ≥ 0 by setting, so that the solutions of (2.18) are given by W = W SN ± σ = − C 2A ± σ.We rewrite (2.15) in terms of (a, A, C, D) Clearly, σ = 0 corresponds to the degenerate saddle-node case in which P +,s 0 and P +,c 0 merge, where we note that w SN 0 satisfies (2.6) for 0 < C 2 < A 2 (independent of a).In fact, we can consider the unfolding of the saddle-node bifurcation by the additional assumption that 0 < σ 1, (j = 1, 2), where the +-case represents the saddle P +,s 0 and the −-case the center P +,c 0 : w +,c 0 < w +,s 0 -see Fig. 4. In this parameter region, the slow reduced system (2.15) features a homoclinic orbit (w hom , q hom ) to P +,s 0 and a family of periodic solutions around the center point P +,c 0 (Fig. 4).
Remark 2.1.We conclude from (2.24) that the reduced slow flow on M + 0 is fully determined by the values of (a, A, B, D).Clearly, the (linear) mapping (Ψ, Φ, Ω, Θ) → (A, B, D) has a kernel: we can vary one of the parameters -for instance Φ -and determine (Ψ, Ω, Θ) such that this does not have an effect on the reduced flow (2.24) on M + 0 (by choosing (Ψ(Φ), Θ(Φ), Ω(Φ)) such that (A, B, D) are kept at a chosen value).We will make use of this possibility extensively in section 3.

Critical points and homogeneous background states
Since the critical points P j = (b j , p j , w j , q j ) of the full ε = 0 system (1.7) must have p j = q j = 0, their (b, w) coordinates are determined by the intersections of the b-and w-nullclines, where we recall that the b-nullcline determines the slow manifolds M 0 0 and M ± 0 -see Fig. 5. Hence, all critical points P j must correspond to critical points of the slow reduced flows on either one of the (unperturbed) slow manifolds This immediately implies that P 1 = P 0 0 = (0, 0, Ψ/Φ, 0) ∈ M 0 0 .The (potential) critical points on M − 0 can be determined completely analogously to P +,j 0 ∈ M + 0 in section 2.2 -the only difference is that the term +CW in (2.18) must be replaced by −CW.Thus, we conclude that there are two additional critical points P 2 and P 3 if C 2 − 4AD ≥ 0 (and that and C > 0, then P 2 = P −,1 , P 3 = P −,2 and both P −,j ∈ M − 0 ; • if D > 0 and C < 0, then P 2 = P +,1 , P 3 = P +,2 and both P +,j ∈ M + 0 . Naturally, the critical points P j correspond to homogeneous background states (B(x, t), W (x, t)) ≡ ( Bj , W j ) of the full PDE (1.5).In this paper, we focus on the existence of patterns in (1.5) and do not consider the stability of these patterns (which is the subject of work in progress).However, there is a strong relation between the local character of critical points P j in the spatial system (1.7) and their (in)stability as homogeneous background pattern in (1.5) -see for instance [17].Therefore, we may immediately conclude, • the bare soil state ( B, W ) = (0, Ψ/Φ) is stable as solution of (1.5) for Ψ/Φ < 1/a, i.e. as long as (0, Ψ/Φ) corresponds to a critical point on the normally hyperbolic part of M 0 0 (2.6); • background states ( B, W ) that correspond to critical points on M − 0 are unstable; • a background state ( B, W ) that corresponds to a center point on M + 0 is unstable; • a background state ( B, W ) that corresponds to a saddle point on M + 0 is stable as solution of (1.5) if one additional (technical) condition on the parameters of (1.5) is satisfied.
Of course this motivates our choice to study homoclinic and heteroclinic connections between the saddle points on M 0 0 and M + 0 in this paper.Remark 2.2.The singular perturbation point of view also immediately provides insight in the possible occurrence of a Turing bifurcation in (1.5).In the setting of (1.7) -with c = 0 -a Turing bifurcation corresponds to a reversible 1 : 1 resonance Hopf bifurcation [38], i.e. the case of a critical point with 2 colliding pairs of purely imaginary eigenvalues.By the slow/fast nature of the flow of (1.7), such a critical point cannot lay inside one of the 3 possible reduced slow manifolds M 0 0 , M − 0 or M + 0 (critical points not asymptotically close to the boundaries must have 2 O(ε) and 2 O(1) eigenvalues).Thus, critical points that may undergo a Turing/reversible 1 : 1 Hopf bifurcation have to be asymptotically close to the edge of M + ε where it approaches M − ε (where we note that we a priori do not claim that M − ε persists).Indeed, the bifurcation appears in that region -although we refrain from going into the details.See Fig. 18a for a thus found spatially periodic Turing pattern in (1.5).
Remark 2.3.By directly focusing on (2.27) -and thus by not following the path indicated by the singularly perturbed structure of (1.7) -the uniform vegetation background states can also be computed in a more straightforward way:

2.4
The slow flows of the ε = 0 system Condition (2.6) was chosen such that the points (0, 0, w, q) ∈ M 0 0 and (b + (w), 0, w, q) ∈ M + 0 are saddles for the fast reduced limit problem (2.3) (so that the associated background states may be stable as trivial, homogeneous, patterns of (1.5) -section 2.3).Thus, where (2.6) holds, M 0 0 and M + 0 are normally hyperbolic and they thus persist as M 0 ε and M + ε for ε = 0 [39,41].Clearly, M 0 0 is also invariant under the flow of the full system (1.7):M 0 ε = M 0 0 .Moreover, the flow on M 0 ε is only a slight -O(ε) -(linear) perturbation of the unperturbed flow (2.13) on M 0 0 -due to the (asymmetric) −εcq term.As a consequence, only the orientation of the (un)stable manifolds W s,u (P 0 The situation is very different for M + ε .A direct perturbation analysis yields, with  and b + (w) as defined in (2.4).Since we only consider situations in which there are critical points (of the full flow) P +,j on M + 0 , and thus on M + ε , we know (and use) that M + ε is determined uniquely.The slow flow on M + ε is given by (cf. (2.15)), with Thus, for c = 0 the flow on M + ε is a perturbed integrable planar system with 'nonlinear friction term' εcqρ 1 (w).In the case that there is only one critical point P +,s of saddle type on M + 0 -and thus on M + 0 -the impact of this term is asymptotically small.The situation is comparable to that of the flow on M 0 ε w.r.t. the flow on M 0 0 .The stable and unstable manifolds of P +,s restricted to the slow manifolds remain close: W u,s (P +,s )| M + ε is O(ε) close to W u,s (P +,s )| M + 0 (for O(1) values of (w, q)) and the span W u,s (P +,s ) ∪ W u,s (P +,s )| M + ε has becomes slightly asymmetric -cf.(2.22).This is drastically different in the case that there are 2 critical points P +,c -the centerand P +,s -the saddle -on M + ε .We deduce by classical dynamical system techniques -such as the Melnikov method (see for instance [34]) -the following (bifurcation) properties of (2.30), and thus of (2.3).
It is possible to (locally) get full analytical control over the set S per and its boundary manifolds R Hopf and R hom in (a, Ψ, Φ, Ω, Θ)-space by considering the unfolding of the saddle-node bifurcation on M + ε as in section 2.2 (cf.(2.26)).
Remark 2.8.A higher order perturbation analysis yields that the O(ε) corrections to R Hopf and R hom -and thus to S per -explicitly depend on c.
Remark 2.9.Of course one can also establish the persistence of periodic orbits of the slow reduced flow -as in Theorem 2.4 -under the assumption that there is only one critical point P +,c of center type on M + 0 , instead of focusing on the present case in which the reduced slow flow (2.12) has a homoclinic orbit on M + 0 (Theorem 2.4).Since we decided to focus on situations in which there is a saddle point on M + 0 -that is potentially stable as homogeneous background state in (1.5) (section 2.3) -we do not consider this possibility here.Note however that the analysis of this case is essentially the same as presented here.See also Remark 3.7.

Localized front patterns
In this section we use the slow-fast geometry of the phase space associated to (1.7) to establish a remarkably rich variety of localized vegetation patterns (potentially) exhibited by model (1.5).First, we consider various kinds of traveling and stationary 'invasion fronts' that connect the bare soil state to a uniform or an 'oscillating' vegetation state and their associated bifurcation structures (sections 3.2, 3.3 and 3.4), next we study stationary homoclinic 2-front spot and gap patterns (section 3.5) and finally spatially periodic multi-front (spot/gap) patterns (section 3.6).As starting point, we need to control the intersection of W u (P 0 ) and W s (M + ε ).
Remark 3.1.We start by considering localized patterns that correspond to orbits in W u (P 0 ), i.e. patterns that approach the bare soil state ( B, W ) = (0, Ψ/Φ) of (1.5) as x → −∞.In fact, the upcoming results on 1-fronts are all on orbits in (1.7) that connect P 0 ∈ M 0 ε either to a critical point or to a persisting periodic orbit in M + ε (Theorem 2.4): all constructed 1-fronts originate from the uniform bare soil state.The existence of 1-front patterns that approach ( B, W ) = (0, Ψ/Φ) as t → +∞ is embedded in these results through the application of the symmetry (2.2).
We know by (2.10) that W u (M 0 0 ) and W s (M + 0 ) intersect transversely -and thus that ) is 2-dimensional and at leading order (in ε) given by (2.10).Since W u (P ) exponentially close until its w-component reaches w + h (c) (2.9) at which it 'takes off' from M 0 ε to follow the fast flow along the 'parabolic' manifold given by (2.10), all at leading order in ε -see sections 2.1, 2.4.Since w, q only vary slowly (1.7), the (w, q)-components of the orbit W u (P 0 ) ∩ W s (M + ε ) remain constant at leading order during its fast jump: it 'touches down' on M + ε with (at leading order) the same (w, q)-coordinates (Remark 3.3).Therefore, we define the touch down curve T down (c) ⊂ M + ε as the set of touch down points of the orbits W u (P 0 ) ∩ W s (M + ε ) that take off from M 0 ε exponentially close to the intersection u ε ∩ {w = w + h (c)} (2.9), parameterized by c; it is at leading order (in ε) given by, In terms of the projected (w, q)-coordinates by which the dynamics on M + ε are described (2.30), T down (c) describes a smooth 1-dimensional manifold I down = {(w down (c), q down (c))} parameterized by c with boundaries (its endpoints): the family of base points of the Fenichel fibers of W u (P 0 ) ∩ W s (M + ε ) on M + ε -Remark 3.3; at leading order in ε, I down is a straight interval with endpoints determined by the bounds (2.6) on w = w + h (c).
Lemma 3.2.At leading order in ε, , 1 √ 2a ] → I down is bijective and Expression (2.9) a priori does not exclude the possibility that w + h has several extremums as function of c, in ) = 0.The proof -derivation -of this lemma thus requires some careful, but straightforward, analysis.We refrain from going into the details here.
We conclude this section by noticing that heteroclinic connections γ h (ξ) between P 0 ∈ M 0 ε and P +,s ∈ M + ε directly correspond to intersections I down ∩ W s (P +,s )| M + ε (Remark 3.3).However, the coordinates of this intersection determine c (through I down ), while W s (P +,s )| M + ε also varies as function of c.Moreover, by the perturbed integrable nature of the flow on M + ε (2.30), there can a priori be (countably) many intersections Thus, the analysis is more subtle and richer than (perhaps) expected -as we shall see in the upcoming sections.
Remark 3.3.We (for instance) refer to [18] for a more careful treatment of 'take off ' and 'touch down' points/manifolds.In fact, these points/manifolds correspond to base points of Fenichel fibers (that persist under perturbation by Fenichel's third Theorem [39,41]).By construction/definition, an orbit that touches down at a certain (touch down) point on a slow manifold is asymptotic to the orbit of the slow flow that has this point as initial condition.Therefore, if an orbit touches down on a stable manifold of a critical point on the slow manifold, it necessarily is asymptotic to this critical point.

Traveling 1-front patterns -primary orbits
Our first result -on the existence of primary heteroclinic orbits -can be described in terms of the slow reduced flow on M + 0 , or more precise, on intersections of the touch down manifold I down and the restricted stable manifold ) of the reduced slow flow (2.15) on M + 0 .However, it is a priori unclear whether such intersections may exist and how many of such intersections may occur: the many parameters of system (1.7) have a 'nontrivial' effect on I down and W s (P +,s )| M + 0 and thus on their relative positions.To obtain a better insight in this, we 'freeze' the flow of (2.15) by fixing a, A, C, D at certain values.Since B + aΘ = D + (a + 1 4 )A (2.17), this indeed fixes all coefficients of the reduced slow flow (2.15) on M + 0 .At the same time, this leaves a 1-parameter freedom in the parameters Φ, Ψ, Ω, Θ. Defining, we see that for all Φ, the choices yield identical slow reduced flows (2.15).On the other hand, the (leading order) interval I down clearly varies as function of Φ, Note that for χ > 0, the intersection of I down with the w-axis can be varied between the critical w values 4/(1 + 4a) and 1/a by increasing Φ from a(1 + 4a)χ to ∞.In fact, χ > 0 necessarily holds in case there are 2 critical points on M + ε (since in that case A, D > 0, C < 0), while χ can also chosen to be positive in the case that there is only 1 critical point on M + ε .Thus, by choosing Ψ, Ω, Θ as in (3.3) and varying Φ we can control Theorem 3.4.Let P +,s = (b + (w +,s ), 0, w +,s , 0) ∈ M + 0 be a critical point of (1.7) that is a saddle point for the slow reduced flow (2.15) on M + 0 , and consider the touch down manifold I down at leading order given in Lemma 3.2 and the restricted stable manifold W s (P +,s )| M + 0 of the reduced slow flow (2.15).If there is a non-degenerate intersection point ( wprim,0 , qprim,0 ) ∈ I down ∩ W s (P +,s )| M + 0 , then, for ε sufficiently small, there exists for c = c prim a primary heteroclinic orbit γ prim (ξ) = (w prim (ξ), p prim (ξ), b prim (ξ), q prim (ξ)) ⊂ W u (P 0 )∩W s (P +,s ) of (1.7) connecting P 0 ∈ M 0 ε to P +,s ∈ M + ε -where c prim = c prim,0 + O(ε) and c prim,0 is the unique solution of w + h (c) = wprim,0 (2.9).Departing from P 0 (and at leading order in ε), γ prim (ξ) first follows u 0 ⊂ M 0 0 (2.14) until it reaches the take off point (0, 0, wprim , qprim ) from which it jumps off from M 0 0 and follows the fast flow along ; from there, it follows W s (P +,s )| M + 0 towards P +,s .Moreover, there is an open region S 1 s−prim in (a, Ψ, Φ, Ω, Θ) parameter space for which I down and W s (P +,s )| M + 0 intersect transversely; however, there is at most one intersection ( wprim , qprim ) ∈ I down ∩ W s (P +,s )| M + 0 and thus at most one primary heteroclinic orbit γ prim (ξ); in fact, this is the only possible heteroclinic orbit between P 0 and P +,s ; • if there are two critical points on M + ε , the center P +,c and saddle P +,s , i.e. if C 2 − 4AD > 0, C < 0, D > 0, then there are open regions S 1 cs−prim , respectively S 2 cs−prim , in (a, Ψ, Φ, Ω, Θ) parameter space for which I down and W s (P +,s )| M + 0 have 1, resp.2, (transversal) intersections, so that there can be (up to) 2 distinct primary heteroclinic orbit γ j prim (ξ) that travel with different speeds, i.e.
In the case of 2 critical points on M + ε , we shall see that the primary orbits may only be the first of many 'higher order' heteroclinic orbits -see section 3.3.We refer to Fig. 6 for sketches of the constructions in M 0 ε that yield the primary heteroclinic orbits γ prim (ξ) and to Figs. 1a, 13 and 14b for the associated -numerically obtainedprimary 1-front patterns in (1.5)-see especially Fig. 13b in which the the slow-fast-slow structure of a (numerically obtained) heteroclinic front solutions of (1.5) is exhibited by its projection in the 3-dimensional (b, w, q)-subspace of the 4-dimensional phase space associated to (1.7).
Proof.The existence of the heteroclinic orbit γ prim (ξ) follows by construction -Remark 3.3 -from an intersection of I down and W s (P +,s )| M + ε .Thus, we first need to show that a (non-degenerate) intersection By the assumption that ( wprim,0 , qprim,0 ) ∈ I down ∩ W s (P +,s )| M + 0 is a non-degenerate intersection point, we know that the intersection is transversal, and thus that W s (P ) with c = c -also intersects I down transversally as c is varied around c prim,0 in an O(1) fashion.Thus, for c sufficiently (but Since the flows of (2.15) and (2.30) are O(ε) close, we know that ( wprim,0 , qprim,0 ) − ( wdown (c i ), qdown (c i )) = O(ε), which implies that c i (c) = c prim,0 + O(ε).Hence, the O(1) variation of c through c prim,0 yields at leading order (in ε) a horizontal line c i (c) ≡ c prim,0 : there must be a unique intersection c i (c * ) = c * , and thus, by construction, we freeze the flow of (2.15) with A, C, D such that χ > 0 (3.2) and define Φ = Φ +,s such that Ψ(Φ)/Φ = 1/a − χ/Φ = w +,s , the w-coordinate of the saddle P +,s on M + 0 -see (3.3), (3.4).Since q is an increasing function of w on I down and W s (P +,s )| M + ε is decreasing near P +,s -see ) is given by a (strictly) decreasing function q +,s | M + 0 (w) for all w ∈ (4/(1 + 4a), 1/a) since it cannot have extremums: ) and thus to critical points of (2.15).By assumption, there are no critical points besides P +,s , which yields that there indeed can be maximally one intersection To control the case with a center P +,c and saddle P +,s on M + 0 , we again consider the unfolded saddle-node case of Lemma 2.6 and define Φ = Φ +,c such that Ψ(Φ)/Φ = 1/a − χ/Φ = w +,c , the w-coordinate of the center P +,c .The level set {H + 0 (w, q) = H +,s 0 } forms a small (w.r.t. the unfolding parameter σ) homoclinic loop around P +,c that intersects I down (transversally) in two points ( wj,0 prim , qj,0 prim ), j = 1, 2 -see Fig. for values of Φ > Φ +,s (as defined above). 2
The homoclinic orbit (w hom,0 (X), q hom,0 (X)) of (2.15) typically breaks open under the perturbed flow of (2.30), and W s (P +,s )| M + ε either spirals inwards in backwards 'time', i.e. as ξ → −∞, or not.In the former case, there will be (typically many) further intersections (at leading order in ε), we may expect further heteroclinic connections γ h,j in (1.7) connecting P 0 ∈ M 0 ε to P +,s ∈ M + ε beyond the primary orbits γ prim (ξ) established in Theorem 3.4.In fact, it follows directly that γ 1 prim (ξ) and γ 2 prim (ξ) are the only heteroclinic orbits between P 0 and P +,s if (3.5) does not hold.If (3.5) does hold, the (spiraling part of) W s (P +,s )| M + 0 clearly must limit -for ξ → −∞ -on either the center P +,c or, if (a, Ψ, Φ, Ω, Θ) ∈ S per , on the persistent periodic solution (b p,ε (X), p p,ε (X), w p,ε (X), q p,ε (X)) ⊂ M + ε (Theorem 2.4).Therefore, we first formulate a result on the existence of heteroclinic connections between P 0 ∈ M 0 ε and (b p,ε (X), p p,ε (X), w p,ε (X), q p,ε (X)) ⊂ M + ε .Like in Theorem 3.4, this can be done in terms of the unperturbed flow in M + 0 .Theorem 3.5.Assume that (2.32) holds and that (a, Ψ, Φ, Ω, Θ) ∈ S per .Let (w p,0 (X), q p,0 (X)) 19) be the periodic solution of (2.15) that persists (on cs−prim defined in Theorem 3.4 -such that there are 2 (non-degenerate) intersection points ( wj Notice that this result is independent of condition (3.5), i.e.Theorem 3.5 holds independent of the sign of ∆H hom .Moreover, we could formulate similar limiting result concerning heteroclinic 1-front connections between P 0 ∈ M 0 ε and P +,c ∈ M + ε for (a, Ψ, Φ, Ω, Θ) on a certain co-dimension 1 manifold.Since the background state associated to P +,c cannot be stable -section 2.3 -we refrain from going into the details.'wraps around' the persistent periodic solution (b p,ε (X), p p,ε (X), w p,ε (X), q p,ε (X)) (Theorem 3.5) -in backwards time -and since I down intersects this orbit in 2 points (by assumption), there are two countable sets of intersections Proof.The proof goes exactly along the lines of that of Theorem 3.4.2 Theorem 3.5 provides the foundation for a result on the existence of multiple -in fact countably many -distinct traveling 1-front connections between P 0 ∈ M 0 ε and P +,j ∈ M 0 ε for an open set in parameter space -see also the sketches in Figs.2a and 2b.Theorem 3.6.Assume that the conditions of Theorem 3.5 hold and let (a, Ψ, Φ, Ω, Θ) ∈ S h−p .If c j h−p and c j prim have the same sign (for either j = 1 or 2) and if (3.5) holds for c of this sign, then -for ε sufficiently small -there are countably many distinct heteroclinic orbits , and, As in the proofs of Theorems 2.4 and 3.4, we can verify that there indeed are open regions in (a, Ψ, Φ, Ω, Θ)-space for which c j h−p and c j prim have the same sign (for either j = 1, 2 or for both) and such that (3.5) holds, by considering the unfolded saddle-node case of Lemma 2.6.In fact, we know from Lemma 3.2 that c changes sign as the w-coordinate of the intersection point on I down ∩ W s (P +,s )| M + 0 passes through 9/(2 + 9a).Thus (and for instance), all 4 values c j h−p and c j prim , j = 1, 2, must have the same sign as the entire homoclinic orbit spanned by W s (P +,s )| M + 0 either is to the left or to the right of w = 9/(2 + 9a) -more precise, if either (w h,0 , w +,s ) ⊂ (4/(1 + 4a), 9/(2 + 9a)) or (w h,0 , w +,s ) ⊂ (9/(2 + 9a), 1/a) (cf.(2.32)).Note that it follows from (2.25) that w SN 0 = 9/(2 + 9a) implies that C 2 = A /9 (independent of a), so that we can indeed move the homoclinic loop associated to the unfolded saddle-node -i.e.σ 1 as in (2.26) -through w = 9/(2 + 9a) by increasing C 2 ∈ (0, A 2 ) through A 2 /9.On the other hand, it is certainly also possible that c j h−p and c j prim do not have the same sign.Hence, apart from the PDE point of view -from which it is natural to consider stationary patterns -this gives us an additional motivation to study the sign-changing stationary case c = 0 in more detail, as we will briefly do in Remark 3.10 in section 3.4.
Proof.We only consider the case j = 1, i.e. we assume that c 1 h−p and c 1 prim have the same sign and that (3.5) holds for c = c 1 h−p , c prim .The proof for j = 2 goes exactly along the same lines.
For c = c 1 prim , W s (P +,s )| M + ε by assumption spirals inwards in backwards 'time' and 'wraps around' the (perturbed) periodic orbit (w p,ε (X), q p,ε (X)) on (the projection of) M + ε -see Fig. 7. Since (a, Ψ, Φ, Ω, Θ) ∈ S h−p , W s (P +,s )| M + ε must intersect I down countably many times.We define ( w1,1 i , q1,1 i ) as the next intersection of W s (P +,s )| M + ε (c 1 prim ) with I down beyond the 2 primary intersection points: it is the first non-primary intersection point and has q1,1 i > 0. As before, ( w1,1 i , q1,1 i ) ∈ I down = {( wdown (c), qdown (c))} determines the value c 1,1 i through wdown (c) = w1,1 i -where we know that c 1,1 i < c 1 prim since the w-component of I down is a monotonically increasing function of q (Lemma 3.2).Since the perturbation term in (2.30) ) change in the flow of (2.30), hence for all c O(ε) close to c 1 prim , the first non-primary intersection of W s (P +,s )| M + ε (c) and I down -denoted by ( w1,1 ).Thus, the speed c 1,1 i (c) associated to this intersectionby wdown (c) = w1,1 i (c) -must also be O(ε 2 ) close to c 1,1 i .The situation is therefore similar to that in the proof of Theorem 3.4: an O(ε) variation of c around c 1,1 i in (2.30) yields only an O(ε 2 ) change in the c-coordinate associated the first non-primary intersection I down ∩ W s (P +,s )| M + ε (c) so that there must be an unique c = c 1,1 h such that . This establishes the existence of the first non-primary heteroclinic 1-front orbit We can now iteratively consider the first intersection in backwards 'time' -denoted by ( w1,2 h .Completely analogous to the above arguments, we deduce the existence of an unique c = c 1,2 h such that )), which establishes the existence of the next non-primary Theorem 2.4 holds independent of c, which implies that W s (P +,s )| M + ε (c) wraps around the periodic orbit periodic orbit (w p,ε (X), q p,ε (X)) of (2.30) (in backwards 'time') for all c with the same sign as c 1 h−p and c 1 prim (cf.Fig. 7).Thus, there must be countably many heteroclinic orbits γ 1,k h (ξ) -every W s (P +,s )| M + ε (c) intersects I down countably many times -and the associated speeds c 1,k h must all be between c 1 h−p and c 1 prim .Moreover, the decreasing sequence {c 1,k h } ∞ k=1 must have a limit that cannot differ from By establishing the existence of countably many distinct heteroclinic connections between P 0 and P +,s , Theorem 3.6 in a sense considers (one of) the most complex case(s), which is quite far removed from situations in which there are no such connections.To obtain insight in the bifurcations that occur 'in between', we can again freeze the reduced slow flow and vary I down = I down (Φ) by increasing Φ (from a(1 + 4a)χ to ∞ (3.2), (3.3), (3.3)).We consider the most simple case and assume that the homoclinic orbit of the frozen flow lies entirely in the w-region (4/(1+4a), 9/(2+9a)) -so that all c j,k h 's of Theorem 3.6 are positive -and that I down (Φ) ∩ W s (P +,s )| M + 0 = ∅ at Φ = a(1 + 4a)χ (this can easily be achieved by the unfolded saddle-node approach).As Φ increases, I down (Φ) becomes steeper and the intersection I down (Φ) ∩ {q = 0} moves over the entire interval determined by (2.6), i.e. from 4/(1 + 4a) to 1/a.Thus, I down (Φ) moves through the homoclinic loop spanned by W s (P +,s )| M + and through the enclosed persistent periodic orbit established by Theorem 2.4.We tune the parameters such that during the passage of the latter, (3.5) holds and (a, Ψ, Φ, Ω, Θ) ∈ S h−p , i.e. that Theorem 3.6 can be applied -which is also possible.It should be noticed that although the reduced flow (2.15) is frozen, this is not the case for the perturbed flow (2.30), since ρ 1 (w) varies with Θ (2.31) and Θ = Θ(Φ) (3.3), thus the persistent periodic orbit established in Theorem 2.4 is not frozen, but also varies with Φ -this is represented in the sketches of Fig. 8 by the decreasing size of the limiting periodic orbit on M + ε .
Fig. 8 exhibits sketches of 3 configurations of I down (Φ) and W s (P +,s )| M + 0 for increasing Φ, the associated bifurcation scenario is sketched in Fig. 9.In Fig. 8a, Φ has already passed through the first bifurcation value Φ prim at which the first 2 primary heteroclinic orbits γ j prim (ξ), j = 1, 2 of Theorem 3.4 are created, and through a second one, Φ b SN,1 -O(ε) close Φ prim -at which the first 2 secondary orbits appear.This bifurcation is followed by countably subsequent saddle-node bifurcations until Φ reaches Φ b per at which the 2 limiting heteroclinic orbits between P 0 and the persistent periodic solution of Theorem 3.5 appear and we enter into the realm of Theorem 3.6.These orbits next disappear at Φ e per , Fig. 8(b) is similar to Fig. 7 and represents the 2 countable families of heteroclinic orbits that exist for Φ ∈ (Φ b per , Φ e per ) (Theorem 3.6).All these orbits step-by-step disappear in pairs as Φ is increased further: Fig. 8c shows the situation with only 5 left -4 of these will disappear just before Φ reaches Φ +,s at which I down (Φ) passes through P +,s .
e e r J e r b d Z a 8 7 K Z n b h F 6 z 3 L x r o j 4 s = < / l a t e x i t > P +,s < l a t e x i t s h a 1 _ b a s e 6 4 = " Z p p j m f R R L w j T B V t t P e q W i n G R G i 0 We refrain from giving all rigorous details on which the above sketched scenario is based -this is in essence a matter of following the lines set out in the proofs of the preceding results.Moreover, we also refrain from working out all possible alternative bifurcation scenarios that may occur -there are many (sub)cases to consider, some more simple, others more complex than that of Fig. 9. Nevertheless, we do briefly come back to this in the upcoming section -where we consider stationary, sign-changing, case c = 0 case.Remark 3.7.As in Remark 2.9, we note that a result like Theorem 3.5 on the existence of heteroclinic connections between P 0 and a periodic orbit on M + ε can also be established under the assumption that there is only one critical point P +,c of center type on M + 0 .Similar remarks can be made about the upcoming Theorems 3.9, 3.12 and Corollary 3.16.We note -also as in Remark 2.9 -that the analysis of these additional cases is essentially the same as already presented.

Stationary 1-front patterns
In this section, we construct stationary heteroclinic 1-front patterns that are similar to those constructed in Theorems 3.4 and 3.5.We immediately note that if the reduced flow on M + 0 has an unperturbed homoclinic loop W s (P +,s )| M + 0 ∩ W s (P +,s )| M + 0 -as in Fig. 4 -that it persists as homoclinic solution of (1.7) on M + ε for ε = 0since (1.7) with c = 0 is a reversible system (see also Theorem 2.4).Thus, we a priori deduce that there cannot be any further non-primary heteroclinic 1-front connections between P 0 ∈ M 0 ε and P +,s ∈ M + ε as those of Theorem 3.6 for c = 0 (see however also Remark 3.10 for a result similar to Theorem 3.6).In the subsequent sections, we will proceed to construct homoclinic and periodic multi-front patterns -i.e.solutions of (1.7) that jump up and down between M 0 ε and M + ε -and show that there is a richness in these kinds of patterns similar to that of Theorem 3.6.
As in the previous sections, we approach the bifurcation analysis by freezing the flow on M + ε .Thus, we choose Ψ, Ω, Θ as in (3.3) and vary Φ.We know by Lemma 3.2 and (3.4) that for c = 0, the touch down point of W u (P 0 ) ∩ W s (M + ε ) is represented by a vertical line/half line J s−d in the (w, q)-plane (at leading order in ε).Clearly, a 1-front connection between P 0 and P +,s corresponds to those values of Φ for which Theorem 3.8.Let c = 0 and let ε be sufficiently small.Then, there is a co-dimension 1 set R s−1f in (a, Φ, Ψ, Ω, Θ)space for which (stationary) 1-front heteroclinic orbits γ s−1f (ξ) ⊂ W u (P 0 ) ∩ W s (P +,s ) exists in (1.7).More precise, let Ψ, Ω, Θ as in (3.3), then: (A) If P +,s = (w +,s , 0) is the only critical point on (the projection of ) • if χ > 0, then there is a unique value Φ s−1f such that there is a 1-front heteroclinic orbit γ s−1f (ξ) ⊂ W u (P 0 ) ∩ W s (P +,s ) in (1.7); • if χ < 0, w +,s < w s−d = 9/(2 + 9a) (3.6) and there are Φ such that H + 0 (w s−d (Φ), q s−d (Φ)) < H +,s 0 (2.11), (2.22), then there are 2 values Φ j s−1f , j = 1, 2 for which 1-front heteroclinic orbits γ s−1f (ξ) ⊂ W u (P 0 ) ∩ W s (P +,s ) exist in (1.7); • if χ < 0 and either one of the above additional conditions does not hold, then there is no such stationary 1-front orbit.(B) If there are two critical points on M + ε , the center P +,c and saddle < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 q B 0 6 < l a t e x i t s h a 1 _ b a s e 6 4 = " b 9 c o r < l a t e x i t s h a 1 _ b a s e 6 4 = " j l d n X P T / 4 B p N P t 6 b 6 < l a t e x i t s h a 1 _ b a s e 6 4 = " U / J 1 7 w P s s 4 s p D n 7 5 U x a s  < l a t e x i t s h a 1 _ b a s e 6 4 = " Z T 4 e v F c V W n z P 9 z l i + L o s 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " c Z 7 p n h M + 5 r 6 2 K P S j r O 1 C T r U y g i 5 g d D l K 5 S p P y d i 7 E k 5 8 x z d 6 W E 1 k X + 9 V P z P 6 0 f K P b d j 5 o e R o j 6 Z L 3 I j j l S A 0 h z Q i A l K F J 9 p g o l g + q + I T L D A R O m 0 S l k I F y l O v 0 9 e J J 3 j W r 1 R a 1 y f V J o 3 e R x F 2 I c D q E I d z q A J V 9 C C N h C 4 h 0 d 4 h h f j w X g y X o 2 3 e W v B y G f 2 4 B e M 9 y + Q k Z R h < / l a t e x i t > +,s < l a t e x i t s h a 1 _ b a s e 6 4 = " X p U V p F q A J 8 u Figure 9: A sketch of the bifurcation scenario as function of Φ representing the appearance in a saddle-node bifurcation of the 2 primary heteroclinic 1-front orbits γ j prim (ξ), j = 1, 2, of Theorem 3.4, followed by further saddle-node bifurcations leading to the situation governed by Theorem 3.6 in which countably many 1-front orbits exist; these orbits subsequently disappear in another cascade of saddle-node bifurcations eventually leaving only one (primary) 1-front orbit behind (Theorem 3.4).The relative configurations of I down (Φ) and W s (P +,s )| M + 0 sketched in Fig. 8 occur at the Φ-values indicated by the vertical (a), (b) and (c) lines.
We refer to Fig. 14a for an example of a numerical simulation of (1.5) exhibiting a stationary 1-front pattern.Moreover, we notice that -by symmetry (2.2) of (1.7) with c = 0 -the heteroclinic orbit γ s−1f (ξ) ⊂ W u (P 0 ) ∩ W s (P +,s ) has a counterpart ⊂ W u (P +,s ) ∩ W s (P 0 ), i.e. an orbit from P +,s to P .Together, these obits form a heteroclinic cycle between the saddles P +,s to P 0 .

< l a t e x i t s h a 1 _ b a s e 6 4 = " B 9 f z 2 0 T U R y b K y w g T t v 9 I 5 d n I / B M = " >
w m e H s + + V 5 0 j y p O q d V 9 8 a t 1 N w 8 j i L s w T 4 c g g P n U I N r q E M D C H B 4 h G d 4 s e 6 t J + v V e p u 1 F q x 8 Z h d + w X r / A h I 0 j 2 4 = < / l a t e x i t > 1,2 s-1f < l a t e x i t s h a 1 _ b a s e 6 4 = " V n j K / q S 9 a m / T 1 p J W z G y D X 9 D e v w A X d J a / < / l a t e x i t > w h,0 < l a t e x i t s h a 1 _ b a s e 6 4 = " X 9 c d P w l y 5 t Z u j 5 k 3 e e r J e r b d Z a 8 7 K Z n b h F 6 z 3 L x r o j 4 s = < / l a t e x i t > P +,c < l a t e x i t s h a 1 _ b a s e 6 4 = " j E x n v + T R l U h 0 p j 7 4 l 0 e e r J e r b d Z a 8 7 K Z n b h F 6 z 3 L w K I j 3 s = < / l a t e x i t > 1 s-1f < l a t e x i t s h a 1 _ b a s e 6 4 = " c n k z x s q x s s L u a k u f E g p a g z o e 9 c g = "

< l a t e x i t s h a 1 _ b a s e 6 4 = " g W y k r E t w 8 c D i B u j t f M O C T t H + g r 0 = " >
2 q P 2 r L 1 p 7 + P W g p b P b I I J a B / f R G q a O g = = < / l a t e x i t > 9 2 + 9a < l a t e x i t s h a 1 _ b a s e 6 4 = " C b l 0 u 8 f w E p y 5 S e < / l a t e x i t > P +,s < l a t e x i t s h a 1 _ b a s e 6 4 = " B 9 f z 2 0 T U R y b K y w g T t v 9 I 5
Remark 3.10.Together, Theorems 3.8 and 3.9 provide the possibility to establish a result similar to that of Theorem 3.6 in the case that c j h−p and c j prim do not have the same signs.Assume we have -for a certain parameter combination (a, Ψ, Φ, Ω, Θ) such that (2.32) holds -that c 1 h−p < 0 < c 1 prim .This implies that the point J s−d (Φ) on T down -Lemma 3.2 -must lie between the intersections of T down with the unperturbed homoclinic orbit that determines c 1 prim > 0 (Theorem 3.8) and the persisting periodic orbit that determines c 1 h−p < 0 (Theorem 3.9).Thus, J s−d (Φ) determines a level set H + 0 (w, q) = H + 0 (w s−d (Φ), q s−d (Φ)) = H ∈ (H +,c 0 , H +,c 0 ) and we know by Theorem 3.9 that Φ = Φ 2 p−2f ( H): (w s−d (Φ), q s−d (Φ)) is the touchdown point of a heteroclinic orbit between the bare soil state and the persisting (stationary) periodic orbit determined by the level set H + 0 (w, q) = H.It then follows by arguments similar to those in the proof of Theorem 3.6 that there are countably many c-values 0 < ... < c 1,k h < ... < c 1,1 h < c 1,0 h = c 1 prim with c 1,k h ↓ 0 for k → ∞ for which non-primary heteroclinic connections between P 0 and P +,s exist -as in Theorem 3.6.The main difference with Theorem 3.6 is that c = 0 determines a stationary orbit and not an attracting one: for c slightly above c = 0, the unstable manifold W s (P +,s )| M + ε only spirals inwards very weakly (in backward time).As a consequence, the number of intersections W s (P +,s )| M + ε ∩ T down with H + 0 (w, q) > H increases (without bound) as c ↓ 0.

Stationary homoclinic 2-front patterns: vegetation spots and gaps
In this section we construct stationary 2-front patterns that correspond to vegetation spots or vegetation gaps -the latter sometimes also interpreted as fairy circles.These patterns are observed in nature and appear as stable patterns in simulations of (1.1)/(1.5)-see [77] and Figs.14c and 15a.The patterns/orbits to be constructed are symmetric with respect to the reversibility symmetry in (1.5) that persists as (2.2) into (1.7)-with c = 0.As a consequence, we may expect that these patterns are generic, in the sense that they exist in open regions within parameter spacesee for instance [17].Notice that this is unlike the stationary -but non-symmetric -1-front patterns of the previous section, that only exist for (a, Φ, Ψ, Ω, Θ) ∈ R s−1f , an explicitly determined co-dimension 1 manifold (Theorem 3.8).
We first consider the (localized) spots: localized vegetated regions embedded within bare soil.Thus, these spots correspond to solutions of (1.7) that are homoclinic to the bare soil state P 0 .Singularly perturbed models of the type (1.7), can have homoclinic (pulse) solutions of various types.The localized vegetation (spot) patterns constructed in [60,3] in the context of the extended Klausmeier model -also called generalized Klausmeier-Gray-Scott model [69] -make a fast excursion away from the slow manifold that contains the critical point associated to the bare soil state following a homoclinic solution of the fast reduced system.As a consequence these spots are 'narrow', their size scales with ε.Although such pulses are also exhibited by the present model -see Fig. 15b and Remark 3.14we focus here on 2-front patterns, i.e. orbits homoclinic to P 0 that 'jump' from M 0 ε to M + ε , follow the slow flow on M + ε over an O(1) distance and jump back again -by its second fast reduced heteroclinic front -to M 0 ε .Since the extended Klausmeier model only has one slow manifold [60,3], such orbits cannot exist in that model.Moreover, these patterns have a well-defined width that does not decrease to 0 as ε ↓ 0, a property that appears to be natural in observed ecosystems [77].
The construction of the most simple -primary, cf.section 3.2 -singular 'skeleton structure' Γ s−2f ⊂ R 4 -the phase space of (1.7) -of a stationary homoclinic 2-front orbit γ s−2f (ξ) to P 0 of (1.7) is relatively straightforward (but somewhat involved/technical).Since the homoclinic orbit γ s−2f (ξ) ⊂ W u (P 0 ), it follows u along M 0 ε , takes off from M 0 ε by following the fast reduced flow and touches down on M + ε near an element of I down -see section 3.1.In fact, since we consider stationary spots, the touch down point is near a point J s−d ∈ J s−d (3.6) -where we note that γ s−2f (ξ) cannot exactly touch down on I down /J s−d since γ s−2f (ξ) / ∈ W s (M + ε ) -see the proof of (upcoming) Theorem 3.11.The touch down point J s−d determines (at leading order) a level set H s−d of the Hamiltonian H + 0 (w, q) (2.19) of the slow reduced flow (2.15) on M + ε : . As long as it remains (exponentially) close to M + ε , the homoclinic orbit-to-be-constructed remains asymptotically close to the level set H + 0 (w, q) = H s−d .This construction provides the first half of skeleton Γ s−2f , the second part follows by the symmetry (2.2) -with c = 0. Completely analogous to I down , one can define I off as the points on M + ε that determine the evolution of orbits in W s (P 0 ) ∩ W u (M + ε ) after their jump from s ⊂ M 0 ε through the fast field in backwards time.In fact, it follows by the symmetry (2.2) that I off and its stationary counterpart J s−o correspond exactly to the reflections of I down and J s−d with respect to the q-axis -where we consider I down /J s−d and I down /J s−d within the (projected) 2-dimensional representation of M + ε as in Lemma 3.2 and in (3.6).Thus, for a given Φ, the take off point J s−d is given by J s−o (Φ) = (w s−d (Φ), −q s−d (Φ)); this point also lies on the level set H s−d -since H + 0 (w, q) (of course) also is symmetric in q → −q.
We define the region S s−2f in (a, Ψ, Φ, Ω, Θ)-space for which the point J s−d (and thus J s−o ) can be constructed (as above) and there is a solution of the slow reduced flow (2.15) on M + ε that connects J s−d to J s−o -so that Γ s−2f indeed exists as closed singular 'loop'.Obviously, S s−2f = ∅ -see also Fig. 11 -however, the fact that both J s−d , J s−o ∈ {H + (w, q) = H s−d } does not necessarily imply that (a, Ψ, Φ, Ω, Θ) ∈ S s−2f for all values for which these points exist on M + ε .For instance, in the case that the saddle P +,s is the only critical point on M + ε , J s−d and J s−o are not connected by a solution of (2.15) if H s−d > H +,s 0 -the value of H + 0 (w, q) at P +,s -see Fig. 11a.Moreover, if H s−d < H +,s 0 there is an additional condition on Φ that is determined by the relative positions of w +,s (the w-coordinate of P +,s ), 9/(2 + 9a) (the w-coordinate of J s−d /J s−o ) and Ψ/Φ (the w-coordinate of the bare soil state associated to P 0 ).Here, we refrain from working out the full 'bookkeeping' details by which (the boundary of) S s−2f is determined -see also a further brief discussion following Theorem 3.11.We refer to Fig. 11a for a case w +,s < 9/(2 + 9a) < Ψ/Φ (and implicitly χ > 0) for which J s−d can only be connected to 2 a(2 + 9a)χ (3.7) -since we need that q s−d (Φ) < 0. Notice that the sketch in Fig. 11a in principle also covers a (sub)case of the situation with two critical points P +,c and P +,s on M + ε and that Fig. 11b considers the case 9/(2 + 9a) < w +,s < Ψ/Φ for this situation.Clearly, there are no further restrictions on Φ if H s−d < H +,s 0 if there are two critical points P +,c and P +,s on M + ε , since the orbits on the level set associated to H s−d (Φ) are periodic, while one again has to impose Φ < Φ s−t to have a connection between J s−d and J s−o for level sets outside the homoclinic loop, i.e. for H s−d > H +,s 0 .Finally, we note that the skeleton structure Γ s−2f can in principle also be constructed for (a, Ψ, Φ, Ω, Θ) such that there is no critical point on M + ε , or only one critical point that is not a saddle but a center (in the limit ε → 0).Summarizing, the (open) region S s−2f is defined such that for parameter combinations (a, Ψ, Φ, Ω, Θ) ∈ S s−2f , a singular skeleton Γ s−2f ⊂ R 4 can be constructed as above.In the limit ε → 0, Γ s−2f is spanned by a piece of u ⊂ M 0 0 from P up to the (ε → 0 limit of the) take off point from M 0 0 (that has the same (w, q)-coordinates as J s−d in the limit ε → 0), the jump through the fast field along (a piece of) W u (M 0 0 ) ∩ W s (M + 0 ) (2.10) towards the ε → 0 limit of the (projected) touch down point J s−d ∈ M + 0 , the connection along M + 0 to J s−d by the slow reduced flow (on the level set {H + 0 (w, q) = H s−d }) up to (the ε → 0 limit of) the take off point J s−o , followed by a fast jump backwards along (a piece of) W s (M 0 0 ) ∩ W u (M + 0 ) to M 0 0 and a final piece of s (up to P 0 ) -see Figs.11a and 11b for 2 sketches of projections of Γ s−2f on M + 0 that skip both jumps through the fast field.The proof of the persistence of Γ s−2f for ε = 0 relies heavily on the reversibility symmetry of (1.7) with c = 0 (2.2).Theorem 3.11.Let (a, Ψ, Φ, Ω, Θ) ∈ S s−2f and Γ s−2f ⊂ R 4 be the singular skeleton constructed above.Then, there is for ε > 0 sufficiently small a symmetric homoclinic 2-front orbit γ s−2f (ξ) = (b s−2f (ξ), p s−2f (ξ), w s−2f (ξ), q s−2f (ξ)) ⊂ W u (P 0 ) ∩ W s (P 0 ) of (1.7) with c = 0 that merges with Γ s−2f as ε ↓ 0. The associated stationary pattern (B(x, t), W (x, t)) = (b s−2f (x), w s−2f (x)) in (1.5) represents a stationary localized vegetation spot embedded in bare soil.
We refer to Figs. 1b and 14c for numerical observations of these 2-front spot patterns.In Fig. 16, the projection of the 2-front orbit on the (w, b)-plane is given; it clearly shows the slow-fast-slow-fast-slow nature of the pattern: it first follows M 0 ε (slowly), jumps to M + ε , follows the slow flow on M + ε , jumps back to M 0 ε and slowly returns to P 0 .
To get some insight in the boundaries of S s−2f -and thus is the bifurcations of γ s−2f (ξ) -we can (as usual) 'freeze' the flow on M + 0 by choosing Ψ, Ω, Θ as in (3.3) and vary Φ.We need to be aware though that u,s and P 0 do vary with Φ (i.e. they are not frozen).In the situation sketched in Fig. 11a -thus with w +,s < 9/(2 + 9a) < Ψ/Φ and χ > 0 -we see that the distance between the 2 fronts of γ s−2f (ξ) approaches ∞ as Φ ↓ Φ s−1f < Φ s−t as defined in Theorems 3.8 and 3.9: γ s−2f (ξ) obtains the character of the superposition of the 1-front heteroclinic orbit γ s−1f (ξ) of Theorem 3.8 between P 0 and P +,s and its symmetrical counterpart -(2.2) -connecting P +,s back to P 0 .The other boundary corresponds to Φ ↑ Φ s−t : the distance traveled along M + ε decreases to 0 and γ s−2f (ξ) detaches from M + ε .However, since the angle between u and s is determined by √ Φ (2.14), this can only happen as also Ψ/Φ ↓ 9/(2 + 9a), i.e. as the 'projected triangle' of Fig. 11a that represents the skeleton Γ s−2f entirely contracts to a point -see also Remark 3.14.The bifurcational structure associated to the situation sketched in Fig. 11b is quite different: as Φ decreases towards Φ 1 s−1f < Φ s−t , γ s−2f (ξ) does not merge with a (superposition of 2) 1-front(s) γ s−1f (ξ) of Theorem 3.8.In fact, γ s−2f (ξ) does not bifurcate at all, the distance between the 2 fronts of γ s−2f (ξ) remains bounded as Φ passes through Φ 1 s−1f (the main difference between the cases Φ > Φ 1 s−1f and Φ < Φ 1 s−1f is the sign of H s−d − H +,s 0 : γ s−2f (ξ) follows an orbit of (2.15) outside the homoclinic loop for Φ < Φ 1 s−1f ).Moreover, as Φ increases towards Φ s−t , the distance γ s−1f (ξ) travels along M + ε also does not go to 0, in fact γ s−1f (ξ) (almost) follows the entire periodic orbit of the level set {H + 0 (w, q) = H s−d (Φ s−t ) = H s−t }, the critical/limiting orbit of Theorem 3.9.Again, this can only happen if also the projection of P 0 on M + ε merges with this periodic orbit.We refrain from going any further into the details of these -and other -bifurcations of γ s−2f (ξ).
Proof of Theorem 3.11.The proof follows the geometrical approach developed in [19,18], equivalently the more analytical approach of [40] could be employed.The construction of γ s−2f (ξ) is based on the 'intermediate' orbit , the heteroclinic connection between P 0 and M + ε that touches down on J s−d ∈ M + ε ; γ i−1f (ξ) follows the slow flow along M + ε and is thus asymptotically close to the skeleton Γ s−2f up to the take off point ε ) and thus cannot take off from M + ε ).The homoclinic orbit γ s−2f (ξ) is constructed as a symmetric orbit -i.e. an orbit that passes through the plane {p = q = 0} at its 'midpoint' -that is exponentially close to γ i−1f (ξ) up to the point it takes off from -by definition.We define for some sufficiently small σ (independent of ε), the (bounded) 1-dimensional curve C σ i−1f ⊂ {b = 1/2} as the (first, transversal) intersection of W u (P 0 ) and {b = 1/2} that is at a distance of maximal σ away from P i−1f ; in other words, provides a parametrization of orbits γ(ξ) in W u (P 0 ) near γ i−1f (ξ) .In fact, the saddle structure of the fast flow around M + ε cuts W u (P 0 ), and thus i−1f cross through the plane {p = 0} near M + ε so that their b-coordinate changes direction; the b-coordinates of γ(ξ)'s with γ(0) in C σ i−1f \ C σ,r i−1f do not change direction, these γ(ξ)'s pass along M + ε without the possibility of returning to M 0 ε .Thus, we may uniquely parameterize orbits γ(ξ) ⊂ W u (P 0 ) that pass through {p = 0} near by M + ε that are σ-close to γ i−1f (ξ) by the distance d between their initial point γ(0) ∈ C σ,r i−1f and P i−1f (∈ ∂C σ,r i−1f ) -where we have implicitly used the fact that [39,41].We denote these γ(ξ)'s by γ(ξ; d).
The construction of stationary homoclinic 2-front gap patterns -localized bare soil areas surrounded by vegetation -goes along exactly the same lines as the above construction of localized spot patterns.The main difference is that the homoclinic orbits-to-be-constructed are ⊂ W u (P +,s ) ∩ W s (P +,s ) so that the structure of orbits taking off and touching down now has to start out from the saddle P +,s ∈ M + ε .Nevertheless, the construction of the skeleton structure Γ g−2f is completely similar to that of Γ s−2f .Therefore, we only provide the essence of the construction of Γ g−2f .
First, we need to assume that there is a critical point P +,s ∈ M + 0 of saddle type.The skeleton structure Γ g−2f consists of a piece of W u (P +,s ) ⊂ M + 0 from P +,s up to the (ε → 0 limit of the) take off point J +,0 g−o from M + 0 , followed by (a piece of) W u (M + 0 ) ∩ W s (M 0 0 ) (2.10) up to the (ε → 0 limit of the) touch down point J +,0 g−d ∈ M 0 0 (that has the same (w, q)-coordinates as J +,0 g−o in the limit ε → 0).Note that the take off/touch down points J +,0 g−0 /J +,0 g−d differ essentially from their counterparts as J s−d /J s−0 (3.6) considered so far: while Nevertheless, the coordinates of all take off/touch down points are at leading order determined by their ε → 0 limits (2.10) with w ± h (0) = 9/(2 + 9a) (2.9).The next piece of Γ g−2f consists of a symmetric part of a (cosh-type) orbit along M 0 0 of the (linear) slow reduced flow (2.13) up to the (ε → 0 limit of the) take off point J 0,+ g−o , which is again followed by a fast jump backwards along (a piece of) W s (M + 0 ) ∩ W u (M 0 0 ) to the (ε → 0 limit of the) touch down point J 0,+ g−d ∈ M + 0 .The final piece is the symmetrical counterpart of the first piece: the flow of (2.15) along M + 0 from the final touch down point back towards P +,s -see Fig. 12a for a sketch of a projection of Γ s−2f (without its fast jumps).
The (open) region S g−2f is defined by those (a, Ψ, Φ, Ω, Θ)-combinations for which Γ g−2f can be constructed.We note that there cannot be points in the intersection of S s−2f as defined in Theorem 3.11 and S g−2f : the (projections e s 9 y 8 A 2 p J I < / l a t e x i t > P +,s < l a t e x i t s h a 1 _ b a s e 6 4 = " B 9 f z 2 0 T U R y b K y w g T t v 9 I 5 w m e H s + + V 5 0 j y p O q d V 9 8 a t 1 N w 8 j i L s w T 4 c g g P n U I N r q E M D C H B 4 h G d 4 s e 6 t J + v V e p u 1 F q x 8 Z h d + w X r / A h I 0 j 2 4 = < / l a t e x i t > P 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " F Z 0 8 u A J 0 j M N B X w U o P i I W Z W j r i 0 e s 9 y 8 A 2 p J I < / l a t e x i t > P +,s < l a t e x i t s h a 1 _ b a s e 6 4 = " B 9 f z 2 0 T U R y b K y w g T t v 9 I 5 w m e H s + + V 5 0 j y p O q d V 9 8 a t 1 N w 8 j i L s w T 4 c g g P n U I N r q E M D C H B 4 h G d 4 s e 6 t J + v V e p u 1 F q x 8 Z h d + w X r / A h I 0 j 2 4 = < / l a t e x i t > P 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " F Z 0 8 u A J 0 j M N B X w U o P i I W Z W j r i 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " M 3 q v 0 W G f w q 9 z E Q n D t a T E 8 / L S K j o = " > A A A B 6 n i c b V D L S s N A F L 2 p r 1 p f V Z d u B o v g q i R a q u 4 K b l z W R x / Q x j K Z T t q h k 0 m Y m Q g l 9 B P c u F D E r V / k z r 9 x k g Z R 6 4 E L h 3 P u 5 d 5 7 v I g z p W 3 7 0 y o s L a + s r h X X S x u b W 9 s 7 5 d 2 9 t g p j S W i L h D y U X Q 8 r y p m g L c 0 0 p 9 1 I U h x 4 n H a 8 y W X q d x 6 o V C w U d 3 o a U T f A I 8 F 8 R r A 2 0 i 2 / j w f l i l 2 1 M 6 B F 4 u S k A j m a g / J H f x i S O K B C E 4 6 V 6 j l 2 p N 0 E S 8 0 I p 7 N S P 1 Y 0 w m S C R 7 R n q M A B V W 6 S n T p D R 0 Y Z I j + U p o R G m f p z I s G B U t P A M 5 0 B 1 m P 1 1 0 v F / 7 x e r P 1 z N 2 E i i j U V Z L 7 I j z n S I U r / R k M m K d F 8 a g g m k p l b E R l j i Y k 2 6 Z S y E C 5 S 1 L 9 f X i T t k 6 p z W q 1 d 1 y q N m z y O I h z A I R y D A 2 f Q g C t o Q g s I j O A R n u H F 4 t a T 9 W q 9 z V s L V j 6 z D 7 9 g v X 8 B f N 2 O H w = = < / l a t e x i t > + p < l a t e x i t s h a 1 _ b a s e 6 4 = " m v D U w R j q < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 M B J n n T g N m y 5 G v l g X M n 2 c r p p K 9 0 = " > A A A B 9 X i c b V D J S g N B E K 2 J W 4 x b 1 K O X x i B E l D C j w e U W 8 O I x i l k g m Y S e T k / S p G e h u 0 c J w / y H F w + K e P V f v P k 3 9 k y C q P F B w e O 9 K q r q O S F n U p n m p 5  of the) take off/touch down points J s−d /J s−0 lie on the level set {H + 0 (w, q) = H s−d } for which H s−d < H +,s 0 , the value of H + 0 (w, q) for P +,s and its (un)stable manifolds (2.22).By construction, J +,0 g−0 , J +,0 g−d ⊂ {H + 0 (w, q) = H +,s 0 }compare Fig. 11a to Fig. 12a.This also implies that S s−2f ∩ S g−2f = ∅, in fact, ∂S s−2f ∩ ∂S g−2f ⊃ R s−1f as defined in Theorem 3.8: both the homoclinic spots γ s−2f (ξ) of Theorem 3.11 and the gaps γ g−2f (ξ) of (upcoming) Theorem 3.13 merge with the heteroclinic cycle spanned by the standing 1-front γ s−1f (ξ) of Theorem 3.8 and its symmetrical counterpart as S s−2f ∪ S g−2f approaches R s−1f .Theorem 3.13.Let (a, Ψ, Φ, Ω, Θ) ∈ S g−2f and Γ g−2f ⊂ R 4 be the singular (gap) skeleton constructed above.Then, there is for ε > 0 sufficiently small a symmetric 2-front orbit γ g−2f (ξ) = (b g−2f (ξ), p g−2f (ξ), w g−2f (ξ), q g−2f (ξ)) ⊂ W u (P +,s ) ∩ W s (P +,s ) of (1.7) with c = 0 that merges with Γ g−2f as ε ↓ 0. The associated stationary pattern (B(x, t), W (x, t)) = (b g−2f (x), w g−2f (x)) in (1.5) represents a stationary localized bare soil gap embedded in vegetation.
Of course the proof of this Theorem goes exactly along the lines of the proof of Theorem 3.11.The main difference between the cases of (stationary, symmetric, homoclinic) 2-front spots and (stationary, symmetric, homoclinic) 2front gaps is that there cannot be periodic orbits on M 0 ε -the slow reduced flow (2.13) on M 0 0 is linear -so that there cannot be any higher order localized gap patterns (as in Theorem 3.12 for localized spots).We refer to Figs. 1c and 15a for numerical observations of -(most likely) stable -localized 2-front spot and gap patterns in PDE (1.5).
Remark 3.14.We refer to [43,44] for studies of the process of a homoclinic 2-front orbit between a slow manifold M 1 ε and a second slow manifold M 2 ε detaching from M 2 ε to become a slow-fast homoclinic to M 1 ε that only makes 1 homoclinic excursion through the fast field (instead of 2 fast heteroclinic jumps between M 1 ε and M 2 ε ).The focus of [43] is on the (exchange of ) stability between these 2 types of homoclinic patterns and the associated bifurcations -especially as localized stripes in 2 space dimensions.In the present work the situation is somewhat more involved than in [43,44], since the skeleton structures as sketched in Figs.11 and 12a must become asymptotically small in this transition -which is not necessary in the setting of [43,44].Notice that this implies that here the W -component of the homoclinic pulse becomes 'small' -a certain well-defined magnitude in ε -during this transition, but that this is not the case for the B-component since the orbit still has to make (almost) a full jump between M 0 ε and M + ε .See Figs.14c, 15b and 16 in section 4.1.
Proof of Theorem 3.15.This proof can be set up very much along the lines of the proofs of similar results -the existence of spatially periodic patterns in the (generalized) Gierer-Meinhardt equation -in [21], therefore we restrict ourselves to the essential ingredients of the proof here.
The approach is similar to that of the proof of Theorem 3.11: we construct an orbit that intersects the plane {p = q = 0} and apply the reversibility symmetry (2.2) -with c = 0.However, unlike for the homoclinic orbits in Theorem 3.11, we do not 'start out' -as ξ → −∞ -at the critical point P 0 ∈ M 0 ε , but choose the initial condition of the orbit-to-be-constructed at the -exponentially short -interval As in the case of the homoclinic spots patterns -Theorem 3.12 -we may immediately conclude that there are countably many (families of) higher order periodic patterns if Γ + ρ (as defined in Theorem 3.15) is part of a periodic orbit on M + 0 (determined by the level set {H + 0 (w, q) = H ρ }) -see Fig. be as defined in Theorem 3.15, with Γ + ρ such that it is part of a closed orbit on M + ε determined by the level set {H + 0 (w, q) = H ρ }.Then, for ε > 0 sufficiently small, there is a countable family of symmetric multi-front periodic orbits γ i mf,ρ (ξ) = (b i mf,ρ (ξ), p i mf,ρ (ξ), w i mf,ρ (ξ), q i mf,ρ (ξ)) of (1.7) with c = 0 (i = 1, 2, ...) that merges in the limit ε → 0 with the extended skeleton structure spanned by Γ 0 ρ , the fast jump over W u (M 0 0 )∩W s (M + 0 ) with w + h = 9/(2+9a), q = −ρ (2.10), the full closed orbit of {H + 0 (w, q) = H ρ } that contains Γ + ρ and the fast jump back over W u (M + 0 )∩W s (M 0 0 ) with w − h = 9/(2 + 9a), q = ρ.The associated stationary patterns (B(x, t), W (x, t)) = (b i mf,ρ (x), w i s−2f (x)) in (1.5) are symmetric periodic spot/gap patterns with an increasing number of oscillations in the vegetated areas.
Finally, we note that the families of 'higher order' periodic patterns γ i mf,ρ (ξ) are only the first of further -more complex -families containing periodic (and aperiodic) patterns of increasing complexity.We refer to [21] for the precise settings and proofs, here we only give a sketch of one specific example.However, this sketch provides the main ideas by which all further orbits may be constructed.
ξ) does not close as it passes along I b , instead it keeps on following γ i mf,ρ (ξ) as it makes it second jump towards M + ε .By the approach of the proofs of Theorems 3.11 and 3.15, we can now tune b so that it has its second take off from M + ε precisely and that it passes through {p = q = 0} while following Γ + ρ (without making any further circuits over the periodic orbit on M + ε that contains Γ + ρ ).It follows by the application of the reversibility symmetry (2.2) that for this value of b, γ ρ,b (ξ) is a symmetric periodic orbit that 'starts' at Γ 0 ρ ⊂ M 0 ε , jumps to M + ε to make i circuits along M + ε , jumps back again to M 0 ε , follows Γ 0 ρ ⊂ M 0 ε to return again to M + ε where it follows Γ + ρ and subsequently immediately jumps back again to M 0 ε -from which it repeats the same path, etc..Note that the associated periodic pattern in (1.5) consists of an alternating array of 2 different types of localized vegetation spots.Clearly, this procedure can be further refined to establish the existence of patterns containing arbitrary arrays of arbitrarily many different types of vegetation spots -under the conditions of Corollary 3.16.
Remark 3.17.We decided to focus in this paper on stationary 2-and multi-front patterns.Of course, (1.5) also exhibits traveling multi-front patterns -see for instance Fig. 14d in which a vegetation spot travels towards a stationary, stable (and attracting) spot of the type established by Theorem 3.11.System (1.7) can also have homoclinic orbits to P 0 for (certain specific values of ) c = 0, i.e. vegetation spots may be traveling with constant speed (without changing shape).An approach along the lines of [22] indicates that bifurcations to traveling spots appear when the touch down manifold I down (Lemma 3.2) is tangent to a level set {H + 0 (w, q) = H} of the slow reduced flow on M + 0 with c = 0 at the (non-transversal) intersection I down ∩ {H + 0 (w, q) = H} (recall that I down is parameterized by c).A similar property holds for bifurcations of stationary spatially periodic multi-front patterns into traveling spatially periodic (wave train) patterns.A simple investigation of the relative orientations of I down and the various possible phase plane configurations of the slow reduced flow on M + 0 shows that there indeed are parameter combinations (a, Ψ, Φ, Ω, Θ) at which these bifurcations into traveling 2-/multi-fronts must occur.This bifurcation may have relevant ecological implications, nevertheless, we refrain from going into the details here (and leave this to future work) -see also section 4.2.

Simulations
The motivation for the numerical simulations presented in this section is threefold: 1) to illustrate some of the analytic results of the previous sections (without doing a systematic search for all constructed patterns) 2) to give a brief outlook beyond the worked out analysis to solution types that may be constructed by the geometric set-up developed here and, finally, 3) to give a flavor of the rich dynamics that PDE (1.1)/(1.5)exhibits.All numerical simulations have been carried using MATLAB's 'pdepe' routine.The corresponding parameter settings are specified in the captions of the figures.Almost all figures show a snapshot in time of the spatial profile of the PDE solution/pattern after it converged to a stationary or uniformly traveling solution.
The opening figure of this section, Fig. 13a, can be seen as a binding element between the papers that have a more ecological emphasis and motivated the present work -see [76,77] and the references therein -and the analysis here.It displays a traveling front solution as established by Theorem 3.4 for a parameter regime comparable to the one from [77] (with slight adjustment in the parameters to compensate for the choice of a 1-D model -A).This profile is then shown in Fig. 13b as a projection in (b, w, q)-space to illustrate that it indeed starts on the slow manifold M 0 ε and then jumps to the M + ε slow manifold.As established by the analysis, the solution first follows the unstable manifold associated to the bare soil state ( B0 , W 0 ) (as a solution of (1.7)), makes a fast excursion through the fast field to then touch down on M + ε following the stable manifold associated to the uniform vegetation state ( B+ , W + ).Note, of course, that this figure contains two approximations: first, the manifold M + ε is only accurate up to second order in ε -see section 2.4 -while the flows on M 0 ε and M + ε are computed numerically (using MATLAB routines).As demonstrated in section 3, heteroclinic 1-front orbits can occur both as traveling -Theorem 3.4 -or stationary patterns -Theorem 3.8.We confirm this numerically in Figs.14a and 14b.Note that these fronts may either represent the retreat of vegetation by the invasion of the bare soil state into the homogeneous vegetation statec > 0 -or the expansion of a (homogeneously) vegetated area into the bare soil state -c < 0. In fact, in order to find the stationary 1-front, we need to tune a single parameter -Φ in the statement of Theorem 3.8 and Ψ in Fig. 14 -to the border point between the ranges of left-traveling and right-traveling 1-fronts (the 'Maxwell point' [6] described by the co-dimension 1 set R s−1f in Theorem 3.8).
The existence of the homoclinic stationary 2-front pattern depicted in Fig. 14c was established in Theorem 3.11.Note that the level of vegetation on the plateau that determines the spot remains relatively far away from the value B+ of the uniform vegetation state ( B+ , W + ).This is caused by the fact that the homoclinic orbit associated to the spot pattern follows an orbit on the slow manifold M + ε of (1.7) that does not approach the critical point associated to ( B+ , W + ) on M + ε -see the sketch of the skeleton structure in Fig. 11.From the ecological point of view, a vegetation spot benefits from soil water diffusion from the adjacent water-rich bare soil areas -see the W -profile in Fig. 13a -besides direct rainfall, and therefore has higher biomass density as compared with uniform vegetation.This also explains (in ecological terms) why the biomass density at the edge of a front, a spot or a gap is higher -a property also exhibited by 'fairy circles' [73,77].See also the upcoming discussion below of 2-front vegetation gaps (Fig. 15a).In Fig. 14d, we show the dynamics of the 2 interacting fronts of an evolving 2-front pattern: the distance between the fronts slowly increases while it settles into a stationary standing spot -see Remark 3.17 and the discussion in section 4.2.
One of the original motivations to analyze far-from-equilibrium patterns in the Gilad et al. model, was to gain a fundamental understanding of 'fairy circles' -a somewhat subtle phenomenon (for instance) observed in western Namibia [73,77].The homoclinic stationary 2-front gap patterns of (1.5) established by Theorem 3.13 and shown in Fig. 15a indeed have the strongly localized nature of observed fairy circles.Moreover, the spot/gap patterns of Theorem 3.15 represent the observed (nearly) periodic arrays of fairy circles (see Fig. 17a and notice that the ratio between the lengths of the vegetated state and the bare soil patches typically varies from 0 to ∞ in this family (section 3.6)).As noted, fairy circle gap patterns have an excess of vegetation at the edge of the gap as distinctive feature -see for instance the images in [73,77].In mathematical terms, this means that the connecting fronts are non-monotonous.In the context of the present model, this non-monotonicity is caused by the orientation and curvature of the slow manifold M + ε relative to the path traced by the gap pattern over M + ε -see the projection in Fig. 16a for a representation of this 'geometrical mechanism' for spot patterns.We refer to [25] for a 1-component model in which the non-monotonicity of the fronts originates from nonlocal effects.In section 3.5 -and especially in Remark 3.14 -we discussed the bifurcation of the homoclinic slow-fast-slow-fastslow 2-front spot pattern of Theorem 3.11 into a homoclinic slow-fast-slow pulse pattern as it 'detaches' from M + ε .Such a (numerically stable) 'detached' spot pattern of pulse type is shown Fig. 15b.In Fig. 16, the detachment process is shown by projections of both the 2-front spot pattern of Fig. 14c and the pulse pattern of Fig. 15b on the (w, b)-plane: as the parameter Ψ -which is linearly related to the rainfall P in the original model (1.6) -is decreased below a critical value Ψ * = 1.2952 the vegetated plateau disappears and the 2-front spot solution transforms into a 1-pulse solution.Note that this 1-pulse solution is of the 'classical' Klausmeier-Gray-Scott (and/or Gierer-Meinhardt [18]) type already discussed in the introduction of section 3.5: its existence may be established by the methods of [18,23] and the references therein.We also found that the spatially periodic spot/gap patterns of Theorem 3.15 may have quite a large domain of attraction: Fig. 17b shows the evolution of a traveling vegetation front into the bare soil state that leaves behind a spatially periodic spot/gap pattern -Fig.17a.This behavior may possibly be related to the existence of a Turing bifurcation -see Remark 2.2 -of the uniform vegetation state and calls for further  studies.Finally, we show in Fig. 18a such a numerically obtained, almost sinusoidal, small amplitude Turing pattern that bifurcated from a (destabilized) uniform vegetation state and note that there are paths through parameter space on which this Turing pattern evolves into a (periodic) multi-pulse pattern -built from homoclinic 1-pulses of Gierer-Meinhardt/Gray-Scott type as in Fig. 15b and typically observed in Klausmeier-type models [60,65,69] that subsequently touches down on M + ε like the solitary pulses of Figs.15b and 14c, to indeed take the shape of the periodic fairy circle-type spot/gap pattern of Theorem 3.15 and Fig. 17a.By further tuning parameters it may also happen that the stationary, spatially periodic, Turing pattern undergoes a Hopf bifurcation (in time), resulting in an oscillating pattern that is periodic both in space and in time -see Fig. 18b.(Note, however, that it is not clear whether this may occur for ecologically feasible parameters -see [66].)

Discussion
Of course, the potential relevance of the various singular slow/fast patterns constructed in this paper is ultimately determined by their stability as solutions of PDE (1.1)/(1.5).In general, this is a seriously hard problem to study analytically.However, the singularly perturbed nature of the patterns considered here enables us to explicitly and rigorously analyze the (spectral) stability of the constructed (multi-)front patterns.In fact, the explicit 'control' we established on the slow-fast structure of the (multi-)fronts provides the perfect (and necessary) starting point for a spectral stability analysis along the lines of (for instance) [20,7] and [14] (for the spatially periodic patterns).This is especially the case for the basic front/spot/gap/periodic patterns of Theorems 3.4, 3.8, 3.11, 3.13, 3.15 shown in Fig. 1.
However, the question whether the non-basic, 'higher order' patterns (for instance) sketched in Fig. 2 can be stable also requires novel mathematical insights and methods.All constructed higher order patterns involve the existence of persisting periodic solutions on the slow manifold M + ε -see Theorems 3.5, 3.6, 3.9, 3.12 and Corollary 3.16.Therefore, the structure of the spectrum associated to the stability of the higher order patterns essentially depends on the preliminary question about the spectrum and stability of the persistent periodic solutions on the slow manifold of Theorem 2.4 -and especially of their homoclinic (or heteroclinic) limits also considered in Theorem 2.4.In fact, this issue is not (at all) specific for the explicit model here.We claim that higher order patterns of the type sketched in Fig. 2 will generically appear as potentially stable solutions in a fully general class of singularly perturbed reaction-diffusion models that includes (1.5), By going into a traveling framework -and thus introducing ξ = x − ct, U (x, t) = u(ξ), V (x, t) = v(ξ), p(ξ) = u ξ (ξ), q = v ξ (ξ)/ε as in section 1 -(4.1) is reduced into the 4-dimensional form of (1.7).By taking the ε → 0 limit, we find that the 2-dimensional (reduced) slow manifolds are determined by F (v 0 , u) = 0 (and p = 0, (v 0 , q 0 ) ∈ R 2 ) -see section 2 -which generically determines J ≥ 1 branches, locally given by graphs, M j 0 = {(u, p, v, q) ∈ R 4 : u = f j (v), p = 0}, j = 1, 2, ..., J, with f j (v) such that F (f j (v), v) ≡ 0 (cf.(2.5) and note that J = 3 for (1.5)).For those (parts of) M j 0 that are normally hyperbolic, M j 0 persists as M j ε , that is approximately given by, ∂F ∂u (f j (v), v), p j 1 (v) = f j (v).
(cf. (2.28), (2.29)).Thus, completely analogous to the analysis in section 2.4, we find that the slow flow on a persisting, normally hyperbolic 2-dimensional slow manifold M j ε is given by a planar Hamiltonian system perturbed by a nonlinear friction term, with X = εξ (cf.(2.30)).Typically, the unperturbed ε → 0 limit v XX + G(f j (v), v) -i.e. the reduced slow flow on M j 0 -is nonlinear and has families of periodic solutions and homoclinic or heteroclinc orbits to critical points on M j ε that correspond to (potentially stable [17]) homogeneous background states (U (x, t), V (x, t) ≡ ( Ū , V ) of PDE (4.1) -as is the case for (2.15) on M + 0 .Thus, indeed, the situation is completely similar to that of section 2.4: using Melnikov-type arguments persistence results equivalent to Theorem 2.4 may be deduced, also in the present general setting.The geometric framework of orbits 'jumping up and down' between two (normally hyperbolic) slow manifolds M j ε and M k ε presented in section 3.1 is based on the persistence of both the stable and unstable manifolds W s,u (M j,k ε ) of M j,k ε and thus of the intersections W u (M j ε ) ∩ W s (M k ε ) and W s (M j ε ) ∩ W k (M k ε ).Therefore, we may use the arguments, methods and insights of section 3 to deduce the equivalents of the 'higher order' existence Theorems 3.5, 3.6, 3.9, 3.12 and Corollary 3.16 in the setting of general system (4.1).Moreover, this also implies that bifurcation scenarios as sketched in Fig. 9 appear generically (where we notice that the sketch in Fig. 9 was just a first example -many other scenarios may occur).In fact, the geometrical setting allows us to (for instance) explicitly establish the existence of heteroclinic networks of orbits jumping between various slow manifolds M j ε and (slowly) following periodic orbits on M j ε in between its fast jumps -like the networks considered in [53,54] and the references therein.Thus, the above noted preliminary (and essential) issue of the spectrum associated to the stability of the persisting periodic and homoclinic solutions on M + ε of Theorem 2.4 also has a fully general -and thus fundamental -counterpart, with a similar relevance for the higher order patterns (almost) heteroclinic to these orbits.In other words, insight in the spectrum associated to the stability of the orbits on M j ε established by a generalization of Theorem 2.4 for (4.1) is expected to yield explicit insight in the stability and bifurcations of the higher order patterns in PDE (4.1) established by the generalizations of Theorems 3.5, 3.6, 3.9, 3.12, etc. and the subsequent more complex 'networks'.
This will be the subject of future work, both in the setting of explicit system (1.5) -which will also include a systematic numerical search for the higher order patterns sketched in Fig. 2 -and in the general setting of (4.1).
To optimally embed the present analysis in the ecological context, we first need to obtain insight in the ranges of the (scaled) parameters (a, Ψ, Φ, Ω, Θ) of (1.5) that may correspond to ecologically realistic settings of the unscaled model .Finding such parameter ranges is possible (more than in the Klausmeier model) since (1.1) is directly linked to the more elaborate 3-component model of Gilad et al. [28] - [77], A -and thus to concrete underlying ecological mechanisms [29,47,61].A crucial question for the potential ecological relevance of the above discussed higher order patterns is whether there are realistic values of (Λ, K, E, M, P, N, R, Γ) for which there are 2 critical points on M + ε , i.e. for which C 2 − 4AD > 0, C < 0, D > 0 (section 2.3) -where (A, C, D) is related to (Λ, K, E, M, P, N, R, Γ) in a rather nonlinear fashion by (1.3), (1.5), (1.6) (2.16), (2.17).Naturally, this will be part of our upcoming work on (1.5).
Each of the higher order invasion fronts established by Theorem 3.6 and sketched in Figs.2a and 2b travels with a different speed -in fact, the (discrete) family may even limit on a stationary front pattern (Remark 3.10).Thus, when stable, these invasion fronts may introduce the possibility of slowing down gradual desertification.Moreover, stationary multi-front patterns may bifurcate into traveling patterns with the same structure -see Remark 3.17 for a brief sketch of the underlying geometrical mechanism.When stable, the appearance of such traveling multi-front patterns -either localized spots or spatially periodic wave trains -may have a similar ecological interpretation and relevance: localized vegetated states may even reverse desertification by invading bare soil -see [74,76,78].Together, the various traveling 1-front patterns and traveling multi-spots form an interacting group of invasion patterns within the transition zone between the bare soil state and a homogeneous vegetation state; in principle all entities in this group travel with different speed.Understanding pattern formation in this zone -and especially also understanding the translation and/or expansion of this zone in terms of the parameters in the model -may have a direct ecological significance.In mathematical terms, such a study may also be performed by a front interaction analysis along the lines of [67,11,12] -although the dynamics generated by (1.5) in this ecological transition zone is expected to be richer than that of the generalized FitzHugh-Nagumo model considered there.
To truly obtain ecological relevance, we must consider the model in 2 space dimensions.Clearly, the extension to more than 1 space dimension does pose fundamental challenges, moreover 2-dimensional systems show much richer dynamical behaviors associated with propagating fronts -see for instance [35,36].However, the results obtained here form a foundation upon which aspects of the step from 1 to 2 space dimensions can be taken -see for instance [68,64,71,72] and the references therein.By extending the patterns constructed here trivially in the second spatial direction, the above mentioned stability analysis can be directly extended to include the stability (and bifurcations) of planar (multi-)fronts/interfaces (where we for simplicity neglect the (technical) fact that (1.5) takes a somewhat different form in R 2 , [47,77], A).Unlike in the extended Klausmeier model [60,64], localized stripes are of 2-front type and may thus be expected to possibly be stable -see [1] for a rigorous treatment in a generalized Klausmeier type model (posed on a sloped terrain without a diffusion term for the water component -like the original Klausmeier model [42]).Naturally, the interfaces will evolve and their curvature driven dynamics may be studied analytically along the lines of [51].Especially in the above discussed multi-front transition region between bare soil and homogeneous vegetation, the ecosystem dynamics generated by the model may be very rich and complex -see for instance [35,36,37].As a final direction of possible future research, we note that our results may be used to establish the existence -and later stability -of localized patterns in the original 3-component model of Gilad et al. [28].Since (1.1) and thus (1.5) -is obtained from the nonlocal, 3-component model of Gilad et al. (see (A.1)) in a systematic wayi.e. by taking several limits ( [47,77], A) -it may be expected that it is possible to establish the persistence of patterns constructed here into the nonlocal, 3-component setting, especially since these patterns are constructed geometrically through transversal intersections of invariant manifolds.Once again, this is interesting and relevant both from mathematical and ecological point of view: (asymptotically) small nonlocal and topographical terms may have a significant effects, even on the most simple -'basic' -(vegetation) patterns exhibited by a model [2,24].
A further simplification we make is related to the nonlocal forms of the biomass growth rate, G B , and the water uptake rate, G W , in (A.1).We assume, consistently with the plant species in the Namibian fairy-circles ecosystems, that the roots, described by the root kernel G(X, X , T ), are laterally confined.We employ this assumption by taking the lateral root extension of a seedling, S 0 , to be very small.Using the limit S 0 → 0 in the integrals in (A.2c) we obtain, for a 1-dimensional system, the simpler algebraic expressions Inserting these expressions in (A.1) we obtain the 2-component model (1.1).Finally, we note that in [77] this reduction was performed in 2 space dimensions and that the general n-dimensional situation is considered in [47].
B The derivation of the scaled model which can be brought into the form, Note that our choice to scale the factor of the term W B 2 in the B-equation to +1 implies -together with the (implicit, natural) assumption that B and B have the same signs (1.2) -that we have chosen to consider EK > 1.
Of course, it may happen that 0 < EK ≤ 1.In these cases, either the term W B 2 disappears from the equation -in the 'non-generic' case EK = 1 -or its pre-factor can be scaled to −1.All of the analysis in this work can also be performed for EK ≤ 1, without any conceptual differences.However, we chose to focus of EK > 1 -and thus on a +W B 2 term in (1.5) -to not further complicate the necessary 'algebra'.

C Lemma 2.6 and the Bogdanov-Takens bifurcation scenario
A planar ODE of the form y X = z , z X = β 1 + β 2 y + y 2 + syz + G(y, z)., for β 1 , β 2 ∈ R, s = +1 is known to possess two fixed points -a saddle and a focus -with an unstable periodic orbit (that emerged from the focus in a Hopf bifurcation in the open parameter region) The right border {β 2 < 0, β 1 = 0} marks the Hopf bifurcation, while the left border {β 2 < 0 , β 1 = − 6 25 β 2 2 + o(β 2 2 )} describes the region where a homoclinic orbit emerged from the periodic orbit (whose period tends to infinity towards that border).

Figure 5 :
Figure 5: Various relative configurations of the nullclines (2.27) and the associated critical points for w ∈ U a (2.6).

Figure 6 :
Figure 6: Sketches of the intersections of I down and W s (P +,s )| M + 0 in M + 0 , i.e. the leading order configurations as described by the (integrable) slow reduced flow (2.15), in the 2 cases considered in Theorem 3.4: there is one critical point P +,s of saddle type on M + 0 or there is a center P +,c and a saddle P +,s on M + 0 .
Fig 6a -it follows that there must be a transversal intersection I down ∩ W s (P +,s )| M + 0 for values of Φ in an (open) interval around Φ +,s .Transversality implies that the intersection persists under varying A, C, D around their initially frozen values, which establishes the existence of the open region S 1 s−prim in (a, Ψ, Φ, Ω, Θ)-space for which I down and W s (P +,s )| M + 0 intersect.Moreover, the manifold W s 6(b).By varying Φ around Φ = Φ +,c and A, C, D around their initially frozen values, we find the open region S 2 cs−prim in (a, Ψ, Φ, Ω, Θ)-space for which both elements of the intersection I down ∩ W s (P +,s )| M + 0 persist: for (a, Ψ, Φ, Ω, Θ) ∈ S 2 cs−prim , (1.7) has 2 (distinct) primary heteroclinic orbits γ j prim (ξ), j = 1, 2, that correspond to 1-front patterns traveling with speeds c 1 prim = c 2 prim -where c j,0 prim is the unique solution of w + h (c) = wj,0 prim .Finally, we note that the existence of the open set S 1 cs−prim follows by considering I down ∩ W s (P +,s )| M + 0

1 a
< l a t e x i t s h a 1 _ b a s e 6 4 = " t y 1 o u

Figure 8 :
Figure 8: Sketches of 3 relative configurations of I down (Φ) and W s (P +,s )| M + 0 for increasing Φ that represent 3 distinct stages in the bifurcation scenario of Fig. 9.
e per < l a t e x i t s h a 1 _ b a s e 6 4 = "g p R O 0 k 7 Y E 8 e A A U e j k S I M E 0 q 2 u M Q = " > A A A B 9 X i c b V D L S s N A F J 3 U V 4 2 v q k s 3 g 0 V w V R I r P n Y F N y 6 r 2 A e 0 a Z l M b 9 q h k 0 m Y m S g l 9 D / c u F D E r f / i z r 9 x 0 g Z R 6 4 E L h 3 P u 5 d 5 7 / J g z p R 3 n 0 y o s L a + s r h X X 7 Y 3 N r e 2 d 0 u 5 e U 0 W J p N C g E Y 9 k 2 y c K O B P Q 0 E x z a M c S S O h z a P n j q 8 x v 3 Y N U L B J 3 e h K D F 5 K h Y A G j R B up 1 6 2 P W D + N Q U 5 7 Y N v 9 U t m p O D P g R e L m p I x y 1 P u l j + 4 g o k w A C N J I y < / l a t e x i t > prim < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 wt A x j b b U 8 r X O G F 0 5 N x l 5 S c 0 z 9 o = " > A A A B 9 H i c b V D L S g N B E O y N r 7 i + o h 6 9 D A b B U 9 h V 8 X E L e P E Y x T w g W c L s Z J I M m Z l d Z 2 Y D Y c l 3 e P G g i F c / x p t / 4 2 y y i B o L G o q q b r q 7 w p g z b T z v 0 y k s L a + s r h X X 3 Y 3 N r e 2 d 0 u 5 e Q 0 e J I r R O I h 6 p V o g 1 5 U z S u m G G 0 1 a s K B Y h p 8 1 w d J 3 5 z T F V m k X y 3 k x i G g g 8 k K z P C D Z W C j q 1 Ie u m s W J i 6 r r d U t m r e D O g R e L n p A w 5 a t 3 r e e r F f r b d 5 a s P K Z f f g F 6 / 0 L 4 / 2 N L g = = < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " a W k t 2 s G O p C F h 9 q a w t v j u o B a e h J I = " > AA A B + n i c b V D L S s N A F L 2 p r 1 p f q S 7 d D B a h b k p q x c e u 4 M Z l F f u A N p T J d N I O n T y Y m S g l 5 l P c u F D E r V / i z r 9 x k g Z R 6 4 E L h 3 P u n T v 3 O C F n U l n W p 1 F Y W l 5 Z X S u u l z Y 2 t 7 Z 3 z P J u R w a R I L R N Ah 6 I n o M l 5 c y n b c U U p 7 1 Q U O w 5 n H a d 6 W X q d + + o k C z w b 9 U s p L a H x z 5 z G c F K S 0 O z P M j e i B 0 e 0 S S u k q

r 6 2 K
P S j r O 1 C T r U y g i 5 g d D l K 5 S p P y d i 7 E k 5 8 x z d 6 W E 1 k X + 9 V P z P 6 0 f K P b d j 5 o e R o j 6 Z L 3 I j j l S A 0 h z Q i A l K F J 9 p g o l g + q + I T L D A R O m 0 S l k I F y l O v 0 9 e J J 3 j W r 1 R a 1 y f V J o 3 e R x F 2 I c D q E I d z q A J V 9 C C N h C 4 h 0 d 4 h h f j w X g y X o 2 3 e W v B y G f 2 4 B e M 9 y + T n Z R j < / l a t e x i t > (b) w 3 b t I g a n 0 w 8 H h v h p l 5 b s i Z 0 r b 9 a S 0 t r 6 y u r R c 2 i p t b 2 z u 7 p b 3 9 j g o i S W i b B D y Q P R c r y p m g b c 0 4 + X 5 5 l r R r V e e 4 W r + q V x r 1 I o 4 S 2 A P 7 4 B A 4 4 B Q 0 w C V o g h Z A 4 B 4 8 g m f w Y j w Y T 8 a r 8 T Z t n T O K m R 3 w C 8 b 7 F y 9 O l k 4 = < / l a t e x i t > (w s-d , q s-d ) < l a t e x i t s h a 1 _ b a s e 6 4 = " O G Y I H 8 o U R 7 X a 8 2 m u t B K P P a 0 h 1 b w w 3 b t I g a n 0 w 8 H h v h p l 5 b s i Z 0 r b 9 a S 0 t r 6 y u r R c 2 i p t b 2 z u 7 p b 3 9 j g o i S W i b B D y Q P R c r y p m g b c 0

Figure 10 :
Figure 10: Sketches of 3 configurations of the line J s−d (3.6) and W s (P +,s )| M + 0 as consider in Theorem 3.8: (a) P +,s is the only critical point on M + ε and χ > 0; (b) P +,s is the only critical point on M + ε and χ < 0; (c) a center and a saddle on M + ε with w h,0 < 9/(2 + 9a) < w +,s .

Figure 11 :
Figure 11: Sketches of 3 (projected) 'skeleton structures' of stationary 2-fronts homoclinic to P 0 : (a) (b) Two examples of the skeleton structures Γ s-2f as considered in Theorem 3.11; (c) The extended skeleton Γ exts−2f of Theorem 3.12 and an associated higher order 2-front homoclinic with an additional full extra 'spatial oscillation' during its passage along M + ε .
9 B P c u F D E r V / k z r 9 x k g Z R 6 4 E L h 3 P u 5 d 5 7 v I g z p W 3 7 0 y o s L a + s r h X X S x u b W 9 s 7 5 d 2 9 t g p j S W i L h D y U X Q 8 r y p m g L c 0 0 p 9 1 I U h x 4 n H a 8 y W X q dx 6 o V C w U d 3 o a U T f A I 8 F 8 R r A 2 0 i 2 / V 4 N y x a 7 a G d A i c X J S g R z N Q f m j P w x J H F C h C c d K 9 R w 7 0 m 6 C p W a E 0 1 m p H y s a Y T L B I 9 o z V O C A K j f J T p 2 h I 6 M M k R 9 K U 0 K j T P 0 5 k e B A q W n g m c 4 A 6 7 H 6 6 6 X i f 1 4 v 1 v 6 5 m z A R x Z o K M l / k x x z p E K V / o y G T l G g + N Q Q T y c y t i I y x x E S b d E p Z C B c p 6 t 8 v L 5 L 2 S d U 5 r d a u a 5 X G T R 5 H E Q 7 g E I 7 B g T N o w B U 0 o Q U E R v A I z / B i c e v J er X e 5 q 0 F K 5 / Z h 1 + w 3 r 8 A e d W O H Q = = < / l a t e x i t > l u r 9 1 i 7 k g g X 8 t p y H u M T j y y Z A g K J X U 1 / e 6 5 5 A x e B M f J / 2 4 y 7 z g L g 6T p K 8 X L d O 1 X b v s G J Z Z q l o l p 6 q I b b u O I r Z p Z S i C O e p 9 / b 0 7 C F D E s C 8 R h U J 0 b C u U v R h y S R D F S a E b C R x C N I E j 3 F H U h w y L X p y d n x i H S h k Y w 4 C r 8 q W R q d 8 n Y s i E m D J P d T I o x + K 3 l 4 r / e Z 1 I D k 9 6 M f H D S G I f z R Y N I 2 r I w E i z M A a E Y y T p V B G I O F G 3 G m g M O U R S J V b I Q j h N 4 X y 9 / J c 0 S 6 Z d N i u X l W L t a h 5 H H u y D A 3 A E b O C C Gr g A d d A A C M T g A T y B Z + 1 e e 9 R e t N d Z a 0 6 b z + y C H 9 D e P g F e M 5 a d < / l a t e x i t > 0 p < l a t e x i t s h a 1 _ b a s e 6 4 = " D F s S I U t e x i t s h a 1 _ b a s e 6 4 = " p B k z c I c R / H B / b o U H q b s 2 r m 7 9 D c I = " > AA A B + X i c b V D L S s N A F J 3 U V 6 2 v q E s 3 w S K 4 C k k f i e 4 K b l x W s Q 9 o Y 5 h M J + 3 Q y Y O Z S b W E / I k b F 4 q 4 9 U / c + T d O 0 i K + D l w 4 n H M v 9 9 7 j x Z R w Y R g f S m l l d W 1 9 o 7 x Z 2 d r e 2 d 1 T 9 w + 6 P E o Y w h 0 U 0 Y j 1 P c g x J S H u C C I o 7 s c M w 8 C j u O d N L 3 K / N 8 O M k y i 8 E f M Y O w E c h 8 Q n C A o p u a p 6 d 5 s a m Z s O A y + 6 T + M s c 9 W q o d u m b d Y t z d B r T a N m N S U x T d u S x N S N A l W w R N t V 3 4 e j C C U B D g W i k P O B a c T C S S E T B F G c V Y Y J x z F E U z j G A 0 l D G G D u p M X l m X Y i l Z H m R 0 x W K L R C / T 6 R w o D z e e D J z g C K C f / t 5 e J / 3 i A R / p m T k j B O B A 7 R Y p G f U E 1 E W h 6 D N i I M I 0 H n k k D E i L x V Q x P I I B I y r E o R w n k O 6 + v l v 6 R b 0 8 2 6 3 r h q V F v X y z j K 4 A g c g 1 N g A h u 0 w C V o g w 5 A Y A Y e w B N 4 V l L l U X l R X h e t J W U 5 c w h + Q H n 7 B N j k l K g = < / l a t e x i t > W s (P +,s ) < l a t e x i t s h a 1 _ b a s e 6 4 = " T d a G Z H Y U F k l c q W b e r M y O 0 g R w 3 D 8 = " > A A A B 9 X i c b V D J S g N B E K 2 J W 4 x b 1 K O X x i B E l D C j w e U W 8 O I x i l k g mY S e T k / S p G e h u 0 c J w / y H F w + K e P V f v P k 3 9 k y C q P F B w e O 9 K q r q O S F n U p n m p 5 F b N M R k j I e 0 o 6 m P P S r t O L s 6 Q Q d a G S A 3 E L p 8 h T L 1 5 0 S M P S k n n q M 7 P a x G 8 q + X i v 9 5 n U i 5F 3 b M / D B S 1 C f T R W 7 E k Q p Q G g E a M E G J 4 h N N M B F M 3 4 r I C A t M l A 6 q k I V w m e L s + + V 5 0 j y p W K e V 6 k 2 1 V L u d x Z G H P d i H M l h w D j W 4 h j o 0 g I C A R 3 i G F + P B e D J e j b d p a 8 6 Y z e z C L x j v X 7 i 3 k j E = < / l a t e x i t > J 0,+ < l a t e x i t s h a 1 _ b a s e 6 4 = " R l 8 c B d 4 e r F Q y x 6 J d A I C n A u S 1 8 u k = " > A A A B 7 n i c b V D J S g N B E K 1 x j X G L e v T S G A R B C T M u 6 D H g R T x F M Q s k Y + j p 9 C R N e n q G 7 h o h D P k I L x 4 U 8 e r 3 e P N v 7 C w H T X x Q 8 H i v i q p 6 Q S K F Q d f 9 d h Y W l 5 Z X V n N r + f W N z a 3 t w s 5 u z c S p Z r z K Y h n r R k A N l 0 L x K g q U v J F o T q N A 8 n r Q v x 7 5 9 S e u j Y j V A w 4 S 7 k e 0 q 0 Q o G E U r 1 W 8 f M / f k e N g u F N 2 S O w a Z J 9 6 U F G G K S r v w 1 e r E L I 2 4 Q i a p M U 3 P T d D P q E b B J B / m W 6 n h C W V 9 2 u V N S x W N u P G z 8 b l D c m i V D g l j b U s h G a u / J z I a G T O I A t s Z U e y Z W W 8 k / u c 1 U w y v / E y o J E W u 2 G R R m E q C M R n 9 T j p C c 4 Z y Y A l l W t h b C e t R T R n a h P I 2 B G / 2 5 X l S O y 1 5 Z 6 W L u / N i + X 4 a R w 7 2 4 Q C O w I N L K M M N V K A K D Pr w D K / w 5 i T O i / P u f E x a F 5 z p z B 7 8 g f P 5 A 1 7 M j v 4 = < / l a t e x i t > J +,0 < l a t e x i t s h a 1 _ b a s e 6 4 = " C i U W / d n T S M i 2 G S x N T z s C h v 4 L N c o = " > A A A B 7 n i c b V D J S g N B E K 1 x j X G L e v T S G A R B C T M u 6 D H g R T x F M Q s k Y + j p 9 C R N e n q G 7 h o h D P k I L x 4 U 8 e r 3 e P N v 7

7 5 9 Figure 12 :
Figure 12: Sketches of 2 skeleton structures of multi-front patterns: (a) The stationary homoclinic 2-front gap pattern of Theorem 3.13.(b) The spatially periodic multi-front spot/gap pattern of Theorem 3.15 as (projected) closed orbit.

Figure 14 :
Figure 14: (a) A heteroclinic stationary 1-front pattern of (1.5);(b) Two traveling 1-front patterns connecting the bare soil state to a homogeneous vegetation state, one invading the bare soil (c < 0), the other the vegetation state (c > 0); (c) A homoclinic stationary 2-front spot pattern; (d) Evolution of the middles of the 2 interacting fronts of an evolving 2-front pattern.

Figure 18 :
Figure 18: (a) A small amplitude stationary spatially periodic solution generated by a Turing bifurcation (see Remark 2.2.(b) A pattern that is periodic in space and time that appeared by decreasing Ψ from the Turing pattern of Fig. 18a through a Hopf bifurcation (in time).

(B. 2 ) 4 )
By choosing δ and γ as in(1.3),we arrive at, use our freedom in α and β to simplify the B-equation and scale the factors of the B-and W B 2 -terms (to −1 and to +1 respectively) -which is achieved by the choices in(1.3),       B t = KE (KE − 1) 2 W − 1 B + W B 2 − W B 3 + B xx ,This is equivalent to (1.5) by definitions (1.4) and (1.6).
S per is given by (2.35) and its boundaries R Hopf and R hom by the upper, respectively lower, boundary of(2.35).