# Mean resolvent operator of a statistically steady flow

Colin Leclercq1<EMAIL_ADDRESS>Denis Sipp1  1ONERA, DAAA, Université
Paris Saclay, 8 rue des Vertugadins, 92190 Meudon, France

###### Abstract

This paper introduces a new operator relevant to input-output analysis of
flows in a statistically steady regime far from the steady base flow: the mean
resolvent $\mathsfbi{R}_{0}$. It is defined as the operator predicting, in the
frequency domain, the mean linear response to forcing of the time-varying base
flow. As such, it provides the statistically optimal linear time-invariant
approximation of the input-output dynamics, which may be useful, for instance,
in flow control applications. Theory is developed for the periodic case. The
poles of the operator are shown to correspond to the Floquet exponents of the
system, including purely imaginary poles at multiples of the fundamental
frequency. In general, evaluating mean transfer functions from data requires
averaging the response to many realizations of the same input. However, in the
specific case of harmonic forcings, we show that the mean transfer functions
may be identified without averaging: an observation referred to as ‘dynamic
linearity’ in the literature (Dahan et al., 2012). For incompressible flows in
the weakly unsteady limit, i.e. when amplification of perturbations by the
unsteady part of the periodic Jacobian is small compared to amplification by
the mean Jacobian, the mean resolvent $\mathsfbi{R}_{0}$ is well-approximated
by the well-known resolvent operator about the mean-flow. Although the theory
presented in this paper only extends to quasiperiodic flows, the definition of
$\mathsfbi{R}_{0}$ remains meaningful for flows with continuous or mixed
spectra, including turbulent flows. Numerical evidence supports the close
connection between the two resolvent operators in quasiperiodic, chaotic and
stochastic two-dimensional incompressible flows.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see Keyword PDF for the full list). Other
classifications will be added at the same time.

MSC Codes  (Optional) Please enter your MSC Codes here

##  1 Introduction

###  1.1 Existence of transfer functions for time-varying base flows ?

In the early stages of ‘natural’ (i.e. caused by small-amplitude two-
dimensional Tollmien–Schlichting waves) laminar-turbulent transition of an
incompressible flow on a smooth flat plate, the flow may be safely assumed to
respond to perturbations in a linear time-invariant (LTI) fashion. Otherwise,
the response of open shear flows to external excitations, even infinitesimal,
is generally more complex. Turbulent flows certainly are not time-invariant
indeed, but neither are idealistic two-dimensional flows like the ones behind
a backward-facing step or in the wake of one or a few cylinders at moderate
Reynolds numbers. Even at very low Reynolds number, unsteadiness may set in as
a result of noise amplification or intrinsic instability. Hence, infinitesimal
extrinsic perturbations do not interact with a steady base flow but a time-
evolving one.

However, for feedback control using modern robust control techniques, an LTI
model of the input-output dynamics is generally required. In the frequency
domain, LTI dynamics indicates the existence of transfer functions, hence the
question: how can we define transfer functions for time-varying, statistically
steady base flows ?

###  1.2 Option A: ‘dynamic linearity’, i.e. using small-amplitude harmonic
forcings

Some authors have used harmonic forcings to identify transfer functions on
simulations of statistically steady flows (Dahan et al., 2012; Dalla Longa et
al., 2017; Evstafyeva et al., 2017) which are either laminar or turbulent, and
possess either peaked or broadband spectra. For sufficiently low forcing
amplitude, the authors evaluate a frequency response by taking the ratio of
the complex amplitudes of the output and input at the forcing frequency (the
concept of nonlinear transfer function (Noiray et al., 2008) is not relevant
here due to the choice of low amplitude of the forcing). This observation
hints at the existence of a transfer function in situations where it should
not theoretically exist, given the time-varying nature of the base flow. The
authors use the term ‘dynamic linearity’ to characterize this interesting
behaviour but no theoretical arguments seem currently available to understand
it.

###  1.3 Option B: Linearizing about the mean flow

A second option consists in assuming that the linear response $\mathbf{u}$ to
infinitesimal forcing $\mathbf{f}$ is characterized by the Jacobian
matrix111With the exception of appendix C, we adopt in this paper a spatially-
discrete framework where the incompressible velocity field is discretized over
$N$ degrees of freedom. $\mathsfbi{J}_{\overline{\mathbf{U}}}$ evaluated about
the mean flow $\overline{\mathbf{U}}$, i.e.

$\mathrm{d}_{t}\mathbf{u}=\mathsfbi{J}_{\overline{\mathbf{U}}}\mathbf{u}+\mathbf{f}.$
(1)

The resolvent operator about the mean flow

$\mathsfbi{R}_{\overline{\mathbf{U}}}(s)=(s\mathsfbi{I}-\mathsfbi{J}_{\overline{\mathbf{U}}})^{-1}$
(2)

then describes an LTI input-output dynamics in the frequency domain ($s$ being
the Laplace variable), and has been successfully used for open-loop (Moarref &
Jovanović, 2012; Luhar et al., 2014; Toedtli et al., 2019; Yeh & Taira, 2019;
Liu et al., 2021) and closed-loop (Leclercq et al., 2019) control of
statistically steady flows. However, such an input-output model remains only
empirically validated so far, and there is no clear justification for its use,
since the mean flow is not an invariant solution of the nonlinear unforced
system. The use of such model is however motivated by a large body of
literature devoted to linear analysis about a time-averaged mean flow, with
clear predictive power.

Indeed, modal analysis of $\mathsfbi{J}_{\overline{\mathbf{U}}}$ yields
accurate predictions of both the dominant oscillation frequency in the wake of
a cylinder (Hammond & Redekopp, 1997; Pier, 2002; Barkley, 2006; Sipp &
Lebedev, 2007; Mittal, 2008) and the associated coherent vortical structures.
The growth rate of the leading eigenvalue, whose real part predicts the
frequency, is nearly neutral, a property often called RZIF (Turton et al.,
2015; Bengana & Tuckerman, 2019; Bengana et al., 2019; Bengana & Tuckerman,
2021) for Real Zero Imaginary Frequency, in a context that is now broader than
just cylinder wake flow. This powerful property has been used to build self-
consistent models (Mantič-Lugo et al., 2014, 2015; Meliga, 2017; Bengana &
Tuckerman, 2021) successfully mimicking the nonlinear flow with the simplest
ingredients. The idea of marginal stability of the mean flow is not recent
(Malkus, 1956) but may only be justified under specific conditions, i.e. weak
nonlinearity (Noack et al., 2003; Sipp & Lebedev, 2007), monochromatic
oscillations (Mezić (2013), criterion (b) of §3.5; Turton et al. (2015)) or
weak unsteadiness relative to the mean flow (Mezić (2013); criterion (a) of
§3.5).

In the case of linearly stable flows, or noise-amplifiers, the resolvent
operator about the turbulent mean flow has been extensively used to predict
second-order statistics (Jovanovic & Bamieh, 2001; Jovanovic & Georgiou, 2010;
Moarref et al., 2014; Beneddine et al., 2016; Zare et al., 2017; Towne et al.,
2020) and in particular coherent flow structures characterized by spectral POD
modes (Semeraro et al., 2016a, b; Beneddine et al., 2016; Schmidt et al.,
2018; Towne et al., 2018). The resolvent operator has also been used to obtain
deterministic reduced-order models of unsteady flows from measurements of the
mean flow and a few local probes (Gómez et al., 2016; Beneddine et al., 2017;
Symon et al., 2019). The input-output framework is also particularly fruitful
for elucidating receptivity mechanisms through the analysis of the optimal
forcing modes associated with the resolvent modes (Hwang & Cossu, 2010;
Garnaud et al., 2013; Sartor et al., 2015; Jeun et al., 2016; Semeraro et al.,
2016b).

The intuition that coherent structures must extract their energy from the mean
flow dates back to the contributions of Lee et al. (1990); Butler & Farrell
(1993); Del Álamo & Jimenez (2006); Cossu et al. (2009); Pujals et al. (2009);
McKeon & Sharma (2010) among others. However, despite the unchallenged
efficacy of resolvent and modal analyses about the mean flow for modelling
purposes, a clear justification for these procedures is still lacking
(Beneddine et al., 2016; Jovanović, 2021). As already said, the first
difficulty in terms of interpretation arises from linearizing about a
quantity, the mean flow, which is not an invariant solution of the governing
equations. In the input-output framework of McKeon & Sharma (2010), the
resolvent operator $\mathsfbi{R}_{\overline{\mathbf{U}}}$ is forced by a
finite-amplitude internal forcing arising from nonlinear interactions of the
perturbations. Although this framework is an exact reformulation of the
Navier–Stokes equations, with no approximation involved, the procedure remains
ad hoc in the sense that alternative formulations may be obtained by
decomposing the flow about any alternative reference state, for instance the
steady base flow. One practical reason for the popularity of this particular
formulation is the recurrent observation that for energetic frequencies of the
flow, the matrix $\mathsfbi{R}_{\overline{\mathbf{U}}}(\mathrm{i}\omega)$ is
often nearly rank 1 and the leading spectral POD mode is often well-
approximated by the leading resolvent mode (Beneddine et al., 2016; Towne et
al., 2018). The clear advantage in this situation is that coherent flow
structures may be predicted without having to consider the unknown nonlinear
forcing term.

However, it is also increasingly clear that incorporating information about
the nonlinear forcing into the linear operator, either in the form of an eddy
viscosity (Illingworth et al., 2018; Morra et al., 2019; Symon et al., 2019;
Madhusudanan et al., 2019; Pickering et al., 2021), or through a low-rank
state-feeback operator (Zare et al., 2017), may significantly enhance its
predictive power. For instance, the alignment between the leading resolvent
mode and leading spectral POD mode may increase from 0.4 to nearly 0.9 by
tuning the effective viscosity above the molecular viscosity (Pickering et
al., 2021). All these recent contributions confirm that
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ is a useful, yet suboptimal linear
input-output representation. Worse still, (Karban et al., 2020) recently
showed that resolvent analysis about a mean flow is in general ill-posed as
the poles of the operator depend on an arbitrary choice of formulation for the
governing equations. In the case of a supersonic heated jet, discrepancies of
up to $40\%$ were noted in the optimal gain at energetic frequencies,
depending on the choice between the use of conservative versus primitive
variables.

All these observations motivate the search for a new input-output operator i)
incorporating information about the nonlinear forcings, ii) which should be
optimal in some sense and iii) which definition should be intrinsic. As we
shall see, the new operator will help understand the ‘dynamic linearity’
phenomenon and provide a physical interpretation to
$\mathsfbi{R}_{\overline{\mathbf{U}}}$. Its poles will also exactly verify the
RZIF property in the periodic case, whereas
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ only approximately does so.

###  1.4 Option C: Averaging the time-varying linear response

In this paper, we will not make the ad hoc hypothesis of linearity about the
mean flow. We will only keep the assumption of linearity, which remains valid
as long as the forcing amplitude is low enough, but now the base flow is an
exact time-varying solution of the Navier–Stokes equations. Until §4, the base
flow will be assumed to be periodic with period $T$ corresponding to a
fundamental angular frequency $\omega_{0}:=2\pi/T$. We represent in figure 1,
the response $\mathbf{u}$ to an infinitesimal forcing signal $\mathbf{f}$
started at time $\tau_{0}$ relative to an arbitrary origin of time $\tau=0$ on
the base flow trajectory $\widetilde{\mathbf{U}}(\tau)$.

Figure 1: (a) Periodic base flow $\widetilde{\mathbf{U}}(\tau)$ for an
arbitrarily defined origin of time $\tau=0$. (b) Linear responses
$\mathbf{u}(t;\phi_{i})$ to the same infinitesimal forcing $\mathbf{f}(t)$
started at different instants $\tau=\tau^{i}_{0}$. Like the underlying base
flow $\mathbf{U}(t;\phi_{i})={\widetilde{\mathbf{U}}(\tau=t+\tau^{i}_{0})}$,
the responses are parametrized by the phase $\phi=\phi_{i}$ of the periodic
base flow at $t=0$, such that $\omega_{0}\tau_{0}=\phi\text{ mod }2\pi$. (c)
The response to linear forcing being a distribution, the statistically optimal
LTI model is obtained by computing the mean response
$\langle\mathbf{u}(t)\rangle_{\phi}$.

Because the base flow is not steady but periodic, the linear response to
forcing is not unique but parametrized by the phase $0\leq\phi<2\pi$ of the
limit cycle when the forcing starts. This key parameter for our analysis is
defined such that

$\omega_{0}\tau_{0}=\phi\text{ mod }2\pi.$ (3)

Indeed, the input-output dynamics is now characterized by a linear time-
periodic (LTP) system

$\mathrm{d}_{t}\mathbf{u}=\mathsfbi{J}(t;\phi)\mathbf{u}+\mathbf{f},\quad\mathsfbi{J}(t+T;\phi)=\mathsfbi{J}(t;\phi),$
(4)

with a periodic Jacobian itself parametrized by the phase $\phi$ through the
base flow $\mathbf{U}(t;\phi)={\widetilde{\mathbf{U}}(\tau=t+\tau_{0})}$. For
$\phi$ uniformly distributed in $[0,2\pi)$, there is in fact a distribution of
timeseries $\mathbf{u}(t;\phi)$ corresponding to a given input
$\mathbf{f}(t)$. The statistically optimal LTI approximation of this LTP
system should predict a single output $\mathbf{v}(t)$ minimizing the mean
error with respect to $\mathbf{u}(t;\phi)$, i.e.

$\mathbf{v}(t)=\operatorname*{argmin}_{\tilde{\mathbf{v}}(t)}\lim_{T^{\prime}\to\infty}\dfrac{1}{T^{\prime}}\int_{0}^{T^{\prime}}\langle\|\mathbf{u}(t;\phi)-\tilde{\mathbf{v}}(t)\|^{2}\rangle_{\phi}\,\mathrm{d}t,$
(5)

where $\langle.\rangle_{\phi}$ denotes an ensemble average with respect to
$\phi$. This output is obviously the mean response
$\mathbf{v}(t)=\langle\mathbf{u}(t)\rangle_{\phi}$. We define the mean
resolvent operator $\mathsfbi{R}_{0}$ (the notation will become clear in §C.2)
as the operator predicting the mean response, in the frequency domain, for any
given input

$\boxed{\text{mean resolvent
}\mathsfbi{R}_{0}:\mathbf{f}(s)\mapsto\langle\mathbf{u}(s)\rangle_{\phi}}.$
(6)

This transfer function is well-defined because the relationship between the
input $\mathbf{f}$ and the output $\langle\mathbf{u}\rangle_{\phi}$ is by
construction causal, linear and time-invariant (i.e. independent of $\phi$).

We see that in option C the order of operations is reversed compared to option
B: we linearize then average instead of averaging then linearizing. In option
A, linearization is involved as well, but there is no averaging at all.
Instead, a very specific form of input is chosen, which is harmonic in time.
The three options are summarized in figure 2.

Figure 2: Three ways of defining a LTI operator for a periodic base flow.

###  1.5 Goal and outline of the paper

The aim of this paper is to investigate option C and establish connections
with options A and B. In particular, we will show that option A is a specific
way to evaluate $\mathsfbi{R}_{0}$, while
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ is a (reduced-order) approximation of
$\mathsfbi{R}_{0}$, under appropriate conditions.

In §2 and §3, we will continue to focus on incompressible periodic base flows.
In §2 we will perform a numerical experiment on the fluidic pinball in a
periodic regime of oscillations at $Re=100$, in order to compare options B and
C, namely mean resolvent versus resolvent about the mean flow. In §3, we will
study the properties of the mean resolvent operator, with the goal to
elucidate the observations made in §2 and draw connections with the resolvent
operator about the mean flow and ‘dynamic linearity’. Connections with the
Koopman operator (§3.3), the RZIF property (§3.6) and model reduction (§3.5)
will also be made. In §4 we perform the same numerical experiment as in §2 but
for more complex base flows, i.e. quasiperiodic, chaotic and stochastic. We
will consider the fluidic pinball at $Re=110$ and $Re=120$ as well as the
backward-facing step flow (with exogenous stochastic forcing) at $Re=500$. The
goal is to stress the strong connection between mean resolvent and resolvent
about the mean flow in all these dynamically distinct cases. No computations
are done in a compressible case but in §4.3 we indicate the implications of
our theoretical analysis in §3 to such situations. A theoretical extension of
§3 is proposed for the quasiperiodic case in appendix D. We summarize our
findings in §5.

##  2 Resolvent about the mean flow versus mean resolvent: fluidic pinball at
$Re=100$

In this section, we perform a numerical experiment on the case of the
incompressible two-dimensional fluidic pinball in a periodic regime at
$Re=100$. The goal is to compare transfer functions based on options B and C,
i.e. average then linearize versus linearize then average. The precise
configuration consists in the the numerical setting of Deng et al. (2020). We
define the unit length to be the diameter $D$ of each cylinder and the
(convective) time unit to be $D/U_{\infty}$, where $U_{\infty}$ is the
upstream velocity norm. The Reynolds number is defined as
$Re:=U_{\infty}D/\nu$ where $\nu$ is the kinematic viscosity. The centers of
the ‘front’, ‘top’ and ‘bottom’ triangles are respectively located at
$(x_{f},y_{f})=(-1.5\cos(\pi/6),0)$, $(x_{t},y_{t})=(0,0.75)$ and
$(x_{b},y_{b})=(0,-0.75)$, hence forming an equilateral triangle of side
length $1.5$. Instead of measuring the full input-output relation between
$N$-dimensional inputs and outputs, we will focus in section §2 on the single-
input single-output (SISO) transfer between an actuator signal $u$ and a
measurement $y$ such that

$\mathbf{f}(t)=\mathbf{B}u(t),\quad
y(t;\phi)=\mathbf{C}^{T}\mathbf{u}(t;\phi).$ (7)

The SISO viewpoint is more practical for the numerical experiment we perform,
and it is also relevant in the context of flow control. More specifically, we
will consider three probes $\mathbf{C}$ extracting the $y$-velocity component
in the wake of the ‘top’ cylinder at $y=y_{t}$ and respectively $x=2,4,8$
(white filled circles in figure 3). Measurements of the periodic base flow
will be denoted $Y=\mathbf{C}^{T}\mathbf{U}$. When necessary, the index
‘a’,‘b’ or ‘c’ will be used to specify which sensor we are referring to. The
quantity $\mathbf{B}$ is a discretized Gaussian volume force field
$\mathcal{B}(x,y;x_{0},y_{0},\sigma_{x},\sigma_{y})$ located at the ‘top’ of
the ‘top’ cylinder

$\mathcal{B}(x,y;x_{0},y_{0},\sigma_{x},\sigma_{y})=\left[0,\exp\left({-\dfrac{(x-x_{0})^{2}}{2\sigma_{x}^{2}}-\dfrac{(y-y_{0})^{2}}{2\sigma_{y}^{2}}}\right)\right]^{T},$
(8)

with $\sigma_{x}=0.01$, $\sigma_{y}=0.1$ and the center $(x_{0},y_{0})$ of the
Gaussian (white triangle in figure 3) relative to the center $(x_{t},y_{t})$
of the cylinder being given by $(x_{0}-x_{t},y_{0}-y_{t})=(0,0.55)$.

Numerical discretization details (computational domain, mesh, spatial and
temporal schemes) are provided in §A.

###  2.1 Unsteady base flow

In figure 3(a), we plot the streamlines and isocontours of velocity norm for
the mean flow. We notice it is asymmetric, as noted by Deng et al. (2020). In
figure (b), we plot the phase portrait of the flow, using the vertical
velocity probes ‘b’ and ‘c’: it is a closed curve, which confirms the periodic
nature of the dynamics at $Re=100$. Fourier spectra are shown in figure 5(i)
for the 3 sensors (a,b,c) located in the wake of the cylinders. These spectra
confirm the periodic nature of the flow since they are discrete with a single
fundamental frequency at $\omega_{0}=0.73$ and harmonics. Higher frequencies
are more energetic in the far wake than in the near wake. Indeed,
perturbations amplify while being convected, causing stronger nonlinear energy
transfers downstream.

###  2.2 Estimating mean transfer functions from input-output data

The resolvent operator about the mean flow
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ and the mean resolvent operator
$\mathsfbi{R}_{0}$ are both multiple-input-multiple-output (MIMO) transfer
functions with $N$-dimensional inputs and outputs. SISO transfer functions
between any time-invariant actuator-sensor pair $(\mathbf{C},\mathbf{B})$ may
be derived from these operators using

$G_{\overline{\mathbf{U}}}(s):=\mathbf{C}^{T}\mathsfbi{R}_{\overline{\mathbf{U}}}(s)\mathbf{B},\quad\text{and}\quad\langle
G(s)\rangle_{\phi}:=\mathbf{C}^{T}\mathsfbi{R}_{0}(s)\mathbf{B}.$ (9)

In order to estimate the mean transfer function $\langle G(s)\rangle_{\phi}$
from input-output data, we introduce ‘frequency response realizations’

$G(s;\phi,u):=\dfrac{y(s;\phi)}{u(s)},$ (10)

which are obtained by taking the ratio between the Laplace transforms of the
input $u$ and the output $y$. The quantity $G$ is not directly associated with
a transfer function since the linear response to forcing about a time-varying
base flow is not time-invariant; hence the dependence on the phase $\phi$ and
the input signal $u$. However, upon averaging with respect to $\phi$ (we
recall that $\phi$ is assumed to be uniformly distributed in $[0,2\pi)$), the
linear input-output dynamics becomes time-invariant and the dependence
disappears. We then recover the mean transfer function based on the mean
resolvent

$\langle G(s)\rangle_{\phi}=\langle G(s;\phi,u)\rangle_{\phi}.$ (11)

For any complex frequency $s$ and input signal $u$, the variance
$\mathrm{Var}_{\phi}\,{G}(s;u)$ of the complex random variable $G(s;\phi,u)$
may be computed as

$\displaystyle\mathrm{Var}_{\phi}\,{G}$ $\displaystyle:=\langle|G-\langle
G\rangle_{\phi}|^{2}\rangle_{\phi},$
$\displaystyle=\langle|G|^{2}\rangle_{\phi}-|\langle G\rangle_{\phi}|^{2},$
(12)

and the ratio

$\eta(s;u):=\dfrac{\sqrt{\mathrm{Var}_{\phi}\,{G}}}{|\langle
G\rangle_{\phi}|}$ (13)

may be interpreted as the inverse of a ‘signal-to-noise’ ratio, allowing us to
quantify the validity of the time-invariance hypothesis. The system is nearly
LTI at a given frequency $s$ if $\eta\ll 1$ for any choice of $u$. A geometric
interpretation of $\eta$ will be proposed in figure 4 of §2.4.

The ratio $\eta$ is also useful to evaluate the convergence of the mean
estimate $\langle G\rangle_{\phi,N_{s}}$ to the true mean $\langle
G\rangle_{\phi}$, using $N_{s}\gg 1$ realizations of $G$. Since all
realizations are independent and identically distributed random variables of
variance $\mathrm{Var}_{\phi}\,G$, the mean estimate is, by the central limit
theorem, a normal random variable of variance

$\mathrm{Var}_{\phi}\,\langle
G\rangle_{\phi,N_{s}}=\dfrac{\mathrm{Var}_{\phi}\,G}{N_{s}}.$ (14)

The convergence criterion $\sqrt{\mathrm{Var}_{\phi}\,\langle
G\rangle_{\phi,N_{s}}}\ll|\langle G\rangle_{\phi}|$ for the mean estimate
therefore translates into the condition $\eta\ll\sqrt{N_{s}}$ for the given
choice of $u$ and $s$.

(a)  | (b)
---|---
Figure 3: Fluidic pinball at $Re=100$: periodic behaviour, (a) Mean velocity
norm and streamlines and (b) phase portraits using $y$-velocity sensors
$Y_{\mathrm{b}}$ and $Y_{\mathrm{c}}$.

###  2.3 Procedure

The nonlinear simulation is initialized at $\tau=0$ with the steady base flow
(fixed-point) perturbed by $\mathbf{B}$ (to quickly reach the limit-cycle).
After an initial transient of more than 1000 convective time-units, the flow
converges to a periodic regime of self-sustained oscillations.

This periodic base flow is then perturbed linearly using impulsive forcings
$u(t)=\delta(t)$ by simply initializing the perturbation field with the
actuator, i.e. $\mathbf{u}(0)=\mathbf{B}$. Other choices of input signals may
be made for estimating mean transfer functions, which will be discussed later
in §3.2.2 and §3.2.3. Following impulsive forcing, both nonlinear (for the
periodic base-flow) and linear (for the perturbation) simulations are run for
over more than 1000 convective time units. We compute $N_{s}=140$ impulse
responses, by uniformly distributing the relative initial time $\tau_{0}$ over
a period $T=8.64$ of the underlying limit-cycle, i.e.
$\Delta\tau_{0}=T/N_{s}$. This amounts to uniformly distributing the initial
phase $\phi$ over $[0,2\pi)$, as required.

Note that since impulse responses never decay, we can only compute frequency
responses shifted in the right-half plane, i.e. obtained for
$s=\mathrm{i}\omega+\sigma$ with $\sigma>0$. This is indeed mandatory for
convergence of the Laplace transform. The value of $\sigma=0.01$ is chosen
such that the mean impulse response modulated by $\mathrm{e}^{-\sigma t}$
reaches negligible values towards the end of the time window of 1000
convective time units. The Fourier transform of the exponentially modulated
signal, estimated via the discrete Fourier transform, then yields the Laplace
transform of the impulse over the shifted imaginary axis. The frequency
response over $\mathrm{i}\mathbb{R}$ (minus the singularities) may then be
retrieved by analytic continuation, but this is not done here.

The present analysis may remind the reader of the work by Yeh & Taira (2019),
who performed resolvent analysis about the mean flow on a shifted imaginary
axis with $\sigma>0$.

###  2.4 Results

First of all, a geometric interpretation of the ratio $\eta$ is proposed in
figure 4, where three cases are considered: (a) $\eta\ll 1$, (b)
$1\leq\eta<\sqrt{N_{s}}$ and (c) $\sqrt{N_{s}}\ll\eta$, corresponding to the
$\left\\{\text{probe,frequency}\right\\}$ pairs of table 1. For a given value
of $s$, the locus of $G(s;\phi)$ is represented by a closed blue curve in the
complex plane as $\phi$ varies between 0 and $2\pi$. The red dot indicates the
barycenter of the curve, which is the estimated mean transfer $\langle
G\rangle_{\phi,N_{s}}$. The length of the red arrow approximates $|\langle
G\rangle_{\phi}|$. There are two dashed circles centered on the red dot: a
black one with radius approximating $\sqrt{\mathrm{Var}_{\phi}G}$ (length of
black arrow), and a green one with radius equal to
$\sqrt{\mathrm{Var}_{\phi}\langle G\rangle_{\phi,N_{s}}}$ (length of green
arrow). In panel (a), the radius of the black circle is much smaller compared
to the red arrow and the entire blue distribution may therefore be
approximated by its red barycenter: the dynamics is quasi-time-invariant for
that $\left\\{\text{probe,frequency}\right\\}$ pair. In panel (b), the black
arrow is twice the size of the red arrow, hence replacing the entire blue
distribution by the single red dot is a very crude approximation: the dynamics
is not quasi-time-invariant. However, the green arrow is still much smaller
than the red arrow (see inset), which means that the estimation of the mean
transfer function is converged. But in panel (c), not only is the black arrow
much larger than the red arrow (dynamics not time-invariant), but the green
arrow is also much larger than the red arrow (see inset) so there are not
enough samples to confidently estimate the mean transfer function. We finally
recall that only the exact mean $\langle G\rangle_{\phi}$ is independent from
the forcing signal, while all other quantities (blue curve and its variance,
mean transfer estimate and its variance) depend on the specific choice of
$u(t)$.

probe | frequency $s$ | $\eta$
---|---|---
a | $0.01+8.00\mathrm{i}$ | $0.0554$
a | $0.01+0.813\mathrm{i}$ | $2.06$
c | $0.01+4.00\mathrm{i}$ | $73.6$
Table 1: Values of the ratio $\eta$ evaluated from $N_{s}=140$ impulses for
three $\left\\{\text{probe,frequency}\right\\}$ pairs. Figure 4: Geometric
interpretation of the ratio $\eta$ for the three cases (a) $\eta\ll 1$, (b)
$1<\eta<\sqrt{N_{s}}$ and (c) $\sqrt{N_{s}}\ll\eta$ of table 1; see main text
for description.

We then proceed to a more systematic examination of $\eta$ for
$u(t)=\delta(t)$ in row (ii) of figure 5, for the three probes (a,b,c) over
the frequency range $0\leq\omega\leq 10$. The threshold $\eta=1$ is indicated
with a dashed black line, while the ratio $\eta=\sqrt{N_{s}}$ is indicated
with a green dashed line. The value of $\eta$ with respect to these two
critical values is also indicated in the Bode diagrams of the mean transfer
function $\langle G\rangle_{\phi}$ represented in rows (iii) and (iv), for the
gain and phase respectively. No shading corresponds to a frequency range where
$\eta<1$, i.e. the dynamics is quasi-LTI, light-grey shading corresponds to
$1\leq\eta<\sqrt{Ns}$, i.e. dynamics not LTI but mean estimate converged, and
dark-grey shading corresponds to $\sqrt{N_{s}}\leq\eta$, i.e. mean estimate
not even converged.

Figure 5: Fluidic pinball in the periodic regime at $Re=100$. The labels
(a,b,c) correspond to the three sensor positions. (i) Fourier spectrum of
measurement $Y$ of the periodic base flow. The spectrum is discrete and no
averaging is done so a single Hann window is applied over the entire signal of
more than 1000 time units. (ii) Ratio $\eta$ (defined in (13)) for $N_{s}=140$
impulsive forcings $u(t)=\delta(t)$. The black dashed line indicates the
threshold $\eta=1$ far below which the system is nearly LTI (with respect to
$u=\delta$). The green dashed line indicates the threshold $\eta=\sqrt{N_{s}}$
far below which the estimate of the mean transfer function is converged. Gain
(iii) and phase (iv) of mean frequency response $\langle
G\rangle_{\phi,N_{s}}$ (solid blue) and frequency response about the mean flow
$G_{\overline{\mathbf{U}}}$ (dashed red). Light-grey shading indicated
frequency ranges where $1<\eta<\sqrt{Ns}$ and dark-grey shading corresponds to
$\eta\geq\sqrt{N_{s}}$. In the phase plots, the two black-dotted curves
indicate a shift of $+/-\pi$ with respect to the phase of
$G_{\overline{\mathbf{U}}}$. For rows (ii)-(iii)-(iv), all quantities are
evaluated on the shifted imaginary axis $\sigma+\mathrm{i}\omega$, with
$\sigma=0.01$.

The first remark is that the ratio $\eta$ becomes larger as the probe moves
downstream. As a result, converging the mean transfer estimate requires more
samples downstream than upstream (see greater extent of dark grey regions).
Correlatively, the frequency range of validity of the quasi-time-invariant
region shrinks as the probe moves downstream. Both observations seem
correlated with the fact that the unsteady fluctuations are more energetic
downstream than upstream (see power spectra in row (i)). The ratio $\eta$ is
also greater near resonance frequencies for the same reason: the base flow
fluctuations are greater near these frequencies. It is interesting to note
however that for all probe positions, the LTI approximation is valid in some
frequency range near the maximum gain (but excluding resonances).

Whether $\eta$ is small compared to 1 or not, we always observe a qualitative
agreement between the Bode diagrams of $G_{\overline{\mathbf{U}}}$ and
$\langle G\rangle_{\phi}$, as long as the mean estimate is converged. For
probe (a), the two diagrams are nearly undistinguisable from one another, and
the deviation remains very small for probe (b) as well, except near resonant
frequencies $k\omega_{0}$. For probe (c) though, the deviation becomes
noticeable, even when $\eta<1$.

Only the fundamental frequency leads to a resonance peak in
$G_{\overline{\mathbf{U}}}$, whereas higher-order harmonics are also visible
in $\langle G\rangle_{\phi}$. These extra peaks are indicative of the presence
of poles at $s=\mathrm{i}k\omega_{0}$ in the mean transfer function, but the
amplitude of the peaks rapidly decreases with $k$ and as the probe moves
upstream.

In the theoretical section, we will attempt to elucidate the aforementioned
observations: i) why are $\langle G\rangle_{\phi}$ and
$G_{\overline{\mathbf{U}}}$ generally so similar, ii) why does the quality of
the LTI approximation deteriorate downstream, iii) why are there poles at
$s=\mathrm{i}k\omega_{0}$ in the mean transfer function and iv) why are
higher-order resonances weaker and only noticeable downstream ?

##  3 Theory for incompressible periodic base flows

In this theoretical section, we come back to the more general MIMO viewpoint
between a $N$-dimensional forcing $\mathbf{f}(t)$ and a $N$-dimensional
response $\mathbf{u}(t;\phi)$.

###  3.1 Phase-dependent frequency response of LTP system

Denote $\mathbf{\Phi}(t,t^{\prime};\phi)$ the propagator from $t^{\prime}$ to
$t$ of the linear time-periodic system (4), such that

$\mathbf{u}(t;\phi)=\mathbf{\Phi}(t,0;\phi)\mathbf{u}_{0}+\int_{0}^{t}\mathbf{\Phi}(t,t^{\prime};\phi)\mathbf{f}(t^{\prime})\,\mathrm{d}t^{\prime},$
(15)

where $\mathbf{u}_{0}$ is the initial condition on $\mathbf{u}$, regardless of
$\phi$. By Floquet’s theorem (Magruder et al., 2018), we know that there
exists a real $T$-periodic matrix $\mathsfbi{P}(t;\phi)$, invertible at all
times $t$, and a constant complex matrix $\mathsfbi{K}$ such that
$\mathbf{\Phi}(t,t^{\prime};\phi)=\mathsfbi{P}(t;\phi)\mathrm{e}^{\mathsfbi{K}(t-t^{\prime})}\mathsfbi{P}^{-1}(t^{\prime};\phi)$.
There are multiple possible determinations of the matrix
$\mathsfbi{K}=(1/T)\mathrm{log}(\mathbf{\Phi}(T,0))$ depending on the
definition of the logarithm. Assuming $\mathsfbi{K}$ to be diagonalizable
$\mathsfbi{K}=\mathsfbi{S}\mathbf{\Lambda}\mathsfbi{S}^{-1}$ and choosing the
principal determination for the logarithm, all the eigenvalues fall into the
fundamental strip
$-\omega_{0}/2<\mathrm{Im}(\lambda_{i})\leq\omega_{0}/2,~{}i=1,\dots,N$. These
are called the principal Floquet exponents and they may be ranked in
decreasing order of growth rate, i.e.
$\mathrm{Re}(\lambda_{1})\geq\mathrm{Re}(\lambda_{2})\geq\dots\geq\mathrm{Re}(\lambda_{N})$.
Other determinations of the logarithm lead to eigenvalues
$s_{ij}=\lambda_{i}+\mathrm{i}j\omega_{0}$ which are Floquet exponents as
well, located in complementary strips
$\omega_{0}(j-1/2)<\mathrm{Im}(s_{ij})\leq\omega_{0}(j+1/2)$ in the complex
plane. The direct Floquet modes $\mathbf{v}^{i}$ are the columns of the
$T$-periodic complex matrix
$\mathsfbi{V}(t;\phi):=\mathsfbi{P}(t;\phi)\mathsfbi{S}$. The adjoint Floquet
modes $\mathbf{w}^{i}$, such that
$(\mathbf{w}^{i},\mathbf{v}^{j})=\delta_{ij}$ for the canonical inner product
$(\mathbf{a},\mathbf{b})=1/T\int_{0}^{T}\mathbf{a}^{H}\mathbf{b}\,\mathrm{d}t$
on the space of periodic functions in $\mathbb{C}^{N}$ (where H denotes the
complex transpose), are the columns of the matrix $\mathsfbi{W}(t;\phi)$ such
that $\mathsfbi{W}^{H}\mathsfbi{V}=\mathsfbi{I}$ at all times (and phases
$\phi$). Using this eigendecomposition, the propagator reads

$\mathbf{\Phi}(t,t^{\prime};\phi)=\mathsfbi{V}(t;\phi)\mathrm{e}^{\mathbf{\Lambda}(t-t^{\prime})}\mathsfbi{W}^{H}(t^{\prime};\phi).$
(16)

We may now expand the direct and adjoint Floquet modes as Fourier series. If
the $j\textsuperscript{th}$ Fourier coefficient of $\mathsfbi{V}$ is
proportional to $\mathrm{e}^{\mathrm{i}j\phi}$, then so is the
$j\textsuperscript{th}$ Fourier coefficient of $\mathsfbi{W}$ since
$\mathsfbi{W}^{H}\mathsfbi{V}=\mathsfbi{I}$ at all times. Therefore, we can
write the following Fourier series

$\mathsfbi{V}(t;\phi)=\sum_{j}\widehat{\mathsfbi{V}}_{j}\mathrm{e}^{\mathrm{i}j(\omega_{0}t+\phi)},\quad\mathsfbi{W}(t^{\prime};\phi)=\sum_{j}\widehat{\mathsfbi{W}}_{j}\mathrm{e}^{\mathrm{i}j(\omega_{0}t^{\prime}+\phi)}.$
(17)

Injecting (16) and (17) in (15) and assuming $\mathbf{u}_{0}=\boldsymbol{0}$,
we have

$\displaystyle\mathbf{u}(t;\phi)=$
$\displaystyle\int_{0}^{t}\mathsfbi{V}(t;\phi)\mathrm{e}^{\mathbf{\Lambda}(t-t^{\prime})}\mathsfbi{W}^{H}(t^{\prime};\phi)\mathbf{f}(t^{\prime})\,\mathrm{d}t^{\prime},$
$\displaystyle=$
$\displaystyle\int_{0}^{t}\left[\sum_{j}\widehat{\mathsfbi{V}}_{j}\mathrm{e}^{\mathrm{i}j(\omega_{0}t+\phi)}\right]\mathrm{e}^{\boldsymbol{\Lambda}(t-t^{\prime})}\left[\sum_{l}\widehat{\mathsfbi{W}}^{H}_{l}\mathrm{e}^{-\mathrm{i}l(\omega_{0}t^{\prime}+\phi)}\right]\mathbf{f}(t^{\prime})\,\mathrm{d}t^{\prime},$
$\displaystyle=$
$\displaystyle\sum_{j,l}\mathrm{e}^{\mathrm{i}(j-l)\phi}\int_{0}^{t}\left[\widehat{\mathsfbi{V}}_{j}\mathrm{e}^{(\boldsymbol{\Lambda}+\mathrm{i}j\omega_{0}\mathsfbi{I})(t-t^{\prime})}\widehat{\mathsfbi{W}}^{H}_{l}\right]\left[\mathrm{e}^{\mathrm{i}(j-l)\omega_{0}t^{\prime}}\mathbf{f}(t^{\prime})\right]\,\mathrm{d}t^{\prime},$
$\displaystyle=$
$\displaystyle\sum_{n}\mathrm{e}^{\mathrm{i}n\phi}\int_{0}^{t}\left[\sum_{j}\widehat{\mathsfbi{V}}_{j}\mathrm{e}^{(\boldsymbol{\Lambda}+\mathrm{i}j\omega_{0}\mathsfbi{I})(t-t^{\prime})}\widehat{\mathsfbi{W}}^{H}_{j-n}\right]\left[\mathrm{e}^{\mathrm{i}n\omega_{0}t^{\prime}}\mathbf{f}(t^{\prime})\right]\,\mathrm{d}t^{\prime}.$
(18)

We recognize convolution products of causal functions on the right-hand-side
hence it is easy to take the Laplace transform, yielding the linear time-
periodic input-output relation

$\boxed{\mathbf{u}(s;\phi)=\sum_{n}\mathrm{e}^{\mathrm{i}n\phi}\mathsfbi{R}_{n}(s)\mathbf{f}(s-\mathrm{i}n\omega_{0})},$
(19)

with the operators

$\boxed{\mathsfbi{R}_{n}(s)=\sum_{j}\sum_{k=1}^{N}\dfrac{\widehat{\mathbf{v}}^{k}_{j}{(\widehat{\mathbf{w}}^{k}_{j-n})^{H}}}{s-(\lambda_{k}+\mathrm{i}j\omega_{0})}},$
(20)

where $\widehat{\mathbf{v}}^{k}_{j}$ and $\widehat{\mathbf{w}}^{k}_{j}$
designate the $k\textsuperscript{th}$ column of $\widehat{\mathsfbi{V}}_{j}$
and $\widehat{\mathsfbi{W}}_{j}$. Unlike in LTI systems, the output at
frequency $s$ depends on the input at an infinite number of frequencies
$s-\mathrm{i}n\omega_{0}$. The output depends explicitly on the phase $\phi$
through the cross-frequency transfers
$\mathrm{e}^{\mathrm{i}n\phi}\mathsfbi{R}_{n}(s)$ for $n\neq 0$. The poles of
the $\mathsfbi{R}_{n}$ are exactly the Floquet exponents of the LTP system.
The adjoint Floquet modes characterize receptivity to forcing of the
corresponding direct Floquet modes.

###  3.2 Mean resolvent operator

We now seek the operator which predicts the mean output
$\langle\mathbf{u}(s)\rangle_{\phi}$ for a given input $\mathbf{f}(s)$, by
averaging (19) with respect to $\phi$. By doing so, all the cross-frequency
transfers vanish and the only remaining transfer is $\mathsfbi{R}_{0}$

$\boxed{\langle\mathbf{u}(s)\rangle_{\phi}=\mathsfbi{R}_{0}(s)\mathbf{f}(s)},$
(21)

which is consistent with our initial choice of notation for the mean resolvent
in the introduction.

####  3.2.1 Resonances at natural frequencies

Next, we note that since the periodic base flow $\mathbf{U}(t)$ is self-
sustained (no external forcing is required), there is a zero Floquet exponent
in the LTP system with associated Floquet mode equal to
$\mathrm{d}_{t}\mathbf{U}$ (see appendix B). If we restrict our attention to
the case where the periodic base flow is linearly stable, then the zero
Floquet exponent is the leading one, i.e. $\lambda_{1}=0$ and
$\mathbf{v}^{1}=\mathrm{d}_{t}\mathbf{U}$222In case of instability
$\lambda_{1}>0$ but there is still a zero exponent $\lambda_{j}=0$ for some
$j>1$., so that
$\widehat{\mathbf{v}}^{1}_{j}=\mathrm{i}j\omega_{0}\widehat{\mathbf{U}}_{j}$
where $\widehat{\mathbf{U}}_{j}$ denotes the $j^{\textsuperscript{th}}$
harmonic of the Fourier decomposition of $\mathbf{U}(t)$. Therefore, using
expression (20), the mean resolvent reads

$\boxed{\mathsfbi{R}_{0}(s)=\sum_{j\neq
0}\dfrac{\mathrm{i}j\omega_{0}{\widehat{\mathbf{U}}}_{j}(\widehat{\mathbf{w}}^{1}_{j})^{H}}{s-\mathrm{i}j\omega_{0}}+\sum_{j}\sum_{k=2}^{N}\dfrac{\widehat{\mathbf{v}}^{k}_{j}(\widehat{\mathbf{w}}^{k}_{j})^{H}}{s-(\lambda_{k}+\mathrm{i}j\omega_{0})}}.$
(22)

We immediately see that there are purely imaginary poles at
$s=\mathrm{i}j\omega_{0}$ as expected from the numerical experiment on the
fluidic pinball. The coefficients associated to these poles in the expansion,
also called the residuals

$\mathrm{Res}_{\mathrm{i}j\omega_{0}}\mathsfbi{R}_{0}(s)=\mathrm{i}j\omega_{0}{\widehat{\mathbf{U}}}_{j}(\widehat{\mathbf{w}}^{1}_{j})^{H},$
(23)

tend very rapidly to 0 in the operator norm as $|j|\to\infty$, because the
$j^{\textsuperscript{th}}$ Fourier harmonic of both $\mathbf{U}$ and
$\mathbf{w}^{1}$ are involved (the decay is $o(1/j^{m})$ for any integer $m$
if we consider $\mathcal{C}^{\infty}$ periodic structures). This explains why
resonance peaks caused by imaginary poles in figure 5(ii), become decreasingly
visible as $|j|$ increases. The fact that the residuals depend on
$\widehat{\mathbf{U}}_{j}$ also explains why resonances at high frequencies
are however more visible downstream than upstream. Indeed, larger oscillation
amplitude in the far wake causes more pronounced nonlinear energy transfers to
higher harmonics (see base flow spectra in row (i) of figure 5), hence a
better observability $|\mathbf{C}^{T}\widehat{\mathbf{U}}_{j}|$ of high-order
harmonics $\widehat{\mathbf{U}}_{j}$ as the sensor moves downstream.

####  3.2.2 System identification using frequency-rich inputs

The upside of using frequency-rich signals for system identification is that
it allows identifying the dynamics for multiple frequencies at once. The
downside though, is that many realizations of the same input signal may be
necessary to converge the mean response prior to identification, as we have
seen in §2 (in particular row (iv) in figure 5). The number of realizations
necessary to converge the mean is very dependent on the input signal chosen.
This dependence can be made explicit, using (19) and the definition of the
variance:

$\mathrm{Var}_{\phi}\mathbf{u}(s)=\sum_{n\neq
0}\|\mathbf{f}(s-\mathrm{i}n\omega_{0})\|_{\mathsfbi{R}_{n}^{H}\mathsfbi{R}_{n}}^{2},$
(24)

where
$\|\mathbf{a}\|_{\mathsfbi{M}}=(\mathbf{a}^{H}\mathsfbi{M}\mathbf{a})^{1/2}$
is the norm induced by the positive-definite matrix $\mathsfbi{M}$. A
convenient choice of input signal is one which minimizes the variance,
allowing for a minimal amount of realizations. Clearly, an input signal with a
broadband spectrum like an impulse $\mathbf{f}(t)=\mathbf{B}\delta(t)$ (for
which $\delta(s)=1$ for all $s$) is likely to generate a lot of output
variance and is not necessarily a favourable choice for identification. An
alternative to frequency-rich inputs is to use non-resonant harmonic forcings
for identification, as we shall see in the next section.

Another advantage of using a broadband input signal for identification is that
if the ratio $\eta$ is very small compared to 1 for such an input signal, then
the ratio is likely to be very small for any input signal, hence the LTI
approximation based on the mean response is likely to be meaningful.
Expression (24) also explains why the quality of the LTI approximation
deteriorates as the sensor moves downstream in figure 5(iv) for
$u(t)=\delta(t)$. Indeed,

$\mathrm{Var}_{\phi}G(s;u=\delta)=\sum_{n\neq
0}|\mathbf{C}^{T}\mathsfbi{R}_{n}(s)\mathbf{B}|^{2},$ (25)

which, according to (20), increases downstream since the observability
$|\mathbf{C}^{T}\widehat{\mathbf{v}}_{j}^{k}|$ of the Floquet modes probably
grows downstream (this is obvious for the first Floquet mode since
$\widehat{\mathbf{v}}^{1}_{j}=\mathrm{i}j\omega_{0}\widehat{\mathbf{U}}_{j}$).

####  3.2.3 System identification using harmonic inputs: ‘dynamic linearity’

The Laplace transform of harmonic forcings of the form
$\mathbf{f}(t)=\mathbf{B}\mathrm{e}^{\mathrm{i}\omega t}$ is given by

$\mathbf{f}(s)=\dfrac{\mathbf{B}}{s-\mathrm{i}\omega}.$ (26)

Plugging this expression into the input-output relation (19) leads to

$\mathbf{u}(s;\phi)=\sum_{n}\mathrm{e}^{\mathrm{i}n\phi}\dfrac{\mathsfbi{R}_{n}(s)\mathbf{B}}{s-\mathrm{i}(\omega+n\omega_{0})}.$
(27)

The temporal response may now be evaluated by inverting the Laplace transform.
For the linearly stable limit-cycle considered here, the decaying Floquet
exponents lead to a transient response while all the purely imaginary poles in
the right-hand-side of (27), either in the transfer operators or the forcing,
lead to a permanent contribution. For non-resonant forcing frequencies
$\omega\neq n\omega_{0}$, we have

$\displaystyle\mathbf{u}(t;\phi)$ $\displaystyle=\mathrm{transient}$
$\displaystyle+\mathsfbi{R}_{0}(\mathrm{i}\omega)\mathbf{B}\mathrm{e}^{\mathrm{i}\omega
t}+\sum_{j\neq
0}\mathrm{i}j\omega_{0}\widehat{\mathbf{U}}_{j}(\widehat{\mathbf{w}}_{j}^{1})^{H}\mathbf{B}\mathrm{e}^{\mathrm{i}j\omega_{0}t}$
$\displaystyle+\sum_{n\neq
0}\mathrm{e}^{\mathrm{i}n\phi}\left[\mathsfbi{R}_{n}(\mathrm{i}(\omega+n\omega_{0}))\mathbf{B}\mathrm{e}^{\mathrm{i}(\omega+n\omega_{0})t}+\sum_{j\neq
0}\mathrm{i}j\omega_{0}\widehat{\mathbf{U}}_{j}(\widehat{\mathbf{w}}_{j-n}^{1})^{H}\mathbf{B}\mathrm{e}^{\mathrm{i}j\omega_{0}t}\right].$
(28)

The first line of (28) collects the transient contributions, the second line
gathers terms which are phase-independent, while the third line is phase-
dependent and vanishes upon averaging. However, the point here is not to take
an average but notice that the Fourier coefficient at the forcing frequency,
which may be obtained using an harmonic average

$\widehat{\mathbf{u}}(\omega;\phi)=\lim_{T^{\prime}\to\infty}\dfrac{1}{T^{\prime}}\int_{0}^{T^{\prime}}\mathbf{u}(t;\phi)\mathrm{e}^{-\mathrm{i}\omega
t}\,\mathrm{d}t$ (29)

is phase-independent and exactly given by the product of the input amplitude
$\widehat{\mathbf{f}}(\omega)=\mathbf{B}$ by the mean resolvent operator
evaluated at $s=\mathrm{i}\omega$:

$\boxed{\widehat{\mathbf{u}}(\omega)=\mathsfbi{R}_{0}(\mathrm{i}\omega)\widehat{\mathbf{f}}(\omega)}.$
(30)

In other terms, using small-amplitude harmonic forcings, it is possible to
obtain frequency samples of the mean transfer function $\langle
G\rangle_{\phi}=\mathbf{C}^{T}\mathsfbi{R}_{0}\mathbf{B}$ without the need to
carry out an ensemble average over several realizations. This is in essence
the ‘dynamic linearity’ phenomenon described by (Dahan et al., 2012; Dalla
Longa et al., 2017; Evstafyeva et al., 2017) and used therein for LTI system
identification on time-varying base flows.

In principle, using a superposition of harmonic forcings

$\mathbf{f}(t)=\mathbf{B}(A_{1}\mathrm{e}^{\mathrm{i}\omega_{1}t}+A_{2}\mathrm{e}^{\mathrm{i}\omega_{2}t}+\dots),$
(31)

it should even be possible to sample the mean transfer function at all input
frequencies at once, using a single input realization, as long as
$\omega_{i}-\omega_{j}\neq n\omega_{0}$ and $\omega_{i}\neq n\omega_{0}$ for
any $i,j>0$.

###  3.3 Connection with the Koopman operator

The propagator from $t^{\prime}=0$ to $t$

$\mathbf{\Phi}(t,t^{\prime}=0;\phi)=\sum_{j}\sum_{k=1}^{N}\widehat{\mathbf{v}}^{k}_{j}\mathrm{e}^{(\lambda_{k}+\mathrm{i}j\omega_{0})t}\mathrm{e}^{\mathrm{i}j\phi}(\mathbf{w}^{k}(0;\phi))^{H}$
(32)

corresponds to the Koopman operator associated with the full-state observable
$\mathbf{u}$ in the case of linear dynamics about the time-periodic base flow
(Mezić & Surana, 2016). The $\widehat{\mathbf{v}}^{k}_{j}$ are the Koopman
modes, the Floquet exponents $\lambda_{k}+\mathrm{i}j\omega_{0}$ are the
Koopman eigenvalues, and
$\mathbf{u}_{0}\mapsto\mathrm{e}^{\mathrm{i}j\phi}(\mathbf{w}^{k}(0;\phi))^{H}\mathbf{u}_{0}$
are the phase-dependent Koopman eigenfunctions.

Averaging the propagator with respect to $\phi$ leads to the mean propagator

$\langle\mathbf{\Phi}(t,0)\rangle_{\phi}=\sum_{j}\sum_{k=1}^{N}\widehat{\mathbf{v}}^{k}_{j}\mathrm{e}^{(\lambda_{k}+\mathrm{i}j\omega_{0})t}{(\widehat{\mathbf{w}}^{k}_{j})^{H}}$
(33)

based on the phase-averaged Koopman eigenfunctions
$\mathbf{u}_{0}\mapsto\langle\mathrm{e}^{\mathrm{i}j\phi}(\mathbf{w}^{k}(0;\phi))^{H}\mathbf{u}_{0}\rangle_{\phi}={(\widehat{\mathbf{w}}^{k}_{j})}^{H}\mathbf{u}_{0}$.

By comparing equations (33) and (20) for $n=0$, we see that the mean resolvent
operator is the frequency domain representation of the mean propagator from
$0$ to $t$, i.e. of the mean Koopman operator for the full-state observable of
the LTP system

$\boxed{\mathsfbi{R}_{0}(s)=\mathcal{L}[\langle\mathbf{\Phi}(t,0)\rangle_{\phi}]},$
(34)

where $\mathcal{L}[.]$ denotes the Laplace transform.

###  3.4 Connection with LTI dynamics about the mean flow

####  3.4.1 Connection with the harmonic transfer operator

Evaluating (19) at various output frequencies $s+\mathrm{i}k\omega_{0}$, it is
possible to introduce the harmonic transfer operator
$\underline{\mathsfbi{H}}(s;\phi)$ such that

$\displaystyle\underbrace{\begin{pmatrix}\vdots\\\
\mathbf{u}(s-\mathrm{i}\omega_{0})\\\ \mathbf{u}(s)\\\
\mathbf{u}(s+\mathrm{i}\omega_{0})\\\
\vdots\end{pmatrix}}_{\underline{\mathbf{u}}(s;\phi)}=$
$\displaystyle\underbrace{\begin{pmatrix}\ddots&\ddots&\ddots&\ddots&\ddots\\\
\ddots&\mathsfbi{R}_{0}(s-\mathrm{i}\omega_{0})&\mathsfbi{R}_{-1}(s-\mathrm{i}\omega_{0})\mathrm{e}^{-\mathrm{i}\phi}&\mathsfbi{R}_{-2}(s-\mathrm{i}\omega_{0})\mathrm{e}^{-2\mathrm{i}\phi}&\ddots\\\
\ddots&\mathsfbi{R}_{1}(s)\mathrm{e}^{\mathrm{i}\phi}&\mathsfbi{R}_{0}(s)&\mathsfbi{R}_{-1}(s)\mathrm{e}^{-\mathrm{i}\phi}&\ddots\\\
\ddots&\mathsfbi{R}_{2}(s+\mathrm{i}\omega_{0})\mathrm{e}^{2\mathrm{i}\phi}&\mathsfbi{R}_{1}(s+\mathrm{i}\omega_{0})\mathrm{e}^{\mathrm{i}\phi}&\mathsfbi{R}_{0}(s+\mathrm{i}\omega_{0})&\ddots\\\
\ddots&\ddots&\ddots&\ddots&\ddots\end{pmatrix}}_{\underline{\mathsfbi{H}}(s;\phi)}\underbrace{\begin{pmatrix}\vdots\\\
\mathbf{f}(s-\mathrm{i}\omega_{0})\\\ \mathbf{f}(s)\\\
\mathbf{f}(s+\mathrm{i}\omega_{0})\\\
\vdots\end{pmatrix}}_{\underline{\mathbf{f}}(s)}.$ (35)

The block $\underline{\mathsfbi{H}}_{jk}$ of this infinite matrix operator
characterizes the transfer from the frequency $s+\mathrm{i}k\omega_{0}$ of the
input $\mathbf{f}$ to the frequency $s+\mathrm{i}j\omega_{0}$ of the output
$\mathbf{u}$ (Wereley & Hall, 1990, 1991; Zhou & Hagiwara, 2002; Zhou, 2008).
The mean resolvent operator appears along the diagonal of the harmonic
transfer operator, as it characterizes the phase-independent transfer from any
frequency of the input to the same frequency at the output. In particular, we
have

$\boxed{\mathsfbi{R}_{0}(s)=\underline{\mathsfbi{H}}_{00}(s)}.$ (36)

####  3.4.2 The harmonic transfer operator as a feedback loop

Decompose the Jacobian as a mean $\overline{\mathsfbi{J}}$ and a periodic
perturbation $\mathsfbi{J}^{\prime}(t;\phi)$, which we expand as a Fourier
series

$\mathsfbi{J}(t;\phi)=\overline{\mathsfbi{J}}+\underbrace{\sum_{j}\widehat{\mathsfbi{J}}^{\prime}_{j}\mathrm{e}^{\mathrm{i}j(\omega_{0}t+{\phi})}}_{\mathsfbi{J}^{\prime}(t;\phi)}.$
(37)

Because $\mathsfbi{J}^{\prime}(t;\phi)$ is real, the harmonics have Hermitian
symmetry, i.e.
$\widehat{\mathsfbi{J}}^{\prime}_{-j}=\widehat{\mathsfbi{J}}^{\prime*}_{j}$,
where $(.)^{*}$ denotes the complex conjugate (not the conjugate transpose
$(.)^{H}$). Moreover, we also obviously have
$\widehat{\mathsfbi{J}}^{\prime}_{0}=\mathsfbi{0}$. Since we are considering
the incompressible Navier–Stokes equations, the nonlinearity is quadratic
hence the mean Jacobian is equal to the Jacobian operator about the mean flow

$\overline{\mathsfbi{J}}=\mathsfbi{J}_{\overline{\mathbf{U}}}.$ (38)

Assuming $\mathbf{u}^{0}=0$, taking the Laplace transform of (4) and plugging
(37)-(38) yields

$s\mathbf{u}(s;\phi)=\mathsfbi{J}_{\overline{\mathbf{U}}}\mathbf{u}(s;\phi)+\sum_{j}\widehat{\mathsfbi{J}}^{\prime}_{j}\mathrm{e}^{\mathrm{i}j\phi}\mathbf{u}(s-\mathrm{i}j\omega_{0};\phi)+\mathbf{f}(s),$
(39)

which in turn leads, through harmonic balance (Khalil, 2002), to an
alternative form of the harmonic transfer operator, also referred to as the
harmonic resolvent operator (Padovan et al., 2020; Franceschini et al., 2022)

$\underline{\mathsfbi{H}}(s;\phi)=(\underline{\mathsfbi{D}}(s)-\underline{\mathsfbi{T}}(\phi))^{-1}.$
(40)

The infinite matrix $\underline{\mathsfbi{D}}(s)$ is block-diagonal, while
$\underline{\mathsfbi{T}}(\phi)$ is block-Laurent (Krishna Kumar & Kulkarni,
2015) with zero blocks on the diagonal:

$\displaystyle\underline{\mathsfbi{D}}(s)$
$\displaystyle=\begin{pmatrix}\ddots&\ddots&\ddots&\ddots&\ddots\\\
\ddots&(s-\mathrm{i}\omega_{0})\mathsfbi{I}-\mathsfbi{J}_{\overline{\mathbf{U}}}&\mathsfbi{0}&\mathsfbi{0}&\ddots\\\
\ddots&\mathsfbi{0}&s\mathsfbi{I}-\mathsfbi{J}_{\overline{\mathbf{U}}}&\mathsfbi{0}&\ddots\\\
\ddots&\mathsfbi{0}&\mathsfbi{0}&(s+\mathrm{i}\omega_{0})\mathsfbi{I}-\mathsfbi{J}_{\overline{\mathbf{U}}}&\ddots\\\
\ddots&\ddots&\ddots&\ddots&\ddots\end{pmatrix},$ (41)
$\displaystyle\underline{\mathsfbi{T}}(\phi)$
$\displaystyle=\begin{pmatrix}\ddots&\ddots&\ddots&\ddots&\ddots\\\
\ddots&\mathsfbi{0}&\widehat{\mathsfbi{J}}_{1}^{\prime*}\mathrm{e}^{-\mathrm{i}\phi}&\widehat{\mathsfbi{J}}_{2}^{\prime*}\mathrm{e}^{-2\mathrm{i}\phi}&\ddots\\\
\ddots&\widehat{\mathsfbi{J}}^{\prime}_{1}\mathrm{e}^{\mathrm{i}\phi}&\mathsfbi{0}&\widehat{\mathsfbi{J}}_{1}^{\prime*}\mathrm{e}^{-\mathrm{i}\phi}&\ddots\\\
\ddots&\widehat{\mathsfbi{J}}^{\prime}_{2}\mathrm{e}^{2\mathrm{i}\phi}&\widehat{\mathsfbi{J}}^{\prime}_{1}\mathrm{e}^{\mathrm{i}\phi}&\mathsfbi{0}&\ddots\\\
\ddots&\ddots&\ddots&\ddots&\ddots\end{pmatrix}.$ (42)

By recasting (40) in the form

$\underline{\mathsfbi{H}}(s;\phi)=(\underline{\mathsfbi{I}}-\underline{\mathsfbi{D}}^{-1}(s)\underline{\mathsfbi{T}}(\phi))^{-1}\underline{\mathsfbi{D}}^{-1}(s),$
(43)

the operator may be interpreted as a feedback loop between the two blocks
$\underline{\mathsfbi{D}}^{-1}(s)$ and $\underline{\mathsfbi{T}}(\phi)$, as
illustrated in figure 6.

Figure 6: Block representation of the harmonic transfer operator
$\underline{\mathsfbi{H}}(s;\phi)$ as a feedback loop between two blocks:
$\underline{\mathsfbi{D}}^{-1}(s)$ accounts for interactions with the mean
flow and $\underline{\mathsfbi{T}}(\phi)$ accounts for interactions with the
harmonics of the periodic Jacobian.

A parallel may be drawn with the describing function methodology (Gelb &
Vander Velde, 1968), which applies to oscillating systems composed of an LTI
block in feedback loop with a nonlinear time-invariant block, instead of a
linear time-varying one as here. The describing function methodology may be
used to obtain transfer functions parametrized by the forcing amplitude
(Noiray et al., 2008), whereas the present methodology yields transfer
functions parametrized by the phase $\phi$. The two approaches are not
mutually exclusive and one may extend the present framework by adding a
nonlinear time-invariant block as well to consider the effect of forcing
amplitude for harmonic inputs.

The operator $\underline{\mathsfbi{D}}^{-1}(s)$ is block-diagonal hence does
not transfer energy from one frequency to another. Physically, this block
represents interactions of the perturbation with the mean flow, and this
process is independent of the phase $\phi$. This block takes a forcing at the
input and delivers a velocity at the output. The feedback block does the
opposite and corresponds to interactions with the fluctuating part of the
periodic base flow. It contains only off-diagonal terms associated to the
various harmonics of the periodic perturbation flow $\mathbf{U}^{\prime}$,
meaning that it can only transfer energy from one frequency to another, and
this process is phase-dependent. This alternative interpretation of the
harmonic transfer operator is key to understanding the connection of the mean
resolvent operator to the resolvent operator about the mean flow, as we shall
see next.

####  3.4.3 The resolvent about the mean flow approximates the mean resolvent

The inverse
$(\underline{\mathsfbi{I}}-\underline{\mathsfbi{D}}^{-1}(s)\,\underline{\mathsfbi{T}}(\phi))^{-1}$
in (43) may be expanded as a Neumann series $\sum_{k\geq
0}(\underline{\mathsfbi{D}}^{-1}(s)\,\underline{\mathsfbi{T}}(\phi))^{k}$ such
that

$\underline{\mathsfbi{H}}(s;\phi)=\underbrace{\underline{\mathsfbi{D}}^{-1}(s)}_{\underline{\mathsfbi{H}}^{0}(s)}+\underbrace{\underline{\mathsfbi{D}}^{-1}(s)\underline{\mathsfbi{T}}(\phi)\underline{\mathsfbi{D}}^{-1}(s)}_{\underline{\mathsfbi{H}}^{1}(s;\phi)}+\underbrace{\underline{\mathsfbi{D}}^{-1}(s)\underline{\mathsfbi{T}}(\phi)\underline{\mathsfbi{D}}^{-1}(s)\underline{\mathsfbi{T}}(\phi)\underline{\mathsfbi{D}}^{-1}(s)}_{\underline{\mathsfbi{H}}^{2}(s;\phi)}+\dots.$
(44)

The series converges if
$\|\underline{\mathsfbi{D}}^{-1}(s)\,\underline{\mathsfbi{T}}(\phi)\|_{2}<1$333Here
$\|.\|_{2}$ denotes the spectral norm, i.e. maximum singular value, not the
$H_{2}$ norm.. Assume $\mathrm{Im}(s)=:\sigma>\sigma_{\mathrm{max}}$ where
$\sigma_{\mathrm{max}}\geq 0$ is the maximum growth rate among the poles of
$\underline{\mathsfbi{D}}^{-1}$ and $\underline{\mathsfbi{H}}$, then the
series converges if a) base flow unsteadiness is sufficiently weak or b)
$\sigma$ is sufficiently large (see appendix C.1). Each contribution
$\underline{\mathsfbi{H}}^{i}(s;\phi)$ corresponds to going $i$ times around
the loop, i.e. interacting $i$ times with the unsteady part of the base flow.

Expansion (44) on the harmonic transfer operator leads to a similar expansion
on the mean resolvent operator since the latter operator is a sub-block of the
former according to (36), i.e.

$\mathsfbi{R}_{0}=\underline{\mathsfbi{H}}^{0}_{00}+\underline{\mathsfbi{H}}^{1}_{00}+\underline{\mathsfbi{H}}^{2}_{00}+\dots.$
(45)

This expansion may be interpreted in the following way:
$\underline{\mathsfbi{H}}^{i}_{00}$ collects contributions of
$\underline{\mathsfbi{H}}$ which deliver an output at the same frequency than
the input, after interacting $i$ times with the fluctuating part of the base
flow. Using the general terms

$\displaystyle\underline{\mathsfbi{D}}_{jk}^{-1}(s)$
$\displaystyle:=\mathsfbi{R}_{\overline{\mathbf{U}}}(s+\mathrm{i}j\omega_{0})\delta_{jk},$
(46) $\displaystyle\underline{\mathsfbi{T}}_{jk}(\phi)$
$\displaystyle:=\widehat{\mathsfbi{J}}^{\prime}_{j-k}\mathrm{e}^{\mathrm{i}(j-k)\phi},$
(47)

we find

$\displaystyle\underline{\mathsfbi{H}}^{0}_{00}$
$\displaystyle=\mathsfbi{R}_{\overline{\mathbf{U}}},$ (48)
$\displaystyle\underline{\mathsfbi{H}}^{1}_{00}$ $\displaystyle=\mathsfbi{0},$
(49) $\displaystyle\underline{\mathsfbi{H}}^{2}_{00}$
$\displaystyle=\sum_{j}\mathsfbi{R}_{\overline{\mathbf{U}}}(s)\widehat{\mathsfbi{J}}^{\prime*}_{j}\mathsfbi{R}_{\overline{\mathbf{U}}}(s+\mathrm{i}j\omega_{0})\widehat{\mathsfbi{J}}^{\prime}_{j}\mathsfbi{R}_{\overline{\mathbf{U}}}(s),$
(50) $\displaystyle\underline{\mathsfbi{H}}^{3}_{00}$
$\displaystyle=\sum_{j,k}\mathsfbi{R}_{\overline{\mathbf{U}}}(s)\widehat{\mathsfbi{J}}^{\prime*}_{j}\mathsfbi{R}_{\overline{\mathbf{U}}}(s+\mathrm{i}j\omega_{0})\widehat{\mathsfbi{J}}^{\prime}_{j-k}\mathsfbi{R}_{\overline{\mathbf{U}}}(s+\mathrm{i}k\omega_{0})\widehat{\mathsfbi{J}}^{\prime}_{k}\mathsfbi{R}_{\overline{\mathbf{U}}}(s),$
(51) $\displaystyle\vdots$

The order 0 term is equal to the resolvent operator about the mean flow, which
corresponds to interacting 0 times with the fluctuating Jacobian. The order 1
term is exactly 0, because it is impossible to interact just once with the
fluctuating Jacobian and come out of the loop at the same frequency than the
input. The first correction to the resolvent operator about the mean flow then
arises at order 2, when successive interactions with
$\widehat{\mathsfbi{J}}^{\prime}_{k}$ and its complex conjugate
$\widehat{\mathsfbi{J}}^{\prime*}_{j}=\widehat{\mathsfbi{J}}^{\prime}_{-j}$
occur, allowing the output to be at the same frequency than the input.

Expansion (45) may be used to bound the absolute and relative differences
between $\mathsfbi{R}_{0}$ and $\mathsfbi{R}_{\overline{\mathbf{U}}}$.
Introduce the following system norm

$\|\mathsfbi{G}\|_{\infty,\sigma}:=\sup_{\omega\in\mathbb{R}}\|\mathsfbi{G}(\sigma+\mathrm{i}\omega)\|_{2},$
(52)

and, for $\sigma>\sigma_{\mathrm{max}}$, the small parameter

$\epsilon_{\sigma}:=\sum_{k}\|\mathsfbi{R}_{\overline{\mathbf{U}}}\widehat{\mathsfbi{J}}^{\prime}_{k}\|_{\infty,\sigma},$
(53)

characterizing the ($\sigma$-dependent) loop gain associated with the feedback
interconnection of $\underline{\mathsfbi{D}}^{-1}(s)$ and
$\underline{\mathsfbi{T}}(\phi)$, averaged with respect to $\phi$. Using (45)
and the definition (53) of $\epsilon_{\sigma}$, we may write

$\displaystyle\|\mathsfbi{R}_{0}-\mathsfbi{R}_{\overline{\mathbf{U}}}\|_{2}(s)$
$\displaystyle\leq\|\mathsfbi{R}_{\overline{\mathbf{U}}}\|_{2}(s)(\epsilon_{\sigma}^{2}+\epsilon_{\sigma}^{3}+\dots),$
$\displaystyle=\|\mathsfbi{R}_{\overline{\mathbf{U}}}\|_{2}(s)\dfrac{\epsilon_{\sigma}^{2}}{1-\epsilon_{\sigma}}.$
(54)

The geometric series $\sum_{k=0}^{\infty}\epsilon_{\sigma}^{k}$ converges to
$(1-\epsilon_{\sigma})^{-1}$ if and only if $\epsilon_{\sigma}<1$, which
occurs if either a) base flow unsteadiness is sufficiently weak or b) $\sigma$
is large enough (see appendix C.2). Since
$\|\mathsfbi{R}_{\overline{\mathbf{U}}}\|_{2}(s)\to 0$ as either
$\sigma\to\infty$ at fixed $\omega$ or $|\omega|\to\infty$ at fixed $\sigma$,
the absolute difference in (54) also tends to zero in these limits, i.e.
$\|\mathsfbi{R}_{0}-\mathsfbi{R}_{\overline{\mathbf{U}}}\|_{2}(s)\to 0$.
Moreover, the relative error at fixed $s$

$\boxed{\dfrac{\|\mathsfbi{R}_{0}-\mathsfbi{R}_{\overline{\mathbf{U}}}\|_{2}(s)}{\|\mathsfbi{R}_{\overline{\mathbf{U}}}\|_{2}(s)}=O(\epsilon_{\sigma}^{2})},$
(55)

is order 2, hence vanishes very quickly with the small parameter
$\epsilon_{\sigma}$. As already said, there are two independent ways to make
$\epsilon_{\sigma}$ small. Taking large $\sigma$ (option b) allows for low
relative difference between the two operators even in strongly unsteady base
flows. But in this case, using an LTI approximation of the input-output
dynamics becomes insufficient (see ratio $\eta$ introduced in §2.2).
Therefore, we argue that the key reason why
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ is a physically relevant operator in
incompressible flows, is because it approximates $\mathsfbi{R}_{0}$ when base
flow unsteadiness is weak (option a). The fact that the relative error is
order 2 with respect to base flow unsteadiness may explain the robust
agreement observed between the two operators. In the case of the fluidic
pinball (see Bode diagrams in figure5(iii,iv)), agreement was observed for a
low value of $\sigma=0.01$, indicating small amplification of the
perturbations by the unsteady part of the base flow.

It is interesting to note that the small correction to
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ involves the unsteady part of the
Jacobian operator, which indirectly incorporates information about the
endogenous nonlinear forcings of McKeon & Sharma (2010) into the linear
operator $\mathsfbi{R}_{0}$. Indeed, the fluctuating part of the Jacobian is
aware of the harmonic balance between the various Fourier components of the
nonlinear base flow. In this regard, the present work may be seen as an
extension of the recent contributions seeking to incorporate information about
nonlinear forcings into the linear operator through either a turbulent
viscosity (Morra et al., 2019; Pickering et al., 2021) or a state-feedback
operator (Zare et al., 2017).

###  3.5 Reduced-order models of the mean resolvent

Another way to write the mean resolvent is in the form

$\boxed{\mathsfbi{R}_{0}(s)=\underline{\mathsfbi{C}}(s\underline{\mathsfbi{I}}-\underline{\boldsymbol{\Lambda}})^{-1}\underline{\mathsfbi{B}}},$
(56)

where

$\displaystyle\underline{\mathsfbi{C}}$
$\displaystyle=[\dots,\widehat{\mathsfbi{V}}_{-1},\widehat{\mathsfbi{V}}_{0},\widehat{\mathsfbi{V}}_{1}\dots],$
(57) $\displaystyle\underline{\mathsfbi{B}}$
$\displaystyle=[\dots,\widehat{\mathsfbi{W}}_{-1},\widehat{\mathsfbi{W}}_{0},\widehat{\mathsfbi{W}}_{1}\dots]^{H},$
(58) $\displaystyle\underline{\mathsfbi{I}}_{jk}$
$\displaystyle=\mathsfbi{I}\delta_{jk},$ (59)
$\displaystyle\underline{\boldsymbol{\Lambda}}_{jk}$
$\displaystyle=(\boldsymbol{\Lambda}+\mathrm{i}j\omega_{0}\mathsfbi{I})\delta_{jk}.$
(60)

The input (resp. output) matrix $\underline{\mathsfbi{B}}$ (resp.
$\underline{\mathsfbi{C}}$) has $N$ columns (resp. rows) and an infinite
number of rows (resp. columns), while the diagonal state matrix
$\underline{\boldsymbol{\Lambda}}$ is infinite-dimensional. Expression (56) is
associated to a state-space representation

$\displaystyle\mathrm{d}_{t}\underline{\mathbf{x}}$
$\displaystyle=\underline{\boldsymbol{\Lambda}}\,\underline{\mathbf{x}}+\underline{\mathsfbi{B}}\,\mathbf{f},$
(61) $\displaystyle\langle\mathbf{u}\rangle_{\phi}$
$\displaystyle=\underline{\mathsfbi{C}}\,\underline{\mathbf{x}},$ (62)

where the internal state
$\underline{\mathbf{x}}=[\dots,\mathbf{x}_{-1},\mathbf{x}_{0},\mathbf{x}_{1},\dots]^{T}$
is an infinite column vector, even though the input $\mathbf{f}$ and the
output $\langle\mathbf{u}\rangle_{\phi}$ both have a finite dimension $N$.
However, for large enough $|j|$ the residuals of the poles
$\lambda_{k}+\mathrm{i}j\omega_{0}$ in $\underline{\boldsymbol{\Lambda}}$
become negligible, as they involve the Fourier coefficients
$\widehat{\mathsfbi{W}}_{j}$ in $\underline{\mathsfbi{B}}$ and
$\widehat{\mathsfbi{V}}_{j}$ in $\underline{\mathsfbi{C}}$ (for
$\mathcal{C}^{\infty}$ functions, the decay is $o(1/j^{m})$ for any positive
integer $m$). Hence the infinite-dimensional state vector
$\underline{\mathbf{x}}$ may be projected onto a finite-dimensional one by
simply eliminating high-frequency poles based on some cutoff value for the
residual norm.

The resolvent operator about the mean flow approximates $\mathsfbi{R}_{0}$ and
it only has $N$ poles: it may therefore be interpreted as a reduced-order
model of order $N$ of the mean resolvent. Even though the approximation is
good for weakly unsteady flows, the operator
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ does not take into account the
fluctuating part of the Jacobian, hence the model may not be optimal for the
fixed order $N$. Moreover, there is no reason a priori to choose a model order
equal to $N$: in model reduction, the order is fixed by the minimal number of
poles necessary to capture the input-output dynamics up to a given precision.
The appropriate order may therefore be chosen smaller or greater than $N$
depending on the desired precision.

###  3.6 The ‘RZIF’ property

We recall that $\lambda_{1}=0$ (see §3.2.1 and appendix B) therefore
$s=\mathrm{i}j\omega_{0}$ are Koopman eigenvalues of the mean propagator
satisfying exactly the so-called RZIF property originally discussed in the
context of mean flow stability analysis (Turton et al., 2015). Since
$\mathbf{v}^{1}=\mathrm{d}_{t}\mathbf{U}$, the associated Koopman modes
$\widehat{\mathbf{v}}_{j}^{1}=\mathrm{i}j\omega_{0}\widehat{\mathbf{U}}_{j}$
are parallel to the Fourier components $\widehat{\mathbf{U}}_{j}$ of the
periodic base flow. This is only approximately true in the case of
eigenvectors of the mean Jacobian associated with RZIF eigenvalues, unless the
oscillations of the periodic base flow are monochromatic (Mezić, 2013; Turton
et al., 2015), the dynamics is weakly non-linear (Noack et al., 2003; Sipp &
Lebedev, 2007) or weakly unsteady Mezić (2013). Therefore, the proposed
framework appears to be more generic than linear analysis about a mean flow,
as the poles of $\mathsfbi{R}_{0}$ and associated modes exactly satisfy the
RZIF property, while this is only approximately true for
$\mathsfbi{R}_{\overline{\mathbf{U}}}$.

Finally, we add that the RZIF property is also verified for unstable periodic
base flows, but in this case the zero Floquet exponent does not correspond to
the leading one, i.e. $\lambda_{1}\neq 0$, which has a strictly positive
growth rate (see appendix B).

##  4 Towards more complex base-flows

The previous sections were only concerned with periodic base flows, but the
definition (6) of the mean resolvent $\mathsfbi{R}_{0}$ may be easily
generalized to any other statistically steady base flow. Instead of taking the
mean output with respect to the phase $\phi$, we may take more generally the
mean output with respect to the relative time $\tau_{0}$ at which the forcing
signal is released:

$\boxed{\text{mean resolvent
}\mathsfbi{R}_{0}:\mathbf{f}(s)\mapsto\langle\mathbf{u}(s)\rangle_{\tau_{0}}}.$
(63)

For such a definition to hold, we assume that the single trajectory
$\widetilde{\mathbf{U}}(\tau)$ of the unsteady base flow, about which we
linearize the dynamics, covers the entire dynamical attractor. Otherwise, an
average over multiple base flow realizations also needs to be performed.

We start by performing similar numerical experiments as in §2, but for more
complex incompressible two-dimensional base flows. The numerical procedure
explained in §2.3 is repeated, but now varying the time $\tau_{0}$ at which
the impulse is triggered, rather than the phase $\phi$, which is no longer
defined. The definitions introduced in §2.2 are extended by replacing
$\langle.\rangle_{\phi}$ with $\langle.\rangle_{\tau_{0}}$, since now
‘frequency response realizations’ are parametrized by $\tau_{0}$.

In §4.1, we consider the fluidic pinball in the quasiperiodic and chaotic
regimes by increasing the Reynolds number to respectively $Re=110$ and
$Re=120$ (Deng et al., 2020). For the quasiperiodic case, an extension of the
periodic theory is possible and reported in appendix D. Next, in §4.2, we
consider the stochastic flow past a backward-facing step at $Re=500$, as in
Hervé et al. (2012). Two values are considered for the variance $\sigma_{w}$
of the stochastic forcing needed to sustain unsteadiness in the base flow.

By considering these extra cases, we cover all possible signs of the maximum
Lyapunov exponent (MLE), which is the maximal growth rate associated with the
linearized dynamics about a statistically steady base flow. The MLE may be
estimated from the mean impulse response according to

$\mathrm{MLE}=\lim_{t\to\infty}\dfrac{1}{t}\ln\dfrac{\|\langle\mathbf{u}(t)\rangle_{\tau_{0}}\|}{\|\langle\mathbf{u}(0)\rangle_{\tau_{0}}\|}.$
(64)

In the periodic case, it is equivalent to take the ensemble average with
respect to $\tau_{0}$ or $\phi$ and the MLE corresponds to
$\mathrm{Re}(\lambda_{1})=0$. In the quasiperiodic case, we also have
$\mathrm{MLE}=0$, for the chaotic regime $\mathrm{MLE}>0$, while for the
stochastic regime $\mathrm{MLE}<0$ (results of these computations are reported
in appendix E). The goal is to provide numerical evidence of the strong
connection between the mean resolvent and the resolvent about the mean flow in
all these cases, even though our theory only applies to the periodic and
quasi-periodic cases.

Finally, in §4.3, we discuss implications of our previous theoretical analysis
to the case of compressible flows.

###  4.1 Quasi-periodic and chaotic flows: fluidic pinball at $Re=110$ and
$Re=120$

The numerical setup is identical to that of §2 (with more details in appendix
A) for $Re=100$.

####  4.1.1 Unsteady base flow

In figure 7(a), we plot the streamlines and velocity norm isocontours for the
mean flow in the (i) quasiperiodic case and the (ii) chaotic case. We notice
that the mean flow is asymmetric in the quasi-periodic case, but symetric in
the chaotic case. The recirculation zone is more extended downstream in the
latter case.

In figure 7(b), we plot the phase portrait of the flow, using the vertical
velocity probes 2 and 3. The amplitude of the fluctuations increases with the
Reynolds number. The phase portrait is well-structured in the quasiperiodic
case but it is impossible to infer the dimensionality of the torus from this
two-dimensional projection of the attractor. The complex structure of the
phase portrait in figure 7(bii) is typical of chaotic dynamics.

Fourier spectra are shown in row (i) of figures 8 and 9 for each flow at each
probe location (a,b,c). These confirm the qualitatively different dynamics of
the flows: the spectrum is discrete in the quasiperiodic case, while it is
continuous in the chaotic case, despite the apparent presence of peaks. In the
latter case, they possess an intrinsic non-zero bandwidth, unlike in the
periodic and quasi-periodic cases, where the ‘thickness’ of the peak only
arises from the estimation error of the discrete Fourier transform. It is not
directly possible to infer the number of basic incommensurate frequencies in
the quasiperiodic case by examining the spectra, but there clearly is energy
in a much greater number of frequencies compared to the periodic case (see
figure 5(i) for comparison).

(ai)  | (bi)
---|---
(aii)  | (bii)
Figure 7: Fluidic pinball in different regimes: (i) quasiperiodic at $Re=110$
and (ii) chaotic at $Re=120$. (a) Mean velocity norm and streamlines and (b)
phase portraits using $y$-velocity sensors $Y_{\mathrm{b}}$ and
$Y_{\mathrm{c}}$. Figure 8: Fluidic pinball in the quasiperiodic regime at
$Re=110$. The labels (a,b,c) correspond to the three sensor positions. (i)
Fourier spectrum of measurement the unsteady base flow $Y$. The spectrum is
discrete and no averaging is done so a single Hann window is applied over the
entire signal of more than 1000 time units. (ii) Ratio $\eta$ (defined in
(13)) for $N_{s}=204$ impulsive forcings $u(t)=\delta(t)$. The black dashed
line indicates the threshold $\eta=1$ far below which the system is nearly LTI
(with respect to $u=\delta$). The green dashed line indicates the threshold
$\eta=\sqrt{N_{s}}$ far below which the estimate of the mean transfer function
is converged. Gain (iii) and phase (iv) of mean frequency response estimate
$\langle G\rangle_{\tau_{0},N_{s}}$ (solid blue) and frequency response about
the mean flow $G_{\overline{\mathbf{U}}}$ (dashed red). Light-grey shading
indicated frequency ranges where $1<\eta<\sqrt{N_{s}}$ and dark-grey shading
corresponds to $\eta\geq\sqrt{N_{s}}$. In the phase plots, the two black-
dotted curves indicate a shift of $+/-\pi$ with respect to the phase of
$G_{\overline{\mathbf{U}}}$. For rows (ii)-(iii)-(iv), all quantities are
evaluated on the shifted imaginary axis $\sigma+\mathrm{i}\omega$, with
$\sigma=0.01$. Figure 9: Same caption as figure 8 for the chaotic regime at
$Re=120$, with $N_{s}=204$ samples and $\sigma=0.03$.

####  4.1.2 Procedure

In both the quasi-periodic and chaotic cases, $N_{s}=204$ impulse responses
are computed over more than 1000 convective time units, each starting 2
convective time units after the preceeding one, i.e. $\Delta\tau_{0}=2$. We
assumed that this procedure allowed to sample a representative part of these
attractors with a single realization of $\widetilde{\mathbf{U}}(\tau)$. We
recall that $\sigma>\mathrm{MLE}$ is necessary for convergence of the Laplace
transform of the linear impulse response. In the quasi-periodic case, the
shifted frequency response is evaluated for $\sigma=0.01$ since
$\mathrm{MLE}=0$. For the chaotic case, we choose $\sigma=0.03$ since
$\mathrm{MLE}=0.02$ (see appendix E).

####  4.1.3 Results

Results are reported respectively in figures 8 and 9 for $Re=110$ and 120, and
they are qualitatively very similar to the periodic case in figure 5. The
ratio $\eta$ increases as the probe moves downstream, and as a consequence the
assumption of time-invariant input-output dynamics is quite poor for probes
(b) and (c), whereas it is quite good for probe (a). Unlike for $Re=100$, the
ratio $\eta$ is never smaller than 1 for probe (c) at $Re=110$ and 120, even
when the mean estimate is converged. The convergence of the mean estimate is
also slower downstream than upstream, and the number of samples $N_{s}=204$ is
insufficient at almost all frequencies for probe (c) using impulsive forcings
$u(t)=\delta(t)$.

Good overall agreement between the Bode diagrams of
$G_{\overline{\mathbf{U}}}$ and $\langle G\rangle_{\tau_{0}}$ is observed for
the two flow regimes. For probe (a), the Bode diagrams of the two transfer
functions are nearly undistinguishable from one another in both cases, except
in the vicinity of some energetic peaks in the power spectrum of the base flow
(see row (i)). For some high enough frequencies $\omega>4$, there are also
noticeable differences between the two transfer functions at probes (b) and
(c), even when $\eta<\sqrt{N_{s}}$. However, mean estimate convergence is only
ensured when $\eta\ll\sqrt{N_{s}}$ so it is difficult to conclude on these
deviations.

In the quasiperiodic case, we notice the presence of resonance peaks in the
mean transfer functions, at natural frequencies of the unsteady base flow,
which are not visible in $G_{\overline{\mathbf{U}}}$. Resonances at high
frequencies are not visible in the near wake and only manifest at sensors (b)
and (c) in the far wake. The observations made for the quasi-periodic case are
similar to the periodic case and the same justification may be provided in
both cases (see appendix D for an extension of the theory developed in §3 to
the quasiperiodic case). It is interesting to note that the close agreement
between $\langle G\rangle_{\tau_{0}}$ and $G_{\overline{\mathbf{U}}}$ carries
over to the chaotic regime, even though we do not have a theory for base flows
with a continuous power spectrum.

###  4.2 Stochastic flow: backward-facing step at $Re=500$

The numerical setup and code for the nonlinear flow are identical to that of
Hervé et al. (2012); Sipp & Schmid (2016). The mesh has $63,902$ triangles,
$32,792$ vertices and $N=258,970$ velocity degrees of freedom. The length unit
corresponds to the height $h$ of the step, while the (convective) time unit
corresponds to $h/U_{\infty}$. The time step is $\mathrm{d}t=0.002$. The main
difference with the fluidic pinball is that the steady base flow is linearly
stable hence requires constant forcing in order to reach a statistically
stationary unsteady state. The forcing is chosen to be of the form
$\mathbf{F}(t)=\mathbf{B}_{w}w(t)$ with $w$ a Gaussian-white noise and
$\mathbf{B}_{w}$ a discretized Gaussian volume force
$\mathcal{B}(x,y;x_{0},y_{0},\sigma_{x},\sigma_{y})$ with
$(x_{0},y_{0})=(-0.5,0.25)$ and $\sigma_{x}=\sigma_{y}=0.1$ (see equation (8)
for the definition of $\mathcal{B}$ and white squares centered at
$(x_{0},y_{0})$ in figures 10 (a,b)). Two noise variances are considered:
$\sigma_{w}=\sqrt{10}$ and $\sigma_{w}=10$, respectively representing weakly
and strongly nonlinear dynamics as in Hervé et al. (2012).

Linear perturbations about the time-varying base flow are solved in a similar
fashion to the case of the fluidic pinball (see appendix A). Note that the
forcing $\mathbf{F}$ produces the nonlinear base flow $\widetilde{\mathbf{U}}$
and should not be confused with $\mathbf{f}=\mathbf{B}u$, which produces the
linear response $\mathbf{u}$ about $\widetilde{\mathbf{U}}$. The actuator
$\mathbf{B}$ is a discretized version of $\mathcal{B}$ with
$(x_{0},y_{0})=(-0.05,0.01)$ and $(\sigma_{x},\sigma_{y})=(0.01,0.1)$ (see
white triangles in figures 10 (a,b)). We again use three probes to monitor the
flow response: two $y$-velocity probes in the shear layer at respectively
$(2.5,0)$ and $(5,0)$, and a friction sensor at $12.5\leq x\leq 12.9$ on the
lower wall, behind reattachment (see white dots and rectangle in figures
10(a,b)). The perturbation and base flow measurements will also be denoted
$y=\mathbf{C}^{T}\mathbf{u}$ and $Y=\mathbf{C}^{T}\mathbf{U}$ with index
‘a’,‘b’ or ‘c’ when necessary.

####  4.2.1 Unsteady base flow

The effect of the noise variance on the mean flow is quite visible in figure
10(a) for $\sigma_{w}=\sqrt{10}$ and 10(b) for $\sigma_{w}=10$. For stronger
perturbations, the recirculation region shrinks, with a reattachement point at
$x=10.7$ in the first case and $x=8.3$ in the second case. Timeseries of
sensor (a) are also shown in figure 10(c): the signals are clearly stochastic
and the fluctuation amplitude is consistent with $\sigma_{w}$. The Fourier
spectra are shown in row (i) of figures 11 and 12 for the two values of
$\sigma_{w}$ and the three probes (a,b,c). The spectra are broadband and there
is more perturbation energy at higher $\sigma_{w}$, as expected. Also,
perturbation energy increases further downstream, from sensor (a) to (c). At
sensor (a), the spectrum for $\sigma_{w}=10$ is roughly equal to that for
$\sigma_{w}=\sqrt{10}$, multiplied by a frequency-independent factor of
$\sqrt{10}$, which is indicative of quasi-linear behaviour. The signature of
strong nonlinear effects is visible downstream, as the spectra for both values
of $\sigma_{w}$ no longer have the same overall shape: saturation occurs
around $\omega\approx 1$, where linear amplification mechanisms dominate.

(a)
(b)
(c)

Figure 10: Isocontours of mean velocity norm and streamlines for the backward-
facing step flow at $Re=500$: (a) $\sigma_{w}=\sqrt{10}$, (b) $\sigma_{w}=10$,
where $\sigma_{w}$ is the rms of the amplitude $w(t)$ multiplying the external
noise field $\mathbf{B}_{w}$ (white square), such that
$\mathbf{F}(t)=\mathbf{B}_{w}w(t)$. The actuator is $\mathbf{B}$ (white
triangle), $y$-velocity probes are placed at locations ‘a’ and ‘b’ in the
shear layer (white circles), while a friction measurement (white rectangle) is
taken at ‘c’. (c) Timeseries of $Y_{\mathrm{a}}$ for $\sigma_{w}=\sqrt{10}$
(dashed black line) and $\sigma_{w}=10$ (blue solid line).

####  4.2.2 Procedure

Again, we use a single realization of $\widetilde{\mathbf{U}}(\tau)$
corresponding to a single random timeseries $w(\tau)$ but vary the relative
time $\tau_{0}$ of linear impulsive forcing to obtain various impulse
responses. The nonlinear flow is initialized with the steady base solution
(fixed point), and after a transient of more than 150 time units, it is
considered statistically stationary and linear impulses are then performed.
The time lag $\Delta\tau_{0}$ between two impulses is of 2 convective time
units, as in the quasiperiodic and chaotic flows. We run $N_{s}=308$
linearized impulse responses, over 220 time units only. Indeed, while impulse
responses never decay for the fluidic pinball, since $\mathrm{MLE}\geq 0$,
they do for the backward-facing step flow since $\text{MLE}<0$ (see appendix
E) and reach negligible amplitude over that period of time, making shorter
integration possible. Since $\mathrm{MLE}<0$, the Laplace transform may be
evaluated directly on the imaginary axis.

####  4.2.3 Results

Figure 11: Same caption as figure 5 for the backward-facing step flow at
$\sigma_{w}=\sqrt{10}$ (low-noise case), with $N_{s}=308$ samples. Unlike for
the fluidic pinball, the frequency response is on the imaginary axis, with no
shift since $\mathrm{MLE}<0$. The continuous power spectrum is obtained by
averaging over Hann windows of 86 convective times each, with $50\%$ overlap,
for a total signal length of 860 time units.

Results for the stochastic case are presented in figures 11 and 12,
respectively corresponding to low $\sigma_{w}=\sqrt{10}$ and high
$\sigma_{w}=10$ noise variance. Interestingly, the results are qualitatively
very similar to the fluidic pinball, even though the nature of the dynamics is
quite different. The main difference is that there are no longer sharp
variations in the gain and phase at very specific frequencies, in accordance
with the absence of resonance peak in the base flow spectrum (see row (i)).

The value of $\eta$ increases downstream, slowing down the convergence of the
mean estimate and progressively invalidating the time-invariant approximation.
This is correlated to spatial perturbation growth in the streamwise direction.
Similarly, increasing the noise variance $\sigma_{w}$ has a negative impact on
$\eta$, as unsteady perturbations become stronger. For probes (a) and (b),
there is a region where the time-invariant approximation is satisfactory, i.e.
$\eta<1$, which is located around the maximum gain, for both values of
$\sigma_{w}$. The frequency of maximum linear amplification of $u$ coincides
with that of maximum nonlinear amplification of the white noise $w$ (see base
flow power spectrum in row (i)).

Figure 12: Same caption as figure 11 for the high-noise case $\sigma_{w}=10$,
with $N_{s}=308$ samples.

###  4.3 Strongly compressible flows

It was recently pointed out by Karban et al. (2020) that resolvent analysis
about the mean flow is ill-posed in the case of strongly compressible flows,
where large discrepancies in the optimal gains are observed depending on the
formulation of the Navier–Stokes equations, whether in primitive or
conservative variables. The inconsistency was attributed to the fact that the
mean flow eigenvalues are not conserved under a nonlinear coordinate change.
Numerical experiments have not been carried out here for compressible cases,
but we wish to highlight the potential of the mean resolvent to resolve this
issue, since we showed that the poles of the operator correspond to Floquet
exponents (and Koopman eigenvalues) of the system in the LTP case (see §3.1
and §3.3), which are invariant under a nonlinear coordinate change (see
appendix F).

Another point we want to stress is the assumption (38) we made early on in
§3.4 that the mean Jacobian is equal to the Jacobian operator about the mean
flow. This is no longer true for cubic nonlinearities found in the
compressible Navier–Stokes equations written either in primitive or
conservative formulation. However, we may still use expansion (48)–(51) if we
replace $\mathsfbi{R}_{\overline{\mathbf{U}}}$ with

$\mathsfbi{R}_{\overline{\mathsfbi{J}}}(s)=(s\mathsfbi{I}-\overline{\mathsfbi{J}})^{-1},$
(65)

which we may call the resolvent operator about the mean Jacobian.
Alternatively, a quadratic formulation of the compressible Navier–Stokes
equations may be used (Vigo, 1998; Iollo et al., 2000), even though this form
may not easily handle shock discontinuities due to its non-conservative form.

##  5 Conclusions and outlook

This paper is concerned with the definition of a time-invariant operator best
characterizing linear input-output behaviour in statistically steady flows.
Rather than making the ad hoc assumption that the governing equations should
be linearized about a time-invariant mean flow, we force the time-varying
tangent system about unsteady trajectories on the attractor and collect the
responses to the same input signal for various realizations of the unsteady
base flow (this is done by varying the relative time $\tau_{0}$ at which the
input is triggered on a single base flow trajectory). We then consider the
ensemble average of responses produced by the given input to obtain a mean
transfer function. By considering two-dimensional incompressible
configurations, the fluidic pinball and the backward-facing step flow, we
investigated four possible dynamical regimes: periodic, quasi-periodic,
chaotic and stochastic. Inverting the order between the linearization and
averaging steps lead to interesting findings.

First of all, our framework allows to quantify the validity of the time-
invariant hypothesis for a given input signal and probe, frequency by
frequency. The ratio $\eta$ between the standard deviation and the module of
the mean quantifies the uncertainty associated with the mean frequency
response and may be used for robust controller design. In all flows
considered, the frequency response based on
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ appears to approximate the mean
frequency response very well, but the two objects are not identical. In
particular, for the periodic case, extra resonance peaks at exact multiples of
the fundamental frequency(ies) are visible in the mean transfer function.

We were able to explain these observations in the periodic and quasiperiodic
cases, by averaging the harmonic transfer operator (Wereley & Hall, 1990,
1991; Padovan et al., 2020; Franceschini et al., 2022) with respect to the
relative initial time of forcing $\tau_{0}$, or simply phase
$\omega_{0}\tau_{0}=\phi\,\text{mod}\,2\pi$ in the periodic case. Mean
transfer functions based on the resulting mean resolvent operator
$\mathsfbi{R}_{0}$, may be identified from input-output data without the need
for averaging over several input realizations if harmonic forcings are used,
even though the linearized system is not time-invariant. This justifies the
principle of ‘dynamic linearity’ put forward by Dahan et al. (2012); Dalla
Longa et al. (2017); Evstafyeva et al. (2017) in the context of flow control.

We then showed that the poles of $\mathsfbi{R}_{0}$ correspond to
(quasi-)Floquet exponents of the LT(Q)P system, which may also be interpreted
as Koopman eigenvalues. As a result, these poles are invariant under a
nonlinear coordinate change, i.e. do not depend on an arbitrary choice of
formulation for the Navier–Stokes equation. This strong property may help
solve the ambiguity of classical resolvent analysis about a mean flow recently
noted by Karban et al. (2020), which is particularly noticeable in strongly
compressible flows. There is also a set of marginally stable poles at
multiples of the fundamental frequency, which causes the resonance peaks in
the mean transfer functions. This so-called ‘RZIF’ property (Turton et al.,
2015), which is only approximate when considering the poles of
$\mathsfbi{R}_{\overline{\mathbf{U}}}$, appears to be exact when considering
the poles of $\mathsfbi{R}_{0}$.

Next, we further investigated the connection between
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ and $\mathsfbi{R}_{0}$ and found that
the former operator approximates the latter in the incompressible case (the
link is lost for cubic nonlinearities) within the weakly unsteady limit, where
amplification by the unsteady part of the base flow is small compared to
amplification due to the mean flow. The relative difference between the two
operators also vanishes for large positive growth rates (and the absolute
difference vanishes for large frequencies). There are however two key
differences between the two operators. One is that $\mathsfbi{R}_{0}$
indirectly incorporates information about the endogenous nonlinear forcing
term of McKeon & Sharma (2010), which are not taken into account by
$\mathsfbi{R}_{\overline{\mathbf{U}}}$, unless turbulent viscosity (Morra et
al., 2019; Pickering et al., 2021) or state-feedback (Zare et al., 2017) is
introduced. The missing information is embedded in the fluctuating part of the
Jacobian, which is aware of the nonlinear equilibrium between the various
Fourier components of the base flow. The other difference is that
$\mathsfbi{R}_{0}$ has an infinite-dimensional internal state, while
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ has a finite-dimensional one. However,
only a finite number of poles have a significant contribution to
$\mathsfbi{R}_{0}$, explaining why a finite-dimensional approximation of this
operator, in terms of internal state, is possible.

The definition of a new resolvent operator $\mathsfbi{R}_{0}$ is intended to
extend that of the usual resolvent operator
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ about the mean flow, therefore, any
analysis which may be done with $\mathsfbi{R}_{\overline{\mathbf{U}}}$
(reduced-order modelling, data-assimilation of second-order statistics, input-
output analysis, flow control, etc.) may, in principle, also be done with
(improved) $\mathsfbi{R}_{0}$, using yet to be defined numerical methods. For
instance, a compelling prospect would be to compute the singular modes of
$\mathsfbi{R}_{0}$ and compare them with the spectral POD modes to see if they
align better than the singular modes of
$\mathsfbi{R}_{\overline{\mathbf{U}}}$. It is hoped that the new operator may
help understand and overcome the observed shortcomings of the classical
resolvent operator, i.e. modelling of nonlinear forcings and compressibility
effects. However, it is important to note that the present paper does not
directly address turbulent nor compressible flows, and future work needs to be
done to confirm the relevance of $\mathsfbi{R}_{0}$ to these flows.

From a theoretical standpoint, an interesting perspective would be to extend
the present formalism to the case of continuous or mixed spectra
characteristic of chaotic, stochastic and turbulent flows. Koopman operator
theory may be key to this endeavour (Črnjarić-Žic et al., 2019). From a
methodological standpoint, a key question would be to find an optimal way to
compress the infinite state of the mean resolvent operator into a finite one.
Since the mean resolvent is related to the mean Koopman operator, this problem
may perhaps be tackled using extensions of the DMD method (Schmid, 2010;
Williams et al., 2015; Herrmann et al., 2021). Recent papers enhancing the
predictive power of $\mathsfbi{R}_{\overline{\mathbf{U}}}$ using optimization
techniques (Zare et al., 2017; Pickering et al., 2021) may possibly bear
connections with this compression problem.

[Declaration of interests]The authors report no conflict of interest.

[Author ORCIDs]Colin Leclercq, https://orcid.org/0000-0002-8262-0697; Denis
Sipp, https://orcid.org/0000-0002-2808-3886

##  Appendix A Fluidic pinball: discretization and numerical methods

As in Deng et al. (2020), the computational domain extends over $-6\leq x\leq
20$ and $-6\leq y\leq 6$, with the same Dirichlet boundary condition
$\widetilde{\mathbf{U}}=(1,0)$ at the inflow $x=-6$, ‘top’ $y=6$ and ‘bottom’
$y=-6$ boundary conditions. A standard outflow boundary condition
$\widetilde{P}\mathbf{n}-Re^{-1}\nabla\widetilde{\mathbf{U}}\cdot\mathbf{n}=\mathbf{0}$
is used at $x=20$, which is similar to the stress-free boundary condition of
Deng et al. (2020). The boundary condition on the cylinders is obviously no-
slip and impermeability, i.e. $\widetilde{\mathbf{U}}=\mathbf{0}$. The freely
available software FreeFem++ (Hecht, 2012) is used to time-march the
incompressible Navier–Stokes equations discretized on Taylor–Hood finite
elements; $P2$ for velocity components and forcings, and $P1$ for pressure.
The mesh comprises $66,668$ triangles, $33,739$ vertices and the total number
of velocity degrees of freedom, including both components, is $N=268,296$.

In order to validate spatial resolution, we produce 4 mesh $M1-M4$ of variable
refinement, and compute the lift and drag coefficients for the asymmetric
steady base flow at $Re=75$ and $Re=100$. The coefficients are based on the
total force by unit length exerted on all three cylinders normalized by
$1/2\rho U_{\infty}^{2}D$, where $\rho$ denotes density. The results are given
in table 2 and justify our choice of mesh M2 for the present study. Note that
the values provided here do not match with those of Deng et al. (2020) as the
authors did not include the viscous contribution to the forces (private
communication). For completeness, we provide in figure 13 the temporal
evolution of the lift coefficient in the periodic regimes at $Re=68,75,100$
analyzed by Deng et al. (2020), using the usual definition of the
coefficients. This may be useful for validation purposes in future studies on
this configuration.

The solvers are based on a sequential open-source code without time-splitting
($\mathrm{https://github.com/denissipp/AMR\\_Sipp\\_Schmid\\_2016}$; Sipp &
Schmid (2016)). Code parallelization was found unnecessary for this study,
which requires multiple realizations of impulsive forcings which can be run in
parallel indeed. The nonlinear problem is solved in perturbative form with
respect to the time-invariant base flow $\widetilde{\mathbf{U}}_{b}$ (i.e.
fixed-point), after computing the latter using a Newton method. The time-
stepping is semi-implicit: the diffusion term and interaction of the
perturbation
$\widetilde{\mathbf{U}}^{\prime\prime}:=\widetilde{\mathbf{U}}-\widetilde{\mathbf{U}}_{b}$
with $\widetilde{\mathbf{U}}_{b}$ are implicit, while the interaction of
$\widetilde{\mathbf{U}}^{\prime\prime}$ with itself is explicit
(Adams–Bashforth). The linear non-autonomous problem about the time-varying
base flow is also solved semi-implicitely: the diffusion term and interaction
of the perturbation $\mathbf{u}$ with $\widetilde{\mathbf{U}}_{b}$ are
implicit, while the interaction of $\mathbf{u}$ with
$\widetilde{\mathbf{U}}^{\prime\prime}$ is explicit (Adams–Bashforth as well).
The time-stepping scheme is second-order for both problems, however we found
the latter to be more sensitive to step size than the former. Temporal
resolution tests were carried out for linear impulses about the quasiperiodic
flow at $Re=110$. Figure 14 shows a portion of an impulse response for three
values of $\mathrm{d}t=0.01,0.05,0.0025$, indicating convergence for the lower
value of $\mathrm{d}t=0.0025$ on all three probes. Finally, the linearized
time-stepping code was validated by computing the impulse response about the
time-averaged mean flow. The corresponding ‘frequency response realization’
(shifted in the right-half plane, as specified in §2.3) was obtained by taking
the Laplace transform of the impulse response and then compared with
frequential samples computed using the resolvent operator about the mean flow
(on the same shifted axis).

| $M1$ | $M2$ | $M3$ | $M4$
---|---|---|---|---
Nb. of triangles | 19628 | 66668 | 78538 | 102965
Nb. of vertices | 10041 | 33739 | 39725 | 52005
$C_{L}$ asymmetric base flow $Re=75$ | 0.0383305 | 0.0381093 | 0.0381902 | 0.0381794
$C_{L}$ asymmetric base flow $Re=100$ | 0.0610811 | 0.0608825 | 0.0609652 | 0.0609484
$C_{D}$ asymmetric base flow $Re=75$ | 3.91464 | 3.91629 | 3.91598 | 3.91634
$C_{D}$ asymmetric base flow $Re=100$ | 3.66636 | 3.66776 | 3.66721 | 3.66758
Table 2: Resolution tests on lift and drag coefficients $C_{L}$ and $C_{D}$
for asymmetric base flow at $Re=75$ and $Re=100$. Figure 13: Time-delay
embedded plot of lift coefficient for the periodic regimes at (a) $Re=68$, (b)
$Re=75$ and (c) $Re=100$. Figure 14: Realization of linear impulse response
about quasiperiodic flow at $Re=110$ for $\mathrm{d}t=0.01$ (solid red),
$\mathrm{d}t=0.005$ (dashed blue) and $\mathrm{d}t=0.0025$ (dotted black), on
all 3 probes (a,b,c).

##  Appendix B Marginal Floquet exponents for self-sustained periodic base
flow

We find it useful to include a proof of this classical result. The periodic
base flow $\mathbf{U}$ is a solution of the unforced discretized
incompressible Navier–Stokes equations projected onto the space of solenoidal
velocity fields

$\mathrm{d}_{t}\mathbf{U}=\mathbf{F}(\mathbf{U}),$ (66)

and by differentiating with respect to time

$\mathrm{d}_{t}(\mathrm{d}_{t}\mathbf{U})=\underbrace{\mathrm{\nabla}_{\mathbf{U}(t)}\mathbf{F}}_{\mathsfbi{J}(t)}\mathrm{d}_{t}\mathbf{U}.$
(67)

Since $\mathrm{d}_{t}\mathbf{U}$ is a solution of the unforced LTP system, its
evolution from $t^{\prime}$ to $t$ is governed by the propagator

$\displaystyle\mathrm{d}_{t}\mathbf{U}(t)$
$\displaystyle=\mathbf{\Phi}(t,t^{\prime})\mathrm{d}_{t}\mathbf{U}(t^{\prime}),$
$\displaystyle=\mathsfbi{V}(t)\mathrm{e}^{\mathbf{\Lambda}(t-t^{\prime})}\mathsfbi{W}^{H}(t^{\prime})\mathrm{d}_{t}\mathbf{U}(t^{\prime}).$
(68)

For $t^{\prime}=t-T$, using the periodicity of $\mathrm{d}_{t}\mathbf{U}(t)$
and the Floquet modes, we have

$(\mathsfbi{I}-\mathrm{e}^{\mathbf{\Lambda}t})\mathsfbi{W}^{H}(t)\mathrm{d}_{t}\mathbf{U}(t)=0$
(69)

Clearly, $\mathrm{d}_{t}\mathbf{U}(t)$ is a Floquet mode associated with the
fundamental Floquet exponent $0$ (and copies $\mathrm{i}j\omega_{0}$ in
complementary strips). Since the periodic base flow is stable and Floquet
exponents are ranked in decreasing order of growth rate, we have
$\lambda_{1}=0$ and $\mathbf{v}^{1}=\mathrm{d}_{t}\mathbf{U}$. Note that if
the base flow is an unstable periodic orbit of the unforced system, then there
is still a zero Floquet exponent, but it will no longer be the leading one
since the maximum Lyapunov exponent has to be strictly positive.

##  Appendix C Convergence of operator expansions

In this appendix, we consider operators acting on continuous functions of
space and replace fonts accordingly; for instance, the Jacobian matrix
$\widehat{\mathsfbi{J}}_{i}^{\prime}$ is replaced by the operator
$\widehat{\mathcal{J}}_{i}^{\prime}$. For simplicity, we will consider the
specific case of spatially periodic $\mathcal{C}^{\infty}$ functions in an
infinite one-dimensional domain with period $L$, such that

$\displaystyle\mathcal{R}_{\overline{U}}(s)$
$\displaystyle=\left[s\mathcal{I}-(-\overline{U}\partial_{x}+Re^{-1}\partial^{2}_{xx})\right]^{-1},$
(70) $\displaystyle\widehat{\mathcal{J}}^{\prime}_{i}$
$\displaystyle=-\widehat{U}_{i}^{\prime}\partial_{x},$ (71)

where $\overline{U}$ and $\widehat{U}_{i}^{\prime}$ are scalar constants
(homogeneous base flow along $x$) and $\mathcal{I}$ is the identity. The
fundamental wavenumber $2\pi/L$ is denoted $k_{L}$. In this functional space,
Fourier modes
$\\{e_{n}:x\mapsto\mathrm{e}^{\mathrm{i}nk_{L}x},n\in\mathbb{Z}\\}$ form an
orthonormal basis for the canonical inner product
$(a,b)_{\mathrm{H}}:=\frac{1}{L}\int_{0}^{L}a^{*}(x)b(x)\,\mathrm{d}x$.

###  C.1 Harmonic transfer operator

The continuous-in-space-version of series (44) converges if and only if the
spectral radius $\rho(\underline{\mathcal{L}})$ of the operator
$\underline{\mathcal{L}}(s;\phi):=\underline{\mathcal{D}}^{-1}(s)\underline{\mathcal{T}}(\phi)$
is strictly less than one (Suzuki, 1976). The operator
$\underline{\mathcal{L}}$ is an infinite-matrix of operators
$\underline{\mathcal{L}}=(\underline{\mathcal{L}}_{jl})_{j,l\in\mathbb{Z}}$,
where
$\underline{\mathcal{L}}_{jl}=\mathcal{R}_{\overline{U}}(s+\mathrm{i}j\omega_{0})\widehat{\mathcal{J}}^{\prime}_{j-l}\mathrm{e}^{\mathrm{i}(j-l)\phi}$.
Introduce $\underline{e}^{n}$, the infinite-diagonal-matrix of operators such
that $\underline{e}_{j,l}^{n}=e_{n}\delta_{jl}$. The columns of this matrix
are orthonormal with respect to the inner product
$[\underline{a},\underline{b}]=\sum_{j}(\underline{a}_{j},\underline{b}_{j})_{\mathrm{H}}$.
The operator $\underline{\mathcal{L}}$ is said to be block-diagonal because it
leaves invariant the subspaces $V_{n}$ generated by the columns of
$\underline{e}^{n}$ (and these subspaces are in a direct sum). More
specifically, for every $n$, we may represent the operator
$\underline{\mathcal{L}}|_{V_{n}}$ as an infinite-matrix
$\underline{\mathsfbi{L}}^{n}$, i.e.
$\underline{\mathcal{L}}\underline{e}^{n}=\underline{e}^{n}\underline{\mathsfbi{L}}^{n}$,
where (for values of $s$ which do not cancel the denominator; it suffices that
$\sigma=\mathrm{Im}(s)>0$)

$\underline{\mathsfbi{L}}^{n}_{jl}=-\dfrac{\widehat{U}_{j-l}^{\prime}\mathrm{e}^{\mathrm{i}(j-l)\phi}\mathrm{i}nk_{L}}{s+\mathrm{i}j\omega_{0}+\overline{U}\mathrm{i}nk_{L}+Re^{-1}n^{2}k^{2}_{L}},$
(72)

therefore the operator norm

$\|\underline{\mathcal{L}}\|=\sup_{n}\|\underline{\mathsfbi{L}}^{n}\|$ (73)

(see p. 30 of Conway (1990)). Moreover
$\rho(\underline{\mathcal{L}})\leq\|\underline{\mathcal{L}}\|$ so it is
sufficient to show that $\sup_{n}\|\underline{\mathsfbi{L}}^{n}\|<1$ in
appropriate conditions. For any given $n$, we have
$\underline{\mathsfbi{L}}^{n}=\mathrm{i}nk_{L}\,\underline{\mathsfbi{r}}^{n}(s)\underline{\mathsfbi{u}}^{\prime}$,
where

$\displaystyle\underline{\mathsfbi{r}}^{n}_{jl}(s)$
$\displaystyle=(s+\mathrm{i}j\omega_{0}+\overline{U}\mathrm{i}nk_{L}+Re^{-1}n^{2}k^{2}_{L})^{-1}\delta_{jl},$
(74) $\displaystyle\underline{\mathsfbi{u}}^{\prime}_{jl}$
$\displaystyle=\widehat{U}^{\prime}_{j-l}\mathrm{e}^{\mathrm{i}(j-l)\phi},$
(75)

are respectively a diagonal matrix and a Laurent matrix of Fourier
coefficients of $U^{\prime}(t;\phi)$. The operator norm being
submultiplicative, we have
$\|\mathsfbi{L}^{n}\|\leq|n|k_{L}\|\underline{\mathsfbi{r}}^{n}\|\|\underline{\mathsfbi{u}}^{\prime}\|$.
For $\sigma>0$, we have
$\|\underline{\mathsfbi{r}}^{n}(s)\|\leq(\sigma+Re^{-1}n^{2}k_{L}^{2})^{-1}$.
Moreover,
$\|\underline{\mathsfbi{u}}^{\prime}\|\leq\|U^{\prime}\|_{\infty}=\max_{t}|U^{\prime}(t)|$
(see chapter III in Gohberg et al. (2003)), therefore

$\displaystyle\|{\underline{\mathcal{L}}(s)}\|$ $\displaystyle\leq
k_{L}\|U^{\prime}\|_{\infty}\max_{n}\dfrac{n}{\sigma+Re^{-1}n^{2}k_{L}^{2}},$
(76)
$\displaystyle\leq\dfrac{\|U^{\prime}\|_{\infty}}{2k_{L}}\sqrt{\dfrac{Re}{\sigma}},$
(77)

(evaluating the maximum over the set of real numbers). Interestingly, there
are two independent ways to make the norm of $\mathcal{L}(s)$ small: either a)
weak base flow unsteadiness $\|U^{\prime}\|_{\infty}$ or b) large $\sigma$.

To check the bound (77) numerically, we choose
$U^{\prime}(t;\phi):=2\sum_{j=1}^{\infty}\mathrm{e}^{-j}\cos[j(\omega_{0}t+\phi)]$
such that $\widehat{U}^{\prime}_{j}=\mathrm{e}^{-|j|}$ and
$\|U^{\prime}\|_{\infty}=2/(\mathrm{e}-1)\approx 1.16$. We also choose
$Re=k_{L}=\overline{U}=1$ and $\omega_{0}=\pi$. For each value of $n$, we
define the finite-dimensional approximation
$\mathsfbi{L}^{n,m}:=(\underline{\mathsfbi{L}}^{n})_{-m\leq i,j\leq m}$ of
$\underline{\mathsfbi{L}}^{n}$. We may then approximate the operator
$\underline{\mathcal{L}}$ by the finite block-diagonal matrix
$\mathsfbi{L}^{m}:=\mathrm{blkdiag}[\mathsfbi{L}^{-m,m},\dots,\mathsfbi{L}^{0,m},\dots,\mathsfbi{L}^{m,m}]$.
For such finite-dimensional matrices, we simply have
$\rho(\mathsfbi{L}^{m})=\max_{n}\rho(\mathsfbi{L}^{n,m})$ and
$\|\mathsfbi{L}^{m}\|_{2}=\max_{n}\|\mathsfbi{L}^{n,m}\|_{2}$. The results are
plotted in figure 15: in panel (a), we fix $s=0.1+\mathrm{1i}$ and vary $0\leq
m\leq 250$ while in panel (b) we fix $m=250$ and $\omega=1$ and vary
$-1\leq\log_{10}(\sigma)\leq 4$. Panel (a) confirms that $m=250$ is large
enough to converge the spectral radius and norm of
$\underline{\mathcal{L}}(s)$ for the $s$ considered (we verified that this the
case for all values of $\sigma$ in panel (b) as well). Panel (b) confirms
expression (77) for the bound and the decrease of the spectral norm as
$O(\sigma^{-1/2})$, which guarantees that the Neumann series (44) converges at
any $\omega$ for large enough $\sigma>0$.

(a) (b)

Figure 15: (a) Evolution of the spectral norm/radius of the finite-dimensional
approximation $\mathsfbi{L}^{m}$ of $\underline{\mathcal{L}}(s)$ with $m$ for
$s=0.1+\mathrm{i}$. (b) Evolution of the spectral norm/radius of
$\mathsfbi{L}^{250}$ for $s=\sigma+\mathrm{i}$ versus $\sigma$, with bound
from equation (77).

###  C.2 Mean resolvent operator

Since the mean resolvent is a sub-block of the harmonic transfer operator, the
convergence of the Neumann series (44) implies the convergence of expansion
(45) for the mean resolvent. Moreover, similarly to (77), it may be shown that
for any integer $j$

$\|\mathcal{R}_{\overline{U}}\mathcal{J}^{\prime}_{j}\|_{\infty,\sigma}\leq\dfrac{|\widehat{U}_{j}^{\prime}|}{2k_{L}}\sqrt{\dfrac{Re}{\sigma}},$
(78)

where $\|.\|_{\infty,\sigma}$ is the maximum value of the operator norm over
all real frequencies $\omega$, for a given value of $\sigma>0$. Since
$U^{\prime}(t)$ is $\mathcal{C}^{\infty}$,
$\widehat{U}^{\prime}_{j}=o(1/j^{m})$ for any integer $m$, and in particular
for $m\geq 2$, therefore the quantity $\sum_{j}|\widehat{U}^{\prime}_{j}|$ is
well-defined and so is
$\epsilon_{\sigma}:=\sum_{j}\|\mathcal{R}_{\overline{U}}\widehat{\mathcal{J}}^{\prime}_{j}\|_{\infty,\sigma}$
(continuous version of definition (53).

Again, it is interesting to note that there are two independent ways to make
$\epsilon_{\sigma}$ small. One is a) by reducing
$\sum_{j}|\widehat{U}^{\prime}_{j}|$, which is a measure of base flow
unsteadiness, the other b) is by increasing $\sigma$.

##  Appendix D Extension of theory to quasiperiodic flows

The theory presented in §3 appears to extend to quasi-periodic flows. Indeed,
$m$-torii possess $m$ basic incommensurate frequencies
$\boldsymbol{\Omega}=[\omega_{0},\omega_{1},\dots,\omega_{m-1}]^{T}$, hence
their discrete Fourier spectrum is indexed by $\mathbb{Z}^{m}$ instead of
$\mathbb{Z}$. By simply introducing a bijection $\mathbf{b}$ from $\mathbb{Z}$
to $\mathbb{Z}^{m}$ such that $\mathbf{b}(0)=(0,\dots,0)$, it seems that we
can reuse the same formalism as in §3.4.2. Using this bijection, the Fourier
expansion of the quasi-periodic Jacobian reads

$\mathsfbi{J}(t;\tau_{0})=\mathsfbi{J}_{\overline{\mathbf{U}}}+\sum_{k}\widehat{\mathsfbi{J}}^{\prime}_{k}\mathrm{e}^{\mathrm{i}\mathbf{b}(k)\cdot\boldsymbol{\Omega}(t+\tau_{0})},$
(79)

and the random phase $\phi\in[0,2\pi)$ is now replaced with a random relative
initial time $\tau_{0}$ at which point the linear forcing is switched on. The
hatted quantities are now defined as harmonic averages (Arbabi & Mezić, 2017)

$\widehat{\mathsfbi{J}}^{\prime}_{k}:=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\mathsfbi{J}^{\prime}(t;\tau_{0})\mathrm{e}^{-\mathrm{i}\mathbf{b}(k)\cdot\boldsymbol{\Omega}(t+\tau_{0})}\,\mathrm{d}t.$
(80)

In general, $\mathbf{b}$ is not odd, i.e. $\mathbf{b}(k)\neq-\mathbf{b}(-k)$
for $k\neq 0$, hence we do not have Hermitian symmetry, i.e.
$\widehat{\mathsfbi{J}}^{\prime}_{-k}\neq\widehat{\mathsfbi{J}}^{\prime*}_{k}$.
However, since $\mathsfbi{J}^{\prime}$ is real, for all $k\neq 0$ there must
exist $p$ such that
$\widehat{\mathsfbi{J}}^{\prime}_{p}=\widehat{\mathsfbi{J}}^{{}^{\prime}*}_{k}$.
We now provide an example of a bijection from $\mathbb{Z}$ to
$\mathbb{Z}^{m}$. Cantor’s pairing function
$\pi(\alpha_{1},\alpha_{2}):=\frac{1}{2}(\alpha_{1}+\alpha_{2})(1+\alpha_{1}+\alpha_{2})+\alpha_{2}$
defines a well-known bijection from $\mathbb{N}^{2}$ to $\mathbb{N}$. It may
be generalized by recurrence over $m\geq 2$, using
$\pi_{m}(\alpha_{1},\dots,\alpha_{m}):=\pi(\pi_{m-1}(\alpha_{1},\dots,\alpha_{m-1}),\alpha_{m})$
and $\pi_{1}:=\mathrm{Id}$; the $\pi_{m}$ are called Cantor tuple functions
and represent bijections from $\mathbb{N}^{m}$ to $\mathbb{N}$. A well-known
bijection from $\mathbb{Z}$ to $\mathbb{N}$ is provided by the function $f$
such that $f(\alpha)=2\alpha$ if $\alpha\geq 0$ and $f(\alpha)=-2\alpha-1$
otherwise. By composing $f$ and $\pi_{m}$, a bijection $\mathbf{b}$ from
$\mathbb{Z}$ to $\mathbb{Z}^{m}$ may be formed, which satisfies the condition
$\mathbf{b}(0)=(0,\dots,0)$.

Using this bijection, we recover an equivalent of the harmonic transfer
operator in the quasi-periodic regime. The block
$\underline{\mathsfbi{H}}_{jk}$ maps the input $\mathbf{f}$ at frequency
$s+\mathrm{i}\mathbf{b}(k)\cdot\boldsymbol{\Omega}$ to the output $\mathbf{u}$
at frequency $s+\mathrm{i}\mathbf{b}(j)\cdot\boldsymbol{\Omega}$. For
$j,k\in\mathbb{Z}$, the general terms of the two infinite block-matrices
$\underline{\mathsfbi{D}}(s)$ and $\underline{\mathsfbi{T}}(\phi)$ simply
change to

$\displaystyle\underline{\mathsfbi{D}}_{jk}(s)$
$\displaystyle:=[(s+\mathrm{i}\mathbf{b}(j)\cdot\boldsymbol{\Omega})\mathsfbi{I}-\mathsfbi{J}_{\overline{\mathbf{U}}}]\delta_{jk},$
(81) $\displaystyle\underline{\mathsfbi{T}}_{jk}(\tau_{0})$
$\displaystyle:=\widehat{\mathsfbi{J}}^{\prime}_{j-k}\mathrm{e}^{\mathrm{i}(\mathbf{b}(j)-\mathbf{b}(k))\cdot\boldsymbol{\Omega}\tau_{0}}.$
(82)

We conclude, as in the periodic case, that the resolvent operator
$\mathsfbi{R}_{\overline{\mathbf{U}}}$ about the mean Jacobian approximates
the mean resolvent operator $\mathsfbi{R}_{0}$ and that corrective terms are
of order 2 with respect to $\|\mathsfbi{J}^{\prime}\|_{\mathcal{W}}$.

Following Mezić (2020), an equivalent to Floquet’s theorem in the quasi-
periodic regime, due to Sell (1978), may also be formulated: under appropriate
conditions, there exists a quasi-periodic transformation
$\mathsfbi{P}(t;\tau_{0})$ and a constant matrix $\mathsfbi{K}$ such that the
system (4) admits a propagator of the form
$\mathbf{\Phi}(t,t^{\prime};\tau_{0})=\mathsfbi{P}(t;\tau_{0})\mathrm{e}^{\mathsfbi{K}(t-t^{\prime})}\mathsfbi{P}^{-1}(t^{\prime};\tau_{0})$.
Direct and adjoint quasi-Floquet modes $\mathsfbi{V}(t;\tau_{0})$ and
$\mathsfbi{W}(t;\tau_{0})$ may be defined in the same way as before, but they
are now quasi-periodic instead of periodic, and adjoint with respect to the
inner product
$(\mathbf{a},\mathbf{b}):=\lim_{T\to\infty}(1/T)\int_{0}^{T}\mathbf{a}^{H}\cdot\mathbf{b}\,\mathrm{d}t$.
The propagator may then be written as

$\mathbf{\Phi}(t,t^{\prime};\tau_{0})=\mathsfbi{V}(t;\tau_{0})\mathrm{e}^{\mathbf{\Lambda}(t-t^{\prime})}\mathsfbi{W}^{H}(t^{\prime};\tau_{0}).$
(83)

with

$\mathsfbi{V}(t;\tau_{0})=\sum_{k}\widehat{\mathsfbi{V}}_{k}\mathrm{e}^{\mathrm{i}\mathbf{b}(k)\cdot\boldsymbol{\Omega}(t+\tau_{0})},\quad\mathsfbi{W}(t;\tau_{0})=\sum_{k}\widehat{\mathsfbi{W}}_{k}\mathrm{e}^{\mathrm{i}\mathbf{b}(k)\cdot\boldsymbol{\Omega}(t+\tau_{0})}$
(84)

and $\boldsymbol{\Lambda}$ a diagonal matrix of quasi-Floquet exponents. We
therefore arrive at an expression similar to (19)

$\mathbf{u}(s;\phi)=\sum_{n}\mathrm{e}^{\mathrm{i}\mathbf{b}(n)\cdot\boldsymbol{\Omega}\tau_{0}}\mathsfbi{R}_{n}(s)\mathbf{f}(s-\mathrm{i}\mathbf{b}(n)\cdot\boldsymbol{\Omega}),$
(85)

with the operators

$\mathsfbi{R}_{n}(s)=\sum_{j}\sum_{k=1}^{N}\dfrac{\widehat{\mathbf{v}}^{k}_{j}(\widehat{\mathbf{w}}^{k}_{j-n})^{H}}{s-(\lambda_{k}+\mathrm{i}\mathbf{b}(j)\cdot\boldsymbol{\Omega})}.$
(86)

Just as before

$\langle\mathbf{u}(s)\rangle_{\tau_{0}}=\mathsfbi{R}_{0}(s)\mathbf{f}(s).$
(87)

##  Appendix E Maximum Lyapunov exponent in various configurations

The maximum Lyapunov (64) is evaluated from the partial-state measurement
$\mathbf{y}=[y_{\mathrm{a}},y_{\mathrm{b}},y_{\mathrm{c}}]^{T}$ as

$\mathrm{MLE}=\lim_{t\to\infty}\dfrac{1}{t}\ln\dfrac{\|\langle\mathbf{y}(t)\rangle_{\tau_{0}}\|}{\|\langle\mathbf{y}(0)\rangle_{\tau_{0}}\|}.$
(88)

In the periodic case, taking an ensemble average with respect to $\tau_{0}$ is
equivalent to taking an ensemble average with respect to $\phi$. Results are
shown in figure 16 for the five cases considered. Even though the timeseries
are not long enough to reach full convergence of the maximum Lyapunov exponent
($\mathrm{MLE}$), they clearly demonstrate the three sought behaviours: (a)
$\mathrm{MLE}<0$, (b) $\mathrm{MLE}=0$, (c) $\mathrm{MLE}>0$. The noise
amplifier has negative MLE which depends on the noise intensity: the larger
$\sigma_{w}$, the larger $|\mathrm{MLE}|$. Both the periodic and quasiperiodic
regimes have null MLE. In the quasiperiodic case, the number of zero Lyapunov
exponents provides the dimensionality of the torus (Oteski et al., 2015), i.e.
the number of basic incommensurate frequencies, but here we only computed the
maximum one. In the chaotic case (fluidic pinball at $Re=120$), we find
$\mathrm{MLE}\approx 0.02$.

Figure 16: Estimation of maximum Lyapunov exponent (MLE) for the 5
simulations, based on the mean impulse response; (a) backward-facing step flow
for (ai) $\sigma_{w}=\sqrt{10}$ and (aii) $\sigma_{w}=10$, fluidic pinball at
(bi) $Re=100$, (bii) $Re=110$ and (cii) $Re=120$.

##  Appendix F Invariance of the Floquet exponents to a coordinate change

Consider the bijective nonlinear coordinate change $\mathbf{h}$ on the base
flow

$\mathbf{U}_{\mathbf{h}}=\mathbf{h}(\mathbf{U}),$ (89)

then the invertible periodic matrix

$\mathsfbi{N}(t;\phi)=\nabla_{\mathbf{U}(t;\phi)}\mathbf{h}$ (90)

characterizes the change of coordinate of the linear variables

$\mathbf{u}_{\mathbf{h}}=\mathsfbi{N}(t;\phi)\mathbf{u},\quad\mathbf{f}_{\mathbf{h}}=\mathsfbi{N}(t;\phi)\mathbf{f}.$
(91)

The propagator of the LTP system for the new input/output variables is simply
given by

${\mathbf{\Phi}_{\mathbf{h}}(t,t^{\prime};\phi)=\mathsfbi{V}_{\mathbf{h}}(t;\phi)\mathrm{e}^{\mathbf{\Lambda}(t-t^{\prime})}\mathsfbi{W}_{\mathbf{h}}^{H}(t^{\prime};\phi)}$
(92)

where direct/adjoint Floquet modes change basis

$\mathsfbi{V}_{\mathbf{h}}(t;\phi)=\mathsfbi{N}(t;\phi)\mathsfbi{V}(t;\phi),\quad\mathsfbi{W}_{\mathbf{h}}(t^{\prime};\phi)=\mathsfbi{N}(t^{\prime};\phi)\mathsfbi{W}(t^{\prime};\phi)$
(93)

but Floquet exponents in the fundamental strip remain characterized by the
same matrix $\mathbf{\Lambda}$.

## References

  * Arbabi & Mezić (2017) Arbabi, H. & Mezić, I. 2017 Study of dynamics in post-transient flows using Koopman mode decomposition. Phys. Rev. Fluids 2, 124402.
  * Barkley (2006) Barkley, D. 2006 Linear analysis of the cylinder wake mean flow. Europhys. Lett. 75, 750.
  * Beneddine et al. (2016) Beneddine, S., Sipp, D., Arnault, A., Dandois, J. & Lesshafft, L. 2016 Conditions for validity of mean flow stability analysis. J. Fluid Mech. 798, 485–504.
  * Beneddine et al. (2017) Beneddine, S., Yegavian, R., Sipp, D. & Leclaire, B. 2017 Unsteady flow dynamics reconstruction from mean flow and point sensors: an experimental study. J. Fluid Mech. 824, 174–201.
  * Bengana et al. (2019) Bengana, Y., Loiseau, J.-Ch., Robinet, J.-Ch. & Tuckerman, L. S. 2019 Bifurcation analysis and frequency prediction in shear-driven cavity flow. J. Fluid Mech. 875, 725–757.
  * Bengana & Tuckerman (2019) Bengana, Y. & Tuckerman, L. S. 2019 Spirals and ribbons in counter-rotating Taylor–Couette flow: Frequencies from mean flows and heteroclinic orbits. Phys. Rev. Fluids 4, 044402\.
  * Bengana & Tuckerman (2021) Bengana, Y. & Tuckerman, L. S 2021 Frequency prediction from exact or self-consistent mean flows. Phys. Rev. Fluids 6, 063901.
  * Butler & Farrell (1993) Butler, K. M. & Farrell, B. F. 1993 Optimal perturbations and streak spacing in wall-bounded turbulent shear flow. Phys. Fluids A: Fluid Dynamics 5, 774–777.
  * Conway (1990) Conway, J. B. 1990 A course in functional analysis, 2nd edition, , vol. 96. Springer.
  * Cossu et al. (2009) Cossu, C., Pujals, G. & Depardon, S. 2009 Optimal transient growth and very large–scale structures in turbulent boundary layers. J. Fluid Mech. 619, 79–94.
  * Črnjarić-Žic et al. (2019) Črnjarić-Žic, N., Maćešić, S. & Mezić, I. 2019 Koopman operator spectrum for random dynamical systems. Journal of Nonlinear Science pp. 1–50.
  * Dahan et al. (2012) Dahan, J. A., Morgans, A. S. & Lardeau, S. 2012 Feedback control for form-drag reduction on a bluff body with a blunt trailing edge. J. Fluid Mech. 704, 360–387.
  * Dalla Longa et al. (2017) Dalla Longa, L., Morgans, A. S. & Dahan, J. A. 2017 Reducing the pressure drag of a d-shaped bluff body using linear feedback control. Theo. Comput. Fluid Dyn. pp. 1–11.
  * Del Álamo & Jimenez (2006) Del Álamo, J. C. & Jimenez, J. 2006 Linear energy amplification in turbulent channels. J. Fluid Mech. 559, 205–213.
  * Deng et al. (2020) Deng, N., Noack, B. R., Morzyński, M. & Pastur, L. R. 2020 Low-order model for successive bifurcations of the fluidic pinball. J. Fluid Mech. 884.
  * Evstafyeva et al. (2017) Evstafyeva, O., Morgans, A. S. & Dalla Longa, L. 2017 Simulation and feedback control of the ahmed body flow exhibiting symmetry breaking behaviour. J. Fluid Mech. 817.
  * Franceschini et al. (2022) Franceschini, L., Sipp, D., Marquet, O., Moulin, J. & Dandois, J. 2022 Identification and reconstruction of high-frequency fluctuations evolving on a low-frequency periodic limit cycle: application to turbulent cylinder flow. J. Fluid Mech. 942.
  * Garnaud et al. (2013) Garnaud, X., Lesshafft, L., Schmid, P. J. & Huerre, P. 2013 The preferred mode of incompressible jets: linear frequency response analysis. J. Fluid Mech. 716, 189–202.
  * Gelb & Vander Velde (1968) Gelb, A. & Vander Velde, W. E. 1968 Multiple-input describing functions and nonlinear system design. McGraw-Hill.
  * Gohberg et al. (2003) Gohberg, I., Kaashoek, M.s A & Spitkovsky, I. M. 2003 An overview of matrix factorization theory and operator applications. Factorization and integrable systems pp. 1–102.
  * Gómez et al. (2016) Gómez, F., Blackburn, H. M., Rudman, M., Sharma, A. S. & McKeon, B. J. 2016 A reduced-order model of three-dimensional unsteady flow in a cavity based on the resolvent operator. J. Fluid Mech. 798, R2.
  * Hammond & Redekopp (1997) Hammond, D. A. & Redekopp, L. G. 1997 Global dynamics of symmetric and asymmetric wakes. J. Fluid Mech. 331, 231–260.
  * Hecht (2012) Hecht, F. 2012 New development in FreeFem++. J. Numer. Math. 20 (3-4), 251–265.
  * Herrmann et al. (2021) Herrmann, B., Baddoo, P. J., Semaan, R., Brunton, S. L. & McKeon, B. J. 2021 Data-driven resolvent analysis. J. Fluid Mechanics 918, A10.
  * Hervé et al. (2012) Hervé, A., Sipp, D., Schmid, P. J. & Samuelides, M. 2012 A physics-based approach to flow control using system identification. J. Fluid Mech. 702, 26–58.
  * Hwang & Cossu (2010) Hwang, Y. & Cossu, C. 2010 Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow. J. Fluid Mech. 664, 51–73.
  * Illingworth et al. (2018) Illingworth, S. J., Monty, J. P. & Marusic, I. 2018 Estimating large-scale structures in wall turbulence using linear models. J. Fluid Mech. 842, 146–162.
  * Iollo et al. (2000) Iollo, A., Lanteri, S. & Désidéri, J.-A. 2000 Stability properties of pod–galerkin approximations for the compressible navier–stokes equations. Theo. Comput. Fluid Dyn. 13 (6), 377–396.
  * Jeun et al. (2016) Jeun, J., Nichols, J. W. & Jovanović, M. R. 2016 Input-output analysis of high-speed axisymmetric isothermal jet noise. Phys. Fluids 28 (4), 047101.
  * Jovanovic & Bamieh (2001) Jovanovic, M. & Bamieh, B. 2001 Modeling flow statistics using the linearized Navier–Stokes equations. In Proceedings of the 40th IEEE Conference on Decision and Control, , vol. 5, pp. 4944–4949. IEEE.
  * Jovanovic & Georgiou (2010) Jovanovic, M. & Georgiou, T. 2010 Reproducing second order statistics of turbulent flows using linearized Navier–Stokes equations with forcing. In 63rd Annual Meeting of the APS Division of Fluid Dynamics.
  * Jovanović (2021) Jovanović, M. R. 2021 From bypass transition to flow control and data-driven turbulence modeling: an input–output viewpoint. Annu. Rev. Fluid Mech. 53, 311–345.
  * Karban et al. (2020) Karban, U., Bugeat, B., Martini, E., Towne, A., Cavalieri, A.V.G., Lesshafft, L., Agarwal, A., Jordan, P. & Colonius, T. 2020 Ambiguity in mean-flow-based linear analysis. J. Fluid Mech. 900, R5.
  * Khalil (2002) Khalil, H. K. 2002 Nonlinear systems; 3rd ed.. Prentice-Hall.
  * Krishna Kumar & Kulkarni (2015) Krishna Kumar, G. & Kulkarni, S. H. 2015 Ann. Funct. Anal. 6, 148–169.
  * Leclercq et al. (2019) Leclercq, C., Demourant, F., Poussot-Vassal, C. & Sipp, D. 2019 Linear iterative method for closed-loop control of quasiperiodic flows. J. Fluid Mech. 868, 26–65.
  * Lee et al. (1990) Lee, M. J., Kim, J. & Moin, P. 1990 Structure of turbulence at high shear rate. J. Fluid Mech. 216, 561–583.
  * Liu et al. (2021) Liu, Q., Sun, Y., Yeh, C.-A., Ukeiley, L. S., Cattafesta, L. N. & Taira, K. 2021 Unsteady control of supersonic turbulent cavity flow based on resolvent analysis. J. Fluid Mech. 925, A5.
  * Luhar et al. (2014) Luhar, M., Sharma, A. S. & McKeon, B. J 2014 Opposition control within the resolvent analysis framework. J. Fluid Mech. 749, 597–626.
  * Madhusudanan et al. (2019) Madhusudanan, A., Illingworth, S. J. & Marusic, I. 2019 Coherent large-scale structures from the linearized navier–stokes equations. J. Fluid Mech. 873, 89–109.
  * Magruder et al. (2018) Magruder, C. C., Gugercin, S. & Beattie, C. A 2018 Linear time-periodic dynamical systems: an ${H}_{2}$ analysis and a model reduction framework. Mathematical and Computer Modelling of Dynamical Systems 24, 119–142.
  * Malkus (1956) Malkus, W.V.R. 1956 Outline of a theory of turbulent shear flow. J. Fluid Mech. 1, 521–539.
  * Mantič-Lugo et al. (2014) Mantič-Lugo, V., Arratia, C. & Gallaire, F. 2014 Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake. Phys. Rev. Lett. 113, 084501.
  * Mantič-Lugo et al. (2015) Mantič-Lugo, V., Arratia, C. & Gallaire, F. 2015 A self-consistent model for the saturation dynamics of the vortex shedding around the mean flow in the unstable cylinder wake. Phys. Fluids 27, 074103.
  * McKeon & Sharma (2010) McKeon, B. J. & Sharma, A. S. 2010 A critical-layer framework for turbulent pipe flow. J. Fluid Mech. 658, 336–382.
  * Meliga (2017) Meliga, P. 2017 Harmonics generation and the mechanics of saturation in flow over an open cavity: a second-order self-consistent description. J. Fluid Mech. 826, 503–521.
  * Mezić (2013) Mezić, I. 2013 Analysis of fluid flows via spectral properties of the Koopman operator. Annu. Rev. Fluid Mech. 45, 357–378.
  * Mezić (2020) Mezić, I. 2020 Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry. J. Nonlin. Sci. 30, 2091–2145.
  * Mezić & Surana (2016) Mezić, I. & Surana, A. 2016 Koopman mode decomposition for periodic/quasi-periodic time dependence. IFAC-PapersOnLine 49 (18), 690–697.
  * Mittal (2008) Mittal, S. 2008 Global linear stability analysis of time-averaged flows. Int. J. Num. Meth. Fluids 58, 111–118.
  * Moarref & Jovanović (2012) Moarref, R. & Jovanović, M. R. 2012 Model-based design of transverse wall oscillations for turbulent drag reduction. J. Fluid Mech. 707, 205–240.
  * Moarref et al. (2014) Moarref, R., Jovanović, M. R., Tropp, J. A., Sharma, A. S. & McKeon, B. J. 2014 A low-order decomposition of turbulent channel flow via resolvent analysis and convex optimization. Phys. Fluids 26 (5), 051701.
  * Morra et al. (2019) Morra, P., Semeraro, O., Henningson, D. S. & Cossu, C. 2019 On the relevance of reynolds stresses in resolvent analyses of turbulent wall-bounded flows. J. Fluid Mech. 867, 969–984.
  * Noack et al. (2003) Noack, B. R., Afanasiev, K., Morzyński, M., Tadmor, G. & Thiele, F. 2003 A hierarchy of low-dimensional models for the transient and post-transient cylinder wake. J. Fluid Mech. 497, 335–363.
  * Noiray et al. (2008) Noiray, N., Durox, D., Schuller, T. & Candel, S. 2008 A unified framework for nonlinear combustion instability analysis based on the flame describing function. J. Fluid Mech. 615, 139–167.
  * Oteski et al. (2015) Oteski, L., Duguet, Y., Pastur, L. & Le Quéré, P. 2015 Quasiperiodic routes to chaos in confined two-dimensional differential convection. Phys. Rev. E 92, 043020.
  * Padovan et al. (2020) Padovan, A., Otto, S. E. & Rowley, C. W. 2020 Analysis of amplification mechanisms and cross-frequency interactions in nonlinear flows via the harmonic resolvent. J. Fluid Mech. 900, A14.
  * Pickering et al. (2021) Pickering, E., Rigas, G., Schmidt, O. T., Sipp, D. & Colonius, T. 2021 Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets. J. Fluid Mech. 917.
  * Pier (2002) Pier, B. 2002 On the frequency selection of finite-amplitude vortex shedding in the cylinder wake. J. Fluid Mech. 458, 407–417.
  * Pujals et al. (2009) Pujals, G., García-Villalba, M., Cossu, C. & Depardon, S. 2009 A note on optimal transient growth in turbulent channel flows. Phys. Fluids 21, 015109.
  * Sartor et al. (2015) Sartor, F., Mettot, C. & Sipp, D. 2015 Stability, receptivity, and sensitivity analyses of buffeting transonic flow over a profile. AIAA Journal 53, 1980–1993.
  * Schmid (2010) Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–28.
  * Schmidt et al. (2018) Schmidt, O. T., Towne, A., Rigas, G., Colonius, T. & Brès, G. A. 2018 Spectral analysis of jet turbulence. J. Fluid Mech. 855, 953–982.
  * Sell (1978) Sell, G. R. 1978 The structure of a flow in the vicinity of an almost periodic motion. Journal of Differential Equations 27 (3), 359–393.
  * Semeraro et al. (2016a) Semeraro, O., Jaunet, V., Jordan, P., Cavalieri, A. V. & Lesshafft, L. 2016a Stochastic and harmonic optimal forcing in subsonic jets. p. 2935.
  * Semeraro et al. (2016b) Semeraro, O., Lesshafft, L., Jaunet, V. & Jordan, P. 2016b Modeling of coherent structures in a turbulent jet as global linear instability wavepackets: Theory and experiment. Int. J. Heat Fluid Flow 62, 24–32.
  * Sipp & Lebedev (2007) Sipp, D. & Lebedev, A. 2007 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.