Pattern formation in the damped Nikolaevskiy equation

The Nikolaevskiy equation has been proposed as a model for seismic waves, electroconvection and weak turbulence; we show that it can also be used to model transverse instabilities of fronts. This equation possesses a large-scale"Goldstone"mode that significantly influences the stability of spatially periodic steady solutions; indeed, all such solutions are unstable at onset, and the equation exhibits so-called soft-mode turbulence. In many applications, a weak damping of this neutral mode will be present, and we study the influence of this damping on solutions to the Nikolaevskiy equation. We examine the transition to the usual Eckhaus instability as the damping of the large-scale mode is increased, through numerical calculation and weakly nonlinear analysis. The latter is accomplished using asymptotically consistent systems of coupled amplitude equations. We find that there is a critical value of the damping below which (for a given value of the supercriticality parameter) all periodic steady states are unstable. The last solutions to lose stability lie in a cusp close to the left-hand side of the marginal stability curve.


I. INTRODUCTION
The Nikolaevskiy equation has been widely studied, because of its application to several physical systems and its interesting nonlinear dynamics. It can also be written in the alternative form after introducing u = φ x . The equation (1) was first derived, in an extended form, by Nikolaevskiy as a model for seismic waves in the earth's crust [1,2]. More recently, other applications of (1) or (2) have been proposed. Fujisaka and Yamada [3] and Tanaka [4] have derived (2) as a possible phase equation arising in reaction-diffusion systems. Turbulence in electroconvection was studied experimentally by Hidaka et al. [5], who proposed that the chaotic dynamics was of a type described qualitatively by (1). More generally, the Nikolaevskiy equation can be regarded as a canonical model for a patternforming system with the additional feature of symmetry under the transformation φ → φ + constant, in the form (2) [6], or with Galilean symmetry, in the form (1) [7]. Therefore, (1) is a suitable model for pattern-forming systems with these additional symmetries. Finally, as we show in Section II below, (2) can be used to describe finite-wavelength instabilities of travelling fronts.
The Nikolaevskiy equation exemplifies what is known as "soft-mode turbulence". This form of chaotic dynamics arises from the interaction between a pattern of finite wavenumber, which appears for r > 0, and a long-wave neutral (or "Goldstone") mode. The neutral mode arises as a direct consequence of the additional symmetries discussed above. In a sufficiently large domain, all spatially periodic steady states are unstable, as was shown by Tribelsky and Velarde [8], and numerical investigations of (1) show irregular chaotic patterns [6,7,9]. In Ref. [7] we showed that in this turbulent regime, (1) exhibits an unusual scaling, with the amplitude of u proportional to the 3/4 power of the supercriticality parameter r, whereas it is more usually proportional to the 1/2 power in other pattern-forming systems. Fujisaka et al. [10] extended (2) to two spatial dimensions and found that this same scaling law holds.
In this paper, we consider the effect of adding a weak damping term to the neutral mode in (1). The motivation for this is that in real experimental systems, the symmetry that gives rise to the neutral mode will often be weakly broken. For example, true Galilean symmetry is broken by the existence of distant boundaries. In the case of electroconvection, a weak magnetic field damps the neutral mode [11]. A second motivation for including the damping term is to improve our understanding of the appearance of chaotic dynamics in the undamped equation (1). If the neutral mode is strongly damped, then the dynamics of patterns is governed by the Ginzburg-Landau equation and stable, steady patterns are observed. As the damping is reduced, a transition to Nikolaevskiy chaos can be observed and investigated. Our work extends the recent study by Tribelsky [12].
In Sections II and III, we motivate and discuss the undamped and damped versions of the Nikolaevskiy equation. Then in Section IV we compute numerically the steady spatially periodic "roll" solutions of the damped equation and compute their secondary stability. This numerical calculation shows that there is an additional pocket of stable rolls, beyond those considered by Tribelsky [12]; indeed this pocket is significant in containing the last rolls to remain stable as the damping coefficient is reduced. Most of the stability boundary of the rolls corresponds to large-scale modulational instabilities; we then analyse these instabilities in three regimes, depending on the relative sizes of the supercriticality and damping parameters. In the strong-damping regime, the relevant weakly nonlinear description is a modified Ginzburg-Landau equation. In the two other regimes, where the damping is either moderate or weak, our analyses are based on three coupled amplitude equations, for the amplitude and phase of the rolls and for an associated large-scale field; in each case, these amplitude equations may be reduced to a single, third-order evolution equation for the phase. In all three regimes, we examine the secondary stability of rolls, and in certain cases we are also able to determine the effects of the leading nonlinear term and comment on the direction of the bifurcation. We demonstrate the subcritical onset of instability in some cases, including the no-damping case, consistent with the observed sudden onset of the instability in numerical simulations of (1).

II. THE NIKOLAEVSKIY EQUATION AS A MODEL FOR TRANSVERSE INSTABILITY OF FRONTS
In this section we explain how the Nikolaevskiy equation may arise as a model for finite-wavenumber transverse instabilities of planar fronts. Consider a planar travelling front in a homogeneous, isotropic medium, arising from, for example, a combustion problem or a system of reaction-diffusion equations. Suppose that the front is travelling at speed c in the z direction and that the x coordinate is directed along the front. Consider now small perturbations to the uniform front. Let the position of the front under perturbation be denoted by z = ct + φ(x, t). Now the evolution equation for φ must have the property that it does depend on the value of φ, because of translational invariance; the equation can depend only on x-derivatives of φ. Furthermore, reflection symmetry in x means that the only permissible linear terms involve even derivatives of φ. Therefore, the linearized equation for φ must have the form where the a j are constants. The corresponding dispersion relation is for the growth of modes φ ∼ exp(λt+ıkx). From (4) it is apparent that there are two possible types of instability of the front. If a 2 < 0 and a 4 < 0, then the planar front is unstable to a band of wavenumbers near k = 0 ( Figure 1, dashed line). This case has been widely studied, and leads, in the weakly nonlinear regime, to the Kuramoto-Sivashinsky equation [13,14], which, with suitable rescalings of x and t, may be written in the form But if a 2 , a 4 and a 6 are all positive, then long-wave modes are stable and a finite-wavenumber instability is possible, as illustrated by the solid line in Figure 1. After a rescaling of x and t, (3) then corresponds to the linear terms in (2). These two possible types of instability are sometimes referred to as "type II" and "type I", respectively [15,16]. In either case, a nonlinear term must be added to the equation to limit the growth of the instability; its derivation is the same for either type of instability, and, as with the linear terms, only x-derivatives of φ may appear. The nonlinear term can be derived by a simple geometrical argument, considering the propagation of a tilted front [13]. Such an argument leads to a nonlinear term which in our scaling takes the form (∂φ/∂x) 2 /2, as in (2) and (5).
We have shown that (2) may be regarded as a model equation for the transverse instability of a planar front at finite wavenumber. More generally, it is apparent that any physical system that gives rise to the Kuramoto-Sivashinsky equation might also lead to the Nikolaevskiy equation, depending on which of the two possible forms shown in Figure 1 is taken by the dispersion relation. Another example [3,4,17] is the derivation of nonlinear phase equations; this leads to the Nikolaevskiy and Kuramoto-Sivashinsky equations, in different parameter regimes, as the dispersion relation changes from type I to type II. The argument in the case of phase equations is essentially identical to that given above for fronts, following from the invariance under addition of a constant to the phase, and under reflection in x.
There are many studies of the transverse linear stability of fronts in the literature, arising from a variety of different physical and chemical systems. condensed matter heated by a laser, Anisimov et al. [18] found a linear spectrum of type I for transverse instability. Both types of instability were found in an experimental and theoretical study of a reaction front between two chemicals in a Hele-Shaw cell [15,19,20].

III. THE DAMPED NIKOLAEVSKIY EQUATION
The complex chaotic behavior of the Nikolaevskiy equation at onset arises from the interaction between the unstable modes near k = 1 and the neutral mode at k = 0. However, in real physical systems it is likely that the symmetry giving rise to the neutral mode (such as translation or Galilean symmetry) will be weakly broken. For example, in the case of a propagating front, considered in Section II, the translation symmetry may be broken by the presence of boundaries, or by a slight variation in the basic state with position. With such a weak symmetry breaking, the k = 0 mode will not be truly neutral, but will be weakly damped instead.
A further motivation for including a damping term in the Nikolaevskiy equation is that it allows us to study the transition from ordered stable patterns to soft turbulence. Such a probing of the origins of the turbulent state is not possible in the original Nikolaevskiy equation (1), since chaos is observed directly at onset (provided the domain is sufficiently large). Introduction of a damping term gives another parameter, which can be used to control and investigate the onset of chaos. A similar approach has been used to examine the dynamics of the Kuramoto-Sivashinsky equation [21,22].
We thus consider in this paper the "damped Nikolaevskiy" equation with damping coefficient ν ≥ 0; in particular we shall be concerned with the dynamics of (6) near the onset of pattern formation. In examining the stability of roll solutions, we shall find it useful to consider various cases for the relative sizes of the supercriticality parameter r = 2 and the damping coefficient ν, and hence we write where µ = O(1); we shall examine below a variety of pertinent values for s. Of course there are many ways in which one might add damping to the Nikolaevskiy equation. The damping term in (6) is chosen so that all modes are linearly damped (with the exception of the mode with wavenumber k = 1); large-scale modes decay at a rate −ν. Furthermore, significant analytical simplifications follow from the fact that in (6) the onset of linear instability is independent of ν (i.e., the critical value of r is r c = 0 and the critical wavenumber is k c = 1, for any ν ≥ 0). Note that our choice of damping term differs from that of Tribelsky [12], who considers instead a damping term proportional to −u − ∂ 2 u/∂x 2 . Consequently, in his formulation both r c and k c become functions of ν, leading to additional algebraic complications. However, the two formulations are equivalent under appropriate rescalings of x and t, and appropriate mappings between the two parameter sets.

IV. SECONDARY STABILITY OF ROLLS: NUMERICAL RESULTS
We begin our investigation of (6) by numerically computing spatially periodic steady "roll" solutions, and determining their secondary stability. To do so, first we fix and compute the roll solutionū(x) for a given choice of k and ν by integrating (6) forwards in time in a box of length 2π/k. In our pseudospectral numerical code to achieve this, u is approximated by the truncated Fourier series Disturbances u (x) to this roll solution then satisfy, from (6), and these disturbances are sought in the form (where the sum is again truncated for numerical purposes). Equation (8) is thus reduced to a matrix eigenvalue problem for the growth rate σ. By computing σ(p; k, , ν) for p in the range −k/2 ≤ p ≤ k/2, we are then able to determine the stability of the roll solutionū(x). Results for = 0.1 and = 0.05 are summarized in Figure 2, where rolls are stable inside the regions bounded by solid curves. For large enough damping (i.e., towards the top of each plot), the familiar Eckhaus stability boundaries are recovered. For sufficiently small damping, by contrast, all rolls are unstable -this shows that the known instability of all rolls in the undamped Nikolaevskiy equation extends to the case of (sufficiently small) finite damping, for a fixed value of the supercriticality parameter. Intermediate numerical results show the transition between the two limits.
Several features of Figure 2 merit comment. First, the region of stable rolls, which is roughly symmetrical about k = 1 at large damping, becomes distinctly asymmetrical as the damping is reduced; rolls with k > 1 tend to be less stable than those with k < 1. Second, stability is finally lost, as ν is decreased, in a small cusp near the left-hand side of the marginal stability curve. For = 0.1, this corner is at ν ≈ 0.02 and k ≈ 0.950 (the roll existence boundary is, for this value of ν, at k ≈ 0.9493); for = 0.05, we have found it much harder to compute the very thin cusp right down to its tip. The left and right sides of the cusp correspond, respectively, to monotonic and oscillatory long-wavelength instabilities. In fact the entire boundary corresponds to long-wavelength instabilities (|p| 1), with the exception of the segment marked "f ", where the instability has finite wavelength (p = O(1)). For even smaller values of than shown in the figure, this finite-wavelength segment is absent, a point to which we shall return in Section V C 1.
A more conventional picture of the stability balloon for rolls is shown in Figure 3, where the damping parameter ν is fixed to be 0.02 and the stability of rolls is shown in the (k, r)-plane. Again rolls are stable inside the solid curves. The region of stable rolls lies almost entirely in the region k < 1. For r > 0.003, all rolls are unstable, except those exceedingly close to the left-hand part of the marginal stability curve; these rolls are last to lose stability as r is increased. Finally, for r > 0.01, all rolls are unstable. A corresponding picture for smaller damping coefficient ν would show stable rolls in a similar, but smaller region, confined closer to the point of onset (r = 0, k = 1). In the next section, we shall argue on the basis of asymptotic arguments that in the limit ν → 0 (approaching the undamped Nikolaevskiy equation), the size of the stable region shrinks to zero. Crucially, however, we expect a small, but finite, region of stable rolls in the (k, r)-plane for any nonzero value of the damping coefficient ν. It is only for the undamped case ν = 0 that all rolls are unstable at onset.  (6), computed numerically for (a) = 0.1 and (b) = 0.05. Rolls with wavenumber k exist between the dashed curves, and are stable inside the region enclosed by solid curves. The annotations on the stability boundary indicate the type of instability suffered by the rolls upon crossing that part of the boundary: modulational steady ("s"), modulational Hopf ("H ") or finite-wavelength ("f ") instability. The crosses correspond to numerical simulations of the initial-value problem (6) in a large computational box and indicate the most extreme values of the wavenumber for which rolls are found to be stable to small perturbations (necessarily restricted to those that fit into the computational box); that all crosses lie inside the roll stability boundary provides a consistency check on our results. (6), computed numerically for ν = 0.02. Rolls with wavenumber k exist between the dashed curves, and are stable inside the region enclosed by solid curves.

V. SECONDARY STABILITY OF ROLLS: ANALYTICAL RESULTS
The numerical results in the previous section reveal a number of different instabilities of rolls in the damped Nikolaevskiy equation. In this section, we consider in detail three scalings for ν in (6) which shed particular light on the numerical stability results illustrated above: these correspond to choosing in turn s = 1, s = 3/2 and s = 2 in (7). The case s = 1 represents strong damping, in the sense that the decay rate ν = µ of the large-scale mode is O( ), which is much greater than the O( 2 ) growth rate of the pattern-forming mode. The effect of this strong damping is simply to modify the usual Ginzburg-Landau amplitude equation and thus modify the Eckhaus stability boundary, as we now demonstrate. Our conclusions below are consistent with those obtained by Tribelsky [12], although his analysis proceeds directly from a substitution of (9) into (8), rather than the amplitude-equation framework developed below.
We apply to (6) the usual weakly nonlinear scalings and expand the solution as where X = x and T = 2 t, and where c.c. denotes the complex conjugate of the preceding term.
, an equation arises for the mean mode f , which represents a balance between the driving term (|A| 2 ) X and the damping term −µf : Note that the linear terms f T and f XX do not appear at this order. The governing equation for A is also obtained at O( 3 ), by applying the usual solvability condition. It turns out to be from which f may be eliminated using (11) to give a single amplitude equation for A in the form FIG. 4: Stability boundaries of rolls according to (19), in (q, µ)-space, where the damping coefficient ν = µ. Rolls exist for −1/2 < q < 1/2 and are stable in the region between the two curves, where indicated. In the limit µ → 0, the stable rolls are those with q < 0. In the limit as µ → ∞, the usual Eckhaus stability boundary is recovered, so that rolls are stable for This is the usual Ginzburg-Landau equation for A, but with an additional nonlinear derivative term arising from the damping. Extensions to the Ginzburg-Landau equation including terms such as ı(|A| 2 ) X A are well known [23,24], but the different scalings involved here have the consequence that no further terms need be included to ensure that all terms of a given asymptotic order are present: with our scalings, (13) is complete. Note that in the limit µ → ∞, we recover the usual Ginzburg-Landau equation for A.
It is straightforward to analyse the stability of patterns in (13) -see Mancebo and Vega [25]. Steady patterns with A = a 0 exp ıqX exist for q 2 < q 2 c = 1/4, with the real amplitude a 0 = 6(1 − 4q 2 ) 1/2 . The stability of these solutions may be determined by writing A = (a 0 + a(X, T )) exp ıqX and linearizing in the perturbation a, to give After writing a = b + ıc, and separating (14) into real and imaginary parts, we find that b and c obey Then by supposing b and c each to be proportional to exp(λT + ılX), we find that the growth rate λ satisfies only monotonic instability is possible, and there can be instability if and only if It is clear from (18) that the most "dangerous" modes are those with small l, so the pattern is stable if The corresponding stability boundaries are shown in Figure 4.
In the limit µ → ∞, the usual Eckhaus stability condition (that rolls are stable for q 2 < q 2 e ≡ 1 3 q 2 c = 1 12 ) is recovered from (19). As µ is decreased, however, an asymmetry develops in the stability region: patterns are stable for f (µ) < q < g(µ), where f (µ) ∈ [− 1 2 , −q e ) and g(µ) ∈ [0, q e ) are each increasing functions of µ. When µ is small, the stability condition is roughly q(1 − 4q 2 ) < 0, so that patterns with q > 0 are unstable and patterns with q < 0 are stable. A corresponding conclusion is reached by Tribelsky (whose results are plotted in his Figure 1(a)). If we compare the results of this analysis with the numerical stability calculations described above in Section IV, we see that the stability boundaries in Figure 4 correspond, roughly, to those in Figure 2 above ν ≈ 0.2 (for = 0.1) and above ν ≈ 0.07 (for = 0.05).
The above analysis shows that for small µ, patterns with q < 0 are stable in (13). However, numerical simulations of (13) in a periodic domain, started from a random, small amplitude initial condition do not always lead to stable patterns. A typical simulation, for µ = 0.1, in a domain of size 100, is shown in Figure 5 (left). The solution is dominated by a state with 8 rolls in the domain, with a wavenumber q = −8 × 2π/100 ≈ −0.5026, which is just outside the range of wavenumbers for which steady patterns exist. This mode must therefore decay, until a bursting event occurs, involving other wavenumbers, from which the 8-roll state again emerges. Simulations with different domain sizes show that (13) seems to exhibit a preference for modes near to the left-hand limit of the region of existence of steady states, q = −1/2. Note that this behavior is shared by the damped Nikolaevskiy equation (6), as shown in Figure 2. Another qualitative similarity between the behavior of (13) and the Nikolaevskiy equation is that both exhibit irregular bursting behavior (albeit with rather different detailed structures). For comparison, a simulation of the undamped Nikolaevskiy equation (1) in a domain of size 200, for r = 0.01, is shown in Figure 5 (right).
In the limit µ → 0, it is clear from (17) that the eigenvalues λ become large, of order µ −1/2 . This indicates that a different scaling is required in which the damping ν is smaller than O( ) and the growth rate of the instability is greater than O( 2 ); this next scaling is considered in the following section.
B. Intermediate damping: s = 3/2 The analysis above corresponds to damping that is sufficiently strong to enslave the large-scale mode to gradients of the pattern amplitude. More subtle effects of damping can be explored if it is rather weaker: if s = 3/2 in (6). Our starting point is now the general roll solution of (6) with ν = 3/2 µ, which may be written in the form where a 0 = 6(1 − 4q 2 ) 1/2 .
A proper treatment of the stability of rolls requires consideration of the evolution of both amplitude and phase of the rolls, together with a large-scale mode. All three modes couple together, and their relative scalings, together with appropriate length and time scales for their evolution, are chosen below to enable a balance in the linearized perturbation problem; this is essentially the balance considered in [7]. However, by making appropriate absolute scalings for the three perturbation quantities, we are able to extend the linear results in [7] and accommodate the first nonlinear term in the evolution equations that follow. We shall see below that this enables us to describe the sub-or supercritical nature of the onset of instability of the rolls. It turns out that the appropriate scaling for the three perturbation quantities is accomplished by writing u = (a 0 + 1/2 a(X, T ))e ıx e ıq x e ı 1/4 φ(X,T ) + c.c. + · · · + 7/4 f (X, T ) + · · · , where now Note that in order to simplify notation we use the same symbols (X and T ) in different sections to represent different length and time scales. Notation within any one section is consistent, so no confusion should arise. Then after a systematic substitution of (20) in (6) and a consideration of terms at successive powers of 1/4 we eventually find The linear terms in this equation may be identified with the linear system considered in [7], with the addition here of the new damping term −µf in (22) (the quantities b, c and f in [7] correspond, respectively, to our a, a 0 φ and f ); the nonlinear term was not computed in [7]. These three equations may then be reduced to the nonlinear phase equation In analysing (24), it is useful to consider first the linearized problem, for which the right-hand side is simply −16a 2 0 qφ XX . This term represents the coupling between φ, a and f : in its absence, no instability can result in (24). Thus in Sections V B 1 and V B 2 below we assume a 2 0 q = O(1). However, we should note that the analysis described below must be reconsidered when |a 2 0 q| is small (cf. [8]) -we see immediately that there are two circumstances in which a separate consideration is warranted (a 0 1 and |q| 1); we shall return to this important point below, in Section V C.

No damping
We begin by summarizing relevant results in the absence of damping. When µ = 0 we recover from (24) the secondary stability results obtained by Tribelsky and Velarde [8] and elsewhere by ourselves [7]: there is monotonic instability of the rolls with q > 0, and such rolls are unstable to disturbances with wavenumbers l in the range 0 < |l| < l cm ≡ (a 2 0 q) 1/4 ; there is oscillatory instability for q < 0, and here the instability strikes for 0 < |l| < l co ≡ (−2a 2 0 q/25) 1/4 . Thus all rolls are unstable at onset in the undamped problem. It is instructive to extend the linearized analysis by considering the weakly nonlinear development of these instabilities, according to (24). We shall focus on computing the coefficients of the nonlinear terms in Landau equations for the amplitudes of the disturbances; the signs of the real parts of these coefficients will indicate the direction of bifurcation to the perturbed state. We consider first the monotonic instability for q > 0, and examine the evolution of a disturbance whose wavenumber satisfies |l − l cm | = O(δ), where 0 < δ 1. Expanding where φ 1 = A(T 2 ) sin lX and T 2 = δ 2 T , we find at O(δ 3 ) that Thus the bifurcation is always subcritical and hence leads to the observed explosive growth of the instability in simulations of the undamped Nikolaevskiy equation [7,8]. For the oscillatory instability for q < 0, we instead suppose that |l − l co | = O(δ); again expanding φ as in (25), but now with where ω = (−48a 2 0 q/25) Since q < 0, both travelling waves and standing waves branch supercritically, with the travelling wave stable. In practice we expect the explosive development of the monotonic instability to dominate the smooth onset of the oscillatory instability in simulations of the undamped Nikolaevskiy equation, except under extremely controlled (and contrived) conditions.

Nonzero damping
For µ > 0, we again begin by examining the linearized behavior of (24). A monotonic instability again occurs for all rolls with q > 0, with the stability margin being given by l 4 cm + µl 2 cm − qa 2 0 = 0. It thus follows that l 4 cm (µ) = qa 2 0 − µl 2 cm < qa 2 0 = l 4 cm (0), and hence the rolls are unstable to a smaller range of disturbance wavenumbers than in the undamped case. An oscillatory instability arises for rolls with q < q 0 , where q 0 = −µ 2 /(2a 2 0 ); such rolls are unstable to disturbances with 0 < |l| < l co , where now 5l 2 co = −µ + (−2qa 2 0 ) 1/2 . We may rewrite the condition q < q 0 for instability as Then since the cubic function in (28) has a minimum value of −8 √ 3 at q = − √ 3/6, this oscillatory instability can occur only for µ 2 < 8 √ 3, that is, for sufficiently weak damping (this is consistent with no oscillatory instability being observed in the previous scaling of stronger damping in Section V A). In the limit of small µ the stability boundaries in (28) approach q = 0 and q = −1/2, so that all rolls with q < 0 are predicted to be unstable to an oscillatory disturbance, just as in the undamped case. The stability boundaries for this scaling are illustrated in Figure 6. Note that the curved line of the Hopf bifurcation is consistent with the numerical stability results for the full damped Nikolaevskiy equation (6) shown in Figure 2 (except, of course, that the analysis above does not capture the finite-wavelength instability denoted by "f " in Figure 2). Figure 6 suggests that two regions of stable rolls, near q = 0 and q = −1/2, extend right down to zero damping. But in each of these cusped regions, |a 2 0 q| = |6q(1 − 4q 2 ) 1/2 | 1, and we recall that the present scaling is not valid in this limit. Thus we should disregard the two cusps in Figure 6 at small damping and instead examine more closely the stability of rolls in this limit; to resolve the correct behavior of the stability boundary at small damping, we consider in the next section a further (and final) scaling for ν.

C. Weak damping: s = 2
The final small-asymptotic description of the roll stability boundary may be described by setting s = 2, so that ν = 2 µ. In this case, the damping rate of the mean mode is of the same order as the growth rate of the patternforming mode. There are two cases to consider where the scaling of the previous section breaks down, corresponding to |q| 1 and a 0 1.

Central corner (|q| 1)
As Tribelsky and Velarde [8] have pointed out for the case of no damping, the scalings of Section V B break down when q is small. To resolve the small-q behavior, we adopt Tribelsky and Velarde's scalings (but in an amplitudeequation framework) and consider u ∼ (a 0 + 2 a(X, T ))e ıx e ı 2 qx e ı φ(X,T ) + c.c.
where now X = x and T = 2 t and a 0 = 6. Note that whereas (20) concerns basic roll wavenumbers k = 1 + O( ), (29), by contrast, concerns the much narrow wavenumber band k = 1 + O( 2 ). We find, after a substitution of (29) in (6), These three equations may then be reduced to the single nonlinear phase equation We now consider the linearized version of (33) and examine the dispersion relation for infinitesimal disturbances proportional to exp(σT + ılX). It is helpful first to recall the stability results of Tribelsky and Velarde [8] for the case µ = 0, which are summarized in Figure 7, and which are a special case of the more general analysis to be presented below. As can be seen in the figure, all rolls are unstable for µ = 0.
Since q m (l, µ) > q m (l, 0), it follows that the stability margin tends to move to the right as the damping is increased, tending to stabilize the rolls (see Figure 8). Rolls undergo an oscillatory instability (σ = ıω) when with onset frequency given by ω 2 = 24l 4 + 2(4µ + 41)l 2 + 2µ.
The oscillatory stability boundary is plotted in Figure 8 for various values of µ. Damping is seen to shift the oscillatory stability curve to the left, again stabilizing the rolls. A further effect of nonzero damping is to stabilize rolls to oscillatory disturbances in the limit that the perturbation wavenumber l → 0. As can be seen from Figure 8, with increasing damping the two stability boundaries separate and eventually part, allowing a small band of stable rolls. We find that stable rolls exist for µ > µ c ≈ 8.445 (the corresponding critical value of ν is ν c = µ c 2 ); thus, for a given value of the supercriticality parameter, there is a threshold value of the damping coefficient to allow stable rolls near k = k c , as previously shown by Tribelsky [12]. Figure 9 shows the region of stable rolls in (µ, q)-space (thus the apparent small-q cusp in Figure 6, which appears to extend down to zero damping, is in fact a wedge terminating at finite damping).
The results above seem at odds with the numerical results presented in Figure 2, since there is no central wedge in the numerical stability boundary. However, it turns out that for the values of used in Figure 2, the small-q wedge is masked by an additional finite-wavenumber instability. For a smaller value of , the theory in this section does indeed correctly predict the shape of the stability boundary, as shown in Figure 10 for = 0.01. (It becomes increasingly challenging to compute the entire stability boundary numerically, so only the relevant part is presented for this value of .) Rolls are stable in the wedge down to ν ≈ 8.445 2 , and the last rolls to destabilize have wavenumber approximately 1 − 4.049 2 . These results agree with corresponding results of Tribelsky [12].   (6), computed numerically for = 0.01, near the central wedge region corresponding to Figure 9. Crosses indicate numerical results and the solid line corresponds to our analytical results (cf. Figure 9). Note that for this value of (and presumably also for even smaller values) the numerical stability boundary contains this central wedge rather than a section corresponding to that marked "f " in Figure 2.

Left-hand corner (a0 1)
We now examine the second case in which the asymptotic results of Section V B break down: a 0 1. This regime proves to include the last rolls to become destabilized as the damping is reduced. We thus investigate the stability of rolls close to the marginal stability boundary. Recall from Section V B and in particular Figure 6 that we expect stable rolls to persist to small values of the damping near the left-hand side of the marginal stability curve.
We introduce the notation k = k − m , k + m for the left-and right-hand marginal stability boundaries, respectively (that is, for given values of r and ν, rolls exist for wavenumbers k in the range k − m < k < k + m ). We shall begin by considering both left-and right-hand parts of the marginal stability curve, since both lead to a 0 1. Our analysis will then show that only the left-hand part of the marginal curve is relevant, being capable of supporting stable rolls. We find, from substitution in the equation It turns out that the correct scaling to resolve the stability boundary of the rolls is obtained by setting |k−k ± m | = O( 2 ), in which case |ū| = O( 3/2 ), with the relevant space and time scales for the evolution of perturbations being X = x and T = 2 t. We consider here only the linearized secondary stability problem and we are thus at liberty to scale arbitrarily the absolute sizes of the three perturbation quantities, keeping their relative sizes fixed. Thus we suppose that u ∼ 3/2 (a 0 + a(X, T ))e ıkx e ıφ(X,T ) + c.c.
where the wavenumber of the rolls under consideration is just inside the marginal curve, so that A weakly nonlinear calculation of the roll solution itself shows that the amplitude and wavenumber are related by Then a consideration of the terms in (6) which are linear in the perturbation quantities φ, a and f shows that their evolution is governed by These three equations may then be combined to yield a single linearized phase equation in the form For modes proportional to exp(σT + ılX), this gives the dispersion relation In the limit a 0 → 0, (41) gives the growth rates σ = −(µ + l 2 ) < 0, σ = 4l(1 − l) and σ = 4l(−1 − l), and hence all rolls are unstable (to a monotonic disturbance) sufficiently close to the marginal curve. We gain some insight into the possible existence of stable rolls near the marginal curve by next considering the limit of small perturbation wavenumbers, l → 0. For fixed a 0 and µ, in the limit l → 0, the leading-order balance in (41) is Thus near k = k + m rolls are unstable, to monotonic disturbances, since (2 + a 2 0 /µ) > 0. No further analysis of this case is necessary: there is no pocket of stable rolls possible near the right-hand marginal stability curve; this conclusion is consistent with Figure 2.
Near k = k − m , by contrast, the conclusion of monotonic instability follows from (42) only if a 2 0 < 2µ (i.e., sufficiently close to the marginal curve). When a 2 0 > 2µ, (42) instead indicates oscillations, and we must go further in our consideration of the small-l problem in order to determine the stability of the rolls. By carrying out a small-l expansion of the growth rate in the form σ = σ 1 l + σ 2 l 2 + · · ·, we find from (41) that σ 2 1 = 8(2 − a 2 0 /µ) and σ 2 = 4(a 2 0 /µ 2 − 1). Thus when rolls are stable in the small-l limit (the lower threshold for a 2 0 corresponding to a monotonic instability, the upper threshold to an oscillatory instability). The relation (43) indicates that rolls enjoy this region of stability only for sufficiently large damping: µ > µ c ≡ 2. Note that since 2 < 8.445 this cusp contains the last stable rolls as the damping is decreased.
In terms of the original variables in the damped Nikolaevskiy equation (6), we thus predict (for small ν) the cusp to be at ν ∼ 2r = 2 2 . We may readily check these analytical results against the numerical secondary stability calculations of Section IV. For = 0.1, as in Figure 2, we thus predict the apex of the cusp to lie at ν = 0.02 (which we do indeed find numerically). Similarly, for = 0.05 we predict the apex to lie at ν = 0.005, and the smallest value of ν for which we have been able to compute the cusp is ν = 0.006, which provides reasonable agreement with the theory. (As → 0 it becomes exceedingly difficult to compute accurately this part of the numerical secondary stability boundary.) For completeness, we note that we have also considered numerically the full dispersion relation (41) for arbitrary values of l: we find that (at least for values of a 0 and µ up to 4, for which we have carried out calculations) the stability boundary is indeed determined by the small-l expansion.

VI. NUMERICAL SIMULATIONS
In this section we present full numerical simulations of the damped Nikolaevskiy equation (6), illustrating the transition from regular patterns to soft-mode turbulence as the damping is reduced. The simulations use periodic boundary conditions and employ a Fourier spectral method for the spatial discretization. The time-stepping is an explicit second-order Runge-Kutta form of the exponential time differencing method [26] that computes the stiff linear part of the equation exactly, permitting O(1) time steps. The size of the domain is 200π, allowing a resolution of 0.01 in wavenumber.
In the simulations shown in Figure 11, the initial condition comprises rolls with wavenumber k = 0.96, plus a small random perturbation. In each case, the simulation was run for a long time to remove transient features, and then plotted over a subsequent relatively short time interval. Note that these time intervals are not the same in each case. The driving parameter was fixed at r = 0.01 and the sequence of simulations shows the effect of reducing ν. For these parameters, all rolls that fit into the computational domain are predicted to be unstable for ν < ν c where ν c ≈ 0.057. The bifurcation at ν = ν c is the Hopf bifurcation whose stability boundary is shown in Figures 2 and 6.
The simulations confirm the predicted bifurcation type, since for ν < ν c a small oscillatory modulation to the rolls grows. This bifurcation seems to be supercritical since the oscillations equilibrate to a small periodic modulation of the rolls. For ν = 0.04 this oscillation takes the form of a standing wave, as shown in Figure 11 (top). At ν = 0.03, the solution is no longer periodic in time and shows occasional irregular bursts of instability. In this case, the dominant wavenumber is k = 0.95, at the left-hand boundary of the roll-existence region (i.e., near k = k − m ), consistent with the analysis of the preceding section and the behavior of (13). Also, counter-intuitively, the amplitude of the solution becomes significantly lower as the damping is decreased. This is because of the transition from the O( ) scaling in the strongly damped case to the O( 3/2 ) scaling [7] in the undamped, unstable regime. For ν = 0.02, the bursts are more frequent and a broad spectrum of wavenumbers is present in the solution. Finally, the lowest panel of Figure 11 shows the case of zero damping. This is significantly different from the weakly damped case, showing grain boundaries between regions of travelling rolls.

VII. DISCUSSION AND CONCLUSIONS
The Nikolaevskiy equation is an important model of a wide range of physical systems, including certain convection problems, phase instabilities and transverse instabilities of fronts. It arises naturally in systems with a finitewavenumber instability and a translation symmetry for the dependent variable. However, in many applications it is likely that this symmetry will be weakly broken.
The damped Nikolaevskiy equation (6) allows us to investigate the effects of weak symmetry breaking, and describe the transition between the appearance of soft-mode turbulence at onset when ν = 0 and the more gradual development of complex dynamics typical of damped systems when the damping coefficient ν > 0. Our numerical investigation of the secondary stability problem for steady roll solutions shows how the instability of all roll states in the undamped case evolves into the more common Eckhaus scenario, whereby rolls with wavenumbers sufficiently close to critical are stable, when the damping is sufficiently great. If we fix the supercriticality parameter r, then it follows from the asymptotic results of Section V C 2 that there is a critical value of the damping coefficient ν = ν c (given by ν c ∼ 2r when r is small) below which all rolls are unstable; for ν > ν c some rolls are stable.
However, it is more common in applications to fix parameters such as the damping coefficient ν and vary the supercriticality parameter r. Then the results of Section V C 2 may be interpreted as indicating that, at least for small values of ν, there are stable rolls provided 2r < ν, i.e., close enough to onset. Thus the damped Nikolaevskiy equation differs fundamentally from the undamped version in that there is not a direct transition to soft-mode turbulence at onset when ν > 0; instead, there is a sequence of bifurcations beginning with some initially stable roll state. Of course, when ν is small this bifurcation sequence may occur very close to onset (i.e., for r = O(ν)).
Our work complements a recent study by Tribelsky [12], but differs in some significant respects. The first and most important difference is our discovery and analysis of the cusp of small-amplitude stable rolls, as described in Section V C 2. These rolls are the last to be destabilized as the damping is reduced, and so are particularly significant. Another difference is that our numerical study of the secondary stability problem has shown that rolls in the central wedge of apparent stability, which is predicted both by ourselves and by Tribelsky using asymptotic arguments, are in fact unstable to short-wave modes, except for exceedingly small values of ν. So this central wedge (see Section V C 1) in fact forms part of the stability boundary only for very small damping.
Finally, we reiterate that our choice of damping term in (6) is by no means unique: indeed, Tribelsky [12] made a different choice, as discussed in Section III. We conclude by noting that we have also computed numerically the secondary stability boundaries corresponding to those shown in Figure 2, but for Tribelsky's version of the damped Nikolaevskiy equation, and find similar results.