|Home | About | Journals | Submit | Contact Us | Français|
The development of spontaneous stationary vegetative patterns in an arid isotropic homogeneous environment is investigated by means of various weakly nonlinear stability analyses applied to the appropriate governing equation for this phenomenon. In particular, that process can be represented by a fourth-order partial differential time-evolution logistic equation for the total plant biomass per unit area divided by the carrying capacity of its territory and defined on an unbounded flat spatial domain. Those patterns that consist of parallel stripes, labyrinth-like mazes, rhombic arrays of rectangular patches, and hexagonal distributions of spots or gaps are generated by the balance between the effects of short-range facilitation and long-range competition. Then those theoretical predictions are compared with both relevant observational evidence and existing numerical simulations as well as placed in the context of the results from some recent nonlinear pattern formation studies.
Rayleigh-Bénard buoyancy-driven convection has to date provided perhaps the best-studied example of nonlinear pattern selection (reviewed by Koschmieder ). One of the methods traditionally used to predict such pattern selection is a weakly nonlinear stability analysis that, although incorporating nonlinearities of the relevant model system, basically pivots a perturbation procedure about the critical point of linear stability theory (reviewed by Wollkind et al. ). The advantage of such an approach over strictly numerical procedures is that it allows one to deduce quantitative relationships between system parameters and stable patterns, which are valuable for experimental design and difficult to accomplish by using simulation alone. Recently, there has been considerable interest generated in pattern formation and selection during normal-incidence ion-sputtered erosion of metallic or semiconductor solid surfaces and during diffraction of radiation induced by the injection of a laser pump field into a ring cavity containing a purely absorptive two-level atomic sodium vapour medium. In order to predict the sequence of solid surface morphologies and nonlinear optical patterns actually observed during such ion-sputtered erosion and laser-induced diffraction, respectively, Pansuwan et al.  and Wollkind et al.  performed both rhombic- and hexagonal-planform weakly nonlinear stability analyses on the governing damped Kuramoto-Sivashinsky and modified Swift-Hohenberg equations, respectively, appropriate for modelling these phenomena.
We wish to continue our examination of nonlinear phenomena by investigating the development of spontaneous stationary equilibrium vegetative patterns in an arid isotropic homogeneous environment (reviewed by Rietkerk et al. ). Pattern formation in this phenomenon can be modelled for the isotropic case by a partial differential time-evolution logistic equation describing the total plant biomass per unit area divided by its carrying capacity that under a weak gradient-low density truncation is reduced to the fourth order in its spatial variables defined on an unbounded two-dimensional flat domain . Under appropriate conditions, those patterns that consist of parallel stripes, labyrinth-like mazes, rhombic arrays of rectangular patches, and hexagonal distributions of spots or gaps can be generated in an initially homogeneous environment by the balance between the effects of short-range facilitation and long-range competition. We shall perform these same rhombic- and hexagonal-planform weakly nonlinear stability analyses on that truncated fourth-order partial differential time-evolution logistic equation, and then compare the theoretical predictions obtained with both relevant observational evidence and existing numerical simulations as well as place them in the context of some recent pattern formation studies. We begin below with a brief sketch of the reduction procedure required to obtain the governing truncated fourth-order partial-differential time-evolution logistic equation.
in order to model spatio-temporal patterning for botanical ecosystems. Here ρ(x, y, t) is the total plant biomass per unit area divided by its carrying capacity where (x, y) a two-dimensional spatial coordinate system and t the time while F1 · F2 and F3 represent a nonlocal birth term and a local death term, respectively. In particular, F1, which represents the actual propagation distribution function and involves the growth effects of reproduction, seed dissemination, germination, and maturation of young plants, depends mainly on vegetative aerial structures like the canopy in the case of trees or bushes, while F2, which represents the inhibition distribution function and involves effects that inhibit the development of vegetation such as root systems and the uptake of scarce resources like water and nutrients, depends mainly on the availability of free soil for vegetative colonization and on soil vegetative resources of limited supply. Specifically, for the isotropic case employing Gaussian weighting functions in conjunction with a Taylor series expansion for ρ, these distribution functions take the following form :
where the gradient operator (/x, /y), 0 1, and 2 · while L and Λ are the facilitation-to-competition range and interaction ratios, respectively. Finally,
where the mortality-to-growth rate ratio μ serves as a measure of the environment's aridity. We note that in Equation (1) time and space have been scaled with the inverse of the vegetation growth rate, which corresponds roughly to the period spent by the plants to achieve adult size and with the interplant competition range, which can be approximated by the radius of the superficial root system of the dominant plants, respectively.
It is now possible to truncate the partial derivatives appearing in Equation (1) by introducing scaling laws such that
and hence to third order
Observe that this development is consistent with the scaling laws
for which each term in Equation (3a) is of O(ε3). Here, the parameter ε can be identified with the mean biomass density α, to be defined below in Equation (8), which is relatively small for arid or semi-arid environments. These scaling laws differ from those introduced by Lefever et al. , and that discrepancy will be discussed in more detail in the last section of this paper. Therefore as observed by Lejeune and Tlidi , Couteron and Lejeune , Lejeune et al. , and Lejeune et al. , Equation (3) has been derived in a low density and weak gradient limit from a generalized logistic equation describing nonlocal interactions. This truncated equation applies to water-poor ecosystems where average plant phytomass density is low with respect to carrying capacity (close-packing density) of unstressed vegetation and for which pattern wavelengths are large in comparison with the average size of the dominant plant form but very small when compared with the territorial length-scale characteristic of that arid environment. Since under this latter condition, the actual territorial boundary of that arid environment does not significantly influence the patterns , it then seems reasonable as a first approximation for us to consider Equation (3a) on an unbounded two-dimensional flat spatial domain.
Rewriting that equation in the form
we first seek a uniform stationary solution, ρ = ρ0, of it satisfying
and find that
Note that ρ0 = 0 which represents bare ground always exists while when ρ±0 are real and positive, these other two solutions represent spatially homogeneous plant distributions. If 0 ≤ Λ, μ ≤ 1 then ρ+0 ≥ 0, ρ−0 ≤ 0 and hence only ρ+0 is biologically meaningful. If Λ > 1, then ρ+0 > 0 for 0 ≤ μ ≤ μ = 1 + (Λ − 1)2/4 while ρ−0 > 0 for 1 ≤ μ ≤ μ. Performing a linear stability analysis of these distributions to homogeneous disturbances by considering a solution of (4) of the form
yields the perturbation equation
and we have linear stability in Equation (6) when f′ (ρ0) < 0 and instability when f′ (ρ0) > 0, it can be concluded from Equation (7) that ρ0 = 0 is linearly stable for μ > 1 and unstable for 0 ≤ μ < 1, ρ0 = ρ+0 is linearly stable everywhere it exists, and ρ0 = ρ−0 is linearly unstable everywhere it exists.
We are interested in determining conditions under which nonzero homogeneous distributions of vegetation that are linearly stable in the absence of spatial effects can become unstable to heterogeneous perturbations. This is reminiscent of the diffusive instabilities that occur during chemical Turing pattern formation . Thus we need only consider the solution ρ+0 which will now be designated by
In order to eliminate the possibility of bistability between that homogeneous solution and the zero state, in what follows after Lejeune et al.  we shall restrict our attention to
In this context, we adopt the far-field condition that
which for the special representative case for Λ = 1 reduces to
It is the weakly nonlinear stability of this homogeneous solution to one-dimensional and two-dimensional rhombic- or hexagonal-planform heterogeneous perturbations with which we shall be concerned in Sections 2–4, respectively.
As a necessary prelude to our two-dimensional rhombic- and hexagonal-planform nonlinear stability analyses of the truncated fourth-order logistic equation (4) to be presented in the next two sections, we perform a weakly nonlinear stability analysis of its homogeneous state in this section by seeking a real one-dimensional Stuart –Watson  type solution of it through third-order terms of the form 
where the amplitude function A1(t) satisfies the Landau equation
and q = 2π/λ, λ being the wavelength of the class of spatially periodic perturbations under investigation. Substituting the solution of Equation (12) into Equation (4), we obtain a sequence of problems, one for each pair of values m and n which corresponds to a term of the form Λn1 (t) cos(mqx) appearing in Equation (12a). Then the n = m = 0 problem is satisfied identically by virtue of the fact that
while the linear stability problem for n = m = 1 yields the secular equation
which is a parabola in the q2 − σ plane. Note that if β ≥ α, then σ < 0 for the parameter range of Equation (9) while if β < α then this parabola has a maximum value at its vertex (q2c, σc) with
where the discriminant of the quadratic in q2 of Equation (14a) is given by
If, in addition, D(α, β) > 0, 0 < σ < σc(α, β) for q21 < q2 < q22 where
This scenario is plotted in Figure 1 for the special case of Λ = 1 when α = 0.2 and β = 0.02. Hence, in what follows we shall equate the q and σ of Equation (12) to
where βc = βc(α) is defined implicitly by D(α, βc) = 0 or explicitly by
Therefore, the homogeneous solution is linearly stable for β > βc(α), neutrally stable for β = βc(α), and unstable for 0 < β < βc(α). Thus β = βc(α) serves as the marginal stability curve in the α − β plane. This situation is plotted in Figure 2 for the special case of Λ = 1. In that instance
while, since βc′;(α) = 1 − 3α1/2, βc′(α0) = 0 implies
Here we are restricting our analysis to the critical wavenumber alone. In Section 5, we shall discuss the results of the so-called side-band instabilities analysis originally introduced contemporaneously by Segel  and Newell and Whitehead  that examines the consequences of investigating other wavenumbers in the side band centred about that critical wavenumber as well.
A plot in the α − β plane of the marginal stability curve β = βc (α) of Equation (16) for Λ = 1 on which the growth rate σc(α, β) = 0 denoted by the solid line. Here the dashed ...
Continuing our description of the results of this one-dimensional expansion procedure, the second-order problems corresponding to n = 2 and m = 0 or 2 can be solved in a straightforward manner to yield
Although there are also two third-order problems, it is permissible to concentrate our attention exclusively on the one corresponding to n = 3 and m = 1 which contains the Landau constant a1 for the Fredholm-type method of solvability we shall use. That problem may be represented as
Here we have only denoted the β dependence of the quantities in question for ease of exposition. Now taking the limit of Equation (18) as β → βc = βc(α), employing Equations (14b), (14c), (15), and (17), and assuming the requisite continuity at β = βc, we obtain the relevant solvability condition represented symbolically as
or given explicitly by
A main feature of the weakly nonlinear stability theory implicit to the formulation of our expansion of Equation (12) with q qc = qc(α, β) is a phenomenological interpretation of the problem under examination based upon the long-time behaviour of the solution A1 (t) to the truncated amplitude Equation (12b). In general, the dynamics of such a Landau equation with real coefficients can be divided into the four qualitatively different cases represented by the possibility of σ, its growth rate, and a1, its Landau constant, being either positive or negative. These cases can be catalogued as follows:
From the description of the four cases just summarized, we can see that the stability behaviour of the Landau equation is crucially dependent upon the sign of a1. Hence in order to determine this behaviour, it is necessary that we next examine the formula for a1 given by Equation (18d). Towards that end we plot a1 of Equation (18d) versus α in Figure 3 with Λ = 1. From this figure, we observe that a1 > 0 whenever
In order to investigate the possibility of occurrence for the truncated fourth-order logistic equation of the type of two-dimensional vegetative patterns mentioned earlier, we shall first consider a rhombic planform solution of Equation (4) of the form [19,30]
Here we are employing the notation ρnjmk for the coefficient of each term in Equation (22a)–(22b) of the form An1(t)Bj1(t) cos[qc(mx + kz)] where ρ1010 = ρ0101 = 1, implicitly. Then substituting this rhombic planform solution of Equation (22) into Equation (4), we again obtain a sequence of problems, each of which corresponds to one of these terms. Solving those problems we find that
as defined in the previous section while
and b1, in particular, satisfies
Having developed these formulae for its growth rate and Landau constants, we turn our attention to the rhombic planform amplitude Equation (23c), which possess the following equivalence classes of critical points (A0, B0):
Assuming that a1, a1 + b1 > 0 and investigating the stability of the critical points of Equation (24) by seeking a solution of our amplitude equations of the form
one finds that p has the associated roots 
which yield the stability criteria that I is stable for σ < 0; II, for σ > 0, b1 > a1; and V, for σ > 0, a1 > b1. Note that I and II, as in the one-dimensional analysis of the previous section, represent the uniform homogeneous and supercritical striped states, respectively, while V can be identified with a rhombic pattern possessing the characteristic angle ϕ (see below and ). We next use these criteria to refine our one-dimensional predictions relevant to the former states due to the presence of the latter. Towards that end, we examine the signs of a1 ± b1 for Λ = 1, α (α1, α2), and ϕ (0, π/2) with ϕ = π/2 (or equivalently 90°) representing a square planform. We first illustrate this procedure by defining the ratio of Landau constants 
and plotting that quantity versus α in Figure 6 for a fixed value of ϕ, namely ϕ = 90°. Here, there exist two α-intervals of stable square rhombic patterns, where a1 ± b1 > 0 or equivalently − 1 < γ < 1, given by
provided in addition that σ > 0 or 0 < β < βc(α). These intervals have been indicated by shading in this figure and observe that for α values between them γ > 1or b1 > a1 > 0. Hence for α (0.1176, 0.2094) and 0 < β < βc(α) there exist stable supercritical striped patterns. Now, repeating the process used to produce Figure 6 but for other values of ϕ, we find the same generic behaviour as for ϕ = π/2 in that there are two intervals of stable rhombic patterns for α (αm, αl) or α (αr, αM) with stable stripes between them for α (αl, αr) when 0 < β < βc(α) where
and summarize these results in Table 1.
α-range for stable rhombic patterns versus ϕ.
which implies that only stripes can be stable with respect to rhombic patterns of such small characteristic angle. Next, we plot our ratio of Landau constants versus ϕ [0, π] for the fixed values of Λ = 1 and α = 0.12 in Figure 7. Restricting ourselves to the interval of interest ϕ [0, π/2], we see that there are two bands of stable rhombic patterns flanking ϕ = π/3 for ϕ (ϕm, ϕl) or ϕ (ϕr, ϕM) when 0 < β < βc(α) where
which have been designated by shading, with no stable patterns between these bands for ϕ (ϕl, ϕr) and stable stripes outside of them for ϕ [0, ϕm) or ϕ (ϕM, π/2] when 0 < β < β c(α). This figure has been drawn for the extended interval ϕ [π/2, π] in order to demonstrate graphically the symmetry about ϕ = π/2 characteristic of rhombic patterns since
and thus, for α = αc = 0.16 and 0 < β < βc = 0.032, stripes are stable versus rhombic patterns. Again, repeating the process used to produce Figure 7 but for other α ≠ αc, we find the same generic behaviour as for α = 0.12 and summarize these results for selected values in Table 2.
Angle range for stable rhombic patterns versus α.
The occurrence of this vertical asymptote is a consequence of b1 of Equation (29) possessing a negative numerator and a denominator having a double root at ϕ = π/3.
We close this section with a morphological instability interpretation of the potentially stable critical points II and V of our amplitude equations when σ > 0 or 0 < β < βc (α) relative to the vegetative patterns under investigation. Then, to the lowest order the deviation from its homogeneous configuration of the function associated with these critical points is given by
The contour plot in the x − y plane for this deviation function with A0 > 0 and B0 = 0 relevant to critical point II is identical to that in Figure 4, of the previous section. Here, as indicated earlier, regions of higher density appear dark and those of lower density, light. Clearly, such alternating dark and light parallel bands should be identified with a striped vegetative pattern as anticipated above. In order to make an analogous interpretation of critical point V, we consider the deviation function (34) with B0 = A0 > 0 and allow ϕ to take on the values 90° and 58°, sequentially. We represent the contour plot for that function with ϕ = 90° in Figure 8. From the checkerboard appearance of the array, it is equally clear that this critical point should be identified with a vegetative pattern of square planform. Finally, in a similar manner, we generate the contour plot of Figure 9 for 58°, which forms a family of rectangles. Here, ϕ = 58° (or equivalently 1.01) is a representative rhombic angle associated with stable patterns for α = 0.15 (see Table 2). To demonstrate that the same thing happens with any ϕ (0, π/2], we first put Equation (34) for critical point V in the form
From Equations (35) and (36), we can then deduce that the intersecting level curves in the associated contour plots are two families of straight lines possessing slopes of
respectively, which intersect at right angles since
Given that, in general, this is a rectangular planform we need to explain in what sense such a critical point can be identified with a rhombic pattern. To do so, we form the quadrilateral depicted by solid lines in Figure 9. Its sides are each composed of two half-diagonals, collectively contained in the four light rectangles surrounding a dark one. Thus, each side of the quadrilateral has the length of one of these diagonals
and hence, the quadrilateral is a rhombus of characteristic angle ϕ. Further, ϕ also plays a role in characterizing the family of rectangles. Each member of that family has aspect ratio
where w and l are its width and length, respectively, while ϕ/2 serves as its angle of inclination as well, an assertion most easily verified by the relation
Therefore, in what follows we shall refer to such a vegetative pattern as a rhombic array of rectangular patches.
The weakly nonlinear stability behaviour of the amplitude-phase equations (39c), (39d) to be described below depends only on the values of their growth rate and Landau constants. We can determine the solvability conditions for these Landau constants most easily by introducing the transformation
which reduces Equation (39) to
Now we are employing the notation ρnjMk for the coefficient of each higher-order term in Equation (42) of the form Λn1(t)Bj1(t) cos(Mqcx/2) cos(k√3qcy2). Then substituting the solution of Equations (39)–(42) into Equation (4) and solving the resultant sequence of problems, we find specifically that
and, in particular, the other two Landau constants satisfy
Note from Equations (44a) and (44b) that the expression a2(α) of Equation (44d) does not contain the component limβ→βc(α) ρ1111 since the coefficient of that quantity, namely limβ→βc(α) (4α0 + Λ − 1 − 3α + q2c/2 − q4c/8), is identically equal to zero by virtue of Equation (44a). As pointed out by Wollkind et al.  such independence may be expected as a general rule, thus eliminating the necessity of determining ρ1111.
Having developed formulae for their growth rate and Landau constants, we return to the six-disturbance hexagonal-planform amplitude-phase equations (39c) and (39d). In cataloguing the critical points of these equations and summarizing their orbital stability behaviour, it is necessary to employ the quantities
Then, those critical points can be represented as equivalence classes of the form (A0, B0, B0, 0, 0, 0) corresponding to ϕ1 = ϕ2 = ϕ3 = 0 with
Using similar but slightly more complicated considerations than those employed in the last section, one can deduce orbital stability criteria for these critical points that are posed in terms of σ . Thus critical point I is stable in this sense for σ < 0 while the stability behaviour of II and III±, which depends on the signs of a0 and 2a2 − a1 as well, has been summarized in Table 3 under the assumption that a1 + a2 > 0.
Orbital stability behaviour of critical points II and III±.
In this parameter range, A+0 > 0 and A−0 < 0. Finally, critical point IV, which reduces to II for σ = σ1 and to III± for σ = σ2 and is hence called a generalized cell , is not stable for any value of σ. Here, we use the term orbital stability of pattern formation to mean a family of solutions in the plane that may interchange with each other but do not grow or decay into a solution type from a different family. Such an interpretation depends upon the translational and rotational symmetries inherent to the amplitude-phase equations, these invariancies also limiting each equivalence class of critical points to a single member that must be considered explicitly. The same reasoning was implicit to our cataloguing of the critical points of the rhombic planform amplitude equations in the previous section.
We next offer a morphological interpretation of the potentially stable critical points catalogued above relative to the vegetative patterns under investigation. Then, critical points I and II represent the uniform homogeneous and supercritical striped states, respectively, described in both of the last two sections. Observing that to lowest order, the deviation function associated with critical points III± satisfies
and noting that the function g0(x, y) exhibits hexagonal symmetry , we can deduce that these critical points represent hexagonal lattices possessing individual hexagons with dark circular regions at their centres and light boundaries for III+ and with light circular regions at their centres and dark boundaries for III−. The contour plots relevant to these critical points are depicted in Figures 10 and and1111 where Figure 10 is for III+ and Figure 11, for III−. Hence, in what follows, we shall identify hexagonal close-packed distributions of vegetative spots or gaps with critical points III+ or III−, respectively.
Having summarized those stability criteria and morphological identifications, we now return to our expressions for the Landau constants of Equations (39c) and (39d) given by Equations (18d), (44c), and (44d). First, we examine the signs of a1 + 4a2, a1 + a2, a0, and 2a2 − a1 versus α in Figures 12 and and1313 for Λ = 1. From the results of this examination, we observe that besides α1, α2, and αc defined in Equation (19) and Equation (28b) there exist the following other significant values of α
These values of α are compiled in Table 4 for Λ = 0.825 and 1.
The α values of Equation (49).
In order to compare our theoretical predictions with numerical simulations and relevant observational evidence, we must represent the stability results of Table 3 graphically in the α − β plane. To do so, it is now necessary for us to generalize the approach which produced the marginal stability σc = 0 curve β = βc(α) of Equation (15) from Equation (14) in order that we may generate the analogous loci associated with σ = σj for j = −1, 1, and 2, respectively. That is, employing Equations (14), (18d), (44c), (44d), and (45), we solve
for β to obtain
respectively. Since all the quantities required for the identification of the vegetative patterns of Table 3 have been evaluated, we can represent graphically the regions corresponding to these patterns in the α − β plane of Figure 14 with Λ = 1, where the loci of Equations (15) and (50b) are denoted by σc and σj with j = −1, 1, and 2, respectively, in that figure which has been plotted for 0.0600 ≤ α ≤ 0.2108. Then we observe from Figure 14 in conjunction with Tables 3 and and44 that for a1 + 4a2 > 0 all (when 2a2 − a1 < 0) or part (when 2a2 − a1 > 0) of the region (σc, a1 > 0), where the one-dimensional analysis of Section 2 predicted striped vegetation patterns, is further divided into two subregions characterized by hexagonal patterns consisting of either vegetative spots (when a0 < 0) or gaps (when a0 > 0), respectively. In the overlap regions satisfying σ1 < σ < σ2 or
where stripes and spots (α−1 < α < αc) or stripes and gaps (αc < α < α2) are predicted, the initial conditions determine which stable equilibrium pattern of each pair will be selected. In this context, we define
which from Figure 4 implies
There also exists a region of bistability corresponding to σ−1 < σ < 0 or
the homogeneous distribution being stable for σ < 0 or β > βc (α) and hexagons for σ−1 < σ < σ2. Given that σ−1 < 0 for a0 ≠ 0 or α ≠ αc = 0.16, the hexagons persisting in this overlap region would be subcritical in nature. Finally, to justify the truncation procedure inherent to the asymptotic representation of Equations (39c) and (39d), it is necessary that its coefficients satisfy the size constraint 
which, as can be seen in Figure 15, is valid in the parameter range of interest. Hence, we may conclude that such a truncation procedure is justified for our hexagonal planform weakly nonlinear stability analysis of Equation (4). Further, given the hexagonal close-packed nature of the distributions associated with III±, we have referred to them collectively as hexagons above.
We are now ready to compare these theoretical predictions developed in Sections 3 and 4 with relevant numerical simulations and observational evidence as well as place them in the context of some recent pattern formation studies. We begin naturally by comparing our predictions with the numerical simulations of Lejeune and his co-workers: namely with those contained in Lejeune and Tlidi , Couteron and Lejeune , Lejeune et al. , and Lejeune et al. . In particular, Lejeune et al.  numerically integrated Equation (3) on a square grid with periodic boundary conditions for Λ = 0.825 and L = 0.1 after constructing the corresponding bifurcation diagram of associated vegetative states as a function of μ to select appropriate values of this aridity control parameter favouring hexagonal pattern formation. These authors reported the following sequence of spatially periodic vegetation states obtained from that bifurcation diagram as aridity was increased: first, a pattern of spots of sparser vegetation or gaps appeared forming a hexagonal lattice, which then transformed into an alternation of stripes of sparser and thicker vegetation, and finally into a hexagonal pattern of vegetation spots separated by bare ground. Recalling that β = L2 and α is inversely related to μ by Equation (11b), we can see that the sequence of morphologies obtained from Figure 14 when the horizontal line
is traversed in the direction of increasing α through the patterned region
is in qualitative agreement with the sequence just catalogued. Here after Borckmans et al. , we employ the notation B/C to indicate bistability between patterns B and C. Next, Lejeune et al.  performed the numerical simulations for the fixed values of Λ and L as described above and generated vegetation patterns constituted of hexagonal gaps, alternating stripes, or hexagonal spots when μ was assigned the values of 0.965, 0.970, or 0.980, sequentially. Although Λ = 1 chosen in Sections 2–4 for ease of exposition is a special case since such a choice eliminates the quadratic term of f (ρ), it is a representative value of that parameter for our purposes. The qualitative behaviour of our stability analyses is the same for the other values of 0 < Λ < 1 defined by Equation (9) as it was for Λ = 1. This is demonstrated conclusively relevant to the simulation results of  by Table 4, which lists the special α values of Equation (49) for Λ = 0.825, which was the value employed during that simulation, as well as those for Λ = 1. To compare our theoretical predictions with these simulation results, we must first construct a stability diagram analogous to Figure 14 but for Λ = 0.825 rather than Λ = 1. Figure 16 represents such a diagram plotted in the α − β plane. Note in this context that αc as defined implicitly by Equation (49e) has the explicit formulation as a function of Λ given by
This critical value of density has been denoted by the vertical line α = αc in both these figures and separates the two types of vegetative hexagonal patterns from each other since spots can only occur for α < αc while gaps can only occur for α > αc. Then using Equation (11a) with Λ = 0.825, it is possible for us to convert Figure 16 into Figure 17, which represents the corresponding stability diagram plotted in the μ − β plane. Again note that the vertical line in this figure is given by
where the latter quantity corresponds to the αc of Equation (57) when Λ = 0.825. Finally, in order to compare the simulation results just described with our theoretical predictions, we now select β = L2 = 0.01 consistent with the value of L = 0.1 employed by Lejeune et al. . Upon examination of Figure 17, we can deduce the predicted morphological sequence of vegetative patterns summarized in Table 5 as the horizontal transit curve β = 0.01 is traversed in the direction of increasing μ. Observe from Table 5 that the complete patterned region occurs for the interval μ (μc−, μc+) where
are the two critical points of bifurcation as identified by Lejeune et al.  while the three μ values of 0.965, 0.970, or 0.980 chosen by the latter authors for their simulation purposes lie within the subintervals for which only gaps, only stripes, or only spots are stable, respectively. Hence our theoretical predictions are in quantitative agreement with these simulation results. Further, note from Figures 16 and and1717 that the term homogeneous appearing in the stable pattern column of the first two entries of Table 5 refers to uniform vegetative distributions of relatively high density while the same term appearing in the corresponding column of its last two entries refers to similar distributions of relatively low density. The results of our weakly nonlinear analyses are only strictly valid in the vicinity of the marginal stability curve. Often, however, in situations of this sort it is possible to extend these results to regions of the relevant parameter space substantially removed from the marginal curve. In order to determine the validity of such an extrapolation, it is standard operating procedure to perform numerical simulations of the original governing partial differential equation(s) in those extended regions and compare these simulation results with those theoretical predictions. This is precisely what was done by Firth and Scroggie  in their nonlinear optical pattern formation problem involving an atomic sodium vapour ring cavity and more recently by Golovin et al.  in their chemical Turing pattern formation problem of the Brusselator model with superdiffusion. In the latter instance, the theoretical predictions given in their Figure 3 were compared with the corresponding numerical simulations appearing in their Figure 5. Indeed Lejeune et al.  also compared their theoretical arid environmental vegetative pattern formation predictions and numerical simulations by means of their Table II which included the maximum and minimum phytomass density values for the simulated patterns of their Figure 5 and the corresponding values for the analytical bifurcation diagram given in their Figure 4 that was a plot of the stationary states of their amplitude equations versus μ. We reproduce those results in our Table 6. Considering the amplitude of the difference of these maximum and minimum values the table shows the relative deviation between the analytical and numerical values to be less than 10% in each of those cases, as pointed out by Lejeune et al. . These cases, summarized earlier upon the introduction of Table 5, were for Λ = 0.825 and β = 0.01 with the aridity parameter μ = 0.965, 0.970, and 0.980 corresponding to gaps, stripes, and spots, respectively, all of which, as can be seen from Figure 17, lie relatively deep within the patterned region.
Morphological stability predictions versus μ for β = 0.01 and Λ = 0.825.
So far, we have been concerned with the bifurcation behaviour of Equation (4) when 0 < Λ ≤ 1 only for β in the range defined by Equation (55). We shall now examine what transpires should we consider Figure 14 in conjunction with the horizontal transit curve
Here βcrit = 0.015 corresponds to the local minimum value of the curve β = β1 (α) for α > αc. Then we obtain the sequence of morphologies
as α is increased along the transit curve of Equation (60). The simulation results of Lejeune and Tlidi  and Couteron and Lejeune  are in a parameter range that yields bifurcation behaviour of this generic sort. In particular, these authors numerically integrated Equation (3) on a square grid with periodic boundary conditions for Λ = 1.2 and L = 0.2 obtaining a hexagonal pattern of gaps for μ = 0.98, alternating stripes for μ = 1.0, and a hexagonal pattern of spots for μ = 1.01. To compare our theoretical predictions with these simulation results, we again construct stability diagrams analogous to Figures 16 and and1717 but for Λ = 1.2 in both the α − β plane of Figure 18 and the μ − β plane of Figure 19. In doing so, we first point out that this requires a relaxation of our original constraint of 0 ≤ Λ, μ ≤ 1 and second, that neither Lejeune and Tlidi  nor Couteron and Lejeune  constructed a bifurcation diagram relevant to their two-dimensional simulations although Lejeune et al.  did generate a so-called one-dimensional bifurcation diagram for this case (see below). Observe that Figures 16 and and1717 generated for Λ = 0.825 are in qualitative agreement with Figures 18 and and19,19, respectively, generated for Λ = 1.2. We also note that the vertical lines in those figures are given by
respectively. Next we select β = L2 = 0.04 consistent with the value of L = 0.2 employed by Lejeune and Tlidi  and Couteron and Lejeune . Then upon examination of Figure 19, we can deduce the predicted morphological sequence of vegetative patterns summarized in Table 7 as the horizontal transit curve β = 0.04 is traversed in the direction of increasing μ.
Morphological stability predictions versus μ for β = 0.04 and Λ = 1.2.
Observe from Table 7 that this morphological sequence is both in qualitative agreement with the sequence of Equation (61) and quantitatively compatible with the simulation results under examination once the μ values in that table have been rounded off to the two places of accuracy reported by Lejeune and Tlidi  and Couteron and Lejeune . Upon examination of Figure 19, we can see that the parameter values of Table 7 lie even deeper within the patterned region of this figure than the parameter values of Table 6 did in Figure 17. In this context, Lejeune et al.  reported their simulation μ values to three places of accuracy. We note that for all the simulation results of Lejeune and his co-workers described in this section which generated striped patterns,
to the relevant degree of accuracy. In order for us to demonstrate that this occurrence was not coincidental, it is necessary that we consider our rhombic planform results of Section 3 in conjunction with the morphological stability predictions summarized in Table 3.
Towards this end, we begin by noting that, although these rhombic results of Section 3 were developed explicitly for the special case of Λ = 1, they are qualitatively valid for any other value of Λ satisfying the condition in Equation (3b) – e.g. specifically for Λ = 0.825 or 1.2. Consistent with Equation (52), we implicitly define α±2 by
From the results of Section 3, we can conclude that for any
all such stripes are unstable with respect to rhombic patterns having characteristic angle
In particular, for α in these intervals rhombic arrays of this sort are difficult to distinguish from hexagonal distributions since the allowable bands of their characteristic angles closely flank π/3 1.047 (see Table 2 and Figure 9). Hence, if we substitute this type of rhombic pattern as the stable morphology for the intervals
then hexagonal or nearly hexagonal arrays can be anticipated for all α ≠ αc over the instability region. In this context, we note that all our pattern-formation figures have used the homogeneous solution value of α as the threshold to trigger the colour change from light to dark. Thus, all spatial regions characterized by ρ ≥ α appear dark and those by ρ < α light. If in the α-intervals of Equation (64e) we adopt the protocol that αc represents this critical threshold instead, then in the high threshold interval where α < αc light regions predominate and hence the rhombic arrays in that interval would be characterized by dark rectangular patches with rounded corners completely surrounded by the light background or pseudo-‘spots’, whereas in the low threshold where α > αc dark regions predominate and hence those arrays in that interval would be characterized by light rectangular patches completely surrounded by the dark background or pseudo-‘gaps’. Hence, after Wollkind and Stephenson  using such a high threshold–low threshold type of argument, it is plausible to assume that these pseudo-‘spots’ are the preferred patterns for α−2 < α < αc while pseudo-‘gaps’ are preferred for αc < α < α+2. Finally, from Figures 14 and and1616–19, in conjunction with the morphological stability results of Section 3, we can conclude that the occurrence of stripes requires α = αc, or equivalently μ = μc since such patterns are only stable for both planforms at this critical value of α or μ and when β < βc(αc). Here where our analysis predicts a stable parallel striped pattern the equivalence class for the critical point designated as II in the last section actually contains the three solutions 
These represent a family of stripes aligned parallel to the y axis as per our original identification plus two similar families of stripes making angles of ± 60° with them, for which stable coexistence with a member of either the original family or one another is impossible . Then for α = αc (or equivalently μ = μc) and β < βc(αc) as initial conditions vary from point to point over a flat environment, such families of stripes can give rise to polygonal arcs the boundaries of which would appear quite random in orientation. Indeed, the simulation results classified as striped patterns by Lejeune and Tlidi  and Couteron and Lejeune  have the appearance of such curved elongated stripes or interconnected bicontinuous patterns in the relevant reproductions contained therein.
So far we have limited our discussion to analyses for which the wavenumber was restricted to the critical wavenumber alone. In order to investigate the consequences of considering other wavenumbers in the instability side band centred about this critical wavenumber, it would be necessary to convert the Landau-type amplitude ordinary differential equations in time to Ginzburg-Landau type partial differential equations by adding the appropriate spatial derivative terms to them. As reviewed in detail by Wollkind et al. , such an analysis, although beyond the scope of our present investigation always yields two additional instabilities besides those discussed above: namely, zig-zag and cross-band relevant to the interaction of oblique and perpendicular modes, respectively. Specifically, it should be pointed out that the stripes simulated by Lejeune et al.  for β = 0.01, Λ = 0.825, and μ = 0.97 were of the zig-zag variety.
Up until now, we have been concerned with a1 > 0 type behaviour for the simulation results of Lejeune and his co-workers when Λ = 1.2 and L = 0.2. Lejeune et al.  constructed a one-dimensional bifurcation diagram corresponding to this problem for Λ = 1.2 and L = 0.2. They found isolated vegetation patches as a simulation outcome occurring outside a patterned aridity interval when μ > μ = 1.01 and interpreted such localized structures to be a spatial compromise between the periodic patchy vegetative and bare stable states. We plot a1 versus α in Figure 20 for Λ = 1.2 relevant to this situation and note that α = 0.1 is equivalent to μ = 1.01. Note that Figure 4 for Λ = 1 is in qualitative agreement with Figure 20 for Λ = 1.2. Here the latter case yields the values α1 = 0.1028 and α2 = 0.2821. Observe from Figure 19 that for the μ range included in Table 7 the transit curve β = 0.04 lies totally within the region where σ, a1 > 0. Should μ be extended just beyond the bounds of that table, however, a1 < 0 although σ > 0 (see Figures 18–20). This behaviour is highly reminiscent of the morphological phase separation instabilities (σ > 0) in thin liquid films analysed by Tian and Wollkind . The main result of their analysis was that a1 > 0 type patterns of this sort consisting of uniform distributions of nanodroplets, bicontinuous ridges, and polygonal cells separated by ultrathin flat films could occur sequentially as mean layer thickness was increased over a particular interval with rupture occurring outside of that interval where a1 < 0. This dewetting-type rupture occurred by isolated drop formation in relatively thin layers, but by isolated hole formation in thicker ones. Analogous behaviour was discovered by Golovin et al.  in their examination of chemical Turing pattern formation for the Brusselator problem with superdiffusion. Given the similarity of behaviour among all these phenomena, we conjecture that outside the bounds of Table 7 where σ > 0 and a1 < 0 there will be desertification characterized by the formation of isolated vegetation patches for μ > 1.0098 and by the formation of isolated patches of bare ground for μ < 0.9769.
It remains for us to compare our theoretical predictions with observational evidence. Striking periodic vegetation patterns covering widespread areas of arid or semi-arid regions of Africa, Australia, the Americas, and the Near East became noticeable through aerial photography nearly 70 years ago [2,13,21]. In this instance, an arid region refers to an environment that is characterized by an extended dry season where yearly potential evaporation exceeds yearly rainfall and water availability is a limiting factor of plant growth . The patterns reported consisted of gaps (‘pearled or spotted bush’), labyrinths, stripes (‘tiger bush’), and spots (‘leopard bush’). These included bushy vegetation punctuated by bare spots in Niger, labyrinths of perennial grasses in Israel, striped patterns of bushy vegetation in Niger, and spots of trees or shrubs in Ivory Coast and French Guiana . From our one-dimensional morphological interpretation of the equivalence class for critical point II of Equation (65) it follows that the patterns usually described by ecologists as tiger bush, lanes, groves, and bands should be associated with what we have termed parallel stripes while those described as labyrinths, mazes, and arcs should be associated with what we have termed bicontinuous patterns. Similarly leopard or pearled bush should be associated with critical points or III+ or III−, respectively. We wish to re-examine further these hexagonal patterns in order to determine sufficient conditions under which they can give rise to regions characterized by bare ground in agreement with relevant observational evidence. Towards that end, we return to our hexagonal symmetry function g0 of Equation (48) and consider it in more detail. In particular, focusing our attention upon a single hexagonal cell we note that each such individual g0 cell has an elevated central region with a maximum elevation of 3 at its centre which is bounded by a level curve of zero elevation appearing as the nearly circular loci found in Figures 10 and and11.11. The peripheral portion of each cell exterior to that central region is depressed, with the hexagonal cellular boundary of variable depth ranging from − 3/2 at its vertices to − 1 at the midpoint of its edges. Specifically, at the centre of these cells, the density function characteristic of III− has a value
while on its boundaries the density function characteristic of III+ ranges from
we see that there is about the same tendency for these density values of Equation (66) to be equal to, or less than, zero, and hence represent bare ground in agreement with the actual appearance of leopard or pearled bush, respectively. Finally, we shall investigate the predicted wavelength of these vegetation patterns. From Equations (14b) and (20), we can deduce that
which for the special case of Λ = 1 reduces to
these critical wavelengths being generated by the interplay between short-range facilitative and long-range competitive plant interactions. Specifically for α = α±0 of Equation (66) we can see from Equation (67c) that
and hence conclude that the vegetative distributions in leopard bush have a tendency to be more widely spaced than the bare patches which regularly punctuate the vegetation cover in pearled bush. Further the characteristic wavelength of tiger bush which corresponds to α = αc = 0.16 for this case satisfies
consistent with field observations.
We close by placing our theoretical predictions in the context of some recent vegetative pattern formation studies. Rietkerk et al.  developed an interaction-diffusion model system consisting of a set of three partial differential equations describing the spatio-temporal behaviour of plant density, soil water, and surface water, respectively. They performed numerical simulations on their model system with reflecting boundary conditions. Typical spatial patterns on flat ground generated in this manner were represented over a two-dimensional domain by Rietkerk et al.  for different amounts of rainfall R measured in mm/day. For R = 0.75 a hexagonal spotted pattern formed, changing into labyrinths for R = 1; while a hexagonal gap pattern was generated for R = 1.25. Von Hardenberg et al.  proposed a similar interaction–diffusion model system that only considered plant biomass and ground water but introduced more complicated plant growth and ground water flux terms. They studied this model by integrating their system numerically at different precipitation values ρ. For a plain (flat) landscape, von Hardenberg et al.  found that vegetation patterns of various forms evolved for p in a particular interval with hexagonal spots occurring at a relatively low p value, holes in a uniform coverage at a relatively high p value, and stripes (or labyrinths) at a p value intermediate to these. In both instances, the average plant biomass density increased with increasing R or p, respectively, although the scale of the patterns was in the range of metres or tens of metres for the Rietkerk et al.  model but only on the order of centimeters for the von Hardenberg et al.  model due to the fact that surface water diffusion which served as the driving mechanism for the former is much greater than ground water diffusion, the driving mechanism for the latter . Klausmeier  devised a related interaction–advection–diffusion model for vegetation and surface water that included a spatial gradient effect for its water rather than a diffusive one. Numerical integration of that system yielded no periodic vegetative pattern formation for flat ground and stripes when the downhill flow of rainwater was allowed for hillsides. Sherratt  presented a detailed theoretical analysis of pattern formation in the Klausmeier  model, derived formulae for the wavelength and migration speed of the predicted patterns, and systematically investigated how these depended on model parameters. Observe that the morphological sequence results described above that were produced by numerical simulation of the Rietkerk et al.  and von Hardenberg et al.  model systems are in qualitative agreement with our theoretical predictions for the propagator-inhibitor Lefever et al.  model equation once one realizes that the aridity parameter of the latter is inversely related to the rainfall and precipitation parameters of those systems. Also, we note that our systematic nonlinear stability analyses of vegetative pattern formation for this model equation is presented in the same spirit as was Sherratt's  detailed theoretical analysis of the Klausmeier  model system.
We conclude this discussion with some additional comments about our analyses in relation to those of Lejeune and his co-workers. As pointed out in the sketch of the development of the propagator-inhibitor model equation (3a) presented in Section 1, Lefever et al.  introduced scaling laws at variance to those catalogued in Equation (3b). Specifically, they took
while retaining our other scales of Equation (3b). Then, the terms
while all the other terms in Equation (3a) remain of O(ε3). Within the framework of our nonlinear stability analyses, let us examine the validity of the scaling laws of Equation (3b) where they differ from Equation (68a): namely
The first condition of Equation (69a) is consistent with all asymptotic scalings for the Stuart–Watson method of weakly nonlinear stability theory  while, using the marginal stability curve of Equation (15b) since our analysis is expected to be asymptotically valid as β → βc(α), we find that
which is consistent with the second condition of Equation (69a). Observe as noted earlier that with the scaling laws inherent to Equation (69) all the terms included in Equation (3a) are now of O(ε3). Indeed, in their seminal work Lefever and Lejeune  actually retained some terms of higher order appearing in the intermediate steps of Equation (2) which is the reason we did not attempt to compare our theoretical predictions with the simulation results contained therein. Also, Lejeune and Tlidi  claimed that they never observed rhombic patterns in their numerical simulations because such patterns are unstable with respect to stripes and hexagons. Our rhombic planform results of Section 3 show this assertion to be in error with respect to stripes except when α = αc. Actually, their simulation results were presented for both small- and large-scale windows. For the extended domains the striped patterns in the small windows were converted to zig-zag patterns while the hexagonal patterns become much less regular. All of these developments are consistent with our earlier discussion concerning the potential bistability of various types of stripes for α = αc and between rhombic and hexagonal patterns for α ≠ αc. In fact, Golovin et al.  found bistability between squares and hexagons in certain regions of their gravity number-capillary number parameter space for a Bénard–Marangoni surface tension-driven convection problem with poorly conducting boundaries.
In summary, we have performed systematic rhombic- and hexagonal-planform weakly nonlinear stability analyses on a truncated fourth-order partial differential time-evolution logistic equation describing the total plant biomass per unit area scaled with its carrying capacity over an unbounded flat spatial domain where vegetative pattern formation is induced by the intrinsic aridity of the region. Our main result was the production of plots of closed-form conditions in the relevant parameter space under which different vegetative patterns could be found. These patterns consisted of rhombic arrays of rectangular patches, parallel stripes or labyrinth-like mazes, and hexagonal distributions of spots or gaps. For a fixed value of the facilitation-to-competition interaction ratio, Λ, the stability behaviour of those vegetative patterns depended on two positive parameters α and β, which represented the mean plant biomass density and the square of the facilitation-to-competition range ratio, respectively. We showed that if β was greater than a certain threshold value only a uniform homogeneous distribution was possible, while should that not be the case, then the vegetative pattern selected depended on the size of α as well. In particular, one-dimensional parallel striped, labyrinth, or bicontinuous patterns were only possible if α = αc and β < βc; spots if α < αc and β2 < β < β−1; and gaps if α > αc and β2 < β < β−1; while if α−2 < α < αc or αc < α < α+2 and β < βc, rhombic patterns were stable for characteristic angles ϕ in two narrow bands flanking π/3. In addition, the same value of α, namely αc, playing a central role with respect to stripe formation in both rhombic- and hexagonal-planform anlayses served as a partial but independent check on these analyses. Further, for Λ held constant, α was determined by μ, the measure of aridity, with α a decreasing function of μ within the patterned range. Such vegetative patterning marked the transition between homogeneous Savannas or deserts when βcrit < β < βc and between isolated patches of bare ground or vegetation when β < βcrit. Hence this propagation–inhibition model was unified in the sense that it could be used to generate a wide variety of vegetative patterns which are observable in arid or semi-arid environments. We then compared these theoretical predictions with relevant numerical simulations and that observational evidence, finding consistency in the former case and very good agreement when parameter values were chosen appropriately in the latter case. Finally, we placed those results in the context of some recent nonlinear vegetation pattern formation studies. In conclusion, this aridity-driven vegetative distributions generation problem, involving a single evolution propagator–inhibitor logistic equation defined over a flat environment, is compatible with our long-range aim of employing the simplest reasonable natural science models that preserve the essential features of pattern formation and are still consistent with observation.
We would like to acknowledge Oliver Lejeune (as well as an anonymous reviewer) for his constructive comments that greatly contributed to this revised version of our manuscript.