2.1. Nonlinear forcing in the NSE

The nonlinear (quadratic) terms $\boldsymbol{f}$ are defined in physical space in (2.3) through a pointwise multiplication, while in (discrete) Fourier space the nonlinear forcing at a wavenumber–frequency triplet $\boldsymbol{k}_3$ can be written in terms of a convolution of the velocity fields and velocity gradients at $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ :

(2.6) \begin{align} \boldsymbol{f}(\boldsymbol{k}_3,y) = - \sum _{ \boldsymbol{k}_{1}+\boldsymbol{k}_{2} = \boldsymbol{k}_3} \boldsymbol{u}(\boldsymbol{k}_1,y) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(\boldsymbol{k}_2,y). \end{align}

The requirement of $\boldsymbol{k}_1 + \boldsymbol{k}_2 = \boldsymbol{k}_3$ is the triadic compatibility or resonance constraint, a result of the quadratic nature of the nonlinearity.

In the resolvent formulation, the Fourier-transformed NSE are written in an input–output form, where $\boldsymbol{f}(\boldsymbol{k},y)$ is considered an input forcing to the resolvent operator $\mathcal{H}(\boldsymbol{k},y)$ as shown in (2.5). It can be seen that the linear operator does not modify the scale of the input, such that each wavenumber–frequency triplet operates independently from all others. On the other hand, (2.6) shows that the nonlinear forcing is responsible for the coupling of different scales, and therefore the distribution of energy among scales.

The triadic interactions are visually depicted in figure 1. First, the velocity fields at $\boldsymbol{k}_1$ nonlinearly interact with the velocity gradients at $\boldsymbol{k}_2$ , generating part of the forcing at $\boldsymbol{k}_3= \boldsymbol{k}_1 + \boldsymbol{k}_2$ . This triadic contribution to the forcing is studied using the forcing coefficients $P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ defined later in § 2.5. The forcing at $\boldsymbol{k}_3$ is then passed through the linear resolvent operator to generate the velocity response at the same wavenumber–frequency triplet $\boldsymbol{k}_3$ . This triadic contribution to the velocity response, which involves both the nonlinear convolution and the linear resolvent operator, is studied using the response coefficients $R(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ defined later in § 2.6.

Figure 1.

Diagram for the triadic interactions. The velocity and velocity gradient at $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ interact nonlinearly to generate part of the forcing at $\boldsymbol{k}_3 = \boldsymbol{k}_1 + \boldsymbol{k}_2$ . The full forcing is a convolution sum of all pairs of $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ that are triadically compatible with $\boldsymbol{k}_3$ , which forces the resolvent operator to generate the response. The red inner box contains the triadic contributions to the nonlinear forcing, studied using the coefficients $P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ , and the black outer box contains the triadic contributions to the response, studied using the coefficients $R(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ .

2.2. From time to frequency

To perform the spatiotemporal analyses, the time domain data needs to be transformed into the Fourier domain. For this purpose, the Welch method is applied where the temporal snapshots are segmented, with the Hann window function applied to each temporal segment before taking the temporal Fourier transform. Finally, the results are averaged across all temporal segments. To analyse the effect of the window function, similar to Nogueira et al. (Reference Nogueira, Morra, Martini, Cavalieri and Henningson2021) and Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021), the momentum equation is multiplied by the window function $w(t)$ and the spatial dimensions are temporarily ignored as they do not play a role in this analysis of the temporal window function,

(2.7) \begin{align} w(t) \frac {\partial }{\partial t} \boldsymbol{{u}}(t) + w(t)\mathcal{L}\boldsymbol{{u}}(t) = w(t) \boldsymbol{{f}}(t). \end{align}

In this equation, $\mathcal{L}$ is the linear part of the Navier–Stokes operator, and the equation can then be rewritten as

(2.8) \begin{align} \frac {\partial }{\partial t} \left [ w(t)\boldsymbol{{u}}(t) \right ] + \mathcal{L} \left [ w(t)\boldsymbol{{u}}(t) \right ] - \boldsymbol{{u}}(t)\frac{\text{d}}{\text{d}t}w(t) = [ w(t)\boldsymbol{{f}}(t) ]. \end{align}

Performing the Fourier transform in time (denoted by $\mathcal{F} [ \boldsymbol{\cdot }]$ ) on the signals with the window function applied, $\boldsymbol{{u}}(\omega ) = \mathcal{F} [ w(t)\boldsymbol{{u}}(t) ], \boldsymbol{f}(\omega ) = \mathcal{F} [ w(t)\boldsymbol{{f}}(t) ]$ , and defining the spurious forcing, $\boldsymbol{{s}}(\omega ) = \mathcal{F} [ \boldsymbol{{u}}(t) ({\text{d}}/{\text{d}t})w(t) ]$ , the above equation can be rewritten using the resolvent operator $\mathcal{H}$ ,

(2.9) \begin{align} \boldsymbol{{u}}(\omega ) = \mathcal{H} \left ( \boldsymbol{f}(\omega ) +\boldsymbol{{s}}(\omega ) \right ), \end{align}

which can be reorganized as

(2.10) \begin{align} \boldsymbol{\tilde {u}}(\omega ) = \boldsymbol{{u}}(\omega ) - \mathcal{H}\boldsymbol{{s}}(\omega ) = \mathcal{H} \boldsymbol{f}(\omega ). \end{align}

Here $\boldsymbol{\tilde {u}}(\omega )$ are the velocity Fourier modes corrected for the spurious forcing due to the window function, which satisfies the input–output form of the resolvent analysis introduced in (2.5).

Additionally, the forcing Fourier modes, $\boldsymbol{f}(\omega ) = \mathcal{F} [ w(t)\boldsymbol{{f}}(t) ]$ , can be computed in two equivalent ways,

(2.11) \begin{align} \boldsymbol{f}(\omega ) = \mathcal{F} \left [ w(t) \boldsymbol{{u}}(t) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{{u}}(t) \right ] = \mathcal{F} \big [ \sqrt {w(t)} \boldsymbol{{u}}(t) \big ] * \mathcal{F} \big [ \sqrt {w(t)} \boldsymbol{\nabla }\boldsymbol{{u}}(t) \big ], \end{align}

where $*$ denotes the convolution operator. The first method computes $\boldsymbol{{f}}(t)$ in physical time, then applies the window function and takes the Fourier transform, while the second method computes $\boldsymbol{f}(\omega )$ in the Fourier space through a convolution. Note that to ensure the equivalence between the two methods, the second method requires the use of the $\sqrt {w(t)}$ as the window function applied to the velocity and velocity gradients. For this work, the periodic Hann window function is utilized in the analyses, and the choice of the window function is not expected to impact the overall structure of the resulting coefficients.

2.3. Singular value decomposition of the resolvent operator

The discrete resolvent operator $\mathcal{H}(\boldsymbol{k},y)$ is constructed using a compact finite difference in $y$ matching the numerical scheme of the DNS described later in § 3. A singular value decomposition can then be performed on the discrete operator,

(2.12) \begin{align} \mathcal{H}(\boldsymbol{k},y) = \sum _{q} \psi _q(\boldsymbol{k},y) \sigma _q(\boldsymbol{k}) \phi _q^*(\boldsymbol{k},y), \end{align}

where $\psi _q$ are the singular response modes (henceforth referred to as resolvent modes), $\sigma _q$ are the (ordered) singular values and $\phi _q$ are the singular forcing modes. Superscript $(\boldsymbol{\cdot })^*$ denotes a complex conjugate. Velocity Fourier modes can then be expressed as

(2.13) \begin{align} \boldsymbol{{u}}(\boldsymbol{k},y) = \sum _{q} \chi _q(\boldsymbol{k}) \sigma _q(\boldsymbol{k}) \psi _q(\boldsymbol{k},y), \end{align}

where $\chi _q(\boldsymbol{k}) = \phi ^*_q(\boldsymbol{k},y) \boldsymbol{{f}}(\boldsymbol{k},y)$ are the nonlinear weights obtained by projecting the nonlinear forcing (if known) onto the singular forcing modes.

2.4. Energy transfer by the nonlinear forcing

We start by defining the spectral TKE as $e(\boldsymbol{k},y) = \left | u(\boldsymbol{k},y) \right |^2 + \left | v(\boldsymbol{k},y) \right |^2 + \left | w(\boldsymbol{k},y) \right |^2$ , which is the energy of Fourier modes at given $\boldsymbol{k}$ . An equation for the spectral TKE can be written as

(2.14) \begin{align} & \underbrace { \textit{Re} \left \{ u^*(\boldsymbol{k},y) v(\boldsymbol{k},y) \overline {U}'(y) \right \} \vphantom {\frac{\text{d}^2}{\text{d}y^2}} }_{\text{Production}} + \underbrace { \frac {k^2 }{Re} e(\boldsymbol{k},y) + \frac {1}{Re} \frac{\text{d}}{\text{d}y}u_i^*(\boldsymbol{k},y) \frac{\text{d}}{\text{d}y} u_i(\boldsymbol{k},y) }_{\text{Viscous dissipation}} \nonumber \\[0.5em] +& \underbrace { \textit{Re} \left \{ \frac{\text{d}}{\text{d}y} \left [ v^*(\boldsymbol{k},y) p(\boldsymbol{k},y) \right ]\right \} }_{\text{Pressure transport}} - \underbrace { \frac {1}{2} \frac {1}{Re} \frac{\text{d}^2}{\text{d}y^2} e(\boldsymbol{k},y) \vphantom {\textit{Re} \left \{ \frac{\text{d}}{\text{d}y}\right \} } }_{\text{Viscous transport}} = \underbrace { \textit{Re} \left \{ u_i^*(\boldsymbol{k},y) f_i(\boldsymbol{k},y) \vphantom {\frac{\text{d}}{\text{d}y}} \right \} }_{\text{Turbulent transport}}, \end{align}

where $\textit{Re} \left \{ \boldsymbol{\cdot }\right \}$ indicates the real part, and the summation notation is used with the subscripts $i$ . To derive this, the three momentum equations are first multiplied with the complex conjugate of the three velocity components, summed and then combined with the continuity equation. It is then added to the equation for $e(-\boldsymbol{k},y)$ , utilizing $e(\boldsymbol{k},y) = e(-\boldsymbol{k},y)$ and the symmetry of Fourier modes to obtain the final form. Details are given in Appendix A.1 of Huang (Reference Huang2025).

The turbulent transport term can be alternatively written as

(2.15) \begin{align} u_i^*(\boldsymbol{k},y) f_i(\boldsymbol{k},y) =& \ \boldsymbol{u}^*(\boldsymbol{k},y) \boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{k},y) \nonumber \\[0.5em] =& - \boldsymbol{u}^*(\boldsymbol{k},y) \boldsymbol{\cdot }\sum _{\boldsymbol{k}_1 + \boldsymbol{k}_2 = \boldsymbol{k} } \boldsymbol{u}(\boldsymbol{k}_1,y) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(\boldsymbol{k}_2,y), \end{align}

which shows that the turbulent transport is a term that involves the nonlinear forcing, and therefore transports energy between different wavenumber–frequency triplets $\boldsymbol{k}$ . Additionally, Schmid & Henningson (Reference Schmid and Henningson2001) and Barthel (Reference Barthel2022) have demonstrated that the nonlinear turbulent transport is energy conserving on a triad by triad basis.

Alternatively, the resolvent formulation $\boldsymbol{u}(\boldsymbol{k},y) = \mathcal{H}(\boldsymbol{k},y) \, \boldsymbol{{f}}(\boldsymbol{k},y)$ can be utilized to obtain

(2.16) \begin{align} e(\boldsymbol{k},y) =& \ \boldsymbol{u}^*(\boldsymbol{k},y) \, \mathcal{H}(\boldsymbol{k},y) \boldsymbol{{f}}(\boldsymbol{k},y) \nonumber \\[0.5em] =& -\boldsymbol{u}^*(\boldsymbol{k},y) \, \mathcal{H}(\boldsymbol{k},y) \sum _{\boldsymbol{k}_1 + \boldsymbol{k}_2 = \boldsymbol{k} } \boldsymbol{u}(\boldsymbol{k}_1,y) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(\boldsymbol{k}_2,y). \end{align}

The formulation of (2.14) and (2.16), both derived from the NSE, are mathematically equivalent, yet they have slightly different interpretations. Equation (2.14) is an energy balance equation. The energy transported into or out of $\boldsymbol{k}$ through nonlinear interactions (the turbulent transport term), is balanced by four other mechanisms, and is therefore not directly correlated to the increase or decrease of the spectral TKE at this $\boldsymbol{k}$ . Equation (2.16) on the other hand, considers the nonlinearity as a forcing that drives the turbulent perturbations and activates the linear mechanisms of production, pressure transport and viscous transport and dissipation (all contained in the linear resolvent operator $\mathcal{H}$ ), therefore providing a direct link between the nonlinearity and the spectral TKE.

2.5. Triadic contributions to the forcing

The contribution from the interaction between a pair of triplets $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ to the resulting forcing at $\boldsymbol{k}_3 = \boldsymbol{k}_1 + \boldsymbol{k}_2$ , can be quantified through a forcing coefficient $P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ ,

(2.17) \begin{align} P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}) = - \mathbb{E}\left \{ \left \langle \boldsymbol{f}(\boldsymbol{k}_1 + \boldsymbol{k}_2,y) , \ \boldsymbol{u}(\boldsymbol{k}_1,y) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(\boldsymbol{k}_2,y) \right \rangle _y \right \}, \end{align}

which is the inner product between $-\boldsymbol{u}(\boldsymbol{k}_1,y) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(\boldsymbol{k}_2,y)$ and $\boldsymbol{f}(\boldsymbol{k}_3,y)=\boldsymbol{f}(\boldsymbol{k}_1 + \boldsymbol{k}_2,y)$ . Here $\mathbb{E}\left \{\boldsymbol{\cdot }\right \}$ is the expected value operator, indicating an average over all temporal segments (or different realizations) and the weighted inner product $\left \langle a(y), \ b(y) \right \rangle _y$ is defined as the integral over a certain $y$ range in the continuous domain and approximated in the discrete space using a weight matrix $W$ with the appropriate integration coefficients on the diagonal:

(2.18) \begin{align} \left \langle a(y), \ b(y) \right \rangle _y = \int _y a^*(y) b(y) {\:\mathrm{d}} y \approx a^* W b. \end{align}

In this work, three separate wall-normal ( $y$ ) ranges will be used, loosely corresponding to the near-wall $y^+ \in (0, 30)$ , overlap $y^+ \in (30, 200)$ and wake regions $y^+ \in (200, 550)$ , to study the wall-normal variations of these coefficients.

The intentionally un-normalized coefficients take the forcing magnitude into consideration, e.g. a large fractional contribution to a small magnitude forcing is treated as unimportant. Note that these coefficients, defined for DNS Fourier modes, differ from those, e.g. in McKeon (Reference McKeon2017), which are defined for the interactions between resolvent modes.

The results span the (large) parameter space spanned by $(k_x, k_z, \omega , y)$ . To facilitate computation and visualization of $P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ , we define $P_{k_x}$ , $P_{k_z}$ and $P_{\omega }$ by summation in four of the six scale dimensions:

(2.19a) \begin{align} P_{k_x}(k_{x1}, k_{x2}) &= \sum _{k_{z1}}\sum _{k_{z2}}\sum _{\omega _1}\sum _{\omega _2} P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}), \\[-12pt] \nonumber \end{align}

(2.19b) \begin{align} P_{k_z}(k_{z1}, k_{z2}) &= \sum _{k_{x1}}\sum _{k_{x2}}\sum _{\omega _1}\sum _{\omega _2} P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}), \\[-12pt] \nonumber \end{align}

(2.19c) \begin{align} P_{\omega }(\omega _1, \omega _2) &= \sum _{k_{x1}}\sum _{k_{x2}}\sum _{k_{z1}}\sum _{k_{z2}} P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}). \\[10pt] \nonumber \end{align}

Here $P_{k_x}$ defined in (2.19a ) describes the interaction in the streamwise direction between $k_{x1}$ and $k_{x2}$ , summed over all possible $k_z$ and $\omega$ interactions. Similarly, $P_{k_z}(k_{z1}, k_{z2})$ describes the interaction between $k_{z1}$ and $k_{z2}$ , and $P_{\omega }$ describes the interaction between $\omega _1$ and $\omega _2$ .

Using the Hermitian symmetry of the velocity Fourier modes, it can be shown that the forcing coefficients are also Hermitian functions. However, they are asymmetric about their two arguments, due to the action of the velocity gradient tensor. The forcing coefficient satisfies

(2.20) \begin{align} P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}) = P^*(-\boldsymbol{k}_{1}, -\boldsymbol{k}_{2})\neq P(\boldsymbol{k}_{2}, \boldsymbol{k}_{1}), \end{align}

with $P_{k_x}$ , $P_{k_z}$ and $P_{\omega }$ satisfying the same property. We retain the two separate coefficients associated with each combination of $(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ to maximize the information about the dominant interactions within each triad, which will be lost under a combined symmetric coefficient.

Utilizing (2.6), we can obtain the summation property of the forcing coefficient,

(2.21) \begin{align} \sum _{ \boldsymbol{k}_1 } P(\boldsymbol{k}_{1}, \boldsymbol{k}_{3} - \boldsymbol{k}_{1}) = \mathbb{E}\left \{\left \langle \boldsymbol{f}(\boldsymbol{k}_3,y) , \ \boldsymbol{f}(\boldsymbol{k}_3,y) \right \rangle _y\right \}, \end{align}

where the right-hand side of the equation is the spectral energy of the forcing at $\boldsymbol{k}_3$ , a real positive quantity. As a result, the forcing coefficients can be interpreted as the triadic contributions to the forcing spectral energy. The positive and negative real parts of the coefficients represent energy injection and extraction, respectively, while the complex parts of the coefficients cancel out in an integral sense due to the Hermitian symmetry described in (2.20), providing an additional constraint on the ensemble of nonlinear interactions and resulting in a real positive sum. Similar properties can be obtained for $P_{k_x}$ , $P_{k_z}$ and $P_{\omega }$ :

(2.22a) \begin{align} \sum _{k_{x1}} P_{k_x}(k_{x1}, k_{x3}-k_{x1}) &= \sum _{k_{z3}}\sum _{\omega _3} \mathbb{E}\left \{ \left \langle \boldsymbol{f}(\boldsymbol{k}_3,y) , \ \boldsymbol{f}(\boldsymbol{k}_3,y) \right \rangle _y \right \}, \\[-12pt] \nonumber \end{align}

(2.22b) \begin{align} \sum _{k_{z1}} P_{k_z}(k_{z1}, k_{z3}-k_{z1}) &= \sum _{k_{x3}}\sum _{\omega _3} \mathbb{E}\left \{ \left \langle \boldsymbol{f}(\boldsymbol{k}_3,y) , \ \boldsymbol{f}(\boldsymbol{k}_3,y) \right \rangle _y \right \}, \\[-12pt] \nonumber \end{align}

(2.22c) \begin{align} \sum _{\omega _1} P_{\omega }(\omega _{1}, \omega _{3}-\omega _{1}) &= \sum _{k_{x3}}\sum _{k_{z3}} \mathbb{E}\left \{ \left \langle \boldsymbol{f}(\boldsymbol{k}_3,y) , \ \boldsymbol{f}(\boldsymbol{k}_3,y) \right \rangle _y \right \}. \\[10pt] \nonumber \end{align}

2.6. Triadic contributions to the velocity response

The forcing contributions at all $\boldsymbol{k}$ are not equally amplified. To study the triadic contribution from the interaction between a pair of triplets $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ to the resulting velocity response at $\boldsymbol{k}_3 = \boldsymbol{k}_1 + \boldsymbol{k}_2$ , including the effect of the linear resolvent operator, we will pass $-\boldsymbol{u}(\boldsymbol{k}_1,y) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(\boldsymbol{k}_2,y)$ through $\mathcal{H}(\boldsymbol{k}_3,y)$ and then take the inner product with the resulting velocity response $\boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y)$ to define the response coefficient,

(2.23) \begin{align} R(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}) = - \mathbb{E}\left \{ \left \langle \boldsymbol{\tilde {u}}(\boldsymbol{k}_1+\boldsymbol{k}_2,y), { \mathcal{H}(\boldsymbol{k}_1+\boldsymbol{k}_2,y) \left [ \boldsymbol{u}(\boldsymbol{k}_1,y) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(\boldsymbol{k}_2,y) \right ] } \right \rangle _y \right \}, \end{align}

where $\boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y)$ , introduced previously in (2.10), are the velocity Fourier modes with a correction to remove the effect of the window function in the temporal Fourier transform. The coefficients are again intentionally un-normalized to take the response magnitude into consideration.

Similar to the previous section, the two-dimensional coefficients $R_{k_x}$ , $R_{k_z}$ and $R_{\omega }$ are defined by summation in four of the six dimensions:

(2.24a) \begin{align} R_{k_x}(k_{x1}, k_{x2}) &= \sum _{k_{z1}}\sum _{k_{z2}}\sum _{\omega _1}\sum _{\omega _2} R(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}), \\[-12pt] \nonumber \end{align}

(2.24b) \begin{align} R_{k_z}(k_{z1}, k_{z2}) &= \sum _{k_{x1}}\sum _{k_{x2}}\sum _{\omega _1}\sum _{\omega _2} R(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}), \\[-12pt] \nonumber \end{align}

(2.24c) \begin{align} R_{\omega }(\omega _1, \omega _2) &= \sum _{k_{x1}}\sum _{k_{x2}}\sum _{k_{z1}}\sum _{k_{z2}} R(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}). \\[10pt] \nonumber \end{align}

The Hermitian symmetry and the asymmetry about the two arguments are also satisfied by the response coefficient,

(2.25) \begin{align} R(\boldsymbol{k}_{1}, \boldsymbol{k}_{2}) = R^*(-\boldsymbol{k}_{1}, -\boldsymbol{k}_{2})\neq R(\boldsymbol{k}_{2}, \boldsymbol{k}_{1}), \end{align}

with $R_{k_x}$ , $R_{k_z}$ and $R_{\omega }$ satisfying the same property.

Utilizing (2.6) and (2.10), the summation properties of the response coefficient can also be obtained,

(2.26) \begin{align} \sum _{ \boldsymbol{k}_1 } R(\boldsymbol{k}_{1}, \boldsymbol{k}_{3} - \boldsymbol{k}_{1}) = \mathbb{E}\left \{ \left \langle \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) , \ \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) \right \rangle _y \right \}, \end{align}

where the right-hand side of the equation is the spectral TKE of the velocity response at $\boldsymbol{k}_3$ , and is again a real positive quantity. As a result, the response coefficients can be interpreted as the triadic contributions to the spectral TKE. The positive and negative real parts of the coefficients represent the injection and extraction of spectral TKE at a given wavelength or frequency, while the complex parts of the coefficients cancel out due to the Hermitian symmetry described in (2.25). Similar properties can be obtained for $R_{k_x}$ , $R_{k_z}$ and $R_{\omega }$ :

(2.27a) \begin{align} \sum _{k_{x1}} R_{k_x}(k_{x1}, k_{x3}-k_{x1}) &= \sum _{k_{z3}}\sum _{\omega _3} \mathbb{E}\left \{ \left \langle \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) , \ \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) \right \rangle _y \right \}, \\[-12pt] \nonumber \end{align}

(2.27b) \begin{align} \sum _{k_{z1}} R_{k_z}(k_{z1}, k_{z3}-k_{z1}) &= \sum _{k_{x3}}\sum _{\omega _3} \mathbb{E}\left \{ \left \langle \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) , \ \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) \right \rangle _y \right \}, \\[-12pt] \nonumber \end{align}

(2.27c) \begin{align} \sum _{\omega _1} R_{\omega }(\omega _{1}, \omega _{3}-\omega _{1}) &= \sum _{k_{x3}}\sum _{k_{z3}} \mathbb{E}\left \{ \left \langle \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) , \ \boldsymbol{\tilde {u}}(\boldsymbol{k}_3,y) \right \rangle _y \right \}. \\[10pt] \nonumber \end{align}

The calculation of these coefficients is computationally intensive in the wavenumber–frequency domain. Thus, for efficiency, they are performed in the physical domain for all variables except the one parameterizing the coefficient. For example, $P_{k_x}$ is calculated in $(z,t)$ ,

(2.28) \begin{align} & P_{k_x}(k_{x1}, k_{x2}) \nonumber \\ &= -\frac {1}{{N_z N_t}} \sum _{z} \sum _{t} \mathbb{E}\left \{ \left \langle \boldsymbol{f}(k_{x1}+k_{x2},y,z,t), { \boldsymbol{u}(k_{x1},y,z,t) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}(k_{x2},y,z,t) } \right \rangle _y \right \}. \end{align}

This form provides the alternative physical interpretations of the coefficients as quantifying the importance of interactions between $k_{x1}$ and $k_{x2}$ averaged over all spanwise locations and time instead of a summation over wavenumbers and frequencies in the Fourier domain.

A similar procedure can be followed for $P_{k_z}$ , $P_{\omega }$ , $R_{k_x}$ , $R_{k_z}$ and $R_\omega$ , ensuring that the full six-dimensional $P(\boldsymbol{k}_{1}, \boldsymbol{k}_{2})$ is never computed, reducing computation cost, memory and storage requirements. Further details may be found in Huang (Reference Huang2025).

2.7. Relation to triple correlation, bispectrum and trispectrum

The forcing and response coefficients studied here may be related to the more well-known triple correlation and bispectrum. Following Lii, Rosenblatt & Van Atta (Reference Lii, Rosenblatt and Van Atta1976), the three-point spatial triple correlation for three state variables $q_l$ , $q_m$ and $q_n$ can be defined as

(2.29) \begin{align} R_{\textit{lmn}}(\boldsymbol{r},\boldsymbol{r'}) = \langle q_l(\boldsymbol{x})q_m(\boldsymbol{x}+\boldsymbol{r})q_n(\boldsymbol{x}+\boldsymbol{r'}) \rangle _{\boldsymbol{x}}, \end{align}

where $\langle \boldsymbol{\cdot }\rangle _{\boldsymbol{x}}$ represents a spatial average. Two triple spatial Fourier transforms of $R_{\textit{lmn}}(\boldsymbol{r},\boldsymbol{r'})$ in $\boldsymbol{r}, \boldsymbol{r'}$ lead to the three-dimensional spatial bispectrum:

(2.30) \begin{align} B_{\textit{lmn}}(\boldsymbol{\hat {k}}, \boldsymbol{\hat {k}'}) = q_{l}^{*}(\boldsymbol{\hat {k}}+\boldsymbol{\hat {k}'})q_m(\boldsymbol{\hat {k}})q_n(\boldsymbol{\hat {k}'}). \end{align}

The forcing coefficients proposed in this work in (2.17) can be seen as a spatiotemporal extension of the bispectrum, considering the spatiotemporal wavenumbers, $\boldsymbol{k}$ , that is suitable for the analysis of wall-bounded flows with the three homogeneous coordinates of $x$ , $z$ and $t$ . An additional difference is that the forcing coefficients in this work are specifically designed to study the contribution of an interacting pair on the resulting forcing, which to our knowledge has not been studied before. The three terms contributing to the interaction coefficients are the velocity $\boldsymbol{u}$ , velocity gradient $\boldsymbol{\nabla }\boldsymbol{u}$ and forcing $\boldsymbol{f}$ in contrast to previous studies that generally focused on the three terms being the same flow quantity: for example, $\partial u/\partial x$ in Lii et al. (Reference Lii, Rosenblatt and Van Atta1976) to study the nonlinear energy transfer between scales and velocity $\boldsymbol{u}$ in Schmidt (Reference Schmidt2020), where the BMD is introduced to compute modes associated with triadic interactions through maximization of the integral bispectral density.

Both the BMD (Schmidt Reference Schmidt2020) and the TOD (Yeung et al. Reference Yeung, Chu and Schmidt2024) consider third-order statistics. The BMD focuses on the two-point autobispectral density matrix: $\boldsymbol{B}(\boldsymbol{x}, \boldsymbol{x'}, \omega _1, \omega _2) = \boldsymbol{u}(\boldsymbol{x},\omega _1) \boldsymbol{u}(\boldsymbol{x},\omega _2) \boldsymbol{u}^*(\boldsymbol{x}',\omega _1 + \omega _2)$ , while TOD focuses on the two-point cross-bispectral covariance tensor, $\boldsymbol{S}(\boldsymbol{x}, \boldsymbol{x'}, \omega _1, \omega _2) = [ \boldsymbol{u}(\boldsymbol{x},\omega _1) \boldsymbol{\cdot }\boldsymbol{\nabla }] \boldsymbol{u}(\boldsymbol{x},\omega _2) \boldsymbol{u}^*(\boldsymbol{x}',\omega _1 + \omega _2)$ . Both methods have implied spatial cross-correlation and aim to extract optimal modes with respect to certain metrics. However, without consideration of the linear operator, both BMD and TOD are directly related to the turbulent transport term in (2.14), and characterize the transfer of energy between different frequencies. This work, on the other hand, does not intend to extract optimal modes, but aims to further the understanding of the combined effect of the linear and nonlinear mechanisms. The response coefficients proposed in (2.23) can be considered as a weighted version of the spatiotemporal extension to the bispectrum. This spatiotemporal extension considers spatial wavenumbers, which allows for the study of nonlinear triadic interactions between structures of different spatial sizes. More importantly, the response coefficients include the effect of the linear resolvent operator, acting as the weight matrix. This inclusion of the linear resolvent operator takes into account the effects of other terms in the spectral TKE (2.14), such as viscosity, pressure and production by the mean shear. For example, energy transfer from one pair of $(\boldsymbol{k}_1,\boldsymbol{k}_2)$ to $\boldsymbol{k}_3 = \boldsymbol{k}_1 + \boldsymbol{k}_2$ with energy distributed mainly at a wall-normal location with strong mean shear could potentially increase the mode amplitude at this location and induce a stronger production, resulting in some energy amplification, therefore contributing more significantly to the spectral TKE at $\boldsymbol{k}_3$ , while another pair transferring energy to the same $\boldsymbol{k}_3$ but with energy distributed mainly at another $y$ location with weak mean shear will not have this effect. The inclusion of the linear resolvent operator as a weight matrix accounts for these mechanisms, and therefore provides a new and more complete picture of the effects of the triadic energy transfer on the resulting spectral TKE. Additionally, the nonlinear forcing can be decomposed into irrotational and solenoidal parts using the Helmholtz decomposition (Rosenberg Reference Rosenberg2018; Morra et al. Reference Morra, Nogueira, Cavalieri and Henningson2021). The irrotational part of the forcing has no effect on the resulting velocity fields, and is naturally eliminated when multiplied by the resolvent operator. Therefore, the response coefficient defined in (2.23) naturally removes the effect of the inactive irrotational forcing, which differs from both the forcing coefficients defined in (2.17) and the bispectrum analyses.

As a point of clarification, in the space–time Fourier domain, amplitude modulation appears in a manner similar to ‘beating’, i.e. by the linear superposition of triadically consistent (resonant) Fourier modes, as discussed in detail in, e.g. Duvvuri & McKeon (Reference Duvvuri and McKeon2015). It may be considered a combination of linear and nonlinear effects. Consider a large scale at $\boldsymbol{k}_1$ and a small scale at $\boldsymbol{k}_2$ . Here $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ nonlinearly interact to generate a mode at $\boldsymbol{k}_1 + \boldsymbol{k}_2$ ; $\boldsymbol{k}_1$ and $-\boldsymbol{k}_2$ nonlinearly interact to generate a mode at $\boldsymbol{k}_1 - \boldsymbol{k}_2$ . When linearly superposing the two modes at $\boldsymbol{k}_1 + \boldsymbol{k}_2$ and $\boldsymbol{k}_1 - \boldsymbol{k}_2$ , we can obtain a ‘beating’ mode with frequency $\boldsymbol{k}_2$ and amplitude modulation at frequency $\boldsymbol{k}_1$ , which is the amplitude modulation mechanism. A different operation on resonant triads leads to the energy transfer discussed here, where we focused on the nonlinear interactions and not the linear superposition. In particular, the amplitude modulation is related to one component of the nonlinear forcing (corresponding to the small-scale stress, typically in $u$ but without considering the gradient operation) for low $k_x$ . Therefore, amplitude modulation has a more direct relationship to the skewness and bispectrum. With the inclusion of the gradient operator and linear resolvent operator in this study, the relationship to the amplitude modulation coefficient becomes less direct.

In summary, the proposed response coefficients in this work are novel in two primary respects. First of all, the response coefficients utilize spatial–temporal triplets, which allows for the study of nonlinear interactions between structures of different sizes and frequencies. More importantly, the response coefficients incorporate the linear resolvent operator, which accounts for the effect of linear mechanisms when studying the nonlinear triadic energy transfer. This study of nonlinear energy transfer while taking into account of the effects of the linear mechanism has not been performed previously to the best of our knowledge.

4.1. Three-dimensional analyses of forcing and response coefficients

For fixed $\boldsymbol{k}_3$ , the six-dimensional forcing and response coefficients $P(\boldsymbol{k}_1,\boldsymbol{k}_2)$ and $R(\boldsymbol{k}_1,\boldsymbol{k}_2)$ defined in (2.17) and (2.23) collapse into three-dimensional quantities, which can be computed and stored with reasonable resources, and can reveal information about the spatiotemporal nature of dominant nonlinear interactions. The magnitude of the forcing and response coefficients are plotted in figure 8 as functions of the contributing scale $\boldsymbol{k}_1$ , noting the (strong) constraint that $\boldsymbol{k}_2 = \boldsymbol{k}_3 - \boldsymbol{k}_1$ .

A first observation is that the majority of the high coefficient values are associated with $k_{x1} \approx 1$ , indicating that the large-scale structures play a dominant role in both forcing and response. The largest coefficient magnitudes in the figures are $2.2{\times 10^{-10}}$ for $P$ , and $2.0{\times 10^{-9}}$ for $R$ , reflecting that, while near-wall modes are known to be energetic in a local event, their limited wall-normal reach restricts the integrated energy of such modes. The energy of the selected mode of $\boldsymbol{k}_3 = [4, \ 28, \ 2.492]$ is $1.3{\times 10^{-8}}$ . As a result, the largest magnitude of $R$ , $2.0{\times 10^{-9}}$ corresponds to 15 % of the total modal energy level, and thus the associated nonlinear interaction is a major contributor towards this representative mode of the near-wall cycle.

Further, it can be observed that most of the high magnitude coefficients (both forcing and response) reside near a single plane with an almost constant $c_1 = \omega _1 / k_{x1}$ . This is especially evident in figure 9, which show the same data summed in $k_{z1}$ and plotted as contour plots in the $k_{x1}-\omega _1$ plane. Three different wave speeds are also shown for reference: $c_1 = \omega _1/k_{x1} = 1$ , where the wave speed matches the centreline velocity; $c_1 = c_3 = 0.63$ , the wave speed of the selected $\boldsymbol{k}_3$ ; and $c_1 = 0.3$ ( $c_1^+ = 6$ ). It can be observed that most of the energetic regions are bounded between $c_1^+ = 6$ and $c_1 = 1$ , centring roughly around $c_1 = c_3 = 0.63$ . Note that $c_1 = c_3$ if and only if $c_2 = c_3$ from the triadic compatibility constraint (except for $k_x = 0$ where $c$ is no longer well defined).

Figure 9.

Magnitude of (a) $\sum _{k_{z1}} P(\boldsymbol{k}_1, \boldsymbol{k}_3 - \boldsymbol{k}_1)$ , the forcing coefficients summed in $k_{z1}$ and plotted as functions of $k_{x1}$ and $\omega _1$ for $\boldsymbol{k}_3 = [4, \ 28, \ 2.492]$ , and (b) $\sum _{k_{z1}} R(\boldsymbol{k}_1, \boldsymbol{k}_3 - \boldsymbol{k}_1)$ , the response coefficients. The black dash–dotted lines in both figures mark the wave speed $c_1 = \omega /k_x = 1$ ; the black dashed lines for $c_1 = c_2 = c_3 = 0.63$ , the wave speed of the selected $\boldsymbol{k}_3$ ; and the black dotted lines for $c_1 = 0.3$ ( $c_1^+ = 6$ ). The + markers in both figures denote the locations of $\boldsymbol{k}_3$ . Panel (a) includes an additional set of red lines, with the red dash–dotted line for $c_2 = 1$ and red dotted line for $c_2 = 0.3$ ( $c_2^+ = 6$ ).

To understand the reason behind the dominance of triadic interactions with $c_1 \approx c_2 \approx c_3$ , we examine the energy (rather than the projections) of the forcing and response generated by the triads $\boldsymbol{k}_1$ and $\boldsymbol{k}_2 = \boldsymbol{k}_3 - \boldsymbol{k}_1$ . In figure 10(a), the energy $E_{\kern-1pt f}(\boldsymbol{k}_1, \boldsymbol{k}_2)$ and figure 10(b), the energy $E_u(\boldsymbol{k}_1, \boldsymbol{k}_2)$ are plotted in the $k_{x1}-\omega _1$ plane similar to figure 9 with the same wave speeds marked. Here $E_u(\boldsymbol{k})$ , the $y$ -integrated TKE of the DNS velocity Fourier modes, as defined in (3.2), is shown in figure 10(c) to provide context of the energy distribution across all interactions. Comparing the three figures, it can be observed that the energetic regions of the triadically generated forcing $E_{\kern-1pt f}(\boldsymbol{k}_1, \boldsymbol{k}_2)$ are located in the vicinity of $c_1 = c_2 = c_3$ , while the triadically generated response $E_u(\boldsymbol{k}_1, \boldsymbol{k}_2)$ lie even closer to $c_1 = c_2 = c_3$ . By contrast, the most energetic regions of full DNS spectrum in figure 10(c) are located much closer to $c = 1$ , reflecting the energetic dominance of fast-moving modes that have a large wall-normal extent and therefore a large integral contribution to the energy (see, e.g. Bourguignon et al. Reference Bourguignon, Sharma, Tropp and McKeon2014).

Figure 10.

Contour plots of (a) $\sum _{k_{z1}} E_{\kern-1pt f}(\boldsymbol{k}_1, \boldsymbol{k}_3 - \boldsymbol{k}_1)$ , the forcing energy, and (b) $\sum _{k_{z1}} E_u(\boldsymbol{k}_1, \boldsymbol{k}_3 - \boldsymbol{k}_1)$ , the response energy, generated by the interactions between $\boldsymbol{k}_1$ and $\boldsymbol{k}_3 - \boldsymbol{k}_1$ summed in $k_{z1}$ and plotted as a function of $k_{x1}$ and $\omega _1$ for $\boldsymbol{k}_3 = [4, \ 28, \ 2.492]$ . Panel (c) is the contour of $\sum _{k_{z}} E_u(\boldsymbol{k})$ , the $y$ integrated kinetic energy of the DNS Fourier modes summed in $k_z$ and plotted as a function of $k_{x}$ and $\omega$ . The three lines mark the same wave speeds as in figure 9: dash–dotted lines for $c_1 = 1$ , dashed lines for $c_1 = c_2 = c_3 = 0.63$ , and dotted lines for $c_1^+ = 6$ .

The triadically generated forcing $\boldsymbol{f}(\boldsymbol{k}_1,\boldsymbol{k}_2,y)$ is the product between the velocity modes at $\boldsymbol{k}_1$ and the gradient of velocity modes at $\boldsymbol{k}_2$ . The spatial $y$ -localization of the modes around their respective critical layers shown in figure 3 provides one simple insight into the dominant response when modes with $\boldsymbol{k}_1$ , $\boldsymbol{k}_2$ and $\boldsymbol{k}_3$ are colocated in $y$ .

4.2. Relative importance of forcing weights and the action of the resolvent

In this section, we utilize the linear resolvent operator to assess the relative importance of the forcing weights on individual resolvent modes in a triadic interaction and the action of the resolvent in generating the velocity response.

Following § 2, the nonlinear weights $\chi _q(\boldsymbol{k}_3)$ obtained by projecting the forcing Fourier modes $\boldsymbol{f}(\boldsymbol{k}_3,y)$ onto the resolvent forcing modes $\phi _q(\boldsymbol{k}_3,y)$ are defined as

(4.4) \begin{align} \chi _q(\boldsymbol{k}_3) = \left \langle \phi _q(\boldsymbol{k}_3,y), \ \boldsymbol{f}(\boldsymbol{k}_3,y) \right \rangle _y. \end{align}

Similarly, the contribution to the weight on the $q$ th resolvent forcing mode from individual triadic interactions, $\boldsymbol{f}(\boldsymbol{k}_1,\boldsymbol{k}_2,y)$ , can also be defined as

(4.5) \begin{align} \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_2) = \left \langle \phi _q(\boldsymbol{k}_1+\boldsymbol{k}_2,y), \ \boldsymbol{f}(\boldsymbol{k}_1,\boldsymbol{k}_2,y) \right \rangle _y. \end{align}

Thence, forcing and response associated with an individual triadic interaction can be written in terms of weighted resolvent forcing modes as

(4.6a) \begin{align} \boldsymbol{f}(\boldsymbol{k}_1,\boldsymbol{k}_2,y) &= \sum _q \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_2) \phi _q(\boldsymbol{k}_1+\boldsymbol{k}_2,y), \\[-10pt] \nonumber \end{align}

(4.6b) \begin{align} \boldsymbol{\tilde {u}}(\boldsymbol{k}_1,\boldsymbol{k}_2,y) &= \sum _q \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_2) \sigma _q(\boldsymbol{k}_1+\boldsymbol{k}_2)\psi _q(\boldsymbol{k}_1+\boldsymbol{k}_2,y). \\[10pt] \nonumber \end{align}

Finally, the following property can also be obtained from (4.3a ):

(4.7) \begin{align} \chi _q(\boldsymbol{k}_3) = \sum _{\boldsymbol{k}_1} \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_3-\boldsymbol{k}_1). \end{align}

Thus, $\chi _q(\boldsymbol{k}_3)$ describes the amount of the $q$ th resolvent forcing mode contained in the DNS forcing Fourier mode at $\boldsymbol{k}_3$ , while $\chi _q(\boldsymbol{k}_1,\boldsymbol{k}_3-\boldsymbol{k}_1)$ is the amount contributed towards the $q$ th resolvent forcing mode by the interaction between $\boldsymbol{k}_1$ and $\boldsymbol{k}_2 = \boldsymbol{k}_3 - \boldsymbol{k}_1$ . Similarly $\sigma _q(\boldsymbol{k}_3)\chi _q(\boldsymbol{k}_3)$ describes the amount of the $q$ th resolvent response mode contained in the DNS velocity Fourier mode at $\boldsymbol{k}_3$ , while $\sigma _q(\boldsymbol{k}_3)\chi _q(\boldsymbol{k}_1,\boldsymbol{k}_3-\boldsymbol{k}_1)$ is the amount contributed towards the $q$ th resolvent response mode by the interaction between $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ .

Applying these definitions in $P(\boldsymbol{k}_1, \boldsymbol{k}_2)$ in (2.17) and $R(\boldsymbol{k}_1, \boldsymbol{k}_2)$ in (2.23) and for the energy in (4.2a ) and (4.2b ) yields

(4.8a) \begin{align} P(\boldsymbol{k}_1, \boldsymbol{k}_2) &= \sum _q \mathbb{E}\big \{\chi _q^*(\boldsymbol{k}_1+\boldsymbol{k}_2) \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_2)\big \}, \\[-10pt] \nonumber \end{align}

(4.8b) \begin{align} R(\boldsymbol{k}_1, \boldsymbol{k}_2) &= \sum _q \sigma _q^2(\boldsymbol{k}_1+\boldsymbol{k}_2) \mathbb{E}\big \{\chi _q^*(\boldsymbol{k}_1+\boldsymbol{k}_2) \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_2)\big \} \\[10pt] \nonumber \end{align}

and

(4.9a) \begin{align} E_{\kern-1pt f}(\boldsymbol{k}_1, \boldsymbol{k}_2) &= \sum _q \mathbb{E}\big \{ \big | \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_2) \big |^2\big \}, \\[-10pt] \nonumber \end{align}

(4.9b) \begin{align} E_u(\boldsymbol{k}_1, \boldsymbol{k}_2) &= \sum _q \sigma _q^2(\boldsymbol{k}_1+\boldsymbol{k}_2) \mathbb{E}\big \{ \big | \chi _q(\boldsymbol{k}_1,\boldsymbol{k}_2) \big |^2\big \}. \\[10pt] \nonumber \end{align}

Note that $\sigma _q$ can be moved outside of the expected value operator as these are quantities computed from the deterministic resolvent operator (cf. Towne et al. Reference Towne, Schmidt and Colonius2018). The linear resolvent operator amplifies the different forcing modes differently, and the effect is captured in the response coefficients $R(\boldsymbol{k}_1, \boldsymbol{k}_2)$ and response energy $E_u(\boldsymbol{k}_1, \boldsymbol{k}_2)$ through the magnitudes of the $\sigma _q$ .

For the near-wall cycle mode, $\boldsymbol{k}_3 = [4, \ 28, \ 2.492]$ , it can be observed from figure 11(a) that the first two singular values are larger than the rest, showing that the linear resolvent operator predominately amplifies the first two modes and is low-rank at this $\boldsymbol{k}_3$ . Figure 11(b) shows that the nonlinear weights $\chi _q(\boldsymbol{k}_3)$ for the first few modes are the smallest and thus contributing only a small portion towards the forcing Fourier mode at this $\boldsymbol{k}_3$ . However, due to the strong amplification of the linear resolvent operator, these first few modes contribute significantly towards the velocity Fourier modes, as evident in figure 11(c), where it is shown that the first two modes have the largest $\left | \sigma _q\chi _q \right |$ . The forcing Fourier mode receives contributions from a wide range of resolvent forcing modes, while the velocity Fourier mode is dominated mainly by the first few resolvent response modes due to the low-rank nature of the resolvent.

These results are consistent with the work of Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021), where it is shown that the forcing has significant projection onto the suboptimal resolvent forcing modes. In fact, the authors showed that the projections onto the first few optimal forcing modes are smaller than the rest, as also observed here. On the other hand, the bulk of the velocity responses is shown to be well approximated using a rank-two approximation with the first two resolvent modes, which can also be observed in this study.

Rosenberg & McKeon (Reference Rosenberg and McKeon2019) exploited the Helmholtz decomposition to demonstrate that only the solenoidal component of the forcing effects a velocity response (since the dilatational part can be absorbed into a modified pressure). A large forcing weight, $\chi _q$ , can be associated with small solenoidal content of high $q$ forcing modes, leading to a small velocity response as seen here.

Table 2.

Three selected triads that contribute to $\boldsymbol{k}_3 = \boldsymbol{k}_1+\boldsymbol{k}_2 = [4, \ 28, \ 2.492]$ , with the magnitude of the forcing and response coefficients $| P(\boldsymbol{k}_1, \boldsymbol{k}_2) |$ , $| R(\boldsymbol{k}_1, \boldsymbol{k}_2) |$ and the energy of the triadically generated forcing and response $E_{\kern-1pt f}(\boldsymbol{k}_1, \boldsymbol{k}_2)$ , $E_u(\boldsymbol{k}_1, \boldsymbol{k}_2)$ .

Figure 11.

Plots of (a) the resolvent singular values $\sigma _q(\boldsymbol{k}_3)$ , (b) the power spectral density of the nonlinear weights $\chi _q(\boldsymbol{k}_3)$ , obtained by taking the inner product between the resolvent forcing modes and the DNS forcing Fourier mode, and (c) the product of the two as a function or the resolvent mode number $q$ .

Consider the three example pairs of $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ that are triadically compatible with $\boldsymbol{k}_3 = [4, \ 28, \ 2.492]$ listed in table 2 and the associated forcing and response magnitudes in figure 12. Triad 1 has strong forcing and response magnitudes. Triad 2 generates a relatively large response from a small forcing, and triad 3 generates a small response even though the forcing magnitude is large.

Figure 12.

Plots of (a) nonlinear weights $\mathbb{E}\{| \chi _q(\boldsymbol{k}_1, \boldsymbol{k}_2) |^2\}$ of triadic interactions between $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ , and (b) the nonlinear weights multiplied by the resolvent singular values $\mathbb{E}\{| \sigma _q(\boldsymbol{k}_1+\boldsymbol{k}_2)\chi _q(\boldsymbol{k}_1, \boldsymbol{k}_2) |^2\}$ as a function of the resolvent mode number $q$ . Blue markers are for triad 1 in table 2, red markers for triad 2 and yellow markers for triad 3.

For this near-wall mode, the full forcing is shown to have significant projection onto suboptimal modes, while the full response is relatively low-rank. This is also true on an interaction-by-interaction basis and thus may provide insight into the modelling of Reynolds stress gradients.

The results presented in this section could potentially benefit from additional averaging with the results for negative $k_z$ , as the positive and negative $k_z$ pairs correspond to left and right oblique waves and are expected to be statistically symmetric. However, care must be taken for results concerning individual short temporal segments, as the symmetry might only hold statistically. Finally, averaging with the results for negative $k_x$ or negative $\omega$ is not expected to be valid as those correspond to upstream travelling waves, which are different from the downstream travelling waves studied here.

5.1. Triadic interactions in the streamwise direction

Consider first the streamwise forcing and response coefficients $\left | P_{k_x}(k_{x1}, k_{x2}) \right |$ and $\left | R_{k_x}(k_{x1}, k_{x2}) \right |$ in figure 14(a,e). A dominant horizontal band is observed for the forcing. This horizontal band, centred around $k_{x1} = 0$ , corresponds to the interaction between the streamwise large-scale velocity modes with velocity gradients at all scales, depicted in the Feynman diagram in figure 13(a). By comparison, the same band is present in the response, figure 14(e), but with lower amplitude relative to other types of interaction. This is attributable to the action of the resolvent on the forcing from figure 14(a).

Figure 14.

Heatmaps of the magnitude of (a–d) the streamwise forcing coefficients $| P_{k_x}(k_{x1}, k_{x2}) |$ , and (e–h) the streamwise response coefficients $| R_{k_x}(k_{x1}, k_{x2}) |$ . Panels (a)–(d) except the insert share the same logarithmic colourbar, same for (e)–(h). The insert in (a) and (e) corresponds to representations of the rectangular regions enclosed in the dashed white lines using linear scale colourbars. The $y$ -integration limits for the inner product in (2.18) are: (a,e) all $y^+$ ; (b,f) $y^+ \in (0, 30)$ ; (c,g) $y^+ \in (30, 200)$ ; and (d,h) $y^+ \in (200, 550)$ . The streamwise wavenumber for the velocity fields, $k_{x1}$ , is on the vertical axis; $k_{x2}$ for the velocity gradient is on the horizontal axis; $k_{x3} = k_{x1} + k_{x2}$ for the resulting forcing and response is constant along lines with slopes of $-1$ . The diagonal white dotted lines in (a) and (e) and the black dashed lines in the inserts mark the location of $k_{x3} = 4$ ( $\lambda _x^+ = 865$ ).

For the forcing coefficient in figure 14(a–d), the high values in the horizontal band reflect the concentration of energy in the streamwise large-scale structures at $k_{x1} = 0$ and $\pm 0.5$ . These modes are also tall, i.e. they have a large wall-normal extent as shown in figure 3, and therefore are capable of interacting with $k_{x2}$ modes of any size, centred at any $y$ location. The combined effect leads to a significant amount of forcing energy generated across a range of scales by these large-scale modes, manifesting as the energetic horizontal band. Additionally, the values in this horizontal band decay relatively slowly as $\left | k_{x2} \right |$ increases. The decreasing energy associated with velocity modes with increasing $\left | k_{x2} \right |$ is partially compensated by increasing magnitude of their spatial gradients.

Finally, comparing the wall-normal variation of the horizontal, $k_{x1}=0$ , band across figure 14(b–d), it can be observed that $P_{k_x}$ exhibits a significant presence across all three $y$ ranges due to the tall large-scale modes. The strength of the interactions in the wake region decays, however, for large $\left | k_{x2} \right |$ , (the interaction type shown in figure 13 f), reflecting the reduction in energy at small scales far from the wall.

The response coefficients in figure 14(e–h) show that the horizontal, $k_{x1} = 0$ band does not correspond to the most energetic response. The strong forcing present for $k_{x1} = 0$ at large $k_{x3}$ results in little response energy. This showcases the effect of the linear resolvent operator, which amplifies the small scales less than the large scales and is also consistent with the cascade picture where energy is transferred from the large scales to the small scales through nonlinear interactions and dissipated at the small scales through viscosity (Jiménez Reference Jiménez2012).

Focusing instead on the different contributors towards the small scales (large $k_{x3}$ ), it can be observed that the nonlinear interactions between large-scale velocity modes (small $\left | k_{x1} \right |$ ) and small-scale velocity gradients (large $\left | k_{x2} \right |$ ) (figure 13 f) dominate the response. These long-range interactions (in $k_x$ ) indicate coherence between the large and small scales, consistent with the superposition and modulation mechanisms observed in Marusic et al. (Reference Marusic, Mathis and Hutchins2010) and Andreolli et al. (Reference Andreolli, Gatti, Vinuesa, Örlü and Schlatter2023).

Lower values of $\left | P_{k_x} \right |$ and $\left | R_{k_x} \right |$ are observed for $k_{x2}=0$ (vertical band, representing interactions of the type shown in figure 13 b). The modes contributing to the velocity and velocity gradient are reversed in (2.19a ) and (2.24a ) relative to the $k_{x1}=0$ , horizontal band. Spatial derivatives for the large scale are weaker and the energy of the small scales is lower than for the horizontal band, underscoring the asymmetry of interactions within a given triad.

For the diagonal band with $k_{x3} \approx 0$ (interactions of type figure 13 c), a weak signature of interactions between two similar size modes generating a streamwise large-scale forcing can be observed for the forcing coefficients in figure 14(a–d). For the response coefficients in figure 14(e–h), this diagonal band is stronger, and is the most dominant band of $R_{k_x}$ , due to the preferential amplification of the large-scale modes by the linear resolvent operator.

Within this diagonal region, forcing coefficient $P_{k_x}$ for interactions between two large scales in the central region (interactions shown in figure 13 d) is stronger than the interactions between two small scales in the corner regions (especially evident in the linear scale inserts). This result is consistent with the findings of Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021), where it is demonstrated that large-scale forcing modes are mostly the result of interactions by large-scale structures, although with a relatively weak influence from the interactions between the smaller scales. In addition, the forcing components generated from small-scale interactions are almost non-existent in the outer region (shown in figure 14 d), due to the small scales having little energy presence far from the wall.

On the contrary, the response coefficient $R_{k_x}$ , with an energetic diagonal band extending all the way to the small scales in figure 14(e–h), shows that this type of interaction between small scales affecting the large scales (figure 13 e), has a non-negligible albeit relatively weaker influence on the velocity response. This phenomenon is also demonstrated in the work of Illingworth, Monty & Marusic (Reference Illingworth, Monty and Marusic2018), where it is shown that including an eddy viscosity into the resolvent framework, which is intended to model the effect of this type of interaction, improves the performance of the resolvent analysis for the large scales. It should be pointed out that figure 14 examines only the magnitude of the interaction coefficients, where large values show strong importance. The direction of energy transfer, whether energy injection or extraction will be examined later using the phase of the coefficients, and it will be shown that this type of interaction between small scales weakens the spectral TKE of the the large scales, consistent with the energy cascade (Jiménez Reference Jiménez2012), and QL/GQL dynamics, where the Reynolds stress divergence from the small scales modifies the streamwise-mean flow in such a way as to render the large scales more stable and likely less energetic.

Finally, it should be noted that the fact that the three bands have peak values at $k_{x1}, k_{x2}, k_{x3} = 0$ is an artefact of the current numerical simulation. With the current $x$ domain length of $4\pi$ , energy from the unresolved large scales manifests as streamwise constant structures at $k_x = 0$ , making it the most energetic wavenumber. It should also be pointed out that these are still the mean-subtracted velocity fluctuations, that are streamwise constant, yet non-constant in spanwise direction ( $k_z \neq 0$ ) and/or non-constant in time ( $\omega \neq 0$ ). As the box size approaches infinity, providing higher and higher resolution for the streamwise wavenumber, these bands are expected to display some width in $k_x$ with multiple wavenumbers very close to $\pm k_{xl}$ , the streamwise wavenumber of the most energetic large-scale streak, presumably the very large-scale motion, displaying high levels of importance.

Figure 15.

Heatmaps of the magnitude of (a–d) the spanwise forcing coefficients $| P_{k_z}(k_{z1}, k_{z2}) |$ , and (e–h) the spanwise response coefficients $| R_{k_z}(k_{z1}, k_{z2}) |$ in the same format as figure 14. The diagonal white dotted lines in (a) and (e) and the black dashed lines in the inserts mark the location of $k_{z3} = 28$ ( $\lambda _z^+ = 124$ ).

5.2. Triadic interactions in the spanwise direction

The spanwise forcing and response coefficients $\left | P_{k_z}(k_{z1}, k_{z2}) \right |$ and $\left | R_{k_z}(k_{z1}, k_{z2}) \right |$ in figure 15 behave very similarly to the streamwise coefficients, with the exception that the three bands no longer show single-banded structures located around $k_{x}\approx 0$ , but instead show dual-band structures located at $k_{z} \approx \pm 3$ . This is as expected since the large-scale structures, which are the most energetic modes in the flow field, with $k_x = 0$ and $k_x = 0.5$ are also most energetic at $k_z \approx \pm 3$ . The dual-banded structures are prominent for the energetic horizontal band of $P_{k_z}$ and diagonal band of $R_{k_z}$ , while less obvious for the less energetic bands.

The teardrop shaped region at the highest $k_{z2}$ and $k_{z3}$ values observed in the forcing coefficient $P_{k_z}$ in figure 15(a–d) is likely an artefact of convergence, although note that it is not observed in the response coefficient $R_{k_z}$ . The forcing coefficients in this region have magnitudes that are orders of magnitude smaller than the peak values, making them sensitive to numerical errors. Although the corner areas of the results involving extreme $k_z$ values are less robust quantitatively, the overall structure of the results is expected to remain consistent.

5.3. Triadic interactions of the temporal frequencies

The overall structure of the frequency forcing and response coefficients $\left | P_{\omega }(\omega _{1}, \omega _{2}) \right |$ and $\left | R_{\omega }(\omega _{1}, \omega _{2}) \right |$ plotted in figure 16 is again similar to that of the streamwise coefficients, with the exception that the single-banded structures in the streamwise direction now becomes multibanded. These prominent discrete high-value lines located in the horizontal, vertical and diagonal bands with $\omega _1, \omega _2, \omega _3 \approx 0$ are due to the discreteness in the streamwise wavenumber $k_x$ , an artefact of the finite simulation domain length, as similarly demonstrated in Gómez et al. (Reference Gómez, Blackburn, Rudman, McKeon, Luhar, Moarref and Sharma2014). The frequency can be related to the phase speed for a given $\boldsymbol{k}$ via $\omega = ck_x$ ; at a given wave speed, $c$ , an increase in $k_x$ to the next discrete wavenumber, with an increment of $0.5$ fixed by the simulation domain length, will cause $\omega$ to increase by $0.5c$ . Therefore, increments in $\omega$ are largest for large scales with high wave speeds, appearing as discrete lines, and smaller for small scales with low wave speeds, reflected as the smooth varying background.

Figure 16.

Heatmaps of the magnitude of (a–d) the temporal forcing coefficients $| P_{\omega }(\omega _{1}, \omega _{2}) |$ , and (e–h) the temporal response coefficients $| R_{\omega }(\omega _{1}, \omega _{2}) |$ in the same format as figure 14. The diagonal white dotted lines in (a) and (e) and the black dashed lines in the inserts mark the location of $\omega _{3} = 2.492$ ( $\omega _3^+ = 0.0953$ ).

A more detailed analysis reveals that the five most prominent lines observed in figure 16 for both $P_{\omega }$ and $R_{\omega }$ are located at $\omega \approx 0$ , $\pm 0.4$ and $\pm 0.8$ . These lines correspond well with the energetic modes at $k_x=0$ , $k_x=\pm 0.5$ (with wave speeds of $c = \omega /k_x \approx 0.8$ shown in figure 3 a) and $k_x= \pm 1$ (with $c \approx 0.8$ shown in figure 3 b), respectively. Furthermore, by observing the relative intensities between the smooth background and discrete lines across different $y$ locations, it can be seen that the former is more prominent near the wall in figures 16(b) and 16(f) where the discreteness is barely visible, while away from the wall in figures 16(d) and 16(h) the opposite is true. Combined with the tall $y$ extent for large scales and the concentration of energy near the wall for small scales, as demonstrated in figure 3, these plots confirm that the smooth background mostly shows the triadic interactions between the small scales, while the discrete lines are mostly the result of interactions involving the large scales. Although the discreteness is an artefact of the finite simulation domain length, interactions involving large scales are expected to be important regardless of the domain length.

5.4. Constructive and destructive triadic interactions

Both the forcing and response coefficients are complex numbers, with the magnitude providing information about the importance of a triadic interaction, while the phase provides information about constructive and destructive interference. From (2.27a ) to (2.27c ), it can be observed that the response coefficients sum to yield the spectral TKE of the response modes, which are real positive quantities. Therefore, the imaginary parts of the coefficients cancel upon summation, while the real parts provide information about the constructive and destructive contributions to the spectral TKE. A positive real part of $R(\boldsymbol{k}_1,\boldsymbol{k}_2)$ indicates that the interaction between $\boldsymbol{k}_1$ and $\boldsymbol{k}_2$ causes an increase in spectral TKE of the response mode at $\boldsymbol{k}_3$ , while a negative real part indicates a decrease of spectral TKE. To analyse the constructive and destructive interference, we utilized the phase angle of the response coefficients: a phase angle within $(-\pi /2, +\pi /2)$ or the right half-plane of the complex plane indicates a positive real part with both positive or negative imaginary part; while a phase angle within $(+\pi /2,+\pi ) \cup (-\pi ,-\pi /2)$ or the left half-plane indicates a negative real part. Since the sign of the imaginary part is of less interest, we utilize the absolute value of the phase angles, where $\left | \angle R \right | \in [0,\pi /2)$ indicates constructive interference or in other words, increase of spectral TKE, and $\left | \angle R \right | \in (\pi /2,\pi ]$ indicates destructive interference or in other words, decrease of spectral TKE. In figure 17, the absolute value of the phase of $R_{k_x}$ , $R_{k_z}$ and $R_{\omega }$ are plotted, with phase angles close to 0 in red indicating constructive interference and phase angles close to $\pi$ in blue indicating destructive interference.

Figure 17.

Heatmaps of the absolute values of the phase angles of the response coefficients: (a) $| \angle R_{k_x}(k_{x1}, k_{x2})|$ , (b) $| \angle R_{k_z}(k_{z1}, k_{z2}) |$ and (c) $| \angle R_{\omega }(\omega _{1}, \omega _{2}) |$ . Phase angles close to 0 (red) indicate constructive interference, while phase angles close to $\pi$ (blue) indicate destructive interference.

It should be noted that the constructive and destructive interference studied here are different from typical studies of energy transfer focusing on the turbulent transport term. As noted in (2.14), a nonlinear energy transfer into a mode at $\boldsymbol{k}_3$ by turbulent transport does not necessarily result in an increase in spectral TKE due to the presence of other linear mechanisms. However, in this study, with both the nonlinear energy transfer and linear energy amplification mechanisms studied together, constructive or destructive interference directly indicates an increase or decrease of the spectral TKE, providing a different perspective compared with other energy transfer studies.

The overall structure of the phase plots can be seen as a red hour-glass structure spanning from the bottom-left to the top-right, which is obvious across all three panels, and a blue hour-glass structure spanning from the top-left to bottom right, which is obvious for $R_{k_x}$ and $R_{\omega }$ , but exhibits more complex behaviour for $R_{k_z}$ . To analyse this structure, we will focus on lines of constant $\boldsymbol{k}_3$ , which are lines with slopes of $-1$ , parallel to the dashed lines marking the extreme values of $\boldsymbol{k}_3$ . Equation (2.27a ) indicates that summing along this line of constant $k_{x3}$ gives the energy of all modes at $k_{x3}$ . Within this line, the central region with small $k_{x1}$ and $k_{x2}$ , representing the interactions between the large scales, generally contributes positively to the energy, while the corner regions, representing the interactions between small scales, generally reduce energy. This is consistent with the energy cascade, where TKE is being generated at the large scales, transferred to small scales through triadic interactions and dissipated at the small scales by viscosity. For the phase of $R_{\omega }$ , the general structure behaves similarly. However, the central region is fuzzier while generally remaining red, contributing positively towards the energy. This is likely due to the fact that, unlike $k_x$ where only large scales contribute towards low $k_x$ , for $\omega$ , both large scales with high wave speeds and small scale with low wave speeds can contribute to low $\omega$ . This phenomenon results in the coexistence of the high-value discrete lines and the smooth varying background in the magnitude of $R_{\omega }$ discussed in the previous section, and most likely contributes to the fuzziness in the central region. For the phase of $R_{k_z}$ , the general structure remains similar. It is observed that the overall structure of these results is robust and sufficiently converged; however, the details of the transitional regions between constructive and destructive interference require more data to reduce the variations for future studies focusing on these regions.