Criteria for the (in)stability of planar interfaces in singularly perturbed 2-component reaction-diffusion equations

We consider a class of singularly perturbed 2-component reaction-diffusion equations which admit bistable traveling front solutions, manifesting as sharp, slow-fast-slow, interfaces between stable homogeneous rest states. In many example systems, such as models of desertification fronts in dryland ecosystems, such fronts can exhibit an instability by which the interface destabilizes into fingering patterns. Motivated by the appearance of such patterns, we propose two versions of a 2D stability criterion for (transversal) long wavelength perturbations along the interface of these traveling slow-fast-slow fronts. The fronts are constructed using geometric singular perturbation techniques by connecting slow orbits on two distinct normally hyperbolic slow manifolds through a heteroclinic orbit in the fast problem. The associated stability criteria are expressed in terms of the nonlinearities of the system and the slow-fast-slow structure of the fronts. We illustrate and further elaborate the general set-up by explicitly working out the existence and transversal (in)stability of traveling fronts in a number of example systems/models. We analytically establish the instability of invading bare soil/vegetation interfaces against transversal long wavelength perturbations in several dryland ecosystem models and numerically recover fingering vegetation patterns counter-invading an invading desertification front.


Introduction
In this paper we consider a general system of 2-component singularly perturbed reaction-diffusion equations, for (x, y) ∈ R 2 , i.e. on the full unbounded 2-dimensional plane, with U (x, y, t), V (x, y, t) : R 2 × R + → R, F (U, V ; µ) and G(U, V ; µ) sufficiently smooth, and τ > 0, µ ∈ R m , parameters. Originally motivated by the fact that fronts connecting stable homogeneous states, also called coexistence states in ecology, have a direct ecological interpretation as invasion fronts -see for instance [16,19] -we investigate the existence and stability of planar interfaces between stable states in systems of reaction-diffusion equations. We focus on model (1.1) and assume throughout this paper that system (1.1) has (at least) two stable homogeneous background states, (U (x, y, t), V (x, y, t)) ≡ (Ū ± ,V ± ). Moreover, we assume that it is singularly perturbed: we consider 0 < ε 1 so that the diffusive spreading speed of the V -component in (1.1) is much larger than that of the U -component. Since pattern formation in ecosystems is typically driven by counteracting feedback mechanisms on widely different spatial scales [39], this is a very natural assumption in the setting of ecosystem models. Although ecosystem models certainly are not restricted to 2-component models (see section 5), there is a large variety of ecological models that are covered by (1.1) -see [1,2,11,16,17,19,27,31,33,42,43,48,50] and the references therein. Moreover, models of type (1.1) occur throughout the literature in a wide variety of (non-ecological) settings -see [10,12,20,23,24,25,35,36,44,45,46,49] and the references therein.
In the case of a scalar reaction-diffusion equation, W t = ∆W + H(W ; µ), (1.2) the assumption that (1.2) has two stable homogeneous rest states, W (x, y, t) ≡W ± immediately implies by phase plane techniques that there must be a monotonic front solution W (x, y, t) = W h (x − c * t) that travels with a well-defined, unique, critical speed c * in the (longitudinal) direction -which we identify with x -and that can be seen as a connection between the stable statesW ± : lim ξ→±∞ W h (ξ) =W ± (with ξ = x − c * t). Its longitudinal (spectral) stability, i.e. its stability with respect to perturbations that only depend on x, is determined by the linearized Sturm-Liouville problem, L w (ξ)w =w ξξ + c * w + H w (W h (ξ))w = λw.
It follows by the monotonicity of W h (ξ) that the translational eigenvalue is the critical eigenvalue, i.e. λ c = 0, and thus that the front is stable [29]. The stability of the front of (1.2) as a traveling planar interface in R 2 is determined by allowing for transversal perturbations through the inclusion of Fourier modes e i y in the linearization Ansatz, which yields the 1-parameter family of spectral problems that can be reduced to the longitudinal problem by replacing λ by λ + 2 . Since the parabolic curve λ c ( ) = − 2 has its maximum at = 0 it follows that the planar interface is also stable with respect to transversal perturbations and thus as traveling planar interfacial solution of (1.2) on R 2 . In this paper we consider the question whether similar general insights on the existence and stability of planar interfaces between stable states can be deduced for non-scalar reaction-diffusion equations.
In the 2-component setting of (1.1), introducing traveling coordinate ξ = x − ct and u(ξ) = U (x, y, t), v(ξ) = V (x, y, t), reduces (1.1) to a system of coupled second order (ordinary differential) equations in ξ, which has, by taking ε → 0, the 2-parameter family as fast reduced limit. The critical points of (1.4) are determined by the manifold {(u, v) : F (u, v) = 0}. Away from the (degenerate) transition points at which F u (u, v) = 0, F (u, v) = 0 determines J ≥ 1 branches that can be written as graphs {u = f j (v)} (j = 1, 2, ..., J) with f j (v) such that F (f j (v), v) ≡ 0. The background states (Ū ± ,V ± ) of (1.1) must correspond to points on one of these branches (since F (Ū ,V ) = G(Ū ,V ) = 0). We define two specific branches {u = f ± (v)} by selecting j's such that U ± = f ± (V ± ). Next to our assumption that the background states (Ū ± ,V ± ) are stable as solutions of (1.1), we impose a second assumption that underlies our analysis: we assume that (Ū ± ,V ± ) are on different branches of {(u, v) : F (u, v) = 0}, i.e. that f + (v) ≡ f − (v). In the setting of the singularly perturbed 4-dimensional dynamical system associated to (1.3) these assumptions together imply that the background states correspond to saddle points for the slow flows on two different 2-dimensional normally hyperbolic slow manifolds M ± that are at leading order in ε given by M ± 0 = {(u, p) : u = f ± (v), p = 0} -see section 2.1. The traveling interfaces (U h (ξ), V h (ξ)) we study here correspond to heteroclinic orbits between the saddles associated to (Ū ± ,V ± ) on these slow manifolds. The fact that M − = M + implies that these interfaces/orbits must consist of 3 parts: a slow part along M − , followed by a fast jump from  [19] (a) The manifold {F (u, v) = 0} (in red) with normally hyperbolic submanifolds M ± 0 indicated by continuous red curves (and non-normally hyperbolic parts dashed), the stable background states (Ū ± ,V ± ) as intersections of the normally hyperbolic slow manifolds and the manifold {G(u, v) = 0} (in blue) and the projection of the singular slow-fast-slow front (in green) that makes a fast jump between two well-defined points (u − * , v * ) -with here u − * = 0 -and (u + * , v * ) along the heteroclinic solution u * (ξ) of (1.4). (b) The components U h (ξ) and V h (ξ) of the front. Note that the fast jump u * (ξ) is monotonic between its endpoints u ± * , but that U h (ξ) typically is not monotonic and that V h (ξ) only varies slowly (in X = εξ).
M − = M + , concluded by a final slow part along M + (see Fig. 1).
The details of the construction of these slow-fast-slow orbits are given in section 2.1 for τ = O(1) (and in section 3.1 for τ = O(ε) -see below). For the formulation and interpretation of the upcoming (in)stability criteria we need to provide somewhat more information on the (leading order) fast jump. The analysis of the (slow) flows on M ± 0 yields an explicit critical value v * (that does not depend on c) at which the jump must occur -see Fig. 1. Planar fast reduced system (1.4) -with v 0 = v * -has (by construction) two saddle points (u ± * , 0) with u ± * = f ± (v * ), these are connected by a monotonic heteroclinic orbit u * (ξ) -see again Fig. 1 -for a uniquely determined value (cτ ) * , which thus determines the speed c = c * = (cτ ) * /τ : u * (ξ) is the leading order approximation of the fast jump and c * determines the (leading order) speed of the full interface (U h (ξ), V h (ξ)). Note that thus neither u * (ξ), u ± * nor the product c * τ depends directly on τ (at leading order in ε): these expressions are solely determined by v * .
In this paper we do not consider the analytical details of the stability of the front (U h (ξ), V h (ξ)) with respect to longitudinal perturbations, i.e. perturbations that only depend on x ∈ R (or ξ ∈ R). Although this may be a nontrivial technical endeavour for a given system, we note that the stability can be established by the methods developed in [1,6,12,14,35,36,49] and the references therein. Thus, we assume that the front is longitudinally stable and thus that the translational eigenvalue is the most critical eigenvalue: λ c = 0 (as in the scalar case). The central question considered in this paper is: Under what condition(s) are longitudinally stable slow-fast-slow fronts of system (1.1) (un)stable as planar traveling interfaces in R 2 , i.e. when are these fronts (un)stable with respect to transversal perturbations?
We do not pursue this question in its full generality, although this again can also be done for a given explicit system: we focus on obtaining general (in)stability conditions against long wavelength transversal perturbations. More precise: given a longitudinally stable interface, we introduce ∈ R by the factor e i y in the spectral stability Ansatz as the wavenumber of transversal perturbations so that λ c becomes a function of and consider | | 1. Our main general (in)stability criteria concern the sign of λ 2,c , the planar desertification front between a bare soil state (in yellow) and a homogeneously vegetated state (in green) in ecosystem model (4.11) of [19] (with µ ∈ R 5 and µ = (3.5, 1.1, 3.2, 1.0, 0.4), ε = 0.04). The simulations employed finite differences for spatial discretization and MATLAB's ode15s routine for time integration. The interface is unstable against transversal long wavelength perturbations by criterion (1.6) for all ecologically realistic parameter combinations (and ε sufficiently small) -see section 4.1.2. As first observed in [19]: the interface is counter-invaded by fingering vegetation patterns.
leading order coefficient of the expansion of λ c ( ) for | | 1: λ c ( ) = λ 2,c 2 + O( 4 ) (since λ c is a function of 2 by the reflection symmetry of (1.1) in y). Our approach is based on the explicit construction of the λ = 0-eigenfunction of the adjoint of the (matrix) operator associated to the spectral stability of the fronts (U h (ξ), V h (ξ)) -which is based on the singular slow-fast-slow structure of the fronts (see section 2.2).
In the so far considered case τ = O(1) we define, Our first (in)stability criterion reads, where we note that this condition does not depend on the value of τ (as long as it is O(1) with respect to ε). Based on the monotonicity of u * (ξ) and the fact that one typically knows the signs of v * and u * (ξ), the simplicity of this criterion allows us to immediately draw conclusions on the instability of traveling planar interfaces in the dryland ecosystem models considered in [1,16,19,27], and thus on the counter-invasion of fingering vegetation patterns against an invading desertification front [27] -see Fig. 2 and section 4.1. In fact, if (1.1) is of activator-inhibitor type [34] with (for instance) V (x, y, t) as activator -i.e. F v (u, v) > 0 -and U (x, y, t) as inhibitor -G u (u, v) < 0 -it follows from (1.5) and (1.6) that the (in)stability of (invading) planar interfaces against transversal long wavelength perturbations typically is decided by the signs ofŪ Clearly, λ 2,c = −1 < 0 in the scalar case of (1. below -that λ 2,c > 0 occurs as natural as λ 2,c < 0: a longitudinally stable front is unstable against transversal long wavelength perturbations under a simple, explicitly verifiable, condition on the vector field (F, G) of (1.1) -see section 4. In this paper we do not go into the nonlinear nature of this instability. Typically, it gives rise to 'fingering patterns' [19,23,45]. This type of growth pattern originates from fluid dynamics and is associated to the sideband mechanism [7], i.e. the instability with respect to perturbations with wave numbers |k| 1, while k = 0 itself is (marginally) stable [8,10] -as is the case here (with in the role of k). Although understanding the onset of counter-invading fingers - [19] and Fig. 2 -formed the original motivation for the research presented in this paper, we note that λ 2,c > 0 does not necessarily induce the formation of fingers: see the simulation of Fig. 6 in which a transversally unstable invading vegetation front in an ecosystem model based on [1,16] exhibits evolving bounded cusps reminiscent of Kuramoto-Sivashinsky dynamics associated to a sideband instability [13,26] (however, in this model invading desertification fronts also exhibit counter-invading vegetation fingers as in [19] and By rescaling space -(x,ỹ) = (εx, εy) -and defining D = ε 2 , we may interpret the parameters (τ, D) as measures for the relative magnitudes of the evolution rates of U and V in time (τ ) and space (D). The choices sofar made -τ = O(1) and D 1 -are based on the ecosystem models that originally motivated our investigations. However, in the non-ecological literature on models of the type (1.1) often the scaling [20,23,24,25,36,44,45] and the references therein. In fact, in [23,24,25] the parameters (ν, η) are introduced, with ν = ε 2 , η = τ /ε in the notation of (1.1). In the notation of [23,24,25] we have so far thus considered 0 < ν 1, η 1 , while η = O(1) (with respect to ε or ν) is the critical -richer -scaling. In other words, τ = O( √ D) is a 'significant degeneration' in the terminology of [15]. Therefore, we introduceτ by setting τ = ετ and set up an existence/stability analysis similar to the τ = O(1) case considered so far.
For τ = O(ε), the existence of the planar fronts can be approached by the same geometrical set-up as for τ = O(1). However, the 'algebra' becomes more involved: for τ = O(1) one can first determine v * by only considering the slow flows on the manifolds M ± 0 ; c * is determined by the fast field, using v * as (given) input. As was already noted, the fast heteroclinic jump u * (ξ) in (1.4) occurs for a given value (τ c) * = τ c * ; thus as τ decreases, c * increases: c must be scaled asc/ε for τ = ετ . However, in that casec will appear as O(1) term in the slow reduced flows: as for τ = O(1),ṽ * can still be determined, however, it becomes a function ofc. The same (geometric) procedure by which c * is determined for given v * , now determinesc as function of v * (c):ṽ * andc are determined simultaneously by algebraic relations (see section 3.1). Taking τ = O(ε) has a similar impact on the stability analysis: the essence of the method is not affected, but working out the details becomes more involved. As a consequence, the criterion for the sign ofλ 2,c , i.e. for the (in)stability of the associated planar interface against transversal long wavelength perturbations, changes from (1.6) into, sign(λ 2,c ( µ;τ )) = − sign(F * ( µ;τ )) × sign(M * ( µ;τ )) (1.7) whereF * is completely similar to F * (cf. (1.5) and (3.8)), but whereM * ( µ;τ ) is a somewhat more involved Melnikov-type expression (see (3.10)).
There is a significant difference between criteria (1.6) and (1.7): unlike (1.6), (1.7) depends explicitly onτ . To better understand the impact ofτ , we deduce a geometrical interpretation of the implication of a change in sign of Melnikov expressionM * ( µ;τ ). To do so, we restrict ourselves for simplicity to stationary fronts -as is common in studies of (1.1) (see for instance [23,35,45]). We show in section 3.2 that a sign change ofM * ( µ;τ ) initiates a bifurcation into traveling waves: ifτ * is such thatM * ( µ;τ * ) = 0, then (generically) a pair of counter-propagating traveling fronts bifurcates off from the stationary front asτ passes throughτ * . Note that this implies that an expression obtained in the context of the stability analysis of planar interfaces against transversal perturbations has an interpretation in terms of the existence analysis of one-dimensional fronts -see section 5.
In sections 4.2 and 4.3, we consider versions of the example systems formulated and studied in [23,24,25,44,45]. We show how criteria (1.6) and (1.7) immediately provide analytical insight in the instability of planar interfaces for (x, y) ∈ R 2 . We test the outcome of our analysis against numerical evaluations of the spectral curve λ c ( ) in the context of a FitzHugh-Nagumo type model based on [23,24,25] in section 4.2 and compare our results in section 4.3.1 to those of [44,45] in which cylindrical spatial domains are considered. Moreover, in section 4.3.2 we apply (1.7) to the counter-propagating fronts initiated by the bifurcation into traveling waves in a special case of the model in 4.2 to show that these interfaces inherit the instability against long wavelength transversal perturbations from the original stationary interface, whileM * ( µ;τ ) has changed sign for the stationary front. We conclude the paper with a section in which we discuss both the ecological and mathematical implications of our findings.
Finally, we briefly remark on the writing style of the paper. Based on the present state of the art literature on geometric singularly perturbation theory (for the existence problem) and the analysis of singularly perturbed spectral operators (for the stability), all our findings can be established rigorously. However, for transparency of presentation we have chosen not to formulate our main results as theorems, and thus not to provide all details of the proofs of these theorems. Instead, we have chosen to present our work in the present hands-on setting, with a central role for section 4 in which we explicitly show how our methods work in the context of a number of diverse and explicit problems. In a way, our paper can be seen as a 'user guide': we have written the paper such that it can directly be used for the analysis of (the stability of) planar interfaces in explicit systems of the type (1.1).
Remark 1.1. We only consider F (U, V ), G(U, V ) that do not explicitly depend on ε. This is a natural assumption from the ecological point of view [11], however, there is a large body of literature on singularly perturbed reaction-diffusion equations in which F (U, V ) and G(U, V ) do depend on ε -such as the classical Gray-Scott/Gierer-Meinhardt-type models (cf. [11,49] and the references therein).

The case τ = O(1)
In this section we first set up the (geometrical) construction of the travelling slow-fast-slow front solutions of (1.1)/(1.3) and next consider its stability as 2-dimensional planar interface against against transversal long wavelength perturbations. For clarity, we list the assumptions that will be imposed on the existence problem in this section (and the next): • (E-I) There are two homogeneous background states (Ū − ,V − ) and (Ū + ,V + ) that are stable as solutions of (1.1).
• We refer to the additional assumption (E-IV) in upcoming section 2.1. We note that similar front patterns as the ones (to be) constructed here -and in section 3.1 for τ = O(ε) -have been established in [20,35,44,45,46], also for general systems of the type (1.1), but typically under stronger/more explicit conditions on the reaction terms (F (U, V ), G(U, V )) of (1.1) than in the present setting -see also section 4.3.1 for some more details. However, the precise nature of the conditions imposed in the literature are for a large part for practical reasons, the methods used (and developed) can typically also be applied to more general systems (see again section 4.3.1). We also refer to the concluding paragraph of (upcoming) section 2.1 and to Remark 2.1 for an extended discussion on the assumptions under which slow-fast-slow fronts can be constructed.
There is one crucial additional assumption for the subsequent stability analysis: • (S-I) The constructed slow-fast-slow heteroclinic traveling front is stable against longitudinal perturbations, i.e. perturbations that only depend on x.
We refer to section 2.2 for a more precise statement and a brief discussion.

The construction of slow-fast-slow interfaces
There are 2 equivalent ways to write (1.3) as a 4-dimensional spatial dynamical system, with X = εξ the slow (spatial) coordinate. The fronts/interfaces of interest correspond to heteroclinic orbits γ h (ξ) in (2.1) that connect its critical points (Ū − , 0,V − , 0) and (Ū + , 0,V + , 0) for a critical value c * of the speed c. The orbits γ h (ξ) can be constructed explicitly by the methods of geometric singular perturbation theory. To do so, we reintroduce the two slow manifolds M ± 0 , , v ± ± ) = 0 (where we note that M ± 0 may be unbounded, i.e. that v ± + may be ∞ and/or v ± − may be −∞). As stated in (E-III), it is necessary for the construction of γ h (ξ) to assume that the projections of M ± 0 on the v-axis overlap, i.e. that Naturally, (Ū ± , 0,V ± , 0) ∈ M ± 0 are critical points of the full system (2.1) and where we note that these (planar) systems are integrable (Hamiltonian) under the implicit assumption that The assumption that both (Ū ± ,V ± ) are stable as homogeneous backgrounds states of (1.1) has a direct impact on the nature of the critical points (Ū ± , 0,V ± , 0) of (2.1), of the critical points (V ± , 0) of (2.3) and of the slow manifolds M ± 0 . To see this, we note that the states (Ū ± ,V ± ) are stable if the eigenvalues λ ± 1,2 (k, ) of (2.4) satisfy Re λ ± 1,2 (k, ) < 0 for all k, ∈ R (whereF ± u = ∂F ∂U (Ū ± ,V ± ), etc., and k, are the (Fourier) wavenumbers in the x, y-directions) [10]. For 0 < ε 1, this is the case if (and only if), [11], where we note that the last 2 conditions come from the observation that the λ-independent term in (2.4) must be positive for all k, ∈ R (consider k = = 0 and Moreover, the critical point (Ū ± , 0) of the fast reduced system (1.4) at the level v 0 =V ± also must be a saddle (F ± u < 0). Since 0 < ε 1 this implies that the eigenvalues λ ± j , j = 1, ..., 4 associated to the critical points (Ū ± , 0,V ± , 0) of (2.1) satisfy 3)). Moreover, the various branches u = f j (v) of the manifold F (u, v) = 0 are separated from each other by transition points at which F u (u, v) = 0 ( Fig. 1(a) , v ± ± ) = 0 in the case of the curves f ± (v). As a consequence, the sign of F u (f ± (v), v) does not change along M ± 0 : the stability of (Ū − , 0,V − , 0) implies that F u (f ± (v), v) < 0 on M ± 0 . This means that M ± 0 are both normally hyperbolic, so that they persist as invariant manifolds M ± for the flow of (2.1) (for ε sufficiently small [18,28]) -see also Remark 2.1.
(V − , 0) on M − 0 (Fig. 3). Likewise, it will follow the stable manifold W s ((V + , 0)) of the saddle (V + , 0) on M + 0 during its second slow part. Neither the v-nor the v X = q-component of γ h (ξ) will change (at leading order in ε) during the fast jump (1.4). Thus, the γ h -skeleton can be constructed if the unstable manifold 0)), follows the fast field, and 'touches down' on M + 0 again at (v * , q * ). Or more precisely, γ h (ξ) follows the fast flow from the point (u The speed c * of the traveling interface is determined by the passage through the fast field (at v = v * , q = q * ). In general, i.e. for general c, the unstable manifold of the saddle (u − * , 0) of (1.4) will not be connected to the stable manifold of (u − * , 0) -where we recall that the assumed stability of background states (Ū ± ,V ± ) implies that the critical points of (1.4) are saddles. The situation is now completely similar to that of constructing heteroclinic traveling fronts for scalar equation (1.2) -in fact, by (re)defining F (u, v * ) as H(u), (1.4) is identical to the (spatial) traveling wave planar dynamical system associated to (1.2). Thus, for a given v * we can (by phase plane techniques) determine a unique value c * of c for which there is a (planar) heteroclinic connection (u * (ξ), p * (ξ)) between the saddles (u ± * , 0). By construction, p * (ξ) = u * ,ξ (ξ) is either strictly positive or strictly negative: (u * , p * ) forms the (monotonic) fast component of the above skeleton structure (Fig. 1). (It should be noted that although the fast jump u * (ξ) is monotonic, the u-component u h (ξ) of the full slow-fast-slow connecting front typically is not monotonicsee Fig. 1

.)
In the setting of a given model it may be a technical undertaking to explicitly set up the construction of the skeleton of the heteroclinic orbit Here, in the general setting, we assume that this has been done. By definition, the piecewise approximations ) we obtain similar leading order approximations of u h (ξ) in the slow fields. Naturally, u h (ξ) is approximated by u * (ξ) during the fast jump -see Figs. 1 and 3. Note that p h (ξ) = u h,ξ (ξ) and q h (X) = v h,X (X) by (2.1). From the geometrical point of view, the construction of the full slow-fast-slow front is now based on a transversal intersection of the (3-dimensional) stable and unstable manifolds W u (M − 0 ) and W s (M + 0 ) (in the ε → 0 limit) and it thus follows by the methods of geometric singular perturbation theory that asymptotically close to this skeleton there indeed is an orbit -the orbit γ h (ξ) -that is a heteroclinic connection between the critical points (Ū − , 0,V − , 0) and (Ū + , 0,V + , 0). This orbit corresponds by construction to the (planar) traveling interface Wrapping up the existence analysis, we return to the conditions under which the traveling interface (U h (ξ), V h (ξ)) can be constructed. Additional to (E-I) -(E-III) stated at the beginning of this section, that did not consider the slow flows on M ± 0 , we now add the fourth condition: Note that this certainly is not 'automatically' the case, on the other hand, W u ((V − , 0)) ∩ W s ((V + , 0)) may consist of several points (v * ,i , q * ,i ) -see Fig. 3 -so that there may exist several distinct connecting fronts for a given choice of parameters µ. For a given explicit system, there also is a 'hidden' condition under which the skeleton structure can be built: the take off/touch down points (u ± * , 0, v * , q * ) must be on . Although this must be the case by construction, in a given system a manifold M ± 0 may extend beyond one of its end/transition points (f ± (v ± ± ), v ± ± ) as non-normally hyperbolic manifold -see for instance (4.13) in section 4.1.2 -and then it may happen that (the projections of) In that case at least one of the points (u ± * , 0, v * , q * ) will not be a saddle for the fast flow (1.4) so that the fast jump between M ± 0 cannot be made (by a monotonic fast connection u * (ξ) -see Remark 2.1).
We conclude that the existence problem for (1.1) differs from that of (1.2) in the sense that there always is a uniquely determined traveling connecting (monotonic, planar) interface in (1.2) (when there are two stable background states), while we found that this certainly is not the case for singularly perturbed 2-component model (1.1). Nevertheless, the conditions on F (U, V ) and G(U, V ) are quite mild and the slow-fast-slow interfaces (U h (ξ), U h (ξ)) can be constructed for many models -as we shall see in section 4.
Remark 2.1. As in [11], we consider a slow manifold of (2.1) to be normally hyperbolic if it is normally hyperbolic for all c with cτ = O(1). The dashed parts of {F (U, V ) = 0} in Fig. 1(a) are normally hyperbolic for c = 0, but elliptic for c = 0, and thus not considered as normally hyperbolic in the present analysis.
Orbits that (for c = 0) either take off or touch down from branches of {F (U, V ) = 0} that are not normally hyperbolic in the present sense, i.e. that are elliptic for c = 0, will oscillate and thus correspond to fronts with non-monotonic fast jumps. Such fronts are expected to be (longitudinally) unstable. Note that by the above arguments, the critical points associated to the stable background states (Ū ± ,V ± ) necessarily must be on invariant manifolds that are normally hyperbolic for all c = O(1).
Remark 2.2. The savanna ecosystem model of [48] only has one (ecologically relevant) normally hyperbolic slow manifold that typically contains critical points (for the slow flow) associated to several different stable homogeneous background states [11]. Thus, the model of [48] violates assumption (E-II) and cannot exhibit the slow-fast-slow interfaces considered here. However, if (E-II) does not hold, there will be heteroclinic connections between saddle points on the slow manifold that do not depart from the slow manifold. The existence, stability and bifurcations of such 'fully slow' localized patterns is studied in [11].

Transversal instability of the interface
To determine the (spectral) stability of the traveling interface ( ) as solution of (1.1), we consider perturbations of the form, and derive the linearized system, which can equivalently be written as where the 2 × 2 matrix operator L(ξ) is given by with Sturm-Liouville operators We know by assumption (E-I) on the stability of the endstates (Ū ± ,V ± ) of the interface that the essential spectrum of L(ξ) is in the stable complex half-plane, moreover we know by the translational invariance of (1.1) that λ(0) = 0 is an eigenvalue of (2.8) As stated in the beginning of this section (and in the Introduction), we also impose assumption (S-I): the interface is stable against 1-dimensional perturbations, i.e. longitudinal perturbations that do not depend on y. In other words, we assume that all elements of the discrete spectrum of L(ξ) have Re(λ(0)) < 0, except for the -under this assumption -critical translational eigenvalue λ c (0) = 0. Moreover, we additionally impose the (generic) assumption that the dimension of the kernel of L is equal to 1, i.e. that the parameters µ of (1.1) are chosen away from their (co-dimension 1) bifurcation sets -see also section 5. Naturally, these assumptions need to be validated in the setting of an explicit model -something that typically can be done by exploiting the singularly perturbed nature of (1.1) by the methods of [1,6,12,14,35,36,49] and the references therein (although one should not underestimate the technicalities involved).
To analyze the impact of the y-direction, and thus the possibility of transversal instabilities, we consider the -dependence of the critical eigenvalue λ c ( ). As explained in the Introduction, we only consider long wavelength perturbations along the interface, i.e. 0 < 1. Recalling that the spectral problem (2.7)/(2.8) varies in 2 instead of , we expand By substitution of (2.11) into (2.8) we find at O( 2 ) the inhomogeneous equation, Since L(ξ) has a nontrivial kernel, the inhomogeneous term needs to satisfy a solvability/Fredholm condition [29]. Clearly L is not self-adjoint (and neither are its diagonal elements L u and L v (2.10)). We therefore define the adjoint operator L A (ξ), which yields an explicit expression for λ 2,c , (2.14) The planar interface (U h (ξ), V h (ξ)) is unstable with respond to long wavelength perturbations along the interface if λ 2,c > 0. It follows from (2.14) that we can explicitly characterize this potential destabilization mechanism by determining (a leading order approximation in ε of) (ū A c (ξ),ū A c (ξ)).
To do so, we rewrite L A (ū,v) = (0, 0) in a form similar to (2.7) and split this equation up into two slow parts and one fast part. The slow parts are at leading order given by, where we have approximated the slow parts (u h (X), v h (X) by (u ± (X), v ± (X)) and suppressed the u-, s-subscripts (i.e. u ± (X) is either = u − u (X) or u + s (X), etc.). The fast part is (at leading order) given by ≡ v * as leading order approximations during the fast jump. By taking ε → 0 in theū-equation of (2.15) we see thatū is determined byv in the slow fields, in the slow fields (as ε → 0), which implies by (2.3) thatv ± (X) = α ± v ± X (X) for some constants α ± . Taking ε → 0 in the fast field yields thatv(ξ) ≡v * , a constant, and that Since the kernel of L u is spanned by u * ,ξ this implies that (at leading order in ε). Thus, under the (generic) assumption that As a consequence, L A uū = 0 (at leading order), which implies thatū for some α * . By matching over the fast field, we can eliminate 2 of the 3 free factorsᾱ ± and α * . First we note thatv must be continuous over the fast field so thatᾱ 3) it also follows thatv * = εᾱq * . Recalling thatv = O(ε), we measure the accumulated change inv ξ over the fast field by integratingv ξξ , (2.18), which must bridge the jump in thev ξ 's coming from the slow flows, ). Now we can evaluate (2.14), although we need to be a bit more precise about our choice for the fast field, In contrast, the integral over u h,ξū A c is dominated by its fast part and is of O(1), (2.18). Thus, we find by (2.14) and (2.19), (at leading order in ε), which indeed yields (1.6) (with F * and G * as in (1.5)).
3 Taking τ asymptotically small: transitions at τ = ετ In the introduction, the impact of having τ = O(ε) was briefly sketched. This sketch will be worked out in more detail in the beginning of section 3.1. Now that we have derived the leading order approximation of λ 2,c (2.22), we can also sketch how the stability analysis of section 2.2 is affected by setting τ = ετ . A crucial ingredient in the derivation of (2.22) is the observation thatv A c must be O(ε), which follows from the application of the leading order solvability condition (2.17). However, as τ becomes O(ε), this condition is satisfied forv A c = O(1). As a consequence, the integral over v h,ξv A c is no longer O(ε), which has a significant impact on (2.22) and thus on (1.6), as we shall see in the upcoming section.
Unlike for τ = O(1), i.e. criterion (1.6), instability criterion (1.7) -withM * (τ ) to be derived in the upcoming subsection -has become an explicit function ofτ . Thus, the sign of λ 2,c may change asτ is varied. In section 3.2 we show that such a change of sign generates a bifurcation.
Remark 3.1. The three subconditions of condition (2.5) that establish the stability of the background states (Ū ± ,V ± ) of (1.1) reduces to two for τ = ετ , : These conditions are equivalent to the condition that the associated critical points of (Ū ± , 0,V ± , 0) of (2.1) correspond to saddle points for the slow reduced flows (2.3) of the slow manifolds M ± 0 . In other words, if τ = ετ then each saddle point for a slow flow on (normally hyperbolic slow manifold) M ± 0 necessarily corresponds to a stable homogeneous background state (which is not the case for τ = O(1)).

Existence and transversal (in)stability
As argued in the Introduction, we follow τ = ετ by setting c =c/ε, so that the 'friction term' in (1.4) remains O(1), For a given v 0 there thus still is a well-defined value ofc = C(v 0 )/τ for which there is an heteroclinic orbit between the saddles (f ± (v 0 ), 0) of (3.1), i.e. for which the fast system enables a jump between the slow manifolds M − 0 and M + 0 . However, the slow reduced problem is no longer integrable in this case, For the stability analysis we again assume that the front (Ũ h (x −c * t/ε),Ṽ h (x −c * t/ε)) is longitudinally stable. To investigate the stability with respect to y-dependent perturbations, we definẽ and derive the equivalent of (2.14), (cf. (2.15)) in the slow fields and of, v ± =αṽ ± X ec * X ,ū =α * ũ * ,ξ (ξ) ec * τ ξ in the slow fields and confirm that indeedα = O(1) (and not O(ε) as for τ = O(1)). The matching procedure also proceeds as for τ = O(1) and yields an outcome almost identical to (2.19) To further simplify notation, we introducẽ 9) and the Melnikov-type expression introduced in the Introduction, However, in a given system it may not be possible to takeτ sufficiently small, since a bifurcation may take place asτ is decreased. The origin of this bifurcation may be geometrical: it follows from (3.1) that |c| must become large asτ becomes small to have a jump through the fast field. However, the intersection of unstable manifold W u ((V − , 0)) and stable manifold W s ((V + , 0)) of the projected flows on M ± 0 may only persist up to a certain maximal value of |c| -see Fig. 4(b): the interface may for instance have 'disappeared' in a saddle-node bifurcation. Moreover, in the case that we know that the interface is unstable with respect to transversal perturbations for τ = O(1), i.e. if sign(F * ) = sign(G * ) forτ sufficiently large (1.6), then the sign ofλ 2,c may also change at the pole of leading order approximation (3.11) associated to M * =F * Ĩs +τG * Ĩf crossing through 0 for decreasingτ (3.11). We shall consider this in the upcoming section. (Here it should be noted that a pole in the leading order approximation (3.11) does not imply that λ 2,c itself is non-smooth: one needs to perform a higher order analysis to resolve the leading order character of λ 2,c near the zero of the denominator in (3.11).)

Bifurcation into traveling waves
As was already noted in the Introduction, the expressionM * ( µ;τ ) (3.10) has a Melnikov function-like character, especially forc * = 0 (see (3.8), (3.9)). Therefore, it is natural to expect that a zero ofM * ( µ;τ ) may be related to a 'structural geometrical change' in the existence analysis. In this subsection we consider the special case of stationary fronts, i.e. fronts for which the existence problem does not depend onτ . We show that in this case a change of the sign ofM * ( µ;τ ) indeed has a direct impact on the existence problem: it initiates a bifurcation into traveling waves, i.e. asM * ( µ;τ ) crosses through 0 -for instance as function ofτ -a pair of counter-moving traveling fronts branches off from the stationary fronts.
Existence problem (2.1) with τ = ετ and c =c/ε has a stationary front solution if (3.1) has a heteroclinic orbit forc = 0 for a value of v 0 =ṽ * which also is the v-coordinate of the intersection of the unstable manifold W u ((V − , 0)) and the stable manifold W s ((V + , 0)) of the projected flows (3.2) on M ± 0 withc = 0. In other words, C(ṽ) and V(c) of (3.3) must be so that C (V(0)) = 0, which imposes a co-dimension 1 condition on the parameters µ of (1.1) -see below. We assume that µ indeed has values so that C (V(0)) = 0 and thus use the fact that (1.1) has a stationary front solution for any τ = ετ > 0. Given this set-up, we ask the question: under which additional condition can there be slowly traveling fronts, i.e. fronts traveling with speedc =ĉδ with (artificial) small parameter 0 < ε δ 1 (that is independent of ε), such that these fronts merge with the stationary front as δ → 0?
We first consider the fast problem (3.1). We have assumed that there is a heteroclinic orbitũ * (ξ) between the critical points (ũ − * , 0) and (ũ + * , 0) -withũ ± * = f ± (ṽ * ) -in the integrable caseĉ = 0. Thus, we define the family of Hamiltonians,H parameterized by v 0 . The assumption on the existence of a stationary front implies see also [35,45] and note that (3.13) is the explicit version of the expression C(ṽ * ) = 0 (3.3). Naturally, one expects that if a traveling front appears from the stationary front, it will not make its fast jump at exactly the same value ofṽ * , therefore we introducev 1 byṽ * (c) =ṽ * + δv 1 (withv 1 =v 1 (ĉ) and recall that c = δĉ). As a consequence, the u-coordinates of the 2 saddle points (ũ ± * (c), 0) = (f ± (ṽ * (c)), 0) of (3.1) also undergo an O(δ) shift, where we note thatṽ * (0) =ṽ * andũ ± * (0) =ũ ± * . By (3.12) and (3.13) we have, Thec = 0-heteroclinic jump between (f − (ṽ * (c)), 0) and (f + (ṽ * (c)), 0) thus needs to accumulate a total change in 'energy' of Naturally, at v 0 =ṽ * (v),H f,ξ = −δĉτ u 2 ξ (3.1), (3.12), so that whereũ * (ξ) is the heteroclinic orbit atṽ * =ṽ * (0). Thus, we conclude by (3.8), (3.9) that, (3.14) Since the slow reduced flows on manifold M ± 0 are also integrable as δ → 0, we define the slow Hamiltonians, so thatH ± s (v, q) = 0 at the saddle points (V ± , 0) and at the intersection (ṽ * ,q * ). Note thatṽ * is determined by, which thus yields the explicit version of V(0) (3.3): together, (3.13) and (3.16) determine the co-dimension 1 condition on the parameters µ, C(V(0)) = 0. Since againH ± s,X = −ĉδq 2 (3.2), we observe by a similar Melnikov argument as above that the unstable manifold W u ((V − , 0)) on M − 0 and the stable manifold W s ((V + , 0)) on M + 0 of the (δ-)perturbedĉ = 0-flows intersect the lines v =ṽ * (⊂ M ± 0 ) at the level sets Thus, the O(δ) corrections to the q-coordinates of the intersections of W u ((V − , 0)) and W s ((V + , 0)) with the vertical {v =ṽ * } are determined by,q As a consequence, we know that O(δ) close to (ṽ * ,q * ), W u ((V − , 0)) and W s ((V + , 0)) are up to O(δ 2 ) corrections given by the linear approximations, Naturally, these two lines intersect and their intersection gives the leading order approximation of the intersection of (the projections of) W u ((V − , 0)) and W s ((V + , 0)) for the perturbed flows (3.2) withc = δĉ. By the above definition, this intersection has v-coordinateṽ * (ĉ) =ṽ * + δv 1 and we conclude by (3.8), (3.9) that,v Note that the analysis so far only provides a first approximation: the lineτ =τ * = −F * Ĩs /G * Ĩf in the (τ ,c)-plane only gives a leading order vertical approximation of the (parabolic) curve of bifurcating traveling waves (again in the (τ ,c)-plane) that is expected for a non-degenerate bifurcation. More specifically, the orientation of this curve determines whether the bifurcation is sub-or supercritical: with the present analysis we do not know yet the precise nature of the bifurcation (and under which condition(s) it is non-degenerate). This more detailed information can be determined explicitly by a higher order asymptotic analysis (see for instance [11]), however, we refrain from going into this here. Moreover, it is natural to ask whether the bifurcating traveling interfaces will be stable against transversal (long wavelength) perturbations or not: clearly the sign ofM * changes for the basic stationary pattern asτ passes through τ * , but the sign ofM * is a priori not clear for the bifurcating traveling fronts. Such a higher order analysis can also be performed as we show in section 4.3.2 in the context of a specific model.

Ecosystem models and other example systems
The analysis of this paper was inspired by the insights of [37] on the positive impact of spatial patterns on the resilience of ecosystems -see section 5 -and of [19] on potentially reversing desertification through a fingering instability of an invasion front. Therefore, we exemplify the general approach by first considering two explicit ecosystem models in sections 4.1.1 and 4.1.2. However, the ecosystem models considered typically do not have a (small) τ -pre-factor so that we next consider a FitzHugh-Nagumo-type model as studied in [23,24,25]: first in section 4.2 for τ = O(1), next for τ = ετ in section 4.3.2. In section 4.3.1, we consider the relation between the present paper and [44,45], in which the stability of planar interfaces on cylindrical domains is studied for models of the type (1.1), especially by analyzing the explicit prototypical example model considered in [44,45].

Invasion fronts in a Klausmeier-type dryland ecosystem model
On flat terrains, the ecosystem model proposed and studied in [1] coincides with the model of [16] reduced to one (instead of two) species, In this model, U (x, y, t) represents vegetation or biomass and V (x, y, t) water (or 'the limiting resource' [16]) -that diffuses much faster (i.e. the assumption 0 < ε 1 is natural in this model [39]). As in the Klausmeier model [31,41] and its reaction-diffusion modification introduced in [47], µ 1 > 1 represents the mortality rate of the vegetation and µ 3 > 0 the (yearly) precipitation rate. Through its zero U = 1/µ 2 (i.e. µ 2 > 0), the additional factor (1−µ 2 U ) in the U -equation represents the carrier capacity of the vegetation (in absence of mortality).
It is possible to write down an explicit expression for v * . Since this expression plays a crucial role in the companion paper [5] -see also section 5 -we also briefly formulate it here. By the linearity of the flow on M − 0 , we know that (v * , q * ) ∈ W u ((V − , 0)) is given by (v * , v * − µ 3 ). By the Hamiltonian associated to integrable system (4.5) on M + 0 we thus find, which can be simplified into, with v 0 = v * that exists for a specific choice c * of c. By the cubic nature of (4.2), both c * and u * (ξ) can be expressed explicitly in terms of v * . Since we will also encounter cubic fast reduced systems in several of the upcoming sections, we present a brief derivation of the formulae for c * and u * (ξ) in a general setting. Let u * (ξ) be a heteroclinic orbit of between the critical points (β − , 0) and (β + , 0) and assume that u * (ξ) also satisfies, (4.8) so that indeed lim ξ→±∞ u * (ξ) = β ± . By taking the derivative of (4.8) and substitution of (4.8) into (4.7), we obtain two polynomial expressions (in u) for u ξξ that coincide for, By (4.8) we also have an explicit heteroclinic solution of (4.7), 2 , which implies that the interface is stationary for v * = v * (µ 1 , µ 2 , µ 3 ) = 9 2 µ 1 µ 2 and that the desert invades the savanna for µ ∈ R 3 such that v * < 9 2 µ 1 µ 2 (and vice versa for v * > 9 2 µ 1 µ 2 ). The unstable interface does not finger but generates a dynamic interplay of bounded cusps similar to the Kuramoto-Sivashinsky dynamics associated to a sideband instability [13,26].
We conclude by the methods of geometric singular perturbation theory that a traveling planar invasion front (U h (ξ), V h (ξ)) between the bare soil state (Ū − ,V − ) and savanna state (Ū + ,V + ) indeed exists -see Fig. 5(c). Here, we do not go into the details of the spectral stability of (U h (ξ), V h (ξ)) with respect to perturbations that only depend on x (or ξ) -as usual we note that the stability can be established by the methods developed in [1,6,12,14,35,36,49] and the references therein. Thus, as in section 2.2 we assume that the front is stable with respect to perturbations that do not depend on y and that the critical eigenvalue curve λ c ( ) has λ c (0) = 0. By section 2.2 we know that the local parabolic character of λ c ( ) is given by (2.11) and (2.22) and that the (in)stability of the invasion front (U h (ξ), V h (ξ)) with respect to long wavelength transversal perturbations is determined by (1.6).
With F (U, V ; µ) and G(U, V ; µ) as in (4.1), we immediately find by the above analysis that, and u * ,ξ (ξ) > 0 by construction. Moreover, from which we conclude that the invasion fronts in dryland ecosystem model (4.1) cannot be stable with respect to transversal perturbations. Like in model (4.11), the invading desertification front is typically counter-invaded by fingering vegetation patterns. However, the transversally unstable invasion fronts do not necessarily finger: an invading vegetation front may generate dynamically evolving bounded cusps typical for Kuramoto-Sivashinsky-type dynamics near a sideband instability [13,26] -see Fig. 6.

Reversing desertification by front instabilities
In [19] a novel mechanism is proposed through which desertification may be reversed. It is shown that an invasion front of a bare soil state into a homogeneously vegetated that is stable with respect to longitudinal perturbations may be unstable with respect to transversal perturbations as planar interface. The vegetation fingers driven by this instability grow back into the bare soil state. Thus, instead of the collapse of the savanna state into a desert state that is predicted by the 1-dimensional point of view, (realistically) allowing for 2-dimensional perturbations along the interface predicts the counter-invasion of vegetation patterns into the desert, with as endstate a fully vegetated patterned domain (see [19] and Fig. 2).
The model considered in [19] reads, (cf. [50]). Note that model (4.1) can be seen as a simplification of this (somewhat) more realistic model. The evaporation rate is scaled to −1 in (4.1), while here it is unscaled (parameter µ 3 ) and takes into account the effect vegetation has on evaporation (µ 4 ). Moreover, following [31], the nonlinear mechanism by which plants increase water infiltration is modeled (and scaled) as U 2 V in (4.1), while in (4.11) the fact that plants extend spatially through their roots is included in the infiltration term, which yields the extended term U (1 + µ 1 U ) 2 V [33]. Apart from that the systems are only scaled in slightly different ways. As we shall see, the upcoming analysis is also quite similar -but slightly more technical and less explicit -to that of the simplified Klausmeier-type model (4.1). Moreover, the outcome is in essence identical.
The fast reduced system for traveling waves is now given by, and the associated two normally hyperbolic slow manifolds of interest are thus given by, (4.14) that intersects the {u = 0}-axis at v = 1 and has a minimum at, Note that in order to have heteroclinic front solutions jumping from M − 0 to M + 0 we need to impose that u m > 0 and v m < 1, i.e. that µ 1 > 1 2 -see Fig. 1(a). However, we need more: we also need that the critical points of the reduced slow flows associated on the homogeneous background states (Ū ± ,V ± ) are on M ± 0 , i.e. we need that 0 <V − < 1 andV + < v m . Moreover, since these background states are assumed to be stable, we also need that the associated critical points on M ± 0 are saddles -see section 2.1. Naturally, the homogeneous background states are determined by the intersections of either with u = 0 or with v = v F (u) (4.14). The former yields the bare soil state (Ū − ,V − ) = (µ 2 /µ 3 , 0), and we thus impose the condition µ 2 /µ 3 < 1. Since for u > 0 v G (u) must have a unique zero and lim u→∞ v G (u) = 0, we conclude that v G (u) has a maximum at u = u M > 0 ( Fig. 1(a)). Hence, the intersection {v = v F (u)} ∩ {v = v G (u)} has at most 2 elements that both correspond to homogeneously vegetated savanna states (see also Remark 4.1). Although the soil moisture level below a bare soil state can both be higher or lower than that below a vegetated state [30,38], model (4.11) has been derived from a more extended model under assumptions under which the water level in the bare soil state is necessarily higher than that in the vegetated state (cf. [33]). Therefore, we need to impose the conditionV + <V − , i.e. v G (Ū + ) < µ 2 /µ 3 or It now follows that only the largest element of {v = v F (u)} ∩ {v = v G (u)} can correspond to a critical point on M + 0 , and thus to (Ū + ,V + ): the other one is an unstable state. Note that this situation is similar to that of the (simplified) Klausmeier model (4.1).
As for (4.1), the slow flow on M − 0 is linear, v XX + 1 − v = 0, while the flow on M + 0 is given by where we have used (4.14) to simplify the expression. Since we assume that the background state (Ū + ,V + ) is stable as solution of (4.11), we know that the associated critical point (V + , 0) of (4.18) -its only critical point -must be a saddle (section 2.1). Thus it follows -by arguments completely similar to those in the previous section -that the (linear) unstable manifold W u ((V − , 0)) on M − 0 must intersect the stable manifold W s ((V + , 0)) on M + 0 (in their combined projection) -exactly as for model (4.1, see Fig. 5(b) -thereby defining the fast jumping point (v * , q * ) ∈ W u ((V − , 0)) ∩ W s ((V + , 0)) (with v m <V + < v * <V − = µ 2 /µ 3 and q * < 0).
To complete the skeleton structure, we only need to establish that there is a c * for which there is a heteroclinic connection u * (ξ) in (4.12) between the critical points (u − * , 0) = (0, 0) and (u + * , 0) = (f + (v * ), 0). Although the nonlinearity of (4.12) is quartic and we thus do not have explicit expression for c * and u * (ξ) as in the cubic case of (4.2), it indeed follows by phase plane arguments based on the integrable nature of (4.12) at c = 0 that c * and u * (ξ) exist and are uniquely determined. (In fact, an explicit condition on the value of v * at the Maxwell point, i.e. the parameter combination for which the front is stationary, can still be obtained (cf. section 3.2).) Therefore, we again conclude by the methods of geometric singular perturbation theory that a traveling planar invasion front (U h (ξ), V h (ξ)) between the desert state (Ū − ,V − ) and savanna state (Ū + ,V + ) exists. Moreover, we assume that the front is stable with respect to longitudinal perturbations and thus that the critical eigenvalue curve λ c ( ) has λ c (0) = 0 -as has become usual by now. By comparing (4.11) to (1.1) we see that, F v (u * (ξ), v * )u * ,ξ (ξ) = u * (ξ)(1 + µ 1 u * (ξ)) 2 (1 − u * (ξ)) u * ,ξ (ξ) > 0, (since 0 < u * (ξ) < 1 (cf. Fig. 1(a))), and that, SinceŪ + < u + * (see again Fig. 1(a)), it follows by (4.17), i.e. the fact that model (4.11) has been derived under assumptions that imply that the water level below the bare soil state is higher than below the vegetated state, that G(u + * , v * ) − G(u − * , v * ) < 0. Hence, under the overall assumption underlying our approach that ε is sufficiently small, we may conclude by (1.6) that all planar invasion fronts covered by model (4.11) are also unstable with respect transversal perturbations: an invading planar desertification front will typically be counter-invaded by a patterned, labyrinthine, vegetated state -see the simulations in [19] and Fig. 2.
Finally we note that it follows by the same arguments as above that a planar invasion front in the (a priori) un-ecological case thatV + >V − , i.e. a setting for which (4.11) cannot be derived, is not unstable with respect to long wavelength transversal perturbations, in the sense that λ 2,c < 0 (2.22), (1.6), since u + * <Ū + in that case. Thus, it could be stable (although one first needs to analyze the spectral stability of the interface for all ∈ R to arrive at that conclusion). Remark 4.1. In [19], the existence and stability of 'weak fronts' is studied by an amplitude equation approach near the co-dimension 2 point at which all 3 background states collide (i.e. asymptotically close to (µ 1,c , µ 2,c ) = ((µ 3 + µ 5 )/(2µ 3 ), µ 3 ) for which (Ū ,V ) = (1, 0) is a triple zero of F (U, V ) = G(U, V ) = 0 -if µ 4 = 0, which is the case considered in [19]). By deriving an expression for the curvature driven normal velocity of the interface, an explicit leading order upperbound ε c is obtained below which the front is unstable with respect to transversal perturbations -thus, diffusion coefficient ε is considered to be an O(1) parameter in this analysis. Since the assumption that ε is sufficiently small is underlying all our (asymptotic) analysis, the result of [19] agrees with our analysis. On the other hand, in [19] a simulation is shown of a 'linearly stable desertification front' that is 'unstable to finite-amplitude transverse modulations' (Fig. 4 in [19] in which µ = (3.5, 1.15, 3.2, 1.0, 0.5) ∈ R 5 ). This precise parameter combination is at the boundary of the applicability of our existence analysis, nevertheless, our stability analysis predicts that this interface must be spectrally (and thus linearly) unstable with respect to transversal long wavelength perturbations -if ε is sufficiently small. A careful simulation that is allowed to run for a sufficiently long time with ε = 1/ √ 150 ≈ 0.08 as in [19] shows that this value of ε indeed is sufficiently small: also this front is unstable, i.e. it develops fingers (like all other simulations shown in [19]).
Remark 4.2. The model considered in [27] can be seen as 'intermediate' between models (4.1) and (4.11) -like (4.11), it is also based on [50]. Application of our methods to this model goes along the very same lines as here and yields completely similar results.

A FitzHugh-Nagumo model
We consider a modified version of the FitzHugh-Nagumo models studied in [23,24,25], (for µ 1 ≥ 0). In this section we consider τ = O(1), in section 4.3.2 we explicitly study the impact of the pre-factor τ in (1.1) in the context of specific choices of the functions f ±,c (v). To mimic the systems studied in [23,24,25], we restrict ourselves to the symmetric case f − (−v) = −f + (v) and we note that (4.19) corresponds to the systems in [23,24,25] with Note that these flows are by construction symmetrical under v → −v, X → −X. It thus follows (by the symmetry) that if the slow reduced flow (4.20) on M + 0 has a saddle (V + , 0) withV + such that −f + (−V + ) < f c (V + ) < f + (V + ) and such that its stable manifold W s ((V + , 0)) intersects the v X -(= q)axis (on M + 0 ), then the unstable manifold W u ((V − , 0)) of the slow flow on M − 0 necessarily intersects the q-axis at the same point. Thus, by the general approach of section 2.1, there exists a traveling slow- where we note that necessarily µ 1 ≥ 0 (2.5)). Due to the cubic character (in U ) of F (U, V ) in (4.19), we have explicit expressions both for c * and u * (ξ) = u * (x − c * t), (4.9), (4.10). In the special (stationary) case that f c (0) = 0, that includes the models of [23,24,25] we obtain in a straightforward manner, (1.5), which yields by (2.22), Under the assumption that the front is longitudinally stable -which can be checked analytically -we conclude that the interface will be destabilized by transversal perturbations if (f + ) (0) < (f c ) (0) -so that for these f +,c (v) fingering may typically be expected (for τ = O(1)). In Fig. 7 we validate our asymptotic analysis by two numerical evaluations of the critical curve λ c ( ) for FitzHugh-Nagumo(-type) model (4.19): one with λ 2,c < 0 and the other with λ 2,c > 0. We conclude from this figure that our analysis indeed has the expected accuracy. We refer to [5] for similar numerical validations -especially for radially symmetric (multi-)front patterns -in the context of ecosystem model (4.1).
Lastly, we note that in the (general) case f c (0) = 0, taking f c (0) ∈ (−f + (0), f + (0)) as parameter yields by a slightly more involved calculus that λ 2,c (f c (0)) changes sign at a certain value of

On cylindrical domains
The stability of planar interfaces in a general class of singularly perturbed 2-component reaction-diffusion systems has also been considered in [44,45]. The class of models considered in [44,45] is more restricted than the models considered here. Several simplifying assumptions are formulated on the functions F (U, V ) and G(U, V ), it is for instance assumed that the slow flows on M ± 0 do not have additional critical points between the saddles at (V ± , 0) and the intermediate 'jumping line' {v = v * }. Especially in the setting of ecological models this indeed is a restriction [11,27]. Moreover, ecological models also typically have τ = O(1) -in fact τ = 1 -while in [44,45] only the case τ = ετ is considered. However, we have seen that τ = O(1) is more simple than τ = O(ε) and the nature and strength of the methods of [44,45] -that are based on [35,36] -strongly suggest that the more analytical approach of these works (compared to the present more geometrical set-up) can also be successfully applied beyond the class of models considered in [44,45]. In fact, we expect that any model that can be studied along the lines of this paper can also be analyzed by the methods of [35,36,44,45].
A much more important difference between our present study and [44,45] is that in the latter works the spatial domain is considered to be cylindrical: instead of (x, y) ∈ R 2 considered here, (x, y) ∈ R × Ω y with Ω y = (0, L) bounded (in fact, in [45], (x, y) ∈ Ω x × Ω y with also Ω x bounded (as in [35]), while Ω y ⊂ R N −1 in [44] -these differences can be overcome though). The (in)stability criteria obtained here are against long wavelength transversal perturbations: if an interface is unstable in this sense, it will only be unstable on a cylindrical domain R × (0, L) for L sufficiently large. So, in a way [44,45] consider the more complex issue of the destabilization of the interface for a given size L of the domain Ω y -in fact, a large part of the analytical focus of [44,45] is on the nature of the fastest growing unstable 'eigenmode', a topic we do not consider here.
A clear advantage of the more simple setting of the present work is that it enables us to also obtain simple instability results for the models considered in [44,45] under the assumption that L is sufficiently large. Instead of considering the models of [44,45] in their full generality, we focus on an example model also considered in [35,44,45], Clearly, the normally hyperbolic slow manifolds M ± 0 are given by, and the reduced slow flows on M ± 0 by, Both flows have (0, 0) as critical point and we want this to be a center in both cases: µ 2 > 0 and µ 1 + µ 2 > 0. We also want the slow flows to both have exactly one additional critical point on M ± 0 , i.e. µ 3 > max{0, 2(µ 1 + 2µ 2 )}, which then must be saddles: (V ± , 0) withV ± > 0. Since one of these saddles always must be inside the homiclinic loop attached to the other it follows that W u ((V − , 0)) must intersect W s ((V + , 0)) at (v * , q * ) with v * > min{V ± } > 0. Since v * < max{V ± } < 1 4 for µ 1 > 0 -the case considered [44,45] -the existence of the traveling front follows immediately; for µ 1 < 0, there may be 3 critical points on M − 0 in which case a fast heteroclinic jump may not always be possible so that additional (technical) conditions need to be added to ensure the existence of a traveling front. Since by construction u * (ξ) > 0, u * ,ξ (ξ) > 0 and, it follows immediately by (1.6) that the planar traveling interface must be unstable in the case µ 1 > 0 considered in [44,45]: even if it is stable with respect to longitudinal perturbations, it is unstable against transversal long wavelength perturbations.
However, we so far only considered τ = O(1), while τ = ετ = O(ε) in [44,45]. Nevertheless, this does impact neither the instability result itself nor its simplicity (for µ 1 > 0). As in [45] and section 3.2, we focus on the case of a stationary planar interface, i.e. we assume that µ ∈ R 3 (in (4.24)) is such that H f (ũ − * , 0) =H f (ũ + * , 0) = 0 (cf. (3.12) and [35,45]). Since the fast field associated to (4.24) is cubic, we noteṽ * = (3.9), (4.8), (4.9), (4.25)), which implies that the planar interface is unstable with respect to transversal perturbations for,τ at which it undergoes a bifurcation into traveling waves (by section 3.2). Thus, also for τ = O(ε), the straightforward observations (4.25) immediately yield that a longitudinally stable stationary front on R 1 must be unstable as planar interface on R 2 up to the critical valueτ * at which it (also) destabilizes as 1dimensional front. In the setting of [44,45], this implies that the interface is always unstable on cylindrical domains R × (0, L) for L sufficiently large. For a given L, the instability condition is much more involved and needs to be evaluated numerically, see [45] (where we note that in [45] also the x-domain is bounded, which may impact the nature of the front).

The (in)stability of bifurcating traveling planar interfaces
Finally, we consider the case τ = ετ in the setting of FitzHugh-Nagumo model (4.19). Naturally, bifurcations may take place due to the geometry of the now asymmetrical slow reduced flows (forc = 0, cf. Fig. 4(b)). However, here we focus on the bifurcation into traveling waves as discussed in section 3.2. In this section we consider the question of the (in)stability of the traveling interfaces that bifurcate off the stationary front with respect to transversal perturbations. Note that it follows from the analysis in section 3.2 thatλ 2,c associated to the primary stationary interface has changed sign in the bifurcation (3.11), however, the existence/Melnikov analysis of section 3.2 does not provide insight in the transversal stability of the bifurcating waves, and more specifically not in the associated sign ofλ 2,c .
As in 3.2, this issue can be studied in the general setting of equation (1.1). Especially since this is a somewhat technical affair, we for simplicity restrict ourselves to a simple case for which explicit analysis is possible: we consider the FitzHugn-Nagumo model (4.19) in the most simple case in which the slow reduced flows on M ± 0 are linear (cf. the approximations in [23]). Thus, we set f + (v) = 1 + µ 2 v and thus (4.20)). Naturally, we need to assume that µ 1 − µ 2 > 0 (so that the critical point of (4.26) is a saddle). It follows that the slow components of the skeleton structure that determines the front are given by, 27) and thus that,ṽ * =ṽ ± (0) =c (µ 1 − µ 2 ) c 2 + 4(µ 1 − µ 2 ) (4.28) (at leading order in ε). By the cubic nature of F (U, V ), the relation C(v 0 ) that determines for which value of 'friction term'cτ there is a heteroclinic orbit in (3.1) can again be determined explicitly, 9). By (4.28), (3.3) thus takes the explicit form, The further simplifying assumption f c (v) = µ 3 v yields, with the additional condition µ 3 −µ 2 > 0 (sinceτ > 0). Thus, the situation is as in section 3.2: independent of the value ofτ there is a stationary front that undergoes a bifurcation into (counter-propagating) travelling waves/planar interfaces asτ decreases through In fact, we have already obtained more than in section 3.2: we know that the bifurcation is supercritical (i.e. the bifurcating waves appear for values ofτ at which the stationary wave is unstable). Note however that there is also a difference with the general setting of section 3.2: due to the assumed symmetry f − (−v) = −f + (v) the stationary interface exists in (4.19) without an additional condition on the parameters µ (∈ R 3 in (4.19) with the present choices of f + (v) and f c (v)). Note also that we can copy (4.22) for the stationary fronts and that by (4.27),Ĩ Thus, assuming that the stationary front remains stable against longitudinal perturbations -which can be established analytically -we conclude that the associated planar interface is unstable against (long wavelength) transversal perturbations forτ >τ * , whereτ * of (4.30) indeed coincides with a zero in the denominator of (4.32) -as was already established by the general analysis of section 3.2.
Under the assumption that the bifurcating travelling interfaces are stable as solutions of x ∈ R, which is natural by the above obtained supercritical nature of the bifurcation (and can be established analytically), we can now consider the question of the stability of the interfaces with respect to long wavelength transversal perturbations. Since we know from section 3.1 that this is determined by (3.11), we can do this by evaluatingF * ,G * (3.8),Ĩ f andĨ s (3.9) forc = 0.
As in section 3.2, we consider the system near the bifurcation: we introduceτ > 0 and second small parameter δ, i.e. 0 < ε δ 1 and independent of ε -by settingτ =τ * − δ 2τ , so that by (4.28) and (4.29), , (4.33) withc,ĉ andv 1 as defined in section 3.2. It follows by the derivation in appendix A that M * (τ ) = F * Ĩs +τG * Ĩf = 8 3 for δ sufficiently small. By (3.11) and (A.6) we now have for the bifurcating traveling waves, (at leading order in ε) and conclude that (in this specific case) the counter-propagating planar interfaces that appear from the bifurcation into traveling waves are unstable with respect to long wavelength transversal perturbations -like the longitudinally stable stationary interface they originate from (4.32).
Note that although the stationary front has become longitudinally unstable by the bifurcation into traveling waves, one can still apply (4.32) to conclude that the curve λ t ( ) through the translational eigenvalue -which thus no longer is the critical curve λ c ( ) -has a local maximum at λ t = 0, i.e. that the orientation of λ t ( ) flips asτ decreases throughτ * . However, it follows from (4.35) -and the analysis in the appendix -that the orientation of the critical curve λ c ( ) does not change by the transition from the stationary interface to the bifurcating traveling interfaces. Especially the fact that the sign of (4.34) does not depend on the parameters in the problem suggests that there may be an underlying mechanism, i.e that there is a direct relation between the sign of λ 2,c for the (longitudinally stable) standing interface before the bifurcation and that of the bifurcating traveling interfaces. Such a general relation may be uncovered by a similar analysis as in appendix A in the general setting of section 3.2. Here, we refrain from going further into this -see also the upcoming discussion.

Discussion
In [2], predictions of model studies of vegetation patterns in dryland ecosystems have been validated by observational evidence. These model studies go back to [42,43] in which the persistence of Turing patterns in a (generalized Klausmeier) model of type (1.1) under slowly varying parameters was investigated. The main conclusion of [2] was that the multistability induced by the richness of spatial (Turing) patterns increases the resilience of an ecosystem: instead of a catastrophic collapse, a patterned state may adapt under worsening external circumstances to a nearby adjacent (vegetated and patterned) state. It is postulated in [37] that this is not special for dryland ecosystems: any ecosystem that exhibits spatial patterns is expected to have an increased resilience. In fact, it is argued in [37] that spatial patterning may provide the ecosystem with a mechanism to evade tipping, and thus to circumvent catastrophic collapse.
Although the patterns of [2,42,43] that initiated the insights of [37] are all of Turing pattern type, it is observed in [37] that there are other types of patterns -with backgrounds that are not of Turing typethat will have the same positive impact on the resilience of an ecosystem. An important role is given to the patterns that originate from the interfaces between homogeneous states. These coexistence states occur naturally in various kinds of ecosystems, moreover, these ecosystems are typically modeled by systems of reaction-diffusion equations like (1.1) -see [3,16,17,19,48,50] and the references therein. Thus, the conditions obtained here on the instability of planar interfaces and the associated onset of fingering patterns -and especially the insight of section 4.1 on the natural occurrence of this instability in ecosystem models -have, by [37], a direct interpretation in the context of the resilience of ecosystems.
Nevertheless, the present results mostly only are a first step towards understanding the relevance of transversal instabilities in ecosystems.
• Non-planar fronts. The present focus on straight, flat interfaces orthogonal to the longitudinal direction of a 1-dimensional front is quite restricted -especially from the ecological point of view. In the companion paper [5], radially symmetric curved fronts are considered -in the explicit context of ecosystem model (4.1). Using geometric blow-up methods, first a spot pattern -i.e. a localized (circular) vegetated area surrounded by bare soil -is constructed rigorously. Gap patterns can be constructed in a very similar fashion and we find that parameter condition (4.6) on the existence of stationary planar fronts determines the transition between the existence of spots and gaps (in parameter space). It is additionally indicated in [5] how the same geometric approach yields the existence of ring and target patterns. Moreover, the stability of the gap and spot patterns is considered. It is shown that spots and gaps with a sufficiently large radius are unstable with respect to the same (fingering) instability as encountered here: simulations show that the evolving fingers merge the vegetated and bare soil patches into a fully mixed labyrinthine state.
• Multiple fronts and non-flat terrains. Existence conditions (E-I) -(E-IV) allow for the construction of double-front homoclinic stripe patterns -either of gap or of spot type. Moreover, these localized stripes are the 'endpoints' of a 1-parameter family of periodic stripe patterns -see [27] for the construction of these patterns in a model of intermediate complexity between (4.1) and (4.11) (Remark 4.2). Stripe patterns are typically observed on sloped terrains [9,21], therefore it is natural to study these patterns in an extended version of model (1.1) that includes an advection term '+µ m+1 V x ' in the V -equation -which models the downslope flow of water in an ecosystem setting [1,31,40,41]. Central questions to consider are: will the multiple-front patterns (also) be unstable with respect to transversal long wavelength perturbations, and what is the impact of parameter µ m+1 = µ m+1 (ε)? (This is the subject of work in progress.) • Multi-component systems. Ecosystem models typically contain more than one type of vegetation and often also separate equations for soil water and overland water -see for instance [3,16,30,33,37] and the references therein. The 2-component ecosystem models of [1,2,17,19,27,31,48,50] are either obtained by reduction from a more extended model or postulated as conceptual, simplified, model. To achieve a next level in ecological validity, the present analysis needs to be extended to multi-component systems. The observation that many ecosystem models have various ordered small parameters opens up the possibility for setting up a geometric multi-scale analysis -as initial explorations indicate.
The present analysis also triggers a number of mathematical research questions, of which we mention a few.
• Bifurcations of fronts and interfaces. In section 2.2, stability assumption (S-I) was formulated more precisely: apart from assuming that the interface is stable against longitudinal perturbations, it was also imposed that the dimension of the kernel of operator L (2.9) is equal to 1. In other words, it is assumed that the parameters (τ, µ) of (1.1) are chosen away from potential (local) bifurcation sets. This assumption underlies the approach of this section -and the paper -since the leading order approximations (2.22) and (3.11) -that underlie criteria (1.6) and (1.7) -can now be derived by approximating the 1-dimensional kernel of the adjoint operator L A (2.12). As (τ, µ) or (τ , µ) passes through a bifurcational (co-dimension 1) set (τ * , µ * )/(τ * , µ * ) -for instance as in section 3.2 -then there is a second eigenvalue curve λ b ( ) that passes through λ c ( ). Naturally, λ b (0) = λ c (0) = 0: the dimension of the kernel of L is (at least) 2 at the bifurcation set. Thus, one has to extend the approach that led to (2.22) and (3.11) to determine the local character -as function of | | 1 -of λ b ( ) and λ c ( ) as the bifurcational set is crossed. This is an interesting line of research, especially since the results of section 4.3.2 suggest that there may be a relation between the signs of λ 2,c of the longitudinally stable interfaces before and after the bifurcation.
• Breakdown of existence conditions. Although we write at the end of section 2.1 'the conditions on F (U, V ) and G(U, V ) are quite mild and the slow-fast-slow interfaces (U h (ξ), U h (ξ)) can be constructed for many models', it is certainly not always possible to construct a traveling front solution that corresponds to a heteroclinic orbit of (2.1) if only condition (E-I) holds -unlike in the scalar case (1.2). Conditions (E-III) and (E-IV) may be violated for generic classes within the full family of systems (1.1) (see Remark 2.2 and [11] for (E-II)). It is natural to expect, under the fundamental assumption (E-I), that (1.1) 'must' exhibit some kind of localized patterns that approach (Ū ± ,V ± ) for x → ±∞. We are not aware of any systematic studies of this issue. We note that such a study could start with considering the implications of varying µ such that the intersection W u ((V − , 0)) ∩ W s ((V + , 0)) vanishes (i.e. the parameters are varied beyond the validity (E-IV)). The example of [12] shows that this may lead to (spatially localized) finitetime blow-up behavior.