2.1. Equations of motion

We start by introducing a decomposition of an instantaneous velocity field, $\boldsymbol{u}(\boldsymbol{x},t)$ , as follows:

(2.1) \begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \boldsymbol{U}(\boldsymbol{x},t) + \boldsymbol{u'}(\boldsymbol{x},t), \end{equation}

where $\boldsymbol{U}(\boldsymbol{x},t)(=(U(y,t),0,W(y,t)))$ and $\boldsymbol{u'}(\boldsymbol{x},t)$ represent the velocity averaged over $x$ and $z$ and the fluctuating velocity, respectively. Further, following Doohan et al. (Reference Doohan, Bengana, Yang, Willis and Hwang2022), where detailed interaction energetics of the self-sustaining process are studied at different length scales, we decompose the fluctuating velocity into

(2.2) \begin{equation} \boldsymbol{u'}(\boldsymbol{x},t) = \boldsymbol{u}_1( y,z,t) + \boldsymbol{u}_2(\boldsymbol{x},t), \end{equation}

with $\boldsymbol{u}_1(y,z,t) = \langle \boldsymbol{u'}(\boldsymbol{x},t)\rangle _{x}$ , where $\langle \boldsymbol{\cdot }\rangle _{x}$ denotes the average in the $x$ -direction. We note that the fluctuating flow averaged in the streamwise direction, $\boldsymbol{u}_1$ , characterises the streamwise elongated flow structures (streaks), while the rest of the fluctuating velocity represents the streamwise varying structures (meandering motions of streaks due to streak instability or transient growth and quasi-streamwise vortices). The pressure can be decomposed similarly, i.e. $P(\boldsymbol{x},t) = P + p_1+p_2$ . The resulting dimensionless equations of motion are given by

(2.3) \begin{align} \frac {\partial \boldsymbol{U}}{\partial t} + (\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U} + \langle (\boldsymbol{u'} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u'} \rangle _{x,z} & = -\boldsymbol{\nabla }P + \frac {1}{Re}{\nabla} ^2 \boldsymbol{U}, \end{align}

(2.4) \begin{align} \frac {\partial \boldsymbol{u}_1}{\partial t} + (\boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U} + (\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_1 + (\boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_1 &+ \langle (\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_2 \rangle _{x} - \langle (\boldsymbol{u'} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u'} \rangle _{x,z}\nonumber\\ & = -\boldsymbol{\nabla }p_1 + \frac {1}{Re} {\nabla} ^2 \boldsymbol{u}_1, \end{align}

(2.5) \begin{align} \frac {\partial \boldsymbol{u}_2}{\partial t} + [(\boldsymbol{U} + \boldsymbol{u}_1) \boldsymbol{\cdot }\boldsymbol{\nabla }] \boldsymbol{u}_2 + (\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) & [\boldsymbol{U} + \boldsymbol{u}_1] + (\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_2 - \langle (\boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_2\rangle _{x}\nonumber\\ & = -\boldsymbol{\nabla }p_2 + \frac {1}{Re} {\nabla} ^2 \boldsymbol{u}_2, \end{align}

where $\langle \boldsymbol{\cdot }\rangle _{x,z}$ denotes the average in the $x$ - and $z$ -directions. We emphasise that the velocity decomposition introduced here is intended specifically for minimal flow units. In such domains, the restriction on admissible streamwise length scales ensures that the streamwise averaged component $\boldsymbol{u}_1$ provides a meaningful characterisation of streak dynamics, while the remaining component $\boldsymbol{u}_2$ captures their streamwise meandering motions caused by the instability or transient growth. In very large domains, where streaks of widely varying streamwise extent and multiple meandering or instability wavelengths coexist simultaneously, this decomposition and the associated global diagnostics derived from it are not, in general, expected to provide a complete description of the dynamics.

2.2. Perturbation energy budget

We now consider a perturbation applied to a given reference flow field from a statistically steady turbulent state:

(2.6) \begin{equation} \Delta \boldsymbol{u} = \boldsymbol{u}^{(2)} - \boldsymbol{u}^{(1)}, \end{equation}

where $\boldsymbol{u}^{(1)}$ is the ‘unperturbed’ reference flow field and $\boldsymbol{u}^{(2)}$ is the corresponding perturbed one with $\Delta \boldsymbol{u}$ . Using (2.6), and the equations of motion given in (2.3), (2.4) and (2.5), the equations for each componet of the perturbation field (i.e. $\Delta \boldsymbol{U}$ , $\Delta \boldsymbol{u}_1$ and $\Delta \boldsymbol{u}_2$ ) are obtained, as documented in Appendix A: see (A1), (A2) and (A3). Using the corresponding energy equations for each perturbation field (see (A4), (A5) and (A6) in Appendix A), the time evolution of the volume-averaged perturbation energy for $\Delta \boldsymbol{U}$ , $\Delta \boldsymbol{u}_1$ and $\Delta \boldsymbol{u}_2$ is given by

(2.7) \begin{align} \frac {\partial E_{\Delta \boldsymbol{U}}}{\partial t} &= P_{\Delta \boldsymbol{U}}+T_{\Delta \boldsymbol{U}}-\epsilon _{\Delta \boldsymbol{U}}, \end{align}

(2.8) \begin{align} \frac {\partial E_{\Delta \boldsymbol{u}_1}}{\partial t} &= P_{\Delta \boldsymbol{u}_1}+T_{\Delta \boldsymbol{u}_1}-\epsilon _{\Delta \boldsymbol{u}_1}, \end{align}

(2.9) \begin{align} \frac {\partial E_{\Delta \boldsymbol{u}_2}}{\partial t} &= P_{\Delta \boldsymbol{u}_2}+T_{\Delta \boldsymbol{u}_2}-\epsilon _{\Delta \boldsymbol{u}_2}, \end{align}

where

(2.10) \begin{equation} E_{\Delta \boldsymbol{U}}= \frac {1}{2}\langle \Delta \boldsymbol{U} \boldsymbol{\cdot }\Delta \boldsymbol{U} \rangle , \quad E_{\Delta \boldsymbol{u}_1}= \frac {1}{2}\langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\Delta \boldsymbol{u}_1 \rangle \quad \text{and} \quad E_{\Delta \boldsymbol{u}_2}= \frac {1}{2}\langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\Delta \boldsymbol{u}_2 \rangle \end{equation}

are the volume-averaged kinetic energies of the three flow components ( $\langle \boldsymbol{\cdot }\rangle \equiv 1/V \int _\varOmega (\boldsymbol{\cdot })\, \textrm {d}V,$ where $\varOmega$ is the given flow domain and $V$ is its volume). Here, the three production terms are written as

(2.11) \begin{align} P_{\Delta \boldsymbol{U}} = & -\biggl (\left \langle \Delta \boldsymbol{U} \boldsymbol{\cdot }(\Delta \boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U}^{(1)} \right \rangle +\left \langle \Delta \boldsymbol{U} \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_1^{(1)} \right \rangle \nonumber\\ &+ \left \langle \Delta \boldsymbol{U} \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_2^{(1)} \right \rangle \biggr ), \end{align}

(2.12) \begin{align} P_{\Delta \boldsymbol{u}_1} = & -\biggl ( \left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{U}^{(1)} \right \rangle + \left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_1^{(1)} \right \rangle \nonumber\\ &+\left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_1^{(1)} \right \rangle + \left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_2^{(1)} \right \rangle \biggr ), \end{align}

(2.13) \begin{align} P_{\Delta \boldsymbol{u}_2}= & -\biggl ( \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_2^{(1)} \right \rangle + \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_2^{(1)} \right \rangle \nonumber\\ &+\left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{U}^{(1)} \right \rangle + \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_1^{(1)} \right \rangle \nonumber\\ &+ \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}_2^{(1)} \right \rangle \biggr ). \end{align}

Similarly, the transport terms are given by

(2.14) \begin{align} T_{\Delta \boldsymbol{U}}= & -\biggl ( \left \langle \Delta \boldsymbol{U} \boldsymbol{\cdot }(\boldsymbol{u}_1^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_1 \right \rangle + \left \langle \Delta \boldsymbol{U} \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_1 \right \rangle \nonumber\\ &+ \left \langle \Delta \boldsymbol{U} \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_2 \right \rangle + \left \langle \Delta \boldsymbol{U} \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_2 \right \rangle \biggr ), \end{align}

(2.15) \begin{align} T_{\Delta \boldsymbol{u}_1}= & -\biggl ( \left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\boldsymbol{u}_1^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{U} \right \rangle + \left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{U} \right \rangle \nonumber\\ &+\left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_2 \right \rangle + \left \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_2 \right \rangle \biggr ), \end{align}

(2.16) \begin{align} T_{\Delta \boldsymbol{u}_2}= &-\biggl ( \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{U} \right \rangle + \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{U} \right \rangle \nonumber\\ &+ \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_1 \right \rangle + \left \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla })\Delta \boldsymbol{u}_1 \right \rangle \biggr ). \end{align}

Finally, the dissipation and diffusion terms are

(2.17) \begin{equation} \begin{aligned} \epsilon _{\Delta \boldsymbol{U}} &= -\left \langle \frac {1}{Re} \Delta \boldsymbol{U} \boldsymbol{\cdot }{\nabla} ^2 \Delta \boldsymbol{U} \right \rangle , \quad\; \epsilon _{\Delta \boldsymbol{u}_1} = -\left \langle \frac {1}{Re} \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }{\nabla} ^2 \Delta \boldsymbol{u}_1 \right \rangle , \quad\; \nonumber\\ \epsilon _{\Delta \boldsymbol{u}_2} &= -\left \langle \frac {1}{Re} \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }{\nabla} ^2 \Delta \boldsymbol{u}_2 \right \rangle . \end{aligned} \end{equation}

It can be seen that (2.7), (2.8) and (2.9) are made up of a term expressing the time evolution of the energy, a set of production terms, a set of transport terms and a viscous term, respectively. In general, the production terms contribute to the growth of the perturbation energy and the viscous terms take energy away from the perturbation energy. The transport terms do not add or subtract any energy to the ‘total’ energy budget (see (2.18)), as they cancel each other when added together, but rather redistribute energy amongst the three different components ( $\Delta \boldsymbol{U}$ , $\Delta \boldsymbol{u}_1$ and $\Delta \boldsymbol{u}_2$ ). Therefore, in this work, we use the term ‘transport’ in the sense that they add/subtract energy to one component and subtract/add it to another, without affecting the total budget. This is indeed the behaviour observed for these terms and will be analysed in § 3.

When adding (2.7), (2.8) and (2.9), the transport terms do not add or subtract energy from the system as a whole and cancel with each other. Therefore, adding the three equations leads to the following perturbation energy balance equation:

(2.18) \begin{equation} \frac {\mathrm{d} E_\Delta }{\mathrm{d}t} = P_\Delta - \epsilon _\Delta \end{equation}

with

(2.19) \begin{equation} E_\Delta = \frac {1}{2}\langle \Delta \boldsymbol{u} \boldsymbol{\cdot }\Delta \boldsymbol{u} \rangle , \end{equation}

where production and dissipation, $P_\Delta$ and $\epsilon _\Delta$ , are defined as

(2.20) \begin{equation} P_\Delta = -\left \langle \Delta \boldsymbol{u} \boldsymbol{\cdot }(\Delta \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}^{(1)} \right \rangle , \quad \epsilon _\Delta = \frac {2}{Re}\langle s_{\Delta ,ij}s_{\Delta ,ij}\rangle \end{equation}

with

(2.21) \begin{equation} s_{\Delta ,ij} = \frac {1}{2} \left ( \frac {\partial \Delta u_i}{\partial x_j} + \frac {\partial \Delta u_j}{\partial x_i} \right ). \end{equation}

Here, the dissipation term is written in tensor notation for simplicity.

2.3. Numerical set-ups

Direct numerical simulations in the present study are carried out using the open-source code Channelflow 2.0 ( https://www.channelflow.ch/ ), a software to analyse incompressible channel flows. The solver is spectral with the Fourier–Galerkin method used in the wall-parallel dimensions and the Chebyshev-tau method used in the wall-normal dimension. In time, the Navier–Stokes equations are integrated using a third-order semi-implicit backward differentiation scheme. The domain is discretised into a grid of $16 \times 33 \times 16$ cells in $x$ , $y$ and $z$ , respectively (after dealiasing), the same discretisation as the one used by Hamilton et al. (Reference Hamilton, Kim and Waleffe1995). The friction Reynolds number is obtained as $Re_\tau$ = 34 from the simulations, where $Re_\tau =u_\tau h/\nu$ and $u_\tau$ is the friction velocity.

3.1. Short-term response: Lyapunov evolution

We start by evaluating the time evolution of the three energy components, $E_{\Delta \boldsymbol{U}}$ , $ E_{\Delta \boldsymbol{u}_1}$ and $ E_{\Delta \boldsymbol{u}_2}$ , for $t\in [0,300]$ , together with the corresponding production and dissipation. The time evolution of these energy components is shown in figure 2(a). After the initial short transient period, the perturbation energy grows exponentially, governed by the leading Lyapunov exponent. The dominant energy term is $ E_{\Delta \boldsymbol{u}_2}$ , with the second largest contributor being $ E_{\Delta \boldsymbol{u}_1}$ followed by $ E_{\Delta \boldsymbol{U}}$ . This order is in direct opposition with the energy distribution in the reference field, where we found the energy of the reference field $\boldsymbol{U}^{(1)}$ , $\boldsymbol{u}_1^{(1)}$ and $\boldsymbol{u}_2^{(1)}$ components to amount to 68 %, 28 % and 4 % of the total energy (defined as the sum of the three component energies), respectively. This distribution follows the expected trends, given that the main flow contains the largest amount of energy and the (larger) streaky structures in Couette flow are more energetic than the small scales. This shows how, despite the reference field $\boldsymbol{u}_2^{(1)}$ component only containing a small fraction of the reference field energy, it clearly contains the largest amount of perturbation energy.

The exponential growth trend of the three components is evident, although small fluctuations are noticeable, likely arising from the ensemble-averaging process. Given the dominant role played by the smallest chaotic scales in the perturbation energy growth during this period and seeing the simultaneous growth of the three components observed in figure 2(a), this suggests that the growth of $\Delta \boldsymbol{u}_2$ contributes to the growth of the other components through the coupling terms in the budget equations (see (2.7), (2.8) and (2.9) and further discussion later). The production and dissipation of each of the three energy components are shown in figure 2(b). Here, the production and dissipation of the perturbation energy are normalised by $E_{\Delta }$ , so that their evolution is visualised by removing the contribution of the leading Lyapunov exponent. The dominant production and dissipation is that of $\Delta \boldsymbol{u}_2$ , which, together with the dominance of $ E_{\Delta \boldsymbol{u}_2}$ , confirms the previous findings of Inubushi et al. (Reference Inubushi, Takehiro and Yamada2015) and Nikitin (Reference Nikitin2018). Furthermore, in line with the previous discussion, the productions for $\Delta \boldsymbol{u}_1$ and $\Delta \boldsymbol{U}$ are positive and slightly higher than the corresponding dissipations, indicating that the growth dominated by $\Delta \boldsymbol{u}_2$ indeed affects that of $\Delta \boldsymbol{u}_1$ and $\Delta \boldsymbol{U}$ .

Figure 3.

Time evolution of production terms normalised by the rate of change of perturbation energy from (2.9). Asterisks ( $\ast$ ) mark the individual terms. Thicker lines are used to represent the dominant terms for ease of visualisation.

The contributions of the component-wise production and transport terms outlined in § 2.2 are now evaluated. Figure 3 shows the production terms from (2.9), which are found to contain the leading mechanisms. All the production and transport terms are further shown in figure 13, in Appendix D, but are omitted from the main text to ease the visualisation. As seen on figure 3 (corresponding to figure 13 e in Appendix D), the largest contributions to the production of perturbation energy are made by two comparable terms, $-\langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_1^{(1)} \rangle$ and $-\langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_2^{(1)} \rangle$ from $P_{\Delta \boldsymbol{u}_2}$ in (2.13). The first term captures the energy production of $\Delta \boldsymbol{u}_2$ caused by the shear of the streamwise-averaged reference field fluctuations (i.e. the shear created by streaks) and the second term captures the production by the shear of the streamwise dependent flow, which has been understood to be crucial for turbulent dissipation in the self-sustaining process (Hernandez & Hwang Reference Hernández and Hwang2020; Doohan et al. Reference Doohan, Willis and Hwang2021, Reference Doohan, Bengana, Yang, Willis and Hwang2022; Hernandez, Yang & Hwang Reference Hernández, Yang and Hwang2022). The dominance of these terms suggests that during the early stage of the growth of perturbations, their production is primarily related to turbulent dissipation processes, consistent with observations in HIT (Mohan et al. Reference Mohan, Fitzsimmons and Moser2017; Ge et al. Reference Ge, Rolland and Vassilicos2023).

We move on to the analysis of the transport terms, in figure 13(b,d,f). When adding the three spatially averaged perturbation energy evolution equations (i.e. (2.7), (2.8) and (2.9)), the net effect of the transport terms in the global perturbation energy evolution is zero. The primary transport terms are $- \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }) \Delta \boldsymbol{u}_1 \rangle$ and its counterpart, $- \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }) \Delta \boldsymbol{u}_2 \rangle$ . These terms characterise the transfers of energy between $E_{\Delta \boldsymbol{u}_1}$ and $ E_{\Delta \boldsymbol{u}_2}$ , which respectively represent the streaky structures of the perturbation field, and the streamwise dependent flow of the perturbation field, occurring via the advection of the shear of one of the two components by $\boldsymbol{u}^{(2)}$ . However, they are much smaller than the leading production terms shown in figure 13(e) and fluctuate around zero, suggesting that their roles in perturbation amplification are minor. We also note that the interactions represented by some of these transport terms are expected to be negligible during the short-term response (which is to a large extent linear), since they correspond to triple interactions of the perturbations with themselves. These interactions are of much lesser importance than those with the reference field, which agrees with the observations made from figure 13. For these reasons, and to ease the visualisation of the dominant terms, these terms have not been shown in the main text.

Figure 4.

Time evolution of perturbation energy.

3.2. Long-term response: beyond the Lyapunov time

Now, we explore the long-term response phase, which comes once the asymptotic phase response, characterised by an exponential evolution through the leading Lyapunov exponent (i.e. the short-term response), saturates. The time evolution of the perturbation energy in this phase is shown in figure 4. It can be seen that after $t\approx 300$ , the perturbation energy keeps growing until $t\approx 3000$ , when it saturates. The nature of this growth will be discussed in detail (see § 4) after the key contributing terms have been identified (see later).

We first evaluate the time evolution of $ E_{\Delta \boldsymbol{U}}$ , $E_{\Delta \boldsymbol{u}_1}$ and $ E_{\Delta \boldsymbol{u}_2}$ , and the corresponding production and dissipation, as shown in figure 5. Starting by analysing the evolution of the energy components, it can be seen that up to $t\approx 750$ , the dominant energy term is $ E_{\Delta \boldsymbol{u}_2}$ . This is a continuation of the behaviour and trends observed in the previous phase. However, the value of $ E_{\Delta \boldsymbol{u}_2}$ saturates at $t\approx 600$ . From $t\approx 750$ onwards, $ E_{\Delta \boldsymbol{u}_1}$ becomes the dominant term, exhibiting a continuous growth trend until $t\approx 3000$ , at which point it reaches its saturation value. This is consistent with previous observations of propagating the perturbation to progressively larger scales observed in HIT (Boffetta & Musacchio Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018), as $\Delta \boldsymbol{u}_2$ may be regarded as representing small scales, given its direct involvement in dissipation associated with the turbulent energy cascade (Hernandez & Hwang Reference Hernández and Hwang2020; Doohan, Willis & Hwang Reference Doohan, Willis and Hwang2021; Hernandez et al. Reference Hernández, Yang and Hwang2022). Therefore, from this point on, the main contributor to perturbation growth becomes the streaky field $\Delta \boldsymbol{u}_1$ , replacing $\Delta \boldsymbol{u}_2$ that dominated growth throughout the short-term exponential phase. In a similar fashion to the short-term behaviour analysis, presented in § 3.1, we can draw a comparison between the behaviour of the reference and perturbation energy distributions. We see that, once again, the term with the largest amount of perturbation energy, in this case, $E_{\Delta \boldsymbol{u}_1}$ , is in direct contrast with the reference field energy distribution, where the $\boldsymbol{U}^{(1)}$ component is the dominant one.

Figure 5.

Time evolution of (a) perturbation energy for the three components, $ E_{\Delta \boldsymbol{U}}$ , $ E_{\Delta \boldsymbol{u}_1}$ and $ E_{\Delta \boldsymbol{u}_2}$ , and (b) production and dissipation of each component.

Turning the focus now to the evolution of the production and dissipation of each of the energy components, shown in figure 5(b), the general trends of the three different components show the same saturation behaviour previously observed in the time evolution of the perturbation energy components. For the case of $ E_{\Delta \boldsymbol{u}_1}$ , the production and dissipation grow, in accordance with the observations from the growth of $ E_{\Delta \boldsymbol{u}_1}$ . An interesting behaviour is seen in the evolution of $ P_{\Delta \boldsymbol{u}_2}$ and $ \epsilon _{\Delta \boldsymbol{u}_2}$ , with dissipation overpowering production shortly after the saturated stage kicks in. Since the energy remains constant during this period (see figure 5a), other mechanisms must be feeding energy into $ E_{\Delta \boldsymbol{u}_2}$ , which will be seen from the evolution of the production and transport terms (see the subsequent discussion, with figure 6). The same argument is true for the difference in the saturation values between $ P_{\Delta \boldsymbol{u}_1}$ and $\epsilon _{\Delta \boldsymbol{u}_1}$ during the saturated stage.

Figure 6.

Time evolution of (a,c) production and (b,d) transport terms. Panels (a,b) are terms from (2.8) and panels (c,d) from (2.9). Thicker lines are used to represent the dominant terms for easiness of visualisation.

To understand the physical mechanisms that characterise the perturbation growth, the time evolution of the production and transport terms for each of the three energy components is now evaluated. Following a similar pattern to that used to analyse these terms in § 3.1, the terms containing the dominant and most relevant mechanisms are shown in figure 6 (in this case, the production and transport terms from (2.8) and (2.9)), whilst all the terms are shown in figure 14 in Appendix D. Starting by the analysis of the production terms, shown in figure 6(a, c), the trends seen are consistent with those seen when evaluating the time evolution of production. Starting with figure 6(a), corresponding to the production terms of $ E_{\Delta \boldsymbol{u}_1}$ , it is observed that there is one clearly dominant term that contributes to the growth in time of the perturbation energy: $- \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U}^{(1)} \rangle$ . This term represents the contribution of the interactions between the difference field streaks $\Delta \boldsymbol{u}_1$ and the shear of the reference field mean flow $\boldsymbol{U}^{(1)}$ . Therefore, the growth of $ E_{\Delta \boldsymbol{u}_1}$ previously reported is primarily driven by shear in the reference field mean flow. Similarly to this case, there is production generated by shear of the reference field streaks, $- \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_1^{(1)} \rangle$ , although this is much smaller than that by shear of the reference field mean flow. It is also worth noting that $- \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_2^{(1)} \rangle$ seems to follow a slightly decreasing trend, having a considerable region of negative values for $t\geq 1000$ . Therefore, this term acts as a sink in the late stage of perturbation growth, unlike its role in the early stage (see figure 13 c).

Subsequently, looking at figure 6(c), showing the production terms of the energy budget equation of $ E_{\Delta \boldsymbol{u}_2}$ , the two clearly dominant terms are $- \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U}^{(1)} \rangle$ and $-\langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_1^{(1)} \rangle$ . The latter was found to be dominant together with $-\langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_2 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_2^{(1)} \rangle$ during the exponential phase (i.e. the short-term response). This follows the same pattern observed in the evolution of the production terms corresponding to $ E_{\Delta \boldsymbol{u}_1}$ , where the two dominant terms are associated with shear of the mean and streaky flows in the reference field.

Finally, briefly analysing figure 14(a) in Appendix D, corresponding to the terms of the energy production of $ E_{\Delta \boldsymbol{U}}$ , it can be seen that the dominant term, on a much lesser scale than those of the other energy balance equations, is $- \langle \Delta \boldsymbol{ U} \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}_1^{(1)} \rangle$ , which characterises the contribution of the interactions between the shear of the reference field streaks advected by the perturbation field streaks and the perturbation field mean.

Moving on to the transport terms, shown in figure 6(b,d), it can be observed that there is one clear transport mechanism that occurs between $ E_{\Delta \boldsymbol{u}_1}$ and $ E_{\Delta \boldsymbol{u}_2}$ . The two leading terms involved in this transport are $- \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }) \Delta \boldsymbol{u}_2 \rangle$ and $- \langle \Delta \boldsymbol{u}_2 \boldsymbol{\cdot }(\boldsymbol{u}_2^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }) \Delta \boldsymbol{u}_1 \rangle$ (which are shown in figure 6 b,d and in figure 14 d,f, respectively). The former has a negative value and the latter is positive, while their magnitudes are approximately in balance. This indicates that some energy in $\Delta \boldsymbol{u}_1$ is transferred to $\Delta \boldsymbol{u}_2$ and eventually dissipated in $\Delta \boldsymbol{u}_2$ .

4.1. Algebraic growth of perturbations

In figure 7(a), we examine how $ E_{\Delta \boldsymbol{u}_1}$ evolves over time from $t\approx 1000$ to $t\approx 2500$ , after the short-term exponential growth phase. The growth appears to favour a linear growth over time, as found by Boffetta & Musacchio (Reference Boffetta and Musacchio2017) and Berera & Ho (Reference Berera and Ho2018) for HIT (see also the later analysis). However, it must be pointed out that the flow in the present study is at a very low Reynolds number, where there is no clear separation between the length scales of production and dissipation. It is also highly anisotropic and inhomogeneous. Furthermore, the time required for $\Delta \boldsymbol{u}_1$ to be decorrelated is $O(10^3)h/{U_w}$ , an order of magnitude higher than the integral time scale of $O(10^2)h/{U_w}$ in this flow, since the typical time period of the self-sustaining process is only $T U_w/h\approx 70{-}80$ (Hamilton et al. Reference Hamilton, Kim and Waleffe1995). It remains unclear whether this exceptionally long decorrelation time is a particular feature of the minimal flow unit or of wall-bounded turbulence. However, it is evident that the theoretical framework proposed by early studies (Lorenz Reference Lorenz1969; Leith & Kraichnan Reference Leith and Kraichnan1972; Boffetta & Musacchio Reference Boffetta and Musacchio2017; Berera & Ho Reference Berera and Ho2018; Ge et al. Reference Ge, Rolland and Vassilicos2025) is not applicable to the present observations. The large standard deviation of the data shown in the figure suggests that large fluctuations around the mean are an essential feature of this system. This agrees with the chaotic nature of the system, which is characterised by the large differences in trajectory resulting from small differences in initial conditions. We note that, despite the standard deviation contour covering some negative values, the energy cannot become negative, and this results from simply adding and subtracting the value of the standard deviation to the ensemble average at each point.

To better understand the underlying physical process of the linear growth of the perturbation energy over time, the terms present in the energy balance (2.8) for $\Delta \boldsymbol{u}_1$ are shown in figure 7(b). Similarly to $ E_{\Delta \boldsymbol{u}_1}$ , both production and dissipation are seen to grow linearly in time from $t\approx 1000$ to $t\approx 2500$ , whilst transport decreases linearly, acting as a sink feeding energy into $ E_{\Delta \boldsymbol{u}_2}$ (see also figure 6 b,d for the related energetics). From figure 6(a), it is clear that the production of the perturbation energy is dominated by $- \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U}^{(1)} \rangle$ . Therefore, since $E_{\Delta \boldsymbol{u}_1} \sim \mathcal{O}(\Delta \boldsymbol{u}_1)^2$ , for $t \in [1000,2500]$ , a dimensional form of production for $\Delta \boldsymbol{u}_1$ is estimated to be

(4.1) \begin{equation} P_{\Delta \boldsymbol{u}_1} \sim \mathcal{O}\left (\frac {E_{\Delta \boldsymbol{u}_1} u_\tau }{h} \right ), \end{equation}

where the shear of $\boldsymbol{U}^{(1)}$ , expected to be dominated by $\partial U^{(1)}/\partial y$ , is set to $O(u_\tau /h)$ due to the poor separation between the viscous inner and inertial outer length scales at the present very low Reynolds number ( $u_\tau$ is the friction velocity). Further, we assume $T_{\Delta \boldsymbol{u}_1} \sim C_tt$ from figure 7(b), where $C_t$ is a negative constant characterising the amount of transport, although it is small. We also set the dimensional form of dissipation to be

(4.2) \begin{equation} \epsilon _{\Delta \boldsymbol{u}_1} \sim \mathcal{O}\left (\nu \frac {E_{\Delta \boldsymbol{u}_1} }{l_\Delta ^2} \right ), \end{equation}

where $l_\Delta$ is a representative length scale of $\Delta \boldsymbol{u}_1$ in the $y$ – $z$ plane and is expected to change in time without loss of generality. Then, for $t \in [1000,2500]$ , the resulting energy balance of (2.8) is expected to be

(4.3) \begin{equation} \frac {\partial E_{\Delta \boldsymbol{u}_1}}{\partial t} = \mathcal{O}\left (t \frac {u_\tau }{h} \right ) + \mathcal{O}(C_tt)-\mathcal{O}\left ( \frac {\nu t}{ l_\Delta ^2} \right ) \sim \text{constant}. \end{equation}

The mentioned energy balance now indicates that $l_\Delta$ will need to be approximately constant for $E_{\Delta \boldsymbol{u}_1} \sim t$ , so that the linear time dependences in the terms of production, transport and dissipation are balanced out. Only in this case does a constant value of $\partial E_{\Delta \boldsymbol{u}_1}/\partial t$ appear to be possible. Since the length scale of the perturbation has been understood to grow in time in general (Lorenz Reference Lorenz1969; Leith & Kraichnan Reference Leith and Kraichnan1972), this is presumably achieved when the length scale of the perturbation is grown enough, so that it can no longer grow due to the length scale posed by the boundary: i.e. $l_\Delta \sim O(h)$ .

Figure 8.

Time evolution of $ E_{\Delta \boldsymbol{u}_1}$ for a single realisation. The points shown are located at $t=200,425,1250,2225,3650$ , corresponding to the snapshot times shown in figure 9.

To confirm the above-mentioned scenario, we visualise the spatial distribution of the perturbation field at different time instances in a single realisation. Figure 8 shows the time evolution of $E_{\Delta \boldsymbol{u}_1}$ for a single realisation. It includes a set of times ( $t=200, 425, 1250, 2225, 3650$ ), at which the velocity fields are visualised in figure 9. In figure 9(a,c,e,g,i), $\Delta u_1$ is shown, while panels (b,d,f,h,j) display $\Delta v_1$ at the corresponding times. Figure 9(a–d) correspond to the exponential growth and posterior transition period before the linear growth. It is seen that both $\Delta u_1$ and $\Delta v_1$ have large magnitudes at certain localised positions in space. As time progresses, the components spread and single patches occupy large spaces in the box. Figure 9(e–h) correspond to the linear growth period and there does not seem to be any further substantial growth on the length scale in the components. A similar behaviour is observed in the remaining time step, during the saturated period (figure 9 i, j).

Figure 9.

Contours of the perturbation velocity components (a,c,e,g,i) $\Delta u_1$ and (b,d,f,h,j) $\Delta v_1$ for a single realisation shown in the $y$ – $z$ plane. The times of each snapshot are (a,b) $t=200$ , (c,d) $t=425$ , (e,f) $t=1250$ , (g, h) $t=2225$ and (i,j) $t=3650$ . See figure 8.

Figure 10.

Time evolution of the integral length scale in the (a) wall-normal and (b) spanwise directions for the three different velocity components of $\Delta \boldsymbol{u}_1$ . Here, the correlation for the integral length scale is taken at $y=0$ . The integral length scales have been calculated using: $l = \int _{0}^{\infty }f(r)\,{\rm d}r$ , where $f(r)$ is the autocorrelation function.

To statistically characterise the above-mentioned observation, figure 10 shows the time evolution of the integral length scales of $\Delta \boldsymbol{u}_1$ obtained at the channel centre. Since $\Delta \boldsymbol{u}_1$ has no streamwise dependence, the evolution of the wall-normal and spanwise integral length scales is shown for each of the velocity components. The results shown confirm the hypothesis set for the previous theoretical development. The length scale of the perturbation field grows during the initial exponential, Lyapunov-dominated, growth period of the perturbation and during the transition period before the linear growth kicks in. Once the length scale can no longer grow and reach an approximately constant, the growth of the perturbation energy slows down and is approximately proportional to $t$ : i.e. $l_\Delta \sim O(h)$ for $E_{\Delta \boldsymbol{u}_1} \sim t$ .

4.2. Lift-up effect as the mechanism of long-term perturbation growth

The late stage evolution of a perturbation field visualised in figure 9(g–j) depicts patches of large (absolute) magnitude streamwise and vertical velocity in the top/bottom and central regions of the channel. This is an interesting analogue to the structures observed in other turbulent Couette flow fields: for example, compare figure 9(i) with figure 7(a) from Zong et al. (Reference Zong, Maia, Cavalieri and Hwang2025). Furthermore, in figure 11, it can be seen that the energy of $\Delta \boldsymbol{u}_1$ is clearly dominated by its streamwise component for $t\gtrsim 400$ and, therefore, the streamwise component of $E_{\Delta \boldsymbol{u_1}}$ is the primary driver of the perturbation growth.

Figure 11.

Time evolution of $E_{\Delta u_1}$ , $ E_{\Delta v_1}$ and $ E_{\Delta w_1}$ .

These observations motivate a deeper analysis of the nature of the dominant production term in figure 6(a): $- \langle \Delta \boldsymbol{u}_1 \boldsymbol{\cdot }(\Delta \boldsymbol{u}_1 \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{U}^{(1)} \rangle$ . Given the flow setting considered, this term can be expanded into $- \left \langle \Delta u_1 \Delta v_1 ( {\partial U^{(1)}}/{\partial y}) \right \rangle - \left \langle \Delta w_1 \Delta v_1 ({\partial W^{(1)}}/{\partial y} )\right \rangle$ , where, as it can be deduced from figure 11, the former is clearly dominant over the latter. Since the dominant term is a production term, dividing by $\Delta u_1$ gives a velocity evolution expression:

(4.4) \begin{equation} \frac {\partial \Delta u_1}{\partial t} \propto -\Delta v_1 \frac {\partial U^{(1)}}{\partial y}. \end{equation}

The form of (4.4) resembles the mathematical definition of the well-known lift-up effect. Much work has been dedicated to understanding the contributions from this physical wall turbulence mechanism (e.g. Ellingsen & Palm Reference Ellingsen and Palm1975; Landahl Reference Landahl1980, Reference Landahl1990; Butler & Farrell Reference Butler and Farrell1993; Hwang & Cossu Reference Hwang and Cossu2010a ). As explained by Landahl (Reference Landahl1980), considering infinitesimal disturbances ( $u^*$ , $v^*$ , $w^*$ ) applied to an inviscid shear flow (with mean velocity $\overline {U}(y)$ ) and integrating its streamwise momentum equation over the streamwise direction gives

(4.5) \begin{equation} \frac {\partial \langle u^* \rangle _x}{\partial t} = -\langle v^* \rangle _x \frac {{\rm d} \overline {U}(y)}{{\rm d}y}. \end{equation}

Further, by following the same process for the vertical momentum equation, it can be shown that $\langle v^* \rangle _x$ is independent of time. Integrating (4.5) in time results in

(4.6) \begin{equation} \langle u^* \rangle _x = \langle u^*_0 \rangle _x - t\langle v^*_0 \rangle _x \frac {{\rm d} \overline {U}(y)}{{\rm d}y}, \end{equation}

where $\langle u^*_0 \rangle _x$ and $\langle v^*_0 \rangle _x$ are initial conditions. Equation (4.6) demonstrates that, in the inviscid limit, any parallel shear flow will exhibit an algebraic linear growth of kinetic energy, irrespective of the existence of an inflection point, given a finite streamwise-independent perturbation. This mechanism was later referred to as the lift-up effect, and is the dominant process in the generation of streamwise elongated streaks from low-energy streamwise vortices and their subsequent amplification. In viscous flows, the growth is transient in time and said growth is related to the non-normality of the Navier–Stokes operator (Schmid & Henningson Reference Schmid and Henningson2001).

The clear resemblance between (4.4) and (4.5) suggests that the perturbation energy grows primarily via the lift-up effect. Therefore, the perturbation field component $\Delta u_1$ , which is the main responsible for the growth of the overall field, grows via the algebraic instability described by the lift-up effect. Interestingly, this instability mechanism, which is known to be key in the transition to turbulence and indeed in the dynamics of turbulence, can be seen to also be the main component responsible for the observed long-term algebraic growth of perturbations in wall turbulence.