Beyond optimal disturbances: a statistical framework for transient growth

Abstract The theory of transient growth describes how linear mechanisms can cause temporary amplification of disturbances even when the linearized system is asymptotically stable as defined by its eigenvalues. This growth is traditionally quantified by finding the initial disturbance that generates the maximum response at the peak time of its evolution. However, this can vastly overstate the growth of a real disturbance. In this paper, we introduce a statistical perspective on transient growth that models statistics of the energy amplification of the disturbances. We derive a formula for the mean energy amplification and spatial correlation of the growing disturbance in terms of the spatial correlation of the initial disturbance. The eigendecomposition of the correlation provides the most prevalent structures, which are the statistical analogue of the standard left singular vectors of the evolution operator. We also derive accurate confidence bounds on the growth by approximating the probability density function of the energy. Applying our analysis to Poiseuille flow yields a number of observations. First, the mean energy amplification is often drastically smaller than the maximum. In these cases, it is exceedingly unlikely to achieve near-optimal growth due to the exponential behaviour observed in the probability density function. Second, the characteristic length scale of the initial disturbances has a significant impact on the expected growth, with large-scale initial disturbances growing orders of magnitude more than small-scale ones. Finally, while the optimal growth scales quadratically with Reynolds number, the mean energy amplification scales only linearly for certain reasonable choices of the initial correlation.


Introduction
A natural approach for analyzing the stability of a steady fluid flow is to linearize and calculate the eigenvalues of the linearized Navier-Stokes operator.Underlying this analysis is the assumption that there will be disturbances to the steady flow, and though their magnitude is difficult to know a priori, it is likely a value small enough that nonlinear mechanisms are not relevant.This approach, known as modal stability theory, is agnostic to the shape of any particular disturbance -if there is a positive eigenvalue, any disturbance will grow, otherwise, any disturbance will decay asymptotically.However, the modal approach predicts stability when experiments tell us otherwise.Famously, Reynolds found that at high velocities, pipe flow transitions to turbulence [25].Efforts to ground this instability in modal theory floundered: pipe flow has all stable eigenvalues.The same is true for Couette flow as well as plane Poisseuille flow at low Reynolds numbers; these flows have only stable eigenvalues, but are observed to transition [32].
The key to their instability can be, in fact, a linear mechanism [28].Perhaps counterintuitively, a linearized Navier-Stokes operator with all stable eigenvalues can lead to short-term growth in the magnitude of disturbances before they decay at the rate prescribed by the least stable eigenvalue.This transient growth is possible only when the linearized Navier-Stokes operator is non-normal, i.e., its eigenvectors are not orthogonal.This permits one eigenvector to initially subtract from another, but this cancellation can cease if one eigenvector vanishes faster than the other, leading to growth.The magnitude of this growth can be remarkable -often more than one-thousand-fold at its peak [37].While the initial disturbances are assumed to be too small for nonlinear effects to be important, when they are amplified by three orders of magnitude, the assumption of linearity may no longer hold, and nonlinearities may bring the flow away from the laminar steady state.Rather than leading directly to transition, the nonlinearities activated by the amplified disturbance might bring the flow to a new state.Instability in this state, known as secondary instability, is more likely than the primary growth to lead to turbulence [27].
The metric reported in the literature to quantify transient growth is the ratio of kinetic energy of the maximally amplified disturbance to its initial kinetic energy.This metric is usually referred to as G(t), though in this paper we call it G opt (t) to distinguish it from suboptimal and mean growths.Significant effort has been devoted to studying G opt (t) both analytically and numerically.
In channel flow, it can be shown to have quadratic dependence on the Reynolds number Re when the product of the streamwise wavenumber and Reynolds number is small, αRe << 1 [7].Under the same conditions, the time at which the maximum occurs increases linearly with Re.Indeed, numerical experiments show that there is quadratic scaling in the optimal growth and linear scaling in the optimal time for plane-Poisseuille [37], Couette [37], Blasius boundary layer [2,9], and pipe [29] flows.In all of these cases, the optimal streamwise wavenumber α is zero or very small, and the optimal spanwise wavenumber β is order unity [27].
Minimal seed theory [14] provides a nonlinear analog of optimal transient growth analysis.At each initial energy, it identifies the disturbance that achieves the greatest growth when evolved according to the full nonlinear Navier-Stokes equations.When the initial energy is just large enough that the optimal disturbance leads to sustained turbulence, the disturbance is called the minimal seed -the smallest disturbance leading to transition.Though minimal seeds are initially amplified by linear mechanisms [22], they can differ substantially in shape from the optimal disturbances in linear transient growth [23].This gives a lower bound for the energy level that disturbances must achieve to spark transition.
In either the linear or nonlinear context, considering optimal disturbances gives an upper bound on the growth experienced by disturbances, but we propose that a more complete picture of the possible growth is needed.In linear transient growth, only the optimal initial disturbance experiences G opt growth.Indeed, if the initial disturbance were one of the eigenvectors of the linearized Navier-Stokes operator, it would decay monotonically.Of course, real disturbances to the flow will not exactly coincide with the optimally amplified disturbance, so in order to quantify their growth, one needs to explore the space of suboptimal disturbances.Is most of this space inhabited by disturbances that decay or by ones that grow?Is the growth of real disturbances on the order of G opt , on average?What is the probability that a random disturbance will come close to G opt ?Motivated by these questions, we investigate transient growth from a statistical perspective in this paper.A statistical view serves both to model the uncertainty and variation in the spatial form of initial disturbances and to fully explore the high dimensional space that these disturbances occupy.We derive an equation for the mean energy of the amplified random disturbances, and dividing this by the mean initial energy gives a metric for the mean energy amplification, which we term G mean .This depends on the statistics of the incoming disturbances, and the formula we report for G mean involves the correlation matrix of the initial disturbances.The correlation matrix at time t can also be derived in terms of the initial correlations.Its eigendecomposition can be viewed as a particular variant of proper orthogonal decomposition and provides the most statistically prevalent structures, which serve as the statistical analog of the left singular vectors of the evolution matrix.
Quantifying the likelihood that a disturbance grows beyond a particular level requires knowledge of the probability density function (PDF) of the energy amplification.Whereas the mean energy amplification depends only on the correlation matrix of the incoming disturbances, the entire distribution of incoming disturbances is needed to calculate the PDF of the growth.Moreover, there is no general formula relating the two.However, we observe empirically that the PDF is nearly exponential, and this leads to an approximation strategy for it.We use the approximate PDF to derive accurate confidence bounds on the growth, i.e., energy levels which p% of the disturbances do not exceed, for some desired p.The exponential behavior of the PDF also means that if G mean is significantly below G opt , it is extremely unlikely for an initial disturbance to achieve near-G opt growth.
Throughout the paper, we demonstrate the statistical framework using plane-Poisseuille flow.Equipped with a statistical lens, numerous observations readily emerge.At each wavenumber pair (α, β), the correlation length in the wall-normal direction has a dramatic impact on G mean , with correlation lengths on the order of the channel half-height growing to nearly half of G opt .If, however, the correlation length is short compared to the channel half-height, G mean can be orders of magnitude smaller than G opt .G opt and G mean achieve their maximum values at similar locations in wavenumber space, but the peak is substantially narrower in α for G mean .This indicates that 3-dimensional disturbances, ones which contain a range of wavenumbers, further undershoot G opt .In the three-dimensional case, G mean is a function of the three-dimensional correlation matrix.We observe that when this correlation is isotropic, G mean is roughly 2% of G opt at Re = 1000.Surprisingly, we find that G mean scales nearly linearly with Re, so the gap between it and G opt widens with increasing Reynolds number.Therefore, G opt increasingly overstates the growth of random disturbances.
Even considering disturbances near the optimal wavenumber pair (α = 0, β = 2), the probability of exceeding certain levels of growth can be extremely low.We show that the distribution of energy is nearly exponential, i.e., the probability of exceeding a particular energy level decays exponentially.Therefore, if G mean is relatively small relative to G opt , there is little chance of observing growth on the order of the optimal value.For a correlation length of one-fourth the channel half-height, fewer than 0.01% of disturbances achieve G opt /2 growth for Re = 1000.
The combined effects of the non-normality of the linearized Navier-Stokes operator and randomness have been analyzed before.In particular, Farrell and Ioannou [3] considered the linearized Navier-Stokes equations forced continuously by white-in-time noise with some spatial correlation.They showed that the expected energy, once statistical stationarity is reached, can be obtained by solving a Lyapunov equation involving the linearized Navier-Stokes operator.Our study is distinct from this work in that we instead consider the physical model of transient growth -impulsive disturbances that evolve unforced under the action of the linearized Navier-Stokes equations.This leads to a dependence on the length scales present in the initial disturbances, which we investigate extensively.The present work is also different from what has been called statistical stability [18,19].That work is concerned with the stability of the statistical state of turbulent flow, whereas our study investigates statistics of transient growth.
The remainder of the paper is organized as follows.In § 2, plane-Poisseuille flow and the numerics used to perform the calculations are described.§ 3 gives a review of transient growth.In § 4, we derive a formula for the mean energy amplification and compare it to the optimal growth for Poisseuille flow, first for disturbances at one pair of wavenumbers, then for disturbances containing a range of wavenumbers.We investigate the probability density function of the growth and detail an accurate approximation strategy for it in § 5. Finally, in § 6, we conclude the paper and offer some closing remarks.

Flow description and numerics
Plane-Poissseuille flow is the steady, laminar flow between two plates separated by 2h in the y direction.The flow is in the x direction, and the plates are infinite in both the streamwise (x) and spanwise (z) directions.It is driven by a constant pressure gradient and the streamwise velocity field is given by The flow is then non-dimensionalized by the channel half-height and the centerline velocity.Because the governing equations and base-flow are homogenous in x and z, it is convenient to take the Fourier transform of disturbances to the base flow in these directions.The associated wavenumbers in the streamwise and spanwise directions are denoted α and β, respectively.For example, the transformed wall-normal velocity is Employing the usual velocity-vorticity formulation of the linearized Navier-Stokes equations yields the following equations for the evolution of disturbances [24], 3) The Orr-Sommerfeld, cross-term, and Squire operators are (2.4b) (2.4c) Above, all quantities are non-dimensionalized, and U = U (y) is the base-flow, η is the transformed wall-normal vorticity, k 2 = α 2 + β 2 is the squared wavevector magnitude, and (•) ′ indicates a wallnormal derivative ∂ ∂y .We use the code provided in [27], which uses a Chebyshev discretization of the linearized Navier-Stokes equations (2.3) [11,24].All norms presented in our numerical results are based on the kinetic energy of a disturbance.It can be shown, by using incompressibility and Parseval's theorem, that the energy of a disturbance in the transformed velocity-vorticity coordinates is [ (2.5)

Optimal transient growth
Here, we review the linear effects responsible for transient growth in a system with all negative eigenvalues.For a more thorough review, see Ref. [28].Expressing the Navier-Stokes equations as q(x, t) = N (q(x, t)), a steady solution q(x) satisfies N (q(x)) = 0. Though their size is likely small, disturbances to the base flow are inevitable.Denoting these disturbances as q(x, t), their dynamics are analyzed by linearizing around the base flow, where A is the Jacobian around the base flow, The problem is discretized as q(t) = Aq(t), where q(t) ∈ R N is the discretized state vector describing the disturbance.
The solution to (3.5) is q(t) = M t q(0), where the evolution operator is the matrix exponential If all of the eigenvalues of the linear operator A have a negative real part, then the linear system is stable in the sense that the norm of any disturbance will eventually decay, i.e., lim t→∞ q(t) = 0.This sense of stability, usually referred to as modal stability, is mathematically powerful -it is a property of the system, not of any particular disturbance.If the eigenvalues are negative, any disturbance decays eventually, but if there is a positive eigenvalue, any disturbance arising in a physical scenario will have a non-zero projection onto the associated eigenvector, and will thus grow exponentially.
The theory of transient growth offers the additional insight that even if all the eigenvalues are stable, if A is non-normal, i.e., its eigenvectors are non-orthognoal, the decay need not be monotonic.The eigenvectors summed together to construct an initial disturbance may mostly cancel each other initially, but because they vanish at different rates, after some time, there may no longer be cancellation, which leads to growth of the disturbance.The linear operators arising in fluids systems, especially in shear flows, can be highly non-normal [37].The ability for these systems to produce growth is quantified in the literature by the maximal amplification that a disturbance may undergo, G opt (t) ≡ max This quantity is usually referred to simply as G. Here, we have termed it G opt to specify that it is the optimal growth among all possible initial disturbances and to distinguish it from G mean , which will arise later in the paper.Its peak in time is referred to in this paper as G opt max (usually referred to simply as G max ).The norm • is based on the kinetic energy of the disturbance and can be written e(q) = q 2 = q * Wq.
W is a weight matrix (required to be Hermetian and positive-definite), and we make frequent use of the decomposition L * L = W.For later use, the inner product that induces the norm is q 1 , q 2 = q * 2 Wq 1 .It can be shown that the optimal growth may be written [24] G opt (t where σ 2 1 (•) returns the first (squared) singular value of the argument.The structures that undergo the most growth up to time t and the structures resulting from the amplification may also be obtained via the singular value decomposition of the weighted evolution operator, (3.10) The optimal output and input modes are recovered as U = L −1 Ũ and V = L −1 Ṽ, respectively.The first column of V is the initial disturbance that grows by G opt (t), and the first column of U is the structure that results.
The largest initial growth rate experienced by any disturbance can be expressed in terms of the optimal growth as .11)By expanding the matrix exponential to first order terms in t, it is easily shown that this optimal growth rate is given by the numerical abscissa [36] a opt = κ 1 (LAL −1 + (LAL −1 ) * ), (3.12)where κ 1 (•) returns the first eigenvalue of the argument.
(b) So long as the disturbance remains small enough, the linear approximation (3.2) remains valid, and the disturbance will decay to zero.However, if the growth is large enough, it can elevate a disturbance from the regime where linearity governs to one where nonlinear effects are relevant.These nonlinear effects can in turn lead the flow away from the base state, eventually causing transition.The growth can indeed be quite large, owing to the severe non-normality in the linearized Navier-Stokes operator in shear flows.Figure 2(a) shows G opt (t) for various streamwise and spanwise wavenumbers in plane Poisseuille flow at Re = 1000.For α = 0, β = 2, G opt max is nearly 200. Figure 2(b) shows G opt max for a range of wavenumbers.Streamwise-elongated structures (small α) are capable of larger growth than shorter structures (larger α).The peak in wavenumber space is at α = 0, β = 2.04, so structures of finite spanwise (z) length experience the most growth.
To motivate the remainder of this paper, we show 1000 random trajectories along with G opt (t) at Re = 1000, α = 0, β = 2 in Figure 3. G opt indeed bounds the trajectories, but, notably, they all substantially undershoot it.The details of the distribution used to generate Figure 3

Expected energy amplification
In light of Figure 3, an obvious question is: how much energy, on average, do the amplified disturbances achieve?We derive a formula for the mean energy of the amplified disturbances in terms of the correlation matrix of the initial disturbances.The expected energy divided by the expected initial energy is termed G mean .We elaborate on the difference between this and the expected value of the ratio of these energies at the end of the following subsection, but, in short, G mean is more physically meaningful, produces a simpler mathematical result, and requires less a priori knowledge of the initial disturbances.
Just as in the standard treatment of transient growth, the physical model that we consider consists of the discretized base flow q, which is impulsively perturbed at t = 0 by q(0).As before, the disturbance q(0) may represent the entire (three-dimensional) flow in space, the flow at a particular pair of wavenumbers, or the flow at a particular location in the streamwise direction.The evolution of the disturbances is governed by the Navier-Stokes equations linearized around the base flow, so (3.5) holds.However, our statistical framework differs from the standard treatment of transient growth in that the disturbance q(0) is now a random variable with some distribution, and we study the statistics of the disturbance after some time q(t).In particular, we are interested in the energy of the growing disturbance in comparison to that of the initial disturbance.
Experimenting with various choices of the initial correlation for plane-Poisseuille reveals that the expected energy can be substantially smaller than G opt max .This is especially true when the correlation length is short relative to the channel half-height.Furthermore, G mean max drops off more rapidly with larger α than does G opt max , which causes the mean energy amplification for three-dimensional disturbances to be quite small unless their energy is focused sharply at α = 0. Surprisingly, we observe that for isotropically correlated initial disturbances, the mean energy amplification scales near-linearly with Re, in contrast to the quadratic scaling of G opt max .

G mean
For simplicity, we omit the weight matrix in the derivations (by setting it to the identity), reporting the formulae with it at the end, so e(q(t)) = q * (t)q(t).The energy may alternatively be written as the trace of the outer product, because the diagonals of q(t)q * (t) are the terms summed in the inner product.In terms of the evolution operator, (4.1) becomes The expected value of this expression gives the expected energy of the amplified disturbances, The expectation commutes with the trace and evolution matrices, giving The expectation of the outer product of the initial disturbances is their correlation matrix, so the expected energy of the growing disturbance is expressed in terms of the correlations of the initial disturbances, A metric for the expected growth of the disturbances, which we term G mean , is provided by the ratio of the expected energy and initial energy, This quantity is not the same as the expected value of the growth; this difference is discussed at the end of this subsection.If a weight matrix W = L * L is used to define the energy, then (4.7) becomes G opt (t) is given in terms of the SVD of the evolution operator.To express G mean (t) in a similar manner, we make use of the fact that the trace of a matrix is the sum of its eigenvalues and that the eigenvalues of BB * are the squared-singular values of B for any matrix B. Using these two facts, (4.7) can be written where B is defined by the factorization C 00 = BB * .In the case of a weight matrix, (4.9) becomes Upper and lower bounds for G mean (t) for any possible initial correlation can be obtained by setting C 00 to the outer product of the first input mode with itself and last input mode with itself, i.e., v 1 v 1 * and v N v N * , respectively, yielding the bounds σ 2 1 (LM t L −1 ) and σ 2 N (LM t L −1 ).Notably, the upper bound is G opt (t).In the case that the disturbances are white in space, i.e, C 00 = W −1 , the resulting G mean (t) is the mean-squared singular value of the weighted evolution operator LM t L −1 .
G mean , defined in (4.7), is the ratio of the expected energy of the disturbance at time t to its expected initial energy.This is distinct from the expected ratio of energy, E e(q(t)) e(q(0)) .Physically, the ratio of expected energies is the more salient quantity because whether a particular disturbance leads to transition depends on its final energy (and shape), not on the growth it underwent.In Figure 3, this ratio of expected energies is the mean of the gray curves (at each time).The expected ratio of energies would come from dividing each disturbance by its initial energy, then taking the average, but this inappropriately weights the growth of smaller initial disturbances equal to that of larger ones.Mathematically, the ratio of expected energies is the easier quantity to work with because it depends only on the correlation matrix of the initial disturbances, as shown in (4.7), while the expected ratio of energies depends on the entire distribution of the initial disturbances.If there is no variation in the size of initial disturbances, i.e., if they live on an N -dimensional sphere, the two quantities are the same.More generally, the quantities are the same in the case that the distribution of initial disturbances is separable in radius and direction, as is proven in Appendix A.
Analogous to the numerical abscissa a opt , we define the mean initial growth rate, This derivative can be calculated by expanding the evolution operator to first order, where I is the identity.Dropping the quadratic term and evaluating the derivative gives Finally, leveraging the Hermicity of the correlation matrix, where Re(•) returns the real part of the argument.The upper bound for this quantity is a opt , which is positive if (and only if) G opt > 1, but we have never observed a opt to be positive in our numerical experiments.Indeed, we have never observed a randomly chosen disturbance initially grow.

Correlation and dominant structures
The statistics of the initial disturbances can also be used to augment prediction of the structures that arise from the linear amplification by the evolution operator.Removing the trace from (4.6) gives a formula for the correlation matrix of the disturbance at time t, The dominant flow structures at time t are the eigenvectors of this correlation matrix (multiplied by a weight if desired), The columns φ k t of Φ t are orthogonal in the weighted inner product, i.e., φ i t , φ j t = δ ij .This can be thought of as a particular variant of proper orthogonal decomposition (POD) [17,16,30] in which the data consists of an ensemble of realizations of the disturbances at a specific time t rather than a single time series.The eigenvalues are non-negative, owing to the semi-positive definiteness of the correlation matrix, and represent the expected energy of each structure.More precisely, the k-th eigenvalue is the expected energy of the projection of the disturbance onto the k-th mode φ k t .The eigenvalues sum to the total expected energy, so Therefore, the eigenvalues quantify the expected contribution of each mode to the growth of the disturbance.
The average energy of the disturbance captured by any structure can be quantified [5] by, The first POD mode maximizes this quantity (over normalized modes), and the latter modes maximize it with the constraint that they are orthogonal to all previous ones.For a more thorough review of POD, see Refs.[26,31,35].
The POD modes offer an alternative to the output modes of the evolution matrix for describing the structures that emerge from the linear amplification.The POD modes are the most energetic structures, while the output modes are the modes resulting from the greatest amplification by the evolution operator.In the case that the initial correlation C 00 is white with respect to the weight, the POD modes are equivalent to the output modes, i.e., This result is analogous to the relationship between resolvent modes and spectral POD modes established in Ref. [35].Of course, the initial correlation is unlikely to be white in a real flow, so it is advantageous to use knowledge of the incoming statistics to augment the prediction of these structures.
In the remainder of this section, we experiment with different choices of C 00 for Poisseuille flow and record our observations.We are not aware of any previous studies on the correlations of disturbances within Poisseuille flow.The nature of the disturbances and quantities such as their correlations are certainly sensitive to the specifics of the flow setup.For example, the disturbances generated by vibrations of the boundary are likely substantially different from those caused by surface roughness.Providing a model for the correlations of the initial disturbance is not the topic of this paper, and we do not claim that the choices made are necessarily reflective of the physics in Poisseuille flow.However, trends that emerge, e.g., that longer correlation lengths lead to more growth and that G mean is substantially smaller than G opt , are not specific to our choice of the correlation, and therefore give physical insight despite the current lack of an accurate model for the correlations.

Numerical experiments with disturbances at a single wavenumber pair
Most studies of transient growth in flows with homogeneous directions take the Fourier transform in these directions and calculate the transient growth for disturbances consisting of a single pair of streamwise and spanwise wavenumbers.Here we perform the analogous analysis for G mean in Poisseuille flow.The correlation at a particular α and β can be written Ĉ00 (y 1 , y 2 ; α, β) = Ĉvv 00 (y 1 , y 2 ; α, β) Ĉvη 00 (y 1 , y 2 ; α, β) Ĉηv 00 (y 1 , y 2 ; α, β) Ĉηη 00 (y 1 , y 2 ; α, β) where the diagonal terms are the autocorrelations of wall-normal velocity and wall-normal vorticity, and the off diagonal terms are the cross correlations between these two variables.It can be shown analytically that for a disturbance to experience large growth, its initial energy should be concentrated in its wall-normal velocity rather than wall-normal vorticity [7].Therefore, we choose only the vertical velocity autocorrelation to be non-zero and take it to be Gaussian in the wall-normal direction with correlation length λ, i.e., Ĉvv 00 (y 1 , y

.22)
The normalization A has no impact on G mean because this constant affects the expected energy of the amplified disturbances and that of the initial ones equally.In our numerics, it is chosen so that when the initial correlation is discretized in y, its trace is unity.

G mean for a single wavenumber pair
Figure 4 shows G mean (t) (solid) for various wavenumbers and G opt (t) (dashed) for the same wavenumbers, both as functions of time for Re = 1000.Whether the mean is on the same order as the maximum depends on the characteristics of the correlations of the initial disturbances.We refer to the peak of G mean (t) in time as G mean max .For the relatively long correlations in (a), G mean max is roughly half G opt max for the most amplified wavenumbers, while for the shorter correlation length (b), the ratio is closer to one-tenth.Figure 5 shows the first time unit of G mean (t) using the same parameters as Figure 4. Despite the fact that G mean grows to be relatively large, it initially decays sharply for all wavenumbers.The initial decay rate can be calculated with (4.14).
Figure 6(a) shows G mean max for a range of λ −1 at Re = 1000.The correlation length λ greatly impacts the mean energy amplification, with longer correlation lengths corresponding to more growth and shorter ones to less growth.It is likely this trend is explained by the fact that short-wavelength (in y) disturbances are quickly dissipated by viscosity before they can extract energy from the mean shear [21].Figure 6(b) shows the time at which G mean is maximized.This time is relatively independent of the correlation length, but changes substantially with the wavenumber pair.
The wavenumber dependance is further explored in Figure 7. Figure 7 (a) shows this dependence for λ −1 = 1, which is near the peak for the maximally amplified wavenumber in Figure 6 (a).The location of the peak in wavenumber space is near that of G opt max seen in Figure 2; however, G mean max decays much more rapidly with α than does G opt max .This indicates that to achieve large- Figure 7: Dependence of G mean on the streamwise and spanwise wavenumbers at Re = 1000 for different correlation lengths.The shape is similar to the contour of G opt max at the same Reynolds number (Figure 2), but, notably, the support in α is substantially narrower for G mean .This indicates that the energy of the disturbance must be quite concentrated at the large-growth wavenumbers to achieve significant growth.quadratically with Reynolds number (for small values of αRe) [7].In Figure 9, we show the scaling of G mean max with Reynolds number for a variety of correlation lengths at α = 0, β = 2.These appear to obey the same scaling.

Dominant structures for a single wavenumber pair
Now we examine the structures that emerge in Poisseuille flow at a single wavenumber pair, as described in § 4.1.2.The key question is: to what extent do the output modes resemble the principal components of the correlation matrix, i.e., the POD modes?The former are the structures resulting from the largest amplification by the evolution operator (see (3.10)), whereas the latter use the statistics of the initial disturbances to inform which structures are most energetic (see (4.16)).

Sfrag replacements
Figure 10: Evolution of correlations and their POD modes for Re = 1000, α = 0, β = 2, and λ = 1.The initial energy is imposed to be in the velocity (top), but it is quickly shifted to vorticity (middle).The first POD mode and first output mode quickly become similar (bottom).For these parameters, G mean (t) peaks near t = 80 (see Figure 4) .
Figure 10 shows the evolution of the correlations and the vorticity component of their POD modes for α = 0, β = 2.We impose the initial correlation to be of the form in (4.22) with λ = 1, so all the energy is initially in the velocity.The evolution operator rapidly shifts this energy to the vorticity, and two counter-rotating vortices emerge.The first POD mode (blue) reflects this with two peaks of opposite sign at the peaks in the vorticity correlation.The velocity component of the POD mode is not plotted because it rapidly decays to 0. Notably, the first output mode quickly becomes nearly identical to the first POD mode despite the former not depending on the initial correlation.Indeed, for this wavenumber pair, the first few POD modes from different initial correlation matrices quickly become similar to one another and to the first few left singular values of the evolution matrix.There is only moderate gain separation in the singular values, so the similarity between the modes is surprising.
Figure 11 compares the modes more thoroughly.For each wavenumber pair, Figure 11(a) shows the average energy captured (see (4.19)) by the first output modes at the peak time of G mean as a fraction of that captured by the first POD mode at the same time.The energy captured by the POD mode is the maximum possible, so a value near unity indicates that the output mode is very effective, while a value near zero indicates the opposite.Whether or not a structure is visible in the flow depends on the energy it captures.Near α = 0, the output mode captures nearly as much as is possible, while at higher α, it captures substantially less.Capturing energy at low α is more important as the growth is the greatest here, so the leading output mode does a good job of predicting the structures observed in Poisseuille flow [1].The large discontinuity in Figure 11(a) in the lower right of the plot occurs because the peak time for these wavenumbers is t = 0, so the structures at this time only depend on the initial correlation, not on the linearized Navier-Stokes operator, so the output modes here cannot hope to capture any structure.
Figure 11(b) compares the output modes to the POD modes via the inner product, again at the peak time in G mean .At each wavenumber pair, the plot shows the matrix of square inner products G ij = | u i , φ j | 2 up to three modes in each basis.For α = 0, this matrix is nearly diagonal, indicating, once again, that the output modes and POD modes are very similar, even for the subleading modes.At higher α, the modes become less similar as shown by the off-diagonal terms in the matrix of inner products.
We also experimented with correlation matrices that do not respect the symmetry about y = 0 in the channel (arising, e.g., due to vibrations in one of the plates, but not the other).These resulted in less similarity in the modes, and the square inner products were in the range 0.5 − 0.9, even at α = 0.However, so long as the correlation was symmetric, the modes were quite similar for α = 0.The output modes are therefore a good model for the POD modes under these conditions, and leveraging the statistics may not present much advantage in predicting the structures.We stress, however, that the energy of each structure is highly dependant on the statistics, so the SVD of the evolution operator does not provide the associated energies accurately.Whether the POD modes and output modes for other flows coincide to the extent that they do in Poisseuille flow may be an interesting topic for future investigation.
It can be shown that the POD modes in flows with homogeneous directions are still delta functions in wavenumber space in those directions [16].Therefore, despite the fact that the disturbances will not themselves be delta functions in wavenumber, the similarity in the POD modes and output modes observed in this subsection still applies in the three-dimensional case.However, the behavior of G mean for three-dimensional disturbances can be markedly different, as shown next.

Three-dimensional disturbances: inclusion of multiple wavenumbers
Clearly, just as real initial disturbances will not identically match the maximally amplified one, real disturbances do not exist at just one pair of streamwise and spanwise wavenumbers.Parallel flow offers an analytical simplification to an analysis of transient growth -each streamwise and spanwise wavenumber pair may be considered separately in its ability to produce growth.However, this tempts further exaggeration of the possibility for large-scale transient growth.In modal stability, one need not add these wavenumbers back together to get an answer as to the long-term behavior -if any wavenumber pair grows exponentially, so will the entire disturbance.However, if one wavenumber pair experiences large transient growth, this only implies a large gain for the entire disturbance to the extent that its initial energy is concentrated at that wavenumber pair.In this subsection, we incorporate a range of wavenumbers and show that the weight for each pair is determined by the Fourier transform in x and z of the three-dimensional correlation.When these three-dimensional correlations are incorporated, substantially less growth is observed.We also observe a linear scaling with Reynolds number for an isotropic correlation in contrast with the quadratic scaling observed for G opt max and for G mean max at a particular α and β.
One way to calculate G mean for disturbances containing multiple wavenumbers would be to define a large domain in x and z (to approximate the desired infinite directions), calculate A, and define a discrete correlation matrix for the full three-dimensional problem.The mean energy amplification would then be given by (4.7) and the correlation matrix by (4.15).However, this strategy is needlessly computationally intensive because it does not take advantage of the analytic simplification possible in parallel flow.Instead, we can add the expected energies at each wavenumber together, modulated by the energy of the incoming disturbances at each wavenumber.Denoting Figure 12: G mean max for a three-dimensional isotropic correlation with correlation length λ at Re = 1000.Even at the optimal correlation length, the inclusion of all wavenumbers causes G mean max to be roughly 2% of G opt max at α = 0, β = 2.
The exponential term can be interpreted as the expected energy at each wavenumber pair implied by the correlation.For Poisseuille flow at Re = 1000, the most amplified wavenumbers are near α = 0, β = 2 (see Figures 2 and 7), and the amplification drops off rapidly as α moves away from zero.To concentrate energy near α = 0, the correlation length must be quite long.This longer correlation also promotes growth because, as we described in the previous section, the longer the correlation length in y, the more growth is observed.However, with a long correlation length, energy is concentrated near β = 0, which does not experience much growth (see Figures 2 and 7).Also detracting from G mean is the fact that the maxima occur at significantly different times for different wavenumbers (see Figure 6).The combined effect of these factors can be seen in Figure 12.At Re = 1000, even with the correlation length that promotes the most growth (λ −1 = 1.2),G mean max is only 2.5% of G opt max at the optimal wavenumbers.
Figure 13 shows contours of G mean max for a range of correlation lengths and Reynolds numbers.The vertical dashed line divides the asymptotically stable and unstable regions.At Reynolds numbers higher than this, G mean max is technically infinite, but there is an initial peak in G mean (t) long before the instability dominates.In this figure, and all subsequent ones which show G mean max above the critical Reynolds number, we plot the magnitude of the initial peak in G mean .The effect of the correlation length is similar across Reynolds numbers.Comparing Figure 13 with Figure 8 (its single-wavenumber analog), we see that the isotropic correlation matrix severely limits growth at all Reynolds numbers.Indeed, the difference becomes greater as the Reynolds number increases.
The Reynolds number scaling is shown in Figure 14.Unlike G opt max or G mean max at a particular wavenumber pair, G mean max for an isotropic three-dimensional correlation scales nearly linearly with Reynolds number.This is a surprising result -the three-dimensional G mean (t) is obtained in (4.26) by integrating over G mean (t) at particular wavenumbers.However, critically, these singlewavenumber values for G mean peak at different times, as can be seen in Figure 6.The difference in the scaling means that the difference between G mean max and G opt max becomes larger with Reynolds num- Figure 13: G mean max for an isotropic initial correlation for a range of correlation length and Reynolds number.Similar dependence on correlation length λ is observed at all Re, and all values are substantially lower than their single-wavenumber counterparts in Figure 8. ber, i.e., G opt max increasingly overpredicts the mean energy amplification with increasing Reynolds number.

Non-isotropic correlation
As a generalization of the isotropic correlation investigated above, we next consider the ellipsoid where λ x , λ y , and λ z are the correlation lengths in the streamwise, spanwise, and wall-normal directions, respectively.For example, this allows for the correlations to persist longer in x than in y or z, as may result from the advective nature of the flow [10].With this extra freedom relative to the isotropic case, we may ask whether the Reynolds number scaling remains linear, as it is in that case, or becomes quadratic, as it is for G opt .The answer depends on the correlation lengths chosen, but we find that if we fix λ x at some non-zero value and vary Reynolds number, the scaling is linear.
When discretized in y, the ellipsoid correlation becomes where, again, C y 00 is the discretized y-dependant part.Upon taking the Fourier transform, the correlation in wavenumber space is and various λ x .The scaling is initially quadratic but becomes linear at higher Re.Longer correlations in x remain quadratic to higher Re.
quadratic up to higher Reynolds numbers.

Estimating the PDF
As Figures 6 and 12 show, the expected value of the energy may be orders of magnitude smaller than G opt max for certain reasonable incoming correlations.With a formula for the expected energy (4.7), one may still ask whether large gains are possible or negligibly unlikely.In this section, we discuss methods to estimate the probability density function (PDF) of the energy.Whereas the mean energy depends only on the initial correlations, the PDF depends on the entire distribution of the incoming disturbances.We first describe a basic Monte Carlo approach for estimating the PDF and apply this to two candidate distributions for the incoming disturbances, noticing that the PDF of the energy drops nearly exponentially.In the two subsequent subsections, two distributions for the initial disturbances are considered, a multivariate Gaussian and a transformation of the uniform distribution on the N -sphere.For these distributions, it is possible to analytically approximate the PDF of the energy using the exactly calculable moments of the energy distribution.With an accurate estimate of the PDF, we can calculate percentile curves for the trajectories.
We denote the probability density function of a random variable X as f X (x) ≡ lim dx→0 P r{X ∈ [x, x + dx]}/dx, where P r{•} denotes probability.If X ∈ R N is a vector, then the PDF is defined (5.1)The incoming disturbances follow some distribution q(0) ∼ f q(0) (q 0 ), ( and this implies a distribution q(t) ∼ f q(t) (q t ) (5.3) for the disturbances q(t) = M t q(0) at time t.These distributions are the most descriptive statistical information about the disturbances; any statistic of the disturbances is implied by the full distribution.For example, the correlation matrix C 00 for the initial disturbances is implied by the distribution of the initial disturbance f q(0) .The converse, however, is not true -there are many distributions with the same correlation matrix.In the last section, we showed that G mean only depends on the correlation matrix, so there was no need to consider the form of the underlying distribution.However, to calculate the PDF of the energy of the disturbance at some time, e(t) ∼ f e(t) (e t ), (5.4) the full distribution of initial disturbances is needed.

Monte Carlo
With a means of sampling initial disturbances from f q(0) (q 0 ), samples of the growing disturbances can be generated by multiplying the initial ones by M t , and samples of e(t) are finally obtained by computing their norm.An estimation of the PDF can be obtained using standard methods, such as ksdensity in Matlab.Figure 17 shows the empirical PDF resulting from performing this Monte Carlo with two different distributions of the initial disturbances (described later), both with the same correlation.The Monte Carlo is performed at t = 100, Re = 1000, α = 0, β = 2 with a correlation of the form (4. in the following subsections.Two observations are apparent.First, both distributions result in a very similar PDF for the energy of the amplified disturbance.That the PDFs are similar indicates that while the PDF of the energy is a function of the full distribution of incoming disturbances, reasonable distributions of initial disturbances with the same correlation will give similar PDFs for the energy.Second, the PDFs decay nearly exponentially.The exponential decay indicates that it is very unlikely that the energy of an amplified disturbance substantially exceeds G mean (t).The exponential decay also allows for accurate a priori approximation of the PDF, which we discuss in the following subsections for the two distributions of initial disturbances.

Multivariate Gaussian
First, we assume the initial disturbances follow the multivariate Gaussian with mean 0 and correlation C 00 , q(0) ∼ N (0, C 00 ).
(5.5) Also, we assume the correlation to have unit trace, so G mean (t) = E[e(q(t))].More explicitly, this distribution, as represented by its PDF, is If C 00 is rank-deficient, the inverse C −1 00 and determinant |C 00 | are modified to the pseudoinverse and pseudodeterminant.
Any linear function of q(0) also follows a multivariate Gaussian distribution [33], so the disturbance q(t) = M t q(0) some time in the future is distributed as q(t) ∼ N (0, C tt ), (5.7)where the correlation C tt = M t C 00 M * t comes from (4.15).The energy is e(t) = q(t) 2 and we seek to estimate its PDF f e(t) (e t ).This distribution is one of a well-studied class -quadratic forms in multivariate Gaussian variables.The moments of these distributions are known [20], and for the case at hand, the r-th moment µ r can be calculated recursively as where and µ 0 = 1.Note that µ 1 = G mean , i.e., the first moment recovers the expected energy.
Our goal is to estimate the right tail of the PDF of the energy in order to approximate the probability of exceeding a particular energy.As shown in Figure 17, the right tail of the empirical PDF displays nearly exponential decay.Therefore we assume its form to be To find the decay rate, we find γ such that the r-th moment of the exponential ansatz matches the r th moment of the true distribution, given in (5.8).The true distribution, estimated via the Monte Carlo, is near-exponential for high energies but not low ones, so to find the correct exponential parameter we equate a relatively high moment, since this weights the high-energy tail of the distribution heavily.Denoting the moment equated as r, the r-th moment of the exponential distribution (5.9) is µ r = r! γ r .Equating this to the true r-th moment given in (5.8) and solving for the exponential decay rate gives This is an analytical approximation; the only role of the previous Monte Carlo was to suggest the exponential form (5.9). Figure 18(a) shows this approximation strategy using r = 4.The confidence bounds are derived by integrating the approximate PDF.The Reynolds number, wavenumbers, and correlation length are Re = 1000, α = 0, β = 2, λ −1 = 5.Out of the 10 7 trajectories used to generate the empirical distribution, 99.005% were below the 99% confidence bound.

Transformation of a uniform distribution on the N-sphere
Second, we assume that the disturbances are distributed as some transformation of a uniform distribution on the surface of the N -sphere, where u = 1 with probability 1, and u ∼ Ru for any rotation matrix R. Samples from this distribution can be easily generated by normalizing samples from an i.i.d.multivariate Gaussian to be the same radius R, u = R x x with x ∼ N (0, I).The transformation Ψ can be chosen so that the initial disturbances have some desired correlation C 00 .By choosing Ψ such that  e.g., setting it to the Cholesky decomposition of N C 00 , where N is the dimension of the state, the correlation matrix of the initial disturbances is The first equality holds because E[uu * ] = 1 N I.The disturbance energy some time later is given by e(t) = q(t) 2 = Tu 2 .The final equality expresses the energy as the square norm of the uniform distribution acted on by a matrix T = M t Ψ.
The moments for this transformation of the uniform distribution were derived by von Neumann [38].The result is [12] E where the ξ are defined as the power series coefficients of and ζ j is defined in terms of the singular values σ i (LT), as Using these moments, we can approximate the true PDF using the technique described in the previous subsection.Figure 18 shows the result along with confidence bounds derived by integrating the approximated PDF.Just as before, the Reynolds number, wavenumbers, and correlation length are Re = 1000, α = 0, β = 2, λ −1 = 5.Of the 10 Figure 19: The one-thousand trajectories in Figure 3 overlayed with the calculated mean and percentile curves.The Reynolds number is 1000, the wavenumbers are α = 0, β = 2, and λ −1 = 5 with the correlation given in (4.22).The mean drastically undershoots the optimum, and the energy is exponentially distributed, so it is unlikely for a disturbance to grow by near G opt .

Conclusions
Standard transient growth analyses are based on the maximum growth experienced by any initial disturbance.While this is a useful upper bound on linear growth, it can vastly overpredict the growth experienced by real disturbances.We have developed a statistical framework to explore the space of real disturbances and quantify their growth.We demonstrated the framework and its ability to extract insights on Poisseuille flow.
The framework can be summarized using Figure 19.G opt , the quantity used in the literature to quantify transient growth, far overshoots all one thousand random trajectories, which are the same ones shown in Figure 3. G mean (t), gives the mean energy divided by the expected initial energy, i.e., the mean energy amplification.It is a function of the correlation matrix of the initial disturbances but does not depend explicitly on the particular form of the distribution of initial disturbances.As Figure 19 shows, G mean (t) can be significantly lower than G opt (t).For threedimensional disturbances, the three-dimensional correlation matrix determines the average energy at each wavenumber.With a realistic correlation matrix, suboptimal wavenumbers will account for a significant portion of the energy, and this, coupled with the fact that different wavenumbers peak at different times, further widens the gap between G mean and G opt .The confidence bounds in Figure 19 give the energy levels that p% of the trajectories undershoot.The levels are calculated analytically by integrating the PDF of the energy, not by performing a Monte Carlo.The energy PDF cannot be calculated exactly for a general distribution of initial disturbances; however, because the energy PDF is nearly exponential, it can be approximated accurately.
Applied to Poisseuille flow, this statistical view reveals a number of insights.For a single wavenumber pair, the correlation length in the wall-normal direction emerges as an important parameter in determining G mean .For long correlation lengths, G mean is on the same order as G opt (as much as half for certain wavenumbers), but for short correlation lengths, G mean is orders of magnitude smaller than G opt .The dependence of G mean on the streamwise and spanwise wavenumbers is different than that of G opt : while the peak is still near α = 0, β = 2, it is substantially narrower in α.Three-dimensional disturbances contain energy at all wavenumber pairs, so the narrower peak of G mean in α leads to substantially less growth.At Re = 1000, with an isotropic correlation of correlation length λ, G mean max is only 2% of G opt even at the most growth-promoting λ.Furthermore, for this form of the three-dimensional correlation, G mean max grows near-linearly, while G opt max grows quadratically, leading the latter to increasingly overpredict the mean energy amplification as the Reynolds number increases.
The formulae derived depend on the correlation matrix of the initial disturbances.We are not aware of previous studies on the statistics of these disturbances in Poisseuille flow, but they are likely to be dependent on the source of disturbances.The figures reported, e.g., for the ratio of G mean /G opt , are not meant to be taken as quantitative predictions of what would be observed in an experiment.Rather, they are meant to show trends and to emphasize that, when various factors are accounted for, the growth of real disturbances can be significantly smaller than G max .Though G mean depends on the full correlation matrix, the correlation length is a particularly important feature in determining G mean .The observation that long correlations lead to more growth may serve as a practical guide when an accurate model for correlations is unavailable.
We have also described a statistical approach to determining the structures that emerge, in the form of modes from a particular POD problem.These structures also depend on the initial correlations but were observed only to differ slightly from the standard output modes of the evolution operator, which do not depend on the correlations.This indicates that the output modes do a good job of capturing the energy of the growing disturbances regardless of the correlation matrix.Nevertheless, they do not accurately predict the energy of each structure, whereas the POD eigenvalues do.
We have discussed the statistical framework in the context of temporal stability, wherein an initial disturbance at a particular time is assumed, then evolved forward in time without further forcing to the linear dynamics.In spatial stability, a disturbance is introduced at a particular streamwise location, and its growth is then calculated as it evolves downstream as a function of the streamwise coordinate.Transient growth has been investigated in the context of spatial stability [8,4], and the framework developed in this paper applies equally to spatial stability by exchanging t for x and the linearized Navier-Stokes operator for a spatial evolution operator [34].Spatial stability may be an excellent application for two reasons.First, applied to spatial stability, G mean (x 1 ) represents the ratio of turbulent intensities at x = x 1 and x = 0. Consider a laminar boundary layer excited at x = 0 by free-stream turbulence.The turbulent intensity at various points along the boundary layer is critical for determining where transition occurs, and its evolution in space is given by G mean (x).In this context, we have shown that the turbulent intensity depends not only on the input turbulent intensity, encapsulated by the diagonal terms of the initial correlation matrix, but on the off-diagonal terms as well.Second, the correlation matrix of the initial disturbances may be easier to model physically in spatial stability than in temporal stability.In this case, incoming turbulence is often the source of the disturbances and turbulent correlations have been studied thoroughly.In the example of a boundary layer, the free-stream turbulence that excites the leading edge may be modeled by the von Kármán spectrum [13], which implies the correlations within the initial disturbances.Also, with such a model for the initial correlation, the POD modes of the space-evolved correlation may provide a better basis for the structures that arise downstream than the output modes of the matrix exponential.

Figure 1 :
Figure 1: Schematic of Poisseuille flow.The red waves represent disturbances with particular wavenumbers.

Figure 2 :
Figure 2: Optimal gains for plane Poisseuille flow at Re = 1000: (a) The maximal gain over all initial disturbances G opt (t) for various choices of wavenumbers.(b) Maximal gain, also maximized over time for a range of streamwise and spanwise wavenumbers α, β.

Figure 4 :Figure 5 :
Figure4: G mean (t) (solid) and G opt (t) (dashed) for various wavenumber at Re = 1000.The mean energy amplification is substantially higher for the longer correlation length λ −1 = 1 than for the shorter one λ −1 = 5. g replacements

Figure 6 :
Figure 6: The effect of correlation length on G mean max for Poisseuille flow at Re = 1000.(a) G mean (t) maximized over time vs. inverse correlation length for various streamwise and spanwise wavenumbers.More coherent disturbances (large λ) tend to grow more, but there is a non-infinite optimum.(b) The time at which G mean is maximized vs. inverse correlation length.The maximum time does not vary much with λ but does with α and β, with shorter wavenumbers corresponding to an earlier maximization time.The maximization time drops to zero when the correlation becomes such that G mean (t) never exceeds 1.

Figure 8 :
Figure 8: G mean max at α = 0, β = 2 as a function of Reynolds number and inverse correlation length λ −1 .Similar behavior is observed over the range of Re.

Figure 8 5 Figure 9 :
Figure8shows G mean opt at α = 0, β = 2 for a range of Reynolds numbers and correlation lengths.Similar dependence on correlation length is observed at all Reynolds numbers with the peak occuring when the correlation length is roughly the channel half-height.G opt max is known to scale

Figure 11 :
Figure 11: Comparison of the POD and output modes.(a) The ratio of energy captured by the first output mode to that of the first POD mode for a range of wavenumbers.The modes are comparable for low α.(b) The square inner products of the first three POD modes and the first three output modes.In both plots, the modes at each wavenumber are compared at the peak time in G mean .

Figure 15 :
Figure15: G mean max optimized over the three correlation lengths (left axis) and the optimal correlation lengths (right axis) vs. Re.The scaling of G mean max is quadratic here, and the optimal correlations do not change substantially with Re.

Figure 16 :
Figure 16: Scaling of G mean max with Re for[λ −1 y , λ −1 z ] = [1, 1.7]and various λ x .The scaling is initially quadratic but becomes linear at higher Re.Longer correlations in x remain quadratic to higher Re.

1 Figure 17 :
Figure17: Empirical PDF of energy at t = 100 resulting from initial disturbances distributed as a transformation of the uniform distribution and a multivariate Gaussian, both with the same correlation with correlation length λ −1 = 5.Note the similarity in the two PDFs and the nearexponential decay.

1 Figure 18 :
Figure18: Empirical PDF for the disturbance energy and the approximation thereof.The approximation is the exponential distribution with the same fourth moment as the true distribution and is calculated without the Monte Carlo.
7trajectories used to calculate the empirical distribution, 98.981% were less than the 99% confidence bound.