# Large-scale streak instabilities of transitional channel flow

Nicola Ciola, Yohann Duguet, Jean-Christophe Robinet, Pietro de Palma, Stefania Cherubini

► To cite this version:

Nicola Ciola, Yohann Duguet, Jean-Christophe Robinet, Pietro de Palma, Stefania Cherubini. Large-scale streak instabilities of transitional channel flow. Journal of Fluid Mechanics, 2025, 1025, pp.A45. 10.1017/jfm.2025.10974. hal-05435618

HAL Id: hal-05435618

https://hal.science/hal-05435618v1

Submitted on 30 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Large-scale streak instabilities of transitional channel flow

Nicola Ciola<fnref n="1" />,<fnref n="3" /> ID, Yohann Duguet<fnref n="2" /> ID, Jean-Christophe Robinet<fnref n="1" /> ID, Pietro De Palma<fnref n="3" /> and Stefania Cherubini<fnref n="3" /> ID

**Corresponding author:** Nicola Ciola, n.ciola@phd.poliba.it

(Received 19 June 2025; revised 15 October 2025; accepted 17 November 2025)

The emergence of large-scale spatial modulations of turbulent channel flow, as the Reynolds number is decreased, is addressed numerically using the framework of linear stability analysis. Such modulations are known as the precursors of laminar–turbulent patterns found near the onset of relaminarisation. A synthetic two-dimensional base flow is constructed by adding finite-amplitude streaks to the turbulent mean flow. The streak mode is chosen as the leading resolvent mode from linear response theory. In addition, turbulent fluctuations can be taken into account or not by using a simple Cess eddy viscosity model. The linear stability of the base flow is considered by searching for unstable eigenmodes with wavelengths larger than the base flow streaks. As the streak amplitude is increased in the presence of the turbulent closure, the base flow loses its stability to a large-scale modulation below a critical Reynolds-number value. The structure of the corresponding eigenmode, its critical Reynolds number, its critical angle and its wavelengths are all fully consistent with the onset of turbulent modulations from the literature. The existence of a threshold value of the Reynolds number is directly related to the presence of an eddy viscosity, and is justified using an energy budget. The values of the critical streak amplitudes are discussed in relation with those relevant to turbulent flows.

**Key words:** shear-flow instability, transition to turbulence, pattern formation

## 1. Introduction

Any description of a turbulent flow starts with the identification of the most relevant coherent structures and of their dynamics. For shear flows in the vicinity of solid walls,

<footnote n="1">¹DynFluid, Arts et Métiers Paris /CNAM, 151 Bd de l'Hôpital, 75013 Paris, France</footnote>

<footnote n="2">²LISN-CNRS, Université Paris-Saclay, 507 Rue du Belvédère, 91405 Orsay, France</footnote>

<footnote n="3">³DMMM, Politecnico di Bari, Via Re David, 200, 70125 Bari, Italy</footnote>
streamwise streaks, i.e. spanwise modulations of the streamwise velocity, together with the associated streamwise vortices, are by far the most reported coherent structures (Kline et al. 1967). Streamwise streaks do not emerge from the laminar flow as an instability but rather as a manifestation of the non-normal amplification of perturbations induced by streamwise vortices that could already be present in the associated laminar shear flow (Schmid, Henningson & Jankowski 2002). The distribution of the spanwise widths of the streaks peaks at a wavelength close to 100 wall units (Kline et al. 1967). This scale is associated with the immediate vicinity of the wall, however, in this study, we are mostly concerned with marginally low-Reynolds-number values for which streaks and vortices extend across the channel gap, making the distinction between the near-wall and outer regions insignificant. Streamwise streaks are themselves linearly unstable to various types of velocity perturbations with different possible length scales (Schoppa & Hussain 2002) and spatial symmetries (Andersson et al. 2001); see Lozano-Durán et al. (2021) for a review of such instabilities in the literature. At the scale of the streaks themselves, streak instabilities generate additional vorticity that feeds back on the streamwise rolls and prevents them from viscously decaying (Hamilton, Kim & Waleffe 1995; Waleffe 1997). This nonlinear mechanism, known as the self-sustaining process (SSP), is local in nature and explains the existence of non-laminar flow, including turbulence, finite-amplitude unsteady perturbations and edge states (Waleffe 1997). Streak instabilities at larger scales have also been reported as one possible explanation for the occurrence of large-scale motions (LSMs) frequently reported in high-Reynolds-number experiments and simulations (Jiménez 2013; de Giovanetti et al. 2017). Such an instability scenario coexists with a transient growth mechanism, justifying LSMs as the flow disturbances most amplified by non-normal effects (Del Alamo & Jimenez 2006; Pujals et al. 2009). Interestingly, high-Reynolds-number flows have always been a topic of investigation in the shear flow literature, as justified by numerous real world applications. However, at the lowest Reynolds numbers where turbulence exists, the structure of the turbulent flow features a different class of large-scale coherent structures, which is so far poorly understood. Focusing on the example of channel flow between two infinitely large parallel plates and driven by a pressure gradient, as the friction Reynolds number ($Re_τ$) is decreased to just below $95 \pm 1$, the turbulent flow begins to display modulations with an oblique wavevector (Shimizu & Manneville 2019; Kashyap, Duguet & Dauchot 2020), interpreted recently as a linear instability of the turbulent regime (Kashyap, Duguet & Dauchot 2022). For lower $Re_τ$, the amplitude of the modulations increases and the weaker zones relaminarise, leading to laminar–turbulent patterns (Tsukahara et al. 2005; Tuckerman et al. 2014, 2020). For even lower $Re_τ$, the patterns break into solitary bands of turbulence (Shimizu & Manneville 2019). A similar phenomenology was reported in most wall-bounded flows, yet sometimes without a clear evidence for the modulational stage as e.g. in plane Couette flow (Gomé et al. 2023). With the hope that the simplicity of the spatial modulations is likely to hide a simpler phenomenology, we choose to investigate the channel case; hence, we focus in this article on the modulational regime of turbulent channel flow and aim at predicting quantitatively such modulations using direct numerical simulation and tools from linear stability analysis (LSA).

Applying hydrodynamic stability concepts to a turbulent base flow poses serious conceptual challenges. In principle, in the deterministic context of the Navier–Stokes equations, linearisation is carried out in the neighbourhood of some relevant equilibrium flow state characterised by a temporal symmetry, e.g. a steady state, a travelling wave or a relative periodic orbit. Such flow states should be exact solutions of the governing equations. In shear flows, such self-sustained states always feature three-dimensional streamwise streaks that are linearly unstable (Waleffe 2001; Kawahara et al. 2012).
A few recent attempts to link the occurrence of laminar–turbulent patterning to bifurcations of such three-dimensional solutions of the Navier–Stokes equations were reported in the context of pipe flow (Chantry, Willis & Kerswell 2014) and plane Couette flow (Reetz, Kreilos & Schneider 2019). In the channel flow geometry, several solutions akin to patterning were also reported recently (Paranjape et al. 2020, 2023) but they have so far not been linked to spatially unmodulated solutions. Given that the exact dynamical relevance of such solutions to the turbulent dynamics remains subtle (Cvitanović 2013), we find it here preferable to focus on a different approach that actively takes into account the presence of turbulent fluctuations. Linearising around a chaotic attractor using concepts such as Lyapunov exponents and covariant Lyapunov vectors (Inubushi, Takehiro & Yamada 2015) is such a possibility, albeit a computationally heavy one. Taking into account the very fine details of the base flow, both spatially and temporally, results anyway in an increasingly complex analysis from which no simple physical interpretation might emerge. A recent attempt to extract the dominant large-scale instability using this approach in plane Couette flow did not lead to a prediction of modulations characterised by an oblique wavevector (Ishikawa, Takehiro & Yamada 2018). Instead, we adopt here a different theoretical approach, based on the choice of a simpler base flow, where statistical symmetries are as much as possible taken into account.

In the last decade, hydrodynamic stability calculations have begun to become relevant for the generation of coherent structures within a turbulent context, once it was realised that LSA could be carried out, in practice, in the vicinity of the turbulent mean flow rather than the laminar base flow. Because of the invariance of the problem with respect to translations in time and in the two planar directions, the mean flow, understood in practice as the time-averaged velocity field, depends only on the wall-normal coordinate $y$, but not on the streamwise coordinate $x$, the spanwise coordinate $z$ or the time $t$. Intriguingly, this approach fails systematically in channel flow: the mean flow of turbulent channel flow was reported to be linearly stable to all perturbations, both for high (Reynolds & Tiederman 1967) and low values of $Re_τ$ (Kashyap, Duguet & Dauchot 2024). A glance at any velocity field in the turbulent regime (see e.g. figure 1 in Kashyap et al. (2022)) reveals the ubiquity of long-lasting streamwise streaks. The relatively rapid emergence time of the large-scale modulations, comparable to the viscous time scale, suggests that the $y$-dependent mean flow is poorly relevant as a base flow for LSA in the present context. Indeed, from most visualisations it is clear that the large-scale modulations grow on a streak field rather than on a bare, structure-free base flow. We hence revisit the linearisation of Kashyap et al. (2025) by adding well-controlled streamwise streaks to the laminar base flow, and by allowing for linear instabilities with length scales larger than the characteristic size of streaks. One of the motivations for including streamwise streaks is their coherence, i.e. their presence in the mean flow if the time average holds over a finite time only. Finite-time averaging is, however, not rigorous proper Reynolds averaging, and linear LSA around such a base flow is ambiguous to interpret theoretically. Because of the turbulent context, fluctuations and their nonlinear self-interaction can also not be ignored. Reynolds stresses must be accounted for and modelled. The quest for an optimal parametrisation of turbulent closures in channel flow is outside the scope of this study. We focus here on a simple and popular eddy viscosity closure, already used in former mean flow stability computations (Del Alamo & Jimenez 2006; Pujals et al. 2009; Morra et al. 2019), which we can decide to either include or not in the analysis. Although this eddy viscosity model is simplistic, as we shall see, it already leads to meaningful conclusions regarding the possibility of large-scale instabilities as precursors of turbulent modulations and laminar–turbulent patterning.

This paper is structured as follows. Section 2 contains the flow governing equations and the primary resolvent analysis from which synthetic streaks are computed. The results of
the LSA of the streaky base flow are shown in § 3. Section 4 contains an energy budget analysis applied to the unstable mode found. The model assumptions are discussed in light of the results in § 5. Conclusions and outlooks are eventually given in § 6.

## 2. Flow set-up and resolvent modes

### 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 \le y \le 2$) and the (imposed) channel bulk velocity $U_b$ ($U_b = \int_0^2 \int_0^{L_z} u \, dydz/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 = v/u_\tau$ and the friction velocity $u_\tau = \sqrt{\tau_w/\rho}$, where $\rho$ is the fluid density, $v$ 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

$$ \left\{ \begin{aligned} \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} &= -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} + \mathbf{f}_b, \\ \nabla \cdot \mathbf{u} &= 0, \end{aligned} \right. \qquad (2.1) $$

where $p$ is the hydrodynamic pressure, $Re = U_b h/v$ is the Reynolds number and $f_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 channel flow 2.0 code (Gibson et al. 2021). 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 2019; Kashyap et al. 2020; Kashyap et al. 2022). 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 $Re_\tau = u_\tau h/v$, which can be computed a posteriori, to ease comparison with previous works.

At the lowest Reynolds number considered here ($Re_\tau = 71$), the turbulent flow appears as a pattern of oblique bands (figure 1c). The alternance of extended blue and red regions indicates the presence of the large-scale flow which coexists with the pattern (Duguet & Schlatter 2013). At a higher value of $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 1b). Instead, for $Re_\tau$ above 95, turbulence appears in DNS to be unmodulated by large-scale oblique structures (figure 1a), consistently with Kashyap et al. (2022), and will be qualified as featureless after Shimizu & Manneville (2019).
<figure><img src="image_3.png" /><figcaption>Figure 1. Contours of streamwise velocity fluctuations in the plane $y^+ \approx 35$ from DNSs at (a) $Re_\tau = 98$, (b) $Re_\tau = 92$ and (c) $Re_\tau = 71$ ($y \approx 0.36$ for $Re_\tau = 98$, 92 and $y \approx 0.49$ for $Re_\tau = 71$). The black vertical segment at the bottom left corner of each panel indicates a spanwise length of 1000 wall units.</figcaption></figure>
## 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. (2022), 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 2010b; McKeon & Sharma 2010; Moarref et al. 2013; Sharma & McKeon 2013) or the initial value problem (Del Alamo & Jimenez 2006; Pujals et al. 2009). The results of this analysis are briefly revised here at low Re.

Using the Reynolds decomposition, the solution of (2.1) can be partitioned as the sum of the mean flow $\overline{U} = \overline{U}(y)e_x$ plus fluctuations. Then, following Reynolds & Hussain (1972), 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:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{u}'}{\partial t} + \overline{\mathbf{U}} \cdot \nabla \mathbf{u}' + \mathbf{u}' \cdot \nabla \overline{\mathbf{U}} + \nabla p' - \frac{1}{Re} \nabla^2 \mathbf{u}' - \nabla \cdot [\nu_t (\nabla \mathbf{u}' + (\nabla \mathbf{u}')^T)] = \mathbf{f}, \quad (2.2) \\[1.5em] \displaystyle \nabla \cdot \mathbf{u}' = 0, \end{array} \right. $$

for the coherent velocity fluctuation $\mathbf{u}'$ (respectively $p'$ for the coherent pressure fluctuation), where $\mathbf{f}$ is a generic forcing term (cf. Morra et al. (2019)). 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 2006; Pujals et al. 2009; Hwang & Cossu 2010b; Park, Hwang & Cossu 2011; Alizard 2015; Ciola et al. 2024), the Cess (1958) eddy viscosity formula, as reported by Reynolds & Tiederman (1967), is employed

$$ \nu_t(y) = \frac{1}{2Re} \left\{ \left[ 1 + \left( \frac{Re_\tau k}{3} (2y - y^2)(3 - 4y + 2y^2) \left(1 - e^{(|y-1|-1)Re_\tau/A}\right) \right)^2 \right]^{1/2} - 1 \right\}, \quad (2.3) $$

with $k = 0.426$, $A = 25.4$ (Del Alamo & Jimenez 2006) 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 2018; Morra et al. 2019; Symon et al. 2023). In this study, we will compare stability results with and without eddy viscosity. Borrowing the terminology from Morra et al. (2019), (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 (1967).

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 $\hat{\mathbf{u}}(y; k_x, k_z, \omega) = \mathcal{H}(y; k_x, k_z, \omega) \hat{\mathbf{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 $\hat{f}(y)$ that maximises the amplitude ratio between response and forcing (cf. Hwang & Cossu 2010b)

$$R(k_x, k_z, \omega) = \max_{\hat{f} \neq 0} \frac{\|\hat{u}\|^2}{\|\hat{f}\|^2} = \|\mathcal{H}(y; k_x, k_z, \omega)\|^2, \quad (2.4)$$

where the standard $L^2$ norm in $y$ is used both for $\hat{u}$ and $\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 2014; Symon et al. 2018). After Fourier transform in $x$, $z$ and $t$ and discretisation in $y$, the system (2.2) can be written as $(i\omega\mathbf{B} - \mathbf{L})\hat{q} = \hat{f}$, where $\hat{q}$ contains the four primitive variables $\hat{u}$, $\hat{v}$, $\hat{w}$ and $\hat{p}$ for each point of the computational grid, $\mathbf{B}$ is the prolongation matrix (Cerqueira & Sipp 2014), $\mathbf{L}$ is the discretised linear Navier-Stokes operator and $i$ the imaginary unit. The discretised resolvent operator can be computed by a matrix inversion as $\mathbf{H} = (i\omega\mathbf{B} - \mathbf{L})^{-1}$ and decomposed such that $\mathbf{H} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^T$. The most amplified harmonic forcing is given by the first column of $\mathbf{V}$, the corresponding response mode by the first column of $\mathbf{U}$ and the energy amplification factor $R$ by the square of the leading singular value of $\mathbf{H}$, which is the first diagonal element of $\mathbf{\Sigma}$. For the discretisation of differential operators, $N_y = 65$ Chebyshev collocation points have proven sufficient at the low values of $Re$ considered. The code was validated by reproducing the results of Hwang & Cossu (2010b).

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 1997) or data-driven techniques where streaks are directly extracted from DNSs (Hack & Zaki 2014) or experimental data (Liu et al. 2024). 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. 2024). Streamwise streaks are routinely explained by the advection of streamwise vortices by the lift-up effect (Ellingsen & Palm 1975; Landahl 1980), 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 $f$ represents the vortices and the response $u'$ 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 (2010b). The problem is solved for a range of $Re$ relevant to the modulational regime of turbulent patterns in the channel, in particular $Re_\tau = [71, 78, 84, 91, 95, 98, 105]$.

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
<figure><img src="image_4.png" /><figcaption>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) $Re_{\tau} = 71$ and (c-d-g-h) $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.</figcaption></figure>

shaded contours and the transverse components as arrows in the $y-z$ plane for two values of $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. Stability analysis

### 3.1. Formulation

The lift-up effect is only one of the several mechanisms involved in the SSP of wall turbulence (Hamilton *et al.* 1995; Waleffe 1997). 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 2002; Park *et al.* 2011; Alizard 2015). In the classical SSP picture, that stage precedes the feedback onto the streamwise vortices. More recently, de Giovanetti *et al.* (2017) and Ciola *et al.* (2024) 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 $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 $\mathbf{u}_s$. Due to the linearity of (2.2), $\mathbf{u}_s$ is defined up to a multiplicative constant. Following our previous work (Ciola et al. 2024), $\mathbf{u}_s$ is rescaled such that

$$ \frac{\max_{y,z} u_s - \min_{y,z} u_s}{2\bar{U}_c} = 1, \qquad (3.1) $$

where $u_s$ (as a scalar field) is the streamwise velocity component of the streaks field and $\bar{U}_c$ is the centreline velocity of the turbulent mean profile.

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

$$ U = \bar{U} + A_s \mathbf{u}_s, \qquad (3.2) $$

with $A_s > 0$ interpreted as the streak amplitude. Consistently with the previous section, $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 $U$ verifies an equation of the generic form

$$ \left\{ \begin{array}{l} U \cdot \nabla U + \nabla P - \frac{1}{Re} \nabla^2 U - \nabla \cdot [v_t(\nabla U + (\nabla U)^T)] - f_b = \mathbf{g}, \\ \nabla \cdot U = 0, \end{array} \right. \qquad (3.3) $$

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

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

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p - \frac{1}{Re} \nabla^2 \mathbf{u} - \nabla \cdot [v_t(\nabla \mathbf{u} + (\nabla \mathbf{u})^T)] - f_b = \mathbf{g}, \\ \nabla \cdot \mathbf{u} = 0. \end{array} \right. \qquad (3.4) $$

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

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

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{u}''}{\partial t} = -U \cdot \nabla \mathbf{u}'' - \mathbf{u}'' \cdot \nabla U - \nabla p'' + \frac{1}{Re} \nabla^2 \mathbf{u}'' + \nabla \cdot [v_t(\nabla \mathbf{u}'' + (\nabla \mathbf{u}'')^T)], \\ \nabla \cdot \mathbf{u}'' = 0. \end{array} \right. \qquad (3.5) $$

The forcing term $\mathbf{g}$ cancels out in the secondary perturbation equation because it is constant, i.e. it depends only on $U$ and not on $\mathbf{u}''$. Similarly, $f_b$ cancels out because we assume that $\mathbf{u}''$ 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 $U$ and the same viscosity $v_t$, are expected to yield the same stability results for $u''$.

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

$$ \mathbf{u}'' = \tilde{\mathbf{u}}(y, z)e^{\sigma t + i k_x x} + \text{complex conjugate}, \quad \sigma \in \mathbb{C} \text{ and } k_x \in \mathbb{R}, \quad (3.6) $$

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

$$ \left\{ \begin{array}{l} \sigma \tilde{u} = -i k_x U \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} - i k_x \tilde{p} + \left( \frac{1}{Re} + v_t \right) \nabla^2 \tilde{u} + \frac{dv_t}{dy} \left( \frac{\partial \tilde{u}}{\partial y} + i k_x \tilde{v} \right), \\ \\ \sigma \tilde{v} = -i k_x U \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}{Re} + v_t \right) \nabla^2 \tilde{v} + 2 \frac{dv_t}{dy} \frac{\partial \tilde{v}}{\partial y}, \\ \\ \sigma \tilde{w} = -i k_x U \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}{Re} + v_t \right) \nabla^2 \tilde{w} + \frac{dv_t}{dy} \left( \frac{\partial \tilde{w}}{\partial y} + \frac{\partial \tilde{v}}{\partial z} \right), \\ \\ 0 = i k_x \tilde{u} + \frac{\partial \tilde{v}}{\partial y} + \frac{\partial \tilde{w}}{\partial z}, \end{array} \right. \quad (3.7) $$

$$ \text{with } \nabla^2 = -k_x^2 + \partial_y^2 + \partial_z^2. $$

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

$$ \sigma \mathbf{B}_0 \tilde{\mathbf{q}}_0 = \mathbf{A}_0 \tilde{\mathbf{q}}_0, \quad (3.8) $$

with $\tilde{\mathbf{q}}_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 $\mathbf{A}_0$ and $\mathbf{B}_0$ are matrices of size $M \times M$. The matrix $\mathbf{B}_0$ is the so-called prolongation matrix that arises when using a velocity-pressure formulation (Cerqueira & Sipp 2014).

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. (2017). 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

$$ \sigma \mathbf{B} \tilde{\mathbf{q}} = \mathbf{A} \tilde{\mathbf{q}}, \quad (3.9) $$

where $\tilde{\mathbf{q}}$ has $N_u$ times the size of the original $\tilde{\mathbf{q}}_0$ vector in (3.8), i.e. $MN_u$, and the matrices $\mathbf{A}$ and $\mathbf{B}$ have now size $MN_u \times 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 $\mathbf{A}$ has a block-circulant structure. One can hence diagonalise it into a block-diagonal matrix $\hat{\mathbf{A}} = \mathbf{P}^H \mathbf{A} \mathbf{P}$
with $N_u$ blocks $\hat{\mathbf{A}}_0, \hat{\mathbf{A}}_1, \dots, \hat{\mathbf{A}}_{N_u}$, thereby reducing one large eigenvalue problem of size $N_u M \times N_u M$ down to $N_u$ manageable eigenproblems, each of size $M \times M$. The matrix $\mathbf{P} \in \mathbb{C}^{N_u M \times N_u M}$ is given by the Kronecker product $\mathbf{J} \otimes \mathbf{I}_M$ with $\mathbf{J} \in \mathbb{C}^{N_u \times N_u}$, $J_{j+1,k+1} = \rho_j^k / \sqrt{N_u}$ and $\mathbf{I}_M$ the identity matrix of size $M \times M$ (Schmid et al. 2017). Each of these spectral subproblems involves a matrix of size $M$, namely $\hat{\mathbf{A}}_j = \mathbf{A}_0 + \rho_j \mathbf{A}_1 + \rho_j^2 \mathbf{A}_2 + \dots + \rho_j^{N_u-1} \mathbf{A}_{N_u-1}$ for each $j = 0, \dots, N_u - 1$, featuring $\mathbf{A}_0, \mathbf{A}_1, \dots$ which are the blocks of $\mathbf{A}$. The general eigensolutions of (3.9) can be written (Schmid et al. 2017) as

$$ \tilde{\mathbf{q}} = (\tilde{\mathbf{q}}_j, \rho_j \tilde{\mathbf{q}}_j, \dots, \rho_j^{N_u-1} \tilde{\mathbf{q}}_j)^T, \quad (3.10) $$

based on the known eigenvectors $\tilde{\mathbf{q}}_j$ of the matrices $\hat{\mathbf{A}}_j$.

In (3.10), the complex numbers $\rho_j$, $j = 1, \dots, 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{\mathbf{q}}$ 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 t)$. The $\gamma$ values are defined in [0, 1). They take only discrete values $0, 1/N_u, 2/N_u, \dots$. 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. 2024). 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 $\hat{\mathbf{A}}_j$ yields the spectrum of the original matrix $\mathbf{A}$. In practice, the latter is found by looping over $j = 0, \dots, 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 1988), 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. 2024).

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. 2024).

### 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. (2024), based on
<figure><img src="image_5.png" /><figcaption>Figure 3. Eigenvalues for the streak stability problem (only a subset of the computed spectra is shown) for $Re_t = 71$ and $k_x = 0.18$. The eigenvalues (● for stable modes, ★ for unstable modes) are coloured with the corresponding root of unity factor $\gamma = j/N_u$ for $j = 0, \dots, 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.</figcaption></figure>

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 2002). 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. 2001), 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. (2022).

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/ik_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
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 4. Variation of the leading eigenvalues with $Re$ and streak amplitude $A_s$ for $k_x = 0.18$ (● for stable modes, ★ for unstable modes). (a–b) Growth rate; (c–d) phase velocity. (a–c) Quasi-laminar model; (b–d) eddy viscosity model.</figcaption></figure>

from those in (0, 0.5) by spanwise reflection (see § 3.2). Focusing on $Re_τ = 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 ≈0.16, whereas for the eddy viscosity model it is larger, ≈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 $γ > 0$.

The $Re_τ$ 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 $Re$. A critical amplitude can still be defined, independently of $Re_τ$, 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 $Re_τ$. As a result, for $A_s$ large enough, as $Re_τ$ 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 $Re$ with the eddy viscosity model is a qualitative improvement with respect to the mean flow analysis in Kashyap et al. (2024).

The bottom panels of figure 4 show the phase velocity of the leading modes with respect to $Re_τ$ and $A_s$. The values globally overestimate by roughly 5%–10% the advection velocities of turbulent bands measured from DNS (Tuckerman et al. 2014), however, they display a similar decrease with $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><img src="image_10.png" /><figcaption>Figure 5. Eigenvalues coloured by the large-scale spanwise energy ratio $r_{LS}$ (only a subset of the computed spectra is shown) for $Re_τ = 71$ and $k_x = 0.18$ (• for stable modes, ★ for unstable modes). (a) Quasi-laminar model and $A_s = 0.2$; (b) eddy viscosity model and $A_s = 0.5$.</figcaption></figure>

### 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 $λ_{z,s}$, yet the question arises whether some eigenmodes have a particularly prominent large-scale component. In this respect, the detuning factor $γ$ yields only partial information. Therefore, an additional quantification of the large-scale property is introduced by

$$r_{LS} = \frac{\int_{0}^{2} \sum_{|k_z| < k_z^c} |\check{u}(y, k_z)|^2 dy}{\int_{0}^{2} \sum_{k_z} |\check{u}(y, k_z)|^2 dy}, \quad (3.11)$$

where $\check{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 $Re$. The new quantity $r_{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.

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_{LS}$ was designed to give the ratio of the
<figure><img src="image_11.png" /><figcaption>Figure 6. Leading eigenmode for the eddy viscosity model with $Re_τ = 71$, $k_x = 0.18$ and $A_s = 0.50$ ($λ_{z,s} \approx 1.41$ at this value of $Re_τ$). The modes are normalised to $max_{y,z}|Re(\tilde{u})| = 1$. Real part of (a) streamwise velocity component, (b) wall-normal velocity component and (c) spanwise velocity component.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>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π/λ_{z,s}$.</figcaption></figure>

energy of this large-scale component with respect to the total energy of the eigenmode. The other waves have wavelengths near $λ_{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°$ with respect to the streamwise direction. This angle is not far from the 23° reported by Kashyap et al. (2022) and the 22.5° found by Benavides & Barkley (2025).

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
N. Ciola, Y. Duguet, J.-C. Robinet, P. De Palma and S. Cherubini

<figure><img src="image_13.png" /><figcaption>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.</figcaption></figure>

velocity components forming the large-scale flow. The displayed velocity field is fully consistent with most observations in the patterning regime (Duguet & Schlatter 2013; Tuckerman et al. 2014). 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. 1995; Waleffe 1997). 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. (2022), namely $k_x = 0.18$, the eddy viscosity model gives a group of eigenmodes characterised by a large-scale modulation that becomes unstable as $Re_τ$ 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 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. (2024) 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 $Re_τ$ 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 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.

## 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. (2022), a critical *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 *Re* for a given *A*<sub>s</sub> (note that the critical *Re* is defined only for sufficiently high *A*<sub>s</sub>, figure 4b) or define a critical *A*<sub>s</sub> for each *Re*, so that a neutral curve *Re*−*A*<sub>s</sub> can be computed. For each point
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Figure 9. Variation of the leading growth rate $\sigma_r$ with $Re$ and streamwise wavenumber $k_x$ (● for stable modes, ★ 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.</figcaption></figure>

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 $Re_τ$ 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 $Re$. This means that, for fixed $A_s$, the instability is obtained by decreasing $Re$, as explained in § 3.3 and consistently with Kashyap et al. (2022). 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 $Re$ in the considered range.

Panels (b) and (c) of figure 10 show the critical wavenumbers. They increase with $Re$, indicating that the instability tends towards smaller wavelengths at high $Re$. This is also
<figure><img src="image_20.png" /><figcaption>Figure 10. Critical $A_s$, $k_x$ and large-scale $k_z$ as a function of $Re_τ$ for the eddy viscosity model.</figcaption></figure>

consistent with DNS observations (Kashyap et al. 2020). At $Re_τ = 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. (2022). 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 $Re$ from $25.5^\circ$ at $Re_τ = 71$ to $20.0^\circ$ at $Re_τ = 106$ (not shown).

## 4. Energy budget analysis

In order to interpret physically the results of the stability analysis of the streaky base flow, an equation equivalent to the Reynolds–Orr equation for the most unstable eigenvector is derived in this section (Albensoeder, Kuhlmann & Rath 2001; Schoppa & Hussain 2002; Schmid et al. 2002; Alizard 2015). For a simpler derivation an index notation is introduced, where $x_1, x_2$ and $x_3$ correspond, respectively, to $x$, $y$ and $z$ and $u_1, u_2$ and $u_3$ to $u$, $v$ and $w$. Summation over repeated indices is implicit, and $i$ takes the values 1, 2 and 3, whereas $j$ takes only the values 2 and 3. Lastly, in this section an asterisk denotes the complex conjugate.

### 4.1. Derivation

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

$$ \sigma \tilde{u}_i^* \tilde{u}_i = - \iota k_x U \tilde{u}_i^* \tilde{u}_i - U_j \tilde{u}_i^* \frac{\partial \tilde{u}_i}{\partial x_j} - \tilde{u}_i^* \tilde{u}_j \frac{\partial U_i}{\partial x_j} - \iota k_x \tilde{u}^* \tilde{p} - \tilde{u}_j^* \frac{\partial \tilde{p}}{\partial x_j} - k_x^2 \left( \frac{1}{Re} + v_t \right) \tilde{u}_i^* \tilde{u}_i + \left( \frac{1}{Re} + v_t \right) \tilde{u}_i^* \frac{\partial^2 \tilde{u}_i}{\partial x_j \partial x_j} + \tilde{u}_i^* \frac{dv_t}{dy} \frac{\partial \tilde{u}_i}{\partial y} + \iota k_x \tilde{u}^* \tilde{v} \frac{dv_t}{dy} + \tilde{u}_j^* \frac{dv_t}{dy} \frac{\partial \tilde{v}}{\partial x_j}. \quad (4.1) $$

Moreover, the following identities are considered:

$$ -\iota k_x \tilde{u}^* \tilde{p} - \tilde{u}_j^* \frac{\partial \tilde{p}}{\partial x_j} = -\frac{\partial}{\partial x_j} \left( \tilde{u}_j^* \tilde{p} \right) \quad (\text{using the continuity equation}); \quad (4.2) $$

$$ \frac{1}{Re} \tilde{u}_i^* \frac{\partial^2 \tilde{u}_i}{\partial x_j \partial x_j} = \frac{1}{Re} \frac{\partial}{\partial x_j} \left( \tilde{u}_i^* \frac{\partial \tilde{u}_i}{\partial x_j} \right) - \frac{1}{Re} \frac{\partial \tilde{u}_i^*}{\partial x_j} \frac{\partial \tilde{u}_i}{\partial x_j}; \quad (4.3) $$

$$ v_t \tilde{u}_i^* \frac{\partial^2 \tilde{u}_i}{\partial x_j \partial x_j} = \frac{\partial}{\partial x_j} \left( v_t \tilde{u}_i^* \frac{\partial \tilde{u}_i}{\partial x_j} \right) - v_t \frac{\partial \tilde{u}_i^*}{\partial x_j} \frac{\partial \tilde{u}_i}{\partial x_j} - \tilde{u}_i^* \frac{dv_t}{dy} \frac{\partial \tilde{u}_i}{\partial y}. \quad (4.4) $$

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

$$ \sigma = T_a + \mathcal{P} - \mathcal{D} - \mathcal{D}^c + C, \quad (4.5) $$

where

$$ T_a = -\iota k_x \underbrace{\int_{\Omega} U \tilde{u}_i^* \tilde{u}_i \, d\Omega}_{T_{a1}} - \underbrace{\int_{\Omega} U_j \tilde{u}_i^* \frac{\partial \tilde{u}_i}{\partial x_j} \, d\Omega}_{T_{a2}}; \quad (4.6) $$

$$ \mathcal{P} = - \underbrace{\int_{\Omega} \tilde{u}^* \tilde{v} \frac{\partial U}{\partial y} \, d\Omega}_{\mathcal{P}_{uy}} - \underbrace{\int_{\Omega} \tilde{u}^* \tilde{w} \frac{\partial U}{\partial z} \, d\Omega}_{\mathcal{P}_{uz}} - \underbrace{\int_{\Omega} \tilde{v}^* \tilde{v} \frac{\partial V}{\partial y} \, d\Omega}_{\mathcal{P}_{vy}} - \underbrace{\int_{\Omega} \tilde{v}^* \tilde{w} \frac{\partial V}{\partial z} \, d\Omega}_{\mathcal{P}_{vz}}; \quad (4.7) $$

$$ - \underbrace{\int_{\Omega} \tilde{w}^* \tilde{v} \frac{\partial W}{\partial y} \, d\Omega}_{\mathcal{P}_{wy}} - \underbrace{\int_{\Omega} \tilde{w}^* \tilde{w} \frac{\partial W}{\partial z} \, d\Omega}_{\mathcal{P}_{wz}}; \quad (4.8) $$

$$ \mathcal{D} = \frac{k_x^2}{Re} \int_{\Omega} \tilde{u}_i^* \tilde{u}_i \, d\Omega + \frac{1}{Re} \int_{\Omega} \frac{\partial \tilde{u}_i^*}{\partial x_j} \frac{\partial \tilde{u}_i}{\partial x_j} \, d\Omega; \quad (4.9) $$

$$ \mathcal{D}^c = k_x^2 \int_{\Omega} v_t \tilde{u}_i^* \tilde{u}_i \, d\Omega + \int_{\Omega} v_t \frac{\partial \tilde{u}_i^*}{\partial x_j} \frac{\partial \tilde{u}_i}{\partial x_j} \, d\Omega; \quad (4.10) $$

$$ C = \iota k_x \underbrace{\int_{\Omega} \tilde{u}^* \tilde{v} \frac{dv_t}{dy} \, d\Omega}_{C_1} + \underbrace{\int_{\Omega} \frac{dv_t}{dy} \tilde{u}_j^* \frac{\partial \tilde{v}}{\partial x_j} \, d\Omega}_{C_2}. \quad (4.11) $$

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:

$$ T_p = - \int_{\Omega} \frac{\partial}{\partial x_j} (\tilde{u}_j^* \tilde{p}) \, d\Omega \equiv 0; \quad (4.12) $$

$$ T_d = \frac{1}{Re} \int_{\Omega} \frac{\partial}{\partial x_j} \left( \tilde{u}_i^* \frac{\partial \tilde{u}_i}{\partial x_j} \right) d\Omega \equiv 0; \quad (4.13) $$

$$ T_d^c = \int_{\Omega} \frac{\partial}{\partial x_j} \left( v_t \tilde{u}_i^* \frac{\partial \tilde{u}_i}{\partial x_j} \right) d\Omega \equiv 0. \quad (4.14) $$

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 $\mathrm{Im}(\sigma)$ but does not contribute to the growth rate $\mathrm{Re}(\sigma)$ (as expected for a transport term); (ii) $\mathcal{D}$, $\mathcal{D}^c$ and $C$ are purely real with $\mathcal{D} \ge 0$ and $\mathcal{D}^c \ge 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><img src="image_21.png" /><figcaption>Figure 11. Energy budget for the leading eigenmode at $Re_τ = 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.</figcaption></figure>

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 $Re_τ = 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 $P_{uy}$, $P_{uz}$ and $D$. The remaining production terms are several orders of magnitude smaller. For the eddy viscosity case, the eddy viscosity dissipation $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 $v_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.

## 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, $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 $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 $Re_τ$. Conversely, for the eddy viscosity model (panel (d)), the eddy viscosity dissipation increases faster than the production when the $Re_τ$ is increased. Consequently, at some critical value of $Re_τ$ the eddy viscosity dissipation takes over the production and stabilises the eigenmode. This suggests a $Re$ dependence for the growth rates encoded in the $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, $P_{uy}$ and $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 > 0.4$.
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Figure 12. Energy budget for the leading eigenmode with $k_x = 0.18$. (a–b) Variation with streak amplitude $A_s$ for fixed $Re_τ = 71$. (c–d) Variation with $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.</figcaption></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>Figure 13. Energy budget for the leading eigenmode as $k_x$ is varied. $Re_τ = 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.</figcaption></figure>

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 $D^c$.

## 5. Discussion

The streak field selected in this study corresponds to the most amplified structure, computed using resolvent analysis after linearising around the mean flow. An eigenvalue problem was formulated by linearising around a synthetic base flow, constructed as the mean flow, to which the optimal streaks have been added with a tuneable amplitude. It was solved numerically by allowing for detuned eigenmodes, i.e. spatially subharmonic modes
whose wavelength is not a simple multiple of the spanwise wavelength of the streaks. Visual comparison between the most unstable large-scale eigenmode and the structure of laminar–turbulent patterns, as observed in numerical simulations, suggests that this instability is indeed relevant to the appearance of the pattern.

The starting hypothesis for this approach differs greatly from that of Liu & Gayme (2021), who have considered linearisation around the analytical laminar base flow solution. The streak instability considered here takes place in a turbulent environment whose effect needs to be modelled. In fact, this work shows that the turbulent fluctuations contribute to the sought instability in two ways. First, the coherent part of the fluctuations (the streaks) make the flow unstable to various wavelengths. Secondly, the turbulent fluctuations, via the eddy viscosity term, damp high wavelengths selecting specific large-scale modes, thereby introducing a cutoff in wavenumber. As a result of this competition, large-scale modes can be destabilised at low enough *Re* provided *A*<sub>s</sub> is large enough. This result highlights the need to correctly model unresolved turbulent motions in linear studies of turbulence, in agreement with other recent works (Illingworth *et al.* 2018; Morra *et al.* 2019; Symon *et al.* 2023). The present model is an improvement compared with the approach in Kashyap *et al.* (2024), which featured the eddy viscosity in the governing equations but not the streak mode in the base flow.

Another important assumption of the model is the high degree of symmetry of the optimal streaks obtained from resolvent analysis. By construction, the resolvent modes are harmonic in *z* and uniform in *x*. As a consequence of the harmonicity in *z* they do not modify the mean flow. In particular, the typical 100 wall units spacing of the streaks was chosen as the spanwise wavelength for the resolvent analysis. In Appendix C, it is shown that this assumption can be improved in future studies by choosing also a non-zero streamwise wavenumber.

The values of the critical amplitude reported earlier deserve a discussion in connection with the streak amplitude found either in the literature on streak instability or in DNSs of turbulent flows. The analysis of the energy budgets in § 4 shows that the desired Reynolds dependence is induced by the presence of an eddy viscosity dissipation. However, the eddy viscosity dissipation implies also an increase of the critical amplitude of the instability. This fact was already observed by Alizard (2015), who compared his critical amplitudes with those of Schoppa & Hussain (2002), after conversion of his values to the strength factor introduced by the latter. Alizard (2015) found significantly greater critical strength factors with respect to those from Schoppa & Hussain (2002) and attributed this difference to his choice of eddy viscosity. Despite the very different *Re* values, we can directly compare our critical values with those of Alizard (2015) since we use the same definition of streak amplitude. We find critical amplitudes that are correspondingly larger (≈0.40 versus ≈0.18). Therefore we further investigate this difference using DNS data. Direct numerical simulations in moderate-size periodic domains (*L*<sub>*x*</sub> × *L*<sub>*y*</sub> × *L*<sub>*z*</sub> = 35 × 2 × 15) have been performed and the following quantities introduced:

$$A_{s,xav} = \frac{\max_{y,z} \langle u' \rangle_x - \min_{y,z} \langle u' \rangle_x}{2\bar{U}_c}, \qquad (5.1)$$

$$A_{s,xmax} = \max_x \frac{\max_{y,z} u' - \min_{y,z} u'}{2\bar{U}_c}, \qquad (5.2)$$

where *u'* is the streamwise velocity fluctuation, $\langle \cdot \rangle_x$ denotes the streamwise average and $\bar{U}_c$ is the centreline velocity of the turbulent mean profile. Also, $A_{s,xav}$ measures the amplitude of streamwise-uniform streaks whereas $A_{s,xmax}$ measures the maximum local amplitude of generic streaks. The smoothing effect of the streamwise average implies
<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>Figure 14. Histograms of streak amplitudes measured from DNS in a domain of size $L_x \times L_y \times L_z = 35 \times 2 \times 15$ for (a–b) $Re_τ = 71$ and (c–d) $Re_τ = 106$. (a–c) Amplitudes from streamwise-averaged fields (5.1); (b–d) streamwise maximum amplitudes of non-averaged velocity fields (5.2).</figcaption></figure>

$A_{s,xav} < A_{s,xmax}$. The distribution of 50 000 temporal samples of these quantities is shown in figure 14 for two $Re$. Here, $A_{s,xav}$ takes values between 0.10 and 0.30 whereas $A_{s,xmax}$ takes values between 0.35 and 0.50. We conclude that, for the streamwise-uniform streaks of our base flow, a critical amplitude 0.4 is too high in average. We do not exclude that such an amplitude could be reached locally within a turbulent flow, thereby making a connection with extreme events (Hack & Schmidt 2021). The need for large $A_s$ may be a consequence of the various approximations made in the construction of the base flow, notably its over-symmetrisation. It is not excluded that, if the base flow streaks were three-dimensional or characterised by more than one spanwise wavelength, the critical amplitude would decrease, yet this remains to be verified. Appendix D documents the results obtained using a slightly modified base flow. It is shown that the instability is retrieved by keeping $A_s$ constant and increasing just the amplitude of the transverse velocity components. This suggests that improving the modelling of the base flow can lower the critical streak amplitude.

Notwithstanding the approximations made, the results of the eddy viscosity model support the initial expectation of this work, namely that large-scale structures can emerge out of the instability of a base flow featuring smaller-scale structures. This idea is illustrated in figure 15, which shows the time-averaged pre-multiplied energy spectrum of the streamwise velocity component, integrated along the wall-normal direction and obtained from a DNS at $Re_τ = 71$ in a large domain with $L_x \times L_z = 250 \times 125$. The flow for $Re$ below critical is characterised by two scales: the streaks and the modulations. There is no marked scale separation between the two, but the two peaks are clearly discernible in figure 15. The peak corresponding to the streaks is centred around the spanwise wavelength $\lambda_{z,s}^+ = 100$. Streaks also have a characteristic streamwise wavelength yet we have supposed $k_x = 0$ as a first approximation. The modulations are characterised by $k_x \approx 0.1$ and $k_z \approx 0.5$, consistently with previous DNS studies (Kashyap et al. 2020, 2022). In figure 15,
<figure><img src="image_32.png" /><figcaption>Figure 15. Wall-normal integrated pre-multiplied energy spectrum for the streamwise velocity component, time-averaged from DNS at $Re_τ = 71$. The simulation domain is $L_x \times L_y \times L_z = 250 \times 2 \times 150$. Black line denotes the structures with spanwise wavelength $λ_z^+ = 100$ (as for the base flow streaks considered in this work). Black stars denote the large-scale wavelengths of all the unstable modes obtained using the eddy viscosity model, $A_s = 0.5$, and having large-scale energy ratio $r_{LS} \ge 2\%$ (see text for detail).</figcaption></figure>

the wavelengths of the unstable modes obtained with the eddy viscosity model at $Re_τ = 71$ are superimposed on the DNS spectrum, all unstable modes with $r_{LS} > 2\%$ being selected. The range $k_x \le 0.5$ was considered, but unstable modes with $r_{LS} > 2\%$ were found only for $k_x \in [0.15, 0.4]$. The spanwise wavelength of the eigenmode was extracted using spanwise Fourier transform. The unstable modes are found in a region of the spectrum close to where the large-scale turbulent modulations lie. Closer inspection reveals that they cover very well the spanwise wavenumber range of the modulations but less accurately the streamwise wavenumber range. Note that, by assumption, the base flow is invariant in $x$, therefore the unstable mode is monochromatic in $x$. As discussed above, this is a first approximation. The results may be improved provided characteristic streamwise wavelengths of the streaks are taken into account in the base flow. This improvement is not straightforward from the computational point of view but remains a stimulating perspective for future studies.

The present modelling effort is directed specifically at the analysis of the linear instability of the turbulent flow, in line with the studies of Kashyap *et al.* (2022, 2024). As such, the model is representative of the physics only in the incipient modulational regime, when the solution does not depart sensibly from the turbulent attractor. However, the link between the oblique structure shown in figure 8 and genuine laminar-turbulent patterns, although visually stimulating, is not trivial. Both the eigenmode and the base flow are three-dimensional vector fields and the result of adding them together is not simple. In particular, active parts of the eigenmode can either tame the turbulence or reinvigorate it.

Fully developed laminar-turbulent patterns remain thus outside the scope of the present model. Indeed, laminar flow is not a solution of the nonlinear model (3.4). Therefore, nonlinear simulations of the instability following (3.4) are not expected to be
physically relevant to laminar–turbulent patterning, given that the flow would not have
the possibility to develop the laminar holes observed in DNS (Kashyap et al. 2020).
Nonlinear simulations of (3.4) have been performed anyway and they confirm this
expectation: the oblique stripe saturates in amplitude but no proper laminar hole can
develop because the streaks are constantly forced everywhere in the domain (not shown).
An improved model capturing the linear instability while allowing for both laminar and
turbulent solutions is the object of current investigation. The recent efforts of Kashyap
et al. (2025) and Benavides & Barkley (2025) are a useful starting point in this sense.

**6. Conclusions**

This work is motivated by the recent observation that laminar–turbulent patterns in channel flow emerge, as the Reynolds number is decreased, from spatial modulations that themselves arise as a linear instability of the underlying turbulent flow (Kashyap et al. 2022). This top-to-bottom point of view with decreasing *Re* differs from the bottom-up description with increasing *Re* (Duguet 2024), although there is no incompatibility between the two. Our work reveals that the modulations observed in turbulent channel flow can be modelled as a detuned instability of a mean flow superimposed with streamwise streaks.

For consistency with the observations from DNS, the effect of the turbulent background
needs to be taken into account through an eddy viscosity closure. The energy budget
analysis shows that a critical Reynolds number is predicted by the model owing to
the Reynolds dependence of the eddy viscosity. This suggests that the pattern-forming
instability is inherent to the turbulent flow, in accordance with the ideas in Kashyap
et al. (2022), the difference being that some knowledge of the coherent structures of
the turbulent flow (the streaky part) is also incorporated. Modelling linear instabilities in
turbulent environments is never trivial and several assumptions need to be made, including
the choice of a closure. We have discussed some of the most delicate points in our present
formulation and have suggested directions for future improvements. Beyond the need for
better modelling, future studies should address the question of generalisation to other
wall-parallel flows like plane Couette flow, Taylor–Couette flow, Couette–Poiseuille flow
and the asymptotic suction boundary layer. At present, it is not clear whether the linear
modulational instability discovered for the channel flow is equally relevant for these other
configurations. The subject of pattern formation in turbulent shear flows remains full of
opportunities for future studies.

**Acknowledgements.** The authors would like to thank C. Cossu for useful discussions. N.C. was funded by the Italian Ministry of University and Research. Y.D. thanks the Dynfluid laboratory for its hospitality, while N.C. thanks A. Franchini for useful discussions.

**Funding.** This work was granted access to the HPC resources of TGCC under the allocation 2024-[A0152A06362] and 2025-[A0172A06362] made by GENCI.

**Declaration of interests.** The authors report no conflict of interest.

**Appendix A. Computation and validation of the streak forcing term**

The forcing **g** in the system (3.4) is numerically computed using the nonlinear time-
stepping code **channelflow** (Gibson *et al.* 2021). To detail the procedure, let us denote
the left-hand side of the first equation in (3.4) as ∂**u**/∂t + **N**(**u**). We seek **g** such that, by
construction, **U** verifies (3.3), i.e. **N**(**U**) = **g**. To approximate **N**(**U**) using a time stepper,
<figure><img src="image_33.png" /><figcaption>Figure 16. (a–b) Forcing term **g** in ((3.3)–(3.4)) (shaded contours) compared with the forcing term **f** in (2.2) (black lines: solid for positive values, dashed for negative values) for $A_s = 1 \times 10^{-4}$. (a) Wall-normal component; (b) spanwise component. The line contours have the same iso-levels as the shaded contours to make the comparison meaningful. (c) Comparison of the streak forcing amplitude ($A_{f,i} = (\max_{y,z} f_i - \min_{y,z} f_i)/2$) as a function of the streak amplitude ($A_s$). Here, **f** and **g** are defined as in panels (a–b).</figcaption></figure>

and only for this purpose, we consider the initial value problem

$$ egin{cases} 
\displaystyle \frac{\partial \mathbf{u}}{\partial t} = -\mathcal{N}(\mathbf{u}), \\ 
\displaystyle \nabla \cdot \mathbf{u} = 0, \\ 
\mathbf{u}(t=0) = \mathbf{U}. 
\end{cases} 
\tag{A1} $$

Advancing this system of one time step $\Delta t$, one obtains a velocity field $U^1 \neq U$. Then, the forcing term $g = \mathcal{N}(U)$ can be approximated by the finite difference $-(U^1 - U)/\Delta t$.
<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>Figure 17. Convergence of the eigenvalues with the number $N_u$ of base flow units (eddy viscosity model). (a) Eigenvalues for $Re_τ = 71$, $k_x = 0.18$ and $A_s = 0.5$. (b) Leading growth rate as a function of Reynolds number for $k_x = 0.18$ and $A_s = 0.5$.</figcaption></figure>

This approach comes with a numerical error due to the finite-difference approximation. In practice, however, it turns out accurate enough. Indeed, advancing the system (3.4) with the time-stepping code starting from $U$, this field remains steady up to a small numerical error, as expected by construction of the forcing. Moreover, simulations of the linear instability of the streaky base flow have been performed using the nonlinear code and the growth rate predicted by the stability analysis was retrieved up to a <5% error (not shown). Therefore, the linear and nonlinear codes are consistent and the numerical procedures are accurate enough for the purposes of this work.

Another way to validate the forcing computed with the time-stepper code is to consider small $A_s$. In the limit $A_s \to 0$, (2.2) and (3.3) coincide with $u' = u_s$ and $g \to f$ because the difference between the equations $A_s^2 u_s \cdot \nabla u_s$ becomes negligible. This is shown in panels (a) and (b) of figure 16 for $A_s = 1 \times 10^{-4}$. As $A_s$ increases, the two forcings differ, especially on the streamwise component (figure 16c). This difference at finite $A_s$ motivates the need for a specific procedure to compute $g$.

## Appendix B. Convergence of the eigenvalues

In the main text, the dependence of the stability results on key parameters like the streak amplitude $A_s$, the Reynolds number $Re_τ$ and the streamwise wavenumber $k_x$ was thoroughly analysed. The stability analysis, however, involves also numerical parameters like the number of collocation points for the discretisation of the problem ($N_y$ in the wall-normal direction and $N_z$ in the spanwise direction) and the number of repeated base flow units which effectively imposes the spanwise size of the domain ($N_u$). This appendix shows that, as expected, the results are independent of these parameters when their values are sufficiently high.

Figure 17(a) shows the effect of increasing the number of independent base flow units on the spectrum: the branches of detuned modes tend towards continuous branches. However, as demonstrated by figure 17(b), the leading eigenvalue (and corresponding mode) is already well captured with $N_u = 50$, which means that the branch is sufficiently well discretised for this number of units.

In a similar way, figure 18(a) shows that doubling the number of collocation points in both spatial directions affects minimally the eigenvalues. The same is true for the spatial structure of the eigenmodes (not shown). Moreover, figure 18(b) shows that the $Re$ dependence of the leading growth rate is not affected either.
<figure><img src="image_36.png" /><figcaption>Figure 18. Convergence of the eigenvalues with the number of collocation points (eddy viscosity model). (a) Eigenvalues for $Re_τ = 71$, $k_x = 0.18$ and $A_s = 0.5$. (b) Leading growth rate as a function of $Re$ for $k_x = 0.18$ and $A_s = 0.5$.</figcaption></figure>

Therefore, the choice of $N_y = 65$, $N_z = 20$ and $N_u = 50$ employed for all the computations analysed in the article does not affect the conclusions of this work.

## Appendix C. Three-dimensional streaks

In the main text, streaks are computed by maximising the resolvent norm for $k_x = 0$ and $\lambda_{z,s}^+ = 100$. However, if the maximum amplification factor (cf. (2.4))

$$R_{max}(k_x, k_z) = \max_{\omega} R(k_x, k_z, \omega), \qquad (C1)$$

is considered as a function of $\lambda_z = 2\pi/k_z$ for $k_x = 0$, then the maximum amplification does not correspond to $\lambda_z^+ = 100$ (see figure 19a). Indeed, the curve is characterised by only one peak whose wavelength scales in outer rather than inner units. In contrast, at higher $Re_τ$, there are two peaks, one of which scales in wall units and corresponds to $\lambda_z^+ \approx 100$ (see the results of Hwang & Cossu (2010b), reproduced successfully, not shown here). As $Re_τ$ is decreased the peak in wall units disappears and the curve displays only one peak which does not match 100 wall units. This behaviour was also documented in other studies of optimal harmonic forcing at relatively low $Re$ (Hwang & Cossu 2010a; Willis, Hwang & Cossu 2010; Pujals, Cossu & Depardon 2010). One possible explanation is the strong hypothesis $k_x = 0$. Streaks are elongated in the streamwise direction, hence $k_x = 0$ is a decent first approximation. Streaks in turbulent flow, however, are characterised by a distribution of (finite) streamwise wavelengths which, according to several observations, is centred around $\lambda_x^+ \approx 1000$ wall units (Smits, McKeon & Marusic 2011). Figure 19(b) shows that, when the optimal amplification factor curves are shown for $k_x^+ = 2\pi/1000$, the peak at $\lambda_z^+ = 100$ is recovered. For the purpose of this work, this implies that the choice of a finite streamwise wavelength for the streaks would be more relevant. Moreover, the value of 100 inner units represents in real flows only the mean spanwise wavelength (Kline et al. 1967). An alternative would be the generalisation to a three-dimensional base flow featuring three-dimensional streaks. Such a choice is technically possible, however, it would make the stability computations more expensive and would require different numerical methods. It is therefore outside the scope of the present work.

## Appendix D. Instability of a modified base flow

The leading resolvent modes have energy predominantly in the streamwise velocity component and almost no energy in the transverse components (Hwang & Cossu 2010b).
<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /><figcaption>Figure 19. Maximum energy amplification factor as a function of the spanwise wavelength of the forcing (in wall units, $\lambda_z^+$) for streamwise wavenumber (a) $k_x = 0$ and (b) $k_x = 2\pi/1000$.</figcaption></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /><figcaption>Figure 20. Effect of increasing the transverse components amplitude in the base flow. (a) Leading growth rate as a function of $Re$ for $k_x = 0.18$. Here, $A_s$ is kept constant at 0.2 but the transverse velocity components of the base flow $V$ and $W$ are multiplied by a factor $B_{vw}$. (b) Variation of the energy budgets with the prefactor $B_{vw}$ ($A_s = 0.2$, $Re_\tau = 71$, $k_x = 0.18$). For the meaning of the labels see (4.5)–(4.11).</figcaption></figure>

However, in actual turbulent flows, the velocity field possesses three genuine components, at least because in the SSP picture streamwise streaks are accompanied by streamwise vortices. In order to analyse the influence of this three-dimensionality on critical thresholds, a numerical experiment was performed (using the eddy viscosity model) by augmenting the amplitude of transverse components of the vector field $\mathbf{u}_s$ while $A_s$ was kept constant. This is easily implemented since the streamwise invariance of the field $\mathbf{u}_s$ implies that streamwise and transverse components are decoupled in the continuity equation. Therefore if, after the rescaling of $\mathbf{u}_s$, $v_s$ and $w_s$ are multiplied by a factor $B_{vw}$, the resulting velocity field is still divergence free. The dependence of the leading growth rate on $Re_\tau$ is shown in figure 20(a) for different values of $B_{vw} \ge 1$. When $B_{vw}$ is large enough, a critical Reynolds number appears exactly as in figure 4(b). Thus, with a slightly modified base flow, the right $Re$ dependence of the growth rates can be retrieved also for $A_s = 0.2$, which is comparable to other streak instabilities critical amplitudes (Park et al. 2011; Alizard 2015; Ciola et al. 2024). Importantly, the augmented amplitudes of the transverse components are of the order of the mean amplitudes of the transverse fluctuations observed in DNS (not shown). Therefore, the modified base flow does not contain any extreme fluctuations in contrast to the previous base flow with $A_s = 0.5$. The energy budgets obtained using the new base flow (figure 20b) show
similar dissipation trends with respect to those in figure 12(b) but different trends on the production contributions, with the production term $P_{vz}$ playing a non-negligible role for large $B_{vw}$. This numerical experiment shows that certain quantitative aspects of the stability analysis depend on explicit choices in the modelling of the base flow. By contrast, the results highlighted in the main text are qualitatively robust to minor modifications of the base flow.

REFERENCES

ALBENSOEDER, S., KUHLMANN, H.C. & RATH, H.J. 2001 Three-dimensional centrifugal-flow instabilities in the lid-driven-cavity problem. *Phys. Fluids* **13** (1), 121–135.

ALIZARD, F. 2015 Linear stability of optimal streaks in the log-layer of turbulent channel flows. *Phys. Fluids* **27** (10).

ANDERSSON, P., BRANDT, L., BOTTARO, A. & HENNINGSON, D.S. 2001 On the breakdown of boundary layer streaks. *J. Fluid Mech.* **428**, 29–60.

BENAVIDES, S.J. & BARKLEY, D. 2025 Model for transitional turbulence in a planar shear flow. *Proc. Royal Soc. A* **481** (2323), 20250391.

CERQUEIRA, S. & SIPP, D. 2014 Eigenvalue sensitivity, singular values and discrete frequency selection mechanism in noise amplifiers: the case of flow induced by radial wall injection. *J. Fluid Mech.* **757**, 770–799.

CESS, R.D. 1958 A survey of the literature on heat transfer in turbulent tube flow. *Res. Rep.* 8–0529.

CHANTRY, M., WILLIS, A.P. & KERSWELL, R.R. 2014 Genesis of streamwise-localized solutions from globally periodic traveling waves in pipe flow. *Phys. Rev. Lett.* **112** (16), 164501.

CIOLA, N., DE PALMA, P., ROBINET, J.-C. & CHERUBINI, S. 2024 Large-scale coherent structures in turbulent channel flow: a detuned instability of wall streaks. *J. Fluid Mech.* **997**, A18.

CVITANOVİĆ, P. 2013 Recurrent flows: the clockwork behind turbulence. *J. Fluid Mech.* **726**, 1–4.

DEL ALAMO, J.C. & JIMENEZ, J. 2006 Linear energy amplification in turbulent channels. *J. Fluid Mech.* **559**, 205–213.

DUGUET, Y. 2024 Puffing along. *Nat. Phys.* **20** (8), 1227–1227.

DUGUET, Y. & SCHLATTER, P. 2013 Oblique laminar-turbulent interfaces in plane shear flows. *Phys. Rev. Lett.* **110** (3), 034502.

ELLINGSEN, T. & PALM, E. 1975 Stability of linear flow. *Phys. Fluids* **18** (4), 487–488.

GIBSON, J. et al. et al. 2021 Channelflow2.0. Available at https://www channelflow ch/.

DE, G., MATTEO, S., JIN, H. & HWANG, Y. 2017 Streak instability in turbulent channel flow: the seeding mechanism of large-scale motions. *J. Fluid Mech.* **832**, 483–513.

GOMÉ, S., TUCKERMAN, L.S. & BARKLEY, D. 2023 Patterns in transitional shear turbulence. Part 2. Emergence and optimal wavelength. *J. Fluid Mech.* **964**, A17.

HACK, M.J.P. & ZAKI, T.A. 2014 Streak instabilities in boundary layers beneath free-stream turbulence. *J. Fluid Mech.* **741**, 280–315.

HACK, M.J.P. & SCHMIDT, O.T. 2021 Extreme events in wall turbulence. *J. Fluid Mech.* **907**, A9.

HAMILTON, J.M., KIM, J. & WALEFFE, F. 1995 Regeneration mechanisms of near-wall turbulence structures. *J. Fluid Mech.* **287**, 317–348.

HERBERT, T. 1988 Secondary instability of boundary layers. *Annu. Rev. Fluid Mech.* **20**, 487–526.

HWANG, Y. & COSSU, C. 2010a Amplification of coherent streaks in the turbulent couette flow: an input-output analysis at low Reynolds number. *J. Fluid Mech.* **643**, 333–348.

HWANG, Y. & COSSU, C. 2010b Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow. *J. Fluid Mech.* **664**, 51–73.

ILLINGWORTH, S.J., MONTY, J.P. & MARUSIC, I. 2018 Estimating large-scale structures in wall turbulence using linear models. *J. Fluid Mech.* **842**, 146–162.

INUBUSHI, M., TAKEHIRO, S.-I. & YAMADA, M. 2015 Regeneration cycle and the covariant lyapunov vectors in a minimal wall turbulence. *Phys. Rev. E* **92** (2), 023022.

ISHIKAWA, T., TAKEHIRO, S.-I. & YAMADA, M. 2018 An orbital instability of minimal plane couette turbulence. *Phys. Fluids* **30** (3).

JIMÉNEZ, J. 2013 Near-wall turbulence. *Phys. Fluids* **25** (10).

JOUIN, A., CIOLA, N., CHERUBINI, S. & ROBINET, J.-C. 2024 Detuned secondary instabilities in three-dimensional boundary-layer flow. *Phys. Rev. Fluids* **9** (4), 043901.

KASHYAP, P.V., DUGUET, Y. & DAUCHOT, O. 2024 Linear stability of turbulent channel flow with one-point closure. *Phys. Rev. Fluids* **9** (6), 063906.
KASHYAP, P.V., DUGUET, Y. & DAUCHOT, O. 2020 Flow statistics in the transitional regime of plane channel flow. *Entropy* **22** (9), 1001.

KASHYAP, P.V., DUGUET, Y. & DAUCHOT, O. 2022 Linear instability of turbulent channel flow. *Phys. Rev. Lett.* **129** (24), 244501.

KASHYAP, P.V., MARÍN, J.F., DUGUET, Y. & DAUCHOT, O. 2025 Laminar-turbulent patterns in shear flows: evasion of tipping, saddle-loop bifurcation, and log scaling of the turbulent fraction. *Phys. Rev. Lett.* **134** (15), 154001.

KAWAHARA, G., UHLMANN, M. & VAN VEEN, L. 2012 The significance of simple invariant solutions in turbulent flows. *Annu. Rev. Fluid Mech.* **44** (1), 203–225.

KLINE, S.J., REYNOLDS, W.C., SCHRAUB, F.A. & RUNSTADLER, P.W. 1967 The structure of turbulent boundary layers. *J. Fluid Mech.* **30** (4), 741–773.

LANDAHL, M.T. 1980 A note on an algebraic instability of inviscid parallel shear flows. *J. Fluid Mech.* **98** (2), 243–251.

LIU, C. & GAYME, D.F. 2021 Structured input-output analysis of transitional wall-bounded flows. *J. Fluid Mech.* **927**, A25.

LIU, T., SEMIN, B., GODOY-DIANA, R. & WESFREID, J.E. 2024 Lift-up and streak waviness drive the self-sustained process in wall-bounded transition to turbulence. *Phys. Rev. Fluids* **9** (3), 033901.

LOZANO-DURÁN, A., CONSTANTINOU, N.C., NIKOLAIDIS, M.-A. & KARP, M. 2021 Cause-and-effect of linear mechanisms sustaining wall turbulence. *J. Fluid Mech.* **914**, A8.

MCKEON, B.J. & SHARMA, A.S. 2010 A critical-layer framework for turbulent pipe flow. *J. Fluid Mech.* **658**, 336–382.

MOARREF, R., SHARMA, A.S., TROPP, J.A. & MCKEON, B.J. 2013 Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels. *J. Fluid Mech.* **734**, 275–316.

MORRA, P., SEMERARO, O., HENNINGSON, D.S. & COSSU, C. 2019 On the relevance of Reynolds stresses in resolvent analyses of turbulent wall-bounded flows. *J. Fluid Mech.* **867**, 969–984.

PARANJAPE, C.S., DUGUET, Y. & HOF, B. 2020 Oblique stripe solutions of channel flow. *J. Fluid Mech.* **897**, A7.

PARANJAPE, C.S., YALNIZ, G., DUGUET, Y., BUDANUR, N.B. & HOF, B. 2023 Direct path from turbulence to time-periodic solutions. *Phys. Rev. Lett.* **131** (3), 034002.

PARK, J., HWANG, Y. & COSSU, C. 2011 On the stability of large-scale streaks in turbulent couette and poiseuille flows. *Comptes Rendus Mécanique* **339** (1), 1–5.

PUJALS, G., COSSU, C. & DEPARDON, S. 2010 Forcing large-scale coherent streaks in a zero-pressure-gradient turbulent boundary layer. *J. Turbul.* **11**, N25.

PUJALS, G., GARCÍA-VILLALBA, M., COSSU, C. & DEPARDON, S. 2009 A note on optimal transient growth in turbulent channel flows. *Phys. Fluids* **21** (1).

REETZ, F., KREILOS, T. & SCHNEIDER, T.M. 2019 Exact invariant solution reveals the origin of self-organized oblique turbulent-laminar stripes. *Nat. Commun.* **10** (1), 2277.

REYNOLDS, W.C. & TIEDERMAN, W.G. 1967 Stability of turbulent channel flow, with application to malkus's theory. *J. Fluid Mech.* **27** (2), 253–272.

REYNOLDS, W.C. & HUSSAIN, A.K.M.F. 1972 The mechanics of an organized wave in turbulent shear flow. Part 3. Theoretical models and comparisons with experiments. *J. Fluid Mech.* **54** (2), 263–288.

SCHMID, P.J. & BRANDT, L. 2014 Analysis of fluid systems: stability, receptivity, sensitivity: lecture notes from the flow-nordita summer school on advanced instability methods for complex flows, stockholm, sweden, 2013. *Appl. Mech. Rev.* **66** (2), 024803.

SCHMID, P.J., DE, P., FOSAS, M. & PEAKE, N. 2017 Stability analysis for n-periodic arrays of fluid systems. *Phys. Rev. Fluids* **2** (11), 113902.

SCHMID, P.J., HENNINGSON, D.S. & JANKOWSKI, D.F. 2002 Stability and transition in shear flows. *Appl. Mathé. Sci. Appl. Mech. Rev.* **142-55** (3), B57–B59.

SCHOPPA, W. & HUSSAIN, F. 2002 Coherent structure generation in near-wall turbulence. *J. Fluid Mech.* **453**, 57–108.

SHARMA, A.S. & MCKEON, B.J. 2013 On coherent structure in wall turbulence. *J. Fluid Mech.* **728**, 196–238.

SHIMIZU, M. & MANNEVILLE, P. 2019 Bifurcations to turbulence in transitional channel flow. *Phys. Rev. Fluids* **4** (11), 113903.

SMITS, A.J., MCKEON, B.J. & MARUSIC, I. 2011 High-Reynolds number wall turbulence. *Annu. Rev. Fluid Mech.* **43** (1), 353–375.

SYMON, S., MADHUSUDANAN, A., ILLINGWORTH, S.J. & MARUSIC, I. 2023 Use of eddy viscosity in resolvent analysis of turbulent channel flow. *Phys. Rev. Fluids* **8** (6), 064601.

SYMON, S., ROSENBERG, K., DAWSON, S.T.M. & MCKEON, B.J. 2018 Non-normality and classification of amplification mechanisms in stability and resolvent analysis. *Phys. Rev. Fluids* **3** (5), 053902.
TSUKAHARA, T., SEKI, Y., KAWAMURA, H. & TOCHIO, D. 2005 DNS of turbulent channel flow at very low reynolds numbers. In Fourth International Symposium on Turbulence and Shear Flow Phenomena, Begel House Inc.

TUCKERMAN, L.S., CHANTRY, M. & BARKLEY, D. 2020 Patterns in wall-bounded shear flows. *Annu. Rev. Fluid Mech.* **52** (1), 343–367.

TUCKERMAN, L.S., KREILOS, T., SCHROBSDORFF, H., SCHNEIDER, T.M. & GIBSON, J.F. 2014 Turbulent-laminar patterns in plane poiseuille flow. *Phys. Fluids* **26** (11).

WALEFFE, F. 1997 On a self-sustaining process in shear flows. *Phys. Fluids* **9** (4), 883–900.

WALEFFE, F. 2001 Exact coherent structures in channel flow. *J. Fluid Mech.* **435**, 93–102.

WILLIS, A.P., HWANG, Y. & COSSU, C. 2010 Optimally amplified large-scale streaks and drag reduction in turbulent pipe flow. *Phys. Rev. E* **82** (3), 036321.