{'fraction_non_alphanumeric': 0.04540492774716706, 'fraction_numerical': 0.018557022559517623, 'mean_word_length': 4.243526846552194, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simple and Superlattice Turing Patterns in Reaction-Diffusion Systems: Bifurcation, Bistability, and Parameter Collapse

This paper investigates the competition between both simple (e.g. stripes, hexagons) and ``superlattice'' (super squares, super hexagons) Turing patterns in two-component reaction-diffusion systems. ``Superlattice'' patterns are formed from eight or twelve Fourier modes, and feature structure at two different length scales. Using perturbation theory, we derive simple analytical expressions for the bifurcation equation coefficients on both rhombic and hexagonal lattices. These expressions show that, no matter how complicated the reaction kinectics, the nonlinear reaction terms reduce to just four effective terms within the bifurcation equation coefficients. Moreover, at the hexagonal degeneracy -- when the quadratic term in the hexagonal bifurcation equation disappears -- the number of effective system parameters drops to two, allowing a complete characterization of the possible bifurcation results at this degeneracy. The general results are then applied to specific model equations, to investigate the stability of different patterns within those models.


Introduction
Regular patterns arise in a wide variety of physical, chemical and biological systems that are driven from equilibrium. The origin of these spatial patterns can often be traced to a symmetry-breaking instability of a spatially-homogeneous state. Such an instability mechanism was proposed by Alan Turing in 1952 [1] for pattern formation in chemical systems featuring both reaction and diffusion. The key insight of Turing was that diffusion, which is usually thought of as a stabilizing mechanism, can actually act as a destabilizing mechanism of the spatially-uniform state. This is the case, for example, in a class of two-component activator-inhibitor systems, provided the inhibitor diffuses more rapidly than the activator. This competition between long-range inhibition and short-range activation can lead to spatially inhomogeneous steady states called Turing patterns [2]. Experimental evidence of Turing patterns was discovered only recently [3], and a number of experiments on the Turing instability have now been carried out on the chlorite-iodide-malonic acid (CIMA) reaction in a gel [3,4,5]. These experiments have demonstrated the existence and stability of a variety of regular patterns, including stripes, hexagons, and rhombs. At the same time, a two variable model of the CIMA reaction has been proposed by Lengyel and Epstein [6].
Equivariant bifurcation theory [7] is a powerful tool for analyzing the nonlinear evolution of symmetry-breaking instabilities in pattern-forming systems. This approach determines the generic behaviors associated with bifurcation problems that are equivariant with respect to a given group of linear transformations. In particular, the form of the bifurcation problem is restricted by symmetry considerations, limiting the behavior to some finite set of possibilities. The problem of determining which type of behavior occurs in a given physical system then reduces to one of computing the coefficients of specific nonlinear terms in the bifurcation problem. In this paper we perform this computation for a general two-component reaction-diffusion system in a neighborhood of a Turing bifurcation. We consider two different symmetry groups for the problem, one associated with Turing patterns that tile the plane in a square lattice, and the other associated with patterns that tile the plane in a hexagonal lattice. The symmetry groups are D 4+ T 2 and D 6+ T 2 , respectively, where D n (n = 4, 6) is the dihedral group of symmetries of the fundamental tile, and T 2 is the two-torus of translation symmetries associated with the doubly-periodic solutions. The fundamental irreducible representations of these groups, which are four-and six-dimensional, lead respectively to bifurcation problems that address competition between stripes and simple squares, and between stripes and simple hexagons. The higher dimensional irreducible representations, which are eight-and twelve-dimensional, lead to bifurcation problems that address competition between these simple Turing patterns and a class of patterns that are spatially-periodic on a larger scale than the simple patterns, but still arise in a primary bifurcation from the homogeneous state at the Turing bifurcation point [8,9,10,11]. An example of such a "superlattice" pattern was recently observed in an experiment on parametrically excited surface waves [12]. Finally, the higherdimensional group representations also allow us to investigate the relative stability properties of rhombs and squares or rhombs and hexagons.
A number of earlier studies have used bifurcation theory to investigate Turing pattern formation for reaction-diffusion systems. For instance, Othmer [13] described a group theoretic approach for analyzing the Turing bifurcation nearly twenty years ago, in 1979. Recently Callahan and Knobloch [14] derived and analyzed the general equivariant bifurcation problems for three-dimensional Turing patterns having the spatial periodicity of the simple cubic, face-centered cubic and body-centered cubic lattices. They applied their general results to the Brusselator and the Lengyel-Epstein CIMA reaction-diffusion models [15]. Gunaratne, et al. [16] developed an approach to pattern formation problems based on a model equation akin to Newell-Whitehead-Segel, but with full rotational symmetry. They used this model to investigate bistability of hexagons and rhombs, and tested their predictions against laboratory experiments.
In addition to these general bifurcation studies of Turing pattern formation, a number of specific reaction-diffusion models have been investigated analytically and numerically. Dufiet and Boissonade [17] investigated numerically the formation of two-dimensional Turing patterns for the Schnackenberg reaction-diffusion system. The formation of three-dimensional Turing structures has been investigated for the Brusselator by DeWit, et al. [18]. Turing patterns for the two-variable reaction-kinetics based model of the CIMA reaction, proposed by Lengyel and Epstein [6], have been investigated using bifurcation theory [19], as well as numerically [20]. This paper extends existing bifurcation studies of Turing patterns in reactiondiffusion systems in two primary ways. The first, as discussed above, is that we consider the higher dimensional (8 and 12, respectively) irreducible group representations associated with the square and hexagonal lattice tilings. This corresponds to choosing a fundamental tile for the periodic patterns large enough to admit many copies of the simple Turing patterns -stripes, hexagons, squares, and rhombs -as well as new solutions in the form of "superlattice" patterns. In this way we extend the standard stripes vs. hexagons stability analysis, while also determining conditions under which other more exotic Turing patterns might arise.
The second way in which our analysis differs from earlier studies is that it does not restrict attention to a specific model (e.g. Schnackenberg, Brusselator, Lengyel-Epstein CIMA model, etc.). The coefficients in the bifurcation equations are derived for a general two-component reaction diffusion system. This allows us to make a number of general statements about Turing patterns in these systems. For example, we show that the bifurcation results depend on four effective parameters that may be simply computed from the nonlinear terms in the reaction-kinetics. Moreover, in the case of the standard degenerate hexagonal bifurcation problem, we show that the details of the reaction-kinetics enter the analysis through just two parameters. Even more surprising is the result that at this degenerate point all angle-dependence, as associated with rhombic patterns, drops out of the corresponding bifurcation problems; in this case we show that all rhombic patterns are unstable to rolls, even when hexagons are the preferred planform. This is in contrast to the general stability analysis of Gunaratne et al. [16], which suggest that rhombs with angle close to π/3 should be stable when hexagons are stable. Thus our analysis indicates that two-component reaction-diffusion models may have a more rigid mathematical structure than anticipated.
Our paper is organized as follows: Section 2 reviews the necessary conditions for a Turing instability in a two-component reaction-diffusion system. It also provides the necessary background on the mathematical framework for the bifurcation analysis. Section 3 calculates the bifurcation equation coefficients using perturbation methods, with the final results summarized in an Appendix. Section 4 describes the concept of parameter collapse, where multiple system parameters collapse into a small number of effective parameters. Special attention is given to the degenerate bifurcation problem on the hexagonal lattice, where we show that the details of the reaction-kinetics enter the analysis in this case through only two effective parameters. The resulting predictions are then then tested numerically; specifically, we present an example of tri-stability between stripes, simple hexagons, and super hexagons in a neighborhood of the degenerate bifurcation problem. Section 5 applies the general results to a model of the CIMA reaction proposed by Lengyel and Epstein [6] and a model demonstrating super squares. Finally, we suggest some directions for future work and summarize the key points of the paper in the conclusions section.

Linearized Problem
We begin by summarizing the conditions for a Turing instability to occur as some parameter is varied in the reaction-diffusion system We assume that the species diffuse at different rates and that v measures the concentration of the more rapidly diffusing species. Thus the ratio of diffusion coefficients is K > 1. The Turing instability first occurs as a symmetry-breaking steady state bifurcation of a spatially-uniform steady state of (1); associated with the instability is a critical wave number q c = 0, which is determined below. Without loss of generality, we take the spatially-uniform state to be u = v = 0. We expand f (u, v) and g(u, v) about this state as f (u, v) = au + bv + . . . g(u, v) = cu + dv + . . . yielding the linearized reaction-diffusion system The stability of the spatially uniform solution is determined by substituting (2) to obtain the eigenvalue problem The eigenvalues σ 1 (q), σ 2 (q) depend on the parameters a, b, c, d, K, as well as the wavenumber of the perturbation q. A Turing bifurcation occurs for parameter values where there is a zero eigenvalue for some q = q c = 0 with all other modes being damped, i.e., σ 1 (q c ) < σ 2 (q c ) = 0, Re σ j (q) < 0, j = 1, 2, for all q = q c .
These inequalities imply that at the bifurcation point. The critical wave number q c is given by Since Det(L(q)) is quadratic in q 2 and concave up, the second of conditions (4) is assured if q c is a minimum of Det(L(q)): ∂Det(L(q)) ∂q Thus the parameters at the Turing bifurcation (a = a and must also satisfy the inequalities to ensure that q 2 c > 0 and that perturbations at q = 0 decay. From the initial assumption that K > 1, it follows that a c > 0 and d c < 0, which for b c < 0, c c > 0 is consistent with the interpretation of u as the concentration of the activator and v as the concentration of the inhibitor in (1).

Spatially-Periodic Solutions of the Nonlinear Problem
For the two-dimensional spatially-unbounded problem there are infinitely many neutrally stable Fourier modes at the bifurcation point: all Fourier modes associated with wave vectors on the critical circle |q| = q c are neutrally stable. In the remainder of the paper we focus on spatially doubly-periodic solutions of (1): where r ≡ (x, y) and d 1 , d 2 are linearly independent vectors in R 2 . We further restrict to the cases where |d 1 | = |d 2 |. We then write each field u, v in terms of its Fourier series, e.g.
where k 1 and k 2 are unit vectors that satisfy for an appropriate choice of q. Since the Fourier transform of spatially-periodic solutions is discrete, only a finite number of Fourier modes can have wave vectors q that lie on the critical circle |q| = q c at the bifurcation point; all other modes will be linearly damped. Thus, in this setting, finite-dimensional bifurcation theory may be applied in a neighborhood of the Turing bifurcation point. For Turing patterns that tile the plane in a square lattice the unit vectors k 1 and k 2 are separated by an angle of π/2, whereas in the hexagonal case they are separated by 2π/3; all angles that are not a multiple of π/3 or π/2 correspond to rhombic lattices. The parameter q in (10), which is determined by the imposed period, also dictates the the spacing of the lattice points in the Fourier transform, and may in turn be used to determine which lattice points (if any) lie on the critical circle. Of particular interest here is the ratio of this scale to the natural scale of the problem q c . For instance, if q/q c = 1 then on the square lattice there are four lattice points which lie on the critical circle: (j, k) = (±1, 0), (0, ±1) in (10), corresponding to Fourier modes e iqcx , e −iqcx , e iqcy , e −iqcy .
In this paper we consider a discrete set of values of q, parameterized by an integer pair (m, n), which admit wave vectors in (10) with length q c , i.e., for which |q(m k 1 + n k 2 )| = q c . Specifically, on the square lattice the values of q are related to (m, n) by and on the hexagonal lattice by  In a neighborhood of the bifurcation point, , we apply our analysis separately to the cases where there are eight critical Fourier modes on the square lattice and twelve on the hexagonal lattice (see Figure 1). Specifically, on the square lattice, we investigate primary steady solution branches of the form where z 1 , . . . , z 4 are the complex-valued amplitudes of the critical modes, · · · represents the linearly-damped Fourier modes that are slaved to the critical modes, and K 1 = (m, n), K 2 = (−n, m), K 3 = (n, m), and K 4 = (−m, n). Note that, while K 1 · K 2 = K 3 ·K 4 = 0, the angle between K 1 and K 3 depends on the integer pair (m, n). Similarly, on the hexagonal lattice, where . K 2 and K 3 are obtained by rotating K 1 by ±2π/3, and K 5 and K 6 are obtained by rotating K 4 by ±2π/3 (see Figure 1). Note that the angle between K 1 and K 4 again depends on (m, n).
Dionne and Golubitsky [8] used the equivariant branching lemma [7] to prove that for the square lattice cases (13), the following types of steady solutions arise in pitchfork bifurcations from the spatially-uniform state: where in each case x is a real amplitude. (Examples of super and anti-squares are given in Section 5.2, Figure 9). While the striped and simple square patterns are essentially the same for each m and n, the solutions in the form of rhombs, super squares and anti-squares differ for each co-prime m and n since the angle between K 1 and K 3 (and K 1 and K 4 ) depends on m and n. The general form of the amplitude equations describing the branching and relative stability of the square lattice patterns is derived in [9]. To cubic order the bifurcation equations take the simple forṁ with similar expressions forż 2 ,ż 3 ,ż 4 and their conjugates. Here µ is the bifurcation parameter and a 1 , . . . , a 4 are real. Linearizing (15) about a particular solution generates a set of eigenvalues which determines the stability of that solution with respect to perturbations on the lattice. The signs of these eigenvalues are given in Table 1 in terms of a 1 , . . . , a 4 ; a positive eigenvalue indicates instability. Note that higher order terms are necessary to determine the relative stability of the super square and antisquare patterns [9]. In the next section we compute the coefficients a 1 , . . . , a 4 needed to evaluate the signs of the eigenvalues in terms of the angle between K 1 and K 3 and in terms of parameters in the general reaction diffusion system (1). Dionne and Golubitsky [8] also considered the hexagonal lattice problem and used group theoretic methods to determine that there are primary branches of the form • Stripes: z ≡ (z 1 , z 2 , z 3 , z 4 , z 5 , z 6 ) = (x, 0, 0, 0, 0, 0) • Simple hexagons: z = (x, x, x, 0, 0, 0) • Rhombs: z = (x, 0, 0, x, 0, 0), (x, 0, 0, 0, x, 0), and (x, 0, 0, 0, 0, x) Recently, Silber and Proctor [11] showed that there is an additional primary solution branch that has triangular symmetry: where, near the bifurcation point, the argument of the complex amplitude z is ∼ π 3 , 2π 3 . The stripes and simple hexagons are the same for each integer pair (m, n), while the other solutions change with the values of m and n. The general form of the bifurcation equations describing the evolution of the critical modes z 1 , . . . , z 6 is derived in [9]. The cubic truncation iṡ Anti-Squares same as super squares, except with similar expressions forż 2 , . . . ,ż 6 . While the stripes and rhombs solutions arise from the uniform state in a pitchfork bifurcation, all other solutions bifurcate transcritically. The relative stability of the primary patterns are summarized from [9] in Table 2. Note that unless the quadratic coefficient γ is zero, all branches bifurcate unstably. Thus we will consider a degenerate bifurcation problem for which γ ≈ 0. Also note that, as in the case of super squares and anti-squares, the relative stability of the super hexagons and super triangles is determined by terms higher than cubic in (16) [11].

Reduction to Bifurcation Equations
Although symmetry determines the form of the bifurcation equations, the particular system determines a n (and b n ) in equations (15) (and (16)). Each coefficient a n is associated with a particular amplitude z n in (15), and each of these amplitudes is in turn associated with a Fourier wave vector K n . In this section we use perturbation theory to derive the coefficients a n , γ, and b n from the full system (1).
The coefficients a 1 , . . . , a 4 and b 1 (= a 1 ), b 4 , b 5 , b 6 may all be calculated by considering a single bifurcation problem involving just two critical Fourier modes; the remaining coefficients γ and b 2 are calculated by considering a bifurcation to hexagons. To understand the need for only two modes in the computation of a 1 , . . . , a 4 , b 4 , b 5 , b 6 consider, for example, the square lattice, where the K 2 , K 3 , and K 4 vectors are all related to the K 1 vector by an angle θ. Let θ j denote the angle between K 1 and K j ; then θ 2 = π/2 and θ 4 = π/2 + θ 3 where θ 3 is related to the integer pair (m, n) by (see Figure 1). The coefficients a j are a function of the angle θ j : Similarly, on the hexagonal lattice, where θ 5 = θ 4 + 2π/3 and θ 6 = θ 4 + 4π/3, with θ 4 determined by the integers (m, n). From symmetry considerations we know that h(θ) = h(π − θ). A calculation with two critical modes will determine h(θ), and hence determine all of the angle-dependent coefficients in the bifurcation equations.
To compute the bifurcation equations (19) and (21) from equation (1) we first Taylor-expand f (u, v) and g(u, v) through cubic order in u and v: Here F 2 (u, v), G 2 (u, v), F 3 (u, v), and G 3 (u, v) denote the quadratic and cubic nonlinear terms: where the subscripts on f and g denote partial derivatives evaluated at u = v = 0. We use perturbation theory to compute the small amplitude, slow evolution on the center manifold. Let where w 1 takes the form of either (18) or (20) for the respective rhombic or hexagonal calculations, and W 2 and W 3 represent higher order terms to be determined. The vector (u 1 , v 1 ) T is a null-vector of L(q c ), at the bifurcation point; see equation (3). We further assume that the linear parameters a, b, c, d, K are within O(ǫ 2 ) of the critical bifurcation values, e.g. a = a c + ǫ 2 λ a , etc. where a c , b c , c c , d c , K c satisfy (8) and determine q c via substitution into (7). With these assumptions, equation (1) may be written as plus higher order terms, where
There are no terms on the right hand side of (24) with wave number q c , so the solvability condition is ∂ ∂t 1 w rh = 0.
Since there is no dependence on the t 1 -time scale it will be omitted from the remaining analysis. Solving equation (24) determines where u n and v n (n = 2, . . . , 5) are constant and satisfy where c 2 = 1 and c 3 = c 4 = c 5 = 2. Explicit expressions for u n and v n are given in the Appendix. The equation at O(ǫ 3 ) is where (U 2 , V 2 ) are given by (25). There are now terms on the right hand side of the equation with spatial dependence e ik1·r and e ik θ ·r ; the operator L 0 is not invertible for terms with this Fourier dependence. To ensure that these terms will lie in the range of L 0 we apply the Fredholm alternative theorem to obtain the solvability condition. In this computation we use the fact that the nullspace of the adjoint linear operator L † 0 is spanned by (ũ,ṽ)e ik1·r , (ũ,ṽ)e ik θ ·r , where (ũ,ṽ) = c, −a c + q 2 c .
The inner product of equation (27) with these left null-vectors yields the solvability condition, which takes the expected form of the bifurcation equations (19). After rescaling time to absorb the factor u 1ũ + v 1ṽ , which is always positive, we obtain the rhombic bifurcation equation coefficients and where The expressions for u n and v n , n = 2, . . . , 5, are given in the Appendix.

Hexagonal Lattice Computation
The hexagonal calculation proceeds in a fashion similar to the rhombic case, but requires some additional work due to the presence of resonant terms with wave number q c which are generated by the quadratic nonlinearities.
To simplify the presentation we let z 1 = z 2 = z 3 = z h in (20). Since the form of the amplitude equation is known from equation (16), and b 1 (= a 1 ) is known from the rhombic lattice calculation, the remaining coefficients γ and b 2 may be extracted from the quadratic and cubic coefficient expressions, respectively.
At O(ǫ 2 ) the equation to be solved is where w hex is given by (20). The term w 2 hex generates four terms with distinct wave numbers: From w 7 we see that there are now terms on the right hand side of (32) with wave number q c , for which L 0 is not invertible. Unlike in the rhombic case the t 1 time scale is necessary here, reflecting that the small-amplitude dynamics are dominated by the quadratic nonlinearities. As in the rhombic calculation the Fredholm alternative theorem is applied, giving the solvability condition The factor (u 1ũ +v 1ṽ ) is positive (from conditions (8) and (9)) and will later be absorbed by a rescaling of time. Note that if γ is O(ǫ) then the t 1 -time scale is unnecessary. This is the situation of interest here, since we know from Table 2 that when γ is O(1) there are no stable small-amplitude steady states. However, in the following computation we must retain both time scales to compute the general form of the cubic coefficient b 2 .
In our subsequent analysis we will focus on the degenerate problem γ = 0. Continuing as in the rhombic case, the solution of (32) is where u 2 , v 2 , u 5 , and v 5 are known from the rhombic lattice computation, and u 6 and v 6 are computed in a similar fashion. Equations (32) and (33) lead to the following solvable equation for (u 7 , v 7 ): We choose the particular solution of this singular problem to be Again the Fredholm alternative theorem must be applied, giving a solvability condition of the form By letting z(t) = ǫz h (t 1 , t 2 ), µ = ǫ 2 λ, and rescaling time to absorb the factor u 1ũ + v 1ṽ on the left hand side, we obtain the reconstituted bifurcation equatioṅ z = µz + γz * 2 + (b 1 + 2b 2 )|z| 2 z.
The coefficients µ, γ, b 1 , and b 2 are given in the Appendix.

Parameter Collapse
Naively one might expect that the cubic coefficients a n and b n in the bifurcation equations may be adjusted to take on any desired value, by varying the many nonlinear coefficients in the Taylor expansions of f (u, v) and g(u, v). However, the expressions for a n and b n , given in Appendix A, clearly show that all eight cubic coefficients contained within F 3 (u 1 , v 1 ) and G 3 (u 1 , v 1 ) collapse into the single expression β ≡ũF 3 (u 1 , v 1 ) +ṽG 3 (u 1 , v 1 ). Similarly, the six quadratic coefficients in F 2 (u, v) and G 2 (u, v) are contained in the coefficient expressions solely in the four terms F 2 (u 1 , v 1 ), G 2 (u 1 , v 1 ), η 1 and η 2 , given by (31). Moreover, these last four quantities satisfy the relation Therefore the six quadratic coefficients collapse into three independent parameters. In summary, the fourteen nonlinear coefficients in the Taylor expansions of f (u, v) and g(u, v) enter our analysis in just four independent combinations.

Hexagonal degeneracy
The parameter collapse is even more dramatic at the degeneracy γ = 0 in the hexagonal bifurcation problem (16). An additional striking phenomena occurs at this point, which results in all θ-dependence of the coefficients dropping out, thereby severely restricting the types of patterns which may bifurcate stably. Below we show that the bifurcation equations depend on just two quantities determined from the reaction kinetics of (1). Thus, we may simply and completely characterize the possible bifurcation scenarios for all reaction-diffusion systems of the form (1) restricted to the hexagonal bifurcation problem at the degeneracy. The degeneracy condition is which implies by substitution into (35). These two equations may be solved for F 2 and η 1 in terms of G 2 and η 2 , and therefore these four parameters collapse into the single parameter combination η 2 G 2 (u 1 , v 1 ) in the coefficient expressions for a n and b n . The bifurcation coefficients b 1 = a 1 , h(θ), and b 2 contain the terms u 2 , v 2 , . . . , u 6 , v 6 , as given in the Appendix (see also equations (29) and (30)). These terms are determined by solving equations of the form (cf. (26)) Here L 0 (q c ) is the singular matrix (3) from the linear calculation, c n is a constant, and f n (θ) is nonzero; for example, f 3 (θ) = 1 + 2 cos θ. Multiplying the above equation by the left null-vector (ũ,ṽ) of L 0 (q c ) gives uu n + Kṽv n = 0 at the degeneracy. Using this relation and relation (37), we see that From the definitions of the null-vectors (u 1 , v 1 ) and (ũ,ṽ) it immediately follows that u n η 1 + v n η 2 = 0, n = 2, . . . , 6.
These expressions may then be substituted into the bifurcation coefficient expressions from the Appendix, giving

Hexagons 3β+α=0
Stripes α β 15β+2α=0 Figure 2: Stability of patterns at the hexagonal degeneracy γ = 0. The details of f (u, v) and g(u, v) in (1) enter the stability analysis via two effective parameters: α, which depends on the quadratic nonlinearities, and β, which depends on the cubic nonlinearities. Rhombs and super hexagons are always unstable, as shown in Section 4.1.
All patterns bifurcate unstably for β > 0 and/or 15β + 2α > 0. For the region labeled "hexagons" we have determined only that hexagons are neutrally stable, since their stability for γ = 0 depends on terms higher than cubic in (16).
for γ = 0. In short, the quadratic nonlinearities from F 2 (u, v) and G 2 (u, v) completely disappear from the bifurcation coefficients b 1 and h(θ), and appear in b 2 solely through the term This further implies that all θ-dependence disappears from the bifurcation equations, and that rhombs and super hexagons are always unstable at the bifurcation point, for the degenerate problem. The signs of the eigenvalues for rhombs and super hexagons are given in Table 2. The first two eigenvalues for rhombs depend on the signs of b 1 + b 4 and b 1 − b 4 , which cannot both be negative. Hence rhombs are always unstable. The signs of the first three super hexagon eigenvalues are determined by the signs of All three cannot be made negative, since e 1 = −(3e 2 + 8e 3 ). Stripes and hexagons may bifurcate stably; Figure 2 gives their stability assignments in the (α, β)-parameter plane. Note that, for the cubic truncation of the bifurcation problem, hexagons are only known to be neutrally stable, since the eigenvalue −γx + O(x 3 ) in Table 2 depends on higher-order terms, which are not computed here. For the corresponding square lattice computation, stripes are stable for β < 0; all other solutions are unstable for γ = 0. This picture can be made even simpler for systems with F 2 (u, v) ∝ G 2 (u, v), such as the CIMA model considered in the next section. In these systems, γ = 0 typically implies F 2 (u 1 , v 1 ) = G 2 (u 1 , v 1 ) = 0, leading to α = 0 in Figure 2. Stripes are then stable provided β < 0. We next unfold the degenerate problem by considering 0 < |γ| ≪ 1 in (20), keeping b 1 = 3β, b 2 = α + 6β, and b 4 = b 5 = b 6 = 6β. There is now a small range of values µ ∈ (µ 1 , µ 2 ) for which super hexagons (or triangles) are stable, provided that α < 0, β ∈ (α/3, −α/21). For example, Figure 3 shows the unfolded bifurcation diagram for β = −1/3, α = −2 (0 < γ ≪ 1). It is straightforward to construct a system of the form (1) with these values. Let g(u, v) = 4.05u − 2.8v − 0.04u 2 + uv.
This system gives γ ≈ 0.03, α ≈ −2, β ≈ −1/3, and according to the unfolding calculation gives stable super hexagons (or triangles) in the range 3.07 × 10 −5 < λ a < 1.38 × 10 −4 . A value of λ a = 7.6 × 10 −5 corresponds to µ = 5.54 × 10 −4 in Figure 3, at which point rolls, hexagons, and super hexagons are all stable according to our analysis. Figure 4 shows a numerical integration of this system with λ a = 7.6 × 10 −5 , using a  with reaction-kinetics given by (39). The sole difference between these numerical runs was in the initial condition. Compare with the bifurcation diagram given in Figure 3. Numerical simulations were performed with a pseudo-spectral Crank-Nicholson scheme, using second-order Adams-Bashforth on the nonlinear terms. The rectangular domain has aspect ratio √ 3, with a 128×128 grid.
pseudo-spectral Crank-Nicholson scheme on a 128×128 rectangular grid of aspect ratio √ 3. All three pictures are at the same parameter values, and differ solely in the initial condition. The amplitude of the fundamental mode in each case is within two percent of the steady-state values predicted by the amplitude equations.

Lengyel-Epstein model
The results of the previous sections will now be applied to a reaction-diffusion system proposed by Lengyel and Epstein [6] as a reduced model for the CIMA reaction: Note that the nonlinear terms differ solely by a multiplicative constant, which will further simplify the hexagonal lattice computations. Our analysis assumes a steady-state solution of u = v = 0. Thus we let After expanding the nonlinearities to cubic order the equation may be written as after rescaling time and space to absorb a factor of 1 + x 2 0 . The coefficients in (41) are related to those in (40) by The overall factor of δb in the v t -equation will later be removed from the bifurcation problem by rescaling time. That is, although the overall factor of δb enters the linear calculations -in particular the necessary conditions (9) for Turing pattern formationit simply scales out of the final bifurcation problem. In addition to the two remaining equation parameters,â and K, the lattice angle θ will enter our bifurcation analysis.  Table 1 and evaluating at different points in theâ − θ plane. The region near lattice angle θ = π/3 has been removed, since this angle gives a hexagonal interaction,. Note also that atâ ≈ 9.4 rolls are stable for all θ. This value ofâ yields the degeneracy γ = 0 in the hexagonal bifurcation problem (16).
Choosing K to fix the system at the Turing bifurcation point, the resulting bifurcation scenario is then determined byâ and θ, where, for example, θ is related to the choice of square lattice by (17).
We first consider the relative stability of patterns which lie on a square lattice. Figure 5 is a plot in the (θ,â) plane, generated by calculating the bifurcation coefficients a 1 , a 2 = h(π/2), a 3 = h(θ), and a 4 = h(π/2 − θ), as given in the Appendix. These coefficients are then substituted into the eigenvalue expressions in Table 1, generating the stability assignments in the figure. The region surrounding θ = π/3 has been removed since h(θ) diverges as θ → π/3; the hexagonal lattice analysis is required at this point. Also, in the vicinity of θ = π/3 the domain of validity of the bifurcation results is very small; certain slaved modes are only weakly damped, leading to a smalldivisor problem in the computation of h(θ) (see Figure 6 for an example). Finally, we emphasize that anti-and super square patterns only exist for a discrete set of θ values, as discussed in Section 2. Figure 6: Critical modes (2, 1) and (−2, 1) on a (2, 1)-square lattice interact nonlinearly to generate a mode (0,2) that lies very near to the critical circle. We expect that this proximity will greatly reduce the domain of validity of the bifurcation calculation.
A number of general conclusions may be drawn from Figure 5. At all values ofâ, simple squares are unstable for some range of θ-values. We therefore conclude that they are unstable, near onset, in the unbounded domain. Rolls are similarly unstable for allâ with the exception ofâ ≈ 9.4, where rolls are stable to perturbations in any direction θ. Below we discuss the significance of this particularâ-value. Figure 7 shows numerical simulations at two points in Figure 5. A square box size of L = √ 13(2π/q c ) gives critical lattice vectors of (±3, 2)(2π/L) and (±2, 3)(2π/L) (cf. Figure 1). A (3, 2)-lattice generates mode angles of θ 3 = cos −1 12/13 ≈ π/8 and θ 4 = π/2 − θ 3 ≈ 3π/8. Figure 5 suggests that for θ = 3π/8 rolls become unstable to rhombs atâ ≈ 9.8. The rolls in Figure 7 were generated from random initial conditions atâ = 9.74; the rhombs were generated from random initial conditions atâ = 9.8. For these simulations δ = 1,b = 5, and the numerics were performed quite close to the bifurcation point, using a value of K = K c + 2.9 × 10 −4 . These numerical results are in fact sensitive to the size and shape of the computational domain; the choice of the computational domain damps modes that might otherwise destabilize the rhombs.
In Figure 5, rolls appear to become stable for all θ at the pointâ ≈ 9.4. This is the point at which the hexagonal degeneracy condition γ = 0 is satisfied; since the nonlinear  terms differ by a constant, the only way to satisfy γ ≡ũF 2 (u 1 , v 1 ) +ṽG 2 (u 1 , v 1 ) = 0 in this model is to take F 2 (u 1 , v 1 ) = G 2 (u 1 , v 1 ) = 0, which occurs at At this point, hexagons are unstable to rolls (cf. Figure 2 for α = 0, β < 0). The results summarized by Table 2 may be used to determine certain aspects of the hexagonal bifurcation problems near the degeneracy γ = 0. In section 4.1 it was pointed out that the unfolding scenario is fairly simple; in this model, the unfolding computation is even simpler. Since F 2 (u 1 , v 1 ) = 0, the parameter α is zero. And since we know from section 4.1 that there is no θ-dependence, the unfolding is described by the single bifurcation diagram given in Figure 8.

An Example of Super Squares
Section 2 described two types of superlattice patterns for square lattices: super squares and anti-squares. Figure 9 gives an example of each of these patterns for (m, n) = (5, 2), corresponding to θ 3 = cos −1 (20/29). Each pattern has a π/2 rotational symmetry. Super squares are invariant under reflections of the square domain; anti-squares have symmetries involving both translations and reflections (see [9] for details).
In Section 4 a system was constructed which supported stable super hexagons; by using the Appendix coefficient expressions, coupled with the stability information in Table 1, we may similarly construct a system which supports either super or anti-  Figure 10 shows a numerical integration of this system at a = 1.88, just above onset, showing both the initial condition and final steady state.

Conclusions
From the general two-component reaction-diffusion system (1), we have derived analytic expressions for the coefficients of the leading nonlinear terms in the bifurcation equations for rhombic, square and hexagonal lattices. These coefficients allow us to calculate the relative stability of patterns which are periodic on either square lattices (stripes, squares, rhombs, and super squares or anti-squares) or hexagonal lattices (stripes, simple hexagons, rhombs, super hexagons or super triangles) at the onset of the Turing instability. In particular, the dependence of stability upon system parameters may be calculated explicitly for specific reaction-diffusion models, as we demonstrated for several model equations, including the Lengyel-Epstein CIMA model.
One of the surprising results of our general analysis is that details of the reactionkinetics enter the computation of the bifurcation coefficients in a very limited fashion. In the case of the degenerate hexagonal bifurcation problem we show that the coefficients of the cubic terms in the bifurcation problem depend on the reaction-kinetics through just two effective parameters, that are simple to compute. Moreover, we find that at this degeneracy the angle dependence, which is expected for the rhombic lattice problems, drops out completely, and that rhombs and super hexagons never bifurcate stably at the Turing bifurcation point of two-component reaction-diffusion models. However, in the unfolding of this degenerate bifurcation problem we find that super hexagons may co-exist stably with simple hexagons and stripes for a limited range of the bifurcation parameter. This tri-stability property is demonstrated by numerical simulation of a reaction-diffusion system in the vicinity of the degenerate bifurcation point.
We plan to carry out more extensive numerical investigations that test some of the predictions of our bifurcation analysis. This will enable us to investigate the domain of validity of our bifurcation analysis, which is restricted to small-amplitude spatiallyperiodic Turing patterns. We also intend to determine the nature of the curious behavior at the hexagonal degeneracy, and to determine to what extent parameter collapse is exhibited by systems with three or more components, e.g. for Turing patterns in activator-inhibitor-immobilizer models [21].