1. Introduction
Energy transfer plays a key role in the organisation and evolution of turbulent flows. It is responsible for the multi-scale nature of turbulence through the Richardson–Kolmogorov turbulent energy cascade (Kolmogorov Reference Kolmogorov1941) and lends insight into the self-sustaining process (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997). Energy transfer for an individual scale is described by the spectral turbulent kinetic energy (TKE) equation, which contains a nonlinear term sometimes referred to as turbulent transport. As noted by Domaradzki et al. (Reference Domaradzki, Liu, Härtel and Kleiser1994), the nonlinearity poses considerable theoretical difficulties by permitting inter-scale energy exchange. It is not possible, for example, to study a scale in isolation without a closure model and, in the context of a large eddy simulation (LES), subgrid models need to account for the influence of small scales on the large scales of interest. An improved understanding of nonlinear interactions in turbulent flows, therefore, is essential to improve turbulence modelling and simulation.
It is also known that linear mechanisms are important in energy transfer. These are described well by the linear operator obtained after linearising the Navier–Stokes equations around a suitable base flow (Schmid & Henningson Reference Schmid and Henningson2001). This operator is highly non-normal due to the mean shear found in wall-bounded flows (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). As a result, infinitesimal disturbances may experience significant transient growth by extracting energy from the mean shear (Butler & Farrell Reference Butler and Farrell1992; Reddy & Henningson Reference Reddy and Henningson1993). Linear mechanisms have also been identified in mean (time-averaged) flows by the resolvent analysis of McKeon & Sharma (Reference McKeon and Sharma2010). In this framework, the equations are linearised around the mean flow to obtain the resolvent operator that maps the nonlinear terms, treated as an intrinsic forcing, to the velocity in the frequency domain. Energy production by the linear system is excited by the nonlinear forcing which, when large enough, also excites dissipative modes to dissipate energy (Sharma Reference Sharma2009). The nonlinear forcing is itself composed of quadratic interactions between various outputs of the linear amplification process to complete a feedback loop (McKeon, Sharma & Jacobi Reference McKeon, Sharma and Jacobi2013). When the flow is dominated by a single (fundamental) Fourier mode, the nonlinear forcing is dominated by interactions between the fundamental mode and its harmonics (Rosenberg, Symon & McKeon Reference Rosenberg, Symon and McKeon2019). In a turbulent channel flow, however, it is less tractable to isolate the principal nonlinear interactions that make up the nonlinear forcing since many Fourier modes are energetic and interact with each other.
An objective of this paper, therefore, is to investigate the extent to which energy transfer is correctly modelled by resolvent analysis. To address this question, we first examine how energy is produced, dissipated and transferred among various scales in turbulent channel flow at low Reynolds numbers. Similar to other studies (Dar, Verma & Eswaran Reference Dar, Verma and Eswaran2001; Mizuno Reference Mizuno2016; Cho, Hwang & Choi Reference Cho, Hwang and Choi2018; Lee & Moser Reference Lee and Moser2019), we calculate these terms in spectral space and integrate them over the wall-normal direction. The true energy transfer from DNS is compared to predictions from the optimal resolvent mode, which is often representative of the true velocity field observed in DNS or experiments (McKeon Reference McKeon2017). The agreement can be improved by adding the Cess (Reference Cess1958) eddy viscosity profile to the resolvent operator as has been done in many studies (Hwang & Cossu Reference Hwang and Cossu2010; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019; Symon, Illingworth & Marusic Reference Symon, Illingworth and Marusic2020). This has been likened to a crude model for the energy cascade by Hwang (Reference Hwang2016), suggesting that it assumes the role of turbulent transport in resolvent analysis. To provide insight into the matter, we quantify the contribution of eddy viscosity to the energy balance for each scale and compare it to the true nonlinear transfer. Finally, we investigate the role of non-normality, quantified by the inner product between the optimal resolvent forcing and response modes, in the energy balance. The influence of non-normality is reduced by eddy viscosity leading to an improvement in predictions from resolvent analysis.
The flows selected for this study are the P4U exact coherent state (ECS) of Park & Graham (Reference Park and Graham2015) and turbulent flow in a minimal channel at low Reynolds number. The former is a nonlinear travelling wave whose mean properties resemble those of near-wall turbulence. It is a particularly appealing choice to study energy transfer since it is low-dimensional (Sharma et al. Reference Sharma, Moarref, McKeon, Park, Graham and Willis2016; Rosenberg & McKeon Reference Rosenberg and McKeon2019) and travels at a fixed convection velocity. As such, all computations can be performed on a standard personal computer, and no integration in time is necessary to obtain each term in the energy budget. To verify that the transfer mechanisms are similar in a time-evolving flow, we compare the results for the P4U ECS to those of more standard turbulence in a ‘minimal flow unit’ (Jiménez & Moin Reference Jiménez and Moin1991).
The paper is organised as follows. In § 2, the relevant equations for resolvent analysis, energy transfer and the eddy viscosity model are derived. The simulation parameters for the P4U ECS and minimal channel flows are described in § 3. The energy balances computed from DNS and resolvent analysis are compared for the ECS in § 4 and for the minimal channel in § 5. In § 6, we examine the influence of non-normality on the ability of the first resolvent mode to describe energy transfer processes. We also analyse the role of eddy viscosity on the efficiency of resolvent modes as a basis for the velocity fluctuations. This leads to a discussion on the typical scales for which an eddy viscosity leads to an improvement before we conclude in § 7.
2. Methods
In § 2.1, we describe the governing equations for plane Poiseuille flow and their non-dimensionalisation. A brief overview of resolvent analysis is provided in § 2.2. The energy balance for each scale is then derived from the fluctuation equations in § 2.3, and we show that this balance is maintained for each resolvent mode. Finally, we describe the eddy viscosity model in § 2.4.
2.1. Plane Poiseuille flow equations
The non-dimensional Navier–Stokes equations for statistically steady, turbulent plane Poiseuille flow are
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} = - \boldsymbol{\nabla} p + \frac{1}{Re_{\tau}}{\nabla}^2\boldsymbol{u}, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
where 
$\boldsymbol {u}(\boldsymbol {x},t) = [u,v,w]^{\textrm {T}}$ is the velocity in the 
$x$ (streamwise), 
$y$ (spanwise) and 
$z$ (wall-normal) directions, 
$p(\boldsymbol {x},t)$ is the pressure and 
$\boldsymbol {\nabla } = [\partial /\partial x,\partial /\partial y, \partial /\partial z]^\textrm {T}$. The friction Reynolds number 
$Re_{\tau } = u_\tau h /\nu$ is defined in terms of the friction velocity 
$u_{\tau }$, channel half height 
$h$ and kinematic viscosity 
$\nu$. No-slip boundary conditions are applied at the walls, and periodic boundary conditions are imposed in the streamwise and spanwise directions. The density of the fluid is 
$\rho$ and the velocities are non-dimensionalised by 
$u_\tau$, the spatial variables by 
$h$ and the pressure by 
$\rho u_{\tau }^2$. A ‘
$+$’ superscript denotes spatial variables that have been normalised by the viscous length scale 
$\nu /u_{\tau }$.
2.2. Resolvent analysis
We begin by a Reynolds-decomposition of (2.1a), which leads to the following equations for the fluctuations:
\begin{equation} \frac{\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}' + \boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} + \boldsymbol{\nabla} p' - \frac{1}{Re_{\tau}}\boldsymbol{\nabla}^2 \boldsymbol{u}' = -\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}' + \overline{\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'} = \boldsymbol{f}', \end{equation}
where 
$(\bar {\cdot })$ and 
$(\cdot )'$ denote a time-average and fluctuation, respectively, and 
$\boldsymbol {U} = [U(z),0,0]^\textrm {T}$ is the time-averaged velocity field. Equation (2.2) is written such that all linear terms appear on the left-hand side while all nonlinear terms appear on the right-hand side. Equation (2.2) is then Laplace-transformed in time and Fourier-transformed in the homogeneous directions 
$x$ and 
$y$:
\begin{equation} \hat{\boldsymbol{u}}(k_x,k_y,s) = \frac{1}{(2 {\rm \pi})^3} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \boldsymbol{u}'(x,y,z,t) \exp({st-\textrm{i}k_xx-\textrm{i}k_yy})\,\textrm{d}x\,\textrm{d}y\,\textrm{d}t. \end{equation}
Upon integration of (2.3), we set 
$s=\textrm {i}\omega$ to consider the frequency response 
$\hat {\boldsymbol {u}}(\boldsymbol {k})$, where 
$(\hat {\cdot })$ denotes the Fourier-transformed coefficient and the wavenumber triplet 
$\boldsymbol {k} = (k_x,k_y,\omega )$ consists of streamwise wavenumber 
$k_x$, spanwise wavenumber 
$k_y$ and temporal frequency 
$\omega$. The equivalent wavelengths in the streamwise and spanwise directions are 
$\lambda _x = 2{\rm \pi} /k_x$ and 
$\lambda _y = 2{\rm \pi} /k_y$.
The equations are arranged into state-space form (Jovanović & Bamieh Reference Jovanović and Bamieh2005) after substituting (2.3) into (2.2),
\begin{gather} \textrm{i}\omega \hat{\boldsymbol{q}}(\boldsymbol{k}) = \boldsymbol{A}(k_x,k_y)\hat{\boldsymbol{q}}(\boldsymbol{k}) + \boldsymbol{B}(k_x,k_y)\hat{\boldsymbol{f}}(\boldsymbol{k}) , \end{gather}
\begin{gather}\hat{\boldsymbol{u}}(\boldsymbol{k}) = \boldsymbol{C}(k_x,k_y)\hat{\boldsymbol{q}}(\boldsymbol{k}) , \end{gather}
where 
$\hat {\boldsymbol {q}}$ consists of the wall-normal velocity and vorticity 
$\hat {\eta } = \textrm {i}k_y\hat {u} - \textrm {i}k_x\hat {v}$. The matrices 
$\boldsymbol {A}$, 
$\boldsymbol {B}$ and 
$\boldsymbol {C}$ represent the discretised forms of the linearised Navier–Stokes operator, the forcing operator and the output operator, respectively, and are defined in appendix A. These operators are independent of 
$\omega$ but are functions of the wavenumber pair 
$(k_x,k_y)$ under consideration. In the interest of readability, this dependence is omitted for the rest of the paper.
Once (2.4) is recast into input-output form, i.e.
\begin{equation} \hat{\boldsymbol{u}}(\boldsymbol{k}) = \boldsymbol{C}(\textrm{i}\omega \boldsymbol{I}- \boldsymbol{A})^{-1} \hat{\boldsymbol{f}}(\boldsymbol{k}) = \mathcal{H}(\boldsymbol{k}) \hat{\boldsymbol{f}}(\boldsymbol{k}), \end{equation}
a linear operator called the resolvent 
$\mathcal {H}(\boldsymbol {k})$ relates the input forcing 
$\skew4\hat {\boldsymbol {f}}(\boldsymbol {k})$ to the output velocity 
$\hat {\boldsymbol {u}}(\boldsymbol {k})$. Even if the nonlinear forcing is unknown, the resolvent operator can be characterised by the singular value decomposition
\begin{equation} \mathcal{H}(\boldsymbol{k}) = \hat{\boldsymbol{\varPsi}} (\boldsymbol{k}) \boldsymbol{\varSigma}(\boldsymbol{k}) \hat{\boldsymbol{\varPhi}}^*(\boldsymbol{k}), \end{equation}
where 
$\hat {\boldsymbol {\varPsi }}(\boldsymbol {k}) = [\hat {\boldsymbol {\psi }}_1(\boldsymbol {k}), \hat {\boldsymbol {\psi }}_2(\boldsymbol {k}), \ldots , \hat {\boldsymbol {\psi }}_p(\boldsymbol {k})]$ and 
$\hat {\boldsymbol {\varPhi }}(\boldsymbol {k}) = [\hat {\boldsymbol {\phi }}_1(\boldsymbol {k}),\hat {\boldsymbol {\phi }}_2(\boldsymbol {k}), \ldots , \hat {\boldsymbol {\phi }}_p(\boldsymbol {k})]$ are orthogonal basis functions for the velocity and nonlinear forcing, respectively. The diagonal matrix 
$\boldsymbol {\varSigma }(\boldsymbol {k})$ ranks the 
$p$th structure by its gain 
$\sigma _p(\boldsymbol {k})$ using an inner product that is proportional to its kinetic energy, i.e. 
$\left \langle \hat {\boldsymbol {\psi }},\hat {\boldsymbol {\psi }} \right \rangle = \int _{-h}^h \hat {\boldsymbol {\psi }}^* \boldsymbol {\cdot } \hat {\boldsymbol {\psi }}\,\textrm {d}z$. Consequently, the structure 
$\hat {\boldsymbol {\psi }}_1(\boldsymbol {k})$, referred to as the optimal or first resolvent mode, is the most amplified response by the linear dynamics contained in the operator. The true velocity field is the weighted sum of resolvent modes, i.e.
\begin{equation} \hat{\boldsymbol{u}}(\boldsymbol{k}) = \sum_{p = 1}^N \hat{\boldsymbol{\psi}}_p(\boldsymbol{k}) \sigma_p(\boldsymbol{k}) \chi_p(\boldsymbol{k}), \end{equation}
where 
$\chi _p(\boldsymbol {k})$ is the projection of 
$\hat {\boldsymbol {\phi }}_p(\boldsymbol {k})$ onto 
$\skew4\hat {\boldsymbol {f}}(\boldsymbol {k})$.
2.3. Energy balance
We now derive the energy balance that must be satisfied by the velocity field and individual resolvent modes. Equation (2.2) is rewritten in index notation
\begin{equation} \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} = - u'_j\frac{\partial u_i'}{\partial x_j} + \overline{u'_j\frac{\partial u_i'}{\partial x_j}} . \end{equation}
Similar to (2.2), all nonlinear terms appear on the right-hand side although they are not treated as an unknown forcing. The indices 
$i,j =1,2,3$ and 
$U_i = (U(z),0,0)$ is the mean velocity, which is a function of the wall-normal direction only. It can be noted, therefore, that 
$U_1 = U$, 
$U_j = 0$ if 
$j=2,3$ and 
$\partial U_i/\partial x_j \neq 0$ for 
$i = 1$ and 
$j = 3$ only. The kinetic energy of the full system is characterised by the inner product between (2.8) and 
$u_i'$ integrated over the volume 
$V$:
\begin{gather} \underbrace{ \frac{1}{2} \int_V \frac{\partial u_i^{'2}}{\partial t} \,\textrm{d}V }_{\skew4\dot{E}(t)} = \underbrace{- \int_V u'_i u_j' \frac{\partial U_i}{\partial x_j} \,\textrm{d}V }_{P(t)} \underbrace{ - \frac{1}{Re} \int_V \frac{\partial u_i'}{\partial x_j} \frac{\partial u_i'}{\partial x_j} \,\textrm{d}V}_{D(t)}, \end{gather}
\begin{gather}\int_V u_i'\,f_i'\,\textrm{d}V = \int_V u_i'u_j'\frac{\partial u_i'}{\partial x_j}\,\textrm{d}V = 0. \end{gather}
Equation (2.9) is the Reynolds–Orr equation (Schmid & Henningson Reference Schmid and Henningson2001), where the evolution of kinetic energy in the system is a balance between production and dissipation, which must be negative. Due to the conservative nature of the nonlinear terms, their contribution to the Reynolds–Orr equation sums to zero when integrated over the volume as expressed in (2.10). For a statistically stationary flow, a time average of (2.9) implies that production balances dissipation since 
$\overline {\textrm {d}E/\textrm {d}t} = 0$.
The kinetic energy for a specific spatial scale is obtained after multiplying (2.8) by 
$u_i^{'*}$ and Fourier-transforming in 
$x$ and 
$y$. The result is integrated over the wall-normal direction and time-averaged to arrive at the spectral turbulent kinetic energy (TKE) equation:
\begin{align} \overline{\frac{\partial \skew4\hat{E}(k_x,k_y)}{\partial t}} &= \underbrace{- \overline{\left\langle \frac{\textrm{d}U}{\textrm{d}z}\hat{u}(k_x,k_y),\hat{w}(k_x,k_y)\right\rangle}}_{\hat{P}(k_x,k_y)} \underbrace{-\frac{1}{Re} \overline{ \left\langle \frac{\partial \hat{u}_i(k_x,k_y)}{\partial x_j}, \frac{\partial \hat{u}_i(k_x,k_y)}{\partial x_j} \right\rangle}}_{\hat{D}(k_x,k_y)} \nonumber\\ &\quad \underbrace{ - \overline{ \left\langle \hat{u}_i(k_x,k_y), \frac{\partial}{\partial x_j} \widehat{u_iu_j}(k_x,k_y) \right\rangle } }_{\hat{N}(k_x,k_y)} = 0. \end{align}
The pressure terms vanish in (2.11) after integrating over the channel height (Aubry et al. Reference Aubry, Holmes, Lumley and Stone1988). Following Muralidhar et al. (Reference Muralidhar, Podvin, Mathelin and Fraigneau2019), we consider the real part of (2.11), which consists of three terms: production, viscous dissipation and nonlinear transfer. In general, production 
$\hat {P}$ is positive for a given scale as perturbations extract energy from the mean flow. Viscous dissipation 
$\hat {D}$, on the other hand, is guaranteed to be real and negative according to (2.11) as it is the mechanism through which kinetic energy is removed from the system and converted into heat. Nonlinear transfer 
$\hat {N}$ may be positive or negative depending on the scale selected. If 
$\hat {P} > \hat {D}$, for example, then 
$\hat {N} < 0$ in order to achieve a balance. In a similar fashion, if 
$\hat {P} < \hat {D}$, then 
$\hat {N} > 0$. The integral of 
$\hat {N}$ over all 
$k_x$ and 
$k_y$, nevertheless, is zero as stated in (2.10).
To obtain the energy balance for resolvent modes, which are defined for a wavenumber triplet 
$\boldsymbol {k}$, (2.8) is Fourier-transformed in 
$x$, 
$y$ and 
$t$ and multiplied by 
$\hat {u}_i^{'*}$. The result is integrated over the wall-normal direction:
\begin{equation} -\left\langle \frac{\textrm{d}U}{\textrm{d}z} \hat{u} (\boldsymbol{k}), \hat{w} (\boldsymbol{k}) \right\rangle -\frac{1}{Re} \left\langle \frac{\partial \hat{u}_i(\boldsymbol{k})}{\partial x_j}, \frac{\partial \hat{u}_i(\boldsymbol{k})}{\partial x_j} \right\rangle - \left\langle \hat{u}_i(\boldsymbol{k}),\hat{f}_i(\boldsymbol{k}) \right\rangle = 0. \end{equation}
In this form, the nonlinear forcing 
$\skew4\hat {\boldsymbol {f}}(\boldsymbol {k})$ appears explicitly in the energy balance. We can now express the velocity field in terms of resolvent modes. In the special case where 
$\skew4\hat {\boldsymbol {f}}(\boldsymbol {k})$ is white noise, the velocity field can be written as
\begin{equation} \hat{\boldsymbol{u}}(\boldsymbol{k}) = \sum_{p = 1}^N \hat{\boldsymbol{\psi}}_p(\boldsymbol{k}) \sigma_p(\boldsymbol{k}). \end{equation}Substituting (2.13) into (2.12) yields
\begin{align} & \sum_p \sigma_p(\boldsymbol{k}) \left( \left\langle \frac{\textrm{d}U}{\textrm{d}z} \hat{\boldsymbol{\psi}}_p^{i=1}(\boldsymbol{k}), \hat{\boldsymbol{\psi}}_p^{j=3} (\boldsymbol{k}) \right\rangle + \frac{1}{Re} \left\langle \frac{\partial \hat{\boldsymbol{\psi}}_p^{i}(\boldsymbol{k})} {\partial \boldsymbol{x}_j}, \frac{\partial \hat{\boldsymbol{\psi}}_p^i (\boldsymbol{k})}{\partial \boldsymbol{x}_j} \right\rangle \right) \nonumber\\ &\quad +\sum_p \left\langle \hat{\boldsymbol{\psi}}_p(\boldsymbol{k}),\hat{\boldsymbol{\phi}}_p(\boldsymbol{k})\right\rangle = 0. \end{align}
Each term in the sum can be decoupled since the basis functions 
$\hat {\boldsymbol {\psi }}_p$ are orthogonal. This means that production, dissipation and nonlinear transfer must be balanced across each resolvent mode. If 
$\sigma _1(\boldsymbol {k}) \gg \sigma _2(\boldsymbol {k})$, this would imply that the sum over all 
$p$ is dominated by the first resolvent mode, or
\begin{align} \sigma_1(\boldsymbol{k}) \left( \left\langle \frac{\textrm{d}U}{\textrm{d}z} \hat{\boldsymbol{\psi}}_1^{i=1}(\boldsymbol{k}), \hat{\boldsymbol{\psi}}_1^{j=3} (\boldsymbol{k}) \right\rangle + \frac{1}{Re} \left\langle \frac{\partial \hat{\boldsymbol{\psi}}_1^{i}(\boldsymbol{k})}{\partial \boldsymbol{x}_j}, \frac{\partial \hat{\boldsymbol{\psi}}_1^i (\boldsymbol{k})}{\partial \boldsymbol{x}_j} \right\rangle \right) +\left\langle \hat{\boldsymbol{\psi}}_1(\boldsymbol{k}),\hat{\boldsymbol{\phi}}_1(\boldsymbol{k})\right\rangle = 0. \end{align}
The bulk of production, dissipation and nonlinear transfer for a particular scale 
$\boldsymbol {k}$, therefore, could also be accounted for by the first resolvent mode. It should be noted that (2.15) assumes the nonlinear forcing is white in space even though in the true flow it is not. Nevertheless, (2.15) is a reasonable assumption as long as 
$\sigma _1(\boldsymbol {k}) \gg \sigma _2(\boldsymbol {k})$ and the nonlinear forcing does not preferentially project onto a suboptimal resolvent forcing mode, i.e. 
$\left \langle \hat {\boldsymbol {\phi }}_1(\boldsymbol {k}),\skew4\hat {\boldsymbol {f}}(\boldsymbol {k}) \right \rangle \approx \left \langle \hat {\boldsymbol {\phi }}_{p>1}(\boldsymbol {k}),\skew4\hat {\boldsymbol {f}}(\boldsymbol {k}) \right \rangle$ (Beneddine et al. Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016).
2.4. Eddy viscosity model
If the nonlinear forcing is not white noise, then (2.14) is not applicable since it does not take into account the complex amplitude of each mode. One method to model the nonlinear forcing is to add an eddy viscosity to the linearised equations after performing a triple decomposition of the velocity field 
$\tilde {\boldsymbol {u}}$ into a mean component 
$\boldsymbol {U}$, coherent motions 
$\boldsymbol {u}$ and incoherent turbulent fluctuations 
$\boldsymbol {u}'$ (Reynolds & Hussain Reference Reynolds and Hussain1972). The equations governing the coherent velocity and pressure are
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} + \boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ \nu_T ( \boldsymbol{\nabla} {\boldsymbol{u}} + \boldsymbol{\nabla} {\boldsymbol{u}}^\textrm{T}) \right] = \boldsymbol{d}, \end{equation}
where 
$\nu _T(z)$ is the total effective viscosity, and 
$\boldsymbol {d} = -\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u} + \overline {\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}}$ is the forcing term. It should be noted that 
$\boldsymbol {d}$ is different from 
$\boldsymbol {f}'$ in (2.2). Following Reynolds & Tiederman (Reference Reynolds and Tiederman1967) and Hwang & Cossu (Reference Hwang and Cossu2010), we use the Cess (Reference Cess1958) eddy viscosity profile of the form
\begin{equation} \nu_T(z) = \frac{\nu}{2}\left(1 + \left[ \frac{\kappa}{3}(1-z^2)(1+2z^2)\left(1-\exp({|z-1| {Re_{\tau}}/{A}})\right) \right]^2 \right)^{1/2} + \frac{\nu}{2}, \end{equation}
where 
$\kappa = 0.426$ and 
$A = 25.4$ are chosen based on a least-squares fit to experimentally obtained mean velocity profiles at 
$Re_{\tau } = 2000$ (del Álamo & Jiménez Reference del Álamo and Jiménez2006). Although the Reynolds numbers studied here are lower, we have verified that the results are not sensitive to the values of these constants.
Fourier-transforming (2.16) in time and the homogeneous directions and rearranging it into the following input-output form yields
\begin{equation} \hat{\boldsymbol{u}}(\boldsymbol{k}) = \mathcal{H}^e(\boldsymbol{k})\hat{\boldsymbol{d}}(\boldsymbol{k}), \end{equation}
where 
$\mathcal {H}^e(\boldsymbol {k})$ is a new linear operator that relates the forcing 
$\hat {\boldsymbol {d}} (\boldsymbol {k})$ to the velocity field 
$\hat {\boldsymbol {u}}(\boldsymbol {k})$. Similar to the resolvent operator, we can analyse structures that are preferentially amplified by performing a singular value decomposition
\begin{equation} \mathcal{H}^e(\boldsymbol{k}) = \hat{\boldsymbol{\varPsi}}^e (\boldsymbol{k}) \boldsymbol{\varSigma}^e (\boldsymbol{k}) \hat{\boldsymbol{\varPhi}}^{*,e}(\boldsymbol{k}), \end{equation}
although the individual modes 
$\hat {\boldsymbol {\psi }}_p^e(\boldsymbol {k})$ do not satisfy the energy balance in (2.14). Instead, the addition of eddy viscosity in (2.16) introduces two terms, the first of which is
\begin{equation} \hat{V}(k_x,k_y) = -\overline{ \left\langle(\nu_T(z)-\nu) \frac{\partial \hat{u}_i(k_x,k_y)}{\partial x_j}, \frac{\partial \hat{u}_i(k_x,k_y)}{\partial x_j} \right\rangle}, \end{equation}
where the kinematic viscosity 
$\nu$ has been subtracted in order to remove the contribution of viscous dissipation 
$\hat {D}(k_x,k_y)$. The remainder 
$\hat {V}(k_x,k_y)$ therefore represents the additional dissipation provided by the wall-normal varying portion of 
$\nu _T$. Similar to 
$\hat {D}(k_x,k_y)$, this term is real and negative, signifying that it removes energy. The second term is related to the wall-normal gradient of 
$\nu _T$
\begin{equation} \hat{G}(k_x,k_y) = - \overline{ \left\langle \frac{d \nu_T}{dz} \hat{u}_i(k_x,k_y), \frac{\partial \hat{u}_i(k_x,k_y)}{\partial z} + \frac{\partial \hat{w}(k_x,k_y)}{\partial x_i} \right\rangle} . \end{equation}
Unlike 
$\hat {V}(k_x,k_y)$, the sign of 
$\hat {G}(k_x,k_y)$ cannot be determined a priori.
The combined effect of 
$\hat {V}(k_x,k_y)$ and 
$\hat {G}(k_x,k_y)$ is referred to as eddy dissipation 
$\widehat {Edd}(k_x,k_y)$, i.e.,
\begin{equation} \widehat{Edd}(k_x,k_y) = \hat{V}(k_x,k_y) + \hat{G}(k_x,k_y). \end{equation}
Eddy dissipation is computed in §§ 4 and 5 using the true velocity field to determine its accuracy in modelling the effect of nonlinear transfer in (2.11). If 
$\widehat {Edd}(k_x,k_y) \approx \hat {N}(k_x,k_y)$, then it is expected that eddy viscosity will lead to an improvement in the structures predicted by resolvent analysis.
3. Flow descriptions
In this section, we describe the two flows that are analysed from an energy transfer perspective. These are the P4U ECS computed by Park & Graham (Reference Park and Graham2015) and turbulent channel flow in the minimal unit (Jiménez & Moin Reference Jiménez and Moin1991), which are discussed in §§ 3.1 and 3.2, respectively.
3.1. P4U ECS
The P4U ECS, henceforth referred to as P4U, is a nonlinear travelling wave with a friction Reynolds number of 
$Re_{\tau } = 85$ and fixed wave speed of 
$c^+ = 14.2$. As seen in table 1, P4U is solved in a computational domain with 24 equally spaced grid points in the streamwise and spanwise directions, which have lengths of 
${\rm \pi}$ and 
${\rm \pi} /2$, respectively. There are 81 points in the wall-normal direction on a Chebyshev grid. The spatial structure of P4U is in the form of low-speed streaks, which are wavy in the streamwise direction, straddled by counter-rotating vortices. As mentioned by Park & Graham (Reference Park and Graham2015), its structure is qualitatively similar to near-wall turbulence, and its mean velocity profile closely resembles a standard turbulent mean. This is seen more clearly in figure 1(a) where the mean profile for P4U in blue is compared to the Cess model in red at 
$Re_{\tau } = 85$. Despite good overall agreement, the P4U ECS profile has a more wavy nature since the structure has a single convection velocity. Additional similarities between P4U and standard turbulence, such as the agreement between the Reynolds shear stress profiles, can be found in Park & Graham (Reference Park and Graham2015).
Table 1. Relevant parameters for the flows under consideration.
Figure 1. (a) Mean velocity profiles for P4U (blue) and Cess model (red) at 
$Re_{\tau } = 85$. (b) Mean velocity profiles for the minimal channel (blue) and the DNS of Lee & Moser (Reference Lee and Moser2015) (black).
Even though the simulation size is small, there are still many wavenumber pairs which may participate in the transfer of energy. We begin by computing the kinetic energy of each wavenumber pair,
\begin{equation} \skew4\hat{E}(k_x,k_y) = \tfrac{1}{2}\overline{\left\langle\hat{\boldsymbol{u}}(k_x,k_y), \hat{\boldsymbol{u}}(k_x,k_y)\right\rangle}, \end{equation}and plot the most energetic pairs in figure 2(a). The area and colour intensity of the square marker at the centre of each tile are directly proportional to the kinetic energy. The most energetic scale is streamwise-constant with a spanwise width of approximately 100 wall units, which matches the near-wall streak spacing of Smith & Metzler (Reference Smith and Metzler1983). Furthermore, most of the kinetic energy is concentrated in structures with small streamwise wavenumbers.
Figure 2. The kinetic energy 
$\hat {E}(k_x,k_y)$ of the most energetic wavenumber pairs in (a) P4U and the (b) minimal channel. The area and colour intensity of the square marker at the centre of each tile are directly proportional to the kinetic energy.
3.2. Minimal channel
The minimal channel flow is computed for 
$Re_{\tau } = 180$ using an unstructured finite difference solver (see Chung, Monty & Ooi (Reference Chung, Monty and Ooi2014), for details) on a domain with dimensions 
${\rm \pi} \times {\rm \pi}/4 \times 2h$ in the streamwise, spanwise and wall-normal directions. There are 96 and 48 equally spaced points in the streamwise and spanwise directions, respectively, and 128 points in the wall-normal direction on a Chebyshev grid. The mean profile for the minimal channel in figure 1(b) is smoother than the P4U mean profile since there exists a distribution of wave speeds which range between 
$0 < c^+ < 19.4$ as seen in table 1. The minimal channel mean profile is in good agreement with the mean profile of Lee & Moser (Reference Lee and Moser2015) for most areas of the flow other than the wake region, where the minimal channel mean profile overshoots the one from Lee & Moser (Reference Lee and Moser2015). This phenomenon has been observed by Jiménez & Moin (Reference Jiménez and Moin1991) and stems from the fact that the minimal domain is too small to accommodate the largest structures which reside in the outer region. Despite this disagreement, there is no impact on near-wall turbulence in the buffer and viscous regions where the bulk of energy resides (Jiménez & Moin Reference Jiménez and Moin1991; Jiménez & Pinelli Reference Jiménez and Pinelli1999).
Similar to P4U, the kinetic energy for the most energetic wavenumber pairs is plotted in figure 2(b). Although there are more energetic scales in the minimal channel since the friction Reynolds number is higher, the relative distribution of energy among the scales is quite similar to P4U. The most energetic scale is also streamwise-constant with a spanwise width of approximately 100 wall units. This supports the notion that P4U is a simple model for turbulent channel flow at very low Reynolds number. Therefore, to facilitate visualisation later in the paper, we choose to plot only those wavenumber pairs that appear in figure 2 although the energy balance will be computed across all of them.
4. Results: P4U ECS
In this section, we analyse energy transfer for P4U. We begin with a comparison of production, dissipation and nonlinear transfer across the most energetic scales in § 4.1. These results are compared to the resolvent predictions in § 4.2. Finally, the additional dissipation introduced by eddy viscosity is quantified for each scale and compared to nonlinear transfer in § 4.3.
4.1. Energy balance
Production, dissipation and nonlinear transfer are computed for P4U and are illustrated in figure 3 for the subset of wavenumber pairs discussed in the previous section. A square marker appears at the centre of each tile. Both colour intensity and area of the square marker indicate each term's magnitude. The colours red and blue denote positive and negative quantities, respectively. In order to satisfy (2.11), the sum across tiles which appear in the same position in each of the three figure panels must be zero. Additionally, the sum over all tiles in figure 3(c) is approximately zero since the nonlinear terms are conservative when summed over all scales (this sum would be exactly zero if all wavenumber pairs were displayed in the figure).
Figure 3. Contributions of (a) production, (b) dissipation and (c) nonlinear transfer to the energy balance of each Fourier mode for P4U. Both colour intensity and area of the square marker indicate each term's magnitude. The colours red and blue denote positive and negative quantities, respectively.
Half of the production terms in figure 3(a) are positive, with the largest energy-producing modes being the streamwise-constant modes. The maximum production occurs at 
$(k_x,k_y) = (0,4)$, which is also the most energetic mode in the flow (see figure 2a). Production is negative for some scales. Of particular note is that production is negative for all of the spanwise-constant modes. Even though 
$\hat {P} \approx 0$ for most of these spanwise-constant modes, the same cannot be said for 
$(2,0)$, for which the production is negative and of large amplitude. In fact, its magnitude is comparable to that of 
$(0,8)$ even though it is less energetic – i.e. 
$|\hat {P}(2,0)| \approx |\hat {P}(0,8)|$ even though 
$\hat {E}(2,0) < \hat {E}(0,8)$. As expected, all dissipation terms in figure 3(b) are negative.
The nonlinear transfer in figure 3(c) contains both positive and negative terms as the sum over all scales must equal zero. Consistent with the turbulent cascade, most values are positive, indicating that they are receiving energy from nonlinear transfer. The most notable exception is the 
$(0,4)$ mode, which must redistribute energy to other scales since dissipation offsets less than half of production. The additional scales that lose energy due to nonlinear transfer all have low streamwise wavenumbers. The 
$(0,8)$ mode is one that receives energy from nonlinear transfer since 
$\hat {D}(0,8) > \hat {P}(0,8)$. Perhaps surprisingly, the spanwise-constant modes receive a considerable share of the nonlinearly-transferred energy. In particular the 
$(2,0)$ mode receives more energy than any other mode. The 
$(4,0)$ and 
$(6,0)$ modes also receive rather than donate energy. Therefore, in addition to a cascade of energy from large scales to small scales, there is also a significant transfer from scales that are streamwise-constant to scales that are spanwise-constant. Indeed, the 
$(2,0)$ mode (the largest recipient) is in fact larger in scale than the 
$(0,4)$ mode (the largest donor). Thus, in addition to a cascade, there also exists a transfer to scales of a similar scale but with a different orientation of their wavenumber vector.
4.2. Resolvent predictions
Having considered the true energy balance from (2.11), we now focus on its counterpart for the first resolvent mode in (2.15). To do so, it is necessary to set 
$\omega = c^+k_x$ since the wave speed is fixed at 
$c^+ = 14.2$. The mean profile is computed directly from P4U although we use the Cess eddy viscosity profile in the following section to calculate eddy dissipation. Figure 4 illustrates the production, dissipation and nonlinear transfer in a manner analogous to that of figure 3. The resolvent prediction for production in figure 4(a) is positive for every scale, and the largest value occurs when 
$(k_x,k_y) = (0,4)$. The predictions for the largest scales are similar to the true values of production in figure 3(a) and reflect the resolvent operator's ability to identify linear amplification mechanisms.
Figure 4. Contributions of (a) production, (b) dissipation and (c) nonlinear transfer from the first P4U resolvent mode.
The dissipation and nonlinear transfer from the first resolvent mode in figures 4(b) and 4(c), respectively, are less similar to the true values in figures 3(b) and 4(c). For all scales, the dissipation is nearly equal and opposite to production, resulting in very small values for nonlinear transfer in figure 4(c). Most tiles appear empty since 
$\hat {N} \ll \hat {P},\hat {D}$ for all wavenumber pairs. A similar phenomenon is observed by Jin, Symon & Illingworth (Reference Jin, Symon and Illingworth2020) for the first resolvent mode in low Reynolds number cylinder flow. It can be concluded that suboptimal resolvent modes are necessary to correctly model nonlinear transfer.
4.3. Eddy dissipation
As discussed in § 2.4, one way to model nonlinear transfer is through the use of an eddy viscosity. Figure 5(a) presents the eddy dissipation from (2.20), which is negative for all wavenumber pairs considered. Unlike nonlinear transfer, therefore, eddy dissipation is not conservative and contributes net energy loss to every scale. Ideally the eddy dissipation would resemble 
$\hat {N}(k_x,k_y)$ in figure 3(c), so its error 
$\epsilon (k_x,k_y)$ with respect to nonlinear transfer is computed in figure 5(b) using the expression
\begin{equation} \epsilon(k_x,k_y) = \frac{\widehat{Edd}(k_x,k_y) - \hat{N}(k_x,k_y)}{|\hat{N}(k_x,k_y)|}. \end{equation}
The error for all wavenumber pairs exceeds 1 other than 
$(0,4)$ where 
$\epsilon \approx 0.28$. The size of the square marker in figure 5(b) reflects the magnitude of the error, and the smallest marker coincides with the tile belonging to 
$(0,4)$. The fact that 
$\epsilon$ is lowest for this scale indicates that the eddy viscosity is most effective for highly amplified linear mechanisms. In other words, the eddy viscosity works best for scales where viscous dissipation is not sufficient to balance production. This is as expected since the motivation for including an eddy viscosity is to model turbulent transfer of energy from the large scales.
Figure 5. (a) Nonlinear transfer modelled by eddy viscosity for P4U and (b) its error compared to the true nonlinear transfer in figure 3(c).
As noted in (2.22), the eddy dissipation consists of two terms. The first, 
$\hat {V}(k_x,k_y)$, stems from 
$\nu _T(z)$ and can be interpreted as a wall-normal-varying effective Reynolds number. The second, 
$\hat {G}(k_x,k_y)$, arises due to the gradient of the eddy viscosity profile 
$\nu _T'(z)$. In figure 6, the contribution of both terms illustrates that the majority of eddy dissipation is due to 
$\hat {V}(k_x,k_y)$ in figure 6(a). As expected, 
$\hat {V}(k_x,k_y)$ is negative for all scales whereas 
$\hat {G}(k_x,k_y)$ in figure 6(b) is mostly positive albeit small. Therefore, the success behind eddy viscosity lies in its locally modifying the Reynolds number via 
$\hat {V}(k_x,k_y)$.
Figure 6. Contributions from (a) 
$\hat {V}(k_x,k_y)$ and (b) 
$\hat {G}(k_x,k_y)$ to the eddy dissipation 
$\widehat {Edd}(k_x,k_y)$ for P4U.
5. Results: minimal channel
Having analysed the energy transfer for the P4U ECS, this section examines the same quantities for the minimal channel. Since each wavenumber pair has a distribution of temporal frequencies, all terms in the energy balance are time-averaged.
5.1. DNS and resolvent energy balances
Production, dissipation and nonlinear transfer for the minimal channel are illustrated in figure 7. For almost all wavenumber pairs shown, production is positive as seen in figure 7(a) with the maximum occurring for 
$(0,8)$. The only scales where production is negative are spanwise-constant, i.e. 
$k_y = 0$. The dissipation in figure 7(b) is negative for all scales, as expected. Even though the largest dissipation occurs for 
$(0,8)$, its value is comparable to that for other wavenumbers.
Figure 7. Contributions of (a) production, (b) dissipation and (c) nonlinear transfer to the energy balance of each Fourier mode for the minimal channel.
The nonlinear transfer in figure 7(c) illustrates that the surfeit of energy not dissipated by viscosity from 
$(0,8)$ is redistributed to other scales. Moreover, all scales which lose energy due to nonlinear transfer are clustered around low streamwise and spanwise wavenumbers. This is consistent with the turbulence cascade in that energy from the large scales trickles down to smaller scales which are more effective at dissipating energy. All but five scales in figure 7(c) receive energy from nonlinear transfer with the largest amounts going to spanwise-constant structures. For laminar flows, these Tollmien–Schlichting-type waves (Tollmien Reference Tollmien1930; Schlichting Reference Schlichting1933) are the first to become unstable. In the minimal channel, however, they play a damping role.
The energy balance for the first resolvent mode is presented in figure 8. Similar to P4U, we use the DNS mean in resolvent analysis and the Cess eddy viscosity profile to calculate eddy dissipation. Since each wavenumber pair has a distribution of energetic temporal frequencies, we compute the singular values across a discretisation of 
$\omega$ and choose the 
$\omega$ that results in the largest amplification. For example, 
$\omega = 0$ leads to the largest amplification for 
$(0,8)$ and 
$\omega = 26$ for 
$(2,8)$. It can be remarked that viscous dissipation in figure 8(b) is sufficient to completely counteract production for the majority of scales considered. Since the sum of all three terms must be zero for each resolvent mode, it follows that nonlinear transfer is negligible for nearly all scales, as seen in figure 8(c).
Figure 8. Contributions of (a) production, (b) dissipation and (c) nonlinear transfer from the first resolvent mode in the case of the minimal channel.
To identify for which scales eddy viscosity can model nonlinear transfer, the eddy dissipation is computed and displayed in figure 9(a). As expected, it is negative for all scales even though nonlinear transfer tends to be positive outside the cluster around 
$(0,8)$. The error 
$\epsilon$, as defined in (4.1), is thus large for the majority of scales as seen in figure 9(b). The only scale where 
$\epsilon < 1$ is 
$(0,8)$. Although 
$\epsilon > 1$ for every other scale, those where nonlinear transfer is negative, such as 
$(0,16)$ or 
$(2,8)$, have lower values of 
$\epsilon$ than scales where nonlinear transfer is positive.
Figure 9. (a) Nonlinear transfer modelled by eddy viscosity for the minimal channel and (b) its error compared to the true nonlinear transfer in figure 7(c).
In figure 10, the eddy dissipation is split into contributions from 
$\hat {V}(k_x,k_y)$ and 
$\hat {G}(k_x,k_y)$. Similar to P4U, the dominant term in figure 10(a) is 
$\hat {V}(k_x,k_y)$, and it is always negative. 
$\hat {G}(k_x,k_y)$ in figure 10(b), on the other hand, is negative primarily for 
$k_x = 0$ modes although it is positive for most scales. Even though the role of 
$\hat {G}(k_x,k_y)$ seems negligible for the minimal channel, it should be noted that the Reynolds number is only a factor of two greater than that of P4U. It is possible that this term is important at higher Reynolds numbers where 
$\nu _T'(z)$ becomes large near the wall.
Figure 10. Contributions from (a) 
$\hat {V}(k_x,k_y)$ and (b) 
$\hat {G}(k_x,k_y)$ to the eddy dissipation 
$\widehat {Edd}(k_x,k_y)$ for the minimal channel.
5.2. Comparison of P4U and minimal channel
The energy transfer processes in P4U and the minimal channel are similar. Production is positive for the majority of scales, and its maximum occurs for 
$(\lambda _x^+,\lambda _y^+) \approx (\infty ,100)$, which corresponds to 
$(k_x,k_y) = (0,4)$ in P4U and 
$(k_x,k_y) = (0,8)$ in the minimal channel. The production for spanwise-constant scales, on the other hand, is mostly negative. Dissipation is always negative, but it is insufficiently large to counterbalance production for the largest streamwise-constant scales. Nonlinear transfer contains both positive and negative terms as the sum over all scales must be zero according to (2.10). Energy is primarily removed from the largest structures and redistributed to smaller ones.
A polarisation effect can be observed where nonlinear transfer reallocates energy from streamwise-constant to spanwise-constant modes. This transverse cascade is a feature of the self-sustaining process (Hamilton et al. Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997) and has been observed by Mamatsashvili et al. (Reference Mamatsashvili, Gogichaishvili, Chagelishvili and Horton2014, Reference Mamatsashvili, Khujadze, Chagelishvili, Dong, Jiménez and Foysi2016) in magnetohydrodynamic turbulence and homogeneous shear turbulence, respectively. For turbulent channel flow, it occurs since production is large for structures where 
$k_x < k_y$ as they most efficiently leverage the lift-up mechanism. Dissipation, on the other hand, is proportional to the size of a structure and has no preferred direction. Finally, our choice to integrate over the wall-normal domain and analyse energy transfer in 
$(k_x,k_y)$ space permits quantification of energy loss due to eddy dissipation. For both flows, the only scale where eddy dissipation can quantitatively predict energy loss due to nonlinear transfer is the most energetic scale 
$(\lambda _x^+,\lambda _y^+) \approx (\infty ,100)$.
The only notable difference between the two flows is that nonlinear transfer has a clearer pattern for the minimal channel. Scales where this term is negative are localised in one cluster of large structures. For P4U, there is more scatter primarily due to the 
$(0,8)$ mode, which is able to dissipate more energy than it produces.
6. Non-normality in the energy balance
The results in §§ 4 and 5 indicate that the optimal resolvent mode does not accurately account for energy transfer between scales. As we will explain in § 6.1, the root of this discrepancy is non-normality induced by the mean shear. In § 6.2, we show that the influence of non-normality can be weakened by eddy viscosity. We analyse the efficiency of two sets of modes in reconstructing the energy balance in § 6.3. The first set is resolvent modes (see (2.6)) obtained from 
$\mathcal {H}(\boldsymbol {k})$, which has no eddy viscosity. The second set is eddy modes (see (2.19)) obtained from 
$\mathcal {H}_e(\boldsymbol {k})$, which does have eddy viscosity. In § 6.4, it is demonstrated that the eddy viscosity is most effective for high aspect ratio modes where the influence of non-normality is most pronounced. Finally, we discuss implications for resolvent modelling in § 6.5.
6.1. Competition between production and nonlinear transfer
In this section, we demonstrate that non-normality leads to a competition between the production and nonlinear-transfer terms in the energy balance. To simplify the discussion, we will only consider the first resolvent mode although similar arguments can be made for suboptimal modes. Equation (2.15) is rearranged into the following form:
\begin{equation} \underbrace{ \sigma_1(\boldsymbol{k}) \left\langle \frac{\textrm{d}U}{\textrm{d}z} \hat{\boldsymbol{\psi}}_1^{i=1}(\boldsymbol{k}), \hat{\boldsymbol{\psi}}_1^{j=3} (\boldsymbol{k}) \right\rangle}_{\hat{P}(\boldsymbol{k})} + \underbrace{\sigma_1(\boldsymbol{k}) \frac{1}{Re} \left\langle \frac{\partial \hat{\boldsymbol{\psi}}_1^{i}(\boldsymbol{k})}{\partial \boldsymbol{x}_j}, \frac{\partial \hat{\boldsymbol{\psi}}_1^i (\boldsymbol{k})}{\partial \boldsymbol{x}_j} \right\rangle}_{\hat{D}(\boldsymbol{k})} + \underbrace{\alpha (\boldsymbol{k} )}_{\hat{N}(\boldsymbol{k})}= 0, \end{equation}
where 
$\alpha (\boldsymbol {k}) = \left \langle \hat {\boldsymbol {\psi }}_1(\boldsymbol {k}), \hat {\boldsymbol {\phi }}_1(\boldsymbol {k}) \right \rangle$. It should be noted from (6.1) that if 
$\sigma _1(\boldsymbol {k}) \gg \alpha (\boldsymbol {k})$, then nonlinear transfer is small in comparison to production. This is consistent with the results in figures 4 and 8 where production was balanced primarily by dissipation. The first singular value 
$\sigma _1 (\boldsymbol {k})$ can be related to 
$\alpha (\boldsymbol {k})$ through
\begin{equation} \sigma_1(\boldsymbol{k}) \approx \underbrace{\frac{1}{(\textrm{i}\omega - \lambda_{min}(\boldsymbol{k}))}}_{C_1} \cdot \underbrace{\frac{1}{\left\langle \hat{\boldsymbol{\psi}}_1(\boldsymbol{k}),\hat{\boldsymbol{\phi}}_1(\boldsymbol{k}) \right\rangle}}_{1/\alpha(\boldsymbol{k})} . \end{equation}
We refer the reader to Symon et al. (Reference Symon, Rosenberg, Dawson and McKeon2018) for a derivation of (6.2) and note that it is applicable for resonant mechanisms such as 
$k_x = 0$ modes in channel flow or the shedding mode in cylinder flow. In (6.2), the first singular value is rewritten as the product of two terms. The first is 
$C_1$, the inverse distance between the imaginary axis and the least stable eigenvalue 
$\lambda _{min}(\boldsymbol {k})$ of the linearised Navier–Stokes operator 
$\boldsymbol {A}$. The second is a metric of non-normality (Chomaz Reference Chomaz2005), which is inversely proportional to 
$\alpha (\boldsymbol {k})$.
Equation (6.2) is substituted into (6.1)
\begin{equation} \underbrace{\frac{C_1}{\alpha (\boldsymbol{k})} \left\langle \frac{\textrm{d}U}{\textrm{d}z} \hat{\boldsymbol{\psi}}_1^{i=1}(\boldsymbol{k}), \hat{\boldsymbol{\psi}}_1^{j=3} (\boldsymbol{k}) \right\rangle}_{\hat{P}(\boldsymbol{k})} + \underbrace{\frac{C_1}{\alpha (\boldsymbol{k}) Re} \left\langle \frac{\partial \hat{\boldsymbol{\psi}}_1^{i}(\boldsymbol{k})}{\partial \boldsymbol{x}_j}, \frac{\partial \hat{\boldsymbol{\psi}}_1^i (\boldsymbol{k})}{\partial \boldsymbol{x}_j} \right\rangle}_{\hat{D}(\boldsymbol{k})} + \underbrace{\alpha (\boldsymbol{k} )}_{\hat{N}(\boldsymbol{k})}= 0, \end{equation}
so that each term in the energy balance is related to 
$\alpha (\boldsymbol {k})$. If we consider 
$\hat {P}(\boldsymbol {k})$ and 
$\hat {N}(\boldsymbol {k})$ in (6.3) relative to 
$\alpha (\boldsymbol {k})$, then
\begin{equation} \hat{P}(\boldsymbol{k}) \sim \frac{1}{\alpha (\boldsymbol{k})},\quad \hat{N}(\boldsymbol{k}) \sim \alpha (\boldsymbol{k}),\quad \frac{\hat{P}(\boldsymbol{k})}{\hat{N}(\boldsymbol{k})} \sim \frac{1}{\alpha (\boldsymbol{k})^2}. \end{equation}
Equations (6.4a–c) illustrate that production is larger than nonlinear transfer by a factor of 
$\alpha (\boldsymbol {k})^2$. In parallel shear flows, the lift-up mechanism (Landahl Reference Landahl1980) aligns 
$\hat {\boldsymbol {\psi }}_1(\boldsymbol {k})$ in the 
$x$-component and 
$\hat {\boldsymbol {\phi }}_1(\boldsymbol {k})$ in the 
$y$- and 
$z$-components. This component-type non-normality (Marquet et al. Reference Marquet, Lombardi, Chomaz, Sipp and Jacquin2009) results in small values of 
$\alpha (\boldsymbol {k})$. The optimal resolvent mode, consequently, does not reflect the true physics of the flow for which, even for the low-Reynolds-number flows considered here, production and nonlinear transfer are of comparable magnitude.
As an aside, we note that the same arguments are applicable to cylinder flow, for which Jin et al. (Reference Jin, Symon and Illingworth2020) found that nonlinear transfer for the first resolvent mode was nearly zero. The primary difference for spatially-developing flows is that convective non-normality (Marquet et al. Reference Marquet, Lombardi, Chomaz, Sipp and Jacquin2009) is responsible for minimising 
$\alpha (\boldsymbol {k})$. Mean advection, instead of mean shear, localises 
$\hat {\boldsymbol {\psi }}_1(\boldsymbol {k})$ downstream of the cylinder and 
$\hat {\boldsymbol {\phi }}_1(\boldsymbol {k})$ upstream of the cylinder.
6.2. Counteracting non-normality with eddy viscosity
Returning to the channel flow case, the root of non-normality is the coupling term 
$-\textrm {i}k_yU'$ of the resolvent operator. Since only the spanwise wavenumber appears in this term, higher aspect ratio structures where 
$k_x < k_y$ are more prone to amplification. The simplest way to explain the weakening of non-normality by eddy viscosity is that it adds additional viscosity to counteract mean shear in the operator, resulting in a smaller disparity between 
$\hat {P}(\boldsymbol {k})$ and 
$\hat {N}(\boldsymbol {k})$ for the first resolvent mode.
As a representative example, we plot the first resolvent and eddy mode in figure 11 for 
$(k_x,k_y) = (0,4)$ and compare it to the DNS – i.e. the true Fourier mode for P4U. In terms of the streamwise velocity component, the first eddy mode is in better agreement with the DNS than the first resolvent mode. In terms of the spanwise and wall-normal components, however, neither the eddy nor the resolvent mode are in close agreement with the DNS, which has less smooth behaviour that can be attributed to the fixed convection velocity of the structure. Nevertheless, the eddy mode has more energy in these velocity components than its resolvent counterpart. In table 2, it can be seen that the inner product 
$\left \langle \hat {\boldsymbol {\psi }}_1(\boldsymbol {k}),\hat {\boldsymbol {\phi }}_1(\boldsymbol {k}) \right \rangle$, which is inversely proportional to 
$\alpha (\boldsymbol {k})$, has increased by a factor of nearly three by adding eddy viscosity. In other words, non-normality is weaker, and this can also be observed for the 
$(2,0)$ and 
$(2,4)$ wavenumber pairs.
Figure 11. The (a) streamwise, (b) spanwise and (c) wall-normal velocity components of 
$(k_x,k_y) = (0,4)$. The first resolvent mode is in blue, the first eddy mode is in orange and the DNS is in yellow. The modes are normalised by the peak value of the streamwise velocity component.
Table 2. Measures of non-normality for energetic wavenumbers in P4U.
It is worth mentioning that eddy viscosity is not the only way to counteract mean shear. In Rosenberg & McKeon (Reference Rosenberg and McKeon2019), a componentwise analysis of the resolvent operator yielded two distinct families of modes which, when correctly weighted, destructively interfere to reduce bias towards the streamwise velocity component. The approach has been applied in McMullen, Rosenberg & McKeon (Reference McMullen, Rosenberg and McKeon2020) to higher Reynolds number flows where the destructive interference is more pronounced due to stronger non-normality and higher mean shear.
6.3. Resolvent and eddy reconstructions of the energy budget
Since the first resolvent mode cannot provide the nonlinear transfer that was observed in the DNS, we hypothesise that the role of suboptimal modes in energy transfer is important. We test this hypothesis by reconstructing the energy budget from resolvent modes for P4U since there is a unique wave speed for every wavenumber pair (see appendix B for the minimal channel case). The weights of resolvent modes 
$\tilde {\chi }_p(\boldsymbol {k})$ are determined by projecting them onto the velocity field as done in Sharma et al. (Reference Sharma, Moarref, McKeon, Park, Graham and Willis2016):
\begin{equation} \tilde{\chi}_p(\boldsymbol{k}) = \sigma_p(\boldsymbol{k}) \chi_p (\boldsymbol{k})= \left\langle\hat{\boldsymbol{\psi}}_p (\boldsymbol{k}), \hat{\boldsymbol{u}}(\boldsymbol{k}) \right\rangle. \end{equation}
The approximate velocity field 
$\hat {\boldsymbol {u}}_a(\boldsymbol {k})$ can be written as
\begin{equation} \hat{\boldsymbol{u}}_a(\boldsymbol{k}) = \sum_{p=1}^n \tilde{\chi}_p(\boldsymbol{k}) \hat{\boldsymbol{\psi}}_p(\boldsymbol{k}), \end{equation}
where 
$n$ is the number of resolvent modes used in the approximation. The reconstructed energy budget is evaluated as a function of 
$n$ by replacing 
$\hat {\boldsymbol {u}}(\boldsymbol {k})$ with 
$\hat {\boldsymbol {u}}_a(\boldsymbol {k})$ in (2.12). Since eddy viscosity can lead to improvements, we will also do this for eddy modes by replacing 
$\hat {\boldsymbol {\psi }}_p(\boldsymbol {k})$ with 
$\boldsymbol {\psi }^e_p(\boldsymbol {k})$ in (6.6).
In figure 12(a), the reconstructions are plotted for 
$(k_x,k_y) = (0,4)$, which is both the most energetic structure in the flow and most amplified by the resolvent. Solid and dotted lines denote resolvent and eddy modes, respectively. Production, which appears as red, requires almost 60 resolvent modes to be adequately captured. Moreover, modes 10 to 30 contribute almost zero net production, and therefore modes 30 to 60 are needed. This might seem at odds with the success of the resolvent in identifying sources of production, but it is actually consistent. Because the first resolvent mode is strongly biased towards production, the suboptimal modes cannot balance it with nonlinear transfer unless 
$\tilde {\chi }_1/\sigma _1 \ll 1$, i.e. the projection of 
$\skew4\hat {\boldsymbol {f}}$ onto 
$\hat {\boldsymbol {\phi }}_1$ is small. As figure 12(a) confirms, the optimal resolvent mode contributes very little production to the energy balance. Therefore, the nonlinear forcing has to be heavily biased towards suboptimal modes in order to maintain an energy balance.
Figure 12. (a) Production (red), dissipation (blue) and nonlinear transfer (black) for (a) 
$(k_x,k_y) = (0,4)$, (b) 
$(k_x,k_y) = (2,0)$ and 
$(k_x,k_y) = (2,4)$ in P4U reconstructed from 
$n$ modes. Solid and dotted lines denote resolvent modes and eddy modes, respectively.
The reconstruction using eddy modes, on the other hand, performs reasonably well with a single mode. Although 30 eddy modes are required to converge to the true value, this is considerably less than the 60 resolvent modes needed. Dissipation (blue) and nonlinear transfer (black) converge more quickly to their true values for both sets of modes. We note that the addition of more modes occasionally results in more reconstruction error since the modes are ranked by a kinetic energy norm that is distinct from 
$\hat {P}$, 
$\hat {D}$ and 
$\hat {N}$.
Resolvent and eddy modes perform equally well for 
$(k_x,k_y) = (2,0)$, as seen in figure 12(b). In contrast to 
$(0,4)$, the spanwise-constant structure has negative production and positive nonlinear transfer. The first eddy mode, which was nearly able to capture all terms in the energy budget, is no longer sufficient to reconstruct any term. Although they converge slightly faster than the resolvent modes, both sets require 
$n \approx 20$ in order to converge to within 5 % of their true values.
We finally consider the oblique wave 
$(k_x,k_y) = (2,4)$ in figure 12(c). This structure is less straightforward than the previous two cases. The eddy modes reconstruct production with just 10 modes whereas almost 30 resolvent modes are needed. Dissipation and nonlinear transfer, on the other hand, require more than 30 eddy modes to converge to the true values while they are reconstructed with roughly the same number of resolvent modes needed to reconstruct production (
$n \approx 28$). It can be concluded from figure 12 that the primary benefit of adding eddy viscosity to the operator is in the reconstruction of production. Eddy viscosity may impair reconstruction of dissipation and nonlinear transfer although its impact on these quantities is minor.
6.4. Role of aspect ratio
For a more rigorous comparison between eddy and resolvent modes, we choose an error threshold 
$\mathcal {T}$. The number of modes 
$n$ is gradually increased until every term in the reconstructed energy balance is within 
$\mathcal {T}\%$ of the true value. We then compute the difference 
${\rm \Delta} n = n_{res} - n_{edd}$ between the number of resolvent 
$n_{res}$ and eddy modes 
$n_{edd}$ required and present the results in figure 13. The colour red indicates that eddy modes are more efficient, or 
${\rm \Delta} n > 0$, and thus more resolvent modes must be included to reconstruct the energy budget. Alternatively, the colour blue designates scales for which resolvent modes are more efficient, or 
${\rm \Delta} n < 0$. We also assess the impact of 
$\mathcal {T}$ on 
${\rm \Delta} n$ by setting the threshold to 
$\mathcal {T} = 25\,\%$ in figure 13(a) and reducing it to 
$\mathcal {T} = 1\,\%$ in figure 13(b).
Figure 13. 
${\rm \Delta} n$ for (a) 
$\mathcal {T} = 25\,\%$ and (b) 
$\mathcal {T} = 1\,\%$. Red denotes positive values and blue negative values.
The first observation that can be made regarding figure 13 is that the choice of 
$\mathcal {T}$ has a relatively minor influence on the results. When 
$\mathcal {T}$ is reduced from 25 % to 1 %, only five tiles change sign. If we recall that most of the kinetic energy and transfer processes are among 
$k_x = 0$ and 
$k_x = 2$ scales, then the only important scale impacted by 
$\mathcal {T}$ is 
$(k_x,k_y) = (2,4)$. The reason it is positive for 
$\mathcal {T} = 25\,\%$ is that eddy modes more quickly reconstruct the dominant production term, as seen in figure 12(c). It is negative for 
$\mathcal {T} = 1\,\%$ because the eddy modes take longer to reconstruct dissipation and nonlinear transfer.
The most striking trend for both values of 
$\mathcal {T}$ is that the eddy modes are a more efficient basis when 
$k_x < k_y$. Streamwise-constant structures, in particular, are represented with more than 20 fewer modes if the velocity field is projected onto eddy modes. When we define aspect ratio as 
, then 
is a conservative estimate for when eddy modes are a superior basis to resolvent modes. The pattern is less clear for structures with 
since the sign of 
$\varDelta n$ is sensitive to the choice of 
$\mathcal {T}$. Excluding spanwise-constant structures, those with 
are slightly better represented with resolvent modes although they tend to be weak energetically.
While we only consider P4U in this section, parallel analyses for the minimal channel in appendix B are consistent with the above. We opt to focus only on P4U for simplicity since the Fourier modes for a particular wavenumber pair are deterministic and travel at a single convection velocity. The details for handling statistical variability and multiple convection velocities are discussed in appendix B.
6.5. Implications for resolvent modelling
Despite its shortcomings, resolvent analysis without eddy viscosity successfully identifies linear mechanisms that extract energy from the mean flow. As such, it predicts with remarkable accuracy the energy spectra of turbulent flows as shown by Moarref et al. (Reference Moarref, Sharma, Tropp and McKeon2013). Resolvent analysis has also proven to be a valuable tool in designing flow control strategies in turbulent flows, e.g. Luhar, Sharma & McKeon (Reference Luhar, Sharma and McKeon2014), and spatially-developing flows such as cavities, e.g. Leclercq et al. (Reference Leclercq, Demourant, Puossot-Vassal and Sipp2019).
Resolvent analysis without eddy viscosity is less applicable for estimation. It was shown in Illingworth, Monty & Marusic (Reference Illingworth, Monty and Marusic2018), for example, that adding eddy viscosity improved linear estimation of turbulent channel flow at 
$Re_{\tau } = 1000$. Towne, Lozano-Durán & Yang (Reference Towne, Lozano-Durán and Yang2020) estimated velocity statistics from limited known measurements using resolvent analysis with eddy viscosity. Subsequent work in jets (Pickering et al. Reference Pickering, Towne, Jordan and Colonius2020) highlights the potential of eddy viscosity to provide improved estimates of spatially-developing flows where the energy of resolvent modes extends further downstream than the true fluctuations, e.g. Rosenberg et al. (Reference Rosenberg, Symon and McKeon2019).
7. Conclusions
We have investigated energy transfer for the P4U ECS and low Reynolds number turbulent flow in a minimal channel. For every wavenumber pair, a balance must be achieved across production, dissipation and nonlinear transfer. For both flows, production for the energetic scales is generally positive, and the largest contribution is generated by the near-wall streaks with a spanwise spacing of 
$\lambda _y^+ \approx 100$. Production is negative for some scales such as those that are spanwise-constant. Dissipation is negative for every scale although it is not sufficiently large to counteract production from the most energetic structures. As such, nonlinear transfer redistributes energy to ensure each scale achieves an energy balance. The net effect of nonlinear transfer across all scales is zero, but it is negative for the largest scales and positive for the rest. Spanwise-constant structures are among the largest recipients of energy through nonlinear transfer since both their production and dissipation are negative. These Tollmien–Schlichting-type waves are the first to become unstable for the laminar flow at 
$Re = 5772$ (Orszag Reference Orszag1971), but they play a damping role in the minimal channel.
Energy transfer in the DNS was compared to predictions from resolvent analysis. Similar to DNS, each mode has to satisfy a balance across production, dissipation and nonlinear transfer. Since the first resolvent mode is often representative of the true velocity field, we computed its energy balance and compared it to DNS. The first resolvent mode was successful in identifying the main production mechanisms in the flow. These are the most-amplified structures by resolvent analysis and highlight the role of linear mechanisms in the sustenance of wall-bounded turbulence. For nearly all scales, production was counteracted primarily by dissipation. Nonlinear transfer was nearly zero for the optimal resolvent mode even though it played a major role in redistributing energy in the DNS.
We quantified the additional dissipation introduced by eddy viscosity to model nonlinear transfer for each scale. Its quantitative accuracy was limited to the most energetic mode, and we noted that it primarily removes energy, suggesting it is less applicable for the many scales that receive energy through nonlinear transfer. The gradient of eddy viscosity added energy for some scales although its contribution to the eddy dissipation was negligible. The addition of eddy viscosity had an impact on the number of resolvent modes required to reconstruct the energy budget. We only considered the P4U ECS since there was a unique wave speed and deterministic structure for each wavenumber pair to simplify the analysis. We determined that the eddy basis performed better for high-aspect-ratio structures, particularly with respect to reconstructing production. In some cases, as many as 20 fewer eddy modes than resolvent modes were needed to properly reconstruct all terms in the energy balance. The resolvent basis was slightly better for high aspect ratio structures although these were energetically less significant. A threshold of 
was a conservative estimate for when eddy modes were more efficient.
We showed that eddy viscosity improved the basis for high aspect ratio structures by counteracting non-normality, which was quantified by the inner product (or overlap) between the optimal resolvent forcing and response modes. Non-normality resulted in a trade-off between production and nonlinear transfer in the energy balance equation for the first resolvent mode. Higher non-normality resulted in more production but less nonlinear transfer. The energy budgets of the most amplified scales, i.e those with high singular values, were therefore poorly reconstructed by the optimal resolvent mode since it had almost no overlap with the optimal forcing mode. The damping introduced by eddy viscosity mitigated this trade-off and resulted in a better basis for the velocity field. Its refinement for higher Reynolds number flows could improve the potential of linear models for estimation.
Acknowledgements
The authors wish to thank J. S. Park and M. D. Graham for providing the P4U solution analysed in this article. The authors are also very grateful to M. Xie for providing the DNS data in the case of the minimal channel. Finally, the authors would like to thank the referees for their suggestions, which improved the clarity of the paper.
Funding
The authors acknowledge the financial support of the Australian Research Council.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Linear operators
After elimination of the pressure, the linearised Navier–Stokes equations can be rewritten for the wall-normal velocity 
$\hat {w}$ and wall-normal vorticity 
$\hat {\eta } = \textrm {i}k_y\hat {u}-\textrm {i}k_x\hat {v}$. The matrices 
$\boldsymbol {A}$, 
$\boldsymbol {B}$ and 
$\boldsymbol {C}$ that appear in (2.4) are
\begin{gather} \boldsymbol{A} = \boldsymbol{M}\left[\begin{array}{cc} \mathcal{L}_{OS} & 0 \\ -\textrm{i}k_y U' & \mathcal{L}_{SQ} \end{array} \right], \end{gather}
\begin{gather}\boldsymbol{B} = \boldsymbol{M} \left[\begin{array}{ccc} -\textrm{i} k_x \mathcal{D} & -\textrm{i} k_y \mathcal{D} & -k^2 \\ \textrm{i}k_y & -\textrm{i}k_x & 0 \end{array} \right], \end{gather}
\begin{gather}\boldsymbol{C} = \frac{1}{k^2}\left[\begin{array}{cc} \textrm{i}k_x \mathcal{D} & -\textrm{i}k_y \\ \textrm{i}k_y \mathcal{D} & \textrm{i}k_x \\ k^2 & 0 \end{array} \right]. \end{gather}
Both 
$\mathcal {D}$ and 
$'$ represent differentiation in the wall-normal direction, and 
$k^2 = k_x^2+ k_y^2$. The mass matrix 
$\boldsymbol {M}$ is defined as
\begin{equation} \boldsymbol{M}(k_x,k_y) = \left[\begin{array}{cc} \varDelta^{-1} & 0 \\ 0 & \boldsymbol{I} \end{array} \right], \end{equation}
where 
$\varDelta = \mathcal {D}^2-k^2$ and 
$\boldsymbol {I}$ is the identity matrix. The Orr–Sommerfeld 
$\mathcal {L}_{OS}$ and Squire 
$\mathcal {L}_{SQ}$ operators are
\begin{gather} \mathcal{L}_{OS} = -\textrm{i}k_xU\varDelta + \textrm{i}k_xU'' + (1/Re_{\tau})\varDelta^2, \end{gather}
\begin{gather}\mathcal{L}_{SQ} = -\textrm{i}k_xU + (1/Re_{\tau}) \varDelta. \end{gather}With the addition of eddy viscosity, they become
\begin{gather} \mathcal{L}_{OS} = -\textrm{i}k_xU\varDelta + \textrm{i}k_xU'' + \nu_T \varDelta^2 + 2 \nu_T'\mathcal{D}\varDelta + \nu_T''(\mathcal{D}^2 + k^2), \end{gather}
\begin{gather}\mathcal{L}_{SQ} = -\textrm{i}k_xU + \nu_T \varDelta + \nu_T'\mathcal{D}. \end{gather}Appendix B. Non-normality in the minimal channel case
As alluded to in § 6.4, analysis of the minimal channel is complicated by the fact that each wavenumber pair has a distribution of energetic temporal frequencies. For the sake of simplicity, we choose to analyse the most energetic temporal frequency for each wavenumber pair only. The second complication is that, unlike the P4U case, there is no deterministic mode to describe the velocity field or nonlinear forcing at each 
$\boldsymbol {k}$ in the minimal channel. As discussed by Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018), this stems from the statistical variability in turbulent flows. Instead, we can replace 
$\hat {\boldsymbol {u}}(\boldsymbol {k})$ in (6.5) with the most energetic mode from spectral proper orthogonal decomposition (SPOD) (Karhunen Reference Karhunen1946; Loève Reference Loève1955; Lumley Reference Lumley1970; Picard & Delville Reference Picard and Delville2000), which optimally represents the space-time flow statistics for any 
$\boldsymbol {k}$.
The SPOD modes are computed using the procedure described in Towne et al. (Reference Towne, Schmidt and Colonius2018) and Muralidhar et al. (Reference Muralidhar, Podvin, Mathelin and Fraigneau2019); a brief summary is provided here. Following Welch's method (Welch Reference Welch1967), the DNS data for a particular 
$(k_x,k_y)$ are split into overlapping segments containing 128 snapshots with 50 % overlap and are Fourier-transformed in time. The data are arranged into the matrix
\begin{equation} \hat{\boldsymbol{Q}}(\boldsymbol{k}) = \big[\hat{\boldsymbol{q}}_{\omega}^{(1)} \quad \hat{\boldsymbol{q}}_{\omega}^{(2)} \quad \cdots \quad \hat{\boldsymbol{q}}_{\omega}^{(s)} \big], \end{equation}
where each row represents a different temporal frequency contained in the 
$s$th segment. The cross-spectral density matrix for a specific wavenumber triplet 
$\hat {\boldsymbol {S}}(\boldsymbol {k})$ is
\begin{equation} \hat{\boldsymbol{S}}(\boldsymbol{k}) = \hat{\boldsymbol{Q}}(\boldsymbol{k})\hat{\boldsymbol{Q}}^*(\boldsymbol{k}) . \end{equation}
The SPOD eigenvectors 
$\hat {\boldsymbol {V}}(\boldsymbol {k})$ and eigenvalues 
$\boldsymbol {\varLambda }(\boldsymbol {k})$ can be obtained via an eigenvalue decomposition of the cross-spectral density matrix
\begin{equation} \hat{\boldsymbol{S}}(\boldsymbol{k})\hat{\boldsymbol{V}}(\boldsymbol{k}) = \hat{\boldsymbol{V}}(\boldsymbol{k})\boldsymbol{\varLambda}(\boldsymbol{k}). \end{equation} In figure 14, we compare the most energetic SPOD, resolvent and eddy modes for 
$\boldsymbol {k} = (0,8,0)$. The streamwise component of SPOD and the eddy mode are nearly identical whereas the resolvent mode predicts a much wider structure. The key improvement, however, is observed for the spanwise and wall-normal components. While the SPOD and eddy modes are not equivalent, the eddy mode is much closer to the DNS than the resolvent mode in terms of magnitude. Without the presence of eddy viscosity, the first resolvent mode is more biased towards the streamwise velocity component. In table 3, it is observed that eddy viscosity increases the overlap of the optimal resolvent response and forcing modes by a factor of approximately three. Similar to P4U, eddy viscosity reduces the influence of the mean shear (and hence non-normality) such that the disparity between production and nonlinear transfer is not large. Eddy viscosity also weakens non-normality for the most energetic spanwise-constant structure 
$(2,0)$ and oblique wave 
$(2,8)$.
Figure 14. The (a) streamwise, (b) spanwise and (c) wall-normal velocity components of 
$\boldsymbol {k} = (0,8,0)$ in the minimal channel. The first resolvent mode is in blue, the first eddy mode is in orange and the DNS is in yellow. The modes are normalised by the peak value of the streamwise velocity component.
Table 3. Measures of non-normality for energetic wavenumbers in the minimal channel.
Figure 15 shows the resolvent reconstructions for production and dissipation only. Nonlinear transfer is not included since there is statistical variability in 
$\skew4\hat {\boldsymbol {f}}(\boldsymbol {k})$ depending on the time series that is Fourier-transformed to obtain the 
$\omega = 0$ structure. As such, it would be necessary to perform SPOD for the nonlinear forcing in addition to the velocity field. Production is the term of interest since it is slowest to converge for P4U. The wavenumber pairs considered in figure 15 reinforce that this observation is not specific to P4U. The resolvent reconstruction of the near wall streaks in figure 15(a) is particularly slow for production, requiring over 100 modes. Similar to the P4U results, a plateau region emerges around 
$20 < n < 50$ where the addition of resolvent modes has virtually no impact on production. Dissipation, on the other hand, converges around 
$n = 50$. Only 40 eddy modes, meanwhile, are needed to reconstruct the energy budget. In fact, a substantial portion is captured by the first pair of eddy modes alone, which is consistent with their close resemblance to the SPOD modes in figure 14.
Figure 15. Production (red) and dissipation (blue) for (a) 
$(k_x,k_y) = (0,8)$, (b) 
$(k_x,k_y) = (2,0)$ and (c) 
$(k_x,k_y) = (2,8)$ in the minimal channel reconstructed from 
$n$ modes. Solid and dotted lines denote resolvent modes and eddy modes, respectively.
The reconstructions for 
$(k_x,k_y) = (2,0)$ and 
$(k_x,k_y) = (2,8)$ in figures 15(b) and 15(c), respectively, reinforce the efficiency of eddy modes in reconstructing the energy budget. Nevertheless, the first pairs of eddy modes for these two scales are much less effective than they were for 
$(k_x,k_y) = (0,8)$. Nearly 10 eddy modes are needed for 
$\hat {P}(2,0)$ and 
$\hat {P}(2,8)$ to be close to the true value whereas the first pair of eddy modes alone for 
$\hat {P}(0,8)$ were able to account for nearly 70 % of production.
Finally, we consider 
${\rm \Delta} n$ for the minimal channel case in figure 16. Eddy modes outperform resolvent modes for nearly every wavenumber pair regardless of the threshold 
$\mathcal {T}$ chosen. Due to the dimensions of the computational domain, however, 
for nearly every scale that appears in figure 16. Based on the P4U results, therefore, it is expected that eddy modes will be more efficient.
Figure 16. 
${\rm \Delta} n$ for (a) 
$\mathcal {T} = 25\,\%$ and (b) 
$\mathcal {T} = 1\,\%$. Red denotes positive values and blue negative values.
