2.1. Fourier-based resolvent analysis

The non-dimensional incompressible Navier–Stokes equations are given by

(2.1a,b)\begin{equation} \frac{\partial\bar{u}_i}{\partial{t}} + \bar{u}_j\frac{\partial\bar{u}_i}{\partial{x_j}} ={-}\frac{\partial\bar{p}}{\partial x_i} + \frac{1}{{{Re}}}\frac{\partial^2\bar{u}_i}{\partial{x}_j\partial{x}_j},\quad \frac{\partial\bar{u}_i}{\partial{x_i}} = 0, \end{equation}

where $\bar {u}_i$ is the total velocity (including the mean and the fluctuating component) in the $x_i$ direction and $\bar {p}$ is the total pressure. The Reynolds number is given by ${{Re}} = u^* L^*/\nu$, where $\nu$ is the kinematic viscosity, and $u^*$ and $L^*$ are, respectively, a reference velocity and length scale used to non-dimensionalise $\bar u_i$, $x_i$ and $t$. Likewise, $\bar p$ is non-dimensionalised by a reference density $\rho ^*$ and $u^*$. The non-dimensionalisations for each of the cases studied in this work are given in table 1. The total velocity can be split into $\bar {u}_i = U_i + u_i$. Here, $U_i := \langle \bar u_i \rangle$ represents the average over ensembles and homogeneous directions, with $\langle {\cdot } \rangle$ denoting the averaging operation, and $u_i$ is the fluctuating component. Similarly, pressure can be decomposed as $\bar {p} = P + p := \langle \bar p \rangle + p$.

Table 1. Length and velocity used to non-dimensionalise the Navier–Stokes equations for each case considered in this paper.

We can split (2.1a,b) into equations for the mean and the fluctuating components of the flow

(2.2a,b) \begin{gather} \frac{\partial{U_i}}{\partial t} + \left\langle \bar{u}_j \frac{\partial\bar{u}_i}{\partial{x_j}}\right\rangle ={-}\frac{\partial{P}}{\partial{x_i}} + \frac{1}{{{Re}}} \frac{\partial^2U_i}{\partial{x_j}\partial{x_j}},\quad \frac{\partial{U_i}}{\partial{x_i}}= 0, \end{gather}

(2.3a,b)\begin{gather}\frac{\partial{u_i}}{\partial t} + U_j\frac{\partial{u_i}}{\partial{x_j}} + u_j\frac{\partial{U_i}}{\partial{x_j}} ={-}\frac{\partial{p}}{\partial{x_i}} + \frac{1}{{{Re}}} \frac{\partial^2u_i}{\partial{x_j}\partial{x_j}} + f_i,\quad \frac{\partial{u_i}}{\partial{x_i}}= 0, \end{gather}

where $f_i$ is the remaining nonlinear terms in the fluctuating equations. The equations above do not have an analytic solution unless in very particular situations and are most commonly solved numerically. Note that some of the terms in the fluctuating equations may be zero depending on the flow configuration. In traditional resolvent analysis and in this work, (2.3a,b) are Fourier-transformed in the homogeneous spatial directions. For the cases considered in this work, these are the streamwise and spanwise directions. Homogeneity in these directions also implies $\partial U_i /\partial x_1 = \partial U_i /\partial x_3 = 0$. Further simplifications of (2.3a,b) are also possible in some cases, for example when certain components of the mean flow are zero. Moreover, in traditional resolvent analysis, the mean flow is statistically stationary (and therefore independent of time), and the fluctuation equations are assumed to be periodic in time and further Fourier transformed in time.

After applying a Fourier transformation in time and in one or more spatial dimensions in which the system is homogeneous, and discretising the equations over the remaining spatial dimensions of the grid, we obtain the following:

(2.4a,b) \begin{equation} \widehat{\boldsymbol{\mathsf{D}}_t} \hat{\boldsymbol{\mathsf{u}}}_i + \widehat{\boldsymbol{\mathsf{u}}_j} \boldsymbol{\mathsf{dU}}_{i,j} + \boldsymbol{\mathsf{U}}_j \widehat{\boldsymbol{\mathsf{D}}_j} \hat{\boldsymbol{\mathsf{u}}}_i ={-}\widehat{\boldsymbol{\mathsf{D}}}_i \hat{\boldsymbol{\mathsf{p}}} + \frac{1}{{{Re}}} \hat{\boldsymbol{\mathsf{L}}} \hat{\boldsymbol{\mathsf{u}}}_i + \hat{\boldsymbol{\mathsf{f}}}_i,\quad \hat{\boldsymbol{\mathsf{D}}}_i \hat{\boldsymbol{\mathsf{u}}}_i = 0, \end{equation}

where $\widehat {\boldsymbol{\mathsf{u}}_i}$, $\widehat {\boldsymbol{\mathsf{f}}_i}$ and $\hat {\boldsymbol{\mathsf{p}}}$, respectively, denote the transformed and discretised $u_i$, $f_i$ and $p$; $\widehat {\boldsymbol{\mathsf{D}}_t}$ and $\widehat {\boldsymbol{\mathsf{D}}_i}$ are the transformed discrete differentiation operators in time and the $x_i$-directions, respectively; $\hat {\boldsymbol{\mathsf{L}}}$ is the transformed Laplacian. For systems that are homogeneous in the streamwise and spanwise directions and periodic in time, we can write $\widehat {\boldsymbol{\mathsf{D}}_t} = -\mathrm {i} \omega$, $\widehat {\boldsymbol{\mathsf{D}}_1} = \mathrm {i}k_1$, $\widehat {\boldsymbol{\mathsf{D}}_3} = \mathrm {i}k_3$, $\widehat {\boldsymbol{\mathsf{D}}_2} = \boldsymbol{\mathsf{D}}_2$ and $\hat {\boldsymbol{\mathsf{L}}} = -k_1^2 + \boldsymbol{\mathsf{D}}_2^2 -k_3^2$, where $\omega$ is the chosen frequency for the temporal Fourier transform, $(k_1, k_3)$ are the streamwise and spanwise wavenumber pair, and $\boldsymbol{\mathsf{D}}_2$ is the discrete derivative in the $x_2$ direction. We use $\boldsymbol{\mathsf{U}}_{i}$ and $\boldsymbol{\mathsf{dU}}_{i,j}$ to denote the diagonal matrices whose diagonal terms are, respectively, $U_i$ and $\partial U_i/\partial x_j$ evaluated at the grid points, and we note that since the system is homogeneous in $x_1$ and $x_3$, $\boldsymbol{\mathsf{U}}_2 = \boldsymbol{\mathsf{dU}}_{i,1} = \boldsymbol{\mathsf{dU}}_{i,3} = \boldsymbol {0}$. Each of the discretised momentum and continuity equations is an $N_{2}$-dimensional system, where $N_2$ is the spatial resolution in the $x_2$ direction. These linearised equations can then be cast as

(2.5) \begin{equation} \left[\begin{array}{@{}c@{}} \hat{\boldsymbol{\mathsf{u}}}_1(x_2)\\ \hat{\boldsymbol{\mathsf{u}}}_2(x_2)\\ \hat{\boldsymbol{\mathsf{u}}}_3(x_2)\\ \hat{\boldsymbol{\mathsf{p}}}(x_2) \end{array}\right] = \hat{\boldsymbol{\mathsf{H}}}^{(k_1,k_3,\omega)} \left[\begin{array}{@{}c@{}} \hat{\boldsymbol{\mathsf{f}}}_1(x_2)\\ \hat{\boldsymbol{\mathsf{f}}}_2(x_2)\\ \hat{\boldsymbol{\mathsf{f}}}_3(x_2)\\ \boldsymbol{0} \end{array}\right], \end{equation}

where $\hat {\boldsymbol{\mathsf{H}}}^{(k_1,k_3,\omega )}$ denotes the traditional Fourier-based resolvent operator

(2.6) \begin{align} & \hat{\boldsymbol{\mathsf{H}}}^{(k_1, k_3, \omega)} \nonumber\\ &\quad =\left[\left( -\mathrm{i} \omega \boldsymbol{\mathsf{I}} - \frac{1}{{{Re}}} \widehat{\boldsymbol{\mathsf{L}}} + \mathrm{i} k_1 \boldsymbol{\mathsf{U}}_1 + \boldsymbol{\mathsf{U}}_2 \boldsymbol{\mathsf{D}}_2 + \mathrm{i} k_3 \boldsymbol{\mathsf{U}}_3 \right) \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}} \boldsymbol{\mathsf{I}} & \boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol{\mathsf{I}} & \boldsymbol 0 & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol 0 & \boldsymbol{\mathsf{I}} & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 \end{array}\right)\right. \nonumber\\ &\qquad +\left.\left(\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}} \boldsymbol{0} & {\boldsymbol{\mathsf{dU}}}_{1,2} & \boldsymbol{0} & \mathrm{i} k_1 \boldsymbol{\mathsf{I}} \\ \boldsymbol{0} & {\boldsymbol{\mathsf{dU}}}_{2,2} & \boldsymbol{0} & \boldsymbol{\mathsf{D}}_2 \\ \boldsymbol{0} & {\boldsymbol{\mathsf{dU}}}_{3,2} & \boldsymbol{0} & \mathrm{i} k_3 \boldsymbol{\mathsf{I}} \\ \mathrm{i} k_1 \boldsymbol{\mathsf{I}} & \boldsymbol{\mathsf{D}}_2 & \mathrm{i} k_3 \boldsymbol{\mathsf{I}} & \boldsymbol 0 \end{array}\right)\right]^{{-}1}, \end{align}

and the superscript $(k_1,k_3,\omega )$ indicates the choice of streamwise and spanwise wavenumbers $k_1$ and $k_3$ and frequency $\omega$ used in the Fourier transforms. The functional dependence on $x_2$ indicates the discretisation over the wall-normal spatial dimension. Typically, the SVD of the linear resolvent operator $\hat {\boldsymbol{\mathsf{H}}}^{(k_1,k_3,\omega )} \in \mathbb{C}^{4N_2} \times \mathbb{C}^{4N_2}$ is taken to study the left and right singular vectors as response and forcing modes, and the singular values as amplification factors or gains. We denote the principal forcing and response modes by $\boldsymbol {\hat \phi }(x_2) = [ \boldsymbol {\hat \phi }_1^{\rm T}, \boldsymbol {\hat \phi }_2^{\rm T}, \boldsymbol {\hat \phi }_3^{\rm T}, \boldsymbol {0}^{\rm T}]^{\rm T}$ and $\boldsymbol {\hat \psi }(x_2) = [ \boldsymbol {\hat \psi }_1^{\rm T}, \boldsymbol {\hat \psi }_2^{\rm T}, \boldsymbol {\hat \psi }_3^{\rm T}, \boldsymbol {\hat \psi }_p^{\rm T}]^{\rm T}$, respectively. For a wall-normal spatial domain $[0,L_2]$, the modes are normalised such that their integrated kinetic energy satisfies

(2.7)\begin{gather}{[}\hat \phi_1] + [\hat \phi_2] + [\hat \phi_3] = 1, \end{gather}

(2.8)\begin{gather}{[}\hat \psi_1] + [\hat \psi_2] + [\hat \psi_3] = 1, \end{gather}

where we use $[{\cdot }] =({1}/{L_2}) \int _0^{L_2} |{\cdot }|^2\,{\rm d}\kern0.7pt x_2$ to denote the $x_2$-integrated energy.

2.2. Wavelet-based resolvent analysis

2.2.1. Formulation

To extend resolvent analysis to statistically non-stationary flows by accounting for transient behaviour in the mean flow or the fluctuations, we introduce the wavelet-based resolvent analysis. The benefit of the wavelet transform in time is that it preserves both time and frequency information. The wavelet transform projects a function onto a wavelet basis composed of scaled and shifted versions of a mother function $\eta (t)$. The transformed function depends on the scale $\alpha$ and shift $\beta$ parameters, respectively, linked to frequency and time information, whereas the Fourier transform is a function of only frequency. We propose using a wavelet transform in time to the left of the spatially Fourier-transformed equations (2.3a,b) while keeping the Fourier transform in homogeneous spatial directions. In practice, we discretise the spatially Fourier-transformed equations in both time and the wall-normal dimension, and apply a discrete wavelet transform $\boldsymbol{\mathsf{W}}$ to the left of the discrete equations. For a system that is homogeneous in the streamwise and spanwise directions, and that has been wavelet-transformed in time and Fourier-transformed in $x_1$ and $x_3$ for a streamwise and spanwise wavenumber pair $(k_1, k_3)$, we obtain

(2.9a,b)\begin{equation} \widetilde{\boldsymbol{\mathsf{D}}_t} \tilde{\boldsymbol{\mathsf{u}}}_i + \widetilde{\boldsymbol{\mathsf{u}}_j} \widetilde{\boldsymbol{\mathsf{dU}}_{i,j}} + \widetilde{\boldsymbol{\mathsf{U}}_j} \widetilde{\boldsymbol{\mathsf{D}}_j} \tilde{\boldsymbol{\mathsf{u}}}_i ={-}\tilde{\boldsymbol{\mathsf{D}}}_i \tilde{\boldsymbol{\mathsf{p}}} + \frac{1}{{{Re}}} \tilde{\boldsymbol{\mathsf{L}}} \tilde{\boldsymbol{\mathsf{u}}}_i + \tilde{\boldsymbol{\mathsf{f}}}_i,\quad \tilde{\boldsymbol{\mathsf{D}}}_i \tilde{\boldsymbol{\mathsf{u}}}_i = 0. \end{equation}

Here, $\widetilde {\boldsymbol{\mathsf{u}}_i}$, $\widetilde {\boldsymbol{\mathsf{f}}_i}$ and $\tilde {\boldsymbol{\mathsf{p}}}$, respectively, denote the transformed $u_i$, $f_i$ and $p$ discretised over $x_2$, $\alpha$ and $\beta$. We introduce $\widetilde {\boldsymbol{\mathsf{U}}_i} := \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{U}}_i \boldsymbol{\mathsf{W}}^{-1}$ and $\widetilde {\boldsymbol{\mathsf{dU}}_{i, j}} := \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{dU}}_{i, j} \boldsymbol{\mathsf{W}}^{-1}$, where $\boldsymbol{\mathsf{U}}_i$ and $\boldsymbol{\mathsf{dU}}_{i, j}$ are defined in §2.1. As before, $\boldsymbol{\mathsf{U}}_2 = \boldsymbol{\mathsf{dU}}_{i,1} = \boldsymbol{\mathsf{dU}}_{i,3} = \boldsymbol {0}$. Moreover, $\tilde {\boldsymbol{\mathsf{D}}}_1 = \hat {\boldsymbol{\mathsf{D}}}_1 = \mathrm i k_1$, $\tilde {\boldsymbol{\mathsf{D}}}_2 = \hat {\boldsymbol{\mathsf{D}}}_2 = \boldsymbol{\mathsf{D}}_2$, $\tilde {\boldsymbol{\mathsf{D}}}_3 = \hat {\boldsymbol{\mathsf{D}}}_1 = \mathrm i k_3$, $\tilde {\boldsymbol{\mathsf{L}}} = \hat {\boldsymbol{\mathsf{L}}} = -k_1^2 \boldsymbol{\mathsf{D}}_2^2 -k_3^2$ and $\widetilde {\boldsymbol{\mathsf{D}}_t} := \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{D}}_t \boldsymbol{\mathsf{W}}^{-1}$, where $\boldsymbol{\mathsf{D}}_t$ refers to the discrete time differentiation matrix. We denote the temporal resolution of our discretised system by $N_t$; the wavelet-transformed vectors will also have a resolution of $N_t$ in the space of wavelet parameters $\alpha$ and $\beta$. Similarly to the Fourier-based resolvent analysis, the equations can be written in matrix form as

(2.10) \begin{equation} \left[\begin{array}{@{}c@{}} \tilde{\boldsymbol{\mathsf{u}}}_1(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{u}}}_2(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{u}}}_3(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{p}}}(x_2,\alpha, \beta) \end{array}\right] = \widetilde{{\boldsymbol{\mathsf{H}}}}^{(k_1,k_3)} \left[\begin{array}{@{}c@{}} \tilde{\boldsymbol{\mathsf{f}}}_1(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{f}}}_2(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{f}}}_3(x_2,\alpha, \beta)\\ \boldsymbol 0 \end{array}\right], \end{equation}

where the functional dependence on $x_2$, $\alpha$ and $\beta$ represents the discretisation over the wall-normal spatial dimension, and the wavelet shifts and scales, and the wavelet-based resolvent operator $\tilde {\boldsymbol{\mathsf{H}}} \in \mathbb{C}^{4N_t\times N_2} \times \mathbb{C}^{4N_t\times N_2}$ is defined as

(2.11) \begin{align} & \tilde{\boldsymbol{\mathsf{H}}}^{(k_1, k_3)} \nonumber\\ &\quad =\left[\left(\widetilde{\boldsymbol{\mathsf{D}}_t} - \frac{1}{{{Re}}} \widehat{\boldsymbol{\mathsf{L}}} + \mathrm{i} k_1 \widetilde{\boldsymbol{\mathsf{U}}_1} + \widetilde{\boldsymbol{\mathsf{U}}_2} \boldsymbol{\mathsf{D}}_2 + \mathrm{i} k_3 \widetilde{\boldsymbol{\mathsf{U}}_3} \right) \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}} \boldsymbol{\mathsf{I}} & \boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol{\mathsf{I}} & \boldsymbol 0 & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol 0 & \boldsymbol{\mathsf{I}} & \boldsymbol 0 \\ \boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 & \boldsymbol 0 \end{array}\right)\right. \nonumber\\ &\qquad +\left.\left(\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}} \boldsymbol{0} & \widetilde{\boldsymbol{\mathsf{dU}}}_{1,2} & \boldsymbol{0} & \mathrm{i} k_1 \boldsymbol{\mathsf{I}} \\ \boldsymbol{0} & \widetilde{\boldsymbol{\mathsf{dU}}}_{2,2} & \boldsymbol{0} & \boldsymbol{\mathsf{D}}_2 \\ \boldsymbol{0} & \widetilde{\boldsymbol{\mathsf{dU}}}_{3,2} & \boldsymbol{0} & \mathrm{i} k_3 \boldsymbol{\mathsf{I}}\\ \mathrm{i} k_1 \boldsymbol{\mathsf{I}} & \boldsymbol{\mathsf{D}}_2 & \mathrm{i} k_3 \boldsymbol{\mathsf{I}} & \boldsymbol 0 \end{array}\right)\right]^{{-}1}. \end{align}

This formulation allows us to study transient flows using resolvent analysis. We denote the principal forcing and response modes obtained under this formulation by $\boldsymbol {\tilde \phi }(x_2, \alpha, \beta ) = [ \boldsymbol {\tilde \phi }_1^{\rm T}, \boldsymbol {\tilde \phi }_2^{\rm T}, \boldsymbol {\tilde \phi }_3^{\rm T}, \boldsymbol 0^{\rm T}]^{\rm T}$ and $\boldsymbol {\tilde \psi }(x_2, \alpha, \beta ) = [ \boldsymbol {\tilde \psi }_1^{\rm T}, \boldsymbol {\tilde \psi }_2^{\rm T}, \boldsymbol {\tilde \psi }_3^{\rm T}, \boldsymbol {\tilde \psi }_p^{\rm T}]^{\rm T}$, respectively. We denote their respective inverse wavelet-transforms by $\boldsymbol {\breve \phi }(x_2, t) := \boldsymbol{\mathsf{W}}^{-1} \boldsymbol {\tilde \phi }$ and $\boldsymbol {\breve \psi }(x_2, t):= \boldsymbol{\mathsf{W}}^{-1} \boldsymbol {\tilde \psi }$. These are normalised such that integrated kinetic energy satisfies

(2.12)\begin{gather} \frac{1}{T}\frac{1}{L_2} \int_0^{\rm T} \int_0^{L_2} |\breve \phi_1 |^2 + |\breve \phi_1|^2 + |\breve \phi_3|^2 \,{\rm d}\kern0.7pt x_2 \,{\rm d} t =1, \end{gather}

(2.13)\begin{gather}\frac{1}{T} \frac{1}{L_2}\int_0^{\rm T} \int_0^{L_2} |\breve \psi_1|^2 + |\breve \psi_2|^2 + |\breve \psi_3|^2 \,{\rm d}\kern0.7pt x_2 \,{\rm d} t = 1, \end{gather}

where $[0, T)$ represents the temporal domain. We denote the inverse Fourier transforms of the modes to the physical domain by $\boldsymbol {\phi }(x_1, x_2, x_3, t)$ and $\boldsymbol {\psi }(x_1, x_2, x_3, t)$, respectively.

2.2.2. Wavelet-based resolvent analysis with windowing

We can formulate a resolvent map between forcing and response at specific time shifts and scales by defining a windowed resolvent operator

(2.14) \begin{equation} \left[\begin{array}{@{}c@{}} \tilde{\boldsymbol{\mathsf{u}}}_1(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{u}}}_2(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{u}}}_3(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{p}}}(x_2,\alpha, \beta) \end{array}\right] = \boldsymbol{\mathsf{C}}\tilde{\boldsymbol{\mathsf{H}}}^{(k_1,k_3)} \boldsymbol{\mathsf{B}} \left[\begin{array}{@{}c@{}} \tilde{\boldsymbol{\mathsf{f}}}_1(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{f}}}_2(x_2,\alpha, \beta)\\ \tilde{\boldsymbol{\mathsf{f}}}_3(x_2,\alpha, \beta)\\ 0 \end{array}\right], \end{equation}

where $\boldsymbol{\mathsf{B}}$ and $\boldsymbol{\mathsf{C}}$ are windowing matrices on the forcing and response modes, respectively (Jeun, Nichols & Jovanović Reference Jeun, Nichols and Jovanović2016; Kojima et al. Reference Kojima, Yeh, Taira and Kameda2020). The windowing matrices select a subset of the full forcing and response states. For example, to select a particular scale and shift parameter $(\alpha _s,\beta _s)$ for the forcing mode, we set

(2.15)\begin{equation} \boldsymbol{\mathsf{B}} = \text{diag}(\mathbb{1}(\alpha = \alpha_s) \mathbb{1}(\beta = \beta_s)), \end{equation}

where $\mathbb{1}({\cdot })$ is an indicator function. Here, $\boldsymbol{\mathsf{B}}$ selects the relevant columns of $\boldsymbol{\mathsf{W}}^{-1}$ corresponding to the wavelet scales and shifts to which we wish to restrict our forcing. Analogously, $\boldsymbol{\mathsf{C}}$ selects the rows of $\boldsymbol{\mathsf{W}}$ corresponding to the wavelet scales and shifts to which we restrict our response modes. The SVD of the windowed resolvent operator, $\boldsymbol{\mathsf{C}}\tilde {\boldsymbol{\mathsf{H}}}^{(k_1,k_3)} \boldsymbol{\mathsf{B}}$, allows us to identify forcing and response modes restricted to a limited frequency and time interval determined by the selected wavelets.

2.2.3. Choice of wavelet basis

Wavelet transforms are not unique and are determined by the choice of the mother wavelet $\eta (t)$. The translations and dilations of a real mother wavelet are given by

(2.16)\begin{equation} \eta_{\alpha,\beta}(t) \equiv \frac{1}{\sqrt{\alpha}}\eta \left(\frac{t-\beta}{\alpha} \right), \end{equation}

where $\alpha$ and $\beta$ correspond, respectively, to the scale and shift parameters, and, respectively, represent location in the frequency and time domains. The dilations of the wavelet capture information at varying scales, and its translations capture information at different time intervals.

Consider an arbitrary square-integrable function $f(t)$. Its Fourier and wavelet transform are

(2.17)\begin{gather} \hat f(\omega) = \int_{-\infty}^{+\infty}f(t)\exp({- \mathrm{i} \omega t}) \,{\rm d} t, \end{gather}

(2.18)\begin{gather}\tilde f^{(w)} (\alpha, \beta)= \int_{-\infty}^{+\infty}f(t)\eta_{\alpha, \beta}(t)\,{\rm d} t, \end{gather}

where $\mathrm {i} = \sqrt {-1}$. In practice, we define the wavelet transform on a dyadic grid, i.e. $\alpha = 2^\ell$, $\beta = k2^\ell$ for $k,\ \ell \in \mathbb{Z}$. However, we usually do not dilate the mother wavelet indefinitely; the dilations and translations of a scaling function $\zeta (t)$ are used to capture the residual left-over from a scale-truncated wavelet expansion. We define the projection onto these functions as

(2.19a,b)\begin{equation} \tilde f^{(s)} (\alpha, \beta)= \int_{-\infty}^{+\infty}f(t) \zeta_{\alpha, \beta}(t)\,{\rm d} t,\quad\zeta_{\alpha, \beta}(t) \equiv \frac{1}{\sqrt{\alpha}} \zeta \left(\frac{t-\beta}{\alpha} \right). \end{equation}

The wavelet expansion of an arbitrary function $f(t)$ at dyadic scales is thus given by

(2.20)\begin{equation} f(t) = \sum_{\ell ={-}\infty}^L\sum_{k ={-}\infty}^{+\infty} \tilde f^{(w)} (2^\ell, 2^k ) \eta \left (\frac{t}{2^\ell}-k \right) + \sum_{k ={-}\infty}^{+\infty} \tilde f^{(s)}(2^L, 2^k) \zeta \left (\frac{t}{2^L} - k \right ), \end{equation}

where $L \in \mathbb{Z}$ represents the largest scale captured by the wavelet expansion, the $\tilde f^{(w)} (2^\ell, 2^k )$ terms approximate $f(t)$ at scales $-\infty < 2^\ell \leq 2^L$, and the $\tilde f^{(s)}(2^L, 2^k)$ terms capture the residual at scales $2^\ell > 2^L$. For each $\ell$, the projection onto $\eta ( t/2^\ell -k)$ will roughly capture a portion of the frequency content of $f$, centred in a time interval determined by $k$. Larger $\ell$ corresponds to a narrower band of frequencies closer to zero, while larger $k$ corresponds to later times. In a discretised setting, we use the finite resolution wavelet expansion, which approximates (2.20) for a discrete signal and where $2 \leq 2^\ell \leq 2^L \leq N_t$. The choice of the largest scale $L$ will depend on the band of frequencies we wish to isolate with our wavelet-based resolvent operator: narrower bands closer to zero will require larger $L$.

The wavelet and scaling function coefficients are produced by a premultiplication by the discrete wavelet transform $\boldsymbol{\mathsf{W}}$, which approximates the convolution against the wavelets and scaling functions. The choice of the wavelet and scaling function pair determines the properties of $\boldsymbol{\mathsf{W}}$. Wavelets/scaling functions of compact support in time result in banded $\boldsymbol{\mathsf{W}}$, since the convolution with these functions will also have compact support. Orthonormal wavelets/scaling functions result in a unitary $\boldsymbol{\mathsf{W}}$. Each wavelet or scaling function captures a portion of the temporal and frequency domains. There is a trade-off between precision in frequency and precision in time, i.e. one cannot find a function $\eta (t)$ that is well localised in both time and frequency (Mallat Reference Mallat2001). As two extreme examples, consider the Dirac delta centred at $t=1$, which is perfectly localised in time but with an infinite spread in frequency space, and the Fourier mode ${\rm e}^{\mathrm {i}t}$, which is perfectly localised in frequency space at $\omega = 1$ but has infinite spread in time. In the context of windowed resolvent analysis, we may wish to highlight specific bands of the frequency spectrum, or conversely, narrow bands in time, which will inform the choice of wavelet transform.

For this study, we work with the Shannon and Daubechies-16 wavelets. The Shannon wavelet is notable because it acts as a perfect bandpass filter (Mallat Reference Mallat2001; Najmi Reference Najmi2012) and covers a frequency band $N_t / 2^{\ell }([-2{\rm \pi}, -{\rm \pi} ] \cup [{\rm \pi}, 2{\rm \pi} ])$ (figure 1). Though the Shannon wavelet does not have the perfect frequency localisation provided by the Fourier transform, it allows the separation of the frequency content into distinct non-overlapping bands for different scales. One disadvantage of the Shannon wavelet is that it does not have a compact support in time and its corresponding discrete wavelet transform is dense, thus increasing the computational cost of the inversion of $\tilde {\boldsymbol{\mathsf{H}}}$ and SVD of its inverse. For problems where the sparsity of the wavelet transform is important, we use the Daubechies wavelets, which trade the perfect bandpass property in frequency domain for a compact support in time. Higher index Daubechies wavelets will have larger temporal supports and will behave closer to perfect bandpass filters. For both the wavelets described, the discrete transform matrix is unitary (Mallat Reference Mallat2001; Najmi Reference Najmi2012).

Figure 1. Shannon wavelets and scaling function in time (a) and frequency (b) domains. Daubechies-16 wavelets and scaling function in time (c) and frequency (d) domains. The functions shown are two wavelets for $\alpha = 2^L = N_t$ (black —) and $\alpha = 2^{L-3}$ (black $--$), and a scaling function for $\alpha = 2^L$ (red —) for arbitrary shift parameters.

In addition to numerical cost, we also consider the physical interpretation implied by the wavelet transform to inform our choice of wavelets. For example, in the channel flow case where critical layer dynamics are important, the Shannon wavelet would allow us to study wave speeds of interest with more precision due to its being a perfect bandpass filter (Ballouz et al. Reference Ballouz, López-Doriga, Dawson and Bae2023). Constraining the resolvent modes to a particular dilation of the Shannon wavelet would highlight waves of frequencies contained in the frequency band of the chosen wavelet. If more precision is required in the time domain to study time-localised phenomena, a wavelet that is compactly supported in time is more pertinent. This is the case, for example, in Ballouz, Dawson & Bae (Reference Ballouz, Dawson and Bae2024), which studies how a compactly supported resolvent forcing function affects a turbulent channel beyond its time support. In some cases, both properties can be combined through the use of wavelets that are compactly supported in time but whose frequency Fourier spectra decay quickly outside a certain frequency band, making them quasibandpass filters. The Daubechies or the Fejér–Korovkin wavelets satisfy these properties, and in the case of turbulent channel flow, produce similar results as when a Shannon wavelet-transform is used.

2.2.5. Choice of time differentiation matrix

The choice of the discrete time differentiation operator $\boldsymbol{\mathsf{D}}_t$ has a significant impact on the computation of resolvent modes. The sparsity of $\boldsymbol{\mathsf{D}}_t$ controls the sparsity of the resolvent operator, which heavily affects the memory and complexity requirements of the computation of the resolvent modes. However, though a sparse $\boldsymbol{\mathsf{D}}_t$ seems beneficial, it also distorts the time differentiation for high-frequency waves and can falsify the results of the SVD.

To illustrate this, we study the spectra of two time-derivative matrices, $\boldsymbol{\mathsf{D}}_{t,2}$, a second-order centred finite difference matrix, and $\boldsymbol{\mathsf{D}}_{t,F}$, a Fourier derivative matrix. The eigenvectors of both operators are the discrete Fourier modes. The eigenvalues of the Fourier derivative matrix $\boldsymbol{\mathsf{D}}_{t,F}$ are simply the wavenumbers $\mathrm {i} 2{\rm \pi} k /T$ for $k = 1,2\cdots \lceil {N_t/2}\rceil -1$, while those of the second-order centred difference matrix $\boldsymbol{\mathsf{D}}_{t, 2}$ are the modified wavenumbers $\mathrm {i} \sin (2{\rm \pi} k / N_t) N_t/T$. For a fixed wavenumber $k$, we note that

(2.21)\begin{equation} \lim_{N_t \to + \infty} \sin\left (\frac{2{\rm \pi} k}{N_t} \right ) \frac{N_t}{T} = \frac{2{\rm \pi} k} {T}, \end{equation}

and our modified wavenumber converges to the correct value. Now consider the maximum wavenumber $k = \lceil {N_t/2}\rceil - 1$. Suppose without loss of generality that $N_t$ is even so that $\lceil {N_t/2}\rceil -1 = N_t/2 -1$. The limit as the resolution in time increases becomes

(2.22)\begin{equation} \lim_{N_t \to +\infty} \sin \left ( \frac{2{\rm \pi}}{N_t} \left(\frac{N_t}{2}-1 \right ) \right ) \frac{N_t}{T} = \frac{2{\rm \pi}}{T}, \end{equation}

which does not converge to the correct wavenumber $2{\rm \pi} (N_t/2-1)/T$. Indeed, the error grows as $O(N_t)$.

To study how the numerical properties of the resolvent operators introduced in §§ 2.1 and 2.2 vary with the choice of discrete time derivative, we consider an arbitrary matrix $\boldsymbol{\mathsf{A}}$, which can represent the spatial derivatives of the linearised Navier–Stokes equations, and write the following approximation:

(2.23)\begin{align} (\boldsymbol{\mathsf{D}}_{t,2} + \boldsymbol{\mathsf{A}})^{{-}1} &= (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}} + \boldsymbol{\mathsf{D}}_{t,2} - \boldsymbol{\mathsf{D}}_{t,F})^{{-}1} \nonumber\\ &\approx (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}})^{{-}1} - (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}})^{{-}1} (\boldsymbol{\mathsf{D}}_{t,2} - \boldsymbol{\mathsf{D}}_{t,F}) (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}})^{{-}1} + \cdots. \end{align}

Thus,

(2.24)\begin{equation} \Vert (\boldsymbol{\mathsf{D}}_{t,2} + \boldsymbol{\mathsf{A}})^{{-}1} - (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}})^{{-}1}\Vert_2 \leq O(\sqrt{N_t}) \Vert (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}})^{{-}1} \Vert_2^2 + O(N_t). \end{equation}

The lack of convergence as $N_t$ increases suggests that the use of a finite difference operator rather than a Fourier derivative can significantly distort the SVD of the resolvent operator. To benefit from the advantages of a sparse temporal finite difference operator while avoiding spurious SVD modes, we propose using the windowing procedure described in § 2.2.2 to filter-out the wavelet scales associated with the high-frequency wavenumbers more susceptible to distortion. Specifically, rather than choose the windowing matrices $\boldsymbol{\mathsf{B}}$ and $\boldsymbol{\mathsf{C}}$ to highlight a physically interesting range of the frequency spectrum, we use them to exclude the frequencies above a threshold $k_{max} < N_t/2$. The maximum error between the eigenvalues of $\boldsymbol{\mathsf{D}}_{t, 2}$ and $\boldsymbol{\mathsf{D}}_{t, F}$ is given by the Taylor expansion

(2.25)\begin{equation} \sin \left(\frac{2{\rm \pi} k_{max}}{N_t}\right) \frac{N_t}{T} - \frac{2{\rm \pi} k_{max}}{T} = \frac{4{\rm \pi}^3 k_{max}^3} {3 N_t^2 T} + O \left(\frac{k_{max}^5}{N_t^4}\right). \end{equation}

Thus, assuming the chosen wavelet transform $\boldsymbol{\mathsf{W}}$ is unitary,

(2.26)\begin{equation} \Vert \boldsymbol{\mathsf{B}} \boldsymbol{\mathsf{W}} (\boldsymbol{\mathsf{D}}_{t,2} + \boldsymbol{\mathsf{A}})^{{-}1} - (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}})^{{-}1}) \boldsymbol{\mathsf{W}}^{{-}1}\boldsymbol{\mathsf{C}} \Vert_2 \leq O \left(\frac{k_{max}^{3/2}}{N_t}\right) \Vert (\boldsymbol{\mathsf{D}}_{t, F} + \boldsymbol{\mathsf{A}})^{{-}1} \Vert_2^2 . \end{equation}

The error between the SVD of the two operators will decrease as $1/N_t$ provided $k_{max}$ remains fixed. In § 3, we employ this filtering approach for cases using a finite difference time derivative operator.

3.1. Turbulent channel flow

The mean profile of turbulent channel flow at friction Reynolds number ${{Re}}_\tau \approx 186$ is obtained from Bae & Lee (Reference Bae and Lee2021). We non-dimensionalise using the channel half-height $\delta ^*$ and the friction velocity $u_\tau ^*$ as shown in table 1, so that ${{Re}} = {{Re}} _\tau$. We note that $\boldsymbol{\mathsf{U}}_{3}$, $\boldsymbol{\mathsf{U}}_{2}$, $\boldsymbol{\mathsf{dU}}_{2,2}$ and $\boldsymbol{\mathsf{dU}}_{3,2}$ are zero due to the absence of a spanwise and wall-normal contribution to the mean velocity profile.

For resolvent analysis, the wall-normal direction is discretised using a Chebyshev collocation method using $N_2 = 128$, and the mean streamwise velocity profile and its wall-normal derivative from the direct numerical simulation (DNS) are interpolated to the Chebyshev collocation points. We uniformly discretise the temporal domain, $[0, T)$, where $T = 5.5$ ($T^+ = 1023$), with a temporal resolution of $N_t = 128$. The superscript $+$ denotes wall units, which are defined to be $({\cdot })^+ := ({\cdot }) {{Re}}_\tau$ for length and time scales and $({\cdot })^+ := ({\cdot })$ for velocity scales. We impose periodic boundary conditions at the edges of the time window. The temporal boundary conditions are encoded in the choice of time differentiation matrix $\boldsymbol{\mathsf{D}}_t$, which we choose to be a Fourier differentiation matrix. For spatial derivatives in the wall-normal direction, we use first- and second-order Chebyshev differentiation matrices, and impose no-slip and no-penetration boundary conditions at the wall. In choosing $k_1$ and $k_3$, we target spanwise and streamwise wavelengths of $\lambda _1^+ \approx 1000$ and $\lambda _3^+ \approx 100$ in wall units, which are the typical length scales for near-wall streaks and correspond to the peaks in the streamwise spectrum. For our wavelet transform, we choose a two-stage ($L = 2$) Shannon wavelet transform, i.e. the state vectors in the resulting system described by (2.10) will contain terms covering the following three intervals in frequency domain: $N_t/4 [-{\rm \pi}, {\rm \pi}]$ and $N_t/2^\ell [-2{\rm \pi}, -{\rm \pi} ] \cup [-2{\rm \pi}, -{\rm \pi} ]$ for $\ell = 1, 2$. We note that, for this application, the resolvent modes converge despite the relatively low dimension of the resolvent operator. This permits us to use the aforementioned dense differentiation matrices. Sparse finite difference matrices may be used in higher-dimensional problems to improve efficiency. The wavelet- and Fourier-based cases, described by (2.9a,b) and (2.4a,b), respectively, only differ by their time differentiation matrices. These satisfy $\boldsymbol{\mathsf{W}}^{-1}\widetilde {\boldsymbol{\mathsf{D}}_t} \boldsymbol{\mathsf{W}} = \boldsymbol{\mathsf{F}}^{-1} \widehat {\boldsymbol{\mathsf{D}}}_{t, F} \boldsymbol{\mathsf{F}} = \boldsymbol{\mathsf{D}}_t$, where $\widehat {\boldsymbol{\mathsf{D}}}_{t,F}$ is the diagonal matrix with diagonal terms $(\widehat {\boldsymbol{\mathsf{D}}}_{t, F})_{\ell \ell } = -\mathrm {i} \omega _\ell = \mathrm {i} (2{\rm \pi} \ell )/T$ for $\ell = -{N_t}/{2}, \ldots, {N_t}/{2}-1$, and $\boldsymbol{\mathsf{F}}$ is the discrete Fourier transform in time. The Shannon wavelet transform is unitary, and since the SVD is unique up to multiplication by a unitary matrix, we expect the singular values of $\tilde {\boldsymbol{\mathsf{H}}}^{(k_1, k_3)}$ to be the same as those for

(3.1) \begin{equation} \hat{\boldsymbol{\mathsf{H}}}^{(k_1,k_3)} = \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}} \hat{\boldsymbol{\mathsf{H}}}^{(k_1, k_3, \omega_1)} & & & \\ & \hat{\boldsymbol{\mathsf{H}}}^{(k_1, k_3, \omega_2)} & & \\ & & \hat{\boldsymbol{\mathsf{H}}}^{(k_1, k_3, \omega_3)} & \\ & & & \ddots \\ \end{array}\right), \end{equation}

where each $\hat {\boldsymbol{\mathsf{H}}}^{(k_1, k_2, \omega _\ell )}$ is defined according to (2.6). Moreover, we expect the response and forcing modes of both systems to be related by the unitary transform given by the Fourier and inverse-wavelet transform in time, $\boldsymbol{\mathsf{FW}}^{-1}$.

The results for the Fourier-based cases were computed by applying traditional resolvent analysis at each $\omega _i$ captured by our temporal grid, while the wavelet-based resolvent modes were computed by solving the full space–time system at once. We consider the results from wavelet-based resolvent analysis to be converged: the singular values obtained for $(N_t, N_2)= (128, 128)$, and for $(N_t, N_2)= (64, 128)$ produced the same leading singular values. These are not shown in this work. As implied by (3.1), a single wavelet-based resolvent analysis will yield the modes corresponding to all time scales captured by the temporal grid. In order to associate each singular value from the wavelet-based resolvent analysis to a frequency, we Fourier-transform the corresponding response mode in time, and identify the index of the non-zero Fourier component.

In figure 2(a), we show the 10 leading singular values of the wavelet-based resolvent operator, along with the first singular value for the Fourier-based operator for different temporal Fourier parameters $\omega _\ell$. In order to associate the singular values obtained from the wavelet-based method with their frequencies, we Fourier-transform each wavelet-based response mode in time, and identify the index of the non-zero component. The first 10 singular values obtained differ by at most $3\,\%$, with the largest differing by $1.3\,\%$, which matches our expectation. The discrepancy can be explained by numerical and truncation errors. Though Shannon wavelet transforms are unitary in the continuous setting, Shannon wavelets do not have compact support in time. The discrete Shannon transform is thus not a unitary matrix due to the truncation of the wavelet in time, and exhibits a condition number of approximately $1.6$ in this case. Using wavelets that are compactly supported in time, such as the Daubechies wavelets, reduces the discrepancy, as the resulting wavelet transform matrix $\boldsymbol{\mathsf{W}}$ is a unitary operator with a condition number of 1 that better preserves the singular values of $\hat {\boldsymbol{\mathsf{H}}}^{(k_1, k_3)}.$ Increasing the time resolution also reduces the gap between the singular values. We note that due to the symmetry of the problem about the centreline, the singular values appear in equal pairs (McKeon & Sharma Reference McKeon and Sharma2010). This is only visible in figure 2(a) for $\omega = 20.56$, where the pair of singular values deviate slightly for each other due to numerical error. The modes corresponding to the pair of equivalent singular values are reflections of each other about the channel centreline.

Figure 2. (a) First 10 singular values for the wavelet-based resolvent (red) and the largest singular value of the Fourier-based resolvent operator (black) computed for each $\omega _\ell$. The vertical grey lines indicate the frequencies resolved by the chosen temporal grid, and the vertical black lines delimit the frequency band covered by each of the wavelet scales used in the wavelet expansion. (b) Magnitude of the streamwise component of the principal resolvent mode corresponding to $\omega \approx 17.14$. The red line represents $\boldsymbol{\mathsf{F}} \boldsymbol {\breve \psi }_1$, the mode obtained from wavelet-based resolvent analysis after a Fourier transformation in time, and the black line the mode obtained from traditional Fourier-based resolvent analysis $\boldsymbol {\hat \psi }_1$. In both cases, $\lambda _1^+ = 1000$ and $\lambda _3^+ = 100$.

In figure 2(b), we compare the two methods further by plotting the streamwise component of the most amplified resolvent response mode that they each produce. For the Fourier-based method, this corresponds to the frequency $\omega \approx 17.14$. For the wavelet-based method, we first Fourier-transform the principal mode in time, and though not shown, observe that the Fourier mode associated with $\omega \approx 17.14$ is the only non-zero component. This matches our expectation that the wavelet-based resolvent analysis in this case is equivalent to performing the traditional Fourier-based resolvent analysis for each frequency $\omega _\ell$ individually. Figure 2(b) shows that the modes from the two methods match. Despite the slight discrepancy in the singular values, both methods yield the same resolvent modes associated with the maximum singular value, indicating that both methods are equivalent for this stationary case.

Although not shown, the streamwise component of the modes dominate for the principal modes computed with the two methods, and the modes form alternating low- and high-speed streamwise streaks. The principal forcing mode is in the form of streamwise rolls, with a negligible streamwise component. The shape of the modes is thus in line with previous analysis of the self-sustaining process of wall turbulence (Jiménez & Moin Reference Jiménez and Moin1991; Hamilton et al. Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997; Jiménez & Pinelli Reference Jiménez and Pinelli1999; Schoppa & Hussain Reference Schoppa and Hussain2002; Farrell, Gayme & Ioannou Reference Farrell, Gayme and Ioannou2017; Bae et al. Reference Bae, Lozano-Durán and McKeon2021). We do, however, note that the response modes peak at $x_2 \approx 0.21$, corresponding to $x_2^+ \approx 40$, higher than the location of $x_2^+ \approx 15$ preferred by near-wall streaks. This is also in line with the previous literature on Fourier-based resolvent analysis (Schmid, Henningson & Jankowski Reference Schmid, Henningson and Jankowski2002; McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017, Reference McKeon2019): resolvent response modes, both formulated traditionally and using the wavelet basis in time, peak at the critical layer located where $U(x_2) = \omega / k_1$. Using a version of the linearised Navier–Stokes equations that includes eddy viscosity (Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023) causes the modes to peak closer to $x_2^+ = 15$ so that they match observations of the near-wall cycle.

3.2. Transient growth mechanism of turbulent channel flow

The added advantage of the wavelet-based method lies in its ability to preserve temporal localisation. The states in (2.14) encode time and frequency information, which allows us to study transient problems even when the mean profile is statistically stationary. One such transient phenomenon is the Orr mechanism, a linear mechanism first described by Orr (Reference Orr1907) that has been proposed to explain transient energy amplification in shear flows (Landahl Reference Landahl1975; Jiménez Reference Jiménez2013, Reference Jiménez2015, Reference Jiménez2018; Encinar & Jiménez Reference Encinar and Jiménez2020). A two-dimensional physical description is given in Jiménez (Reference Jiménez2013, Reference Jiménez2018): the mean shear profile rotates backward-tilting velocity structures forward (in the positive $x_1-$ direction), effectively extending the wall-normal distances between structures; to compensate, continuity will impose larger wall-normal fluxes, i.e. larger wall-normal velocity perturbations. This effect amplifies the velocity perturbations until the velocity structures are tilted past the normal to the wall, after which the mechanism is reversed and the perturbations are attenuated. The tilting angle of the velocity perturbations does not affect the streamwise or spanwise velocities directly, and after those components are amplified by the growth of the wall-normal component through lift-up, they decay slowly through viscous effects (Jiménez Reference Jiménez2013). Thus, due to the long-lasting effects on the streamwise component of a brief amplification of the wall-normal velocity perturbation, it has been proposed that these Orr bursts are involved in the generation of persistent near-wall streaks (Jiménez Reference Jiménez2013, Reference Jiménez2015; Encinar & Jiménez Reference Encinar and Jiménez2020).

Because of its possible role in regenerating near-wall streaks, the Orr mechanism in linearised wall-bounded flows has been examined in Jiménez (Reference Jiménez2013, Reference Jiménez2015, Reference Jiménez2018) and Encinar & Jiménez (Reference Encinar and Jiménez2020). These studies rely on computing optimal growth trajectories for the linearised equations with viscosity as well as solutions to the inviscid linearised equations The optimal trajectories are defined as those emanating from the optimal initial condition which maximises the growth of kinetic energy under the linearised dynamics (Butler & Farrell Reference Butler and Farrell1993; Schmid et al. Reference Schmid, Henningson and Jankowski2002). The linear trajectories exhibit the characteristic forward tilting of velocity structures in conjunction with the transient amplification of velocity perturbations, suggesting that the Orr-mechanism is a dominant energy amplification mechanism in the linearised system.

Optimal growth trajectories compute the singular modes for the linearised flow map between an initial condition and velocity perturbations at a later time; the question we wish to answer is whether optimal external forcing upon the linearised system, which could originate from nonlinear interactions, also exploits the Orr mechanism. For this, we use the wavelet-based resolvent analysis formulation. We note that traditional resolvent analysis has been used to reveal some evidence of the Orr mechanism in turbulent jets, where it is identified by the tilt of the optimal forcing structures against the jet shear (Tissot et al. Reference Tissot, Zhang, Lajús, Cavalieri and Jordan2017; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020). In an attempt to capture the Orr mechanism in channel flow at $Re_\tau = 2000$ as in Encinar & Jiménez (Reference Encinar and Jiménez2020), we use the mean profile for channel flow at $Re_\tau = 2000$ (Hoyas & Jiménez Reference Hoyas and Jiménez2008), the same grid in the wall-normal direction as § 3.1 with $N_2 = 128$ collocation points, and a uniform temporal grid with $N_t = 256$. As argued in Encinar & Jiménez (Reference Encinar and Jiménez2020), we choose spatial wavelengths $\lambda _1 = \lambda _3 = 0.7$, and we choose a total time of $T = 2.15$.

The choice of $T$ depends on the critical layer dynamics that we wish to highlight. From traditional resolvent analysis applied to turbulent stationary channel flow, we know that the resolvent Fourier modes tend to peak in magnitude at the critical layer, i.e. where $U(x_2) = \omega / k_1$ (Schmid et al. Reference Schmid, Henningson and Jankowski2002; McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017, Reference McKeon2019). In this section, we study the effect of a time-localised forcing on a region $x_2^+ \in [0, 500]$ which contains the inner region of the boundary layer (Hoyas & Jiménez Reference Hoyas and Jiménez2006). This maps to a frequency interval $[0 ,U(x_2^+ = 500)] k_1 = [0, \omega _{max}] \approx [0, 185]$ according to the mean profile for turbulent channel flow (figure 3). Thus, using a total time $T = 2.15$ and a Shannon wavelet transform with $L = 1$, the Shannon scaling functions $\zeta (t/2^L - k)$ would cover the frequency interval $N_t/(2T)[-{\rm \pi}, {\rm \pi}] \approx [-187, 187]$. To restrict our forcing term to this frequency band while localising it time time, we use the windowed wavelet-based resolvent analysis framework from § 2.2.2. We set $\boldsymbol{\mathsf{C}}$ to the identity matrix, allowing the response modes to cover the entire time and frequency range. We choose $\boldsymbol{\mathsf{B}}$ to select the forcing terms corresponding to the relevant the Shannon scaling function. Without loss of generality, we select the shift parameter $\beta = 0$ so that the forcing term is concentrated at a time interval centred at $t=0$.

Figure 3. Mean streamwise velocity profile for channel flow at ${{Re}}_\tau = 2000$. The dashed line corresponds to the upper bound of the frequency band covered by the chosen wavelet, mapped to streamwise velocities using $U_1 = \omega /k_1$.

The principal resolvent response mode obtained from the SVD of $\tilde {\boldsymbol{\mathsf{H}}}\boldsymbol{\mathsf{B}}$ represents the maximally amplified response to a transient forcing term aligned with the selected wavelet, under the dynamics of the linearised Navier–Stokes equations. The resulting principal response mode is confined to the frequency band determined by the forcing, as shown in figure 4(a), which is expected since the time scales are decoupled in resolvent analysis for statistically stationary flows (3.1). The velocity fluctuations can only extract energy from the mean flow in the linearised formulation of resolvent analysis, and since the mean profile is solely constituted by a steady-state component, the response mode can only contain the frequencies injected by the forcing mode. We observe that the spatially integrated energy of the response mode first grows transiently, peaks at $t = 0.72$, then decays, as shown in figure 4(b). This transient growth can be explained by the non-normality of the linearised system (Schmid et al. Reference Schmid, Henningson and Jankowski2002). The response modes at three different times are shown in figure 5. We note that the modes are concentrated in a region $x_2 < 0.4$ with a corresponding Corrsin shear parameter of approximately 10, which justifies the use of a linearised method to analyse the full system (Jiménez Reference Jiménez2013).

Figure 4. (a) Inverse wavelet-transformed principal response mode for channel flow at ${{Re}}_\tau = 2000$ and $k_1 = k_3 \approx 8.98$ under windowed forcing, in the frequency–time plane. The forcing window is highlighted in red. (b) The $x_2$-integrated total energy of the inverse wavelet-transformed principal response mode.

Figure 5. Principal response mode for the channel flow at ${{Re}}_\tau =2000$ and $k_1 = k_3 \approx 8.98$ under transient forcing. The modes are shown at (a) $t = 0.25$, (b) $t = 0.71$ and (c) $t = 1.56$. The panels (i) and (ii), respectively, correspond to half-wavelength locations of $x_3/\delta = {\rm \pi}/k_3 \delta$ and $x_1/\delta = {\rm \pi}/k_1 \delta$.

The transient behaviour of the modes displays characteristics of the Orr mechanism, mainly a synchronisation between the amplification of the wall-normal component of the response mode and its forward tilting. To study the forward tilting of velocity structures in the response mode, we define the tilt angles as in Jiménez (Reference Jiménez2015), i.e.

(3.2)\begin{equation} \theta_i(x_2, t) ={-}\tan^{{-}1} \left( \frac{\partial_{x_2} \angle \boldsymbol{\breve \psi}_i (x_2, t)}{k_1} \right), \end{equation}

where $\angle ({\cdot })$ represents the complex angle. In Jiménez (Reference Jiménez2015), the angle defined above is averaged over a region of interest. We define the energy-weighted average tilt as

(3.3)\begin{equation} \theta_i^{[y^+_a,y^+_b]}(t) = \frac{\displaystyle \int_{y_a}^{y_b} \Vert \boldsymbol{\breve \psi}_i \Vert^2 \theta_i \,{\rm d}\kern0.7pt x_2}{\displaystyle \int_{y_a}^{y_b} \Vert \boldsymbol{\breve \psi}_i \Vert^2 \,{\rm d}\kern0.7pt x_2}, \end{equation}

and the amplitude as

(3.4)\begin{equation} A_i^{[y^+_a,y^+_b]} = \left(\int_{y_a}^{y_b} \Vert \boldsymbol{\breve \psi}_i \Vert^2 \,{\rm d}\kern0.7pt x_2\right)^{{1}/{2}}, \end{equation}

and pick $y_a^+ = 0$ and $y_b^+ = 2000$ to capture the half-channel. The results (figure 6) show that the amplitude of the wall-normal velocity component of $\boldsymbol {\breve \psi }_2$ indeed peaks roughly when $\theta ^{[0, 2000]}_2 \approx 0$ at $t \approx 0.45$, and decays as $\theta _2^{[0,2000]}$ tilts past zero until it vanishes for $\theta _2^{[0,2000]} = {\rm \pi}/4$.

Figure 6. Mode angle (a) and magnitude (b) of the principal response mode as a function of time. The streamwise component of the mode is represented by $-\bullet$ (blue), the wall-normal component by $-\blacksquare$ (green) and the spanwise component by $- \blacktriangle$ (red).

We note that the spanwise component tilts forward much faster than the streamwise and wall-normal components (figure 6). This faster tilt has been attributed to the spanwise component's placement closer to the wall where the shear is stronger (Encinar & Jiménez Reference Encinar and Jiménez2020), and we do indeed find that $\boldsymbol {\breve \psi _3}$ is closer to the wall compared with $\boldsymbol {\breve \psi _2}$. However, we also find that $\boldsymbol {\breve \psi _1}$ is located at the same wall-normal height as $\boldsymbol {\breve \psi _3}$, and hypothesise that it tilts forward at the same rate as $\boldsymbol {\breve \psi _2}$ because of the more direct coupling between the two components via the $u_2 \,{\rm d} U_1/{\rm d}\kern0.7pt x_2$ term in the linearised momentum equations.

Moreover, both the streamwise and spanwise components of the mode peak at $t \approx 0.75$, i.e. after $\theta _1^{[0,2000]}$ and $\theta _3^{[0,2000]}$ cross the zero threshold and approximately $t \approx 0.35$ later than the amplitude peak for the wall-normal component. The streamwise and spanwise components also decay more slowly than the wall-normal component, even as their tilt is relatively constant at $\theta _1^{[0, 2000]} = \theta _3^{[0,2000]} = {\rm \pi}/4$ and after the wall-normal component vanishes. Jiménez (Reference Jiménez2013) explains this using the Squire equation, the linearised advection–diffusion equation for the wall-normal vorticity in which the wall-normal velocity perturbation acts as a forcing term. Even if this forcing disappeared, the streamwise and spanwise velocity components would continue advecting downstream and decay only due to viscous dissipation, which acts at a slower time scale than the mean shear in the near-wall region for the length scales considered. The non-normality of the Squire system further explains the delayed growth of the spanwise component after being forced by the wall-normal component. Note that the amplitude streamwise component rises in tandem with the wall-normal component, possibly due their tighter coupling via the presence of the $u_2 \,{\rm d} U_1/{\rm d}\kern0.7pt x_2$ term in the linearised streamwise momentum equation. The tilt of the streamwise and spanwise components does not directly drive or suppress their amplitudes, since forward-tilting only affects their wall-normal gradients which do not appear in the equations of motion. This allows them to grow even after they attain their maximum positive tilt.

The lifetime of the wall-normal component $\boldsymbol {\breve \psi }_2$ of the principal resolvent response mode differs significantly from the bursting time scales in Jiménez (Reference Jiménez2013, Reference Jiménez2015) and Encinar & Jiménez (Reference Encinar and Jiménez2020), which argue that the Orr-mechanism is a linear inviscid process whose period scales with the local mean shear. Defining the location of the wall-normal velocity fluctuation as its energy-weighted centre of gravity

(3.5)\begin{equation} x_{2, g} := \frac{\displaystyle\int_{0}^{1} \Vert u_2 (t_{max})\Vert^2 x_2 \,{\rm d}\kern0.7pt x_2}{\displaystyle \int_0^{1}\Vert u_2 (t_{max}) \Vert^2 \,{\rm d}\kern0.7pt x_2}, \end{equation}

where $t_{max}$ is the time of maximum amplitude for $u_2$, we find that $x_{2, g} \approx 0.21$, farther away from the wall than the inviscid and optimal growth solutions in Jiménez (Reference Jiménez2013) and Encinar & Jiménez (Reference Encinar and Jiménez2020) located at $x_{2, g} \approx 0.16 \lambda _1 \approx 0.11$. Defining a shear time scale $S^{-1}$ as

(3.6)\begin{equation} S := \frac{{\rm d} U_1}{{\rm d}\kern0.7pt x_2}(x_2 = x_{2, g}), \end{equation}

we find that the mode grows and decays in $tS \approx 13$ local shear units, much more slowly than the inviscid and optimal growth solutions, which grow and decay in $tS \approx 1$ and $tS \approx 3.5$, respectively. Moreover, the Orr mechanism studied in Encinar & Jiménez (Reference Encinar and Jiménez2020) is linked more broadly to energetic bursting events, which have been considered in Jiménez (Reference Jiménez2015), and which exhibit a bursting period of $t = 0.2$ for the wall-normal mode corresponding to $(\lambda _x, \lambda _z) = (0.7, 0.7)$, again faster than the growth time scale for $\boldsymbol {\breve \psi }_2$. In an attempt to explain the discrepancies between the optimal growth solution in Encinar & Jiménez (Reference Encinar and Jiménez2020) and the principal resolvent mode, we note that the measures of optimality differ between the optimal growth and wavelet-based resolvent analysis frameworks; while the optimal growth framework measures energy amplification as a ratio between the initial condition and the solution at a given time, wavelet-based resolvent analysis maximises the integrated kinetic energy of the entire time interval considered. This may better capture energetic structures that persist in time. Moreover, resolvent analysis additionally provides the optimal forcing that produces the response in the velocity field, which is helpful in studying mechanisms that drive linear transient growth in turbulent flows.

However, the bursting time scales of the principal resolvent mode agree with those found using turbulent channel flow data in Encinar & Jiménez (Reference Encinar and Jiménez2020). Using an average shear time scale defined as

(3.7)\begin{equation} S_\varLambda := \frac{(U_1(x_2 = y_b) - U_1(x_2 = y_a))}{y_b - y_a},\quad y_a = 0.7,\enspace y_b = 0.25, \end{equation}

where $\varLambda = [y_a, y_b]$ captures a section of the logarithmic region, we find that $\boldsymbol {\breve \psi _2}$ takes $t S_\varLambda \approx 1.87$ average shear units to grow to half its maximum amplitude, compared with $t S_\varLambda \approx 1$ in Encinar & Jiménez (Reference Encinar and Jiménez2020). The amplitude and time scale of the process captured by $\boldsymbol {\breve \psi }$ can be further compared with other studies of bursting in the near-wall region of turbulent channel flow at moderate ${{Re}}_\tau$. In Flores & Jiménez (Reference Flores and Jiménez2010), turbulent energy is found to peak at time intervals of approximately $6 x_2$, and given that our resolvent modes are centred at $x_2 \approx 0.2$, the trend would predict a bursting period of $t = 1.2$. This roughly matches the time scale $t \approx 1.5$ for the growth and decay of $\boldsymbol {\breve \psi _2}$. Likewise, in Hwang & Bengana (Reference Hwang and Bengana2016), which studies bursting in the logarithmic region of the channel, the relation $t = 2 \lambda _3$ is found to describe the bursting period. For our case where $\lambda _3 = 0.7$, this would correspond to a bursting period of $t = 1.4$, again matching the lifetime of $\boldsymbol {\breve \psi _2}$. Thus, the Orr bursting period predicted by $\boldsymbol {\breve \psi _2}$ roughly agrees with the bursting time scale detected in turbulent channel flow, though this time scale tends to be significantly longer than the shear time scale found to govern the Orr mechanism in linear inviscid and linear optimal growth analyses.

4.1. Turbulent Stokes boundary layer

The Stokes boundary layer is simulated through a channel flow with the lower and upper walls oscillating in tandem at a velocity of $U^*_w(t) = U^*_{max} \cos (\varOmega ^* t^*)$ with no imposed pressure gradient. We non-dimensionalise velocities by $U^*_{max}$ and lengths by $\delta ^*_\varOmega := \sqrt {2\nu /\varOmega ^*}$, which denotes the laminar Stokes boundary layer thickness. Though time and frequency are both non-dimensionalised with $U^*_{max}$ and $\delta ^*_\varOmega$, we will use $t\varOmega$ as our preferred time variable for a clearer comparison with the period, and $\omega /\varOmega$ as our preferred frequency variable as it represents temporal wavenumber in this case. The relevant non-dimensional number is ${{Re}}_\varOmega = U^*_{max} \delta _\varOmega ^* / \nu$. For the current case, we consider ${{Re}}_\varOmega = 1500$, which lies within the intermittently turbulent regime (Hino, Sawamoto & Takasu Reference Hino, Sawamoto and Takasu1976; Akhavan, Kamm & Shapiro Reference Akhavan, Kamm and Shapiro1991; Verzicco & Vittori Reference Verzicco and Vittori1996; Vittori & Verzicco Reference Vittori and Verzicco1998; Costamagna, Vittori & Blondeaux Reference Costamagna, Vittori and Blondeaux2003). This problem has been well-studied numerically and experimentally in the literature (Von Kerczek & Davis Reference Von Kerczek and Davis1974; Hino et al. Reference Hino, Sawamoto and Takasu1976; Jensen, Sumer & Fredsøe Reference Jensen, Sumer and Fredsøe1989; Spalart & Baldwin Reference Spalart and Baldwin1989; Akhavan et al. Reference Akhavan, Kamm and Shapiro1991; Sarpkaya Reference Sarpkaya1993; Blondeaux & Vittori Reference Blondeaux and Vittori1994; Verzicco & Vittori Reference Verzicco and Vittori1996; Vittori & Verzicco Reference Vittori and Verzicco1998; Costamagna et al. Reference Costamagna, Vittori and Blondeaux2003; Carstensen, Sumer & Fredsøe Reference Carstensen, Sumer and Fredsøe2010; Ozdemir, Hsu & Balachandar Reference Ozdemir, Hsu and Balachandar2014).

To generate the mean profile and second-order statistics, we run a DNS using a second-order staggered finite-difference (Orlandi Reference Orlandi2000) and a fractional-step method (Kim & Moin Reference Kim and Moin1985) with a third-order Runge–Kutta time-advancing scheme (Wray Reference Wray1990). Periodic boundary conditions are imposed in the streamwise and spanwise directions and the no-slip and no-penetration boundary conditions are used at the top and bottom walls. The code has been validated in previous studies in turbulent channel flows (Bae et al. Reference Bae, Lozano-Durán, Bose and Moin2018, Reference Bae, Lozano-Durán, Bose and Moin2019; Lozano-Durán & Bae Reference Lozano-Durán and Bae2019) and flat-plate boundary layers (Lozano-Durán, Hack & Moin Reference Lozano-Durán, Hack and Moin2018), though we note that, for this problem, we modify the boundary conditions to accommodate the oscillating walls. The domain size of the channel for the DNS is given by $6{\rm \pi} \times 80 \times 3{\rm \pi}$. The domain is discretised uniformly in the $x_1$- and $x_3$-directions using $64$ points, which corresponds to non-dimensionalised spacings of $\Delta x \approx 0.29$ and $\Delta x_3 \approx 0.15$. For the $x_2$-direction, a hyperbolic tangent grid with $385$ points is used, resulting in $\min (\Delta x_2) \approx 0.01$ and $\max (\Delta x_2) \approx 0.91$. We compute the mean velocity profiles by averaging in homogeneous directions and phase. Figure 7 shows the mean and the streamwise root-mean-square (r.m.s.) velocity profiles at different times. Using $\varOmega$ to denote the non-dimensionalised wall oscillation frequency, we note that $U_1(t\varOmega + {\rm \pi}) = -U(t\varOmega )$ and $U_{i,rms}(t\varOmega + {\rm \pi}) = U_{i,rms}(t\varOmega )$. We observe that the turbulent energy peak occurs near the wall at $x_2 = 1.43$ and $t\varOmega = 2.65$, and propagates away from the wall thereafter.

Figure 7. (a) Mean streamwise velocity profile and (b) streamwise r.m.s. velocity from $t\varOmega = 0$ (blue) to $t\varOmega = {\rm \pi}$ (red). The profiles shown are at $t\varOmega = n {\rm \pi}/8$, $n = 0, 1,\ldots 8$.

To construct the resolvent operator, we first choose the spatial scales for the homogeneous directions. Using the DNS data, we calculate the streamwise energy spectrum at $x_2= 1.43$ and $t\varOmega = 2.65$, the wall-normal location and phase of the peak $U_{1,rms}$. The most energetic streamwise and spanwise scales at that location are $k_1 = 0.67$ and $k_3 = 4.22$, which we choose as the streamwise and spanwise scales for the resolvent operator. To solve the discrete system, we use a Chebyshev grid in the wall-normal direction, with $N_2 = 80$, and a uniform temporal discretisation over one period $T\varOmega = 2{\rm \pi}$, with $N_t = 1600$. As in the previous section, $\boldsymbol{\mathsf{U}}_{3}$, $\boldsymbol{\mathsf{U}}_{2}$, $\boldsymbol{\mathsf{dU}}_{2,2}$ and $\boldsymbol{\mathsf{dU}}_{3,2}$ are zero.

We choose the first-order derivative matrix in time $\boldsymbol{\mathsf{D}}_t$ and the wall-normal spatial derivative matrices to be second-order-accurate centred finite difference matrices. We additionally choose $\boldsymbol{\mathsf{D}}_t$ to be circulant to enforce periodicity in time. We compute the modes for the half-channel and enforce a no-slip and no-penetration boundary condition at the wall, and a free-slip and no-penetration boundary condition at the centreline. Because $\boldsymbol{\mathsf{D}}_t$ is a finite difference matrix rather than a Fourier differentiation matrix, we must implement a filtering step, detailed in § 2.2.5, to exclude the high temporal wavenumbers. To apply this filtering step, we must assume that the high-frequency waves are not physically significant for the turbulent Stokes boundary layer problem. We use a two-stage Daubechies-16 wavelet transform, which is a sparse unitary operator. We note that the Daubechies-16 operator is not a perfect bandpass filter, and the numerical filtering operation simply attenuates the high-frequency waves that produce spurious SVD modes instead of excluding them outright. Nevertheless, due to the high dimensionality of the problem, it remains advantageous to use sparse transforms. We choose to constrain the forcing and response modes to the scaling functions and their shifts, which roughly cover the first quarter of all temporal wavenumbers $k_t = 0\ldots, N_t/8$.

We compare the results obtained with the wavelet-based resolvent modes with the results from harmonic resolvent analysis (Padovan et al. Reference Padovan, Otto and Rowley2020). The latter computes a Fourier-based resolvent analysis simultaneously for multiple temporal wavenumbers and includes the interactions between them as they are coupled by the temporally evolving mean profile. For the harmonic resolvent analysis, we use the same Chebyshev grid as in the wavelet-based method, with $N_2 = 80$. For the sake of comparing with the wavelet-based method and to account for the filtering step, we choose a frequency resolution of $N_t = 1600/4 = 400$. We expect the two methods to produce similar singular values and modes. The singular values and modes would be equivalent in both cases if we use a Fourier differentiation operator for the wavelet-based method as in § 3.1.

The modes obtained from harmonic resolvent analysis agree well with those obtained from wavelet-based resolvent analysis. They occur at the same $x_2$ location, and time (figure 8a,b), and exhibit roughly the same frequency content (figure 8c). Moreover, the SVD of the wavelet-based and harmonic resolvent operators yield very similar singular values. The first twenty singular values are shown in figure 8(d). Despite Daubechies-16 wavelets being imperfect bandpass filters, filtering-out high-frequency waves using the sparse wavelet transform succeeds in producing resolvent modes that match the leading modes from harmonic resolvent analysis. We consider that the results for harmonic resolvent analysis shown in figure 8 are converged: the ones that use a coarser grid with $N_2 = 80$ and $N_t = 300$ produced the same 10 leading singular values up to $10^{-4}$. Moreover, the windowed wavelet-based resolvent operator exhibits significant sparsity and can be analysed efficiently, despite the larger dimension of the system. Indeed, the harmonic and wavelet-based resolvent analyses run in similar wall times.

Figure 8. Magnitude contours ($25\,\%$, $50\,\%$, $75\,\%$ and $90\,\%$ of the maximum value) of (a) the wall-normal component of the principal resolvent forcing mode and (b) the streamwise component of the principal resolvent response mode for the turbulent Stokes boundary layer; (c) the $x_2$-integrated Fourier spectrum in time for the principal response modes; (d) singular values from the SVD of the resolvent operators. Results from harmonic resolvent analysis are shown in red, and those from wavelet-based resolvent analysis in black.

The principal input and output modes corresponding to the chosen spatial scales and boundary conditions in time are shown in figure 9(a). We observe that the modes are located at roughly the same wall normal height as the peaks in $U_{1, rms}$, i.e. $x_2 \approx 0.5$. We also observe that the principal input and output modes are synchronised with the peaks in $U_{1, rms}$, though the modes tend to peak slightly earlier than $U_{1, rms}$. This suggests that linear amplification might provoke the transition to turbulence. The preliminary results presented in figure 9 suggest that the study of energy amplification in the Stokes boundary layer is a good candidate for the use of linearised methods, similarly to the turbulent channel flow (Jiménez Reference Jiménez2013).

Figure 9. Real part of (a) $\breve \phi _1$, and (b) $\breve \psi _1$ for the turbulent Stokes boundary layer. The black contour lines are $U_{1, rms}$ with the levels indicating $50\,\%, 75\,\%, 95\,\%$ of its maximum value. The vertical dashed lines show the times of the amplitude peak for the input mode ($t\varOmega = 2.11$) and output mode ($t\varOmega = 2.18$).

We also observe that the principal input mode precedes the principal output mode in time, with the peak of the former occurring $\Delta t\varOmega \approx 0.063$ before the peak of the latter. Wavelet-based resolvent analysis is able to capture the natural response time between forcing and response terms under the dynamics of the linearised Navier–Stokes equations. This time delay is also in line with a physical interpretation of the modes in which the input modes cause the output modes and must thus occur earlier. The extent to which this captures important causal mechanisms within the full nonlinear system is yet to be determined. In future works, it would be interesting to project flow fields onto these time-separated resolvent forcing and response modes to test whether better correlations can be obtained between them in the transformed bases.

Additionally, we notice that the only non-zero wavelet coefficients of both the principal forcing and response modes are those corresponding to the bottom quarter of the set of resolved frequencies i.e. the lowest four bands in the scalograms shown in figure 10(a). This validates the windowing step described above. We also see in figure 10(a) that the frequency content of the principal modes varies with time. The principal forcing mode is initially composed of lower-frequency waves, whose frequencies are centred in a band $[0, 25\varOmega ]$; these waves are gradually shifted up to frequencies centred in $[25\varOmega, 50\varOmega ]$. Likewise, the waves composing the principal response mode, initially at frequencies centred in $[25\varOmega, 50\varOmega ]$ are also shifted up to higher frequencies. We propose that this frequency shift is due to the time-varying mean streamwise velocity $U_1$, which acts as a convection velocity and accelerates the resolvent forcing and response waves. Thus, we expect the frequency content of the forcing and response modes to vary in tandem with the changing mean streamwise velocity profile. We define the average location of the streamwise modes as

(4.1)\begin{equation} x^{resp.}_{2, {avg}} := \frac{\displaystyle\int_0^{\rm T} \int_0^1 x_2 |\breve \psi_1|^2 \,{\rm d}\kern0.7pt x_2 \,{\rm d} t }{ \displaystyle\int_0^{\rm T} \int_0^1 |\breve{\psi}_1|^2 \,{\rm d}\kern0.7pt x_2 \,{\rm d} t}, \end{equation}

and

(4.2)\begin{equation} x^{forc.}_{2, {avg}} := \frac{\displaystyle\int_0^{\rm T} \int_0^1 x_2 |\breve \phi_1|^2 \,{\rm d}\kern0.7pt x_2 \,{\rm d} t}{ \displaystyle\int_0^{\rm T} \int_0^1 |\breve{\phi}_1|^2 \,{\rm d}\kern0.7pt x_2 \,{\rm d} t}, \end{equation}

and plot the frequency shift due to the mean convection $U_1k_1/\varOmega$ at the average mode locations in figure 10(a). We observe a good correlation between the shift in the frequency content of the forcing and response modes and the change in the mean velocity.

Figure 10. (a) Principal (i) forcing and (ii) response modes in the time-frequency plane; the blue line indicates $U_1k_1/\varOmega$ at $x_{2, {avg}}^{resp.} = x_{2, {avg}}^{forc.} \approx 0.46$; (b) the $x_2$-integrated frequency content in the streamwise component of the forcing (red) and response (black) modes.

We can also use the changing mean velocity profile to explain the difference in the frequency content between the forcing and response modes. We propose that this difference in frequency content is due to the different peaking times of the forcing and response modes: since the modes occur at difference phases of the oscillating mean profile, they will be convected at different velocities. To verify this, we first Fourier-transform $\boldsymbol {\breve \phi }$ and $\boldsymbol {\breve \psi }$ in time to extract their frequency content with better precision, and observe in figure 10(b) that the average frequency shift between the forcing and response modes is

(4.3)\begin{equation} \Delta \omega := \frac{\displaystyle\int_0^{\varOmega Nt/2} \omega[\hat \psi_1] \,{\rm d}\omega}{ \displaystyle\int_0^{\varOmega Nt/2} [\hat \psi_1] \,{\rm d}\omega} - \frac{\displaystyle\int_0^{\varOmega Nt/2} \omega [\hat \phi_1]\,{\rm d}\omega}{\displaystyle\int_0^{\varOmega Nt/2} [\hat \phi_1]\,{\rm d}\omega} \approx 10.2 \varOmega. \end{equation}

We then define the average temporal location of the modes as

(4.4)\begin{equation} t^{resp.}_{avg} := \frac{\displaystyle\int_0^{\rm T} t [\breve \psi_1] \,{\rm d} t}{\displaystyle\int_0^{\rm T} [\breve{\psi}_1] \,{\rm d} t}, \end{equation}

and

(4.5)\begin{equation} t^{forc.}_{avg} := \frac{\displaystyle\int_0^{\rm T} t [\breve \phi_1] \,{\rm d} t}{\displaystyle \int_0^{\rm T} \int_0^1 [\breve{\phi}_1] \,{\rm d} t}. \end{equation}

Assuming that both the optimal forcing and response prefer the same natural frequency $\omega _0$ with a corresponding streamwise wave speed $c_0 = \omega _0/k_1$, we estimate the shift with

(4.6)\begin{equation} (|U(x^{resp.}_{2, {avg}}, t^{resp.}_{avg})| - |U(x^{forc.}_{2, {avg}}, t^{forc.}_{avg}|)) k_1/\varOmega \approx 10.5, \end{equation}

which roughly matches the observed shift. We repeat this analysis for two other sets of spatial parameters, $(k_1, k_3) = (0.67, 2.67)$ and $(k_1, k_3) = (1.33, 2.67)$, though the plots are not shown. For the first set, the estimated frequency shift is found to be $\Delta \omega / \varOmega \approx 27$ using (4.6), roughly matching the measured mean frequency shift of $\Delta \omega / \varOmega \approx 26$ computed using a Fourier transform of the response and forcing modes. Similarly, for the second set of these length scales, the estimated frequency shift due to convection by the mean flow is $\Delta \omega / \varOmega \approx 14.6$ and the measured mean frequency shift is $\Delta \omega / \varOmega \approx 13.4$.

Wavelet-based resolvent analysis for this non-stationary flow thus reveals how a time-varying mean profile affects the linear amplification of perturbations. The mean velocity profile not only determines the spatial structure of the modes like in § 3.1, but also their transient behaviour, and in this case, acts as a convection velocity that modulates their frequency content and wave speeds.

4.2. Channel flow with sudden lateral pressure gradient

Finally, we study a fully developed turbulent channel flow at ${{Re}}_\tau = 186$ that is subjected to a sudden lateral pressure gradient ${\rm d} P/{\rm d}\kern0.7pt x_3 = \varPi \,{\rm d} P/{\rm d}\kern0.7pt x_1$ at $t=0$ with $\varPi = 30$ (Moin et al. Reference Moin, Shih, Driver and Mansour1990; Lozano-Durán et al. Reference Lozano-Durán, Giometto, Park and Moin2021). This flow, commonly referred to as a three-dimensional (3-D) channel flow, has an initial transient period dominated by 3-D non-equilibrium effects. Eventually, the flow will reach a new statistically steady state with the mean flow in the $({\rm d} P/{\rm d}\kern0.7pt x_1, {\rm d} P/{\rm d}\kern0.7pt x_3)$ direction parallel to the wall. In the transient period, the tangential Reynolds stress initially decreases before increasing linearly, with depletion and increase rate that scales as $\varPi x_2/\delta$ (Lozano-Durán et al. Reference Lozano-Durán, Giometto, Park and Moin2021).

The mean flow profiles are obtained from Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2021) and have non-zero streamwise and spanwise components $U_1$ and $U_3$ (figure 11) as well as non-zero wall-normal gradients of streamwise and spanwise components $\boldsymbol{\mathsf{dU}}_{1, 2}$ and $\boldsymbol{\mathsf{dU}}_{3, 2}$. In this section, we non-dimensionalise velocity by initial friction velocity $u^*_{\tau,0}$, lengths by the channel half-height $\delta ^*$ and time by $\delta ^* / u^*_{\tau,0}$. The Reynolds number for this problem is ${{Re}} = {{Re}}_{\tau, 0} := u^*_{\tau, 0} \delta ^* / \nu$. The time domain of the simulation is $T = 2.34$. To construct the discrete resolvent operator, we use a Chebyshev grid of size $N_2 = 65$ in the $x_2$-direction extending from $x_2 = 0$ to $x_2 = 1$. For the spatial derivatives in the $x_2$-direction, we choose second-order-accurate finite difference matrices. We enforce a no-slip and no-penetration boundary condition at the wall and a free-slip and no-penetration condition at the centreline. The boundary condition for the temporal finite difference operator $\boldsymbol{\mathsf{D}}_t$ is chosen to enforce a Neumann-type condition, $\partial _t({\cdot })|_{t=0} = \partial _t({\cdot })|_{t=T} = 0$. To reduce the impact of the boundary condition on the modes at $t=0$ we extend $U_1$ and $U_3$ to the time interval $t \in [-0.58, 2.34]$ and assume $U_i(t \leq 0, y) = U_i(t = 0, y)$ and $dP/dx_3(t < 0) = 0$. When the modes are plotted, we only show the original time domain $t \in [0, 2.34]$ and exclude the contribution from negative times. We use a temporal resolution of $N_t = 1000$ for the extended time frame. In this case, we note that we do not obtain spurious modes due to the distortion of high frequency waves, and that filtering-out those waves as in § 4.1 has little effect on the results.

Figure 11. Mean (a) streamwise and (b) spanwise velocity profile from $t = 0$ (blue) to $t = 2.34$ (red). The times shown are $t = 0, 0.58, 1.17, 1.76, 2.34$. Data taken from Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2021). Time $t$ is non-dimensionalised with $u_{\tau, 0}/\delta$.

Regarding the spatial scales for the homogeneous directions, we choose them to capture near-wall streaks at three different times: $t = 0$, $1.3$ and $1.94$. We thus tune them to represent the aspect ratio characteristic of near wall streaks, i.e. $\lambda _1^+ \approx 10 \lambda _3^+$ for a mean flow with a dominant streamwise component, and $\lambda _3^+ \approx 10 \lambda _1^+$ for a mean flow with a dominant spanwise component. Here, $({\cdot })^+$ indicates the wall scaling with ${{Re}}_{\tau,0}$, before the lateral pressure gradient is applied. To capture near-wall streaks at $t = 0$, we choose $(\lambda _1^+, \lambda _3^+) = (\lambda _{1,0}^+, \lambda _{3,0}^+) := (1000, 100)$ as in § 3.1, which corresponds to the spatial scales preferred by the near-wall streaks at ${{Re}}_\tau = 186$ prior to the lateral pressure gradient. Under the shear conditions at $t = 1.3, 1.94$, we must take into account the stronger mean shear in the spanwise direction (Lozano-Durán et al. Reference Lozano-Durán, Giometto, Park and Moin2021) by multiplying by a factor of ${{Re}}_{\tau }(t) / {{Re}}_{\tau, 0} = u^*_\tau (t) / u^*_{\tau, 0} = u_\tau$, plotted in figure 12(a). We also take into account the new orientation of the streaks by applying a rotation by the wall-shear stress angle $\gamma (t) = \tan ^{-1}(\tau _3/\tau _1)$, where $\tau _i$ is the instantaneous wall-shear stress in the $x_i$ direction (see figure 12b). Using the expressions $\lambda ^+_1 \approx (\lambda ^+_{1, 0} \cos (\gamma (t)) - \lambda ^+_{3,0} \sin (\gamma (t))) u_\tau (t)$ and $\lambda ^+_3 \approx (\lambda ^+_{1, 0} \sin (\gamma (t)) + \lambda ^+_{3,0} \cos (\gamma (t))) u_\tau (t)$, we obtain spatial parameters $(\lambda _1^+, \lambda _3^+) = (264, 1827)$ corresponding to $t = 1.3$, and $(\lambda _1^+, \lambda _3^+) = (329, 2898)$ corresponding to $t = 1.94$.

Figure 12. (a) Friction velocity $u_{\tau }$ and (b) wall shear-stress angle $\gamma = \tan ^{-1}(\tau _3 / \tau _1)$ as a function of time. Data taken from Lozano-Durán et al. (Reference Lozano-Durán, Giometto, Park and Moin2021). The vertical dashed lines are at $t = 0, 1.38, 1.94$, and correspond to the choice of $\lambda _1^+$ and $\lambda _3^+$ for the modes plotted in figure 13.

The resolvent modes for the first wavelength pair $(\lambda _1^+,\lambda _3^+) = (1000, 100)$ are shown in figure 13(a,b). The magnitude of the modes in frequency-time space is also plotted in figure 14(a,b). The resolvent modes are temporally centred around $t = 0$ and exhibit a predominant streamwise component. The modes are located in a region $x_2 < 0.25$, which corresponds to $x_2^+ < 45$, i.e. the buffer region. Thus, at $t = 0$, the modes capture the highly energetic near-wall streaks. The subsequent temporal decay of these modes can be explained by the changing flow conditions, notably the growth of the spanwise wall-shear stress $\tau _3$, and consequently $u_\tau$ (see figure 12). Under these conditions, the spatial scales preferred by the near-wall streaks stretch as $u_\tau$ increases and the wall-shear stress tensor rotates towards the $x_3$ direction.

Figure 13. Real part of (a,c,e) $\breve \psi _1$ and (b,d,f) $\breve \psi _3$ for the turbulent channel subject to a spanwise pressure gradient. The chosen spatial scales are (a,b) $\lambda _1^+ = 1000$, $\lambda _3^+ = 100$, (c,d) $\lambda _1^+ = 264$, $\lambda _3^+ = 1827$ and (e,f) $\lambda _1^+ = 329$, $\lambda _3^+ = 2898$. The vertical dashed lines mark (a,b) $t = 0$, (c,d) $t=1.38$ and (e,f) $t=1.94$.

Figure 14. Principal (a,c,e) forcing (b,d,f) response modes in the frequency-time plane for the turbulent channel subject to a spanwise pressure gradient, for (a,b) $\lambda _1^+ = 1000$, $\lambda _3^+ = 100$, (c,d) $\lambda _1^+ = 264$, $\lambda _3^+ = 1827$ and (e,f) $\lambda _1^+ = 329$, $\lambda _3^+ = 2898$.

The response mode for the second pair of spatial scales, $(\lambda _1^+,\lambda _3^+) = (264, 1827)$, tuned to conditions at $t = 1.38$, are plotted in figure 13(c,d). The frequency-time map of the modes is shown in figure 14(c,d). Similar to the first case, the modes are centred around $t = 1.38$, indicating that the wavelet-based resolvent analysis is able to identify the non-equilibrium effects of the non-stationary flow. We note that the spanwise component of the response mode is much more dominant than the streamwise component, which reflects the new wall-shear angle $\gamma = 75.7^\circ$. Finally, for the third case, $(\lambda _1^+,\lambda _3^+) = (329, 2898)$, tuned to conditions at $t = 1.94$, we observe that the modes (figures 13e,f and 14e,f) are not centred around the target time. We speculate that this is due to the temporal boundary condition at $t = 2.34$. As the flow is not at a statistically steady state at this time, a Neumann boundary condition may not be the most suitable boundary condition. The modes cannot grow beyond the boundary due to the boundary condition and are artificially damped near the end of the temporal domain.

Figure 15 shows the principal response mode in the physical domain for the three target length scale pairs, in the reference frame rotated by $\gamma$ about the $x_2$-axis. We denote the streamwise and spanwise directions in the rotated reference frame by $y_1$ and $y_3$, respectively, with $y_1 = x_1 \cos (\gamma ) - x_3 \sin (\gamma )$ and $y_3 = x_1 \sin (\gamma ) + x_3 \cos (\gamma )$ The modes resemble each other qualitatively, and capture elongated near-wall streaks in the direction of the rotated flow. The length and spanwise spacing of the streaks increases with wall shear stress, as expected. Moreover, the response modes for $(\lambda _1^+,\lambda _3^+) = (264, 1827)$ and $(\lambda _1^+,\lambda _3^+) = (329, 2898)$ (figure 15b) are concentrated closer to the wall than for $(\lambda _1^+,\lambda _3^+) = (1000, 100)$ (figure 15a), which indicates that the region of high-sensitivity to forcing moves closer to the wall as $Re_\tau$ increases. This is in line with the behaviour of near-wall turbulence: for higher $Re_\tau$, the buffer and logarithmic layers, which contain the bulk of turbulent energy in channel flow, are closer to the wall. We note, however, that the mode corresponding to $(\lambda _1^+,\lambda _3^+) = (329, 2898)$ is expected to be located at an even lower wall-normal height than the mode corresponding to $(\lambda _1^+,\lambda _3^+) = (264, 1827)$, and attribute its higher location to the effect of the temporal boundary condition at the end of the temporal domain.

Figure 15. Principal response mode in the reference frame rotated about $x_2$ by an angle $\gamma$, at (a) $t = 0$, (b) $t=1.38$ and (c) $t =1.94$. The contours represent the streamwise component; the arrows represent the velocity field in the wall-parallel plane.

We can compute flow and shear angles from the principal resolvent mode and compare them to the wall shear stress angle $\gamma$ extracted from DNS. We define the mode shear angle as

(4.7)\begin{equation} \breve \gamma(x_2, t) = \left|\tan^{{-}1} \left(\frac{{\rm d} \psi_3/ {\rm d}\kern0.7pt x_2}{{\rm d} \psi_1/{\rm d}\kern0.7pt x_2}\right)\right|_{x_1 = x_{1, {max}}, x_3 = x_{3, {max}}}, \end{equation}

where $\psi _i$ denotes the inverse Fourier and wavelet transform of the $i$th mode component, and $(x_{1, {max}}, x_{3, {max}})$ denotes the location in the $x_1$–$x_3$ plane where the amplitude of the response mode $\boldsymbol {\psi }$ is maximal at each $x_2$. We plot the results for the three spatial parameters at their respective target times (figure 16). For the mode corresponding to $(\lambda _1^+, \lambda _3^+) = (1000, 100)$, the mode shear angle $\breve \gamma$ matches the wall shear angle of zero from DNS at a wall-normal location slightly farther from the wall than the location of peak amplitude. For $(\lambda _1^+, \lambda _3^+) = (264, 1827)$ and $(\lambda _1^+, \lambda _3^+) = (329, 2898)$, we observe that $\breve \gamma$ matches $\gamma$ well at the amplitude peaks. For all three sets of spatial parameters, the mode angles at the wall itself differ significantly from the wall shear angle from DNS, though this is expected since the modes obtained are lifted from the wall and decay to match the wall boundary condition.

Figure 16. Shear angle of the principal response mode for (a) $(\lambda _1^+, \lambda _3^+) = (1000, 100)$ at $t = 0$, (b) $(\lambda _1^+, \lambda _3^+) = (264, 1827)$ at $t = 1.38$, and (c) $(\lambda _1^+, \lambda _3^+) = (329, 2898)$ at $t = 1.94$. The red vertical line (- -, red) indicates the wall shear-stress angle at the chosen times; the horizontal black line indicates the location of the peak of kinetic energy $|\breve \psi _1|^2 + |\breve \psi _2|^2 + |\breve \psi _3|^2$ at the target times, and the dashed horizontal lines indicate where the mode energy is at $10\,\%$ of the peak.

The amplification of the response mode given by the leading singular value of the resolvent operator differs across the chosen spatial parameters. The energy of the principal response mode corresponding to $(\lambda _1^+, \lambda _3^+) = (1000, 100)$, which peaks for $t < 0$, is amplified by a factor of $\sigma _1^2 \approx 8$, while the energy of the modes corresponding to $(\lambda _1^+, \lambda _3^+) = (264, 1827)$ and $(\lambda _1^+, \lambda _3^+) = (329, 2898)$, respectively, peaking at $t = 1.34$ and $t = 1.94$, are, respectively, amplified by a factor of $\sigma _1^2 \approx 0.4$ and $\sigma _1^2 \approx 0.6$, indicating an effective energy suppression. A decrease in Reynolds stresses was similarly observed in the fully turbulent system simulated during the development of the spanwise mean flow (Lozano-Durán et al. Reference Lozano-Durán, Constantinou, Nikolaidis and Karp2021), and it was proposed that the developing spanwise mean flow generates smaller transverse structures close to the wall which disrupt the coherence of the dominant streamwise rolls and inhibit them from vertically transporting momentum upward from the near-wall region.

The behaviour of the wavelet-based resolvent modes allows us to expand this proposed explanation. The optimal forcing for the mean flow conditions at a time $t$ is in the form of rolls approximately pointing in the direction $\gamma (t)$, and its corresponding response mode will eventually align itself with the forcing angle after a period of transient growth as shown in figure 17(a,b). Due to the non-normality of the linearised system, the optimal forcing peaks before the target time $t$, and thus the optimal forcing may be instantaneously misaligned with the rotating mean flow during the transient growth period of the response. Indeed, the forcing mode is not computed to maximise the instantaneous kinetic energy amplification, but maximises instead the time-integrated kinetic energy of the response over the entire temporal domain.

Figure 17. Flow angle of the principal forcing (red) and response (black) modes for (a) $(\lambda _1^+, \lambda _3^+) = (1000, 100)$, (b) $(\lambda _1^+, \lambda _3^+) = (264, 1827)$ and (c) $(\lambda _1^+, \lambda _3^+) = (329, 2898)$. The vertical lines indicates the amplitude peaks for the forcing (red) and response (black) modes, and the dotted black line represents $\gamma (t)$.

More insight can be gained by considering the linearised equation governing the wall-normal vorticity $v_2$,

(4.8)\begin{align} \left(\partial_t + U_1\partial_{x_1} + U_3 \partial_{x_3} - \frac{1}{{{Re}}}(\partial_{x_1}^2 + \partial_{x_2}^2 + \partial_{x_3}^2)\right) v_2 = \left(-\frac{{\rm d} U_1}{{\rm d}\kern0.7pt x_2} \partial_{x_3} + \frac{{\rm d} U_3}{{\rm d}\kern0.7pt x_2} \partial_{x_1}\right) u_2 + g_2, \end{align}

where $g_2 := \partial _{x_3} f_1 - \partial _{x_1} f_3$ denotes the external forcing. We note that in this problem, the mean wall-normal flow $U_2$, and the wall-normal gradients ${\rm d} U_1/ {\rm d}\kern0.7pt x_2$ and ${\rm d} U_3/ {\rm d}\kern0.7pt x_2$ are zero. Introducing a flow angle $\mu := \tan ^{-1}(U_3/U_1)$ and a shear angle $\upsilon := \tan ^{-1} (({\rm d} U_3 /{\rm d}\kern0.7pt x_2)/ ({\rm d} U_1 /{\rm d}\kern0.7pt x_2) )$, we can rewrite (4.8) in the reference frame locally rotated by $\mu$ as

(4.9)\begin{align} & \left(\partial_t + (U_1^2 + U_3^2)^{{1}/{2}} \partial_{x'_1} - \frac{1}{{{Re}}}(\partial_{x'_1}^2 + \partial_{x_2}^2 + \partial_{x'_3}^2 ) \right) v_2 \nonumber\\ &\quad =\left( \frac{{\rm d} U_1^2}{{\rm d}\kern0.7pt x_2} + \frac{{\rm d} U_3^2}{{\rm d}\kern0.7pt x_2} \right)^{{1}/{2}} \left(- \cos(\mu-\upsilon) \partial_{x_3'} + \sin(\mu -\upsilon)\partial_{x_1'}\right) u_2 + g_2, \end{align}

where $(x_1', x_2, x_3')$ denote the coordinates rotated anticlockwise by angle $\mu$ in the $x_1$–$x_3$ plane. The left-hand side of (4.9) is identical to the classical Squire equation for perturbations about a one-dimensional streamwise mean flow, with an effective advection velocity of $(U_1^2 + U_3^2)^{1/2}$. The right-hand side includes a lift-up term modified by the misalignment between the mean shear and velocity profiles. If the mean flow is at equilibrium and $\mu = \upsilon$, we obtain an identical lift-up term to the classical Squire equation, with an effective mean shear of $( ({\rm d} U_1/{\rm d}\kern0.7pt x_2)^2 + ({\rm d} U_3/{\rm d}\kern0.7pt x_2)^2 )^{1/2}$. The lift-up term in the Squire equation constitutes a way for the wall-normal velocity perturbations to force the streamwise and spanwise components, as discussed in § 3.2 and Jiménez (Reference Jiménez2013, Reference Jiménez2018). Through this coupling, the streamwise and spanwise velocity components can be efficiently forced by $f_2$ along with $f_1$ and $f_3$. However, in a non-equilibrium rotating flow, the mean shear profile lags behind the mean velocity, i.e. $\mu \neq \upsilon$, causing a misalignment between the two profiles which is particularly pronounced for regions farther away from the wall. The term pertaining to the misalignment in (4.9) damps the effect of the wall-normal velocity perturbation $u_2$ on the wall-normal vorticity, and consequently, the streamwise and spanwise velocity perturbations. This is illustrated by figure 18: for the initial one-dimensional mean flow configuration in figure 18(a), the optimal forcing mode exploits the lift-up mechanism via a strong wall-normal component; for the rotated configurations in figure 18(b,c), the wall-normal component of the optimal forcing is significantly attenuated.

Figure 18. Amplitudes for the streamwise (blue $\bullet$), spanwise (red $\blacktriangle$) and wall-normal (green $\blacksquare$) components of the principal forcing (dotted line, hollow markers) and response (solid line, filled markers) modes, for (a) $(\lambda _1^+, \lambda _3^+) = (1000, 100)$, (b) $(\lambda _1^+, \lambda _3^+) = (264, 1827)$ and (c) $(\lambda _1^+, \lambda _3^+) = (329, 2898)$.

Thus, the lag of the mean shear profile disrupts the coupling between the wall-normal, and the streamwise and spanwise velocity perturbations. In a non-equilibrium rotating mean flow, wall-normal velocity perturbations are less efficient at producing slow-decaying streaks in the flow direction (than in § 3.2, for example), and the mean flow selects rolls with a much weaker wall-normal component.