2.1. Governing equations

We consider the incompressible flow in a channel between two infinitely extended parallel no-slip walls. The flow is driven by a time-varying forcing (streamwise pressure gradient) such that the streamwise flow rate is held constant. The spanwise flow rate is also imposed constant and equal to zero. The streamwise, wall-normal and spanwise directions are denoted, respectively, by $x$ , $y$ and $z$ . Velocity components in these directions are denoted, respectively, by $u$ , $v$ and $w$ . Quantities without any superscript have been non-dimensionalised by the channel half-gap $h$ (such that $0\leqslant y \leqslant 2$ ) and the (imposed) channel bulk velocity $U_b$ ( $U_b=\int _0^2\int _0^{L_z} u \; {\textrm{d}}y{\textrm{d}}z / 2L_z$ ), where $L_z$ is the spanwise width of the channel. Quantities with a $+$ superscript have been non-dimensionalised by the viscous length scale $\delta _{v} = \nu /u_{\tau }$ and the friction velocity $u_{\tau }=\sqrt {\tau _w/\rho }$ , where $\rho$ is the fluid density, $\nu$ is the kinematic viscosity and $\tau _w$ is the (measured) mean wall shear stress.

The instantaneous flow is governed by the Navier–Stokes equations for incompressible flows

(2.1) \begin{equation} \begin{cases} \displaystyle \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u} = - \boldsymbol{\nabla}\kern-1pt p + \frac{1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u} + \boldsymbol{f}_{\kern-2pt b}, \\[1ex] \displaystyle \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}=0, \end{cases} \end{equation}

where $p$ is the hydrodynamic pressure, $ \textit{Re}=U_bh/\nu$ is the Reynolds number and $\boldsymbol{f}_{\kern-2pt b}$ is the time-varying forcing that keeps the flow rate constant.

Direct numerical simulation (DNS) of (2.1) performed in this study makes use of the channelflow 2.0 code (Gibson et al. Reference Gibson2021). It is based on Fourier decomposition in the streamwise ( $x$ ) and spanwise ( $z$ ) directions and Chebyshev collocation in the wall-normal direction ( $y$ ). Dealiasing in the wall-parallel directions is performed using the $2/3$ rule. A third-order-accurate backward finite-difference scheme is used for time discretisation. The simulations illustrated in figure 1 are performed in a domain of size $L_x \times L_y \times L_z = 250\times 2\times 150$ , in the streamwise, wall-normal and spanwise directions, respectively, in units of the channel half-gap. The collocation points (before dealiasing) are $1500\times 65\times 1800$ , respectively, leading to resolutions comparable to previous works (Shimizu & Manneville Reference Shimizu and Manneville2019; Kashyap et al. Reference Kashyap, Duguet and Dauchot2020; Kashyap et al. Reference Kashyap, Duguet and Dauchot2022). Equivalent resolutions are used for the other domains mentioned in the paper. The control parameter in our simulations is the Reynolds number based on the bulk velocity. However, results are presented as a function of the friction Reynolds number $ \textit{Re}_{\tau }=u_{\tau }h/\nu$ , which can be computed a posteriori, to ease comparison with previous works.

Figure 1. Contours of streamwise velocity fluctuations in the plane $y^+\approx 35$ from DNSs at (a) $ \textit{Re}_{\tau } = 98$ , (b) $ \textit{Re}_{\tau } = 92$ and (c) $ \textit{Re}_{\tau } = 71$ ( $y\approx 0.36$ for $ \textit{Re}_{\tau }=98,92$ and $y\approx 0.49$ for $ \textit{Re}_{\tau }=71$ ). The black vertical segment at the bottom left corner of each panel indicates a spanwise length of $1000$ wall units.

At the lowest Reynolds number considered here ( $ \textit{Re}_{\tau }=71$ ), the turbulent flow appears as a pattern of oblique bands (figure 1 c). The alternance of extended blue and red regions indicates the presence of the large-scale flow which coexists with the pattern (Duguet & Schlatter Reference Duguet and Schlatter2013). At a higher value of $ \textit{Re}_{\tau }=92$ , such a pattern of alternatively laminar and turbulent regions is replaced in DNS by a weaker modulation of a more standard streaky turbulent flow (figure 1 b). Instead, for $ \textit{Re}_{\tau }$ above 95, turbulence appears in DNS to be unmodulated by large-scale oblique structures (figure 1 a), consistently with Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022), and will be qualified as featureless after Shimizu & Manneville (Reference Shimizu and Manneville2019).

2.2. Streaks computation by resolvent analysis

As mentioned in the introduction, DNS visualisations suggest that the modulation of the turbulent channel flow, reported by Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022), exists on top of a field of predominantly streamwise streaks, see figure 1. The objective of this work is to investigate this phenomenology by performing a modal stability analysis of streaks superimposed on the mean shear flow. A wealth of previous studies have demonstrated that typical near-wall structures can be computed using the Navier–Stokes operator linearised around the mean flow either by considering the harmonic forcing problem (Hwang & Cossu Reference Hwang and Cossu2010b ; McKeon & Sharma Reference McKeon and Sharma2010; Moarref et al. Reference Moarref, Sharma, Tropp and McKeon2013; Sharma & McKeon Reference Sharma and McKeon2013) or the initial value problem (Del Alamo & Jimenez Reference Del Alamo and Jimenez2006; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009). The results of this analysis are briefly revised here at low $ \textit{Re}$ .

Using the Reynolds decomposition, the solution of (2.1) can be partitioned as the sum of the mean flow $\boldsymbol{\overline {U}} = \overline {U}(y)\boldsymbol{e}_x$ plus fluctuations. Then, following Reynolds & Hussain (Reference Reynolds and Hussain1972), the turbulent fluctuations are further partitioned into a coherent part and an incoherent part. Neglecting the interactions between the two parts, the nonlinear interaction of the coherent part with itself and modelling the incoherent fluctuations using a classical eddy viscosity model, one obtains the following linear system:

(2.2) \begin{align} \begin{cases} \displaystyle \frac{\partial \boldsymbol{u}^{\prime }}{\partial t} + \boldsymbol{\overline {U}}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}^{\prime } + \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\overline {U}} + \boldsymbol{\nabla}\kern-1pt p^{\prime } - \frac{1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u}^{\prime } - \boldsymbol{\nabla } \boldsymbol{\cdot }\big [ \nu _t\big ( \boldsymbol{\nabla } \boldsymbol{u}^{\prime } + (\boldsymbol{\nabla } \boldsymbol{u}^{\prime })^T \big ) \big ] = \boldsymbol{f}, \\[1ex] \displaystyle \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0, \end{cases} \end{align}

for the coherent velocity fluctuation $\boldsymbol{u}^{\prime }$ (respectively $p'$ for the coherent pressure fluctuation), where $\boldsymbol{f}$ is a generic forcing term (cf. Morra et al. (Reference Morra, Semeraro, Henningson and Cossu2019)). Here, $T$ denotes the transpose of the velocity gradient tensor. The coherent fluctuation is interpreted mainly as the large scales, but it also includes the streaks because of their temporal and spatial coherence. Smaller scales, as well as the remaining other flow structures, that by default we consider as ‘incoherent’, are not explicitly resolved, instead they are modelled as an eddy viscosity.

In this work, as in previous works on channel flow (Del Alamo & Jimenez Reference Del Alamo and Jimenez2006; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009; Hwang & Cossu Reference Hwang and Cossu2010b ; Park, Hwang & Cossu Reference Park, Hwang and Cossu2011; Alizard Reference Alizard2015; Ciola et al. Reference Ciola, De Palma, Robinet and Cherubini2024), the Cess (Reference Cess1958) eddy viscosity formula, as reported by Reynolds & Tiederman (Reference Reynolds and Tiederman1967), is employed

(2.3) \begin{align} \nu _t(y) = \frac{1}{2\textit{Re}}\left \{ \left [ 1 + \left ( \frac{\textit{Re}_{\tau }k}{3} \big ( 2y-y^2 \big ) \big ( 3-4y+2y^2 \big ) \big ( 1-e^{\left ( \lvert y-1 \rvert -1 \right )\textit{Re}_{\tau }/A} \big ) \right )^2 \right ]^{1/2} - 1 \! \right \}\!\!, \end{align}

with $k=0.426$ , $A=25.4$ (Del Alamo & Jimenez Reference Del Alamo and Jimenez2006) and $y \in [0,2]$ . Recent comparative studies have shown that including an eddy viscosity term in the linear operator is crucial for quantitatively improved predictions of the linear model (Illingworth, Monty & Marusic Reference Illingworth, Monty and Marusic2018; Morra et al. Reference Morra, Semeraro, Henningson and Cossu2019; Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023). In this study, we will compare stability results with and without eddy viscosity. Borrowing the terminology from Morra et al. (Reference Morra, Semeraro, Henningson and Cossu2019), (2.2) without the eddy viscosity term will be referred to as the quasi-laminar model, whereas the full equation will be referred to as eddy viscosity model. Note, however, that, in both cases, the eddy viscosity given by (2.3) is used to compute the mean profile solving the one-dimensional boundary value problem of Reynolds & Tiederman (Reference Reynolds and Tiederman1967).

The linear nature of the system (2.2) greatly simplifies its analysis: by considering the Fourier transform of the equation in the wall-parallel directions and in time, the harmonic response is computed as $\boldsymbol{\hat {u}}(y;k_x,k_z,\omega )=\mathcal{H}(y;k_x,k_z,\omega ) \boldsymbol{\hat {f}}(y;k_x,k_z,\omega )$ , where $\mathcal{H}(y;k_x,k_z,\omega )$ is the resolvent operator, the hats denote Fourier transforms in $x$ , $z$ and $t$ , $\omega$ is the angular frequency of both the forcing and the response and $k_x$ and $k_z$ are the streamwise and spanwise wavenumbers. For a given $\{\omega ,k_x,k_z\}$ triplet. The optimal harmonic forcing is defined as the function $\boldsymbol{\hat {f}}(y)$ that maximises the amplitude ratio between response and forcing (cf. Hwang & Cossu Reference Hwang and Cossu2010b )

(2.4) \begin{equation} R(k_x,k_z,\omega ) = \max _{\boldsymbol{\hat {f}}\neq \boldsymbol{0}} \frac{ \lVert \boldsymbol{\hat {u}} \rVert ^2} {\lVert \boldsymbol{\hat {f}} \rVert ^2}=\left \lVert \mathcal{H}(y;k_x,k_z,\omega ) \right \rVert ^2, \end{equation}

where the standard $L^2$ norm in $y$ is used both for $\boldsymbol{\hat {u}}$ and $\boldsymbol{\hat {f}}$ . In this study, the optimal harmonic forcing problem is solved using an in-house code which performs the singular value decomposition of the discretised resolvent operator (Schmid & Brandt Reference Schmid and Brandt2014; Symon et al. Reference Symon, Rosenberg, Dawson and McKeon2018). After Fourier transform in $x$ , $z$ and $t$ and discretisation in $y$ , the system (2.2) can be written as $ ( \iota \omega \unicode{x1D63D} - \unicode{x1D647} ){\hat {\unicode{x1D666}}}={\hat {\unicode{x1D65B} }}$ , where $\hat {\unicode{x1D666}}$ contains the four primitive variables $\hat {u}$ , $\hat {v}$ , $\hat {w}$ and $\hat {p}$ for each point of the computational grid, $ \unicode{x1D63D}$ is the prolongation matrix (Cerqueira & Sipp Reference Cerqueira and Sipp2014), $ \unicode{x1D647}$ is the discretised linear Navier–Stokes operator and $\iota$ the imaginary unit. The discretised resolvent operator can be computed by a matrix inversion as $ \unicode{x1D643}= ( \iota \omega \unicode{x1D63D} - \unicode{x1D647} )^{-1}$ and decomposed such that $ \unicode{x1D643}= \unicode{x1D650}\boldsymbol{\mathsf{\varSigma }} \unicode{x1D651}^T$ . The most amplified harmonic forcing is given by the first column of $ \unicode{x1D651}$ , the corresponding response mode by the first column of $ \unicode{x1D650}$ and the energy amplification factor $R$ by the square of the leading singular value of $\unicode{x1D643}$ , which is the first diagonal element of $\boldsymbol{\mathsf{\varSigma }}$ . For the discretisation of differential operators, $N_y=65$ Chebyshev collocation points have proven sufficient at the low values of $ \textit{Re}$ considered. The code was validated by reproducing the results of Hwang & Cossu (Reference Hwang and Cossu2010b ).

Resolvent analysis is not the main focus of this work, it is used here as a means of defining the streaks that are relevant for the stability analysis. Note that there are other alternative ways to define two-dimensional streak modes, including eigenmodes from the associated Stokes operator (Waleffe Reference Waleffe1997) or data-driven techniques where streaks are directly extracted from DNSs (Hack & Zaki Reference Hack and Zaki2014) or experimental data (Liu et al. Reference Liu, Semin, Godoy-Diana and Wesfreid2024). The present choice of the resolvent modes was motivated by the available machinery together with the need to make the results reproducible. Therefore, at this stage, we do not explore systematically the optimal harmonic response. It is well known that typical near-wall structures in the considered flow are streaks characterised by a dominant spanwise wavelength $\lambda ^+_{z,s} \approx 100$ , defined in wall units. Moreover, in a first approximation, they can be seen as streamwise-independent structures (straight streaks following Liu et al. Reference Liu, Semin, Godoy-Diana and Wesfreid2024). Streamwise streaks are routinely explained by the advection of streamwise vortices by the lift-up effect (Ellingsen & Palm Reference Ellingsen and Palm1975; Landahl Reference Landahl1980), known as a linear mechanism, although explaining their persistence in turbulent flows requires nonlinear concepts such as the SSP. As it happens, resolvent analysis with the operator linearised around the mean flow is well suited to computing such structures, including the one optimally amplified by linear mechanisms. In the resolvent framework, the forcing $\boldsymbol{f}$ represents the vortices and the response $\boldsymbol{u}^{\prime }$ the streak field. The optimal harmonic forcing problem is solved here for $k_x = 0$ and $k^+_z = 2\pi /100$ and using $20$ equi-spaced $\omega$ values in $[0,2\pi ]$ . Optimal amplification is found for steady forcing ( $\omega =0$ ) as in Hwang & Cossu (Reference Hwang and Cossu2010b ). The problem is solved for a range of $ \textit{Re}$ relevant to the modulational regime of turbulent patterns in the channel, in particular $ \textit{Re}_{\tau }=[71, 78, 84, 91, 95, 98, 105]$ .

Figure 2. Optimal harmonic forcing (a–b–c–d) and response velocity field (e–f–g–h) obtained from the resolvent analysis of the mean flow for (a–b–e–f) $ \textit{Re}_{\tau }=71$ and (c–d–g–h) $ \textit{Re}_{\tau }=105$ . In all panel contours denote the streamwise component while quivers stand for transverse components. The scale of the arrows in the top row is ten times larger than in the bottom row. Streamwise uniform case ( $k_x=0$ ). (a–c–e–g) Quasi-laminar model; (b–d–f–h) eddy viscosity model.

As anticipated, the forcing is found to be mainly directed in the wall-normal and spanwise directions, while the response is primarily dominated by the streamwise velocity component. Figure 2 shows the streamwise component of the forcing and response as shaded contours and the transverse components as arrows in the $y{-}z$ plane for two values of $ \textit{Re}_{\tau }$ . The typical picture of streamwise rolls and induced streaks is recognised here. Comparing the quasi-laminar model (panels a–c–e–g) with the eddy viscosity model (panels b–d–f–h), it is found that the eddy viscosity squeezes both streaks and vortices towards the wall. This effect increases with the Reynolds number. This is likely due to the wall-normal dependence of the eddy viscosity (2.3). We now turn towards the LSA of the streaks.

3.1. Formulation

The lift-up effect is only one of the several mechanisms involved in the SSP of wall turbulence (Hamilton et al. Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997). Another step in this cyclic process is the so-called secondary instability of the streaks that have been amplified by the lift-up effect (Schoppa & Hussain Reference Schoppa and Hussain2002; Park et al. Reference Park, Hwang and Cossu2011; Alizard Reference Alizard2015). In the classical SSP picture, that stage precedes the feedback onto the streamwise vortices. More recently, de Giovanetti et al. (Reference de, Matteo, Jin and Hwang2017) and Ciola et al. (Reference Ciola, De Palma, Robinet and Cherubini2024) have shown that streak instabilities can also contribute to the generation of larger-scale structures, thereby allowing for an escape from the one-scale-only process described by the original SSP theory. At lower $ \textit{Re}$ typical of the transitional regime, a similar analysis is of interest to investigate whether the secondary instability of streaks can be related to the emergence of large-scale modulations of the turbulent flow suggested from DNS observations.

As illustrated in the previous section, the resolvent analysis identifies an optimal response with zero frequency. The velocity field associated with this response is hereafter denoted by $\boldsymbol{u}_s$ . Due to the linearity of (2.2), $\boldsymbol{u}_s$ is defined up to a multiplicative constant. Following our previous work (Ciola et al. Reference Ciola, De Palma, Robinet and Cherubini2024), $\boldsymbol{u}_s$ is rescaled such that

(3.1) \begin{equation} \frac{\max _{y,z} u_s - \min _{y,z} u_s}{2\overline {U}_{c}} = 1, \end{equation}

where $u_s$ (as a scalar field) is the streamwise velocity component of the streaks field and $\overline{U}_{\kern-1pt c}$ is the centreline velocity of the turbulent mean profile.

A new base flow $\boldsymbol{U}$ is constructed from the knowledge of the mean flow $\boldsymbol{\overline {U}}$ and the streak field $\boldsymbol{u}_s$ as

(3.2) \begin{equation} \boldsymbol{U} = \boldsymbol{\overline {U}} + A_s \boldsymbol{u}_s, \end{equation}

with $A_s\gt 0$ interpreted as the streak amplitude. Consistently with the previous section, $\boldsymbol{U}=(U(y,z),V(y,z),W(y,z))$ is uniform in the streamwise direction, which simplifies and reduces the computational cost of LSA compared with the case of a genuinely three-dimensional base flow.

Since both the mean flow and the optimal response are steady, we expect that the base flow $\boldsymbol{U}$ verifies an equation of the generic form

(3.3) \begin{equation} \begin{cases} \displaystyle \boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{U} + \boldsymbol{\nabla } P - \frac{1}{\textit{Re}} {\nabla} ^2 \boldsymbol{U} - \boldsymbol{\nabla } \boldsymbol{\cdot }\big [ \nu _t\big ( \boldsymbol{\nabla } \boldsymbol{U} + (\boldsymbol{\nabla } \boldsymbol{U})^T \big ) \big ] - \boldsymbol{f}_{\kern-2pt b} = \boldsymbol{g}, \\[1ex] \displaystyle \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{U}=0, \end{cases} \end{equation}

with some forcing term $\boldsymbol{g}$ to be determined. Such a forcing term is necessary since $\boldsymbol{U}$ is not a steady-state solution to the unforced Navier–Stokes equations. This forcing term  $\boldsymbol{g}$ differs from the forcing term $\boldsymbol{f}$ in (2.2) because (2.2) is linear with respect to $\boldsymbol{u}_s$ whereas (3.3) implicitly contains the nonlinear term $A_s^2\boldsymbol{u}_s\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}_s$ ; with $A_s$ not vanishing.  $\boldsymbol{g}$ can be computed explicitly using a nonlinear time-stepping algorithm as detailed in Appendix A.

Once $\boldsymbol{g}$ is known, the following nonlinear system is fully specified by:

(3.4) \begin{equation} \begin{cases} \displaystyle \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u} + \boldsymbol{\nabla}\kern-1pt p - \frac{1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u} - \boldsymbol{\nabla } \boldsymbol{\cdot }\big [ \nu _t\big ( \boldsymbol{\nabla } \boldsymbol{u} + (\boldsymbol{\nabla } \boldsymbol{u})^T \big ) \big ] - \boldsymbol{f}_{\kern-2pt b} = \boldsymbol{g}, \\[1ex] \displaystyle \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}=0. \end{cases} \end{equation}

By construction, the base flow $\boldsymbol{U}$ is a steady-state solution of this system. Moreover, when the amplitude of the streaks contained in $\boldsymbol{U}$ is sufficiently small, the forcing $\boldsymbol{g}$ converges towards the optimal harmonic forcing obtained from resolvent analysis $\boldsymbol{f}$ , as verified in Appendix A.

Since $\boldsymbol{U}$ is now a steady solution of (3.4), LSA can be used rigorously. By perturbing the base flow $\boldsymbol{U}$ with a small perturbation $\boldsymbol{u}^{\prime \prime }$ such that $\boldsymbol{u}=\boldsymbol{U}+\boldsymbol{u}^{\prime \prime }$ and neglecting quadratic perturbation terms, one gets the following linear system for $\boldsymbol{u}^{\prime \prime }$ and the corresponding pressure perturbation $p^{\prime \prime }$ :

(3.5) \begin{equation} \begin{cases} \displaystyle \frac{\partial \boldsymbol{u}^{\prime \prime }}{\partial t} = -\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}^{\prime \prime } - \boldsymbol{u}^{\prime \prime }\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{U} - \boldsymbol{\nabla}\kern-1pt p^{\prime \prime } + \frac{1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u}^{\prime \prime } + \boldsymbol{\nabla } \boldsymbol{\cdot }\big [ \nu _t\big ( \boldsymbol{\nabla } \boldsymbol{u}^{\prime \prime } + (\boldsymbol{\nabla } \boldsymbol{u}^{\prime \prime })^T \big ) \big ], \\[1ex] \displaystyle \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}^{\prime \prime }=0. \end{cases} \end{equation}

The forcing term $\boldsymbol{g}$ cancels out in the secondary perturbation equation because it is constant, i.e. it depends only on $\boldsymbol{U}$ and not on $\boldsymbol{u}^{\prime \prime }$ . Similarly, $\boldsymbol{f}_{\kern-2pt b}$ cancels out because we assume that $\boldsymbol{u}^{\prime \prime }$ does not modify the flow rate. Again, one can either choose to include the eddy viscosity term or to neglect it. We will compare the two approaches, consistently with the analysis in § 2. Note that two different frozen eddy viscosity models, provided they predict exactly the same mean flow $\boldsymbol{U}$ and the same viscosity $\nu _t$ , are expected to yield the same stability results for $\boldsymbol{u}^{\prime \prime }$ .

Modal stability analysis assumes the following perturbation ansatz ( $\iota$ denotes the imaginary unit):

(3.6) \begin{equation} \displaystyle \boldsymbol{u}^{\prime \prime } = \tilde {\boldsymbol{u}}(y,z) e^{\sigma t + \iota k_x x} + \text{complex conjugate}, \quad \sigma \in \mathbb{C} \text{ and } k_x \in \mathbb{R}, \end{equation}

and an equivalent one for the pressure $p^{\prime \prime }$ . Inserting the ansatz (3.6) into (3.5), the following eigenvalue problem is obtained:

(3.7) \begin{equation} \begin{cases} \displaystyle \sigma \tilde {u} = - \iota k_xU\tilde {u} - V\frac{\partial \tilde {u}}{\partial y} - W\frac{\partial \tilde {u}}{\partial z} - \tilde {v}\frac{\partial U}{\partial y} - \tilde {w}\frac{\partial U}{\partial z} - \iota k_x\tilde {p} + \left ( \frac{1}{\textit{Re}}+\nu _t \right ){\nabla} ^2 \tilde {u} + \frac{{\textrm{d}}\nu _t}{{\textrm{d}}y}\left ( \frac{\partial \tilde {u}}{\partial y} + \iota k_x\tilde {v} \right )\!, \\[2ex]\displaystyle \sigma \tilde {v} = - \iota k_xU\tilde {v} - V\frac{\partial \tilde {v}}{\partial y} - W\frac{\partial \tilde {v}}{\partial z} - \tilde {v}\frac{\partial V}{\partial y} - \tilde {w}\frac{\partial V}{\partial z} - \frac{\partial \tilde {p}}{\partial y} + \left ( \frac{1}{\textit{Re}}+\nu _t \right ){\nabla} ^2 \tilde {v} + 2\frac{{\textrm{d}}\nu _t}{{\textrm{d}}y}\frac{\partial \tilde {v}}{\partial y}, \\[2ex]\displaystyle \sigma \tilde {w} = - \iota k_xU\tilde {w} - V\frac{\partial \tilde {w}}{\partial y} - W\frac{\partial \tilde {w}}{\partial z} - \tilde {v}\frac{\partial W}{\partial y} - \tilde {w}\frac{\partial W}{\partial z} - \frac{\partial \tilde {p}}{\partial z} + \left ( \frac{1}{\textit{Re}}+\nu _t \right ){\nabla} ^2 \tilde {w} + \frac{{\textrm{d}}\nu _t}{{\textrm{d}}y}\left ( \frac{\partial \tilde {w}}{\partial y} + \frac{\partial \tilde {v}}{\partial z} \right )\!, \\[2ex] 0 = \displaystyle \iota k_x \tilde {u} + \frac{\partial \tilde {v}}{\partial y} + \frac{\partial \tilde {w}}{\partial z}, \end{cases} \end{equation}

with ${\nabla} ^2 =-k_x^2 + \partial ^2_y + \partial ^2_z$ .

The secondary stability problem is solved in a two-dimensional $y{-}z$ domain using Fourier collocation in the spanwise direction and Chebyshev collocation in the wall-normal direction.

3.2. Exploiting the spanwise periodicity of the streaks

Standard LSA is usually performed by assuming eigenvectors with a wavelength matching the fundamental wavelength of the base flow. This refers here to the spanwise wavelength $\lambda _{z,s}$ of the streaks. By contrast, the wavelength can be freely chosen in the $x$ -direction. In that case, the linearised system (3.7), once discretised, leads to a generalised eigenvalue problem of the form

(3.8) \begin{equation} \sigma \unicode{x1D63D}_0{\tilde {\unicode{x1D666}}}_0 = \unicode{x1D63C}_0{\tilde {\unicode{x1D666}}}_0, \end{equation}

with ${\tilde {\unicode{x1D666}}}_0$ a vector of size $M$ , containing the four primitive variables $\tilde {u}$ , $\tilde {v}$ , $\tilde {w}$ and $\tilde {p}$ for each point of the computational grid, and $ \unicode{x1D63C}_0$ and $ \unicode{x1D63D}_0$ are matrices of size $M\times M$ . The matrix $ \unicode{x1D63D}_0$ is the so-called prolongation matrix that arises when using a velocity–pressure formulation (Cerqueira & Sipp Reference Cerqueira and Sipp2014).

We are interested in this study in predicting instabilities of the base flow occurring over wavelengths larger than that of the base flow. This can be achieved by simply considering a numerical domain containing several concatenated copies of the spatially periodic base flow. However, solving (3.8) becomes computationally more costly as the domain size increases. This cost can be reduced exploiting the $z$ -periodicity of the base flow by following the study of Schmid et al. (Reference Schmid, De, Fosas and Peake2017). We therefore consider a given (integer) number $N_u$ of base flow units concatenated in the $z$ -direction, such that the total spanwise size of the domain used for LSA is now $N_u\lambda _{z,s}$ . After spatial discretisation, the stability problem (3.7) becomes a new generalised eigenvalue problem of the form

(3.9) \begin{equation} \sigma \unicode{x1D63D}{\tilde {\unicode{x1D666}}} = \unicode{x1D63C}{\tilde {\unicode{x1D666}}}, \end{equation}

where $\tilde {\unicode{x1D666}}$ has $N_u$ times the size of the original ${\tilde {\unicode{x1D666}}}_0$ vector in (3.8), i.e. $\textit{MN}_u$ , and the matrices $A$ and $B$ have now size $\textit{MN}_u \times \textit{MN}_u$ . Owing to the periodicity of the base flow and to the assumed periodicity of the eigenvector over $N_u$ unit cells, the matrix $ \unicode{x1D63C}$ has a block-circulant structure. One can hence diagonalise it into a block-diagonal matrix ${\kern2pt \hat{\kern-2pt\unicode{x1D63C}}}={ \unicode{x1D64B}\,}^H \unicode{x1D63C} \unicode{x1D64B}$ with $N_u$ blocks ${\kern2pt \hat{\kern-2pt\unicode{x1D63C}}}_0,{\kern2pt \hat{\kern-2pt\unicode{x1D63C}}}_1,\ldots ,{\kern2pt \hat{\kern-2pt\unicode{x1D63C}}}_{N_u}$ , thereby reducing one large eigenvalue problem of size $N_u M \times N_uM$ down to $N_u$ manageable eigenproblems, each of size $M \times M$ . The matrix $ \unicode{x1D64B}\in \mathbb{C}^{N_uM\times N_uM}$ is given by the Kronecker product $ \unicode{x1D645}\otimes \unicode{x1D644}_{M}$ with $ \unicode{x1D645}\in \mathbb{C}^{N_u\times N_u}$ , $J_{j+1,k+1}=\rho _j^k/\sqrt {N_u}$ and $ \unicode{x1D644}_M$ the identity matrix of size $M\times M$ (Schmid et al. Reference Schmid, De, Fosas and Peake2017). Each of these spectral subproblems involves a matrix of size $M$ , namely ${\kern2pt \hat{\kern-2pt\unicode{x1D63C}}}_j= \unicode{x1D63C}_0+\rho _j \unicode{x1D63C}_1 + \rho _j^2 \unicode{x1D63C}_2 + \cdots + \rho _j^{N_u-1} \unicode{x1D63C}_{N_u-1}$ for each $j=0,\ldots ,N_u-1$ , featuring $ \unicode{x1D63C}_0, \unicode{x1D63C}_1,\ldots$ which are the blocks of $ \unicode{x1D63C}$ . The general eigensolutions of (3.9) can be written (Schmid et al. Reference Schmid, De, Fosas and Peake2017) as

(3.10) \begin{equation} {\tilde {\unicode{x1D666}}}=\left ({\tilde {\unicode{x1D666}}}_j,\rho _j{\tilde {\unicode{x1D666}}}_j,\ldots ,\rho _j^{N_u-1}{\tilde {\unicode{x1D666}}}_j\right )^T \! , \end{equation}

based on the known eigenvectors ${\tilde {\unicode{x1D666}}}_j$ of the matrices ${\kern2pt \hat{\kern-2pt\unicode{x1D63C}}}_j$ .

In (3.10), the complex numbers $\rho _j$ , $j=1,\ldots ,N_u$ , are the $N_u$ different roots of unity solutions of the equation $\rho ^{N_u}=1$ . This defining property for the $\rho _j$ values ensures that the generalised eigenvector $\tilde {\unicode{x1D666}}$ is periodic in $z$ over the distance $N_u\lambda _{z,s}$ . This allows one to use spectral methods based on periodic Fourier functions. Each $\rho _j$ has modulus one and can be rewritten as $\rho _j=\exp (2\pi \gamma _j \iota )$ . The $\gamma$ values are defined in $[0,1)$ . They take only discrete values $0,1/N_u,2/N_u,\ldots$ These exponents are ideally interpreted, for a given eigenvector, as the ratio between the wavelength of the base flow and that of the eigenvector. In the particular case where $\gamma =0$ ( $\rho _1=1$ ), the eigenvector has the same periodicity as the base flow, i.e. $N_u=1$ . One recovers thus the classical LSA results, since $\rho _1=1$ and $\gamma =0$ . In the special case where $\gamma$ is of the form $1/Q$ with $Q$ a non-zero integer, then the fundamental wavelength of the eigenvector is exactly $Q\lambda _{z,s}$ , and it can be labelled sub-harmonic. In the general case, for finite $N_u$ , only rational values of $\gamma$ can be tackled using this method and irrational values are excluded. The physical range of meaningful values of  $\gamma$ is continuous and involves irrational values as well. In practice, it is advised to consider a value of $N_u$ as large as possible in order to have access to a well-discretised range of values of $\gamma$ in $[0,1)$ . As a consequence, $\gamma$ can be interpreted as a detuning factor (Jouin et al. Reference Jouin, Ciola, Cherubini and Robinet2024). Note moreover that the symmetry $z \leftarrow -z$ inherent to the base flow implies that $\gamma$ and $1-\gamma$ are associated with the same eigenvalue and eigenvector. This restricts the study of the eigenspectrum to the range of values $\gamma \in [0,0.5]$ . Finally, the union of the spectra of the sub-matrices ${\kern2pt \hat{\kern-2pt\unicode{x1D63C}}}_j$ yields the spectrum of the original matrix $ \unicode{x1D63C}$ . In practice, the latter is found by looping over $j=0,\ldots ,N_u-1$ , which is indirectly a loop over the values of $\gamma$ in $[0,0.5]$ . This is the basis of the computation shown in figure 3.

The exponent $\gamma$ plays the same role of the Floquet modulation parameter in spatial Floquet theory (Herbert Reference Herbert1988), despite the fact that the block-circulant matrix method and spatial Floquet analysis (also called Bloch formalism) are formally different. However, a posteriori the two methods give identical results (Jouin et al. Reference Jouin, Ciola, Cherubini and Robinet2024).

The main results of the paper are obtained using $N_u=50$ base flow units. Each unit was discretised with $65$ points in the wall-normal direction and $20$ collocation points in the spanwise direction. The number of points in the spanwise direction refers to the discretisation of a single streak wavelength. Therefore, it turns out to be sufficient for the relatively narrow spanwise width of 100 wall units considered here. The convergence of the results with the number of collocation points and with the number of base flow units is assessed in Appendix B. The code was validated in a previous study (Ciola et al. Reference Ciola, De Palma, Robinet and Cherubini2024).

3.3. Critical Reynolds number

The first question addressed in this work is whether a critical Reynolds number exists below which spatial modulations appear. In the work of Kashyap et al. (Reference Kashyap, Duguet and Dauchot2024), based on a linearisation around the asymptotic mean flow (e.g. without streaks), no such critical threshold was identified. It is well established that streaks become unstable for sufficiently large amplitude (Schoppa & Hussain Reference Schoppa and Hussain2002). However, classical streak instability is characterised by wavelengths which are either identical to ( $\gamma =1$ ) or double ( $\gamma =0.5$ ) the wavelength of the streak (Andersson et al. Reference Andersson, Brandt, Bottaro and Henningson2001), whereas the focus of this work is on instabilities that modulate the streaky flow over much larger scales. Therefore, the second point to be addressed is whether the unstable mode at criticality features an important large-scale component. In this subsection we focus specifically on $k_x=0.18$ , the critical streamwise wavenumber reported by Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022).

Figure 3. Eigenvalues for the streak stability problem (only a subset of the computed spectra is shown) for $ \textit{Re}_{\tau }=71$ and $k_x=0.18$ . The eigenvalues ( $\bullet$ for stable modes, $\star$ for unstable modes) are coloured with the corresponding root of unity factor $\gamma = j/N_u$ for $j=0,\ldots ,N_u/2$ (the eigenvalues for $\gamma \in (0.5,1.0)$ are equal to those for $\gamma \in (0,0.5)$ and the corresponding modes are the same up to a reflection in the spanwise direction). (a–c–e) Quasi-laminar model; (b–d–f) eddy viscosity model. Streak amplitudes increase from top to bottom and are (a) $0.10$ , (c) $0.16$ and (e) $0.20$ for the quasi-laminar model and (b) $0.30$ , (d) $0.40$ and (f) $0.50$ for the eddy viscosity model.

The effect of the streak amplitude $A_s$ on the eigenvalues associated with (3.9) is documented in figure 3. This figure shows a close-up view in the complex plane of the leading branch of eigenvalues, i.e. the least stable (or the most unstable) one. It shows the eigenvalues in the form of complex phase velocities $c=-\sigma /\iota k_x$ . The figure displays only a small part of the computed eigenspectrum. There is no instability elsewhere in the range of parameters investigated. As explained in the previous section, these spectra are the union of $N_u$ sub-spectra, each associated with a different root of unity. Therefore, in figure 3 the eigenvalues are coloured according to the respective factor $\gamma$ , which identifies the relevant root of unity and takes values in $[0,1)$ . In practice, only the range $[0,0.5]$ needs to be shown since the eigenmodes corresponding to $\gamma$ in $(0.5,1)$ are obtained from those in $(0,0.5)$ by spanwise reflection (see § 3.2). Focusing on $ \textit{Re}_{\tau }=71$ and $k_x=0.18$ , the figure shows that the branch becomes unstable as the $A_s$ is increased. The critical amplitude for the quasi-laminar model (left panels) is ${\approx} 0.16$ , whereas for the eddy viscosity model it is larger, ${\approx} 0.4$ (right panels). After inspection of the eigenvalues for several $A_s$ , we report that, for both models, the first mode that becomes unstable at $k_x=0.18$ has $\gamma \gt 0$ .

The $ \textit{Re}_{\tau }$ dependence of the leading eigenvalue is documented in figure 4 for varying $A_s$ . The first row shows the leading growth rates for the two models considered, namely quasi-laminar or using an eddy viscosity. In the quasi-laminar model the growth rates become positive without any apparent dependence on $ \textit{Re}$ . A critical amplitude can still be defined, independently of $ \textit{Re}_\tau$ , but no critical Reynolds number exists. For the eddy viscosity model, the growth rate increases with $A_s$ and depends also on the value of $ \textit{Re}_{\tau }$ . As a result, for $A_s$ large enough, as $ \textit{Re}_{\tau }$ is decreased, the growth rate of the most unstable mode becomes positive at a well-defined critical Reynolds number. A striking point is the qualitative difference between the results of the two models. The presence of a critical $ \textit{Re}$ with the eddy viscosity model is a qualitative improvement with respect to the mean flow analysis in Kashyap et al. (Reference Kashyap, Duguet and Dauchot2024).

The bottom panels of figure 4 show the phase velocity of the leading modes with respect to $ \textit{Re}_\tau$ and $A_s$ . The values globally overestimate by roughly $5\,\%{-}10\,\%$ the advection velocities of turbulent bands measured from DNS (Tuckerman et al. Reference Tuckerman, Kreilos, Schrobsdorff, Schneider and Gibson2014), however, they display a similar decrease with $ \textit{Re}$ .

Besides the fact that the critical value itself depends on the value of $A_s$ , another point of concern is the fact that the instability and the Reynolds-number dependence are observed only for large amplitudes. The issue whether these amplitudes are realistic will be addressed in the discussion section.

Figure 4. Variation of the leading eigenvalues with $ \textit{Re}$ and streak amplitude $A_s$ for $k_x=0.18$ ( $\bullet$ for stable modes, $\star$ for unstable modes). (a–b) Growth rate; (c–d) phase velocity. (a–c) Quasi-laminar model; (b–d) eddy viscosity model.

3.4. Unstable (large-scale) modes

The above results are encouraging and call for an examination of the unstable eigenmodes. Notably, no spanwise wavelength of the modes is explicitly imposed in the ansatz in (3.6). All detuned eigenmodes contain energy in wavelengths larger than the base flow wavelength $\lambda _{z,s}$ , yet the question arises whether some eigenmodes have a particularly prominent large-scale component. In this respect, the detuning factor $\gamma$ yields only partial information. Therefore, an additional quantification of the large-scale property is introduced by

(3.11) \begin{align} r_{\textit{LS}} = \dfrac {\int _0^2 \; \sum \limits _{\lvert k_z \rvert \lt {k_z^c}} \lvert \check {\boldsymbol{u}}(y,k_z) \rvert ^2 \; {\textrm{d}}y}{\int _0^2 \; \sum \limits _{k_z} \lvert \check {\boldsymbol{u}}(y,k_z) \rvert ^2 \; {\textrm{d}}y}, \end{align}

where $\check {\boldsymbol{u}}$ denotes the Fourier transform of the mode in $z$ and $k_z^c$ is a cutoff wavenumber that represents a higher limit for a Fourier mode to be considered large scale. The value $k_z^c=0.5$ was chosen after inspection of the DNS and it corresponds to a wavelength which is roughly $10$ times greater than the wavelength of the streaks depending on $ \textit{Re}$ . The new quantity $r_{\textit{LS}}$ is the ratio of the eigenmode’s energy at large wavelengths to the total eigenmode energy. It is used to colour the eigenvalues in e.g. figure 5, where darker points correspond to eigenmodes with a pronounced large-scale character. Again, the qualitative difference between the quasi-laminar model and the eddy viscosity model is significant. For the eddy viscosity model the leading mode has a large-scale factor near $10\,\%$ , to be compared with much weaker values for the quasi-laminar model.

Figure 5. Eigenvalues coloured by the large-scale spanwise energy ratio $r_{\textit{LS}}$ (only a subset of the computed spectra is shown) for $ \textit{Re}_{\tau }=71$ and $k_x=0.18$ ( $\bullet$ for stable modes, $\star$ for unstable modes). (a) Quasi-laminar model and $A_s=0.2$ ; (b) eddy viscosity model and $A_s=0.5$ .

The spatial structure of this interesting eigenmode for the eddy viscosity model is shown in figure 6 in the $y{-}z$ plane. The eigenmode is characterised by a small-scale modulation with the same spanwise periodicity as the base flow streaks, and by a large-scale modulation. In particular, the large-scale modulation of the wall-normal velocity component appears as an amplitude modulation. On the spanwise component, it appears as a large-scale flow. As for the streamwise component, the large-scale character is less clear but the large-scale modulating feature is still clearly observed. These qualitative observations are confirmed looking at the Fourier decomposition of the modes along $z$ . The real part of each eigenmode component $\hat {u}_i(y,z)$ is decomposed in Fourier series. The wall-normal integrated squared modulus of the Fourier coefficients gives energies as a function of spanwise wavenumber $E_i(k_z)$ , which are plotted in figure 7 for the three velocity components. The figure shows that the eigenmode can be seen as a superposition of waves. The first wave has a small wavenumber and is responsible for the large-scale character of the eigenmode. We note that $r_{\textit{LS}}$ was designed to give the ratio of the energy of this large-scale component with respect to the total energy of the eigenmode. The other waves have wavelengths near $\lambda _{z,s}$ (dashed line in the figure) or smaller. Moreover, figure 7(c) shows that the large-scale wave is dominant in the spanwise velocity component, as suggested by figure 7(c). Lastly, we note that the large-scale spanwise wavenumber of this mode is $k_z \approx 0.36$ . Given that $k_x=0.18$ for this eigenmode, it means that the modulation forms an angle of $\tan ^{-1}(k_x/k_z)\approx 26.6^{\circ }$ with respect to the streamwise direction. This angle is not far from the $23^{\circ }$ reported by Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022) and the $22.5^{\circ }$ found by Benavides & Barkley (Reference Benavides and Barkley2025).

Figure 6. Leading eigenmode for the eddy viscosity model with $ \textit{Re}_{\tau }=71$ , $k_x=0.18$ and $A_s=0.50$ ( $\lambda _{z,s}\approx 1.41$ at this value of $ \textit{Re}_{\tau }$ ). The modes are normalised to $\max _{y,z} \lvert \textit{Re}(\tilde {u}) \rvert = 1$ . Real part of (a) streamwise velocity component, (b) wall-normal velocity component and (c) spanwise velocity component.

Figure 7. Spanwise Fourier decomposition of the leading eigenmode for the same parameters as figure 6. The figure shows the wall-normal integrated squared Fourier coefficient of the real part of the (a) streamwise, (b) wall-normal and (c) spanwise velocity components. The black dashed line denotes the wavenumber of the base flow $2\pi /\lambda _{z,s}$ .

The above analysis is consistent with DNS observations interpreted as a modulation of the small scales together with a wall-parallel large-scale flow. This ideal picture is confirmed by figure 8, which shows the eigenmode integrated along the wall-normal direction $y$ in the horizontal plane $x{-}z$ . The arrows portray the streamwise and spanwise velocity components forming the large-scale flow. The displayed velocity field is fully consistent with most observations in the patterning regime (Duguet & Schlatter Reference Duguet and Schlatter2013; Tuckerman et al. Reference Tuckerman, Kreilos, Schrobsdorff, Schneider and Gibson2014). Similar observations also apply to other unstable eigenmodes that belong to the branch and are detuned. There are also unstable modes which do not show any large-scale modulation, related to classical sinuous/varicose streak instabilities which are involved in the self-sustaining cycle (Hamilton et al. Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997). Their presence is expected from the classical literature on streak instabilities and they are hence not the focus of this work.

When the streamwise wavenumber is fixed to a value consistent with the value from Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022), namely $k_x=0.18$ , the eddy viscosity model gives a group of eigenmodes characterised by a large-scale modulation that becomes unstable as $ \textit{Re}_{\tau }$ is lowered. It is now necessary to assess which streamwise wavenumber is selected by the system. Figure 9 shows the dependence of the leading growth rate on the streamwise wavenumber for several values of $A_s$ and for all $ \textit{Re}$ under consideration. For low $A_s$ , the growth rate monotonically decreases with $k_x$ both for the quasi-laminar model and for the eddy viscosity model. The same behaviour was reported in Kashyap et al. (Reference Kashyap, Duguet and Dauchot2024) for the stability analysis of the mean flow, i.e. the case $A_s=0$ . As the streak amplitude is increased, the growth rate is maximal for a non-zero streamwise wavenumber. For the quasi-laminar model, the maximum is unique and the curves for the different $ \textit{Re}_{\tau }$ overlap. Moreover, the maximum occurs for $k_x \approx 1$ , which implies $\lambda _x \approx $ 3–6. These wavelengths are too short to be representative of the large-scale patterns of interest. Instead, with the eddy viscosity model, the instability is limited to $k_x \approx \mathcal{O}(0.1)$ , i.e. $\lambda _x \approx 30{-}40$ , which is consistent with expectations. In conclusion, only when the Reynolds stresses are taken into account does stability analysis select a relevant large-scale streamwise wavenumber.

Figure 8. Leading eigenmode (real part of the full ansatz) for the same parameters as figure 6. Wall-normal velocity component (shaded contours) at midplane and wall-parallel large-scale flow (black arrows) in the $x{-}z$ plane. The large-scale flow is obtained by integrating the wall-parallel velocity components in the wall-normal direction.

Figure 9. Variation of the leading growth rate $\sigma _r$ with $ \textit{Re}$ and streamwise wavenumber $k_x$ ( $\bullet$ for stable modes, $\star$ for unstable modes). (a–c–e) Quasi-laminar model; (b–d–f) eddy viscosity model. Note the different scale in $k_x$ between left and right panels. Streak amplitudes $A_s$ increase from top to bottom and are (a) $0.10$ , (c) $0.20$ and (e) $0.30$ for the no-closure cases and (b) $0.30$ , (d) $0.40$ and (f) $0.50$ for the cases with eddy viscosity.

3.5. Critical amplitude and wavelengths

This subsection is devoted to a deeper analysis of the eddy viscosity model at critical conditions. The critical parameters are defined as the parameters for which the leading mode has zero growth rate. Usually in linear stability studies, as in Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022), a critical $ \textit{Re}$ and critical wavenumbers are reported. In our model, we have a degree of freedom given by the streak amplitude. One can either define a critical $ \textit{Re}$ for a given $A_s$ (note that the critical $ \textit{Re}$ is defined only for sufficiently high $A_s$ , figure 4 b) or define a critical $A_s$ for each $ \textit{Re}$ , so that a neutral curve $ \textit{Re}{-}A_s$ can be computed. For each point on this neutral curve there is a critical eigenmode with a streamwise wavenumber $k_x^c$ . Moreover, a critical spanwise wavenumber $k_z^c$ can also be defined since the eigenmodes are characterised by one single large-scale wavelength (see figure 7).

The neutral curve is computed for the considered value of $ \textit{Re}_{\tau }$ by coupling a line search algorithm in $k_x$ with a bisection algorithm for $A_s$ . In the inner loop we look for the $k_x$ giving the maximum growth rate for a given $A_s$ . Then, $A_s$ is bisected until the absolute value of the maximum growth rate falls below a tolerance of $10^{-5}$ . The result is presented in figure 10(a). It can be seen that the critical amplitude grows almost linearly with $ \textit{Re}$ . This means that, for fixed $A_s$ , the instability is obtained by decreasing $ \textit{Re}$ , as explained in § 3.3 and consistently with Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022). The neutral curve has not been computed for the quasi-laminar model since the results presented in figure 4 already suggest that the critical amplitude is independent of $ \textit{Re}$ in the considered range.

Panels (b) and (c) of figure 10 show the critical wavenumbers. They increase with $ \textit{Re}$ , indicating that the instability tends towards smaller wavelengths at high $ \textit{Re}$ . This is also consistent with DNS observations (Kashyap et al. Reference Kashyap, Duguet and Dauchot2020). At $ \textit{Re}_{\tau }=95$ , we obtain $k_x^c\approx 0.3$ and $k_z^c\approx 0.71$ , which are slightly greater than the $k_x^c\approx 0.18$ and $k_z^c\approx 0.42$ of Kashyap et al. (Reference Kashyap, Duguet and Dauchot2022). However, the ratio between the wavenumbers, hence the angle of the eigenmode with the streamwise direction, is approximately the same ( ${\approx} 23^{\circ }$ ). We report that our critical angle slowly decreases with $ \textit{Re}$ from $25.5^{\circ }$ at $ \textit{Re}_{\tau }=71$ to $20.0^{\circ }$ at $ \textit{Re}_{\tau }=106$ (not shown).

Figure 10. Critical $A_s$ , $k_x$ and large-scale $k_z$ as a function of $ \textit{Re}_{\tau }$ for the eddy viscosity model.

4.1. Derivation

By taking the scalar product of $\tilde {\boldsymbol{u}}^*$ with the momentum equation in (3.7), one obtains the following scalar equation:

(4.1) \begin{equation} \begin{aligned} \sigma \tilde {u}_i^*\tilde {u}_i = & -\iota k_x U \tilde {u}_i^*\tilde {u}_i - U_{\kern-1pt j}\tilde {u}_i^*\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} - \tilde {u}_i^*\tilde {u}_{\kern-1pt j}\frac{\partial U_i}{\partial x_{\kern-1pt j}} -\iota k_x\tilde {u}^*\tilde {p} - \tilde {u}_{\kern-1pt j}^*\frac{\partial \tilde {p}}{\partial x_{\kern-1pt j}} - k_x^2 \left ( \frac{1}{\textit{Re}} + \nu _t \right ) \tilde {u}_i^*\tilde {u}_i \\[1ex] & + \left ( \frac{1}{\textit{Re}} + \nu _t \right )\tilde {u}_i^*\frac{\partial ^2 \tilde {u}_i}{\partial x_{\kern-1pt j} \partial x_{\kern-1pt j}} + \tilde {u}_i^*\frac{{\textrm{d}}\nu _t}{{\textrm{d}}y}\frac{\partial \tilde {u}_i}{\partial y} + \iota k_x \tilde {u}^*\tilde {v}\frac{{\textrm{d}}\nu _t}{{\textrm{d}}y} + \tilde {u}_{\kern-1pt j}^*\frac{{\textrm{d}}\nu _t}{{\textrm{d}}y}\frac{\partial \tilde {v}}{\partial x_{\kern-1pt j}}. \end{aligned} \end{equation}

Moreover, the following identities are considered:

(4.2) \begin{align} -\iota k_x\tilde {u}^*\tilde {p} - \tilde {u}_{\kern-1pt j}^*\frac{\partial \tilde {p}}{\partial x_{\kern-1pt j}} & = - \frac{\partial }{\partial x_{\kern-1pt j}}\left ( \tilde {u}_{\kern-1pt j}^*\tilde {p} \right ) \qquad \text{(using the continuity equation)}; \end{align}

(4.3) \begin{align} \frac{1}{\textit{Re}}\tilde {u}_i^*\frac{\partial ^2 \tilde {u}_i}{\partial x_{\kern-1pt j} \partial x_{\kern-1pt j}} & = \frac{1}{\textit{Re}}\frac{\partial }{\partial x_{\kern-1pt j}}\left ( \tilde {u}_i^*\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} \right ) - \frac{1}{\textit{Re}}\frac{\partial \tilde {u}_i^*}{\partial x_{\kern-1pt j}}\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}}; \end{align}

(4.4) \begin{align} \nu _t\tilde {u}_i^*\frac{\partial ^2 \tilde {u}_i}{\partial x_{\kern-1pt j} \partial x_{\kern-1pt j}} & = \frac{\partial }{\partial x_{\kern-1pt j}}\left ( \nu _t\tilde {u}_i^*\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} \right ) - \nu _t\frac{\partial \tilde {u}_i^*}{\partial x_{\kern-1pt j}}\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} - \tilde {u}_i^*\frac{{\textrm{d}}\nu _t}{{\textrm{d}}y}\frac{\partial \tilde {u}_i}{\partial y}. \end{align}

Substituting the above identities in (4.1), then integrating on the secondary-stability $y{-}z$ domain ( $\varOmega$ ) and normalising the eigenmode such that $\int _{\varOmega }\tilde {u}_i^*\tilde {u}_i \; {\textrm{d}}\varOmega = 1$ , a decomposition for the complex eigenvalue $\sigma$ is obtained

(4.5) \begin{equation} \sigma = T_a + \mathcal{P} - \mathcal{D} - \mathcal{D}^c + C, \end{equation}

where

(4.6) \begin{align} T_a &= \underbrace {-\iota k_x \int _{\varOmega } U \tilde {u}_i^*\tilde {u}_i \; {\textrm{d}}\varOmega }_{\displaystyle T_{a1}} \underbrace {- \int _{\varOmega } U_{\kern-1pt j}\tilde {u}_i^*\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} \; {\textrm{d}}\varOmega }_{\displaystyle T_{a2}}; \end{align}

(4.7) \begin{align} \mathcal{P} \; &= \underbrace {-\int _{\varOmega } \tilde {u}^*\tilde {v}\frac{\partial U}{\partial y} \; {\textrm{d}}\varOmega }_{\displaystyle \mathcal{P}_{uy}} \underbrace {-\int _{\varOmega } \tilde {u}^*\tilde {w}\frac{\partial U}{\partial z} \; {\textrm{d}}\varOmega }_{\displaystyle \mathcal{P}_{uz}} \underbrace {-\int _{\varOmega } \tilde {v}^*\tilde {v}\frac{\partial V}{\partial y} \; {\textrm{d}}\varOmega }_{\displaystyle \mathcal{P}_{vy}} \underbrace {-\int _{\varOmega } \tilde {v}^*\tilde {w}\frac{\partial V}{\partial z} \; {\textrm{d}}\varOmega }_{\displaystyle \mathcal{P}_{vz}}; \end{align}

(4.8) \begin{align} & \quad \underbrace {-\int _{\varOmega } \tilde {w}^*\tilde {v}\frac{\partial W}{\partial y} \; {\textrm{d}}\varOmega }_{\displaystyle \mathcal{P}_{wy}} \underbrace {-\int _{\varOmega } \tilde {w}^*\tilde {w}\frac{\partial W}{\partial z} \; {\textrm{d}}\varOmega }_{\displaystyle \mathcal{P}_{wz}}; \end{align}

(4.9) \begin{align} \mathcal{D} \; &= \; \frac{k_x^2}{\textit{Re}}\int _{\varOmega }\tilde {u}_i^*\tilde {u}_i \; {\textrm{d}}\varOmega + \frac{1}{\textit{Re}}\int _{\varOmega } \frac{\partial \tilde {u}_i^*}{\partial x_{\kern-1pt j}}\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} {\textrm{d}}; {\textrm{d}}\varOmega ; \end{align}

(4.10) \begin{align} \mathcal{D}^c & = \; k_x^2\int _{\varOmega }\nu _t\tilde {u}_i^*\tilde {u}_i \; {\textrm{d}}\varOmega + \int _{\varOmega } \nu _t\frac{\partial \tilde {u}_i^*}{\partial x_{\kern-1pt j}}\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} \; {\textrm{d}}\varOmega ; \end{align}

(4.11) \begin{align} C \hspace {0.15cm} & = \; \underbrace {\iota k_x \int _{\varOmega } \tilde {u}^*\tilde {v} \frac{{\textrm{d}}\nu _t}{{\textrm{d}}y} \; {\textrm{d}}\varOmega }_{\displaystyle C_1} + \underbrace {\int _{\varOmega } \frac{{\textrm{d}}\nu _t}{{\textrm{d}}y}\tilde {u}_{\kern-1pt j}^*\frac{\partial \tilde {v}}{\partial x_{\kern-1pt j}} \; {\textrm{d}}\varOmega }_{\displaystyle C_2}.\end{align}

Equation (4.5) represents a Reynolds–Orr-type energy budget for an eigenvector perturbation to the two-dimensional base flow, extended to the turbulent regime modelled by a turbulent eddy viscosity. We should note that the following terms integrate to zero because of the boundary conditions in $y$ and $z$ :

(4.12) \begin{align} T_p = & \; -\int _{\varOmega }\frac{\partial }{\partial x_{\kern-1pt j}}\big ( \tilde {u}_{\kern-1pt j}^*\tilde {p} \big )\,{\textrm{d}}\varOmega \equiv 0; \end{align}

(4.13) \begin{align} T_d = & \; \frac{1}{\textit{Re}} \int _{\varOmega }\frac{\partial }{\partial x_{\kern-1pt j}}\left ( \tilde {u}_i^*\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} \right )\,{\textrm{d}}\varOmega \equiv 0; \end{align}

(4.14) \begin{align} T_d^c = & \; \int _{\varOmega } \frac{\partial }{\partial x_{\kern-1pt j}}\left ( \nu _t\tilde {u}_i^*\frac{\partial \tilde {u}_i}{\partial x_{\kern-1pt j}} \right ) \; {\textrm{d}}\varOmega \equiv 0. \end{align}

The relation (4.5) shows that the complex-valued eigenvalue $\sigma$ associated with one given eigenmode results from several contributions: advective transport ( $T_a$ ), production due to base flow gradients ( $\mathcal{P}$ ), dissipation due to molecular viscosity ( $\mathcal{D}$ ), dissipation due to eddy viscosity ( $\mathcal{D}^c$ ) and production due to eddy viscosity gradients ( $C$ ). Of course the last two terms are present only when the eddy viscosity model is used. Some remarks need to be made about these quantities: (i) $T_a$ is purely imaginary, hence it contributes to the angular frequency $\textrm{Im}(\sigma )$ but does not contribute to the growth rate $ \textit{Re}(\sigma )$ (as expected for a transport term); (ii) $\mathcal{D}$ , $\mathcal{D}^c$ and $C$ are purely real with $\mathcal{D} \geqslant 0$ and $\mathcal{D}^c \geqslant 0$ and thus contribute to the damping of the mode, whereas $C$ generally contributes positively to the growth rate; (iii) $\mathcal{P}$ is complex but contributes mainly to the real part of the eigenvalue with a positive amount. The main contribution to the imaginary part of $\sigma$ (hence to the phase velocity of the mode) comes from $T_{a1}$ , i.e. from the mean flow advection.

Figure 11. Energy budget for the leading eigenmode at $ \textit{Re}_{\tau }=71$ and $k_x=0.18$ for (a) the quasi-laminar model ( $A_s=0.2$ ) and (b) the eddy viscosity model ( $A_s=0.5$ ). Blue bars denote positive contributions to the growth rate while orange bars denote negative contributions. For the meaning of the labels see ((4.5)–(4.11)) in the text.

The terms contributing to the real part of the eigenvalues (the modal growth rate) are analysed in more detail on specific examples. We focus on the parameter values $ \textit{Re}_{\tau }=71$ , $k_x=0.18$ , $A_s=0.2$ for the quasi-laminar model and $A_s=0.5$ for the eddy viscosity model. Figure 11 shows the various terms from (4.5) for the two cases. For the quasi-laminar case, the important terms are $\mathcal{P}_{uy}$ , $\mathcal{P}_{uz}$ and $\mathcal{D}$ . The remaining production terms are several orders of magnitude smaller. For the eddy viscosity case, the eddy viscosity dissipation $\mathcal{D}^c$ is another important term. On the contrary, the production due to eddy viscosity gradients $C$ appears as negligible. Nonetheless, this does not imply that the wall-normal variation of the eddy viscosity is not important. All the quantities plotted are integrals that involve the eigenmodes and, in general, it can be expected that the wall-normal variation of $\nu _t$ contributes to the structure of the eigenmodes. It was verified by inspection that similar considerations apply to all the cases considered in the following analysis.

Figure 12. Energy budget for the leading eigenmode with $k_x=0.18$ . (a–b) Variation with streak amplitude $A_s$ for fixed $ \textit{Re}_{\tau }=71$ . (c–d) Variation with $ \textit{Re}$ for fixed streak amplitude $A_s=0.2$ in (c) and $A_s=0.5$ in (d). (a–c) Quasi-laminar model; (b–d) eddy viscosity model. For the meaning of the labels see ((4.5)–(4.11)) in the text.

4.2. Results

Having identified the dominant terms, their dependence on the stability parameters is analysed in figure 12. Concerning the effect of $A_s$ (panels (a) and (b)), for both models the production term inducing instability is the one associated with spanwise gradients of the base flow, $\mathcal{P}_{uz}$ . For the eddy viscosity model, a larger amplitude is needed to compensate for the dissipation due to the eddy viscosity.

The panels (c) and (d) in figure 12 illustrate why a critical $ \textit{Re}$ is observed only within the eddy viscosity model. For the quasi-laminar model (panel (c)), production and dissipation follow essentially the same trend, resulting in a growth rate almost independent of $ \textit{Re}_{\tau }$ . Conversely, for the eddy viscosity model (panel (d)), the eddy viscosity dissipation increases faster than the production when the $ \textit{Re}_{\tau }$ is increased. Consequently, at some critical value of $ \textit{Re}_{\tau }$ the eddy viscosity dissipation takes over the production and stabilises the eigenmode. This suggests a $ \textit{Re}$ dependence for the growth rates encoded in the $ \textit{Re}$ dependence of the eddy viscosity, an interesting result pointing at the role of the turbulent fluctuations in this instability.

The two panels of figure 13 illustrate how the energy budget changes with the streamwise wavenumber $k_x$ . For both models, $\mathcal{P}_{uy}$ and $\mathcal{P}_{uz}$ depart from each other, the former decreasing while the latter increases. For the quasi-laminar model in panel (a), this behaviour changes drastically beyond $k_x=1.5$ , indicating a change in the nature of the leading eigenmode. However, the main difference between the two models is again the presence of the eddy viscosity dissipation that stabilises the modes for $k_x \gt 0.4$ .

In conclusion, the energy budget analysis shows that the inclusion of the eddy viscosity induces the observed parameter dependence mainly via the additional damping term $\mathcal{D}^c$ .

Figure 13. Energy budget for the leading eigenmode as $k_x$ is varied. $ \textit{Re}_{\tau }=71$ . (a) Quasi-laminar model with $A_s=0.2$ ; (b) eddy viscosity model with $A_s=0.5$ . Note the different scale of the abscissa in the two panels. For the meaning of the labels see ((4.5)–(4.11)) in the text.