Spatiotemporal chaos and quasipatterns in coupled reaction-diffusion systems

In coupled reaction-diffusion systems, modes with two different length scales can interact to produce a wide variety of spatiotemporal patterns. Three-wave interactions between these modes can explain the occurrence of spatially complex steady patterns and time-varying states including spatiotemporal chaos. The interactions can take the form of two short waves with different orientations interacting with one long wave, or vice-versa. We investigate the role of such three-wave interactions in a coupled Brusselator system. As well as finding simple steady patterns when the waves reinforce each other, we can also find spatially complex but steady patterns, including quasipatterns. When the waves compete with each other, time varying states such as spatiotemporal chaos are also possible. The signs of the quadratic coefficients in three-wave interaction equations distinguish between these two cases. By manipulating parameters of the chemical model, the formation of these various states can be encouraged, as we confirm through extensive numerical simulation. Our arguments allow us to predict when spatiotemporal chaos might be found: standard nonlinear methods fail in this case. The arguments are quite general and apply to a wide class of pattern-forming systems, including the Faraday wave experiment.


Introduction
Two substances that react and diffuse can form patterns, an insight first highlighted in the work of Alan Turing [1,2]. Motivated by an interest in embryonic morphogenesis, Turing studied discrete and continuum models for the spontaneous emergence of structure in a ring of cells. Depending on the details of the reaction, a Turing-type system may have a stable, spatially-uniform steady state in the absence of diffusion. Two fundamental instabilities may occur. One possibility is a Hopf bifurcation leading to temporal oscillations with a preferred wavenumber of zero. The other possibility, driven by diffusion, is a bifurcation to a steady spatial pattern (Turing pattern) with non-zero wavenumber, typically stripes or hexagons.
The first laboratory experiment to produce a Turing pattern came nearly 40 years after Turing's original work: Ref. [3] reports the observation of patterns in the chlorite-iodide-malonic acid (CIMA) chemical reaction. Since the seminal discoveries of [1,3], there has been a vast literature on reaction-diffusion patterns and their applications, which include animal skin pigmentation [4], the cerebral cortex [5], vegetation ecology [6], plankton colonies [7], and many others [8,9].
A variation on the classic reaction-diffusion system is the so-called coupled (or multilayered) system, in which two or more reaction-diffusion systems are connected together so that they may influence each other. Because of the additional degrees of freedom in these coupled systems, they are an amenable setting in which to investigate how competing instabilities affect pattern formation. Coupled systems can produce a variety of states including simple Turing patterns, standing waves, mixes of Turing patterns and spiral waves, square and hexagonal superlattice patterns, and many more [10,11]. Coupled systems are important in biology, especially in neural, ecological and developmental contexts; see [12] for examples and for an overview of selected results.
Unfortunately, it is difficult to manipulate experiments on the aforementioned biological systems. A common approach to studying coupled reactiondiffusion systems, then, is to study a paradigmatic chemical experiment in the laboratory, as in [13] (see Figure 1). In this work, the experimentalists set up two thin gels, within each of which the chlorine dioxide-iodine-malonic acid (CDIMA) reaction takes place. They put the gels in contact and controlled the strength of coupling between the two layers by modifying the properties of a membrane placed at the interface, resulting in different patterns. To complement laboratory experiments, investigators have studied a host of nonlinear partial differential equation (PDE) models, including the Lengyel-Epstein model of the CIMA reaction [14] and Brusselator model of a generic trimolecular reaction [15].
The work of [16] included a theoretical study of two reaction-diffusion systems coupled together in a parameter regime near a codimension-two Turing-Turing bifurcation point. This work demonstrated that by changing the interlayer coupling strength, one can manipulate the ratio of the length scales associated with two resonantly interacting Turing instabilities and encourage the formation of certain complex patterns in the Brusselator model. In our present work, we will also carry out a theoretical investigation of coupled reactiondiffusion systems and we will also focus on resonant mode interactions. However, in contrast to the set-up in [16], we will use the within-layer diffusion constants as control parameters. In experiments, one could manipulate these diffusion constants by changing properties of the medium of each layer.
Our work complements a robust literature that has examined the role of three-wave interactions in the Faraday system, in which a layer of fluid is vertically vibrated in a time-periodic fashion, potentially producing standing wave patterns. Patterns with two dominant length scales, including quasipatterns and superlattice patterns, have been observed in many Faraday wave experiments [17][18][19][20][21]. The theory of Faraday three-wave interactions was developed in [22][23][24][25][26], among other sources. Much of this body of work took the following approach. Based on symmetry considerations, one can write down amplitude equations describing the slow-time evolution of modes close to a codimensiontwo point where all waves associated with two different length scales are neutrally linearly stable. By detuning from that point and assuming that one of the sets of waves is weakly damped, one can perform a centre manifold reduction and assess the role that the weakly damped mode has on the dynamics of the other modes. At a granular level, this influence is seen as a (potential) contribution to coefficients of cubic terms in the amplitude equations for the primary pattern modes. The leading order influence is determined by quadratic terms in the original amplitude equations.
Our present study focuses on the role of three-mode or three-wave interactions and, pivotally, builds on, clarifies and extends the main ideas of [27].
When there are two (nearly) critical length scales that are not too disparate, two of the shorter wavelength modes with different orientations can interact with one of the longer ones, or two of the longer wavelength modes can interact with one of the shorter ones. In each case, the orientations of the modes are determined by the requirement that two longer wavevectors add up to a shorter one, or that two shorter wavevectors add up to a longer one. Pattern formation can be strongly dominated by these interactions. Rather than slaving away one set of critical modes and studying cubic terms, as described above, we instead see how much understanding may be gleaned by restricting our attention to quadratic terms near the codimension-two point. This approach, namely, studying the effect of three-wave interactions on spatiotemporal pattern formation in reaction-diffusion systems by looking at quadratic coefficients, has proven successful in the past [27,28]. Our present work develops a more exhaustive investigation in the context of layered Turing systems, though the ideas are applicable wherever a pattern-forming system can have two unstable length scales, including the Faraday wave experiment.
The rest of this paper is organised as follows. In Section 2 we outline the basic nonlinear three-wave interactions in the case of pattern formation with two competing wavelengths, and in Section 3 discuss the role of the quadratic coefficients (and in particular their signs) in influencing the resulting patterns. Section 4 presents the two-layer Brusselator model, and Sections 5 and 6 describe the linear and weakly nonlinear theory of the model. Numerical results appear in Section 7, and we conclude in Section 8.

Nonlinear three-wave interactions
We first consider patterns in the variations of a real scalar field U (x, y, t).
Assume the system forms patterns with two distinct length scales. More specifically, and without loss of generality, we assume that waves with wavenumbers k = 1 and k = q (q < 1) become unstable and have growth rates r 1 and r q respectively. At onset, the pattern U (x, y, t) will contain a combination of Fourier modes e ik·x , with |k| = q or |k| = 1. We write, close to onset, where q j are wavevectors on the circle |k| = q, with mode amplitudes w j (t), and k j are wavevectors on the circle |k| = 1, with mode amplitudes z j (t). The overall pattern U is real, so waves come in equal and opposite pairs with complex conjugate amplitudes. The time evolution of the complex mode amplitudes is influenced by nonlinear combinations of other mode amplitudes. The particular combinations that arise are determined by the lengths and orientations of the wavevectors, in a manner that can be explained by focusing on one mode on each circle and examining the lowest-order combinations that influence the chosen mode.
w 1 w 9 w 8ż 1 = Q ww w 6 w 7 w 1 = Q wz (w 8 z 8 + w 9 z 9 ) Figure 2: Three-wave interactions with two wavenumbers k = 1 (outer circle) and k = q (inner circle) that influence the evolution of z 1 (left column) and w 1 (centre column). Each vector is labelled with the amplitude (z 1 , w 1 , . . . ) of the corresponding mode. First row: three wave vectors of the same length (three long or three short). Middle row: two long wave vectors and one short, defining an angle θz = 2 arccos(q/2). Bottom row: one long wave vector and two short, defining an angle θw = 2 arccos(1/2q). This last case only occurs when q > 1 2 . In all cases, the right column gives the quadratic terms in the amplitude equation that result from the three-wave interactions depicted to the left.
The two modes we choose are z 1 (t)e ik1·x and w 1 (t)e iq1·x , as well as their complex conjugates, illustrated in Figure 2. We will develop an ordinary differential equation (ODE) for each mode amplitude and express it as a truncated Taylor series. The linear terms in the evolution equation for z 1 and w 1 are simply r 1 z 1 and r q w 1 , respectively, and the starting point for the ODEs describing the evolution of each mode amplitude iṡ z 1 = r 1 z 1 + nonlinear terms, w 1 = r q w 1 + nonlinear terms. ( Nonlinear functions of U , as written in (1), will involve products of modes and therefore sums of wavevectors. The combinations of modes that influence z 1 and w 1 will be those whose wavevectors add up to k 1 and q 1 respectively. The lowest-order nonlinear terms are quadratic, arising when two vectors (of length 1 or q) add up to k 1 or q 1 . The simplest interactions involve modes at 60 • . The wave vectors in these so-called hexagonal states can be arranged in an equilateral triangle; see Figure 2, top row. If k 1 = k 2 + k 3 (all of length 1), and q 1 = q 2 + q 3 (all of length q) then the equations forż 1 andẇ 1 will have the terms Q zh z 2 z 3 and Q wh w 2 w 3 , where z 2 , z 3 , w 2 and w 3 are the amplitudes of modes with wavevectors k 2 , k 3 , q 2 and q 3 respectively, and Q zh and Q wh are coefficients.
As well as equilateral triangles, one may have isosceles triangles with one short and two long sides ( Figure 2, middle row) and triangles with one long and two short sides ( Figure 2, bottom row). The latter case can only happen if q > 1 2 . The two isosceles triangles define related angles θ z = 2 arccos(q/2), θ w = 2 arccos(1/2q), as seen in Figure 2 [27]. These triangles lead, in different combinations, to contributions indicated in the right column of Figure 2, where Q zw , Q zz , Q ww and Q wz are further coefficients. The mode amplitudes are numbered in order of appearance in Figure 2. The end result is that, at quadratic order, there are 8 modes that couple to each of z 1 and w 1 : These 16 additional modes, 8 with wavenumber 1 and 8 with wavenumber q, will each couple to up to 8 further modes, and each of these further modes will couple to up to 8 more, as so on, as outlined in [27]. One might ask, "where does it all end?" The answer depends on q, as explained in [27]. For q < 1 2 , two short vectors added together do not extend to the outer circle, and so the interactions in Figure 2 (bottom row) do not exist, and the end result is, for example, six modes on the inner circle and 12 on the outer, as in Figure 3(a). This case can lead to superlattice patterns [18,21,23] or quasipatterns [29]. For q = 2 sin π 12 = 2 − √ 3 = 0.5176, the angles θ z and θ w (defined in Figure 2) are 150 • and 30 • respectively, and so all possible threewave interactions can be accommodated within a set of 12 vectors of length 1 interleaved with 12 vectors of length q, as in Figure 3(b). This special value of q is the only one in the range 1 2 < q < 1 where three-wave interactions generate a finite number of modes [27]. For all other 1 2 < q < 1, an infinite number of modes is generated, as illustrated in Figure 3(c) for a generic choice of q.
Of course, equations (2) and (4) go only up to quadratic order. At cubic order, every wave on the two circles couples to every other wave, since k 1 = k 1 + k j − k j = k 1 + q j − q j , for any vectors k j and q j . Given a finite set of modes as in Figure 3(a,b), one can work out the amplitude equations, calculate quadratic and cubic (and higher if needed) coefficients, and analyse which solutions are possible and stable. Doing this for complex periodic patterns is challenging because dozens of modes are involved. For quasipatterns, there is the additional complication that this process, where the small-amplitude pattern is expressed as a power series in a small parameter, leads to divergent series [30,31], though existence of quasipatterns has been proved in the Swift-Hohenberg equation [32] and in Rayleigh-Bénard convection [33]. The case of a potentially infinite set of modes ( Figure 3c) is challenging.
The purpose of our present work is to see how far considerations from just the quadratic level can help understand the outcome when both circles in Fourier space are (potentially) fully occupied. The example we use to illustrate the ideas is the two-layer Brusselator model, in Section 4. Before discussing this model, we review what is known about the role of three-wave interactions in the formation of complex spatiotemporal patterns and outline our hypotheses regarding how quadratic coefficients would influence observed patterns. Figure 3 gives the three qualitatively different possible cases. If q < 1 2 , all three-wave interactions can be accommodated within a set of 6 vectors of length q and 12 vectors of length 1, as in Figure 3(a). This leads to 9 coupled complex amplitude equations. The resulting patterns are spatially periodic when cos θ z and √ 3 sin θ z are both rational, which happens for a dense but measure zero set of q. Otherwise, the resulting patterns are quasiperiodic [29]. The second case, as in Figure 3(b), with q = 2 − √ 3 = 0.5176, leads to 12 vectors of length q and 12 vectors of length 1, necessitating 12 complex amplitude equations. In these two cases, with a finite number of amplitude equations (which will be explored in more detail elsewhere), standard nonlinear methods can be employed to obtain equilibrium points, and (to some extent) their stability, bifurcations, and so forth. In the third case, with 1 2 < q = 2 − √ 3 < 1 as in Figure 3(c), three-wave interactions lead to coupling between an infinite number of modes, and so there is the possibility of an infinite number of amplitude equations. In this scenario, it is not clear that standard nonlinear methods will yield useful information.

Role of the quadratic coefficients
All three cases involve sets of interacting waves, with the strongest interactions happening between groups of three. We illustrate a single set of three interacting waves by taking two outer vectors coupling to an inner one, with wavevectors k 6 + k 7 = q 1 and amplitudes z 6 , z 7 , w 1 , as in Figure 2 (middle row, middle column). The amplitude equations in this case are of the form: Here, A z , A zz , A zw , A w and A wz are cubic coefficients that depend on the details of the problem and that can in principle be calculated from governing equations.
Porter and Silber [34] investigated (5) in detail and found that the dynamics depends on the product of quadratic coefficients Q zw Q zz , as well as the linear and cubic coefficients. Typically, when Q zw Q zz is positive, there are stable equilibria and no time-dependent states. On the other hand, when Q zw Q zz is negative, in addition to stable equilibria, time-periodic solutions and chaotic solutions are possible via Hopf and global bifurcations. In the positive case, the z and w modes can act to reinforce each other, while in the negative case, there can be time-dependent competition between z and w modes. The same conclusion applies equally to the three-wave interaction between two w and one z mode ( Figure 2, bottom row, left column). Here, the relevant combination of quadratic coefficients is Q wz Q ww .
These considerations led Rucklidge et al. [27] to hypothesise how the combinations of quadratic coefficients would influence patterns, essentially supposing that the qualitative conclusion of [34] applies also when there are many sets of interacting waves, even though each individual wave has three-wave interactions with several combinations of modes, as discussed above. Steady patterns should be expected when Q zw Q zz > 0 and Q wz Q ww > 0, and time-dependent patterns should be possible when one or both pairs of quadratic coefficients are of opposite sign. When q > 1 2 , complex patterns, with modes at many different orientations, as in Figure 3(c), may be possible.
Before developing these ideas further, for the purposes of this paper, we distinguish between different types of patterns. Simple patterns are stripes, hexagons (or symmetry-broken hexagons) of either critical wavelength. There are also rhombs, here taken to mean patterns with two modes of equal amplitude on one circle coupled to a third mode on the other circle. Superlattice patterns are dominated by 12 modes at one wavenumber and 6 at the other; here we blur the distinction between spatially periodic superlattice patterns and quasipatterns [29]. With q < 1 2 , there is only one type of superlattice pattern, while with 1 2 < q < 1, there are two types, with six modes on one circle and twelve on the other, either way around. Regular twelve-fold quasipatterns have 12 modes (equally spaced) at each wavenumber, as illustrated in Figure 3(b). The patterns discussed so far may have defects, which can evolve over long timescales [35]. Complex patterns have large numbers of modes, at both wavenumbers, coupled through three-wave interactions, as illustrated in Figure 3(c), but are not simple patterns with defects as defined here. Time dependent (periodic, chaotic) versions of each of these types of patterns are also possible, evolving over shorter timescales. We reserve the term spatiotemporal chaos for the case when complex patterns have persistent chaotic dynamics with many positive Lyapunov exponents as in [36].
With this classification in mind, we extend the hypotheses of [27] as detailed in the points below, and as summarised in Table 1. Here by finding a pattern, we mean that there are combinations of r 1 and r q where that pattern is an asymptotic state obtained when starting from random initial conditions in a domain large enough to accommodate a wide range of wavevector orientations on the two critical circles.
• In all cases, we expect to find steady simple patterns such as stripes and hexagons, possibly with broken symmstry or with defects.
• In addition, with q < 1 2 , we expect to find steady superlattice patterns with wavevectors as in Figure 3(a). We may also find rhombs. If Q zw Q zz is negative, we expect to find time-dependent superlattice patterns (and rhombs) with the same wavevectors, and also spatiotemporal chaos, with all wavevectors on the two circles being active. We do not expect to find steady complex patterns.
• With q > 1 2 , we expect to find both types of steady superlattice patterns. We also expect to find steady complex patterns, with large numbers of wavevectors on both circles. The combinations of quadratic coefficients relevant to the two superlattice cases are Q zw Q zz and Q wz Q ww respectively: if the relevant combination is negative, we expect to find time- dependent superlattice patterns with the same wavevectors. If either or both combination is negative, we expect to find spatiotemporal chaos, with all wavevectors on the two circles being active. If Q zw Q zz and Q wz Q zz are both negative, we expect to find time dependence more readily. In general, we expect to find spatiotemporal chaos more readily than in the q < 1 2 case.
• For the special value q = 2 − √ 3 = 0.5176, we expect to find steady twelve-fold quasipatterns with wavevectors as in Figure 3(b). If one or both of Q zw Q zz or Q wz Q ww is negative, we expect to find time-dependent quasipatterns and spatiotemporal chaos.
These considerations neglect the roles that the hexagonal quadratic coefficients Q zh and Q wh might play.

Two-layer Brusselator model
The Brusselator [15,37] is a canonical model of a reaction-diffusion system. More specifically, it describes an autocatalytic chemical reaction, The products D, E are generally not of interest because they do not enter into the autocatalysis. Therefore, we restrict attention to the reactants X, Y, A, B.
In the Brusselator, it is assumed that A, B are present in great excess, and thus can be treated as constants. Allowing for spatial diffusion, and using the standard theories of reaction kinetics, we write down rate laws for X, Y as the differential equations Through abuse of notation, we have now let A, B, X, Y ≥ 0 represent concentrations of these chemicals rather than symbolising the chemicals themselves. Here, X, Y are time and space dependent chemical concentrations and, as assumed, A, B are constant. Eq. (7) has a spatially homogeneous steady state solution, namely X = A, Y = B/A. We adopt shifted coordinates for the dependent variables, letting X = A + U , Y = B/A + V so that the equilibrium becomes the trivial one, U = 0, V = 0. In these coordinates, (7) is The chemical concentrations are U (x, t) and V (x, t), where x is the planar spatial coordinate x = (x, y). The diffusion constants have been relabelled for clarity of notation, that is, D U = D X and D V = D Y . As in [10,16], we consider a two-layer Brusselator model. The layers are coupled together "diffusively," manifesting as linear terms with coefficients α, β ≥ 0: Here, U 1,2 (x, t) and V 1,2 (x, t) are chemical concentrations in each layer. For convenience, we have used shorthand to represent the nonlinear terms, We have assumed that A and B do not vary across layers, meaning that each excess reactant is present in the same amount in each layer. For all calculations in the remainder of this paper, we take A = 3 and B = 9 as our standard parameter values, as chosen in [10] to model the CIMA reaction.

Linear theory
If we drop the nonlinear terms in (9), we can solve the resulting linear PDE in terms of modes e ik·x that grow as e σt . Here, σ is a growth rate that depends on the wavenumber k = |k|. The linear problem is represented by a 4 × 4 Jacobian matrix J, The growth rates σ are the eigenvalues of J and satisfy the characteristic equation where the coefficients C 0 , . . . , C 3 are (cumbersome) polynomial functions of nine parameters: A, B, D U1 , D V1 , D U2 , D V2 , α, β and the wavenumber k.
If we were to choose D U1 = D U2 and D V1 = D V2 , as done in [16], then the 4 × 4 matrix J can be decomposed in to two 2 × 2 parts. However, in this case, it turns out that two of the quadratic coefficients vanish in the weakly nonlinear theory. To avoid this degeneracy, we take an alternative approach, allowing the diffusion constants to be different in the two layers. As mentioned in Section 4, the experimental context is that we take the chemistry to be identical in the two layers (meaning A, B do not depend on layer) but take the substrates to be different, so their diffusion properties will be different.
Rather than fix the values of the parameters, we are aiming to explore the range of outcomes close to the codimension-two point where patterns with two length scales are simultaneously unstable, for a range of values of the wavenumber ratio. Therefore, we seek parameter values for which σ, when viewed as a function of k, takes on certain values at local maxima. For example, if σ has a local maximum at k = 1 and k = q, for some choice of q, then four conditions must be satisfied: σ = r 1 at k = 1, σ = r q at k = q and dσ dk = 0 at k = 1, q. The derivative dσ dk can be obtained from (12) by differentiating with respect to k: The four conditions result in four equations for the nine parameters listed above (but with k replaced by q), with two additional parameters σ(1) = r 1 and σ(q) = r q . This means that seven of the parameters can be specified, and four found by solving the equations. For example, we take our base parameter values A = 3 and B = 9 [10] and choose q = 2 − √ 3 = 0.5176, appropriate for twelve-fold quasipatterns [38,39]. Additionally, we choose r 1 = r q = 0 in order to be at the codimension-two point, and we choose α = β = 1. Recall that α and β control the diffusion of the two chemicals between the two layers, while D U1 , D U2 , D V1 and D V2 control the diffusion of the chemicals within each layer.   Table 2 and in full in [40].
For this choice of seven parameters, the resulting four polynomial equations for the four remaining unknowns (D U1 , D U2 , D V1 and D V2 ) can be worked out; the simplest (shortest) of these is The coefficients in this (and the other three equations) depend on the choice that we made for A, B, α, β, q, r 1 and r q .
We solve the four polynomial equations numerically, using Bertini [41]. by relabelling) for 0.32 ≤ α ≤ 7.99. We plot the four diffusion coefficients as functions of α for β = 1 in Figure 4 in black for q = 2 − √ 3 = 0.5176, with vertical lines indicating how the diffusion coefficients vary with q, keeping A, B, r 1 and r q fixed. Sample numerical values of the diffusion constants for different choices of α and q are given in Table 2 (below) and in full in [40].
We observe from Figure 4 that the range from the smallest to the largest values of the diffusion constants appears to diverge as α approaches 8. The same happens for other choices of q. In addition, the ordering of the diffusion constants changes around α = 7: for α < 7, we have D U1 < D V1 and D U2 < D V2 , which seems experimentally reasonable, in that one chemical diffuses slower than the other in either substrate, while for α > 7, this is not true. For later calculations, we will choose α to vary between 1 and 7. The experimental relevance of the larger values of α should be treated with caution.
We conclude our discussion of the linear theory with a sample dispersion relation (σ(k) plotted as a function of wavenumber k) in Figure 5, for q = 2 − √ 3 and for a range of α, at the codimension-two point r 1 = r q = 0. The eigenvalue is maximal at k = q and k = 1, but the minimum for q < k < 1 is about −0.1 for α = 1, and only about −0.015 for α = 7.
In this dispersion relation, we control the separation and heights of the growth-rate maxima by varying q, r 1 and r 1 , and solving for the four diffusion coefficients. For smaller q, the peaks are well separated and reasonably sharp, while for larger q the peaks are closer together, broader and the depth of the minimum is less. For this reason, we limit ourselves to q ≤ 0.66. Since we want to keep an interval of negative growth-rate between k = q and k = 1, we (mostly) limit ourselves to r 1 ≤ 0.01 and r q ≤ 0.01.

Weakly nonlinear theory
Once the uniform state U 1 = V 1 = U 2 = V 2 = 0 becomes linearly unstable, solutions will grow exponentially until nonlinear effects become important. The first of these are the three-wave interactions. The weakly nonlinear theory is standard [28,42,43], though made more complicated here because of the codimension-two bifurcation and because there are four scalar fields in (9). We are concerned only with the leading order effect of three-wave interactions, and so we need to compute only up to second order in the weakly nonlinear theory.
We write and write the PDE (9) as where L is a linear operator representing the linear terms, and all the nonlinear terms in (9) are in NLT(u). To explore the properties of solutions close to u = 0, we introduce a small parameter 1, and we expand u in powers of : Recall that in Section 5, we computed values of D U1 , D V1 , D U2 and D V2 such that the linear operator L had zero eigenvalues (r 1 = r q = 0) at two wavenumbers, k = 1 and k = q, at given values of A, B, α and β. We now suppose that the linear operator is perturbed by an order amount so that the growth rates r 1 at k = 1 and r q at k = q are order . In practice we perturb D U1 , D V1 , D U2 and D V2 and do it in such a way that there are local maxima in the growth rate remain at k = 1 and k = q. We can scale r 1 → r 1 and r q → r q and write the linear operator L as where L 0 is a singular linear operator, and L 1 is the largest part of the perturbation of the linear operator from L 0 . Finally, we scale time so that d/dt → d/dt. With these choices of scaling, the time derivative, the linear terms, and the lowest-order nonlinear terms all appear at the same order. Substituting into (17), we have where NLT 2 represents the quadratic nonlinear terms.
The operator L 0 is singular: L 0 e ik·x v 1 = 0 whenever |k| = 1, and L 0 e iq·x v q = 0 whenever |q| = q, where v 1 and v q are the eigenvectors of the zero eigenvalues of the Jacobian matrix (11), with k replaced by 1 and q respectively. We normalise the eigenvectors so that v 1 · v 1 = 1 and v q · v q = 1. Following the example of Section 5 and (15), with A = 3, B = 9, α = 1 and β = 1, we find With these eigenvectors, the general solution to L 0 u 1 = 0 is similar to the expression in (1): where {q j } and {k j } are arbitrary sets of vectors on the two circles |q j | = q and |k j | = 1. Writing u 1 in this way solves the O( ) part of (20). (20) is Recall that L 0 is singular and so cannot simply be inverted to find u 2 as a function of u 1 . Thus, before solving for u 2 , a solvability condition must be imposed. The standard method is to define an inner product between vectorvalued functions f (x) and g(x) on the domain Ω of the problem: wheref is the complex conjugate of f and |Ω| is the area of the domain. We define L † 0 , the adjoint of L 0 , by requiring that for all f and g. We restrict to functions on Ω that satisfy periodic boundary conditions. In this case, the adjoint operator L † 0 is just the transpose of L 0 . Having defined L † 0 , we solve L † 0 e ik·x v † 1 = 0 and L † 0 e iq·x v † q = 0 to find the normalised adjoint eigenvectors v † 1 and v † q , with |k| = 1 and |q| = q. For our example, these are Then, for any u 2 , where k 1 and q 1 represent any vectors on the two critical circles. Thus, taking the inner products of e ik1·x v † 1 and e iq1·x v † q with (23) results in the solvability conditions Taking u 1 to be made up of waves with wavevectors from all the combinations of wavevectors in Figure 2 results in amplitude equations (including the values of the coefficients) up to quadratic order, as written in (4). We will take two specific examples, focusing only on the quadratic coefficients, and compute Q zh , Q zz and Q zw . For these, we need NLT 2 (u 1 ), which is (for A = 3 and B = 9): where (U 1 , V 1 , U 2 , V 2 ) are the four entries in u 1 .
To calculate the various quadratic coefficients described in Section 2, we take the combinations of wavevectors appropriate for each coefficient. For Q zh , we write where k 1 = k 2 + k 3 as in the top left panel of Figure 2, v 1 is the eigenvector as in (21) and c.c. stands for the complex conjugate. In this case, we have where we have highlighted the e ik1·x term, and v (1) 1 and v (2) 1 are the first and second entries in the vector v 1 in (21). There are similar expressions for U 2 2 and U 2 V 2 , involving v We have used the fact that v † 1 is real. Dividing by v † 1 · v 1 and matching to (4) results in an expression for Q zh . For the example set of parameters, Q zh = −0.7018. With r 1 defined to be the growth rate (on the slow time scale) of the wavenumber |k| = 1 modes, the linear term above is r 1 v † 1 · v 1 z 1 . Similar calculations but for different choices of wavevectors yield Q wh , Q zw and Q zz , and Q wz and Q ww . We illustrate with the calculation for Q zz and Q zw , and write where q 1 = k 6 + k 7 as in the middle row centre panel of Figure 2. In this case, we need the e iq1·x and e ik6·x components of U 2 1 and U 1 V 1 : again with similar expressions for U 2 2 and U 2 V 2 . The inner product with e ik6·x v † 1 in the first line of the solvability condition in (28) picks out the Q zwz7 w 1 term, while the inner product with e iq1·x v † q in the second line of the solvability condition picks out the Q zz z 6 z 7 term. These result in equations for the two quadratic Quadratic coefficients Q zz Q zw < 0 Q wz Q ww < 0 q = 0.5176 Figure 6: Weakly nonlinear theory for the two-layer Brusselator model (9), up to quadratic order. We take β = 1 in (9) and q = 0.5176 in (1). The six quadratic coefficients are Q zh (green), Q wh (black), Qzz (red), Qzw (magenta), Qwz (cyan), Qww (blue). In this case, the coefficients Qzw and Qzz have opposite sign for 2.33 < α < 5.00, and Qwz and Qww have opposite sign for 4.76 < α < 6.87. The values of diffusion coefficients are as in Figure 4. The data for this figure is available in full at [40], with selected values in Table 2.
For the example set of parameters, Q zw = −0.8974 and Q zz = −0.1997. Similar calculations yield Q wh = −0.5610, and Q wz = −0.2623 and Q ww = −0.9263 (only available since q > 1 2 ). Examples of the six quadratic coefficients as functions of α are shown in Figure 6, for q = 2 − √ 3 = 0.5176, and for β = 1, with numerical values for  Table 2: Sample values of the diffusion coefficients in (9) and the resulting quadratic coefficients in (4), with A = 3, B = 9, β = 1, r 1 = 0 and rq = 0, for different choices of q and α. The data are illustrated in Figures 4 and 6. A fuller version of this table (for 0.25 ≤ q ≤ 0.66) is available in full in [40]. this and other choices of q given in Table 2 and in [40].
With this choice of parameters, the coefficients Q zw and Q zz have opposite sign for 2.33 < α < 5.00, and Q wz and Q ww have opposite sign for 4.76 < α < 6.87. The behaviour of the quadratic coefficients for other values of q in the range 0.25 ≤ q ≤ 0.66 is similar: there is a range of α for which Q zw Q zz < 0, and (provided q > 1 2 ) there is a range of α for which Q wz Q ww < 0, where the ordering is the same throughout. The two ranges overlap over a limited range of α, centred on α ≈ 4.8 for all q.

Numerical results
Based on the linear and weakly nonlinear calculations in the previous sections, we have carried out a series of numerical simulations of the PDEs in (9). Our main goal is to explore the effect of varying the ratio of length scales, q, in regimes where we can control the signs of the quadratic coefficients. Our choice is to fix the diffusive coupling coefficient β = 1 and vary α with 1 ≤ α ≤ 7 (in steps of 1). With different choices of α, the two pairs of quadratic coefficients can have the same or opposite signs (see Figure 6), though the range where both pairs had opposite sign was very limited. We chose some special values of q, some less than and some greater than 1 2 : q = 1/ √ 7 = 0.3780, to encourage superlattice patterns [44]; q = 2 − √ 3 = 0.5176, to encourage twelve-fold quasipatterns [38,39]; q = 1/ √ 3 = 0.5774, to allow quadratic interactions between six modes on each circle; and q = 1 2 (−1 + √ 5) = 0.6180, to encourage ten-fold quasipatterns [45]. We also chose more "generic" values of q: 0.25, 0.33, 0.44, 0.55 and 0.66. The values of q and the corresponding angles θ z and θ w are listed in Table 3. All chemical properties are frozen with the choice of A = 3 and B = 9 as in [10].
The values of the diffusion coefficients at the codimension-two point r 1 = r q = 0 are given in Figure 4 and in [40]. For each selected case of q and α, we vary the diffusion coefficients to explore small positive and negative values of the two growth rates r 1 and r q . Specifically, setting (r 1 , r q ) = (r cos θ, r sin θ), we choose r = 0.01 (apart from data in Figure 11), with θ varying from 5 • to 355 • in steps of 10 • . For smaller q and α, these choices lead to growth rates σ(k) that are sharply peaked at k = q and k = 1, with a relatively deep negative minimum in between (see Figure 5). However, for larger q and α, the minimum between the two maxima is quite shallow, which means that, even with a small value of r = 0.01, there can be wide bands of unstable wavenumbers.
We start all simulations from small-amplitude random initial conditions in 16π×16π (8×8 of the shorter wavelengths) domains, except in the case when q = 2 − √ 3 where we also start simulations from a small-amplitude quasipattern initial condition. For parameter choices that do not result in a simple pattern, we explore the effect of a larger domain by re-running calculations in 60π × 60π (30 × 30 wavelengths) domains. Both 8 × 8 and 30 × 30 domains are appropriate for twelve-fold quasipatterns [46]. Time simulations are for at least 10000 time units: this is 100 growth times (for r = 0.01) and approximately three diffusion times for the larger domain when considering the smallest values of the diffusion coefficients.
We use 128 × 128 Fourier modes (using FFTW [47], the fastest Fourier transform in the West) in each direction for the 8 × 8 domains, and 512 × 512 Fourier modes for the larger 30 × 30 domains. We use the second-order exponential time differencing (ETD2) [48] scheme for timestepping, with a fixed timestep of 0.01. For this matrix exponential method, we split the linear part of the PDE (9) into diagonal and off-diagonal parts, and we treat the off-diagonal parts as nonlinear terms.
In all, we carried out over 4000 simulations, and the results we present below are an overview of the range of patterns we find. For α ≤ 3, we find a wide range of different patterns, but for α ≥ 4 we find simple patterns (hexagons) almost exclusively. Therefore, we focus on the cases with α = 1, 2 and 3. When α = 1, all quadratic coefficients are negative for all q (see Figure 6, Table 2 and [40]). Therefore, from Table 1, we expect to find only steady patterns. For α = 2, Q zw and Q zz are of opposite sign for q ∈ {0.2500, 0.3300, 0.3780} and are of the same sign for q ∈ {0.4400, 0.5176, 0.5774, 0.6180, 0.6600}, although Q zz is very close to zero for q = 0.4400. For α = 3, Q zw and Q zz are of opposite sign for all q apart from q = 0.6600. For 1 ≤ α ≤ 3, Q wz and Q ww are both negative. We connect some of the observed steady patterns in this section to the three cases of nonlinear wave-vector interactions described in Figure 3 and relate these to our expectations in Table 1. 7.1. Steady patterns with varying q: α = 1 First, we explore steady patterns with α = 1 at fixed r = 0.01 and θ = 45 • , so r 1 = r q = 0.00707, but for varying q (see Figure 7). For q < 1 2 , we see strong hexagonal motifs on a scale set by the smaller wavenumber q, inset with stripes on a scale of wavenumber 1, resembling patterns found by [10]. For q = 0.3780 the pattern is exactly hexagonal, with six equally spaced modes on the inner circle and twelve unequally spaced on the outer, as in Figure 3(a) -this is the simplest example of a superlattice pattern. As q increases beyond 0.5176, the (a) q = 0.2500 (e) q = 0.5176  patterns continue as essentially hexagonal on the scale of the smaller wavenumber, but defects and grain boundaries become more common for larger q.

Quasipatterns
We take q = 2 − √ 3 = 0.5176 and start with small amplitudes for twelve Fourier modes on the circle k = 1 as initial condition to encourage twelve-fold quasipatterns, finding stable examples as in Figure 8(a). This is a periodic approximant to a true quasipattern, but the approximation is particularly accurate in the 30 × 30 domain [46]. There are twelve peaks on the inner and outer circles, interleaved as in Figure 3(b). This kind of quasipattern has been seen in many similar kinds of calculations going back to [38,39].
We also obtain an eight-fold quasipattern, in Figure 8(b). This is surprising since neither θ z nor θ w is a multiple of 45 • (Table 3). In addition, in our 30 × 30 domain, the approximation to a true eight-fold quasipattern is not particularly accurate. Nonetheless, there are eight reasonably clear peaks on the inner circle, with sixteen diffuse peaks on the outer and an additional eight peaks just outside the outer circle, giving the impression of a regular octagon. It may be significant that θ z = 165.6 • (see Table 3), which is close to 15 • less than 180 • , as the twenty-four peaks on and just off the outer circle are spaced roughly 15 • apart.
The third common two-dimensional quasipattern has ten-fold symmetry. We have not found examples of such a quasipattern, but there are hints of a ten-fold motif in calculations with q = 0.6180 (see Figure 9b).

Steady complex patterns
We find many examples of hexagonal patterns with defects as in Figure 7(e-h). In Figure 9, we show two examples of steady complex patterns that are not just straightforward patches of hexagons (as in Figure 7e-h). In Figure 9(a), with q = 0.5176, the pattern has a "swirly" appearance with regions of distorted hexagons in between patches of more regular hexagons. The patches are rotated with respect to each other, leading to twelve broad peaks in the outer circle of the power spectrum. In Figure 9(b), with q = 0.6180, the complex structure of the pattern is more uniformly distributed, both in space and around the two circles in the power spectrum. There are several examples of a ten-fold motif, not surprising given that q is the inverse of the golden ratio.
Both examples are not steady but continue to evolve on timescales longer than 10000 time units. The example in Figure 9(a) eventually anneals to hexagons. The example in Figure 9(b) persists for at least 50000 time units, and is the closest we have found to an example of a steady complex pattern with the infinite set of wavevectors implied by Figure 3(c).

Spatiotemporal chaos
Finally, we show three examples of spatiotemporal chaos in Figure 10(a) and (b), with q = 0.4400 and q = 0.6180 respectively, and in Figure 11, with q = 0.3780 but with larger linear parameters than the others calculations (r = 0.03). The spatiotemporal chaos examples in Figure 10 evolve quite slowly: a frame spacing of 2000 time units is needed to show appreciable differences between the frames. The frame spacing in Figure 11 is 100 times less. Videos of all three examples are available in [40].
In the q < 1 2 example in Figure 10(a), there are evolving patches of elongated hexagons, and in some frames, the power spectrum has twelve peaks on the outer circle. In contrast, the q > 1 2 example in Figure 10(b) has a much more axisymmetric power spectrum and the complexity of the pattern is more uniformly spread across the domain. In the third example in Figure 11, the system alternates between episodes dominated by small hexagons and episodes dominated by larger structures.

Summary and Discussion
One main finding is that we only find persistent time dependence (as opposed to slow healing of defects and coarsening of grain boundaries) when Q zw and Q zz had opposite sign, as in Figures 10 and 11. This is consistent with our a priori expectations outlined in Table 1. With appropriate initial conditions, we find twelve-fold quasipatterns only in the case with q = 2 − √ 3, although we find eight-fold and hints of ten-fold quasipatterns in other cases. We did find examples of steady (or persistent) complex patterns, as in Figure 9. It was noticeable that having two interacting wavelengths encourages patterns with defects.  Figure 10: the time interval between frames is 20 time units. A video is available in [40].
With quadratic coefficients Q zw and Q zz of opposite sign, the preliminary 8 × 8 calculations often yield time-dependent patterns, with relatively simple oscillatory, chaotic or heteroclinic cycle dynamics. When we extend these into 30×30 domains, the simple time dependence is often replaced by spatiotemporal chaos, supporting the infinite set of wavevectors picture implied by Figure 3(c). We suspect that the reason for this is that in 8 × 8 domains, there are relatively few modes available close enough to each circle to participate in the dynamics. In contrast, with 30 × 30 domains, the density of modes in Fourier space is higher, and so modes are more likely to be able to participate in multiple threewave interactions, as in Figure 2. Considering a single set of modes coupled by a three-wave interaction, the modes may be oscillatory. When two (or more) sets of modes, also coupled within themselves by three-wave interaction, have modes in common, the common modes will be torn in different directions by their partners in the different sets, resulting in spatiotemporal chaos.
All interesting cases of time dependence have Q zw and Q zz of opposite sign, and time dependence can happen for all values of q. Having Q wz and Q ww of opposite sign (relevant only for q > 1 2 ) did not lead to persistent time dependence. Having q > 1 2 did appear to help when seeking steady complex patterns (Figure 9b): we find no examples of such patterns with q < 1 2 . This is in contrast to the hypotheses of [27], who argued that having q > 1 2 and having Q wz and Q ww of opposite sign should encourage complex patterns. In fact, we find that q < 1 2 is more interesting than anticipated from the results of [27], especially when Q zw and Q zz have opposite sign: there are many more states possible, including spatiotemporal chaos, going well beyond the steady superlattice example associated with q = 1/ √ 7 = 0.3780, in Figure 7(c). We also had not anticipated finding quasipatterns in the case q < 1 2 (as in Figure 8b), but recent work [29] suggests this warrants more exploration.
In summary, our main numerical findings described above are broadly in line with the a priori expectations in Table 1. In particular, time dependence, and with it complex spatial structure, requires Q zw and Q zz to have opposite sign. The other pair of quadratic coefficients (Q wz and Q ww ) can also have opposite sign for α = 5 and 6, but for these values, we mainly find domains of hexagons. At this point, it is not clear why this happens, nor whether other systems would behave differently. We should emphasise that the time-dependence we have found is not associated with a primary Hopf bifurcation to spiral waves, common in many Turing systems.
It would be interesting to explore these complex patterns in more detail, in the context of coupled reaction-diffusion systems, in the context of amplitude equations, especially in the case q < 1 2 , as initiated in [49], and in the context of simpler model PDEs, such as the ones proposed by [27,38,39,46] as models for Faraday waves. A related PDE is known to produce three-dimensional icosahedral quasipatterns [50] and localised quasipatterns [51], and such structures may also be possible in Turing systems. We plan to undertake further investigations in the future.
Of course, one has to ask whether these kinds of patterns can be found in experiments, and indeed if the mechanisms for forming them are as outlined in this paper. As explained in [13], manipulating the strength of the coupling, and indeed the diffusion constants, as we have done here, is difficult. Nonetheless, the spatially complex experimental patterns reported in [13] (see Figure 1) resemble, at least qualitatively, the images in Figures 10 and 11.