Turing patterns in parabolic systems of conservation laws and numerically observed stability of periodic waves

TL;DR

This paper extends the study of Turing patterns to conservation laws, deriving instability conditions, identifying bifurcating periodic solutions, and numerically analyzing their stability, thus confirming the existence of stable periodic waves in such systems.

Contribution

It introduces the first conditions for Turing instability in conservation laws and demonstrates the emergence and stability of periodic solutions through numerical continuation.

Findings

Stable periodic waves can arise from supercritical Turing bifurcations.

Secondary bifurcations can produce stable waves from sub-critical bifurcations.

Numerical analysis suggests stable periodic solutions are possible in conservation laws.

Abstract

Turing patterns on unbounded domains have been widely studied in systems of reaction-diffusion equations. However, up to now, they have not been studied for systems of conservation laws. Here, we (i) derive conditions for Turing instability in conservation laws and (ii) use these conditions to find families of periodic solutions bifurcating from uniform states, numerically continuing these families into the large-amplitude regime. For the examples studied, numerical stability analysis suggests that stable periodic waves can emerge either from supercritical Turing bifurcations or, via secondary bifurcation as amplitude is increased, from sub-critical Turing bifurcations. This answers in the affirmative a question of Oh-Zumbrun whether stable periodic solutions of conservation laws can occur. Determination of a full small-amplitude stability diagram-- specifically, determination of rigorous Eckhaus-type stability conditions-- remains an interesting open problem.

Figures

Equations

$$u_t+f(u;\epsilon {)}_x=(D^{\epsilon }{u}_x{)}_x,$$

$$u_t=L(\epsilon )u:=-A^{\epsilon }{u}_x+D^{\epsilon }{u}_{xx}$$

$$\Re \sigma (-i\xi A^{\epsilon }-{\xi }^2{D}^{\epsilon })\le -\theta |\xi {|}^2,\theta >0,$$

$$\underset{\xi \ne 0}{\max }\Re \sigma (-i\xi A^{{\epsilon }_{\ast }}-{\xi }^2{D}^{{\epsilon }_{\ast }})=0$$

$${\lambda }_j(\xi )=-i{a}_j\xi -D_{jj}{\xi }^2+O({\xi }^3),$$

$$\sigma (-i\xi A-{\xi }^{2}D)=-{\xi }^2\sigma (D)+O(\xi ),$$

$$\pm i{\xi }_{\ast }\in \sigma (D^{-1}\tilde{A}).$$

$$D^{\epsilon }{u}^{\prime }=f(u;\epsilon )-cu+q,$$

$$A^{\epsilon }=\left(\begin{array}{ccc}1& 0& 0\\ 0& \epsilon & 0\\ 0& 0& 1\end{array}\right)\text{and}D=\left(\begin{array}{ccc}1& 0& 2\\ 0& 1& 1\\ 1& -2& 1\end{array}\right).$$

$$\tilde{A}=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\text{and}\tilde{D}=\left(\begin{array}{cc}1& 2\\ 1& 1\end{array}\right).$$

$${\lambda }_j(\xi )=-i{a}_j\xi -b_j{\xi }^2+O({\xi }^3),j=1,\dots ,n+1.$$

$$u_t+A^{\epsilon }{u}_x+N(u{)}_x=D{u}_{xx},$$

$$A^{\epsilon }:=\left(\begin{array}{ccc}1& 0& 0\\ 0& a_{22}^0+\epsilon & 0\\ 0& 0& 3\end{array}\right),D:=\left(\begin{array}{ccc}1& 0& 2\\ 0& 1& 1\\ 1& -2& 1\end{array}\right),\text{and}N(u):=\beta \left(\begin{array}{c}{u}_1^2\\ 0\\ 0\end{array}\right),$$

$$u_t+A^{\epsilon }{u}_x=D{u}_{xx}.$$

$$-cu+A^{\epsilon }u+N(u)=D{u}^{\prime }+q$$

$$u_t+A^{\epsilon }{u}_x+N(u{)}_x=D{u}_{xx},$$

$$A^{\epsilon }:=\left(\begin{array}{ccc}1& 0& 0\\ 0& a_{22}^0+\epsilon & 0\\ 0& 0& 3\end{array}\right),D:=\left(\begin{array}{ccc}1& 0& 2\\ 0& 1& 1\\ 1& -2& 1\end{array}\right),\text{and}N(u):=\beta \left(\begin{array}{c}{u}_1^3\\ 0\\ 0\end{array}\right),$$

Full text

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. Introduction

The study of periodic solutions of conservation laws and their stability, initiated in [OZ03a, OZ03b] and continued in [Ser05, JZ10], etc., has led to a number of interesting developments, 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 primary results of [OZ03a, PSZ13] was that for the fundamental example of planar viscoelasticity, stable periodic waves do not exist, due to a special variational structure of this particular system; it was cited as a basic open problem whether stable periodic waves could arise for any system of conservation laws, either physically motivated: or artificially contrived.

In the more standard context of reaction diffusion systems and classical pattern formation theory, by contrast, stable periodic solutions are abundant and well-understood, through the mechanism of Turing instability, or bifurcation of small-amplitude, approximately-constant period, periodic solutions from a uniform state. For such waves, stability is completely determined by an associated Eckhaus stability diagram, as derived formally in [Eck65] and verified rigorously in [Mie95, Mie97, Sch96, SZJV16], essentially by perturbation from constant-coefficient linearized behavior. By contrast, the small-amplitude waves investigated up to now (see Remark 3.2) come through more complicated zero-wave number bifurcations in which period goes to infinity as amplitude 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, explaining the absence of Turing instabilities for $2\times 2$ systems considered previously. Finally, we perform a numerical existence/stability study for $3\times 3$ example systems exhibiting Turing instability, 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. Delicacy of numerical approximation as amplitude goes to zero, however, prevents us from obtaining 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 diffusion equations with an associated conservation law may offer guidance in such an investigation.

2. Turing instability for conservation laws

We begin by defining a notion of Turing instability for systems of conservation laws

[TABLE]

$u\in{\mathbb{R}}^{n}$, where $\varepsilon$ is a bifurcation parameter and $D$ for simplicity is taken constant. Linearizing (2.1) about a uniform state $u(x,t)\equiv u_{0}$ yields the family of constant-coefficient equations

[TABLE]

with dispersion relations $\lambda_{j}(\xi)\in\sigma(-i\xi A^{\varepsilon}-\xi^{2}D^{\varepsilon})$, $\xi\in\mathbb{R}$, where $\sigma(\cdot)$ here and elsewhere denotes spectrum of a matrix or linear operator. The state $u_{0}$ is spectrally (hence nonlinearly) stable if

[TABLE]

for all $\xi\in\mathbb{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 $|\xi|\to 0,\infty$– but allowing instability at finite frequencies $|\xi|\neq 0,\infty$. Should this be possible, then performing a homotopy in $\varepsilon$ between stable and unstable states, we may expect generically to arrive at a special bifurcation point $\varepsilon=\varepsilon_{*}$, without loss of generality $\varepsilon_{*}=0$, for which (2.3) holds uniformly away from special points $\xi=\pm\xi_{*}$, at which

[TABLE]

is achieved (note, by complex conjugate symmetry, that extrema appear in $\pm$ pairs) and for which (2.3) fails strictly as $\varepsilon$ is further increased. We may then conclude, by standard bifurcation theory applied to the domain of periodic functions with period $X:=2\pi/\xi_{*}$ the appearance of nontrivial spatially periodic solutions with periods near $X$, similarly as in the reaction diffusion case [Mie95, Mie97, Sch96, SZJV16].

At $\xi=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}$, $j=1,\dots,n$. In the simplest case that $A$ is strictly hyperbolic, in the sense that these $a_{j}$ are distinct, we find by spectral perturbation expansion about $\xi=0$ [Kaw83] that the corresponding eigenvalue expansions are

[TABLE]

so that (2.3) ($\xi\ll 1$) is equivalent to the condition that $D$ have positive diagonal entries $D_{jj}$. Similarly, by spectral expansion about $\xi=\infty$,

[TABLE]

so that (2.3)($\xi=\infty$) is equivalent to the condition that $D$ be unstable, i.e., have eigenvalues with strictly positive real part. Collecting, our conditions are ($\mathcal{C}$) :

$\bullet$ $A^{\varepsilon}$ is diagonal with distinct entries, and

$\bullet$ $D^{\varepsilon}$ 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_{t}=Du_{xx}+g(u)$ 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 $\lambda=\pm i\tau\in\sigma(-i\xi A-\xi^{2}D)$ for $\xi=\pm\xi_{*}$, $\xi_{*}\neq 0$. Then, changing to the moving coordinate frame $x\to\tilde{x}:=x-ct$, for $c:=\tau/\xi_{*},$ or, equivalently, under the the change of coordinates $A\to\tilde{A}:=A-cI$, we have $\lambda=0\in\sigma(-i\xi\tilde{A}-\xi^{2}D)$ for $\xi=\pm\xi_{*}$, i.e., $\det(-i\xi\tilde{A}-\xi^{2}D)=0$ at $\xi=\xi_{*}$, or

[TABLE]

Condition (2.5) may be recognized as the condition for Hopf bifurcation of an equilibrium $u(x,t)\equiv\text{\rm constant}$ of the traveling-wave ODE

[TABLE]

where $q$ is a constant of integration, for which the linearized equation is $u^{\prime}=D^{-1}\tilde{A}u,$ $\tilde{A}$ again diagonal. Thus, we recover by finite-dimensional bifurcation theory the previously-remarked appearance of nontrivial periodic solutions with period near $X=2\pi/\xi_{*}$. We also obtain the alternative bifurcation 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 from the constant solution, generically parametrized nonsingularly by period $X$. Thus, fixing $q=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^{\varepsilon}$ and $D^{\varepsilon}$ satisfying ($\mathcal{C}$) , $\varepsilon\in\mathbb{R}$ a bifurcation parameter, such that (2.4) is violated at $\varepsilon=0$ (instability), but (2.3) is satisfied for all $\xi$ at $\varepsilon=1$ (stability), for example if $A^{1}={\rm Id}$ or $D^{1}={\rm Id}$. For, in this case, the conditions ($\mathcal{C}$) on $A^{\varepsilon}$, $D^{\varepsilon}$ insure that at the largest value $\varepsilon_{*}$ of $\varepsilon$ for which (2.4) is satisfied, the maximum (2.4) is achieved at some $\xi=\xi_{*}\neq 0$, while for $\varepsilon>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 bifurcation, notice first that (2.4) cannot occur for $D=I$, in which case the spectra of $(-i\xi A-\xi^{2}D)$ are simply $\lambda_{j}(\xi)=-i\xi a_{j}-\xi^{2}$; nor can (2.5), since $\sigma(\tilde{A})$ is by assumption real. Thus, we suggest, first, finding examples $\check{A}$, $\check{D}$ satisfying (2.5) either analytically or by checking random matrices, then, setting up a homotopy $D^{\epsilon}:=\epsilon\check{D}+(1-\epsilon)I$ from the identity to $\check{D}$. Since, as just observed, $\sigma(-i\xi\check{A}-\xi^{2}D^{\epsilon})$ is stable for $\varepsilon=0$, while for $\varepsilon=1$ it is at most neutrally stable, having zero eigenvalues at $\xi=\pm\xi_{*}\neq 0$, we find that for some $\epsilon\in(0,1]$, $\sigma(-i\xi\check{A}-\xi^{2}D^{\epsilon})$ is exactly neutral, i.e., a Turing instability, with eigenvalues $\pm i\tau$ at $\xi=\pm\hat{\xi}_{*}$ (note: different from the original $\xi_{*}$ in general!). As described above, this corresponds to a Hopf bifurcation in the traveling-wave ODE for speed $c_{*}:=\tau/\hat{\xi}_{*}$, with limiting wave number $\hat{\xi}_{*}$ and period $X_{*}:=2\pi/\hat{\xi}_{*}$.

3. Negative results

We next describe situations in which Turing instability cannot occur, narrowing our search.

3.1. The $2\times 2$ case

We have the following result for $n=2$, strikingly different from the situation of the reaction diffusion case.

Proposition 3.1**.**

Assuming ($\mathcal{C}$) , there exist no Turing-type instabilities of (2.1) for $n=2$.

Proof.

Take by assumption $A$ diagonal. Since $D^{-1}A$ is real, appearance of a pure imaginary eigenvalue $i\tau$ implies the appearance also of its complex conjugate $-i\tau$, hence trace is zero and determinant is positive. By a scaling transformation $S=\begin{pmatrix}\alpha&0\\ 0&\beta\end{pmatrix}$ not affecting diagonal form of $A$, we may arrange therefore that $D^{-1}A=\begin{pmatrix}c&1\\ -1&-c\end{pmatrix}=:J,$ for some $c^{2}<1$. Noting that $J^{2}=(c^{2}-1)I,$ we may solve to obtain $D=\frac{1}{c^{2}-1}AJ=\frac{1}{c^{2}-1}\begin{pmatrix}a_{1}c&a_{1}\\ -a_{2}&-a_{2}c\end{pmatrix}.$ The requirement that $D$ have positive diagonal implies, with $c^{2}<1$, that $a_{1}c<0$ and $a_{2}c>0$, so that $a_{1}$ and $a_{2}$ have opposite sign. But, $\det D=(c^{2}-1)^{-2}a_{1}a_{2}(1-c^{2})>0$ implies that $a_{1}$ and $a_{2}$ have the same sign, hence these two conditions cannot hold at once. ∎

Example 3.2**.**

The viscoelasticity model $\tau_{t}-u_{x}=d_{11}\tau_{xx}$, $u_{t}+p(\tau)_{x}=d_{22}u_{xx}$ studied by Oh-Zumbrun [OZ03a] falls into the above framework, hence does not admit Turing instabilities. In fact, 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^{-1}A$ is similar to the symmetric matrix $D^{1/2}D^{-1}AD^{-1/2}=D^{-1/2}AD^{-1/2}$, hence has real eigenvalues. More generally, it is easy to see that Turing instability does not occur for $A$ symmetric and $\Re D:=(1/2)(D+D^{T})>0$, since $D^{-1}Av=i\tau v$ would imply $0=\Re i\tau\langle v,Av\rangle=Re\langle v,Dv\rangle=\langle v,\Re Dv\rangle>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 symmetry conditions on the linearized equations. Thus, taking $A$ without loss of generality diagonal, we must specifically seek $D$ nonsymmetric, $D+D^{T}$ 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^{\varepsilon}$ is necessary in ($\mathcal{C}$) . Consider the matrices

[TABLE]

Here, $\sigma(D)=\{1\}$; so $-i\xi A-\xi^{2}D$ is stable for $|\xi|\rightarrow+\infty$. For $|\xi|\rightarrow 0$, we look at $2\times 2$ blocks corresponding to the 1 and 3 entries of $A$ and $D$,

[TABLE]

Then, the two eigenvalues of $-i\xi A-\xi^{2}D$ close to $i\xi$ for $\xi\ll 1$ are by standard spectral perturbation theory $\lambda_{j}(\xi)=-i\xi-\xi^{2}\tilde{d}_{j}$ , where $\tilde{d}_{j}$ are eigenvalues of $\tilde{D}$. We easily see that $\tilde{D}$ has two real eigenvalues with opposite sign because $det(\tilde{D})=-1<0$. Thus, (2.3) is not satisfied for $|\xi|\rightarrow 0$.

Remark 3.3*.*

Though example (3.1), failing ($\mathcal{C}$) , does not itself yield Turing instability, it is quite useful in finding nearby systems that do. For, note perturbation in $\varepsilon$ generates matrices $D^{-1}A$ with nonstable eigenvalues despite $A>0$. Perturbing first $\varepsilon$ to obtain instability, then $A$ still more slightly to recover strict hyperbolicity, we thus obtain an example satisfying ($\mathcal{C}$) with unstable $D^{-1}A$, which yields a Turing bifurcation upon homotopy $D\to 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 (guaranteed in this case by the Hopf bifurcation scenario, for sufficiently small-amplitude waves), nonlinear stability with respect to localized perturbations of the periodic wave considered as a solution on the whole line is determined (up to mild nondegeneracy conditions) by conditions of diffusive spectral stability, as we now describe.

By Floquet theory, the $L^{2}(\mathbb{R})$ spectrum of the linearized operator $L$ about a periodic wave of period $X$ is entirely essential spectrum, corresponding to values $\lambda\in\mathbb{C}$ for which there exist generalized eigenfunction solutions $v(x)=e^{i\xi x}w(x)$, $\xi\in\mathbb{R}$, of the associated eigenvalue equation $(L-\lambda)v=0$ with $w$ periodic, period $X$. The dissipative stability conditions are that this spectrum have real part $\leq-\eta\xi^{2}$, $\eta>0$, for all $\xi\in\mathbb{R}$, and strictly negative for $(\xi,\lambda)\neq(0,0)$.

For transversal orbits with $\varepsilon$ bounded away from $\varepsilon_{*}$, the spectra near $(\xi,\lambda)=(0,0)$ consists of the union of $(n+1)$ smooth spectral curves $\lambda_{j}(\xi)=-ia_{j}\xi+o(\xi)$ through the origin $\lambda=0$, which, under the nondegeneracy condition that $a_{j}$ be distinct, are analytic in $\xi$, admitting second-order expansions

[TABLE]

Moreover, the functions $\lambda_{j}(\xi)$ correspond to the linearized dispersion relations for the associated second-order Whitham system, an associated second-order $(n+1)\times(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_{j}$ (hyperbolicity) and positivity of $\Re b_{j}$ (parabolicity) in (4.1), with high-frequency spectral stability given by $\Re\lambda\leq-\eta<0$ for $|\xi|\geq\eta$, $\eta>0$.

In the case of Turing instability, choosing the period $X_{*}$ such that the wave-numbers $\pm\xi_{*}$ at $\varepsilon=\varepsilon_{*}$ are equal to zero modulo $2\pi/X_{*}$, we find by direct Fourier transform calculation that the constant solution at $\varepsilon=\varepsilon_{*}$ has low-frequency spectrum consisting of $(n+2)$* spectral curves passing through the origin,* with all other spectra satisfying $\Re\lambda\leq-\eta<0$ for some $\eta>0$. The spectra of the bifurcating periodic waves perturbs smoothly from these values as $\varepsilon$ 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 $5$th (initially) neutral spectral curve passing near the origin, with all $5$ of these passing through the origin at the bifurcation point $\varepsilon=\varepsilon_{*}$.

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

[TABLE]

with

[TABLE]

where $a_{22}^{0}=2.605173614560316$. Here, $\varepsilon$ is a bifurcation parameter that we will vary and $u\equiv 0$ is a constant solution of (5.1). By linearization of (5.1) about $u=0$, we have

[TABLE]

We first check Turing-type instability conditions for $u\equiv 0$ in (5.3). Notice that $A^{\varepsilon}$ is strictly hyperbolic and $D$ has positive diagonal entries with $\sigma(D)=\{1\}$, which means that $-i\xi A^{\varepsilon}-\xi^{2}D$ is stable near $\xi=0$ or $\xi=\pm\infty$. We examine numerically stability of $u\equiv 0$ as $\varepsilon$ changes. In Figure 1, we plot the spectrum of $-i\xi A^{\varepsilon}-\xi^{2}D$ with $\varepsilon=-0.2$, $\varepsilon=0$, and $\varepsilon=0.2$. It is seen that the constant solution $u\equiv 0$ is stable for $\varepsilon<0$ and unstable for $\varepsilon>0$. Thus, Turing instability occurs at $\varepsilon=0$, that is, (2.4) is satisfied with $\pm i\tau\in\sigma(-i\xi A^{0}-\xi^{2}D)$ for $\tau\approx 1.5$ and $\xi_{*}\approx\pm 1.16$. As we observed in the previous section, $\pm i\xi_{*}$ are eigenvalues of $D^{-1}(A^{0}-c_{*}I)$ for $c_{*}=\frac{\tau}{\xi_{*}}\approx 1.30$. So the condition for Hopf bifurcation of a constant solution $u\equiv 0$ of the profile equation

[TABLE]

is satisfied at the bifurcating point $\varepsilon=0$ and $c=c_{*}$. Here $q\in{\mathbb{R}}^{3}$ is an integration constant and we fix $q=0$ from now on. In Figure 2, we plot the spectrum of $-i\xi(A^{\varepsilon}-c_{*}I)-\xi^{2}D$ for the same $\varepsilon$ as in Figure 1, showing how this moves the neutral spectrum from $\lambda=\pm i\tau$ to $\lambda=0$.

The Hopf bifurcation leads to periodic profiles bifurcating from the uniform state $u\equiv 0$. In order to solve for these profiles, we let $\varepsilon$ be a free variable and vary the period $X$ and wave speed $c$, approximating associated solutions using the periodic profile solver built into STABLAB, which uses MATLAB’s Newton-based boundary-value problem solver bvp5c. In addition to periodic boundary conditions, the profile solver specifies a phase condition $w\cdot f(y(0))=0$ where $y^{\prime}(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\cdot\dot{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)=\sqrt{\varepsilon}\Re(e^{2\pi ix}v)/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 $\varepsilon=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 wave speed $c^{0}=c-c_{*}$ and period $X$. The bifurcation diagram shows that there is a family of stable 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, corresponding to negative diffusion or ill-posedness of the associated formal slow modulation 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 $\varepsilon$ 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 $\beta=-10$, $\varepsilon=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 through the origin, with second-order Taylor expansions related to the linearized dispersion relation 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 $\varepsilon$ nearby $\varepsilon_{*}$. It explains why the the spectrum of stable periodic waves bifurcating from Turing bifurcation in Figure 3 (d) has an additional 5th curve which is very close to 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 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 $\beta$, then $-u(x)$ is a profile solution for $-\beta$, with the same value of $\varepsilon$. Thus, we are not able to produce a corresponding sub-critical Hopf bifurcation by reversing the sign of $\beta$, 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 $\beta$.

5.2. Cubic nonlinearity

We consider next the system of conservation laws

[TABLE]

with

[TABLE]

where $a_{22}^{0}=2.605173614560316$. Similarly as the quadratic example, we vary $\varepsilon$ as a bifurcation parameter. The stability of $u\equiv 0$ as $\varepsilon$ 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\leq h\leq 1$ via the nonlinearity $N(U)=[\beta(hy_{1}^{3}+(1-h)y_{1}^{2}),0,0]^{T}$. To obtain a sub-critical profile solution for the cubic nonlinearity, we use the approximate symmetry $(\beta,c,\varepsilon)\to(-\beta,-c,-\varepsilon)$, 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 bifurcation. Since $\varepsilon>0$ for the constant solution to be unstable, as seen in Figure 2, the periodic profile shown in Figure 4 (c) exists through a super-critical Hopf bifurcation. Figure 4 (b) shows the stable spectrum of the periodic profile shown in (c). Here $\beta=10$, $c^{0}=0.5$, $X=6$ , and $\varepsilon=8.74e-1$. In Figure 4 (a), we plot a stability diagram in the coordinates of shifted wave speed $c^{0}=c-c_{*}$ and period $X$. We do not find a family of stable waves bifurcating from the Turing instability.

By changing the sign of $\beta$, we find the stable periodic solutions through a sub-critical Hopf bifurcation as demonstrated in Figure 5. Since $\varepsilon<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 $\beta=-10$, $c^{0}=-0.3$, $X=4.5$ , and $\varepsilon=-3.5e-3$. In Figure 5 (a), we plot a stability diagram in the coordinates of shifted wave speed $c^{0}=c-c_{*}$ and period $X$. We do not find a family of stable waves bifurcating from the Turing instability.

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=\lambda v$ where the linear operator $L$ takes the form $L_{j,k}=\sum_{q=1}^{m_{jk}}f_{j,k,q}(x)\frac{\partial^{q}}{\partial x^{q}}$. The coefficients $f_{j,k,q}(x)$ are $X$ periodic. As in [DKCK07], we use a Fourier series to represent the coefficient functions $f_{j,k,q}$, $f_{j,k,q}(x)=\sum_{j=-\infty}^{\infty}\hat{\phi}_{j,k,q}e^{i2\pi jx/X}$, and write the generalized eigenfunctions as $v(x)=e^{i\xi x}\sum_{j=-\infty}^{\infty}\hat{v}_{j}e^{i\pi jx/X}$, where $\xi\in(-\pi/2X,\pi/2X]$ is the Floquet exponent. Substituting these quantities into the eigenvalue problem and equating coefficients gives an infinite dimensional eigenvalue problem for each fixed $\xi$. By truncating the Fourier series at $N$ terms and using MatLabs FFT function to determine the coefficients $\hat{\phi}_{j,k,q}$, we arrive at a finite dimensional eigenvalue problem $L_{N}^{\xi}\hat{v}=\lambda\hat{v}$, which we solve with MATLAB’s eigenvalue solver. All computations were done using STABLAB [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 number of modes used. We typically used 101 Floquet parameters and 41 or 81 Fourier modes when using Hill’s method. Each stability diagram took less then 24 hours to compute on the Leopard WS desktop.

6. Discussion and open problems

We have identified an analog of Turing instability occurring for $n\times n$ systems of conservation laws of dimension $n\geq 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 stable, answering the question posed in [OZ03a, PSZ13] whether there can exist stable periodic solutions of conservation laws.

Moreover, the same numerical investigations indicate that at least for some model parameters, the bifurcation diagram near Turing instability/Hopf bifurcation includes an open region of instability. This opens the possibility for rigorous proof of existence of stable periodic waves through a small-amplitude bifurcation analysis as carried out in [Mie95, Mie97, Sch96, SZJV16] for the reaction diffusion case. Such an analysis we consider an extremely interesting open problem. Note, however, that it is inherently more complicated than the reaction diffusion version, involving $n+2$ bifurcation parameters $(X,c,q)$, $X,c\in\mathbb{R}^{1}$, $q\in\mathbb{R}^{n}$ rather than the two parameters of the reaction diffusion case. For an example of intermediate complexity, we point to the recent analyses [MC00, Suk16] of reaction diffusion equations with a single conserved quantity, featuring a three-parameter bifurcation.