Pattern formation in the damped Nikolaevskiy equation

# Pattern formation in the damped Nikolaevskiy equation

## Abstract

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

 ∂u∂t+u∂u∂x=−∂2∂x2(ru−(1+∂2∂x2)2u) (1)

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

 ∂ϕ∂t+12(∂ϕ∂x)2=−∂2∂x2(rϕ−(1+∂2∂x2)2ϕ) (2)

after introducing .

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 pattern-forming 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 , 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 (7); (6); (9). In Ref. (7) we showed that in this turbulent regime, (1) exhibits an unusual scaling, with the amplitude of proportional to the power of the supercriticality parameter , whereas it is more usually proportional to the 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 in the direction and that the 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 . 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 -derivatives of . Furthermore, reflection symmetry in means that the only permissible linear terms involve even derivatives of . Therefore, the linearized equation for must have the form

 ∂ϕ∂t=a2∂2ϕ∂x2+a4∂4ϕ∂x4+a6∂6ϕ∂x6+⋯, (3)

where the are constants. The corresponding dispersion relation is

 λ=−a2k2+a4k4−a6k6+⋯, (4)

for the growth of modes . From (4) it is apparent that there are two possible types of instability of the front. If and , then the planar front is unstable to a band of wavenumbers near (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 and , may be written in the form

 ∂ϕ∂t+12(∂ϕ∂x)2=−∂2∂x2(rϕ+∂2ϕ∂x2). (5)

But if , and 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 and , (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 -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 , 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 .

There are many studies of the transverse linear stability of fronts in the literature, arising from a variety of different physical and chemical systems. A review of the literature reveals several examples of a finite-wavenumber (type I) instability, consistent with the Nikolaevskiy equation. For example, in their study of an evaporation front for 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 and the neutral mode at . 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 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

 ∂u∂t+u∂u∂x=−∂2∂x2(ru−(1+∂2∂x2)2u)−ν(1+∂2∂x2)2u, (6)

with damping coefficient ; 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

 ν=ϵsμ, (7)

where ; we shall examine below a variety of pertinent values for .

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 ); 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 is and the critical wavenumber is , for any ). Note that our choice of damping term differs from that of Tribelsky (12), who considers instead a damping term proportional to . Consequently, in his formulation both and become functions of , leading to additional algebraic complications. However, the two formulations are equivalent under appropriate rescalings of and , 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 for a given choice of and by integrating (6) forwards in time in a box of length . In our pseudospectral numerical code to achieve this, is approximated by the truncated Fourier series

 ¯u(x)=N∑−N¯unenikx.

Disturbances to this roll solution then satisfy, from (6),

 ∂u′∂t+u′∂¯u∂x+¯u∂u′∂x=−∂2∂x2(ru′−(1+∂2∂x2)2u′)−ν(1+∂2∂x2)2u′, (8)

and these disturbances are sought in the form

 u′(x)=eσt+ipxN∑−Nu′nenikx (9)

(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 for in the range , we are then able to determine the stability of the roll solution .

Results for and 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 at large damping, becomes distinctly asymmetrical as the damping is reduced; rolls with tend to be less stable than those with . Second, stability is finally lost, as is decreased, in a small cusp near the left-hand side of the marginal stability curve. For , this corner is at and (the roll existence boundary is, for this value of , at ); for , 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 (), with the exception of the segment marked “f”, where the instability has finite wavelength (). 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.3.1.

A more conventional picture of the stability balloon for rolls is shown in Figure 3, where the damping parameter is fixed to be and the stability of rolls is shown in the -plane. Again rolls are stable inside the solid curves. The region of stable rolls lies almost entirely in the region . For , 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 is increased. Finally, for , 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 (, ). In the next section, we shall argue on the basis of asymptotic arguments that in the limit (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 -plane for any nonzero value of the damping coefficient . It is only for the undamped case that all rolls are unstable at onset.

## 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 , and in (7).

### v.1 Strong damping: s=1

The case represents strong damping, in the sense that the decay rate of the large-scale mode is , which is much greater than the 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

 u=ϵ(A(X,T)eix+c.c.)+ϵ2u2+ϵ3u3+⋯, (10)

where and , and where c.c. denotes the complex conjugate of the preceding term. At , (6) is satisfied. At , we find

 u2=(−136iA2e2ix+c.c.)+f(X,T),

where the large-scale mode is arbitrary. At , an equation arises for the mean mode , which represents a balance between the driving term and the damping term :

 0=−(|A|2)X−μf. (11)

Note that the linear terms and do not appear at this order. The governing equation for is also obtained at , by applying the usual solvability condition. It turns out to be

 AT=A+4AXX−136|A|2A−ifA, (12)

from which may be eliminated using (11) to give a single amplitude equation for in the form

 AT=A+4AXX−136|A|2A+i(|A|2)XA/μ. (13)

This is the usual Ginzburg–Landau equation for , but with an additional nonlinear derivative term arising from the damping. Extensions to the Ginzburg–Landau equation including terms such as 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 .

It is straightforward to analyse the stability of patterns in (13) – see Mancebo and Vega (25). Steady patterns with exist for , with the real amplitude . The stability of these solutions may be determined by writing and linearizing in the perturbation , to give

 aT=−136a20(a+a∗)+8iqaX+4aXX+ia20(aX+a∗X)/μ. (14)

After writing , and separating (14) into real and imaginary parts, we find that and obey

 bT = −118a20b−8qcX+4bXX, (15) cT = 8qbX+4cXX+2a20bX/μ. (16)

Then by supposing and each to be proportional to , we find that the growth rate satisfies

 λ2+(118a20+8l2)λ+16l4+29a20l2−64q2l2−16ql2a20/μ=0; (17)

only monotonic instability is possible, and there can be instability if and only if

 l2(72l2+a20−288q2−72qa20/μ)<0. (18)

It is clear from (18) that the most “dangerous” modes are those with small , so the pattern is stable if

 (1−12q2)>72q(1−4q2)/μ. (19)

The corresponding stability boundaries are shown in Figure 4.

In the limit , the usual Eckhaus stability condition (that rolls are stable for ) is recovered from (19). As is decreased, however, an asymmetry develops in the stability region: patterns are stable for , where and are each increasing functions of . When is small, the stability condition is roughly , so that patterns with are unstable and patterns with 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 (for ) and above (for ).

The above analysis shows that for small , patterns with 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 , 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 , 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, . 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 , is shown in Figure 5 (right).

In the limit , it is clear from (17) that the eigenvalues become large, of order . This indicates that a different scaling is required in which the damping is smaller than and the growth rate of the instability is greater than ; this next scaling is considered in the following section.

### v.2 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 in (6). Our starting point is now the general roll solution of (6) with , which may be written in the form

 u=ϵa0e(1+ϵq)ix+c.c.+O(ϵ2),

where .

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=ϵ(a0+ϵ1/2a(X,T))eixeiqϵxeiϵ1/4ϕ(X,T)+c.c.+⋯+ϵ7/4f(X,T)+⋯, (20)

where now

 X=ϵ3/4x,T=ϵ3/2t.

Note that in order to simplify notation we use the same symbols ( and ) 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 we eventually find

 ∂ϕ∂T = 4∂2ϕ∂X2−f, (21) ∂f∂T = ∂2f∂X2−μf−2a0∂a∂X, (22) ∂a∂T = 4∂2a∂X2−4a0(∂ϕ∂X)2−8qa0∂ϕ∂X. (23)

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 in (22) (the quantities , and in (7) correspond, respectively, to our , and ); the nonlinear term was not computed in (7). These three equations may then be reduced to the nonlinear phase equation

 (∂∂T−4∂2∂X2)2(∂∂T−∂2∂X2+μ)ϕ=−16a20(q+∂ϕ∂X)∂2ϕ∂X2. (24)

In analysing (24), it is useful to consider first the linearized problem, for which the right-hand side is simply . This term represents the coupling between , and : in its absence, no instability can result in (24). Thus in Sections V.2.1 and V.2.2 below we assume . However, we should note that the analysis described below must be reconsidered when is small (cf. (8)) – we see immediately that there are two circumstances in which a separate consideration is warranted ( and ); we shall return to this important point below, in Section V.3.

#### No damping

We begin by summarizing relevant results in the absence of damping. When 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 , and such rolls are unstable to disturbances with wavenumbers in the range ; there is oscillatory instability for , and here the instability strikes for . 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 , and examine the evolution of a disturbance whose wavenumber satisfies , where . Expanding

 ϕ=δϕ1+δ2ϕ2+⋯, (25)

where

 ϕ1=A(T2)sinlX

and , we find at 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 , we instead suppose that ; again expanding as in (25), but now with

 ϕ1=eiωT(B(T2)eilX+C(T2)e−ilX)+c.c.,

where , we find at the Landau equations

 −qa20BdBdT2 = (−74072639+55√68071i)|B|2+(−401197−10√6399i)|C|2, (26) −qa20CdCdT2 = (−74072639+55√68071i)|C|2+(−401197−10√6399i)|B|2. (27)

Since , 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 , we again begin by examining the linearized behavior of (24). A monotonic instability again occurs for all rolls with , with the stability margin being given by . It thus follows that

 l4cm(μ)=qa20−μl2cm

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 , where ; such rolls are unstable to disturbances with , where now . We may rewrite the condition for instability as

 72q(1−4q2)<−μ2. (28)

Then since the cubic function in (28) has a minimum value of at , this oscillatory instability can occur only for , 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.1). In the limit of small the stability boundaries in (28) approach and , so that all rolls with 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 and , extend right down to zero damping. But in each of these cusped regions, , 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 .

### v.3 Weak damping: s=2

The final small- asymptotic description of the roll stability boundary may be described by setting , so that . In this case, the damping rate of the mean mode is of the same order as the growth rate of the pattern-forming mode. There are two cases to consider where the scaling of the previous section breaks down, corresponding to and .

#### Central corner (|q|≪1)

As Tribelsky and Velarde (8) have pointed out for the case of no damping, the scalings of Section V.2 break down when is small. To resolve the small- behavior, we adopt Tribelsky and Velarde’s scalings (but in an amplitude-equation framework) and consider

 u∼ϵ(a0+ϵ2a(X,T))eixeiϵ2qxeiϵϕ(X,T)+c.c.+⋯+ϵ3f(X,T)+⋯, (29)

where now and and . Note that whereas (20) concerns basic roll wavenumbers , (29), by contrast, concerns the much narrow wavenumber band . We find, after a substitution of (29) in (6),

 ∂ϕ∂T = 4∂2ϕ∂X2−f, (30) ∂f∂T = ∂2f∂X2−μf−2a0∂a∂X, (31) ∂a∂T = 4∂2a∂X2−118a20a+a0(1154a20−8q+12∂2∂X2)∂ϕ∂X−4a0(∂ϕ∂X)2−a0∂f∂X. (32)

These three equations may then be reduced to the single nonlinear phase equation

 (∂∂T−4∂2∂X2)(∂∂T−4∂2∂X2+118a20)(∂∂T−∂2∂X2+μ)ϕ= 2a20(1154a20−8q+∂∂T+8∂2∂X2)∂2ϕ∂X2−16a20∂ϕ∂X∂2ϕ∂X2. (33)

We now consider the linearized version of (33) and examine the dispersion relation for infinitesimal disturbances proportional to . It is helpful first to recall the stability results of Tribelsky and Velarde (8) for the case , 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 .

For general values of , rolls undergo a monotonic instability () when , where

 72qm(l,μ)=2l4+(2μ−71)l2+66+μ.

Since , 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

 q=−μ(2+μ)288l−2−(μ272+3μ16−91144)−(5μ36+677288)l2−2572l4, (34)

with onset frequency given by

 ω2=24l4+2(4μ+41)l2+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 .

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 (the corresponding critical value of is ); thus, for a given value of the supercriticality parameter, there is a threshold value of the damping coefficient to allow stable rolls near , as previously shown by Tribelsky (12). Figure 9 shows the region of stable rolls in -space (thus the apparent small- 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- 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 . (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 , and the last rolls to destabilize have wavenumber approximately . These results agree with corresponding results of Tribelsky (12).

#### Left-hand corner (a0≪1)

We now examine the second case in which the asymptotic results of Section V.2 break down: . 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.2 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 for the left- and right-hand marginal stability boundaries, respectively (that is, for given values of and , rolls exist for wavenumbers in the range ). We shall begin by considering both left- and right-hand parts of the marginal stability curve, since both lead to . 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 for the marginal curve, that

 k±m=1±12ϵ−18ϵ2±(116−14μ)ϵ3+⋯.

It turns out that the correct scaling to resolve the stability boundary of the rolls is obtained by setting , in which case , with the relevant space and time scales for the evolution of perturbations being and . 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(a0+a(X,T))eikxeiϕ(X,T)+c.c.+⋯+ϵ2f(X,T)+⋯, (35)

where the wavenumber of the rolls under consideration is just inside the marginal curve, so that

 k=k±m∓ϵ2κ.

A weakly nonlinear calculation of the roll solution itself shows that the amplitude and wavenumber are related by

 a20=144κ. (36)

Then a consideration of the terms in (6) which are linear in the perturbation quantities , and shows that their evolution is governed by

 a0∂ϕ∂T = 4a0∂2ϕ∂X2−a0f±4∂a∂X, (37) ∂f∂T = ∂2f∂X2−μf−2a0∂a∂X, (38) ∂a∂T = 4∂2a∂X2∓4a0∂ϕ∂X. (39)

These three equations may then be combined to yield a single linearized phase equation in the form

 Missing or unrecognized delimiter for \right (40)

For modes proportional to , this gives the dispersion relation

 (σ+l2+μ)[(σ+4l2)2−16l2]=±8a20l2. (41)

In the limit , (41) gives the growth rates , and , 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, . For fixed and , in the limit , the leading-order balance in (41) is

 σ2∼8l2(2±a20/μ). (42)

Thus near rolls are unstable, to monotonic disturbances, since . 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 , by contrast, the conclusion of monotonic instability follows from (42) only if (i.e., sufficiently close to the marginal curve). When , (42) instead indicates oscillations, and we must go further in our consideration of the small- problem in order to determine the stability of the rolls. By carrying out a small- expansion of the growth rate in the form , we find from (41) that and . Thus when

 2μ

rolls are stable in the small- limit (the lower threshold for 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: . Note that since 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 . We may readily check these analytical results against the numerical secondary stability calculations of Section IV. For , as in Figure 2, we thus predict the apex of the cusp to lie at (which we do indeed find numerically). Similarly, for we predict the apex to lie at , and the smallest value of for which we have been able to compute the cusp is , which provides reasonable agreement with the theory. (As 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 : we find that (at least for values of and up to , for which we have carried out calculations) the stability boundary is indeed determined by the small- 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 time steps. The size of the domain is , allowing a resolution of in wavenumber.

In the simulations shown in Figure 11, the initial condition comprises rolls with wavenumber , 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 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 where . The bifurcation at is the Hopf bifurcation whose stability boundary is shown in Figures 2 and 6.

The simulations confirm the predicted bifurcation type, since for 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 this oscillation takes the form of a standing wave, as shown in Figure 11 (top). At , the solution is no longer periodic in time and shows occasional irregular bursts of instability. In this case, the dominant wavenumber is , at the left-hand boundary of the roll-existence region (i.e., near ), 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 scaling in the strongly damped case to the scaling (7) in the undamped, unstable regime. For , 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 finite-wavenumber 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 and the more gradual development of complex dynamics typical of damped systems when the damping coefficient . 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 , then it follows from the asymptotic results of Section V.3.2 that there is a critical value of the damping coefficient (given by when is small) below which all rolls are unstable; for some rolls are stable.

However, it is more common in applications to fix parameters such as the damping coefficient and vary the supercriticality parameter . Then the results of Section V.3.2 may be interpreted as indicating that, at least for small values of , there are stable rolls provided , 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 ; 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 ).

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

### References

1. I. A. Beresnev and V. N. Nikolaevskiy, Physica D 66, 1–6 (1993).
2. V. N. Nikolaevskiy, Springer Lecture Notes in Engineering, 39, 210–221 (1989)
3. H. Fujisaka and T. Yamada, Prog. Theor. Phys. 106, 315–322 (2001).
4. D. Tanaka, Phys. Rev. E 70, 015202(R) (2004).
5. Y. Hidaka, J.-H. Huh, K.-i. Hayashi, S. Kai and M. I. Tribelsky, Phys. Rev. E 56, R6256–R6259 (1997).
6. M. I. Tribelsky and K. Tsuboi, Phys. Rev. Lett. 76, 1631–1634 (1996).
7. P. C. Matthews and S. M. Cox, Phys. Rev. E 62, R1473–R1476 (2000).
8. M. I. Tribelsky and M. G. Velarde, Phys. Rev. E 54, 4973–4981 (1992).
9. H. Sakaguchi and D. Tanaka, Phys. Rev. E 76, 025201(R) (2007).
10. H. Fujisaka, T. Honkawa and T. Yamada, Prog. Theor. Phys. 109, 911–918 (2003).
11. Y. Hidaka, K. Tamura, S. Kai, Prog. Theor. Phys. Supp. 161, 1–11 (2006).
12. M. I. Tribelsky, Preprint (2006).
13. Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence. Springer, Berlin (1984).
14. G. I. Sivashinsky, 4 1177–1206 (1977).
15. T. Bánsági, D. Horvath, A. Toth, J. Yang, S. Kalliadasis and A. De Wit, Phys. Rev. E 68, 055301(R) (2003).
16. M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851–1112 (1993).
17. D. Tanaka and Y. Kuramoto, Phys. Rev. E 68, 026219 (2003).
18. S. I. Anisimov, M. I. Tribel’skii, Ya. G. Épel’baum, Sov. Phys. JETP 51, 802–806 (1980).
19. J. D’Hernoncourt, S. Kalliadasis and A. De Wit, J. Chem. Phys. 123, 234503 (2005).
20. S. Kalliadasis, J. Yang and A. De Wit, Phys. Fluids 16, 1395–1409 (2004).
21. P. Brunet, Phys. Rev. E 76, 017204 (2007).
22. C. Misbah and A. Valance, Phys. Rev. E 49, 166–183 (1994).
23. A Doelman, PhD thesis, University of Utrecht (1990).
24. R. B. Hoyle, Phys. Rev. E 58, 7315–7318 (1998).
25. F. J. Mancebo and J. M. Vega, J. Fluid Mech. 560, 369–393 (2006).
26. S. M. Cox and P. C. Matthews, J. Comp. Phys. 176, 430–455 (2002).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters