Wavenumber Selection in Pattern Forming Systems

TL;DR

This paper demonstrates through extensive numerical simulations that noise acts as a robust mechanism for selecting a unique wavenumber in pattern forming systems, independent of initial conditions or noise strength.

Contribution

It provides the first large-scale numerical evidence supporting noise-induced wavenumber selection in the stabilized Kuramoto-Sivashinsky equation.

Findings

Unique wavenumber state has highest probability at long times

Wavenumber selection is independent of noise strength

State is robust against initial conditions

Abstract

Wavenumber selection in pattern forming systems remains a long standing puzzle in physics. Previous studies have shown that external noise is a possible mechanism for wavenumber selection. We conduct an extensive numerical study of the noisy stabilized Kuramoto Sivashinsky equation. We use a fast spectral method of integration, which enables us to investigate long time behavior for large system sizes that could not be investigated by earlier work. We find that a state with a unique wavenumber has the highest probability of occurring at very long times. We also find that this state is independent of the strength of the noise and initial conditions, thus making a convincing case for the role of noise as a mechanism of state selection.

Tables

Table 1. Table 1 : Most visited wavenumbers for various sizes, α = 0.22 𝛼 0.22 \alpha=0.22 .

Number of lattice points $N$	System Length $\left(L=Nh\right)$	Most visited state
1024	512	0.6749 $\left(n=55\right)$
1600	800	0.6754 $\left(n=86\right)$
2200	1100	0.6740 $\left(n=118\right)$
3000	1500	0.6744 $\left(n=161\right)$
4000	2000	0.6754 $\left(n=215\right)$

Table 2. Table 2 : Most visited wavenumbers for various sizes, α = 0.20 𝛼 0.20 \alpha=0.20 .

Number of lattice points $N$	System Length $\left(L=Nh\right)$	Most visited state
1024	512	0.6627 $\left(n=54\right)$
1600	800	0.6597 $\left(n=84\right)$
2200	1100	0.6569 $\left(n=115\right)$
3000	1500	0.6576 $\left(n=157\right)$
4000	2000	0.6566 $\left(n=209\right)$

Table 3. Table 3 : Most visited wavenumbers for various sizes, α = 0.17 𝛼 0.17 \alpha=0.17 .

Number of lattice points $N$	System Length $\left(L=Nh\right)$	Most visited state
1024	512	0.6381 $\left(n=52\right)$
1600	800	0.6361 $\left(n=81\right)$
2200	1100	0.6340 $\left(n=111\right)$
3000	1500	0.6367 $\left(n=152\right)$
4000	2000	0.6377 $\left(n=203\right)$

Table 4. Table 4 : Most visited wavenumbers for various sizes, α = 0.24 𝛼 0.24 \alpha=0.24 .

Number of lattice points $N$	System Length $\left(L=Nh\right)$	Most visited state
1024	512	0.6995 $\left(n=57\right)$
1600	800	0.6990 $\left(n=89\right)$
2200	1100	0.6969 $\left(n=122\right)$
3000	1500	0.6953 $\left(n=166\right)$
4000	2000	0.6974 $\left(n=222\right)$

Table 5. Table 5 : Most visited or selected wavenumbers q s subscript 𝑞 𝑠 q_{s} in the thermodynamic limit.

Control Parameter $\alpha$	$q_s$
0.22	0.6748
0.20	0.6567
0.17	0.6384

Equations

$${\partial }_{t}u(x,t)={\hat{L}}_{p}u(x,t)+\hat{N}[u(x,t)]$$

$${\partial }_{t}u(x,t)=-\alpha u(x,t)-{\partial }_x^{2}u(x,t)-{\partial }_x^{4}u(x,t)+{({\partial }_{x}u(x,t))}^2$$

$$\sigma =-\alpha +q^2-q^4$$

$$1/2-\sqrt{1/4-\alpha }\le q^2\le 1/2+\sqrt{1/4-\alpha },\alpha \le {\alpha }_c$$

$${\partial }_{t}u(x,t)=-\alpha u-{\partial }_x^{2}u(x,t)-{\partial }_x^{4}u(x,t)+{({\partial }_{x}u(x,t))}^2+\zeta (x,t)$$

$$⟨\zeta (x,t)⟩=0$$

$$⟨\zeta (x,t)\zeta (x^{\prime },t^{\prime })⟩=2\epsilon \delta (x-x^{\prime })\delta (t-t^{\prime })$$

$$\dot{\Phi }(x,t)=f[\Phi (x,t)]+\zeta (x,t)$$

$$P_T[\varphi (x,t)]=\mathcal{N}\text{exp}[-S_T[\varphi (x,t)]/\epsilon ]$$

$$⟨\zeta (x,t)\zeta (x^{\prime },t^{\prime })⟩=2\epsilon \delta (x-x^{\prime })\delta (t-t^{\prime })$$

$$S_T[\varphi (x,t)]=\frac{1}{2}{\int }_0^T\text{d}t{\int }_0^L\text{d}x{[\dot{\varphi }(x,t)-f[\varphi (x,t)]]}^2$$

$$\frac{\partial u(x,t)}{\partial t}=\hat{L}u(x,t)+\hat{N}[u(x,t),{\partial }_{x}u(x,t),{\partial }_x^{2}u(x,t),\dots ]$$

$$\frac{d{u}_i}{dt}={(\hat{L}u)}_i+\hat{N}[u(x,t),\dots {]}_{i}i=0,1,\dots ,N-1.$$

$$\frac{d{\tilde{u}}_k}{dt}={\widetilde{(\hat{L}u)}}_k+{\tilde{\hat{N}}}_{k}k=0,1,\dots ,\frac{N}{2}-1,-N/2,\dots ,-1$$

$$\hat{L}u(x,t)=-\alpha u(x,t)-{\partial }_x^{2}u(x,t)-{\partial }_x^{4}u(x,t)$$

$$\hat{N}[u(x,t),{\partial }_{x}u(x,t)]={({\partial }_{x}u(x,t))}^2$$

$${\widetilde{(\hat{L}u)}}_k=-\alpha {\tilde{u}}_k+{(2\pi k/Nh)}^2{\tilde{u}}_k-{(2\pi k/Nh)}^4{\tilde{u}}_k$$

$$\frac{d{\tilde{u}}_k}{dt}=-\alpha {\tilde{u}}_k+{(2\pi k/Nh)}^2{\tilde{u}}_k-{(2\pi k/Nh)}^4{\tilde{u}}_k+{\tilde{\hat{N}}}_k[u(x,t),\dots ]$$

$$u_k^{j+1}=\frac{u_k^j+\Delta t{\tilde{N}}_k^j}{1-\Delta t\left[-\alpha +{(\frac{2\pi k}{L})}^2-{(\frac{2\pi k}{L})}^4\right]}$$

$$u_k^{j+1}=\frac{u_k^j+\Delta t{\tilde{N}}_k^j+\sqrt{\frac{2\epsilon N\Delta t}{h}}{\zeta }_k^j}{1-\Delta t\left[-\alpha +{(\frac{2\pi k}{L})}^2-{(\frac{2\pi k}{L})}^4\right]}$$

$$M_T(q)=\frac{1}{T}{\int }_0^T{1}_q(t)\text{dt}$$

$$1_q(t)=\{\begin{array}{ll}1& k=q\text{at time}t\\ 0& \text{otherwise}\end{array}$$

$$\underset{T\to \infty }{\lim }{M}_T(q)=P_{st}(q)$$

$$q_n=\frac{2\pi n}{Nh}$$

Full text

I Introduction

This work addresses the question of wavenumber selection in pattern forming systems. Pattern forming systems are characterized by the emergence of a band of spatially periodic steady states, as a certain quantity called the control parameter is varied Cross and Hohenberg (1993). Examples of pattern forming systems in physics are Rayleigh Bénard convection Cross and Hohenberg (1993) and directional solidification Grossmann et al. (1993). Although a large number of periodicities are mathematically allowed for such systems, experiments and simulations of realistic physical systems have repeatedly shown that only a narrow range of periodicities is realized in practice. The tendency of a system to prefer a narrow set of states out of many possible states is known as wavenumber selection.

Many mechanisms have been proposed to explain this phenomenon. These include evolution from random initial conditions having a power spectrum centered at a given wavenumber Cross and Greenside (2009); Schober et al. (1986), and control parameter ramps Cannell et al. (1983); Kramer et al. (1982); Hohenberg et al. (1985) (see Section III for a brief discussion). Various computational studies Kerszberg (1983a, b); Filho et al. (2005); Obeid et al. (2010) have also investigated the role of additive stochastic noise in wavenumber selection. The idea that noise can trigger wavenumber selection can be justified with the help of a simple dynamical system which evolves in a relaxational manner, i.e. it minimizes a potential energy. If the potential has several local minima, the deterministic system will evolve to one of these minima, depending on the initial conditions. However, in the presence of noise, the system will escape from any local minima and eventually reach the global minimum of the potential energy, where it will then spend most of its time, provided the noise strength is small. This has been investigated in Viñals et al. (1991) for the Swift-Hohenberg model Swift and Hohenberg (1977), which is a potential system. While this example of a potential system illustrates wavenumber selection, the systems that are described in Kerszberg (1983a, b); Filho et al. (2005); Obeid et al. (2010) do not have a potential energy function. The role of noise in inducing wavenumber selection in non-potential systems is not well understood.

In this paper, we further explore the role of noise as a mechanism for wavenumber selection in a model known as the stabilized Kuramoto-Sivashinsky (SKS) equation Misbah and Valance (1994). There are two main reasons for this choice. Firstly, we chose this model because it exhibits rich nonlinear behavior and a band of spatially periodic stationary states, while being relatively simple and one dimensional, which makes simulating it easy. The second reason is that it is non-variational, i.e. the deterministic driving force cannot be expressed as the gradient of a potential. Thus, studying this system would shed light on wavenumber selection in systems where minimization of a potential cannot explain selection of a particular state.

Noise induced wavenumber selection in this model has been studied before using direct numerical simulation Obeid et al. (2010), but only for small systems and a limited range of control parameters. More recently, Qiao et al. Qiao et al. (2016) used path integral methods to find the selected wavenumber for a range of control parameter values. It is of considerable interest to extend the work in Obeid et al. (2010) to larger system sizes, since true selection of a unique final state can only occur in the thermodynamic limit. Therefore, our aim in this work is to determine the noise selected wave number for large system sizes and for a wide range of control parameter values through direct integration and to determine if the results of Qiao et al. (2016) can be reproduced in this way.

This paper is organized as follows: Section II introduces the mathematical formalism used to study pattern forming systems and describes the SKS model. Section III describes in detail the results of Obeid et al. (2010) and Qiao et al. (2016). Section IV describes our computational method and Section V shows some of our results and their interpretation. Finally, Section VI discusses some of the drawbacks of our computational methods and touches on potential improvements that will be the focus of future work.

II Pattern Formation and the SKS Model

Pattern forming systems such as the ones mentioned above are represented by nonlinear partial differential equations that govern the evolution of a given physical quantity. A typical equation for such a system is of the form,

[TABLE]

where $u(x,t)$ is a field representing the quantity of interest, $\hat{L}_{p}$ is a linear differential operator acting on $u(x,t)$ and $\hat{N}[u(x,t)]$ is a nonlinear operator. The subscript $p$ on the linear operator indicates that it depends on the control parameter $p$. As an example, in the case of directional solidification, the quantity of interest $u(x,t)$ is the position of the interface between the liquid and solid phases. The trivial solution (or base state) $u_{b}(x,t)=0$ is generally a stationary state of these equations. However, this solution is stable only for a certain range of values of $p$. After $p$ crosses a critical value $p_{c}$, the trivial, spatially uniform state becomes unstable to periodic perturbations and a band of stable, periodic steady states emerges. To determine when the uniform state becomes unstable, we imagine adding to the uniform base state $u_{b}(x,t)=0$ the perturbation $\delta u\sim e^{iqx+\sigma(q)t}$. This perturbation is periodic in space with wavenumber $q$ and grows with time at a rate $\sigma$. We then substitute $u(x,t)=u_{b}(x,t)+\delta u$ in Eq. (1), retain only terms linear in $\delta u$ and derive an expression for the growth rate $\sigma$ as a function of $q$. If, for a given value of $p$, $Re(\sigma(q))$ is negative, that means that the perturbations with those wavenumbers decay exponentially with time and hence the uniform base state $u_{b}$ is stable to those perturbations. However, if $Re(\sigma(q))$ is positive, then perturbations with wavenumber $q$ grow exponentially with time, implying that the base state is unstable to them. Of course, in practice, the unbounded exponential growth of these perturbations is balanced by the nonlinear terms in the original equation.

To make these ideas more concrete, we illustrate the above steps for the deterministic SKS model. The SKS equation is given by,

[TABLE]

Here, $\alpha$ plays the role of the control parameter, and $u(x,t)$ is a dimensionless field of dimensionless space-time variables. This equation is used to describe directional solidification Misbah and Valance (1994) and the Burton-Cabrera-Frank model of terrace growth Bena et al. (1993). The trivial solution $u_{b}(x,t)=0$ is one of the steady states of this equation. We now add to this solution a small perturbation of the form $\delta u\sim e^{iqx+\sigma(q)t}$. (We restrict ourselves to the case where $\sigma$ is purely real.) Substituting into Eq. (2) and linearizing about the base state $u_{b}(x,t)=0$ state gives us an expression for the growth rate,

[TABLE]

From this, we see that the growth rate is non-negative for $1/2-\sqrt{1/4-\alpha}\leq q^{2}\leq 1/2+\sqrt{1/4-\alpha}$ and $\alpha\leq 1/4$, as shown in Fig. 1. For these modes, the state $u_{b}(x,t)=0$ is unstable. By solving the equation $\frac{\partial\sigma}{\partial q}=0$, we see that the growth rate is maximum for $q=q_{c}=1/\sqrt{2}$, for all $\alpha$. In summary, for $\alpha\leq 1/4=\alpha_{c}$, a band of periodic steady states exists, with wavenumbers given by,

[TABLE]

These periodic states may themselves be unstable to long wavelength periodic perturbations. Hence, in practice, one observes a band of periodic states that is narrower than suggested by Eq. (4) and is called the Eckhaus stable band Cross and Greenside (2009). Wavenumbers within this band are stable to long wavelength periodic perturbations.

III Previous Studies of Wavenumber Selection

As mentioned in the introduction, many possible explanations for wavenumber selection have been suggested. Schober et al. Schober et al. (1986) studied wavenumber selection in the context of the one dimensional Swift-Hohenberg equation Swift and Hohenberg (1977). The deterministic equation was integrated numerically, starting from random initial conditions with a power spectrum centered at a particular wavenumber $\bar{k}$. It was observed that the width of the noise averaged power spectrum (or the structure function) decreased with time until a sharp power spectrum centered about a final wavenumber $k_{\infty}$ was obtained. Their simulations showed that the final wavenumber depended on the value of $\bar{k}$.

Another proposed mechanism is the slow spatial variation of the control parameter, so that it changes from a value below threshold in one part of the experimental apparatus to a value above threshold in another part of the apparatus. This has been investigated experimentally Cannell et al. (1983) for the case of rotating Couette-Taylor flow. The experimental apparatus consists of two coaxial cylinders, separated by a gap that contains fluid. The inner cylinder is allowed to rotate with a given angular velocity. The control parameter, which is a function of the gap, is varied gradually by reducing the gap linearly towards the bottom of the apparatus. The variation in the gap is characterized by the angle made by outer wall with the vertical. The wavelength of the resulting flow pattern was measured in the straight section for different bulk values of the control parameter and was observed to lie in a narrow band, whose width decreased as the control parameter was increased beyond threshold. The observed wavelength was also found to be a periodic function of the aspect ratio (length of straight section/gap). The role of ramps has also been studied analytically Kramer et al. (1982); Hohenberg et al. (1985). Ramps provide an efficient way to obtain a desired, well defined periodic state. However, it must be noted that ref. Hohenberg et al. (1985) shows that varying different quantities which influence the control parameter leads to different selected wavenumbers. This implies that the selected wavenumber depends either on the initial conditions, or in case of control parameter ramps, on the chosen ramping protocol. In other words, the selected wavenumber in a deterministic system is not an intrinsic property of the system itself, but depends on the measuring protocol.

The role of noise in inducing wavenumber selection is the subject of some debate. While thermal noise is small enough that its effect on macroscopic patterns can be neglected Cross and Hohenberg (1993), the same cannot be said about noise in experimental apparatus Kerszberg (1983a, b); Filho et al. (2005); Obeid et al. (2010); García-Ojalvo et al. (1993); Parrondo et al. (1996). As argued in Kerszberg (1983b), it is hard to precisely control experimental conditions in situations such as directional solidification, leading to uncertainties which are most easily modeled as additive noise. Studies have also shown that multiplicative noise can have non-trivial consequences on the dynamics of bifurcations, such as shifting the primary bifurcation point in the Swift-Hohenberg equation García-Ojalvo et al. (1993) and inducing new patterns which are absent in the deterministic system Parrondo et al. (1996). In addition, Kerszberg (1983a, b); Filho et al. (2005); Obeid et al. (2010) have shown evidence of additive noise induced wavenumber selection. In our view, the other mechanisms mentioned above are too deterministic because they only sample a restricted set of perturbations. In contrast, stochastic noise is the most unbiased mechanism possible because it perturbs all periodic states of the deterministic system equally and does not favor a particular wavenumber over the others. It also samples a much larger set of perturbations. Hence, it is our view that more work needs to be done to gain a deeper understanding of the effects of noise, whether additive or multiplicative.

For the particular case of the SKS equation, Obeid et al. Obeid et al. (2010) and Qiao et al. Qiao et al. (2016) have recently investigated noise induced wavenumber selection. We present here a brief summary of their results, which are relevant to this work. Obeid et al. carried out direct numerical simulations of the noisy SKS equation,

[TABLE]

where $\zeta(x,t)$ is an additive Gaussian noise satisfying,

[TABLE]

and

[TABLE]

Here $\varepsilon$ is the noise strength. Eq. (5) was discretized using a simple finite difference scheme and the time integration was performed using the explicit forward Euler method. Various values of $\alpha$ between 0.17 and 0.24 and noise strengths ranging from $10^{-5}$ to $10^{-3}$ were studied. For $\alpha=0.24$, it was found that the state with $q=0.6995$ was most stable, in the sense that the system could not be knocked out of this state for $10^{8}$ time steps, with noise strengths up to $\varepsilon=5\times 10^{-4}$.

For $\alpha$ values farther from the critical value of 0.25, no unique state could be identified as being most stable. For $\alpha=0.2$, the allowed wavenumbers are $0.589\leq q\leq 0.767$. It was found that states with $0.650\leq q\leq 0.712$ all remained stable up to $10^{8}$ time steps and noise strengths up to $\varepsilon=2\times 10^{-3}$. Similarly, for $\alpha=0.17$, states with $0.650\leq q\leq 0.699$ were found to remain stable at $10^{8}$ time steps and noise strengths up to $\varepsilon=2\times 10^{-3}$. Increasing the noise strength made it hard to discern any periodic states and it was concluded that very long computational times would be required to destroy the stability of the states. In summary, the simulations in Obeid et al. (2010) hinted that the noise does indeed select a narrrow band of wavenumbers over others. However, selection of a unique state could only be demonstrated for $\alpha=0.24$. A conclusive numerical demonstration of state selection therefore requires integration over large system sizes, long times and a wider range of control parameters.

A more recent study of state selection in the SKS equation was conducted by Qiao et al. Qiao et al. (2016). They used the least action principle of Freidlin-Wentzell theory Freidlin and Wentzell (2012) to calculate transition probabilities between pairs of periodic steady states for the SKS equation. For the following stochastic process defined on a spatial domain [0, L],

[TABLE]

the probability of a particular trajectory $\Phi(x,t)=\phi(x,t)$ defined over a time interval [0, T] is,

[TABLE]

where $\mathcal{N}$ is a normalization constant independent of $\phi(x,t)$, $\varepsilon$ is the noise strength defined by,

[TABLE]

and $S_{T}[\phi(x,t)]$ is the action, given by,

[TABLE]

Eq. (9) implies that the most probable trajectory connecting two states of the system is the one which minimizes the action. The most likely paths entering and leaving successive periodic states of the SKS equation were computed by finding the minimum action for transitions between those states (for example, $S^{*}_{q_{j}\rightarrow q_{j+1}}$ is the minimum action to go from a periodic state with wavenumber $q_{j}$ to one with wavenumber $q_{j+1}$, and $S^{*}_{q_{j+1}\rightarrow q_{j}}$ is the minimum action for the reverse transition). These values were then used to determine the net direction of transitions between two adjacent states. By implementing this procedure for all pairs of successive steady states, the wavenumber corresponding to the selected state was found. In order to minimize the actions $S_{q_{j}\rightarrow q_{j+1}}$ and $S_{q_{j+1}\rightarrow q_{j}}$, the saddle state between the states $q_{j}$ and $q_{j+1}$ was found numerically. The saddle state of the lowest order amplitude equation Cross and Hohenberg (1993); Kramer and Zimmermann (1985) was used as the initial guess for this procedure. Having found the saddle state, the time reversed deterministic paths between the saddle state and the states $q_{j}$ and $q_{j+1}$ were used as initial guesses to find the minimum action paths (see sec. III B of Qiao et al. (2016) for a detailed discussion).

IV Calculating the Empirical Probability Distribution of Final States

Although Qiao et al. have devised a way to calculate the selected wavenumber, we wish to know if the same results can be obtained by direct integration of the equation of motion. There are two hurdles that must be overcome in order to do so. The first is that farther from threshold, several neighboring states in the middle of the Eckhaus band become very stable to noise Obeid et al. (2010). In order to determine which one, if any, is the most stable, one has to induce transitions between these states, so that the state in which the system spends most of its time is the most stable state. Observing such transitions between highly stable states requires extremely long integration times, as noted in Obeid et al. (2010). The second hurdle is that true state selection would only occur in the thermodynamic limit, which necessitates simulation of large systems. These two hurdles suggest that one looks for a fast and efficient integration algorithm. Here, we use a semi-implicit, Fourier spectral integration method Chen and Shen (1998). Using a semi-implicit time integration scheme instead of explicit time integration allows one to use a significantly larger time step without compromising on accuracy, thus yielding a much higher speed of integration. At the same time, using a Fourier spectral method to approximate spatial derivatives gives much higher accuracy than finite difference methods Trefethen (2000). Thus, the use of a semi-implicit Fourier spectral method enables us to integrate Eq. (5) for long times and large system sizes. The general idea behind the semi-implicit Fourier method is as follows. Consider the following partial differential equation,

[TABLE]

where $\hat{L}$ is a linear differential operator and $\hat{N}$ is a nonlinear functional of $u$ and its spatial derivatives. If we discretize space into $N$ (not to be confused with $\hat{N}$) grid points with lattice spacing $h$ and denote the value of the field $u(x,t)$ at each grid point by $u_{i}(t)$, we obtain a system of ordinary differential equations,

[TABLE]

Taking the discrete Fourier transform of this equation gives,

[TABLE]

For the (deterministic) SKS equation,

[TABLE]

and

[TABLE]

The Fourier transform of the discretized field $u_{i}$ can be found numerically using the Fast Fourier Transform (FFT). If $\widetilde{u}_{k}$ denotes the $k$-th wavenumber component of the discrete Fourier transform of the field $u$, then the discrete Fourier transform of Eq. (15) is Trefethen (2000),

[TABLE]

Evaluating the nonlinear term $\widetilde{\hat{N}_{k}}[u(x,t),\dots]$ in Fourier space directly is computationally expensive, since the Fourier transform of a product of functions is a convolution involving $O(N^{2})$ computations. Therefore, we do the integration by evaluating the nonlinear term in position space and then transforming it to Fourier space. Putting all this together, we get,

[TABLE]

To solve the issue of small time step, we use a semi-implicit integration scheme. We integrate forward in time by treating the linear terms implicitly and the nonlinear term explicitly Chen and Shen (1998). Let $u_{k}^{j}$ be the value of $u_{k}$ at time $t_{j}$. Then we approximate the time derivative by $du_{k}^{j}/dt=\frac{u_{k}^{j+1}-u_{k}^{j}}{\Delta t}$ and Eq. (18) becomes,

[TABLE]

With this semi-implicit scheme, it is possible to increase the time step by a factor of about 50 compared to that in an explicit time integration scheme, without causing instabilities and without loss of accuracy. This results in a significant speed up of the algorithm. Finally, in order to incorporate noise, we let $\zeta_{k}^{j}$ denote the value of the noise term at time $t_{j}$ and Eq. (19) is modified to read,

[TABLE]

where $\langle\zeta_{k}^{j}\rangle=0$ and $\langle\zeta_{k}^{j}\zeta_{k^{\prime}}^{j^{\prime}}\rangle=\delta_{k,-k^{\prime}}\delta_{jj^{\prime}}$. This is simply the Ito formula for the solution of a stochastic differential equation Garcia-Ojalvo and Sancho (1999). The procedure for generating noise in Fourier space which satisfies equations Eqs. (6) and (7) is also given in Garcia-Ojalvo and Sancho (1999), Appendix B.

V Results

Our aim is to compute the empirical probability distribution for the allowed periodic states and determine if the distribution has a peak at a particular wavenumber. If such a peak is present, it would support the hypothesis that noise is a possible mechanism of wavenumber selection. To do this, we use noise strengths that are around 0.1 % of the amplitude of the field $u(x,t)$. This noise strength is quite large, but it allows the system to explore a large number of allowed states. We start the simulation with the system initially placed in one of the periodic states and then perturb it with noise. The noise causes the system to visit other steady states. We then count the number of times each periodic state is visited and compute the fraction of the time spent in each state, for long runs of $\sim 10^{8}$ time steps. The idea is that the system will spend the greatest fraction of time in the most stable or selected state. It is important to clarify what we mean by “visiting a given periodic state”. Since we are using relatively large noise, the state of the system at any given time has a broad power spectrum, with non-zero components for several wavenumbers and it is, strictly speaking, inaccurate to say that the system is in a given periodic state. Instead, we consider the system to be in the neighborhood of a periodic state, if the power spectrum has a peak at that periodicity, and if the Fourier component at this periodicity is at least twice as large as the other Fourier components. Hence, what we actually calculate is the empirical probability for being near a periodic state. To do so, we calculate the time average of the indicator function for a state with wavenumber $q$,

[TABLE]

If the system is in a state with a power spectrum peaked at wavenumber $k$ (and this peak is at least twice as large as other peaks), then the indicator function for wavenumber $q$ is defined by,

[TABLE]

This quantity gives the fraction of time spent near the state with wavenumber $q$ and approaches the stationary probability distribution at very long times.

[TABLE]

Since wavenumber selection has already been demonstrated for a small system with control parameter $\alpha=0.24$ in Obeid et al. (2010), we attempted to carry out the procedure described above for $\alpha=0.22$. We restrict our simulations to the interval $0.17\leq\alpha\leq 0.24$, since for $\alpha<0.16$ the second harmonic in the stationary solution becomes active, leading to complicated instabilities Misbah and Valance (1994). For $\alpha=0.22$, the wave numbers that are stable to perturbations lie within the range $0.6136\leq q\leq 0.7486$, Misbah and Valance (1994). We used periodic boundary conditions, which implies that out of the wave numbers in the above range, only those appear in the simulation that satisfy,

[TABLE]

where $n$ is an integer $(1\leq n\leq N)$. Thus, the imposition of periodic boundary conditions limits the allowed wavenumbers to a discrete set given by Eq. (24), with $n$ representing the number of cells in the solution. For all our simulations, the initial condition was of the form $u_{in}\sim\text{sin}(2\pi n_{in}x/Nh)$ with $n_{in}$ being an integer. For the range $0.6136\leq q\leq 0.7486$, and $N=1024$, $h=0.5$, we have $50\leq n\leq 61$. Our time step was $\Delta t=0.3$.

We first used a small noise strength $\varepsilon=10^{-4}$ to determine which states were the most stable. We found that the states with $n=54,55,56$ were all very stable and the system could not be knocked out of them with this noise strength. So now the problem was to determine which one of these three states was the most stable. We turned up the noise strength gradually and found that at $\varepsilon=1.8\times 10^{-3}$, the $n=54$ and $n=56$ states became unstable, while the $n=55\ (q=0.6749)$ state persisted. When we plotted the fraction of time spent near each state, the resulting histogram developed a peak at this state. This peak persists until the end of the integration (about $10^{8}$ time steps), indicating that the system spends most of its time in that state. The histogram also gets narrower with time, again indicating that the system spends most of its time close to the $n=55$ state. We ran a couple of simulations at $\varepsilon=1.9\times 10^{-3}$ and $\varepsilon=3.9\times 10^{-3}$ and again observed that after spending some time around $n=54$ or $n=56$, the system transitioned to a neighborhood of the $n=55$ state and spent the most amount of time there (see Fig. 2 (a), left). In all three cases, the histogram became nearly stationary after around $10^{8}$ time steps. The observation that the most visited state remains the same in spite of increasing the noise strength is promising evidence of the selection of a unique state by noise. The most visited state was also found to be independent of initial condition.

We repeated these simulations for the same value of $\alpha$, i.e. $\alpha=0.22$, but different system sizes and noise strengths. In each case, we observed that after initially visiting a narrow band of states, the system settled down close to one of them. This state was independent of the noise strength and initial conditions. The largest system size we tried was $N=4000$ lattice points, which gives us the most precise estimate of the selected wavenumber (see Fig. 2(a), right). For this size, the difference between two successive wavenumbers is the smallest and is equal to $\Delta q=2\pi/Nh=0.003$. The most visited states for $\alpha=0.22$ are presented in Table. 1. Thus, we see that our best estimate for the selected wavenumber for $\alpha=0.22$ is $q_{s}=0.6754\pm\Delta q/2$ or $q_{s}=0.6754\pm 0.0015$, corresponding to the case $N=4000$.

Next, we repeated the same process for values of $\alpha$ farther from the critical value $\alpha_{c}=0.25$. As explained in Obeid et al. (2010), states in the middle of the Eckhaus band become more and more stable away from threshold. This made it necessary to use slightly larger noise strengths to destabilize the states and obtain the required probability distributions. For $\alpha=0.20$ and $N=1024$, the Eckhaus stable wavenumbers are $0.5890\leq q\leq 0.7608$ or $48\leq n\leq 62$. For a range of noise strengths between $3\times 10^{-3}$ and $6\times 10^{-3}$ and different initial conditions, we found $n=54\ (q=0.6627)$ to be the most visited state (Fig. 2(b), left). For $N=4000$ lattice points, we found $n=209$ to be the most visited state, see Fig. 2(b), right. Results for various system sizes are shown in Table 2. Finally, results for $\alpha=0.17$ are shown in Table 3 and Fig. 2(c).

Figs. 2(a), (b) and (c) show that the width of the probability distributions shrinks as the system size increases, indicating that the probability distributions approach Dirac delta functions peaked at a unique wavenumber in the thermodynamic ($N\rightarrow\infty$) limit. The probability distributions for each size also become narrower as the noise strength is decreased, but the maximum of the distributions remains at the same wavenumber.

For the sake of completeness, we have also reproduced Obeid’s result for $\alpha=0.24$ and extended it to larger sizes, as shown in Table 4.

The power spectra at the end of integration for $\alpha=0.22$ and low and high noise strengths are shown in Fig. 3. These figures show that the final power spectra have prominent maxima at the most visited wavenumber, in spite of the large noise strengths we have used. Small but visible second harmonic peaks are also observed, indicating that the system is close to a periodic steady state of the deterministic system. The power spectra for $\alpha=0.20$ and $0.17$ are similar. Note that we have plotted the number of cells $n$ on the x-axis, instead of the wavenumber.

We also show plots of the field $u(x,t)$ (for $\alpha=0.22$ and 1024 lattice points) at very long times in Fig. 4. The left plot is for $\varepsilon=0.0018$, while the right plot is for $\varepsilon=0.0039$. Finally, a plot of the selected wavenumbers for 4000 lattice points against $\alpha$ is shown in Fig. 5, with error bars representing the discretization error. It also shows the results obtained in Qiao et al. (2016) for the same values of $\alpha$. Both curves show that as the value of $\alpha$ is decreased below threshold, the selected wavenumber is shifted to the left of the critical wavenumber $q_{c}$ (horizontal line in Fig. 5). However, our values for the selected wavenumber do not agree with theirs far from $\alpha=\alpha_{c}$. We believe that the reason for the disagreement is likely numerical. Since the disagreement is particularly large away from threshold, it is our view that the use of the amplitude equation may be inappropriate. It is possible, for example, that the saddle solution for the amplitude equation is not a sufficiently accurate initial guess for the true saddle state for $\alpha$ values far from 0.25. Secondly, the time reversed deterministic path is the minimum action path only for potential systems Qiao et al. (2016). While it is reasonable to assume that it would be a good initial guess for non-potential systems, it is also possible that the actual minimum action path is far more complicated away from threshold, where the non-potential term in Eq. (5) is not small. The discretized action is a function of the discretized field at all lattice points and all times, and hence could be an extremely complicated function, with possibly multiple minima. It is not clear what kind of numerical uncertainties arise in minimizing this action, and whether the minimum action paths found in Qiao et al. (2016) are local minima or global minima. A thorough analysis of the optimization techniques is needed to explore these issues.

V.1 Extension to thermodynamic limit

In order to obtain an estimate for the selected wavenumber for the thermodynamic limit (in which case a continuous band of wavenumbers is allowed), we used cubic spline interpolation to get a smooth function passing through the discrete points obtained from our simulations. We implemented this procedure for $\alpha=0.22,0.20$ and $0.17$ and $N=4000$ lattice points. The continuous curves thus obtained show which wavenumber would be most probable in the continuum case. The results are shown in Figs. 6 through 8. Our continuum estimates for the selected wavenumbers are given in Table 5.

VI Conclusions

We have performed a detailed numerical study of noise induced wavenumber selection in the one dimensional SKS equation. We have shown that even in the presence of large noise, long time power spectra of the field are most likely to be peaked at a unique wavenumber that does not depend on the initial state or the noise strength. We note again that due to the large noise strengths we had to use, there were non-zero Fourier components at other wavenumbers, which means that for finite system sizes and large noise strengths, there is no true selection, as expected. However, since the position of the peak of the power spectrum is independent of noise strength and since the probability distributions become sharper with increasing system size and decreasing noise strength, it is plausible to conclude that the same wavenumber would be selected for noise strengths much less than the ones used here, except that it would require prohibitively large integration times.

It is interesting to compare our findings with other selection mechanisms proposed in Cannell et al. (1983); Kramer et al. (1982); Hohenberg et al. (1985). As discussed there, control parameter ramps are an efficient way to select a unique, well defined wavenumber. However, Hohenberg et al. (1985) also shows that varying different quantities that lead to the same final control parameter results in different selected wavenumbers. Ref. Schober et al. (1986) has shown that starting from random initial conditions and simulating deterministically, the final wavenumber depends on the initial state. In contrast, our work shows that at long times, the system is most likely to be found in a state with a dominant wavenumber, and this dominant wavenumber does not depend on initial conditions or noise strength, but is an intrinsic property of the SKS system. We have been able to show this for a wide range of control parameters and large system sizes, in contrast to previous work. Our work shows that noise induced wavenumber selection can occur in non-potential systems. Much remains to be done however; a future direction could be to devise efficient importance sampling techniques which would enable us to study transitions between states without using large noise strengths. It would also be interesting to see if it is possible to predict the selected wavenumber analytically, although to our knowledge no such analytical theory exists.