Formation of square lattices in coupled pattern-forming systems

A wide variety of natural and laboratory systems can produce patterns of ripples, hexagons, or squares. The formation of stable square patterns from partial differential equation models requires specific cubic nonlinearities involving higher-order derivatives. Motivated by plant phyllotaxis, we demonstrate that the coupling of more than one pattern-forming system can produce square patterns without these special nonlinearities. Keywords-pattern formation; phyllotaxis; nanoscale structures

The steady-state solutions are approximately given by the sum of only a few Fourier modes.That is, for the case of ripples, u(x, y, •) Ae i k•(x,y) + c.c., (2) where A is a complex amplitude, and c.c. denotes the complex conjugate.For the case of hexagons, there is a triplet k 1 , k 2 , k 3 of wavevectors that are of the same modulus k c = | k i |, j = 1, 2, 3 and satisfy the condition k 1 + k 2 = k 3 (see Fig. 3 (a)) such that u(x, y, •) Σ 3 j=1 A j e i kj•(x,y) + c.c., (3) and for the case of squares, there is a pair of orthogonal wavevectors k 1 , k 2 of the same modulus A j e i kj•(x,y) + c.c.
Mathematical analysis of Equation (1) proceeds by first performing a linear stability analysis of the homogeneous steady-state solution u = 0.This determines the modulus k c of wavevectors that will be present in the pattern.All Fourier modes A(t)e i k•(x,y) with wavevector k on a circle of radius k c are linearly unstable so that their amplitudes A(t) grow in time for t 0. Once the amplitudes become large enough, nonlinear functions of these amplitudes become large enough to dampen the growth and to allow for interactions between the modes that determine the resulting steady-state pattern [5], [14].
While there is an abundance of examples of ripple and hexagon patterns in laboratory and natural systems, square patterns are relatively rare.As one example, in [13] it is shown that a cubic nonlinearity of the form δ∇ • (|∇u| 2 ∇u) in the equations of motion describing the surface evolution of a crystalline material being bombarded by a broad ion beam results in a pattern of squares.We propose in this paper an alternative to cubic terms involving higher-order derivatives, namely the coupling of two pattern-forming systems, that can also result in a pattern of squares.
The idea that coupled mechanisms may result in square patterns comes to us from observations of patterns on plants.Section 2 describes these plant patterns and a model for plant pattern formation proposed in [23] that couples biochemical and biomechanical mechanisms.Although we suggest in [23] that this model can produce square patterns, this paper gives the first numerical and analyti-cal evidence that this is the case by numerical simulations and analysis of a simplified system given in Section 3 that captures the key features of the original model.Numerical simulations of the simplified system produce square patterns for certain choices of parameters.The results agree with linear stability analysis and weakly nonlinear analysis of the simplified system, as given in the Appendix.

II. PHYLLOTAXIS AND A SYSTEM OF COUPLED PATTERN-FORMING EQUATIONS
Phyllotaxis refers to the arrangement (taxis) on plants of leaves (phylla) or their analogs such as bracts on a pine cone, seeds in a seedhead, or spines on a cactus.We describe in Section II.A how a pattern of squares underlies many of these phyllotactic patterns.In Section II.B, we review a systems of coupled PDEs, proposed in [23] as a model for the formation of phyllotactic patterns on plants.The square bracts on the pinecone shown in Fig. 4 (a) have been numbered in sequence of their distance from the center of the cone.Connecting bracts that have adjacent sides results in the eight (yellow) counterclockwise spirals and thirteen (red) clockwise spirals.The numbers of spirals in the clockwise and counterclockwise spiral families are called the parastichy numbers for the pattern.The pair of parastichy numbers 8 and 13 is not unusual.Indeed, the spiral numbers observed on plants are typically consecutive members of the Fibonacci sequence {F j } = 1, 1, 2, 3, 5, 8, 13, 21, . . . .The parastichy numbers may vary even within one plant's pattern.This is illustrated by the sunflower seed heads of Fig. 4 (b) and Fig. 5.The transition between spiral numbers in the Fibonacci sequence can, in fact, be continuous.To understand this, consider the function
The key observation is that the square patterns that are evident at radii r = r 1 and r = r 5 are formed not by only two Fourier modes with wavevectors as in Fig. 3 (b), but by overlapping triads of modes that satisfy summation relations similar to the wavevectors of Fig. 3 (a) that produce hexagons.As depicted in Fig. 3 (c), in order for this to occur, two of the vectors are larger in modulus than the other two.These are the wavevectors corresponding to the smaller amplitudes.For example, at r = r 1 , the wavevectors k 6 , k 9 corresponding to the smaller amplitudes A 8 = A 34 are longer than the wavevectors k 7 , k 8 corresponding to the larger amplitudes A 13 = A 21 .

B. A mechanistic model for phyllotactic patterns
In this section, we review a system of partial differential equation that has been proposed as a model for the formation of phyllotactic patterns and that incorporates biochemical and biophysical mechanisms, each of which could produce a pattern on its own.Growth of a shoot tip or formation of flowers occurs in regions of active cell growth and division at apical meristems.A schematic diagram of a shoot apical meristem (SAM) is shown in Fig. 7. Small bumps called primordia on the plant surface, which will become leaves, form not at the very center of the SAM, but in an annular region which we call the generative region and which is marked as Region 2 in Fig. 7.
What mechanisms lead to the formation of primordia in the generative region?There is evidence for both biomechanical and biochemical mechanisms which may interact with each other.
The idea of the biomechanical mechanism is as follows: If the outer skin (the tunica) of the plant is growing more quickly than the inner tissue, then a compressive stress will build up in the tunica.If that compressive stress increases above a large enough threshold, then the tunica will become Biomath 5 (2016), 1612181, http://dx.doi.org/10.11145/j.biomath.2016.12.181 unstable and buckle under the stress.Biologist Paul Green proposed in the 2000's that primordia are the result of this buckling [9]- [11].
The plant hormone auxin influences cell growth, and other experiments suggest that auxin itself may be spatially patterned in the generative region, with primordia forming where there is a higher auxin concentration [15], [16].Auxin is produced uniformly throughout the generative region, but the key idea of the groups of Kuhlemeier and Meyerowitz [16], [30] is that PIN1 proteins in cell walls transport auxin from cells with lower concentrations of auxin to cells with higher auxin concentrations.This produces an instability that allows for pattern formation.
The biomechanical and biochemical mechanisms may interact in that stress states can impact the action of PIN1 proteins.In [23], we incorporate both mechanisms into a mathematical model of three partial differential equations for the tunica surface deformation w(x, y, t), a potential F (x, y, t) for the stresses in the tunica, and the auxin concentration difference g(x, y, t) from a mean auxin concentration.
We refer the reader to [23] for a complete description of this model, which reads, in nondimensionalized parameters, where [f, g] = f xx g yy + f yy g xx − 2f xy g xy .Equations (6a) and (6b) are the Föppl-von Kárman equations for a thin elastic to which a compressive stress has been applied.Equation (6c) is a continuum approximation of a spatially discrete model proposed in [16] for auxin transport.Both the mechanical system (6a, 6b) and the auxin system (6c) may produce an instability of the homogeneous state w = F = g = 0 to a pattern.For the mechanical system, this instability Fig. 7: A schematic of the plant shoot apical meristem (SAM).Cells form but primorida do not form at the center of the SAM (Region 1).Region 2 is the annular generative region where primoria form.In Region 3, no new primordia form, but there is active cell generation and differentiation.occurs if the compressive stress, expressed in the nondimensional parameter P , exceeds a critical value.For the auxin system, this instability occurs if the relative magnitude of auxin transport compared to auxin diffusion, expressed in the nondimensional parameter H, exceeds a critical value.If both the mechanical instability and the auxin instability are active, the possibly different natural wavelengths of the patterns that would result from either instability alone allow for differences in phyllotactic configurations (the underlying lattice) and the surface deformation.In [23] we analyze a variety of scenarios in which the elastic and auxin instabilities may cooperate or compete.

III. SQUARE PATTERN FORMATION IN COUPLED
EQUATIONS OF SWIFT-HOHENBERG TYPE Equations (6a) and (6c) are both of Swift-Hohenberg type (1) with the addition of nonlinear terms that couple the equations.We expect therefore for the following simpler and more analyti-Biomath 5 (2016), 1612181, http://dx.doi.org/10.11145/j.biomath.2016.12.181 cally tractable system of coupled Swift-Hohenberg equations for two fields u(x, y, t) and v(x, y, t) to produce similar steady-state patterns: This system has the uniform steady-state solution u(x, y, t) = v(x, y, t) = 0. Linear stability analysis, given in the Appendix, reveals that in the absence of linear coupling (α 1 = α 2 = 0), u(x, y, t) = 0 is stable for P ≤ P c = 1, but unstable to Fourier modes with wavevectors of modulus close to |k| = 1 for The modified conditions for instabilities in the presence of linear coupling are given the the Appendix, but the key point is that the parameter L determines how the the wavelength of the pattern favored by the equation for v compares to that for u.If L = 1, then both equations would yield patterns of the same natural wavelength, but if L is larger (smaller) than 1, then the wavelength of the pattern favored by the equation for v will be smaller (larger) than that favored by the equation for u.
If P = P c + χ and H = H c + χ, where χ ∼ O(1), are slightly above their respective bifurcation thresholds (as measured by the small parameter ), Fourier modes with moduli close to 1 and 4 √ L grow in amplitude with time and interact via the nonlinear terms in the equations.The wavevectors for these modes are depicted as the two circles in Fig. 3 (c).In the Appendix, we demonstrate a weakly nonlinear asymptotic analysis that allows us to derive a system of ordinary differential equations for the amplitudes of these excited Fourier modes.This analysis begins with an Ansatz for the form of the solution, namely Here, N is the number of interacting modes in the Fourier expansion of the order-term, A j = A j (T = εt), B j = B j (T = εt), and ε 2 u 1 + ε 3 u 2 + ..., and ε 2 v 1 + ε 3 v 2 + ... are the correction terms.A condition for solvability of the correction terms leads to a set of ordinary differential equations for the time-evolution of the complex amplitudes A j (t) and B j (t).
The result of a numerical simulation of Equation (7) for a choice of parameter values that includes nonlinear coupling terms but not linear coupling (α 1 = α 2 = 0) is shown in Fig. 8.The initial conditions are low-amplitude white noise.We employ a Fourier spectral method with periodic boundary conditions and a fourth-order exponential time differencing Runge-Kutta method for the time stepping as the numerical technique, and the spatial grid is 256 × 256.
The Fourier transform of the pattern shown in Fig. 8 (a) is shown in Fig. 8 (b).There are two circles of excited wavevectors, as marked in Fig. 8 (c).The reason for this is the parameter choice L = 4.7 in (7), which allows for the wavevectors of modulus 4 √ 4.7 1.47 to be excited by the field v, while wavevectors of length 1 are excited by the field u (see the linear stability analysis in the Appendix).These wavevectors interact via nonlinear terms in Equations ( 7), and a discrete set of Fourier modes with wavevectors k 1 , . . ., k 4 such that k 1 + k 2 = k 3 and k 2 + k 3 = k 4 , as marked in Fig. 8 (c) (which may be compared to Fig. 3 (c)), dominates the pattern.
Motivated by the results of the numerical simulation, we carry out in the Appendix the weakly nonlinear asymptotic analysis, choosing in Equation (8) N = 4 and the overlapping triad conditions k 1 + k 2 = k 3 and k 2 + k 3 = k 4 .This results in a system of eight ordinary differential equations for the amplitudes A j , B j , j = 1, . . ., 4 of the Fourier modes in (8).This is the system (15).A numerical simulation of (15) for the parameter values of the simulation in Fig. 8 is shown in Fig. 9.The amplitudes reach a steady state in which This is consistent with the Fourier spectrum shown in Fig. 8 (b,c), in which the modes on the circle of larger radius have smaller amplitude.It is also consistent with the motivation given in Section II.A of patterns observed in plant phyllotaxis.
Other possible patterns and their stability would be found by a bifurcation analysis of the amplitude equations (15).The relevance of this analysis to plant phyllotaxis would require that the length scale of the patterns produced by the system (6) are in accord with those observed on plants.In [23], we provide suggestions on experiments to  determine some of the parameter values relevant to plants.Although motivated by phyllotaxis in this paper, the framework of coupled patternforming systems is relevant to other phenomena.This includes nanoscale pattern formation induced by bombarding a binary alloy by a broad ion beam [1], [2], [7], [29].In these experiments, collision cascades induced by the ions hitting the surface of the alloy result in the sputtering of material from the surface.Curvature-dependent sputter yield and phase separation are two independent pattern-forming mechanisms that may both contribute to the observed pattern.Squares can also appear in Turing patterns: Recently, Li and colleauges form patterns of squares in simulations of coupled reaction-diffusion equations [20].
(the derivation of the nonlinear amplitude equations) for the system (7).

Linear stability analysis
The system (7) has the uniform steady-state solution u = u s = 0, v = v s = 0. We will examine the stability of this solution by introducing a small perturbation of this solution and determining if the perturbation has linear growth or decay.Write the perturbation as where k = (k x , k y ), and û and v are constants.Inserting (9) into the linearization of Equation ( 7), we obtain where k 2 = k 2 x + k 2 y .This vector equation implies that σ is an eigenvalue of the matrix in Equation (10).The two solutions for σ are where σ u (k) .= −k 4 + 2P k 2 − 1 and σ v (k) .= −k 4 + 2Hk 2 − L. We will use σ + ( k) to denote the larger of the two eigenvalues.
Note that if v = 0, we have only the equation The coupling terms alter these values, and for positive Since σ is the growth rate for the system and instabilities occur when σ > 0, the coupling terms actually help create an instability.The choice of σ + is non-trivial for L = 1 even in the case where there is no linear coupling (α 1 = α 2 = 0); σ + = σ u ( k) for some values of k and σ + = σ v ( k) for others.By definition, we always choose the larger of the two values.
The parameters P and H thus serve as bifurcation parameters with respective bifurcation values In the absence of linear coupling, the uniform steadystate solution is linearly stable for P < P cr and H < H cr .Perturbations with certain wavevectors are unstable for P > P cr and H > H cr .Recall that for general values of α 1 and α 2 (with α 1 α 2 > 0), the modified value of σ + ( k) will be greater than either σ u or σ v (say σ * for the general case).
We identify the set of active wavevectors, k, to be those for which σ + (k) > −σ * for some small −σ * < 0. The goal is to analyze the cases where σ + is just above zero, when the active modes begin to interact through the nonlinear terms in the equations.

Nonlinear Analysis
For σ + above the critical value of σ + = 0, there are active modes, and these will interact through the nonlinear terms in (7).We derive evolution equations for the amplitudes of the active modes by asymptotically expanding the solutions u and v with respect to a small parameter ε measuring how far σ + is above 0.That is, we assume that the bifurcation parameters are close to their respective critical values: for χ ∼ O(1), We also assume that the coefficients of the cubic terms are of order −1 : As an Ansatz for the form of the solution to (7), we assume that u and v are an order-ε perturbation from the steady state solution plus correction terms that are of higher order in ε: Here, N is the number of interacting modes in the Fourier expansion of the order-term, A j = A j (T = εt), B j = B j (T = εt), and ε 2 u 1 +ε 3 u 2 + ..., and ε 2 v 1 + ε 3 v 2 + ... are the correction terms.Inserting (13) into the system (7) and collecting coefficients of powers of ε yields, at order ε simply the expression 0 = 0, but at order ε 2 , we obtain the relations The first of these equations has the form (∆ 2 + 2P cr ∆ + 1)u 1 = Ce i k x where P cr = 1 and | k| = 1 (recall that this value was found to be critical for instabilities in the linear case with no coupling when P ≥ 0).This equation has the form u 1 = De i k x .Thus, if the coefficient of e i kj x on the right-hand side of ( 14) is nonzero, resonance allows the solution to grow without bound, and our asymptotic expansion is invalid since the correction terms are no longer small.Our solvability condition is therefore that the coefficient of e i kj x on the right-hand side of (14) be zero.We now examine the case of ( 14) for N = 4, where k 1 + k 2 = k 3 , and k 2 + k 3 = k 4 .Requiring that the coefficients of the e i kj x terms, j = 1 . . . 4 sum to zero results in a set of eight differential equations for the time evolution of the amplitudes A j and B j , j = 1 . . . 4 which read Now observe that for k 2 = 1 We obtain a similar relation for σ v ( k).This allows us to rewrite the linear terms in the amplitude equations as 2χA i = ε −1 σ u ( k)A i and Recalling that T = εt, we rescale A i and B i as Using these and the earlier rescaling γ1,2 = εγ 1,2 , we obtain, after multiplying through by ε 2 , the following set of differential equations (the hat notation has been suppressed on all A i and B i ): state that is a pattern of ripples (Fig 1 (a)), hexagons (Fig 1 (b)), or squares (Fig 1 (c)).

Fig. 4 :
Fig. 4: (a) A pinecone with the bracts numbered in order of their distance to the center.Also marked are eight counterclockwise spirals (in yellow) and thirteen clockwise spirals (in red) formed by connecting adjacent bracts.(b) From [3], a sunflower seed head with clockwise and counterclockwise spirals marked.The spiral families are different near the outer boundary of the seed head compared to the center.