Predicting the emergence of localised dihedral patterns in models for dryland vegetation

Localised patterns are often observed in models for dryland vegetation, both as peaks of vegetation in a desert state and as gaps within a vegetated state, known as `fairy circles'. Recent results from radial spatial dynamics show that approximations of localised patterns with dihedral symmetry emerge from a Turing instability in general reaction--diffusion systems, which we apply to several vegetation models. We present a systematic guide for finding such patterns in a given reaction--diffusion model, during which we obtain four key quantities that allow us to predict the qualitative properties of our solutions with minimal analysis. We consider four well-established vegetation models and compute their key predictive quantities, observing that models which possess similar values exhibit qualitatively similar localised patterns; we then complement our results with numerical simulations of various localised states in each model. Here, localised vegetation patches emerge generically from Turing instabilities and act as transient states between uniform and patterned environments, displaying complex dynamics as they evolve over time.


Introduction
In semi-arid environments, the formation of distinct vegetation patterns has been well documented since their first discovery in sub-Saharan Africa in the 1950s [42,43].Sloped terrains exhibit stripes (or bands) of vegetation [13,14,54,58], while flat terrains are home to periodic arrangements of spots and gaps, as well as labyrinthine patterns [4,22,23,41,51]; see [46] for a detailed review of the mechanisms that cause vegetation patterns.The environmental role of these structures is still debated: are they a sign of desertification, providing early indications of environmental decline [52], or are they nature's safety net, allowing an ecosystem to avoid critical tipping points and stave off extinction [50]?Such questions have generated significant interest in vegetation patterns over several decades, and this looks set to only increase in the face of our changing climate.
One of the key questions when first studying vegetation patterns is what kind of model to consider-the dynamics of an ecosystem are highly complex, and any attempts to model a particular environment will need to reflect the specific qualities of that environment.Vegetation with vertical roots and strong water uptake are subject to a pattern-forming feedback driven by diffusion, while vegetation with widely spread roots are also subject to nonlocal effects, due to water uptake by the roots [46].Any attempts at modelling these types of environments must also reflect these differences; see [32,35,64,66] for diffusive PDE models, and [15,20,48,61] for nonlocal models.
Even with these environment-specific modelling terms there are certain phenomena that appear across many vegetation models, which we call universal phenomena; we can better understand these universal phenomena by identifying the mechanisms that cause them.For example, periodic and labyrinthine patterns are known to emerge from bifurcations with non-zero wavenumber, which we refer to as Turing bifurcations, and so most Figure 1: Example of a localised vegetation pattern with hexagonal symmetry, with vegetation density plotted in the top row and water density plotted in the bottom row.The left-hand column provides a top-down view of the pattern, while the right-hand column provides a three-dimensional view.Here the vegetation has a negative polarity and so consists of gaps, whereas the water has a positive polarity and is thus made up of peaks.These images were produced using VisualPDE [67].
flat terrain vegetation models will undergo such bifurcations for a critical choice of parameters.The presence of a universal phenomenon (such as periodic patterns) implies the presence of a universal mechanism (such as a Turing bifurcation) within a phenomenological model.
A particular phenomenon that we are interested in is the emergence of fully localised patterns, where a patterned state is completely surrounded by a uniform state.A notable example of these are so called 'fairy circles', consisting of roughly circular barren patches surrounded by vegetation [19].Fairy circles have been found to exhibit a consistent wavelength between neighbours, and can be found in fully localised patches or as constituents of a larger periodic structure [17].Originally observed in the Namib desert [65], fairy circles have since been discovered in Western Australia [18], suggesting that they may be more universal than originally supposed.
While localised patterns are well understood mathematically in one spatial dimension [10,69], this is not the case in two dimensions.There has been recent progress on axisymmetric patterns both in vegetation models [7,25] and other contexts [29,37,44,45], mostly using radial spatial dynamics theory pioneered by Scheel [55].However, there is very little theory regarding fully localised non-axisymmetric patterns outside of variational methods [6] and exponential asymptotics [34].We focus on fully localised patterns with dihedral symmetry-such that they are invariant under rotations of 2π/m, for some m ∈ N, and one line of reflection-which includes localised hexagons and squares observed in prototypical pattern-forming systems [38,53].
We consider a simple class of two-component reaction-diffusion models that cover a wide range of vegetation models, where pattern formation is driven by diffusion and domain-covering patterns emerge from Turing bifurcations.For this general class of models, we show that fully localised patterns with dihedral symmetry emerge universally from a Turing bifurcation.To do this, we first present a step-by-step procedure to determine the local expression of a reaction-diffusion system near a Turing bifurcation, which then allows us to leverage recent rigorous results in [27,28] regarding the existence of approximate fully localised dihedral patterns.We present theorems for the existence of two classes of localised dihedral patterns: spot A-type patterns, which can appear as either peaks or gaps, and ring-type patterns that emerge from a pitchfork bifurcation.If a pattern mostly consists of peaks we say it has a positive polarity, and a negative polarity if it mostly consists of gaps; see Figure 1.
This process also reveals four key quantities P 1 , . . ., P 4 , which we call qualitative predictors, that establish the qualitative behaviour of these localised patterns.In particular, (i) P 1 determines the direction of bifurcation from the Turing point; see Figure 2(a).If P 1 > 0 then localised solutions bifurcate as the precipitation increases and if P 1 < 0 then localised solutions bifurcate as the precipitation decreases.
(ii) P 2 determines the relationship between vegetation and water densities of spot A-type patterns; see Figure 2(b).If P 2 > 0 then the vegetation and water share the same polarity and if P 2 < 0 then they have opposing polarity.Such vegetation patterns are respectively called in-phase or anti-phase, both of which have been observed in various semi-arid climates; see [18], for example.This work represents the first analytic evidence for fully localised non-axisymmetric patterns in vegetation models, to the author's knowledge, and establishes a guide for quickly analysing such patterns in other two-component reaction-diffusion systems.We emphasise that our theoretical results are regarding the existence of localised steady-state patterns for reaction-diffusion systems, and so they do not determine the stability of such solutions.As such, these patterns instead represent transient states as a uniform vegetated state transitions to a patterned state; we later observe the nonlinear dynamics of this transition in Section 3.
Since the time scales associated with vegetation models are often larger than other pattern-forming systems, these transient states are more likely to be observed in the field.
We find that spot A-type patterns are highly universal, emerging generically from a Turing bifurcation, whereas ring-type patterns only emerge in narrow parameter regions.In each of the vegetation models we consider spot A-type patterns grow in width over time, invading the surrounding uniform state.As such, these localised patterns can represent a degradation of a uniform vegetated state resulting in a significant change to the environment, such as desertification.
The rest of the paper is structured as follows: in Section 2.1 we derive the local expression for a general two-component reaction-diffusion system near a Turing instability, presenting a step-by-step process that can be easily automated.In Section 2.2 we then introduce four notable reaction-diffusion models for dryland vegetation, and compute their qualitative predictors in each case.Finally, in Section 3 we use exponential time-steppers to numerically simulate three examples of localised dihedral patterns in each vegetation model, comparing the similarities and differences between different models.

Theory
We begin this section by taking a general two-component reaction-diffusion model and deriving its local representation in the neighbourhood of a Turing instability.This allows us to utilise abstract theorems to predict the existence of solutions in various vegetation models.Throughout this procedure we identify several notable quantities which describe the qualitative behaviour of our localised solutions; computing these key quantities then allows us to make heuristic predictions regarding the emergence of localised patterns in each specific reaction-diffusion model.Each stage of the procedure is reduced to solving an algebraic condition, thus the entire process can be easily automated for rapid predictions regarding the emergence of localised patterns in two-component reaction-diffusion models.
Following this, we introduce four vegetation models that belong to our general class of two-component reaction-diffusion equations: the Klausmeier-Gray-Scott, logistic Klausmeier, NFC-Gilad, and von Hardenberg models.For each model we compute our key quantities and predict the behaviour of localised solutions bifurcating from each Turing point, we will then compare these predictions with the numerical results presented in Section 3.

Derivation
In general, a reaction-diffusion equation has the form Here, u(t, x) ∈ R N denotes N coupled variables, with temporal and spatial coordinates t ∈ R and x ∈ R 2 , respectively, and µ ∈ R p denotes p parameters of the system, for some N, p ∈ N. The two-dimensional Laplacian ∆ := ∂ 2 x + ∂ 2 y encodes the spatial diffusion of u, and the N × N dimensional square matrix D ∈ R N ×N defines the diffusion coupling of the system (i.e., how the spatial diffusion of each variable affects the overall temporal growth of u).Finally, the function F ∈ R N accounts for all local reactions in the model.
For this present work we assume that u(t, x) is made up of two variables u(t, x), v(t, x) ≥ 0 which represent the densities of the vegetation biomass and soil water, respectively, and that the diffusion matrix D ∈ R 2×2 is lower triangular, such that the rate of growth of the biomass is independent of the diffusive rate of the water.This assumption models the effect of vegetation on water diffusion through root suction.We also assume that p − 1 parameters are fixed, such that any bifurcations are driven by a single parameter µ ∈ R, which we call the bifurcation parameter (in vegetation models µ is usually the precipitation rate).
As such, the remainder of this work will focus on general reaction-diffusion systems of the form We have defined the system (2.2) so that subsequent analysis will be slightly easier; a more general system like could equivalently be written in the form of (2.2) by rescaling the coordinate system (t, x) with respect to D uu and then defining parameters D v , β and functions f , ĝ accordingly.We are interested in steady states of the system (2.2) and so, for the rest of this section, we will consider the following equations where we have defined In Section 3, we will then return to the full time-dependent PDE system (2.2) for our numerical investigations.
The aim for this section is to use abstract theory to predict the emergence of fully localised dihedral patterns in (2.3) near a Turing instability.By Turing instability, we mean a parameter value at which a uniform steady state loses stability with respect to perturbations with nonzero wave number k > 0; sometimes the term Turing instability is reserved for instabilities caused by changes in the diffusion coefficient D, but we will not use that definition in this paper.To do this, we present the following procedure: ( §2.1.1)Identify uniform steady states of the system (2.3).
( §2.1.2)Compute the spatial eigenvalues of the system (2.3) at the uniform steady state and identify a critical parameter value µ * (k) such that (2.3) undergoes a Turing bifurcation with wave number k > 0.
( §2.1.3)Express the system (2.3) as a Taylor expansion about the uniform steady state at the Turing point µ = µ * (k).
Once we have transformed the system (2.3) into local coordinates in the neighbourhood of a Turing instability, we identify four key quantities P 1 , . . ., P 4 that help us predict various properties of localised dihedral patterns bifurcating from the Turing point ( §2.1.5).

Finding uniform steady states
The first task is to identify uniform solutions of the steady state equation (2.3).That is, we look for solutions that are independent of the spatial coordinate x, such that ∆u = ∆v = 0.Then, the system (2.3) becomes for which we want to find solutions u = u * (µ), v = v * (µ).In general, while these equations can be difficult to solve by hand, especially with additional fixed parameters in a given model, they are very tractable for numerical solvers.For the vegetation models presented later in this section, the problem reduces to solving a cubic-order polynomial, and so one would typically expect to find between one and three uniform solutions to (2.4) for any fixed µ ∈ R.
Having identified some uniform steady state (u, v; µ) = U * (µ) := (u * (µ), v * (µ); µ), we can restrict our analysis to a small neighbourhood of U * (µ) by taking the Taylor series of f and g.We note that, in order to perform weakly nonlinear analysis in Section 2.1.4,we require explicit formulae for the quadratic and cubic nonlinear terms in the Taylor expansion of f and g.Thus, we define u(x) = u * (µ) + U (x) and v(x where U := (U, V ) T and (2.6) Here, the matrix M(µ) is the linear operator of (2.5), F 2 contains all quadratic (in U) nonlinearities of (2.5), and F 3 contains all cubic (in U) nonlinearities of (2.5).Note that we do not truncate the Taylor series at cubic order-the O(|U| 4 ) terms remain in our problem but do not affect our existence results, and so we choose not to write them down explicitly.

Spatial eigenvalue analysis
We now wish to compute the spatial eigenvalues of the linear operator for the system (2.5); in particular, for fixed µ ∈ R we wish to find λ ∈ C such that This is equivalent to solving the eigenvalue problem of M(µ) for λ and so we need to solve σ(λ; µ) = 0, where We briefly comment on how σ(λ; µ) connects to the standard dispersion relation found in one-dimensional Turing analysis.Taking the original time-dependent system (2.2) for u ∈ R 2 , restricted to one spatial dimension, and introducing the ansatz u(t, x) = u * + û e ωt e ikx + c.c. with u * a uniform steady state, we obtain the following dispersion relation Then the condition σ(λ; µ) = 0 corresponds to the critical value ω = 0, such that the uniform state u * changes stability with respect to perturbations with spatial eigenvalues ± √ λ.Using the above definition of M(µ), we find The point µ = µ * ∈ R defines a Turing instability if σ(−k 2 ; µ * ) = 0 for some k ∈ R, σ(λ; µ * + ε) has no real roots, and σ(λ; µ * − ε) has two distinct negative roots for 0 < |ε| ≪ 1; see Figure 3.Note that we have λ σ(λ;μ) adopted the convention that σ(λ; µ * + ε) has no real roots but we have not assumed the sign of ε, and so ε could be positive or negative.Thus, µ = µ * is a Turing point if possesses repeated real roots, or equivalently Supposing that µ * satisfies the above condition, the wave number k ∈ R is then determined to be We note that k 2 < 0 would imply that µ * is not a Turing point, and so we proceed under the assumption that the right-hand side of (2.9) is positive.

Local coordinates near a Turing instability
We centre our system (2.5) at the point µ = µ * where the uniform state U * (µ) undergoes a Turing instability.
To do this, we define µ = µ * + ε, for 0 < |ε| ≪ 1, such that (2.5) becomes (2.10) Here, we have introduced new objects M 1 , M 2 , Q, C to bring our system in line with those considered in [27,28].In particular, we define M 1 := M(µ * ) and M 2 := M ′ (µ * ) as the first terms in the Taylor expansion of M(µ * + ε).For the nonlinear terms F 2 , F 3 we have defined , and Q and C are respective symmetric bilinear and trilinear functions, i.e. they satisfy There are a number of quantities in (2.10) that play a part in the existence theorems of [27,28]; these are mostly not unique to the existence of localised dihedral patterns and so we present them here before starting the next part of our procedure.We begin by defining the following vectors We now introduce the following constants where we have written The constant c 0 determines the direction of bifurcation from the Turing point-we will see this explicitly in the upcoming theorems-while the constant c 3 determines whether domain-covering stripes (otherwise known as rolls) emerge via a subcritical or supercritical pitchfork bifurcation [30].This second definition may seem irrelevant in our current study; however, it was previously determined that one dimensional localised patterns only bifurcate from the Turing point if rolls bifurcate subcritically [69].Finally, we introduce the constant γ, where In the study of axisymmetric patterns, the sign of γ determines whether localised spot solutions have an elevation (a peak) or a depression (a gap) at their centre [29,37].We note that the quantities c 0 , c 3 , γ and the eigenvector Û0 will be important when defining our qualitative predictors P 1 , . . ., P 4 later in this section.
We are now equipped to study the existence theorems for localised approximate dihedral patterns seen in [27,28].

Approximate localised dihedral patterns
We consider fully localised patterns with D m dihedral symmetry, so that where r, θ are the standard polar coordinates for R 2 .We employ a Galerkin scheme for approximating D m solutions of the system (2.10), which we now detail.Taking a truncated Fourier expansion for a fixed choice of N ≥ 0 and projecting onto each cos(mnθ), (2.10) becomes for each n ∈ [0, N ], where h.o.t.s denotes the higher order terms seen in (2.10) and ∆ n := (∂ 2 r + 1 r ∂ r − (mn) 2 r 2 ) is the Laplacian operator ∆ when applied to the n th Fourier mode.Note that our introduction of Q and C has provided a much cleaner expression of the nonlinear terms in (2.17) than if we used F 2 and F 3 .
We now write down existence theorems for localised solutions to (2.17) originally proven in [27,28].We shorten our statement of the theorems for the sake of brevity and readability, presenting only the key points relevant to this work.First, we consider spot A-type solutions, whose existence was the focus of [27]: then there exists some r 0 > 0 such that the Galerkin system (2.17) has an approximate localised dihedral solution for r ∈ [0, r 0 ] and |U A (r, θ)| → 0 exponentially fast as r → ∞.Here J n is the n th order Bessel function of the first kind and the constants {a n } N n=0 are nondegenerate solutions of the quadratic matching condition ) Setting N = 0 recovers the axisymmetric spot A solution found in [37] for the Swift-Hohenberg equation, which is why we refer to these solutions as spot A-type patterns.The proof of Theorem 1 requires tools from radial spatial dynamics, pioneered by Scheel [55], and can be found in Section 4 of [27].Solving the algebraic matching condition (2.19) then provides excellent initial approximations of fully localised dihedral patterns for numerical study; we will use an initial guess of (2.18) with solutions from (2.19) as the starting point for our numerical simulations in Section 3.
Beyond the spot A-type dihedral patterns found in [27], we can also find localised dihedral ring patterns for the system (2.17); the existence of these patterns was the focus of [28], resulting in the following theorem.
Theorem 2 ( [28]).Fix m ∈ N, N ∈ N 0 .If c 3 < 0, then there exists some r 0 > 0 such that the Galerkin system (2.17) has an approximate localised dihedral solution U R (r, θ) ) Setting N = 0 recovers the axisymmetric ring solution found in [37] for the Swift-Hohenberg equation, and so we refer to these solutions as ring-type patterns.The proof of Theorem 2 follows in a similar way to Theorem 1, but with more delicate analysis regarding the exponential decay of solutions for large values of r.
As discussed earlier, these ring-type patterns only emerge when bifurcating stripes are unstable; in particular, their far-field profile is approximated by the localised solution of the cubic nonautonomous Ginzburg-Landau where s := (c 0 ε) r is a rescaled radial coordinate in the far-field.An exponentially decaying solution to (2.22) only exists when c 3 < 0, as proven in [44,63].The algebraic matching condition (2.21) possesses Z 2 symmetry, so −U R (r, θ) is also a solution of (2.17) and ring-type patterns emerge in a pitchfork bifurcation.A numerical study of these ring-type patterns first requires solving (2.22) in order to construct an effective initial guess; we will restrict our numerical simulations in Section 3 to the simpler spot A-type patterns.

Qualitative predictors
We conclude this section by highlighting the following quantities, where the notation [v] i denotes the ith element of a vector v.Each quantity can be computed directly from the original model (2.2) and describes a qualitative property of our localised solutions.We recall, (i) P 1 determines the direction of bifurcation from the Turing point.If P 1 > 0 then localised solutions bifurcate for ε > 0 and if P 1 < 0 then localised solutions bifurcate for ε < 0.
(ii) P 2 determines the relationship between u and v.If P 2 > 0 then u and v share the same polarity and if P 2 < 0 then u and v have opposing polarity.
(iii) P 3 determines the polarity of u.If P 3 > 0 then spot A-type patterns are made up of peaks and if P 3 < 0 then spot A-type patterns are made up of gaps.
(iv) P 4 determines whether or not ring-type patterns exist.If P 4 < 0 then ring-type patterns bifurcate from the Turing point, and if P 4 > 0 then no ring-type patterns emerge.
Let us briefly discuss some intuition for each predictor.First, the sign of P 1 fixes the sign of ε via the condition that σ(λ; µ * + ε) has no real roots for |ε| ≪ 1.One can formally derive this by a linear transformation of σ(λ; µ * + ε) into { Û0 , Û1 } coordinates, so that Then, we see that λ becomes complex when P 1 ε > 0. The quantities P 2 , P 3 both appear in the profile of our spot A type solutions (2.18), which we can write as and so the role of both P 2 , P 3 can be seen explicitly.The final predictor P 4 has already been discussed previously, where we observed that it corresponds to a focussing/defocussing condition in the non-autonomous Ginzburg-Landau equation (2.22).We will see that these predictors provide a quick guide by which to compare and categorise different reaction-diffusion systems with regards to the formation of localised dihedral patterns.

Vegetation models
We now introduce several models for dryland vegetation and compute the quantities P j for j = 1, . . ., 4 in each model.This will then allow us to predict qualitative properties of the types of localised patterns emerging in each model.We again recall that our qualitative predictions are regarding the existence of localised steady-states, which may represent transient states between uniform and patterned environments.
We then support our analysis with a numerical exploration of some spot A-type patterns in Section 3, where we verify our qualitative predictions and observe how localised patterns evolve in each of the following vegetation models.

Klausmeier-Gray-Scott model
We begin with the system introduced by van der Stelt et al. in [64] and commonly referred to as the 'extended Klausmeier' model [58,59], or the 'Klausmeier-Gray-Scott' model [16,36,68,71].The first name is in reference to the Klausmeier model for vegetation patterns on sloped domains, where m > 0 models an effective death rate of vegetation and ν models the advection of water down a uniform hillslope [32].The Klausmeier model can be extended to flat land by setting ν = 0 and including a diffusive term for the water density, thus resulting in (2.24).The name 'Klausmeier-Gray-Scott' refers to the similarities between (2.24) and the Gray-Scott model for chemical reactions [24].We will presently introduce another model derived as an extension of the Klausmeier model, and so we will call (2.24) the Klausmeier-Gray-Scott model to reduce confusion.
The Klausmeier-Gray-Scott model (2.24) is described by the general reaction-diffusion model (2.2) with For this particular model, the functions f , ĝ are simple enough to analyse explicitly, which is done in Appendix A. We discover that (2.24) undergoes a Turing bifurcation if δ v m > 2, and The sign of P 4 varies for different values of δ v , m, and so we numerically plot a diagram for sgn (P 4 ) in Figure 4(a).
The benefit of finding explicit expressions for each P j is that we can judge the robustness of each qualitative property; i.e., how sensitive any numerical solutions are with respect to our choice of parameters.What we find is that the direction of bifurcation (P 1 ), whether solutions are in-phase or anti-phase (P 2 ), and whether the solutions are peaks or gaps (P 3 ) remain consistent across all choices of δ v , m such that δ v m > 2.
Figure 4(a) shows that the subcriticality condition P 4 < 0 only holds for sufficiently small values of δ v , given a fixed m, and so the ring-type patterns found in Theorem 2 emerge in a much smaller parameter region than the spot A-type patterns of Theorem 1.
For our subsequent numerical study in Section 3, we choose the following parameter values taken from [71], where localised periodic states have been observed previously; we will use these parameter values for the remainder of this work when considering (2.24).Then, we identify the Turing point We note that the direction of bifurcation, as predicted by P 1 , does not tell the whole story; the bifurcation curve in Figure 4(b) could undergo a fold and continue into the region where µ < µ * .This is the case for certain localised one dimensional solutions in the Klausmeier-Gray-Scott model [71], and so we might expect the same to be true in this case.

Logistic Klausmeier model
Our second example is the logistic Klausmeier model introduced by Bastiaansen et al. [3], which is an extension of the previous Klausmeier-Gray-Scott model (2.24), where the unrestricted vegetation growth term vu 2 is replaced by a logistic growth term (1 − bu)vu 2 with carrying capacity b −1 > 0.
This adjustment limits the growth of vegetation at any given point, which is motivated as follows.For extremely low precipitation levels it is reasonable to assume that vegetation growth is constrained by the amount of soil water, as modelled by (2.24); however, in environments that can support larger quantities The additional nonlinear term makes any explicit calculations significantly more cumbersome, and so we just present our numerical calculations of P j .We fix b = 1 and choose the following parameter values used to model grass patterns in [32], where the value of the diffusion constant δ v is significantly larger than in the previous example.Then, we find the Turing point where w(t, x) denotes the density of water on the surface of the soil and G u (t, x), G v (t, x) are non-local terms that describe vegetation growth and soil water consumption, given by 0 (1 + ηu(t, x)) 2 u(t, x ′ ) dx ′ , respectively.One of the aims of this model was to describe the behaviour of plants with widely spread out roots, such that the interaction between vegetation and water is no longer strictly local, resulting in a three-component non-local PDE system.As such, this model is not covered by our results unless some simplifications are first performed.
We consider a simplification of the Gilad model introduced by Zelnik et al. in [72], where the authors focused on studying fairy circles on flat terrains in Namibia.For this particular environment, the sandy soil results in very fast changes to the surface water w(t, x)-which converges to an equilibrium state-and the remaining slow dynamics of the system can be reduced to the two-component reaction-diffusion system (2.29) We refer to this model as the 'Namibian fairy circle Gilad' model, or 'NFC-Gilad' model for short.It is exactly of the form in (2.2), with D v = δ v , β = 0, and (2.30) Again, the functions f , ĝ are complicated enough that we only present numerical calculations here.We choose the following parameter values Λ = 16  35 , η =

von Hardenberg model
The final model we consider is the von Hardenberg model introduced in [66], which reproduces many symmetry-breaking transitions found in field observations.In particular, simulations of domain-covering patterns undergo 'gaps → labyrinth → spots' and 'spots → stripes → gaps' transitions [22,23], and simulations of localised patterns undergo 'spot → ring' and 'ring → spots' transitions [47].We highlight some of the terms in (2.32) not seen in the previous models.The coefficient γv 1+σv models vegetation growth due to water intake, which is linear for low levels of water and tends to a constant for larger levels of water; this is instead of the constant term m in (2.24) and (2.26) or the linear term Λv in (2.29).The ∆(v − βu) term models the transport of water via Darcy's law; this includes the suction of water by plant roots at a rate of β > 0, which is not covered by the previous models.

The von Hardenberg model (2.32) is of the form in (2.2), with
and we again only present numerical calculations for P 1 , P 2 , P 3 , P 4 here.We choose the following parameter values γ = 1.6, σ = 1.6, ν = 0.2, ρ = 1.5, δ v = 100, β = 3, (2.34) taken from [66].Notably, we obtain two Turing points for the von  15].As discussed earlier, these one dimensional patterns require c 3 < 0 in order to emerge, and so we expect dihedral ring-type patterns to also emerge in those parameter regimes.

Numerical Results
In the previous section we demonstrated how fully localised patterns emerge in general two-component reaction-diffusion systems near a Turing instability, which we now support with a brief numerical study on the evolution of spot A-type dihedral vegetation patterns in each of the vegetation models introduced previously.The goal of this section is to numerically explore spot A-type vegetation patterns in the four vegetation models introduced previously, in order to (1) verify our qualitative predictions from Section 2.2; (2) explore the role of localised patterns as transient states of vegetation models; and (3) observe how different models with similar initial conditions can result in qualitatively different evolution dynamics.
The numerical codes for producing the results of this section are all available at [26], along with codes for calculating the predictive quantities P 1 , . . ., P 4 in each vegetation model considered here.Furthermore, we provide others with the option to define their own reaction-diffusion models of the form (2.2) and determine the qualitative properties of patterns emerging from a given Turing point.

Implementation
Our numerical approach is as follows: for a system of the form (2.2) we provide the functions f (u, v, µ), ĝ(u, v, µ), including the values of any other parameters, and an initial guess for the Turing point (u * , v * , µ * ).We numerically solve the algebraic equations detailed in Section 2 in order to find a Turing point close to our initial guess, and compute k, Û0 , Û1 , γ, c 0 , c 3 .Then the the values of P 1 , P 2 , P 3 , P 4 can be computed, as in the previous section for each vegetation model.
Next, we choose what type of dihedral pattern we want to find by fixing m, N ∈ N in Theorem 1 and providing an initial guess for the algebraic matching condition (2.19) for spot A-type patterns.This results in an initial profile (2.18) for a localised dihedral pattern bifurcating from the Turing instability; we then input this solution into our time-stepping codes in order to investigate how these solutions evolve over time.
We employ a second-order exponential time-differencing method from Asanté-Asamani et al. [2] for a reaction-diffusion system of the form with homogeneous Neumann boundary conditions.In order to convert our general system (2.2) into this form, we apply the transformation so that our general system of the form In particular, we apply our time-stepping codes to the system Each profile is an approximate localised dihedral pattern obtained in [27] and together they form the initial guesses for our simulations.
and then recover u, v from the transformation u = û, v = βDv Dv−1 û + v.We compute our simulations on a square domain of length 20λ, where λ := 2π k is the wavelength related to each Turing bifurcation.This k−dependent choice of domain means that any patterns found across different models should be of a similar size when plotted.We present our simulations through snapshots of our solutions in time; videos of each simulation can be found at [26].

Simulations
We present three examples of localised dihedral patterns for our numerical study: hexagons (D 6 ), squares (D 4 ) and pentagons (D 5 ).In we provide an initial guess  See Figure 7 for the profile of u − u * in each case.Note that, while we refer to each pattern by the polygon associated with its dihedral symmetry group-i.e., hexagons, squares and pentagons-the patterns do not necessarily appear similar to these polygons.For example, the 'square' pattern is related to a twelve-fold quasipattern, rather than a square lattice.
These three examples of dihedral patterns each possess a different level of rarity in pattern formation; hexagons are very common in experiments and field observations, square lattices are common but our specific square example is not, and pentagons are not very common at all.However, Theorem 1 tells us that each example emerges from the same Turing point in the same parameter regime, and so it would be interesting to compare the time evolution for each example and observe whether the square and pentagon patterns converge to the more common hexagons over time.
Since each of the models we consider is nondimensionalised, we cannot readily compare the amplitude of patterns across different models.However, we observe that, for a given model, each of our initial dihedral patterns converge to roughly the same amplitude, and so we provide an approximate range of values for u in each model.The videos available at [26] show an individual colorbar for each pattern, thus providing more details for enquiring readers., where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at [26].
The Turing bifurcation has an associated wave number k ≈ 0.3177, and so we simulate our patterns in a square box of length 396.The simulations of our three examples-hexagons, square, and pentagons-are presented in Figure 8.We note that the patterns found for (2.24) are all anti-phase, as predicted in Section 2 and shown in the final column of Figure 8; we do not otherwise plot the water density v since the spatial structure is identical to the vegetation density u.
For the hexagon example, solutions form a regular hexagonal lattice of gaps.As time evolves, new gaps emerge on the outer layer of the patch such that the width of the pattern grows continuously.This behaviour is what we would expect for localised patterns that undergo homoclinic snaking (for example [5,27,38] for snaking solutions and [39,40] for growing patterns away from the snaking region) and is also observed in the square and pentagon examples, but with very different structures being formed.
In the square example, the initial pattern contains regions of the uniform state between individual gaps; , where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at [26].
emerging gaps first fill these regions until the resulting pattern is compact, before then emerging on the outer layers.Whereas the hexagon example only consists of circular gaps, the square pattern is also comprised of long thin gaps.For later times, the square pattern is reminiscent of the penta-hepta defects and grain boundaries studied in [60]; this is when a point in an almost-hexagonal lattice has five or seven neighbours, rather than the standard six neighbours.
For the pentagon example, the story is much the same as the square example.Gaps first fill out the initial patterned region before then emerging on the outer layer.Notably, the five gaps at the core of the pattern coalesce and appear to converge to a ring, while the outer layer of the pattern traces out a regular pentagon as times increases.
The Turing bifurcation has an associated wave number k ≈ 0.1612, and so we simulate our patterns in a square box of length 780.The simulations are presented in Figure 9, where we again mostly plot the vegetation density u since the water density v exhibits the same spatial structures in anti-phase., where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at [26].
The logistic term in (2.26) restricts the height of individual gaps, causing the gaps to grow in width instead; this was also observed in [3] for one dimensional patterns.As a result, nearby gaps expand and coalesce, causing a collapse from a patterned state to a compact region of bare soil.New gaps emerge on the periphery of this compact region, in the same positions as in the Klausmeier-Gray-Scott model (2.24), but are soon absorbed into the central gap as time evolves.
Over time these patterns resemble radial fronts connecting the bare state to the vegetated state, similar to the axisymmetric and one dimensional fronts studied in [3,7,9].However, the geometry of the front interface is highly nontrivial, and so it is unclear whether these structure could be analysed using similar techniques.
The Turing bifurcation has an associated wave number k ≈ 0.333, and so we compute our simulations in a square box of length 378.The simulations for our three examples are presented in Figure 10 for the vegetation density u; solutions are again found to be anti-phase, and so we mostly do not plot the water density v since its spatial structure is identical to u. , where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at [26].
Much like with the logistic Klausmeier model (2.26), individual gaps in the NFC-Gilad model grow in width as time evolves.However, unlike the logistic Klausmeier model, these gaps do not coalesce with their neighbours.As a result, we observe patterns with non-uniform elliptic gaps that retain a mosaic-like structure.We note that the interface between the pattern and the surrounding vegetated state is much more circular for (2.29) than in our previous models, even for the square and pentagon examples.
We note that all three examples exhibit a pseudo-hexagonal lattice structure away from their core, which might hint towards the dominance of hexagonal patterns in field observations and experiments.One could then think of the resultant square and pentagon patterns as versions of the hexagonal example with defects at their core.However, it is worth noting that the hexagonal packing observed in the square and pentagon examples is not symmetry breaking; the dihedral symmetry of each example (D 4 for squares, D 5 for pentagons) is preserved throughout, except for when the pentagon example becomes large enough to experience boundary effects.
This bifurcation has an associated wave number k ≈ 0.106, and so we compute our simulations in a square box of length 1190.However, the width of the patterns is considerably smaller than this domain, and so we present our results in the smaller square box of length 595.The water density v is found to be in-phase, as predicted in Section 2, and shares an identical spatial structure to the vegetation density u.Hence, we present our simulations of u in Figure 11.The evolution of these patterns occurs at a slower rate than in previous models, and so we present our simulations at later times, as detailed in Figure 11.
The initial guess (3.6) is very small for this Turing bifurcation, and so it can be difficult to capture numerically.As such, in order to observe localised dihedral patterns in this model, we take a larger choice for the constant C in (3.6).Then, for each example solutions converge to localised peaks spread out in dihedral arrangements.The emergence of new peaks is much slower than in the previous models; individual spots instead grow in width until their core destabilises, causing a transition into a ring or two smaller spots.This behaviour is reminiscent of localised spikes in singularly perturbed reaction-diffusion systems, such as [33], rather than the invading patterns seen in our other models.We note that we observe these spot-to-ring transitions occur even though localised ring-type patterns do not emerge from the uniform state.
Other than the central spot in the hexagon example, every spot is observed transitioning into two smaller spots.The way that each spot splits appears to be very structured; the outer ring of spots appears to split in the azimuthal (i.e.angular) direction, whereas the inner ring of spots appears to split in the radial direction.
Since the outer layer of spots also appears to split symmetrically, we note that the overall dihedral symmetry of the system seems to be preserved by these splittings.
This bifurcation has an associated wave number k ≈ 0.201, and so we compute our simulations in a square box of length 612.Again, we find the water density v to be in-phase, as predicted in Section 2, with an identical spatial structure to the vegetation density u; we present our simulations of u in Figure 12.Again, the time evolution of our patterns is slower than the other models, and so we present our simulations for slightly later times.
In each example solutions are comprised of circular gaps with a uniform size.The hexagon example results in a perfect hexagonal lattice, exhibiting the same behaviour as in the Klausmeier-Gray-Scott model and prototypical pattern-forming systems like the Swift-Hohenberg equation [38].The square and pentagon examples grow similarly to the Klausmeier-Gray-Scott (2.24), but new gaps exhibit hexagonal packing like in the NFC-Gilad models (2.29).Then, each example evolves into a perfect hexagonal lattice with possible defects; these defects then form grain boundaries throughout the localised pattern, such as the horizontal and vertical lines in the square example.As with the previous models, even though the solutions are all evolving into distorted hexagonal lattices, the dihedral symmetry of the solution is preserved.

Discussion
In this work, we have provided a step-by-step guide on how to convert a two-component reaction-diffusion system into its local representation near a Turing instability for a class of vegetation models.Then one can find localised dihedral patterns bifurcating from the Turing instability and obtain four quantities that encode the qualitative behaviour of these localised solutions and how they bifurcate.We found that localised spot A-type patterns emerge generically from a Turing instability, not requiring any subcriticality conditions like localised one dimensional and ring-type patterns.
We then considered four established models for dryland vegetation and predicted the emergence of localised dihedral patterns.For the simple Klausmeier-Gray-Scott model (2.24), we were able to derive explicit forms for each qualitative predictor in terms of the parameters of the model, and hence better understand the qualitative robustness of localised dihedral solutions with regards to our choice of parameters.We found that each model exhibits spot A-type patterns bifurcating from a Turing point, but ring-type patterns do not emerge for the given choices of parameters from the literature.However, in the Klausmeier-Gray-Scott model we saw that there exists a small parameter region in which ring-type patterns emerge from the uniform state, which is also observed for the von Hardenberg model in [12,Figure 15].We then supported our predictions through numerical simulations of each vegetation model, where we observed that the generic spot A-type patterns evolve into dramatically different structures depending on the choice of model.
We note that the dihedral symmetry of each solution is very robust in our simulations, even though the Cartesian discretisation of our spatial domain should actively disfavour the rotational symmetry of our patterns.One possible reason for this might be that the emergence of each peak or gap is strongly influenced by its neighbours.Then, if we begin with a collection of peaks or gaps arranged with dihedral symmetry, the neighbourly effects on a new peak or gap also possess dihedral symmetry, thus enforcing the dihedral structure on the new peaks and gaps.Another possible reason might be that our initial guess imparts a dihedral symmetry on the whole domain which the localisation just makes extremely small away from the origin.Then, as the uniform state destabilises, the already present dihedral oscillations grow in amplitude, resulting in a pattern with dihedral symmetry.
Although our solutions maintain their dihedral symmetry as time evolves, we observe significant differences in the resultant patterns for different vegetation models.For example, if you were looking to model localised in-phase gaps (such as fairy circles observed in Australia [18]), then any prospective reaction-diffusion model must have a Turing instability with P 2 > 0, P 3 < 0. Such a model would not be appropriate for the fairy circles observed in Namibia, however, which manifest as anti-phase gaps and thus require a model with P 2 < 0, P 3 < 0. Furthermore, we note that any vegetation model that exhibits a transition from spots to gaps (or vice versa), such as the von Hardenberg model [23], must thus possess two distinct Turing points with different signs for P 3 .While these properties alone are not sufficient to successfully model complex phenomena such as vegetation patterns, they can provide helpful conditions to be used in conjunction with ecological modelling.We also note that the spot A-type patterns emerge generically from a Turing bifurcation, and so the mere existence of localised patterns is not necessarily a good indicator of whether or not a model is appropriate, since if a model can exhibit two-dimensional periodic patterns it can also exhibit localised dihedral patterns.
Beyond the design of new models, our work also provides some insight into the formation of vegetation patterns.In our simulations, localised patterns represent a transition between two states; a stable patterned state and an unstable uniform state.The patterned state emerges at the origin and begins to invade the uniform state in a radial direction, as seen prototypical pattern-forming systems [39,40]; this is particularly serious in the logistic Klausmeier model (2.26) where the invading patterned state collapses to the bare state, resulting in the desertification of the local environment.Since localised dihedral patterns lead to the further destabilisation of the uniform vegetated state, they may serve the same role as periodic patterns in predicting or mitigating tipping points and desertification effects in semi-arid environments [50,52].
We briefly comment on the popular phenomenon of fairy circles and their relation to the localised gaps discussed in this work.In [19], Getzin et al. define fairy circles as satisfying the following conditions; (a) they must be gaps in a vegetated state, that (b) can form periodic structures, and (c) are not rings surrounded by a bare state.We have found that, for any Turing bifurcation with P 3 < 0, there are localised dihedral patterns that satisfy these criteria, as seen in Figures 8,9,10 & 12. Furthermore, we have found that these gap solutions are distinct from ring-type patterns found in Theorem 2, which emerge in much more restrictive parameter regimes.
We note that there are numerous possible directions in which to numerically explore localised vegetation patterns, many of which we do not consider here.The theory in Section 2 describes the bifurcation of localised patterns from a Turing instability and so a natural first step would be to use numerical continuation schemes to characterise the bifurcation curves of these patterns in more detail.One could adapt the pseudo-spectral codes of [27] to each vegetation model, using a truncated angular Fourier expansion and finite difference to solve and continue solutions with a given D m dihedral symmetry.However, determining the nonlinearities of each model projected onto angular Fourier modes remains a difficult task, and each bifurcation curve would be restricted to the dihedral symmetry group D m , thus preventing symmetry-breaking instabilities.
One could instead use a PDE solver such as pde2path [62], but the large domains required for vegetation models would make any continuation study computationally intensive.It would be interesting to consider the bifurcation structure of these localised patterns in a future numerical study.Furthermore, it would also be interesting to explore the sensitivity of our localised vegetation patterns to the domain boundary.The size of our spatial domain and choice of boundary conditions can both cause nonlinear effects in the time evolution of vegetation patterns; for localised patterns in sufficiently large domains these effects are negligible, but become more important as the vegetation patches grow in width.We again leave this topic for future study.
A possible extension of this work would be to consider nonlocal equations or reaction-diffusion systems with more than 2 components; there are several examples of vegetation models of these forms, including in [20] and [51].However, such problems would require some kind of model reduction, either formally or rigorously, before one could apply the techniques presented in this work.It would also be interesting to try and investigate more complicated two-dimensional localised patterns, such as localised labyrinthine patterns observed numerically [11].Of course, the theory presented in this work goes beyond the study of vegetation patterns, and could be applied to any two-component reaction-diffusion system that undergoes a Turing bifurcation.Some interesting examples would include the Lugiato-Lefever equation in nonlinear optics, and the Gross-Pitaevskii (GP) equation for Bose-Einstein solitons in quantum mechanics.
This work also highlights the need for more topological measures in order to compare localised planar patterns.We discussed the qualitative differences between our numerical simulations, but we are not currently equipped with the tools to easily compare them quantitatively.There have been recent developments in this direction for studying nearly hexagonal lattices [49,60] which could provide further insight into the dynamics observed in our simulations.so that .
Then, we see that , and so P 3 < 0 for all parameter values, since m − 2k 2 = 1 δv (1 + u 2 * ) > 0. Finally, we obtain It is not clear for each m > 0, δ v > 2 m whether P 4 > 0 or P 4 < 0, and so we numerically plot the sign of P 4 in Figure 4(a).

Figure 2 :
Figure 2: (a) Possible bifurcation curves for localised dihedral patterns.The sign of P 1 determines the direction of bifurcation for spot A-type patterns (black, solid) and, if P 4 < 0, ring-type patterns (red, dashed).(b) The different profiles of localised patterns; P 2 determines the phase of the solutions, while P 3 determines the polarity.Here, u denotes vegetation and v denotes water density.

(iii) P 3
determines the polarity of spot A-type patterns; see Figure 2(b).If P 3 > 0 then spot A-type patterns are made up of peaks and if P 3 < 0 then spot A-type patterns are made up of gaps.(iv)P 4 determines whether or not ring-type patterns exist.If P 4 < 0 then ring-type patterns bifurcate from the Turing point, as seen in Figure2(a), and if P 4 > 0 then no ring-type patterns emerge.
10) into a canonical form for a Turing bifurcation.The vectors Û * 0 , Û * 1 form a dual basis for Û0 , Û1 , allowing us to project (2.10) onto the respective coordinates A, B.

Figure 4 :
Figure 4: (a) Plot of the sign of P 4 in (δ v , m)-parameter space; P 4 < 0 in the left-most shaded region (blue) and P 4 > 0 in the right-most shaded region (green).The left-most boundary is the curve δ v m = 2, and so no Turing bifurcation occurs in the unshaded region.The red dot represents the parameter values (2.25) chosen for our numerical study.(b) Schematic bifurcation diagram of the Klausmeier-Gray-Scott model (2.24) with parameter values (2.25), where anti-phase spot A-type patterns emerge as gaps from the Turing point as µ increases.

Figure 5 :
Figure 5: Schematic bifurcation diagrams of (a) the logistic Klausmeier model (2.26) with parameter values (2.28), and (b) the NFC-Gilad model (2.29) with parameter values (2.31).In both cases, anti-phase spot A-type patterns emerge as gaps from the Turing point as µ increases.

2 Figure 6 :
Figure 6: Schematic bifurcation diagram of the von Hardenberg model (2.32) with parameter values (2.34), with in-phase spot A-type patterns emerging as peaks from the lower Turing point (2.35) and as gaps from the upper Turing point (2.36).

Figure 8 :
Figure 8: Time evolution of vegetation density u in the Klausmeier-Gray-Scott model (2.24) with parameter values (2.25).The dark green represents a vegetated state ( u ≈ 1.5), while the lighter yellow represents a bare state ( u ≈ 0).Each row corresponds to a different initial perturbation (3.6) for-from top to bottomhexagons, squares, and pentagons.From left to right, the columns correspond to different points in time t = 100, 200, 300, 400 and 500.The final column is divided into vegetation density (left) and water density (right), where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at[26].

Figure 9 :
Figure 9: Time evolution of vegetation density u in the logistic Klausmeier model (2.26) with parameter values(2.28).The dark green represents a vegetated state ( u ≈ 0.6), while the lighter yellow represents a bare state ( u ≈ 0).The rows and columns are identical to Figure8.The final column is divided into vegetation density (left) and water density (right), where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at[26].

Figure 10 :
Figure 10: Time evolution of vegetation density u in the NFC-Gilad model (2.29) with parameter values (2.31).The dark green represents a vegetated state ( u ≈ 0.6), while the lighter yellow represents a bare state ( u ≈ 0).The rows and columns are identical to Figure 8.The final column is divided into vegetation density (left) and water density (right), where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at[26].

Figure 11 :
Figure 11: Time evolution of vegetation density u in the von Hardenberg model (2.32) with parameter values (2.34) near its first Turing point.The dark green represents a vegetated state ( u ≈ 0.35), while the lighter yellow represents a bare state ( u ≈ 0).The rows are the same as Figure8, the columns correspond to the points in time t = 200, 300, 400, 500 and 600.The final column is divided into vegetation density (left) and water density (right), where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at[26].

Figure 12 :
Figure 12: Time evolution of vegetation density u in the von Hardenberg model (2.32) with parameter values (2.34) near its second Turing point.The dark green represents a vegetated state ( u ≈ 0.35), while the lighter yellow represents a bare state ( u ≈ 0).The rows and columns are the same as Figure 11.The final column is divided into vegetation density (left) and water density (right), where darker and lighter blue represents higher and lower water density, respectively.Videos of each simulation can be found at[26].
The logistic Klausmeier model (2.26) is again covered by the general reaction-diffusion model (2.2) with Hence, the logistic Klausmeier model (2.26) with parameter values (2.28) exhibits anti-phase spot A-type localised gaps, which bifurcate as the precipitation increases beyond the Turing point.Comparing Figure5(a)with Figure4(b), we see that the localised structures emerging in the Klausmeier-Gray-Scott model(2.24)and the logistic Klausmeier model (2.26) are qualitatively the same.

Table 1 :
The Klausmeier-Gray-Scott (2.24), NFC-Gilad (2.29) and von Hardenberg models (2.32) all exhibit similar growth behaviour inFigures 8, 10 & 12, but with significantly different final patterns.For the Klausmeier-Gray-Scott model in Figure8there is no clear packing structure as each pattern grows larger, whereas solutions the NFC-Gilad model in Figure10exhibit nonuniform hexagonal packing and solutions for the von Hardenberg model in Figure12exhibit uniform hexagonal packing.We note that we cannot reasonably compare the speed of growth of the patterns in each model, as each model considered in this work is nondimensionalised with a rescaled time coordinate.For example, for the fixed parameter values given, t = 500 in our simulations corresponds to approximately 50 years in the NFC-Gilad model but 125 years in the logistic Klausmeier model.We have provided novel tools for the design of phenomenological reaction-diffusion models; by computing the quantities P 1 , P 2 , P 3 , P 4 , one is able to quickly check the qualitative behaviour of solutions of a prospective model and compare them against the observed phenomena.While these predictors are derived for localised patterns, they also have applications in general pattern formation.Like spot A-type patterns, domaincovering hexagons also possess a fixed phase and polarity-as predicted by P 2 and P 3 , respectively-while P 4 determines whether unstable stripes bifurcate in the direction determined by P 1 .In order to model the emergence of a particular type of pattern, one would then require a particular type of Turing point determined by the signs of each P i ; this is summarised in Table1.A summary of the different qualitative patterns emerging from a Turing point and their required values for P 2 , P 3 and P 4 .