Optimal control of thin liquid films and transverse mode effects

We consider the control of a three-dimensional thin liquid film on a flat substrate, inclined at a non-zero angle to the horizontal. Controls are applied via same-fluid blowing and suction through the substrate surface. We consider both overlying and hanging films, where the liquid lies above or below the substrate, respectively. We study the weakly nonlinear evolution of the system, which is governed by a forced Kuramoto--Sivashinsky equation in two space dimensions. The uncontrolled problem exhibits three ranges of dynamics depending on the incline of the substrate: stable flat film solution, bounded chaotic dynamics, or unbounded exponential growth of unstable transverse modes. We proceed with the assumption that we may actuate at every point on the substrate. The main focus is the optimal control problem, which we first study in the special case that the forcing may only vary in the spanwise direction. The structure of the Kuramoto--Sivashinsky equation allows the explicit construction of optimal controls in this case using the classical theory of linear quadratic regulators. Such controls are employed to prevent the exponential growth of transverse waves in the case of a hanging film, revealing complex dynamics for the streamwise and mixed modes. We then consider the optimal control problem in full generality, and prove the existence of an optimal control. For numerical simulations, we employ an iterative gradient descent algorithm. In the final section, we consider the effects of transverse mode forcing on the chaotic dynamics present in the streamwise and mixed modes for the case of a vertical film flow. Coupling through nonlinearity allows us to reduce the average energy in solutions without directly forcing the linearly unstable dominant modes.

1. Introduction. The dynamics of thin liquid films in various physical situations form a range of classical and important problems in fluid mechanics, attracting a lot of attention from researchers in the field (see [33,17] and the references therein). In this work, we study the optimal control problem for a gravity-driven, thin, viscous liquid film on an inclined flat substrate, where controls are applied at the substrate surface by means of same-fluid blowing and suction. We allow the fluid to be either overlying or hanging, where the film lies above or below the substrate, respectively. Falling liquid films have a wide range of industrial applications, including coating processes [68,79] and heat and mass transfer [20,76,45,48,74].
For coating processes, a stable film with a flat interface is desirable, whereas heat (and mass) transfer is improved by the presence of interfacial deformations; this is due to increased surface area, the presence of recirculation regions in the wave crests, and effective film thinning where the heat transfer becomes almost purely conductive. Thus, it is desirable to know how to control a liquid film to have a desired predetermined interface.
The experimental work of Kapitza and Kapitza [36] displayed the vast range of dynamical behaviors that could be obtained with flows of this kind and promoted analytical and numerical studies of the problem. In the absence of controls, the system possesses an exact flat film solution [53] (known as the Nusselt solution) with a semiparabolic velocity profile in the streamwise direction. For an overlying film, Yih [91,92] and Benjamin [6] considered the linear stability of the exact Nusselt solution and showed that the film is unstable to long waves beyond a critical Reynolds number which depends on the angle of inclination; for vertical arrangements, the critical Reynolds number is zero. Starting with the full Navier-Stokes equations, families of reduced-order models may be constructed to simplify the problem both analytically and numerically, with the aim of capturing the evolution of the film interface. With a long-wave assumption and formal asymptotics, a Benney equation [7,23] for the interface height may be constructed for Reynolds numbers close to critical; such equations usually retain the effects of gravity, viscosity, and surface tension. These highly nonlinear models have been studied extensively and generalized by a number of authors [41,67,49,5]; however, they lack global existence of solutions, and finite-time blow-ups have been observed in numerical simulations [64,30,66]. Although including many physical mechanisms, Benney equations are not effective at modeling thin film flows at higher Reynolds numbers. Coupled systems of evolution equations for the interface height and fluid flux can be derived as an alternative-notably the integral boundary layer formulation of Kapitza [35] and Shkadov [75] and the weighted residual models [69,70,72].
The lowest rung in the hierarchy of models is occupied by the weakly nonlinear evolution equations; for the thin film flow problem under consideration, these models are relatives of the Kuramoto-Sivashinsky equation (KSE), which governs the dynamics of small perturbations to the flat interface solution. The classical KSE in one space dimension is written (1) η t + ηη x + η xx + η xxxx = 0 and is usually studied with periodic boundary conditions on the interval [0, L]. This partial differential equation (PDE) exhibits a range of dynamical behaviors in windows of the bifurcation parameter L, including steady and traveling waves, time-periodic and quasi-periodic solutions, and full spatiotemporal chaos [37,56,78]. It has been observed that (1) possesses a finite-dimensional global attractor [26], and numerical results suggest a finite energy density in the limit as the periodicity L becomes large; the stronger result that the L ∞ -norm of solutions at large times is bounded uniformly as L → ∞ was proved for the steady problem [47] and appears to hold in general. A number of authors have considered this question analytically [52,26,16,29,12,22,54,24], yet the best result is far from agreeing with the optimal numerical bound. An equipartition of energy was observed for solutions [62,90,83], with a flat range in the power spectrum for the long wavelength modes, rising to a peak at the most active mode, and then decaying exponentially for the high frequencies due to strong dissipation on small scales. Additionally, Collet et al. [15] showed that solutions become instantly analytic, as is observed numerically in the decay of the Fourier coefficients as the wavenumber becomes large. A KSE may be obtained in the thin film flow context from Benney models by seeking the evolution of a small perturbation to the flat film solution. Weakly nonlinear models have been derived for the three-dimensional (3D) film problem under consideration, yielding evolution equations in two spatial dimensions [50,51,86,19] similar to (1); such an equation is the main focus of this work. Variants of (1) arise in core-annular flows [55] or fiber coating problems [18].
Modeling thin liquid films with the 2D Navier-Stokes equations yields evolution equations that overlook transverse and mixed-mode phenomena which in some situations dominate the interface dynamics. For the overlying film problem, an initially 2D flow is seen to transition to 3D waves in both experiments and numerical simulations [44,58,39,38]. Hanging film flows (often referred to as film flows on inverted substrates) are less well understood, with few experimental studies performed. In this case, the linear theory from the 3D formulation predicts a Rayleigh-Taylor (RT) instability in the transverse modes, independent of the Reynolds number. This instability leads to the formation of rivulet structures which were observed in experiments by Markides and Charogiannis [46]. They also observed complicated pulse dynamics in the streamwise direction on the crests of the rivulets, although their parameters were not so extreme that dripping took place. In [13], the authors considered Stokes flow on an inverted substrate, deriving a Benney equation for the streamwise interface evolution (no transverse dynamics), and performed experiments. They found that the fluid parameters for which dripping occurred coincided with the parameters for which their Benney model exhibited absolute instability where small perturbations grow locally (as opposed to convective instability, where perturbations are convected with the flow). It can be surmised from these results that the dripping of a hanging film can be broken down into two instabilities. First, the rivulet structures form due to a transverse RT instability. Then dripping occurs for some choices of parameters; this appears to coincide closely with the regime of absolute instability for the streamwise dynamics. In a related study, Lin, Kondic, and Filippov [42] considered the case of a nonwetted substrate with numerical simulations of a Benney equation and observed that fluid fronts were unstable to a transverse fingering instability. Thin rivulets form with approximately equal width in the spanwise direction, with fast moving "drop-like" waves appearing on the rivulets as observed in the wetted case. However, the authors do not attribute the fingering instability to be of RT type. Controls considered in the current work attempt to prevent dripping of liquid films by averting the initial rivulet formation at the weakly nonlinear level.
There are a number of ways to influence the interfacial dynamics of thin liquid films. Spatially varying topography may be utilized to create patterned steady states but does not have a considerable effect on the stability of solutions [63,21,88,27,61,73]. The openloop controls used in such studies are steady since the topography is fixed. The effect of an electric field which is normal or parallel to the substrate has also been studied extensively; the former has a destabilizing effect on the flow [89,85]. The film flow problem has also been considered with the addition of heating at the wall surface which is not necessarily uniform [71,32,77,31] and even in conjunction with substrate topography [10]. Other possibilities for influencing the dynamics of thin films include the introduction of magnetic fields [4], surfactants [11], and substrate microstructure or coatings to induce effective slip [34] the current work, we focus on the control of the film interface using same-fluid blowing and suction controls at the substrate surface. For the 2D simplification of this problem with an interface in one spatial dimension, the weakly nonlinear evolution can be modeled by the 1D KSE (1) with the introduction of some nonzero right-hand side ζ(x, t) (this is not the case for the majority of the above control methodologies). The optimal control problem for this equation (in fact a generalization with dispersive and electric field effects) was considered by Gomes, Papageorgiou, and Pavliotis [25]. For their numerical experiments, the authors consider the optimization of point-actuator locations (given an initial condition) to stabilize unstable traveling wave solutions. This is in contrast to the current work, where our controls are more regular in space and actuation is not restricted to a finite set of points. An alternative to same-fluid controls at the substrate surface is air-blowing and suction controls via actuators on a plate which is parallel to the substrate on which the fluid lies [57].
The main objective of the current paper is to present the theory and numerical experiments for the optimal control of the KSE in two space dimensions: Here η is a zero-average perturbation of the flat film solution, and κ > 0 (< 0) means that the film is overlying (hanging). The inertial and gravitational effects are manifested in the second derivative terms, and the bi-Laplacian term containing mixed derivatives corresponds to surface tension effects. The function ζ is a control which corresponds to blowing and suction through the substrate surface. In this paper, we consider only passive open-loop control systems; research on active control (closed-loop feedback control) problems for (2) is underway by the authors. The current paper is organized as follows: In section 2, we present the physical problem and discuss the hierarchy of models culminating with the 2D KSE (2). In section 3, we consider purely transverse controls, ζ =ζ(y, t). Seeking an optimal control of this form reduces to a linear quadratic regulator (LQR) problem for each transverse Fourier mode. The resulting optimal control is used to suppress the exponential growth of the purely transverse modes for a hanging film set-up. In section 4, we discuss both the analytical and numerical aspects of the full optimal control problem, where ζ ≡ ζ(x, y, t). Existence of an optimal control is proven, and a forward-backward sweeping method is employed for numerical experiments. Finally, in section 5, we study the impact of transverse modes in the dynamics of the streamwise flow evolution. This study cannot be categorized as an optimal control problem; however, our findings are linked to results in section 4 and may lead to new control strategies for flows with a dominant direction. Our conclusions and a discussion are presented in section 6.
2. Physical problem and hierarchy of models.
2.1. Physical problem and governing equations. We consider a Newtonian fluid with constant density ρ, dynamic viscosity µ, and kinematic viscosity ν, flowing under gravity along a flat infinite 2D substrate inclined at a nonzero angle θ to the horizontal; the schematic of the problem is presented in Figure 1. Same-fluid blowing and suction controls are applied at the substrate surface; the forcing is assumed to be purely perpendicular to the substrate and is introduced into the problem through modifying the impermeability condition at the substrate, but retaining the no-slip condition. The actuator locations are taken to be close  together compared to the interfacial waves, and so we assume that blowing or suction can be performed at all locations on the substrate surface (a study of point-actuated control by the authors will be presented elsewhere). Initially, this presents itself as a boundary control problem; however, it manifests itself as a distributed control problem for the dynamics of the fluid-air interface as we will show. We use coordinates (x, y, z) which are fixed in the plane, with x directed in the streamwise direction, y in the spanwise direction, and z perpendicular to the substrate, as shown for the case of an overlying film in the schematic of Figure 1. Note that, as θ increases, the substrate and axes rotate; we have θ ∈ (0, π/2) for overlying films, θ ∈ (π/2, π) for hanging films, and the vertical film flow with θ = π/2. The surface tension coefficient between the liquid and the surrounding hydrodynamically passive medium is denoted by σ (assumed constant), the local film thickness is denoted by h(x, y, t), a function of space and time, with unperturbed thickness , and the acceleration due to gravity is denoted by g = (g sin θ, 0, −g cos θ). In the case when no blowing or suction is applied, the Navier-Stokes equations (with no-slip and impermeability boundary conditions at the substrate, and the kinematic condition and stress balances at the interface) admit an exact Nusselt solution [53,6] with a film of uniform thickness, i.e., h = , and a streamwise velocity profile which is semiparabolic in z. Velocities and controls are rescaled with the base velocity of this exact solution at the free surface, U 0 = g 2 sin θ/2ν, the lengths scale with , and we rescale the pressure with µU 0 / . We use the nondimensional parameters where the Reynolds number Re measures the ratio of inertial to viscous forces, and the capillary number C measures the ratio of surface tension to viscous forces. The nondimensional Navier- where u = (u, v, w) is the velocity field, p is the pressure, ∇ is the 3D spatial gradient operator, ∇ 2 is the 3D Laplacian operator, and g = (1, 0, − cot θ) is the nondimensional gravitational forcing. We have no-slip conditions at the solid substrate surface, u| z=0 = v| z=0 = 0, and the impermeability condition is modified to account for the controls as w| z=0 = f (x, y, t). The nondimensional kinematic condition and balance of stresses in the two tangential directions and normal direction at the interface are omitted for brevity but are given (with an additional electric field term) in [85].
Actuation through flow boundaries (as is considered in this work) for the Navier-Stokes equations has been considered in the optimal control of turbulent channel flows by Bewley and Moin [8] and in two dimensions with wall slip by Chemetov and Cipriano [14]. Existence and uniqueness results for optimal controls in the case of body forcing controls (where f appears on the right-hand side of the Navier-Stokes equations) are given in [9,1]. Additionally, optimization of wall-parallel velocities for cavity-driven flows has been considered, as well as other problems related to optimization of fluid mixing or thermal convection.

Fully nonlinear Benney equation.
In order to reduce the Navier-Stokes formulation to a forced evolution equation for the film thickness, we make a long-wave assumption. The details of the following derivation are omitted since they are similar to those in [85] for the case of a fully 3D model of an electrified thin film, and the derivation for the 2D case is provided in [82]. We assume that the typical interfacial deformation wavelengths λ are large compared to the unperturbed thickness , set δ = /λ 1, and introduce the change of variables where hats denote O(1) quantities, and the capillary number is rescaled in order to retain the effects of surface tension in the leading order dynamics, with C = O(1). We also assume that the Reynolds number Re is an O(1) quantity. The change of variables (5) is substituted into the governing equations and the hats are dropped. Then, with a systematic asymptotics procedure, in which the flow field is substituted into the kinematic equation, we obtain a fully nonlinear Benney equation in two spatial dimensions, with errors of O(δ 2 ), where we have redefined ∇ = (∂ x , ∂ y ), e x = (1, 0), and ∆ ≡ ∂ 2 x + ∂ 2 y is the 2D Laplacian operator. The variables H and F are approximations of the interface height h and given control f , respectively, correct to O(δ). Notice that F is not only on the right-hand side of (6) but appears also in a nonlinear term in the streamwise direction at O(δ). We have assumed that F t = O(1) in the current variables; otherwise a term involving F t is promoted from the O(δ 2 ) error to the O(δ) terms of (6). As observed by a number of authors [64,30,66], the uncontrolled 1D simplification (no spanwise variation) of this Benney equation exhibits finite-time blow-ups in numerical simulations. The existence of an optimal control for (6) is an open problem (a lack of analytical results for the uncontrolled equation is problematic), and we are not aware of any numerical studies of optimal control for the 1D case or related thin film equations. However, feedback control methods for the 1D version of (6) have been considered by Thompson et al. [81].
2.3. Weakly nonlinear evolution: 2D Kuramoto-Sivashinsky equation. We seek to analyze the evolution of a sufficiently small perturbation about the nondimensional exact constant solution to the Benney equation (6) given by H = 1, F = 0. For this, we substitute H = 1+δη and F = 4δ 2 ζ into (6), where η and ζ are O(1) quantities, and also assume that cot θ is O(1). Truncating terms of o(δ), the resulting equation is Rescaling with and dropping hats gives where κ = 5 cot θ/4Re. Note that the rescaling (8) involves a Galilean transformation, so (9) is in a moving reference frame. For our investigation of optimal controls with actuators at every location on the substrate, the advection term makes no difference and the result in the moving frame can be translated back into the "lab" frame without any issue. Similarly, the advection in the streamwise direction has no effect on the purely spanwise forcing we study in section 5. However, the case of point-actuated controls warrants the consideration of traveling actuator grids. We supplement the 2D KSE (9) with periodic boundary conditions on the rectangle Q = [0, L 1 ] × [0, L 2 ]. For this purpose, we denote the wavenumber vectors byk, with components for k ∈ Z 2 , so that η and ζ may be written in terms of their Fourier series as Here η −k and ζ −k are the complex conjugates of η k and ζ k , respectively, since both the solution and control are real-valued. Equation ( for each k ∈ Z 2 . Since (9) governs the perturbation of the flat film state, we consider initial conditions, denoted by v, with zero spatial mean (v 0 = 0). From (12) with k = 0, the evolution of the mean depends on the control ζ as We restrict to controls with ζ 0 ≡ 0, so that the zero spatial mean of the solution is preserved. The controls considered in this work can thus be thought of as moving fluid from one location to another, while still conserving the total fluid volume. From the ODE description (12), it can be seen that the dynamics of the transverse modes (i.e., k 1 = 0) are linear, being governed by (14) d dtη Here we have denoted the transverse Fourier coefficients η (0,k 2 ) byη k 2 , with the same notation for the transverse components of the forcing ζ. We letP denote the projection onto the subspace of transverse modes and defineη(y, t) to be the image of η under this projection, η =P η. This projection satisfies the linear PDE which is equivalent to the ODE system (14), and may be obtained from (9) by averaging over the streamwise direction. From (12), it can be seen that the dynamics of the streamwise and mixed modes are slaved toη; i.e., the transverse modes only decouple partially from the full nonlinear system. The parameter κ encodes the incline of the substrate, and in the absence of controls, we have three distinct dynamical regimes depending on its value. For κ < 0 we obtain hanging films, and, for ζ = 0, a range of transverse modes with 0 < |k 2 | < (−κ) 1/2 are linearly unstable. Since the governing equation (14) is linear, unbounded exponential growth occurs with rate −κk 2 2 −k 4 2 > 0. A large focus of this work is on stabilizing hanging films with growing transverse modes. For κ ≥ 1 the flow is overlying and the Reynolds number is subcritical, with κ = 1 corresponding to the critical Reynolds number Re c = 5 cot θ/4. In this case, simple energy estimates may be used to show that all solutions converge to zero in the absence of blowing or suction; importantly, the quadratic nonlinearity facilitates the transfer of energy, rather than being a sink or source of energy. Overlying flows with supercritical Reynolds numbers are then found for 0 ≤ κ < 1. When κ = 0, the film is vertical, and the canonical equation (9) reduces to a forced version of the thin film equation obtained by Nepomnyashchy [50,51], Without controls, this equation has been studied both analytically and numerically by a number of authors. Pinto [59,60] proved the existence of a finite-dimensional global attractor and instant analyticity of solutions. Numerical studies have been carried out by Akrivis et al. [3] and Tomlin, Kalogirou, and Papageorgiou [84], where it was shown that solutions possess a finite energy density independent of the periodic domain size with chaotic dynamics emerging on sufficiently large domains. Similar Kuramoto-Sivashinsky-type dynamics are found for 0 < κ < 1.
We apply controls over a finite time interval [0, T ] in order to drive the solution to a zero-mean desired state η(x, t). The desired state is not required to be an exact (possibly unstable) solution of the uncontrolled problem, although many control methodologies make this assumption. For our setting, we employ homogeneous spatial norms, with all functions assumed to be Q-periodic, and have zero spatial mean for all times; the spaces H s 0 (where are defined through their respective norms: Note that the Fourier symbol of the operator (−∆) s/2 is |k| s , so the L 2 0 -norm of ∆η equals the H 2 0 -norm of η, for example. Inner products for the H s 0 and L 2 (0, T ; H s 0 ) spaces are denoted by angled brackets with the appropriate subscripts. The space of admissible controls is denoted by F ad , a nonempty, closed, convex subset of L 2 (0, T ; L 2 0 ), where the norm is defined by (17b) with s = 0. We optimize with respect to the cost functional where s ∈ R and γ > 0 (γ = 0 allows infinitely strong controls). We denote the three components of the cost functional by C s,γ , and C s,γ , respectively. The optimal control ζ * and associated optimal state η * (which solves the 2D KSE (9) with ζ = ζ * ) are defined as minimizers of the cost (18) over all controls in F ad ; i.e., if ζ ∈ F ad has associated state η , then C s,γ (η * , ζ * ; η) ≤ C s,γ (η , ζ ; η). The parameter γ in (18) can be thought of as the cost of using the control, relative to the cost of inaccuracy between η and the desired state η. For small γ we may use larger controls to ensure that the solution is very close to the desired state, but for large γ the controls are expensive, and the difference between η and η is less important. Larger values of the Sobolev index s have the effect of increasing the weighting on the higher frequencies with larger wavelength modes in the solution costing relatively less. We note that more general wavenumber weighting functions can be used which do not necessarily correspond to Sobolev norms. We include the "payoff term" C The definition of the spatial norm used here is a measure of energy density; norm (17a) is related to the usual integral definition for the square of the Sobolev norm through multiplication by L 1 L 2 . Thus, the cost functionals for different choices of the domain Q become comparable as they are independent of the underlying periodicities; this makes sense if we think of our problem as being on an infinite domain on which we impose periodicity, rather than on a bounded domain with periodic boundary conditions. Taking this definition is natural for the study of (9), since it has been observed that (16) without forcing (ζ = 0) possesses a finite energy density as the domain size is increased [84]. Furthermore, values of the norms are comparable to the amplitude of the interface.
For the numerical study of (9) on Q-periodic domains, we utilize a family of implicitexplicit backwards differentiation formula (BDF) methods constructed by Akrivis and Crouzeix [2] for a class of nonlinear parabolic equations under appropriate assumptions on the linear and nonlinear terms. They considered evolution equations of the form where A is a positive definite, self-adjoint linear operator, and B is a nonlinear operator which satisfies a local Lipschitz condition. It was shown that these numerical schemes are efficient, convergent, and unconditionally stable. For us, we have the addition of forcing ζ to the right-hand side, and thus the operators are defined as where c is chosen to ensure that A is positive definite. These schemes were utilized for the corresponding 1D optimal control problem by Gomes, Papageorgiou, and Pavliotis [25], and their applicability was checked for similar multidimensional problems in [3] (including a convergence study) and [85] for a nonlocal problem. In order to perform computations, we truncate the Fourier series (11) to |k 1 | ≤ M , |k 2 | ≤ N , which corresponds to a discretization of the spatial domain Q into (2M + 1) × (2N + 1) equidistant points, and carry out timeintegration of the system in Fourier space using the BDF methods.
3. Optimal transverse control for hanging films. We first consider controlling the transverse instabilities present for hanging films (κ < 0) by applying controls of the form ζ =ζ(y, t) (withζ 0 = 0). We work under the assumption that the dynamics of (9) are bounded if the linear growth in the spanwise dimension is controlled. This is a reasonable assumption given the form of the ODE system (12) 1 and is confirmed by our numerical results. The ODE system (14) has the explicit solution Either side of κ = 0, there are unstable streamwise and mixed modes, governed by (12); decreasing κ only increases the strength of the destabilizing terms relative to the nonlinearity and fourth order damping (and increases the number of unstable modes, similar to the effect of increasing L1 and L2), rather than qualitatively changing the structure of the nonlinear system. The critical value of κ = 0 is only important for the transverse system (15), as it is the point at which unstable modes first appear.
whereṽ k 2 are the transverse Fourier coefficients of a given initial condition v. Equation (21) is often referred to as the control-to-state map. LetP Ξ denote the projection onto the subspace spanned by the unstable purely transverse modes, where the wavenumbers belong to the set Ξ = {k 2 ∈ Z | 0 < |k 2 | < (−κ) 1/2 }, and defineη u =P Ξ η. The transverse modes in Z\Ξ either are neutrally stable or decay exponentially without controls, and have no effect on the boundedness of the solution. Thus, in this section, we take F ad to be the image of L 2 (0, T ; L 2 0 ) under the operatorP Ξ . In fact, the linearity of (14) permits the explicit construction of an optimal control which is smooth in both time and space. If no modes are neutrally stable, all uncontrolled transverse modes are damped, and the long time dynamics will converge to those of the full system (9) withη = 0 if the modes in Ξ are controlled to zero. The desired state is taken to be the orthogonal projection η = (I −P Ξ )η = η −η u , with no components of unstable transverse modes. It is important to note that due to the slaving of the streamwise and mixed modes to the transverse modes, η is not the same as the evolution of the initial condition (I −P Ξ )v under the dynamics (in other words, the evolution under the PDE and the projection I −P Ξ do not commute). The desired state is dependent on the solution locally in time-not a prescribed function. The cost functional (18) simplifies to (22) C s,γ (η,ζ) = 1 2 where the summands are costs for each individual mode. We denote the optimal solution and control with respect to the cost functional by η * andζ * , respectively. For finite T , this results in a set of linear ODE control problems, each a linear quadratic regulator (LQR) problem, for which there is a well-developed theory (see [93], for example). The optimal control (derived in Appendix A) is given byζ * k 2 = r k 2η * k 2 (with stars denoting optimality), where r k 2 satisfies a Riccati equation and final time boundary condition By defining the roots we can give the solution to (23) explicitly as (25) Furthermore, the value of the cost functional attained by the optimal control is Since the problem under consideration is linear and finite-dimensional, these optimal controls are unique. For comparison, in the infinite time-horizon case, i.e., T = ∞, we findṙ k 2 = 0 so that r k 2 = λ − (this is the correct root to stabilize the system (23)). Substituting back into (14) yields the solution where the optimal control is defined byζ * k 2 = λ −η * k 2 . Now we test the optimal transverse control constructed above numerically. So that unstable transverse modes are present, we consider two cases of hanging films, taking parameters (i) κ = −1, L 1 = 40, L 2 = 15, and T = 100, and (ii) κ = −0.25, L 1 = 40, L 2 = 27, and T = 400. The final time T differs between the two simulations as the timescale of the dynamics is dependent on κ, which varies between the two cases. For (i) and (ii), the (0, ±1) and (0, ±2) modes are unstable (all other transverse modes are damped), and we take the following initial condition, which contains contributions from these modes: We chose L 1 to be sufficiently large in both cases so that the dynamics of the controlled solutions are nontrivial, i.e., unstable streamwise and mixed modes are present. Figure 2 shows the optimal controls successfully inhibiting the growth of transverse modes for cases (i) and (ii) in numerical simulations for the choices of s = 2 and γ = 1. For case (i), the uncontrolled solution behaves as sin(4πy/L 2 )e λt (the (0, ±2) modes are more unstable than the (0, ±1) modes) with exponential growth rate λ = −κ(4π/L 2 ) 2 − (4π/L 2 ) 4 ≈ 0.209. With the application of the optimal transverse controls, the solution energy, plotted in Figure 2(a), remains bounded as desired, with a modal steady state (dominated by the (m, 2m) modes for m ∈ Z) emerging as shown in panel (e). In case (ii), the (0, ±1) modes dominate the uncontrolled solution with linear growth rate λ ≈ 0.0106, and the optimal controls successfully prevent unbounded growth, revealing a chaotic attractor; this is apparent from the L 2 0 -norm of η * plotted in Figure 2(b). The cost functional takes the values 3.27×10 −3 and 7.24×10 −4 with the optimal controls for the respective cases, in agreement with analytical values computed with (26). We observe in both cases that the H 2 0 -norm of the transverse modes and the L 2 0 -norm of the optimal controls decay exponentially for the majority of the time interval [0, T ]; see Figure 2(c),(d). For midrange times, the controlled transverse modes decay with exponential decay rate proportional to the uncontrolled linear growth rate, and r k 2 ≈ λ − (as in the infinite time-horizon problem). For example, in case (i), the most unstable (0, ±2) modes are controlled to zero with the greatest decay rate. Figure 3 shows how the (analytical) optimal value of the cost functional (26) for case (i) is affected by changes in s and γ; we also validate our numerical results by superimposing the costs obtained in simulations with s ∈ {0, 1, 2, 3} and γ ∈ {0.1, 1, 10}. We observe that C * s,γ is a decreasing function of s and an increasing function of γ; this is not surprising given the definitions of the norms (17a,b) and the cost functional (22). These cases all provide plots similar to those found in Figure  2  It now remains to comment briefly on the eventual dynamics of the system withη = 0 resulting from application of the above controls over a large time interval. Interestingly, we found that as the gravitational instability is strengthened (κ is decreased), the usual chaotic dynamics observed for the vertical film case in [84] are controlled to attractors of steady modal solutions. It would be expected that a forwards Feigenbaum cascade occurs for overlying films as κ decreases to zero (for fixed L 1 and L 2 ) given the results of Smyrlis and Papageorgiou [78] for the 1D problem; thus films close to vertical (controlled if necessary) appear to be the most prone to chaotic dynamics. With the length parameters L 1 = 40 and L 2 = 27, as in case (ii), we find mostly chaotic dynamics as κ is decreased to −0.44, below which a time-periodic attractor emerges. This behavior dominates until κ = −0.58, when steady modal waves dominated by the (m, 3m) modes for m ∈ Z appear; this is similar to the (m, 2m)-modal wave in Figure 2(e).
As shown in this section, only the unstable transverse modes are responsible for the unbounded behavior of the 2D KSE (9) with κ < 0. Nevertheless, this linear control theory may be extended to all transverse modes, where the summations in the previous lines extend from Ξ to Z, in order to drive the stable transverse modes to zero optimally. Furthermore,η u orη does not necessarily need to be controlled to zero (see Appendix A). It is also worthwhile to note that weakly nonlinear theory predicts that constant transverse forcing (as would arise from spanwise substrate corrugation) would not be a successful control method.
4. Full optimal control. In this section, we consider the general optimal control problem for the 2D KSE (9), where the set of admissible controls F ad is a nonempty, closed, convex subset of L 2 (0, T ; L 2 0 ). We consider the existence of an optimizer for our problem and then give the methodology for numerical simulations. In contrast to section 3, we must resort to an iterative algorithm for the full nonlinear problem. In the numerical simulations that follow the analysis, we take F ad = L 2 (0, T ; L 2 0 ). Following the abstract formulation in [43,80] [59] for the case of ζ = 0 and κ = 0). Theorem 4.1. For initial condition v ∈ L 2 0 and ζ ∈ L 2 (0, T ; L 2 0 ), there exists a unique solution η to (9) in X 1 : . We shall use both parts of the above theorem for the proof of existence of an optimal control. Uniqueness of an optimal control is not guaranteed as the optimization problem is not convex; this is due to the ηη x nonlinearity. The following theorem and proof make no deep assumptions (such as analyticity or regularity in higher Sobolev spaces) which would be expected to be true for equations such as (9), only requiring the above existence and uniqueness theorem. Thus, the proof is given in a very general framework which may be applied to similar problems. However, this restricts us in the range of the index s, and the regularity of the initial condition and desired state.
. Define F ad to be a nonempty, closed, convex subset of L 2 (0, T ; L 2 0 ). Then, there exists an optimal control ζ * ∈ F ad for the 2D KSE (9) with initial condition v which minimizes the cost functional C s,γ defined by (18).
Proof. For a control ζ, we write the solution of (9) with the given initial condition v in terms of the control-to-state map as η(ζ; v), through which we define the reduced cost functionalC s,γ (ζ) = C s,γ (η(ζ; v), ζ; η). Since v ∈ H 2 0 ⊂ L 2 0 , the first part of Theorem 4.1 implies that (η(ζ; v), ζ) ∈ X 1 × F ad . The optimal control problem can thus be recast as the problem of finding the minimizer ofC s,γ (ζ) over F ad . It makes sense to check that there exist ζ ∈ F ad which give a finite value ofC s,γ . The boundedness of the first and last terms of the cost functional (18) are consequences of (η(ζ; v), ζ) ∈ X 1 × F ad (recall that s ≤ 2) and the regularity of η. The final time component is the problematic term which requires the additional regularity of the initial condition: By the second part of Theorem 4.1, we know that η| t=T ∈ H 2 0 ⊆ H s 0 , and since η| t=T ∈ H s 0 , the final time component of (18) is also finite. Thus any forcing in the space of admissible controls yields a finite cost. Since C s,γ ≥ 0, the reduced cost has a finite infimum, (29) inf ζ∈F adC s,γ (ζ) = c ≥ 0.
We cannot yet say that this infimum is attained by a control in the admissible space; however, we know that there is a minimizing sequence {ζ (n) } ∞ n=1 ⊆ F ad with associated states defined by η (n) = η(ζ (n) ; v) such that (30) lim n→∞C s,γ (ζ (n) ) = c.
Without loss of generality, we may assume that the corresponding sequence of costs is monotonically decreasing. From the form of the cost functional (18), we know that {ζ (n) } ∞ n=1 is bounded uniformly in L 2 (0, T ; L 2 0 ), i.e., there is some constant r ≥ 0 (r =C s,γ (ζ (1) ) suffices) such that (31) ζ (n) L 2 (0,T ;L 2 We define K to be the intersection of F ad with the closed ball of radius r in L 2 (0, T ; L 2 0 ). K is a closed, convex, and bounded subset of the Hilbert space L 2 (0, T ; L 2 0 ) and thus is weakly sequentially compact (see Theorems 2.10 and 2.11 in [87]). Then, the sequence {ζ (n) } ∞ n=1 ⊆ K has a weakly convergent subsequence ζ (n) ζ * for the topology of L 2 (0, T ; L 2 0 ) (not relabeled for simplicity) with weak limit ζ * ∈ K ⊆ F ad . For more general forms of the cost functional, it may be necessary to assume that F ad is bounded to make this step. The function ζ * is a candidate for the optimal control.
Multiplying the 2D KSE (9) by the solution and taking the spatial average, we obtain the inequality (see Appendix B for the derivation) L 2 (0,T ;L 2 0 ) ), where the constant C depends on κ and T only. With this estimate and (31), the sequence {η (n) } ∞ n=1 is uniformly bounded in X 1 . Since X 1 is the dual of a separable Banach space, 2 by the Banach-Alaoglu theorem, there is a subsequence (not relabeled for simplicity) with It also follows from (31), (32) that the sequence with terms is uniformly bounded in L 2 (0, T ; H −2 0 ); a Poincaré inequality for the continuous embedding of H 2 0 in L 2 0 is also utilized. To deal with the nonlinearity, we use the estimate where the first inequality follows from η (n) η (n) x = ∂ x (η (n) ) 2 /2 and the definition of the H −1 0norm (17a), and the last uses an Agmon inequality proved in Theorem 4.1 of [28] and again a Poincaré inequality for the embedding of H 2 0 in L 2 0 . The constant C involves the constants from both the Agmon and Poincaré inequalities. Squaring and integrating in time leads to (36) η (n) η (n) where we have used (31), (32). Then the sequence {η (n) η (n) x } ∞ n=1 is bounded uniformly in L 2 (0, T ; H −1 0 ), and thus bounded uniformly in L 2 (0, T ; H −2 0 ) again by a Poincaré inequality. With these results and the uniform boundedness of {ζ (n) } ∞ n=1 , we can conclude that {η is bounded uniformly in L 2 (0, T ; H −2 0 ). With this and the weak convergence result (33a), it follows from Theorem 8.1 in [65] (a compactness result) that (37) η (n) → η * strongly in L 2 (0, T ; H s 0 ) for s < 2.
2 The predual of X1 is the direct sum space L 1 (0, T ; (L 2 0 ) * ) + L 2 (0, T ; (H 2 0 ) * ), where stars denote duals; this is the direct sum of two separable Banach spaces and is thus separable Banach space itself (endowed with an appropriate norm). We now check that the weak limits satisfy the control-to-state map, η * = η(ζ * ; v); the pair (η * , ζ * ) must satisfy the governing equation (9) and the initial data, η * | t=0 = v. The strong convergence result (37) with s = 1 implies that the nonlinearity converges weakly in L 2 (0, T ; (H 2 0 ) * ), where the star denotes the dual, as follows: Let w ∈ L 2 (0, T ; H 2 0 ); then whereC is the constant corresponding to the continuous embedding of H 1 0 in L 4 0 in two space dimensions. The right-hand side of (38) converges to zero as n → ∞, and thus Also, the sequence with terms (34) converges weakly to its corresponding limit in the same space. Since {η , it is weakly convergent, and furthermore by a density argument we can conclude that the weak limit is η * t . Since the sequence {ζ (n) } ∞ n=1 converges weakly in L 2 (0, T ; L 2 0 ) ⊂ L 2 (0, T ; (H 2 0 ) * ), by uniqueness of weak limits we may conclude that (40) η The proof that the optimal state has initial condition v follows similarly to the arguments in Theorem 9.3 of [65].
We remark that the above theorem holds for v ∈ L 2 0 and s ≤ 0 by a similar proof, and extension to any s ∈ R would be viable with analyticity results.
We now present the adjoint framework which allows the construction of such minimizing sequences; this forms the basis of the iterative algorithm we will employ for our numerical simulations. The details of the following methodology (including numerical aspects) can be found in [40], for example. The Lagrangian of the optimization problem is where the L 2 (0, T ; L 2 0 )-inner product corresponds to the energy density norm defined by (17b). Here, p is known as the adjoint variable, which is, in essence, a Lagrange multiplier. The firstorder conditions of optimality for a local optimizer (η * , ζ * , p * ) consist of the governing equation where I − P 0 is the projection onto the space of functions with zero mean and we recall that (−∆) s is the fractional Laplacian of order 2s with Fourier symbol |k| 2s . Both (9) and (43) are necessarily satisfied by a local optimizer (η * , ζ * , p * ). The adjoint equation is backwards in time, being supplied with the final time "transversality" condition which arises due to the payoff term in the cost functional. Finally, we obtain a variational inequality by taking the Fréchet derivative of L with respect to ζ, This inequality informs us that the optimal control ζ * is in the direction of −(γζ + p) from the current control ζ, from which we may construct iteration schemes. Such updates move along the local approximation to the curve of steepest descent. The exclusion of the zero mode, k = 0, in the L 2 (0, T ; L 2 0 )-inner product is vital to ensuring that the control remains in the space of zero average functions. The projection I − P 0 appearing in (43) guarantees that the spatial average of the adjoint variable is fixed at zero (ηp x is not a zero mean function in general). Then, the update direction −(γζ + p) also has zero spatial average, and thus any iteration scheme will yield a sequence of controls preserving this property.
Next we detail the numerical algorithm which is employed to approximate a local optimizer (η * , ζ * , p * ); checking that this is the global optimizer for our infinite-dimensional problem is difficult. In simple terms, the forward-backward sweep method is comprised of iterated simulations of (9) and (43), with control updates after each iteration. Note that (43) may be written in the form (46) − p t + Ap = B (p; η, η), where A is defined by (20a). Thus, the BDF methods outlined in subsection 2.3 are applicable to this backwards in time equation with the same value of c, provided that B satisfies the necessary Lipschitz bounds. The pseudocode for the forward-backward sweeping method is given in a simplistic form in Algorithm 4.1. The returned values (η (n) , ζ (n) ) are a minimizing sequence as in Theorem 4.2, with the limit approximating the optimal state and control pair, (η * , ζ * ). For our control update, we take a step of size c (n) /γ in the direction of −(γζ (n) +p (n) ), which, as above, can be rewritten as a convex combination of the control and adjoint variables at the current iteration step. This classical method is found to perform consistently well compared to more complicated methods which may converge "too quickly"-see the discussion in Chapter 4 of [40]. Rather than performing an expensive line search across a range of convex Choose initial state η 0 (x), desired target state η(x, t), and initial control guess ζ (1) (x, t) for t ∈ [0, T ]. Take c (1) ∈ (0, 1), set c (n) = c (1) for n ∈ N >1 , and choose a tolerance τ 1. Initialize with n = 0, C s,γ (0) = 0, d = τ + 1. while d > τ do n = n + 1. Solve (9) with initial state η 0 and control ζ (n) to obtain η (n) (x, t) for t ∈ [0, T ]. Calculate C s,γ (n) := C s,γ (η (n) , ζ (n) ; η) and update d = |C s,γ (n) − C s,γ (n − 1)|. Solve (43) backwards given final time condition (44) and state η (n) to obtain while C s,γ (n + 1) ≥ C s,γ (n) do Solve (9) with initial state η 0 and control ζ (n+1) to obtain combinations, i.e., seeking a minimum as c (n) is varied, we employ an adaptive step-halving scheme (a backtracking line search) starting with steps of size c (1) = 0.1 (see the update formula in Algorithm 4.1); if a control update results in an increased cost, then the step size (and the successive step sizes) are halved until the updated control yields a lower cost than the previous iteration. The sequence of values of the cost functional can be used to indicate a posteriori that we have achieved convergence to a local minimizer. We have checked the results of our numerical schemes with more complicated updating and searching methods with good agreement.
For the first numerical experiment, we take L 1 = L 2 = 21 and cover the three different dynamical regimes and two critical points of the system with κ ∈ {−0.5, 0, 0.5, 1, 1.5}. We again take the initial condition defined by (28) and choose the desired state η to be the "snaking" transverse wave shown in panel (a) of Figure 4; this corresponds to a steady solution of the nonlocal 2D KSE studied in [85], shown in their Figure 6(c). Controls are applied until the final time T = 5, and we also take parameters s = 0 and γ = 1 for the cost functional. The case of zero forcing is used as the initial control guess, ζ (1) = 0, thus allowing a cost comparison between the uncontrolled and optimally controlled systems. The results of the forward-backward optimization procedure (Algorithm 4.1) for κ = −0.5 are shown in panels (b)-(d) of     the majority of the time interval. The evolution of η * , η * − η, and ζ * for this case is presented in M119390 01.mov [local/web 8.87MB], and we note that the observed optimal control is not too far from a proportional control, i.e., −ζ * ∼ η * − η. The sequence of costs attained by the minimizing sequence (η (n) , ζ (n) ) is displayed in Figure 4(d), including a breakdown into its three components C (j) 0,1 (n) for j = 1, 2, 3. The dotted lines with points at fractional values of the iteration number n correspond to the cases where the updated control results in an increased cost. In these cases, the step size is halved and the forward sweep is repeated, as shown in the plot of c (n) (faded dashed line), taking values on the right axis. The cost breakdowns for the uncontrolled and optimally controlled cases are given in Table 1 for each κ. Rows 2-4 correspond to the uncontrolled case, with cost components C 0,1 = 0 for ζ = 0) and total cost C 0,1 ; these are observed to be decreasing functions of κ. The next four rows give the optimal cost and its breakdown into the three components; these asymptote values were obtained by fitting the sequences of costs to functions of the form ae bn + c. The increase in the optimal cost with κ over these examples may be explained as follows: Since we start from small amplitude initial conditions and such a short time interval [0, T ] for optimization, the linear instabilities work in favor of controlling the solution to the desired state; for κ = 1.5, the control must be strong since the zero solution is exponentially stable. However, we expect a turning point in this behavior as κ becomes very negative (dependent on T ), when the cost of controlling the linear instabilities outweighs the control cost saved due to the linear instabilities aiding the solution growth towards η. Row 9 in Table  1 gives the number of iterations required to ensure that the change in the minimizing sequence of costs is below the desired tolerance τ = 5 × 10 −6 ; this becomes very large as κ decreases further beyond −0.5. Finally, the minimum of the L ∞ -norm of η − η * across the whole time interval is given in row 10; the optimal controls approach the desired state more closely in the L ∞ -sense as κ decreases (for the range of numerical simulations we completed). We also performed numerical experiments over a range of s and γ, although not shown here, and we report that the dependence of the optimal cost C * s,γ on these parameters was the same as in the linear case shown in Figure 3, i.e., a decreasing function of s and an increasing function of γ; the same trend was observed for the individual components, C    L 1 = 32 and L 2 = 21. We take three cases, (i), (ii), (iii), of initial condition and timedependent desired state: where δ (ii) and δ (iii) are 1 for their respective cases and zero otherwise. In case (i), both initial condition and desired state are spatially 1D, whereas in cases (ii) and (iii), a mixed-mode term is added to either v or η, respectively. No transverse modes are included in v or η in any of the cases, nor are these modes linearly unstable since κ = 0. We additionally consider We also fix the final time T = 5 and set control parameters s = 0 and γ = 1 as in the previous numerical experiment. In Figure 5(a), we show a portion of the minimizing sequence of costs for the three cases of (v, η), each with the three choices of ζ (1) given in (49). The limits of the sequences are independent of the initial control guess, and we confirm from inspection of the numerical solution that the approximations of ζ * also agree (this is not necessarily implied by the previous statement). This lends credence to the possibility that the obtained optimizers are global minima of the cost functionals. Figure 5(b),(c),(d) plots the L 2 0 -norms of the optimal control and state over [0, T ] for the three cases. For case (i), both ζ * (and η * ) remain 1D for the entire time interval; the optimal control obtained here is the same as what would be obtained from the 1D simplification of the control problem. In cases (ii) and (iii), two-dimensionality enters into the problem via the initial condition and desired state, respectively, through the addition of a mixed-mode term. We find not only that mixed modes appear in ζ * but also that the projection onto the transverse modes is nontrivial. For these two cases, the L 2 0 -norms of the projectionsP (η * − η) andP ζ * are included in Figure 5(c),(d), and visualizations of the time evolution are given in M119390 02.mov [local/web 8.57MB] and M119390 03.mov [local/web 9.23MB], respectively. This can be understood with the adjoint equation (43), where it can be seen that mixed-mode and streamwise mode activity in the solution excites the purely transverse modes in the adjoint variable. With these results, we make the following conjecture about the spatial dimension of the optimal state and control pair: Assuming that the space of admissible controls F ad cannot be spanned by functions of x and t alone (F ad contains functions dependent on y), the optimal control ζ * (and η * ) is independent of y if and only if v and η are independent of y. Moreover, if F ad contains functions which are independent of x, i.e., transverse modes, and ζ * is dependent on y, theñ P ζ * is nonzero. The presence of transverse modes in the optimal control even if none are present in v or η may at first appear unusual. However, the system of transverse modes decouples partially, and the transverse modes in the control influence the dynamics of the streamwise and mixed modes through the nonlinearity. We investigate the effect of transverse mode forcing in more detail in the next section, analyzing the response of the interface energy to spanwise blowing and suction patterns. We also note that, as in the previous numerical experiment, the optimal states and controls appear to be close to proportional for the majority of [0, T ] (see M119390 02.mov [local/web 8.57MB] and M119390 03.mov [local/web 9.23MB]).
In the situation when F ad is a strict subset of L 2 (0, T ; L 2 0 ), the numerical procedure is exactly as outlined above, but with projections of the updated controls onto F ad at each iteration; this is viable due to the linearity of the adjoint equation (43) in p. We note that both Theorems 4.1 and 4.2 are valid without the zero-average restriction, i.e., taking F ad ⊆ L 2 (0, T ; L 2 ), where the latter space is given an appropriate norm. We repeated a number of the above numerical experiments in this setting, allowing controls in L 2 (0, T ; L 2 ). In this case, the projection P 0 is removed from the adjoint equation. We found that the optimal controls in this space caused large drifts in the spatial average of the solution. Physically

Transverse mode effects.
In this section, we briefly examine the extent to which controlled transverse modes can affect the streamwise and mixed modes through the nonlinear coupling; see the ODE system of Fourier coefficients (12). In this way, purely transverse controls (such as those studied in section 3) may be thought of as indirect controls on the full system. This may be useful in physical situations where it is easier to force (and/or observe) the transverse modes than the other components of the flow.
We focus on the vertical film set-up, κ = 0, which has been studied extensively in [84]. In this case, the transverse modes of uncontrolled solutions are damped for any choice of Q with chaotic dynamics emerging for sufficiently large L 1 and L 2 . In [84], the authors found that the time-averaged energy density behaves as (50) η := lim for sufficiently large length scales; the quantity E L,α considered in [84] is related to η 2 by the factor of L −1 1 L −1 2 due to our definition of the L 2 0 -norm (17a). The estimate (50) is additionally found numerically for the 1D KSE (1), although current analytical bounds are not sharp. The goal of this section is to determine if purely transverse controls may be used to decrease η below 1, i.e., make the fluid interface less energetic on average. Another appealing property of a controlled system would be the regularization of chaos.
The appropriate choice of desired state for this study is η = (I −P Σ )η +ψ(y, t), where Σ is the set of transverse Fourier modes which we can force (P Σ is a projection onto these modes) andψ is a real-valued function of the modes in Σ alone, with the property that ψ −k 2 is the complex conjugate ofψ k 2 . Transverse modes which are not included in Σ decay exponentially for the vertical film case under consideration here; thus we do not take them in our initial conditions to prevent any transient effects. The derivation of optimal controls for this problem is provided in Appendix A, generalizing the controls used in section 3. Rather than specifying any particular control, we assume that (suboptimal) controls are applied so that the desired state is achieved exactly for all time, i.e.,η =ψ. This is reasonable since we have full reachability and controllability for the individual transverse modes; the ODEs (14) allow explicit construction of a controlζ for a given stateη which varies continuously fromP v toψ. With this, we may focus on the flow response to transverse modes with a fixed amplitude.
We consider the simple choice ofψ = A sin(k 2 y) for an amplitude A ≥ 0 and transverse wavenumberk 2 , only containing contributions from the (0, ±k 2 ) modes. We fix the domain lengths L 1 = 120, L 2 = 30 for numerical simulations throughout this section, in which case the uncontrolled system evolves chaotically. Additionally, we take random (real-valued) initial conditions with unstable streamwise and mixed modes so that the solution enters the global attractor rapidly. As noted above, we do not include transverse modes other than the (0, ±k 2 ) modes in the initial condition, a choice which is justified by the numerical results themselves. The time-averaged energy density for both η −η and η is plotted in Figure 6 for k 2 = 1, . . . , 5 and a range of A. These time-averages are approximated by averages over a large time interval [T 1 , T 2 ] after the solution has entered the global attractor (see [84] for the details). Panel monotonic in most cases). Furthermore, we find that the chaotic dynamics are regularized for a sufficiently large value of A, and the transverse wave for a given mode eventually becomes nonlinearly stable where the line in Figure 6(a) touches down on the A-axis; these critical values and the behaviors of the energy densities are heavily dependent on k 2 . For k 2 = 1, the chaotic dynamics are regularized for A ≈ 5.8 and the transverse sine-wave becomes nonlinearly stable at A = 10.22. For k 2 = 2, this happens for much lower amplitudes, with chaos being regularized for A just beyond 2, and the trivial solution in the streamwise and mixed modes (η −η = 0) becoming stable at A = 2.86. Panel (b) shows the more surprising result that a small amplitude transverse wave can lower the total average energy of the system (however, this does not account for any control costs). Due to the form of the nonlinear coupling, we postulate that very high frequency transverse waves will have little effect on the system energy and dynamics (as is confirmed by linear results to follow). We find that the choices of k 2 = 3, 4 are particularly successful in decreasing η ; these results suggest the existence of an optimal mode of control and amplitude (ifψ is restricted to a steady modal wave) to minimize the system energy, dependent on κ, L 1 , and L 2 . From our time-dependent simulations, we also observed that larger values of A delayed the onset of chaos below the threshold of regularized dynamics.
The numerical results in Figure 6(a) show that full nonlinear stability of the transverse wave occurs at the lowest values of A for k 2 = 3, 5. We now use linear theory to investigate the relationship between this critical value A c and k 2 .  Figure 7. Critical amplitude Ac which separates the linearly stable and unstable regimes of (52) plotted against k2.
system of Fourier modes, the linearization of (12) about η =η =ψ for a generalψ(y, t) is for m ∈ Z 2 with m 1 = 0. For the case of κ = 0 andψ = A sin(k 2 y) this becomes The stability of the above infinite-dimensional linear system can be determined by considering the stability of sufficiently large finite-dimensional truncation. With this, we may compute the critical value A c , which separates the unstable and stable regimes, as a function of k 2 .
The results are plotted in Figure 7, and the linear theory agrees with the critical values of the amplitude found in our nonlinear time-dependent simulations when the transverse wave becomes attractive to all initial conditions. Although not plotted, A c for k 2 ≥ 10 is monotonically increasing, as suggested by the form of the nonlinearity. Although not as informative as Figure 6, the results shown in Figure 7 are a good predictor of whether the flow has a weak or strong response to a given frequency. The consideration of more generalψ and different parameter choices is beyond the scope of our study. Transverse controlsζ allow us to attain any such transverse wave state (which then may be thought of as a control in its own right); however, the reachability and controllability of the streamwise and mixed-mode system for the controlψ is unclear, in which the control acts multiplicatively through the nonlinear coupling, rather than additively. We believe that these results may have applications to other control problems for multidimensional flows with a dominant direction-for example in aerodynamics.
6. Conclusions. In this paper we considered the distributed control of a Kuramoto-Sivashinsky equation for gravity-driven thin film flow overlying or hanging from a 2D flat substrate. Blowing and suction controls applied at the substrate surface appear as a forcing term in the weakly nonlinear evolution equation. For hanging films (κ < 0), optimal controls which are constant in the streamwise direction were constructed in section 3 to impede the exponential growth of linearly unstable transverse modes; these controls were applied successfully in numerical simulations. This spanwise instability is physical, predicting the formation of rivulets which may be a precursor to dripping for certain parameter regimes; see the experiments performed by Brun et al. [13]. In a nonidealized situation, it may be much more difficult to prevent such an instability from developing. It may also be the case that the path to dripping may take a different route upon the control of this initial instability, bypassing rivulet formation. Will controls merely delay the onset of dripping? It is clear that weakly nonlinear analysis alone will not give an acceptable answer as the processes of rivulet formation and dripping are inherently nonlinear. The use of stabilizing electric fields to prevent the dripping instability (blow-up) of the rivulets is currently under investigation by the authors. Furthermore, it may be viable to use the constructed optimal controls to develop feedback controls for the manifestation of the same spanwise instability in the models higher up the hierarchy-for example the Benney equation (6). This instability appears nonlinearly in these models, and thus the explicit construction of optimal controls is not possible, and iteration methods are computationally expensive.
In section 4, we considered a more general class of controls which varied spatially in both the streamwise and spanwise directions yet restricted to zero spatial average; without this assumption, the spatial average of the solutions was seen to vary greatly. A detailed proof of existence of an optimal control was given, outlining the general strategy that may be used for similar problems. Using the adjoint formulation, we constructed a forward-backward sweeping algorithm for the problem and successfully applied it in numerical experiments. The optimal controls for the problem under consideration are dependent on the initial condition and parameters, and a large number of iterations is required for the numerical algorithm to converge to the optimizer. Furthermore, a large amount of data must be stored unless checkpointing is used, which slows the algorithm. Thus, it is difficult to construct a (near-) optimal control in real time if the problem is nonlinear and multidimensional, as in our case, and so it may be unfeasible to use such controls in applications. However, the computed optimal controls indicate the efficacy of proportional controls, where ζ = −α(η − η); a study of (point-actuated) proportional controls by the authors will be presented elsewhere.
In section 5, we studied the effects of purely spanwise forcing in a nonoptimal setting. We focused on the vertical falling film case, for which controls are not required to ensure bounded solutions. With extensive numerical simulations, we found that small amplitude sinusoidal transverse waves yielded a decrease in the average energy of the fluid interface, with large amplitude waves being nonlinearly stable, regularizing the chaotic dynamics.
A large motivation for this paper and current work by the authors is the construction of controls for other systems in the hierarchy of models, with the ultimate goal of using control methodologies derived for the KSE to control the Navier-Stokes equations-from inducing recirculation regions (to improve heat transfer) to stabilizing the exact Nusselt solution and preventing dripping for hanging film arrangements. This will also serve as a test for the weakly nonlinear models in describing the dynamics of the full system. the general desired state (53) η = (I −P Σ )η +ψ(y, t),ψ(y, t) = k 2 ∈Σψ k 2 (t)e ik 2 y .
Here,P Σ is a projection onto a subset of the transverse Fourier modes, andψ is a real-valued function of the modes in Σ alone, with the property thatψ −k 2 is the complex conjugate of ψ k 2 ; importantly then, Σ must satisfy the property that if k 2 ∈ Σ, then −k 2 ∈ Σ. The term "tracker" is used to indicate that the desired state is possibly nontrivial and time-dependent. In section 3, this theory was employed for Σ = Ξ andψ = 0, and the full generalization gives a viable choice of control to be used in conjunction with the results in section 5. The linearity of the problem allows us to build optimal controls for each mode. From the k 2 th component of the cost functional, we have the Hamiltonian (54) H k 2 = |k 2 | 2s η k 2 (t) −ψ k 2 (t) 2 + γ|ζ k 2 (t)| 2 + p k 2 −(κk 2 2 +k 4 2 )η k 2 +ζ k 2 , where p k 2 is the adjoint variable for the k 2 th mode (with complex conjugate p −k 2 ). Then, from Hamilton's equations, we have a two point boundary value problem, where the boundary conditions are (58)η k 2 (0) = v (0,k 2 ) , p k 2 (T ) = |k 2 | 2s η −k 2 (T ) −ψ −k 2 (T ) .
The final time boundary condition for the adjoint is found by differentiating the cost functional with respect toη k 2 (T ). Taking the complex conjugate of (56) and (57) gives the most useful form of the equations. To solve this two point boundary value problem, we make the ansatz that (59) p −k 2 (t) = −γr k 2 (t)η k 2 (t) + q k 2 (t).