Mean resolvent operator of statistically steady flows

This paper introduces a new operator relevant to input-output analysis and control of flows in a statistically steady regime far from the steady base flow: the mean resolvent $\mathbf{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. Theory is developed for the periodic case. The poles of the operator are shown to correspond to the Floquet exponents of the system, including purely imaginary poles at multiples of the fundamental frequency. In general, evaluating mean transfer functions from data requires averaging the response to many realizations of the same input. However, in the specific case of harmonic forcings, we show that the mean transfer functions may be identified without averaging: an observation referred to as `dynamic linearity' in the literature (Dahan et al., 2012). For incompressible flows in the weakly unsteady limit, i.e. when amplification of external perturbations by the unsteady part of the periodic Jacobian is small compared to amplification by the mean Jacobian, the mean resolvent $\mathbf{R}_{0}$ is well-approximated by the well-known resolvent operator about the mean-flow $\mathbf{R}_{\overline{\mathbf{U}}} $. Although the theory presented in this paper only extends to quasiperiodic flows, the definition of the operator remains meaningful for more complex flows with continuous or mixed spectra, including turbulent flows. Numerical evidence is provided supporting the close connection between $\mathbf{R}_0$ and $\mathbf{R}_{\overline{\mathbf{U}}}$ in the case of quasiperiodic, chaotic and stochastic two-dimensional incompressible flows.


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 flatplate, the flow may be safely assumed to respond to perturbations in a linear time-invariant (LTI) fashion. Otherwise, the response of open shear flows to external excitations, even infinitesimal, is generally more complex. Turbulent flows certainly are not time-invariant indeed, but neither are idealistic two-dimensional flows like the ones behind a backwardfacing step or in the wake of one or a few cylinders at moderate Reynolds numbers. Even at very low Reynolds number, unsteadiness may set in as a result of noise amplification or intrinsic instability. Hence, infinitesimal extrinsic perturbations do not interact with a steady base flow but a time-evolving one.
However, for feedback control using modern robust control techniques, an LTI model of the input-output dynamics is generally required. In the frequency domain, LTI dynamics indicates the existence of transfer functions, hence the question: how can we define transfer functions for time-varying, statistically steady base flows ?
1.2. Option A: 'dynamic linearity', i.e. using small-amplitude harmonic forcings Some authors have used harmonic forcings to identify transfer functions on simulations of statistically steady flows (Dahan et al. 2012;Dalla Longa et al. 2017;Evstafyeva et al. 2017) which are either laminar or turbulent, and possess either peaked or broadband spectra. For sufficiently low forcing amplitude, the authors evaluate a frequency response by taking the ratio of the complex amplitudes of the output and input at the forcing frequency (the concept of nonlinear transfer function (Noiray et al. 2008) is not relevant here due to the low amplitude of the forcing). This observation hints at the existence of a transfer function in situations where it should not theoretically exist, given the time-varying nature of the base flow. The authors use the term 'dynamic linearity' to characterize this interesting behaviour but no theoretical arguments seem currently available to understand it.

Option B: Linearizing about the mean flow
We adopt in this paper a spatially-discrete framework where the incompressible velocity field is discretized over degrees of freedom. A second option consists in assuming that the linear response u to infinitesimal forcing f is characterized by the Jacobian matrix J U evaluated about the mean flow U, i.e. d u = J U u + f. (1.1) The resolvent operator about the mean flow R U ( ) = ( I − J U ) −1 (1.2) then describes an LTI input-output dynamics in the frequency domain ( being the Laplace variable), and has been successfully used for feedback control of a quasiperiodic flow (Leclercq et al. 2019). 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 a fixed-point of the nonlinear unforced system. The use of such model is however motivated by a large body of the literature devoted to linear analysis about a time-averaged mean flow, with clear predictive power. Modal analysis of J 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, which real part predicts the frequency, is nearly neutral, a property often called RZIF (Turton et al. 2015;Bengana & Tuckerman 2019;Bengana et al. 2019;Bengana & Tuckerman 2021) for Real Zero Imaginary Frequency, in a context that is now broader than just cylinder wake flow. This powerful property has been used to build self-consistent models (Mantič-Lugo et al. 2014, 2015Meliga 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;Stewartson & Stuart 1971), but it is not generally understood and theoretical criteria are currently limited to either the weakly nonlinear limit (Noack et al. 2003;Sipp & Lebedev 2007) or the monochromatic oscillation limit (Turton et al. 2015).
In the case of linearly stable flows, or noise-amplifiers, optimal response modes of R U (i ) often match with coherent structures identified by spectral POD modes (Beneddine et al. 2016;Towne et al. 2018) at the frequency , and optimal forcing modes have been analyzed to better understand receptivity mechanisms (Garnaud et al. 2013;Sartor et al. 2015;Luhar et al. 2015;Semeraro et al. 2016). The intuition that coherent structures must extract their energy from the mean flow dates back to the contributions of Butler & Farrell (1993)

Option C: Averaging the time-varying linear response
In this paper, we will not make the ad hoc hypothesis of linearity about the mean flow. We will only keep the assumption of linearity, which remains valid as long as the forcing amplitude is low enough, but now the base flow is time-varying. Until §4, the base flow will be assumed to be periodic with period corresponding to a fundamental angular frequency 0 := 2 / . We represent in figure 1, the response u to an infinitesimal forcing signal f started at time 0 relative to an arbitrary origin of time = 0 on the base flow trajectory. Because the base flow is not steady but periodic, the linear response to forcing is not unique but parametrized by the phase 0 < 2 of the limit cycle when the forcing starts. This key parameter for our analysis is defined such that 0 0 = mod 2 . (1.3) Indeed, the input-output dynamics is now characterized by a linear time-periodic (LTP) system d u = J ( ; )u + f, J ( + ; ) = J ( ; ).
(1.4) with a periodic jacobian itself parametrized by the phase through the base flow U( ; ) = U( − 0 ). For uniformly distributed in [0, 2 ), there is in fact a distribution of timeseries u( ; ) corresponding to a given input f ( ). The statistically optimal LTI approximation of this LTP system should predict a single output v( ) minimizing the mean error with respect to u( ; ), i.e.
where . denotes an ensemble average with respect to . This output is obviously the mean response v( ) = u( ) . We define the mean resolvent operator R 0 (the notation will become , the responses are parametrized by the phase = of the periodic base flow at = 0, such that 0 0 = mod 2 . (c) The response to linear forcing being a distribution, the statistically optimal LTI model is obtained by computing the mean response u( ) . clear in §3.2) as the operator predicting the mean response, in the frequency domain, for any given input mean resolvent R 0 : f ( ) ↦ → u( ) . (1.6) This transfer function is well-defined because the relationship between the input f and the output u is by construction causal, linear and time-invariant (i.e. independent of ). 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.
Focus on Fluids articles must not exceed this page length 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 §2 and §3, we will continue to focus on the case of incompressible periodic base flows. In §2 we will perform a numerical experiment on the case of the fluidic pinball in a periodic regime of oscillations at = 100, in order to compare options B and C, namely mean resolvent versus resolvent about the mean flow. In §3, we will study the properties of the mean resolvent operator, with the goal to elucidate the observations made in §2 and draw connections with the resolvent operator about the mean flow and 'dynamic linearity'. Connections with the Koopman operator ( §3.3), the RZIF property ( §3.6) and model reduction ( §3.5) will also be made. In §4 we perform the same numerical experiment as in §2 but for more complex base flows, i.e. quasiperiodic, chaotic and stochastic. We will consider the fluidic pinball at = 110 and = 120 as well as the backward-facing step flow (with exogenous stochastic forcing) at = 500. The goal is to stress the strong connection between mean resolvent and resolvent about the mean flow in all these dynamically distinct cases. No computations are done in a compressible case but in §4.3 we indicate the implications of our theoretical analysis in §3 to such situations. A theoretical extension of §3 is proposed for the quasiperiodic case in appendix D. We summarize our findings in §5.

Resolvent about the mean flow versus mean resolvent: fluidic pinball at
= 100 In this section, we perform a numerical experiment on the case of the incompressible twodimensional fluidic pinball in a periodic regime at = 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 of each cylinder and the (convective) time unit to be / ∞ , where ∞ is the upstream velocity norm. The Reynolds number is defined as := ∞ / where is the kinematic viscosity. The centers of the 'front', 'top' and 'bottom' triangles are respectively located at ( , ) = (−1.5 cos( /6), 0), ( , ) = (0, 0.75) and ( , ) = (0, −0.75), hence forming an equilateral triangle of side length 1.5. Instead of measuring the full input-output relation between -dimensional inputs and outputs, we will focus in section §2 on the single-input single-output (SISO) transfer between an actuator signal and a measurement such that The SISO viewpoint is more practical for the numerical experiment we perform, and it is also relevant in the context of flow control. More specifically, we will consider three probes C for = 1, 2, 3 extracting the -velocity component in the wake of the 'top' cylinder at = and respectively = 2, 4, 8 (white filled circles in figure 3). Measurements of the periodic 6 base flow will be denoted = C U. When necessary, the notations and will be used to specify which sensor C we are referring to. The quantity B is a discretized Gaussian volume force field B ( , ; 0 , 0 , , ) located at the 'top' of the 'top' cylinder (2.2) with = 0.01, = 0.1 and the center ( 0 , 0 ) of the Gaussian (white triangle in figure 3) relative to the center ( , ) of the cylinder being given by ( 0 − , 0 − ) = (0, 0.55).
Numerical discretization details (computational domain, mesh, spatial and temporal schemes) are provided in §A.

Unsteady base flow
In figure 3(a), we plot the streamlines and isocontours of velocity norm for the mean flow. We notice it is asymetric, as noted by Deng et al. (2020). In figure (b), we plot the phase portrait of the flow, using the vertical velocity probes 2 and 3: it is a closed curve, which confirms the periodic nature of the dynamics at = 100. Fourier spectra are shown in figure 4(i) for the 3 sensors located in the wake of the cylinders, with (a) = 1, (b) = 2, (c) = 3. These spectra confirm the periodic nature of the flow since they are discrete with a single fundamental frequency at 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.

Estimating mean transfer functions from input-output data
The resolvent operator about the mean flow R U and the mean resolvent operator R 0 are both multiple-input-multiple-output (MIMO) transfer functions with -dimensional inputs and outputs. SISO transfer functions between any time-invariant actuator-sensor pair (C, B) may be derived from these operators using In order to estimate the mean transfer function ( ) from input-output data, we introduce 'frequency response realizations' which are obtained by taking the ratio between the Laplace transforms of the input and the output . The quantity is not directly associated with a transfer function since the linear response to forcing about a time-varying base flow is not time-invariant; hence the dependence on the phase and the input signal . However, upon averaging with respect to (we recall that is assumed to be uniformly distributed in [0, 2 )), 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) For any complex frequency and input signal , the variance Var ( ; ) of the complex random variable ( ; , ) may be computed as 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 if 1 for any choice of . The ratio is also useful to evaluate the convergence of the mean estimate , to the true mean , using 1 realizations of . Since all realizations are independent and identically distributed random variables of variance Var , the mean estimate , is, by the central limit theorem, a normal random variable of variance Var . (2.8) The convergence criterion √︁ Var , | | for the mean estimate therefore translates into the condition √ for the given choice of and . In the following, for simplicity, we will use the notation to refer to the mean estimate, whether it is converged or not (which will be clear by comparing with √ ).

Procedure
The nonlinear simulation is initialized at = 0 with the steady base flow (fixed-point) perturbed by 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 ( ) = ( ) by simply initializing the perturbation field with the actuator, i.e. u(0) = 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 baseflow) and linear (for the perturbation) simulations are run for over more than 1000 convective time units. We compute = 140 impulse responses, by uniformly distributing the relative initial time 0 over a period = 8.64 of the underlying limit-cycle, i.e. Δ 0 = / . This amounts to uniformly distributing the initial phase over [0, 2 ), as required.
Note that since impulse responses never decay, we can only compute frequency responses shifted in the right-half plane, i.e. obtained for = i + with > 0. This is indeed mandatory for convergence of the Laplace transform. The value of = 0.01 is chosen such that the mean impulse response modulated by e − 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

Results
In figure 4, the three columns correspond to the estimated mean frequency responses (solid blue lines) from to the three sensors : (a) = 1, (b) = 2, (c) = 3. Comparison with the transfer function about the mean flow U is provided (dashed red curves). Row (ii) corresponds to the gain and the next row (iii) to the phase, while the last row (iv) corresponds to the ratio defined in (2.7), for the particular case of impulsive forcings ( ) = ( ). The two previously identified thresholds of 1 and √ are respectively indicated with a solid-black line and a dashed green line. Vertical magenta lines highlight the first natural frequencies = 0 of the flow for = 1, . . . , 6. The first remark is that the standard deviation of varies a lot with frequency and sensor position, which impacts both the validity of the LTI assumption and the convergence of the sample mean in the same way. The LTI assumption is valid (i.e. blue curve well-below black dashed line) in the near-wake region, close to both the actuator and the separation point, about the frequency ≈ 1 of maximal amplification. It also appears to be generally valid at high enough frequency, even far from the actuator, but not at intermediate frequencies.
At low frequency and in the close vicinity of natural frequencies, the validity is poor. The quality of the LTI assumption deteriorates in the downstream direction, as fluctuations grow in amplitude.
Within the range of acceptable convergence of the mean (i.e. blue curve well-below green dashed line), we notice that for each sensor position, and U closely match both in gain and phase: the transfer function about the mean flow appears to be an excellent approximation of the mean transfer function, regardless of the flow regime. However, appears to possess resonance peaks at multiple natural frequencies (see vertical magenta lines), which are not visible in U . These extra resonant peaks at high frequencies are only visible for downstream sensors and are indicative of the presence of marginally stable poles at = i 0 .
In the theoretical section, we will attempt to elucidate the following observations: i) why are and U generally so similar, ii) why does the quality of the LTI approximation deteriote downstream, iii) why are there poles at = i 0 in the mean transfer function and iv) why are resonances at high frequencies only noticeable downstream ?

Theory for incompressible periodic base flows
In this theoretical section, we come back to the more general MIMO viewpoint between a -dimensional forcing f ( ) and a -dimensional response u( ; ).
3.1. Phase-dependent frequency response of LTP system Denote ( , ; ) the propagator from to of the linear time-periodic system (1.4), such that where u 0 is the initial condition on u, regardless of . By Floquet's theorem (Magruder et al. 2018), we know that there exists a real -periodic matrix P ( ; ), invertible at all times , and a constant complex matrix K such that ( , ; ) = P ( ; )e K ( − ) P −1 ( ; ). There are multiple possible determinations of the matrix K = (1/ )log( ( , 0)) depending on the definition of the logarithm. Assuming K to be diagonalizable K = S S −1 and choosing the principal determination for the logarithm, all the eigenvalues fall into the fundamental strip − 0 /2 < Im( ) 0 /2, = 1, . . . , . These are called the principal Floquet exponents and they may be ranked in decreasing order of growth rate, i.e. Re( 1 ) Re( 2 ) · · · Re( ). Other determinations of the logarithm lead to eigenvalues = + i 0 which are Floquet exponents as well, located in complementary strips 0 ( − 1/2) < Im( ) 0 ( + 1/2) in the complex plane. The direct Floquet modes v are the columns of the -periodic complex matrix V ( ; ) := P ( ; )S. The adjoint Floquet modes w , such that (w , v ) = for the canonical inner product (a, b) = 1/ ∫ 0 a b d on the space of periodic functions in C (where denotes the complex transpose), are the columns of the matrix W ( ; ) such that W V = I at all times (and phases ). Using this eigendecomposition, the propagator reads We may now expand the direct and adjoint Floquet modes as Fourier series. If the th Fourier coefficient of V is proportional to e i , then so is the th Fourier coefficient of W since W V = I at all times. Therefore, we can write the following Fourier series Injecting (3.2) and (3.3) in (3.1) and assuming u 0 = 0, we have 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 with the operators where v and w designate the th column of V and W . Unlike in LTI systems, the output at frequency depends on the input at an infinite number of frequencies − i 0 . The output depends explicitely on the phase through the cross-frequency transfers e i R ( ) for ≠ 0.
The poles of the R are exactly the Floquet exponents of the LTP system. The adjoint Floquet modes characterize receptivity to forcing of the corresponding direct Floquet modes.

Mean resolvent operator
We now seek the operator which predicts the mean output u( ) for a given input f ( ), by averaging (3.8) with respect to . By doing so, all the cross-frequency transfers vanish and the only remaining transfer is R 0 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 U( ) is self-sustained (no external forcing is required), there is a zero Floquet exponent in the LTP system with associated Floquet mode equal to d 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. 1 = 0 and v 1 = d U, so that v 1 = i 0 U where U denotes the th harmonic of the Fourier decomposition of U( ). Therefore, using expression (3.9), the mean resolvent reads We immediately see that there are purely imaginary poles at = i 0 as expected from the numerical experiment on the fluidic pinball. The coefficients associated to these poles in the expansion, also called the residuals tend very rapidly to 0 in the operator norm as | | → ∞, because the th Fourier harmonic of both U and w 1 are involved (the decay is exponential for C ∞ periodic structures). This explains why resonance peaks caused by imaginary poles in figure 4(ii), become decreasingly visible as | | increases. The fact that the residuals depend on U also explains why resonances at high frequencies are however more visible downstream than upstream. Indeed, larger oscillation amplitude in the far wake causes more pronounced nonlinear energy transfers to higher harmonics (see base flow spectra in row (i) of figure 4), hence a better observability |C U | of high-order harmonics U as the sensor moves downstream.

System identification using frequency-rich inputs
The upside of using frequency-rich signals for system identification is that it allows identifying the dynamics for multiple frequencies at once. The downside though, is that many realizations of the same input signal may be necessary to converge the mean response prior to identification, as we have seen in §2 (in particular row (iv) in figure 4). 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.8) and the definition of the variance: where a N = (a Na) 1/2 is the norm induced by the positive-definite matrix N. A convenient choice of input signal is one which minimizes the variance, allowing for a minimal amount of realizations. Clearly, an input signal with a broadband spectrum like an impulse f ( ) = B ( ) (for which ( ) = 1 for all ) is likely to generate a lot of output variance and is not necessarily a favourable choice for identification. An alternative to frequency-rich inputs is to use non-resonant harmonic forcings for identification, as we shall see in the next section.
Another advantage of using a broadband input signal for identification is that if the ratio 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.13) also explains why the quality of the LTI approximation deteriorates as the sensor moves downstream in figure 4(iv) for ( ) = ( ). Indeed, which, according to (3.9), increases downstream since the observability |C v | of the Floquet modes probably grows downstream (this is obvious for the first Floquet mode since v 1 = i 0 U ).

System identification using harmonic inputs: 'dynamic linearity'
The Laplace transform of harmonic forcings of the form f ( ) = Be i is given by Plugging this expression into the input-output relation (3.8) leads to The temporal response may now be evaluated by inverting the Laplace transform. The decaying Floquet exponents lead to a transient response while all the purely imaginary poles in the right-hand-side of (3.16), either in the transfer operators or the forcing, lead to a permanent contribution. For non-resonant forcing frequencies ≠ 0 , we have The first line (3.17) collects the transient contributions, the second line (3.18) gathers terms which are phase-independent, while the third line (3.19) is phase-dependent and vanishes upon averaging. However, the point here is not to take an average but notice that the Fourier coefficient at the forcing frequency, which may be obtained using an harmonic average is phase-independent and exactly given by the product of the input amplitude f ( ) = B by the mean resolvent operator evaluated at = i : In other terms, using small-amplitude harmonic forcings, it is possible to obtain frequency samples of the mean transfer function = C R 0 B without the need to carry out an ensemble average over several realizations. This is in essence the 'dynamic linearity' phenomenon described by (Dahan et al. 2012;Dalla Longa et al. 2017;Evstafyeva et al. 2017) and used therein for LTI system identification on time-varying base flows.
In principle, using a superposition of harmonic forcings ( 3.22) 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 − ≠ 0 and ≠ 0 for any , > 0.

Connection with the Koopman operator
The propagator from = 0 to Averaging the propagator with respect to leads to the mean propagator By comparing equations (3.24) and (3.9) for = 0, we see that the mean resolvent operator is the frequency domain representation of the mean propagator from 0 to , i.e. of the mean Koopman operator for the full-state observable of the LTP system . . .
. . . ( 3.26) The block H of this infinite matrix operator characterizes the transfer from the frequency + i 0 of the input f to the frequency + i 0 of the output u (Wereley & Hall 1990, 1991Zhou & 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.27)

The harmonic transfer operator as a feedback loop
Decompose the Jacobian as a mean J and a periodic perturbation J ( ; ), which we expand as a Fourier series . (3.28) Because J ( ; ) is real, the harmonics have hermitian symmetry, i.e. J − = J * , where (.) * denotes the complex conjugate (not the conjugate transpose (.) ). Moreover, we also obviously have J 0 = 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.29) the operator may be interpreted as a feedback loop between the two blocks D −1 ( ) and T ( ), as illustrated in figure 5. 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 . 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 D −1 ( ) 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 . 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 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 shall see next. This expansion may be interpreted in the following way: H 00 collects contributions of H which deliver an output at the same frequency than the input, after interacting times with the fluctuating part of the base flow. Using the general terms we find . . .
The order 0 term is equal to the resolvent operator about the mean flow, which corresponds to interacting 0 times with the fluctuating Jacobian. The order 1 term is exactly 0, because it is impossible to interact just once with the fluctuating Jacobian and come out of the loop at the same frequency than the input. The first correction to the resolvent operator about the mean flow then arises at order 2, when successive interactions with J and its complex conjugate J * = J − occur, allowing the output to be at the same frequency than the input.
The Neumann expansion of the mean resolvent may be used to bound the relative ∞, norm of the error between the mean resolvent and the resolvent operator about the mean flow, where R 0 ∞, := sup where . 2 denotes the spectral norm (maximum singular value), not the 2 norm. In the process, we will need the Wiener norm (Gohberg et al. 2003) of the matrix-valued periodic function J which is obviously finite given the exponential convergence of the Fourier coefficients of J . Using the two definitions above, we find The coefficient characterizes the loop gain associated with the feedback interconnection of D −1 ( ) and T ( ), averaged with respect to . In the weakly unsteady limit → 0, we find The small parameter depends on because the shifted imaginary axis must lay far enough to the right of the poles of both R 0 and R U (see appendix C). For a given acceptable (so that the Neumann series converges), the relative error vanishes as J 2 W , so very quickly as unsteadiness weakens. The agreement between transfer functions based on R 0 and R U in the fluidic pinball §2 (see figure 4 rows (i,ii) in particular) indicates that amplification by the unsteady part of the periodic base flow is quite weak in this open shear flow.
3.5. Reduced-order models of the mean resolvent Another way to write the mean resolvent is in the form (3.50) I = I , (3.51) = ( + i 0 I) . (3.52) The input (resp. output) matrix B (resp. C) has columns (resp. lines) and an infinite number of rows (resp. columns), while the diagonal state matrix is infinite-dimensional. Expression (3.48) is associated to a state-space representation where the internal state x = [. . . , x −1 , x 0 , x 1 , . . . ] is an infinite column vector, even though the input f and the output u both have a finite dimension . However, for large enough | | the poles + i 0 in become both uncontrollable and unobservable as the Fourier coefficients W in B and V in C vanish exponentially. Since there is a finite number of effectively observable and controllable poles, it should be possible to approximate the input-output dynamics by projecting the infinite-dimensional state vector x onto a finitedimensional subspace.
The resolvent operator about the mean flow approximates R 0 and it only has poles: it may therefore be interpreted as a reduced-order model of order of the mean resolvent. Even though the approximation is good for weakly unsteady flows, the operator R U does not take into account the fluctuating part of the Jacobian, hence the model may not be optimal for the fixed order . Moreover, there is no reason a priori to choose a model order equal to : 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 depending on the desired precision.

The 'RZIF' property
We recall that 1 = 0 (see §3.2.1 and appendix B) therefore = i 0 are Koopman eigenvalues of the mean propagator satisfying exactly the so-called RZIF property originally discussed in the context of mean flow stability analysis (Turton et al. 2015). Since v 1 = d U, the associated Koopman modes v 1 = i 0 U are parallel to the Fourier components U 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 (Turton et al. 2015) or the dynamics is weakly non-linear (Noack et al. 2003;Sipp & Lebedev 2007). Therefore, the proposed framework appears to be more generic than modal analyses of mean flows. We note that a weakly time-varying criterion qualitatively similar to the one proposed in the previous section was already provided to justify the RZIF property in Mezić (2013) (criterion (a) §3.5 therein).
Finally, we add that the RZIF property is also be verified for unstable periodic base flows, but in this case the zero Floquet exponent does not correspond to the leading one, i.e. 1 ≠ 0, which has a strictly positive growth rate (see appendix B).

Towards more complex base-flows
The previous sections were only concerned with periodic base flows, but the definition (1.6) of the mean resolvent R 0 may be easily generalized to any other statistically steady base flow.
Instead of taking the mean output with respect to the phase , we may take more generally the mean output with respect to the relative time 0 at which the forcing signal is released: For such a definition to hold, we assume that the single trajectory U( ) of the unsteady base flow, about which we linearize the dynamics, covers the entire dynamical attractor. Otherwise, an average over multiple base flow realizations also needs to be performed. We start by performing similar numerical experiments as in §2, but for more complex incompressible two-dimensional base flows. The numerical procedure explained in §2.3 is repeated, but now varying the time 0 at which the impulse is triggered, rather than the phase , which is no longer defined. The definitions introduced in §2.2 are extended by replacing . with . 0 , since now 'frequency response realizations' are parametrized by 0 . In §4.1, we consider the fluidic pinball in the quasiperiodic and chaotic regimes by increasing the Reynolds number to respectively = 110 and = 120 (Deng et al. 2020). For the quasiperiodic case, an extension of the periodic theory is possible and reported in appendix D. Next, in §4.2, we consider the stochastic flow past a backward-facing step at = 500, as in Hervé et al. (2012). Two values are considered for the variance 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 In the periodic case, it is equivalent to take the ensemble average with respect to 0 or and the MLE corresponds to Re( 1 ) = 0. In the quasiperiodic case, we also have MLE = 0, for the chaotic regime MLE > 0, while for the stochastic regime MLE < 0 (results of these computations are reported in appendix E). The goal is to provide numerical evidence of the strong connection between the mean resolvent and the resolvent about the mean flow in all these cases, even though our theory only applies to the periodic and quasi-periodic cases. Finally, in §4.3, we discuss implications of our previous theoretical analysis to the case of compressible flows.

Quasi-periodic and chaotic flows: fluidic pinball at
= 110 and = 120 The numerical setup is identical to that of §2 (with more details in appendix A) for = 100.

Unsteady base flow
In figure 6(a), we plot the streamlines and velocity norm isocontours for the mean flow in the (i) quasiperiodic case and the (ii) chaotic case. We notice that the mean flow is asymetric in the quasi-periodic case, but symetric in the chaotic case. The recirculation zone is more extended downstream in the latter case. In figure 6(b), we plot the phase portrait of the flow, using the vertical velocity probes 2 and 3. The amplitude of the fluctuations increases with the Reynolds number. The phase portrait is well-structured in the quasiperiodic case but it is impossible to infer the dimensionality of the torus from this two-dimensional projection of the attractor. The complex structure of the phase portrait in figure 6(bii) is typical of chaotic dynamics.
Fourier spectra are shown in row (i)  the spectrum is discrete in the quasiperiodic case, while it is continuous in the chaotic case, despite the apparent presence of peaks. In the latter case, they possess an intrinsic non-zero bandwidth, unlike in the periodic and quasi-periodic cases, where the 'thickness' of the peak only arises from the estimation error of the discrete Fourier transform. It is not directly possible to infer the number of basic incommensurate frequencies in the quasiperiodic case by examining the spectra, but there clearly is energy in a much greater number of frequencies compared to the periodic case (see figure 4(i) for comparison).

Procedure
In both the quasi-periodic and chaotic cases, = 204 impulse responses are computed over more than 1000 convective time units, each starting 2 convective time units after the preceeding one, i.e. Δ 0 = 2. We assumed that this procedure allowed to sample a representative part of these attractors with a single realization of U( ). We recall that > 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 = 0.01 since MLE = 0. For the chaotic case, we choose = 0.03 since MLE = 0.02 (see appendix E).

Results
Results are reported respectively in figures 7 and 8, and they are qualitatively very similar to the results of the periodic case in figure 4. The quality of the LTI assumption deteriorates as the sensor moves downstream (blue curve more often above black threshold in row (iv)). As a consequence, the convergence of the estimate of the mean transfer function is slower for downstream sensors (blue curve above green threshold). The LTI approximation is however always good at large enough frequencies and poor at low enough frequencies. At intermediate frequencies in the vicinity of the maximum gain, the LTI approximation is reasonably good ( ≈ 1) but degrades locally in the close vicinity of energetic natural frequencies of the base flow. Overall, agreement between 0 and U is excellent on a large frequency range in the near wake, and poorer in the far wake.
In the quasiperiodic case, we notice the presence of resonance peaks in the mean transfer functions, at natural frequencies of the unsteady base flow, which are not visible in U . Resonances at high frequencies are not visible in the near wake and only manifest at sensors 2 and 3 in the far wake. All the observations made for the quasi-periodic case are similar to the periodic case and the same justification may be provided in both cases (see appendix D for an extension of the theory developed in §3 to the quasiperiodic case).
Interestingly, the close agreement between 0 and U carries over to the chaotic regime, even though we do not have a theory for base flows with a continuous power spectrum. with the fluidic pinball is that the 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 F( ) = B ( ) with a Gaussian-white noise and B a discretized Gaussian volume force B ( , ; 0 , 0 , , ) with ( 0 , 0 ) = (−0.5, 0.25) and = = 0.1 (see equation (2.2) for the definition of B and white squares centered at ( 0 , 0 ) in figures 9 (a,b)). Two noise variances are considered: = √ 10 and = 10, respectively representing weakly and strongly nonlinear dynamics as in Hervé et al. (2012).
Note that the forcing F acts on the nonlinear base flow U and should not be confused with f = B , which acts on the linear flow u. The actuator B is a discretized version of B with ( 0 , 0 ) = (−0.05, 0.01) and ( , ) = (0.01, 0.1) (see white triangles in figures  9 (a,b)). We again use three probes to monitor the flow response: two -velocity probes in the shear layer at respectively (2.5, 0) and (5, 0), and a friction sensor at 12.5 12.9 on the lower wall, behind reattachment (see white dots and rectangle in figures 9(a,b)). The perturbation and base flow measurements will also be denoted = C u and = C U with index = 1, 2, 3 (the number increases downstream) when necessary.

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

Procedure
Again, we use a single realization of U( ) corresponding to a single random timeseries ( ) but vary the relative time 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 Δ 0 between two impulses is of 2 convective time units, as in the quasiperiodic and chaotic flows. We run = 308 linearized impulse responses, over 220 time units only. Indeed, while impulse responses never decay for the fluidic pinball, since MLE 0, they do for the backward-facing step flow since MLE < 0 (see appendix E) and reach negligible amplitude over that period of time, making shorter  Figure 10: Same caption as figure 4 for the backward-facing step flow at = √ 10 (low-noise case), with = 308 samples. Unlike for the fluidic pinball, the frequency response is on the imaginary axis, with no shift since 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.
integration possible. Since MLE < 0, the Laplace transform may be evaluated directly on the imaginary axis.

Results
Results for the stochastic case are presented in figures 10 and 11, respectively corresponding to low = √ 10 and high = 10 noise variance. Interestingly, the results are qualitatively very similar to the fluidic pinball, even though the nature of the dynamics is quite different. Unlike for the fluidic pinball though, the LTI assumption and convergence of the mean transfers appear to be poor at high frequency, for downstream sensors 2 and 3 (see row (iv)). Otherwise, whenever the convergence of the mean estimate is satisfactory, the agreement between and U is excellent (see rows (ii) and (iii)). There are no longer sharp variations in the gain and phase at very specific frequencies, which seems in accordance with the absence  Figure 11: Same caption as figure 10 for the high-noise case = 10, with = 308 samples.
of resonance peak in the base flow spectrum (see row (i)). For both values of and all three sensors, best 'signal-to-noise' ratio is achieved around ≈ 1, which also corresponds to the maximum gain. This is consistent with the flow preferentially amplifying the incoming white noise in this frequency band. Increasing the noise variance has an effect similar to moving the sensor downstream: the 'signal-to-noise' ratio strongly deteriorates as fluctuations grow in amplitude (see row (iv)). For = 1 and = √ 10, the mean is converged for all frequencies 10 −2 < < 10 2 whereas for = 3 and = 10, convergence is barely achieved at any point.

Strongly compressible flows
It was recently pointed out by Karban et al. (2020) that resolvent analysis about the mean flow is ill-posed in the case of strongly compressible flows, where large discrepancies in the optimal gains are observed depending on the formulation of the Navier-stokes equations, whether in primitive or conservative variables. The inconsistency was attributed to the fact that the mean flow eigenvalues are not conserved under a nonlinear coordinate change.
Numerical experiments have not been carried out here for compressible cases, but we wish to highlight the potential of the mean resolvent to resolve this issue, since we showed that the poles of the operator correspond to Floquet exponents (and Koopman eigenvalues) of the system in the LTP case (see §3.1 and §3.3), which are invariant under a nonlinear coordinate change (see appendix F).
Another point we want to stress is the assumption (3.29) we made early on in §3.4 that the mean Jacobian is equal to the Jacobian operator about the mean flow. This is no longer true for cubic nonlinearities found in the compressible Navier-Stokes equations written either in primitive or conservative formulation. However, we may still use expansion (3.39)-(3.42) if we replace R U with R J ( ) = ( I − J) −1 , (4.3) which we may call the resolvent operator about the mean Jacobian. Alternatively, a quadratic formulation of the compressible Navier-Stokes equations may be used (Vigo 1998;Iollo et al. 2000), even though this form may not easily handle shock discontinuities due to its non-conservative form.

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 time-march the non-autonomous linear equations about a trajectory of the statistically steady flow, and only then average the results over an ensemble of input realizations triggered at various time instants along the base flow trajectory. By considering two configurations, the fluidic pinball and the backward-facing step flow, we investigated 4 possible dynamical regimes: periodic, quasi-periodic, chaotic and stochastic. Multiple frequency response realizations were obtained by varying the initial time 0 (or phase in the periodic case) of the forcing relative to an unsteady base flow trajectory. This change of perspective lead to interesting findings.
First of all, our framework allows to quantify the validity of the LTI assumption by analyzing the standard deviation of the distribution of frequency response. The LTI approximation is found to be satisfactory in the upstream region of all flows for the frequency range of maximal amplification by the mean flow. The assumption deteriorates as perturbations grow in amplitude, either as the sensor moves downstream or as noise variance increases. In all cases considered though, we found that the transfer function based on R U always approximates to some degree that based on the mean resolvent R 0 , but the two objects are not identical. In the periodic (and quasi-periodic) case(s), extra resonance peaks at multiples of the fundamental frequency(ies) are visible in the mean transfer function.
We were able to explain these observations in the periodic (and quasiperiodic) case(s), by averaging the harmonic transfer operator with respect to (or to the relative initial time of forcing 0 in the quasiperiodic case). The harmonic transfer operator was introduced by Wereley & Hall (1990Hall ( , 1991 and recently brought back to the fore in our community by Padovan et al. (2020);Franceschini et al. (2022) as the harmonic resolvent operator. In the periodic case, it characterizes cross-frequency transfers from any frequency + i 0 at the input, to any frequency + i 0 at the output (a similar definition may be proposed in the quasi-periodic case). We showed that cross-frequency transfers cancel out on average, such that only the transfers for = remain. The remaining LTI operators are all identical and correspond to the same mean resolvent operator evaluated at different frequencies + i 0 . Mean transfer functions, based on this mean resolvent operator, may be identified from input-output data without the need for averaging over several input realizations if harmonic forcings are used, even though the linearized system is not time-invariant. This justifies the principle of 'dynamic linearity' put forward by Dahan et al. (2012);Dalla Longa et al. (2017);Evstafyeva et al. (2017) in the context of flow control. Harmonic identification directly leads to the statistically optimal LTI input-output models for periodic (and quasiperiodic) flows.
Using (quasi-)Floquet theory for the (quasi-)periodic case, we showed that the poles of the mean resolvent operator are exactly the (quasi-)Floquet exponents of the LTP (LTQP) system. 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. Since the mean resolvent is a frequency representation of the mean Koopman operator for the full-state observable of the linear dynamics about the periodic base flow, the poles of R 0 may also be interpreted as Koopman eigenvalues. Among those poles, we find the set i 0 of marginally stable poles at multiples of the fundamental frequency, which causes resonance peaks in the gain of mean transfer functions. The Koopman modes associated with these eigenvalues are parallel to the Fourier coefficient of the periodic base flow U . This is reminiscent of modal analysis of J U for periodic flows, where there is often a nearly marginally stable eigenvalue with frequency almost equal to 0 , and where the corresponding eigenvector is nearly parallel to U 1 . This so-called 'RZIF' property which is only approximate in the case of modal analysis about a mean flow, appears to be exact when considering the Koopman eigenvalues and modes associated with R 0 .
Next, we investigated the connection between the usual resolvent about the mean flow R U , and the mean resolvent R 0 . We found that the former approximates the latter, and that the difference is second order based on a small parameter characterizing the strength of the linear amplification by the unsteady part of the (quasi-)periodic Jacobian. The main difference between the two resolvent operators is that R 0 has an infinite-dimensional internal state, while R U has a finite-dimensional internal state of dimension . Infinite-dimensionality arises from the infinite set of poles corresponding to the (quasi-)Floquet exponents. However, the observability and controllability of these poles + i 0 quickly vanish for large | | explaining why a finite-dimensional (in terms of internal state) approximation of R 0 is possible. The observability of the marginally stable poles at i 0 by a sensor field C is characterized by |C U |, which explains why resonances at higher-order harmonics of 0 are more visible in the far wake than in the near wake of the periodic fluidic pinball. The higher observability of these higher-order harmonics in the far wake also explains the greater variance of the transfer function distribution for downstream sensors.
In our analysis, we assumed a quadratic nonlinearity in the governing equations implying J = J U . While this is true for the incompressible Navier-Stokes equations in primitive variables, it is not for the compressible equations in either primitive or conservative variables. In these cases, the nonlinearity is cubic and no connection can be made between the mean resolvent R 0 and R U . Instead, it is the resolvent operator about the mean Jacobian R J which approximates R 0 . This may help elucidate the ambiguity in resolvent analysis recently noted by Karban et al. (2020) in resolvent analysis about a mean flow, as discrepancies are particularly noticeable in strongly compressible flows. In this context, a compelling prospect would be to compute the optimal gains, forcing and response modes of R 0 and compare them to the corresponding quantities for R U .
An interesting perspective would be to extend the present theory to the case of continuous or mixed spectra characteristic of chaotic, stochastic and turbulent flows. Indeed, numerical results in the two-dimensional chaotic and stochastic cases seem to indicate an excellent match between the mean resolvent operator and the resolvent about the mean flow, which remains unexplained. Koopman operator theory may be key to this endeavour (Črnjarić-Žic et al.

28
2019). Another interesting question would be to find an optimal way to compress the infinite state of the mean resolvent operator into a finite one. Since the mean resolvent is related to the mean Koopman operator, this problem may perhaps be tackled using extensions of the DMD method ( = 6 and 'bottom' = −6 boundary conditions. A standard outflow boundary condition n − −1 ∇ U · n = 0 is used at = 20, which is similar to the stress-free boundary condition of Deng et al. (2020). The boundary condition on the cylinders is obviously no-slip and impermeability, i.e. U = 0. The freely available software FreeFem++ (Hecht 2012) is used to time-march the incompressible Navier-Stokes equations discretized on Taylor-Hood finite elements; 2 for velocity components and forcings, and 1 for pressure. The mesh comprises 66, 668 triangles, 33, 739 vertices and the total number of velocity degrees of freedom, including both components, is = 268, 296.
In order to validate spatial resolution, we produce 4 mesh 1 − 4 of variable refinement, and compute the lift and drag coefficients for the asymetric steady base flow at = 75 and = 100. The coefficients are based on the total force by unit length exerted on all three cylinders normalized by 1/2 2 ∞ , where denotes density. The results are given in table 1 and justify our choice of mesh M2 for the present study. Note that the values provided here do not match with those of Deng et al. (2020) as the authors did not include the viscous contribution to the forces (private communication). For completeness, we provide in figure 12 the temporal evolution of the lift coefficient in the periodic regimes at = 68, 75, 100 analyzed by Deng et al. (2020), using the usual definition of the coefficients. This may be useful for validation purposes in future studies on this configuration.
The solvers are based on a sequential open-source code without time-splitting (https : //github.com/denissipp/AMR_Sipp_Schmid_2016; Sipp & Schmid (2016)). Code parallelization was found unnecessary for this study, which requires multiple realizations of impulsive forcings which can be run in parallel indeed. The nonlinear problem is solved in perturbative form with respect to the time-invariant base flow U (i.e. fixed-point), after computing the latter using a Newton method. The time-stepping is semi-implicit: the diffusion term and interaction of the perturbation U := U − U with U are implicit, while the interaction of U with itself is explicit (Adams-Bashforth). The linear non-autonomous problem about the time-varying base flow is also solved semi-implicitely: the diffusion term and interaction of the perturbation u with U are implicit, while the interaction of u with U is explicit (Adams-Bashforth as well). The time-stepping scheme is second-order for both problems, however we found the latter to be more sensitive to step size than the former. Temporal resolution tests were carried out for linear impulses about the quasiperiodic flow at = 110. Figure 13 shows   (c) the impulse response about the time-averaged mean flow. The corresponding 'frequency response realization' (shifted in the right-half plane, as specified in §2.3) was obtained by taking the Laplace transform of the impulse response and then compared with frequential samples computed using the resolvent operator about the mean flow (on the same shifted axis).

Appendix B. Marginal Floquet exponents for self-sustained periodic base flow
We find it useful to include a proof of this classical result. The periodic base flow U is a solution of the unforced discretized incompressible Navier-Stokes equations projected onto the space of solenoidal velocity fields d U = F(U), (B 1) and by differentiating with respect to time Since d U is a solution of the unforced LTP system, its evolution from to is governed by the propagator Clearly, d U( ) is a Floquet mode associated with the fundamental Floquet exponent 0 (and copies i 0 in complementary strips). Since the periodic base flow is stable and Floquet exponents are ranked in decreasing order of growth rate, we have 1 = 0 and v 1 = d U. Note that if the base flow is an unstable periodic orbit of the unforced system, then there is still a zero Floquet exponent, but it will no longer be the leading one since the maximum Lyapunov exponent has to be strictly positive.

Appendix C. Convergence of the expansion of the harmonic transfer operator
Consider the space 2 (Z, C ) of infinite sequences (x ) ∞ =−∞ of -dimensional complex vectors indexed by Z of finite 2-norm, i.e. is strictly less than one, i.e. D −1 ( )T ( ) < 1 (Suzuki 1976). Since the operator norm is always larger than the spectral radius, a sufficient condition is that D −1 ( )T ( ) 2 < 1. Because operator norms are submultiplicative, D −1 ( )T ( ) 2 D −1 ( ) 2 T ( ) 2 so we may bound the two operators separately.
Each block D −1 ( ) = R U ( + i 0 ) on the diagonal of D −1 corresponds to the resolvent operator about the mean flow evaluated at a different frequency on the shifted imaginary axis + i , where = Re( ) and = Im( ). Because of the block-diagonal structure of D −1 , we have D −1 ( ) 2 = max R U ( + i 0 ) 2 R U ∞, with the ∞, norm defined in (3.43).
(C 1) This condition may be interpreted as an expression of the small-gain theorem (Vidyasagar 1978) for the feedback loop between the two blocks D ( ) and T ( ). Now, since R U ∞, → 0 as → ∞, there exist a critical value + > 0 such that the series converges on the entire half-plane > + . In fact, this critical value is necessarily larger than both the maximum Lyapunov exponent (i.e. the maximum growth rate of all Floquet exponents, equal to 0 if the periodic base flow is linearly stable) and the maximum growth rate of the mean flow eigenvalues (which sometimes is strictly positive ( We conclude, as in the periodic case, that the resolvent operator R U about the mean Jacobian approximates the mean resolvent operator R 0 and that corrective terms are of order 2 with respect to J W . Following (Mezić 2020), an equivalent to Floquet's theorem in the quasi-periodic regime, due to Sell (1978), may also be formulated: under appropriate conditions, there exists a quasi-periodic transformation P ( ; 0 ) and a constant matrix K such that the system (1.4) admits a propagator of the form ( , ; 0 ) = P ( ; 0 )e K ( − ) P −1 ( ; 0 ). Direct and adjoint quasi-Floquet modes V ( ; 0 ) and W ( ; 0 ) may be defined in the same way as before, but they are now quasi-periodic instead of periodic, and adjoint with respect to the inner product

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

Appendix F. Invariance of the Floquet exponents to a coordinate change
Consider the bĳective nonlinear coordinate change h on the base flow