Controlling roughening processes in the stochastic Kuramoto–Sivashinsky equation

We present a novel control methodology to control the roughening processes of semilinear parabolic stochastic partial differential equations in one dimension, which we exemplify with the stochastic Kuramoto-Sivashinsky equation. The original equation is split into a linear stochastic and a nonlinear deterministic equation so that we can apply linear feedback control methods. Our control strategy is then based on two steps: first, stabilize the zero solution of the deterministic part and, second, control the roughness of the stochastic linear equation. We consider both periodic controls and point actuated ones, observing in all cases that the second moment of the solution evolves in time according to a power-law until it saturates at the desired controlled value.


Introduction
Roughening processes arise in nonequilibrium systems due to the presence of different mechanisms acting on multiple time and length scales and are typically characterized by a time-fluctuating ''rough'' interface whose dynamics are described in terms of a stochastic partial differential equation (SPDE). Examples are found in a broad range of different applications, including surface growth dynamics such as e.g. surface erosion by ion sputtering processes [1,2], film deposition in electrochemistry [3,4], or by other methods [5,6], fluid flow in porous media [7][8][9], fracture dynamics [10] and thin film dynamics [11][12][13][14][15], to name but a few. Not surprisingly, understanding the dynamics of the fluctuating interface in terms of its roughening properties, which often exhibit roughness) to evolve towards any desired value. By considering either periodic or point actuated controls, our results show that the second moment of the solution grows in time according to a power-law with a well-defined growth exponent until it saturates to the prescribed value we wish to achieve.
It is important to note that other control strategies have been proposed previously for controlling the surface roughness and other quantities of interest, such as the film porosity and film thickness in various linear dissipative models, including the stochastic heat equation, the linear sKS equation, and the Edwards-Wilkinson (EW) equation; see e.g. [5,6,[19][20][21][22][23][24][25]. However, it should also be emphasized that most of these works involve the use of nonlinear feedback controls which change the dynamics of the system and require knowledge of the nonlinearity at all times, something that may be difficult to achieve. We believe that our framework offers several distinct advantages since the controls we derive and use are linear functions of the solution which do not affect the overall dynamics of the system and also decrease the computational cost. Another recent study is Ref. [26] which considered a deterministic version of the KS equation, and presented a numerical study of the effects of the use of ion bombardment which varies periodically in time on the patterns induced by the ion beams on an amorphous material. In particular, this study found that rocking the material sample about an axis orthogonal to the surface normal and the incident ion beam, which corresponds to making the coefficients of the KS equation periodic in time, can lead to suppression of spatiotemporal chaos.
The work presented in this paper is motivated by earlier research carried out by our group: on one hand, the study of noise induced stabilization for the KS equation [27,28] and, on the other hand, the study of optimal and feedback control methodologies for the KS equation and related equations that are used in the modeling of falling liquid films [17,18,29]. It was shown in [27,28] that an appropriately chosen noise can be used in order to suppress linear instabilities in the KS equation, close to the instability threshold. Furthermore, it was shown in [17,18] that nontrivial steady states and unstable traveling wave solutions of the deterministic KS equation can be stabilized using appropriate optimal and feedback control methodologies. In addition, similar feedback control methodologies can be used in order to stabilize unstable solutions of related PDEs used in the modeling of falling liquid films, such the Benney and weighted-residuals equations.
The rest of the paper is structured as follows. Section 2 introduces the sKS equation and discusses means to characterize the roughening process of its solution. In Section 3 we outline the general linear control methodology which is applied to the case of periodic controls in Section 4, and point actuated controls in Section 5. A summary and conclusions are offered in Section 6.

The stochastic Kuramoto-Sivashinsky (sKS) equation
Consider the sKS equation: normalized to 2π domains (x ∈ [0, 2π ]) with ν = (2π /L) 2 > 0, where L is the size of the system, with periodic boundary conditions (PBCs) and initial condition u(x, 0) = φ(x). ξ (x, t) denotes Gaussian mean-zero spatiotemporal noise, which is taken to be white in time, and whose strength is controlled by the parameter σ : where G(x − x ′ ) represents its spatial correlation function. We can, in principle, consider the control problem for SPDEs of the form (1) driven by noise that is colored in both space and time. Such a noise can be described using a linear SPDE (Ornstein-Uhlenbeck process) [30]. The noise term can be expressed in terms of its Fourier components as: whereẆ k (t) is a Gaussian white noise in time and the coefficients q k are the eigenfunctions of the covariance operator of the noise.
For the noise to be real-valued, we require that the coefficients q k verify q −k = q k .
Proofs of existence and uniqueness of solutions to Eq. (1) can be found in [31,32], for example. The behavior of Eq. (1) as a function of the noise strength, and for particular choices of the coefficients {q k } has been analyzed in detail in [27,28]. In particular, it was shown that sKS solutions undergo several state transitions as the noise strength increases, including critical on-off intermittency and stabilized states.
The quadratic nonlinearity in Eq. (1) is typically referred to as a Burgers nonlinearity. We note that an alternative version of Eq. (1) is found by making the change of variable u = −h x , giving rise to The main effect of this transformation is to change the dynamics of the average u 0 (t) = 1 2π  2π 0 u(x, t) dx of the solution. Indeed, Eq. (1) with PBCs preserves the value of u 0 whereas as a consequence of the nonlinear term (h x ) 2 , Eq. (4) does not conserve the mass h 0 (t) = 1 2π  2π 0 h(x, t) dx. Both equations have received a lot of attention over the last decades, with Eq. (1) more appropriate in mass-conserved systems such as the dynamics of thin liquid films [28,[11][12][13][14][15], and Eq. (4) relevant in modeling surface growth processes such as surface erosion by ion sputtering processes [3,4,1,2,33,22,34]. It is also worth mentioning that the quadratic nonlinearity appearing in Eq. (4) is the same as that in the Kardar-Parisi-Zhang (KPZ) equation [35,36] In fact extensive work indicates that Eqs. (4) and (5) are asymptotically equivalent, something referred to as the ''Yakhot conjecture'' [37][38][39]. Throughout the remainder of this study we will refer to Eq. (1) as the sKS equation with Burgers nonlinearity and Eq. (4) as the sKS equation with KPZ nonlinearity.

Surface roughening
An important feature of systems involving dynamics of rough surfaces is that one often observes the emergence of scale invariance both in time and space, i.e., the statistical properties of quantities of interest are described in terms of algebraic functions of the form f (t) ∼ t β or g(x) ∼ x α , where α and β are referred to as scaling exponents. An example of this is the surface roughness, or variance of u(x, t), which is defined as We remark that u 0 may or may not depend on time, depending on whether we consider the Burgers or the KPZ nonlinearities. Usually the above quantity grows in time until it reaches a saturated regime, in which the fluctuations become statistically independent of time and are scale-invariant up to some typical length scale of the system, say ℓ s . This behavior can be expressed as: where ⟨. . .⟩ denotes average over different realizations, β is the so-called growth exponent [16], and t s and r s are the saturation time and saturated roughness value, respectively, which depend on the length scale ℓ s . In particular, at a given time t < t s , the correlation of these fluctuations is on a spatial length scale which grows in time as ℓ c ∼ t 1/z . Therefore, saturation occurs whenever ℓ c = ℓ s from which we find r s ∼ ℓ α s with α = βz. In this context, the exponents α and z are the roughness and dynamic exponents, respectively, and their particular values determine the type of universality class [40]. For example, it is known that the long-time behavior of the KPZ equation (5), is characterized by the KPZ universality class with α = 1/2 and z = 3/2, while its linear version, which is referred to as the EW equation, is characterized by the EW universality class with α = 1/2 and z = 2 [16,41,42,36].
Alternatively, the solution u(x, t) can also be written in terms of its Fourier representation whereû k (t) are the Fourier components. By making use of Parseval's identity, we can compute the expected value of r(t) 2 as follows: where we have defined the power spectral density S( we propose a control methodology precisely for this purpose.

Linear feedback control methodology
The methodology to control the roughness of the sKS solution consists of two main steps. First, using a standard trick from the theory of semilinear parabolic SPDEs, see e.g. [32], we define w to be the solution of the linear sKS equation: and write the full solution u of Eq. (1) as The important point here is to note that Eq. (11) is now a deterministic PDE with random coefficients and so we are in a position where we can apply the methodology for nonlinear deterministic PDEs we have developed in previous works [17,18], to stabilize its zero solution-something possible as long as w and its first derivative are bounded in an appropriate sense (see Section 4.2 for a justification of this point). We therefore introduce the controlled equation for v: where ) is the number of controls, and b det n (x) are the control actuator functions. Here we use [x] to denote the integer part of x.
Once the zero solution of the equation for v has been stabilized, the second step is to control the roughness of the solution by applying appropriate controls to the linear SPDE (10) for w so that the solution is driven towards the desired surface roughness r d . In the following, we apply this methodology to the sKS equation, Eq. (1) or (4), by choosing two different types of controls, namely periodic controls, when the controls are applied throughout the whole domain and point actuated ones, when the control force is applied in a finite number of positions in the domain.

Derivation of the controlled equation
and take the inner product with the basis functions e ikx to obtaiṅ with k ∈ Z and a dot denoting a time derivative. We have intro- and note that g k are functions of the coefficients of v and w.
Next we define the following vectors and matrices. We denote are the coefficients of the (slow) unstable modes, and z v With these definitions we rewrite the infinite system of ordinary differential equations (ODEs) (14) aṡ The key point now is to note that if there exists a matrix K det such that all the eigenvalues of the matrix A u + B det u K det have negative real part, then the controls given by f det where K det n is the nth row of K det , stabilize the zero solution of Eq. (12) (see [17,18] for previously derived methodologies for deterministic systems). The proof of this follows the same type of Lyapunov argument as for the deterministic KS equation and is justified as long as we have nice bounds on w, something we will demonstrate below. It should be emphasized that in Eq. (15) we have suppressed the influence of the nonlinearity on the SPDE without assuming knowledge of its value at all times and without changing the fundamental dynamics, in contrast to previous work [5,6].
The next step is to control the stochastic linear equation for w such that the value of the second moment evolves towards a desired target. To this end we write and take the inner product with the basis functions to obtain the following infinite system of ODEs for the Fourier coefficientṡ Here The solution to system (18) iŝ and it easily follows that We observe that in this case the expected surface roughness depends only on the eigenvalues of the linear operator L = −ν∂ 4 x −∂ 2 x ; these can be controlled using feedback control to direct the evolution towards the desired value of surface roughness r d .
Hence we introduce the controlled equation for w, (21) where m 2 = 2l 2 is the number of controls (l 2 needs to be larger than or equal to the number of unstable modes and will be specified later), and we choose the functions b rand We also notice that we do not need to control the eigenvalue corresponding to the constant eigenfunction (k = 0), since it does not contribute to the surface roughness. By truncating the system into N modes (with N sufficiently large so that the contribution from higher modes can be neglected) and taking inner products with the basis functions, we arrive aṫ Remark 1. An important point to note is that because of the choice of periodic functions for b rand n , the system (22) is decoupled. In fact, with such a choice of actuator functions, the matrix B rand u is the identity matrix, and B rand s± are zero matrices. As will be shown in Section 5, this is not the case for point actuated controls.
The surface roughness for m 2 = 2l 2 controls is therefore given  .
If we denote the desired surface roughness as r 2 where we have used the fact that the coefficients q k are real with q −k = q k (see Eq. (3)). The chosen eigenvalues for the controlled modes are λ k , and we take them to be λ k = λ for all k to arrive at To control the surface roughness we therefore define the controls f rand k such that the new eigenvalues satisfy the following relation Finally, putting Eqs. (12) and (21) together, yields the controlled equation for the full solution u

Proof of applicability of the control methodology
Our aim here is to prove that the solution v can indeed be controlled to zero even though Eq. (12) has random coefficients, i.e. the terms (vw) x and ww x . We will show that by adopting a similar argument as used for the proof of existence and uniqueness of solutions of the sKS equation [32], we can apply a Lyapunov-type argument as in the deterministic KS equation [18]. We use (24) to write the solution of Eq. (21) as x and F is an operator discretized as We take G to be a trace class operator, so that it satisfies [32, Writing where we used ⟨e k , e j ⟩ = 0 and f k = λ + νk 4 − k 2 . Since we are assuming that the covariance matrix G is such that assumption (3.1) in [32] is satisfied, we have that w(t) ∈L 2 (0, 2π ), the space of mean zero L 2 functions, almost surely, for any time t. This also means that there exists a continuous version of w [30] that we shall consider from now on. Now we define B(u, v) = uv x and b(u, v, w) = ⟨B(u, v), w⟩ =  2π 0 uv x w dx, which satisfy the following relations [32,43]: and [32, Proposition (2.1)]: On the other hand, we notice that the existence of the matrix K det implies that the operator A, such that for some positive constant a, which in turn depends on the eigenvalues we choose for the controlled operator. Therefore, multiplying Eq. (12) by v and integrating by parts yield where we have used Young's inequality and relations (27) and (28). The term c∥Av∥ 2 L 2 can be controlled using sufficiently strong controls and the last term on the right-hand-side is a constant that depends on the desired surface roughness and which can be controlled by choosing large enough eigenvalues. Therefore, by choosing the controls such that a is large enough, ∥v∥ 2 L 2 is a Lyapunov function for this system and the zero solution for the controlled equation for v is stable.

Numerical results
We apply now the methodology outlined above with periodic controls to the sKS with either the Burgers nonlinearity (Eq. (1)) or the KPZ nonlinearity (Eq. (4)). For simplicity, we consider white noise in both space and time (q k = 1). All our numerical experiments are performed using spectral methods in space and a second-order backward differentiation formula scheme in time [44].

Controlling the roughening process
We solved Eqs. (1) and (4) for ν = 0.05 and σ = 0.5, controlling its solutions towards various desired surface roughnesses r d . The results are presented in Fig. 1. We observe that in both cases the solution exhibits a power-law behavior at short times of the form given by Eq. (7) until the solution saturates to the desired value of the roughness. It is interesting to note that the exponent in all cases is the same and with the value β ≈ 0.43, independently of the type of nonlinearity and desired surface roughness (note that the exponent in Fig. 1 is ≈ 0.85 = 2β, since we are plotting ⟨r(t) 2 ⟩). This becomes even clearer if time and surface roughness where x = t/r 1/β d . Fig. 2 shows that all the different cases presented in Fig. 1 collapse into a single curve which is given by (31) with the universal value β ≈ 0.43.
We also study the effect of changing the domain by varying the parameter ν. Fig. 3 shows the numerical results obtained when we fix the target value r d and change the parameter ν. We observe that changing the domain does not change the growth rate (yielding the same growth exponent β ≈ 0.43) but it does slightly affect the final value of the roughness. An important point to note is that since we are controlling the surface roughness of the solution r(t) to be at a specified value r d , the saturated state in which the statistical properties become stationary, is reached whenever r(t) = r d . Therefore, the saturation time, and the long-time roughness value, should not depend on the system size. Figs. 4 and 5 show typical snapshots of the spatiotemporal evolution of the controls (left panels) and the time evolution of their L 2 -norm (right panels) for the controlled sKS solution with Burgers nonlinearity with r 2 d = 20 and two values of the system size, namely ν = 0.03 and ν = 0.05 (see Fig. 3 for the corresponding surface roughness). We observe that neither the norm of the controls nor their amplitudes increase as a function of the system size.

Changing the shape of the solution
It is important to emphasize that in addition to controlling the roughness of the solution of the sKS equation, we can also change its shape, something that could have ramifications in technological applications such as materials processing. We quantify this by considering the surface roughness of the solution to be its distance to the desired state. Ifū(x) is the ultimate desired shape of the solution, then the quantity we are trying to control now becomes Using Parseval's identity we compute the expected value of r(t) 2 To control the shape of the solution, we can therefore control the solution of Eq. (12) for v to the desired shape rather than controlling it to zero. This in turn implies the use of f det We use the steady states of the KS equation for a chosen value of ν to define the desired shapē u. Results are shown in Fig. 6 for ν = 0.5, where we can see that the solution is fluctuating around the imposed shape.

Point actuated controls
We now consider controls that are point actuated and not x , we obtain the following infinite system of linear stochastic ODEṡ where the coefficients b k n are defined from the functions b(x) = δ(x − x n ) as before, b k n =  2π 0 b n (x)e ikx dx. We can see that the difference between the above system and the periodic controls one given by (22), is that now the system is coupled. In fact the coupling matrix is not symmetric, and most importantly, it does not commute with its transpose. Therefore the solution does not follow directly and we cannot easily write the second moment of the coefficients as a function of the eigenvalues as in the previous section. To obtain the controlled equation we thus need to apply a different approach.
Let the controls F = [f 1 , . . . , f m ] be such that F = Kŵ whereŵ is a vector containing the Fourier coefficients of w, and the matrix K is to be determined. Since the equations are not decoupled we cannot multiply by w and integrate to find directly the second moment of the coefficients. However, we can make use of results derived in [45] which provide simplified formulas for the first and second moments of systems analogous to (34). Let Ξ be the vector Ξ k = ξ k and C = A + BK where A = diag −νk 4 + k 2 and B kn = b k n , so that we can write the truncated system (22) aṡ w = Aŵ + BKŵ + Ξ := Cŵ + Ξ .
We also assume without loss of generality that m(0) = E(ŵ(0)) = 0 and P(0) = E(ŵ(0)ŵ(0) T ) = 0. Then Theorem 4 in [45] states  I is an appropriately sized identity matrix and the zeros stand for zero matrices of appropriate size. We compute e Mt and conclude that and Remark 2. In the periodic case, the matrix C is diagonal, so CC T = C T C , and this is exactly the same as In addition, when choosing the eigenvalues of C , we can ensure that it is invertible and therefore C + C T is also invertible, which gives so as t → ∞, P(t) → −σ 2  C + C T  −1 and where λ k are the chosen eigenvalues of C . Hence we recover the same result as before.
It is important to note that the matrix C is not normal, i.e. it does not commute with its transpose, and the eigenvalues of C + C T do not satisfy the useful properties that allow us to obtain (30). However, we are not interested in knowing the full matrix P(t), but only its trace tr(P(t)) = tr (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) By now making use of the linearity of the trace and its continuity to pass it inside the infinite sum, we get We also note that Similarly we can prove that the terms multiplied by t n n! are of the form tr  (C + C T ) n−1  and we finally obtain tr(P(t)) = σ 2  tr(I)t + tr(C + C T ) We proceed by assuming that C + C T is invertible, so that we can multiply by I = (C +C T ) −1 (C +C T ) and add and subtract pertinent terms to obtain tr(P(t)) = −σ 2 tr Remark 3. This does not change the proof provided in Section 4.2, it only changes the formula for the covariance so that the bounds are still valid.

Remark 4.
We emphasize that the following assumptions were made here: (a) C + C T needs to be invertible.
(b) In order for the surface roughness to converge to a finite value, we require all of the eigenvalues of C + C T to be negative, so that the exponential part disappears.

Computation of the matrix K
In Eq. (42) we need to control the trace of D −1 and we can do that by prescribing the eigenvalues of D. Hence we can control the surface roughness by finding a matrix K such that the eigenvalues of are a given set {µ 1 , . . . , µ N }. Since we only wish to prescribe the eigenvalues of D, rather than knowing all of its entries, we can tackle this problem by using the information provided by the characteristic polynomial, χ D , of D. We know that where J is a subset of {1, . . . , N}. Equivalently we can express χ D in terms of the sum over all its diagonal minors, i.e., where η k is the sum over all of the diagonal minors of size k of D.
This translates into a system of N nonlinear algebraic equations, for the m × N entries of the matrix K -see [46] for details on the solution to this problem. For the purposes of our study, we will make use of a nonlinear solver (e.g., Matlab's fsolve) to obtain the matrix K . Given the structure of the matrix B and the fact that the system is underdetermined, convergence is rather slow when solving the problem directly. We overcome this by performing a change of variables: we obtain the SVD decomposition of B by finding matrices X and Y such thatB = XBY T , and multiply Eq. (43) by X T on the left and by X on the right. We then defineK = Y T KX , A = X T AX andD = X T DX , so that we obtain the equatioñ This is of the same form as (43), but where the matrixB is diagonal. We find that this accelerates the convergence of the system (for moderate values of N) and we were able to get satisfactory numerical results, which we now present.

Numerical results
We apply the methodology presented in the previous section with point actuated controls to the sKS with the Burgers nonlinearity (Eq. (1)) (similar results are expected for the KPZ nonlinearity (Eq. (4))). We solved Eq. (1) for ν = 0.4 and σ = 0.5. For this value of ν the linear operator has 3 unstable modes and we apply m = 3 controls. We note that even though we do not need to control the mode corresponding to the first moment of the solution when using periodic controls, we benefit from doing so in this case, since the matrix D would not be invertible if we allowed for a zero eigenvalue. We consider either space-time white noise (q k = 1) or colored noise with q k = |k| −1 , (which is chosen to decay at a fast rate so that the system can be truncated at a smaller value of N) and control the solution towards various desired values r d of the surface roughness.
The results are depicted in Fig. 7 where we observe that the solution still exhibits a power-law behavior with similar exponent as in the periodic case (there we found β ≈ 0.43) until it saturates at the desired value of the surface roughness. It is noteworthy that even though we obtained satisfactory results for the range of values of r 2 d selected in Fig. 7, further increase of r d does not lead to the expected saturated results. This may be due to the relatively small system truncation value N = 21 that was found necessary in order to obtain convergence. Further work is required in this direction but this is beyond the scope of the present study.

Conclusions
We have presented a generic methodology for controlling the surface roughness of nonlinear SPDEs exemplified by the sKS equation with either the Burgers nonlinearity or the KPZ nonlinearity and using periodic or point actuated controls.
We have shown that with the appropriate choice of periodic controls the solutions of these equations can be forced to have a wide range of prescribed surface roughness values, defined to be the distance of the solution to its mean value. We are also able to force the solutions to a prescribed shape e.g. steady state solutions of the deterministic KS equation. We find that the solution to the controlled problem exhibits a power-law behavior with a universal exponent β ≈ 0.43, which is not affected by changes in the length of the domain, and is found to be independent of the type of nonlinearity of the sKS equation.
When using point actuated feedback controls, the problem becomes considerably harder to solve due to the fact that the resulting system of linear ODEs is not decoupled. This leads to the need to solve a new matrix problem which is similar to a matrix Lyapunov equation; to the best of our knowledge such a problem has not been tackled before. The complexity of this problem makes it harder to solve for a large system truncation value N, but we have obtained satisfactory results when controlling towards a range of surface roughness values for moderate N. This is an interesting separate problem and our detailed results and associated algorithm for its solution can be found in [46].
We believe that our framework offers several distinct advantages over other approaches. First, the controls we derived are Fig. 7. Squared value of the surface roughness of the solutions to the sKS equation with Burgers nonlinearity for ν = 0.04, σ = 0.5 and different values of the desired surface roughness, ranging from 2 to 6. Left: using space-time white noise; Right: using colored noise described by the coefficients q k = |k| −1 . We applied m = 3 point actuated controls, which were located at the positions x 1 = π 3 , x 2 = π, x 3 = 5π 3 .
linear functions of the solution u, and this in turn decreases the computational cost of their determination. Second, our splitting methodology allows us to deal with the nonlinear term directly rather than including it in the controls, thus rendering the resulting equation essentially linear and easier to handle. One interesting observation is that feedback control methodologies can be used, in principle, in order to accelerate the convergence of infinite dimensional stochastic systems such as the sKS and the KPZ equations to their steady state. This might prove to be a useful computational tool when analyzing the equilibrium properties of such systems, e.g. calculating critical exponents, studying their universality class, etc. Accelerating convergence to equilibrium and reducing variance by adding appropriate controls that modify the dynamics while preserving the equilibrium states has already been explored for Langevin-type samplers that are used in molecular dynamics [47,48]. In addition, it would be interesting to investigate how our methodology could be used to control the kinetic roughening process of the system. In particular, our results show that the dynamics towards saturation is described in terms of power laws. Whether we can control the values of the associated scaling exponents during such scale-free behavior is something that requires a systematic study of different stochastic models by controlling them to evolve towards large values of the surface roughness. We shall examine these and related issues for the sKS and KPZ equations in future studies.