{'fraction_non_alphanumeric': 0.03338685565583301, 'fraction_numerical': 0.012789989480095233, 'mean_word_length': 4.723067173637516, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Linear instability of turbulent channel flow

Laminar-turbulent pattern formation is a distinctive feature of the intermittency regime in subcritical plane shear flows. By performing extensive numerical simulations of the plane channel flow, we show that the pattern emerges from a spatial modulation of the turbulent flow, due to a linear instability. We sample over many realizations the linear response of the fluctuating turbulent field to a temporal impulse, in the regime where the turbulent flow is stable, just before the onset of the instability. The dispersion relation is constructed from the ensemble-averaged relaxation rates. As the instability threshold is approached, the relaxation rate of the least damped modes eventually reaches zero. The method allows, despite the presence of turbulent fluctuations and without any closure model, for an accurate estimation of the wavevector of the modulation at onset.

(Dated: September 14, 2022) Laminar-turbulent pattern formation is a distinctive feature of the intermittency regime in subcritical plane shear flows. By performing extensive numerical simulations of the plane channel flow, we show that the pattern emerges from a spatial modulation of the turbulent flow, due to a linear instability. We sample over many realizations the linear response of the fluctuating turbulent field to a temporal impulse, in the regime where the turbulent flow is stable, just before the onset of the instability. The dispersion relation is constructed from the ensemble-averaged relaxation rates. As the instability threshold is approached, the relaxation rate of the least damped modes eventually reaches zero. The method allows, despite the presence of turbulent fluctuations and without any closure model, for an accurate estimation of the wavevector of the modulation at onset. Turbulent channel flow is one of the most studied prototypes of inhomogeneous anisotropic turbulence. It has been evidenced, both experimentally and numerically, that at moderate flow rates -quantified by the Reynolds number Re -it exhibits a spatio-temporally intermittent regime featuring robust large-scale turbulent structures amid a laminar background [1][2][3]. The dynamical origin of such patterned turbulence in channel flow, as well as in other shear flows remains however actively debated [4][5][6][7][8][9][10][11][12][13][14].
On the lower end in Re of the coexistence regime, turbulent patches grow and split, or decay, resulting in strongly fluctuating dynamics. It was suggested that the stochastic nature of these processes, which decides whether turbulence will either spread or recede and eventually decay, could be described in the framework of non-equilibrium critical phenomena and specifically of directed percolation (DP), with the laminar state acting as the absorbing phase [15]. A major achievement of the past two decades has been to provide strong experimental and numerical evidences in favour of this scenario in a few shear flows [13,16,17]. On the theoretical side, this regime has been described by an effective one-dimensional model of fronts in an excitable medium [18,19]. Quite remarkably, and yet not theoretically understood, numerical simulations of the stochastic version of that model reproduce the DP scenario.
Increasing Re, individual turbulent patches leave place to a coexistence organized at the flow scale [4,5], in the form of a well organized periodic pattern of alternating laminar and turbulent bands, inclined at a well defined angle to the mean flow [6]. Considering the proliferation of turbulence as a problem of front propagation, it is tempting to view this pattern as packed arrays of individual localized structures. Yet, periodic pattern solutions have not been identified as solutions to the effective excitable dynamics [18,19]. An alternative viewpoint is to consider the pattern as emerging from the featureless turbulence found at larger Re. Pioneering studies have demonstrated experimentally that the pattern developing in pCf and TCf, characterised by two competing orientation of alternate sign, is fully captured by the dynamics of two coupled Ginburg-Landau equations with noise [6,7]. Recent visualizations, obtained in well-resolved numerical simulations of large domains channel flows, unveil small-amplitude harmonic modulations of the turbulent flow for values of Re larger than those at which genuine laminar-turbulent coexistence is reported [3]. Statistical signatures of low-wall-shearrate intermittency have been found, at Re-values usually associated with featureless turbulent flows [20]. Altogether these results suggest the possibility of a largewavelength instability of the turbulent flow itself, as already proposed in [6]. As recently suggested on the basis of a spatiotemporal extension of a classical selfsustained turbulence model [21,22], the instability could be of Turing type [9]. However, there is no theoretical evidence for such a linear instability of the turbulent mean flow obtained from one-point closure models [23]. These different viewpoints reflect global doubts regarding the origins of the patterned state and ambiguity as to whether the starting point for modeling the periodic pattern should be the spatial organization of the isolated turbulent patches when increasing Re, or the linear instability of the turbulent flow, including its fluctuations, when decreasing Re.
Here we bring direct evidence in favour of the linear instability scenario in the case of the channel flow. To do so, we perform extensive numerical simulations and sample the linear response of the turbulent flow to a temporal impulse, in the regime where the large wavelength modulations are damped. The dispersion relation is then constructed from the ensemble-averaged relaxation rates, for decreasing values of Re. The smallest relaxation rate approaches zero for some critical value Re c , pointing at the spatial structure of the modes which grow at the instability onset. The method can be seen as the temporal counterpart of the spatial linear response considered in [24]. It is intrinsically statistical in the sense that it establishes an average dispersion relation for the instability modes, from which the quantitative onset for the spatial modulation can be identified.
The incompressible flow considered in this study is driven in the streamwise direction x by a constant pressure gradient. The other Cartesian coordinates y and z are respectively wall-normal and spanwise. All length scales are nondimensionalised by the channel half-gap h, and velocities u = (u x , u y , u z ) by U cl , the centerline velocity of the classical laminar plane Poiseuille flow U(y) = 1 − y 2 driven by the same pressure gradient. Time is reported in units of h/U cl . The velocity field is decomposed as u = U(y)e x + u where u denotes the perturbation to the laminar base flow. Spatial averages are indicated with • x,y,z where the subscript indicates the direction over which the average is computed. Time averages are indicated by•. Ensemble averages are indicated as • e . Fourier amplitudes are denoted with•. The selected control parameter is the friction Reynolds number Re τ = u τ h/ν, where ν is the kinematic viscosity of the fluid, u τ = τ xz /ρ is the friction velocity, with τ xz the mean shear rate fixed by the pressure gradient, and ρ the fluid density. Turbulent simulations were performed with the spectral solver Channelflow2.0 [25] in a domain of L x = 2L z =250 for times up to t = 4, 000. These simulations are resolved with a resolution of N x = N z =1024 (including dealiasing with the 2/3 rule) and N y = 65 comparable to [3]. The most recent investigations have reported laminar-turbulent patterns for 50 Re τ 90, and independent turbulent bands for lower values of Re τ down to ≈ 36 [3,20,26,27].
Large-scale modulations close to Re τ ≈ 90, as well as genuine laminar-turbulent patterning for Re τ 90, are unambiguous from Fig. 1, which displays the instantaneous kinetic energy in the wall-normal direction both at full spatial resolution and after application of a low-pass filter. It was checked that the modulations and the pattern are robust with respect to the doubling and quadrupling of the numerical domain in both x and z (see Appendix-A of Supplementary material). The exact range of existence of the modulations, and notably their onset, are difficult to judge from visualizations alone because of the turbulent fluctuations, whose standard deviation can exceed the amplitude of the modulation. It is also sensitive to the choice of the visualized quantity. Conversely, the emergence of large-scale patterns, as Re τ decreases, appears clearly as a low-wavenumber signature in the time-averaged two-dimensional energy spectrum of the y-averaged fluctuating streamwise component u x y (Figure 2-a). Apart from the small-scale modes, corresponding to the turbulent fluctuations, one clearly observes a set of large-scale modes excited at Re τ = 92 (but absent at Re τ = 110). We also note an increase of the energy contained in the small-scale modes, and in the modes separating them from the large-scale ones, as Re τ decreases. This highlights the persistent role of nonlinear triadic interactions between different scales. Still, two maxima are readily identified.
Exploiting this scale separation, we define E LSF (Re τ ), and E SSF (Re τ ), the dimensionless amplitude of the large-scale and small-scale flows, as the energy content of the spectral subdomains S LSF = {0 < k z < 1} and S SSF = {1 k z }, respectively : E LSF is dominated by the low-k modes inside the spectral subdomain S highlighted in Fig. 2-a(left). We checked that this scale separation, based solely on k z , appropriately delineates the two energy peaks observed in the spectra of all the turbulent fields we analysed. E LSF and E SSF are shown as functions of Re τ in Fig. 2-c obtained from random initial conditions (R IC ) or during slow ascent, respectively descent, annealing in Re τ . One observes a clear increase of E LSF in contrast with the marginal increase of E SSF as Re τ decreases from the featureless turbulent regime (Re τ 110) to the well-defined pattern one (50 Re τ 90), with no sign of hysteresis. We note that E LSF is never strictly zero even at high Re τ . Whether the above observations result from a true bifurcation or are simply a mere crossover cannot be decided by simply looking at Figure 2. This is what motivates the following analysis where we show that the rise of E LSF is due to a linear instability of the turbulent flow.
Establishing the linear instability of a flow with arbitrary time-dependence can be addressed in different ways. One possibility is to study the linear stability of the mean flow using the Orr-Sommerfeld formalism. This strategy, whether conducted at high [28] or transitional [29] Re, predicts linear stability. At the opposite end, taking into account all temporal fluctuations is in principle possible using Lyapunov analysis. However for turbulent flows the number of positive Lyapunov exponents is prohibitively huge [30] because of the chaoticity at small scales down to the Kolmogorov scale. The turbulent scales where these instabilities dominate are however not the emerging large-scales visible in Fig. 2, which suggests the computation of alternative quantities. The general idea is to study the linear response of the flow to a temporal impulse. If the flow is linearly stable, the disturbance should relax, otherwise it should grow and lead to a bifurcated flow. However, the reference flow being turbulent, the analysis must be conducted at a statistical level. Besides the spatial structure of the temporal impulse should be agnostic to the turbulent spectrum. We therefore proceed as follows. A representative turbulent state in the statistically steady regime at the required value of Re τ , simulated for t < 0, is perturbed at t = 0 using a divergence-free noise field, before the simulations runs further without noise, for t > 0, and we monitor the temporal evolution of the modulus of the Fourier amplitudes of large-scale modes, |τ| (Re τ , k x , k z , t), with (k x , k z ) ∈ S = {0.075 k x 0.22, 0.2 k z 0.5} (the highlighted square area in Fig. 2-a, part of S LSF ). For large enough Re τ , the disturbed flow relaxes back towards the steady turbulent state. The individual time series |τ| however showcase a strongly fluctuating decay. This computational decay experiment is therefore repeated over 40 different realizations of the noise field and the modulus of the spectral amplitude of each largescale Fourier mode, |τ|, is ensemble-averaged over all realizations to yield |τ| e , as illustrated in figure 3-a for Re τ = 120 and (k x , k z ) = (0.12, 0.3). Ensemble-averaging brings clarity into the system's response: past an initially nonlinear decrease of |τ| e , a clear exponential de-cay towards a finite value A s is observed. Given that this exponential decay is relevant only at the late stage of the relaxation, changing the amplitude of the initial noise field is not relevant. This exponential decay captures the averaged linear response of the turbulent state with respect to a temporal impulse. The corresponding growth rate σ is evaluated by estimating first the saturation level A s and then fitting an exponential decay to ( |τ| e − A s ), using a straight line fit in logarithmic scale, as portrayed in figure 3-b (See Appendix B of Supplementary material for a detailed step by step description of the procedure). As a first step, the analysis is carried out along the diagonal of the spectral window S. Fig. 3c,d show the strong dependence of the growth rate on both k = k 2 x + k 2 z and Re τ . More specifically, one observes that, for Re τ = 96, the growth rate of the mode corresponding to k = 0.38 is close to vanishing, suggesting the proximity of a linear instability. In principle one could expect monitoring the average exponential growth of such a large-scale mode beyond the instability threshold. However, not only would the growth rate be hard to measure accurately near onset, one would also need to isolate the featureless turbulent state in a regime where it is unstable.
We therefore concentrate on the decay rates and extract the mean dispersion relation for the linear response of the turbulent flow (Fig. 4). The data is fit with a paraboloid surface (Fig. 4a) of the form : where (k xc , k zc ) is the critical wavevector. The coefficients of Eq 3 obtained for different values of Re τ are reported in the supplementary material Appendix B. The dispersion curves approach the neutral axis as the value of Re τ is decreased and eventually cross it for Re τ = 94. The estimated critical value for the instability is Re τ = 95±1. The critical wavevector (k xc = 0.18 ± 0.025, k zc = 0.42 ± 0.05) is obtained from the above parabolic fit of the decay rates, estimated for all (k x , k z ) ∈ S, as illustrated in Fig. 4-a for Re τ = 110 and Re τ = 96. It perfectly matches the one measured directly at onset and lead to an inclination of the pattern with the streamwise direction of 23±0.5 • , consistently with the measurements reported in [20]. This quantitative agreement validates the proposed methodology, i.e. the statistical analysis of the temporal impulse response can be considered as a new experimental/numerical method to address the linear stability analysis of a steady, but fluctuating dynamics. We emphasize again that the base flow for the analysis is the turbulent flow itself, including all fluctuations [31], not the mean flow.
Altogether our results provide direct evidence for a linear instability of the turbulent state itself, as first conjectured in Ref. [6]. This linear instability leads to a spatial modulation of the turbulent flow, the amplitude of which grows and saturates according to weakly nonlinear contributions [7]. For low enough Re τ , the modulation breaks into a pattern of alternated turbulent and laminar bands. Further decreasing Re τ these bands gain in independence and a proper stochastic front dynamics sets in.
Our work paves the way for future works in two main directions. First, one would like to identify the instability mechanism. A possible candidate, commonly encountered across diverse noisy chemical and biological systems [32] relies on the Turing instability [9,33]. It is based on the competition between an inhibitor and an activator field with different diffusivities [34]. However this approach requires modelling of the turbulent diffusivity using e.g. simple first-moment closures [35,36]. Another possible approach is to consider a generalised stability analysis taking into account higher-order moments of the fluctuations [37]. Both approaches are based on closure assumptions. The instability unveiled in the present work represents an ideal and simple case to test these assumptions. The second future direction of research consists in identifying the strongly nonlinear scenario along which the pattern looses its spatial coherence. It remains a formidable challenge. ing computational resources from IDRIS (Institut du Développement et des Ressources en Informatique Scientifique) and the support of its staff. The developing team of channelflow.ch is also gratefully thanked. The authors also acknowledge constructive discussions with D. Barkley, L.S. Tuckerman, S. Gomé and P. Manneville.

SUPPLEMENTARY INFORMATIONS Appendix A : Illustration of pattern formation for different domain sizes
The robustness of pattern formation has been verified with simulations conducted for different domain sizes with equivalent numerical resolutions. As an example, fig 5 illustrates pattern formation occurring at Re τ = 80, visualized with the help of E v (x, z) for different domain sizes i.e L x = 2L z = 250, L x = 2L z = 500, L x = 2L z = 1000. It demonstrates the robustness of the wavelength with respect to changes of computational domain. Simulations with the inverse aspect ratio, i.e 2L x = L z = 250, confirm this observation.