Large Amplitude Radially Symmetric Spots and Gaps in a Dryland Ecosystem Model

We construct far-from-onset radially symmetric spot and gap solutions in a two-component dryland ecosystem model of vegetation pattern formation on flat terrain, using spatial dynamics and geometric singular perturbation theory. We draw connections between the geometry of the spot and gap solutions with that of traveling and stationary front solutions in the same model. In particular, we demonstrate the instability of spots of large radius by deriving an asymptotic relationship between a critical eigenvalue associated with the spot and a coefficient which encodes the sideband instability of a nearby stationary front. Furthermore, we demonstrate that spots are unstable to a range of perturbations of intermediate wavelength in the angular direction, provided the spot radius is not too small. Our results are accompanied by numerical simulations and spectral computations.


Introduction
The phenomenon of vegetation pattern formation in dryland ecosystems has attracted attention in the last several decades as a mechanism of ecological resilience.While frequently considered to be an early warning sign for desertification [16,17,26,28,32,33,34], the formation of localized vegetation patches or patterns has also been viewed as a means of evading such critical transitions [30].The interaction of infiltration feedback mechanisms and competition for water resources results in the formation of vegetation patches [25,31,37,43,45].On sloped terrain, one observes vegetation stripes, or bands, aligned perpendicular to the slope, while on flat ground spots, gaps, or disorganized labyrinth patterns are prevalent [1,6,7,13,24,40].Spots, gaps, rings or other radially symmetric patterns (sometimes called "fairy circles") have been observed extensively in drylands in Australia and Africa [14,27,29] (and other ecosystems, such as submarine seascapes [35]), and have served as the focus of many studies of self-organization in ecosystems.
Vegetation pattern formation is frequently modeled by multi-component reaction diffusion systems.In such models, there is a well-developed theory of spot and stripe pattern formation near the onset of Turing instabilities [17,15,39].However, far less is known analytically concerning large amplitude or far-from-onset planar vegetation patterns.A number of studies have considered existence and stability properties of banded vegetation [2,3,10,38], as well as desertification fronts [4,12], but less is known concerning far-from-onset radially symmetric vegetation patches.Prior work has considered small amplitude radial solutions [18] and 1D simplifications [19].However, to our knowledge, no rigorous studies exist concerning large amplitude radially symmetric vegetation patches in a dryland ecosystem model.
We consider the model introduced in [2] which is a modification of the dryland ecosystem model originally proposed by Klausmeier [20].This model also coincides with that studied in [11], in the case of a single species.Here U represents water availability, V represents vegetation density, and the parameters a, m, b are positive and represent rainfall, mortality, and inverse of soil carrying capacity, respectively.The diffusion coefficient δ 2 represents the ratio of timescales of diffusion of water vs. vegetation.Here, we assume that water diffuses much faster than vegetation, so that 0 < δ 1 is a small parameter [33].We note that Klausmeier's model originally had b = 0.While a singular perturbation analysis of this case is possible (see, e.g.[3]), this limit is highly singular, requiring several rescalings and blow up techniques to account for passage near a non-hyperbolic slow manifold.Here, we therefore focus on the case b > 0.
We are interested in the formation of radially symmetric vegetation spots (localized vegetation patches surrounded by bare soil) and gaps (localized regions of bare soil in an otherwise uniformly vegetated state) in (1.1).Exploiting the small parameter δ 1, we will use geometric singular perturbation methods to construct radially symmetric solutions through a spatial dynamics approach in the radial coordinate.Our approach is similar to that in [41], in which the authors construct far-from-onset spot solutions in a 3-component FitzHugh-Nagumo model.However, the nonlinearities in (1.1) introduce complications in the analysis due to the fact that the reduced flow on the resulting slow manifolds is no longer linear as in the case in [41].However, we will show that radial solutions can still be constructed in (1.1).
Remark 1.1.In the setting of classical models like Gray-Scott, Gierer-Meinhardt and Schnakenberg, the existence, stability and interactions of radially symmetric localized 'spikes' has been studied -see for instance [5,23,44] and the references therein.However, these patterns differ from the ones considered here since they have a homoclinic nature: unlike the present gaps and spots, the regions in which these spikes are not close to the background states are asymptotically small.The patterns considered here can be seen as two-dimensional (radially symmetric) versions of the one-dimensional 'mesa patterns' studied in [22] and in [19] in the setting of vegetation patterns.Like in the present model, and unlike in [41,42], the 'slow flows' (for the spatial dynamics) considered in these papers are nonlinear -which is typically the case for ecosystem models [9].
We search for stationary solutions, which are radially symmetric, so that (U, V )(x, y, t) = (u, v)(r), r ∈ [0, ∞), and ∆ = ∂ 2 r + r −1 ∂ r under planar radial symmetry, where r = (x 2 + y 2 ) 1/2 .Such solutions satisfy the ordinary differential equation Figure 1: Shown is a schematic of a radial profile for a vegetation spot solution (left) and gap solution (right) of (1.3) as in Theorem 1.2 and 1.3, respectively.The profile contains a sharp transition from the vegetated state to the bare soil state in an interval of width O(δ) of the critical radius r I = O (1).
which can be rewritten as a first order non-autonomous system This system admits up to three equilibria given by the steady states above: the desert state P 0 = (a, 0, 0, 0), and if a m > 2(b + 1 + b 2 ), there are two additional equilibria P 1,2 := (U 1,2 , 0, V 1,2 , 0) corresponding to uniform vegetation.
Spots, gaps, and other localized radially symmetric solutions are constructed as orbits of (1.4) which are asymptotic as r → ∞ to one of these steady states, and which are bounded as r → 0. In order to find solutions of this system, we construct candidate solution orbits in different subsets of the spatial domain r ∈ [0, ∞) and glue these together to build a solution on the entire domain.In particular, we construct the solution over three primary regions: the core r ∈ [0, r c ], where r c = O(δ), the far field r ∈ [r I , ∞), where r I = O(1), and the transition region(s) in between; see Figure 1.
Our main result concerning the existence of spots is the following.
or equivalently, where U 2 is defined as in (1.2).Then there exists V c (a, b, m) > 0 and r I (a, b, m) > 0 such that for all sufficiently small δ > 0, (1.1) admits a stationary, bounded, radially symmetric vegetation spot solution (U, V ) = Similarly, we have the following theorem concerning the existence of gaps.
Theorem 1.3.(Existence of gaps) Fix a, m, b > 0 satisfying (1.5).Suppose or equivalently, where U 2 is defined as in (1.2).Then there exists r I (a, b, m) > 0 and ν > 0 such that for all sufficiently small δ > 0, (1.1) admits a stationary, bounded, radially symmetric vegetation gap solution (U, V ) = (u g , v g )(r; a, b, m, δ) satisfying where V 2 is defined as in (1.2), with a single interface between the vegetated and desert states occurring at the radius r = r I (a, b, m) + O(δ).
See Figure 2 for a visualization of the condition (1.5) necessary for the existence of spots/gaps in Theorems 1.2-1.3.
Remark 1.4.The methods used in this paper to construct radially symmetric spot and gap solutions can be applied in a similar manner for the construction of solutions such as rings, targets, or other radially symmetric profiles, with different, perhaps more complex, conditions on parameters which ensure their existence.We provide some numerical evidence for the existence of such solutions in §5 and describe how one might go about constructing these orbits, but do not provide the lengthy technical details here.
The spot and gap solutions of Theorems 1.2-1.3 will be constructed as heteroclinic orbits for r ∈ [0, ∞) using geometric singular perturbation theory.The main idea, following [41], is to find the orbits as intersections of a core manifold of solutions which remain bounded at r → 0, and a far-field manifold of solutions which decay to one of the states P 0 or P 2 as r → ∞.These manifolds are each three-dimensional, and, unlike in [41], it is not possible to obtain an explicit description of these manifolds due to the fact that the flow on (one of) the slow manifolds of (1.4) is nonlinear.This introduces complications which are handled through the use of an intermediate scaling and careful qualitative analysis of the nonlinear non-autonomous reduced flow on this slow manifold.The non-autonomous nature of the equation (1.4) makes this a somewhat challenging construction.To demonstrate how these orbits are constructed using the slow/fast geometry of (1.4), it is helpful to first consider the simpler construction of traveling or stationary planar front solutions of (1.1), which manifest as heteroclinic orbits in an appropriate traveling equation with a similar geometry to that of (1.4).
It is shown in companion paper [4] that the invasion fronts in (1.1) -that can be seen as spots or gaps with radius r I → ∞ -are unstable with respect to a sideband/finger instability.Therefore, we next consider the stability of spots.Without going into the details of a fully rigorous analysis, we first consider the spectral stability problem for large spots.By a careful asymptotic analysis we recover the sideband instability mechanism and conclude that large spots are unstable and will typically form finger-like patterns -like the planar invasion fronts.By similarity, we conclude that the same is true for large gaps.These observations are confirmed by direct simulations: in Fig. 9 in §5 we show snapshots of a spot and a gap that both evolve towards labyrinthine patterns after undergoing such an instability.Moreover, we subsequently conclude by considering perturbations within a range of intermediate wave numbers that spots (and gaps) with O(1) radius r I must also be unstable: only sufficiently small spots and gaps may possibly be stable -see again §5 and especially Fig. 10 for a brief numerical study.These small gaps correspond to fairy circles [14] and also appear as stable vegetation patterns in numerical simulations of the dryland ecosystem model of [46] -a model that can be seen as a slightly more extended version of (1.1) [4] (and that has been deduced from more extended models to study fairy circles).
The set-up of the paper is as follows.In §2, we consider the construction of traveling fronts in (1.1), while in §3, we treat the radially symmetric case and provide the proofs of Theorems 1.2-1.3.The spectral stability of the radial spot and gap solutions is considered in §4 using formal asymptotic methods, and we include numerical simulations and a brief discussion in §5.

Stationary and traveling planar fronts
To motivate the approach for the existence analysis of radially symmetric solutions of (1.1), we first consider the (simpler) case of constructing stationary and traveling front solutions.We pose the traveling wave ansatz (u, v)(x, y, t) = (u, v)(ξ), where ξ = x − δct is a traveling wave coordinate; here c = 0 corresponds to stationary front solutions, while c = 0 corresponds to solutions which propagate with wave speed δc.
Using a geometric singular perturbation approach, we can construct bistable front solutions as perturbations from slow/fast heteroclinic orbits between the desert state (U, V ) = (a, 0) and the vegetated state (U, V ) = (U 2 , V 2 ) in the traveling wave ODE (2.1) We note that in the case c = 0, the system (2.1) corresponds to (1.3) in the far field limit r → ∞.Thus, it is natural to first consider the geometry of (2.1), as solutions with radial symmetry are constructed by matching a solution which is bounded near the core r = 0, with a stationary solution which (approximately) satisfies (2.1) in the far field limit r → ∞.

Slow/fast analysis
We can write (2.1) as a first order system which we refer to as the slow system, and upon rescaling ξ = δζ, we obtain the equivalent fast system Setting δ = 0 in (2.3) yields the layer problem (2.4) This system admits up to three equilibria depending on u, given by M 0 (u) := (0, 0), and if u > 4bm, M ± (u) := (v ± (u), 0), where When u = 4bm, the two equilibria M + (4bm) = M − (4bm) coincide.Computing the linearization a short computation shows that for m, b > 0 the fixed points M 0 and M + (u), u > 4bm are always of saddle type, while M − (u) is focus or node for c = 0, and a center for c = 0.The equilibrium M + (4bm) = M − (4bm) is nonhyperbolic.Taken together, the set of equilibria of the layer equation (2.4) corresponds to the critical manifold obtained by setting δ = 0 in (2.2).We therefore decompose M 0 into three branches , where and the latter two manifolds M ± 0 meet along the nonhyperbolic fold curve F := {q = 0, v = 1/2b, u = 4bm}.The manifolds M 0 0 and M + 0 are normally hyperbolic and of saddle-type.The reduced flow on each branch of the critical manifold is obtained by setting δ = 0 in (2.2), and is given by where v * (u) = 0, v ± (u) on the branches M 0 0 and M ± 0 , respectively.
Figure 3: A singular heteroclinic orbit between P 2 and P 0 is formed by concatenating slow orbits on the critical manifolds M 0 0 and M + 0 with a fast orbit φ vd of the layer problem (2.4).

Singular orbits
The desert equilibrium state P 0 = (a, 0, 0, 0) lies on the branch M 0 0 , while the vegetated state The state P 2 = (U 2 , 0, V 2 , 0) can lie on either M − 0 or M + 0 ; the latter occurs provided a > 4mb + m/b.In this case, a linear stability analysis shows that both (U 0 , V 0 ) and (U 2 , V 2 ) are temporally stable as homogeneous equilibria of the PDE (1.1) [4], and we can proceed to construct singular bistable fronts connecting P 0 and P 2 (or vice versa); see Figure 3.
To do this we form a singular slow/fast/slow heteroclinic orbit between P 0 , P 2 , by concatenating a slow orbit on M 0 0 with another slow orbit on M + 0 via a fast heteroclinic orbit in the layer problem.Within the layer problem (2.4), for any value of (u, p) satisfying u > 4mb, by adjusting the speed c appropriately it is possible to construct fast heteroclinic orbits φ dv (ζ) = (v dv , q dv )(ζ) from M 0 0 to M + 0 and φ vd (ζ) = (v vd , q vd )(ζ) from M + 0 to M 0 0 , corresponding to fast desert-to-vegetation and vegetation-to-desert fronts, respectively.These orbits can be computed explicitly as follows.Considering (2.4), we search for solutions satisfying the ansatz q = Cv(v − v + (u)) for some value of C = 0. Substituting into (2.4),we obtain the algebraic equation which can be solved to find C = ± bu/2 and wave speeds We then obtain the explicit profiles (up to translation) (2.12) Figure 4: The stationary fronts φ vd and φ dv of (2.4).
Examining the reduced flow on M 0 0 we see that within M 0 0 , the equilibrium (a, 0), corresponding to P 0 , is a saddle-type equilibrium with (un)stable manifolds W u/s (a, 0) given by the lines p = ±(u − a).On M + 0 , the equilibrium (U 2 , 0), corresponding to P 2 , is also saddle-type equilibrium, with (un)stable manifolds W u/s (U 2 , 0) given by the level set E(u, p) = 0 of the conserved quantity noting that E(U 2 , 0) = 0.
To construct a singular heteroclinic orbit from say P 2 to P 0 , we follow the slow unstable manifold W u (U 2 , 0) of P 2 , then a fast jump φ vd , then the slow stable manifold W s (a, 0), given by the line p = −(u − a).Since the slow variables (u, p) are constant across the fast jump, this is only possible if there is an intersection of W u (U 2 , 0) and W s (a, 0) when W u (U 2 , 0) is projected onto M 0 0 , which occurs if there exists u = u * such that E(u * , a − u * ) = 0, or equivalently For an open region in (a, b, m) parameter space, there is a critical value U 2 < u * < a, depending on a, b, m, which satisfies this criterion, and therefore the speed c * is determined so that the fast layer jump φ vd exists for (u, p) = (u * , a−u * ).We denote by (u +,∞ , p +,∞ )(ξ) the solution of the reduced flow on M + 0 corresponding to the unstable manifold W u (U 2 , 0), which satisfies (u +,∞ , p +,∞ )(0) = (u * , a−u * ), and we denote by (u 0,∞ , p 0,∞ )(ξ) the solution of the reduced flow on M 0 0 corresponding to the stable manifold W s (a, 0), which satisfies (u 0,∞ , p 0,∞ )(0) = (u * , a − u * ).
Thus the singular front solution follows the solution (u +,∞ , p +,∞ )(ξ) of the reduced flow on M + 0 , followed by the fast front φ vd (ζ), at finally the solution (u 0,∞ , p 0,∞ )(ξ) of the reduced flow on M 0 0 .For 0 < δ 1, these singular fronts can be shown to perturb to front solutions of (2.1) using standard methods of geometric singular perturbation theory, using the wave speed c as a free bifurcation parameter.The construction of fronts from P 0 to P 2 follows similarly.

Stationary fronts
Stationary fronts correspond to the case c = 0. We describe the geometry of the singular orbit(s) in this case, as it will be useful in the forthcoming construction of spot and gap solutions.Proceeding as in §2.2, we find Figure 5: Shown are singular orbits representing stationary fronts of (2.1) obtained by concatenating slow orbits of (2.9) on the critical manifolds M 0 0 and M + 0 with fast orbits φ vd , φ dv of the layer problem (2.4) for c = 0, u = u f .that (2.4) admits a pair of heteroclinic orbits where (2.17) where again φ dv (ζ) represents the desert-to-vegetation state transition front which jumps from M 0 0 and M + 0 , while φ vd (ζ) represents the vegetation-to-desert state transition front which jumps from M + 0 and M 0 0 ; see Figures 4-5.Note that we assume a > u f , so that the equilibrium P 0 lies 'above' the critical fronts φ dv , φ vd .A lengthy but straightforward computation shows that , so that the equilibrium P 2 lies 'below' the fronts φ dv , φ vd .
While this pair of heteroclinic orbits exists for any a, b, m > 0, in order to construct a singular slow/fast heteroclinic orbit, we still require an intersection of W u (U 2 , 0) and W s (a, 0) in the reduced flow when W u (U 2 , 0) is projected onto M 0 0 .Since the jump height u = u f is fixed by the condition c = 0, this only occurs if which gives an implicit condition on the parameters a, b, m for which a singular stationary front exists.We will show in §3 that this condition describes the boundary (in parameter space) which separates the existence region of radial spots versus gaps, which will be constructed as slow/fast fronts in the non-autonomous system (1.4).The radius of the spot/gap will serve as a free parameter which allows for the construction of a solution for parameters which satisfy (1.6) in the case of spots, or (1.8) in the case of gaps.

Sideband (in)stability of planar fronts
Given a front solution constructed as in §2.2, we briefly examine the stability of the front as a planar interface, which will help motivate our formal stability results in the case of radial spot and gap solutions in §4.We assume that a heteroclinic vegetation-to-desert front solution (u h , v h )(ξ; δ) exists which connects the state P 2 to P 0 (the case of a desert-to-vegetation front is similar), and has been constructed as a perturbation from one of the singular slow/fast fronts in §2.2, with speed c = c h (δ) = c vd (u * ) + O(δ).
We linearize (1.1) about this front solution in a comoving frame using an ansatz (U, V ) = (u h , v h )(ξ; δ) + e λt+i y (u, v)(ξ) for ∈ R, which results in the eigenvalue problem Due to translation invariance, this eigenvalue problem has a solution when λ = = 0, with eigenfunction given by the derivative (u h , v h )(ξ; δ).For the purposes of this discussion, we focus only on this critical, marginal eigenvalue, and we assume that all other spectrum of the front for = 0 (that is, the spectrum corresponding to 1D longitudinal perturbations in the direction of propagation) is bounded away from the imaginary axis in the left half plane.
Focusing on this critical eigenvalue, we consider its continuation for small | |: as the eigenvalue problem only depends on through terms of O( 2), we anticipate that this critical translation eigenvalue expands as This eigenvalue describes the stability of the front to long wavelength perturbations transverse to the front, and the stability is thus determined by the sign of the coefficient λ c,2 .In a companion paper [4], we have developed a procedure to compute this coefficient in a general class of two-component singularly perturbed reaction diffusion systems, which includes (1.1) as an example, and we briefly describe the results here.
We rewrite the stability problem (2.19) in the form where We now expand Substituting into (2.21), at leading order we have the eigenvalue problem

.23)
This leads to the Fredholm solvability condition where denotes the bounded solution to the adjoint equation where From this we obtain an expression for λ 2,c In [4], it is shown that to leading order λ 2,c is given by > 0, so that the traveling planar fronts are always unstable to long wavelength transverse perturbations.This has implications for the (in)stability of spot/gap solutions, as we will show in §4 using formal asymptotic methods that the spectrum for radial spots/gaps of sufficiently large radius is approximated by that of a nearby stationary front solution.

Existence of radially symmetric spots and gaps
With the construction of the front solutions in §2 in mind, we now focus on the construction of a vegetation spot solution, consisting of a single vegetation patch localized near r = 0, with a single interface at some radius r = r I (to be specified), at which the profile transitions from the vegetated state in the core to the desert state in the far-field.The case of gaps is similar, and we will briefly outline the differences in §3. 4.
Throughout the analysis we treat 0 < δ 1 as a singular perturbation parameter.At times, it will also be convenient to consider (1.4), which we refer to as the 'slow' system, with respect to the rescaled radial coordinate s = r/δ, which results in the system We refer to (3.1) as the 'fast' system.
The spot solution will be constructed as a perturbation from the singular limit structure associated with (3.1), and consists of three pieces: The core vegetated and far-field desert states are given as slow orbits which lie near equilibria on saddle-type slow manifolds M 0 δ , M + δ within (3.1) (to be defined below), while the interface between these states is given by a fast layer orbit between these slow manifolds which is inserted at a particular radius r = r I .The construction has very similar geometry as in the construction of traveling fronts in §2, though with some complications due to the nonautonmous nature of the equation and the singularity at r = 0. Additionally, since the spot interface is stationary, the speed is not available as a free parameter; however, the jump value r = r I can be thought of as a free parameter which is chosen in such a way in order to ensure that the stable manifold of the far-field desert equilibrium and the unstable manifold associated with the core vegetated states precisely intersect transversely across the fast jump.

Far field region and fast transition
In the far-field, we consider r ∈ [r, ∞) for arbitrary r > 0 fixed independently of δ > 0.Here we define a new variable k := 1/r and similarly consider the region k ∈ [0, k] where k := 1/r for the corresponding (autonomous) system which is a slow-fast system with two fast variables and three slow variables.

Slow manifolds away from the core
Setting δ = 0, we see that (3.2) admits a three dimensional critical manifold defined by 2) on the fast scale s = r/δ, we have the equivalent system where k = (δs) −1 .Setting δ = 0 in this system yields the layer problem in which the slow variables (u, p, k) act as parameters.As in this case of the layer problem (2.4) associated with the traveling fronts in §2, this system admits up to three equilibria depending on u, given by M 0 (u) := (0, 0), and if u > 4bm, M ± (u) := (v ± (u), 0), where When u = 4bm, the two equilibria M + (4bm) = M − (4bm) coincide.Computing the linearization we find that this corresponds to (2.6) in the stationary case c = 0, so that for m, b > 0, the fixed points M 0 and M + (u), u > 4bm are always of saddle type, while M − (u) is a center for u > 4bm.The equilibrium M + (4bm) = M − (4bm) is nonhyperbolic with a double-zero eigenvalue.Taken together, the set of equilibria of the layer equation (3.5) corresponds to the critical manifold The manifolds M 0 0 and M + 0 are normally hyperbolic, while M − 0 is not.We note that compared with the analysis in §2, these manifolds are actually subsets of a 5-dimensional ambient space due to the additional slow variable k, but we slightly abuse notation and continue to refer to these as M * 0 , for * = 0, ±, as they are defined by the same algebraic conditions.Further, any compact portion of M 0 0 or M + 0 admits local (un)stable manifolds W s,u (M * 0 ), for * = 0, +, comprised of the union of the local (un)stable manifolds W s,u (M * (u)) of the equilibria M * (u), for * = 0, +.
To obtain the reduced dynamics on the critical manifolds, we consider (3.2) for δ = 0, given by where we substitute v = 0 (in the case of M 0 0 ) or v = v ± (u) (in the case of M ± 0 ) into the p-equation.By standard results of geometric singular perturbation theory, (restricting to the region u > 4bm in the case of M + 0 ) for all sufficiently small δ > 0 any compact portions of the normally hyperbolic invariant manifolds M 0 0 and M + 0 perturb to three-dimensional locally invariant manifolds M 0 δ and M + δ , which are C 1 -O(δ)-close to their singular counterparts.The slow flow on M 0 δ and M + δ is an O(δ)-perturbation of the reduced flow (3.8).Similarly, the local (un)stable manifolds W s,u (M * 0 ), * = 0, + perturb to four-dimensional locally invariant manifolds W s,u (M * δ ), for * = 0, + which are again O(δ)-close to their singular counterparts, and consist of the fast (un)stable fibers associated with orbits which lie on the slow manifolds M 0 δ and M + δ .

Fast transition layers
We aim to construct fast transition layers consisting of fast jumps between the critical manifolds M 0 0 and M + 0 .We return to the fast system (3.4) and the associated layer problem (3.5) for values of k ∈ [0, k].We note that the layer problem (3.5) corresponds to (2.4) in the stationary case c = 0. We recall from §2 that for values of u > 4mb, this problem has three equilibria, M 0 (u) and M ± (u), and at the critical value u = u f := 9bm 2 , the layer problem admits a heteroclinic loop between M 0 (u f ) and M + (u f ), while a homoclinic orbit to either M 0 (u) or M + (u) exists for values of u > u f or u < u f , respectively.As in §2.3, we assume that a > u f or equivalently, a m > 9b 2 , so that the equilibrium P 0 lies above the heteroclinic loop in the layer problem, and we further assume , so that the equilibrium P 2 lies on M + 0 and below the heteroclinic loop.The two heteroclinic orbits comprising the heteroclinic loop provide an opportunity for orbits to jump between the manifolds M 0 0 and M + 0 , and effectively transition between the desert/vegetated states.We recall from §2.3 that these two orbits φ dv (s) = (v dv , q dv )(s) and φ vd (s) = (v vd , q vd )(s) are given explicitly as and see Figure 4.As in §2.3, φ dv (s) represents the desert-to-vegetation state transition front which jumps from M 0 0 and M + 0 , while φ vd (s) represents the vegetation-to-desert state transition front which jumps from M + 0 and M 0 0 .These singular orbits serve as candidate interfaces between the desert and vegetated states.
Taking the union over values of the slow variables (p, r), or equivalently values of (p, k), the orbits φ dv form a two-parameter family of orbits lying in the intersection of the four-dimensional manifolds W u (M 0 0 ) and W s (M + 0 ), and likewise the orbits φ vd form a two-parameter family of intersections between W u (M + 0 ) and W s (M 0 0 ).In order to determine that these intersections are non-degenerate and persist in a suitable sense for δ > 0, we show transversality of the intersections in the remaining slow variable u.Lemma 3.1.Consider (3.4) for δ = 0 and b, m > 0, and fix p, k > 0. The following hold: (i) The manifolds W u (M 0 0 ) and W s (M + 0 ) intersect transversely along the three-dimensional manifold (ii) The manifolds W u (M + 0 ) and W s (M 0 0 ) intersect transversely along the three-dimensional manifold Proof.We focus on the statement (i), as the proof of (ii) is nearly identical.To prove transversality, it remains to show that the intersection breaks transversely when varying the remaining slow variable u near u = u f = 9bm 2 .
We accomplish this by computing the splitting of the manifolds W u (M 0 0 ) and ) is a solution of the fast layer equation (3.5) at u = u f .The adjoint equation associated with the linearization of (3.5) about φ dv at u = u f is given by and admits a unique bounded solution (up to multiplication by a constant) ψ dv (s) := (−q dv (s), v dv (s)).To leading order, the splitting distance of the manifolds W u (M 0 0 ) ∩ W s (M + 0 ) is determined to leading order in |u − u f | by the Melnikov integral where F (v, q; u) denotes the right-hand-side of (3.5).We compute from which we determine that the intersection The proof of (ii) proceeds similarly, and transversality is then determined by the Melnikov coefficient where ψ vd (s) := (−q vd (s), v vd (s)).In that case, we similarly find that Figure 6: Shown are the dynamics of (3.20) for δ = 0 on the far-field manifold M 0 0 .The manifold B far 0 is the set of all solutions which converge to the equilibrium p far 0 = (a, 0, 0) and remain bounded in the far field as k → 0.
We now construct the set of orbits which remain bounded in the far-field, and in particular those which converge to the desert state v = 0; these orbits will be asymptotic to the manifold M 0 δ .We first examine the reduced flow (3.8) on M 0 0 , given by (3.20) This system admits a single equilibrium p far 0 = (a, 0, 0) corresponding to the desert steady state of (1.1) in the far-field limit r → ∞.This fixed point is nonhyperbolic when considered as a fixed point of (3.20), but is of saddle-type when restricted to the invariant plane k = 0. Within M 0 0 , this fixed point admits a two-dimensional center-stable manifold B far 0 representing the set of solutions (u, p)(r) of (3.20) which remain bounded as r → ∞.To obtain a more explicit description of this set, we instead express (3.20) as the linear equation the solutions of which are given in terms of the modified Bessel functions I 0 , K 0 of the first and second kind.In particular, the unique solution of (3.21) which is bounded as r → ∞ and satisfies u(r) = ū for given ū ∈ R and r > 0 is given by For any r > 0, we can therefore express B far 0 as where K 1 (z) = −K 0 (z), and we note that since K 0 (z), K 1 (z) → 0 exponentially as |z| → ∞, the quantities K 0 (1/k), K 1 (1/k) are well defined (and converge to zero exponentially) as k → 0.
Within the stable manifold W s (M 0 0 ) of M 0 0 , we can construct the stable manifold of B far 0 as the set W s (B far 0 ) of stable fibers over trajectories in B far 0 .This set comprises the singular limit of all solutions which are bounded in the far-field.The set of orbits B far 0 perturbs within M 0 0 for sufficiently small δ > 0 to a two-dimensional manifold B far δ consisting of orbits within M 0 δ which are bounded as r → ∞ and in particular converge to the equilibrium p far 0 = (a, 0, 0).Likewise, as a subset of W s (M 0 0 ), the manifold W s (B far 0 ) perturbs to a three-dimensional invariant manifold W s (B far δ ) ⊂ W s (M 0 δ ) consisting of stable fibers lying over trajectories within B far δ ⊂ M 0 δ .The manifold W s (B far δ ) thus describes the set of solutions which remain bounded as r → ∞ and in particular those converging to the homogeneous equilibrium P 2 of (3.1).In light of the results of Lemma 3.1 in §3.1.2,the perturbed manifolds W u (M + δ ) ∩ W s (M 0 δ ) intersect transversely in a three-dimensional manifold which lies within O(δ) of the subspace u = u f .Viewed as a three-dimensional submanifold of W s (M 0 δ ), the manifold W s (B far δ ) is an O(δ)-perturbation of the singular manifold W s (B far 0 ) consisting of stable fibers over trajectories within B far 0 ⊂ M 0 0 defined as in (3.23).The manifold W s (B far 0 ) therefore transversely intersects W u (M + 0 ) in a two-dimensional manifold H far 0 ⊂ H vd consisting of the orbits within the subspace {u = u f }.This transverse intersection persists for all sufficiently small δ > 0, with W s (B far δ ) transversely intersecting W u (M + δ ) in a two-dimensional manifold H far δ which lies within O(δ) of H far 0 , and likewise O(δ)-close to the subspace {u = u f }.

The core region
In this section, we construct a three-dimensional manifold of orbits which remain bounded and converge to a set of uniformly vegetated states as r → 0.

3.2.1
The center-unstable core manifold W cu δ (C) For the core region, we consider r ∈ [0, r c ], where r c = δs c , or equivalently s ∈ [0, s c ], for some s c > 0 fixed independently of δ > 0. We use a blow-up rescaling z = log s to obtain the dynamics in the core region as Note that this system admits a two-dimensional manifold of equilibria at s = 0 defined by E = {p = q = s = 0}.When δ = 0, (3.25) becomes and near values of (u, v) satisfying mv − uv 2 (1 − bv) = 0, solutions of (3.26) which are bounded as z → −∞ can be expanded in terms of modified Bessel functions as follows.For given u > 4bm, there are three solutions Fixing one of these choices, we consider v ≈ v * (u), so that v = v * (u) + v, where |v| 1 is assumed small.
In the case of spots, we are interested in solutions near the vegetated state in the core, hence we take v * (u) = v + (u).Substituting into (3.26),we obtain where Our aim is to construct solutions of (3.25) which remain bounded as z → −∞ (s → 0).In view of (3.26), when δ = 0, any such solution must satisfy p ≡ 0, while the equation for (v, q) can be re-expressed in terms of the fast variable s as which, at the linear level, is a zero-order Bessel-type equation whose solutions can be expressed as linear combinations of the modified Bessel functions I 0 ( √ κs), K 0 ( √ κs) of the first and second kind.The function I 0 (ζ) is bounded at ζ = 0, while K 0 diverges logarithmically.Therefore, given u 0 > 4mb and s c > 0, we linearize (3.25) about the solution (u, p, v, q, s) = (u 0 , 0, v + (u 0 ), 0, e z ), integrate, and solve the resulting fixed-point equation in terms of the Bessel functions I 0 (•), K 0 (•).Defining the subset C := {p = q = s = 0, v = v + (u), u > 4bm} ⊂ E, this allows us to construct a local three-dimensional center-unstable manifold W cu δ (C) of solutions which are bounded as s → 0 (z → −∞), and in particular converge to C as s → 0. This manifold admits the expansion for sufficiently small |c 0 | and δ 1.Here I 1 = I 0 is the first-order modified Bessel function of the first kind.The manifold W cu δ (C) is parameterized by u 0 and c, both of which will be selected by intersecting with the three-dimensional far-field manifold W s (B far δ ) to obtain a spot solution in the full five-dimensional phase space.

Core transition region
In the previous subsection, we constructed the core center-unstable manifold W cu δ (C) of solutions which remain bounded as s → 0. Given s c > 0, we obtained expansions for the manifold W cu δ (C) valid for 0 ≤ s ≤ s c , provided δ is taken sufficiently small.We now aim to track this manifold into the region r = δs = O(1).
Hence we consider r ∈ [δs c , r 0 ] where s c is as above, fixed large independently of δ, and r 0 > 0 is fixed independently of δ, and δ is taken sufficiently small.We define the quantity δ := δs c , and we note that since s c 1 will be fixed independently of δ, in the limit δ → 0 we have that δ = O(δ) so that δ can be bounded as small as desired.We return to the fast system (3.1),appending an equation for r, which results in the following system.
We view this as a slow-fast system with timescale separation parameter 1/s c .In the region r ≥ δ, the variables u, p, r are slow, while v, q are fast.Rescaling s = s c ζ, we obtain the corresponding slow system Letting 1/s c → 0, this system admits a three-dimensional critical manifold M c 0 = {q = 0, mv = uv 2 (1 − bv)}, and the branch M c,+ 0 0 is normally hyperbolic, of saddle type.The reduced flow on M c,+ 0 is given by noting that the vector field is uniformly bounded in the region of interest r ∈ [ δ, r 0 ].In order to determine the dynamics in this region, we desingularize the system and rescale the independent variable dζ = r δ d ζ, which results in the system (3.33) The system (3.33)admits an invariant manifold at r = 0 with dynamics and a line of equilibria given by 0 := {(u, p) = (u 0 , 0), u 0 > 4bm} which are attracting within the manifold {r = 0}, each with a one dimensional stable manifold.In the normal (r) direction, this line of equilibria is repelling and admits a unique two-dimensional unstable manifold W u ( 0 ), which satisfies the following proposition, the proof of which follows by standard invariant manifold theory.Proposition 3.2.Consider (3.33).The line of equilibria 0 admits a unique two-dimensional unstable manifold W u ( 0 ), which for all sufficiently small r 0 > 0 can be represented as a graph W u ( 0 ) = {p = h 0 (u, r), 0 ≤ r < r 0 } where However, it does not suffice to restrict our attention to small values of r, and in fact r may need to be taken large.
To understand how solutions originating near r = 0 behave for large values of r, we need further information on the nonlinear vector field (3.33).In general, detailed estimates are not available for W u ( 0 ) when r is not small.However, more can be said near a value of u which satisfies a − u − uv + (u) 2 = 0. Note that this is the condition satisfied by the equilibrium P 2 of the full system (3.1)provided this equilibrium lies on the branch v = v + (u).
In this case, since a−U 2 −U 2 v + (U 2 ) 2 = 0, the line 2 := {u = U 2 , p = 0} is invariant.Examining the linearization of (3.33) about the invariant line 2 reveals a single zero eigenvalue with eigenvector (1, 0, 0), and a negative eigenvalue λ = −1, while the dynamics along 2 are simply r ζ = r.Therefore, there exists a two-dimensional, normally attracting manifold W c ( 2 ) which contains the line 2 , which can be represented as a graph where Moreover, in the region 0 ≤ r 1, the manifold W c ( 2 ) coincides with W u ( 0 ), and hence for simplicity we denote the union of these manifolds by W u ( 0 ).We emphasize that this (combined) manifold W u ( 0 ) is locally invariant and normally attracting.See Figure 7 for a depiction of W u ( 0 ).In the full system (3.30), the manifold M c,+ 0 perturbs to a three-dimensional slow manifold M c,+ sc which is O(1/s c )-close to M c,+ 0 , with slow flow given by an O(1/s c ) perturbation of the reduced flow (3.32).In particular, the manifold W u ( 0 ) perturbs within M c,+ sc to a locally invariant manifold W u sc ( 0 ).Furthermore, the four-dimensional stable/unstable manifolds W s,u (M c,+ 0 ) formed by the union of the stable/unstable fibers of basepoints on M c,+ 0 also perturb to stable/unstable manifolds W s,u (M c,+ sc ).We can identify the subset of these manifolds corresponding to the (un)stable fibers of basepoints on W u sc ( 0 ) as three-dimensional locally invariant manifolds W s,u (W u sc ( 0 )).We now track the core center-unstable manifold W cu δ (C) through the region r ∈ [δs c , r 0 ].We recall that, given any fixed (large) s c > 0, W cu δ (C) admits the expansion (3.29), which is valid up to s = s c for all sufficiently small δ > 0. At s = s c (corresponding to the subspace r = δ), W cu

Dynamics on M + δ
In the region r ≥ r 0 , where r 0 > 0 is taken sufficiently small and fixed independently of δ, we return to the fast system (3.1) and append an equation for r (3.37) In this region, when δ = 0 this system admits a critical manifold M 0 defined by (3.3), which can be decomposed into the branches M 0 = M 0 0 ∪M − 0 ∪F ∪M + 0 as in (2.8).For all sufficiently small δ > 0, (any normally hyperbolic portion of) M + 0 perturbs to a three-dimensional invariant manifold M + δ and its four-dimensional (un)stable manifolds perturb to four dimensional locally invariant manifolds W s,u (M + δ ).As a result of the analysis of the previous section, we know that W cu δ (C) approaches the set r = r 0 aligned along the strong unstable fibers of orbits on M + δ which are O(1/s c + δ)-close to the intersection W u ( 0 ) ∩ {r = r 0 }.By Proposition 3.2, this set is given by the graph and the projection of W cu δ (C) ∩ {r = r 0 } onto M + δ along the unstable fibers within W u (M + δ ) is therefore within O(1/s c + δ) of this graph.
Likewise, we consider the far-field stable manifold W s (B far δ ), which we recall from §3.1.3transversely intersects W u (M + δ ) in a two-dimensional manifold H far δ which lies within O(δ) of the set H far 0 given by where H vd is as in Lemma 3.1 and we recall k = 1/r.In other words, W s (B far δ ) intersects W u (M + δ ) transversely along the unstable fibers of orbits lying within O(δ) of the set where r > 0 is arbitrary.
We aim to show the existence of r = r I such that the manifolds W s (B far δ ) and W cu δ (C) intersect transversely at r = r I near the fast jump in the set {u = u f }.To do this, we will track orbits on W cu δ (C) as they evolve according to the dynamics of (3.37) until reaching the set {u = u f }.
In order to track W cu δ (C) through this region, we examine the reduced flow on M + 0 , given by where we've introduced the variable dξ = dr.At r = r 0 orbits of W cu δ (C) are aligned along the unstable fibers of orbits crossing Γ in .At u = u f , W s (B far δ ) are aligned along orbits crossing Γ out .Hence we aim to show in the reduced flow (3.41) that the forward evolution of trajectories in Γ in transversely intersects the set Γ out within the set {u = u f }.This transverse intersection will then persist under perturbation, thereby obtaining the transverse intersection of W s (B far δ ) and W cu δ (C) for sufficiently small δ > 0. We have the following proposition, which is the main result of this section.Proposition 3.3.Consider the set Γ in of initial conditions at r = r 0 for the system (3.41).The forward evolution of Γ in under the flow of (3.41) traces out a two-dimensional manifold Γ in which intersects the set {u = u f } in a curve Γ f .If the parameters a, b, m, satisfy then, within {u = u f }, there exists r = r I > 0 such that the curve Γ f transversely intersects Γ out at r = r I .
We begin with the following lemma which describes the set Γ out , given by the graph of the function p out (r) in (3.40).
Lemma 3.4.Regarding the function p out (r), the following hold.
(ii) lim Proof.For (i), we recall that p out (r) = a − u f K 0 (r) K 1 (r).Note that a − u f > 0. We have that Using the integral form for products of Bessel functions [8, §10.32.17], we see that since cosh(2t) − 1 ≥ 0 and K 0 (x) − K 2 (x) < 0 for all x ≥ 0, which completes the proof of (i).
We now describe the evolution of the set of initial conditions Γ in to the set u = u f .At r = r 0 , we represent the initial conditions Γ in via (3.38) as the graph p = p in (u) := h 0 (u, r 0 ).For a/m > max {9b/2, 4b + 1/b}, the function f (u) := −(a − u − uv + (u) 2 ) admits a unique zero u = U 2 ∈ (4bm, u f ), coinciding with the uniformly vegetated equilibrium state P 2 .We have the following lemma.
Lemma 3.5.Fix u in ∈ (U 2 , u f ).The forward evolution of the initial condition in Γ in given by p = p in (u in ) at r = r 0 eventually reaches the set {u = u f } at some value of (p, r) = (p f , r f )(u in ).
Proof.We show that for any initial condition u in ∈ (U 2 , u f ), under the flow of (3.41), that u ξ is nondecreasing and therefore eventually the u coordinate will reach u = u f at some r = r f (u in ).Since u ξ = p, we achieve this by showing that p ξ ≥ 0 along such a trajectory, which ensures that u ξ > p in (u in ) > 0. We note that p ξ > 0 initially at r = r 0 (via (3.38) and (3.41)), and at any location where p ξ = 0, we have that p ξξ = p r 2 + f (u)p > 0. Thus, u ξ = p is nondecreasing, which ensures that u will increase towards u in .
Taken over all values of u in ∈ (U 2 , u f ), we obtain a curve (p, r) = (p f , r f )(u in ) parameterized by the initial u-coordinate u in ∈ (U 2 , u f ).In order to prove Proposition 3.3, we show that the curve (p, r) = (p f , r f )(u in ) transversely intersects the curve p = p out (r) (that is, the set Γ out ) within the plane u = u f , for some value of u in and corresponding r = r f (u in ) =: r I .
Proof.We begin with the statement concerning the sign of r f (u in ).Consider two trajectories with two different initial conditions u in = u in,1 and u in = u in,2 , with u in,1 < u in,2 , which trace out solution curves (u, p) = (u(r; u in ), p(r; u in )).At r = r 0 , we have that p in (u in,1 ) = p(r 0 ; u in,1 ) < p(r 0 ; u in,2 ) = p in (u in,2 ) by (3.38).
We claim this implies p(r; u in,1 ) < p(r; u in,2 ) for all r > r 0 .Suppose for contradiction that p(r; u in,2 ) = p(r; u in,2 ), or for some r = r, which represents the first r value where the two trajectories cross.(Note that since p(r 0 ; u in,1 ) < p(r 0 ; u in,2 ), and p(r; u in ) is continuous, we know that there exists this first value r.) Thus p(r; u in,2 ) < p(r; u in,2 ) for all r < r, which implies that u(r; u in,2 ) < u(r; u in,2 ) for all r < r, since u in,1 < u in,2 and u ξ = p.Recall that p ξ = − p r + f (u).Then, since f (u) is an increasing function of u, we have that p ξ (r; u in,1 ) < p ξ (r; u in,2 ).(Note that this is a strict inequality since f (u) > 0.) This contradicts the fact that these solution curves intersect at r since p(r; u in,2 ) < p(r; u in,2 ) for all r < r.
Therefore p(r; u in,2 ) < p(r; u in,2 ) for all r.This also means that u (r; u in,1 ) < u (r; u in,2 ) for all r.Thus, since u in,1 < u in,2 , we also obtain that u(r; u in,1 ) < u(r; u in,2 ) for all r.Recall that u(r; u in,1 ) < u(r; u in,2 ) < u f for r < r f (u in,2 ).Thus, r f (u in,2 ) < r f (u in,1 ).Thus, whenever u in,1 < u in,2 , we have that r f (u in,2 ) < r f (u in,1 ), so that r f (u in ) is a strictly decreasing function.
We now turn to the sign of p f (u in ).Again, we consider two trajectories with initial conditions u in = u in,1 , u in,2 , with u in,1 < u in,2 .Suppose for contradiction that p f (u in,1 ) ≤ p f (u in,2 ).Since u (r; u in,1 ), u (r; u in,2 ) > 0, there exists a least u = ũ < u f and r1 , r2 < r f at which u(r 1 ; u in,1 ) = u(r 2 ; u in,2 ) = ũ and p(r 1 ; u in,1 ) = p(r 2 ; u in,2 ) = p.Since f (u) > 0, and using a similar argument as above, we have that r1 > r2 .
We express the solution curve (u, p) = (u(r; u in,1 ), p(r; u in,1 )) as a graph u = u 1 (p) over p, and similarly the curve (u, p) = (u(r; u in,1 ), p(r; u in,1 )) as u = u 2 (p).Then at p = p, we have from which we see that du 1 dp < du 2 dp , which contradicts the fact that u 1 (p) must cross u 2 (p) from below.Thus, we conclude that for any u in,1 < u in,2 , we have p f (u in,1 ) > p f (u in,2 ), so p f (u in ) is a strictly decreasing function of u in .
Lemma 3.6 above shows that p f (u in ) and r f (u in ) are both strictly decreasing functions of u in ∈ (U 2 , u f ).We now consider the limiting behavior of (p f , r f )(u in ) as u in approaches the limits u in = U 2 and u in = u f .In preparation, we consider (3.41) in the limit r → ∞, resulting in the vector field which admits the conserved quantity Note that E(U 2 , 0) = 0, and define p f,∞ to be the unique positive solution of E(u f , p f,∞ ) = 0, corresponding to the intersection of the unstable manifold of the saddle equilibrium (u, p) = (U 2 , 0) of (3.43) with the set u = u f .We have the following.
(i) lim Proof.The limit (i) follows directly from the definition of (p f (u in ), r f (u in )) and (3.38).
For (ii), we aim to compute the limit lim with the invariant line 2 corresponding to the fixed point (u, p) = (U 2 , 0).Hence to determine the behavior of trajectories lying on Γ in with values of u ≈ U 2 , we can track such trajectories along the invariant line 2 to large values of r.In particular, for any fixed R 0 1, there exists δ R0 such that the forward evolution of Γ in traces out a two dimensional manifold Γ in which contains the invariant line 2 and intersects the plane r = R 0 in a curve which can be represented as a graph p = h We set r = 1/k and arrive at the system The invariant set k = 0 (corresponding to r = ∞) contains the limiting system (3.43) for the variables (u, p), whose solutions lie on level sets of the function E(u, p), with the saddle-type equilibrium (u, p) = (U 2 , 0) satisfying E(U 2 , p) = 0.The two branches of this level set correspond to the one-dimensional stable and unstable manifolds W u/s,∞ (U 2 , 0) of the equilibrium within the invariant set k = 0, which are tangent to the lines p = ± f (U 2 )(u − U 2 ).The manifolds W u/s,∞ (U 2 , 0) extend to two-dimensional center-stable/center-unstable manifolds W cu/cs,∞ (U 2 , 0) for small k 1 which intersect along the invariant line 2 = {u = U 2 , p = 0}.Let k 0 := 1/R 0 1.Recall the manifold Γ in intersects the set k = k 0 in a a curve which can be represented as a graph over |u − U 2 | < δ R0 given by p ) at the linear level and transversely intersects W cs,∞ (U 2 , 0) along the invariant line 2 .Thus as k → 0, Γ in aligns along the branch of the level set E(U 2 , p) = 0 corresponding to W u,∞ (U 2 , 0), and contains the invariant line 2 .From this, we see that as u in → U 2 , r f (u in ) → ∞ and the corresponding solution approaches the manifold W u 0 (U 2 , 0), so that p f (u in ) → p f,∞ as claimed.
Combining this with the results of Lemma 3.4, we are able to complete the proof of Proposition 3.3.
Proof of Proposition 3.3.As described above, the forward evolution of Γ in reaches the set u = u f in a curve (p f , r f )(u in ) parameterized by u in ∈ (U 2 , u f ).By Lemma 3.6, p f (u in ) and r f (u in ) are both strictly decreasing functions of u in ∈ (U 2 , u f ), so we can express the curve (p f , r f )(u in ) as a graph Furthermore, by Lemma 3.4, within the set u = u f , Γ out is given by a graph p = p out (r) which satisfies p out (r) < 0 for r ∈ (0, ∞) with lim r→0 p out (r) = ∞ and lim Therefore, in order for the sets Γ in and Γ out to intersect transversely at some r = r I , it only remains to check whether p f,∞ > a − u f .Equivalently, a transverse intersection occurs provided or equivalently (3.48)

Proof of Theorems 1.2-1.3
The results of the preceding sections §3.1-3.3 allow us to complete the construction of radial spot solutions in (1.1).
Proof of Theorem 1.2.As in §3.3, we track the far-field manifold W s (B far δ ) into a neighborhood of M + δ , where it aligns along the unstable fibers of orbits which cross the set {u = u f } within O(δ) of the set Γ out .Likewise, the manifold W cu δ (C) of solutions bounded at the core can be tracked into a neighborhood of M + δ , where it aligns within O(1/s c + δ) of the unstable fibers of orbits crossing the set {r = r 0 } along the curve Γ in .By Proposition 3.3, the forward evolution of Γ in reaches the set u = u f in the curve Γ f which transversely intersects Γ out at some r = r I , provided Therefore, the manifolds W s (B far δ ) and W cu δ (C) intersect transversely, corresponding to a radial spot solution bounded on r ∈ [0, ∞), with a single sharp interface occurring at r = r I + O(δ).The value V c (a, b, m) is determined by the coordinate v + (u 0 ) in corresponding fiber of W cu δ (C) in the limit δ → 0; see (3.29).Finally, the condition (1.7) can be obtained from a lengthy but straightforward computation by carrying out the integration in (3.49) and using the steady equation satisfied by U 2 .
Regarding gaps, we similarly complete the proof of Theorem 1.3.
Proof of Theorem 1.3.The argument is similar to that of Theorem 1.2.We briefly outline the differences which result in the opposite condition (1.8).
The geometry of the construction is quite similar, except opposite in that the far field manifold consists of solutions asymptotic to the equilibrium (U 2 , V 2 ) on the critical manifold M + δ , and the core manifold consisting of solutions bounded as r → 0 is constructed from orbits originating near the desert state (U 0 , V 0 ) on M 0 δ .In this case building the core manifold is somewhat less involved since the flow on M 0 δ is linear, and the solutions are given explicitly in terms of modified Bessel functions.However, in the far field one must deal with the nonlinear flow on M + δ .In this case, the manifold W cu δ (C) intersects W s (M + δ ) in the set {u = u f } transversely along the unstable fibers of orbits lying within O(δ) of the set where r > 0 is arbitrary.We note that, using asymptotic properties of modified Bessel functions, it can be shown similarly as in §3.3 that p in is an increasing function of r with lim r→0 p in (r) = 0 and lim In the far field, analogously to the case of spots, we can construct the two-dimensional far field manifold B far 0 within M + 0 as the set of solutions of the system (3.41) which remain bounded as r → ∞.In the limit r → ∞, this system approaches the system (3.43), the solutions of which are given by level sets of the conserved quantity (3.44).Let p far (r I ) denote the p coordinate at r = r I of the solution which is bounded as r → ∞ and satisfies u(r I ) = u f .Then using similar arguments as in §3.3, we see that p far is an increasing function of r I which satisfies lim Thus in order to have an intersection of W cu δ (C) and the stable fibers W s δ (B far δ ), and thus a radial gap solution with a single interface at some value of r = r I , a sufficient condition is which is precisely the opposite condition as that which guarantees the existence of spots.The estimate for v g (r) as r → 0 is due to the exponential decay along the fast fibers of W u (M 0 δ ), and the fact that the subspace {v = q = 0} is invariant under the flow of (3.25) for δ > 0.

Rings, targets, and other radially symmetric solutions
The techniques used in §3.1-3.4 to construct spot and gap solution could be used to construct other localized solutions with radial symmetry as follows: The general strategy is the same, in that to construct a solution which is asymptotically constant, and bounded as r → 0, as in the case of spots/gaps, we construct a core manifold W u δ (C) of states originating near one of the desert or vegetated steady states (U 0 , V 0 ) or (U 2 , V 2 ).We similarly construct a far-field stable manifold W s δ (B far δ ) of solutions bounded as r → ∞ consisting of stable fibers of one of the steady states (U 0 , V 0 ) or (U 2 , V 2 ), which could be the same or different from the state near the core.
Then, to construct a radially symmetric profile with a desired number of interfaces, the core manifold W u δ (C) is tracked along a number of fast jumps alternating between M 0 δ and M + δ as a sequence of radii r j , j = 1, 2, 3, . .., in between which W u δ (C) follows the slow flow of the corresponding slow manifold M 0 δ or M + δ , entering along its fast stable fibers, and exiting aligned along its fast unstable fibers according to the exchange lemma.This procedure could be used, in principle, to construct ring or target profiles with any desired (finite) number of interfaces.Some examples obtained numerically are presented in §5.

Spot instabilities
In this section, we examine the stability of the spot solutions from Theorem 1.2, and in particular we demonstrate several instabilities exhibited by these solutions when considering 2D perturbations.As we are primarily interested in demonstrating potential instability mechanisms, we do not take a rigorous approach, but rather employ formal asymptotic arguments; however, we emphasize that rigorous results could be obtained using similar methods as in the existence analysis in §3.
We linearize (1.1) about a radial spot solution (u sp , v sp )(r) = (u sp , v sp )(r; a, b, m, δ) of Theorem 1.2 using an ansatz of the form for ∈ Z, which results in the eigenvalue problem We consider the essential spectrum associated with spot solutions in §4.1.The point spectrum for wave numbers , and is evaluated asymptotically in the limit of spots of large radius in §4.3, where the sideband stability from nearby planar front solutions (see §2.4) is recovered in the limit r I 1.Finally, the point spectrum for large wavenumbers | | 1 is considered in §4.4.

Essential Spectrum
The essential spectrum associated with the radial spot solution (u sp , v sp )(r; a, b, m, δ) is determined by considering the limit r → ∞ in (4.1), and computing the 1D essential spectrum of the asymptotic rest state lim Proof.Letting r → ∞, and writing (4.1) as a first order system, we obtain The essential spectrum Σ ess consists of λ ∈ C for which the matrix A ∞ is not hyperbolic.A short computation shows that this can only occur in the region {λ ∈ C : Re λ ≤ −β}, where β = min{1, m} > 0.

Point spectrum for | | = O(1)
Near the interface r = r I , we change variables to r = r I + δs for s ∈ (−ν| log δ|, ν| log δ|) for some ν 1. Due to the exponential convergence of the front of the fast subsystem between the critical manifolds M 0 0 and M + 0 , for ν chosen sufficiently large, this interval captures the portion of the spot solution which lies outside an O(δ)-neighborhood of the slow manifolds M 0 δ and M + δ .In this region, the eigenvalue problem (4.1) becomes For wavenumbers | | = O(1) with respect to δ, we expand solutions of this eigenvalue problem in terms of the reduced fast system which is an eigenvalue problem of Sturm-Liouville type, obtained by linearizing the fast subsystem about the front φ vd , and which has a solution given by the derivative v vd when λ = 0.For (4.3), we expand the eigenfunction and eigenvalue parameter λ( ) = δλ 1 ( ) + O(δ 2 ), as well as the solution Substituting into (4.3),we obtain to leading order where Considering the first equation of (4.3), we write as a first order system so that ū1 is constant to leading order.We now expand the existence problem across the fast jump near r ≈ r I as and differentiate the second equation with respect to s to obtain to leading order Substituting into the second equation of (4.3), we have to leading order Using the fact that L 0 is self adjoint, we take the inner product of this equation with v vd , and obtain the solvability condition .
Since v vd is strictly negative, the sign of λ 1 is given by the sign of the prefactor ((u 1 ) s − ū1 ).The value of (u 1 ) s is easily determined through the expansion of the existence problem since so that (u 1 ) s = a − u f K 0 (r I ) K 1 (r I ) corresponding to the leading order p-value across the fast jump at r = r I .From this we obtain the solvability condition It remains to determine the constant ū1 in (4.11).To determine ū1 , we recall from (4.7) that ū1 is constant, whilst p1 satisfies to leading order where p0 1 , p+ 1 denote the limiting values of p1 on either side of the fast jump, when the solution approaches the critical manifolds M 0 0 , M + 0 , respectively.To determine ū1 , we construct bounded eigenfunctions (ū 0 , p0 )(r) and (ū + , p+ )(r) in the slow regions near M 0 0 , M + 0 , respectively, such that across the fast jump at r = r I , we have ū0 (r I ) = ū+ (r I ) p0 (r I ) = p+ (r (4.12) We begin by analyzing the linearized equation in the slow variables on M 0 0 .We note that here v sp (r) = 0 to leading order; inspecting (4.1) and recalling λ = δλ 1 , we obtain that v = 0 to leading order on M 0 0 , so that u satisfies the leading order equation which is a modified Bessel's equation.For each , this equation admits a unique solution which is bounded as r → ∞, which is the modified Bessel function of the second kind K (r).We therefore obtain the leading order solution (ū 0 , p0 )(r) = (αK , αK )(r) in the slow region r > r I .
Near M + 0 , the slow reduced equations are not as straightforward, due to the nonlinear reduced flow on M + 0 .In particular inspecting (4.1), to leading order we have that v satisfies where v + (r) := v + (u + (r)) and u + (r) is the solution in the slow region for r < r I .Substituting into (4.1),we obtain the leading order equation for u in the slow region u. (4.15) The solutions of this equation do not appear to have a nice representation in terms of special functions.However, this equation is of the form where f + (r) has a well-defined limit as r → 0. As r → 0, the equation behaves like a Bessel-type equation, and it is possible to show that there is a unique (up to a constant multiple) solution u = u * (r) for each which is bounded as r → 0. To see this, we rewrite (4.16) as where η = log r.The system (4.17) has a fixed point at the origin which admits a two-dimensional unstable manifold corresponding to a one-dimensional space of solutions of the nonautonomous linear system (4.16) which are bounded as η → −∞ (r → 0).This space is spanned by a nontrivial solution, which we denote by u = u * (r), which will serve as a candidate eigenfunction in the slow region r < r I .We therefore obtain the leading order solution (ū + , p+ )(r) = (βu * , β(u * ) )(r) in the slow region r < r I .Using the conditions (4.12), we obtain which we can solve to determine and so In general, we require information about the solution u * to be able to determine the sign of this quantity as a function of .This is nontrivial to do in general, as u * likely does not have a direct representation in terms of special functions.However, in certain limiting cases we can approximate (4.20).For sufficiently large spots r I 1, we argue in §4.3 that such spots inherit instabilities from nearby planar front solutions (see §2.4).We consider the case of large wavenumbers | | 1 in §4.4.

Large spots: recovering the sideband instability
In this section, we consider the critical eigenvalue expression (4.20) in the case of a very large radial spot solution, that is r I 1. Near the core, such a solution is approximately constant, while at the interface, the solution resembles a stationary planar front between the desert and vegetated equilibrium states.In the limit r I → ∞, in the far field the solution approaches the stationary planar front and inherits the (in)stability properties of the front.To see this, in this section we estimate the expression (4.20) in the asymptotic limit r I → ∞.
Remark 4.2.To investigate this limit, one option would be to apply the approach from §2.4 to (4.3) under the assumption r I 1, as for the stability of traveling fronts.However, since we do not intend to repeat the analysis from [4] which results in the expression (2.28), in this section we provide a more direct method by estimating the expression (4.20) in the asymptotic limit r I → ∞.
We estimate the expression (4.20) for finite values of ∈ Z as r I → ∞.The expression is explicit (in terms of special functions) except for the value of u * (r I ), where u * (r) is the unique bounded solution (up to a constant) of (4.15) as r → 0. In the case of large r I 1, we can approximate u * (r I ) as follows.Note that when = ±1, (4.15) reduces to which admits a solution bounded at r = 0, given by the derivative u + (r), where u + (r) is the core solution on . Using reduction of order, we can find another linearly independent solution of this equation, given by and the Wronskian of u 1 := u + and u 2 is given by W (u 1 , u 2 )(r) = r −1 .We now assume r I 1 and attempt to construct the leading-order bounded solution of (4.15) on the interval [0, r I ] as r I → ∞.We split this interval into [0, Since the vector field for the existence problem on M + 0 for large r is given by (3.43), by taking R sufficiently large, we can arrange for (u where to leading order in 2 − 1 r 2

I
. We write the solution of this system using variation of constants as , and using (4.15), we have that w = w * (r) satisfies the equation where Appending an equation for r and rescaling the spatial coordinate, we have the equivalent autonomous system shows that Finally, returning to (4.29), we have that 1 sufficiently large (but O(1) with respect to δ), so that spots of radius r I = O(1) are always unstable.In particular, the spots we have constructed in §3 are unstable to (suitably) large wave numbers .However, for sufficiently small spots (r I = o(1) with respect to δ), the above argument is no longer valid, and it is not possible to rule out stable spots; see §5 for some solutions obtained numerically.However, the regime r I = o(1) lies outside the scope of the existence analysis in §3.
Next, we note that if = ¯ δ −1/2 for some 0 < ¯ = O(1), a similar analysis as in §4.2 results in the solvability condition (in place of (4.20)) so that spots are always unstable to wavenumbers of the form = ¯ δ −1/2 where ¯ is small but O(1).For larger wavenumbers, the first term in (4.33) dominates so that λ 1 ( ) < 0. Thus a switch from unstable to stable wavenumbers occurs at some critical = ¯ δ −1/2 .
Finally, if = ¯ δ for some 0 < ¯ = O(1), to construct a bounded solution of (4.3), to leading order we must have u ≡ 0, and the fast equation in (4.3) 8 (bottom left panel), which develops similar finger-type patterns.Simulations were performed in Matlab using finite differences for spatial discretization with periodic boundary conditions, and Matlab's ode15s routine for time integration.that these spectral computations show agreement with the analysis in §4.4, in that the spots/gaps are unstable for a range of 'large' wave numbers (note here that 1/ √ δ ≈ 4.47), and that λ( ) becomes negative for sufficiently large | |.
A natural question concerns the nature of these linear instabilities in the nonlinear dynamics of the spots/gaps.In the large radius limit, we expect such solutions to inherit the sideband instability of the nearby stationary front; in [4], it was demonstrated that this sideband instabilty can lead to the appearance of finger-like patterns along the front interface, which can in turn lead to labyrinthine patterns which expand spatially into the homogeneous states.By performing direct numerical simulations using the unstable spot and gap solutions from Figure 8 as initial data, we see a similar instability manifest along the (circular) interface; see Figure 9 for snapshots of these simulations.We leave a more detailed study of the appearance of such finger-like patterns, and the relation to the corresponding instabilities in the stationary front interface to future work.
While the spots and gaps of Theorems 1.2-1.3 are unstable for radii r I = O(1) (or larger) with respect to δ, the relation (4.27) is no longer valid when r I is not large, and the analysis in §4.2 is not valid if r I = o(1) as δ → 0. Hence it may be possible to find smaller spots or gaps which are stable.Figure 10 depicts spot and gap solutions of smaller radii (but nearby in parameter space to those in Figure 8), for which we see that λ 1 ( ) is no longer well approximated by (4.27).We see in this case that λ( ) is negative aside from the double zero eigenvalue λ(±1) = 0 due to translation invariance.Figure 11 shows a continuation of the eigenvalues λ 1 ( ), = 0, 1, 2, 3, 4, 5 for decreasing r I for these spot and gap solutions, where we observe that each eigenvalue eventually stabilizes as r I decreases.Additionally, as described in §3.5, we are similarly able to find (seemingly) stable radially symmetric ring and target patterns (see Figure 10), which could be obtained in a similar manner to the spots/gaps of Theorems 1.2-1.3 by constructing solutions with several sharp interfaces at distinct radii r i .
Obtaining the stability of spots of smaller radii rigorously appears to be a challenging problem, and we leave this to future work; in particular our existence analysis does not immediately extend to this regime.Additionally, unlike prior works which considered the stability of radially symmetric solutions using singular perturbation methods [41,42] in a 3-component FitzHugh-Nagumo system, the solutions of the linearized equations (4.15)  do not have explicit representations in terms of special functions, which makes it difficult to determine λ 1 ( ) for smaller values of r I .This is related to the challenges which arise in the existence analysis in §3 due to the nonlinear reduced flow on the slow manifold M + δ .Another natural line of further research concerns investigating whether the present insights obtained on specific model (1.1) can be lifted to the general setting of the 2-component singularly perturbed reaction diffusion systems considered in the companion paper [4].In that paper, we study the (in)stability of planar fronts with respect to longitudinal perturbations, as we did for (1.1) in section §2, and derive two general, and relatively simple, criteria on the emergence of the sideband instability mechanism (that typically leads to finger-like patterns).The present analysis indicates that the same mechanism may drive the (in)stability of large spots and gaps in the general setting of [4], although one should not underestimate the technicalities involved in establishing the counterparts of Theorems 1.2 and 1.3 for general model.It is less from the combination of the present insights and those [4], under conditions spots gaps with a radius r I of O(1) will be unstable (as is the case for (1.1)).The nature of the analysis in section §4.4 suggests that it is possible to derive a general (in)stability result for spots and gaps with radius r I = O(1) against perturbations with | | 1.This suggests that also in the general setting, spots and gaps with radius r I 'sufficiently small' are potentially the 'most stable' (radially symmetric) localized patterns.Thus, the issue of the existence and stability of spots and gaps of sufficiently small radius is a central question and resolving that question may explain the abundance of 'spikes' in the literature on localized patterns in singularly perturbed reaction-diffusion systems -see [5,23,44] and the references therein and Remark 1.1.These spikes have a fully homoclinic nature, in the sense that they are not close to a concatenation of almost heteroclinic orbits: away from the slow manifold M 0 δ they only follow the fast (spatial) dynamics, they do not follow the slow flow on a second slow manifold M + δ -as is the case for the patterns constructed here (with r I = O(1)).Thus, by studying spot and gap patterns with radius r I decreasing from being O(1) to asymptotically small, one needs to zoom in on the subtle process through which a localized pattern detaches from M + δ during its jump away from and back to M 0 δ -see also [21].Lastly, we briefly describe the appearance of far-from-onset spot patterns in (1.1).While the analysis of Theorems 1.2-1.3only applies to the construction of a single spot or gap solution, we anticipate that one could construct periodic spot or gap patterns by tiling the plane with well-separated copies of the primary spot or gap solution.See Figure 12 for results of direct numerical simulations in (1.1) which result in the appearance of spatially periodic spot and gap lattice patterns.While the construction of such patterns is beyond the scope of this work, we note that a spatial dynamics approach, such as that described in [36] could be used to construct such large amplitude spot patterns.However, the question of stability of the resulting patterns is likely very challenging.

Figure 2 :
Figure 2: Plotted are the curves a m = 9b 2 (black), a m = 4b + 1 b (dotted red), and a m = 9b 2 + 2 b (dashed blue).The shaded region corresponds to the region in parameter space where the condition (1.5) is satisfied.

Figure 7 :
Figure 7: Shown is the slow flow on M c,+ sc .The manifold W cu δ (C) aligns along the unstable fibers W u (C c,+ sc ) of the manifold C c,+ sc in the full system (3.30).Within M c,+ sc , C c,+ sc aligns along W u sc ( 0 ) under the forward evolution of (3.30).
δ (C) is aligned along the unstable fibers within W u (M c,+ sc ) of base point orbits on M c,+ sc lying on a curve C c,+ sc = {(u, p, r) : r = δs c , p = O(δ, 1/s c )}.In particular W cu δ (C) transversely intersects the stable fibers of these orbits within W s (M c,+ sc ).Tracking under the forward-flow of (3.30), by the exchange lemma, W cu δ (C) aligns O(e −sc )-close to the unstable fibers within W u (M c,+ sc ) of the forward evolution C c,+ sc of the manifold C c,+ sc within M c,+ sc ; see Figure 7.As the flow on M c,+ sc is an O(1/s c ) perturbation of the reduced flow (3.32), we are thus able to determine how W cu δ (C) emerges at r = r 0 , noting that for larger values of r, we only have detailed estimates on W cu δ (C) near the solution u = U 2 of a − u − uv + (u) 2 = 0.

. 25 ) 1 r 2 I
Recalling u * (r) = u 1 (r)+ 2 − ũ(r), in order to construct a bounded solution as r → 0, we must have (u, u )(r) ≈ (C 3 I , C 3 I )(κ 0 r) at r = r I − R for some constant C 3 in order to match with the bounded solution I on the interval [0, r I − R].Using the fact that u 1 (r I − R) decays exponentially in R as R → ∞, while u 2 (r I − R) grows exponentially as R → ∞ and I ( √ κ 0 (r I − R)) ∼ e − √ κ0R I ( √ κ 0 r I ) for 1 R r I , we see that we must choose C 2 ≈ 0 in order to ensure u(r) remains bounded as r → 0. Thus we obtain the solution u ) which, for large | | 1, admits two invariant manifolds which are defined up to r = 0, given by w ± (r; ) = ±1 + O(| | −1 ).The manifold w = 1 is attracting while w = −1 is repelling.For large | |, in order for u * (r) to be bounded as r → 0, the solution w * (r) must lie on the manifold w + (r; ), and hence w * (r) ∼ 1 for large | |.Therefore, (u * ) (r) u * (r) ∼ | | r as | | → ∞.A similar argument (or using the asymptotic expressions in [8, §10.41])

. 34 ) 2 r 2 I. 2 r 2 I.
which is just the Sturm-Liouville problem (4.4), but with λ shifted by ¯ The problem (4.4) has an eigenvalue at 0 due to translation invariance of the front v vd , while any other eigenvalues are bounded away from the imaginary axis.Hence we have that any eigenvalues of (4.34) lie in the region λ ∈ C : Reλ ≤ − ¯ Thus there are no further instabilities of λ( ) for | | = O(1/δ) or larger.

Figure 8 :
Figure 8: Stationary radial profiles obtained in (1.3) for b = 1.0, m = 0.5, δ = 0.05 corresponding to a spot solution (top panels, a = 2.625) and a gap solution (bottom panels, a = 2.665).The left panels depict the corresponding radial profiles (u, v)(r) (u-profile in blue, v-profile in green), while the right panels depict the corresponding critical eigenvalues (blue dots) λ( ) for −12 ≤ ≤ 12. Also plotted (red) is the critical eigenvalue curve for a (slowly) traveling front found for the same parameter values: for a = 2.625 (top), the corresponding front has speed c = 0.012, while for a = 2.665 (bottom), the front has speed c = −0.013.The radial profiles were obtained by solving the stationary equation (1.3) using Matlab's fsolve routine, where finite differences were employed for the spatial discretization with Neumann boundary conditions.The eigenvalues λ( ) were obtained by linearizing (1.1) about the radial profile and using Matlab's eigs routine.

Figure 9 :
Figure 9: (Upper panels) Snapshots of direct numerical simulation of a spot solution for a = 2.625, b = 1.0, m = 0.5, δ = 0.05, with initial data given by the solution in Figure 8 (top left panel).The spot develops finger-type patterns along the interface which spread throughout the domain.(Lower panels) Snapshots of direct numerical simulation of gap solution for a = 2.665, b = 1.0, m = 0.5, δ = 0.05, with initial data given by the solution in Figure 8 (bottom left panel), which develops similar finger-type patterns.Simulations were performed in Matlab using finite differences for spatial discretization with periodic boundary conditions, and Matlab's ode15s routine for time integration.

Figure 10 :
Figure 10: Radially symmetric solutions obtained for the parameter values = 1.0, m = 0.5, δ = 0.05 for different values of a. From top to bottom: spot (a = 2.55), gap (= 2.765), ring (a = 2.538), and target pattern (a = 2.78).The various solutions survive in direct numerical simulations suggesting they are indeed stable.For each row, the left panel depicts the planar profile obtained by direct numerical simulation in (1.1) using finite differences for spatial discretization with periodic boundary conditions.The middle panel depicts a radial profile obtained by solving (1.3) using finite differences and Neumann boundary conditions, and the right panel depicts the critical eigenvalue λ( ) for −12 ≤ ≤ 12, obtained by linearizing (1.1) about the radial profile and using Matlab's eigs routine.

Figure 11 :
Figure11: Shown is a numerical continuation of the eigenvalues λ 1 ( ) for = 0, 1, 2, 3, 4, 5 for a spot solution (left) and gap solution (right) as a function of the interface radius r I .The curves were obtained by continuing the unstable spot and gap solutions from Figure8to the stable spot and gap solutions of Figure10by adjusting the parameter a for fixed (b, m, δ).The interface location r I was computed at each step by approximating the inflection point of the vegetation profile of the corresponding solution.From the figure, we see that the unstable eigenvalues λ 1 ( ), = 2, 3, 4, 5 eventually stabilize as r I decreases.