Spatio-temporal secondary instabilities near the Turing-Hopf bifurcation

In this work, we provide a framework to understand and quantify the spatiotemporal structures near the codimension-two Turing-Hopf point, resulting from secondary instabilities of Mixed Mode solutions of the Turing-Hopf amplitude equations. These instabilities are responsible for solutions such as (1) patterns which change their effective wavenumber while they oscillate as well as (2) phase instability combined with a spatial pattern. The quantification of these instabilities is based on the solution of the fourth order polynomial for the dispersion relation, which is solved using perturbation techniques. With the proposed methodology, we were able to identify and numerically corroborate that these two kinds of solutions are generalizations of the well known Eckhaus and Benjamin-Feir-Newell instabilities, respectively. Numerical simulations of the coupled system of real and complex Ginzburg-Landau equations are presented in space-time maps, showing quantitative and qualitative agreement with the predicted stability of the solutions. The relation with spatiotemporal intermittency and chaos is also illustrated.


Antecedents
In the vicinity of the CTHP, oscillating patterns, steady spatial patterns and homogeneous oscillations can be found 5 . They are the result of a primary instability, that is, the first transition occurring in a dynamical system as the result of a bifurcation 13 ; the Turing-Hopf bifurcation in this case. These three solutions can be described generally by assuming a Mixed Mode solution x t u( , ), which combines the influence of both a spatial stationary pattern with wavenumber k c and an homogeneous oscillation with angular frequency ω c : where H and T are the amplitudes of the Hopf and Turing modes, respectively, measuring the weight of each mode in the solution, and c.c. denotes the complex conjugate. k c and ω c can be obtained from the analysis of the linearized system, whereas the spatiotemporal evolution of the amplitudes can be obtained by the method of multiscales, as a perturbation series in powers of the distance of two parameters of the system to the CTHP. At first order of approximation, amplitudes T and H obey the coupled system of GL equations where a t , b t , c t are the real coefficients of the Turing mode equation, and a h = α r + iα i , b h = β r + iβ i and c h = γ r + iγ i are the complex coefficients of the equation for the wave mode 1 . μ t and μ h are proportional to the distance of two parameters to the CTHP. The approximate solution derived from the solution of (2) in (1), describes the spatiotemporal solutions of a great variety of dynamical systems near the Turing-Hopf bifurcation. Given this generality of the amplitude equations, the peculiarity of each dynamical system will be reflected in the coefficients of the approximated equations in (2), which depend on the parameters of the original system. The system (2) has three known families of solutions: (1) pure Turing structures, (2) plane waves, and (3) MM solutions, where both the spatial pattern and wave propagation occur 3 . The general form of the solution of (2) is = = .
− Ω iRx i Qx i t where Z, A and Ω depend on the wavenumbers R and Q 3 as: where τ α ≡ c / t r and σ γ ≡ a / r t . When Z ≠ 0 and A ≠ 0, (8) and (4) define the MM solution 3 . The pure Turing mode is recovered in the limit τ → 0 (α → ∞ r ) where www.nature.com/scientificreports www.nature.com/scientificreports/ These limiting cases (τ → 0 and σ → 0) will be important below to identify the weak coupling limit. The conditions for the existence of a MM solution require positive radicands in (4a) and (4b), which at the same time define the marginal conditions for R and Q. In particular, for the homogeneous state, Q = R = 0, from (4) we get that the existence conditions of the MM solution are τδ τσ t t h r 1 respectively. Since μ t and μ h measure the distance to their respective bifurcation, then δ μ μ ≡ / t h is a measure of the distance to both bifurcations since the cases |δ| ≈ 0 and δ | |  1 reflect more proximity of the system to the Turing or to the Hopf bifurcation, respectively. In what follows it will be assumed that existence conditions of the MM solution in (7), are fulfilled.
Equation (3) are exact solutions of (2) at some infinite set of the parameters, provided that the existence conditions stated before are fulfilled. However, near the CTHP, the occurrence of one or other solution among the three available families will depend upon their respective stability, which in turn it is expected to depend upon the proximity of the system to either one or another bifurcation. There are some few examples of RD systems where the stability of the solutions has been numerically studied 5,8,9,12 . However, up to now, the formal analysis is known only for the homogeneous state. As it is clear from substituting = = R Q 0 in Eq. (3), this corresponds to a spatial pattern of constant wavenumber oscillating in time with constant phase 6 . This is the expected solution on very small domains or in systems where no spatial perturbations occur 10 .
However, for most of the relevant systems, it is expected that inhomogeneous perturbations can cause also inhomogeneous states of the solution, i.e. states of the MM solution where R or Q are different from zero 12 . The effect of these inhomogeneous perturbations can be understood by substituting (3) in (1): where it is clear that R ≠ 0 accounts for a change in the effective wavenumber of the final spatial pattern respect to k c , and Q ≠ 0 produces a phase destabilization of the temporal solution and results in a wave behavior of the solution or even in phase chaos. We can consider these two types of solutions as secondary instabilities, that is, the result of a second transition to a new regime where the previous state (in our case, a oscillatory pattern of frequency ω c and wavenumber k c ), is in itself the result of a bifurcation (in our case, the Turing-Hopf bifurcation). As it was expected 5 and numerically corroborated 12 , the change of effective wavenumber and the phase unstabilization of the Mixed Mode, explained in (8), can be understood as generalizations of the Eckhaus and BFN instabilities previously described for the real and complex uncoupled GL equations separately 11 . Since our results on the analysis of (2) must recover these two well-known instabilities in the limit of zero coupling between Turing and Hopf modes, let us briefly consider the secondary instabilities of the real and complex GL equations separately.
If the system is near a Turing bifurcation but away from the Hopf bifurcation, the system obeys (2a) with c t = 0, i.e., the real GL equation. The solution of the system is a pattern with wavenumber k c + R, where the values of R for stability are determined by Eckhaus stability criterion: where λ Eck measures the temporal growth of sideband perturbations with wavenumber k 11 . In the supercritical 17 . Perturbations with larger wavenumbers can be destabilized if the domain size is large enough, since the periodic spatial solution admits only an integer number of wavelengths.
When the system is near a Hopf bifurcation, the system obeys (2b) with c h = 0, i.e., the complex GL equation. In this case, the solution of the system is a plane wave with phase φ = − Ω Qx t. The temporal growth λ BFN of sideband perturbations of wavenumber k is determined by where the first term of the right hand side is proportional to the group velocity of the wave, www.nature.com/scientificreports www.nature.com/scientificreports/ Recently, we have numerically identified these solutions in a RD system close to the CTHP, and shown that, in the Eckhaus instability, a change in the wavenumber of the spatial pattern occurs while the entire solution oscillates in phase and, in the case of the BFN phase instability, a spatial pattern with fixed positions of crests and troughs is combined with the propagation of a wavefront that change their height 12 . We also claimed that it is necessary to take into account inhomogeneous perturbations since the current analysis, which considers only homogeneous states, overestimates the regions of stability of pure Hopf and Turing solutions. This means that oscillatory spatial pattern can reach larger regions of the parameter space as the size of the domain is increased.
Since patterns arising close to the CTHP point have been frequently reported in experimental [18][19][20] and numerical [21][22][23] investigations and some of them could be the result of such secondary instabilities of the MM solution, in this work we provide a methodology to understand and quantify the stability of inhomogeneous states of the MM solution of (1). In the next section we describe a methodology to study these structures in a very general framework and to establish their stability conditions in terms of the system parameters.

Methodology
The purpose of this work is to quantify the secondary instabilities of the Mixed mode solution. Thus, the growth of inhomogeneous perturbations λ will be measured by following the standard procedure described in ref. 11 , where secondary instabilities of the real and complex GL equations are quantified separately. Our purpose is to obtain and solve the dispersion relation in order to identify how each secondary instability of the Turing and Hopf modes, given by (9) and (11) separately, is modified by the presence of the other mode. Since this analysis requires to find the solutions of a fourth order polynomial for λ(k), we will use a standard perturbative method for solving algebraic equations, where a parameter (in our case, the wavenumber of the sideband perturbation k) is small 15,16 . As we will see, with this approach a stability diagram where the kind of stability of the MM solutions with given R and Q in (3) can be theoretically predicted.
To find the secondary instabilities of the MM solution, we consider the sideband modulations of the harmonics 11 . The solutions (3) are perturbed as where p, q, r and s are complex numbers and the overline denotes the complex conjugate. By substituting (12) in (2), and neglecting non linear terms, we obtain the eigenvalue problem (M − λI)x = 0, where = p q r s x ( , , , ) T is the vector of perturbations, I is the 4 × 4 identity matrix and M is the matrix of coefficients with diagonal entries: 2 . The dispersion relation of the eigenvalue problem can be written as where f i (k) are polynomials as functions of the perturbation parameter k. In general, (13) has four possible roots, which represent the different behaviors of the solutions in (12). We will show that two roots are zero when k = 0 and, therefore, constitute bifurcations. It turns out that these two solutions correspond to generalizations of λ Eck and λ BFN given in (9) and (11), respectively. The stability criteria of the MM solution to homogeneous conditions will be obtained from the remaining two solutions. The roots of (13) are difficult to obtain and their exact values are complicated expressions from which no useful conclusions can be derived. In view of this, the idea of the perturbation method is to express the eigenvalue λ and the functions f i in (13) as a power series of the perturbation parameter k (assumed small) and then, solve algebraic equations recurrently. This will allow us to compute the solutions of λ up to second order, and compare with the Eckhaus and BFN criteria calculated in (9) and (11). Therefore, we expand λ and f i as www.nature.com/scientificreports www.nature.com/scientificreports/ By substituting in (13), the result can be arranged in powers of k as λ λ λ is equated to zero, the system can be solved recurrently. The values of the coefficients f ij are given in the Supplementary Material.
Since in our case = = f f 0 00 10 , the zero order solution is obtained from In this case, the two stable eigenvalues at order , which do not correspond to bifurcation eigenvalues, so they will be ignored. However, as the conditions in (16) are determined at zero order in k in (12), they determine the stability to homogeneous perturbations of the MM solution. Actually, (7) and (16) were deduced in ref. 10 as existence and stability conditions of the homogeneous MM solution.
For higher orders in k, we see that nothing new is obtained at first order, since λ λ P( , ) 1 0 20 . This means that and γ γ γ ≡ / i r . The labels Eck and BFN are included since these eigenvalues will be related to such instabilities, when higher orders in k are included. At order k, the traditional Eckhaus eigenvalue in (9) is 0 as in (17); in the same way, the group velocity of the wave solution in (10), is recovered from (17) in the limit σ → 0.
Since f 10 = 0, at order k 3 the polynomial P 3 depends only on λ 0 = 0, λ 1 and λ 2 . Therefore, it should be calculated for each value of λ 1 obtained in the previous paragraph.
For the Eckhaus related eigenvalue, λ λ λ in (5). In this case, we recover Eckhaus stability criterion in (9) for the real GL equation 17 . Notice however, that the Eckhaus criterion depends now on Q, the wavenumber of the wave solution, through Z. Our analysis predicts that this instability depends also on χ in (18), which is related to the complex coefficients of the Hopf mode, α, β and γ.
For the BFN related eigenvalue, λ λ λ λ = = P ( 0, , ) 3 2 / f 20 3 . The complete substitution is given at the Supplementary Material, but at first order in σ, it turns out to be To understand this eigenvalue, consider the limit of null interaction with the Turing mode (σ, → Z 0 and → A A Hopf ). In this case, we recover the BFN criterion (10) of the complex GL equation 2 . However, in our case, due to the proximity with the CTHP, both modes interact and the stability of the spatial pattern and wave solutions depend on each other.
From (20) it is straightforward to deduced that the origin is a critical point where λ This criterion was deduced in ref. 4 for the phase stability of the homogeneous MM solution.
In conclusion, we have proved that the generalization of the Eckhaus instability near the CTHP in (8) with R ≠ 0, represents a roll pattern solution, with an effective wavenumber moving from k c to k c + R when R and Q are inside the stability region given by the condition λ < 0 www.nature.com/scientificreports www.nature.com/scientificreports/ oscillates in time or propagates slowly as a wave 24 .
Eck t 2 , and therefore the destabilization of the system occurs when b t ≤ 0 as in the traditional Eckhaus criterion (9).
We also proved that the second approximated eigenvalue solution of (13) has a real and a pure imaginary part. The imaginary part, given by λ BFN 1 in (17), is related to the group velocity of the wave when Q ≠ 0 in (8). The unstable values of R and Q are defined by the condition λ > 0 BFN 2 . This generalization of the Benjamin-Feir-Newell instability near the CTHP establishes the phase destabilization of the system into wave undulations or phase chaos combined with a spatial pattern. Therefore, this route to chaos has a semi-organized background due to the interaction with the Turing mode [25][26][27][28] .
The solutions near the CTHP are both Eckhaus and BFN stable when λ Eck 2 , λ < 0 BFN 2 , respectively. These approximated criteria where calculated up to order k 2 in (14). However, since f 10 = 0, this process can be continued to higher orders in k. In particular, the fourth order approximation can be important to determine if the instabilities are of the long-wave type 2 . Now, for the regular Eckhaus instability in the real GL equation, the marginal and stability conditions are usually represented in the (μ t , R) plane as two parabolas. With this plot, the number of possible stable and unstable wavenumbers R, as function of the distance to the bifurcation μ t , can be measured. In our case, however, the visualization of the results is complicated by the fact that for the MM solution, there are two distances to the bifurcation and two wavenumbers indexing the solution, R and Q, and therefore, such a plot should be four-dimensional.
However, for a given set of parameters (i.e. for given μ h and μ t ), a level curve of the existence conditions in (4) and of the Eckhaus and BFN stability conditions given in (19) and (20) can be plotted as function of R and Q. An example of this stability diagram is given in Fig. 1a. Each point in this plot represents a MM solution given in (3), in terms of the wavenumbers R and Q. Given the even symmetry of the stability and existence conditions on R and Q, only the first quadrant is plotted. The outer curves R + and Q + are determined by the existence conditions (4), whereas the inner curves R − and Q − depend on the stability conditions in (19) and (20).
These conditions determine four regions, where the MM solution exists: 1) The region of stability, where MM solution with initial conditions in (3) will stay the same even under small random perturbations; 2) The region of Eckhaus unstable solutions, where the MM solution will keep the phase with constant slope Q, but the wavenumber of the Turing mode will decay to a band of stable wavenumbers; 3) The region of BFN unstable solutions, where MM solution in (3) will keep the Turing wavenumber but the phase will lose the wavenumber Q of the initial condition and, finally, 4) The region where both instabilities occur and, therefore, the initial pattern with wavenumbers R and Q in (3) will change when random perturbations are applied. In the next section, we will provide quantitative and qualitative evidence of the validity of (19) and (20) and consequently, of the stability diagrams thus constructed.

numerical Results
Up to now, our understanding of the spatiotemporal solutions resulting from the eigenvalues in (19) and (20) has been given through the analogy with the well-known Eckhaus and BFN instabilities. In this Section, numerical evidence of our theoretical findings is provided. Therefore, the Turing-Hopf amplitude equations in (2) (4), (16), (19) and (20). The regions in the stability diagram are labeled as S (stable), E (Eckhaus unstable), P (Phase instability) and B (both Eckhaus and BFN unstable). Parameters of (2)  www.nature.com/scientificreports www.nature.com/scientificreports/ of uniform rolls and the apparition of strong plateaus in T(x, t). Both behaviors are present when the solution is in the Eckhaus-BFN unstable region as illustrated in Fig. 2b.
In this way, we have shown the agreement of the theoretical and numerical results of the secondary instabilities in Fig. 2a, for a given set of parameters where the distances to the Turing-Hopf bifurcation remain fixed. The full understanding of the validity of the instabilities in (19) and (20) have to be extended to different set of parameters varying such distances. This will allow to measure the validity of the amplitude formalism as function of the proximity to the CTHP. This study requires a detailed analysis of the Fourier modes of the Turing and Hopf solutions for different stability diagrams, which requires future numerical investigations.
In what follows, we focus on showing that the existence and stability conditions given in (4), (19) and (20), produce different arrangements of the stability diagrams. The idea is to explain how the features of these plots can be used to understand diverse qualitative properties of the solutions of the GL equations in (2), in terms of the domain size. In particular, we are interested in showing that the domain sizes where instabilities occur are intrinsically related to the stability diagram, since the critical domain size for the occurrence of any instability is inversely proportional to the unstable values of R or Q.
All subsequent figures (Figs 4-6) show the spacetime solutions maps of the real and complex GL equations (2), solved for the amplitude of the Turing mode, and for the phase and modulus of the Hopf mode. As in Fig. 1, each row corresponds to a different domain size, namely, L = 25, 50, 75 and 100 and the same random initial conditions are used for illustrative purposes. In all these figures, the transient times are displayed in the vertical direction up to a total time t max = 100. In these plots, some spatio-temporal structures absent in the real and complex GL equations, separately, can be distinguished. The parameter values for each figure are given in the legend. The stability diagram of the solution with wavenumbers R and Q is also shown.
In Fig. 4 we contrast two examples where the critical size for the Eckhaus instability is very different. In the case of Fig. 4a, this instability cannot occur independently of the BFN instability and requires larger domain sizes. This is explained by the fact that the critical domain size for Eckhaus is inversely proportional to the unstable values of R. In contrast, the Eckhaus instability in Fig. 4b occurs for relatively small domains. Notice that in both cases, the coupling between the Turing and Hopf modes causes that the modulus |H(x, t)| exhibits interesting behaviors just where T(x, t) = 0, i.e. the vertical lines separating black and white zones in the space-map times of T. Near these zones, a kind of backbone pattern and a foliar structure appear in for |H| in Fig. 4a,b, respectively. All these structures result from the interaction of Turing and Hopf modes.
The qualitative correspondence between the stability diagram and the numerical solution is also exemplified in Fig. 5a,b, where the stability diagram is qualitatively similar in both cases. However, we have proved numerically that the Eckhaus instability occurs in Fig. 5b for shorter domains as predicted in our formalism. It is also confirmed in the numerical simulations that the solution with R = Q = 0 is unstable for any domain size in both cases. This statement was also proved for some parameter values of a reaction-diffusion system, where phase instability can occur even for very short domains (see ref. 12 ). In these Figures, patterns resemble those founded in refs. [29][30][31] for defect chaos and temporal intermittency.
Finally, in Fig. 6 we illustrate the phase chaos. This study suggest the appealing possibility of using the coupled real and complex GL equations in (2) for the study of intermittency and chaos already studied for the sole complex GL equation (see refs 30,31 ) for dynamical systems near the CTHP, where both behaviors appear together with a prevailing spatial pattern. www.nature.com/scientificreports www.nature.com/scientificreports/ Therefore, by means of numerical simulations, we have identified some qualitative relations between the size of the domain and the stability diagrams. We show that sometimes the homogeneous state is always unstable (Fig. 5), that in some systems the Eckhaus instability cannot occur without a phase unstabilization (Figs 4a and 5), and that the sizes of the domain where any instability occurs are directly related to the position of the unstable regions on the stability diagrams (by comparing Fig. 5a with Figs 5a and 1b with 6). Finally, with the numerical simulations we were able to detect that the coupling between Turing and Hopf modes gives place to structures not displayed by the real or complex GL equation separately, as the foliar and backbone patterns of Fig. 4.
Needless to say that these aspects predicted by the stability diagram should be corroborated numerically by the methodology used for Fig. 2a. However, we expect that the qualitative predictions provided by these examples will help to understand general features of the solutions in dynamical system near the CTHP.

Discussion
In this work we have provided analytic expressions for the Eckhaus and BFN instabilities in dynamical systems near the CTHP, which allowed to understand the change of wavenumber of the spatial pattern and the unstabilitization of the phase in the Mixed mode solution. A linear stability analysis of the perturbed oscillatory pattern solutions of the Turing-Hopf amplitude equations, yielded a fourth order polynomial dispersion relation (13) for the growth of the sideband perturbations. This polynomial was solved using standard perturbation techniques, thus the conditions for the stability of Mixed mode solutions with wavenumber R and Q in (3) were stablished. In order to verify the validity of our results, equations (19) and (20) were tested in three different ways. First we   www.nature.com/scientificreports www.nature.com/scientificreports/ have computed (19) and (20) in the limit of no coupling between Turing and Hopf modes, and the well known Eckhaus and BFN criteria of stability were recovered. Second, the predicted stability of the solutions in (3) was compared with the numerical simulations of (2), obtaining pretty good agreement. Finally, we show that the obtained stability diagrams are useful to understand some qualitative features of the numerical solutions as the size of the domain increases.
Our theoretical results are useful to understand diverse spatiotemporal solutions of a system near the CTHP, reported in numerical simulations and experiments. In particular, they corroborate that inhomogeneous perturbations give place to changes in the wavenumber of the spatial pattern and to undulations of the phase recently reported by us in a numerical study of a RD system 12 . Our present work allows us to relate them to generalizations of the Eckhaus and BFN instabilities of the MM solution, valid when the system is near the CTHP. Regarding consistency, we also show how all the known results concerning the existence and stability conditions of the MM solution to homogeneous perturbations 3,5,10 are recovered as limit cases.
More importantly, our analysis of the MM solution allows to construct the stability diagram of the solutions with wavenumbers R and Q using, besides the existence conditions implicit in (4), the stability results of the MM solution to inhomogeneous perturbations. Our methodology to quantify the stability of the solutions is a useful tool to identify the relationship between the domain size and both the phase instability and the change of the wavenumber of the spatial pattern, since these occur until a critical domain size, capable to fit wave modes that are destabilized by the Eckhaus and/or BFN instability, is reached. We validate our stability diagrams by means of numerical simulations of (2) for a set of parameters, giving a support to our theoretical findings, but further numerical work in this line is still pending.
In view of this, future investigations must be directed to establish the validity of (19) and (20) for a wider range of parameters, mainly away from the CTHP. Besides, although the scope of our results is limited to one-dimensional systems, we believe that the perturbative technique used in this work can be extended to study another type of systems of GL type.
In this work, we have shown the presence of the secondary instabilities from two kinds of experiments. We have shown that if the amplitude of the random perturbations is small compared with that of the tested modes, for example, the solution with wavenumbers (R, Q) = (0, 0), the result is a stable and homogeneous solution verified in Fig. 2a. This stability assumes that the amplitudes of T and H in the initial condition are much greater compared with that of the perturbations. However, for Fig. 1, the amplitudes of the mode with (R, Q) = (0, 0) is almost the same that for any other combination (R, Q) chosen randomly. In this case, what the system experiences is a competing of modes 32 which cannot be analyzed with either (12) or (19) and (20). However, as we have discussed in the previous section, these two formulas for predicting secondary instabilities provide a qualitative picture of what are the possible solutions of systems not prepared with any particular choice of initial conditions as the sizes of the domain increases. Our results show that by increasing the sizes of the domain the apparition of secondary instabilities is facilitated. The form of the solutions seems to be related to the intensity of the random perturbations and the sizes of the domain. How these two factors are related to the multiplication of Fourier modes in the solution and what is their relation with multiple modes of the perturbations remain as unsolved and interesting problems that should be addressed for the full understanding of the dynamical system.
It is very important to realize that our results on the stability of the solutions of the Turing-Hopf amplitude equations can be translated to the study of a particular dynamical system, for example a RD system. In this case, the difficulty lies in to determine how a single initial condition in the RD system gives place to two instabilities represented by the combination of two wavenumbers in (3). This is an unsolved problem. One would expect that the selection of the phase wavenumber, Q, and of the spatial modulation of the pattern, R, occurs by the proximity of these two wavenumbers to the bifurcation values, i.e., k = 0 for the limit cycle solution and k = k c for the Turing pattern. Therefore perturbations of a initial condition similar to (8) c.c., should be studied when k c is relatively large, as compared with R and Q, and therefore, each wavenumber can be related to the Turing and Hopf modes, respectively. However, this hypothesis for RD systems have to be tested in future works and constitutes a very interesting research line.
We believe that our work opens the way to future research linking the solutions of a RD system near the CTHP, the solutions of the amplitude equations in (2) and the stability predictions of the MM solution provided here. Finally it should be said that some of our numerical simulations of (2) encouraging to pursue in the study of intermittency and chaos over systems where a spatial pattern prevails, and to study some specific structures not arising in the separate real and complex GL equations.