TURING PATTERNS IN PARABOLIC SYSTEMS OF CONSERVATION LAWS AND NUMERICALLY OBSERVED STABILITY OF PERIODIC WAVES BLAKE BARKER, SOYEUN JUNG, AND KEVIN ZUMBRUN Abstract. Turingpatternsonunboundeddomainshavebeenwidelystudiedinsystemsofreaction- diffusion equations. However, up to now, they have not been studied for systems of conservation 7 laws. Here, we (i) derive conditions for Turing instability in conservation laws and (ii) use these 1 conditionstofindfamiliesofperiodicsolutionsbifurcatingfromuniformstates,numericallycontin- 0 uing these families into the large-amplitude regime. For the examples studied, numerical stability 2 analysis suggests that stable periodic waves can emerge either from supercritical Turing bifurca- tionsor,viasecondarybifurcationasamplitudeisincreased, fromsub-criticalTuringbifurcations. n ThisanswersintheaffirmativeaquestionofOh-Zumbrunwhetherstableperiodicsolutionsofcon- a J servation laws can occur. Determination of a full small-amplitude stability diagram– specifically, 6 determinationofrigorousEckhaus-typestabilityconditions–remainsaninterestingopenproblem. 1 ] P A 1. Introduction . h The study of periodic solutions of conservation laws and their stability, initiated in [OZ03a, t OZ03b] and continued in [Ser05, JZ10], etc., has led to a number of interesting developments, a m particularly in the related study of roll-waves in inclined shallow-water flow. For an account of [ these developments, see, e.g., [JNRZ12] and references therein. However, in the original context of conservation laws, so far no example of a stable periodic wave has been found. Indeed, one of the 1 v primary results of [OZ03a, PSZ13] was that for the fundamental example of planar viscoelasticity, 9 stable periodic waves do not exist, due to a special variational structure of this particular system; 8 it was cited as a basic open problem whether stable periodic waves could arise for any system of 2 conservation laws, either physically motivated: or artificially contrived. 4 0 In the more standard context of reaction diffusion systems and classical pattern formation the- 1. ory, by contrast, stable periodic solutions are abundant and well-understood, through the mech- 0 anism of Turing instability, or bifurcation of small-amplitude, approximately-constant period, pe- 7 riodic solutions from a uniform state. For such waves, stability is completely determined by an 1 associated Eckhaus stability diagram, as derived formally in [Eck65] and verified rigorously in : v [Mie95, Mie97, Sch96, SZJV16], essentially by perturbation from constant-coefficient linearized i X behavior. By contrast, the small-amplitude waves investigated up to now (see Remark 3.2) come r through more complicated zero-wave number bifurcations in which period goes to infinity as am- a plitude goes to zero and the stability analysis is far from constant-coefficient (see, e.g., [Bar14] in the successfully-analyzed case of shallow-water flow). Our simple goal in this paper, therefore, is to seek stable periodic waves via a conservation law analog of Turing instability. In the first part, we find an analog of Turing instability, with which we are able to generate large numbers of examples of spatially periodic solutions of conservation laws. Next, we find an interesting dimensional restriction to systems of three or more coordinates, Date: January 17, 2017. Research of B.B. was partially supported under NSF grant no. DMS-1400872. Research of S.J. was partially supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIP) (No. 2016009978). Research of K.Z. was partially supported under NSF grants no. DMS-0300487 and DMS-0801745. 1 explaining the absence of Turing instabilities for 2×2 systems considered previously. Finally, we performanumericalexistence/stabilitystudyfor3×3examplesystemsexhibitingTuringinstability, answering in the affirmative the fundamental question posed in [OZ03a, PSZ13] whether there can exist stable spatially periodic solutions of systems of conservation laws, at least at the level of numerical approximation. These studies suggest that, for supercritical Turing bifurcation, stable waves can emerge through the small-amplitude limit and persist up to rather large amplitudes. For sub-critical Turing bifurcations, all emerging waves are necessarily initially unstable, but appear in some cases to undergo secondary bifurcation to stability as amplitude is further increased. The numerically observed stability of intermediate-amplitude waves we regard as conclusive. Delicacyofnumericalapproximationasamplitudegoestozero,however,preventsusfromobtaining a detailed stability diagram near the Turing bifurcation or even from making definitive conclusions about stability in that regime. Rigorous spectral stability analysis for conservation laws in this regime, analogous to those of [Mie95, Mie97, Sch96, SZJV16] in the reaction diffusion case, we regard therefore as a very interesting open problem. The studies in [MC00, Suk16] of reaction diffusionequationswithanassociatedconservationlawmayofferguidanceinsuchaninvestigation. 2. Turing instability for conservation laws We begin by defining a notion of Turing instability for systems of conservation laws (2.1) u +f(u;ε) = (Dεu ) , t x x x u ∈ Rn, where ε is a bifurcation parameter and D for simplicity is taken constant. Linearizing (2.1) about a uniform state u(x,t) ≡ u yields the family of constant-coefficient equations 0 (2.2) u = L(ε)u := −Aεu +Dεu t x xx with dispersion relations λ (ξ) ∈ σ(−iξAε −ξ2Dε), ξ ∈ R, where σ(·) here and elsewhere denotes j spectrum of a matrix or linear operator. The state u is spectrally (hence nonlinearly) stable if 0 (2.3) (cid:60)σ(−iξAε−ξ2Dε) ≤ −θ|ξ|2, θ > 0, for all ξ ∈ R [Kaw83]. Following the original philosophy applied by Turing [Tur52] to reaction diffusion systems, we seek a natural set of conditions guaranteeing low- and high-frequency stability– i.e., that (2.3) hold for |ξ| → 0,∞– but allowing instability at finite frequencies |ξ| (cid:54)= 0,∞. Should this be possible, then performing a homotopy in ε between stable and unstable states, we may expect generically to arrive at a special bifurcation point ε = ε , without loss of generality ε = 0, for which (2.3) holds ∗ ∗ uniformly away from special points ξ = ±ξ , at which ∗ (2.4) max(cid:60)σ(−iξAε∗ −ξ2Dε∗) = 0 ξ(cid:54)=0 is achieved (note, by complex conjugate symmetry, that extrema appear in ± pairs) and for which (2.3) fails strictly as ε is further increased. We may then conclude, by standard bifurcation theory applied to the domain of periodic functions with period X := 2π/ξ the appearance of nontrivial ∗ spatially periodic solutions with periods near X, similarly as in the reaction diffusion case [Mie95, Mie97, Sch96, SZJV16]. At ξ = 0, (2.3) yields that A is hyperbolic, in the sense that it has real semisimple eigenvalues. Without loss of generality, therefore, take A diagonal, with entries a , j = 1,...,n. In the sim- j plest case that A is strictly hyperbolic, in the sense that these a are distinct, we find by spectral j perturbation expansion about ξ = 0 [Kaw83] that the corresponding eigenvalue expansions are λ (ξ) = −ia ξ−D ξ2+O(ξ3), j j jj 2 so that (2.3) (ξ (cid:28) 1) is equivalent to the condition that D have positive diagonal entries D . jj Similarly, by spectral expansion about ξ = ∞, σ(−iξA−ξ2D) = −ξ2σ(D)+O(ξ), so that (2.3)(ξ = ∞) is equivalent to the condition that D be unstable, i.e., have eigenvalues with strictly positive real part. Collecting, our conditions are (C) : • Aε is diagonal with distinct entries, and • Dε has positive diagonal entries and eigenvalues with strictly positive real part. These are to be contrasted with Turing’s conditions in the reaction diffusion case u = Du +g(u) t xx that D be symmetric positive and A := dg(u) be symmetric negative definite [Tur52]. 2.1. Turing instability and Hopf bifurcation. Let (2.4) hold, with λ = ±iτ ∈ σ(−iξA−ξ2D) for ξ = ±ξ , ξ (cid:54)= 0. Then, changing to the moving coordinate frame x → x˜ := x − ct, for ∗ ∗ c := τ/ξ , or, equivalently, under the the change of coordinates A → A˜ := A − cI, we have ∗ λ = 0 ∈ σ(−iξA˜−ξ2D) for ξ = ±ξ , i.e., det(−iξA˜−ξ2D) = 0 at ξ = ξ , or ∗ ∗ (2.5) ±iξ ∈ σ(D−1A˜). ∗ Condition (2.5) may be recognized as the condition for Hopf bifurcation of an equilibrium u(x,t) ≡ constant of the traveling-wave ODE (2.6) Dεu(cid:48) = f(u;ε)−cu+q, whereq isaconstantofintegration,forwhichthelinearizedequationisu(cid:48) = D−1A˜u,A˜againdiago- nal. Thus, we recover by finite-dimensional bifurcation theory the previously-remarked appearance of nontrivial periodic solutions with period near X = 2π/ξ . We also obtain the alternative bifur- ∗ cation criterion (2.5). This simplifies the problem a great deal; for one thing, we are now working with real matrices, as occur for symbols in the reaction diffusion case, and not complex ones. 2.1.1. Dimensional count. From the usual Hopf bifurcation theorem for ODE, we find that for each fixed nearby q, c, there exists a one-parameter family of nontrivial periodic solutions bifurcating fromtheconstantsolution,genericallyparametrizednonsingularlybyperiodX. Thus,fixingq = 0, we obtain a 2-parameter family of periodic solutions, generically well-parametrized by c and X. 2.2. Finding Turing instabilities. To find Turing instability, we may seek Aε and Dε satisfying (C) , ε ∈ R a bifurcation parameter, such that (2.4) is violated at ε = 0 (instability), but (2.3) is satisfied for all ξ at ε = 1 (stability), for example if A1 = Id or D1 = Id. For, in this case, the conditions (C) on Aε, Dε insure that at the largest value ε of ε for which (2.4) is satisfied, the ∗ maximum (2.4) is achieved at some ξ = ξ (cid:54)= 0, while for ε > 0 there must be strictly positive real ∗ part eigenvalues, again bounded uniformly away from zero. As another approach, starting from the observation relating Turing instabilities and Hopf bifur- cation, notice first that (2.4) cannot occur for D = I, in which case the spectra of (−iξA−ξ2D) are simply λ (ξ) = −iξa −ξ2; nor can (2.5), since σ(A˜) is by assumption real. Thus, we suggest, j j first, finding examples Aˇ, Dˇ satisfying (2.5) either analytically or by checking random matrices, then, setting up a homotopy D(cid:15) := (cid:15)Dˇ +(1−(cid:15))I from the identity to Dˇ. Since, as just observed, σ(−iξAˇ− ξ2D(cid:15)) is stable for ε = 0, while for ε = 1 it is at most neutrally stable, having zero eigenvalues at ξ = ±ξ (cid:54)= 0, we find that for some (cid:15) ∈ (0,1], σ(−iξAˇ−ξ2D(cid:15)) is exactly neutral, ∗ i.e., a Turing instability, with eigenvalues ±iτ at ξ = ±ξˆ (note: different from the original ξ in ∗ ∗ general!). As described above, this corresponds to a Hopf bifurcation in the traveling-wave ODE for speed c := τ/ξˆ, with limiting wave number ξˆ and period X := 2π/ξˆ. ∗ ∗ ∗ ∗ ∗ 3 3. Negative results We next describe situations in which Turing instability cannot occur, narrowing our search. 3.1. The 2×2 case. We have the following result for n = 2, strikingly different from the situation of the reaction diffusion case. Proposition 3.1. Assuming (C) , there exist no Turing-type instabilities of (2.1) for n = 2. Proof. Take by assumption A diagonal. Since D−1A is real, appearance of a pure imaginary eigen- value iτ implies the appearance also of its complex conjugate −iτ, hence trace is zero and deter- (cid:18) (cid:19) α 0 minant is positive. By a scaling transformation S = not affecting diagonal form of A, we 0 β (cid:18) (cid:19) c 1 may arrange therefore that D−1A = =: J, for some c2 < 1. Noting that J2 = (c2−1)I, −1 −c (cid:18) (cid:19) a c a we may solve to obtain D = 1 AJ = 1 1 1 . The requirement that D have positive c2−1 c2−1 −a −a c 2 2 diagonal implies, with c2 < 1, that a c < 0 and a c > 0, so that a and a have opposite sign. But, 1 2 1 2 detD = (c2 −1)−2a a (1−c2) > 0 implies that a and a have the same sign, hence these two 1 2 1 2 conditions cannot hold at once. (cid:3) Example 3.2. The viscoelasticity model τ −u = d τ , u +p(τ) = d u studied by Oh- t x 11 xx t x 22 xx Zumbrun[OZ03a]fallsintotheaboveframework,hencedoesnotadmitTuringinstabilities. Infact, periodic waves arise in this model through Bogdanov-Takens bifurcation associated with splitting of two or more equilibria, a more complicated bifurcation far from constant-coefficient behavior. 3.2. Simultaneous symmetrizability. Another case in which Turing instabilities do not occur is when A and D are simultaneously symmetrizable, or, equivalently, can be converted by change of coordinates to be both symmetric. For, then, in the new coordinates, D, being symmetric positive definite, has a square root, and so D−1A is similar to the symmetric matrix D1/2D−1AD−1/2 = D−1/2AD−1/2, hence has real eigenvalues. More generally, it is easy to see that Turing instability does not occur for A symmetric and (cid:60)D := (1/2)(D+DT) > 0, since D−1Av = iτv would imply 0 = (cid:60)iτ(cid:104)v,Av(cid:105) = Re(cid:104)v,Dv(cid:105) = (cid:104)v,(cid:60)Dv(cid:105) > 0, a contradiction. This recovers the well-known fact that existence of a viscosity-compatible convex entropy for the system (2.1) implies nonexistence of non-constant stationary solutions, since existence of such an entropy implies the corresponding symmetryconditionsonthelinearizedequations. Thus,takingAwithoutlossofgeneralitydiagonal, we must specifically seek D nonsymmetric, D+DT nonpositive in order to find Turing instability. 3.3. Nonstrict hyperbolicity. Finally, we give a simple example showing that the condition of strict hyperbolicity of Aε is necessary in (C) . Consider the matrices     1 0 0 1 0 2 (3.1) Aε = 0 ε 0 and D = 0 1 1. 0 0 1 1 −2 1 Here, σ(D) = {1}; so −iξA−ξ2D is stable for |ξ| → +∞. For |ξ| → 0, we look at 2×2 blocks corresponding to the 1 and 3 entries of A and D, (cid:18) (cid:19) (cid:18) (cid:19) 1 0 1 2 (3.2) A˜= and D˜ = . 0 1 1 1 Then, thetwoeigenvaluesof−iξA−ξ2D closetoiξ forξ (cid:28) 1arebystandardspectralperturbation theory l (ξ) = −iξ − ξ2d˜ , where d˜ are eigenvalues of D˜. We easily see that D˜ has two real j j j eigenvalues with opposite sign because det(D˜) = −1 < 0. Thus, (2.3) is not satisfied for |ξ| → 0. 4 Remark 3.3. Though example (3.1), failing (C) , does not itself yield Turing instability, it is quite useful in finding nearby systems that do. For, note perturbation in ε generates matrices D−1A with nonstable eigenvalues despite A > 0. Perturbing first ε to obtain instability, then A still more slightly to recover strict hyperbolicity, we thus obtain an example satisfying (C) with unstable D−1A, which yields a Turing bifurcation upon homotopy D → I. We in fact used this method to generate the examples of Section 5. (We have generated other examples in other ways, that were not reported here; all exhibited similar behavior, however.) 4. Spectral and nonlinear stability Before describing our numerical investigations, we briefly recall the abstract stability framework developed in [OZ03a, JZ10, JNRZ12], etc., relevant to stability of the nontrivial periodic waves bifurcating from a constant solution at Turing instability. First, recall [JZ10, JNRZ12] that, under the condition of transversality of the associated periodic orbit of the traveling-wave ODE (guaran- teed in this case by the Hopf bifurcation scenario, for sufficiently small-amplitude waves), nonlinear stabilitywithrespecttolocalizedperturbationsoftheperiodicwaveconsideredasasolutiononthe whole line is determined (up to mild nondegeneracy conditions) by conditions of diffusive spectral stability, as we now describe. By Floquet theory, the L2(R) spectrum of the linearized operator L about a periodic wave of period X is entirely essential spectrum, corresponding to values λ ∈ C for which there exist generalized eigenfunction solutions v(x) = eiξxw(x), ξ ∈ R, of the associated eigenvalue equation (L−λ)v = 0 with w periodic, period X. The dissipative stability conditions are that this spectrum have real part ≤ −ηξ2, η > 0, for all ξ ∈ R, and strictly negative for (ξ,λ) (cid:54)= (0,0). For transversal orbits with ε bounded away from ε , the spectra near (ξ,λ) = (0,0) consists of ∗ the union of (n+1) smooth spectral curves λ (ξ) = −ia ξ+o(ξ) through the origin λ = 0, which, j j under the nondegeneracy condition that a be distinct, are analytic in ξ, admitting second-order j expansions (4.1) λ (ξ) = −ia ξ−b ξ2+O(ξ3), j = 1,...,n+1. j j j Moreover, the functions λ (ξ) correspond to the linearized dispersion relations for the associated j second-order Whitham system, an associated second-order (n+1)×(n+1) system of conservation laws formally governing slow modulational behavior [Whi11, Ser05, JNRZ12]. Thus, low-frequency diffusive spectral stability is equivalent to well-posedness (hyperbolic-parabolicity) of the Whitham system, which is in turn equivalent to reality of a (hyperbolicity) and positivity of (cid:60)b (parabol- j j icity) in (4.1), with high-frequency spectral stability given by (cid:60)λ ≤ −η < 0 for |ξ| ≥ η, η > 0. In the case of Turing instability, choosing the period X such that the wave-numbers ±ξ at ∗ ∗ ε = ε are equal to zero modulo 2π/X , we find by direct Fourier transform calculation that the ∗ ∗ constant solution at ε = ε has low-frequency spectrum consisting of (n+2) spectral curves passing ∗ through the origin, with all other spectra satisfying (cid:60)λ ≤ −η < 0 for some η > 0. The spectra of the bifurcating periodic waves perturbs smoothly from these values as ε is increased, hence high- frequency diffusive stability is guaranteed. However, low-frequency stability is now determined by a possibly complicated bifurcation of (n+2) spectral curves involving the (n+1) “Whitham curves” (4.1) passing through the origin plus an additional curve originating from the constant limit passing close to but not through the origin. These curves are clearly visible in the numerically approximated spectra displayed below in Section 5 for example systems with n = 3: namely, 4 Whitham curves passing through the origin, with a 5th (initially) neutral spectral curve passing near the origin, with all 5 of these passing through the origin at the bifurcation point ε = ε . ∗ 5 5. Numerical investigations Guided by the results of Sections 2, 3, and 4, we now perform the main work of the paper, carrying out numerical existence and stability investigations for periodic solutions of systems of conservation laws arising through Turing bifurcation from the uniform state in dimension n = 3. Numerics are carried out using the MATLAB-based package STABLAB developed for this purpose [BHLZ]. 5.1. Quadratic nonlinearity. We first consider the system (5.1) u +Aεu +N(u) = Du , t x x xx with 1 0 0 1 0 2 u2 1 (5.2) Aε := 0 a0 +ε 0, D := 0 1 1, and N(u) := β0, 22 0 0 3 1 −2 1 0 where a0 = 2.605173614560316. Here, ε is a bifurcation parameter that we will vary and u ≡ 0 is 22 a constant solution of (5.1). By linearization of (5.1) about u = 0, we have (5.3) u +Aεu = Du . t x xx We first check Turing-type instability conditions for u ≡ 0 in (5.3). Notice that Aε is strictly hyperbolic and D has positive diagonal entries with σ(D) = {1}, which means that −iξAε −ξ2D is stable near ξ = 0 or ξ = ±∞. We examine numerically stability of u ≡ 0 as ε changes. In Figure 1, we plot the spectrum of −iξAε −ξ2D with ε = −0.2, ε = 0, and ε = 0.2. It is seen that the constant solution u ≡ 0 is stable for ε < 0 and unstable for ε > 0. Thus, Turing instability occurs at ε = 0, that is, (2.4) is satisfied with ±iτ ∈ σ(−iξA0−ξ2D) for τ ≈ 1.5 and ξ ≈ ±1.16. As we ∗ observed in the previous section, ±iξ are eigenvalues of D−1(A0−c I) for c = τ ≈ 1.30. So the ∗ ∗ ∗ ξ∗ condition for Hopf bifurcation of a constant solution u ≡ 0 of the profile equation (5.4) −cu+Aεu+N(u) = Du(cid:48)+q is satisfied at the bifurcating point ε = 0 and c = c . Here q ∈ R3 is an integration constant and ∗ we fix q = 0 from now on. In Figure 2, we plot the spectrum of −iξ(Aε−c I)−ξ2D for the same ∗ ε as in Figure 1, showing how this moves the neutral spectrum from λ = ±iτ to λ = 0. (a) (b) (c) Figure 1. Plot with dots of a sampling of the spectrum of the constant solution, −iξA−ξ2D, with (a) ε = −0.2, (b) ε = 0, (c) ε = 0.2. The dashed vertical line marks the imaginary axis. The Hopf bifurcation leads to periodic profiles bifurcating from the uniform state u ≡ 0. In order to solve for these profiles, we let ε be a free variable and vary the period X and wave speed c, approximatingassociatedsolutionsusingtheperiodicprofilesolverbuiltintoSTABLAB,whichuses MATLAB’s Newton-based boundary-value problem solver bvp5c. In addition to periodic boundary 6 (a) (b) (c) Figure 2. Plot with dots of a sampling of the spectrum of the constant solution, −iξ(A−c I)−ξ2D, with (a) ε = −0.2, (b) ε = 0, (c) ε = 0.2 and c = c ≈ 1.30. ∗ ∗ The dashed vertical line marks the imaginary axis. conditions, the profile solver specifies a phase condition w·f(y(0)) = 0 where y(cid:48)(x) = f(y(x)) is the profile ODE ((2.6) in the present case) and w is a random vector. Unless w is a degenerate choice, w ·y˙(t) = 0 for some t by periodicity of y and Rolle’s Theorem, so this phase condition chooses a solution (at least locally) uniquely. To numerically solve the profile equation with a quadratic √ nonlinearity, we first obtain a solution by using as an initial guess u(x) = ε(cid:60)(e2πixv)/10, where v is the real part of an eigenvector, whose corresponding eigenvalue has non-zero imaginary part, of the profile Jacobian evaluated at the fixed point (0,0,0)T. That is, we start with an initial guess consisting of a strategically scaled periodic solution of the linearized equations at the bifurcation point ε = 0. Once we have a profile solution via this guess, we use continuation to solve for other profiles with nearby period X and speed c, obtaining thereby a full 2-parameter family of approximate solutions parametrized by (c,X), as described in Section 2.1.1. In Figure 3 (a) and (b), we plot the stability bifurcation diagram in the coordinates of shifted wavespeedc0 = c−c andperiodX. Thebifurcationdiagramshowsthatthereisafamilyofstable ∗ waves bifurcating from the Turing bifurcation. There is a small region of instability occurring from a “parabolic” Whitham instability, or change in curvature of a neutral spectral curve through the origin, correspondingtonegativediffusionorill-posednessoftheassociatedformalslowmodulation Whitham equations, which separates the region of stability near the Turing bifurcation point and the larger stability region. Figures 3 (d)-(f) demonstrate this onset of Whitham-type instability as seen in the spectrum of the bifurcating periodic waves. In Figure 3 (c), we see that the spectrum of the background constant solution becomes unstable as ε increases, so that the periodic profile shown in Figure 3 (g) comes into existence through a super-critical Hopf bifurcation. Finally, in Figure 3 (g), we plot the periodic profile for β = −10, ε = 2.82e−3, c = c +4.06e−3, X = 5.44. ∗ We note that, as described in Section 4, there are generically 4 neutral spectral curves passing throughtheorigin, withsecond-orderTaylorexpansionsrelatedtothelinearizeddispersionrelation for a formal Whitham slow-modulation approximation. This is clearly visible in Figure 3 (d)-(f). However, as seen in Figure 2 (b), the constant solution has 5 spectral curves passing through the origin at the bifurcation point and the spectra of bifurcating periodic waves perturbs from these 5 curves. So, at the bifurcation point, there is a 5th neutral curve passing through the origin, which remains nearby for values of ε nearby ε . It explains why the the spectrum of stable periodic waves ∗ bifurcatingfromTuringbifurcationinFigure3(d)hasanadditional5thcurvewhichisverycloseto the origin but not through the origin. Stability of small-amplitude waves is determined by behavior of these 5 neutral curves, either by movement of the maximum real part of the 5th curve into the unstable or stable half-plane (“co-periodic” stability, corresponding with super- or sub-criticality of the associated Hopf bifurcation), or by a “Whitham-type” instability consisting of loss of tangency 7 (a) (b) (c) (d) (e) (f) 4×10-3 y1 y2 2 y3 y(x)0 -2 -4 (g) 0 1 2 3 4 5 6 x Figure 3. (a)Stabilitybifurcationdiagraminthecoordinatesofshiftedwavespeed c0 = c − c and period X. Pink dots (light dots in grayscale) and black dots ∗ correspond respectively to stable and unstable waves. (b) Zoom in of (a) showing a family of stable waves in parameter space leading to the point of the Turing bifurcation. There is a small region of instability separating the stable waves near the Turing bifurcation point and the large stability region. (c) Plot of the spectrum of the zero constant solution when ε = 2.82e − 3, c = c + 4.06e − 3, and X = ∗ 5.44, indicating that the Turing bifurcation corresponds to a supercritical Hopf bifurcation. (d) Plot of the spectrum of a periodic wave in the family of stable wavesbifurcatingfromtheTuringbifurcation. (e)Plotofthespectrumofaperiodic wave in the family of unstable waves separating the two regions of stability. (f) Plot of the spectrum of a periodic wave in the large stability region. (g) Plot of the bifurcating periodic profile when ε = 2.82e−3, c = c +4.06e−3, and X = 5.44, ∗ with component one marked with a solid line, component two with a dashed line, and component three with a dot-dashed line. Throughout β = −10 and a dashed line marks the imaginary axis. to the imaginary axis (first-order, or “hyperbolic” instability) or change in curvature (2nd order, or “parabolic” instability) of one of the 4 neutral curves through the origin; see Section 4. For the quadratic nonlinearity, if u(x) is a profile solution for a fixed β, then −u(x) is a profile solution for −β, with the same value of ε. Thus, we are not able to produce a corresponding sub-critical Hopf bifurcation by reversing the sign of β, but a mirror super-critical bifurcation. To find examples of stable periodic profiles corresponding to both sub and super-critical Hopf bifurcations, we change the quadratic nonlinearity to a cubic nonlinearity in the next example, removing this symmetry and allowing us to change from super- to sub- by changing the sign of β. 8 5.2. Cubic nonlinearity. We consider next the system of conservation laws (5.5) u +Aεu +N(u) = Du , t x x xx with 1 0 0 1 0 2 u3 1 (5.6) Aε := 0 a0 +ε 0, D := 0 1 1, and N(u) := β0, 22 0 0 3 1 −2 1 0 where a0 = 2.605173614560316. Similarly as the quadratic example, we vary ε as a bifurcation 22 parameter. The stability of u ≡ 0 as ε varies is already shown in Figure 1 and Figure 2. Starting from the super-critical periodic profile solutions found previously for the quadratic nonlinearity, we obtain a solution for the cubic nonlinearity by continuation in a homotopy variable 0 ≤ h ≤ 1 via the nonlinearity N(U) = [β(hy3 +(1−h)y2),0,0]T. To obtain a sub-critical profile 1 1 solution for the cubic nonlinearity, we use the approximate symmetry (β,c,ε) → (−β,−c,−ε), which is valid at the linear periodic level only. Thereafter, we solve for profiles using continuation. In Figure 4, we plot the bifurcating stable periodic solution through a super-critical Hopf bifur- cation. Since ε > 0 for the constant solution to be unstable, as seen in Figure 2, the periodic profile showninFigure4(c)existsthroughasuper-criticalHopfbifurcation. Figure4(b)showsthestable spectrum of the periodic profile shown in (c). Here β = 10, c0 = 0.5, X = 6 , and ε = 8.74e−1. In Figure 4 (a), we plot a stability diagram in the coordinates of shifted wave speed c0 = c−c and ∗ period X. We do not find a family of stable waves bifurcating from the Turing instability. By changing the sign of β, we find the stable periodic solutions through a sub-critical Hopf bifurcation as demonstrated in Figure 5. Since ε < 0 for the constant solution to be stable, as seen in Figure 2, the periodic profile shown in Figure 5 (c) exists through a sub-critical Hopf bifurcation. Figure 5 (b) shows the stable spectrum of the periodic profile shown in (c). Here β = −10, c0 = −0.3, X = 4.5 , and ε = −3.5e−3. In Figure 5 (a), we plot a stability diagram in the coordinates of shifted wave speed c0 = c−c and period X. We do not find a family of stable ∗ waves bifurcating from the Turing instability. 1 0.4 y 1 0.3 y2 0.5 y 3 0.2 λm( ) 0 y(x)0.1 I 0 -0.5 -0.1 -1 -0.2 (a) (b) -0.4 -0.2Re( λ )0 0.2 (c) 0 2 4x 6 8 Figure 4. (a) Stability diagram in the coordinates of shifted wave speed c0 = c−c and period X for β = 10. Pink dots (light dots in grayscale) and black dots ∗ correspond respectively to stable and unstable waves. (b) For a stable wave, we plot in (b) its spectrum and in (c) the wave itself, with β = 10, c0 = 0.5, X = 6 , and ε = 8.74e−1. A dashed line marks the imaginary axis in (b). 5.3. Numerical stability method. To determine the spectrum of the periodic profiles, we used Hill’s method. The associated eigenvalue problem is given by Lv = λv where the linear oper- ator L takes the form L = (cid:80)mjkf (x) ∂q . The coefficients f (x) are X periodic. As j,k q=1 j,k,q ∂xq j,k,q in [DKCK07], we use a Fourier series to represent the coefficient functions f , f (x) = j,k,q j,k,q (cid:80)∞ φˆ ei2πjx/X, and write the generalized eigenfunctions as v(x) = eiξx(cid:80)∞ vˆ eiπjx/X, j=−∞ j,k,q j=−∞ j 9 1 0.5 λ ) m( 0 I -0.5 -1 -0.4 -0.2 0 0.2 (a) (b) Re( λ ) 0.3 y 1 0.2 y 2 y X = 5.4 3 0.1 y(x) 0 0.3 0.2 -0.1 ǫ 0.1 -0.2 0 -0.3 -0.2 -0.15 -0.1 -0.05 0 1 2 3 4 5 (c) x (d) c0 Figure 5. (a)Stabilitydiagraminthecoordinatesofshiftedwavespeedc0 = c−c ∗ and period X for β = −10. Pink dots (light dots in grayscale) and black dots correspond respectively to stable and unstable waves. For a stable wave, we plot in (b) its spectrum and in (c) the wave itself, with β = −10, c0 = −0.3, X = 4.5, and ε = −3.5e−3. A dashed line marks the imaginary axis in (b). In (d) we plot a curve showing existence, up to numerical approximation, of periodic profiles of period X = 5.4 in the parameters c0 and ε when β = −10 and the nonlinearity is cubic. A thin horizontal line marks the axis. where ξ ∈ (−π/2X,π/2X] is the Floquet exponent. Substituting these quantities into the eigen- value problem and equating coefficients gives an infinite dimensional eigenvalue problem for each fixedξ. BytruncatingtheFourierseriesatN termsandusingMatLabsFFTfunctiontodetermine the coefficients φˆ , we arrive at a finite dimensional eigenvalue problem Lξ vˆ = λvˆ, which we j,k,q N solvewithMATLAB’seigenvaluesolver. AllcomputationsweredoneusingSTABLAB[BHLZ]. For further information about Hill’s method and its convergence properties, see [CD10, DK06, JZ12]. 5.4. Computational statistics. All computations were carried out on a Macbook pro quad core or a Leopard WS desktop with 10 cores. Computing a profile took approximately 2 seconds or less, and computing the spectrum via Hill’s method took on average 20-60 seconds depending on the numberofmodesused. Wetypicallyused101Floquetparametersand41or81Fouriermodeswhen using Hill’s method. Each stability diagram took less then 24 hours to compute on the Leopard WS desktop. 6. Discussion and open problems WehaveidentifiedananalogofTuringinstabilityoccurringforn×nsystemsofconservationlaws of dimension n ≥ 3, leading to a large family of spatially periodic traveling waves. Our numerical stability investigations give convincing numerical evidence that at least some of these waves are 10