Mohamadreza Ahmadi

Giorgio Valmorbida

Dennice Gayme

Antonis Papachristodoulou

A Framework for Input-Output Analysis of Wall-Bounded Shear Flows

Keywords: Control theory, Navier-Stokes equations, Transition to turbulence, Nonlinear instability, Channel flow

We propose a new framework to evaluate input-output amplification properties of nonlinear models of wall-bounded shear flows, subject to both square integrable and persistent disturbances. We focus on flows that are spatially invariant in one direction and whose base flow can be described by a polynomial, e.g. streamwise constant channel, Couette and pipe flows. Our methodology is based on the notion of dissipation inequalities in control theory and provides a single unified approach to examining flow properties such as energy growth, worst case disturbance amplification, and stability to persistent excitation (i.e., input-to-state stability). It also enables direct analysis of the nonlinear partial differential equation (PDE) rather than of a discretized form of the equations, thereby removing the possibility of truncation errors. We demonstrate how to numerically compute the input-output properties of the flow as the solution of a (convex) optimization problem. We apply our theoretical and computational tools to plane Couette, channel and pipe flows. Our results demonstrate that the proposed framework leads to results that are consistent with theoretical and experimental amplification scalings obtained in the literature.

Introduction

Literature Review

In the study of hydrodynamic stability, one aims to find out if a given flow state is stable or unstable, and if so, how these instabilities lead to turbulence. Conventional hydrodynamic stability methods involves linearization of the Navier-Stokes equations around a base flow followed by spectral analysis, which aims to indentify the Reynolds number Re at which this solution becomes unstable [START_REF] Drazin | Hydrodynamic Stability[END_REF]). However, the gap between the linear stability limits and the ones observed in experiments is often significant, with flows routinely transitioning to turbulence well below the critical Reynolds number identified in the analysis (a phenomenon refered to as sub-critical transition). This discrepancy has long been attributed to the use of the linearized Navier-Stokes operator [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF]). Other theoretical methods for studying stability of flows are often based on spectral truncation of the Navier-Stokes equations into a system of ordinary differential equations (ODEs). However, this approach can lead to truncation errors thereby creating a mismatch between the dynamics of the truncated model and the partial differential Navier-Stokes equations. Both of these drawbacks were addressed in [START_REF] Goulart | Global stability analysis of fluid flows using sumof-squares[END_REF]; [START_REF] Chernyshenko | Polynomial sum of squares in fluid dynamics: a review with a look ahead[END_REF]), where the authors introduced a method involving keeping a number of modes from the Galerkin expansion of the nonlinear Navier-Stokes equations and bounding the energy of the remaining modes, leading to a nonlinear ODE description of the flow dynamics. The resulting ODE description allowed for techniques from control theory, such as Lyapunov's direct method [START_REF] Khalil | Noninear Systems[END_REF]), to be used to study stability and was implemented computationally based on sum-of-squares programming [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF]). It was shown in [START_REF] Huang | Sum-of-squares of polynomials approach to nonlinear stability of fluid flows: an example of application[END_REF]) that, in the case of rotating Couette flow, this method can find a global stability limit, which is better than the energy method but not as good as the linear stability limit †.

Even in the seminal paper by [START_REF] Reynolds | An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous and the law of resistance in parallel channels[END_REF], it was observed that external excitations and body forces play an important role in flow instabilities. Mechanisms such as energy amplification of external excitations and body forcings have shown to be crucial in understanding transition to turbulence, as highlighted by [START_REF] Joseph | Stability of fluid motions[END_REF]. Therefore, instead of studying stability, researchers began to focus on energy growth and were able to uncover additional flow properties through the new paradigm of input-output analysis. This led to a focus on so-called transient growth as a primary factor for sub-critical transition to turbulence; i.e., although the perturbations to the linearized Navier-Stokes equation are stable (the eigenvalues have negative real parts), they undergo high amplitude transient amplifications that steer the trajectories out of the region of linearization. The root cause of the transient growth phenomenon is the non-normality of the stable Navier-Stokes operator that has been linearized about a base flow. The focus on energy growth and the associated input-output properties of the flow has led to studying the resolvent operator or ε-pseudospectra to uncover when transition occurs, based on the general solution to the linearized Navier-Stokes equations [START_REF] Schmid | Nonmodal stability theory[END_REF]). In particular, (McKeon & Sharma (2010)) used resolvent analysis to study the amplification scalings from an input composed of nonlinear terms and periodic forcings for turbulent pipe flows.

The input-output properties can be characterized based on the class of forcings (stochastic noise vs square integrable signals) and the flow model (linear vs nonlinear and finite-dimensional vs infinite dimensional) one considers. For stochastic forcings (Gaussian noise), energy amplification to the linearized Navier-Stokes equations in wallbounded shear flows was studied by [START_REF] Farrell | Stochastic forcing of the linearized Navier-Stokes equations[END_REF]. In a similar vein, [START_REF] Bamieh | Energy amplification in channel flows with stochastic excitation[END_REF] used the stochastically forced linearized Navier-Stokes equation to show that the input-output energy amplification from streamwise constant excitations to perturbation velocities in channel flows is proportional to Re 3 . The amplification scaling of the linearized Navier-Stokes equation was further characterized in [START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF]) and [START_REF] Jovanović | Modeling, analysis, and control of spatially distributed systems[END_REF]), where the authors studied the influence of each component of the body forces. For square integrable forcings, (Jovanović 2004, Chapter 9) and [START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF]) provided worst-case amplification scalings for incompressible viscous channel flows based on the linearized Navier-Stokes equations. However, an input-output analysis for the nonlinear (PDE) flow model is absent from the literature.

Contribution

Our work extends the input-output analysis paradigm. Specifically, we propose a method based on dissipation inequalities [START_REF] Willems | Dissipative dynamical systems part I: General theory[END_REF]), to study input-output amplification in wall-bounded shear flows. This method extends previous results in the following ways. First, it applies directly to nonlinear PDE models of the flow rather than linearized [START_REF] Vazquez | Control of turbulent and magnetohydrodynamic channel flows: boundary stabilization and state estimation[END_REF]) or nonlinear equations that have been discretized and projected onto finite bases. We extend the class of forcings studied in previous work from square integrable forcings to persistent but bounded disturbances, which allows us to analyze both traditional notions of input-output flow properties (e.g. energy amplification and disturbance amplification) as well as new stability properties such as the control theoretic notion of input-to-state stability. The latter refers to the stability of the flow field under persistent deterministic excitations, which has not been previously analyzed. Our technique has the added advantage of enabling the study of multiple input-output flow properties for a broad class of wall-bounded shear flows within a single framework.

One important observation from previous work on input-output analysis is the importance of streamwise-constant perturbations, see e.g. Gustavsson (1991a); [START_REF] Farrell | Stochastic forcing of the linearized Navier-Stokes equations[END_REF]; [START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF]. The importance of streamwise constant structures is also supported with experimental observations of the ubiquity of streamwise coherent structures, see e.g. [START_REF] Kim | Very large scale motion in the outer layer[END_REF]; [START_REF] Kline | The structure of turbulent boundary layers[END_REF]; [START_REF] Bullock | Structural similarity in radial correlations and spectra of longitudinal velocity fluctuations in pipe flow[END_REF]; [START_REF] Hutchins | Evidence of very long meandering features in the logarithmic region of turbulent boundary layers[END_REF] and analysis showing that streamwise constant models reproduce important flow properties [START_REF] Gayme | A streamwise constant model of turbulence in plane[END_REF][START_REF] Gayme | Amplification and nonlinear mechanisms in plane Couette flow[END_REF]). We therefore focus on this class of problems.

We complement the analysis framework by providing a computationally tractable method to evaluate input-output properties by solving a convex optimization [START_REF] Boyd | Convex optimization[END_REF]) problem. Convex optimization problems can be solved efficiently by interior-point methods [START_REF] Nesterov | Interior-point Polynomial Algorithms in Convex Programming[END_REF]) and there are several tools available for solving them (see e.g. [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF]). The strength of the method is that the results can be directly extended to more complex flow geometries as long as they can be described or approximated [START_REF] Weihs | On the polynomial approximation of boundary-layer flow profiles[END_REF]; [START_REF] Rubin | Polynomial interpolation methods for viscous flow calculations[END_REF]) by semi-algebraic sets, i.e., sets characterized by a number of polynomial equalities and inequalities.

We apply then the theory and computational tools to compute the input-output properties of streamwise constant plane Couette, plane Poiseuille, and Hagen-Poiseuille flows. Our results lead to upper-bounds that are consistent with the transient growth results in the literature. For channel flows, we show the upper-bounds obtained using our method agree with the results in [START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF]) and (Jovanović 2004, Chapter 9) in terms of worst-case disturbance amplification and we show that our framework can be used to study pipe flows as well. Moreover, we observe that the stability bounds to persistent forcings and the experimental Reynolds numbers for transition to turbulence are very similar.

Preliminary mathematical results related to this work were presented in [START_REF] Ahmadi | A convex approach to hydrodynamic analysis[END_REF]). The current paper is different from [START_REF] Ahmadi | A convex approach to hydrodynamic analysis[END_REF]) in several aspects. From a theoretical standpoint, the current paper provides a method for energy growth analysis and extends the formulation to both flows between parallel plates and flows in pipes. In addition, it presents the mathematical proofs of the input-output analysis framework and the computational formulation based on convex optimization. From the examples standpoint, we applied the framework to investigate the input-output properties of plane Couette flow, plane Poiseuille, and the Hagen-Poiseuille flow. Furthermore, the current paper includes a detailed comparison with previous results in the literature and an examination of flow structures corresponding to maximum input-output amplifications.

Organization

In the next section, we briefly describe the flow model studied in the paper. In Section 3, we propose the flow input-output analysis framework based on dissipation inequalities. In Section 4, we show how the input-output analysis can be computationally implemented as the solution to a convex optimization problem. In Section 5, we demonstrate the effectiveness of the proposed framework by applying it to study input-output properties of plane Couette flow, plane Poiseuille flow, and Hagen-Poiseuille flow. Finally, in Section 6, we present some concluding remarks and provide directions for future research.

The Flow Perturbation Model

We consider the dynamics of forced incompressible shear flows in a given domain with time-independent boundary conditions. The instantaneous flow velocity ū(t, x), ū : R 0 × Ω → R 3 , and pressure p(t, x), p : R 0 × Ω → R, are governed by the Navier-Stokes equations

∂ t ū = 1 Re ∇ 2 ū -ū • ∇ū -∇p + d, 0 = ∇ • ū, (2.1)
where t > 0, x ∈ Ω = Ω i × Ω j ⊂ R × R with i = j, i, j ∈ I are spatial coordinates with I being the index set of the coordinates and ∂ s (•) = ∂(•) ∂s . The dependent variable d : R 0 × Ω → R 3 is the input vector representing exogenous excitations or body forces. ∇ 2 is the Laplacian operator, ∇ denotes the gradient, and ∇ • u denotes the divergence of u.

We consider perturbations (u, p) to a laminar/nominal solution (U , P ), which are spatially invariant in one of the directions, for example x m , m ∈ I, i.e., ∂ xm = 0. This leads to a 2 dimensional, 3 component (2D/3C) flow field. Let I 0 = I -{m}. The velocity field can be decomposed as

ū = u + U , p = p + P, (2.2)
where U and P are steady state solutions to (2.1), i.e.,

0 = 1 Re ∇ 2 U -U • ∇U -∇P, 0 = ∇ • U . (2.3) Substituting (2.
2) in (2.1) and using (2.3), we obtain the perturbation dynamics

∂ t u = 1 Re ∇ 2 u -u • ∇u -U • ∇u -u • ∇U -∇p + d, 0 = ∇ • u. (2.4)
In the rest of this paper, we study the properties of (2.4). We concentrate on perturbations with no-slip boundary conditions u| ∂Ω ≡ 0 (in the direction with solid boundaries)

and periodic boundary conditions (in the spatially homogeneous direction). In a similar manner, we extend the results to pipe flows (cylindrical coordinates) as discussed in Appendix C. Throughout the paper, it is assumed that the perturbation dynamics (2.4) admit a solution that is sufficiently differentiable and square integrable. Next, we introduce the input-output analysis methods based on dissipativity theory.

Dissipation Inequalities for Shear Flows

In systems and control theory, dissipativity [START_REF] Willems | Dissipative dynamical systems part I: General theory[END_REF][START_REF] Willems | Dissipative dynamical systems[END_REF]; [START_REF] Hill | Dissipative dynamical systems: basic input-output and state properties[END_REF]) † establishes a relationship between the energy stored in the system represented by a continuous, non-negative functional V (u), known as the storage functional, and the power supplied to the system W (u, d, y), known as the supply rate, with d and y being the inputs and outputs of the system, respectively. In particular, in this study, we are concerned with y = u. This relationship is often given by a dissipation inequality (in differential form) as

dV (u) dt W (u, d, y). (3.1)
A system is called dissipative with respect to the supply rate W (u, d, y), if there is a nonnegative functional V (u) that satisfies (3.1). Dissipativity theory has a close connection with Lyapunov stability theory [START_REF] Khalil | Noninear Systems[END_REF]). In particular, dissipativity theory can be understood as a generalization of the Lyapunov stability theory to systems with inputs and outputs.

Given the dissipation inequality (3.1) with a fixed supply rate, the main challenge is to find a corresponding storage functional that satisfies the dissipation inequality along the solutions of the flow. In fact, kinetic energy was shown to be a candidate storage functional for some input-output properties. In the special case of an irrotational flow under no-slip, stress-free or periodic boundary conditions, if we set V to be the kinetic energy of the perturbations V (u) = 1

2 Ω |u| 2 dΩ and W (u, d) = Ω u • d dΩ as the supply rate, we can show that dissipation inequality (3.1) holds. Following (Doering & Gibbon 1995, p. 31), the total kinetic energy of the perturbations satisfies the following equality

dV (u) dt = - 1 Re ∇u 2 L 2 Ω - Ω u • ∇U • u dΩ + Ω u • d dΩ.
The above equality implies that the kinetic energy of the perturbations in the flow changes according to three effects: the energy dissipated by viscosity, the energy either injected or dissipated depending on the base flow, and the energy expended by the external force.

Since the viscosity term 1 Re ∇u 2 L 2 Ω is always non-negative, we can obtain the following inequality

dV (u) dt - Ω u • ∇U • u dΩ + Ω u • d dΩ.
If the base flow U is such that the term Ω u•∇U •u dΩ is non-negative, we can conclude † Note that the notion of dissipativity used here should not be confused with dissipative operators in semigroup theory [START_REF] Lumer | Dissipative operators in a banach space[END_REF]). The latter is concerned with proving the existence of a contraction semigroups and used to prove well-posedness of solutions to PDEs, [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF]; whereas, the dissipativity notion we use here is concerned with the input-output properties of a dynamical system.

that the following dissipation inequality holds

dV (u) dt Ω u • d dΩ.
This is a well-known dissipation inequality that corresponds to passivity. Passivity has been used to study finite-dimensional linear discretizations of the Navier-Stokes equation with the nonlinearity u•∇u being modeled as an input [START_REF] Sharma | Relaminarisation of reτ = 100 channel flow with globally stabilising linear feedback control[END_REF][START_REF] Heins | Passivity-based output-feedback control of turbulent channel flow[END_REF]).

The general dissipation inequality framework allows us to consider more general energy inequalities rather than only the passivity inequality. In particular, our formulation considers weighted kinetic energy as the storage functional and three different supply rates. As will be shown in Section 4, we present an algorithmic way to find the storage functionals based on convex optimization for the class of fluid flows studied in this paper.

Input-Output Properties

We now define the three types of input-output properties that we study within the dissipativity framework and discuss their relation to common notions in the literature.

The first property that we study is the maximum energy growth due to initial perturbation velocities that are square integrable over the spatial domain for the nonlinear Navier-Stokes equation (2.4). In the context of linear systems, this property corresponds to maximum transient growth [START_REF] Butler | Three-dimensional optimal perturbations in viscous shear flow[END_REF]; [START_REF] Reddy | Energy growth in viscous channel flows[END_REF]; Gustavsson (1991b)). For a vector l, we denote by l the transpose of l. Definition 3.1 (Energy Growth). For any u satisfying (2.4) with d ≡ 0and initial flow perturbations that are square integrable over the spatial domain, i.e., u 0

L 2 Ω = Ω u 0 (θ)u 0 (θ) dθ 1 2 < ∞, if there exists a constant γ > 0 such that u L 2 [0,∞),Ω γ u(0, •) L 2 Ω , (3.2) where u L 2 [0,T ),Ω = T 0 Ω u (τ, θ)u(τ, θ) dτ dθ 1 2
, then we say that the flow perturbations have bounded energy growth with bound γ.

The next property of interest is related to amplifications from square integrable body forces or disturbances (see (Jovanović 2004, Chapter 9) for analogous results pertaining to a linearized model of channel flows). The square integrable forcings are of special interest, because they can be interpreted as finite energy forcings. We refer to this class of amplifications as worst-case disturbance amplification. Definition 3.2 (Worst-Case Disturbance Amplification). If there exist

η i > 0, i ∈ I such that u 2 L 2 [0,∞),Ω i∈I η 2 i d i 2 L 2 [0,∞),Ω , (3.3)
for all u satisfying (2.4) subject to zero initial perturbations u(0, x) ≡ 0, ∀x ∈ Ω and square integrable disturbances over both spatial domain and time, i.e.,

d i 2 L 2 [0,∞),Ω < ∞,
i ∈ I, then we say that the flow has bounded worst-case disturbance amplification with bounds η i , i ∈ I.

The above property is equivalent to the induced L 2 -norm in control theory (Van der Schaft (2017)). In other words, each η i upper-bounds the peak amplification of perturbation velocities from the forcing in the direction i, d i , when the forcings in other directions are set to zero, i.e., d j = 0, j ∈ I, i = j. That is,

η i sup di L 2 [0,∞),Ω =0 u L 2 [0,∞),Ω d i L 2 [0,∞),Ω
.

Note that the above relation provides an upper-bound for the amplification of disturbances of any magnitude. For nonlinear systems, such as the flow dynamics studied here, the amplification may depend on the magnitude of the disturbance d L 2 in a nonlinear fashion (Ahmadi et al. 2016, Example I). In this case, any value of η satisfying (3.3) provides an upper bound to the worst-case disturbance amplification. Therefore, in order to obtain tight bounds to the worst-case amplification, we are interested in the minimum values of η i , i ∈ I for which (3.3) holds.

From a practical perspective, stability of a base flow is often not very meaningful, because small disturbances may cause an unstable behavior. We instead require a notion of stability that relates disturbances to perturbation velocities. In particular, since the definition of the worst-case disturbance amplification requires the forcings to be square integrable, it excludes persistent forcings, e.g. constant and sinusoidal forcings, that are defined for all time. To include these classes of forcings in a nonlinear context †, we employ the concept of input-to-state stability [START_REF] Sontag | Input to state stability: Basic concepts and results[END_REF]) to study the class of upperbounded forcings. We refer to this extended notion of stability, as stability to persistent disturbances. Prominent among the features of this property are that forcings that are bounded, eventually small, integrally small, or convergent lead to perturbation velocities with the respective property. Furthermore, this property quantifies in what manner initial perturbation velocities affect transient behavior. Flows with this property do not have unstable behavior for persistent (nonvanishing) forcings.

To characterize this property, let us introduce a few comparison functions. Let K denote the class of non-negative functions that are strictly increasing and zero for zero argument, and K ∞ denote the class of functions that, in addition, become unbounded as their argument goes to infinity. For a vector l denote by |l|, its absolute value given by √ l l.

Definition 3.3 (Stability to persistent disturbances). For all u satisfying (2.4) with upper-bounded disturbances, i.e., d

L ∞ [0,∞) = sup τ ∈[0,∞) |d(τ )| < ∞, if there exist some scalar ψ > 0, functions β, β, χ ∈ K ∞ , and σ ∈ K, such that u(t, •) L 2 Ω β e -ψt χ u(0, •) L 2 Ω + β sup τ ∈[0,t) Ω σ |d(τ, x)| dΩ , (3.4)
for all t > 0, then we call the flow stable to persistent disturbances.

Property (3.4) implies convergence to the base flow (U , P ) in the L 2 Ω -norm (the norm corresponding to the space of square integrable functions over the spatial domain) when the disturbances are not present (d ≡ 0). Indeed, the β e -ψt χ u(0, •) L 2 Ω term dominates for small t, and this serves to quantify the magnitude of the transient growth as a function of the size of the initial state u(0,

•) L 2 Ω . t u L 2 Ω ≈ β Ω σ( d L ∞ ) dΩ ≈ β • χ u 0 L 2 Ω Figure 1:
The stability to persistent disturbances property combines transient growth (overshoot) and asymptotic behavior.

Moreover, as t → ∞, we obtain

lim t→∞ u(t, •) L 2 Ω β Ω σ(|d(•, x)|) L ∞ [0,∞) dΩ β Ω σ( d(•, x) L ∞ [0,∞) ) dΩ , (3.5)
where, σ, ) is bounded by assumption, σ ∈ K (continuous and bounded), and the integral is over a bounded domain Ω, the quantity on the right-hand side of (3.5) is also bounded. Hence, as long as the external excitations or body forces d are upper-bounded, the perturbation velocities u are bounded in the L 2 Ω -norm, meaning that they remain square integrable over the flow geometry.

β ∈ K. Since d(•, x) L ∞ [0,∞
In fact, by input-to-state superposition theorem (Sontag ( 2013)), we can shows that stability to persistent disturbances is the conjunction of two properties, one of them concerned with asymptotic bounds on the perturbation velocities, in the sense of u(t, •) L 2 Ω , as a function of the magnitude of the forcings, and the other one providing a transient term obtained when we ignore forcings (see Figure 1).

We now demonstrate how the problem of verifying the properties in Definitions 3.1-3.3 can be cast as verifying a set of dissipation inequalities. These results, which can be derived from (Ahmadi et al. 2016, Theorem 6), allows for the extension of well known methods for stability, input/output, and optimal perturbation analysis of linear systems to the full nonlinear Navier-Stokes equations.

We begin by studying energy growth.

Theorem 3.4. Consider the perturbation model (2.4) with d ≡ 0 and let

u(0, •) 2 L 2 Ω < ∞, i.e.
, the initial perturbations be square integrable over the spatial domain. If there exist a positive semidefinite storage functional V (u) and a positive scalar γ such that

V (u) γ 2 u(t, •) 2 L 2 Ω , (3.6) dV (u(t, x)) dt - Ω u (t, x)u(t, x) dΩ, (3.7)
then it has bounded energy growth as given by (3.2).

The next result pertains to worst-case energy amplification.

Theorem 3.5. Consider the perturbation model (2.4) and let

d(t, •) 2 L 2 Ω < ∞, i.e
., the disturbances be square integrable over the spatial domain. If there exist a positive semidefinite storage functional V (u) and positive scalars {η i } i∈I such that

dV (u(t, x)) dt - Ω u (t, x)u(t, x) dΩ + Ω i∈I η 2 i d 2 i (t, x) dΩ, (3.8)
then the perturbation velocities (2.4) have worst-case disturbance amplification upperbounds η i , i ∈ I as in (3.3).

Lastly, we study amplifications to persistent forcings.

Theorem 3.6. Consider the perturbation model (2.4) and let

d L ∞ [0,∞) < ∞, i.e.
, the disturbances have an upper bound. If there exist a positive semidefinite storage functional V (u), a positive scalar ψ, and functions

β 1 , β 2 , σ ∈ K ∞ such that β 1 ( u(t, •) L 2 Ω ) V (u) β 2 ( u(t, •) L 2 Ω ), (3.9) dV (u(t, x)) dt -ψV (u(t, x)) + Ω σ(|d(t, x)|) dΩ, (3.10)
then the perturbation velocities described by (2.4) are stable to persistent disturbances as given by (3.4)

with χ = β 2 , β = β -1 1 • 2 and β = β -1 1 • 2 ψ , where • implies function composition.
In the following, we derive classes of storage functionals V (u) suitable for the analysis of perturbation dynamics (2.4) invariant in one of the three spatial coordinates. We consider two classes of flows, namely, channel flows with perturbations that vary in two spatial dimensions and time discussed in Section 3.2 and pipe flows invariant in the axial direction discussed in Appendix C.

A Lyapunov/Storage Functional for Flows Between Parallel Plates

In Cartesian coordinates, for a scalar function v, ∇v =

i ∂ xi v - → e i and ∇ 2 v = i ∂ 2 xi v,
where -→ e i is the unit vector in the direction x i . For a vector valued function

w = i w i - → e i , and the divergence ∇ • w is given by ∇ • w = i ∂ xi w i .
In the following, {x 1 , x 2 , x 3 } corresponds to {x, y, z} (the streamwise, wall-normal, and spanwise directions) and I = {1, 2, 3}. Additionally, we adopt Einstein's multi-index notation over index j, that is the sum over repeated indices j, e.g.,

v j ∂ xj u i = j v j ∂ xj u i .
The perturbation model (2.4) can be re-written as (3.11) where i, j ∈ I. To simplify the exposition, without loss of generality, we assume that the perturbations are invariant with respect to x 1 . Since x i , i = 1, 2, 3 are arbitrary, this does not affect the formulation.

∂ t u i = 1 Re ∇ 2 u i -u j ∂ xj u i -U j ∂ xj u i -u j ∂ xj U i -∂ xi p + d i , 0 = ∂ xj u j .
The next proposition states that, by choosing a suitable storage functional structure (weighted kinetic energy of the perturbation velocities), the time derivative of the storage functional turns out to be upper-bounded by a quadratic form in the velocity fields u and their spatial derivatives. This property paves the way for a convex optimization based method to check stability and input-output properties. Convex optimization is a subfield of optimization that studies the problem of minimizing convex functions over convex sets. The convexity makes optimization easier than the general case since local minimum must be a global minimum, and first-order conditions are sufficient conditions for optimality [START_REF] Boyd | Convex optimization[END_REF]). Convex optimization problems can be solved efficiently by interior-point methods [START_REF] Nesterov | Interior-point Polynomial Algorithms in Convex Programming[END_REF]). Convex optimization was used by [START_REF] Moarref | A low-order decomposition of turbulent channel flow via resolvent analysis and convex optimization[END_REF] to obtain a low-order decomposition of the Navier-Stokes equations based on resolvent modes.

Proposition 3.7. Consider the perturbation model (3.11) subject to periodic or noslip boundary conditions u| ∂Ω = 0. Assume the velocity perturbations in (3.11) are invariant with respect to x 1 . Let I 0 = {2, 3} and

V (u) = 1 2 Ω u Qu dΩ, (3.12)
where Q = q1 0 0 0 qj 0 0 0 qi > 0, q i = q j for i = j, i, j ∈ I 0 , be a candidate storage functional.

Then, the time derivative of (3.12) along the solutions to (3.11) satisfies

dV (u) dt - i∈I q i Ω C Re u 2 i + U j u i ∂ xj u i + u j u i ∂ xj U i -u i d i dΩ, (3.13)
where C is a positive constant that only depends on the domain Ω.

The proof of this proposition is given in Appendix A.

Remark that a special case of (3.12) was used in [START_REF] Joseph | Contributions to the nonlinear theory of stability of viscous flow in pipes and between rotating cylinders[END_REF]) to study the stability of viscous flows (subject to streamwise constant perturbations) in pipes and between rotating cylinders. The authors referred to this structure as the two energy function. In the formulation presented in this paper, assuming invariant perturbations in the x 1 -direction, and since the inequalities hold if scaled by a positive scalar, we can represent the two energy function as

V (u) = 1 2 Ω u q 0 0 0 1 0 0 0 1 u dΩ,
where q is a positive scalar. The "optimal" value for this constant was then calculated analytically for the pipe Poiseuille and the Taylor-Couette flow by [START_REF] Joseph | Contributions to the nonlinear theory of stability of viscous flow in pipes and between rotating cylinders[END_REF]. Indeed, considering streamwise-constant perturbations provides a degree of freedom in the Lyapunov/storage functional, which can be used to study flow properties at higher Reynolds numbers. When the perturbations are spatially varying in all three directions, the Lyapunov/storage functional converts to the energy functional, which has been used to study the stability of viscous flows [START_REF] Serrin | On the stability of viscous fluid motions[END_REF]). However, since for 3D perturbations there is no degree of freedom in the structure of the Lyapunov/storage functional, the flow properties can only be studied for low Reynolds numbers.

Note that in (3.13) the Poincaré constant, C, appears. There are several estimates for the optimal Poincaré constant. The optimal constant [START_REF] Payne | An optimal Poincare inequality for convex domains[END_REF]) we use in this paper is

C(Ω) = π 2 D(Ω) , (3.14)
where D(Ω) is the diameter of the domain Ω (the longest distance between any two point in the domain Ω). Proposition 3.7 allows us to provide an algorithmic method for input-output analysis of fluid flows based on convex optimization. These convex optimization problems are in terms of linear matrix inequalities and polynomial matrix inequalities. This formulation is delineated in more detail in the next section.

Matrix Inequalities for Input-Output Analysis of Streamwise Constant Perturbations

In this section, we show that the input-output analysis problem outlined in Section 3 for the class of streamwise constant perturbations can be converted into a set of matrix inequalities. These matrix inequalities can be solved by convex optimization, provided that the base flow is a polynomial in the spatial coordinates and the flow geometry is described by a semi-algebraic set †. Examples are laminar base flows that are linear or parabolic, and turbulent flows that can be represented by polynomial fits ‡ (or by piecewise polynomial functions).

To present a convex method for checking the conditions in Theorems 3.4-3.6 (see also Corollary B.1 in Appendix B), we restrict our attention to streamwise constant perturbations in the x 1 -direction with base flow U = U m (x 2 , x 3 ) -→ e 1 , where -→ e 1 denotes the unit vector in the x 1 -direction.

In order to present the procedure, we first need to define the following notation. For a square matrix M , M 0 (M 0) implies that the matrix is positive semidefinite (positive definite), i.e., all the eigenvalues of M are non-negative (positive). Similarly, M 0 (M ≺ 0) signifies that -M 0 (-M 0). By I n×n , we denote the square matrix of dimension n × n with diagonal entries set to 1.

Corollary 4.1. Consider the perturbation dynamics given by (3.11), that are constant in the streamwise direction x 1 and with base flow U = U m (x) -→ e 1 , where

x = (x 2 , x 3 ). Let I 0 = {2, 3}. If there exist positive constants {q l } l∈I with q i = q j , i, j ∈ I 0 , {η l } l∈I , ψ, and functions {σ l } l∈I such that

M (x) =    C Re q 1 q1∂x j Um(x) 2 q1∂x i Um(x) 2 q1∂x j Um(x) 2 C Re q j 0 q1∂x i Um(x) 2 0 C Re q i    , i, j ∈ I 0 , i = j. (4.1) I) when d ≡ 0, M (x) -I 3×3 0, x ∈ Ω, (4.2) † Let R[x]
be the set of polynomials with real coefficients. A set is semi-algebraic if it can be described by a finite number of polynomial equalities and inequalities. That is, S ⊂ R n for some closed field, say R, is defined by a set of polynomial equalities and inequalities as follows

S = {x ∈ R n | pi(x) 0, i = 1, 2, . . . , np, qi(x) = 0, i = 1, 2, . . . , nq} , where {pi} np i=1 , {qi} nq i=1 ∈ R[x],
where R[x] denotes the set of polynomials in the variable x with real coefficients.

‡ Note that by the Weierstrass approximation theorem [START_REF] De Branges | The Stone-Weierstrass theorem[END_REF]), any continuous function defined on a closed interval can be uniformly approximated as closely as desired by a polynomial function.

II)

N (x) =              -q1 2 0 0 M (x) -I 3×3 0 - qj 2 0 0 0 -qi 2 -q1 2 0 0 η 2 1 0 0 0 - qj 2 0 0 η 2 i 0 0 0 -qi 2 0 0 η 2 j              0, (4.3) for i, j ∈ I 0 , i = j and x ∈ Ω, III) σ l (x) 0, x ∈ Ω, l ∈ I and Z(x) =              -q1 2 0 0 M (x) -W 0 - qj 2 0 0 0 -qi 2 -q1 2 0 0 σ 1 (x) 0 0 0 - qj 2 0 0 σ j (x) 0 0 0 -qi 2 0 0 σ i (x)              0, (4.4)
for i, j ∈ I 0 , i = j and x ∈ Ω, where W = ψq1 0 0 0 ψqj 0 0 0 ψqi . Then, it follows that I) the flow energy growth is bounded by γ 2 = max i∈I q i as described by (3.2), II) the worst-case disturbance amplification (induced L 2 norm from the disturbances to perturbation velocities) is bounded by η i , i ∈ I as in (3.3) when the initial perturbations have zero velocity, III) the flow is stable to persistent disturbances in the sense of (3.4) with σ(|d|) =

i∈I σ i (x)d 2 i .
The proof of the above Corollary is given in Appendix B.

When U m (x) is a polynomial function, inequalities (4.1)-(4.4) are polynomial matrix inequalities that should be checked for all x ∈ Ω. If the set Ω is a semi-algebraic set, i.e.,

Ω = x ∈ R 2 | g l (x) = 0, f k (x) > 0, l = 1, 2, . . . , L, k = 1, 2, . . . , K ,
where {g l } L l=1 and {f k } K k=1 are polynomial functions, then these inequalities can be cast as a sum-of-squares program by applying Corollary D.4. We show in the next section that this assumption is indeed the case for several well-known flows. For a brief introduction to sum-of-squares programming refer to Appendix D. Note that once the input-output analysis problem is cast as a sum-of-squares program, it can be checked using available MATLAB toolboxes such as SOSTOOLS [START_REF] Papachristodoulou | SOSTOOLS: Sum of squares optimization toolbox for MATLAB[END_REF]) and YALMIP [START_REF] Löfberg | YALMIP : A toolbox for modeling and optimization in MATLAB[END_REF]).

We can compute the bound on the maximum energy grown described in (3.2) by solving an optimization problem. To this end, we solve minimize

{qi} i∈I max i∈I q i subject to M (x) -I 3×3 0, ∀x ∈ Ω, q i > 0, i ∈ I.
(4.5)

In order to find upper-bounds on the worst-case disturbance amplification (the induced L 2 -norm) from the body forces or disturbances d to the perturbation velocities u as described in (3.3), we solve the following optimization problem minimize

{qi} i∈I ,{ηi} i∈I i∈I η 2 i subject to N (x) 0, ∀x ∈ Ω, η i > 0, i ∈ I, q i > 0, i ∈ I.
(4.6)

Numerical Results

In this section, we apply the proposed method to three cannonical flows, namely, plane Couette flow, plane Poiseuille flow, and Hagen-Poiseuille flow. For each case we compute bounds on energy growth, worst-case disturbance amplification, and stability to persistent forcings (i.e., input-to-state stability). Our analysis of the worst-case disturbance amplification includes a comparative analysis of the influence of each of the disturbance components. For the case of stability to persistent disturbances, we find the maximum Reynolds number for which stability to persistent disturbances holds.

For each result, we provide explicit descriptions of the matrix inequalities that need to be solved to study each of the three input-output properties based on Theorems 3.4-3.6. These matrix inequalities can then be solved using available MATLAB toolboxes such as SOSTOOLS [START_REF] Papachristodoulou | SOSTOOLS: Sum of squares optimization toolbox for MATLAB[END_REF]) and YALMIP [START_REF] Löfberg | YALMIP : A toolbox for modeling and optimization in MATLAB[END_REF] as described in Appendix D.

Plane Couette Flow

We consider the flow of viscous fluid between two parallel plates, where the gap between the plates is much smaller than the length of the plates as illustrated in Figure 2. The base flow is given by U = (y, 0, 0) = y -→ e x and P = P 0 .

We consider no-slip boundary conditions u| 1 y=-1 = 0 in the wall-normal direction and u(t, y, z) = u(t, y, z + L) in the spanwise direction. The Poincaré constant is then given by C

= π 2 √ L 2 +2 2 .
A storage functional of the form (3.12), can be used to compute bounds on all three quantities of interest (i.e., energy growth, worst-case amplification, and stability to persistent forcings) with q y = q z . We assume that spatial invariance in the x-direction, and thus (m = x, j = y, i = z). The M matrix in (4.1) is then given by

V (u) =
M =   qxC Re qx 2 0 qx 2 qyC Re 0 0 0 qyC Re   (5.2)
Given a channel width L = π we can compute the input-output energy growth by solving the optimization problem (4.5) with M given by (5.2). The resulting energy growth γ 2 as a function of Reynolds number is plotted in Figure 3. For small Reynolds numbers γ 2 ∼ O(Re), whereas for larger Reynolds numbers γ 2 ∼ O(Re 3 ). Therefore, it can be inferred that γ 2 ∼ c 0 Re + c 1 Re 3 with c 0 , c 1 > 0. This scaling is consistent with the results obtained by [START_REF] Bobba | Highly optimized transitions to turbulence[END_REF], where the maximum energy growth of steamwise constant (nonlinear) plane Couette flow was calculated analytically based on a streamwise constant (2D/3C) model. It is worth noting that we calculated the scalings by interpolation of the upper-bounds computed from the convex optimization problems in terms of the Reynolds numbers. In the cases wherein a polynomial fit could not be found, we measured the slope of the line in the amplification versus Reynolds number plot, Figure 3.

We next compute the worst-case amplification using the matrix inequality (4.3), which for this particular flow is given by

N =              -qx 2 0 0 M -I 3×3 0 - qy 2 0 0 0 - qy 2 -qx 2 0 0 η 2 x 0 0 0 - qy 2 0 0 η 2 y 0 0 0 - qy 2 0 0 η 2 z              0 with M as in (5.2).
The resulting upper-bounds on the worst-case amplification for 2D/3C plane Couette flow obtained by minimizing η 2

x , η 2 y , and η 2 z , are given in Figure 4. Since the flow is stable for all Reynolds numbers, the worst-case amplifications are increasing monotonically with Reynolds number. The results in Figure 4 imply

η 2 x ∼ a 0 Re 2 + a 1 Re 3 , η 2 y ∼ b 0 Re 2 + b 1 Re 4 and η 2 z ∼ c 0 Re 2 + c 1 Re 4 with a 0 , a 1 , b 0 , b 1 , c 0 , c 1 > 0.
This implies that worst-case amplification in all three components of disturbances grow with a Re 2 ratio for very low Reynolds numbers. For Reynolds numbers approximately greater than 1, the streamwise disturbances are amplified at a rate proportional to Re 3 ; whereas, the wall-normal and spanwise disturbance components disturbance amplification scales as Re 4 .

The upper-bounds depicted in Figure 4 can be compared with Corollary F.2 (see Ap-pendix F), wherein it was demonstrated that η 2 x = f 0 Re 2 , η 2 y = g 0 Re 2 + g 1 Re 4 and η 2 z = h 0 Re 2 + h 1 Re 4 for the linearized plane Couette flow with constants f 0 , g 0 , g 1 , h 0 , h 1 > 0 by solving Lyapunov equations.

Finally, we compute the stability to persistent forcings (input-to-state stability) property using inequality (4.4) from Corollary 4.1. Applying this to the 2D/3C Couette flow of interest leads to

Z =              -qx 2 0 0 M -W 0 - qy 2 0 0 0 - qy 2 -qx 2 0 0 σ x 0 0 0 - qy 2 0 0 σ y 0 0 0 - qy 2 0 0 σ z              0 with M given in (5.2) and W = qxψ 0 0 0 qyψ 0 0 0 qyψ
. We fix ψ = 10 -4 and the channel width to L = 2π. In this case, we obtain Re ISS = 316. The quantity Re ISS = 316 is interestingly close to the empirical transition Reynolds number Re ≈ 350 obtained experimentally by [START_REF] Tillmark | Experiments on transition in plane Couette flow[END_REF]. In this sense, it turns out that Re ISS is surprisingly close to the Reynolds number above which transition occurs. Understanding this connection is an interesting direction for future work.

In order to understand the above result on stability to persistent disturbances, we carried out numerical experiments to investigate the flow structures that receive maximum amplification from persistent disturbances. The experiments were undertaken for the linearized Navier-Stokes equation through the Orr-Somerfield equations. Appendix E discusses the details of these numerical experiments. Notice that these results are based on solving linear matrix inequalities that ensure stability to persistent forcings for the ODE obtained through spatial discretization of the Orr-Somerfield equations. Our simulations were carried out using a 50 × 50 grid on the wave number space k xk z (k x , k z ∈ [0, 150]) and the linear matrix inequalities were computed at each point in the grid. Then, the wave numbers corresponding to the maximum amplification are selected (especially, we are interested to find k x corresponding to maximum amplification, as this is the streamwise direction) and the corresponding flow structure is simulated. As expected the maximum amplification corresponds to the streamwise constant case k x = 0. Figure 5 illustrates the flow structures that receive maximum amplification at Re = 316.

It is also worth mentioning that certificates for stability to persistent disturbances of the linearized Navier-Stokes equation, as discussed in Appendix E, could be constructed for all Reynolds numbers, which is in contrast to the nonlinear case studied in this example. This illustrates that stability to persistent disturbances is a fundamentally nonlinear phenomenon.

Plane Poiseuille Flow

Similar to the plane Couette flow, we consider the flow of viscous fluid between two parallel plates, where the gap between the plates is much smaller than the length of the plates. Here, the plates are stationary and the flow is induced by a pressure gradient in the flow direction. The flow geometry is depicted in Figure 6.

The domain Ω is defined as Ω = {(y, z) | -1 < y < 1, 0 < z < L}. The base flow is given by U = U m (y) -→ e x = (1y 2 ) -→ e x and P = 1 -4x Re . We consider no-slip boundary u(y, z) L conditions u| 1 y=-1 = 0 and u(t, y, z) = u(t, y, z + L). The Poincaré constant is then given by C

-3 -2 -1 0 1 2 3 z -1 -0.5 0 0.5 1 y v(y, z) -3 -2 -1 0 1 2 3 z -1 -0.5 0 0.5 1 y w(y, z) -3 -2 -1 0 1 2 3 z -1 -0.
Ω y z x y = 1 y = -1 U = (1 -y 2 ) - → e x
= π 2 √ L 2 +2 2 .
We can again study the the input-output properties of the flow using the storage (5.

3)

The results are illustrated in Figure 7, which shows that the maximum energy amplification is described by γ 2 ∼ b 0 Re + b 1 Re 2 , with b 0 , b 1 > 0. This result is consistent with transient growth calculations of [START_REF] Reddy | Energy growth in viscous channel flows[END_REF]), in which the authors showed that the transient growth of the linearized plane Poiseuille flow model behaves like O(Re 2 ) for large Reynolds numbers.

For worst-case amplification analysis, we use inequality (4.3) which for this flow is given by the following matrix inequality

N =              -qx 2 0 0 M (y) -I 3×3 0 - qy 2 0 0 0 - qy 2 -qx 2 0 0 η 2 x 0 0 0 - qy 2 0 0 η 2 y 0 0 0 - qy 2 0 0 η 2 z              0, y ∈ (-1, 1),
with M as in (5.3). By minimizing η 2 x , η 2 y , and η 2 z , we compute upper-bounds on the worst-case amplification for the plane Poiseuille flow that are depicted in Figure 8. From Figure 8, it can be inferred that η 2

x ∼ a 0 Re 2 + a 1 Re 3 , η 2 y ∼ b 0 Re 2.2 + b 1 Re 4 and η 2 z ∼ c 0 Re 2 + c 1 Re 4 with a 0 , a 1 , b 0 , b 1 , c 0 , c 1 > 0. From this result, we can infer that the worstcase amplification in all three components of disturbances grow with a Re 2 ratio for low Reynolds numbers. For Reynolds numbers approximately greater than ≈ 5, the streamwise disturbances are amplified proportional to Re 3 ; whereas, the wall-normal and spanwise disturbance components show Re 4 amplification.

We again evaluate stability to persistent disturbances using inequality (4.4) from Corollary 4.1 flow, i.e.,

Z =              -qx 2 0 0 M (y) -W 0 - qy 2 0 0 0 - qy 2 -qx 2 0 0 σ x (y) 0 0 0 - qy 2 0 0 σ y (y) 0 0 0 - qy 2 0 0 σ z (y)              0, y ∈ (-1, 1),
with M given in (5.3) and W = qxψ 0 0 0 qyψ 0 0 0 qyψ . We fix ψ = 10 -4 and L = 2π. In this case, we obtain Re ISS = 1855. The quantity Re ISS = 1855 can be compared with the empirical Reynolds number at the onset of turbulence Re ≈ 2000 as discussed by [START_REF] Grossmann | The onset of shear flow turbulence[END_REF]. Once again, we observe that Re ISS provides a bound that is consistent with the Reynolds number at which transition to turbulence occurs in channel flows; an observation that leads to some interesting future research questions that are beyond the scope of the current work. Analogous to the plane Couette flow, we undertook numerical experiments to find the flow structures subject to maximum amplification from persistent forcings. Again, we found that the maximum amplification corresponds to the streamwise constant case k x = 0. Figure 9 illustrates the flow structures that receive maximum amplification from persistent forcings at Re = 1855.

Hagen-Poiseuille Flow

In this section we examine the input-output properties of Hagen-Poiseuille flow using an extension of the analysis framework that is detailed in Appendix C. In particular, we consider the flow of viscous fluid driven by a pressure gradient, illustrated in Figure 10. The domain Ω is defined as

u(y, z) -3 -2 -1 0 1 2 3 z -1 -0.5 0 0.5 1 y v(y, z) -3 -2 -1 0 1 2 3 z -1 -0.5 0 0.5 1 y w(y, z) -3 -2 -1 0 1 2 3 z -1 -0.5 0
Ω = {(r, θ) | 0 < r < 1, 0 < θ < 2π}.
The base flow is given by U = U m (r) -→ e z = (1r 2 ) -→ e z and P = 1 -4z Re . We again assume that flow is invariant in the streamwise direction z, based on observations, e.g. [START_REF] Schmid | Optimal energy density growth in Hagen-Poiseuille flow[END_REF] that axially constant perturbations are subject to maximum background energy amplification in pipe flow. Then, the perturbation dynamics is given by (C 2) in Appendix C with U m (r) = 1r 2 . Moreover, we assume no-slip boundary conditions u| r=1 = 0.

We consider the storage functional given in (C 3). Then, substituting U m and F , we have

M c (r) =   qzC Re -rq z 0 -rq z qrC Re 0 0 0 q θ C Re   .
(5.4)

Computation of the bounds on maximum energy growth in Hagen-Poiseuille flow as a function of disturbance direction and velocity component is an extension of previous analysis such as [START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF]) that considered channel flows. We obtain these bounds by solving the optimization problem (4.5) with M = M c (r) as (5.4). The results are illustrated in Figure 11. The results imply that the maximum energy growth is described by γ 2 ∼ b 0 Re + b 1 Re 2 , with b 0 , b 1 > 0. This bound is in agreement with the calculations and numerical experiments of [START_REF] Schmid | Optimal energy density growth in Hagen-Poiseuille flow[END_REF], where the authors investigated transient growth based on the linearized Navier-Stokes equations for pipe flow.

Considering M c (r) as in (5.4), inequality (C 12) becomes

N c (r) =              -qz 2 0 0 M c (r) -I 3×3 0 -qr 2 0 0 0 -q θ 2 -qz 2 0 0 η 2 z 0 0 0 -qr 2 0 0 η 2 r 0 0 0 -q θ 2 0 0 η 2 θ              0, r ∈ (0, 1).
(5.5)

Minimizing η 2 z , η 2 r and η 2 θ subject to the above inequality provides upper-bounds on the worst-case disturbance amplification, which are plotted in Figure 12 as a function of Reynolds number. Figure 12 indicates that the perturbations are amplified as η 2 z ∼ a 0 Re 2 +a 1 Re 3 , η 2 θ ∼ b 0 Re 2 +b 1 Re 4 , and η 2 r ∼ c 0 Re 2 +c 1 Re 4 with a 0 , a 1 , b 0 , b 1 , c 0 , c 1 > 0. Thus, similar to channel flows, for low Reynolds numbers, the worst-case amplification from all three disturbance components scale as Re 2 . For Reynolds numbers greater than ≈ 8, the amplification from axial (flow direction) disturbances grow proportional to Re 3 ; whereas, the worst-case amplification from azimuthal and radial disturbances increase as Re 4 .

We compute the bounds on stability to persistent forcings using the following polynomial matrix inequality (5.6) where W c = ψqz 0 0 0 ψqr 0 0 0 ψq θ , which is analogous to inequality (C 13) in Corollary C.2. The maximum Reynolds number for that the system remained stable to persistent forcings (i.e. remained input-to-state stable) was Re ISS = 1614. Obtaining this bound required the use of degree 10 polynomials in σ z (r), σ θ (r) and σ r (r). This result is again consistent with the lower bound on the Reynolds number for which transition to turbulence was observed empirically by [START_REF] Peixinho | Decay of turbulence in pipe flow[END_REF], i.e., Re ≈ 1800 and provides further support that this is a connection worth investigating.

Z c (r) =              -qz 2 0 0 M c (r) -W c 0 -qr 2 0 0 0 -q θ 2 -qz 2 0 0 σ z (r) 0 0 0 -qr 2 0 0 σ r (r) 0 0 0 -q θ 2 0 0 σ θ (r)              0, r ∈ (0, 1),

Flow Energy Growth

Worst-Case Amplification ISS vs. Transition Table 1: Summary of the numerical results using the proposed framework (boldfaced), and results obtained in the literature.

Plane Couette O(Re 3 ), O(Re 3 ) †   O(Re 3 ) O(Re 4 ) O(Re 4 )   ,   O(Re 3 ) O(Re 4 ) O(Re 4 )   ‡ 316, 340 ¶ Plane Poiseuille O(Re 2 ), O(Re 2 )   O(Re 3 ) O(Re 4 ) O(Re 4 )   ,   O(Re 3 ) O(Re 4 ) O(Re 4 )   † † 1855, 2000 ‡ ‡ Hagen-Poiseuille O(Re 2 ), O(Re 2 ) ¶ ¶   O(Re 3 ) O(Re 4 ) O(Re 4 )   , - 1614

Discussion

We studied stability and input-output properties of fluid flows with spatially invariant perturbations in one of the directions using dissipation inequalities. Our framework generalizes certain types of input-output analysis techniques to the nonlinear Navier-Stokes equations, thereby matching more closely with experimental results. The proposed input-output analysis method introduces a unified framework for addressing a broad range of questions related to transition (transient growth and input-output analysis) that can be adapted to a large class of flows, specially, for those in which the base flow is given by a polynomial of spatial coordinates and the flow geometry is described by a semi-algebraic set. We showed how the input-output framework can be computationally implemented based on convex optimization. For illustration purposes, we applied the proposed method to study several examples of flows between parallel plates and pipe flow.

Table 1 lists the numerical results based on the proposed framework for plane Couette flow, plane Poiseuille flow, and the Hagen-Poiseuille flow. For energy growth and worstcase amplification, the table outlines the amplification scalings at high Reynolds numbers. Energy growth results for all three flows are consistent with the theoretical and experimental amplification scalings in the literature. Our worst-case amplification scalings for plane Couette flow and plane Poiseuille flow were similar to the scalings calculated by [START_REF] Bamieh | Energy amplification in channel flows with stochastic excitation[END_REF]; [START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF]. In addition to comparing the scalings we obtained using our framework for channel flows, we carried out numerical experiments to study the worst-case amplification scalings in Hagen-Poiseuille flow. Our results indicate that, similar to channel flows, perturbations in the direction of the base flow are least amplified in Hagen-Poiseuille flow. For transition analysis, we compare the maximum Reynolds numbers for which stability to persistent disturbances could be certified to the Reynolds numbers for which transition to turbulence was observed experimentally. We made an interesting observation from the results that Re ISS is very similar to the Reynolds number for which transition to turbulence occurs in practice.

Future research will focus on applying the framework proposed here to turbulent flows [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF]). In particular, we study time-averaged mechanical energy dissipation. For a channel flow of channel length h, the mechanical energy dissipation per unit mass is given by

ε := ν 3 h 4 ∇u 2 L 2 Ω ,
where ν is the kinematic viscosity. [START_REF] Doering | Variational bounds on energy dissipation in incompressible flows: Shear flow[END_REF] proposed a variational method for bounding this quantity based on the background flow decomposition. The method has been significantly successful in finding the time-averaged mechanical energy dissipation scaling with respect to the root-mean-square velocity U and the longest length scale, i.e., it was shown that

ε c 1 ν U 2 2 + c 2 U 3 .
and bounds on c 1 and c 2 were obtained for different flows [START_REF] Doering | Energy dissipation in body-forced turbulence[END_REF] (2004)). In order to find bounds on the time-averaged mechanical energy dissipation, we can consider the following dissipation inequality

dV (u) dt ν 3 h 4 ∇u 2 L 2 Ω -C, (6.1)
where C > 0 is a constant. Minimizing C while searching over the storage functional V (u) gives upper bounds on the time-averaged mechanical energy dissipation. It is worth pointing out that the solution to the dissipation inequality (6.1) was studied using the background method (Chernyshenko (2017)) for incompressible viscous flows. For the Kuramoto-Sivashinsky PDE, semi-definite programs have been used to study timeaveraged mechanical energy dissipation [START_REF] Goluskin | Bounds on mean energy in the kuramoto-sivashinsky equation computed using semidefinite programming[END_REF]) and a corresponding toolbox was developed (see also the application of the toolbox in analyzing bounds for Rayleigh-Bérnard convection in Fantuzzi ( 2018)).

Another interesting problem for future research is identifying the regions of attraction for different flow configurations. For example, in the case of Taylor-Couette flow, after decomposing the Navier-Stokes equation about different flow states, one can search for estimates of the region of attraction inside which each flow state is stable.

In addition, input-output amplification mechanisms of turbulent flows is also an intriguing prospective research direction. In this regard, del Alamo & Jiménez (2006); [START_REF] Pujals | A note on optimal transient growth in turbulent channel flows[END_REF], consider a non-polynomial model for turbulent mean velocity profiles and turbulent eddy viscosities. Polynomial approximations (of high degrees) of such nonlinear models fit the formulation given in this paper.

Lastly, more general storage functional structures can be considered. More specifically, given the nonlinear dynamics of the Navier-Stokes equations, one can consider the following class of storage functionals

V (u) = Ω u u 2 Q u u 2 dΩ.
However, a convex formulation using the above structure is not clear at the moment. Other structures for the storage functionals and convex methods for finding them were proposed by [START_REF] Gahlawat | A convex sum-of-squares approach to analysis, state feedback and output feedback control of parabolic PDEs[END_REF]; [START_REF] Shivakumar | Computing input-output properties of coupled PDE systems[END_REF] for linear PDEs. Extention of the latter results to nonlinear PDEs can pave the way for input-output analysis of the Navier-Stokes equations using more general storage functional structures. Thus, inequality (3.7) is also satisfied. Applying Item I from Theorem 3.4, we infer that the system has bounded energy growth. II) Re-arranging terms in (B 2) yields

- i∈I q i Ω C(Ω) Re u 2 i + U j u i ∂ xj u i + u j u i ∂ xj U i -u i d i dΩ - i∈I Ω u 2 i dΩ + i∈I Ω η 2 i d 2 i dΩ (B 4)
Then, from (3.13) in Proposition 3.7,

dV (u) dt - i∈I q i Ω C(Ω) Re u 2 i + U j u i ∂ xj u i + u j u i ∂ xj U i dΩ, we deduce that dV (u) dt - i∈I Ω u 2 i dΩ + i∈I Ω η 2 i d 2 i dΩ.
From Theorem 3.5, we infer that, under zero initial conditions, the perturbation velocities satisfy (3.3). III) Adopting (3.12) as a storage functional, we have

λ m (Q) Ω u u dΩ V (u(t, x)) λ M (Q) Ω u u dΩ,
where λ m (Q) denotes the minimum eigenvalue of Q. Hence, (3.9) is satisfied with β 1 (•) = min i∈I q i (•) 2 and β 2 (•) = max i∈I q i (•) 2 . Re-arranging the terms in (B 3), we obtain

- i∈I Ω q i C(Ω) Re u 2 i + q i U j u i ∂ xj u i + q i u j u i ∂ xj U i -q i u i d i dΩ - i∈I ψ i Ω q i u 2 i dΩ + Ω σ(|d|) dΩ (B 5)
From (3.13) in Proposition 3.7, it follows that

dV (u) dt -ψV (u) + Ω σ(|d|) dΩ. (B 6)
Then, from Theorem 3.6, we infer that the perturbation velocities are stable to persistent focings (3.4).

B.1. Proof of Corollary 4.1

The proof is straightforward and follows from computing conditions (B 2)-(B 3) considering perturbations that are constant in x 1 , the base flow U = U m -→ e 1 , and σ(|d|) = i∈I σ i (x)d 2 i . Since the flow perturbations are constant in x 1 and the base flow is given by U = U m -→ e 1 , we have U j u i ∂ xj u i = 0, i ∈ I. The right hand side of (3.13) hence changes to

A = Ω C(Ω) Re q i u 2 i + C(Ω) Re q j u 2 j + C(Ω) Re q 1 u 2 1 + u 1 ∂ xi U m u i + u 1 ∂ xj U m u j dΩ (B 7)
for i, j ∈ I 0 , i = j. Since the integrand in the above expression is quadratic in u i , i ∈ I, it can be rewritten as

A = Ω u1 uj ui M (x) u1 uj ui dΩ. (B 8)
with M (x) given in (4.1). I) Given storage functional structure (3.12), we have

V (u(t, x)) λ M (Q) Ω u u dΩ,
where λ M (Q) denotes the maximum eigenvalue of Q. Since Q is diagonal, we have λ M (Q) = max i∈I q i . Therefore, (3.6) is satisfied with γ 2 = max i∈I q i . Re-arranging terms in (B 1) yields

- i∈I q i Ω C(Ω) Re u 2 i + U j u i ∂ xj u i + u j u i ∂ xj U i dΩ i∈I Ω u 2 i dΩ.
Applying Proposition 3.7 with d ≡ 0, we obtain

dV (u) dt i∈I Ω u 2 i dΩ.
Thus, inequality (3.7) is also satisfied. Applying Theorem 3.4, we infer that the system has bounded energy growth.

II) Inequality (B 2) is changed to

A + Ω (q i u i d i + q j u j d j + q 1 u 1 d 1 ) dΩ - Ω (u 2 i + u 2 j + u 2 1 ) dΩ + Ω (η 2 i d 2 i + η 2 j d 2 j + η 2 1 d 2 1 ) dΩ 0, ( B 9)
for i, j ∈ I 0 , i = j. Moving all the terms under one integral gives

A + Ω q i u i d i + q j u j d j + q 1 u 1 d 1 -u 2 i -u 2 j -u 2 1 + η 2 i d 2 i + η 2 j d 2 j + η 2 1 d 2 1 dΩ 0. (B 10)
Since the integrand in the above expression is quadratic in u i 's and d i 's, it can be rewritten as

Ω   u1 uj ui d1 dj di   N (x)   u1 uj ui d1 dj di   dΩ 0, ( B 11)
where N is defined in (4.3). Consequently, if (4.3) is satisfied for all x ∈ Ω, (B 11) holds and from Item II in Corollary B.1 we infer that, subject to zero initial conditions, the worst-case amplification from disturbances to perturbation velocities is bounded by η i , i ∈ I as in (3.3).

III) The proof follows the same lines as the proof of Item II above.

defined as ∇ c (•) = ∂ r (•) - → e r + 1 r ∂ θ (•) - → e θ + ∂ z (•) - → e z and ∇ 2 c (•) = 1 r ∂ r (r∂ r (•)) + 1 r 2 ∂ 2 θ (•) + ∂ 2 z (•).
The Navier-Stokes equations in cylindrical coordinates are then given by

∂ t ūr = 1 Re ∇ 2 c ūr - ūr r 2 - 2 r 2 ∂ θ ūθ -ū • ∇ c ūr + ū2 θ r -∂ r p + d r ∂ t ūθ = 1 Re ∇ 2 c ūθ - ūθ r 2 + 2 r 2 ∂ θ ūr -ū • ∇ c ūθ - ūθ ūr r - 1 r ∂ θ p + d θ ∂ t ūz = 1 Re ∇ 2 c ūz -ū • ∇ c ūz -∂ z p + d z 0 = 1 r ∂ r (rū r ) + 1 r ∂ θ ūθ + ∂ z ūz , ( C 1)
where ū = (ū r , ūθ , ūz ) .

We consider the flow perturbations that are invariant in the axial direction (zdirection). The base flow is given by U = U m (r, θ) -→ e z and P . For such flows, substituting ū = u + U and p = P + p in (C 1), the perturbation dynamics is obtained as

∂ t u r = 1 Re ∇ 2 c u r -u r ∂ r u r - u θ ∂ θ u r r + u 2 θ r - u r r 2 Re - 2∂ θ u θ r 2 Re -∂ r p + d r , ∂ t u θ = 1 Re ∇ 2 c u θ -u r ∂ r u θ - u θ ∂ θ u θ r - u r u θ r - u θ r 2 Re - 2∂ θ u θ r 2 Re - 1 r ∂ θ p + d θ , ∂ t u z = 1 Re ∇ 2 c u z -u r ∂ r u z -u r ∂ r U m - u θ ∂ θ U m r - u θ ∂ θ u z r + d z , 0 = ∂ r (ru r ) + ∂ θ u θ , ( C 2)
wherein u = (u r , u θ , u z ) .

Proposition C.1. Consider the perturbation dynamics in cylindrical coordinates (C 2) with periodic or no-slip boundary conditions u | ∂Ω = 0. The time derivative of storage functional

V (u) = 1 2 Ω ur u θ uz qr 0 0 0 q θ 0 0 0 qz ur u θ uz rdrdθ, ( C 3)
with q r = q θ , satisfies

dV (u) dt - Ω q r C Re u 2 r + q z ∂ r U m u r u z + q z C Re u 2 z + q z r ∂ θ U m u θ u z + q θ C Re u 2 θ -q r u r d r -q θ u θ d θ -q z u z d z rdrdθ, (C 4)
where C > 0.

Proof: The time derivative of the storage functional (C 3) is given by

dV (u) dt = Ω -rq r u 2 r ∂ r u r -q r u r u θ ∂ θ u r + q r u r u 2 θ -rq r ∂ r pu r + q r Re ru r ∇ 2 c u r - q r u 2 r Rer - 2q r u r ∂ θ u θ rRe + q r ru r d r dθdr + Ω -rq θ u r u θ ∂ r u θ -q θ u r u θ ∂ θ u r -q θ u r u 2 θ -q θ ∂ θ pu θ + q θ Re ru θ ∇ 2 c u θ - q θ u 2 θ rRe + 2q θ ∂ θ u r u θ rRe + q θ ru θ d θ drdθ + Ω -rq z u r u z ∂ r u z -rq z u r u z ∂ r U m -q z ∂ θ U m u θ u z -q z u θ u z ∂ r u z + q z Re ru z ∇ 2 c u z + q z ru z d z drdθ. (C 5)
We first focus on the terms involving pressure p in the right-hand side of (C 5). From the incompressibility condition ∂ r (ru r ) + ∂ θ u θ = 0 and the fact that q r = q θ , we obtain

Ω (-rq r ∂ r pu r -q θ ∂ θ pu θ ) drdθ = Ω (q r ∂ r (ru r )p + q θ ∂ θ up) drdθ = Ω q r p (∂ r (ru r ) + ∂ θ u) drdθ = 0. (C 6)
where, in the first equality above, we used integration by parts and the boundary conditions. At this point, we consider the other terms. Using integration by parts, boundary conditions and the incompressibility condition it can be shown that

Ω -rq r u 2 r ∂ r u r -q r u r u θ ∂ θ u r drdθ = Ω q r u 2 r 2 ∂ r (ru r ) + q r u 2 Ω - 2q r u r ∂ θ u θ rRe - 2q θ ∂ θ u r u θ rRe drdθ = Ω 2q r rRe (-u r ∂ θ u θ + u r ∂ θ u θ ) drdθ = 0.
Then, the time derivative expression (C 5) simplifies to

dV (u) dt = Ω q r Re ru r ∇ 2 c u r - q r u 2 r Rer + q r ru r d r dθdr + Ω q θ Re ru θ ∇ 2 c u θ - q θ u 2 θ rRe + q θ ru θ d θ drdθ + Ω -rq z u r u z ∂ r U m -q z ∂ θ U m u θ u z + q z Re ru z ∇ 2 c u z + q z ru z d z drdθ. (C 7) Factoring out r yields dV (u) dt = Ω q r Re u r ∇ 2 c u r - q r u 2 r r 2 Re + q r u r d r rdθdr + Ω q θ Re u θ ∇ 2 c u θ - q θ u 2 θ r 2 Re + q θ u θ d θ rdrdθ + Ω -q z u r u z ∂ r U m - q z r ∂ θ U m u θ u z + q z Re u z ∇ 2 c u z + q z u z d z rdrdθ. (C 8)
Since the terms qru 2 r r 2 Re and q θ u 2 θ r 2 Re are non-negative, it follows that

dV (u) dt Ω q r Re u r ∇ 2 c u r + q θ Re u θ ∇ 2 c u θ + q z Re u z ∇ 2 c u z -q z u r u z ∂ r U m - q z r ∂ θ U m u θ u z -q r u r d r -q θ u θ d θ -q z u z d z rdrdθ = - Ω q r Re |∇ c u r | 2 + q θ Re |∇ c u θ | 2 + q z Re |∇ c u z | 2 + q z u r u z ∂ r U m + q z r ∂ θ U m u θ u z -q r u r d r -q θ u θ d θ -q z u z d z rdrdθ, (C 9)
where in the last equality above integration by parts and the boundary conditions were used. Applying the Poincaré inequality [START_REF] Payne | An optimal Poincare inequality for convex domains[END_REF]), we obtain (C 4).

C.1. Convex Formulation: Pipe Flows

Similar to the case of channel flows, in the following, we propose a convex formulation for pipe flows. The method relies on inequality (C 4). Note that for cylindrical coordinates I = {r, θ, z} and I 0 = {r, θ}.

Corollary C.2. Consider the perturbation dynamics given by (C 2), streamwise constant in the z-direction with base flow U = U m (r, θ) -→ e z . Suppose that there exist positive constants {q l } l∈I with q r = q θ , ψ, and functions {σ l } l∈I such that

M c (r, θ) = C Re qz 1 2 qz∂rUm qz 2 ∂ θ Um r 1 2 qz∂rUm C Re qr 0 qz 2 ∂ θ Um r 0 C Re q θ , ( C 10)
I) M c (r, θ) -I 3×3 0, (r, θ) ∈ Ω, (C 11) II) N c (r, θ) =              -qz 2 0 0 M c (r, θ) -I 3×3 0 -qr 2 0 0 0 -q θ 2 -qz 2 0 0 η 2 z 0 0 0 -qr 2 0 0 η 2 r 0 0 0 -q θ 2 0 0 η 2 θ              0, (r, θ) ∈ Ω, (C 12) III) σ l (r, θ) 0, (r, θ) ∈ Ω, l ∈ I and Z c (r, θ) =              -qz 2 0 0 M c (r, θ) -W c 0 -qr 2 0 0 0 -q θ 2 -qz 2 0 0 σ z (r, θ) 0 0 0 -qr 2 0 0 σ r (r, θ) 0 0 0 -q θ 2 0 0 σ θ (r, θ)              0, (r, θ) ∈ Ω, ( C 13)
where W c = ψqz 0 0 0 ψqr 0 0 0 ψq θ . Then, it follows that I) the flow has bounded energy growth γ 2 = max(q r , q θ , q z ) as given by (3.2), II) subject to zero initial conditions, the induced L 2 norm from inputs to perturbation velocities is bounded, III) the perturbation velocities are ISS in the sense of (3.4) with σ(|d|) = l∈I σ l (r, θ)d 2 l . Note that, depending on ∂ θ U m , M c and therefore N c and Z c can be functions of 1 r . Then, inequalities (C 10)-(C 13) become intractable. To circumvent this problem, since r is positive, we can multiply (C 10)-(C 13) by positive powers of r making the resulting inequalities solvable by convex optimization methods.

Appendix D. Polynomial Optimization and Sum-of-Squares Programming

Let R[x] denote the set of polynomials in x with real coefficients and Σ[x] ⊂ R[x] the set of such polynomials with a sum-of-squares decomposition. We employ sum-ofsquares programming in our computational formulations whenever ∇U depends on the spatial variables for channel flows (one exception in the plane Couette flow) and for all pipe flows. That is, we convert different input-output analysis problems into a sumof-squares program (SOSP), i.e., an optimization problem involving a linear objective function subject to a set of polynomial constraints as given below minimize

c∈R N w c subject to a 0,j (x) + N i=1 p i (x)a i,j (x) = 0, j = 1, 2, . . . , J, a 0,j (x) + N i=1 p i (x)a i,j (x) ∈ Σ[x], j = J + 1, J + 2, . . . , J, ( D 1)
where w ∈ R N is a vector of weighting coefficients, c ∈ R N is a vector formed of the (unknown) coefficients of

{p i } N i=1 ∈ R[x] and {p i } N i= N +1 ∈ Σ[x], a i,j (x) ∈ R[x] are given scalar constant coefficient polynomials, p i (x) ∈ Σ[x] are sum-of-squares polynomial (SOSP) variables.
The main idea behind sum-of-squares programming is that checking whether a polynomial is positive can be relaxed into checking whether it has a sum-of-squares decomposition. Formally, if there exists a sum-of-squares decomposition for p

(x) ∈ R[x], i.e., if there exist polynomials f 1 (s), . . . , f m (x) ∈ R[x] such that p(x) = m i=1 f 2 i (x),
then it follows that p(x) is non-negative. Unfortunately, the converse does not hold in general † ; that is, there exist non-negative polynomials which do not have an sumof-squares decomposition. An example of this class of non-negative polynomials is the Motzkin's polynomial [START_REF] Motzkin | The arithmetic-geometric inequality[END_REF]) given by

p(x) = 1 -3x 2 1 x 2 2 + x 2 1 x 4 2 + x 4 1 x 2 2 , ( D 2)
which is non-negative for all x ∈ R 2 but is not a SOS. This imposes some degree of conservatism when utilizing sum-of-squares based methods. Generally, determining whether a given polynomial is positive is an NP-hard problem [START_REF] Bovet | Introduction to the Theory of Complexity[END_REF]) (except for degrees less than 4); but, sum-of-squares decompositions provide a conservative, yet computationally feasible method for checking non-negativity. The next lemma gives an intriguing formulation to the sum-of-squares decomposition problem.

Lemma D.1 [START_REF] Choi | Sums of squares of real polynomials[END_REF]). A polynomial p(x) of degree 2d belongs to Σ[x] if and only if there exist a positive semi-definite matrix Q (known as the Gram matrix) and a vector of monomials Z(x) which contains all monomial of x of degree d such that p(x) = Z T (x)QZ(x).

For instance, consider the polynomial

p(x 1 , x 2 ) = 5x 4 1 + 2x 2 2 -x 2 1 x 2 2 -2x 3 1 x 2 -2x 1 x 3 2 . Factoring out the monomials in p as Z(x 1 , x 2 ) = (x 2 1 , x 2 2 , x 1 x 2 ) , we have p(x 1 , x 2 ) =   x 2 1 x 2 2 x 1 x 2    
q 11 q 12 q 13 q 21 q 22 q 23 q 31 q 32 q 33

    x 2 1 x 2 2 x 1 x 2   = Z T (x)QZ(x).
Then, by comparison, we have q 11 = 5, q 22 = 2, 2q 12 + q 33 = -1, q 13 = -1 and q 23 = -1.

If we can find a Q 0, then p(x 1 , x 2 ) is positive semi-definite.

In [START_REF] Chesi | On convexification of some minimum distance problems[END_REF]) and [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF]) it was demonstrated that the answer to the query that whether a given polynomial p(x) is sum-of-squares or not can be investigated via semi-definite programming methodologies.

Lemma D.2 [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF]). Given a finite set

{p i } m i=0 ∈ R[x], the existence of a set of scalars {a i } m i=1 ∈ R such that p 0 + m i=1 a i p i ∈ Σ[x] (D 3)
is a linear matrix inequality (LMI) ‡ feasibility problem.

Furthermore, we often need to verify whether a matrix with polynomial entries is positive (semi)definite. To this end, we use the next lemma from [START_REF] Prajna | Nonlinear control synthesis by sum of squares optimization: A Lyapunov-based approach[END_REF]).

Lemma D.3 [START_REF] Prajna | Nonlinear control synthesis by sum of squares optimization: A Lyapunov-based approach[END_REF]). Denote by ⊗ the Kronecker product. Suppose † Exceptions [START_REF] Reznick | Some concrete aspects of Hilbert's 17th problem[END_REF]): Univariate polynomials of any even degree; Quadratic polynomials in any number of variables; Quartic polynomials in two variables.

‡ An LMI is an expression of the form A0 + x1A1 + x2A2 + • • • + xmAm 0, where x ∈ R m and Ai ∈ S n , i = 1, 2, . . . , m (symmetric n × n matrices). An LMI specifies a convex constraint on the variables x.

bation velocity, η := ∂ z u -∂ x w, the cross-stream component of perturbation vorticity (z denotes the spanwise direction). Velocity has been non-dimensionalized by U 0 , the maximum velocity in the channel; length has been non-dimensionalized by L, the width of the channel. The Reynolds number is defined as Re := U0L ν , where ν is the kinematic viscosity. Considering no-slip boundary conditions at y = ±1, we have v = ∂ y v = η = 0 at y = ±1. Recall that for the plane Couette flow U = y, and for the plane Poiseuille flow U = 1y 2 .

Consider a single Fourier component v = ve ikxx+ikzz , (E 2) η = ηe ikxx+ikzz . (E 3)

Physical variables being identified with the real part of these complex form. The field equations can be written in the compact form

∂ t v η = L 0 C S v η , ( E 4)
in which the Orr-Sommerfield operator L , the Square operator S , and the coupling operator C are defined as in which the linear dynamical operator, A , is the discretized form of L 0 C S . This means that the infinite dimensional dynamical system (E 1), is approximated as a finite dimensional dynamical systems.

L = ∆ -1 -ik x U ∆ + ik x ∂ 2 y U + ∆∆ Re , ( E 5
The discretized operator A was calculated using the codes available in [START_REF] Schmid | Stability and transition in shear flows[END_REF], Appendix A) using Chebyshev discretization. For both flows, we considered N = 50. Then, the state-space form (E 10) is a linear system that has to be studied. In the following, we obtain linear matrix inequality conditions to check inputto-state stability (ISS) of a linear system. Now, consider the following linear dynamical system Applying the bounds in (E 13), we obtain

c 1 ( ζ 2 ) e -ψt c 2 ( ζ 0 2 ) + 1 ψ sup τ ∈[0,t)
S (d(τ )) .

Since c 1 ∈ K, its inverse exists and belongs to K. Thus, Noting that S is positive semidefinite, we have

ζ 2 c -1 1 2e -ψt c 2 ( ζ 0 2 ) + c -1 1 2 ψ S d L ∞ [0,t)
.

The following corollary gives sufficient conditions based on linear matrix inequalities to check the conditions of Theorem E.1.

Corollary E.2. Consider system (E 11). If there exist symmetric matrices P and S, and a positive scalar ψ such that P 0, S 0 (E 15) In order to the find the maximum ISS amplification, we solve the following optimization problem For the wave numbers that correspond to the maximum ISS amplification, we obtain the direction in which maximum amplification is attained. To this end, we carry out a singular-value decomposition of P (since P is symmetric the singular values and eigenvalues coincide) and we obtain the eigenvector in A that corresponds to the maximum singular value.

Appendix F. Induced L 2 [0,∞),Ω -norms for the Linearized 2D/3C Model

In [START_REF] Jovanović | Modeling, analysis, and control of spatially distributed systems[END_REF]), the authors calculated componentwise H ∞ -norms for the linearized 2D/3C model by finding the maximum singular values. This result is described as follows.

Theorem F.1 (Thoerem 11, p. 93 in Jovanović (2004)). For any streamwise constant channel flows with nominal velocity U (y), the H ∞ norms of operators H rs (ω, k z , Re) that maps d s into u r , {r = x, y, z; s = x, y, z}, are given by 

 H xx ∞ (k z ) H xy ∞ (k z ) H xz ∞ (k z ) H yx ∞ (k z ) H yy ∞ (k z ) H yz ∞ (k z ) H zx ∞ (k z ) H zy ∞ (k z ) H zz ∞ (k z )   =   h xx (k z )

Figure 2 :

 2 Figure 2: Schematic of the plane Couette flow geometry.

Figure 3 :Figure 4 :

 34 Figure 3: Upper bounds on the maximum energy growth for plane Couette flow in terms of Reynolds numbers.

Figure 5 :

 5 Figure 5: The perturbation flow structures with maximum amplification from persistent forcings at Re = 316 for plane Couette flow.

Figure 6 :

 6 Figure 6: Schematic of the plane Poiseuille flow geometry.

Figure 7 :

 7 Figure 7: Upper bounds on the maximum energy growth for plane Poiseuille flow in terms of Reynolds numbers.

Figure 8 :

 8 Figure 8: Upper bounds on the worst-case amplification of plane Poiseuille flow for different Reynolds numbers.

Figure 9 :Figure 10 :

 910 Figure 9: The perturbation flow structures with maximum amplification to persistent disturbances Re = 1855 for plane Poiseuille flow.

Figure 11 :

 11 Figure 11: Upper bounds on the maximum energy growth for Hagen-Poiseuille flow in terms of Reynolds numbers.

Re

Figure 12 :

 12 Figure 12: Upper bounds on the worst-case amplification of the Hagen-Poiseuille flow in terms of different Reynolds numbers.

  , 1800 † Bobba et al. (2002) ‡ Jovanović & Bamieh (2005) ¶ Tillmark & Alfredsson (1992) Reddy & Henningson (1993) † † Jovanović & Bamieh (2005) ‡ ‡ Grossmann (2000) ¶ ¶ Schmid & Henningson (1994)[START_REF] Peixinho | Decay of turbulence in pipe flow[END_REF]

  ; Childress et al. (2001); Alexakis & Doering (2006); Rollin et al. (2011); Tang et al.

  k y ηk x ∂ y v) , k x η + k y ∂ y v) . (E 9)For numerical simulations of the Orr-Somerfield equation (E 1), we consider its discrete equivalent for an N -level discretization (over space)ζ = v1 • • • vN η1 • • • ηN ,and the initial value problem (E 1) can be rewritten asζ = A ζ, (E 10)

  ζ = A ζ + Bd, t > 0, (E 11)whereζ(0) = ζ 0 , ζ ∈ R 2N , d ∈ R 2N and B = I 2N ×2N. This is the perturbed version of the discrete system (E 10). We are interested in studying the ISS of (E 10). That is, given d ∈ L ∞ , we have the following inequality for allζ 0 ∈ R 2N ζ(t) 2 β (t, ||ζ 0 || 2 ) + σ d L ∞ β ∈ KL, σ ∈ K and ζ(t) 2 is the Euclidean 2-norm, i.e., ζ(t) 2 = √ ζ ζ.Theorem E.1. Consider system (E 11). If there exists an ISS-storage function V (ζ) and a positive semidefinite function S, c 1 , c 2 ∈ K, and a positive scalar ψ satisfyingc 1 ( ζ 2 ) V (ζ) c 2 ( ζ 2 ) , (E 13)and∂ t V (ζ) -ψV (ζ) + S(d), (E14) then solutions of (E 11) satisfy estimate (E 12) with β(•) = c -1 1 2e -ψt c 2 (•) and σ(•) = both sides of (E 14) by e ψt , gives e ψt ∂ t V (ζ) -e ψt ψV (ζ) + e ψt S(d) which implies d dt e ψt V (ζ) e ψt S(d). Integrating both sides of the latter inequality from 0 to t yields e ψt V (ζ(t)) -V (ζ 0 ) the last inequality, we applied the Hölder inequality. Then, e ψt V (ζ(t)) -V (ζ 0 ) e ψt -Dividing both sides of the last inequality above by the non-zero term e ψt and re-arranging the terms gives V (ζ(t)) e -ψt V (ζ 0 ) + 1 ψ sup τ ∈[0,t) S (d(τ )) .

  c 2 ( ζ 0 2 ) + c -1 1 2 ψ sup τ ∈[0,t)S (d(τ )) .

  to (E 11) satisfy (E 12) with for β(•) = 2λ M (P ) λm(P ) e -ψt (•) is a result of applying Theorem E.1 by considering V (ζ) = ζ P ζ and S(d) = d Sd.

  minimize P,S (λ 1λ 2 ) subject to S λ 1 I, P λ 2 I, (E 15), and (E 16).(E 17) Then, the system satisfies inequality (E 12) with β(•) = 2λ M (P )

  Re h xy (k z )Re 2 h xz (k z )Re 2 0 h yy (k z )Re h yz (k z )Re 0 h

zy (k z )Re h zz (k z )Re   , (F 1)

where k z represent the wavenumber in x z (spanwise direction).

† Email address for correspondence: mrahmadi@caltech.edu

† In the fluids literature, the ensemble average energy density or the H 2 -norm has been used to study amplifications from Gaussian stochastic forcings to the linearized flow dynamics[START_REF] Farrell | Stochastic forcing of the linearized Navier-Stokes equations[END_REF];[START_REF] Jovanović | Componentwise energy amplification in channel flows[END_REF]). The H 2 -norm is equivalent to the (root mean square) RMS-value of the linearized flow response to white noise forcings. However, extension of H 2 analysis to the nonlinear Navier-Stokes equations is an open problem.

Appendix A. Proof of Proposition 3.7

The time derivative of storage functional (3.12) along the solutions of (3.11) can be computed as

We first consider the term Ω q i u j u i ∂ xj u i dΩ. Using the boundary conditions, integration by parts and the incompressibility condition ∂ xj u j = 0, we obtain

Secondly, we consider the pressure terms Ω q i u i ∂ xi p dΩ. Since the perturbations are assumed invariant in x 1 implying ∂ x1 p = 0, we have

where, in the first equality above, we use integration by parts and, in the second inequality, we use the boundary conditions. Then, since q 2 = q 3 , using the incompressibility condition

Finally, from (A 2) and (A 3), we conclude that the time derivative of the storage functional (A 1) is

Integrating by parts the u i ∇ 2 u i term in (A 4) and using the boundary conditions, we get

Applying the Poincaré inequality [START_REF] Payne | An optimal Poincare inequality for convex domains[END_REF]) to (A 5), we obtain (3.13).

Appendix B. Derivation of the Convex Programs for Channel Flows

The next corollary proposes integral inequalities under which properties such as energy growth bounds, worst-case amplification, and stability to persistent forcings can be inferred for the flow described by (3.11).

Corollary B.1. Consider the perturbation dynamics described by (3.11) subject to periodic or no-slip boundary conditions u| ∂Ω = 0. Assume the velocity perturbations are constant with respect to x 1 . Let I 0 = {2, 3}. If there exist positive constants q i , i ∈ I, with q i = q j , i, j ∈ I 0 , positive scalars ψ, {η i } i∈I , and σ ∈ K such that I) when d ≡ 0,

)

Then, I) system (3.11) has bounded energy growth as described by (3.2) with γ 2 = max i∈I q i ; II) under zero initial perturbations u(0, x) ≡ 0, the worst-case amplification from disturbances to perturbation velocities is bounded by η i , i ∈ I as in (3.3); III) the perturbation velocities described by (3.11) are stable to persistent forcings in the sense of (3.4).

Proof: Each item is proven as follows. I) Given storage functional structure (3.12), we have

Therefore, (3.6) is satisfied with γ 2 = max i∈I q i . Re-arranging terms in (B 1) yields

Applying Proposition 3.7 with d ≡ 0, we obtain

Appendix C. Pipe Flows: Cylindrical Coordinates

In this appendix, we extend the proposed method to flows in cylindrical coordinates (r, θ, z). In cylindrical coordinates, the gradient and Laplacian operators are, respectively,

is symmetric and of degree 2d for all x ∈ R n . In addition, let Z(x) ∈ R n×1 [x] be a column vector of monomials of degree no greater than d and consider the following conditions

(C) There exists a positive semi-definite matrix Q such that

for any v ∈ R n . Then (A) ⇐ (B) and (B) ⇔ (C).

We can also check the positivity of a matrix with polynomial entries F (x) ∈ R n×n [x] inside a set Ω ⊂ R n . It turns out that if the set is semi-algebraic, Putinar's Positivstellensatz (Lasserre 2009, Theorem 2.14) can be used.

If the coefficients of F (x) depend affinely in unknown parameters and the degree of N (x) is fixed, checking whether (D 4) holds can be cast as a feasibility test of a convex set of constraints, a semi-definite program (SDP), whose dimension depends on the degree of the polynomial entries of F (x) and N (x).

Algorithms for solving sum-of-squares programs are automated in MATLAB toolboxes such as SOSTOOLS [START_REF] Papachristodoulou | SOSTOOLS: Sum of squares optimization toolbox for MATLAB[END_REF]) and YALMIP [START_REF] Löfberg | YALMIP : A toolbox for modeling and optimization in MATLAB[END_REF]), in which the sum-of-squares problem is parsed into an SDP formulation and the SDPs are solved by LMI solvers such as SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF].

For example, one can solve a SOS optimization problem as described in (D 1) using SOSTOOLS, by following the steps given below

Appendix E. Details of Numerical Experiments for Flow Structures

In the following, we describe the details of the numerical experiments carried out to obtain the flow structures for the plane Couette flow and the plane Poiseuille flow. We begin by describing the linearized Navier-Stokes equation and its corresponding discretization [START_REF] Farrell | Stochastic forcing of the linearized Navier-Stokes equations[END_REF]).

The non-dimensional linearized Navier-Stokes equations governing the evolution of disturbances in steady mean flow with streamwise velocity varying only in the crossstream direction are

Re ∆η,

where U (y) is the mean streamwise velocity component, v is the cross-section pertur-

We are interested in studying the induced L 2 -norms from inputs d x , d y , d z to u = (u x , u y , u z ) . The following corollary provides the induced norms of interest.

Corollary F.2. For any streamwise constant channel flows with nominal velocity U (y), we have

Proof: From (F 1), we infer that   

Thus, we have

Then, multiplying both sides of the above equality by the transpose of vector

. (F 8)

In order to see the influence of each

, we set d y = d z = 0 obtaining

.

It suffices to set f 1 (k z ) = h 2 xx (k z ). Similarly, we have