2.1. Linear-time-varying dynamical system

We begin with a general dynamical system of the form

(2.1)\begin{equation} \frac{{\rm d} \boldsymbol q(t)}{{\rm d} t} = \mathcal{N}(\boldsymbol q(t)),\end{equation}

where $\boldsymbol q(t)\in \mathbb {R}^N$ is the state vector and $\mathcal {N}:\mathbb {R}^N\rightarrow \mathbb {R}^N$ is a nonlinear function that describes the dynamics of the system. We decompose the state $\boldsymbol q(t)$ into a base state $\bar {\boldsymbol q}(t)$ and a perturbation $\boldsymbol q'(t)$ such that $\boldsymbol q(t)= \bar {\boldsymbol q}(t) + \boldsymbol q'(t)$. We substitute the decomposition into (2.1) and adopt the Taylor series expansion. After neglecting terms that are third order or higher, we obtain

(2.2)\begin{equation} \frac{{\rm d} \bar{\boldsymbol q}(t)}{{\rm d} t} + \frac{{\rm d} \boldsymbol q'(t)}{{\rm d} t} = \mathcal{N}(\bar{\boldsymbol q}(t)) + \left. \frac{\partial \mathcal{N}}{\partial \boldsymbol q} \right\vert_{\bar{\boldsymbol q}(t)} \boldsymbol q'(t) + {O}^2(\boldsymbol q'(t)),\end{equation}

which can be cast as

(2.3)\begin{equation} \frac{{\rm d} \boldsymbol q'(t)}{{\rm d} t} = {\boldsymbol{\mathsf{A}}}(t) \boldsymbol q'(t) + \boldsymbol f'(t),\end{equation}

where ${\boldsymbol{\mathsf{A}}}(t) = \partial \mathcal {N}/ \partial \boldsymbol q \vert _{\bar {\boldsymbol q}(t)}\in \mathbb {R}^{N\times N}$ is the Jacobian operator, and

(2.4)\begin{equation} \boldsymbol f'(t) = \mathcal{N}(\bar{\boldsymbol q}(t)) - \frac{{\rm d} \bar{\boldsymbol q}(t)}{{\rm d} t} + {O}^2(\boldsymbol q'(t)), \end{equation}

which contains the residual terms arising from (2.1) when $\bar {\boldsymbol q}(t)$ is not an exact solution and the second-order terms that are nonlinear in $\boldsymbol q'(t)$. From the general expression of the linearized perturbation dynamics governed by (2.3), we will derive the classical and harmonic resolvent analysis frameworks in §§ 2.2 and 2.3, respectively, depending on the nature of the base state (i.e. stationary or time varying).

2.2. Classical resolvent analysis framework

In the classical resolvent analysis we consider the base state to be a steady solution of (2.1) or a statistically stationary time-averaged mean of the solution. Then (2.3) becomes

(2.5)\begin{equation} \frac{{\rm d} \boldsymbol q'(t)}{{\rm d} t} = {\boldsymbol{\mathsf{A}}}\, \boldsymbol q'(t) + \boldsymbol f'(t),\end{equation}

where ${\boldsymbol{\mathsf{A}}}\in \mathbb {R}^{N\times N}$ is the time-invariant Jacobian operator that governs the dynamics of the unsteady perturbations $\boldsymbol q'(t)$ with an unsteady forcing $\boldsymbol f'(t)$. The forcing $\boldsymbol f'(t)$ has different interpretations in the classical resolvent analysis that are worth discussing. If $\bar {\boldsymbol q}$ is an exact solution of (2.1) such that $\mathcal {N}(\bar {\boldsymbol q})=0$, and the perturbations $\boldsymbol q'(t)$ are small enough to discard ${O}^2(\boldsymbol q'(t))$ terms, we can treat $\boldsymbol f'(t)$ solely as an exogenous forcing such as a control input or free-stream disturbance (Jovanović & Bamieh Reference Jovanović and Bamieh2005). Such a scenario can arise for fluid flows at low Reynolds numbers. On the other hand, if the perturbations $\boldsymbol q'(t)$ around the base state are large enough (e.g. in turbulent flows), we retain the ${O}^2(\boldsymbol q'(t))$ terms and $\boldsymbol f'(t)$ has a nonlinear dependence on $\boldsymbol q'(t)$. Then $\boldsymbol f'(t)$ becomes a combination of endogenous forcing originating from the nonlinear interactions of the perturbations and any other presence of external influence in the system (McKeon & Sharma Reference McKeon and Sharma2010). Following previous studies (McKeon & Sharma Reference McKeon and Sharma2010; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018; Liu et al. Reference Liu, Sun, Yeh, Ukeiley, Cattafesta and Taira2021), we disregard the dependence of $\boldsymbol f'(t)$ on the state perturbation $\boldsymbol q'(t)$ and consider $\boldsymbol f'(t)$ as an unknown forcing term in the subsequent analysis. Substituting the Fourier transform of the perturbation $\boldsymbol q'(t)$ and the forcing $\boldsymbol f'(t)$,

(2.6a,b)\begin{equation} \boldsymbol q'(t) = \int_{-\infty}^{\infty} \hat{\boldsymbol q}'_{\omega}\, \mathrm{e}^{\mathrm{i} \omega t}\, \mathrm{d}\omega,\quad \boldsymbol f'(t) = \int_{-\infty}^{\infty} \hat{\boldsymbol{f}}'_{\omega}\, \mathrm{e}^{\mathrm{i} \omega t}\, \mathrm{d}\omega, \end{equation}

into (2.5) yields

(2.7)\begin{equation} \hat{\boldsymbol q}'_{\omega} = {\boldsymbol{\mathsf{R}}}_{\omega} \hat{\boldsymbol{f}}'_{\omega}, \end{equation}

where ${\boldsymbol{\mathsf{R}}}_{\omega } := (\mathrm {i} \omega {\boldsymbol{\mathsf{I}}} - {\boldsymbol{\mathsf{A}}})^{-1}\in \mathbb {C}^{N\times N}$ is the classical resolvent operator (Jovanović & Bamieh Reference Jovanović and Bamieh2005; McKeon & Sharma Reference McKeon and Sharma2010), where ${\boldsymbol{\mathsf{I}}}$ is a identity matrix. The operator ${\boldsymbol{\mathsf{R}}}_{\omega }$ acts as an open-loop transfer function from the input forcing $\hat {\boldsymbol {f}}'_{\omega }$ to the output response $\hat {\boldsymbol q}'_{\omega }$ at the frequency $\omega$.

2.3. Harmonic resolvent analysis framework

The harmonic resolvent analysis (Padovan et al. Reference Padovan, Otto and Rowley2020) extends the classical resolvent analysis for a time-varying system. In the harmonic resolvent analysis, we consider the base state $\bar {\boldsymbol q}(t)$ to be time periodic with a fundamental period $T$ and a fundamental frequency $\omega _p =2{\rm \pi} /T$. Since the Jacobian matrix ${\boldsymbol{\mathsf{A}}}(t)$ is evaluated about the base state $\bar {\boldsymbol q}(t)$, it inherits the time periodicity with the same period $T$ of the base state. In the present work we are interested in understanding the dynamics of a $T$-periodic perturbation $\boldsymbol q'(t)$ developing around the periodic base state. The perturbation dynamics need not be exactly $T$ periodic, and the analysis can be generalized for any $nT$-periodic perturbation, for an integer $n$ (Padovan & Rowley Reference Padovan and Rowley2022). We expand ${\boldsymbol{\mathsf{A}}}(t)$, $\boldsymbol q'(t)$ and $\boldsymbol f'(t)$ in terms of the Fourier series as

(2.8a–c)\begin{equation} {\boldsymbol{\mathsf{A}}}(t) = \sum_{k ={-}\infty}^{k= \infty} \hat{\!{\boldsymbol{\mathsf{A}}}}_k \, \mathrm{e}^{\mathrm{i} k\omega_p t},\quad \boldsymbol q'(t) = \sum_{k={-}\infty}^{k=\infty} \hat{\boldsymbol q}'_k \, \mathrm{e}^{\mathrm{i} k\omega_p t},\quad \boldsymbol f'(t) = \sum_{k={-}\infty}^{k=\infty} \hat{\boldsymbol{f}}'_k \, \mathrm{e}^{\mathrm{i} k\omega_p t}. \end{equation}

Substituting the Fourier series expansions into (2.3) yields

(2.9)\begin{equation} [\boldsymbol T \hat{\boldsymbol q}']_k := \mathrm{i} k\omega_p \hat{\boldsymbol q}'_k - \sum_{j={-}\infty}^{j=\infty} \hat{\!{\boldsymbol{\mathsf{A}}}}_{k-j} \hat{\boldsymbol q}'_j = \hat{\boldsymbol{f}}'_k,\quad \forall \, k,j\in \mathbb{Z}, \end{equation}

which represents a system of infinitely coupled equations, where perturbation $\hat {\boldsymbol q}'_k$ at the frequency $k\omega _p$ is coupled with the perturbation $\hat {\boldsymbol q}'_j$ at frequency $j\omega _p$ through the base state at frequency $(k-j)\omega _p$. In a matrix form, we can express the coupled system of equations as

(2.10)\begin{equation} {\boldsymbol{\mathsf{T}}} \hat{\mathcal{Q}} = \hat{\mathcal{F}},\end{equation}

where ${\boldsymbol{\mathsf{T}}}$ is an infinite-dimensional Toeplitz matrix of the form

(2.11)

with the infinite-dimensional state perturbation vector and the forcing vector of

(2.12a)$$\begin{gather} \hat{\mathcal{Q}} = [ \begin{array}{ccccccc} \dots & \hat{\boldsymbol q}'_{{-}2} & \hat{\boldsymbol q}'_{{-}1} & \hat{\boldsymbol q}'_{0} & \hat{\boldsymbol q}'_{1} & \hat{\boldsymbol q}'_{2} & \dots \end{array} ]^T, \end{gather}$$

(2.12b)$$\begin{gather}\hat{\mathcal{F}} = [ \begin{array}{ccccccc} \dots & \hat{\boldsymbol{f}}'_{{-}2} & \hat{\boldsymbol{f}}'_{{-}1} & \hat{\boldsymbol{f}}'_{0} & \hat{\boldsymbol{f}}'_{1} & \hat{\boldsymbol{f}}'_{2} & \dots \end{array} ]^T, \end{gather}$$

respectively. The diagonal of the matrix ${\boldsymbol{\mathsf{T}}}$ contains the block matrices of the form ${\boldsymbol{\mathsf{R}}}^{-1}_k := (\mathrm {i} k\omega _p {\boldsymbol{\mathsf{I}}}- \,\hat {\!{\boldsymbol{\mathsf{A}}}}_0)\in \mathbb {C}^{N\times N}$. The off-diagonal blocks of ${\boldsymbol{\mathsf{T}}}$ are the Fourier components of the Jacobian matrix $\,\hat {\!{\boldsymbol{\mathsf{A}}}}_j\in \mathbb {C}^{N\times N}$ at the frequency $j\omega _p$ with $j\in \mathbb {Z} \backslash \{0\}$. Since the matrix ${\boldsymbol{\mathsf{A}}}(t)$ is real valued, the Fourier component $\,\hat {\!{\boldsymbol{\mathsf{A}}}}_k$ is the complex conjugate of the component $\,\hat {\!{\boldsymbol{\mathsf{A}}}}_{-k}$ and vice versa. If the base state is steady then $\,\hat {\!{\boldsymbol{\mathsf{A}}}}_j=0,\ \forall \, j\in \mathbb {Z} \backslash \{0\}$, and the off-diagonal blocks of the matrix ${\boldsymbol{\mathsf{T}}}$ become zero. In the resulting system, the state perturbations at different frequencies are decoupled, and the non-zero diagonal elements of ${\boldsymbol{\mathsf{T}}}$ are the inverse of the classical resolvent operators at frequencies $k\omega _p$.

Next, we need to define an input–output relation between the forcing $\hat {\mathcal {F}}$ and state perturbation $\hat {\mathcal {Q}}$ in the frequency domain using the inverse of the operator ${\boldsymbol{\mathsf{T}}}$ in (2.10). However, if $\bar {\boldsymbol q}(t)$ satisfies (2.1) exactly, the operator ${\boldsymbol{\mathsf{T}}}$ is singular and contains a non-zero vector in the right null space (see Appendix A). If, however, the dynamics develop due to external forcing, then the null space becomes trivial. The existence of a singularity when the null space is non-trivial prevents an inversion of the operator ${\boldsymbol{\mathsf{T}}}$. Following the work by Padovan et al. (Reference Padovan, Otto and Rowley2020), we can restrict the domain and range of the operator ${\boldsymbol{\mathsf{T}}}$ to remove the singularity. By defining $\hat {\boldsymbol w}$ as a unit norm vector in the right null space, we can use the elementary orthogonal projector ${\boldsymbol{\mathsf{P}}}_{\mathcal {X}} ={\boldsymbol{\mathsf{I}}} - \hat {\boldsymbol w} \hat {\boldsymbol w}^*$ to project vectors on the subspace $\mathcal {X}$ that is an orthogonal complement to the right null space of ${\boldsymbol{\mathsf{T}}}$. Here, $({\cdot })^*$ denotes the complex conjugate transpose of a variable. Similarly, we can restrict the range of ${\boldsymbol{\mathsf{T}}}$ to a subspace $\mathcal {W}$, which is an orthogonal complement of the left null space of ${\boldsymbol{\mathsf{T}}}$ using the elementary projector ${\boldsymbol{\mathsf{P}}}_{\mathcal {W}} = {\boldsymbol{\mathsf{I}}} - \hat {\boldsymbol u} \hat {\boldsymbol u}^*$, where $\hat {\boldsymbol u}$ is a unit norm vector in the left null space of ${\boldsymbol{\mathsf{T}}}$. The computation of the unit norm vectors and the associated projection operators in practical applications are detailed in the study by Padovan et al. (Reference Padovan, Otto and Rowley2020). The restricted operator ${\boldsymbol{\mathsf{T}}}_w: \mathcal {X}\rightarrow \mathcal {W}$ is invertible and the input–output relation from (2.10) is obtained as

(2.13)\begin{equation} \hat{\mathcal{Q}} = {\boldsymbol{\mathsf{H}}} \hat{\mathcal{F}}, \end{equation}

where ${\boldsymbol{\mathsf{H}}} := {\boldsymbol{\mathsf{T}}}_w^{-1}$ is the harmonic resolvent operator. The operator ${\boldsymbol{\mathsf{H}}}$ maps the Fourier coefficients of the $T$-periodic forcing $\hat {\mathcal {F}}$ to the Fourier coefficients of the output state perturbation $\hat {\mathcal {Q}}$ of the same period $T$.

2.4. Harmonic resolvent formulation for compressible flow

2.4.1. Navier–Stokes equation

In this section we derive the harmonic resolvent formulation for the fluid flow governing equations. We consider the compressible NSE in the conservative formulation. The Cartesian coordinate system $x_i (i=1, 2, 3)$, time $t$, density $\rho$, three components of velocity $u_i$, pressure $p$, temperature $T$ and total energy $E$ are non-dimensionalized as

(2.14) \begin{align} x_i = \frac{\tilde x_i}{L}, \qquad t = \frac{\tilde t}{L/ \tilde u_{\infty}}, \qquad \rho = \frac{\tilde\rho}{\tilde \rho_{\infty}}, \qquad u_i = \frac{\tilde u_i}{\tilde u_{\infty}}, \qquad p = \frac{\tilde p}{\tilde \rho_{\infty} \tilde u_{\infty}^2}, \qquad T =\frac{\tilde T}{\tilde T_{\infty}}, \qquad E = \frac{\tilde E}{\tilde \rho_{\infty} \tilde u_{\infty}^2},\end{align}

where the variables denoted with the symbol $\widetilde {({\cdot })}$ are the dimensional quantities, the variables with the subscript $\infty$ denote free-stream values and $L$ is a dimensional reference length. We introduce three dimensionless numbers, namely the Reynolds number $Re$, Prandtl number $Pr$ and Mach number $Ma$,

(2.15a–c)\begin{equation} Re = \frac{\tilde \rho_{\infty} \tilde u_{\infty} L}{\tilde \mu_{\infty}},\quad Pr = \frac{\tilde \mu_{\infty} \tilde c_p}{\tilde \kappa},\quad Ma = \frac{\tilde u_{\infty}}{\tilde a_{\infty}}, \end{equation}

where $\tilde \mu _{\infty }$ is the free-stream dynamic viscosity, $\tilde a_{\infty }$ is the speed of sound in the free stream, $\tilde c_p$ is the specific heat at constant pressure and $\tilde \kappa$ is the thermal conductivity of the fluid. Then we can compactly write the NSE in a non-dimensional form as

(2.16)\begin{equation} \frac{\partial \boldsymbol q}{\partial t} + \frac{\partial \boldsymbol F_j^e}{\partial x_j} + \frac{\partial \boldsymbol F_j^v}{\partial x_j} =0 , \end{equation}

where $\boldsymbol q = [\rho, m_i, \rho E]^T\in \mathbb {R}^5$ is the vector of the conservative state variables with $m_i:= \rho u_i$ being the three components of the momentum, $\boldsymbol F_j^e$ and $\boldsymbol F_j^v$ represents the Euler flux and viscous flux vector, respectively. For a thermally and calorically perfect gas, total energy $E$ is given by

(2.17)\begin{equation} E = \frac{p}{\rho(\gamma_s -1)} + \frac{1}{2} u_k u_k, \end{equation}

where $\gamma_s$ is the ratio of specific heat. The Euler flux $\boldsymbol F_j^e$ and the viscous flux $\boldsymbol F_j^v$ (Bugeat et al. Reference Bugeat, Chassaing, Robinet and Sagaut2019) are given by

(2.18a,b) \begin{equation} \boldsymbol F_j^e = \left[ \begin{array}{@{}c@{}} m_j\\ m_i u_j + p \delta_{ij}\\ (\rho E+p)u_j \end{array} \right] ,\qquad \boldsymbol F_j^v = \left[ \begin{array}{@{}c@{}} 0\\ \displaystyle -\dfrac{1}{Re} \tau_{ij}\\ \displaystyle -\dfrac{1}{Re} u_i \tau_{ij} - \dfrac{\kappa}{(\gamma_s-1)Re\, Pr\, Ma^2} \dfrac{\partial T}{\partial x_j} \end{array} \right] ,\end{equation}

where

(2.19)\begin{equation} \tau_{ij} = \mu \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} - \frac{2}{3} \frac{\partial u_k}{\partial x_k} \delta_{ij}\right) \end{equation}

is the viscous stress tensor for a Newtonian fluid. Note that the coefficients of the viscous stress term and the temperature gradient term in (2.18b) depend on the reference variables used to non-dimensionalize the governing equations. The dynamic viscosity $\mu$ is a function of the temperature, which is calculated using a power law as $\mu (T) = (T/T_{\infty })^{0.76}$ (Garnier, Adams & Sagaut Reference Garnier, Adams and Sagaut2009). We use the equation of state to relate the pressure, density and temperature as

(2.20)\begin{equation} p = \frac{1}{\gamma_s\, Ma^2} \rho T. \end{equation}

We take the value of the Prandtl number $Pr=0.7$ and the specific heat ratio $\gamma _s=1.4$, which are the standard values for air.

2.4.2. Linearized NSE

We decompose the state variables $\boldsymbol q(t)$ into a periodic base state $\bar {\boldsymbol q}(t)$ and an unsteady perturbation $\boldsymbol q'(t)$ as $\boldsymbol q(t) = \bar {\boldsymbol q}(t) + \boldsymbol q'(t)$. Substituting the decomposition into (2.16) and linearizing around the base state, we obtain

(2.21)\begin{equation} \frac{\partial \boldsymbol q'(t)}{{\rm d} t} + \frac{\partial}{\partial x_j} \mathcal{F}_j^e(\bar{\boldsymbol q}(t),\boldsymbol q'(t)) + \frac{\partial}{\partial x_j} \mathcal{F}_j^v(\bar{\boldsymbol q}(t),\boldsymbol q'(t)) = \boldsymbol f'(t), \end{equation}

where $\boldsymbol q'(t) =[\rho ',m'_i,\rho E']$ is the vector of conservative variable perturbations, and $\boldsymbol f'(t)$ contains the terms that are nonlinear in $\boldsymbol q'(t)$ and the residual terms if $\bar {\boldsymbol q}(t)$ is not an exact solution of (2.16) as

(2.22)\begin{equation} \boldsymbol f'(t) ={-}\frac{\partial \bar{\boldsymbol q}(t)}{{\rm d} t} - \frac{\partial}{\partial x_j} \mathcal{F}_j^e(\bar{\boldsymbol q}(t)) - \frac{\partial}{\partial x_j} \mathcal{F}_j^v(\bar{\boldsymbol q}(t)) -{O}^2(\boldsymbol q'(t)). \end{equation}

The linearized Euler flux $\mathcal {F}_j^e(\bar {\boldsymbol q}(t),\boldsymbol q'(t))$ reads

(2.23) \begin{equation} \mathcal{F}_j^e(\bar{\boldsymbol q}(t),\boldsymbol q'(t)) = \left[ \begin{array}{@{}c@{}} m'_j\\ \displaystyle \dfrac{\bar m_i}{\bar \rho} m'_j + \dfrac{\bar m_j}{\bar \rho} m'_i - \dfrac{\bar m_i \bar m_j}{\bar{\rho}^2} \rho' + p'\delta_{ij}\\ \displaystyle \left(\gamma_s \bar{\rho E} -\dfrac{\gamma_s-1}{2} \dfrac{\bar m_k \bar m_k}{\bar \rho} \right) u'_j + \dfrac{\bar m_j}{\bar \rho} (\rho E)' + \dfrac{\bar m_j}{\bar \rho} p' \end{array} \right] , \end{equation}

and $\mathcal {F}_j^v(\bar {\boldsymbol q}(t),\boldsymbol q'(t))$ is the linearized viscous flux vector with

(2.24) \begin{equation} \mathcal{F}_j^v(\bar{\boldsymbol q}(t),\boldsymbol q'(t)) = \left[ \begin{array}{@{}c@{}} 0\\ \displaystyle - \dfrac{1}{Re} \tau'_{ij}\\ \displaystyle -\dfrac{1}{Re}\, u'_i \bar{\tau}_{ij} -\dfrac{1}{Re}\, \bar u_i \tau'_{ij} - \dfrac{\bar \mu}{(\gamma_s-1)Re\, Pr\, Ma^2}\, \dfrac{\partial T'}{\partial x_j} \end{array} \right] , \end{equation}

where

(2.25a,b)\begin{equation} \bar \tau_{ij} = \bar \mu \left(\frac{\partial \bar u_i}{\partial x_j} + \frac{\partial \bar u_j}{\partial x_i} - \frac{2}{3} \frac{\partial \bar u_k}{\partial x_k} \delta_{ij}\right),\quad \tau'_{ij} = \bar \mu \left(\frac{\partial u'_i}{\partial x_j} + \frac{\partial u'_j}{\partial x_i} - \frac{2}{3} \frac{\partial u'_k}{\partial x_k} \delta_{ij}\right) \end{equation}

and $\bar u_i := \bar m_i/\bar \rho$ denotes the velocity components of the base state. We neglect the terms with viscosity perturbation $\mu '$ by assuming its negligible variation with temperature. The perturbations of primitive variable velocity $u_i'$, pressure $p'$ and temperature $T'$ is calculated respectively as

(2.26a)$$\begin{gather} u_i' = \frac{1}{\bar \rho} m_i' - \frac{\bar m_i}{\bar{\rho}^2} \rho', \end{gather}$$

(2.26b)$$\begin{gather}p' = (\gamma_s-1)\left[(\rho E)' - \frac{\bar m_k}{\bar \rho} m_k' +\frac{1}{2} \frac{\bar m_k \bar m_k}{ \bar{\rho}^2} \rho'\right], \end{gather}$$

(2.26c)$$\begin{gather}T' = \gamma_s (\gamma_s -1) Ma^2 \left[\frac{(\rho E)'}{\bar \rho} - \frac{\bar{\rho E}}{\bar{\rho}^2}\rho' - \frac{\bar m_k}{\bar{\rho}^2} m'_k + \frac{\bar m_k \bar m_k}{\bar{\rho}^3} \rho'\right]. \end{gather}$$

After substituting all the expressions into (2.21), we obtain the governing equation of unsteady perturbations developing over a time-varying compressible fluid flow in the time domain.

2.4.3. Construction of operator ${\boldsymbol{\mathsf{T}}}$

To facilitate the conversion of the LNSE from the time domain to the frequency domain, we rewrite (2.21) as

(2.27)\begin{equation} \frac{{\rm d} \boldsymbol q'(t)}{{\rm d} t} + {\boldsymbol{\mathsf{L}}} \boldsymbol q'(t) + \frac{\partial}{\partial x_j} [\mathcal{G}_j^e(\bar{\boldsymbol q}(t),\boldsymbol q'(t)) + \mathcal{F}_j^v(\bar{\boldsymbol q}(t),\boldsymbol q'(t))] = \boldsymbol f'(t), \end{equation}

where we have grouped the terms of the linearized Euler flux that contain the product between the base state and perturbation state variables in $\mathcal {G}_j^e(\bar {\boldsymbol q}(t),\boldsymbol q'(t))$ and the rest of the terms containing only the perturbation state variables in ${\boldsymbol{\mathsf{L}}}\boldsymbol {q}'(t)$ (expressions are given in Appendix B). Then, we expand the periodic base state and perturbation using the Fourier series as

(2.28a,b)\begin{equation} \bar{\boldsymbol q}(t) = \sum_{k \in \varOmega} \hat{\bar{\boldsymbol q}}_k \, \mathrm{e}^{\mathrm{i} k\omega_p t},\quad \boldsymbol q'(t) = \sum_{k\in \tilde{\varOmega}} \hat{\boldsymbol q}'_k \, \mathrm{e}^{\mathrm{i} k\omega_p t},\end{equation}

where both $\varOmega,\tilde {\varOmega } \subseteq \{k\omega _p\}\ \forall \, k\in \mathbb {Z}$, are sets of integer multiples of the fundamental frequency $\omega _p = 2{\rm \pi} /T$ with $T$ being the fundamental period of the base flow. While one can consider an infinite number of frequencies for the base state and the perturbation, in practical computations, we truncate the number of frequencies in the sets $\varOmega$ and $\tilde {\varOmega }$ to a finite extent. Usually, $\varOmega =\{-m,\ldots,-1,0,1,\ldots,m\}\omega _p$ contains a small number of frequencies associated with the dominant frequency $\omega _p$ and its harmonics present in the base flow. These frequencies approximate the dominant dynamics of the large-scale coherent structures in the fluid flows. Then we seek to study the dynamics of perturbations with frequencies in the set of $\tilde {\varOmega }=\{-n,\ldots,-1,0,1,\ldots,n\}\omega _p$, where $n$ is chosen by the maximum frequency of perturbation that one wishes to resolve and $n\geq m$. Substituting the Fourier expansions into (2.27), we obtain the following system of a finite number of coupled equations:

(2.29)\begin{equation} [{\boldsymbol{\mathsf{T}}}\hat{\boldsymbol q}']_k = \hat{\boldsymbol{f}}'_k \end{equation}

with

(2.30)\begin{equation} [{\boldsymbol{\mathsf{T}}}\hat{\boldsymbol q}']_k = \mathrm{i}k\omega_p \hat{\boldsymbol q}'_k + {\boldsymbol{\mathsf{L}}} \hat{\boldsymbol q}'_k + \frac{\partial}{\partial x_j} \sum_{\substack{l\in \tilde{\varOmega}\\(k-l)\in\varOmega}} [\hat{\mathcal{G}}_j^e(\hat{\bar{\boldsymbol q}}_{k-l},\widehat{\boldsymbol q'}_l) + \hat{\mathcal{F}}_j^v(\hat{\bar{\boldsymbol q}}_{k-l},\widehat{\boldsymbol q'}_l)].\end{equation}

Here the number of equations is linked to the number of perturbation frequencies in the set $\tilde {\varOmega }$. To assemble the matrix ${\boldsymbol{\mathsf{T}}}$, we need to find the expressions for the equations corresponding to the frequencies in $\tilde \varOmega$ one at a time. For simplicity, we consider a set of base flow frequencies $\varOmega = \{-1,0,1\}\omega _p$ along with a set of perturbation frequencies $\tilde {\varOmega } = \{-2,-1,0,1,2\}\omega _p$ to demonstrate the construction of the operator ${\boldsymbol{\mathsf{T}}}$. The equation corresponding to the perturbation frequency $-2\omega _p$ can be obtained as

(2.31)\begin{align} [{\boldsymbol{\mathsf{T}}}\hat{\boldsymbol q}']_{{-}2} &={-}\mathrm{i}2\omega_p \hat{\boldsymbol q}'_{{-}2} + {\boldsymbol{\mathsf{L}}} \hat{\boldsymbol q}'_{{-}2} + \frac{\partial}{\partial x_j} [\hat{\mathcal{G}}_j^e(\hat{\bar{\boldsymbol q}}_{0},\hat{\boldsymbol q'}_{{-}2}) + \hat{\mathcal{F}}_j^v(\hat{\bar{\boldsymbol q}}_{0},\hat{\boldsymbol q'}_{{-}2})]\nonumber\\ &\quad + \frac{\partial}{\partial x_j} [\hat{\mathcal{G}}_j^e(\hat{\bar{\boldsymbol q}}_{{-}1},\widehat{\boldsymbol q'}_{{-}1}) + \hat{\mathcal{F}}_j^v(\hat{\bar{\boldsymbol q}}_{{-}1},\widehat{\boldsymbol q'}_{{-}1})], \end{align}

where we have neglected the terms containing $\hat {\bar {\boldsymbol q}}_{k-l}$ with $(k-l)\not \in \varOmega$ in the expansion of the sum. For brevity, we show how to perform the Fourier expansion of the terms in $\hat {\mathcal {G}}_j^e(\hat {\bar {\boldsymbol q}}_{k-l},\widehat {\boldsymbol q'}_{l})$ for the linearized compressible NSE in Appendix B. Similarly, for other frequencies $k\omega _p$ in the set $\tilde \varOmega$, we can obtain the expression for $[ {\boldsymbol{\mathsf{T}}}\hat {\boldsymbol q}']_k$ using (2.30). Then the system of equations in a matrix form is

(2.32) \begin{equation} {\boldsymbol{\mathsf{T}}} = \left[ \begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{}} {\boldsymbol{\mathsf{R}}}_{{-}2}^{{-}1} & \hat{{\boldsymbol{\mathsf{G}}}}_{{-}1} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} \\ \hat{{\boldsymbol{\mathsf{G}}}}_1 & {\boldsymbol{\mathsf{R}}}_{{-}1}^{{-}1} & \hat{{\boldsymbol{\mathsf{G}}}}_{{-}1} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} \\ {\boldsymbol{\mathsf{0}}} & \hat{{\boldsymbol{\mathsf{G}}}}_1 & {\boldsymbol{\mathsf{R}}}_{0}^{{-}1} & \hat{{\boldsymbol{\mathsf{G}}}}_{{-}1} & {\boldsymbol{\mathsf{0}}}\\ {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & \hat{{\boldsymbol{\mathsf{G}}}}_1 & {\boldsymbol{\mathsf{R}}}_{1}^{{-}1} & \hat{{\boldsymbol{\mathsf{G}}}}_{{-}1}\\ {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & \hat{{\boldsymbol{\mathsf{G}}}}_1 & {\boldsymbol{\mathsf{R}}}_{2}^{{-}1} \end{array} \right] , \end{equation}

where

(2.33a)$$\begin{gather} \hat{{\boldsymbol{\mathsf{G}}}}_k \hat{\boldsymbol q}' := \frac{\partial}{\partial x_j} [\hat{\mathcal{G}}_j^e(\hat{\bar{\boldsymbol q}}_{k},\widehat{\boldsymbol q'}) + \hat{\mathcal{F}}_j^v(\hat{\bar{\boldsymbol q}}_{k},\widehat{\boldsymbol q'})], \end{gather}$$

(2.33b)$$\begin{gather}{\boldsymbol{\mathsf{R}}}^{{-}1}_k := (\mathrm{i}k\omega_p{\boldsymbol{\mathsf{I}}} + {\boldsymbol{\mathsf{L}}} + \hat{{\boldsymbol{\mathsf{G}}}}_0). \end{gather}$$

The operator ${\boldsymbol{\mathsf{T}}}$ has dimension $\mathbb {C}^{5N_f\times 5N_f}$, where $N_f$ is the number of frequencies in the set $\tilde \varOmega$. The number of non-zero blocks in each row of the operator ${\boldsymbol{\mathsf{T}}}$ depends on the number of base flow frequencies in the set $\varOmega$. To numerically solve the system of equations, we discretize (2.30) using a finite volume scheme in the present work and obtain the discrete matrices ${\boldsymbol{\mathsf{L}}}\in \mathbb {R}^{5N_g\times 5N_g}$ and $\hat {{\boldsymbol{\mathsf{G}}}}_k \in \mathbb {C}^{5N_g\times 5N_g}$ that constitutes the blocks of the operator ${\boldsymbol{\mathsf{T}}}$, with $N_g$ being the number of discrete grid points. We note that only the matrices $\hat {{\boldsymbol{\mathsf{G}}}}_k$ need to be constructed, and $\hat {{\boldsymbol{\mathsf{G}}}}_{-k}$ can be obtained by taking the complex conjugate of $\hat {{\boldsymbol{\mathsf{G}}}}_{k}$. After assembling the discrete operator ${\boldsymbol{\mathsf{T}}}$, we can remove its singularity, if necessary, using the method outlined in § 2.3. Then the input–output perturbation dynamics of a time-periodic fluid flow in the frequency space are represented by

(2.34)\begin{equation} \hat{\mathcal{Q}} = {\boldsymbol{\mathsf{H}}} \hat{\mathcal{F}},\end{equation}

where, ${\boldsymbol{\mathsf{H}}} \in \mathbb {C}^{{(5N_g)N_f\times (5N_g)N_f}}$ is the discrete harmonic resolvent operator, $\hat {\mathcal {Q}}\in \mathbb {C}^{(5N_g)N_f\times 1}$ contains the collection of Fourier coefficients of the discrete state variable perturbations and $\hat {\mathcal {F}}\in \mathbb {C}^{(5N_g)N_f\times 1}$ represents the Fourier coefficients of the discrete forcing variables.

2.5. Modal decomposition of the harmonic resolvent operator

We seek to obtain a reduced-order representation of the input–output dynamics using the modal decomposition of the harmonic resolvent operator. In particular, we need to identify the most amplified output perturbation and the corresponding input perturbation characterized by a gain describing the amplification level. We can define the gain as a ratio of the output to input perturbation energy. To measure the perturbation energy, we introduce the discrete inner products $\langle \hat {\boldsymbol q}',\hat {\boldsymbol q}' \rangle _{q} = \hat {\boldsymbol q}^{'*} {\boldsymbol{\mathsf{W}}}_{q}\, \hat {\boldsymbol q}'$ and $\langle\, \boldsymbol {\hat f}',\boldsymbol {\hat f}'\rangle _{f} = \boldsymbol {\hat f}^{'*} {\boldsymbol{\mathsf{W}}}_{f} \boldsymbol {\hat f}'$ in the output and input space, respectively. Then using the inner product, we can define a norm to measure the perturbation energy in both spaces. The positive definite weight matrices ${\boldsymbol{\mathsf{W}}}_{q}$ and ${\boldsymbol{\mathsf{W}}}_{f}$ depend on the choice of energy norm one wishes to optimize. For compressible flows, a widely used measure of the perturbation energy is given by Chu's norm (Chu Reference Chu1965; Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996), which in the non-dimensional form is

(2.35) \begin{align} \|\hat{\boldsymbol q}'_p\|_E^2 = \hat{\boldsymbol q}_p^{' *} {\boldsymbol{\mathsf{W}}}_{C}\hat{\boldsymbol q}'_p = E_{Chu} = \frac{1}{2}\int_{\mathcal{V}} \hat{\boldsymbol q}_p^{' *} \text{diag}\left(\frac{\bar T_0}{\gamma_s Ma^2 \bar \rho_0},\bar \rho_0,\frac{\bar \rho_0}{\gamma_s(\gamma_s-1)Ma^2 \bar T_0}\right) \hat{\boldsymbol q}'_p \, \mathrm{d}\mathcal{V}, \end{align}

where $\boldsymbol q_p = [\rho ',u'_i,T']^T$ is the vector of primitive variable perturbation and $\hat {\boldsymbol q}'_p$ are the Fourier coefficients. The variables denoted with $({\cdot })_0$ represent the time-averaged quantity. Since the harmonic resolvent formulation is derived using the conservative state variable, we need to modify (2.34) to accommodate the primitive variable perturbations before applying Chu's norm. The transformed equation in the primitive state variable reads

(2.36)\begin{equation} \hat{\mathcal{Q}}_p = \underbrace{{\boldsymbol{\mathsf{M}}}^{{-}1} {\boldsymbol{\mathsf{H}}} {\boldsymbol{\mathsf{M}}}}_{{\boldsymbol{\mathsf{H}}}_p} \hat{\mathcal{F}}_p, \end{equation}

where $\hat {\mathcal {Q}}_p$ is the vector of Fourier coefficients of the output perturbation $\boldsymbol q'_p(t)$, $\hat {\mathcal {F}}_p$ contains the Fourier coefficients of $\boldsymbol f'_p(t)$. The operator ${\boldsymbol{\mathsf{H}}}_p$ governs the input–output dynamics of the primitive state variable perturbations and the details of the matrix ${\boldsymbol{\mathsf{M}}}$ are given in Appendix C. We seek to maximize the gain

(2.37)\begin{equation} \varGamma^2 = \max_{\hat{\mathcal{F}}_p} \frac{\|\hat{\mathcal{Q}}_p\|_E^2}{\|\hat{\mathcal{F}}_p\|_E^2} = \max_{\hat{\mathcal{F}}_p} \frac{\langle \hat{\mathcal{Q}}_p,\hat{\mathcal{Q}}_p \rangle_{q}}{\langle \hat{\mathcal{F}}_p,\hat{\mathcal{F}}_p \rangle_{f}} = \max_{\hat{\mathcal{F}}_p} \frac{\hat{\mathcal{Q}}^*_p {\boldsymbol{\mathsf{W}}}_{C}\, \hat{\mathcal{Q}}_p}{\hat{\mathcal{F}}_p^* {\boldsymbol{\mathsf{W}}}_{C}\, \hat{\mathcal{F}}_p}, \end{equation}

where we use the same measure of perturbation energy using Chu's norm (i.e. ${\boldsymbol{\mathsf{W}}}_c$ contain the Chu's norm weights along with the discrete integration weights) in both input and output spaces leading to ${\boldsymbol{\mathsf{W}}}_q = {\boldsymbol{\mathsf{W}}}_f = {\boldsymbol{\mathsf{W}}}_c$. It is not necessary to use the same energy norm in both spaces. Since ${\boldsymbol{\mathsf{W}}}_c$ is a symmetric positive definite matrix, we can perform a Cholesky factorization as ${\boldsymbol{\mathsf{W}}}_c = {\boldsymbol{\mathsf{N}}}^* {\boldsymbol{\mathsf{N}}}$. Using the factorization in (2.37), we can transform Chu's energy norm to a discrete $L_2$ norm as

(2.38) \begin{align} \varGamma^2 = \max_{\hat{\mathcal{F}}_p} \frac{\hat{\mathcal{F}}_p^* {\boldsymbol{\mathsf{H}}}_p^* {\boldsymbol{\mathsf{N}}}^* {\boldsymbol{\mathsf{N}}} {\boldsymbol{\mathsf{H}}}_p \hat{\mathcal{F}}_p}{\hat{\mathcal{F}}_p^* {{\boldsymbol{\mathsf{N}}}}^* {\boldsymbol{\mathsf{N}}} \hat{\mathcal{F}}_p} &= \max_{\hat{\mathcal{U}}_p} \frac{\hat{\mathcal{U}}_p^* {\boldsymbol{\mathsf{N}}}^{{-}1,*} {\boldsymbol{\mathsf{H}}}_p^* {{\boldsymbol{\mathsf{N}}}}^* {\boldsymbol{\mathsf{N}}} {\boldsymbol{\mathsf{H}}}_p {\boldsymbol{\mathsf{N}}}^{{-}1} \hat{\mathcal{U}}_p}{\hat{\mathcal{U}}_p^* \hat{\mathcal{U}}_p }\quad \text{(with}\qquad \hat{\mathcal{U}}_p = {\boldsymbol{\mathsf{N}}} \hat{\mathcal{F}}_p) \nonumber\\ &= \max_{\hat{\mathcal{U}}_p} \frac{\|{\boldsymbol{\mathsf{N}}} {\boldsymbol{\mathsf{H}}}_p {\boldsymbol{\mathsf{N}}}^{{-}1} \hat{\mathcal{U}}_p\|_2^2}{\|\hat{\mathcal{U}}_p\|_2^2} \nonumber\\ &= \|{\boldsymbol{\mathsf{N}}} {\boldsymbol{\mathsf{H}}}_p {\boldsymbol{\mathsf{N}}}^{{-}1}\|_2^2 = \sigma_1^2. \end{align}

So the solution to the optimization problem can be obtained by the SVD of the weighted harmonic resolvent operator

(2.39)\begin{equation} {\boldsymbol{\mathsf{N}}} {\boldsymbol{\mathsf{H}}}_p {\boldsymbol{\mathsf{N}}}^{{-}1} = \tilde{{\boldsymbol{\mathsf{U}}}} \boldsymbol \varSigma \tilde{{\boldsymbol{\mathsf{V}}}}^*, \end{equation}

where $\boldsymbol \varSigma = \text {diag}(\sigma _1,\sigma _2,\ldots )$ contains the ranked singular values of the operator ${\boldsymbol{\mathsf{N}}} {\boldsymbol{\mathsf{H}}}_p {\boldsymbol{\mathsf{N}}}^{-1}$ in descending order and the maximum energy gain $\varGamma ^2$ is given by the leading singular value $\sigma _1^2$. The response modes are the columns of the matrix ${\boldsymbol{\mathsf{U}}} = {\boldsymbol{\mathsf{N}}}^{-1} \tilde {{\boldsymbol{\mathsf{U}}}}$ and the forcing modes are the columns of the matrix ${\boldsymbol{\mathsf{V}}} = {\boldsymbol{\mathsf{N}}}^{-1} \tilde {{\boldsymbol{\mathsf{V}}}}$. The optimal forcing and response modes are given by the first column of the matrix ${\boldsymbol{\mathsf{V}}}$ and ${\boldsymbol{\mathsf{U}}}$, respectively. The forcing and response modes are orthonormal in their respective inner products, that is, ${\boldsymbol{\mathsf{V}}}^* {\boldsymbol{\mathsf{W}}}_c {\boldsymbol{\mathsf{V}}} = {\boldsymbol{\mathsf{I}}}$ and ${\boldsymbol{\mathsf{U}}}^* {\boldsymbol{\mathsf{W}}}_c {\boldsymbol{\mathsf{U}}} = {\boldsymbol{\mathsf{I}}}$. Next, we recover the complete decomposition of the operator as

(2.40)\begin{equation} {\boldsymbol{\mathsf{H}}}_p = {\boldsymbol{\mathsf{U}}} \boldsymbol{\varSigma} {\boldsymbol{\mathsf{V}}}^* {\boldsymbol{\mathsf{W}}}_c, \end{equation}

allowing the output response to be expanded as

(2.41)\begin{equation} \hat{\mathcal{Q}}_p = \sum_{k} \boldsymbol U_k \sigma_k \lambda_k \quad \text{with } \lambda_k = \boldsymbol V_k^* {\boldsymbol{\mathsf{W}}}_c \hat{\mathcal{F}}_p. \end{equation}

If the operator ${\boldsymbol{\mathsf{H}}}_p$ is low rank and $\sigma _1 \gg \sigma _2$, then we can use a rank-1 approximation to get a reduced-order representation of the dynamics

(2.42)\begin{equation} \hat{\mathcal{Q}}_p = \boldsymbol U_1 \sigma_1 \lambda_1 \quad \text{with } \lambda_1 = \boldsymbol V_1^* {\boldsymbol{\mathsf{W}}}_c \hat{\mathcal{F}}_p. \end{equation}

If the projection $\lambda _1$ of the input on the forcing mode $\boldsymbol V_1$ is maximal, the output response will have structures similar to the response mode $\boldsymbol U_1$ scaled by the singular value $\sigma _1$. In other words, if we want to excite the optimal response, our input to the system needs to be aligned as closely as possible to the optimal forcing mode. The reduced-order representation has profound significance in understanding the physics of the time-periodic fluid flows and developing inputs for flow control.

We have seen till now that the SVD of the weighted harmonic resolvent operator sheds light on the global energy amplification mechanism in the time-periodic flow. However, since the time-periodic base flow admits cross-frequency interaction between perturbations, it is possible to study the energy amplification between a pair of input and output perturbations at different frequencies (Padovan et al. Reference Padovan, Otto and Rowley2020). In particular, we want to maximize the gain between the input energy perturbation $\hat {\boldsymbol f}'_j$ at the frequency $j\omega _p\in \tilde \varOmega$ to the output response $\hat {\boldsymbol q}'_{k}$ at the frequency $k\omega _p\in \tilde \varOmega$. As shown before in (2.38), the solution to the optimization leads to the SVD of the weighted operator ${\boldsymbol{\mathsf{H}}}^w_{j,k}\in \mathbb {C}^{5N_g\times 5N_g}$, which is the corresponding block of the operator ${\boldsymbol{\mathsf{N}}} {\boldsymbol{\mathsf{H}}}_p {\boldsymbol{\mathsf{N}}}^{-1}$ that couples the input $\hat {\boldsymbol f}'_j$ to the output $\hat {\boldsymbol q}'_k$ as

(2.43)\begin{equation} \widehat{\boldsymbol q'}_{k} = {\boldsymbol{\mathsf{H}}}^w_{j,k}\, \hat{\boldsymbol f}'_j. \end{equation}

The optimal singular value $\sigma _{1,(j,k)}$ provides a measure of how effectively the forcing $\hat {\boldsymbol f}'_j$ can excite the response at $\hat {\boldsymbol q}'_k$.

2.6. Harmonic resolvent analysis of $nT$-periodic perturbations

The harmonic resolvent analysis framework we have discussed till now considers the set perturbation frequencies with the same period as the fundamental base flow frequency. However, Padovan & Rowley (Reference Padovan and Rowley2022) showed that it is possible to analyse the dynamics of a $nT$-periodic perturbation developing over a $T$-periodic base flow using a modified harmonic resolvent operator ${\boldsymbol{\mathsf{H}}}(\mathrm {i}\gamma ) = (\mathrm {i}\gamma {\boldsymbol{\mathsf{I}}}+{\boldsymbol{\mathsf{T}}})^{-1}$ parameterized by $\gamma \in [-\omega _p/2,\omega _p/2]$, with $\omega _p = 2{\rm \pi} /T$. The operator ${\boldsymbol{\mathsf{H}}}(\mathrm {i}\gamma )$ describes the input–output perturbation dynamics of the set of perturbation frequencies $\tilde {\varOmega }_{\gamma } = \gamma + \{-n,\ldots,-1,0,1,\ldots,n\}\omega _p$. For $\gamma =0$, we obtain the harmonic resolvent operator defined in (2.34). Note that the operator ${\boldsymbol{\mathsf{T}}}$ only depends on the set of base flow frequencies and remains the same for any set of perturbation frequencies defined by $\gamma$. In defining ${\boldsymbol{\mathsf{H}}}(\mathrm {i}\gamma )$, Padovan & Rowley (Reference Padovan and Rowley2022) has used a slightly different projection operator than the one described in § 2.3. We use the projection operator defined at the end of § 2.3 to reduce the computational cost in this work. The modal decomposition of the modified harmonic resolvent operator ${\boldsymbol{\mathsf{H}}}(\mathrm {i}\gamma )$ is performed in the same way as described in § 2.5 to analyse the dynamics of the $nT$-periodic perturbations.

4.1. Cavity flow at $Ma_\infty =0.6$

We discuss the classical and harmonic resolvent analyses of the cavity flow at Mach 0.6 due to its simplicity of containing only one fundamental frequency. In the classical resolvent analysis we use the time-averaged mean flow state as the base flow to linearize the governing equations, equivalent to linearizing the equations about the base flow frequency set $\varOmega =\{0\}$. In the harmonic resolvent analysis we consider a set of truncated base flow frequencies of $\varOmega =\{-1,0,1\}\omega _p$ leading to the number of base flow frequencies $N_b=3$. We vary the set of perturbation frequencies $\tilde \varOmega = \{-m,\ldots,-1,0,1\ldots,m\}\omega _p$ with $m= 1,3,5$, resulting in the number of perturbation frequencies $N_f= 3,7$ and $11$, respectively. We build the linear operators for both analyses using a smaller domain with extent $x/D\in [-3,8]$ and $y/D\in [-1,7]$. Approximately $48\ 000$ grid points are used to discretize the domain. At the cavity surface and the upstream and downstream wall, the velocity perturbations and the wall-normal gradient of the pressure perturbation are specified to be zero. The velocity and density perturbations and gradient of pressure perturbations are specified as zero at the inlet. Sponge layers with an extent of $1D$ are applied near the top and outflow boundaries to dampen the perturbations and prevent any reflection back into the domain. We use $10$ random test vectors to compute the SVD of the harmonic resolvent operator using the randomized algorithm (Ribeiro et al. Reference Ribeiro, Yeh and Taira2020) along with three power iterations. For the classical resolvent operator, we compute the SVD using the Krylov-based Arnoldi-iteration method with the number of singular values of $4$, Krylov space of $128$ vectors and a residual tolerance of $10^{-10}$.

Variation of the first two singular values of the classical resolvent operator as a function of frequency is shown in figure 5(a). The flow is most responsive to perturbation at the frequency $\omega _p$ where the large separation between the optimal ($\sigma _1$) and sub-optimal ($\sigma _2$) singular values indicates a rank-1 behaviour. Since the leading two singular values overlap at higher frequencies ($\omega \geq \omega _p$), the rank-1 feature is no longer valid. In figure 5(b) we plot the first 10 singular values of the harmonic resolvent operator constructed by considering three different numbers of perturbation frequencies ($N_f$) with the same number of base frequencies $N_b=3$. Solving for an increasing number of perturbation frequencies stacked in one singular vector leads to a reduction in the optimal singular value ($\sigma _1$) when $N_f$ increases from 3 to 7, but a further increase of $N_f$ from 7 to 11 implies a minimal effect on capturing the dominant dynamics of the flow as shown in figure 5(b). The effect of increasing the perturbation frequencies in the set $\tilde {\varOmega }$ on the remaining sub-optimal singular values is trivial. Compared with the singular values of the classical resolvent operator at the dominant frequency ($\omega _p$), the separation between the optimal ($k=1$) and sub-optimal ($k=2$) singular values of the harmonic resolvent operator is smaller. Unlike the classical resolvent operator, the singular values of the harmonic resolvent operator do not reveal the energy amplification of an individual frequency but rather a combined effect of all the coupled frequencies in a set. However, since we solve for cross-frequency modes in one stacked vector, the relative amplitude of each perturbation in the harmonic resolvent mode will reveal their relative dominance, which we will see shortly.

Figure 5. (a) The first and second singular values of the classical resolvent operator at the fundamental frequency and its harmonics. (b) The singular values of the harmonic resolvent operator constructed using $N_b =3$ and $N_f =3, 7, 11$ for the cavity flow at $Ma_{\infty }= 0.6$.

To get insight into the coherent structures that get preferentially amplified through linear mechanisms selected by the two transfer functions (i.e. the classical and harmonic resolvent operators), we look into the real component of the streamwise velocity response modes in figure 6 along with the forcing mode that generates those responses. The forcing and response modes for the classical resolvent analysis correspond to the optimal singular values ($\sigma _1$) at different frequencies. These frequencies, and hence the modes, are decoupled, and each mode oscillates individually in time at the corresponding frequency. Moreover, since the modes are decoupled and follow a unit normalization, the amplitude of each mode can be arbitrary compared with each other. The response mode at frequency zero is mostly located inside the cavity, which resembles the centrifugal mode observed at low frequencies in the cavity flow (Bres & Colonius Reference Bres and Colonius2008). The spatial structure of the corresponding forcing mode inside the cavity overlaps with the spatial structures of the response mode indicating the presence of the so-called wavemaker region (Symon et al. Reference Symon, Rosenberg, Dawson and McKeon2018). At frequency $\omega _p$, the modal structure represents the unstable Rossiter mode II, with the Kelvin–Helmholtz (KH) instabilities mainly localized in the shear-layer region over the cavity. The forcing mode structures at frequency $\omega _p$ are mostly concentrated near and upstream of the cavity leading edge revealing the convective nature of the KH instabilities. The optimal classical resolvent forcing and response modes at higher frequencies (${\geq }2\omega _p$) are unphysical and provide no meaningful information. The inability to obtain a meaningful mode using the classical resolvent analysis at frequencies where the resolvent operator is not low rank as shown in figure 5(a) is a limitation that has been observed in other flows (Symon et al. Reference Symon, Sipp and McKeon2019) before. For an oscillatory flow, which is the case here, Symon et al. (Reference Symon, Sipp and McKeon2019) found the high-rank behaviour (i.e. lack of linear amplification mechanism) of the resolvent operator at harmonics of the fundamental frequency of an airfoil flow. After approximating the nonlinear forcing using the triadic interactions of a few highly amplified resolvent response modes, they obtained meaningful structures at those higher harmonics that agreed with the spectral proper orthogonal decomposition modes. In the context of a mode's interaction, we speculate that the unphysical classical resolvent modes obtained at the higher harmonics in figure 6 are due to the absence of modelling interactions among the modes at different frequencies in the classical resolvent formulation. To resolve the physical modes at those frequencies in the classical resolvent analysis framework, one might need to model the forcing, which includes nonlinear perturbation terms, using a few energetic response modes such as those at frequencies 0 and $\omega _p$. However, in the harmonic resolvent analysis, the cross-frequency interaction is embedded in the modelling framework. Indeed, we see next that it is possible to circumvent the limitation very well within the linear framework using the harmonic resolvent analysis without resorting to modelling the nonlinear forcing.

Figure 6. Real component of the streamwise velocity forcing (blue-green-yellow) and response modes (blue-red) for the cavity flow at $Ma_{\infty } = 0.6$ associated with $\sigma _1$ at each frequency in figure 5(a) for the classical resolvent analysis, and corresponding to $\sigma _1$ in figure 5(b) for the harmonic resolvent analysis with $N_b=3$ and $N_f =11$. The contours of the optimal forcing modes of the harmonic resolvent at frequencies $0, \omega _p, 2\omega _p, 3\omega _p$ and $4\omega _p$ are plotted within the values $\pm 0.26,\pm 2.44, \pm 0.76,\pm 0.24$ and $\pm 0.07$, respectively.

The streamwise velocity component of the harmonic resolvent forcing and response modes corresponding to the optimal singular value is shown in figure 6. We found that the qualitative structure of the forcing and response modes are almost identical regardless of the number of perturbation frequencies ($N_f$). Hence, we only show the modes corresponding to $N_f=11$ for brevity. Since the modes are stacked in one singular vector and orthonormal in their respective inner products, the relative amplitude of the modes at each frequency within one vector varies with the number of perturbation frequencies $N_f$. As the modes are temporally coupled, the mode amplitudes at different frequencies reveal their relative weights of being preferentially excited by an input. At frequency zero, the organization of the response mode structures is identical to the zero frequency DFT mode inside the cavity shown in figure 4(b) with a tail extending over the downstream wall. Unlike the classical resolvent mode, the optimal harmonic resolvent forcing and response mode at frequency zero do not overlap inside the cavity implying a difference in the instability mechanism revealed by both analyses. At frequency $\omega _p$, we observe similar KH mode shapes in the shear-layer region with some qualitative difference near the trailing edge and inside the cavity, compared with the classical resolvent mode, in which the instability also propagates into the cavity following the recirculation contour inside the rear cavity region. Due to the convective nature of the instability at frequency $\omega _p$, the optimal forcing is mostly located upstream of the response mode structures.

The most significant difference between the classical and harmonic resolvent response modes emerges at the higher harmonics of $\omega _p$, as evident in figure 6. The response modes at these higher frequencies resemble the compact KH wavepacket structures with smaller spatial wavelengths along the cavity shear layer, and their concentration shifts toward the trailing edge with increasing frequency. Meanwhile, the relative amplitude of each mode decreases as frequency increases. Since these instability modes at high frequencies (${\geq }2\omega _p$) are also convective, the forcing structures are mainly concentrated upstream of the response. With an increase in frequency, the most responsive region to introduce forcing shifts towards the shear-layer region after the cavity leading edge. Generation of harmonics is a nonlinear process, and by incorporating the base flow at frequency $\omega _p$ into the set about which the linearization is done, we successfully recover physical modes at those higher frequencies from the harmonic resolvent analysis. Again we stress that the cross-frequency interaction between perturbations at different frequencies through the base flow is the underlying mechanism for this outcome.

To understand the effect of the base flow frequency truncation on the cross-frequency interaction, we perform the harmonic resolvent analysis by increasing the number of base frequencies in the set $\varOmega$ from $N_b=3$ to $N_b=7$, i.e. considering the base flow frequency set $\varOmega = \{-3,-2,-1,0,1,2,3\}\omega _p$. We keep the number of perturbation frequencies to $N_f=11$ as a constant. The optimal singular value decreases slightly with the increase in $N_b$ from 3 to 5; however, it remains unchanged with a further increase of $N_b$ to 7. The sub-optimal singular values remain the same for all values of $N_b$. We found the effect of increasing $N_b$ on modal structures minimal except at high frequencies (${>}3\omega _p$) that exhibit more small-scale structures near the cavity trailing edge as discussed in Appendix D. To examine the input–output amplifications from different frequency pairs, we look into the cross-frequency interaction in the set $\tilde \varOmega$ through the block singular values of the harmonic resolvent operator (see § 2.5). We reconstruct a low-rank approximation of the harmonic resolvent operator using the singular values and the corresponding singular vectors from $k=1$ to $k=8$. Then, by performing the SVD of individual blocks of the harmonic resolvent operator, we compute the quantity $E_{j,k}$ following (3.1). The result is shown in figure 7(b), where we have used $N_b=5$ and $N_f=11$ to construct the harmonic resolvent operator. The darker colour in figure 7(b) represents the significant coupling between the pair of frequencies ($\,j\omega _p,k\omega _p$). The diagonal blocks reveal strong self-interaction at frequencies $0$ and $\omega _p$. The off-diagonal blocks show the extent of the cross-frequency interactions. The higher harmonics (${>}3\omega _p$) mostly interact with the frequency $\omega _p$. From the input–output perspective, the input at frequency $\omega _p$ can generate output at frequencies up to $5\omega _p$, as revealed by the column corresponding to $\omega _p$.

Figure 7. (a) The singular values of the harmonic resolvent operator constructed using $N_b =3, 5, 7$ and $N_f = 11$ for the cavity flow at $Ma_{\infty }= 0.6$. (b) Fractional variance $E_{j,k}$ of block singular values of the harmonic resolvent operator obtained using $N_b= 5$ and $N_f= 11$.

The analysis that was performed until now reveals the effect of truncation of both base flow and the perturbation frequencies on capturing the dominant dynamics of perturbations. The truncation of perturbation frequencies ($N_f$) affects the relative amplitude of the modes to some extent but minimally affects the spatial structures of the modes. The truncation of the frequencies in the base flow set about which linearization is performed influences the modal structures at relatively high frequencies. Consequently, the choice of truncation of base flow frequency depends on the analysis objective. If one wishes to model the structures accurately at high frequencies, including more frequencies in the base flow set will enhance the accuracy. Otherwise, if the objective is to get insight into the dominant amplified structures in the flow, linearizing about a few energetic frequencies should suffice.

Next, we analyse the dynamics of perturbations that are not harmonic to the fundamental base flow frequency. We construct the modified harmonic resolvent operator $\boldsymbol H(\mathrm {i}\gamma )$ for values of $\gamma$ in the range $(0,\omega _p/2]$ to identify the amplification of perturbations over the set of frequencies $\tilde \varOmega _{\gamma } = \gamma + \{-5,\ldots,-1,0,1,\ldots,5\}\omega _p$. The variation of the optimal singular value as a function of $\gamma$ is shown in figure 8(a). The optimal singular value peaks for $\gamma =0.45\omega _p$, indicating that significant amplification to perturbation for frequencies in the set $\tilde \varOmega _{0.45\omega _p}$ can happen. To identify the dominant frequency in the set that causes the amplification, we look into the cross-frequency amplification map quantified by $E_{j,k}$ following (3.1) in figure 8(b). We observe dominant self-interactions and strong cross-frequency interaction by the input at frequency $-0.55\omega _p$. This frequency is likely to be associated with Rossiter mode I, although it shows deviation from the theoretical prediction of frequency using (4.1). We plot the modal structures at two representative frequencies $-1.55\omega _p$ and $-0.55\omega _p$ in figure 9(a). At both frequencies, upstream-located forcing generates a response along the cavity shear layer and downstream with an overlapping region near the leading edge. It is instructive to examine the time-dependent evolution of the perturbation from the set $\tilde \varOmega _{0.45\omega _p}$ by converting the modes into the time domain according to (2.28b). The spatiotemporal modes at four instants within one base flow period are shown in figure 9(b). The figures exhibit the temporal growth of the modal structures in the shear layer, which convects downstream and hits the trailing edge. Then perturbations travel upstream following the recirculation inside the cavity, and the cycle repeats.

Figure 8. (a) Variation of the optimal singular value of the harmonic resolvent operator constructed using $N_b = 5$ and $N_f = 11$ as a function of $\gamma$. (b) Fractional variance $E_{j,k}$ of block singular values of the harmonic resolvent operator for the set of frequencies corresponding to $\gamma = 0.45\omega _p$.

Figure 9. (a) Real component of the streamwise velocity mode at frequencies $-1.55\omega _p$ and $-0.55\omega _p$ from the set of perturbation frequencies $\tilde \varOmega _{0.45\omega _p}$, and (b) temporal evolution of the streamwise velocity perturbations at four instants within one base flow period.

4.2. Cavity flow at $M_\infty =0.8$

We also apply the harmonic resolvent analysis to understand the cross-frequency interaction in the cavity flow at Mach $0.8$, which contains more than one resonance. We collect post-transient data and perform DFT to obtain the base flow $\hat {\bar {\boldsymbol {q}}}(x,y)$ as before. The frequency spectrum of the streamwise-momentum component of the base flow is shown in figure 10(a). Two coexisting Rossiter mechanisms drive the base flow oscillation where the frequency of Rossiter mode $\text {I}$ is $St_1=\omega _1 L/2{\rm \pi} u_{\infty }=0.345$ and the frequency of Rossiter mode $\text {II}$ corresponds to $St_2=\omega _2 L/2{\rm \pi} u_{\infty }=0.689$, which agrees well with oscillatory frequency based on the semi-empirical formula of Rossiter (Rossiter Reference Rossiter1964). In addition, we observe the presence of the harmonics of $St_1$ and $St_2$ in the spectrum. Unlike the cavity flow at Mach 0.6, the spectrum of the base flow at Mach 0.8 is not monochromatic and, thus, poses several ways to construct the set of the base flow frequency $\varOmega$. Here we consider three sets of base flow frequency, $\varOmega _{\text {I}} = \{-\omega _1,0,\omega _1\}$, $\varOmega _{\text {II}} = \{-\omega _2,0,\omega _2\}$ and $\varOmega _{\text {III}}=\{-\omega _2,-\omega _1,0,\omega _1,\omega _2\}$, to approximate the time-varying base flow and linearize the dynamics about those frequency sets. We remark here that in this particular flow the frequency of Rossiter mode $\text {II}$ ($\omega _2$) is approximately two times the frequency of Rossiter mode $\text {I}$ ($\omega _1$). Thus, we can combine them in the set $\varOmega _{\text {III}}$ without modifying the theoretical formulation. To perform the harmonic resolvent analysis, we choose the set of perturbation frequencies to be $\tilde \varOmega = \{-5,\ldots,\,-1,0,1,\ldots,5\}\omega _{p}$, where $\omega _p = \omega _1$ for the base flow frequency sets $\varOmega _{\text {I}}$ and $\varOmega _{\text {III}}$, and $\omega _p = \omega _2$ for the set $\varOmega _{\text {II}}$. We follow the same steps used for the cavity flow at Mach $0.6$ to construct the linear operators ${\boldsymbol{\mathsf{L}}}$ and $\hat {{\boldsymbol{\mathsf{G}}}}_k$ and to perform the SVD of the harmonic resolvent operator.

Figure 10. (a) Frequency spectrum of the streamwise momentum of cavity flow at Mach 0.8. Here, (- - -, red) indicates the Rossiter mode frequency prediction using (4.1). (b–d) The leading $10$ singular values of the harmonic resolvent operator constructed using the sets of base flow frequencies $\varOmega _{\text {I}}=\{-\omega _1,0,\omega _1\}$, $\varOmega _{\text {II}}=\{-\omega _2,0,\omega _2\}$, $\varOmega _{\text {III}}=\{-\omega _2,-\omega _1,0,\omega _1,\omega _2\}$.

The first $10$ singular values of the harmonic resolvent operator obtained by linearizing about the three base flow frequency sets are shown in figure 10(b–d). For set $\varOmega _{\text {III}}$, the singular values are computed without removing the singularity and the first two modes are associated with phase shift directions. The optimal singular values ($\sigma _1$) in all three cases are greater than the singular value $\sigma _{10}$ by approximately an order of magnitude. To analyse the cross-frequency interactions, we reconstruct the harmonic resolvent operator using singular values (and associated singular vectors) from $k= 1$ to $k= 8$ in figure 10(b,c), and from $k= 3$ to $k= 10$ in figure 10(d). The fractional variance of the block singular values $E_{j,k}$ of the harmonic resolvent operators calculated according to (3.1) is plotted in figure 11. In figure 11(a) the map shows that the self-interactions at frequencies zero and $\omega _2$ dominate the flow by amplification through the diagonal blocks whereas the cross-frequency interactions are weak. Linearizing about the base flow frequency set $\varOmega _{\text {I}}$ does not model the cross-frequency interactions effectively as Rossiter mode $\text {I}$ is not dominant in the nonlinear flow. When we linearize about the set $\varOmega _{\text {II}}$ to perform the harmonic resolvent analysis, the frequency interaction map exhibits a progressive pattern with an increase in off-diagonal blocks in figure 11(b). The map shows the amplification of perturbations at harmonics of Rossiter mode $\text {II}$ through the cascaded cross-frequency interactions. The base flow frequency set $\varOmega _{\text {III}}$ contains both Rossiter mechanisms present in the nonlinear flow, and the frequency interaction map is shown in figure 11(c). The input–output pair of frequencies ($\omega _1, \omega _1$), ($\omega _1,3\omega _1$) and ($\omega _1,5\omega _1$) are strongly coupled such that the interactions affect the perturbation amplification at odd harmonics of Rossiter mode $\text {I}$. Whereas the interactions between the frequency pairs ($\omega _1, \omega _2$) and ($\omega _1,2\omega _2$) are negligible indicating that perturbation at the Rossiter mode $\text {I}$ frequency does not interact significantly with perturbation at the Rossiter mode $\text {II}$ frequency. A similar conclusion can be drawn by observing the frequency $\omega _2$ and $2\omega _2$ rows of the map. The perturbations at the harmonics of both Rossiter modes are amplified through the interactions among corresponding Rossiter modes and the mean. Alternatively, based on the feature of the column of frequency $\omega _1$, we can also find that the input perturbation at the Rossiter mode $\text {I}$ frequency $\omega _1$ will generate a significant response at frequencies $\omega _1$, $3\omega _1$ and $5\omega _1$ but no response at the Rossiter mode $\text {II}$ frequency $\omega _2$, and vice versa.

Figure 11. Fractional variance $E_{j,k}$ of block singular values of the harmonic resolvent operator constructed using the sets of base flow frequencies (a) $\varOmega _{\text {I}}=\{-\omega _1,0,\omega _1\}$, (b) $\varOmega _{\text {II}}=\{-\omega _2,0,\omega _2\}$ and (c) $\varOmega _{\text {III}}=\{-\omega _2,-\omega _1,0,\omega _1,\omega _2\}$.