Mean resolvent operator of a statistically steady flow

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 ${{\boldsymbol{\mathsf{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., J. Fluid Mech., vol. 704, 2012, pp. 360–387). 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 ${{\boldsymbol{\mathsf{R}}}}_0$ is well-approximated by the well-known resolvent operator about the mean flow. Although the theory presented in this paper extends only to quasi-periodic flows, the definition of ${{\boldsymbol{\mathsf{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 quasi-periodic, chaotic and stochastic two-dimensional incompressible flows.

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 ones, is generally more complex. Turbulent flows certainly are not time-invariant, 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 interact not with a steady base flow but with 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, Morgans & Lardeau 2012; Dalla Longa, Morgans & Dahan 2017; Evstafyeva, Morgans & Dalla Longa 2017) that 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 exist theoretically, 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 available currently to understand it.

1.3. Option B: linearizing about the mean flow

A second option consists in assuming that the linear response $\boldsymbol {u}$ to infinitesimal forcing $\boldsymbol {f}$ is characterized by the Jacobian matrix ${{\boldsymbol{\mathsf{J}}}}_{\bar {\boldsymbol {U}}}$ evaluated about the mean flow $\bar {\boldsymbol {U}}$, i.e.

(1.1)\begin{equation} \mathrm{d}_t\boldsymbol{u}={{\boldsymbol{\mathsf{J}}}}_{\bar{\boldsymbol{U}}}\boldsymbol{u}+\boldsymbol{f}.\end{equation}

As should be clear from (1.1), we adopt in this paper (except in Appendix C) a spatially discrete framework. The flow will be assumed incompressible and discretized over $N$ degrees of freedom. The resolvent operator about the mean flow,

(1.2)\begin{equation} {{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s)=(s{{\boldsymbol{\mathsf{I}}}}-{{\boldsymbol{\mathsf{J}}}}_{\bar{\boldsymbol{U}}})^{{-}1}, \end{equation}

then describes an LTI input–output dynamics in the frequency domain ($s$ being the Laplace variable), and has been used successfully for open-loop (Moarref & Jovanović 2012; Luhar, Sharma & McKeon 2014; Toedtli, Luhar & McKeon 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 a 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 ${{\boldsymbol{\mathsf{J}}}}_{\bar {\boldsymbol {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 real zero imaginary frequency (RZIF; Turton, Tuckerman & Barkley 2015; Bengana & Tuckerman 2019, 2021; Bengana et al. 2019), 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, Arratia & Gallaire 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 be justified only 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 used extensively to predict second-order statistics (Jovanovic & Bamieh 2001; Jovanovic & Georgiou 2010; Moarref et al. 2014; Beneddine et al. 2016; Zare, Jovanović & Georgiou 2017; Towne, Lozano-Durán & Yang 2020) and in particular coherent flow structures characterized by spectral proper orthogonal decomposition (POD) modes (Semeraro et al. 2016a,b; Beneddine et al. 2016; Schmidt et al. 2018; Towne, Schmidt & Colonius 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, Sipp & McKeon 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, Mettot & Sipp 2015; Jeun, Nichols & Jovanović 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, Kim & Moin (1990), Butler & Farrell (1993), Del Álamo & Jimenez (2006), Cossu, Pujals & Depardon (2009), Pujals et al. (2009) and 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 said already, the first difficulty in terms of interpretation arises from linearizing about a quantity, the mean flow, that is not an invariant solution of the governing equations. In the input–output framework of McKeon & Sharma (2010), the resolvent operator ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {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, Monty & Marusic 2018; Morra et al. 2019; Madhusudanan, Illingworth & Marusic 2019; Pickering et al. 2021; Symon et al. 2019), or through a low-rank state-feedback 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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ is a useful, yet suboptimal linear input–output representation. Worse still, (Karban et al. 2020) showed recently 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 will see, the new operator will help us to understand the ‘dynamic linearity’ phenomenon and provide a physical interpretation to ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$. Its poles will also exactly verify the RZIF property in the periodic case, whereas ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ does so only approximately.

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 keep only 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{\rm \pi} /T$. We represent in figure 1 the response $\boldsymbol {u}$ to an infinitesimal forcing signal $\boldsymbol {f}$ started at time $\tau _0$ relative to an arbitrary origin of time $\tau =0$ on the base flow trajectory $\tilde {\boldsymbol {U}}(\tau )$. Because the base flow is not steady but periodic, the linear response to forcing is not unique but parametrized by the phase $0\leqslant \phi <2{\rm \pi}$ of the limit cycle when the forcing starts. This key parameter for our analysis is defined such that

(1.3)\begin{equation} \omega_0\tau_0=\phi\text{ mod }2{\rm \pi}. \end{equation}

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

(1.4a,b)\begin{equation} \mathrm{d}_t\boldsymbol{u}={{\boldsymbol{\mathsf{J}}}}(t;\phi)\,\boldsymbol{u}+\boldsymbol{f},\quad {{\boldsymbol{\mathsf{J}}}}(t+T;\phi)={{\boldsymbol{\mathsf{J}}}}(t;\phi),\end{equation}

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

(1.5)\begin{equation} \boldsymbol{v}(t)=\textrm{argmin}_{\tilde{\boldsymbol{v}}(t)}\lim_{T'\to\infty}\dfrac{1}{T'}\int_0^{T'} \langle\|\boldsymbol{u}(t;\phi)-\tilde{\boldsymbol{v}}(t)\|^2\rangle_\phi\,\mathrm{d}t, \end{equation}

where $\langle \cdot \rangle _\phi$ denotes an ensemble average with respect to $\phi$. This output is obviously the mean response $\boldsymbol {v}(t)=\langle \boldsymbol {u}(t)\rangle _\phi$. We define the mean resolvent operator ${{\boldsymbol{\mathsf{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:

(1.6)\begin{equation} \boxed{\text{mean resolvent }{{\boldsymbol{\mathsf{R}}}}_0: \boldsymbol{f}(s)\mapsto\langle \boldsymbol{u}(s)\rangle_\phi}.\end{equation}

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

Figure 1. (a) Periodic base flow $\tilde {\boldsymbol {U}}(\tau )$ for an arbitrarily defined origin of time $\tau =0$. (b) Linear responses $\boldsymbol {u}(t;\phi _i)$ to the same infinitesimal forcing $\boldsymbol {f}(t)$ started at different instants $\tau =\tau ^i_0$. Like the underlying base flow $\boldsymbol {U}(t;\phi _i)={\tilde {\boldsymbol {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{\rm \pi}$. (c) The response to linear forcing being a distribution, the statistically optimal LTI model is obtained by computing the mean response $\langle \boldsymbol {u}(t)\rangle _\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 an 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 ${{\boldsymbol{\mathsf{R}}}}_0$, while ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ is a (reduced-order) approximation of ${{\boldsymbol{\mathsf{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. quasi-periodic, 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 quasi-periodic 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 centres of the ‘front’, ‘top’ and ‘bottom’ triangles are respectively located at $(x_f,y_f)=(-1.5\cos ({\rm \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 this section on the single-input single-output (SISO) transfer between an actuator signal $u$ and a measurement $y$ such that

(2.1a,b)\begin{equation} \boldsymbol{f}(t)=\boldsymbol{B}\,u(t),\quad y(t;\phi)=\boldsymbol{C}^T\,\boldsymbol{u}(t;\phi).\end{equation}

The SISO viewpoint is more practical for the numerical experiment that we perform, and it is also relevant in the context of flow control. More specifically, we will consider three probes $\boldsymbol {C}$ extracting the $y$ velocity component in the wake of the ‘top’ cylinder at $y=y_t$ and respectively $x=2,4,8$ (white dots in figure 3). Measurements of the periodic base flow will be denoted $Y=\boldsymbol {C}^T\boldsymbol {U}$. When necessary, the index $a$, $b$ or $c$ will be used to specify which sensor we are referring to. The quantity $\boldsymbol {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:

(2.2)\begin{equation} \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]^{\rm T},\end{equation}

with $\sigma _x=0.01$, $\sigma _y=0.1$, and the centre $(x_0,y_0)$ of the Gaussian (white triangle in figure 3) relative to the centre $(x_t,y_t)$ of the cylinder being given by $(x_0-x_t,y_0-y_t)=(0,0.55)$.

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

Numerical discretization details (computational domain, mesh, spatial and temporal schemes) are provided in Appendix 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 that it is asymmetric, as noted by Deng et al. (2020). In figure 3(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 row (i) of figure 5 for the three sensors (a–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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ and the mean resolvent operator ${{\boldsymbol{\mathsf{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 $(\boldsymbol {C},\boldsymbol {B})$ may be derived from these operators using

(2.3a,b)\begin{equation} G_{\bar{\boldsymbol{U}}}(s):=\boldsymbol{C}^T{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s)\,\boldsymbol{B}\quad\text{and}\quad\langle G(s)\rangle_\phi:=\boldsymbol{C}^T{{\boldsymbol{\mathsf{R}}}}_0(s)\,\boldsymbol{B}. \end{equation}

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

(2.4)\begin{equation} G(s;\phi,u):=\dfrac{y(s;\phi)}{u(s)}, \end{equation}

which are obtained by taking the ratio between the Laplace transforms of the input $u$ and the output $y$. The quantity $G$ is not associated directly 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{\rm \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

(2.5)\begin{equation} \langle G(s)\rangle_\phi=\langle G(s;\phi,u)\rangle_\phi. \end{equation}

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

(2.6)\begin{align} \mathrm{Var}_\phi \,G&:=\langle |G -\langle G \rangle_\phi|^2\rangle_\phi\nonumber\\ &\phantom{:}= \langle |G|^2 \rangle_\phi - |\langle G \rangle_\phi|^2, \end{align}

and the ratio

(2.7)\begin{equation} \eta(s;u):=\dfrac{\sqrt{\mathrm{Var}_\phi\,{G}}}{|\langle G\rangle_\phi|}\end{equation}

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.

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.

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

(2.8)\begin{equation} \mathrm{Var}_\phi\,\langle G\rangle_{\phi,N_s}=\dfrac{\mathrm{Var}_\phi\,G}{N_s}.\end{equation}

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$.

2.3. Procedure

The nonlinear simulation is initialized at $\tau =0$ with the steady base flow (fixed point) perturbed by $\boldsymbol {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. $\boldsymbol {u}(0)=\boldsymbol {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{\rm \pi} )$, as required.

Note that since impulse responses never decay, we can compute only 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 $\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, a geometric interpretation of the ratio $\eta$ is proposed in figure 4, where three cases are considered: (a) $\eta \ll 1$, (b) $1\leqslant \eta <\sqrt {N_s}$ and (c) $\sqrt {N_s} \ll \eta$, corresponding to the $\{\text {probe},\text {frequency}\}$ 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{\rm \pi}$. The red dot indicates the barycentre 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 centred on the red dot: a black circle with radius approximating $\sqrt {\mathrm {Var}_\phi G}$ (length of black arrow), and a green circle with radius equal to $\sqrt {\mathrm {Var}_\phi \langle G\rangle _{\phi,N_s}}$ (length of green arrow). In figure 4(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 barycentre: the dynamics is quasi-time-invariant for that $\{\text {probe},\text {frequency}\}$ pair. In figure 4(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 figure 4(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 estimate confidently 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)$.

Table 1. Values of the ratio $\eta$ evaluated from $N_s=140$ impulses for three $\{\text {probe},\text {frequency}\}$ pairs.

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–c) over the frequency range $0\leqslant \omega \leqslant 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\leqslant \eta <\sqrt {Ns}$, i.e. dynamics not LTI but mean estimate converged; and dark-grey shading corresponds to $\sqrt {N_s}\leqslant \eta$, i.e. mean estimate not even converged.

Figure 5. Fluidic pinball in the periodic regime at $Re=100$. The labels (a)–(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 (2.7)) 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_{\bar {\boldsymbol {U}}}$ (dashed red). Light-grey shading indicates frequency ranges where $1< \eta <\sqrt {Ns}$, and dark-grey shading corresponds to $\eta \geqslant \sqrt {N_s}$. In the phase plots, the two black-dotted curves indicate a shift of $+/-{\rm \pi}$ with respect to the phase of $G_{\bar {\boldsymbol {U}}}$. For rows (ii), (iii) and (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) of figure 5). 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 or not $\eta$ is small compared to 1, we always observe a qualitative agreement between the Bode diagrams of $G_{\bar {\boldsymbol {U}}}$ and $\langle G\rangle _{\phi }$, as long as the mean estimate is converged. For probe a, the two diagrams are nearly indistinguishable 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_{\bar {\boldsymbol {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 decreases rapidly with $k$ and as the probe moves upstream.

In § 3, we will attempt to elucidate the aforementioned observations. (i) Why are $\langle G\rangle _\phi$ and $G_{\bar {\boldsymbol {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? (iv) Why are higher-order resonances weaker and noticeable only downstream?

3. Theory for incompressible periodic base flows

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

3.1. Phase-dependent frequency response of an LTP system

Denote by $\boldsymbol {\varPhi }(t,t';\phi )$ the propagator from $t'$ to $t$ of the linear time-periodic system (1.4a,b), such that

(3.1)\begin{equation} \boldsymbol{u}(t;\phi)=\boldsymbol{\varPhi}(t,0;\phi)\,\boldsymbol{u}_0+\int_{0}^t \boldsymbol{\varPhi}(t,t';\phi)\,\boldsymbol{f}(t')\,\mathrm{d}t',\end{equation}

where $\boldsymbol {u}_0$ is the initial condition on $\boldsymbol {u}$, regardless of $\phi$. By Floquet's theorem (Magruder, Gugercin & Beattie 2018), we know that there exists a real $T$-periodic matrix ${{\boldsymbol{\mathsf{P}}}}(t;\phi )$, invertible at all times $t$, and a constant complex matrix ${{\boldsymbol{\mathsf{K}}}}$ such that $\boldsymbol {\varPhi }(t,t';\phi )={{\boldsymbol{\mathsf{P}}}}(t;\phi )\,\mathrm {e}^{{{\boldsymbol{\mathsf{K}}}}(t-t')}\,{{\boldsymbol{\mathsf{P}}}}^{-1}(t';\phi )$. There are multiple possible determinations of the matrix ${{\boldsymbol{\mathsf{K}}}}=(1/T)\log (\boldsymbol {\varPhi }(T,0))$ depending on the definition of the logarithm. Assuming ${{\boldsymbol{\mathsf{K}}}}$ to be diagonalizable, ${{\boldsymbol{\mathsf{K}}}}={{\boldsymbol{\mathsf{S}}}}\boldsymbol {\varLambda }{{\boldsymbol{\mathsf{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)\leqslant \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)\geqslant \mathrm {Re}(\lambda _2) \geqslant \cdots \geqslant \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})\leqslant \omega _0(j+1/2)$ in the complex plane. The direct Floquet modes $\boldsymbol {v}^i$ are the columns of the $T$-periodic complex matrix ${{\boldsymbol{\mathsf{V}}}}(t;\phi ):={{\boldsymbol{\mathsf{P}}}}(t;\phi ){{\boldsymbol{\mathsf{S}}}}$. The adjoint Floquet modes $\boldsymbol {w}^i$, such that $(\boldsymbol {w}^i,\boldsymbol {v}^j)=\delta _{ij}$ for the canonical inner product $(\boldsymbol {a},\boldsymbol {b})=1/T\int _0^T\boldsymbol {a}^H\boldsymbol {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 ${{\boldsymbol{\mathsf{W}}}}(t;\phi )$ such that ${{\boldsymbol{\mathsf{W}}}}^H{{\boldsymbol{\mathsf{V}}}}={{\boldsymbol{\mathsf{I}}}}$ at all times (and phases $\phi$). Using this eigendecomposition, the propagator reads

(3.2)\begin{equation} \boldsymbol{\varPhi}(t,t';\phi)={{\boldsymbol{\mathsf{V}}}}(t;\phi)\,\mathrm{e}^{\boldsymbol{\varLambda} (t-t')}\, {{\boldsymbol{\mathsf{W}}}}^H(t';\phi).\end{equation}

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

(3.3a,b)\begin{equation} {{\boldsymbol{\mathsf{V}}}}(t;\phi)=\sum_j \hat{{{\boldsymbol{\mathsf{V}}}}}_j \exp({\mathrm{i}j(\omega_0 t+\phi)}),\quad {{\boldsymbol{\mathsf{W}}}}(t';\phi)=\sum_j \hat{{{\boldsymbol{\mathsf{W}}}}}_j \exp({\mathrm{i}j(\omega_0 t'+\phi)}). \end{equation}

Injecting (3.2) and (3.3a,b) in (3.1) and assuming $\boldsymbol {u}_0=\boldsymbol {0}$, we have

(3.4) \begin{align} &\boldsymbol{u}(t;\phi)\nonumber\\ &\quad=\int_{0}^t {{\boldsymbol{\mathsf{V}}}}(t;\phi)\exp({\boldsymbol{\varLambda} (t-t')})\, {{\boldsymbol{\mathsf{W}}}}^H(t';\phi)\,\boldsymbol{f}(t')\,\mathrm{d}t'\nonumber\\ &\quad=\int_{0}^t \left[\sum_j \hat{{{\boldsymbol{\mathsf{V}}}}}_j\exp({\mathrm{i}j(\omega_0t\,{+}\,\phi)})\right] \exp({\boldsymbol{\varLambda}(t\,{-}\,t')})\left[\sum_l \hat{{{\boldsymbol{\mathsf{W}}}}}^H_l\exp({-\mathrm{i}l(\omega_0t'\,{+}\,\phi)})\right] \boldsymbol{f}(t')\,\mathrm{d}t'\nonumber\\ &\quad=\sum_{j,l}\exp({\mathrm{i}(j-l)\phi})\int_{0}^t\left[ \hat{{{\boldsymbol{\mathsf{V}}}}}_j \exp({(\boldsymbol{\varLambda}+\mathrm{i}j\omega_0{{\boldsymbol{\mathsf{I}}}})(t-t')})\, \hat{{{\boldsymbol{\mathsf{W}}}}}^H_l\right]\left[\exp({\mathrm{i}(j-l)\omega_0 t'})\,\boldsymbol{f}(t')\right]\mathrm{d}t'\nonumber\\ &\quad=\sum_{n}\exp({\mathrm{i}n\phi})\int_{0}^t\left[\sum_{j}\hat{{{\boldsymbol{\mathsf{V}}}}}_j \exp({(\boldsymbol{\varLambda}+\mathrm{i}j\omega_0{{\boldsymbol{\mathsf{I}}}})(t-t')})\, \hat{{{\boldsymbol{\mathsf{W}}}}}^H_{j-n}\right]\left[\exp({\mathrm{i}n\omega_0 t'})\,\boldsymbol{f}(t')\right]\mathrm{d}t'. \end{align}

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

(3.5)\begin{equation} \boxed{ \boldsymbol{u}(s;\phi)=\sum_n \mathrm{e}^{\mathrm{i}n\phi}\,{{\boldsymbol{\mathsf{R}}}}_n(s)\,\boldsymbol{f}(s-\mathrm{i}n\omega_0)},\end{equation}

with the operators

(3.6)\begin{equation} \boxed{ {{\boldsymbol{\mathsf{R}}}}_n(s)=\sum_j \sum_{k=1}^N\dfrac{ \hat{\boldsymbol{v}} ^k_j{(\hat{\boldsymbol{w}}^k_{j-n})^H} }{s-(\lambda_k+\mathrm{i}j\omega_0)}},\end{equation}

where $\hat {\boldsymbol {v}}^k_j$ and $\hat {\boldsymbol {w}}^k_j$ designate the $k\text {th}$ columns of $\hat {{{\boldsymbol{\mathsf{V}}}}}_j$ and $\hat {{{\boldsymbol{\mathsf{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 }\,{{\boldsymbol{\mathsf{R}}}}_n(s)$ for $n\neq 0$. The poles of the ${{\boldsymbol{\mathsf{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 that predicts the mean output $\langle \boldsymbol {u}(s)\rangle _\phi$ for a given input $\boldsymbol {f}(s)$, by averaging (3.5) with respect to $\phi$. By doing so, all the cross-frequency transfers vanish, and the only remaining transfer is ${{\boldsymbol{\mathsf{R}}}}_0$:

(3.7)\begin{equation} \boxed{ \langle \boldsymbol{u}(s)\rangle_\phi={{\boldsymbol{\mathsf{R}}}}_0(s)\,\boldsymbol{f}(s)}, \end{equation}

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 $\boldsymbol {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 \boldsymbol {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 $\boldsymbol {v}^1=\mathrm {d}_t \boldsymbol {U}$ (in case of instability $\lambda _1>0$ but there is still a zero exponent $\lambda _j=0$ for some $j>1$), so that $\hat {\boldsymbol {v}}^1_j=\mathrm {i}j\omega _0\hat {\boldsymbol {U}}_j$, where $\hat {\boldsymbol {U}}_j$ denotes the $j\text {th}$ harmonic of the Fourier decomposition of $\boldsymbol {U}(t)$. Therefore, using expression (3.6), the mean resolvent reads

(3.8)\begin{equation} \boxed{ {{\boldsymbol{\mathsf{R}}}}_0(s)=\sum_{j\neq0}\dfrac{ \mathrm{i}j\omega_0\,{\hat{\boldsymbol{U}}}_j(\hat{\boldsymbol{w}}^1_j)^H }{s-\mathrm{i}j\omega_0}+\sum_j\sum_{k=2}^N\dfrac{ \hat{\boldsymbol{v}} ^k_j(\hat{\boldsymbol{w}}^k_j)^H }{s-(\lambda_k+\mathrm{i}j\omega_0)}}. \end{equation}

We see immediately 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 with these poles in the expansion, also called the residuals

(3.9)\begin{equation} \mathrm{Res}_{\mathrm{i}j\omega_0}\,{{\boldsymbol{\mathsf{R}}}}_0(s)=\mathrm{i}j\omega_0\,{\hat{\boldsymbol{U}}}_j(\hat{\boldsymbol{w}}^1_j)^H, \end{equation}

tend very rapidly to 0 in the operator norm as $|j|\to \infty$, because the $j\text {th}$ Fourier harmonics of both $\boldsymbol {U}$ and $\boldsymbol {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 $\hat {\boldsymbol {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 figure 5(i)), hence a better observability $|\boldsymbol {C}^T\hat {\boldsymbol {U}}_j|$ of high-order harmonics $\hat {\boldsymbol {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 us to identify 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 (3.5) and the definition of the variance:

(3.10)\begin{equation} \mathrm{Var}_\phi\,\boldsymbol{u}(s)=\sum_{n\neq 0}\|\boldsymbol{f}(s-\mathrm{i}n\omega_0)\|_{{{\boldsymbol{\mathsf{R}}}}_n^H{{\boldsymbol{\mathsf{R}}}}_n}^2,\end{equation}

where $\|\boldsymbol {a}\|_{{{\boldsymbol{\mathsf{M}}}}}=(\boldsymbol {a}^H{{\boldsymbol{\mathsf{M}}}}\boldsymbol {a})^{1/2}$ is the norm induced by the positive-definite matrix ${{\boldsymbol{\mathsf{M}}}}$. A convenient choice of input signal is one that minimizes the variance, allowing for a minimal amount of realizations. Clearly, an input signal with a broadband spectrum like an impulse $\boldsymbol {f}(t)=\boldsymbol {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 will see in the next subsubsection.

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 (3.10) 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,

(3.11)\begin{equation} \mathrm{Var}_\phi\,G(s;u=\delta)=\sum_{n\neq 0} | \boldsymbol{C}^T{{\boldsymbol{\mathsf{R}}}}_n(s)\,\boldsymbol{B}|^2, \end{equation}

which, according to (3.6), increases downstream since the observability $|\boldsymbol {C}^T\hat {\boldsymbol {v}}_j^k|$ of the Floquet modes probably grows downstream (this is obvious for the first Floquet mode since $\hat {\boldsymbol {v}}^1_j=\mathrm {i}j\omega _0\hat {\boldsymbol {U}}_j$).

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

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

(3.12)\begin{equation} \boldsymbol{f}(s)=\dfrac{\boldsymbol{B}}{s-\mathrm{i}\omega}. \end{equation}

Plugging this expression into the input–output relation (3.5) leads to

(3.13)\begin{equation} \boldsymbol{u}(s;\phi)=\sum_n \mathrm{e}^{\mathrm{i}n\phi}\,\dfrac{{{\boldsymbol{\mathsf{R}}}}_n(s)\,\boldsymbol{B}}{s-\mathrm{i}(\omega+n\omega_0)}.\end{equation}

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 (3.13), in either the transfer operators or the forcing, lead to a permanent contribution. For non-resonant forcing frequencies $\omega \neq n\omega _0$, we have

(3.14)\begin{align} \boldsymbol{u}(t;\phi)&=\mathrm{transient}\nonumber\\ &\quad+{{\boldsymbol{\mathsf{R}}}}_0(\mathrm{i}\omega)\,\boldsymbol{B}\exp({\mathrm{i}\omega t})+\sum_{j\neq 0}\mathrm{i}j\omega_0\,\hat{\boldsymbol{U}}_j(\hat{\boldsymbol{w}}_j^1)^H\boldsymbol{B} \exp({\mathrm{i}j\omega_0t})\nonumber\\ &\quad+\sum_{n\neq 0}\exp({\mathrm{i}n\phi})\left[\vphantom{\sum_{n\neq 0}} {{\boldsymbol{\mathsf{R}}}}_n(\mathrm{i}(\omega+n\omega_0))\,\boldsymbol{B} \exp({\mathrm{i}(\omega+n\omega_0)t})\right.\nonumber\\ &\quad \left.+\sum_{j\neq 0}\mathrm{i}j\omega_0\, \hat{\boldsymbol{U}}_j(\hat{\boldsymbol{w}}_{j-n}^1)^H\boldsymbol{B} \exp({\mathrm{i}j\omega_0t})\right]. \end{align}

The first line of (3.14) collects the transient contributions, and the second line gathers terms that 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 to notice that the Fourier coefficient at the forcing frequency, which may be obtained using a harmonic average

(3.15)\begin{equation} \hat{\boldsymbol{u}}(\omega;\phi)=\lim_{T'\to \infty}\dfrac{1}{T'}\int_0^{T'} \boldsymbol{u}(t;\phi)\,\mathrm{e}^{-\mathrm{i}\omega t}\,\mathrm{d}t, \end{equation}

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

(3.16)\begin{equation} \boxed{\hat{\boldsymbol{u}}(\omega)={{\boldsymbol{\mathsf{R}}}}_0(\mathrm{i}\omega)\,\hat{\boldsymbol{f}}(\omega) }. \end{equation}

In other terms, using small-amplitude harmonic forcings, it is possible to obtain frequency samples of the mean transfer function $\langle G\rangle _\phi =\boldsymbol {C}^T{{\boldsymbol{\mathsf{R}}}}_0\boldsymbol {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) and Evstafyeva et al. (2017), and used therein for LTI system identification on time-varying base flows.

In principle, using a superposition of harmonic forcings

(3.17)\begin{equation} \boldsymbol{f}(t)=\boldsymbol{B}(A_1\,\mathrm{e}^{\mathrm{i}\omega_1 t}+A_2\,\mathrm{e}^{\mathrm{i}\omega_2 t}+\cdots), \end{equation}

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'=0$ to $t$,

(3.18)\begin{equation} \boldsymbol{\varPhi}(t,t'=0;\phi)=\sum_j\sum_{k=1}^N \hat{\boldsymbol{v}} ^k_j\exp({(\lambda_k+\mathrm{i}j\omega_0)t}) \exp({\mathrm{i}j\phi})\,(\boldsymbol{w}^k(0;\phi))^H, \end{equation}

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

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

(3.19)\begin{equation} \langle\boldsymbol{\varPhi}(t,0)\rangle_\phi=\sum_j\sum_{k=1}^N \hat{\boldsymbol{v}} ^k_j\exp({(\lambda_k+\mathrm{i}j\omega_0)t})\,{(\hat{\boldsymbol{w}}^k_j)^H}\end{equation}

based on the phase-averaged Koopman eigenfunctions $\boldsymbol {u}_0\mapsto \langle \mathrm {e}^{\mathrm {i}j\phi }(\boldsymbol {w}^k(0;\phi ))^H\boldsymbol {u}_0\rangle _\phi ={(\hat {\boldsymbol {w}}^k_j)}^H\boldsymbol {u}_0$.

By comparing (3.19) and (3.6) 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

(3.20)\begin{equation} \boxed{ {{\boldsymbol{\mathsf{R}}}}_0(s)=\mathcal{L}[\langle \boldsymbol{\varPhi}(t,0)\rangle_\phi]}, \end{equation}

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

3.4. Connection with LTI dynamics about the mean flow

3.4.1. Connection with the harmonic transfer operator

Evaluating (3.5) at various output frequencies $s+\mathrm {i}k\omega _0$, it is possible to introduce the harmonic transfer operator $\underline {{{\boldsymbol{\mathsf{H}}}}}(s;\phi )$ such that

(3.21)\begin{align} &\underbrace{\begin{pmatrix}\vdots \\\boldsymbol{u}(s-\mathrm{i}\omega_0)\\\boldsymbol{u}(s)\\ \boldsymbol{u}(s+\mathrm{i}\omega_0)\\\vdots\end{pmatrix}}_{\underline{\boldsymbol{u}}(s;\phi)}\nonumber\\ &\quad =\underbrace{\begin{pmatrix} \ddots & \ddots & \ddots & \ddots & \ddots\\ \ddots & {{\boldsymbol{\mathsf{R}}}}_{0}(s-\mathrm{i}\omega_0) & {{\boldsymbol{\mathsf{R}}}}_{{-}1}(s-\mathrm{i}\omega_0)\,\mathrm{e}^{-\mathrm{i}\phi} & {{\boldsymbol{\mathsf{R}}}}_{{-}2}(s-\mathrm{i}\omega_0)\,\mathrm{e}^{{-}2\mathrm{i}\phi} & \ddots\\ \ddots & {{\boldsymbol{\mathsf{R}}}}_{1}(s)\,\mathrm{e}^{\mathrm{i}\phi} & {{\boldsymbol{\mathsf{R}}}}_{0}(s) & {{\boldsymbol{\mathsf{R}}}}_{{-}1}(s)\,\mathrm{e}^{-\mathrm{i}\phi} & \ddots\\ \ddots & {{\boldsymbol{\mathsf{R}}}}_{2}(s+\mathrm{i}\omega_0)\,\mathrm{e}^{2\mathrm{i}\phi} & {{\boldsymbol{\mathsf{R}}}}_{1}(s+\mathrm{i}\omega_0)\,\mathrm{e}^{\mathrm{i}\phi} & {{\boldsymbol{\mathsf{R}}}}_{0}(s+\mathrm{i}\omega_0) & \ddots\\ \ddots & \ddots & \ddots & \ddots & \ddots \end{pmatrix}}_{\underline{{{\boldsymbol{\mathsf{H}}}}}(s;\phi)} \nonumber\\ &\qquad\times \underbrace{\begin{pmatrix}\vdots \\\boldsymbol{f}(s-\mathrm{i}\omega_0)\\\boldsymbol{f}(s)\\ \boldsymbol{f}(s+\mathrm{i}\omega_0)\\\vdots\end{pmatrix}}_{\underline{\boldsymbol{f}}(s)}. \end{align}

The block $\underline {{{\boldsymbol{\mathsf{H}}}}}_{jk}$ of this infinite matrix operator characterizes the transfer from the frequency $s+\mathrm {i}k\omega _0$ of the input $\boldsymbol {f}$ to the frequency $s+\mathrm {i}j\omega _0$ of the output $\boldsymbol {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

(3.22)\begin{equation} \boxed{{{\boldsymbol{\mathsf{R}}}}_0(s)=\underline{{{\boldsymbol{\mathsf{H}}}}}_{00}(s)}.\end{equation}

3.4.2. The harmonic transfer operator as a feedback loop

Decompose the Jacobian as a mean $\bar{\boldsymbol{\mathsf{J}}}$ and a periodic perturbation ${{\boldsymbol{\mathsf{J}}}}'(t;\phi )$, which we expand as a Fourier series

(3.23)\begin{equation} {{\boldsymbol{\mathsf{J}}}}(t;\phi)=\bar{\boldsymbol{\mathsf{J}}}+ \underbrace{\sum_{j} \hat{\boldsymbol{\mathsf{J}}}'_j \exp({\mathrm{i}j(\omega_0 t+{\phi})})}_{{{\boldsymbol{\mathsf{J}}}}'(t;\phi)}.\end{equation}

Because ${{\boldsymbol{\mathsf{J}}}}'(t;\phi )$ is real, the harmonics have Hermitian symmetry, i.e. $\hat {{{\boldsymbol{\mathsf{J}}}}}'_{-j}=\hat {{{\boldsymbol{\mathsf{J}}}}}'^{*}_{j}$, where $(\cdot )^*$ denotes the complex conjugate (not the conjugate transpose $(\cdot )^H$). Moreover, we also obviously have $\hat {{{\boldsymbol{\mathsf{J}}}}}'_0={{\boldsymbol{\mathsf{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:

(3.24)\begin{equation} \bar{{{\boldsymbol{\mathsf{J}}}}}={{\boldsymbol{\mathsf{J}}}}_{\bar{\boldsymbol{U}}}.\end{equation}

Assuming $\boldsymbol {u}^0=0$, taking the Laplace transform of (1.4a,b) and plugging (3.23)-(3.24) yields

(3.25)\begin{equation} s\,\boldsymbol{u}(s;\phi)={{\boldsymbol{\mathsf{J}}}}_{\bar{\boldsymbol{U}}}\boldsymbol{u}(s;\phi)+\sum_{j}\hat{{{\boldsymbol{\mathsf{J}}}}}'_j\,\mathrm{e}^{\mathrm{i}j\phi}\boldsymbol{u}(s-\mathrm{i}j\omega_0;\phi)+\boldsymbol{f}(s), \end{equation}

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, Otto & Rowley 2020; Franceschini et al. 2022):

(3.26)\begin{equation} \underline{{{\boldsymbol{\mathsf{H}}}}}(s;\phi)=(\underline{{{\boldsymbol{\mathsf{D}}}}}(s)-\underline{{{\boldsymbol{\mathsf{T}}}}}(\phi))^{{-}1}.\end{equation}

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

(3.27)$$\begin{gather} \underline{{{\boldsymbol{\mathsf{D}}}}}(s)= \begin{pmatrix} \ddots & \ddots & \ddots & \ddots & \ddots\\ \ddots & (s-\mathrm{i}\omega_0){{\boldsymbol{\mathsf{I}}}}-{{\boldsymbol{\mathsf{J}}}}_{\bar{\boldsymbol{U}}} & {{\boldsymbol{\mathsf{0}}}} & {{\boldsymbol{\mathsf{0}}}} & \ddots\\ \ddots & {{\boldsymbol{\mathsf{0}}}} & s{{\boldsymbol{\mathsf{I}}}}-{{\boldsymbol{\mathsf{J}}}}_{\bar{\boldsymbol{U}}} & {{\boldsymbol{\mathsf{0}}}} & \ddots\\ \ddots & {{\boldsymbol{\mathsf{0}}}} & {{\boldsymbol{\mathsf{0}}}} & (s+\mathrm{i}\omega_0){{\boldsymbol{\mathsf{I}}}}-{{\boldsymbol{\mathsf{J}}}}_{\bar{\boldsymbol{U}}} & \ddots\\ \ddots & \ddots & \ddots & \ddots & \ddots \end{pmatrix}, \end{gather}$$

(3.28)$$\begin{gather}\underline{{{\boldsymbol{\mathsf{T}}}}}(\phi)= \begin{pmatrix} \ddots & \ddots & \ddots & \ddots & \ddots\\ \ddots & {{\boldsymbol{\mathsf{0}}}} & \hat{{{\boldsymbol{\mathsf{J}}}}}_1'^{*}\,\mathrm{e}^{-\mathrm{i}\phi} & \hat{{{\boldsymbol{\mathsf{J}}}}}_2'^{*}\,\mathrm{e}^{{-}2\mathrm{i}\phi} & \ddots\\ \ddots & \hat{{{\boldsymbol{\mathsf{J}}}}}'_1\,\mathrm{e}^{\mathrm{i}\phi} & {{\boldsymbol{\mathsf{0}}}} & \hat{{{\boldsymbol{\mathsf{J}}}}}_1'^{*}\,\mathrm{e}^{-\mathrm{i}\phi} & \ddots\\ \ddots & \hat{{{\boldsymbol{\mathsf{J}}}}}'_2\,\mathrm{e}^{2\mathrm{i}\phi} & \hat{{{\boldsymbol{\mathsf{J}}}}}'_1\,\mathrm{e}^{\mathrm{i}\phi} & {{\boldsymbol{\mathsf{0}}}} & \ddots\\ \ddots & \ddots & \ddots & \ddots & \ddots \end{pmatrix}. \end{gather}$$

By recasting (3.26) in the form

(3.29)\begin{equation} \underline{{{\boldsymbol{\mathsf{H}}}}}(s;\phi)= (\underline{{{\boldsymbol{\mathsf{I}}}}}-\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s)\,\underline{{{\boldsymbol{\mathsf{T}}}}}(\phi))^{{-}1}\,\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s),\end{equation}

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

Figure 6. Block representation of the harmonic transfer operator $\underline {{{\boldsymbol{\mathsf{H}}}}}(s;\phi )$ as a feedback loop between two blocks: $\underline {{{\boldsymbol{\mathsf{D}}}}}^{-1}(s)$ accounts for interactions with the mean flow, and $\underline {{{\boldsymbol{\mathsf{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 {{{\boldsymbol{\mathsf{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 $\boldsymbol {U}'$, 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 will see next.

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

The inverse $(\underline {{{\boldsymbol{\mathsf{I}}}}}- \underline {{{\boldsymbol{\mathsf{D}}}}}^{-1}(s)\,\underline {{{\boldsymbol{\mathsf{T}}}}}(\phi ))^{-1}$ in (3.29) may be expanded as a Neumann series $\sum _{k\geqslant 0}(\underline {{{\boldsymbol{\mathsf{D}}}}}^{-1}(s)\,\underline {{{\boldsymbol{\mathsf{T}}}}}(\phi ))^k$ such that

(3.30)\begin{align} \underline{{{\boldsymbol{\mathsf{H}}}}}(s;\phi)=\underbrace{\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s)}_{\underline{{{\boldsymbol{\mathsf{H}}}}}^0(s)}+\underbrace{\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s)\,\underline{{{\boldsymbol{\mathsf{T}}}}}(\phi)\,\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s)}_{\underline{{{\boldsymbol{\mathsf{H}}}}}^1(s;\phi)}+\underbrace{\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s)\,\underline{{{\boldsymbol{\mathsf{T}}}}}(\phi)\,\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s)\,\underline{{{\boldsymbol{\mathsf{T}}}}}(\phi)\,\underline{{{\boldsymbol{\mathsf{D}}}}}^{{-}1}(s)}_{\underline{{{\boldsymbol{\mathsf{H}}}}}^2(s;\phi)}+\cdots.\end{align}

The series converges if $\|\underline {{{\boldsymbol{\mathsf{D}}}}}^{-1}(s)\,\underline {{{\boldsymbol{\mathsf{T}}}}}(\phi )\|_2<1$ (here, $\|\cdot \|_2$ denotes the spectral norm, i.e. maximum singular value, not the $H_2$ norm.). Assume $\mathrm {Im}(s)=:\sigma >\sigma _{max}$, where $\sigma _{max}\geqslant 0$ is the maximum growth rate among the poles of $\underline {{{\boldsymbol{\mathsf{D}}}}}^{-1}$ and $\underline {{{\boldsymbol{\mathsf{H}}}}}$. Then the series converges if (a) base flow unsteadiness is sufficiently weak or (b) $\sigma$ is sufficiently large (see § C.1). Each contribution $\underline {{{\boldsymbol{\mathsf{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 (3.30) 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 (3.22), i.e.

(3.31)\begin{equation} {{\boldsymbol{\mathsf{R}}}}_0=\underline{{{\boldsymbol{\mathsf{H}}}}}^0_{00}+\underline{{{\boldsymbol{\mathsf{H}}}}}^1_{00}+\underline{{{\boldsymbol{\mathsf{H}}}}}^2_{00}+\cdots.\end{equation}

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

(3.32)$$\begin{gather} \underline{{{\boldsymbol{\mathsf{D}}}}}_{jk}^{{-}1}(s):={{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s+\mathrm{i}j\omega_0)\delta_{jk}, \end{gather}$$

(3.33)$$\begin{gather}\underline{{{\boldsymbol{\mathsf{T}}}}}_{jk}(\phi):=\hat{{{\boldsymbol{\mathsf{J}}}}}'_{j-k}\,\mathrm{e}^{\mathrm{i}(j-k)\phi}, \end{gather}$$

we find

(3.34)$$\begin{gather} \underline{{{\boldsymbol{\mathsf{H}}}}}^0_{00}={{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}, \end{gather}$$

(3.35)$$\begin{gather}\underline{{{\boldsymbol{\mathsf{H}}}}}^1_{00}={{\boldsymbol{\mathsf{0}}}}, \end{gather}$$

(3.36)$$\begin{gather}\underline{{{\boldsymbol{\mathsf{H}}}}}^2_{00}=\sum_j {{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s)\, \hat{{{\boldsymbol{\mathsf{J}}}}}'^{*}_{j}\,{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s+\mathrm{i}j\omega_0)\,\hat{{{\boldsymbol{\mathsf{J}}}}}'_{j}\,{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s), \end{gather}$$

(3.37)$$\begin{gather}\underline{{{\boldsymbol{\mathsf{H}}}}}^3_{00}=\sum_{j,k}{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s)\,\hat{{{\boldsymbol{\mathsf{J}}}}}'^{*}_{j}\,{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s+\mathrm{i}j\omega_0)\, \hat{{{\boldsymbol{\mathsf{J}}}}}'_{j-k}\,{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s+\mathrm{i}k\omega_0)\,\hat{{{\boldsymbol{\mathsf{J}}}}}'_{k}\,{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}(s),\\ \vdots \nonumber \end{gather}$$

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 as the input. The first correction to the resolvent operator about the mean flow then arises at order 2, when successive interactions with $\hat {{{\boldsymbol{\mathsf{J}}}}}'_k$ and its complex conjugate $\hat {{{\boldsymbol{\mathsf{J}}}}}'^{*}_{k}=\hat {{{\boldsymbol{\mathsf{J}}}}}'_{-k}$ occur, allowing the output to be at the same frequency as the input.

Expansion (3.31) may be used to bound the absolute and relative differences between ${{\boldsymbol{\mathsf{R}}}}_0$ and ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$. Introduce the system norm

(3.38)\begin{equation} \|{{\boldsymbol{\mathsf{G}}}}\|_{\infty,\sigma}:=\sup_{\omega\in\mathbb{R}}\|{{\boldsymbol{\mathsf{G}}}}(\sigma+\mathrm{i}\omega)\|_2,\end{equation}

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

(3.39)\begin{equation} \epsilon_\sigma:=\sum_k\|{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}\hat{{{\boldsymbol{\mathsf{J}}}}}'_k\|_{\infty,\sigma},\end{equation}

characterizing the ($\sigma$-dependent) loop gain associated with the feedback interconnection of $\underline {{{\boldsymbol{\mathsf{D}}}}}^{-1}(s)$ and $\underline {{{\boldsymbol{\mathsf{T}}}}}(\phi )$, averaged with respect to $\phi$. Using (3.31) and the definition (3.39) of $\epsilon _\sigma$, we may write

(3.40)\begin{align} \|{{\boldsymbol{\mathsf{R}}}}_0-{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}\|_2(s)&\leqslant\|{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}\|_2(s)\, (\epsilon_\sigma^2 + \epsilon_\sigma^3+\cdots)\nonumber\\ &= \|{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}\|_2(s)\,\dfrac{\epsilon_\sigma^2}{1- \epsilon_\sigma}. \end{align}

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 § C.2). Since $\|{{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}\|_2(s)\to 0$ as either $\sigma \to \infty$ at fixed $\omega$ or $|\omega |\to \infty$ at fixed $\sigma$, the absolute difference in (3.40) also tends to zero in these limits, i.e. $\|{{\boldsymbol{\mathsf{R}}}}_0-{{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}\|_2(s)\to 0$. Moreover, the relative error at fixed $s$,

(3.41)\begin{equation} \boxed{ \dfrac{\|{{\boldsymbol{\mathsf{R}}}}_0-{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}\|_2(s)}{\|{{\boldsymbol{\mathsf{R}}}}_{\bar{\boldsymbol{U}}}\|_2(s)}=O(\epsilon_\sigma^2)},\end{equation}

is order 2, hence vanishes very quickly with the small parameter $\epsilon _\sigma$. As said already, 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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ is a physically relevant operator in incompressible flows is because it approximates ${{\boldsymbol{\mathsf{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 figures 5(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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {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 ${{\boldsymbol{\mathsf{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

(3.42)\begin{equation} \boxed{ {{\boldsymbol{\mathsf{R}}}}_0(s)=\underline{{{\boldsymbol{\mathsf{C}}}}}(s\underline{{{\boldsymbol{\mathsf{I}}}}}-\underline{\boldsymbol{\varLambda}})^{{-}1}\underline{{{\boldsymbol{\mathsf{B}}}}}},\end{equation}

where

(3.43)$$\begin{gather} \underline{{{\boldsymbol{\mathsf{C}}}}}=[\dots,\hat{{{\boldsymbol{\mathsf{V}}}}}_{{-}1},\hat{{{\boldsymbol{\mathsf{V}}}}}_0,\hat{{{\boldsymbol{\mathsf{V}}}}}_{1}\dots], \end{gather}$$

(3.44)$$\begin{gather}\underline{{{\boldsymbol{\mathsf{B}}}}}=[\dots,\hat{{{\boldsymbol{\mathsf{W}}}}}_{{-}1},\hat{{{\boldsymbol{\mathsf{W}}}}}_0,\hat{{{\boldsymbol{\mathsf{W}}}}}_{1}\dots]^H, \end{gather}$$

(3.45)$$\begin{gather}\underline{{{\boldsymbol{\mathsf{I}}}}}_{jk}={{\boldsymbol{\mathsf{I}}}}\delta_{jk}, \end{gather}$$

(3.46)$$\begin{gather}\underline{\boldsymbol{\varLambda}}_{jk}=(\boldsymbol{\varLambda}+\mathrm{i}j\omega_0{{\boldsymbol{\mathsf{I}}}})\delta_{jk}. \end{gather}$$

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

(3.47)$$\begin{gather} \mathrm{d}_t \underline{\boldsymbol{x}} = \underline{\boldsymbol{\varLambda}}\,\underline{\boldsymbol{x}} + \underline{{{\boldsymbol{\mathsf{B}}}}}\,\boldsymbol{f}, \end{gather}$$

(3.48)$$\begin{gather}\langle\boldsymbol{u}\rangle_\phi = \underline{{{\boldsymbol{\mathsf{C}}}}}\, \underline{\boldsymbol{x}}, \end{gather}$$

where the internal state $\underline {\boldsymbol {x}}=[\dots,\boldsymbol {x}_{-1},\boldsymbol {x}_0,\boldsymbol {x}_1,\dots ]^{\rm T}$ is an infinite column vector, even though the input $\boldsymbol {f}$ and the output $\langle \boldsymbol {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 {\varLambda }}$ become negligible, as they involve the Fourier coefficients $\hat {{{\boldsymbol{\mathsf{W}}}}}_j$ in $\underline {{{\boldsymbol{\mathsf{B}}}}}$, and $\hat {{{\boldsymbol{\mathsf{V}}}}}_j$ in $\underline {{{\boldsymbol{\mathsf{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 {\boldsymbol {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 ${{\boldsymbol{\mathsf{R}}}}_0$ and it has only $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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {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 discussed originally in the context of mean flow stability analysis (Turton et al. 2015). Since $\boldsymbol {v}^1=\mathrm {d}_t\boldsymbol {U}$, the associated Koopman modes $\hat {\boldsymbol {v}}_j^1=\mathrm {i}j\omega _0\hat {\boldsymbol {U}}_j$ are parallel to the Fourier components $\hat {\boldsymbol {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), or the dynamics is weakly nonlinear (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 ${{\boldsymbol{\mathsf{R}}}}_0$ and associated modes satisfy the RZIF property exactly, while this is only approximately true for ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {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 concerned with only periodic base flows, but the definition (1.6) of the mean resolvent ${{\boldsymbol{\mathsf{R}}}}_0$ may be generalized easily to any other statistically steady base flow. Instead of taking the mean output with respect to the phase $\phi$, more generally we may take the mean output with respect to the relative time $\tau _0$ at which the forcing signal is released:

(4.1)\begin{equation} \boxed{ \text{mean resolvent } {{\boldsymbol{\mathsf{R}}}}_0: \boldsymbol{f}(s)\mapsto \langle\boldsymbol{u}(s)\rangle_{\tau_0}}. \end{equation}

For such a definition to hold, we assume that the single trajectory $\tilde {\boldsymbol {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 numerical experiments similar to those 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 \cdot \rangle _\phi$ with $\langle \cdot \rangle _{\tau _0}$, since now ‘frequency response realizations’ are parametrized by $\tau _0$.

In § 4.1, we consider the fluidic pinball in the quasi-periodic and chaotic regimes by increasing the Reynolds number to respectively $Re=110$ and $Re=120$ (Deng et al. 2020). For the quasi-periodic 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

(4.2)\begin{equation} \mathrm{MLE} = \lim_{t\to\infty}\dfrac{1}{t}\ln \dfrac{\|\langle\boldsymbol{u}(t)\rangle_{\tau_0}\|}{\|\langle \boldsymbol{u}(0)\rangle_{\tau_0}\|}.\end{equation}

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 quasi-periodic case, we also have $\mathrm {MLE}=0$, for the chaotic regime $\mathrm {MLE}>0$, and 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 applies to only the periodic and quasi-periodic cases.

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

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

The numerical set-up 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 quasi-periodic case and chaotic cases. We notice that the mean flow is asymmetric in the quasi-periodic case, but symmetric in the chaotic case. The recirculation zone is more extended downstream in the latter case.

Figure 7. Fluidic pinball in different regimes: (i) quasi-periodic at $Re=110$, and (ii) chaotic at $Re=120$.(a) Mean velocity norm and streamlines. (b) Phase portraits using $y$ velocity sensors $Y_{b}$ and $Y_{c}$.

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 quasi-periodic 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(b-ii) 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–c). These confirm the qualitatively different dynamics of the flows: the spectrum is discrete in the quasi-periodic case, while it is continuous in the chaotic case, despite the apparent presence of peaks. In the latter case, the peaks possess an intrinsic non-zero bandwidth, unlike in the periodic and quasi-periodic cases, where the ‘thickness’ of the peak arises only from the estimation error of the discrete Fourier transform. It is not possible to infer directly the number of basic incommensurate frequencies in the quasi-periodic case by examining the spectra, but clearly there is energy in a much greater number of frequencies compared to the periodic case (see figure 5(i) for comparison).

Figure 8. Fluidic pinball in the quasi-periodic regime at $Re=110$. The labels (a)–(c) correspond to the three sensor positions. (i) Fourier spectrum of measurement of 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 (2.7)) 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_{\bar {\boldsymbol {U}}}$ (dashed red). Light-grey shading indicates frequency ranges where $1< \eta <\sqrt {N_s}$, and dark-grey shading corresponds to $\eta \geqslant \sqrt {N_s}$. In the phase plots, the two black-dotted curves indicate a shift of $+/-{\rm \pi}$ with respect to the phase of $G_{\bar {\boldsymbol {U}}}$. For rows (ii), (iii) and (iv), all quantities are evaluated on the shifted imaginary axis $\sigma +\mathrm {i}\omega$, with $\sigma =0.01$.

Figure 9. Same 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 preceding one, i.e. $\Delta \tau _0=2$. We assumed that this procedure allowed us to sample a representative part of these attractors with a single realization of $\tilde {\boldsymbol {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_{\bar {\boldsymbol {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 ensured only when $\eta \ll \sqrt {N_s}$, so it is difficult to conclude on these deviations.

In the quasi-periodic 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_{\bar {\boldsymbol {U}}}$. Resonances at high frequencies are not visible in the near wake and manifest only 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 quasi-periodic case). It is interesting to note that the close agreement between $\langle G\rangle _{\tau _0}$ and $G_{\bar {\boldsymbol {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 set-up and code for the nonlinear flow are identical to that of Hervé et al. (2012) and 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 compared 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 $\boldsymbol {F}(t)=\boldsymbol {B}_w\,w(t)$, with $w$ a Gaussian white noise, and $\boldsymbol {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 (2.2) for the definition of $\mathcal {B}$ and the white squares centred at $(x_0,y_0)$ in figures 10a,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).

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 root-mean-square of the amplitude $w(t)$ multiplying the external noise field $\boldsymbol {B}_w$ (white square), such that $\boldsymbol {F}(t)=\boldsymbol {B}_w w(t)$. The actuator is $\boldsymbol {B}$ (white triangle), and $y$ velocity probes are placed at locations a and b in the shear layer (white dots), while a friction measurement (white rectangle) is taken at c. (c) Time series of $Y_{a}$ for $\sigma _w=\sqrt {10}$ (dashed black line) and $\sigma _w=10$ (blue solid line).

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 $\boldsymbol {F}$ produces the nonlinear base flow $\tilde {\boldsymbol {U}}$ and should not be confused with $\boldsymbol {f}=\boldsymbol {B}u$, which produces the linear response $\boldsymbol {u}$ about $\tilde {\boldsymbol {U}}$. The actuator $\boldsymbol {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 the white triangles in figures 10a,b). We again use three probes to monitor the flow response: two $y$ velocity probes in the shear layer at $(2.5,0)$ and $(5,0)$, and a friction sensor at $12.5\leqslant x\leqslant 12.9$ on the lower wall, behind reattachment (see the white dots and rectangles in figures 10a,b). The perturbation and base flow measurements will also be denoted $y=\boldsymbol {C}^T\boldsymbol {u}$ and $Y=\boldsymbol {C}^T\boldsymbol {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 figure 10(b) for $\sigma _w=10$. For stronger perturbations, the recirculation region shrinks, with a reattachment point at $x=10.7$ in the first case, and $x=8.3$ in the second case. Time series 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–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 sensor 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 at $\omega \approx 1$, where linear amplification mechanisms dominate.

Figure 11. Same 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.

Figure 12. Same as figure 11 for the high-noise case $\sigma _w=10$, with $N_s=308$ samples.

4.2.2. Procedure

Again, we use a single realization of $\tilde {\boldsymbol {U}}(\tau )$ corresponding to a single random time series $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 quasi-periodic 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}\geqslant 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

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 those for 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 the base flow power spectrum in row (i)).

4.3. Strongly compressible flows

It was pointed out recently 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 that we want to stress is the assumption (3.24) that 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 in either primitive or conservative formulation. However, we may still use expansion (3.34)–(3.37) if we replace ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ with

(4.3)\begin{equation} {{\boldsymbol{\mathsf{R}}}}_{\bar{{{\boldsymbol{\mathsf{J}}}}}}(s)=(s{{\boldsymbol{\mathsf{I}}}}-\bar{{{\boldsymbol{\mathsf{J}}}}})^{{-}1}, \end{equation}

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, Lanteri & Désidéri 2000), even though this form may not handle shock discontinuities easily 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 led to interesting findings.

First, our framework allows us 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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {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 (or frequencies) are visible in the mean transfer function.

We were able to explain these observations in the periodic and quasi-periodic 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{\rm \pi}$ in the periodic case. Mean transfer functions based on the resulting mean resolvent operator ${{\boldsymbol{\mathsf{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) and Evstafyeva et al. (2017) in the context of flow control.

We then showed that the poles of ${{\boldsymbol{\mathsf{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 to solve the ambiguity of classical resolvent analysis about a mean flow noted recently 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 ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$, appears to be exact when considering the poles of ${{\boldsymbol{\mathsf{R}}}}_0$.

Next, we further investigated the connection between ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ and ${{\boldsymbol{\mathsf{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 ${{\boldsymbol{\mathsf{R}}}}_0$ indirectly incorporates information about the endogenous nonlinear forcing term of McKeon & Sharma (2010), which is not taken into account by ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {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 ${{\boldsymbol{\mathsf{R}}}}_0$ has an infinite-dimensional internal state, while ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ has a finite-dimensional one. However, only a finite number of poles make a significant contribution to ${{\boldsymbol{\mathsf{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 ${{\boldsymbol{\mathsf{R}}}}_0$ is intended to extend that of the usual resolvent operator ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ about the mean flow, therefore any analysis that may be done with ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ (reduced-order modelling, data assimilation of second-order statistics, input–output analysis, flow control, etc.) may, in principle, also be done with (improved) ${{\boldsymbol{\mathsf{R}}}}_0$, using yet to be defined numerical methods. For instance, a compelling prospect would be to compute the singular modes of ${{\boldsymbol{\mathsf{R}}}}_0$ and compare them with the spectral POD modes to see if they align better than the singular modes of ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$. It is hoped that the new operator may help us to 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 address turbulent or compressible flows directly, and future work needs to be done to confirm the relevance of ${{\boldsymbol{\mathsf{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, Maćešić & Mezić 2020). 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 dynamic mode decomposition method (Schmid 2010; Williams, Kevrekidis & Rowley 2015; Herrmann et al. 2021). Recent papers enhancing the predictive power of ${{\boldsymbol{\mathsf{R}}}}_{\bar {\boldsymbol {U}}}$ using optimization techniques (Zare et al. 2017; Pickering et al. 2021) may possibly bear connections with this compression problem.