# On non-normality and classification of amplification mechanisms in stability
and resolvent analysis

Sean Symon<EMAIL_ADDRESS>Kevin Rosenberg  Scott T. M. Dawson  Beverley
J. McKeon  Graduate Aerospace Laboratories, California Institute of
Technology, Pasadena, California 91125, USA

###### Abstract

We seek to quantify non-normality of the most amplified resolvent modes and
predict their features based on the characteristics of the base or mean
velocity profile. A 2-by-2 model linear Navier-Stokes (LNS) operator
illustrates how non-normality from mean shear distributes perturbation energy
in different velocity components of the forcing and response modes. The
inverse of their inner product, which is unity for a purely normal mechanism,
is proposed as a measure to quantify non-normality. When the operator is
normal but not self-adjoint, a phase difference between the forcing and
response modes is indicated. In flows where there is downstream spatial
dependence of the base/mean, mean flow advection separates the spatial support
of forcing and response modes which impacts the inner product and leads to
non-normal amplification. Success of mean flow (linear) stability analysis for
a particular frequency depends on the normality of amplification. If the
amplification is normal, the resolvent operator can be rewritten in its dyadic
representation to reveal that the adjoint and forward stability modes are
proportional to the forcing and response resolvent modes at that frequency. If
the amplification is non-normal, then resolvent analysis is required to
understand the origin of observed flow structures. Eigenspectra and
pseudospectra are used to characterize these phenomena. Two test cases are
studied: low Reynolds number cylinder flow and turbulent channel flow. The
first deals mainly with normal mechanisms hence the success of both classical
and mean stability analysis with respect to predicting the critical Reynolds
number and global frequency of the saturated flow. Quantification of non-
normality using the inverse inner product of the leading forcing and response
modes agrees well with the product of the resolvent norm and distance between
the imaginary axis and least stable eigenvalue. An approximate wavemaker
computed from the resolvent modes scales with the length of the mean
recirculation bubble. In the case of turbulent channel flow, structures result
from both normal and non-normal mechanisms, suggesting that both are necessary
elements to sustain turbulence. Mean shear is exploited most efficiently by
stationary disturbances while bounds on the pseudospectra illustrate how non-
normality is responsible for the most amplified disturbances at spatial
wavenumbers and temporal frequencies corresponding to well-known turbulent
structures. Some implications for flow control are discussed.

††preprint: APS/123-QED

##  I Introduction

###  I.1 Background

Decomposing unsteady and turbulent flows into low-rank models is an important
step towards the realization of closed-loop flow control, which has proven to
be elusive (at least in part) due to the high degrees of freedom in most flow
systems. Despite the inherent complexity, a time-averaged flow, or mean, that
is statistically stationary can often be defined and leveraged using the
eigenvalue spectrum of the governing Navier-Stokes equations (NSE) to educe
the frequencies, i.e. the imaginary part of the eigenvalues, and shapes of
coherent structures which appear in the flow. Recent studies have demonstrated
the success of mean flow stability analysis for a variety of flows including
thermosolutal convection Turton  _et al._ (2015), turbulent jets Gudmundsson
and Colonius (2011); Oberleithner  _et al._ (2014); Schmidt  _et al._ (2017a),
and flow over a backward facing step Beneddine  _et al._ (2016). There is also
a significant body of work discussing stability analysis of the mean cylinder
wake which was shown by Barkley (2006) to correctly identify the frequency of
the globally unstable flow above the critical Reynolds number of $Re=47$
(Provansal  _et al._ , 1987; Sreenivasan  _et al._ , 1987; Noack and
Eckelmann, 1994). Notably, classical linear stability analysis of the base
flow, which is an equilibrium solution of the NSE, at supercritical Reynolds
numbers does not predict the correct observed frequency. The base (laminar)
and mean (time-average of the fluctuating velocity field) profiles are
differentiated because of the importance of nonlinearity in sustaining the
latter.

Recent work has endeavored to explain why and when mean stability analysis is
valid. Barkley (2006) suggested that success corresponds to cases where the
Reynolds stresses are unperturbed at order $\epsilon$ when considering
infinitesimal perturbations $\epsilon\tilde{\bm{u}}(x,y)\text{exp}(\lambda t)$
to the mean flow solution. This was confirmed by Sipp and Lebedev (2007) who
performed a weakly nonlinear expansion of the cylinder base flow near
criticality and determined that the nonlinear interaction of the leading
global mode with its conjugate, i.e. the contribution to the mean Reynolds
stresses, significantly outweighed the interaction of the mode with itself
leading to additional harmonics of the wavenumber. Mantič-Lugo  _et al._
(2014) provided a self-consistent framework for the cylinder mean flow by
approximating the divergence of the Reynolds stresses with the leading global
mode and its conjugate.

Sipp and Lebedev (2007) used open cavity flow as a counter example to the
validity of mean stability analysis where the predicted frequencies do not
match direct numerical simulation (DNS) of the flow. This discrepancy can be
attributed to the non-normality of the flow which leads to non-orthogonality
of the global modes and sensitivity of the spectrum to perturbation of the
operator Trefethen  _et al._ (1993). In such circumstances, the behavior of
the system can be more accurately characterized by consideration of the
properties of the pseudospectrum of the LNS operator using resolvent analysis
(e.g. Trefethen  _et al._ , 1993; Schmid and Henningson, 2001) rather than the
spectrum itself. Jovanović and Bamieh (2005) formulated the linearized problem
for laminar channel flow in input-output terms, where the resolvent operator
constitutes the transfer function between them, considering the component-wise
transfer from harmonic exogenous disturbance or forcing (input) to velocity
response (output). There is also a broad literature considering stochastic
forcing, (e.g. Farrell and Ioannou, 1993), and the initial condition,
transient growth problem, (e.g Butler and Farrell, 1992).

McKeon and Sharma (2010) and Hwang and Cossu (2010) considered the resolvent
reformulated with respect to the turbulent mean flow for canonical turbulent
wall flows. The latter authors employed an eddy viscosity to account for the
action of the Reynolds stresses, while the former analysis extends the
approach to include the nonlinear terms as the input forcing to the linear
operator, i.e. closing the feedback loop. McKeon and Sharma (2010) performed a
singular value decomposition of the resolvent to identify the inputs giving
rise to the most amplified responses which are ranked by their gain (singular
value). The approach has been extended to non-parallel flows (e.g. Lu and
Papadakis, 2014; Beneddine  _et al._ , 2016; Jeun  _et al._ , 2016; Schmidt
_et al._ , 2017b). Beneddine  _et al._ (2016) concluded that mean stability
analysis was valid when the dominant singular value of the resolvent operator
was significantly greater than the others at a given frequency and that this
condition holds for flows where there is a dominant convective instability
mechanism and an eigenvalue which is nearly marginally stable. In such
circumstances, it was shown that the eigenmodes are proportional to the
resolvent response modes.

###  I.2 Motivation and scope of the study

A broad review of the approaches discussed above can be found in Taira  _et
al._ (2017). Given the lack of a comprehensive understanding of when spectral
and pseudospectral analyses lead to perturbations with characteristics that
are observed in real flows, and when the base or mean flow will lead to
success, our approach is to synthesize the characteristics of the LNS operator
with a view to answering specific questions.

One of the main objectives of this paper is to understand the conditions under
which eigenanalyses and resolvent analyses yield the same results, i.e. when
spectral and pseudospectral analyses reveal the same most amplified
perturbations. Furthermore, we seek to address when each method is likely to
lead to successful prediction of flow features, and whether the properties of
the base or mean flow (which are inputs to form the linear operators for both
eigen- and resolvent analyses) indicate when these analyses coincide. This
leads to deeper questions about the resolvent operator and the types of
amplification mechanisms it admits.

It is well known that the response of a system to harmonic input (forcing) can
occur due to normal and/or non-normal mechanisms. The former corresponds to
excitation in the vicinity of an eigenvalue (with exponential growth if the
eigenvalue is unstable) while the latter is due to nonmodal effects associated
with the sensitivity of the spectrum to perturbation. Subsequent to the work
of Jovanović and Bamieh (2005), Marquet  _et al._ (2009) and Brandt  _et al._
(2011) investigated the distribution of energy and relative phase between
velocity components in analyses about laminar, base flows in recirculation
bubbles and the flat plate boundary layer, respectively. These studies
distinguished between component-type non-normalities which distribute energy
in different velocity components and convective non-normalities, which
separate the spatial support of forcing upstream of the response. The roots of
these non-normalities are the mean shear and mean flow advection terms,
respectively, in the linearized NSE. These terms also result in the Orr
mechanism Orr (1907) which reorients upstream-leaning forcing modes with the
mean shear such that the response modes are leaning downstream Farrell (1987).
In this paper we extend the analysis to include turbulent mean flows and
quantify non-normality based on inner products between the resolvent forcing
and response modes. The special case of limited spatial overlap of the forcing
and response modes, which is known as a wavemaker (e.g. Pier and Huerre, 2001)
and arises when there is sufficiently strong reverse flow Huerre and Monkewitz
(1985), is also addressed.

Non-normality can be investigated through the lens of the pseudospectrum of
the LNS operator. Most previous studies have only investigated the
pseudospectra of parallel base flows (e.g. Trefethen  _et al._ , 1993; Reddy
_et al._ , 1993; Trefethen  _et al._ , 1999; Schmid and Henningson, 2001;
Schmid, 2007; Schmid and Brandt, 2014). In this study, pseudospectral analysis
is extended to mean flows with or without spatial development. The
contribution from normality of the operator, or resonance, can be quantified
from the spectrum alone by computing the inverse distance between the
imaginary axis and the nearest eigenvalue. The contribution from non-normality
of the operator can be quantified by bounds of the pseudospectrum which may be
dramatically different from the spectrum if the operator is non-normal and
thus sensitive to perturbations. For complex flows where a variety of
frequencies and wavenumbers are active, both normal and non-normal mechanisms
may be simultaneously present and dynamically important.

Lastly, we consider the conditions under which the resolvent operator can be
considered to be low-rank. Here we use the term low-rank to describe an
operator whose behavior can be well approximated by a low-rank operator as
opposed to the meaning that the operator is of low numerical rank. McKeon and
Sharma (2010) identified low-rank characteristics of the resolvent for
turbulent wall flows via the rapid decrease of singular value with increasing
order of the singular value decomposition. They attributed this to the
dominance of specific physical mechanisms that are efficiently captured
through representation in terms of spatio-temporal modes. Subsequent work
(e.g. Gómez  _et al._ , 2016; Beneddine  _et al._ , 2016) has exploited the
low-rank approximation of the resolvent to construct reduced-order models of
complex flows. We show that the question of rank of the resolvent at a
particular frequency can be investigated by analyzing the pseudospectrum as
well as the ratio of the largest to smallest terms in the resolvent operator
itself.

###  I.3 Choice of flows and outline of the paper

We consider low Reynolds number cylinder flow and turbulent channel flow to
demonstrate the contributions of this paper. Flow around a circular cylinder
is an example of an oscillator flow (Huerre and Rossi, 1998) which has
intrinsic dynamics that are insensitive to background noise. The flow exhibits
a region of absolute instability Huerre and Monkewitz (1985); Monkewitz (1988)
which can be approximated using the wavemaker (e.g. Giannetti and Luchini,
2007). Using resolvent analysis and plotting bounds of the pseudospectrum for
a given perturbation magnitude, it is evident that resonance accounts for the
bulk of the amplification. The impact of mean flow advection, however, which
separates the spatial support of the forcing and response modes, introduces
convective non-normality and leads to a resolvent norm which is appreciably
larger than the contribution from resonance alone. In the case of turbulent
channel flow, which is a noise-amplifier Huerre and Rossi (1998), the choice
of temporal frequency and spatial wavenumbers has an impact on the influence
of non-normality. We analyze several structures which are either highly
amplified or representative of known turbulent structures in order to
establish the role of normal and non-normal mechanisms as well as generalize
when low-rank behavior can be expected.

The rest of the paper is organized as follows. In Section II, the governing
equations are derived for the linear operators which form the basis of
stability and resolvent analysis. In Section III, the resolvent operator is
rewritten to differentiate normal and non-normal amplification mechanisms. A
model LNS operator and its resolvent without spatial dependence is used to
isolate the effects of individual terms in the NSE on the forcing and response
modes. The effect of spatial dependence is discussed particularly under the
influence of mean shear and reverse flow. Finally, the resolvent operator is
rewritten using its dyadic representation to formally relate stability and
resolvent modes and show when the rank-1 approximation is appropriate for
normal mechanisms. Application of the findings to circular cylinder flow is
considered in Section IV; the numerical methods are outlined and resolvent
analysis is applied to both base and mean flows. The results are compared to
those from a stability analysis. Section V studies the observations from
earlier sections in the context of turbulent channel flow. The influence of
the wall-normal height and spatial and temporal wavenumbers are shown to play
a major factor in the type of amplification mechanisms which dominate as well
as whether or not the resolvent is low-rank. Conclusions are presented in
Section VI along with the implications for reduced-order modeling and control.

##  II Governing Equations

The relevant operators for the analyses that follow are derived from the
incompressible NSE which are non-dimensionalized by the characteristic length
and velocity scales, $L$ and $U$:

$\displaystyle\partial_{t}\bm{u}+\bm{u}\cdot\nabla\bm{u}$
$\displaystyle=-\nabla p+Re^{-1}\nabla^{2}\bm{u}$ (1a)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0.$ (1b)

The states $\bm{u}(\bm{x},t)$ and $p(\bm{x},t)$ are the spatially- and
temporally-varying velocity and pressure fields, respectively (explicit
statement of the spatial and temporal dependences will be dropped hereon for
conciseness), while $Re$ is the Reynolds number based on $L$ and $U$. After
Reynolds-decomposing the states into a stationary temporal mean (denoted by an
overline) and a fluctuating component (denoted by a prime), one obtains the
mean flow equations:

$\displaystyle\overline{\bm{u}}\cdot\nabla\overline{\bm{u}}+\nabla\overline{p}-Re^{-1}\nabla^{2}\overline{\bm{u}}$
$\displaystyle=-\overline{\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}}$ (2a)
$\displaystyle\nabla\cdot\overline{\bm{u}}$ $\displaystyle=0.$ (2b)

For an exact solution of the NSE, i.e. a true base flow, the divergence of the
Reynolds stress tensor, or the right-hand term in Equation 2a, is zero.
Henceforth we identify the temporal mean in such a case by
$\overline{\bm{u}}=\bm{U}_{0}$. If the term is nonzero, however,
$\overline{\bm{u}}$ does not constitute an exact solution and the action of
Reynolds stresses must be taken into account to satisfy the NSE. We will refer
to $\overline{\bm{u}}$ in this case as a (turbulent) mean flow.

Subtracting the mean momentum equations (Equation 2) from the NSE (Equation 1)
yields the following for the fluctuating quantities:

$\displaystyle\partial_{t}\bm{u}^{\prime}+\overline{\bm{u}}\cdot\nabla\bm{u}^{\prime}+\bm{u}^{\prime}\cdot\nabla\overline{\bm{u}}+\nabla
p^{\prime}-Re^{-1}\nabla^{2}\bm{u}^{\prime}=-\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}+\overline{\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}}=\bm{f}^{\prime}$
(3a) $\displaystyle\nabla\cdot\bm{u}^{\prime}=0.$ (3b)

Equation 3a has been written such that all linear terms appear on the left-
hand side. The nonlinear terms on the right-hand side can be lumped together
as a forcing $\bm{f}^{\prime}$ without loss of generality. The literature
devoted to analyzing these linearized equations is broad and covers a range of
flows. We recap here the analyses required for the subsequent development and
refer the reader to recent reviews and contributions (Chomaz, 2005; Schmid,
2007; Sipp  _et al._ , 2010; Schmid and Brandt, 2014; McKeon, 2017) for
further information.

Our development considers global and spatially periodic modes in the context
of bluff body and wall-bounded turbulent flows, respectively. The conclusions,
nevertheless, are applicable when using the parabolized stability equations
(PSE), (e.g. Beneddine  _et al._ , 2016), although PSE is appropriate only for
the treatment of convective instability and not absolute instability (Herbert,
1997). They are also applicable to both base and mean flows although there are
important differences between the results for these profiles which will be
drawn out. The temporal mean velocity profiles are known throughout the domain
from numerical simulation for the cylinder flow while they are obtained via an
eddy viscosity model (e.g. Reynolds and Tiederman, 1967) for the turbulent
channel flow. We begin by considering the general case in which there exists
invariance only in time. In other words, the analysis is performed in the
frequency domain, such that mode shapes may be functions of all three spatial
dimensions. Decomposition into spatial wavenumbers is applied in the case of
wall-bounded turbulent flows with streamwise and spanwise periodicity in
Section V.

###  II.1 Eigenmode decomposition

The classical (temporal) linear stability analysis, detailed, for example, in
Schmid and Henningson (2001), proceeds under the assumption of small
perturbations to the steady state. It has been performed relative to both base
and mean flows in order to determine frequencies subject to exponential growth
(instability). See e.g. Taira  _et al._ (2017) for further details.
Substituting perturbations of the form

$\bm{u}(\bm{x},t)=\overline{\bm{u}}(\bm{x})+\epsilon\tilde{\bm{u}}(\bm{x})\exp(\lambda
t),$ (4)

where $\epsilon\ll 1$, into Equation 1 yields at $\mathcal{O}(\epsilon)$:

$\displaystyle\lambda\tilde{\bm{u}}$
$\displaystyle=-\overline{\bm{u}}\cdot\nabla\tilde{\bm{u}}-\tilde{\bm{u}}\cdot\nabla\overline{\bm{u}}-\nabla\tilde{p}+Re^{-1}\nabla^{2}\tilde{\bm{u}}$
(5a) $\displaystyle\nabla\cdot\tilde{\bm{u}}$ $\displaystyle=0,$ (5b)

where a tilde will be used to denote that a stability analysis has been
performed. Written explicitly for base flows in operator form, one obtains

$\lambda\bm{B}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\
\tilde{p}\end{array}\right)=\bm{A}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\
\tilde{p}\end{array}\right),$ (6)

where $\bm{A}$ is the LNS operator with respect to the base flow,

$\bm{A}=\left(\begin{array}[]{cc}-\bm{U}_{0}\cdot\nabla()-()\cdot\nabla\bm{U}_{0}+Re^{-1}\nabla^{2}()&-\nabla()\\\
\nabla\cdot()&0\end{array}\right),$ (7)

and

$\bm{B}=\left(\begin{array}[]{cc}1&0\\\ 0&0\end{array}\right).$ (8)

Similarly, we use $\bm{L}$ to denote the LNS operator with respect to the mean
flow, such that for mean flow stability analysis

$\lambda\bm{B}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\
\tilde{p}\end{array}\right)=\bm{L}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\
\tilde{p}\end{array}\right),$ (9)

with

$\bm{L}=\left(\begin{array}[]{cc}-\overline{\bm{u}}\cdot\nabla()-()\cdot\nabla\overline{\bm{u}}+Re^{-1}\nabla^{2}()&-\nabla()\\\
\nabla\cdot()&0\end{array}\right).$ (10)

The resulting eigenvalue problems for base and mean flows are summarized in
the upper row of Table 1. Stability analysis is based on the spectrum of the
LNS operator and $\lambda$ is an indicator of the linear stability of a given
profile. The eigenvectors may be used as a basis for modal decomposition.

| Base Flow | Mean Flow
---|---|---
Stability Analysis | $\lambda\bm{B}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\ \tilde{p}\end{array}\right)=\bm{A}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\ \tilde{p}\end{array}\right)$ | $\lambda\bm{B}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\ \tilde{p}\end{array}\right)=\bm{L}\left(\begin{array}[]{c}\tilde{\bm{u}}\\\ \tilde{p}\end{array}\right)$
Resolvent Analysis |  $i\omega\bm{B}\left(\begin{array}[]{c}\hat{\bm{u}}\\\ \hat{p}\end{array}\right)=\bm{A}\left(\begin{array}[]{c}\hat{\bm{u}}\\\ \hat{p}\end{array}\right)+\bm{C}\hat{\bm{f}}$ |  $i\omega\bm{B}\left(\begin{array}[]{c}\hat{\bm{u}}\\\ \hat{p}\end{array}\right)=\bm{L}\left(\begin{array}[]{c}\hat{\bm{u}}\\\ \hat{p}\end{array}\right)+\bm{C}\hat{\bm{f}}$
Table 1: Operator form of the equations for stability and resolvent analyses.
Variables with a tilde correspond to stability analysis while a caret
indicates resolvent analysis.

The adjoint NSE have been derived by, e.g. Luchini and Bottaro (2014), and the
linearized operators can be written for both base and mean flows, $\bm{A}^{*}$
and $\bm{L}^{*}$, respectively, for the adjoint variables $\tilde{\bm{v}}$ and
$\tilde{q}$,

$\displaystyle\bm{A}^{*}=\left(\begin{array}[]{cc}\bm{U}_{0}\cdot\nabla()-()\cdot(\nabla\bm{U}_{0})^{*}+Re^{-1}\nabla^{2}()&\nabla()\\\
\nabla\cdot()&0\end{array}\right),$ (13)
$\displaystyle\bm{L}^{*}=\left(\begin{array}[]{cc}\overline{\bm{u}}\cdot\nabla()-()\cdot(\nabla\overline{\bm{u}})^{*}+Re^{-1}\nabla^{2}()&\nabla()\\\
\nabla\cdot()&0\end{array}\right).$ (16)

The operators satisfy

$\displaystyle\left<\tilde{\bm{u}},\bm{A}\tilde{\bm{v}}\right>=\left<\bm{A}^{*}\tilde{\bm{u}},\tilde{\bm{v}}\right>,$
(17)
$\displaystyle\left<\tilde{\bm{u}},\bm{L}\tilde{\bm{v}}\right>=\left<\bm{L}^{*}\tilde{\bm{u}},\tilde{\bm{v}}\right>,$
(18)

where $<,>$ is the scalar product associated with the energy in the whole
domain.

For a general operator, $\bm{T}$, that is normal, i.e.
$\bm{T}\bm{T}^{*}=\bm{T}^{*}\bm{T}$, the eigenvectors of $\bm{T}$
corresponding to distinct eigenvalues are orthogonal although the eigenvalues
may be complex. Self-adjoint operators $(\bm{T}=\bm{T}^{*})$, on the other
hand, have orthogonal eigenvectors and real eigenvalues. In general, the LNS
operators $\bm{A}$ and $\bm{L}$ are neither self-adjoint nor normal. These
phenomena account for the differences between Equations 7, 10 and Equations
13, 16 and their influence is discussed in Section III.

###  II.2 Resolvent analysis

For the more general case when the perturbation cannot be considered to be
infinitesimal and the nonlinearity $\bm{f}^{\prime}$ is retained, Equation 3a
can be rewritten in terms of a transfer function between the forcing (input)
and response state (output) (e.g. Schmid and Henningson, 2001; Jovanović and
Bamieh, 2005; McKeon and Sharma, 2010). This transfer function is the (linear)
resolvent operator, which can be defined around either the base or mean flow
(see Table 1). Endogenous nonlinear terms in Equation 3b when $\epsilon$ is
not required to be small or exogenous forcing can be treated equally well via
$\bm{f}^{\prime}$, although the interpretation of the resulting system is
different.

For harmonic forcing and response at temporal frequency $\omega$,

$\bm{f}^{\prime}=\hat{\bm{f}}~{}\text{exp}(i\omega
t),\quad\bm{u}^{\prime}=\hat{\bm{u}}~{}\text{exp}(i\omega t)$ (19)

and a base flow profile, $\overline{\bm{u}}=\bm{U}_{0}$,

$\hat{\bm{u}}=\mathcal{H}(\omega)\hat{\bm{f}},$ (20)

where the caret denotes that the perturbation is associated with a resolvent
analysis. The resolvent operator is given by

$\mathcal{H}(\omega)=\bm{C}^{T}(i\omega\bm{B}-\bm{A})^{-1}\bm{C},$ (21)

where

$\bm{C}=\left(\begin{array}[]{c}1\\\ 0\end{array}\right).$ (22)

Similarly, for a mean flow,

$\mathcal{H}(\omega)=\bm{C}^{T}(i\omega\bm{B}-\bm{L})^{-1}\bm{C}.$ (23)

It should be noted that the sense of the imaginary and real parts of $\omega$
is reversed in the resolvent analysis relative to the definition customary to
the stability literature of Equation 4: here the real part of $\omega$ is the
frequency associated with a mode while the imaginary part is set to zero as
only neutral disturbances are considered. In stability analysis, the imaginary
part of $\lambda$ is the frequency of the disturbance and the real part is the
growth rate.

$\mathcal{H}(\omega)$ can be decomposed via a singular value decomposition
(SVD), e.g. McKeon and Sharma (2010):

$\mathcal{H}(\omega)=\bm{\Psi}(\omega)\bm{\Sigma}(\omega)\bm{\Phi}^{*}(\omega),$
(24)

where $\bm{\Psi}$ and $\bm{\Phi}$ are the left and right singular vectors
corresponding to the response and forcing modes, often called resolvent modes
(see McKeon and Sharma, 2010), respectively. Both sets of singular vectors are
guaranteed to be orthonormal bases and are ranked according to their gain, or
singular value, contained in the diagonal matrix $\bm{\Sigma}$. The resolvent
operator can thus be written as the sum of outer products of the left and
right singular vectors

$\mathcal{H}(\omega)=\sum_{j=1}^{\infty}\hat{\bm{\psi}}_{j}(\omega)\sigma_{j}(\omega)\hat{\bm{\phi}}^{*}_{j}(\omega).$
(25)

$\mathcal{H}(\omega)$ is (approximately) low rank if

$\sum_{j=1}^{p}\sigma^{2}_{j}\approx\sum_{j=1}^{\infty}\sigma^{2}_{j},$ (26)

where $\sigma_{p}\gg\sigma_{p+1}$ and $p$ is small Moarref  _et al._ (2013);
McKeon (2017) . If the leading singular value is significantly greater than
all others ($\sigma_{1}\gg\sigma_{2}$) then the rank-1 approximation can be
invoked and the resolvent is approximated by the outer product of the leading
optimal response and forcing modes:

$\mathcal{H}(\omega)\approx\sigma_{1}\hat{\bm{\psi}}_{1}\hat{\bm{\phi}}_{1}^{*}.$
(27)

The physical interpretation of the resolvent response modes is the response to
forcing that results in a neutrally stable response, i.e. with the real
component of frequency equal to zero. The singular value gives the input-
output gain, here associated with the energy norm.

###  II.3 Spectrum and pseudospectrum of the LNS operator

Analyzing the resolvent corresponds to considering the spectrum of the
perturbed LNS operator:

$\Lambda_{\epsilon}(\bm{A})=\\{z\in\mathbb{C}:z\in\Lambda(\bm{A}+\bm{E})~{}\text{where}~{}\|\bm{E}\|\leq\epsilon\\},$
(28)

where $\Lambda_{\epsilon}$ is the pseudospectrum of $\bm{A}$ under a
perturbation magnitude $\epsilon>0$ and $\|\cdot\|$ is the (operator) 2-norm
Trefethen  _et al._ (1993); Reddy  _et al._ (1993); Taira  _et al._ (2017). An
equivalent definition is given by

$\Lambda_{\epsilon}(\bm{A})=\left\\{z\in\mathbb{C}:\|(z\bm{I}-\bm{A})^{-1}\|\geq\epsilon^{-1}\right\\}~{}\cup~{}\Lambda(\bm{A}),$
(29)

where $\Lambda=\Lambda_{0}$ is the spectrum of $\bm{A}$. Throughout the paper,
we will use $\Lambda$ to represent the set of eigenvalues and $\bm{\Lambda}$
the diagonal matrix of eigenvalues. If $\bm{A}$ is normal,
$\Lambda_{\epsilon}$ can be interpreted as the set of points away from
$\Lambda$ by only less than or equal to $\epsilon$ on the complex plane (see
Taira  _et al._ , 2017). If $\bm{A}$ is non-normal, this distance may be
greater than $\epsilon$ signifying that an eigenvalue is sensitive to
perturbation of the LNS operator.

For a given $z$, the resolvent norm $\|(z\bm{I}-\bm{A})^{-1}\|$ is equal to
the largest value of $\epsilon^{-1}$ such that $z$ is contained within
$\Lambda_{\epsilon}$. The resolvent norm is (by definition) the maximum
singular value of the resolvent operator, and quantifies the system’s
sensitivity to temporal forcing. The neutrally stable response of the system
to harmonic forcing is characterized by the value of the resolvent norm along
the imaginary axis. For a stable, normal operator, the largest response occurs
at the frequency corresponding to the imaginary part of the least stable
eigenvalue $\lambda_{\text{ls}}$ and the resolvent norm is
$1/\text{Real}(\lambda_{\text{ls}})$ since $\text{Real}(\lambda_{\text{ls}})$
is the minimum distance between the eigenvalue and the imaginary axis (Chomaz,
2005). If the operator is marginally stable
$(\text{Real}(\lambda_{\text{ls}})\approx 0)$, then $\epsilon\approx 0$ and
the response is dominated by the corresponding eigenmode. For a stable, non-
normal operator, the frequency and gain of the largest response is less
predictable since it is necessary to find the smallest value of $\epsilon$ for
which the pseudospectrum crosses the imaginary axis.

Having reviewed the derivations and interpretations of the two analyses that
will be considered here, we now elaborate on the mathematical and physical
origins of the associated amplification mechanisms and componentwise transfer.
Simple example operators are used to illustrate the concepts, in the spirit of
Gebhardt and Grossmann (1994), before we consider the complete LNS operator.

##  III Mechanisms for amplification

###  III.1 Resonance and pseudoresonance

The origin of the amplification mechanisms characterized by the resolvent norm
can be identified by expanding the resolvent through an eigenvalue
decomposition of the LNS operator,

$\mathcal{H}(\omega)=\bm{C}^{T}(i\omega\bm{B}-\bm{V}\bm{\Lambda
B}\bm{V}^{-1})^{-1}\bm{C}.$ (30)

Here $\bm{V}$ represents the matrix of eigenvectors of the LNS operator for
either the base or mean flow profile and $\bm{\Lambda}$ the diagonal matrix of
eigenvalues. These can be used to find an upper and lower bound for the
resolvent norm (see Schmid and Henningson, 2001):

$\|i\omega\bm{I}-\bm{\Lambda}\|^{-1}\leq\|\mathcal{H}(\omega)\|\leq\underbrace{\|\bm{V}\|\|\bm{V}^{-1}\|}_{\text{pseudoresonance}}\underbrace{\|i\omega\bm{I}-\bm{\Lambda}\|^{-1}}_{\text{resonance}}.$
(31)

Considering first the far-righthand term in Equation 31, it is clear that
forcing in the vicinity of an eigenvalue, i.e. $\omega=\lambda$, is likely to
lead to amplification due to resonance. This is a normal phenomenon,
predictable from an eigenanalysis. Large amplification also arises in the
event of pseudoresonance when the condition number
$\kappa=\|\bm{V}\|\|\bm{V}^{-1}\|$ is large due to non-orthogonality of the
eigenvectors, a consequence of the non-normal nature of $\bm{L}$ and hence
$\mathcal{H}(\omega)$. In the formulation of Equation 31, the resolvent
therefore contains both the amplification mechanisms associated with the
eigenvalue spectrum accessed via eigenanalysis (normal mode linear stability
analysis) and the non-normal amplification that is possible when the
eigenvectors are not orthogonal to each other. In what follows, we will
identify “optimal pseudoresonance” as the condition under which the
amplification in Equation 31 is associated with pseudoresonance, i.e.
$\kappa\to\infty$, to be compared with the resonant case in which $\kappa=1$.

(a)

(b)

(c)

(d)

Figure 1: Comparison of the pseudospectra and resolvent norm for the operators
given in Equations 34 and 37, which have the same eigenvalues. (a, b) normal
operator $\bm{S}$, (c, d) non-normal operator $\bm{P}$. The eigenvalues, i.e.
the eigenspectrum, are marked by red crosses and color contours outline the
bounds of the perturbed spectrum for constant perturbation magnitudes in (a,
c). The dashed contours in (a) reflect that the pseudospectra are circles
centered on the eigenvalues, which is not true for the non-normal operator in
(c). The resolvent norm in each case (b, d) reflects the value of these
contours along the imaginary axis. Red, dashed horizontal lines indicate the
resonant frequencies of the operator, i.e. the frequencies corresponding to
the eigenvalues, while the blue, solid horizontal line represents the most
highly amplified frequency in the non-normal case.

To demonstrate the amplification characteristics of the resolvent for normal
and non-normal operators, we sketch in Figure 1 the pseudospectra for simple
example operators,

$\displaystyle\bm{S}=\left(\begin{array}[]{cc}-1.5+1.1i&0\\\
0&-1.9-2.2i\end{array}\right),$ (34)
$\displaystyle\bm{P}=\left(\begin{array}[]{cc}-1.5+1.1i&5\\\
0&-1.9-2.2i\end{array}\right),$ (37)

where $\bm{S}$ is a normal operator containing only the eigenvalues of the
non-normal operator $\bm{P}$. Level curves of $\epsilon$ for operators
$\bm{S}$ and $\bm{P}$ satisfy

$\Lambda_{\epsilon}(\bm{S})=\\{z\in\mathbb{C}:\|(z\bm{I}-\bm{S})^{-1}\|\geq\epsilon^{-1}\\},$
(38)

and

$\Lambda_{\epsilon}(\bm{P})=\\{z\in\mathbb{C}:\|(z\bm{I}-\bm{P})^{-1}\|\geq\epsilon^{-1}\\},$
(39)

respectively.

For a normal operator such as $\bm{S}$, $\kappa=1$ and the level curves of
$\epsilon$ are proportional to the distance from the closest eigenvalue. The
resolvent norm for a particular $\omega$ is inversely proportional to the
distance from $i\omega$ to the nearest eigenvalue. The spectrum and
pseudospectra of $\bm{S}$ are shown in Figure 1(a); there are two stable
eigenvalues denoted by red crosses, and the pseudospectra consist of circular
contours centered on the two eigenvalues. Since both eigenvalues are
significantly damped, amplification due to resonance is not possible and the
magnitude of the resolvent norm is less than one (Figure 1(b)). Moreover, the
eigenvalue and singular value decompositions of $\bm{S}$ yield parallel basis
functions, and the singular values are simply the magnitude of the
eigenvalues.

Operator $\bm{P}$, however, is non-normal due to the non-zero off-diagonal
term and, with reference to Equation 31, $\kappa>1$. The shifts of the
eigenvalues of the perturbed operator are not proportional to $\epsilon$, as
indicated by the pseudospectrum isocontours in Figure 1(c), and the resolvent
norm of Figure 1(d) is appreciably larger than that of the normal case in
Figure 1(a), with values exceeding one. Furthermore, the maximum value of the
resolvent norm occurs at a non-resonant frequency, $\omega=0.66$.
Amplification is possible for a linearly stable operator due to
pseudoresonance under forcing at any frequency for which
$\|\mathcal{H}(\omega)\|>1$. It is important to add that even when the primary
contribution to amplification is a normal mechanism, non-normality can still
contribute to exacerbate the response. For example, the amplification at the
frequencies of the eigenvalues in Figure 1 is higher for operator $\bm{P}$
than it is for the purely normal operator $\bm{S}$. That being the case, the
right-hand side of Equation 31 may be large due to one or both terms in the
product.

###  III.2 Eigenvectors, singular functions, and non-normality of the operator

We specialize now to operators with similar behavior to the LNS system and
consider the associated features of eigenvector and singular value
decompositions. We choose a model operator with characteristics similar to the
LNS, in the vein of that explored by Gebhardt and Grossmann (1994) but without
specifying numerical values. The LNS for the two-dimensional velocity field
associated with a one-dimensional base or mean flow variation takes a similar
form to $\bm{M}$, where

$\bm{M}=\left(\begin{array}[]{cc}m_{1}&d\\\ 0&m_{2}\end{array}\right).$ (40)

By selecting a one-dimensional operator, we have elected to neglect spatial
(streamwise and spanwise) dependence of the base flow and therefore the modes
themselves here. Nevertheless, the impact of the various types of term in the
LNS operator on the resolvent modes can be modeled. Here
$\text{Real}(m_{j})<0$ is analogous to the stabilizing role of viscosity
through the $Re^{-1}\nabla^{2}()$ term and $\text{Imag}(m_{j})$ represents
mean flow advection through the $-\overline{\bm{u}}\cdot\nabla()$ term. $d$ is
real and is analogous to mean shear $()\cdot\nabla\overline{\bm{u}}$, which
here is equal to the gradient in the 2-direction of mean flow in the
1-direction. Thus in this simple 2-by-2 example, $d$ models the lift-up
mechanism (Landahl, 1980) by coupling forcing in the $n_{2}$-direction (second
component of the vector) with a response in the $n_{1}$-direction (first
component of the vector).

The resolvent of $\bm{M}$ is:

$\mathcal{H}(\omega)=\left(\begin{array}[]{cc}-1/(m_{1}-i\omega)&d/[(m_{1}-i\omega)(m_{2}-i\omega)]\\\
0&-1/(m_{2}-i\omega)\end{array}\right).$ (41)

In order to isolate the effect of non-normality introduced via the off-
diagonal term in Equation 42, the eigenvalues are assumed to be real and we
limit the immediate development to stationary disturbances ($\omega=0$) to
eliminate the remaining imaginary terms, such that

$\mathcal{H}(\omega=0)=\left(\begin{array}[]{cc}-1/(m_{1})&d/[(m_{1}m_{2})]\\\
0&-1/(m_{2})\end{array}\right).$ (42)

For the limiting case of $d=0$, i.e. zero mean shear, $\bm{M}$ and its
resolvent are self-adjoint and therefore normal. If the least stable
eigenvalue has real part close to zero which would occur, say, if $m_{1}\to
0$, the singular value decomposition can be simplified to

$\lim_{m_{1}\to
0}\text{SVD}(\mathcal{H}(\omega=0))=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)\left(\begin{array}[]{cc}\sigma_{1}&0\\\
0&\sigma_{2}\end{array}\right)\left(\begin{array}[]{ccc}1&0\\\
0&1\end{array}\right),$ (43)

where $\sigma_{1}/\sigma_{2}\to\infty$, i.e. the resolvent is low-rank. Thus
the response can be well predicted from the leading singular vectors
$\hat{\bm{\psi}}_{1}$ and $\hat{\bm{\phi}}_{1}$. For a normal operator, these
are identical to each other,
$\hat{\bm{\psi}}_{1}=\hat{\bm{\phi}}_{1}=[1~{}0]^{T}$, and identical to the
corresponding eigenvectors. The inner product
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|$ quantifies the componentwise
correspondence between the forcing and response modes which, in this limit, is
equal to unity. A schematic of the variation in the $n_{1}$ and $n_{2}$
directions of the forcing and response mode, i.e. the singular vectors, in
this case is shown in Figure 2(a-b).

The optimal pseudoresonance case occurs in the limit $d\to\infty$, for which:

$\lim_{d\to\infty}\text{SVD}(\mathcal{H}(\omega=0))=\left(\begin{array}[]{cc}1-\gamma&-\delta\\\
\delta&1-\gamma\end{array}\right)\left(\begin{array}[]{cc}\sigma_{1}&0\\\
0&\sigma_{2}\end{array}\right)\left(\begin{array}[]{cc}\delta&\gamma-1\\\
1-\gamma&\delta\end{array}\right),$ (44)

where $\sigma_{1}/\sigma_{2}\to\infty$ and $\delta,~{}\gamma\to 0$. The
constants $\delta$ and $\gamma$ are real and positive. The resolvent operator
is still low-rank in this limit, but $\hat{\bm{\psi}}_{1}$ and
$\hat{\bm{\phi}}_{1}$ are now orthogonal to each other and thus the inner
product $|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|\rightarrow 0$. The
perturbation energy in the optimal forcing mode is concentrated in the second
component of the vector while the perturbation energy in the optimal response
mode is concentrated in the first component, as sketched in Figure 2(c-d).

The analogous eigenvalue decomposition of (non-normal) $\bm{M}$ is

$\lim_{d\to\infty}\text{EIG}(\bm{M})=\left(\begin{array}[]{cc}1&1\\\
0&\alpha\end{array}\right)\left(\begin{array}[]{cc}m_{1}&0\\\
0&m_{2}\end{array}\right)\left(\begin{array}[]{cc}1&-1/\alpha\\\
0&1/\alpha\end{array}\right),$ (45)

where $\alpha\to 0$ is a positive, real constant. Unlike the resolvent
response modes which are orthogonal to one another, the eigenvectors are non-
orthogonal, such that $\kappa\to\infty$. Thus, in this case the stability and
resolvent modes are different. Since the eigenvectors are nearly parallel,
they both project equally well onto the optimal resolvent response mode. The
same can be said for the projection of the adjoint modes, which are also
nearly parallel, onto the optimal resolvent forcing mode. The optimal response
and forcing modes, therefore, are linear combinations of multiple stability
and adjoint modes, respectively. Furthermore, it is not clear from this
decomposition that the resolvent operator is low-rank.

As outlined in Figure 2, the model operator given in Equation 40 can be used
to elucidate that the presence of mean shear, which introduces non-normality
into the LNS operator, is to concentrate energy in different velocity
components of the resolvent forcing and response modes. The model LNS operator
does not need a base or mean flow with any special characteristics or spatial
dependence other than an off-diagonal term to reveal this. Nevertheless, the
same mechanism will be present in the case of more complex variation of shear,
as is common in real flows. In the normal limit where there is no mean shear,
the forcing and response are in the same velocity component. In the
pseudoresonant limit, the forcing acts in the $n_{2}$ direction while the
response is in the $n_{1}$ direction.

###  III.3 Eigenvectors, singular functions, and self-adjointness of the
operator

In this section, we isolate the effect (or lack thereof) of self-adjointness
on the characteristics of the resolvent modes. To this end, the mean shear
term is set to zero, i.e. $d=0$ in Equation 40, guaranteeing that the LNS and
resolvent operators are normal:

$\displaystyle\bm{M}=\left(\begin{array}[]{cc}m_{1}&0\\\
0&m_{2}\end{array}\right),$ (48)
$\displaystyle\mathcal{H}(\omega)=\left(\begin{array}[]{cc}1/(m_{1}-i\omega)&0\\\
0&1/(m_{2}-i\omega)\end{array}\right).$ (51)

Again we will limit the discussion to stationary disturbances $(\omega=0)$ as
the presence of $i\omega$ in the resolvent operator will always guarantee that
it is not self-adjoint.

If the eigenvalues are real, $\bm{M}$ and $\mathcal{H}(\omega=0)$ are self-
adjoint and the SVD of the resolvent is

$\lim_{\text{Imag}(m_{j})\to
0}\text{SVD}(\mathcal{H}(\omega=0))=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)\left(\begin{array}[]{cc}m_{1}&0\\\
0&m_{2}\end{array}\right)\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right).$ (52)

Equation 53 is similar to the normal limit case of Equation 43 and illustrates
that the resolvent forcing and response modes are identical.

The simplest way to demonstrate what happens when the resolvent is not self-
adjoint is to assume that the eigenvalues are purely imaginary. In this case,
the SVD of the resolvent is

$\lim_{\text{Real}(m_{j})\to
0}\text{SVD}(\mathcal{H}(\omega=0))=\left(\begin{array}[]{cc}i&0\\\
0&i\end{array}\right)\left(\begin{array}[]{cc}\|m_{1}\|&0\\\
0&\|m_{2}\|\end{array}\right)\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right).$ (53)

Since the singular values of the resolvent are required to be real yet the
eigenvalues are imaginary, the resolvent response modes must be $90^{\circ}$
out of phase from the resolvent forcing modes as is portrayed in Figure 3.

###  III.4 Mean flow advection and the Orr mechanism

In the case of an operator that is not self-adjoint and has non-zero mean
shear, i.e. $d\neq 0$, the spatio-temporal manifestation of the Orr mechanism
is observed. This mechanism results in disturbance amplification by the
reorientation of an input tilted against the mean shear into a response which
is aligned with the mean shear, with the maximum amplification in the
transient case occurring when the disturbance is vertical.

In the continuously forced formulation of the present problem, the Orr
mechanism is attributable to the space dependence of $\overline{\bm{u}}$ in
the mean flow advection term. Assuming a parallel flow
$(\overline{\bm{u}}=\overline{u}(y))$ with non-zero mean shear
$(\partial\overline{u}/\partial y\neq 0)$ and traveling at the local mean
velocity for a given $y_{c}$, the flow above $y_{c}$ moves towards the
downstream direction while the flow below moves upstream in a relative sense.
There is, therefore, a phase difference of $\pi$ across $y_{c}$ which
manifests itself in the response modes (e.g. McKeon, 2017). The decrease in
phase results in the response modes leaning downstream, aligned with the mean
shear. When considering the adjoint LNS operator $\bm{L}^{*}$, the mean flow
advection acts in the upstream direction due to the sign change of
$\bm{\overline{u}}\cdot\nabla()$. The corresponding phase jump across $y_{c}$
is now in the opposite sense and this results in the forcing modes leaning
upstream. Artificially removing mean flow advection from the resolvent
operator suppresses the Orr mechanism. Figure 4 is a cartoon illustrating how
the Orr mechanism results in the tilting of the resolvent forcing and response
modes.

###  III.5 Resolvent (approximate) wavemaker

For a spatially-varying base or mean flow, further statements can be made
concerning the effect of the non-self-adjoint nature of the resolvent operator
on the resolvent mode shapes, namely a difference in the spatial support of
forcing and response modes, with the latter being downstream of the former.
For a spatially-developing base or mean flow, normal mechanisms can be
categorized in terms of either convective or absolute instability depending on
the characteristics of the profile. The non-self-adjoint nature of $\bm{L}$
changes the influence of the convective terms in the adjoint operator, per
Equations 13 and 16, or, for the model operator of Equation 40,

$\bm{M}^{*}=\left(\begin{array}[]{cc}m_{1}^{*}&0\\\
d&m_{2}^{*}\end{array}\right).$ (54)

The direction of mean flow advection is reversed since the adjoint of the
derivative operator introduces a negative sign implying that adjoint
perturbations are transported upstream. Direct or forward perturbations, on
the other hand, are transported downstream. For absolutely unstable flows, for
which perturbations grow both upstream and downstream of the source, the
advection term may no longer separate the spatial support of the
forcing/adjoint and response/forward modes leading to regions of overlap at
resonant frequencies. This region, known as the wavemaker, is traditionally
computed from the eigenmodes, and is associated with non-zero values of
$\mathcal{W}$, where

$\mathcal{W}(\bm{x}_{0})=\|\tilde{\bm{u}}(\bm{x}_{0})\|\|\tilde{\bm{v}}(\bm{x}_{0})\|,$
(55)

and $\bm{x}_{0}$ denotes a position in space (see derivation of Giannetti and
Luchini, 2007).

In cases where amplification is due to normal mechanisms, the resolvent modes
can be used to find the wavemaker as long as they are normalized
appropriately. The wavemaker approximates regions of the flow which are
absolutely unstable or self-sustaining since perturbations are prevented from
convecting due to reverse flow Juniper (2012). Huerre and Monkewitz (1985)
have shown that when a mean profile of hyperbolic tangent form exhibits
greater than 13.6% reverse flow with respect to the free stream, the flow is
absolutely unstable. The streamwise extent of absolute instability and the
wavemaker is finite since flow reversal is confined to a certain portion of
the flow. This information is encoded within the advection term
$\overline{\bm{u}}\cdot\nabla()$ through the sign of $\overline{\bm{u}}$. We
note that the overlap of the resolvent forcing and response modes was
identified as a qualitative proxy for sensitivity for base flows by Brandt
_et al._ (2011). Henceforth we will denote the wavemaker of Equation 55 as the
“true” wavemaker and the approximation using resolvent modes as the
“resolvent” wavemaker. If the flow is convectively unstable, there is no
region of reverse flow and so $\overline{\bm{u}}$ is always positive. In this
case, the optimal response or stability mode will be downstream of the optimal
forcing or adjoint mode. The inner product
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|$ decreases when the modes are
separated in space. Chomaz (2005) noted that this is due to the convective-
type non-normality introduced from advection of the base flow.

The cartoon in Figure 5 summarizes the relative positioning of the various
modes for a 2D normal operator. The top row is characteristic of a convective-
type instability where mean flow advection orients the optimal forcing
upstream of a downstream, detached response. The bottom row represents a case
where there is overlap between the forcing and response modes indicating a
resolvent wavemaker where the flow exhibits an absolute-type instability. The
implications of this approximate wavemaker region are discussed in the context
of cylinder flow in Section IV.

(a)

(b)

(c)

(d)

Figure 2: Cartoon of forcing and response modes corresponding to the operator
given in Equation 40, where the resolvent norm is large, for the limiting
values of $d$. Panels (a, b) are the resonant forcing and response,
respectively, for the case $d=0$ where the operator is self-adjoint with
$\kappa=1$ (and hence normal) while panels (c, d) are the optimal
pseudoresonance forcing and response, respectively, for $d\to\infty$ and hence
$\kappa\to\infty$. Positive/negative isocontours are denoted by solid/dotted
lines and blue/red colors indicate streamwise/transverse components, in the
$n_{1},n_{2}$ directions, respectively. Each mode is nonzero in one velocity
component only.

(a)

(b)

Figure 3: Cartoon illustrating the $\pi/2$ phase difference between (a) the
resolvent forcing mode and (b) the resolvent response mode in Equation 53.
Positive/negative isocontours are denoted by solid/dotted lines, respectively.

(a)

(b)

Figure 4: Cartoon illustrating the action of the Orr mechanism in resolvent
modes. (a) The resolvent forcing mode leans upstream against the mean shear
and (b) the resolvent response mode is rotated by the mean shear to lean
downstream. Positive/negative isocontours are denoted by solid/dotted lines,
respectively.

(a)

(b)

(c)

(d)

Figure 5: Cartoon of the leading forcing and response modes corresponding to
the operator given in Equation 40 when the base or mean flow has variation in
both directions. Panels (a, b) are the forcing and response, respectively,
convective-type instability where mean flow advection results in the forcing
(a) being upstream of the response (b). Panels (c, d) show the forcing and
response, respectively, for an absolute-type instability. The vertical dashed-
dotted lines denote the region of overlap in the streamwise direction of
forcing (c) and response (d), i.e. the resolvent wavemaker. Positive/negative
isocontours are denoted by solid/dotted lines, respectively.

###  III.6 Low-rank approximation of the resolvent through a dyad expansion

We present what is believed to be a novel (if logical) expansion of the above
analysis to demonstrate and formalize the conditions under which analysis of
the resolvent is likely to identify stability modes as the most amplified
disturbance. We also highlight when the rank-1 approximation is appropriate
for amplifications which are normal in character.

As is more customary for the eigenvalue problem (Luchini and Bottaro, 2014;
Schmid and Brandt, 2014) a dyad expansion of the resolvent $\bm{R}$ for a
generic, non-singular linear operator $\bm{Q}$ can be performed,

$\bm{R}=(z\bm{I}-\bm{Q})^{-1}=\sum_{j=1}^{n}\frac{1}{z-\lambda_{j}}\tilde{\bm{g}}_{j}\tilde{\bm{h}}_{j}^{*},$
(56)

where $\tilde{\bm{g}}_{j}$ and $\tilde{\bm{h}}_{j}$ are the $j$th left and
right eigenvectors of $\bm{Q}$, respectively. Since the objective of resolvent
analysis is often the identification of the most amplified neutral
disturbance, $i\omega$ is substituted for $z$ and the eigenvectors of the LNS
operator for $\tilde{\bm{g}}_{j}$ and $\tilde{\bm{h}}_{j}$ into Equation 56 to
give

$\mathcal{H}(\omega)=\sum^{n}_{j=1}\frac{1}{i\omega-\lambda_{j}}\tilde{\bm{u}}_{j}\tilde{\bm{v}}_{j}^{*}.$
(57)

Thus if the real part of an eigenvalue $\lambda_{p}$ is sufficiently close to
zero and the forcing frequency $\omega$ is identical to the imaginary part,
then its contribution to the series dominates over the contributions from all
other eigenvalues. The resolvent, furthermore, can be approximated by the
forward and adjoint eigenvectors corresponding to that frequency weighted by
the inverse distance between the eigenvalue and the imaginary axis:

$\mathcal{H}(\omega)\approx\frac{1}{i\omega-\lambda_{p}}\tilde{\bm{u}}_{p}\tilde{\bm{v}}_{p}^{*}.$
(58)

Equation 58 represents a rank-1 approximation of the resolvent operator using
eigenvectors. In the context of base flows, the resolvent is singular at the
critical Reynolds number since the real part of the least stable eigenvalue is
identically zero when it crosses the imaginary axis. In the case of mean
flows, which tend to be marginally stable (e.g. Reynolds and Tiederman, 1967;
Barkley, 2006; Turton  _et al._ , 2015), Equation 58 is applicable for
eigenvalues near the imaginary axis. It is important to note that an
eigenvalue does not have to be marginally stable, but it must be the dominant
contribution to the series in Equation 57. It is possible to obtain a rank-1
approximation of the resolvent even when an eigenvalue is highly damped as
will be seen in Section IV.

This rank-1 approximation fails when there is not sufficient separation of
eigenvalues at the frequency of interest (see Appendix A for an example). If
there are several eigenvalues in the vicinity of the imaginary axis at a
frequency $\omega$, then the resolvent operator can no longer be approximated
by just one outer product in Equation 58. Equating the two low-rank
approximations of the resolvent operator in terms of eigenvectors (Equation
58) and resolvent modes (Equation 27) implies the following:

$\sigma_{1}\hat{\bm{\psi}}_{1}\hat{\bm{\phi}}_{1}^{*}\approx\frac{1}{i\omega-\lambda_{r}}\tilde{\bm{u}}_{r}\tilde{\bm{v}}^{*}_{r}\implies\hat{\bm{\psi}}_{1}\propto\
\tilde{\bm{u}}_{r},~{}~{}~{}\hat{\bm{\phi}}_{1}\propto\tilde{\bm{v}}_{r},$
(59)

since

$\hat{\bm{\psi}}_{1}\approx\frac{1}{\sigma_{1}(i\omega-\lambda_{r})}\tilde{\bm{u}}_{r}\tilde{\bm{v}}_{r}^{*}\hat{\bm{\phi}}_{1}=\beta\tilde{\bm{u}}_{r},$
(60)

where $\beta$ is a complex constant. The leading resolvent response and
forcing modes are proportional to the forward and adjoint eigenmodes,
respectively, and this holds for any base or mean flow as long as only one
eigenvalue leads to amplification. The similarity between the resolvent
forcing and adjoint stability modes draws out how the resolvent operator
contains sensitivity information, as described by, (e.g. Qadri and Schmid,
2017). The development is less amenable to non-normal mechanisms where the
proximity of an eigenvalue to the imaginary axis does not necessarily govern
the behavior of the resolvent (see, e.g. Figure 1). Some guidelines for when a
rank-1 approximation is expected to work for a general resolvent operator are
presented in Section V.

###  III.7 The relationship between spectral radius and spectral norm for
approximately low-rank operators

For a nonsingular linear operator $\bm{Q}$, we now seek to find an explicit
relationship between the spectrum of $\bm{Q}$ and the spectral norm of its
resolvent $\bm{R}(z)=(z\bm{I}-\bm{Q})^{-1}$. We are particularly interested in
cases where $\bm{R}(z)$ is approximately low-rank (i.e., a small number of
leading singular values are much larger than the others).

The spectral radius $\rho$ of an operator $\bm{Q}$ can be defined through the
eigendecomposition $\bm{V\Lambda V}^{-1}$

$\rho(\bm{Q})=\max_{\lambda_{j}\in\bm{\Lambda}}(|\lambda_{j}|).$ (61)

The spectral radius of the corresponding resolvent operator is

$\rho(\bm{R}(z))=\max_{\lambda_{j}\in\bm{\Lambda}}(|z-\lambda_{j}|^{-1})=\left[\min_{\lambda_{j}\in\bm{\Lambda}}(|z-\lambda_{j}|)\right]^{-1}.$
(62)

Note that, with this definition, Equation 31 may be expressed as

$\rho(\bm{R})\leq\sigma_{1}\leq\kappa\rho(\bm{R}).$

For non-normal operators with large condition numbers, the upper and lower
bounds span a large range, and thus do not give much insight into the size of
the resolvent norm, $\sigma_{1}$. To estimate the resolvent norm in terms of
the spectral radius, we will make use of the relationship Gelfand (1941)

$\rho(\bm{R})=\lim_{n\to\infty}\|\bm{R}^{n}\|^{1/n}.$ (63)

Suppose now that the largest singular value of resolvent operator (for a given
$z$) is much larger than the rest, such that

$\bm{R}(z)=\bm{\Psi}\bm{\Sigma}\bm{\Phi}^{*}\approx\bm{\Psi}\bm{\Sigma}_{1}\bm{\Phi}^{*},$

where $\bm{\Sigma}_{1}$ is $\bm{\Sigma}$ with all but the first singular value
set to zero.

Suppose in addition that this truncation is also accurate for powers of
$\bm{R}$, i.e., that we have

$\bm{R}^{n}\approx(\bm{\Psi}\bm{\Sigma}_{1}\bm{\Phi}^{*})^{n}.$ (64)

Defining the quantity

$r_{ij}=\frac{\bm{\phi}_{i}^{*}\bm{\psi}_{j}}{\bm{\phi}_{j}^{*}\bm{\psi}_{j}},$

we then have

$\displaystyle\left(\bm{\Phi}^{*}\bm{\Psi}\bm{\Sigma}_{1}\right)^{n}$
$\displaystyle=\sigma_{1}^{n}(\bm{\phi}_{1}^{*}\bm{\psi}_{1})^{n}\begin{pmatrix}r_{11}&0&\cdots&0\\\
r_{21}&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\end{pmatrix}.$

We may now estimate the norm of powers of the resolvent as

$\|\bm{R}^{n}\|=\|\bm{\Phi}^{*}\bm{R}^{n}\bm{\Phi}\|\\\
\approx\|\left(\bm{\Phi}^{*}\bm{\Psi}\bm{\Sigma}_{1}\right)^{n}\|\\\
=\sigma_{1}^{n}|\bm{\phi}_{1}^{*}\bm{\psi}_{1}|^{n-1}\|\bm{r}\|,$

where $\bm{r}=[r_{11}\ r_{21}\ \cdots]^{T}$, and we have used the fact that
$\bm{\Phi}$ is unitary. Consequently, assuming that Equation 64 holds,
Equation 63 results in the estimate

$\rho(\bm{R})\approx\sigma_{1}|\bm{\phi}_{1}^{*}\bm{\psi}_{1}|.$ (65)

In other words, we may estimate that the resolvent norm is larger than the
lower bound in Equation 31 by a factor of
$|\bm{\phi}_{1}^{*}\bm{\psi}_{1}|^{-1}$. This analysis relied on the rather
restrictive assumption that only the leading singular value was large. If
there is a pair of large singular values, as is often the case in channel
flows (owing to spatial symmetry across the mid-plane of the channel) then we
may generalize the argument as follows. Suppose that $\sigma_{1}$ and
$\sigma_{2}$ are of comparable size, and that all other singular values are
negligibly small. If we further assume that
$|\bm{\phi}_{1}^{*}\bm{\psi}_{2}|,\ |\bm{\phi}_{2}^{*}\bm{\psi}_{1}|\approx
0$, then we find that

$\displaystyle\left(\bm{\Phi}^{*}\bm{\Psi}\bm{\Sigma}\right)^{n}$
$\displaystyle\approx\begin{pmatrix}\sigma_{1}\bm{\phi}_{1}^{*}\bm{\psi}_{1}&\sigma_{2}\bm{\phi}_{1}^{*}\bm{\psi}_{2}&0&\cdots&0\\\
\sigma_{1}\bm{\phi}_{2}^{*}\bm{\psi}_{1}&\sigma_{2}\bm{\phi}_{2}^{*}\bm{\psi}_{2}&0&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\end{pmatrix}^{n}$
$\displaystyle\approx\sigma_{1}^{n}|\bm{\phi}_{1}^{*}\bm{\psi}_{1}|^{n-1}\begin{pmatrix}r_{11}&0&\cdots&0\\\
r_{21}&0&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\end{pmatrix}+\sigma_{2}^{n}|\bm{\phi}_{2}^{*}\bm{\psi}_{2}|^{n-1}\begin{pmatrix}0&r_{12}&0&\cdots&0\\\
0&r_{22}&0&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\end{pmatrix},$

which, following the same approach as before, gives

$\rho(\bm{R})\approx\max\\{\sigma_{1}|\bm{\phi}_{1}^{*}\bm{\psi}_{1}|,\sigma_{2}|\bm{\phi}_{2}^{*}\bm{\psi}_{2}|\\}.$
(66)

Thus the inverse of $|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|$ can be
interpreted as the contribution of non-normality to the resolvent norm. We can
also identify the product $\sigma_{1}|i\omega-\lambda|$ as a quantification of
non-normality since $|i\omega-\lambda|^{-1}$ represents the resonance
contribution to the resolvent norm. However, since highly amplified modes may
occur at non-resonant frequencies, the contribution from $|i\omega-\lambda|$
is typically overestimated as it is likely for a pseudoeigenvalue to reside
much closer to the imaginary axis than the nearest eigenvalue of the
unperturbed spectrum. As will be seen in Sections IV and V, these two
predictions tend to agree in cases where amplification can be attributed to a
single eigenvalue and mean stability analysis is valid.

Having examined the implications of the structure of the operator on
amplification and forcing and response modes, we now consider two real example
flows. Low Reynolds number cylinder flow is used to investigate the choice of
base or mean flow as the linear stability threshold is crossed (Section IV). A
canonical wall turbulence configuration is employed to identify the influence
of the various terms in the resolvent on the resulting SVD (Section V).

##  IV Application to cylinder flow

We apply a global resolvent analysis to the base and mean velocity profiles
for cylinder flows under the critical Reynolds number $Re_{c}\leq 47$
(Provansal  _et al._ , 1987; Sreenivasan  _et al._ , 1987; Noack and
Eckelmann, 1994), as well as mean flows of the 2D laminar vortex shedding
regime where $Re\leq 189$ (Barkley and Henderson, 1996). Cylinder flow is a
particularly suitable choice to investigate trends associated with the
wavemaker since it exhibits a region of absolute instability.

###  IV.1 Numerical methods

The relevant procedures for computing the two-dimensional base and mean flows,
$\bm{U}_{0}$ and $\overline{\bm{u}}$, are detailed here before applying the
analysis tools.

The NSE (Equation 1) are non-dimensionalized by the cylinder diameter $D$ and
inlet velocity $U_{\infty}$ which are both set to unity. For the base flow
calculation, a uniform inlet velocity condition is prescribed while no-slip
Dirichlet boundary conditions are applied to the cylinder surface, symmetric
conditions to the upper and lower boundaries, and advective conditions to the
outlet. The nonlinear equations for $\bm{U}_{0}$ are solved using a Newton
method on a finite-element mesh generated by FreeFem++ (see Hecht (2012)).
Taylor-Hood finite elements (P1b, P1b, P1 for $U_{0}$, $V_{0}$, and $P_{0}$
respectively) are used for spatial discretizations. The computational domain
$\Omega$ spans $-30\leq x/D\leq 60,-25\leq y/D\leq 25$ with the cylinder
centered at the origin and the mesh is made up of 104,214 triangles resulting
in 365,358 degrees of freedom for velocity and pressure.

A DNS of the cylinder flow is also performed to obtain the mean flow profile
using FreeFem++ with the same boundary conditions and mesh. A second-order
semi-implicit time discretization is employed with a non-dimensional time step
$\Delta t=0.02$. Beyond $Re_{c}$, the simulated flow settles into regular
vortex shedding at a fixed amplitude $A$ and temporal frequency $\omega_{s}$
where the subscript $s$ denotes shedding. The mean flow $\overline{\bm{u}}$ is
computed by time-averaging the DNS state vector over 25 complete shedding
cycles. The linear operators are formed in FreeFem++ and the only boundary
condition which differs with respect to the base flow calculation is at the
inlet where homogeneous boundary conditions are enforced so that the
perturbations vanish at infinity. The eigenvalues are computed using a shift-
and-invert strategy, the details of which are discussed in Nayar and Ortega
(1993). The generalized eigenvalue problem is then solved with the Implicitly
Restarted Arnoldi method using the ARPACK library developed by Lehoucq and
Sorensen (1996).

The singular values of the resolvent operator are computed in a manner
outlined by Sipp and Marquet (2013); a brief summary of the procedure is
presented here. The singular value problem is reformulated as the following
eigenvalue problem:

$\mathcal{H}(\omega)^{*}\mathcal{H}(\omega)\hat{\bm{\phi}_{i}}=\sigma_{i}^{2}\hat{\bm{\phi}_{i}},$
(67)

where $\hat{\bm{\phi}_{i}}$ is the $i$th right singular vector corresponding
to the singular value $\sigma_{i}$ of $\mathcal{H}(\omega)$. The largest
eigenvalues of the Hermitian operator
$\mathcal{H}(\omega)^{*}\mathcal{H}(\omega)$ are computed using the ARPACK
library and the parallel MUMPS solver developed by Amestoy  _et al._ (2001).
The response modes are then computed from Equation 20.

###  IV.2 Base flow velocity profile

A resolvent analysis is performed on the cylinder base flow for various
Reynolds numbers over a range of $\omega$. Contours of the pseudospectrum for
$Re=47$ are overlaid onto the spectrum, which is in agreement with Sipp and
Lebedev (2007) to within the sensitivity to the mesh geometry, in Figure 6(a).
The variation of the resolvent norm along the imaginary axis, i.e.
$\sigma_{1}$, is plotted alongside the second singular value of the resolvent,
$\sigma_{2}$, and the inverse of the distance between the nearest eigenvalue
and the imaginary axis in Figure 6(b). There is only one frequency
$\omega_{\text{max}}$ of interest, where the first singular value is several
orders of magnitude larger than all the others (only two are shown for
clarity). For the cylinder base flow at $Re=47$, $\omega_{\text{max}}=0.742$
corresponds to the true $\omega_{s}$ at the onset of vortex shedding. The
amplification is significantly lower for all other frequencies, including
harmonics, and $\omega_{\text{max}}$ is very close to the imaginary part of
the least stable eigenvalue. In short, resolvent analysis of the base flow
yields useful predictions for the frequency of the unsteady flow in the
vicinity of the bifurcation like its stability counterpart.

(a)

(b)

Figure 6: (a) Spectrum (red dots) and pseudospectrum (filled contours,
$\epsilon$ increasing as colors change from dark to light) of the LNS operator
for the cylinder base flow at $Re=47$. (b) The resolvent norm, $\sigma_{1}$
(solid line), i.e. the value of $\epsilon^{-1}$ along the imaginary axis,
second largest singular value $\sigma_{2}$ (dash-dotted line) and inverse
distance from the imaginary axis to the nearest eigenvalue (dotted red line).

The previous observations reinforce why stability analysis about the base flow
can predict $Re_{c}$. Only one structure at the globally most amplified
frequency is prone to significant amplification at subcritical Reynolds
numbers and it is the first to become unstable. This is characteristic of an
absolute instability mechanism in which frequency selection is not influenced
by background noise. The stability modes and resolvent modes are nearly
identical as seen in Figure 7.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7: Comparison of stability modes (left) with resolvent modes (right) at
the critical Reynolds number $Re_{c}=47$ and a temporal frequency of
$\omega=0.742$. Panels (a, b) are the streamwise component of the forward or
response mode, (c, d) are the transverse component of the forward or response
mode, (e, f) are the streamwise component of the adjoint or forcing mode, and
(g, h) are the transverse component of the adjoint or forcing mode. The
eigenmodes and resolvent mode shapes are essentially indistinguishable for
this flow.

The effect of lift-up is weak since the energy is fairly evenly distributed in
the $u-$ and $v-$components of both the forcing and response modes in Figure
7. Mean flow advection, on the other hand, plays a significant role in the
spatial support of the forcing and response modes which are located upstream
and downstream of the cylinder. Chomaz (2005) made an analogous observation
for the forward and adjoint eigenmodes and attributed this to convective non-
normality. Since the resolvent operator is low-rank, computing
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}$ is a good estimate of the
non-normal amplification experienced by the flow. We obtain a value of
$\|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}\|^{-1}=79.4$ which is in good
agreement with
$\sigma_{1}(\omega_{\text{max}})|i\omega_{\text{max}}-\lambda_{\text{ls}}|=79.3$
(see Table 2) in Figure 6(b). There is also a physical overlap of the forcing
and response modes in the region where the flow is absolutely unstable, which
can be attributed to the extent of reverse flow.

The least stable global mode and its adjoint counterpart are computed for
various Reynolds numbers near and below $Re_{c}$ to illustrate the cylinder
transition from convective to absolute instability. Figure 8 juxtaposes the
$v$-component of the adjoint mode, forward mode, and wavemaker. The forward
mode has unit magnitude while the adjoint has been normalized with respect to
the forward mode such that their inner product is unity. A wavemaker first
appears for $Re=25$, the Reynolds number at which Monkewitz (1988) determined
the cylinder wake is absolutely unstable. There is no wavemaker for lower
Reynolds numbers due to the downstream location of the forward eigenmode which
is a consequence of mean flow advection; the strength of the reverse flow is
not sufficient to produce an overlap region. For the lowest two Reynolds
numbers considered in Figure 8, the contour levels of the forward eigenmode
immediately behind the cylinder are three orders of magnitude smaller than
their higher Reynolds number counterparts. The downstream location where the
contour levels are significant does not appear within the plotted domain. As
the Reynolds number increases, the velocity deficit grows and the reverse flow
directly behind the cylinder strengthens. The forward eigenmode gradually
appears closer to the cylinder until there is a nontrivial overlap between it
and its adjoint counterpart.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

Figure 8: Contours of the transverse velocity for the leading adjoint modes
$\tilde{v}^{\dagger}$ (left) and forward modes $\tilde{v}$ (middle) of the
base flow. The Reynolds numbers ($Re=$15; 25; 35; 45; and 50) increase from
top to bottom. The wavemaker $\mathcal{W}$ (right) is computed using the
forward and adjoint modes. Contour levels are not identical for the adjoint
modes which are normalized based on the forward modes. Note that the
streamwise velocity component has not been plotted even though the wavemaker
depends on this quantity.

Beyond the critical Reynolds number, the region of the flow which is
absolutely unstable is sufficiently long for the flow to become globally
unstable. Perturbations grow exponentially in time until they are saturated by
nonlinearities. The resulting velocity fluctuations are dominated by the
vortex shedding. Once the flow has reached a limit cycle, the shedding
frequency is different from that predicted by resolvent analysis of the base
flow since the frequency of the least stable perturbations is altered during
the saturation process. Additionally, the mean recirculation region behind the
cylinder is shorter than its base flow counterpart in the streamwise
direction.

###  IV.3 Mean velocity profile

The focus of resolvent analysis typically shifts when using the mean velocity
profile rather than the base flow. The goal becomes identification of the
energetically important structures and their frequencies in the unsteady flow
rather than prediction of the external forcing and structure which appears
when the flow becomes unsteady.

Contours of the pseudospectrum corresponding to the mean flow in the 2D
laminar shedding regime are overlaid with the spectrum of the mean flow at
$Re=100$ in Figure 9. The resolvent norm along the imaginary axis is also
plotted alongside the second largest singular value. Similar to the base flow
resolvent, there is only one frequency at which there is a resonant peak.

(a)

(b)

Figure 9: (a) Spectrum (red dots) and pseudospectrum (filled contours,
$\epsilon$ increasing as colors change from dark to light) of the LNS operator
for the cylinder mean flow at $Re=100$. (b) The resolvent norm, $\sigma_{1}$
(solid line), i.e. the value of $\epsilon^{-1}$ along the imaginary axis,
second largest singular value $\sigma_{2}$ (dash-dotted line) and inverse
distance from the imaginary axis to the nearest eigenvalue (dotted red line).

Unlike the base flow case, the most amplified frequency at supercritical
Reynolds numbers correctly predicts the shedding frequency as seen in Figure
10. While the resolvent norm always peaks at a distinct frequency for all
cases, the growth rate of the least stable eigenvalue of the base flow
continues to grow while the frequency remains roughly constant. Figure 10
shows that the largest peak occurs at the stability limit, $Re_{c}$. The
maximum amplification, which here is proportional to the inverse distance
between the eigenvalue and the imaginary axis, indicates the progression of
the least stable pole across the complex plane and over the imaginary axis.
The resolvent norm has not been plotted for supercritical base flows since the
resolvent attempts to quantify the size of perturbation necessary for the
spectrum to cross the neutral axis. For the base flow at the critical Reynolds
number $Re_{c}$ and mean flows where $Re>Re_{c}$, the size of this
perturbation is very small leading to very highly amplified disturbances.

(a)

(b)

Figure 10: (a) The resolvent norm for the critical (solid red line),
subcritical (solid black lines), and supercritical (dotted black lines) base
flows. (b) The resolvent norm for supercritical mean flows.

The peak resolvent norm for the mean flows has no discernible pattern in
Figure 10 since the real part of the eigenvalue is approximately zero. It is
very sensitive, therefore, to the spatial resolution and temporal convergence
of the mean flow in addition to the discretization of $\omega$. Nevertheless,
proportionality between the resolvent and stability mode shapes can be
expected. Substituting $i\omega_{s}$ for $\lambda$ into Equation 57, since the
real part of the marginally stable mode is nearly zero, yields

$\mathcal{H}(\omega_{s})\approx\tilde{\bm{u}}_{s}\tilde{\bm{v}}_{s}^{*}\approx\hat{\bm{\psi}}_{s}\hat{\bm{\phi}}_{s}^{*}.$
(68)

The sum in Equation 57 is dominated by the contribution from the marginally
stable mode so the resolvent operator can be approximated by the outer product
of the marginally stable mode and its adjoint counterpart or the optimal
resolvent response and forcing modes at $\omega_{s}$. Similar to the base flow
case, we can quantify $|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}=26.9$
and this agrees fairly well with the ratio of
$\sigma_{1}(\omega_{s})|i\omega_{s}-\lambda_{\text{s}}|=28.9$ (see Table 2).

Rather than comparing the stability and resolvent mode shapes as shown in
Figure 7, the true and approximated wavemakers associated with the mean flow,
computed using stability and resolvent modes, respectively, are compared for
$Re=100$ in Figure 11(a,b). The stability wavemaker is in good agreement with
the one computed by Meliga  _et al._ (2016). A very close match between the
extent of the two wavemaker regions can be observed, implying that the
underlying modes are indeed proportional to each other. Streamlines from the
mean flow are superimposed to observe how the wavemaker is related to the mean
recirculation bubble, the size of which depends on Reynolds number. The length
of the recirculation bubble scales with the streamwise extent of the wavemaker
region for any Reynolds number for either the mean flow as seen in Figure
11(a,b) or base flow as seen in Figure 11(c,d). The advection term in the
operator $\overline{\bm{u}}\cdot\nabla()$ changes sign in the reverse flow
region behind the cylinder which then guarantees overlap between the adjoint
and forward eigenvectors, per Section III.5.

(a)

(b)

(c)

(d)

Figure 11: Wavemakers for mean flow at $Re=100$ computed from stability modes
(a) and resolvent modes (b). Wavemakers for the base flow at $Re=47$ (c) and
$Re=100$ (d). Blue lines superimpose the mean flow streamlines which delineate
the mean recirculation bubble.

The true (stability) wavemaker regions associated with the base profile are
shown in Figures 11(c) and (d) for $Re=47$ and $Re=100$, respectively. Figures
11(a) and (d) compare the mean and base wavemakers at $Re=100$. The main
difference between base and mean flows is that as the Reynolds number
increases, both the mean recirculation bubble Zielinska  _et al._ (1997) and
wavemaker region shrink, whereas for the unstable base flows, increasing the
Reynolds number will also increase the streamwise length of the recirculation
bubble and wavemaker region (Figures 11(c)-(d)). Therefore, the wavemaker
computed from the mean flow is a better indication of the region of
instability in the flow as it is responsible for the selection of the
nonlinear frequency and amplitude of the limit cycle.

###  IV.4 Resolvent modeling of (low Reynolds number) cylinder flow

In the case of flow around a circular cylinder, the resolvent identifies only
stability mechanisms for both the base and mean flow cases. There are no
frequencies at which amplification occurs due to pseudoresonance. This is
consistent with the work of Abdessemed  _et al._ (2009) who investigated
direct transient growth analysis to study its role in the primary and
secondary bifurcations of cylinder flow. Since only one mode becomes unstable,
the effect of the nonlinearities is to saturate the growth mechanism and this
alters the frequency of the structure (Barkley, 2006). While the base flow is
unable to predict the nonlinear frequency in the saturated state, it can be
correctly predicted from the mean flow. Thus in the case of stability
mechanisms leading to unsteadiness (e.g. cylinder flow, Rayleigh-Bénard
convection (Turton  _et al._ , 2015)), mean stability analysis is successful
at predicting the frequency of the unsteady flow.

##  V Application to wall turbulence

We now consider turbulent channel flow in the context of the discussion from
Section III. Unlike the cylinder flow, which is an oscillator with intrinsic
dynamics that are insensitive to background noise, this flow is an example of
a noise-amplifier; as such, pseudoresonance plays a big role and leads to
significant amplification at non-resonant frequencies. Due to its geometric
simplicity, we choose channel flow at $Re_{\tau}=2000$ which has a parallel
mean velocity profile $\overline{\bm{u}}=\overline{u}(y)$.

###  V.1 Numerical Methods

We again consider the non-dimensional, incompressible NSE where the channel
half-height $h$ and the friction velocity $u_{\tau}=\sqrt{\tau_{w}/\rho}$
(where $\tau_{w}$ is the wall shear stress, $\rho$ is the density) are the
characteristic scales to obtain

$\displaystyle\partial_{t}\bm{u}+\bm{u}\cdot\nabla\bm{u}$
$\displaystyle=-\nabla p+{Re_{\tau}}^{-1}\nabla^{2}\bm{u}$ (69a)
$\displaystyle\nabla\cdot\bm{u}$ $\displaystyle=0,$ (69b)

where $Re_{\tau}=hu_{\tau}/\nu$ and $\nu$ is the kinematic viscosity. The
streamwise direction and spanwise directions are periodic, and the wall-normal
domain extends from $y/h=-1$ to $y/h=1$ with no-slip and no-penetration
conditions imposed at the wall. The velocity field is Reynolds decomposed into
the sum of a spatio-temporal mean and fluctuations,

$\bm{u}(x,y,z,t)=\overline{u}(y)+\bm{u^{\prime}}(x,y,z,t).$ (70)

Here we assume the mean velocity profile is known a priori from an eddy
viscosity model Reynolds and Tiederman (1967) as discussed in Moarref  _et
al._ (2013). We express the fluctuations as Fourier modes in the
streamwise/spanwise directions and in time,

$\hat{\bm{u}}(k_{x},k_{z},\omega;y)=\displaystyle\int_{-\infty}^{\infty}\displaystyle\int_{-\infty}^{\infty}\displaystyle\int_{-\infty}^{\infty}\bm{u}^{\prime}(x,y,z,t)e^{-i(k_{x}x+k_{z}z-\omega
t)}\mathrm{d}x\mathrm{d}z\mathrm{d}t,$ (71)

where $k_{x}$ is the streamwise wavenumber, $k_{z}$ is the spanwise
wavenumber, and $\omega$ is the radial frequency. Upon elimination of the
pressure term, we can express the governing equations in terms of the
fluctuating vertical velocity $\hat{v}$ and normal vorticity
$\hat{\eta}=ik_{z}\hat{u}-ik_{x}\hat{w}$,

$-i\omega\left(\begin{array}[]{c}\hat{v}\\\
\hat{\eta}\end{array}\right)+\left(\begin{array}[]{c
c}k^{2}-\mathcal{D}^{2}&0\\\ 0&1\end{array}\right)^{-1}\left(\begin{array}[]{c
c}\mathcal{L}_{OS}&0\\\
ik_{z}\overline{u}^{\prime}&\mathcal{L}_{SQ}\end{array}\right)\left(\begin{array}[]{c}\hat{v}\\\
\hat{\eta}\end{array}\right)=\bm{J}\hat{\bm{f}},$ (72)

where the Orr-Sommerfeld (OS) and Squire (SQ) operators are given by

$\displaystyle\mathcal{L}_{OS}=ik_{x}\overline{u}(k^{2}-\mathcal{D}^{2})+ik_{x}\overline{u}^{\prime\prime}+\frac{1}{Re_{\tau}}(k^{2}-\mathcal{D}^{2})^{2},$
(73)
$\displaystyle\mathcal{L}_{SQ}=ik_{x}\overline{u}+\frac{1}{Re_{\tau}}(k^{2}-\mathcal{D}^{2}),$
(74)

and

$\bm{J}=\left(\begin{array}[]{c c}k^{2}-\mathcal{D}^{2}&0\\\
0&1\end{array}\right)^{-1}\left(\begin{array}[]{c c
c}-ik_{x}\mathcal{D}&-k^{2}&-ik_{z}\mathcal{D}\\\
ik_{z}&0&-ik_{x}\end{array}\right),$ (75)
$\hat{\bm{f}}=\left(\begin{array}[]{c}\hat{f}_{u}\\\ \hat{f}_{v}\\\
\hat{f}_{w}\end{array}\right)=-\langle\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}\rangle_{\bm{k}}.$
(76)

Here $\mathcal{D}=\frac{\partial}{\partial{y}}$, $k^{2}=k_{x}^{2}+k_{z}^{2}$,
and $\langle\hskip 5.69054pt\rangle_{\bm{k}}$ denotes the Fourier component
associated with the wavenumber vector $\bm{k}=(k_{x},k_{z},\omega)$. The wall-
normal operators are discretized numerically with Chebyshev collocation points
using the suite developed by Weideman and Reddy (2000). We can recast Equation
72 into the following input/output form

$\left(\begin{array}[]{c}\hat{u}\\\ \hat{v}\\\
\hat{w}\end{array}\right)=\mathcal{H}(k_{x},k_{z},\omega)\left(\begin{array}[]{c}\hat{f}_{u}\\\
\hat{f}_{v}\\\ \hat{f}_{w}\end{array}\right),$ (77)

where the resolvent operator $\mathcal{H}$ is given by

$\mathcal{H}(k_{x},k_{z},\omega)=\bm{K}(-i\omega+\bm{L})^{-1}\bm{J},$ (78)

where

$\bm{L}=\bm{G}^{-1}\mathcal{L},$ (79) $\bm{G}=\left(\begin{array}[]{c
c}k^{2}-\mathcal{D}^{2}&0\\\ 0&1\end{array}\right),$ (80)
$\mathcal{L}=\left(\begin{array}[]{c c}\mathcal{L}_{OS}&0\\\
ik_{z}\overline{u}^{\prime}&\mathcal{L}_{SQ}\end{array}\right),$ (81)
$\bm{K}=\frac{1}{k^{2}}\left(\begin{array}[]{c c
c}ik_{x}\mathcal{D}&-ik_{z}\\\ k^{2}&0\\\
ik_{z}\mathcal{D}&ik_{x}\end{array}\right).$ (82)

As before, we can decompose the resolvent operator via the SVD as

$\mathcal{H}(k_{x},k_{z},\omega)=\bm{\Psi}(k_{x},k_{z},\omega)\bm{\Sigma}(k_{x},k_{z},\omega)\bm{\Phi}^{*}(k_{x},k_{z},\omega)$
(83)

The off-diagonal term in $\bm{L}$ is proportional to the mean shear
$\overline{u}^{\prime}$ which is maximum at the wall. It remains large in the
inner region before it begins to decline in the log region. Similar to the
model LNS operator in Section III, mean shear is the primary source of non-
normality leading to significant amplification. Its spatial variation is
important since it has been shown McKeon and Sharma (2010) that a critical-
layer mechanism tends to localize activity at the wall-normal location where
the phase speed of the disturbance is equal to the local mean velocity. This
is explored further by considering three particular wavenumber triplets that
are representative of the near-wall cycle, a very large-scale motion (VLSM),
and a stationary disturbance. The roles of normal and non-normal mechanisms
are studied by analyzing the mode shapes of the leading resolvent response
modes, the pseudospectrum of the LNS operator, and the resolvent norm compared
with the inverse distance between the imaginary axis and the nearest
eigenvalue.

(a)

(b)

Figure 12: (a) Turbulent mean velocity profile for channel flow at
$Re_{\tau}=2000$ plotted in inner units alongside (b) the mean shear.

###  V.2 Near-wall cycle

The first wavenumber combination considered is
$(k_{x},k_{z},c^{+})=(4\pi,40\pi,14)$ which is representative of the near-wall
cycle McKeon and Sharma (2010). Here the wavespeed is given by
$c^{+}=\omega/k_{x}$. Figure 12 shows that the mean shear is very large at the
wall-normal height where the wavespeed matches the local mean, resulting in
the off-diagonal terms of the resolvent operator being large. This is similar
to the model LNS operator in Equation 44 where the influence of non-normality
concentrates energy in different velocity components for $\hat{\bm{\psi}}_{1}$
and $\hat{\bm{\phi}}_{1}$. The optimal resolvent forcing and response modes
are plotted in Figure 13 to illustrate that the forcing is primarily
concentrated in $v$ and $w$ while the response is mostly in $u$.

(a)

(b)

Figure 13: Velocity amplitudes for the optimal forcing mode $\hat{\phi}_{1}$
in (a) and optimal response mode $\hat{\psi}_{1}$ in (b) corresponding to the
wavenumber triplet of $(k_{x},k_{z},c^{+})=(4\pi,40\pi,14)$.

The strength of mean shear suggests that pseudoresonance is the primary driver
of the near-wall cycle mode. The spectrum as well as contours of the
pseudospectrum are plotted in Figure 14 for various $\epsilon$. Figure 14 also
includes the resolvent norm and the contribution from resonance which is
quantified by $|i\omega-\lambda|^{-1}$. The ratio of the resolvent norm to the
contribution from resonance is 19.6, which is of the same order of magnitude
as the value predicted by
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}=4.81$ (see Table 2).
Nevertheless, it is clear from this discrepancy that amplification cannot be
attributed to one particular eigenvalue and that there is no eigenvector which
is proportional to the resolvent mode. Using the expression
$\sigma_{1}|i\omega-\lambda|$ to quantify non-normality is problematic since
there is no unique $\lambda$ which is responsible for amplification.

The $\omega$ corresponding to $c^{+}=14$ is indicated by the horizontal,
dashed blue line. At this frequency, the resolvent norm is significantly
larger than the resonance term suggesting that amplification is due to non-
normal mechanisms. This observation is confirmed by the spectrum where the
least damped eigenvalues are clustered around higher frequencies and
contributes to the resolvent norm for $\omega>300$.

Figure 14: The eigenvalues of the operator $\bm{L}(k_{x}=4\pi,k_{z}=40\pi)$ in
red circles overlaid with contours of the pseudospectrum (left). The resolvent
norm is plotted in the solid black line along with the inverse distance from
the imaginary axis to the nearest eigenvalue which is the red dotted line
(right). The spatial wavenumbers correspond to the near-wall cycle, and the
horizontal, dashed blue line represents the $\omega$ with the largest
resolvent norm, which corresponds to a phase speed of $c^{+}\approx 14$.

It is also worth noting that the eigenvalues are significantly damped, which
results in the leading singular values being on the order of unity. While
these are not large, the rank-1 approximation is still valid since the first
pair of singular values is approximately one order of magnitude larger than
the others as seen in Figure 15.

Figure 15: First 20 singular values of the resolvent operator for
$(k_{x},k_{z},c^{+})=(4\pi,40\pi,14)$.

###  V.3 VLSM

Further from the wall, the mean shear drops several orders of magnitude (see
Figure 12) and the effect of viscosity decreases. As a result, amplification
becomes a mix of both normal and non-normal effects. To reinforce this point,
a wavenumber triplet which is representative of a VLSM is considered. The mode
shapes for $(k_{x},k_{z},c^{+})=(\pi/9,2\pi/3,22)$ are plotted in Figure 16.

(a)

(b)

Figure 16: Velocity amplitudes for the optimal forcing mode $\hat{\phi}_{1}$
in (a) and optimal response mode $\hat{\psi}_{1}$ in (b) corresponding to the
wavenumber triplet of $(k_{x},k_{z},c^{+})=(\pi/9,2\pi/3,22)$.

The forcing is dominated by the $w$-component while the response is dominated
by the $u$-component. The $v$-component of the forcing, notably, is less
significant than the near-wall cycle mode implying that the role of lift-up is
not as pronounced for this mode. The spectrum associated with the streamwise
and spanwise wavenumbers of the VLSM is plotted along with the pseudospectrum
of the LNS operator in Figure 17.

Figure 17: The eigenvalues of the operator $\bm{L}(k_{x}=\pi/9,k_{z}=2\pi/3)$
in red circles overlaid with contours of the pseudospectrum (left). The
resolvent norm is plotted in the solid black line along with the inverse
distance from the imaginary axis to the nearest eigenvalue which is the red
dotted line (right). The spatial wavenumbers correspond to the VLSM mode, and
the horizontal, dashed blue line represents the $\omega$ with the largest
resolvent norm, which corresponds to a phase speed of $c^{+}\approx 22$.

The results are drastically different from the near-wall cycle case as the
resonant forcing of eigenvalues is greater than one so amplification is due to
both terms on the right-hand side of Equation 31. Notably, the spectrum in
Figure 17 resembles that of the base flow at $Re=10,000$ for $k_{x}=1$,
$k_{z}=0$ as there are three distinct branches. As observed by Reddy  _et al._
(1993) and Schmid and Henningson (2001), the eigenvalues at the intersection
of the branches are the most sensitive to perturbations and result in very
large non-normal amplification. The product $\sigma_{1}|i\omega-\lambda|=146$
while $\|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}\|^{-1}$ = 34.9 (see Table
2), suggesting that the nonorthogonality of many eigenfunctions leads to high
pseudoresonance. Similar to the near-wall cycle case, there are no eigenvalues
which exactly match the wave speed associated with the VLSM. The eigenvalue
close to the dotted blue line in Figure 17, however, does seem to impact the
resolvent norm which has an extra bump near its maximum value. This is also
reflected in the red dotted line since the eigenvalue protrudes from the
distinct Y-shape of the spectrum. The maximum singular values are on the order
of $10^{3}$ and the resolvent operator is low-rank as seen in Figure 18.

Figure 18: First 20 singular values of the resolvent operator for
$(k_{x},k_{z},c^{+})=(\pi/9,2\pi/3,22)$.

###  V.4 Stationary disturbances

Finally, we consider stationary disturbances which tend to be the most
amplified by the resolvent operator with singular values exceeding $10^{4}$.
The specific wavenumber triplet selected for this study is
$(k_{x},k_{z},\omega)=(0,2\pi/3,0)$. The roots behind such large amplification
can be traced back to the model operator in Equation 44. In this example, the
$\partial\overline{u}/\partial y\to\infty$ resulting in a low-rank system
which concentrated all the forcing energy in the second velocity component and
the response energy in the first velocity component. When $k_{x}=\omega=0$,
all of the diagonal terms of the resolvent become order $\epsilon$ small since
imaginary terms are eliminated and $\mathcal{D}$ scales with $Re^{-1}$. Thus,
when the LNS operator is inverted, the determinant, which is the product of
the diagonal terms, is very small. The energy for the forcing is almost
totally in the wall-normal and spanwise components as seen in Figure 19 while
the response is almost totally in the streamwise component.

(a)

(b)

Figure 19: Velocity amplitudes for the optimal forcing mode $\hat{\phi}_{1}$
in (a) and optimal response mode $\hat{\psi}_{1}$ in (b) corresponding to the
wavenumber triplet of $(k_{x},k_{z},\omega)=(0,2\pi/3,0)$.

Similar to the near-wall cycle and VLSM modes, the spectrum and contours of
the pseudospectrum are presented in Figure 20 alongside the resolvent norm and
contribution from resonance. All of the eigenvalues are real since the
imaginary terms are eliminated from the resolvent operator when $k_{x}=0$.
Another implication of eliminating mean flow advection, as mentioned by Hack
and Moin (2017), is that the Orr mechanism is absent from this mode.

Figure 20: The eigenvalues of the operator $\bm{L}(k_{x}=0,k_{z}=2\pi/3)$ in
red circles overlaid with contours of the pseudospectrum (left). The resolvent
norm is plotted in the solid black line along with the inverse distance from
the imaginary axis to the nearest eigenvalue which is the red dotted line
(right). The horizontal, dashed blue line represents the $\omega$
corresponding to the phase speed of $c^{+}=0$.

Stationary disturbances are highly amplified and the singular values are
plotted in Figure 21. The rank-1 approximation is quite applicable for this
particular mode as the leading pair of singular values is on the order of
$10^{5}$. Moreover, the contribution from non-normality is well approximated
by $|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}=250$ which agrees quite
well with $\sigma_{1}|i\omega-\lambda|$ = 278. Such agreement can be
attributed to the eigenvalue closest to the imaginary axis which has an
imaginary component that agrees with the most amplified frequency. Unlike the
cylinder case where there exists a convective-type non-normality, the
$k_{x}=0$ modes are an example of the component-type non-normality and so
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|$ is small since the velocity for
the forcing mode is almost completely concentrated in the wall-normal plane
while the velocity for the response mode is nearly all in the streamwise
direction.

Figure 21: First 20 singular values of the resolvent operator for
$(k_{x},k_{z},\omega)=(0,2\pi/3,0)$.

###  V.5 Predicting low-rank nature of resolvent

We briefly return to the discussion of predicting when the resolvent operator
may be low-rank. When the LNS operator is inverted, the terms which were
originally small along the diagonal become large leading to a large resolvent
norm. Marginally stable eigenvalues or stationary disturbances have large
amplification since either one term along the diagonal is small when forcing
at a resonant frequency or all the terms along the diagonal are small since
$k_{x}=\omega=0$. The system is low-rank roughly speaking if the largest entry
of the resolvent operator is appreciably larger than the other terms. This can
happen when the off-diagonal are terms are large with respect to the diagonal
as is the case for the near-wall cycle and stationary disturbance modes or
when one term along the diagonal is small as is the case for the cylinder. If
the frequency of the forcing excites two eigenvalues then two diagonal terms
of the resolvent operator are large and the system is not rank-1 even though
there might be significant amplification (see Appendix A).

###  V.6 Influence of spatial wavenumber and wave speed

Based on the findings of this study and observations from Bourguignon (2012),
it is possible to categorize amplification mechanisms in wall-bounded
turbulence as either normal or non-normal depending on the wavenumber vector
$\bm{k}$ selected. When $k_{x}$ is small, the influence of both mean flow
advection and viscosity is diminished resulting in a non-normal system where
there is high amplification. Low-order modes (those corresponding to the
largest singular values), experience high amplification due to both normality
and non-normality. Higher-order modes also experience amplification due to
non-normality. At higher $k_{x}$, only low-order modes are amplified as long
as they are localized near the critical layer. Higher-order modes experience
low amplification which is proportional to viscosity. The wall-normal height,
furthermore, has implications on the type of amplification as it influences
the choice of wave speed $c$, or $\omega$, as well as the influence of mean
shear $\partial\overline{u}/\partial y$. Closer to the wall, the mean shear is
highest while in the log region, mean shear still plays an important role
resulting in preferential amplification of long streamwise structures.

|  $\bm{k}$ or $\omega$ | $\sigma_{1}$ | $|i\omega-\lambda|^{-1}$ | $\sigma_{1}|i\omega-\lambda|$ | $|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}$
---|---|---|---|---|---
Cylinder Base |  $\omega=0.743$ |  5.78e04 |  729 |  79.3 |  79.4
Cylinder Mean |  $\omega=1.02$ |  1.65e04 |  570 |  28.9 |  26.9
Near-wall cycle |  $(k_{x},k_{z},c^{+})=(4\pi,40\pi,14)$ |  0.502 |  2.56e-02 |  19.6 |  4.81
VLSM |  $(k_{x},k_{z},c^{+})=(\pi/9,2\pi/3,22)$ |  479 |  3.28 |  146 |  34.9
Stationary Disturbance |  $(k_{x},k_{z},\omega)=(0,2\pi/3,0)$ |  8.14e04 |  293 |  278 |  250
Table 2: Quantification of non-normality for the most amplified modes in
cylinder and turbulent channel flows.

##  VI Closing remarks

###  VI.1 Conclusions

The preceding development and example operators provide a context in which to
interpret the recent findings concerning analysis of the resolvent in the
literature. Amplification quantified by the resolvent norm incorporates both
normal and non-normal origins. In the limiting case of a normal operator, the
condition number $\kappa=1$ and resonance due to the eigenvalue is captured.
For a non-normal operator, $\kappa\gg 1$ and pseudoresonant effects can lead
to large amplification. The nature of the non-normality in the resolvent
contains information about the transfer mechanisms from inputs to outputs and
may impact the distribution of energy among different velocity components. The
inverse of the inner product $|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|$,
which is unity for a purely normal mechanism and zero for a purely non-normal
mechanism, quantifies the degree to which the amplification is non-normal. In
the presence of mean shear, there is always some degree of non-normality so a
mechanism can be considered normal in character if
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|$ is close to unity.

The model resolvent operator is also non-self-adjoint unless there is no mean
shear, the eigenvalues are real, and the disturbances are stationary
($\omega=0$). If the eigenvalues are purely imaginary and there is no spatial
dependence, the resolvent forcing and response modes are $90^{\circ}$ out of
phase. In a spatially-developing flow, this translates to the spatial support
of the forcing mode being upstream that of the response mode. It is also
possible to observe the Orr mechanism which aligns the resolvent response and
forcing modes with and against the mean shear, respectively. If there is
reverse flow, the modes overlap and form a wavemaker which is approximately
coincident with the region of the flow that is absolute unstable. Non-
normality can still be quantified via
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}$, which tends to infinity
for a purely convective-type instability mechanism. It is concluded that mode
shapes and locations can be generalized depending on the strength of lift-up
and mean flow advection.

These findings are applicable to both base and mean velocity profiles, but
there is an important distinction. For a base flow, the input forcing must be
provided by an external source. Therefore, the most amplified structure is of
interest as it is likely to be present in the flow when it becomes unsteady.
For the cylinder flow (Barkley, 2006) and thermo-solutal convection (Turton
_et al._ , 2015), the operator is normal in character so the most amplified
structure corresponds to an eigenvalue crossing the imaginary axis. This is
not the case for non-normal operators (e.g. Jovanović and Bamieh, 2005) as the
highest amplification occurs for modes which do not correspond to an
eigenvalue. For a mean flow, the nonlinear term is the source of intrinsic
forcing so the goal is to predict the structure of the unsteady flow (Taira
_et al._ , 2017). In many cases, the most amplified response is normal in
character and accounts for a significant amount of the kinetic energy of the
velocity fluctuations (e.g. McKeon and Sharma, 2010; Beneddine  _et al._ ,
2016). There are circumstances, however, where it does not give the complete
picture as the rank-1 approximation is no longer valid and suboptimal modes
need to be considered.

Expressing the resolvent operator as the outer products of eigenmodes
clarifies the connection between stability and resolvent analysis.
Furthermore, it becomes apparent that mean stability analysis is successful
when the amplification is due to forcing at a frequency in the vicinity of a
unique eigenvalue. In such circumstances, the forward and adjoint stability
modes are proportional to the optimal forcing and response modes at that
frequency. Both sets of modes form a rank-1 approximation, which fails if the
eigenvalues are clustered together, of the resolvent operator. Cylinder flow
is chosen to demonstrate the similarity of stability and resolvent modes when
resonance outweighs pseudoresonance at a particular frequency. In both the
base flow and mean flow cases, amplification occurs at a single frequency
corresponding to the imaginary part of the least stable eigenvalue. The
resolvent modes approximate the region of the flow which is absolutely
unstable by identifying the wavemaker. There is no wavemaker at very low
Reynolds numbers when the flow is only convectively unstable. The cylinder
exemplifies a convective-type non-normality where mean flow advection
separates the spatial support of the forcing mode to be upstream of the
response mode as long as $\overline{u}>0$. Non-normality quantified by
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}$ agrees well with
$\sigma_{1}|i\omega-\lambda|$ for both the base and mean flows.

It is important to note that only normal mechanisms are active in cylinder
flow whereas in more complicated flows such as wall-bounded turbulence, both
normal and non-normal mechanisms are relevant. Three different wavenumber
triplets, representative of the near-wall cycle, VLSM’s, and stationary
disturbances, highlight the competing influences of viscous dissipation, mean
flow advection, and mean shear on not only the most amplified modes but also
the spectrum and pseudospectrum. The importance of each term depends
significantly on the wall-normal height where the perturbations are localized.
In the inner region where there is very high mean shear and viscosity is most
important, amplification is primarily due to pseudoresonant mechanisms.
Forcing in the spanwise/wall-normal plane leads to a large response in the
streamwise direction as seen for the near-wall cycle mode. The eigenvalues are
highly damped yet the resolvent norm is on the order of unity due to the high
sensitivity of the spectrum to perturbations. In the log region, mean shear
and consequently lift-up are weaker yet the declining importance of viscosity
results in eigenvalues which are closer to the imaginary axis. The most
amplified disturbance corresponding to the wavespeed of the VLSM,
nevertheless, is primarily due to pseudoresonance. Consequently, there is
poorer agreement between non-normality quantified by
$|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}$ and
$\sigma_{1}|i\omega-\lambda|$ as it is clear that amplification can no longer
be attributed to a single eigenvalue in the spectrum. Mean flow advection
results in the Orr mechanism (see McKeon, 2017) and hence less overlap between
the forcing and response modes.

Stationary disturbances are the globally most amplified disturbances by
effectively leveraging mean shear. The perturbation energy is almost
exclusively concentrated in the $v-$ and $w$-components of the forcing mode
and in the $u$-component of the response mode. Assuming streamwise constant
disturbances eliminates the mean flow advection term from the resolvent
operator and hence suppresses the Orr mechanism. All of the non-normality,
consequently, can be classified as a component-type non-normality, in contrast
to the cylinder flow, and the eigenvalues of the LNS operator are real. Non-
normality quantified by $|\hat{\bm{\phi}}_{1}^{*}\hat{\bm{\psi}}_{1}|^{-1}$
agrees well with $\sigma_{1}|i\omega-\lambda|$ since amplification can be
attributed to the eigenvalue closest to the imaginary axis. For a generic
wavenumber triplet, the applicability of the rank-1 approximation can be
approximated by the ratio of the largest term in the resolvent operator to the
other entries. Finally, the distribution of energy among various velocity
components may be useful when considering how the nonlinear term, which can be
computed from resolvent response modes (McKeon  _et al._ , 2013), projects
onto the optimal resolvent forcing modes.

###  VI.2 Reduced-order modeling and control

In the pseudo-resonance limit case discussed in Section III, it was
demonstrated that a single resolvent forcing/response pair could characterize
the most amplified disturbance. The same could not be said of the eigenvectors
which were nearly parallel to one another. In reduced-order modeling
applications, the resolvent modes could be a superior basis to the stability
modes, since fewer are needed to represent the most important behavior of the
flow. In the context of control, where reduced-order models play a key role in
rendering the problem computationally tractable, stability modes have been
shown by Barbagallo  _et al._ (2009) to be less effective at representing the
input-output behavior of the system. While their study did not consider
resolvent modes as an expansion basis for model reduction, their capability of
capturing the intrinsic dynamics (stability-related phenomena) as well as non-
modal growth could make them an ideal candidate. Resolvent modes can also
provide relevant information in regard to where to place sensors and
actuators.

The resolvent response modes are outputs which are highly amplified by the
linear dynamics of the NSE, and so sensors could be placed where these are
likely to be strong. The resolvent forcing modes are the ‘trigger’ or input
which leads to high amplification, and so the actuators could be placed to
manipulate the flow in such a way that suppresses these disturbances.

Note, however, that in cases where there is large spatial separation between
resolvent forcing and response modes, it is possible that improved performance
could be attained by sensing and actuating within a wavemaker region (e.g.
Chen and Rowley, 2011; Giannetti and Luchini, 2007), which, as discussed in
Section IV, may also be estimated from resolvent analysis. While the resolvent
decomposition shows potential for control applications (Luhar  _et al._ ,
2014), further refinements could seek to balance the observability and
controllability of the reduced-order model (Moore, 1981; Rowley, 2005),
subject to known sensor and actuator locations, and information about the
nature of the nonlinear forcing.

##  VII Acknowledgments

The work in this study has been financially supported by a National Science
Foundation Graduate Fellowship, AFOSR under grant number FA 9550-16-1-0361 and
ARO under grant number W911NF-17-1-0306. The authors would like to thank Denis
Sipp for providing the resolvent code and Andres Goza who assisted setting it
up on a cluster. Finally, the authors would like to acknowledge Theresa
Saxton-Fox and Ryan McMullen for useful feedback on the manuscript.

##  Appendix A Rank-1 example

An important result of this paper is that the rank-1 approximation fails when
there is not sufficient separation of eigenvalues at the frequency of
interest. The two operators $\bm{T}_{1}$ and $\bm{T}_{2}$ in Equation 84,
which are nearly identical except for the frequency of the second eigenvalue,
illustrate this point. In the case of $\bm{T}_{2}$, substituting the
eigenvalues into Equation 57 results in two terms of roughly equal importance
leading to the failure of the rank-1 approximation as seen in Figure 22 where
the singular values of the resolvent operator for both $\bm{T}_{1}$ and
$\bm{T}_{2}$ are plotted. The rank-1 approximation works better for
$\bm{T}_{1}$ since there is a clear separation of singular values at the
resonant frequencies.

$\bm{T}_{1}=\left(\begin{array}[]{cc}-0.2+0.1i&0\\\
0&-0.22+\boxed{1.12i}\end{array}\right),~{}~{}~{}\bm{T}_{2}=\left(\begin{array}[]{cc}-0.2+0.1i&0\\\
0&-0.22+\boxed{0.12i}\end{array}\right).$ (84)

(a)

(b)

Figure 22: The singular values as a function of the forcing frequency for the
operators $\bm{T}_{1}$ (left) and $\bm{T}_{2}$ (right).

## References

  * Turton  _et al._ (2015) S. E. Turton, L. S. Tuckerman, and D. Barkley, “Prediction of frequencies in thermosolutal convection from mean flows,” Phys. Rev. E 91 (2015).
  * Gudmundsson and Colonius (2011) K. Gudmundsson and T. Colonius, “Instability wave models for the near-field fluctuations of turbulent jets,” J. Fluid Mech. 689, 97–128 (2011).
  * Oberleithner  _et al._ (2014) K. Oberleithner, L. Rukes, and J. Soria, “Mean flow stability analysis of oscillating jet experiments,” J. Fluid Mech. 757, 1–32 (2014).
  * Schmidt  _et al._ (2017a) O. T. Schmidt, A. Towne, T. Colonius, A. V. G. Cavalieri, P. Jordan, and G. A. Brès, “Wavepackets and trapped acoustic modes in a turbulent jet: coherent structure eduction and global stability,” J. Fluid Mech. 825, 1153–1181 (2017a).
  * Beneddine  _et al._ (2016) S. Beneddine, D. Sipp, A. Arnault, J. Dandois, and L. Lesshafft, “Conditions for validity of mean flow stability analysis,” J. Fluid Mech. 798, 485–504 (2016).
  * Barkley (2006) D. Barkley, “Linear analysis of the cylinder wake mean flow,” Europhys. Lett. 75, 750–756 (2006).
  * Provansal  _et al._ (1987) M. Provansal, C. Mathis, and L. Boyer, “Bénard-von Kármán instability: transient and forced regimes,” J. Fluid Mech. 182, 1–22 (1987).
  * Sreenivasan  _et al._ (1987) K. R. Sreenivasan, P. J. Strykowski, and D. J. Olinger, “Hopf bifurcation, Landau equation, and vortex shedding behind circular cylinders,” in  _Proc. Forum on Unsteady Flow Separation_ , edited by K. Ghia (1987) pp. 1–13.
  * Noack and Eckelmann (1994) B. R. Noack and H. Eckelmann, “A global stability analysis of the steady and periodic cylinder wake,” J. Fluid Mech. 270, 297–330 (1994).
  * Sipp and Lebedev (2007) D. Sipp and A. Lebedev, “Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows,” J. Fluid Mech. 593, 333–358 (2007).
  * Mantič-Lugo  _et al._ (2014) V. Mantič-Lugo, C. Arratia, and F. Gallaire, “Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake,” Phys. Rev. Lett. 113 (2014).
  * Trefethen  _et al._ (1993) L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll, “Hydrodynamic stability without eigenvalues,” Science 261, 578–584 (1993).
  * Schmid and Henningson (2001) P. J. Schmid and D. S. Henningson,  _Stability and transition in shear flows_ (Springer, 2001).
  * Jovanović and Bamieh (2005) M. R. Jovanović and B. Bamieh, “Componentwise energy amplification in channel flows,” J. Fluid Mech. 534, 145–183 (2005).
  * Farrell and Ioannou (1993) B. F. Farrell and P. J. Ioannou, “Stochastic forcing of the linearized Navier-Stokes equations,” Phys. Fluids 5, 2600–2609 (1993).
  * Butler and Farrell (1992) K. M. Butler and B. F. Farrell, “Three-dimensional optimal perturbations in viscous shear flow,” Phys. Fluids A 4, 1637–1650 (1992).
  * McKeon and Sharma (2010) B. J. McKeon and A. S. Sharma, “A critical-layer framework for turbulent pipe flow,” J. Fluid Mech. 658, 336–382 (2010).
  * Hwang and Cossu (2010) Y. Hwang and C. Cossu, “Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow,” J. Fluid Mech. 664, 51–73 (2010).
  * Lu and Papadakis (2014) L. Lu and G. Papadakis, “An iterative method for the computation of the response of the linearised Navier-Stokes equations to the harmonic forcing and application to forced cylinder wakes,” Int. J. Numer. Methods Fluids 74, 794–817 (2014).
  * Jeun  _et al._ (2016) J. Jeun, J. W. Nichols, and M. R. Jovanović, “Input-output analysis of high-speed axisymmetric isothermal jet noise,” Phys. Fluids 28 (2016).
  * Schmidt  _et al._ (2017b) O. T. Schmidt, A. Towne, G. Rigas, T. Colonius, and G. A. Brès, “Spectral analysis of jet turbulence,” arXiv preprint arXiv:1711.06296 (2017b).
  * Taira  _et al._ (2017) K. Taira, S. L. Brunton, S. T. M. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal analysis of fluid flows: an overview,” AIAA J. 55, 4013–4041 (2017).
  * Marquet  _et al._ (2009) O. Marquet, M. Lombardi, J. M. Chomaz, D. Sipp, and L. Jacquin, “Direct and adjoint global modes of a recirculation bubble: lift-up and convective non-normalities,” J. Fluid Mech. 622, 1–21 (2009).
  * Brandt  _et al._ (2011) L. Brandt, D. Sipp, J. O. Pralits, and O. Marquet, “Effect of base-flow variation in noise amplifiers: the flat-plate boundary layer,” J. Fluid Mech. 687, 503–528 (2011).
  * Orr (1907) W. Mc. F. Orr, “The stability or instability of the steady motions of a perfect liquid and of a viscous liquid,” Proc. R. Irish Acad. 27, 9–138 (1907).
  * Farrell (1987) B. F. Farrell, “Developing disturbances in shear,” J. Atmos. Sci. 44, 2191–2199 (1987).
  * Pier and Huerre (2001) B. Pier and P. Huerre, “Nonlinear self-sustained structures and fronts in spatially developing wake flows,” J. Fluid Mech. 435, 145–174 (2001).
  * Huerre and Monkewitz (1985) P. Huerre and P. A. Monkewitz, “Absolute and convective instabilities in free shear layers,” J. Fluid Mech. 159, 151–168 (1985).
  * Reddy  _et al._ (1993) S. C. Reddy, P. J. Schmid, and D. S. Henningson, “Pseudospectra of the Orr-Sommerfeld operator,” SIAM J. Appl. Math. 53, 15–47 (1993).
  * Trefethen  _et al._ (1999) A. E. Trefethen, L. N. Trefethen, and P. J. Schmid, “Spectra and pseudospectra for pipe Poiseuille flow,” Comput. Methods in Appl. Mech. Eng. 175, 413–420 (1999).
  * Schmid (2007) P. J. Schmid, “Nonmodal stability theory,” Annu. Rev. Fluid Mech. 39, 129–162 (2007).
  * Schmid and Brandt (2014) P. J. Schmid and L. Brandt, “Analysis of fluid systems: stability, receptivity, sensitivity,” Appl. Mech. Rev. 66 (2014).
  * Gómez  _et al._ (2016) F. Gómez, H. M. Blackburn, M. Rudman, A. S. Sharma, and B. J. McKeon, “A reduced-order model of three-dimensional unsteady flow in a cavity based on the resolvent operator,” J. Fluid Mech. 798, R2 (2016).
  * Huerre and Rossi (1998) P. Huerre and M. Rossi, “Hydrodynamic instabilities in open flows,” in  _Hydrodynamics and nonlinear instabilities_ , edited by C. Godrèche and P. Manneville (Cambridge University Press, 1998).
  * Monkewitz (1988) P. A. Monkewitz, “The absolute and convective nature of instability in two-dimensional wakes at low Reynolds numbers,” Phys. Fluids 31 (1988).
  * Giannetti and Luchini (2007) F. Giannetti and P. Luchini, “Structural sensitivity of the first instability of the cylinder wake,” J. Fluid Mech. 581, 167–197 (2007).
  * Chomaz (2005) J. M. Chomaz, “Global instabilities in spatially developing flows: non-normality and nonlinearity,” Annu. Rev. Fluid Mech. 37, 357–392 (2005).
  * Sipp  _et al._ (2010) D. Sipp, O. Marquet, P. Meliga, and A. Barbagallo, “Dynamics and control of global instabilities in open-flows: a linearized approach,” Appl. Mech. Rev. 63 (2010).
  * McKeon (2017) B. J. McKeon, “The engine behind (wall) turbulence: perspectives on scale interactions,” J. Fluid Mech. 817, P1 (2017).
  * Herbert (1997) T. Herbert, “Parabolized stability equations,” Annu. Rev. Fluid Mech. 29, 245–283 (1997).
  * Reynolds and Tiederman (1967) W. C. Reynolds and W. G. Tiederman, “Stability of turbulent channel flow, with application to Malkus’s theory,” J. Fluid Mech. 27, 253–272 (1967).
  * Luchini and Bottaro (2014) P. Luchini and A. Bottaro, “Adjoint equations in stability analysis,” Annu. Rev. Fluid Mech. 46, 493–517 (2014).
  * Moarref  _et al._ (2013) R. Moarref, A. S. Sharma, J. A. Tropp, and B. J. McKeon, “Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels,” J. Fluid Mech. 734, 275–316 (2013).
  * Gebhardt and Grossmann (1994) T. Gebhardt and S. Grossmann, “Chaos transition despite linear stability,” Phys. Rev. E 50, 3705 (1994).
  * Landahl (1980) M. T. Landahl, “A note on an algebraic instability of inviscid parallel shear flows,” J. Fluid Mech. 98, 243–251 (1980).
  * Juniper (2012) M. Juniper, “Absolute and convective instability in gas turbine fuel injectors,” in  _Proceedings of ASME Turbo Expo 2012: Power for Land, Sea and Air, Copenhagen_ (ASME, New York, 2012).
  * Qadri and Schmid (2017) U. A. Qadri and P. J. Schmid, “Frequency selection mechanisms in the flow of a laminar boundary layer over a shallow cavity,” Phys. Rev. Fluids 2 (2017).
  * Gelfand (1941) I. Gelfand, “Zur Theorie der Charaktere der Abelschen topologischen Gruppen,” Rec. Math. N. S. 9, 49–50 (1941).
  * Barkley and Henderson (1996) D. Barkley and R. D. Henderson, “Three-dimensional Floquet stability analysis of the wake of a circular cylinder,” J. Fluid Mech. 322, 215–241 (1996).
  * Hecht (2012) F. Hecht, “New development in FreeFem++,” J. Numer. Math. 20, 251–265 (2012).
  * Nayar and Ortega (1993) M. Nayar and U. Ortega, “Computation of selected eigenvalues of generalized eigenvalue problems,” J. Comput. Phys. 108, 8–14 (1993).
  * Lehoucq and Sorensen (1996) R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly restarted Arnoldi iteration,” SIAM J. Matrix Anal. Appl. 17, 789–821 (1996).
  * Sipp and Marquet (2013) D. Sipp and O. Marquet, “Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer,” Theor. Comput. Fluid Dyn. 27, 617–635 (2013).
  * Amestoy  _et al._ (2001) P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster, “A fully asynchronous multifrontal solver using distributed dynamic shedding,” SIAM J. Matrix Anal. Appl. 23, 15–41 (2001).
  * Meliga  _et al._ (2016) P. Meliga, E. Boujo, and F. Gallaire, “A self-consistent formulation for the sensitivity analysis of finite-amplitude vortex shedding in the cylinder wake,” J. Fluid Mech. 800, 327–357 (2016).
  * Zielinska  _et al._ (1997) B. J. A. Zielinska, S. Goujon-Durand, J. Dušek, and J. E. Wesfried, “Strongly nonlinear effect in unstable wakes,” Phys. Rev. Lett. 79, 3893–3896 (1997).
  * Abdessemed  _et al._ (2009) N. Abdessemed, A. S. Sharma, S. J. Sherwin, and V. Theofilis, “Transient growth analysis of the flow past a circular cylinder,” Phys. Fluids 21 (2009).
  * Weideman and Reddy (2000) J. A. Weideman and S. C. Reddy, “A MATLAB differentiation matrix suite,” ACM Trans. Math. Softw. 26, 465–519 (2000).
  * Hack and Moin (2017) M. J. P. Hack and P. Moin, “Algebraic disturbance growth by interaction of Orr and lift-up mechanisms,” J. Fluid Mech. 829, 112–126 (2017).
  * Bourguignon (2012) J.-L. Bourguignon,  _Models of turbulent pipe flow_ , Ph.D. thesis, California Institute of Technology (2012).
  * McKeon  _et al._ (2013) B. J. McKeon, A. S. Sharma, and I. Jacobi, “Experimental manipulation of wall turbulence: A systems approach,” Phys. Fluids 25 (2013).
  * Barbagallo  _et al._ (2009) A. Barbagallo, D. Sipp, and P. J. Schmid, “Closed-loop control of an open cavity flow using reduced-order models,” J. Fluid Mech. 641, 1–50 (2009).
  * Chen and Rowley (2011) K. K. Chen and C. W. Rowley, “${H}_{2}$ optimal actuator and sensor placement in the linearised Ginzburg-Landau system,” J. Fluid Mech. 681, 241–260 (2011).
  * Luhar  _et al._ (2014) M. Luhar, A. S. Sharma, and B. J. McKeon, “Opposition control within the resolvent analysis framework,” J. Fluid Mech. 749, 597–626 (2014).
  * Moore (1981) B. Moore, “Principal component analysis in linear systems: Controllability, observability, and model reduction,” IEEE Trans. Autom. Control 26, 17–32 (1981).
  * Rowley (2005) C. W. Rowley, “Model reduction for fluids using balanced proper orthogonal decomposition,” Int. J. Bifurc. Chaos 15, 997–1013 (2005).