Laminar-Turbulent Patterns in Shear Flows : Evasion of Tipping,
Saddle-Loop Bifurcation and Log scaling of the Turbulent Fraction.

The transitional regime to laminar flow in wall-bounded shear flows takes the form of coexisting laminar and turbulent domains Manneville (2015); Tuckerman et al. (2020), resulting from a globally subcritical transition scenario Dauchot and Manneville (1997). When increasing the Reynolds number, $R$, the transition from the laminar flow to the regime of spatio-temporally intermittent turbulent spots follows a directed percolation scenario Pomeau (1986); Lemoult et al. (2016); Chantry et al. (2017); Klotz et al. (2022). Away from the transition, large aspect ratio plane Couette (pCf) and Taylor-Couette (TCf) flows, the flows sheared between two parallel planes, respectively two rotating coaxial cylinders, organize in the form of a regular laminar-turbulent pattern, as first experimentally reported in Prigent et al. (2002, 2003). These patterns have now been reported experimentally or numerically in all wall-bounded shear flows Tsukahara et al. (2014); Duguet et al. (2010); Ishida et al. (2016); Kunii et al. (2019); Paranjape (2019); Shimizu and Manneville (2019); Takeda and Tsukahara (2019); Kashyap et al. (2020), except for the circular Poiseuille flow (cPf) Moxey and Barkley (2010). In planar geometries they emerge, when decreasing R from homogeneous turbulent flow, with a wavelength $\lambda\simeq 50h$, with $h$ the thickness of the direction of mean shear, and an orientation of approximately $25^{\circ}$ with respect to the mean flow Prigent et al. (2003); Tsukahara et al. (2014); Shimizu and Manneville (2019); Kashyap et al. (2020).

According to Prigent et al. (2002), the patterns result from a linear instability of the uniform turbulent flow, allowing for a weakly nonlinear Ginzburg–Landau description Prigent et al. (2003). This hypothesis was first reinforced by the study of the statistics of the spatial Fourier component of the pattern Tuckerman et al. (2008). It recently found its confirmation with the study of the ensemble-averaged response of uniform turbulence to large-wavelength perturbations and the finding of a dispersion relation in very good agreement with the observed pattern Kashyap et al. (2022). This linear instability of the homogeneous turbulent flow also takes place in respectively two-fields Manneville (2012) and six fields Benavides and Barkley (2023) models, obtained from the Navier-Stokes equations through specific closures.

Here we tackle the remaining central issue illustrated in Fig. 1, namely how the strongly nonlinear evolution of the pattern leads, in a characteristic tipping evasion scenario Rietkerk et al. (2021), to the spatio-temporally intermittent regime when decreasing $R$. Starting with the simplest two-fields model Manneville (2012), to which advection and noise are added, we reveal rich unexplored dynamics. More specifically, we show that (i) many stable patterns coexist with wave vectors delimited by the boundaries of the so-called Busse balloon Busse (1978), (ii) a specific wavelength is selected when adding a multiplicative noise that mimics the turbulent fluctuations, (iii) a cascade of instabilities lead to a dramatic increase of the wavelength on the lower $R$ side of the balloon, eventually reaching infinite wavelength, (iv) in this limit, the periodic limit cycle associated with the spatial pattern hits the laminar fixed point into a saddle-loop global bifurcation. According to this scenario, the pattern evolves from a regular harmonic pattern to a non-harmonic one, with a logarithmic divergence of its wavelength, while the turbulent regions remain essentially of constant size. The prediction compares well with existing data from the literature, as long as the pattern is well formed, i.e. before entering the spatio-temporal intermittent regime obeying the directed percolation scenario.

The Waleffe model Waleffe (1997); Dauchot and Vioujard (2000) was decisive in understanding the local self-sustaining process along which streamwise vortices ($V$) draw the mean flow ($M$) from high to low-velocity regions leading to streamwise velocity fluctuations called streaks ($U$). The latter undergo an inflectional instability leading to modulations of amplitude ($W$), feeding back the vortices and thereby sustaining turbulence. An important step forward was to propose a spatially extended version of the model including diffusive coupling Manneville (2012). Assuming fast dynamics for the fields $U$ and $V$, the model was reduced to a two-field model, which exhibits a Turing linear instability of the homogeneous turbulent flow, when decreasing $R$. Here we further expand on this first success by adding advection and a multiplicative noise to account for the fluctuating nature of the turbulent field $W$. The model we study in one dimension then reads:

	$\displaystyle\frac{\partial M}{\partial t}+M\frac{\partial M}{\partial x}$	$\displaystyle=f(M,W)+D_{M}\frac{\partial^{2}M}{\partial x^{2}},$		(1a)
	$\displaystyle\frac{\partial W}{\partial t}+M\frac{\partial W}{\partial x}$	$\displaystyle=g(M,W)+D_{W}\frac{\partial^{2}W}{\partial x^{2}}+\sigma W\eta,$		(1b)

where $D_{M},D_{W}$ are the diffusivities of the fields $M,W$, $\sigma$ is the amplitude of the noise and $\eta$ is a Gaussian random field with zero mean and $\left<\eta(x,t)\eta(x^{\prime},t^{\prime})\right>=\delta(t-t^{\prime})\delta(x-x^{\prime})$ (see Supp. Mat. sec. I-A,B for more details on the model ). The functions $f(M,W)$ and $g(M,W)$ are given by:

	$\displaystyle f(M,W)\!$	$\displaystyle=\!-\alpha_{{}_{M}}(M-M_{0})-\beta_{{}_{M}}(M-\bar{M})W^{4}+\sigma_{{}_{M}}W^{2},$		(2a)
	$\displaystyle g(M,W)\!$	$\displaystyle=\!-\alpha_{{}_{W}}W+\left[\beta_{{}_{W}}(M-\bar{M})-\gamma_{{}_{W}}\right]W^{3}-\sigma_{{}_{M}}MW,$		(2b)

where the dependence on $R$ is encoded in the coefficients (see Supp. Mat. sec. I-B). The absorbing nature of the laminar state $(M=M_{0},W=W_{0}=0)$ emerges from the form of $f$ and $g$ and is respected by the choice of a multiplicative noise.

Figure 2 displays the vector field $(f,g)$, together with the nullclines $f(M,W)=0$ and $g(M,W)=0$, the intersection of which provide the spatially homogeneous fixed points. The laminar solution $(M=M_{0},W=W_{0}=0)$ is linearly stable for all $R$. For $R\geq R_{sn}$, two extra solutions $(M=M_{\pm},W=W_{\pm})$ emerge from a saddle node bifurcation. The lower branch $(M_{-},W_{-})$ is linearly unstable for all $R>R_{sn}$ The parameter values are such that the upper branch $(M_{+},W_{+})$ is linearly stable with respect to homogeneous perturbations for $R>R_{sn}=310$, yet destabilizes via a Turing instability, when decreasing $R$ below $R_{T}=415$ (see also Dauchot and Vioujard (2000); Manneville (2012) and Supp. Mat. sec. I.C).

The simulation of Eqs. (1) with periodic boundary conditions in a domain of size $L=1000$ (see Supp. Mat. sec. II-A) display the main features of the transitional regime in plane shear flows. When slowly annealing $R$ from the homogeneous turbulent state (Fig. 3-a), a periodic pattern emerges, the wavelength of which increases with decreasing $R$. For low enough values of $R$, the pattern is replaced by a disordered juxtaposition of isolated pulses, which eventually vanish when further decreasing $R$ below $R_{g}\simeq 200$. For $R<R_{sn}$, the model exhibits an excitable-like behaviour: when initialized with a large enough amplitude Gaussian profile of $W$, a localized pulse nucleates (Fig. 3-b). Single pulses, as well as superposition of multiple pulses, are stable. The amplitude of the initial condition necessary to trigger a pulse increases with decreasing $R$. Nucleating two pulses one after the other either leads to a fusion of the two pulses when they are too close or to two distinct pulses, with the downstream one being pushed away from the upstream one. For $R>R_{sn}$, the initial Gaussian profile expands downstream, increasing the width between a laminar-turbulent front and a turbulent-laminar one. The turbulent state between the two fronts modulates, leading to a well-formed pattern (Fig. 3-c), as also seen in Alavyoon et al. (1986); Duguet et al. (2010).

The pattern wavelength observed when applying the above annealing procedure is selected by noise. To show this, we first run deterministic simulations ($\sigma=0$) with initial conditions in the form of a harmonically modulated turbulent state. We identify the stable patterns as those which keep the same wave number $k=2\pi/\lambda$ for very long times ($T>100\,000$). The region in the $R-k$ space encompassing all the stable patterns is called the Busse balloon (Busse, 1978). It extends for $R$ up to $460>R_{T}$, retaining the subcriticality of the Turing instability reported in Manneville (2012), and down to $\simeq 240$, far below $R_{sn}$, in a characteristic evasion from the tipping point catastrophe Rietkerk et al. (2021). When initiating the noiseless simulation with a well-formed pattern and decreasing the Reynolds number, the wavelength remains constant until the pattern hits the left edge of the Busse balloon (continuous line trajectories on Fig. 4-a). Conversely, in the presence of noise, patterns initiated with two different wave numbers converge towards a narrow range of wavelength, which is also the one obtained when annealing from the homogeneous turbulent state (dashed line trajectories Fig. 4-a and spatio-temporal diagrams Fig. 4-b,c).

When $R\lesssim R_{sn}$, large wavelength patterns coexist with localized pulses in phase space. The connection between these two families of solutions is better described in the system of coordinate $z=x-ct$, where $c$ is the advection speed of the travelling solutions. One thereby obtains a 4-dimensional dynamical system:

	$\displaystyle\dot{M}$	$\displaystyle=P,$		(3a)
	$\displaystyle D_{{}_{M}}\dot{P}$	$\displaystyle=(M-c)P-f(M,W),$		(3b)
	$\displaystyle\dot{W}$	$\displaystyle=Q,$		(3c)
	$\displaystyle D_{{}_{W}}\dot{Q}$	$\displaystyle=(M-c)Q-g(M,W),$		(3d)

with $\dot{A}\equiv dA/dz$. The fixed points of this dynamical system are $(M=M_{0,\pm},P=0,W=W_{0,\pm},Q=0)$ and are all saddles with two positive eigenvalues. The pattern and the pulse solutions of Eqs. (1) respectively map onto periodic orbits around the turbulent fixed point and homoclinic orbits of the laminar one (see Fig. 5-a). When reducing $R$, the periodic orbits grow around the unstable homogeneous turbulent state (see Fig. 5-b), distort and approach the laminar saddle. The central question is whether a global saddle-loop bifurcation takes place, with the periodic orbits turning into homoclinic ones Strogatz (2000); Doelman et al. (2018). If such a scenario were to take place, it would strongly constrain the wavelength dependence on $R$ in the vicinity of the global bifurcation (Gaspard, 1990). This question can however not be answered using simulations with periodic boundary conditions, which formally only produce periodic solutions.

We thus turn to numerical continuation methods to follow the periodic orbits of Eqs. (3) and their stability Sherratt (2012) (see Sup. Mat. sec. II-B). A pattern solution at large $R\simeq 450$, which advection speed $c\simeq 0.2$ is measured, is seeded as a solution for Eqs. (3). We then perform arc-length continuation, both increasing and decreasing $R$, keeping $c$ fixed (black lines of Fig. 5-c,d), until a secondary bifurcation takes place (red and yellow bounds of Fig. 5,c-d). Starting from a stable solution, close to the lower $R$ instability, we then suggest a small increase in $c$ and iterate the procedure until the full Busse balloon, with a highly non-trivial shape at low $k$, is obtained. In particular, we find two families of solutions with arbitrarily large periods located in the two ”legs”, resp. ”tips”, in the $(k,R)$, resp. $(c,R)$, representation of the Busse balloon. Those solutions eventually connect to the homoclinic orbits starting from the laminar fixed point and therefore confirm the global saddle-loop bifurcation scenario.

Fig. 5-(c,d) also displays the trajectory followed by the pattern solution in the $(R,k)$ and $(R,c)$ parameter spaces during three independent annealing processes as the one reported in Fig. 3-(a). The pattern selected by the noise at large $R$ moves across the Busse balloon, maintaining an essentially constant wavelength until it reaches its left boundary and follows it down to larger and larger wavelength. The growth of the wavelength is then controlled by the saddle-loop bifurcation and, therefore, obey the following logarithmic scaling Gaspard (1990):

$$\lambda=-a\log\epsilon+b,$$

(4)

with $a,b>0$ and where $\epsilon=(R-R_{sl})/R_{sl}$, with $R_{sl}\simeq 250$ the critical value of $R$, where the saddle-loop bifurcation takes place.

In real flows, the saddle-loop bifurcation scenario certainly does not strictly take place, first because of finite size effects and second because other mechanisms such as longitudinal band splitting Shimizu and Manneville (2019) may take over. Nevertheless, the proximity of a global bifurcation should control the behaviour of the wavelength in a critical regime for $R\gtrsim R_{sl}$. We conclude this work by testing this hypothesis with existing experimental and numerical measurements of both the pattern wavelength and the turbulent fraction in plane Poiseuille flow (pPf) and TCf flows. The expected scaling of the turbulent fraction, $F_{t}=w_{t}/\lambda$, is obtained on the basis that the width of the turbulent bands $w_{t}$ is essentially constant in the vicinity of $R_{sl}$. In all cases, we find a remarkable agreement between the reported data and the scaling predicted in the proximity of a saddle-loop bifurcation (Eq. (4)). Note the deviation below $R_{sl}$, as the pattern breaks down and the other critical point $R_{g}<R_{sl}$ for directed percolation is approached (see e.g Klotz et al. (2022)).

Altogether, the linear instability of the homogeneous turbulent flow leads to a plethora of coexisting metastable patterns, and the fluctuations select the wavelength of the observed nonlinear pattern. The limit cycle, associated with the selected wavelength, grows, distorts, and finally hits the laminar fixed point in a global saddle-loop bifurcation, leading to a logarithmic divergence of the wavelength and, as a consequence, of the turbulent fraction. In practice, the divergence is regularized by a crossover to richer dynamics involving more than one dimension. Nevertheless, the global bifurcation and the associated singularity control the dynamics in a critical range accessible experimentally and numerically. This calls for refined simulations and experiments in the vicinity of the crossover.

The present model shares some similarities with the two-fields one-dimensional model introduced in Barkley (2011) to describe the transition to turbulence in pipe flow. Both models describe excitable dynamics, allowing for stable localized pulse solutions corresponding to the turbulent spots. However, while in Barkley (2011) the inhibiting role of turbulence comes from an increased dissipation of the turbulent energy, here it arises from reduced turbulent energy production (see Supp. Mat. sec I.D). This is actually responsible for a structurally different organization of the vector field $(f,g)$. In Barkley (2011), the nullclines have slopes of opposite signs when they cross, here the slopes have the same sign. This last property being a necessary condition for a Turing instability of the upper branch, Barkley (2011) does not exhibit laminar turbulent periodic patterns.

Finally, the scenario described above is not unique to the subcritical transition to turbulence. Reaction-diffusion is most common in many other contexts including biology, chemistry or ecology Turing (1952); Pearson (1993); Lengyel and Epstein (1992); Rietkerk et al. (2021) and generically leads to pattern formation via a linear instability, be the instability of strict Turing type or not. Making predictions beyond the weakly non-linear regime in such spatially extended systems is generally challenging. When the pattern wavelength increases away from the instability, the methodology presented here is a promising path to unveil the general scenario of a saddle-loop bifurcation Gaspard (1990); Plaza et al. (1997).