{'fraction_non_alphanumeric': 0.0475238922675934, 'fraction_numerical': 0.018360845641471186, 'mean_word_length': 4.368130586863583, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Existence of localised radial patterns in a model for dryland vegetation

Localised radial patterns have been observed in the vegetation of semi-arid ecosystems, often as localised patches of vegetation or in the form of `fairy circles'. We consider stationary localised radial solutions to a reduced model for dryland vegetation on flat terrain. By considering certain prototypical pattern-forming systems, we prove the existence of three classes of localised radial patterns bifurcating from a Turing instability. We also present evidence for the existence of localised gap solutions close to a homogeneous instability. Additionally, we numerically solve the vegetation model and use continuation methods to study the bifurcation structure and radial stability of localised radial spots and gaps. We conclude by investigating the effect of varying certain parameter values on the existence and stability of these localised radial patterns.


Introduction
The self-organisation of vegetation in semi-arid ecosystems remains one of the most visually striking examples of pattern formation in nature. Non-uniform distributions of vegetation are thought to result from symmetrybreaking instabilities induced by a lack of resources, such as a scarcity of water, combined with a positive feedback between local vegetation growth and water transport towards the growing vegetation; see [38] for a detailed review of the mechanisms involved in the formation of vegetation patterns, and see Figure 1 for a schematic of how vegetation patterns form and are measured. Semi-arid ecosystems are defined to have an annual rainfall of 250-500 mm; these ecosystems are not rare, with 90% of Israel defined to be semi-arid or arid, and 30% of the Earth's land surface estimated to be semi-arid [53]. With the advent of aerial [33,34] and satellite imagery, and an increasing public awareness of ecological issues such as desertification, the study of vegetation patterns has never been more accessible or relevant; this is reflected in the wealth of new research being completed in this area (see [4,19,27,47,50] for recent work in mathematical modelling, [36,58] for experimental and numerical simulations, and [14,15,43] for remote sensing from satellite imagery).
While sloped terrains often produce slowly migrating stripe patterns [44,48,53], flat terrains are known to produce a variety of different steady structures [6,23,24,41]. These patterns include localised patches, where sparse regions of vegetation are surrounded by bare soil [40], and 'fairy circles', where a sparse region of bare soil is surrounded by uniform vegetation (as seen in the Namib desert [50, Figure 1.]). A single patch or fairy circle is locally axisymmetric and spatially localised; hence, we restrict our focus to localised radial patterns in order to better understand these structures.
There has been a number of analytical and numerical studies of localised vegetation patterns; see [4,5,40,50] for nonlocal mean-field models, [4,5,31] for non-variational PDE models, and [16] for a variational nonlocal a) b) c) Figure 1: a) Due to evaporation of the soil water (illustrated by the blue arrows) and insufficient resources to sustain uniform plant growth, semi-arid ecosystems can exhibit localised patches of vegetation. The soil water is also localised, subsisting in the shadow of the vegetation. b) For sufficiently large domains, these patches can also form part of a larger pattern (such as spots, stripes or hexagons, for example). c) By measuring the vegetation and soil water densities, vegetation patterns can be described by continuum models. Here, we present a schematic density plot of the vegetation in panel b), where the darker green indicates a higher intensity of vegetation.
Nagumo-type model. Notably, Jaïbi et al. [27] rigorously analysed a 2-component reaction-diffusion model in order to prove the existence of localised travelling fronts, as well as stationary spot and gap patterns. This work focused on one-dimensional patterns and utilised singular perturbation theory in order to complete the proof; it may be possible to extend these results for the radial case, however the associated equations will explicitly depend on the radial coordinate, which may introduce some extra technical challenges. Instead, we will take a weakly nonlinear approximation of the reaction-diffusion model introduced by von Hardenberg et al. in [55], which we term the 'von Hardenberg' model, and draw connections to prototypical pattern-forming models, which we will discuss presently.
Within the study of localised states, we wish to highlight an important distinction between localised spikes and localised spots; here we are using the terminology set out in [11], related to the spatial eigenvalues of the surrounding homogeneous state at the bifurcation point. We note that, for radial problems, we typically examine the spatial eigenvalues in the limit as r → ∞. A localised spike is commonly the result of spatial eigenvalues bifurcating from the origin onto the real line; these solutions do not have an associated wave number, and so they have monotonic decay in their tails; see Figure 2 a). In contrast, a localised spot is often the result of purely imaginary eigenvalues bifurcating off of the imaginary axis; such solutions have an associated wave number, and so they decay with spatially-oscillating tails; see Figure 2 b). Localised patterns associated with this bifurcation often possess much more complicated bifurcation curves than localised spikes and are known to sometimes undergo homoclinic snaking, where infinite families of localised patterns are connected via a snaking process, as seen in [57]. For more information about homoclinic snaking, see [8] and references therein.
In the context of vegetation patterns, localised spikes are typically found far from any symmetry-breaking instability, while localised spots are found close to such an instability (see [11] for a review of such structures). This instability is often called a 'Turing' bifurcation after Turing's seminal work on pattern formation [51]. Since we will be studying a reaction-diffusion model, such that the instabilities of the system are diffusiondriven, we will stick to this naming convention; we will also refer to the spike-inducing bifurcation, where real eigenvalues are bifurcating from the origin, as a 'homogeneous' bifurcation. A notable difference between localised spikes and spots can be seen in how they interact; see [4,5,49,50]. Berrios-Caro et al. [4] found that, for a nonlocal reaction-diffusion equation, localised spikes can be either purely attractive or purely repulsive; in contrast, the interaction between localised spots alternates between attractive and repulsive depending the distance separating the patterns, allowing for the existence of bound states. Similar results can be found in [49], where localised spots with monotonic tails are seen to be repulsive while localised gaps with oscillatory tails were found to exhibit bound states. Throughout this work, we will attempt to a) b) Figure 2: a) For localised spikes, the top picture shows a pair of real eigenvalues bifurcating from the origin (blue circle with cross) onto the real line (red circles). This results in a localised solution with monotonic tails, as shown in the bottom picture. b) For localised spots, the top picture shows four purely imaginary eigenvalues bifurcating from ±ik (blue circle with cross), for some k ∈ R, into the complex plane (red circles). This results in a localised solution with spatially-oscillatory tails with wave number k, as shown in the bottom picture.
distinguish between localised spikes and localised spots, identifying their bifurcation points and regions of existence for the von Hardenberg model. However, we first require some background theory regarding the existence of localised radial patterns, which we will now discuss in the context of some prototypical models.
Localised Radial Patterns One of the biggest issues in the analytical study of localised planar patterns is when localisation occurs in multiple directions. There has been some progress regarding two-dimensional patterns, such as [7,25,45], but the majority of analytical work focuses on localised radial patterns. In this case the localisation is restricted to just the radial direction, effectively reducing the problem to one independent variables. Localised radial patterns are seen as the first step to better understanding localised planar patterns [29]; however, radial problems tend to have explicit dependence on the spatial coordinate r, which can cause significant problems in the analysis.
There are a number of papers that investigate the existence of localised radial patterns in various systems, such as [26,32,35,37,54]. However, all of these examples utilise the radial invariant manifold and normal form theory introduced by Scheel [46], and so we quickly summarise the relevant results from [46] and their applications in two prototypical examples.
General results for localised radial patterns: We begin with a general stationary axisymmetric reaction-diffusion system, where D is a diffusion matrix, ∆ r := ∂ rr + 1 r ∂ r is the radial Laplacian operator, U ∈ R 2 is the dependent variable, µ ∈ R p is a parameter, and F is a nonlinear function where F(0; 0) = 0. We can also express this as a first-order ODE system where 1 2 denotes the 2-dimensional identity matrix. The key idea is the following; a localised radial pattern remains bounded close to the origin, which we will call the 'core' region, while exhibiting exponential decay sufficiently far from the origin, which we will call the 'far-field' region. The existence proof follows from a singular perturbation argument, and involves local invariant manifold theory. We first construct the set of all possible small-amplitude solutions that remain bounded in the core region. We call this set of solutions the core manifold, denoted by W cu − (µ), which we construct over the finite sub-interval r ∈ [0, r 0 ]. Likewise, we then construct the set of all small-amplitude solutions that exhibit exponential decay in the far-field region. We call this set of solutions the far-field manifold, denoted by W s + (µ), which is defined for sufficiently large values of r. We note that we are borrowing the notation for standard 'centre-unstable' and 'stable' manifolds to denote the core and far-field manifolds, respectively; this is equivalent to the notation used in previous works regarding localised radial patterns [32,35,37,46]. Once both of these manifolds have been constructed, we evaluate their parametrisations at some coincidental point, which we can choose to be at r = r 0 . Any function that lies on the intersection of the two manifolds is then, by definition, a localised solution.
Due to the explicit r-dependence of the system (1.1), the core manifold W cu − (µ) must be constructed explicitly; this is usually done by introducing a variation-of-constants formula and applying a fixed-point argument. In contrast, because the radial coordinate r is taken to be very large in the far-field region, the construction of the far-field manifold W s + (µ) can be seen as an extension of standard local centre-manifold theory. In order to parametrise the far-field manifold W s + (µ), it is helpful to reduce the problem onto the coordinates of a local centre manifold. In [46], it was proven that the extended vector field of (1.1), where r is sufficiently large, can be reduced to the centre eigenspace; then, the reduced equations have a leading order expansion that is partially determined by the type of bifurcation occurring as r → ∞. The two respective bifurcation types, homogeneous and Turing, are examined in [46]; their reduced equations are found in [46,Lemmas 3.9 & 3.10]. We now present one example of each bifurcation type, starting with a summary of the results of McQuighan and Sandstede [37] regarding a homogeneous bifurcation.

Example 1: Homogeneous Bifurcation
In [37], McQuighan and Sandstede investigate the existence of stationary localised radial patterns in the complex Ginzburg-Landau equation with 2 : 1 forcing. By choosing the parameters of the system such that a homogeneous bifurcation occurs, the problem can be reduced to the real reaction-diffusion system, where U ∈ R 2 , ε is the bifurcation parameter, and C i is a square matrix that is determined by the parameters of the original system, for each i = 1, 2, 3. This can also be written as the following first-order ODE system, The linear solutions of (1.2) are associated with the eigenvectors of C 1 , which has eigenvalues λ = 0 (double), −m, and m for some m ∈ R. For ε > 0 the zero eigenvalues split onto the real line, as seen in Figure 2 a), and so we have a homogeneous bifurcation as stated. By determining four linearly independent solutions of the linearisation of (1.2), the core manifold W cu − (ε) was constructed using the standard approach of a fixed-point argument for a variation-of-constants formula; see [37,Lemma 3.2]. In order to construct the far-field manifold W s + (ε) and reduce the equations to some local centre manifold in the far-field region, McQuighan and Sandstede utilised the theory of strong stable foliations. This is presented in great detail in Section 4 of [37], and so we just note the following: there exists coordinates A, B ∈ R on some autonomous local centre manifold W c + (ε) such that (1.2) reduces to where c 0 3 is a fixed constant and R A , R B contain higher-order terms. Here, σ ∝ 1 r is taken to be an unconstrained variable such that the system becomes autonomous; in order to recover the r-dependent system (1.2), we just define σ(r 0 ) = 1 r0 when evaluating the far-field parametrisation. To find exponentially decaying solutions to (1.3), McQuighan and Sandstede introduced geometric blowup coordinates to track solutions backwards in r to r = r 0 . They determined that exponentially decaying solutions to (1.3) were associated with exponentially decaying solutions to the cubic non-autonomous Ginzburg-Landau equation, where u ∈ R, ε is the bifurcation parameter, and ν is a fixed parameter. As with the previous examples, we express this equation as a first-order ODE system,  Examining the spatial eigenvalues at the bifurcation point as r → ∞, we observe a pair of purely imaginary eigenvalues at λ = ±i; these eigenvalues split off of the imaginary axis for ε > 0, and so we have a Turing bifurcation as illustrated in Figure 2 b). The core manifold W cu − (µ) is constructed in the usual way; see [32, Lemma 1]. The system (1.5) is already reduced to the centre eigenspace, and so no extra local centre manifold reduction is required. Introducing some complex coordinates A, B ∈ C, (1.5) can be written as where c 0 , c 3 are fixed constants and R A , R B contain higher-order terms. Again, σ ∝ 1 r is taken to be an unconstrained variable such that the system becomes autonomous; in order to recover the r-dependent system (1.5), we just define σ(r 0 ) = 1 r0 when evaluating the far-field parametrisation. To find exponentially decaying solutions to (1.6), geometric blow-up coordinates are introduced to allow for tracking solutions backwards in r to r = r 0 . Exponentially decaying solutions to (1.6) were found to be associated with exponentially decaying solutions to the cubic non-autonomous Ginzburg-Landau equation, For c 0 > 0, c 3 < 0, there exists a exponentially decaying solution to (1.7); there also exists a sufficiently small solution that decays exponentially independent of the sign of c 3 . Tracking back through the blow-up coordinate charts, Lloyd and Sandstede [32] proved the existence of an elevated localised spot in (1.5), as well as localised rings. Following this, McCalla and Sandstede [35] also proved the existence of a localised depressed spot, termed spot B. We will use these examples as a guide for our analysis of localised radial vegetation patterns.

Results
We summarise our results in Figure 3. We look for two types of localised solutions bifurcating from the homogeneous instability: elevated spikes and depressed spikes, which we now refer to as 'spikes' and 'gaps', respectively. We find that no localised spikes bifurcate from the homogeneous instability of the von Hardenberg model, while localised gaps bifurcate both sub-and super-critically. However, for standard parameter values, localised gaps are seen to have a negative profile; hence, these solutions are unphysical.
Close to a Turing instability, we prove that three classes of localised patterns exist and bifurcate subcritically from a uniform vegetated state. These solutions include spots A and B, with respective maximum and minimum values at the core r = 0, as well as elevated and depressed rings, whose maximum values are away from the core. We expect each of the localised patterns bifurcating from the Turing instability to remain physical for a reasonable choice of parameter values.
Overview The rest of the paper is structured as follows: we begin by introducing the von Hardenberg model [55] for stationary radial patterns, identifying uniform states and their associated bifurcation diagram. We then introduce the weakly nonlinear approximation for the von Hardenberg model derived in [13]; again, we compute uniform states and the bifurcation diagram of the system. We identify three parameter regions of the precipitation p where localised states may emerge; these bifurcations are categorised by their spatial eigenvalues, and so we label each as either 'homogeneous' or 'Turing' bifurcations, appropriately.
Having categorised the bifurcations of the system, we relate each case to the relevant prototypical examples ((1.2) or (1.5)) and apply weakly nonlinear theory to compute the relevant coefficients required in the analysis. Where possible, we quote results from the prototypical examples; if the results are not identical we highlight the extra work necessary to complete the proof, although the technical details are beyond the scope of this paper.
Following this, we arrive at the following results for the reduced model: we prove the existence of three classes of localised radial patterns bifurcating from the uniform vegetated state at some Turing point p = p 1 , where we recall that p denotes precipitation. These localised patterns have spatially-oscillatory tails with an associated wave number k 1 > 0, and are directly associated with the spot A, spot B and ring solutions found in the Swift-Hohenberg equation [32,35]. We also find evidence of the existence of localised gaps, characterised by depressions at their core and monotonic tails, bifurcating from the critical precipitation p = p c where the bare state destabilises. We expect there to be one solution bifurcating sub-critically from the bare state and another bifurcating super-critically from the uniform vegetated state, though these solutions may be unphysical.
We conclude by finding numerical solutions to the full von Hardenberg model corresponding to the localised spot A and localised gap solutions found in our analysis. We explore the effect of changing the shading parameter ρ and plot the various bifurcation curves of these localised radial structures for distinct parameter choices.

The Model
We focus our analysis on the model introduced by von Hardenberg et al. [55] to explain vegetation dynamics via a two-component reaction-diffusion system. This can be seen as an adaptation of the model presented in [28]; we note that there have been numerous studies of the von Hardenberg model (see [39] for an in-depth description of model terms and behaviour, and [24] for an analysis of two-dimensional pattern selection), including the work of Dawes and Williams [13] in which localised one-dimensional patterns are observed. Hence, this model is a natural first choice for the study of localised radial patterns, since we can then directly compare our results with the one-dimensional case. We consider two non-negative dynamical variables n(x, t) and w(x, t), denoting the plant biomass density and ground water density, respectively; then, the dimensionless model takes the form, Here, ∆ is the planar Laplacian, p represents precipitation, and δ denotes the relative diffusivity of the water with respect to the vegetation, which is taken to be a large, positive parameter. In addition, γ, σ, ρ, β, and ν are taken to be O(1), positive parameters. We briefly summarise the ecological role of each term in (2.1). The term ( γw 1+σw )n models plant growth which has a linear rate for low levels of water. The terms −νn and −n 2 account for the plant mortality and saturation due to limited nutrients, respectively, while the term ∆n models the spread of vegetation. The term p denotes the increase of water due to precipitation, while the −(1 − ρn)w term represents a loss of water due to evaporation. Hence, we require n ≤ 1 ρ in order for the model to remain physical, where the parameter ρ > 0 represents the reduction in evaporation caused by the presence of vegetation, mostly through shading effects. The term −w 2 n models the local uptake of water, mostly determined by transpiration. The ∆(w − βn) term represents the transport of water in the soil as modelled by Darcy's law, where the βn term accounts for the suction of water by the roots of the plant. The −c(w − αn) x term models the advection of soil water down a slope in the x-direction, where c is the runoff velocity. We will restrict our model to level ground, setting c = 0, in order to study localised radial patterns. It would be interesting to reintroduce a slope, where c > 0, and investigate how localised radial patterns are affected by this imposed symmetry-breaking; this is left for future study. Setting c = 0 and considering stationary axisymmetric solutions, (2.1) becomes where r is the planar radial coordinate and ∆ r := ∂ rr + 1 r ∂ r is the planar axisymmetric Laplacian. We begin by discussing the uniform states associated with (2.2).
Uniform States Any spatially uniform solution of (2.2), which we call (n * , w * ), must solve the following algebraic system, Hence, we see that there are two types of spatially uniform equilibria: the 'bare' state (n * , w * ) = (0, p), and 'vegetated' states where n * = γw * 1+σw * − ν and w * satisfies The bare state is known to undergo a homogeneous instability at p = p c := ν γ−νσ (as discussed in [13,24,39]) and so we require that γ > νσ in order for p c to be physically relevant. We note that this condition also guarantees that (2.3) has at least one positive solution w * , by Descartes' rule of signs; for our choice of parameters, as seen in [13,24], there is only one uniform vegetated state that bifurcates from the Bare state at p = p c . Turing bifurcations have been found at some values p = p 1 , p 2 , with associated critical wave numbers k = k 1 , k 2 , where periodic patterns emerge from the vegetated state; see [24]. For a bifurcation diagram of the uniform solutions to (2.2), see Figure 4

a).
Weakly Nonlinear Reduction In order to make (2.2) more analytically tractable, we now simplify the model in the case where δ ≫ 1; to do this, we first introduce ε := δ − 1 2 such that we can write (2.2) as This is now reminiscent of problems studied in [3,10,27], in which existence of solutions can be proven via geometric singular perturbation methods. Similar techniques should be applicable here, although the r-dependent terms in (2.4) may cause some extra complications. Instead, we will perform a weakly nonlinear reduction seen in [13] in order to find a simpler radial reaction-diffusion system that we can analyse. Following this, we will then numerically investigate the full system (2.2) for physically relevant parameter values, using the intuition gained from the reduced model.
In [13], Dawes and Williams presented a simplified model of (2.2) for large values of δ; for ε = δ − 1 2 , the authors found that where (n * , w * ) denotes the uniform vegetated state, p 1 is the minimum Turing point at which the vegetated state destabilises, and k 1 is the critical wave number associated with p 1 . Then, by defining p = p c + ε 2 P , n = ε 2 N , w = p + ε 2 W and r = εr, (2.4) can be rewritten, to leading order, as For simplicity, we perform the transformation (N, W, P ) → N, βW, b −1 0 P and drop the tilde notation for r, such that we can write (2.5) as where a = a0 β and b = b 0 β. Again, we see that there are two uniform states: the 'bare' state B : (N, W ) = (0, 0), and a uniform vegetated state V : Then, we require sgn(P ) = sgn(1 − ab) such that N 0 > 0. We note that 1 > ab for the choice of parameters seen in [13,24]; see (4.1), (4.2). Hence, we see that the vegetated state is bifurcating from the bare state super-critically; the bifurcation diagram for the uniform states of (2.6) is presented in Figure 4 b).

Localised Radial Patterns in the Reduced Model
We now wish to identify possible bifurcations of the reduced system (2.6), where localised solutions may emerge as perturbations from a uniform state. To do this, we will linearise about each uniform state in turn, investigating the spatial eigenvalues in the limit as r → ∞.

Bifurcation Analysis and Canonical Forms
We recall that there are two uniform states associated with (2.6): the bare state B : (N, W ) = (0, 0), and the vegetated state V : . We begin by looking for any localised radial solutions bifurcating from the bare state B.
Case 1 We linearise about the bare state B; then, we can express (2.6) as the first-order system We investigate the eigenvalues of the linear operator L B (r) in the limit as r → ∞; any eigenvalue λ ∈ C of lim r→∞ L B (r) must be a root of the characteristic polynomial, We see that solutions of p B (λ) = 0 take the form λ = 1, −1, and ± √ −P . Hence, for P > 0 we have purely imaginary eigenvalues ±i √ P ; these eigenvalues pass through the origin at P = 0 and, for P < 0, lie on the real axis at ± |P |. Thus, we define P = −ε 2 such that the bifurcation point is at ε = 0, and we expect localised solutions to emerge in the region 0 < ε ≪ 1. Taking ε small, we now write (3.1) as This is the only bifurcation for the bare state B, and so we now turn our attention to localised solutions bifurcating off of the vegetated state V. (2.7) and N 1 , W 1 are taken to be small. Then, we can write (2.6) as the first-order system, We investigate the eigenvalues of the linear operator L V (r) in the limit as r → ∞; an eigenvalue λ ∈ C of lim r→∞ L V (r) must be a root of the characteristic polynomial, where we have defined K : For small values of P , we see that Hence, taking P = ε 2 , we again have a bifurcation at ε = 0 such that we expect localised solutions to emerge in the region 0 < ε ≪ 1. We can now write the roots of (3.4) as and express (3.3) as where L 1 (r) is the linear operator defined in (3.2), and Thus, at the bifurcation point P = 0, we have a homogeneous bifurcation close to which localised radial solutions may emerge, bifurcating from either uniform state. Notably, any solutions bifurcating from the bare state B will emerge in a sub-critical region (P < 0), whereas any solutions bifurcating from the vegetated state V will emerge in a super-critical region (P > 0). As previously discussed, we also expect the vegetated state V to undergo a Turing instability, and so we return to our characteristic polynomial p V (λ).
Case 3 For the vegetated state V, we expect a Turing instability to be associated with the repeated roots of p V (λ), as defined in (3.4). Completing the square, we note that and so p V (λ) has repeated roots at some critical P = P c , where P c satisfies Both choices of P c,± have an associated eigenvalue λ ± , where For the parameter choices seen in [13,24], we find that K < 0; hence, λ 2 + < 0 and λ 2 − > 0, and so a Turing instability occurs at P = P 0 := 2−K+2 The critical eigenvalue is defined as λ = ±ik, where We note that P 0 = k 4 and take P = k 4 + P ; then, we see that We recall that the Turing instability occurs as the spatial eigenvalues bifurcate off of the imaginary axis, as seen in Figure 2 b). This would mean that λ 2 has non-zero imaginary part, and so P must be negative; that is, we expect the Turing instability to be a sub-critical bifurcation. We define P := P 0 − ε 2 such that ε = 0 is the bifurcation point and we expect localised solutions to emerge for 0 < ε ≪ 1. Hence, we can now write the eigenvalues {λ i } 4 i=1 as and express (3.3) as where We introduce the following transformation which has a well-defined inverse .
Then, defining V := PV, r := kr and ε := k −2 ε, we can finally write (3.6) as with the fixed parameter ω := k 2 +1 k 2 . Finally, we drop the tilde notation for simplicity. Thus, we have three regions of interest: the sub-and super-critical regions close to a homogeneous bifurcation at P = 0, where localised solutions may bifurcate from the respective bare or vegetated states, and the sub-critical region close to a Turing bifurcation at P = P 0 , where localised solutions may bifurcate from the vegetated state; see Figure 5. We now look to relate these cases to the prototypical equations seen in [32,35,37] associated with each bifurcation type, and then discuss the existence of possible localised radial solutions in each case.

Cases 1 and 2: Homogeneous Bifurcation
We recall that cases 1 and 2 both occur close to a homogeneous bifurcation, and so we expect that we can relate both to the study of oscillons in the complex Ginzburg-Landau equation [37]. To make this more clear, we express our equations (3.2), (3.5) in a more general form; that is, for cases 1 and 2, we write our equations in the following way, for some U, V ∈ R 2 , C, D ∈ R 2×2 , and f : R 2 → R 2 . In both cases, however U, V and D differ for each case. In particular, for case 1, whereas in case 2, we have By writing our equations in this form, we have expressed cases 1 and 2 in an analogous way to the CGL problem in [37]. We note that the general system (3.9) has a quadratic nonlinearity, which is in contrast to Z 2 symmetry seen in [37]. However, the core parametrisation and the foliations in the far-field are entirely dependent on the linear part of the system, and so the results in [37] should follow identically if C has the same properties as the matrix C 1 in [37]. Hence, we begin by analysing the associated linear problem of (3.9).
We note that, at the bifurcation point ε = 0, the linear operator of (3.9) in the limit as r → ∞ has eigenvalues λ c = 0 (with double multiplicity), λ u = 1, and λ s = −1; these eigenvalues have associated generalised eigenspaces where, We note that this has the same structure as in [37] (in the case when m = 1) for different values of U 0 , U 1 . Then, the linearisation of (3.9) at ε = 0, for some U, V ∈ R 2 , has four linearly independent solutions {W i (r)} 4 i=1 , given by where I ν (r), K ν (r) are the ν th -order modified Bessel functions of the first and second kind, respectively; see [1, Section 9.6]. We also note that span . We define P cu − (r) to be the projection on span {W 1 (r), W 3 (r)} along span {W 2 (r), W 4 (r)} for each fixed r ≥ 0.
These solutions have exactly the same structure as [37]; since the explicit values of U 0 and U 1 are not used in the following analysis, we can conclude the following results from [37].
The first result is regarding the existence and parametrisation of the core manifold W cu − (ε), which contains all small-amplitude solutions of (3.9) on the interval r ∈ [0, r 0 ], for some large fixed r 0 > 0, that remain bounded as r → 0, Lemma 3.1 Fix r 0 > 0; then, there are constants δ 0 , δ 1 , ε 0 > 0 such that the set W cu − (ε) of solutions U(r) of (3.9) for which sup 0≤r≤r0 |U| < δ 0 is a smooth two-dimensional manifold for each ε < ε 0 . Furthermore, Proof. This follows from Lemma 3.2 and Remark 3.1 in [37] and involves a standard application of the variation-of-constants formula on a bounded interval r ∈ [0, r 0 ]; see [32] for a similar example.
Next, we establish the existence and parametrisation of the far-field manifold W s + (ε), which contains all small-amplitude solutions of (3.9), for sufficiently large r, that decay exponentially as r → ∞. To do this, we replace any 1 r terms by the unconstrained variable σ, defined such that (3.9) becomes an extended autonomous system of the form,  (3.14) In order to recover the non-autonomous system (3.9), we would later implement the restriction σ(r 0 ) = r −1 0 , thus ensuring σ(r) = r −1 for all r ≥ r 0 . The parametrisation of W s + (ε) relies on the existence of autonomous local centre-stable and centre manifolds for (3.14), as well as the existence of strong stable foliations; this is tackled in great detail in Section 4 of [37], which we summarise in the following result, Lemma 3.2 There exist constants ε 0 , δ 0 , δ 1 , ρ 0 > 0 such that the set W cs + (ε) of solutions (U (r), σ(r)) of (3.14) for which sup r≥ρ0 {|U (r)| , |σ(r)|} < δ 0 is a smooth four-dimensional manifold for each 0 < ε ≤ ε 0 . Furthermore, (U, σ) ∈ W cs + (ε) has the expansion  [37]. The parametrisation of the local centre-stable manifold is then derived in Section 4.4 of [37], resulting in (4.9).
We note that the functions A and B are defined in [37] as the projections of a point (p, σ) on a local centre manifold W c + (ε) onto W 0 1 and W 0 2 , respectively. In order to parametrise the far-field manifold W s + (ε), we need to determine an expression for (U, σ) ∈ W cs + (ε) such that |U (r)| decays exponentially as r → ∞. Assuming that we can determine the conditions such that A and B are exponentially decaying, and so |U (r)| is also exponentially decaying in (3.15), we then evaluate our parametrisation at σ = 1 r in order to recover solutions U (r) ∈ W s + (ε) for the non-autonomous system (3.9). Following the parametrisation of the far-field manifold W s + (ε) by exponentially-decaying coordinates A, B on some local centre manifold W c + (ε), we now reduce (3.14) to a system of equations for the centre coordinates (A, B, σ), Lemma 3.3 Using the coordinates P c + U (r) = A(r)W 0 1 +B(r)W 0 2 , the vector field (3.14) restricted to W c + (ε) can be written as where c 0 = 1+a 1+a 2 and c 2 = (1 − ab). The remainder terms R A and R B satisfy Proof. This follows identically to Lemma 4.6 in [37] except that we no longer have Z 2 -symmetry, due to the quadratic nonlinearity in (3.9). Instead, we retain the more general normal form seen in [46, (58) on page 42], where [46, Lemma 3.9] implies that for some fixed constants c 0 , c 0 2 . In order to obtain the normal form (3.16), we introduce a new coordinate B such that B = B + R A (A, B, σ, ε). Then, because R A is higher order in B, no new terms are introduced to leading order and we obtain the leading-order expansion seen in (3.16). Finally, the leading order expansions of each coefficient can be computed via weakly nonlinear analysis seen in [9]; see Appendix A.1. We note that these coefficient values can also be found rigorously via normal form transformations, as seen in [46].
In Section 5 of [37], localised solutions were found to take the form A(r) = εA 2 (εr), where A 2 (s) satisfies where c 0 3 is the cubic coefficient in the radial normal form, such that the full system (1.2) has a localised solution if there exists a localised solution to the real non-autonomous Ginzburg-Landau equation (3.17). Similarly, we can take an ansatz of the form A(r) = ε 2 A 2 (εr), such that the quadratic nonlinearity scales at the same rate as the linear part of the normal form (3.16). Then, we see that localised solutions to (3.16) are associated with localised solutions to the real non-autonomous Ginzburg-Landau equation We note that (3.18) can be rescaled such that where we have defined s := √ c 0 s, A 2 := −c 0 2 A 2 . By [30], we know that there exists a unique positive localised solution to (3.19); we note that sgn (A) = −sgn (c 0 2 ), and so we expect to find localised depressed spikes. However there remains some significant technical analysis required to complete this existence proof: including setting up appropriate geometric blow-up charts, proving the existence of a connecting orbit between the two charts related to localised solutions of (3.18), and tracking solutions through the various charts. Hence, the full existence proof is beyond the scope of this paper. We will, however, use this as motivation for our numerical investigation in Section 4, where we observe localised depressed spikes bifurcating from the homogeneous bifurcation point, both sub-and super-critically. In the sub-critical region the spikes decay to the bare state as r → ∞, whereas the spikes in the super-critical region decay to the uniform vegetated state as r → ∞. This is all in agreement with our analysis, and so we expect that these localised solutions exist and bifurcate from the homogeneous bifurcation point P = 0 as previously discussed.
For both cases 1 and 2, we have followed the work of [37] to produce a (currently incomplete) proof of the existence of localised radial spikes bifurcating from the homogeneous bifurcation at P = 0. The general equation (3.9) is sufficiently similar to [37] that much of the analysis follows identically up to the centremanifold reduction in the far-field. We expect that the existence of localised solutions to (3.9) is dependent on the existence of localised solutions to the Ginzburg-Landau equation (3.18), although this analysis has been left for future work. Finally, we applied the result seen in [30] to conclude that, when the uniform vegetated state bifurcates super-critically at P = 0, localised depressed spikes may emerge in both cases 1 and 2. This is supported by our numerical results in which we find localised depressed spikes bifurcating both sub-and super-critically from the homogeneous instability for the full von Hardenberg model (2.1); see Section 4.
We note that in case 1 solutions are bifurcating from the bare state B, and so localised depressed spikes are always unphysical in our vegetation model. In case 2 solutions are bifurcating from the uniform vegetated state V, and so it is possible that these solutions may be physical if the amplitude of the spike is sufficiently small such that N (r) and W (r) remain non-negative. Our results also suggest that models where the uniform vegetated state bifurcates sub-critically, such as [22,40], may allow for localised elevated spikes to emerge close to the critical precipitation p = p c . However, a closer analysis of the model would first be required in order to substantiate such a prediction.

Case 3: Turing Bifurcation
We recall that case 3 occurs close to a Turing bifurcation, and so we expect that we can relate this case to the study of localised radial patterns in the Swift-Hohenberg equation [32]. In fact, the system (3.8) has been transformed in such a way that the linear part is identical to the Swift-Hohenberg equation seen in [32], and so the results of [32,35] will follow almost identically. We first note that the linearisation of (3.8) at the bifurcation point ε = 0 has four linearly independent solutions, as seen in [32, (2.4)], of the form where J ν , Y ν denote the ν th Bessel functions of the first and second kind, respectively. We denote the projection onto the space span{V 1 (r), V 2 (r)} as P cu − (r), with null space given by span{V 3 (r), V 4 (r)}, for each fixed r ≥ 0. We now conclude the following results from [32,35]. The first result is regarding the existence and parametrisation of the core manifold W cu − (ε), which contains all small-amplitude solutions of (3.8) on the interval r ∈ [0, r 0 ], for some large fixed r 0 > 0, that remain bounded as r → 0, Lemma 3.4 Fix r 0 > 0; then, there are constants δ 0 , δ 1 , ε 0 > 0 such that the set W cu − (ε) of solutions U(r) of (3.8) for which sup 0≤r≤r0 |U| < δ 0 is a smooth two-dimensional manifold for each ε < ε 0 . Furthermore,
Proof. This follows from [32, Lemma 1] and involves a standard application of the variation-of-constants formula on a bounded interval r ∈ [0, r 0 ]. As seen in [32], the quadratic coefficient ν is found from a Taylor expansion to be where [X] 1 denotes the first element of a vector X. Hence, we see that where we have used [56, (3) on page 411] to compute the final integral.
Next, we establish the existence and parametrisation of the far-field manifold W s + (ε), which contains all small-amplitude solutions of (3.8), for sufficiently large r, that decay exponentially as r → ∞. To do this, we replace any 1 r terms by the unconstrained variable σ, defined such that (3.8) becomes an extended autonomous system of the form,  In order to recover the non-autonomous system (3.8), we would later implement the restriction σ(r 0 ) = r −1 0 , thus ensuring σ(r) = r −1 for all r ≥ r 0 . The far-field normal form equations are constructed identically to [32]; we introduce the coordinates A, B ∈ C, defined in [32, (3.5

)] as
and present the following result, Lemma 3.5 Fix 0 < m < ∞; then, there exists a change of coordinates

36
, and the remainder terms satisfy Proof. This follows identically from [32, Lemma 2], where we compute the leading order expansions of c 0 , c 3 via weakly nonlinear analysis seen in [9]; see Appendix A.2. Again, we note that these coefficient values can also be found rigorously via normal form transformations, as seen in [46].
Localised solutions of (3.25) are found to take the form A(r) = ε It was determined in [32,35] that there are two types of localised solutions to (3.26) for c 0 > 0: one that satisfies the full nonlinear equation (3.26), and another where |A 2 | ≪ 1 and so solutions follow the linear flow of (3.26). The linear flow of (3.26) has a unique localised solution that can be found explicitly and is independent of c 3 ; it was proven in [52] that the cubic equation , where q 0 > 0 and q + = 0 are fixed constants. By introducing geometric blow-up coordinates, as seen in [35], one can track small-amplitude exponentially decaying solutions back through extended phase spaces by following close the localised solutions of (3.26). This is entirely identical to the work in [32,35], which allows us to state the following theorems: uniformly on bounded intervals [0, r 0 ] as ε → 0, where J 0 (r) is a Bessel function of the first kind of order 0.
Each of these theorems follow identically from the work of [32,35], where the corresponding parameter values are computed in Lemmas 3.4 and 3.5; see Figure 6 for the profile of each solution. We note that ω := k 2 +1 k 2 > 1; hence, spot A solutions exist for all choices of ω, whereas spot B and ring solutions require that ω > ω * ≈ 1.3 in order to emerge. water and vegetation densities (N, W )(r). Then, in the core region 0 ≤ r ≤ k −1 r 0 , we find that Hence, we have found three classes of localised radial patterns bifurcating from the uniform vegetated state at the Turing point P = P 0 . These are all physical solutions for ε sufficiently small and have oscillatory tails with an associated wave number k := 1+ If these solutions persist away from the Turing point, these oscillatory tails may be damped by the stronger localisation, resulting in solutions that appear to have monotonic tails even though they bifurcate from a Turing instability.

Numerical Results for the von Hardenberg Model
We numerically investigate the existence of localised spot A and spike solutions in the full von Hardenberg model (2.2). We employ finite-difference and numerical continuation methods in order to find localised solutions to (2.2) and determine their bifurcation structure. While we expect to be able to also find localised rings and spot B solutions, as proven to exist in the reduced model in Theorems 2 & 3, these solutions typically exist in more restrictive parameter regions than spot A solutions and are harder to find numerically. Hence, while we do expect to be able to find such patterns, we restrict our attention to localised spike and spot A solutions. We also investigate the role of the shading parameter ρ in the emergence of these localised radial patterns and their bifurcation structure. In [24], the default choice of parameters is the following, γ = 1.6, σ = 1.6, ν = 0.2, β = 3, δ = 30, ρ = 1.5. (4.1) In this parameter regime, peaks arranged in a hexagonal lattice bifurcate sub-critically from p = p 1 whereas stripes bifurcate super-critically; see Figure 8 in [24]. In [13], localised one-dimensional solutions do not emerge for this parameter region; in order to numerically find such solutions, the authors introduce a second choice of parameters, {γ, σ, ν, β, δ} as defined in (4.1), ρ = 2.5, (4.2) which ensures that stripes bifurcate sub-critically, and so localised one-dimensional solutions emerge from the Turing point p = p 1 . We wish to see how the choice of ρ affects localised radial patterns; we fix the other parameters {γ, σ, ν, β, δ} as in (4.1) and present our results for four distinct cases: These cases include the choices of ρ seen in [24] and [13] (cases (i) and (iii), respectively), as well as an intermediate value (case (ii)) and an exceeding value (case (iv)) for ρ. We note that our analysis holds as long as a < 1 b , where which, for {γ, σ, ν, β, δ} fixed as in (4.1), is equivalent to the condition ρ < (1+σpc) 2 γpc + p c ≈ 6.41. Thus, each case in (4.3) is valid for our weakly nonlinear analysis.
Numerical Implementation We briefly describe our numerical scheme for finding localised solutions to (2.2) via central finite-difference methods. We recall that the von Hardenberg model (2.2) can be written in the form where U := (n, w) T . We define our finite radial domain 0 ≤ r ≤ r * := 400 and impose the following Neumann boundary conditions lim r→r * ∂ r U = 0, and lim We are looking for solutions that decay exponentially to some uniform state as the radius r increases, which is why we use homogeneous Neumann boundary conditions at the outer boundary; the inner boundary condition is defined such that (4.4) is well-defined in the limit as r → 0. We discretise our radial domain into T := 2000 grid points {r i } T i=1 and define U := (n, w) ∈ R 2T , where n i = n(r i ) and w i = w(r i ) for each i ∈ [1, T ], such that we can write (4.4) as Here, 1 T , 0 T and e T are the identity matrix, zero matrix and vector of ones for R T , respectively. Furthermore, the • notation indicates a Hadamard (or element-wise) product and D r := D 2 +RD 1 ∈ R 2T ×2T approximates the operator ∆ r . The matrices D 1 and D 2 are the respective central finite-difference matrices for first and second derivatives, R := diag(1, r −1 2 , r −1 3 , . . . , r −1 T ) approximates 1/r, and the first and last rows of each matrix are chosen such that the boundary conditions (4.5) are satisfied.
We use Matlab's trust-region solver fsolve to solve the discretised problem (4.6) for a localised initial guess, and then employ path-following codes from [2] in order to continue our solutions in parameter space. At each localised solution U * := (n * , w * ) of (4.6), with p = p * , we calculate the radial stability of our solution by computing the eigenvalues of the Jacobian of the discretised system evaluated at the localised state (n * , w * ). That is, we solve the eigenvalue problem λV = δ U F (U * ; p * )V for some V ∈ R 2T , where the Jacobian δ U F (U * ; p * ) is the following, If we find an eigenvalue λ of (4.7) such that Re(λ) > 0, we conclude that the solution U * is unstable with respect to radial perturbations; likewise, if every eigenvalue λ is such that Re(λ) < 0, we conclude that the solution U * is stable with respect to radial perturbations. We emphasise that radial stability is a necessary, but not sufficient, condition for solutions to be stable in the two-dimensional plane; for example, a solution might be stable with respect to radial perturbations but unstable with respect to perturbations in the x-direction. Hence, from these calculations we can only definitively say when solutions are unstable, not stable.
We note that the size of our radial domain r * must be chosen sufficiently large in order to find localised solutions. If the domain is too small, then solutions will solve the boundary value problem (4.4)-(4.5) but will not exhibit exponential decay as r increases. Our choice of r * = 400 is acceptable in our numerical study, since the observed solutions decay to their uniform state well before reaching the outer boundary. In fact, the solution profiles shown in Figures 7-10 are only displayed for r ∈ [0, 200]; each solution remains flat for the remainder of the radial domain, and so we omit this from the figures. Then, for each localised solution, the profile in the core remains the same for different domain sizes as long as the domain is sufficiently large such that boundary effects are avoided. We now present our numerical findings for cases (i)-(iv).
Case (i) The bifurcation diagram for ρ = 1.5 is plotted in Figure 7, where the L 2 -norm of the vegetation density n(r) is plotted against the shifted precipitation p − p c . We recall that, in [13], this choice of ρ does not exhibit one-dimensional localised patterns bifurcating from the Turing point p = p 1 . As predicted in our earlier analysis, however, we observe three localised radial patterns: localised gap solutions bifurcating suband super-critically from the homogeneous point p = p c , respectively, as well as a localised spot A solution bifurcating sub-critically from p = p 1 .
The localised spot A branch is unstable to radial perturbations as it bifurcates from the uniform vegetated state at p = p 1 ; as p decreases, the spot undergoes a fold bifurcation where it gains stability with respect to radial perturbations. The changing radial profile of the spot can be seen in panels a) and b) of Figure 7.
As p re-approaches the critical value p 1 , the localisation of the spot becomes weaker and its oscillatory tails become more prominent. Eventually the localisation becomes so weak that the oscillating tails reach the outer boundary, and the bifurcation curve of the localised radial spot connects to a domain-covering target pattern; the bifurcation curve for this domain-covering pattern is omitted from Figure 7.
The super-critical localised gap solution remains unstable with respect to radial perturbations for all p > p c . The solution emerges as a monotonic negative perturbation from the uniform vegetated state at p = p c , as seen in Figure 7 c), and increases in amplitude as p increases. The branch also undergoes a fold bifurcation, although the solution remains unstable to radial perturbations. For the upper part of the solution branch, the localised gap gains a positive peak, as seen in Figure 7 d). We note that this peak remains much smaller than the depression at r = 0, however. The sub-critical localised gap solution bifurcates from the bare state at p = p c and remains unstable with respect to radial perturbations. An example of this solution's radial profile is plotted in Figure 7 e). As mentioned previously, the sub-critical localised gap is a negative perturbation from the bare state, and so will always be unphysical in the context of our vegetation model. We also note that, for this choice of ρ as well as the remaining cases (ii)-(iv), the super-critical localised gap is a sufficiently large perturbation from the uniform vegetated state that it is also unphysical. As the precipitation increases, this perturbation appears to increase in amplitude faster than the uniform vegetated state, and so we expect the entire branch to be unphysical. Although the localised gap solutions exist in the mathematical model, they are both unphysical for each of our choices of ρ; hence, for the remaining cases, we will focus on the behaviour the localised spot A solution.
Case (ii) The bifurcation diagram for ρ = 2 is plotted in Figure 8, with the same measures as in case (i). We note that, in [13, Figure 15], this choice of ρ may exhibit one-dimensional localised patterns bifurcating from the Turing point p = p 1 , but only for a very small parameter region of p. As predicted in our earlier analysis, we observe sub-and super-critical localised gap solutions and a localised spot A solution; see Figure  8.
Similarly to case (i), the localised spot A branch is unstable as it bifurcates and then undergoes a fold bifurcation where it gains stability with respect to radial perturbations. The first difference to note between cases (i) and (ii) is the increased region of existence for the localised spot A solution. By increasing ρ from (i) 1.5 to (ii) 2, the length of the spot A branch, in terms of the precipitation p, has grown from (i) 0.01 to (ii) 0.032. Hence, there is a much larger range of precipitation values for which stable localised spot A solutions exist.
Close to the Turing point p = p 1 , the solution has spatially-oscillating tails with an associated wave number k 1 ; see Figure 8 a). However, as the solution branch moves further away from p = p 1 , the localisation becomes sufficiently strong that the solution appears to be monotonic; see Figure 8 b). We emphasise that oscillations with wavenumber k 1 are still present in these solutions, but they have become sufficiently small such that they are not observed qualitatively. This means that a localised spot, typically associated with a Turing bifurcation, can take on an appearance similar to that of a localised spike, typically associated with a homogeneous bifurcation. We note that, following the fold bifurcation, the spot branch returns to p = p 1 , now with a larger amplitude. The localisation of the solution becomes weaker as the precipitation approaches the critical value p = p 1 , and so spatially-oscillating tails with wave number k 1 become more prominent again; see Figure 8 c). As with case (i), the solution branch for the localised spot terminates when the localisation becomes sufficiently weak that its oscillations reach the outer boundary; here, the solution coincides with a radial target pattern whose bifurcation curve we omit from Figure 8.
Case (iii) The bifurcation diagram for ρ = 2.5 is plotted in Figure 9, where the L 2 -norm of the vegetation density n(r) is plotted against the shifted precipitation p − p c . We recall that this is the parameter choice in [13] such that localised one-dimensional solutions exist and engage in homoclinic snaking; see [13, Figure  12]. Again, we observe sub-and super-critical localised gap solutions and a localised spot A solution; see Figure 9.
As with the previous cases, the localised spot A branch is unstable as it bifurcates and then undergoes a fold bifurcation where it gains stability with respect to radial perturbations. Again, compared to the previous cases, we see that there is a larger range of precipitation values for which stable localised spot A solutions exist. A notable change for ρ = 2.5 is that the localised spot A solution now exists and is stable for some part of the region p < p c , where one might expect localised peaks to bifurcate from the homogeneous point p = p c as opposed to the Turing point p = p 1 . As discussed previously, the localisation of the solution is sufficiently strong in this region, such that the spot A may appear to be a monotonic spike even though it has an implicit wave number k 1 . This highlights a potential issue with analysing localised structures experimentally, as their qualitative behaviour may belie the intrinsic properties of their bifurcation.
The behaviour of the localised spot A solution, illustrated in Figure  localised spot A destabilises close to p = p 1 in the bifurcation diagram; see Figure 9. At this point, the core value decreases such that the maximum value of n(r) is no longer at r = 0; see Figure 9 c). This type of behaviour has been seen in [37, Figure 15], in which the radial solution grows outwards and leaves a plateau in the core region. It is unclear whether the same phenomenon would occur in this problem, however, and we leave this question for future study.
Case (iv) Finally, we plot the bifurcation diagram for ρ = 2.7 in Figure 10. We note that this choice of ρ exceeds those used in [13], but we would expect localised one-dimensional solutions to exist and snake within the entire region between the homogeneous point p c and the Turing point p 1 . As with the previous cases, we observe sub-and super-critical localised gap solutions and a localised spot A solution; see Figure 10.
We note that the existence region for stable localised spot A solutions is smaller than in case (iii). The stable part of the spot A branch is now almost entirely contained in the region p < p c , and so the stable localised spot only appears in the form of a localised monotonic spike; see Figure 10 b).
The behaviour of the localised spot A solution, illustrated in Figure 10 a), b) & c), is identical to case (iii) for most p sufficiently far from the Turing point p 1 . However, for ρ = 2.7 we note that instability point seen in case (iii) has moved further from the Turing point; see Figure 10. At this instability, the core value decreases in a similar way to Figure 9 c), except with damped spatially-oscillating tails such that they appear monotonic; see Figure 10 c). Following this unstable branch, the solution decreases in amplitude as it approaches the Turing point and spatially-oscillating tails become more prominent again; see Figure 10 d). This solution appears similar to a localised ring, suggesting that the spot A branch may connect to other localised radial patterns near the Turing instability, but we leave this question for future study.

Conclusion
In this paper, we have investigated the existence of stationary localised radial patterns in the von Hardenberg model for dryland vegetation [55]. We began by applying a weakly nonlinear reduction, introduced by [13], to the stationary axisymmetric von Hardenberg model on flat soil (2.2), resulting in the stationary radial reaction-diffusion system (2.6). By comparing this system to prototypical pattern-forming equations seen in [32,37], we were able to analyse the existence of localised radial patterns close to a homogeneous instability and a Turing instability, respectively.
In a sub-(super-) critical region close to the homogeneous instability, we presented evidence that localised gap solutions to (2.6) bifurcate from the uniform bare state B (vegetated state V). We also note that we expect these gap solutions to be unphysical, at least for some range of parameter values, due to the vegetation density becoming negative. The existence proof for these solutions is almost identical to the work presented in [37]; however, due to a difference in the nonlinearity, further analysis is required to track exponentially decaying solutions in the far-field. This analysis would require tracking solutions through the geometric blow-up charts presented in [37] for our nonlinear solutions, which is beyond the scope of this paper, and so the proof is left incomplete. Also, in a sub-critical region close to the Turing instability, we proved that three classes of localised radial patterns exist and solve (2.6). These patterns are small perturbations from the uniform vegetated state V, and are equivalent to the patterns found in the quadratic-cubic Swift-Hohenberg equation [32,35]; see Figure 6.
Following this, we used our analytical results as a guide to numerically find localised radial patterns in the full von Hardenberg model (2.2). For simplicity, we restricted our investigation to the existence and radial stability of localised spot A and gap solutions, which have been predicted by our weakly nonlinear analysis. We used numerical path following routines to continue our solutions to larger amplitudes, and investigated the effect of changing the parameter ρ, related to shading effects in (2.2), on the bifurcation structure of the system.
Through this work, we have found that no physically-relevant localised radial patterns bifurcate from the homogeneous instability, at least for a range of appropriate parameter values. Therefore, in this parameter region, any observed localised vegetation patches are predicted to be the result of the uniform vegetated state undergoing a Turing instability. As the precipitation decreases, the oscillatory tails of a localised spot become damped such that they may appear monotonic; in this way, a localised spot can have the appearance of a localised spike far from the Turing instability.
Localised Spots We now comment on how our results relate to previous literature regarding localised spots in vegetation models. We first note that the majority of studies of localised spots use models where the uniform vegetated state bifurcates subcritically, and are focused on multi-state solutions, i.e. where there is more than one spot [21,40,59]. However, despite these differences, we find that our results remain in qualitative agreement with the wider literature. Localised spots have been observed in the field (see [40,Fig. 3 & 7], [31, Fig. 1]) and in many cases appear to be have monotonic tails [21,31,40,59]. However, we note that these solutions also tend to exhibit bound states [31,40], which suggests that they are localised patterns emerging from a Turing bifurcation [4,49]. This is in harmony with our results, as we only find localised spots bifurcating from the Turing bifurcation and observe that a localised spot can appear monotonic sufficiently far from its associated Turing point; see Figure 8 b).
We also find that localised spots can exist for every choice of shading parameter ρ in (4.3) as well as for very low precipitation values, where uniform vegetated states are no longer viable, which is in agreement with results for the Gilad et al. model in [40]. This is in contrast to localised one-dimensional patterns [13], otherwise known as stripes, which do not emerge for weak shading effects (ρ < 2) and do not persist for low precipitation values (p < p c ). Furthermore, whereas stripes undergo homoclinic snaking, where additional peaks are added to the solution as p oscillates [13], radial spots do not appear to undergo any snaking behaviour and instead terminate at a periodic-pattern branch, as seen in [59, Fig. 7] and Figures 7-10.
Another phenomenon observed in the literature is the radial destabilisation of localised spots, resulting in a transition to ring-like structures seen in [40]. For ρ = 2.5, 2.7, we observe a radial instability for the localised spot, as seen in Figures 9 & 10, such that the solution transitions towards a more ring-like structure as the precipitation increases. This is consistent with the behaviour described in [40, Fig. 10], but we leave further study for a follow-up investigation.  Fig. 3], which supports our hypothesis that fairy circles are related to the localised spot B solution.
Finally, it may also be possible that the von Hardenberg model is insufficient to find fairy circles, and so a different model may be required. There have been various models used to better study fairy circles; these include the Gilad et al. model [22] which, under certain assumptions, can be reduced to a reaction-diffusion system similar to the model studied in this paper [20,59]. There exist similar reaction-diffusion models for studying fairy circles, as seen in [60], as well as more complex equations such as the non-local models seen in [4,18,49]. Some of these models may prove difficult to study analytically, especially those with non-local terms, but any progress could be a major step in developing our understanding of the emergence of fairy circles.
Ecological Implications We will briefly discuss some of the ecological implications of this work. We observe that localised spots emerge from the uniform vegetated state V at a Turing point p = p 1 ; whereas domain-filling patterns mostly exist in the super-critical region p > p 1 , apart from a possible hysteresis region, localised radial patterns exist in the sub-critical region p < p 1 . Hence, we expect localised patterns to emerge when the precipitation levels are not sufficient to sustain the domain-covering patterns seen in [24], for example. Furthermore, we expect that localised radial patterns can be found for a much wider range of environments than localised stripes, which are associated with the one-dimensional patterns found in [13].
Localised stripes are not found for sufficiently weak shading effects (i.e. for ρ < 2), suggesting that some species of vegetation may only be able to maintain localised spots and not localised stripes.
For vegetation with weaker shading effects, for example cases (i) and (ii) in (4.3), localised spots appear as a local concentration of vegetation surrounded by a sparser vegetated state. An example of this can be seen in French Guiana, where isolated patches of woody vegetation have been found surrounded by grassland; see [31, Fig. 1]. For vegetation with stronger shading effects, for example cases (iii) and (iv) in (4.3), we instead find that localised spots persist as the precipitation decreases until they appear as isolated patches of vegetation surrounded by bare soil; examples of this have been found in the Negev region of Israel, such as in [40,Fig. 3 & 7]. Hence, changes to the bifurcation structure caused by the value of the shading parameter ρ can result in the formation of qualitatively different localised spot solutions in the field.
The saddle node of spot A solutions, as seen in Figures 7-10, defines a tipping point in climate dynamics; when the precipitation decreases beyond this point, solutions can undergo a transition to the stable uniform state. This transition is hysteretic because the uniform state remains stable as the precipitation returns to its original value. For plant species with sufficiently strong shading effects the spot A saddle node point extends into a region of low precipitation such that the uniform vegetated state is no longer viable; then, this tipping point represents an ecological collapse to the bare state B. This transition is sometimes referred to as a 'catastrophic regime shift' [42], and provides an example of desertification induced by a change in environment.
Further Study A possible extension of this work would be to use the singular perturbation methods seen in [10,27] to attempt to prove the existence of stationary localised radial patterns in the full von Hardenberg model (2.2). There is also an interesting question regarding the reintroduction of a sloping term into the original model. Instead of setting c = 0 in (2.1), what if we set c =ĉδ 1 2 such that there is some small symmetry-breaking that favours stripes over spots? Are we still able to find localised radial solutions for sufficiently smallĉ in the reduced model (2.6) and, if so, can we say anything about their stability? There is also the question of whether we can study the existence of more complicated localised patterns, such as solutions with dihedral symmetry [25] or labyrinthine patterns [12], in the von Hardenberg model. These are just a couple possibilities for future work; there has also been recent progress on centre manifold reduction techniques for nonlocal equations [17], which we may be able to adapt for radial solutions to nonlocal vegetation models seen in [40,50], for example. weakly nonlinear analysis, as seen in [8], in order to determine certain coefficients for each of the normal forms. In particular, for each case we will derive a nonlinear non-autonomous Ginzburg-Landau equation that determines the behaviour of small-amplitude solutions for large values of r.

A.1 Cases 1 and 2
For cases 1 and 2, (A.1) has the following form, where, respectively. From the normal form for this bifurcation, we expect to find a quadratic Ginzburg-Landau equation that defines the behaviour of solutions in the far-field, after some rescaling. Hence, we define R = εr and Then, we can write L = −C + ε 2 L 1 , N (u, v) = ε 4 N 1 (u 1 , v 1 ) + O(ε 6 ), where L 1 = ∂ RR + 1 R ∂ R 1 2 − D i , and so (A.1) can be expressed as The first equation implies that (u 1 , v 1 ) ⊺ must take the form of an eigenvector of C with eigenvalue 0; hence, we define (u 1 , v 1 ) ⊺ := A(R) U 0 , where U 0 is defined in (3.10). Then, we choose (u 2 , v 2 ) ⊺ = B(R) U 0 , such that the second equation becomes Then, for both cases 1 and 2, pre-multiplying by P results in the following equation for A(R), Hence, we have found that the coefficients c 0 and c 2 in (3.16) are defined such that c 0 = 1+a 1+a 2 and c 2 = (1 − ab).

A.2 Case 3
For case 3, we write the explicit form of (A.1) as From the normal form for this bifurcation, we expect to find a cubic Ginzburg-Landau equation that defines the far-field behaviour, and so we define R = εr and u(r) v(r) = ε u 1 (r, R) v 1 (r, R) + ε 2 u 2 (r, R) v 2 (r, R) + ε 3 u 3 (r, R) v 3 (r, R) + . . .