# Subcritical transition to turbulence in shear flows

Yohann Duguet

► To cite this version:

Yohann Duguet. Subcritical transition to turbulence in shear flows. Science non linéaire [physics]. Université Paris-Saclay, 2021. tel-04415087

HAL Id: tel-04415087

https://hal.science/tel-04415087

Submitted on 24 Jan 2024

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
Subcritical transition to turbulence in
shear flows

HDR

soutenue le 10.12.2021

pour l'obtention d'une

Habilitation de l'Université Paris-Saclay

(mention mécanique)

par

Yohann Duguet

Composition du jury

*Rapporteurs :* Laurette S. Tuckerman
François Charru
Peter J. Schmid

*Examinateurs :* François Daviaud
Jean-Christophe Robinet
Mis en page avec la classe thesul.
# Sommaire

<table><tr><td>Chapitre 1 Transition to turbulence in canonical shear flows</td><td>3</td></tr><tr><td>1.1 Introduction and Motivation</td><td>4</td></tr><tr><td>1.2 Governing equations</td><td>4</td></tr><tr><td>1.3 Canonical base flows</td><td>4</td></tr><tr><td>1.3.1 Planar flows</td><td>4</td></tr><tr><td>1.3.2 Cylindrical flows</td><td>6</td></tr><tr><td>1.3.3 Notion of forcing protocol</td><td>6</td></tr><tr><td>1.4 Wall turbulence at onset</td><td>7</td></tr><tr><td>1.4.1 High-Reynolds number wall turbulence</td><td>7</td></tr><tr><td>1.4.2 Estimating the onset of wall turbulence</td><td>8</td></tr><tr><td>1.4.3 Relevant interrogations regarding transition to turbulence in shear flows</td><td>10</td></tr><tr><td>1.5 Linear stability of canonical base flows</td><td>14</td></tr><tr><td>1.5.1 Transient growth</td><td>16</td></tr><tr><td>1.5.2 Nonlinear extension</td><td>19</td></tr><tr><td>1.6 Nonlinear concepts</td><td>20</td></tr><tr><td>1.6.1 Self-sustaining process</td><td>20</td></tr><tr><td>1.6.2 Finite-amplitude solutions</td><td>25</td></tr><tr><td>1.6.3 Cycle Expansions</td><td>26</td></tr><tr><td>1.6.4 Chaotic saddles</td><td>29</td></tr><tr><td>1.6.5 Other invariant sets</td><td>29</td></tr></table>

<table><tr><td>Chapitre 2 Dynamical systems description</td><td>33</td></tr><tr><td>2.1 The edge of chaos</td><td>34</td></tr><tr><td>2.1.1 The edge as a smooth separatrix</td><td>34</td></tr><tr><td>2.1.2 Concept of edge state in chaotic saddles</td><td>38</td></tr><tr><td>2.1.3 Edge tracking techniques</td><td>41</td></tr><tr><td>2.2 Edge state dynamics in terms of finite-amplitude solutions.</td><td>45</td></tr><tr><td>2.2.1 Recurrences</td><td>46</td></tr><tr><td>2.2.2 Chaoticity of edge trajectories</td><td>50</td></tr><tr><td>2.3 Minimal states</td><td>52</td></tr><tr><td>2.3.1 Theoretical concept and limitations</td><td>52</td></tr><tr><td>2.3.2 Optimal oblique transition.</td><td>53</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>2.3.3</td>
      <td>Variational technique</td>
      <td>55</td>
    </tr>
    <tr>
      <td>2.3.4</td>
      <td>Dependence on Re</td>
      <td>59</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Spatial localisation of edge states</td>
      <td>59</td>
    </tr>
    <tr>
      <td>2.4.1</td>
      <td>Motivation</td>
      <td>59</td>
    </tr>
    <tr>
      <td>2.4.2</td>
      <td>One-dimensional localisation</td>
      <td>61</td>
    </tr>
    <tr>
      <td>2.4.3</td>
      <td>Two-dimensional localisation</td>
      <td>71</td>
    </tr>
    <tr>
      <td>2.4.4</td>
      <td>plane Couette flow</td>
      <td>71</td>
    </tr>
    <tr>
      <td>2.4.5</td>
      <td>ASBL</td>
      <td>71</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Conclusions</td>
      <td>75</td>
    </tr>
    <tr>
      <td colspan="2">Chapitre 3 Spatiotemporal aspects of subcritical transition</td>
      <td>77</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Transition as a phase transition</td>
      <td>78</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>One-dimensional shear flows</td>
      <td>78</td>
    </tr>
    <tr>
      <td>3.2.1</td>
      <td>Localised puffs</td>
      <td>79</td>
    </tr>
    <tr>
      <td>3.2.2</td>
      <td>Turbulence proliferation</td>
      <td>80</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Two-dimensional shear flows</td>
      <td>86</td>
    </tr>
    <tr>
      <td>3.3.1</td>
      <td>Experimental approach</td>
      <td>86</td>
    </tr>
    <tr>
      <td>3.3.2</td>
      <td>Numerical approach</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.3.3</td>
      <td>Bifurcation diagrams</td>
      <td>92</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>The directed percolation hypothesis</td>
      <td>99</td>
    </tr>
    <tr>
      <td>3.4.1</td>
      <td>Concept and motivation</td>
      <td>99</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Large-scale flows</td>
      <td>105</td>
    </tr>
    <tr>
      <td>3.5.1</td>
      <td>Motivation</td>
      <td>105</td>
    </tr>
    <tr>
      <td>3.5.2</td>
      <td>Large/small scale coupling</td>
      <td>105</td>
    </tr>
    <tr>
      <td>3.5.3</td>
      <td>Spot structure</td>
      <td>107</td>
    </tr>
    <tr>
      <td>3.5.4</td>
      <td>Stripe development from spots</td>
      <td>109</td>
    </tr>
    <tr>
      <td>3.5.5</td>
      <td>Helical spots in annular flows</td>
      <td>113</td>
    </tr>
    <tr>
      <td>3.5.6</td>
      <td>Boundary layer flows</td>
      <td>115</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Multiphysics extensions</td>
      <td>118</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Conclusion</td>
      <td>122</td>
    </tr>
    <tr>
      <td colspan="2">Chapitre 4 A transition theory for spatially developing boundary layer flows</td>
      <td>125</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Phenomenology of Blasius boundary layer transition</td>
      <td>126</td>
    </tr>
    <tr>
      <td>4.1.1</td>
      <td>Flow set-up</td>
      <td>126</td>
    </tr>
    <tr>
      <td>4.1.2</td>
      <td>Bypass scenario</td>
      <td>127</td>
    </tr>
    <tr>
      <td>4.1.3</td>
      <td>Difficulties associated with state space picture</td>
      <td>132</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Competition between two transition scenarios</td>
      <td>133</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>Impact on edge tracking</td>
      <td>133</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Application to wind tunnel experiments</td>
      <td>137</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Low-order modelling</td>
      <td>137</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Low-order modelling using cellular automata</td>
      <td>143</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Conclusion</td>
      <td>147</td>
    </tr>
  </tbody>
</table>
This scientific subject represents my largest scientific investment since my PhD thesis and I thus chose to develop it in the rest of the present habilitation thesis. As a short informal introduction to the topic, I would like to mention how I got acquainted with the world of subcritical transition. I simply got entrained by other people who underwent that transition from rotating flows instabilities to shear flow transition, namely Rich Kerswell from Bristol University, who was apparently following the scheme proposed by Fabian Waleffe before him and got interested in the application of dynamical systems theories to these flows. I got seduced initially by the apparent mathematical formalism used in the related studies as well as by the appealing simplicity and importance of the flow configurations under study. After working on pipe flow in Bristol with R.R. Kerswell together with A.P. Willis and C.C.T. Pringle, I continued along (hopefully) the same lines at the Department of Mechanics of KTH Stockholm together with Dan S. Henningson, Philipp Schlatter and Luca Brandt, who were experts of transitional shear flows but somehow new to the nonlinear interpretation schemes. My interest, influenced directly by numerical feasibility, switched from turbulence itself to transition (via the study of "edge states") for which the dynamical systems paradigm semed more relevant... at least to me. At about the same time I got acquainted with Paul Manneville from LadHyx who convinced me that a spatiotemporal look at such flows is not an option, and suggested instead a spatiotemporal if not statistical physics-oriented perspective. Such a view gets challenged periodically by some of my colleagues from the field, including the late Bruno Eckhardt (whose scientific influence I am hardly able to quantify as of now), Predrag Cvitanović, Olivier Dauchot, Yves Pomeau, Björn Hof, John F. Gibson, Marc Avila, Laurette Tuckerman or Dwight Barkley, to name just the most productive of them. My current understanding of transition results from the interplay between these stimulating as well as often conflicting points of views, but more importantly from the scientific optimism of the people mentioned just above.
# Preamble

Here are a few words to explain the way I structured my habilitation thesis. First of all I chose to write in English because my priority targets are my own students, those gone as well as those to come. These many students are all different, yet the one thing they have in common is that they understand English. English is not my mother tongue although I use it more than French for Science, nevertheless I apologise in advance for all the pain endured by genuine Anglophones upon reading my manuscript. I chose to tell a story here, or (hopefully) give a course at an advanced scientific level. I hence chose a coherent storyline whose elements were subtle from my past published papers. Choosing to dissert on the transition to turbulence was a natural choice because that is so far the most coherent story line I can offer, and because it covers all the time period since my first post-doctoral position. The story told is almost free from technicalities. In particular, although my own everyday work (and that of my students) has long consisted of implementing numerical schemes and post-processing numerical data, I intended to zoom out and tell a physicist story. There won't be any detail about the numerical codes used, about how "details" were implemented. There are several reasons behind this choice. To start with, this is my second thesis and I wanted it to be as different as possible from the first one which was very numerics-oriented. More importantly perhaps, these details can be found in the published articles. Crucially too, including them would shadow or obscure the thinking underlying all these studies. Instead I would ideally like the potential reader to grasp the subtile connections between the mathematician's approach to the Navier-Stokes equations, the engineer's concern about keeping his pipe flow energetically interesting, and the physicist's wish to understand complexity and summarise it using a few modular concepts. Almost all images included in this document have been published somewhere and the bibliographic reference is always included in the caption or in the related text. A small number of images come from as yet unpublished studies, in which case the collaborator(s) for this image is/are also mentioned. One of the main difficulties in writing a research course, or something close to it, is that each result mentioned, be it from me or someone else, is dated. In the relatively recent field of hydrodynamics, only a very small number of results are universally accepted as being part of a basic knowledge package. This becomes especially relevant when transition to turbulence is concerned, namely a field with almost no established mathematical formulas. As a consequence I had to give in the first chapter an overview of the status of the field as I found it when I began working in it. Obviously, and I am in fact glad to have witnessed so much, during the 15 years in which I have worked on transition the field itself has mutated. The reasons are multiple: huge progress in numerical as well as experimental facilities, understanding of certains concepts evolving with time, newborn scientific analogies, new fashionable directions in this field as in other fields... It was difficult to untangle all the different contributions and to extract mine out of the context of the work of a community (since this was the initially requested homework). This let to a tendency to rewrite history in a different order, rather multiple back and forth in time. And last but not least I hope to convince you that I was among the people responsible for that mutation of the field, although most of the people mentioned in the introductory chapter bring so much more than I could ever dream of. With a certain taste for provocation or intellectual stimulation I took a guilty pleasure in revising some of the well-accepted concepts and to re-define them either more formally, or from a more general or exotic perspective. This is hence a highly subjective report, and I apologise in advance for not citing every single contribution.
1

Transition to turbulence in canonical
shear flows
## 1.1 Introduction and Motivation

Hydrodynamic stability lies at the crossroads between hydrodynamics and applied mathematics. The main reason for this intersection is two-fold: it is due first to the highly mathematical formalism developed in the pioneering studies of G.I. Taylor and co-workers, but also to the historical impact that hydrodynamic stability has had later on other fields of applied mathematics such as pattern formation or chaos theory. The fundamental mathematical property of the governing equations, here the incompressible Navier-Stokes equations, is the multiplicity of their solutions. For the same parameters, some solutions appear stable in time while others happen to be unstable in time. The dynamics of the system results from the influence of all these solutions together, and not only from the stable ones. Even stable solutions have, most often, finite-size basins of attraction and do not necessarily act as global attractors, especially in the presence of noise, fluctuations,... or anything that can provoke dynamical transitions. The goal of the present chapter is to introduce the basic concepts necessary for the understanding of the next chapters, together with some historical/bibliographical remarks.

## 1.2 Governing equations

The fluid system of interest in this thesis is, almost without exception, governed by the incompressible Navier-Stokes equations, in their unsteady form. Before any non-dimensionalisation this system of partial differential equations reads

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{f}(t) \quad (1.1) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (1.2) $$

The space is by default assumed three-dimensional. In Eq. 1.1, $\rho$ is the fluid density (usually taken as unity) and $\nu$ the kinematic viscosity of the fluid, which we will assume to be known, constant and homogeneous. The scalar field $p$ denotes the pressure. It can in certain cases incorporate additional potential energies from which additional forces derive, for instance gravity or centrifugal forces. The external force $\mathbf{f}(t)$ does not need to be specified at this point. Most flows with the same Newtonian rheology obey the set of equations above. What distinguishes one flow case from another one is, rather than the governing equations, the domain of study, called $V$, and the boundary conditions that are associated with its boundary $\partial V$. Close to a solid boundary, consistently with the viscous hypothesis, we will usually consider that the fluid obeys no slip. As a consequence, the fluid velocity close to the wall equals the wall velocity:

$$ \mathbf{u}(\mathbf{x}_{wall}) = U_{wall}. \quad (1.3) $$

The other simple boundary conditions that are sometimes considered are free slip (used for mathematical simplicity), stress-free (at free interfaces) or other types. In particular, for "open" flows, in contrast to "closed" flows, $\partial V$ does not consist only of wall-like boundaries where the velocity is prescribed: the problem needs to consist also of entrance and exit conditions. We are left with the impossible task of modelling the state of the flow as it enters $V$, whereas we assume that the flow is not known outside $V$. Several strategies can be adopted: the entrance flow can be modelled, for instance by noise, by incomming flow stemming from another problem... the main hypothesis in this thesis consists of choosing to make the flow periodic, i.e. assuming that outflow is recycled into inflow. This hypothesis corresponds to the imposition of a given fundamental wavelength to the solutions of the system. This approach, highly debatable from a physical point of view, ensures in general that the mathematical problem is well-posed.

## 1.3 Canonical base flows

### 1.3.1 Planar flows

The most trivial solutions to the usual planar wall-bounded flows form a single parametric family of solutions to the incompressible Navier-Stokes equations. For a planar infinite geometry along the variables
x and z, let y stand for the direction normal to the solid plate(s). If a 2D steady solution is sought for, independently of the third variable z, it is written a priori $U = (U_x(x, y), U_y(x, y))$. Incompressibility then leads to $\partial_x U_x + \partial_y U_y = 0$. Searching for an x-independent flow leads to $\partial_y U_y = 0$ i.e. $U_y$ independent of y. If $U_y$ vanishes at the wall it is identically zero everywhere. Yet that does not necessarily need to be the case, for instance when uniform suction of blowing is applied at the wall, then $U_y = v$ at the wall (with v a constant). In all cases the nonlinear term $(U \cdot \nabla)U$ vanishes altogether too, and the streamwise momentum equation reduces to

$$ \nu \frac{\partial^2 U_x}{\partial y^2} = \underbrace{v \frac{\partial U_x}{\partial y}}_{\text{vertical transport}} + \underbrace{\frac{1}{\rho} \frac{\partial P}{\partial x}}_{\text{pressure term}}. \qquad (1.4) $$

Let us label the main terms $A = U_y \partial_y U_x$ and $B = \rho^{-1} \partial_x P$. The classification of the solutions depends on the values (zero or non-zero) of A and B :

- $A = B = 0 \rightarrow$ plane Couette flow (pCf),

- $A = 0$ and $B \neq 0 \rightarrow$ plane Poiseuille flow (pPf) or Couette-Poiseuille flow (CPf),

- $A \neq 0$ and $B = 0 \rightarrow$ Asymptotic Suction boundary layer flow (ASBL),

- $A \neq 0$ and $B \neq 0 \rightarrow$ generic Couette Poiseuille flow with suction.

We list below the usual velocity profiles for the canonical base flows that will be analysed here, together with the usual conventions for dimensional and non-dimensional variables.

### plane Couette flow

Plane Couette flow (pCf) is the simplest example of shear flow, forced by the difference in velocity between the two walls. In a frame travelling at the average velocity, the boundary conditions can be expressed as $U_x(y = \pm h) = \pm U_w$, together with $U_i(y = \pm h) = 0$ for $i = y, z$. The resulting laminar profile is hence given by

$$ U_x(y) = \frac{U_w y}{h}. \qquad (1.5) $$

The shear is $U_w/h$ and is spatially homogeneous. One frequent convention is to non-dimensionalise distances by the half-gap $h$ between the walls and velocities by the wall velocity $U_w > 0$. The resulting velocity profile is simply $U_x = y$. The Reynolds number is hence $Re = U_w h/\nu$.

### plane Poiseuille flow

Plane Poiseuille flow is the flow between two parallel plates, forced by a given axial pressure gradient. Meanwhile, in the frame of the walls, the flow at the walls is at rest, hence the boundary conditions reduce to $U_i(y = \pm h) = 0$ for alis $i$. For simplicity we denote $G = -\partial_x P$ the pressure gradient imposed in the x direction. The resulting laminar velocity profile is given by

$$ U_x(y) = G\left(1 - \left(\frac{y}{h}\right)^2\right). \qquad (1.6) $$

Note that $U_x \ge 0$ when $G \ge 0$, i.e. when the pressure $P$ decreases with $x$. The shear now depends on $y$, it is maximal in magnitude at the walls and it vanishes at the centerplane. Many different non-dimensionalisations coexist for plane Poiseuille flow, depending on the culture of the author and the exact forcing protocol selected. For distances, as for pCf the half-gap $h$ is the usual reference length. As for velocities, the most classical is perhaps non-dimensionalisation by the centerline velocity of the laminar regime $U_{cl} = G$. It is related to the bulk velocity $U_b = Q/2h$, where $Q$ is the flow rate by unit width.
The relation for a laminar plane channel flow is simply $U_b = 2U_{cl}/3$. At equilibrium the wall shear rate $\tau_w = \mu\partial u_x/\partial y$ evaluated and averaged at the walls is constant and is given by $\tau_w/h = G$. It is common practice in pressure-driven flows to use the Reynolds number $Re_\tau$ based on the friction velocity $u_\tau$ defined by $\rho u_\tau^2 = \tau_w$. This results in $Re_\tau = u_\tau h/\nu \sim \sqrt{\tau_w}$ so that $Re_\tau^2 = 3Re_b$ in the laminar regime.

## Asymptotic Suction boundary layer flow

The less known Asymptotic Suction boundary layer (ASBL) profile occurs in the presence of only one wall, which by convention is located at $y = 0$. The flow occurs only in the half-space $y > 0$ while the other half-space $y \le 0$ is ideally a dense porous medium. The system assumes Couette-like boundary conditions for $U_x : U_x(y = 0) = 0$, $U_x(y \to \infty) = U_\infty$. Simultaneously there is constant suction at the wall and $U_y$ is subject to $U_y(y = 0) = -V_S$ at the wall. The resulting base flow profile has two non-zero components rather than one, and is given by

$$U_x(y) = U_\infty \left(1 - \exp\left(-\frac{y}{\delta^*}\right)\right), \qquad (1.7)$$

$$U_y = -V_S. \qquad (1.8)$$

The convention is to non-dimensionalise velocities by $U_\infty$ and distances by the boundary layer thickness $\delta^*$, given analytically by $\delta^* = \nu/V_S$. This yields a definition for the Reynolds number $Re = U_\infty\delta^*/\nu = \frac{U_\infty}{V_S}$. The shear is maximal at the wall and the wall-normal velocity is constant and non-zero everywhere. In any regime this flow verifies the property that the mean friction velocity $u_\tau$ is fixed by the Reynolds number by $u_\tau^2 = U_\infty^2/Re$ [1]. The Reynolds number $Re_\tau$ based on the friction velocity is hence also constant on average when $Re$ is imposed.

### 1.3.2 Cylindrical flows

The same problem can also be solved inside a cylindrical geometry of radius $r_0$ when the mean flow direction is aligned with the cylinder axis (called by convention $z$). When the flow is forced by an axial pressure gradient, the laminar profile is parabolic and corresponds to the **Hagen-Poiseuille solution (HPf)**

$$U_x(r) = U_z(r) = \frac{G}{4\nu} (1 - (r/r_0)^2). \qquad (1.9)$$

A common modification of this geometry consists in adding a coaxial cylindrical wall inside the primary pipe at radius $r = r_i$ whereas the outer wall still lies at $r = r_o$. The flow in the gap between the two cylinders resulting from an applied axial pressure gradient $G$ is referred to as annular Poiseuille flow (aPf). Finally, another canonical annular flow can be suggested where the fluid is entrained, not by an axial pressure gradient, but by the axial translation of the inner rod; it is labelled **annular Couette flow**. If $\eta = r_i/r_o$, $d = r_0 - r_i$, $y = r - r_i$ and $U(r = r_i) = U_w$, it admits the simple analytical logarithmic solution

$$U_x(r) = U_w \frac{\ln\left(\frac{y(1-\eta)+\eta}{d}\right)}{\ln(\eta)}. \qquad (1.10)$$

### 1.3.3 Notion of forcing protocol

In practice, depending on the type of boundary conditions applied, two main classes emerge from the list of simple parallel shear flows above. In cases like plane Couette flow, the fluid flow is *entrained*, whereas for plane Poiseuille flow or pipe flow it is *pressure-driven*. There are several ways to achieve a given pressure-driven flow in practice, depending on which physical quantity is strictly imposed and which one is left free to adapt. The simplest conceptual way to achieve such a flow is to impose a fixed pressure gradient, this is the case for example when the fluid inside a straight vertical pipe is subject to a constant gravity field. Another less direct way to achieve a similar flow is to impose a given flow rate,
for instance by considering that the fluid is pushed or pulled by a piston or a syringe, whose own velocity
is known. Other protocols can be suggested, for instance imposing a fixed energy input. It is not obvious
which protocol characterises a priori a given hydrodynamic device, for instance standard pump engines
are rather believed to operate in the constant power input protocol [2]. Even entrained flow cases are sub-
ject to such a distinction, if one considers that the shear is imposed at the wall, rather than the velocity.
This is the case for instance in standard models of gas-liquid interfaces although this has recently been
questioned in liquids [3]. In the strictly steady regime laminar regime there is no mathematical distinction
between these protocols. This is not true in other circumstances. For instance in the steady turbulent
regime, the mean values of the kinetic energy, dissipation, etc... are independent of the protocol, however
the fluctuations around these mean values are not.

The standard convention in wall turbulence studies (but not always in transition studies!) is to use
different governing parameters depending on the forcing protocol. For flows driven by an imposed flow
rate, the Reynolds number (which necessitates a velocity scale) is usually based on the so-called bulk
velocity $U_b$ (defined classically by the flow rate divided by the equivalent cross-sectional area). The ad-
vantage of using $U_b$ in that case is that $U_b$ stays the same no matter whether the flow is laminar or
turbulent, whereas the value of the corresponding mean pressure gradient is strongly dependent on the
regime. For flows driven by an imposed pressure gradient $G$, one can exploit the equivalence between the
pressure gradient $G$ and the mean wall shear stress $\tau_w$. Such a relation only involves average values and is
hence not valid instantaneously, unlike its fixed flow rate counterpart. A velocity scale $u_\tau$ can be defined
from the knowledge of $\tau_w$ via the relation $\tau_w = \rho u_\tau^2 \sim G$. This allows one to define the friction Reynolds
number $Re_\tau$ based on the velocity scale $u_\tau$ rather than $U_b$. By construction $Re_\tau$ is proportional to $\sqrt{G}$
independently of whether the flow is in a laminar or turbulent regime.

1.4 Wall turbulence at onset

Writing a whole thesis about transition to turbulence without ever documenting the turbulent state
itself would look incomplete. This section does not pretend to contain all the knowledge on wall turbulence.
However it aims at presenting briefly the minimum knowledge necessary for a decent theory of transition,
including a justification of the present parameter range of investigation.

1.4.1 High-Reynolds number wall turbulence

Anyone who has glanced at the beginning of a textbook about turbulence might have faced the same
question: although the meaning of the word is not so ambiguous in everyday life, is there a solid, uni-
versally accepted definition of what *turbulent* means in hydrodynamics? Although laminar flows are
perfectly well defined, turbulent flows are not, although we expect them to have all the worst possible
complications: finite unsteadiness, correlations in time and space, multiscale structure, impredictability...
Etymology hardly helps here, at the very least it indicates (from "turbinis" in Latin) that the flow has to
contain... vortices, which does not really come as a surprise. Everyday intuition might be more useful here.
We all have a mental image of what a *very turbulent* flow looks like from watching rivers, cascades, weather
maps. Hence it sounds wiser to consider first these conditions which make fluid flows *very turbulent* and to
test how far they can be pushed when that same Reynolds number is reduced. In that spirit we should look
at published works on huge values of the Reynolds number $Re$. We emphasize that, no matter how large,
$Re$ should nevertheless be finite since it is well known that the infinite Reynolds number/zero viscosity
limit is singular and might lead to very wrong speculations. This looks preferable to an abstract theory
from scratch à la Bourbaki, because it is more likely to reconnect properly with the real engineering world.

Let us have a look at a set of recent well-resolved turbulent channel flow simulations. Such simulations
are considered in periodic domains and run over long enough times to gather converged flow statistics.
Figure 1.1 displays a rendition of vortex clusters for a high Reynolds number of $Re_\tau \approx 4200$ in a domain
of size $2\pi h \times 2h \times \pi h$. It is confirmed that vortices abound in the flow, but also that many scales
<figure><img src="image_2.png" /><figcaption>FIGURE 1.1 – Three-dimensional visualisation of vortex clusters in a simulation of channel at $Re_τ \approx 4200$, from Ref. [4] (only the bottom of the channel is shown).</figcaption></figure>

coexist on a single snapshot, and that clusters of small scales form larger scales. Figure 1.2, based on a simulation of the same flow at a slightly lower value of $Re_τ$, displays the turbulent kinetic energy of the streamwise velocity component plotted in a $(\lambda_z^+, y^+)$ diagram. $y^+$ represents the distance from the lower wall expressed in units of $\delta_\nu = h/Re_τ$, the classical units of wall turbulence, whereas $\lambda_z^+$ represents a wavelength in the spanwise ($z$) direction, also normalised by $\delta_\nu$. Such a figure is less sexy yet more useful than visualisations. It appears clearly that a whole set of scales ranging from $\lambda_z^+ \approx 40$ up to $\lambda_z^+ \approx 10^4$ (and more for higher $Re_τ$) is present in the flow. Worse : the distribution of scales involved energetically depends strongly on the distance from the wall, with energetic features present from as low as $y^+ \approx 10$. The whole distribution of energy apparently lies around a line $y \sim 0.1\lambda_z$. This echoes the *Attached Eddy Hypothesis* of A. Townsend [5, 6] sketched in figure 1.3. This theory states that the size of each eddy is proportional to its distance from the wall, as opposed to dynamically inert *detached eddies*. In this sketch the wall-parallel space variable can be either $x$ or $z$. Although the multiscale nature of turbulence is also not surprising, what matters here are the numbers. The smallest eddy in this hierarchy, equivalent to the Kolmogorov scale in isotropic settings, turns out to be slightly less than $100\delta_\nu$ no matter how high $Re_τ$. This validates the importance of inner units $O(\delta_\nu)$ for scaling laws. The largest eddies in the hierarchy are simply bounded by the channel geometry, which can be estimated as the channel height $2h$. The multiscale nature of wall turbulence can be attributed grossly to the fact that these two orders of magnitude, $100\delta_\nu$ on one hand and $2h$ on the other hand, *do not coincide*. It is precisely in this scale interval that universal features, such as the law of the wall $\bar{U} \approx \frac{1}{\kappa}\log(y^+) + B$, are valid [8]. This whole picture holds of course as long as $100\delta_\nu \ll 2h$, which supposes $Re_τ$ large enough.

### 1.4.2 Estimating the onset of wall turbulence

The previous estimates imply a simple yet powerful reverse idea : multiscale turbulence cannot sustain if $100\delta_\nu > 2h$. Besides, if $100\delta_\nu = 2h$, the resulting turbulent regime is characterised by one lengthscale only. It does feature any log law or $k^{-5/3}$ energy spectrum typical of turbulent cascade, and the turbulent flow is limited to the buffer layer only, together with a viscous sublayer for $y^+ < 5$. Given the definition of $\delta_\nu = h/Re_τ$, this suggests as a gross estimation that the turbulent state should collapse, should $\delta_\nu = h/Re_τ < 50$ be verified. Such an estimate is as universal as the scaling laws and is independent of the flow case. It has been advanced or re-discovered several times in the literature [9, 10] and proves correct as a rule of thumb.
<figure><img src="image_3.png" /><figcaption>FIGURE 1.2 – Streamwise Energy as a function of $(\lambda_z^+, y^+)$ for a DNS of channel flow at $Re_\tau = 2003$ (from Ref. [5]).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 1.3 – Sketch of the Attached Eddy Hypothesis (from Ref. [5]).</figcaption></figure>
What this estimation does not reveal is what such a collapse should look like if one starts from a pre-
existing turbulent velocity field. There is ample evidence that in almost all wall-contained shear flows,
an adiabatic reduction of the Reynolds number yields a spatially and temporally intermittent state of
turbulence : where it is turbulent, the flow field looks similar to developed turbulence, but it is not tur-
bulent at every location at the same time.

The way turbulence collapses is illustrated in figure 1.4, taken from Ref. [11] in pipe flow driven by a fixed mass flux. The Reynolds number *Re* is progressively reduced in small steps starting from a high enough value at which a turbulent flow field is present, no matter how it got generated. The quantity visualised here is the streamwise velocity fluctuation A similar diagram can be found in Ref. [12] for other velocity components and is reproducible, provided the boundary conditions are periodic in the axial direction and the length of the computational domain is large enough. The Reynolds number is here based on the bulk velocity and the pipe diameter. For *Re* = 2800 and above the flow looks unambiguously turbulent, moreover it appears uniformly turbulent, in the mild sense that the large-scale turbulent properties do not depend on the axial coordinate (note that "homogeneous turbulence" is a different, more demanding concept in turbulence theory). For *Re* = 2600 the flow looks similar except for a "hole" located around *x* = 150 where the turbulent activity is less intense, at least at the moment displayed here. More such holes are present at *Re* = 2400. At *Re* = 2200, something has changed qualitatively : rather than laminar-like holes inside a turbulent sea, it is better described by a sequence of localised turbulent-like structures in an otherwise laminar-like environment. These coherent structures look all alike, they are called *puffs* and they form an apparent pattern with defects inside. Whether or not an exact wavelength emerges from a well-defined relation dispersion, like in Turing instabilities, is today still a matter of debate. As *Re* decreases the number of puffs goes down until extinction. At *Re* = 1680 a last puff survives and it is known from the literature that this puff has a finite lifetime. Once the last puff has died the flow is said to have *relaminarised*. There is a complex dynamics associated with puffs and trains of puffs [13] that cannot be grasped easily from still pictures, but will be explored in Chapter 4.

The collapse of turbulent flows in planar geometries obeys a similar story. It is illustrated in figure 1.5 in the case of plane Poiseuille flow driven by a fixed pressure gradient. The advantage of this kind of forcing over the former fixed mass flow protocol, is precisely that here $Re_τ$ is the control parameter, whereas the traditional $Re$ should be regarded as an (unsteady) output of the simulation. This represents an opportunity to test the above $100\delta_ν = 2h$ criterion. The quantity visualised in the diagram is the cross-flow energy $E_{cf}$ based on the transverse velocity components only, and the spatial variable is the streamwise coordinate $x$ in a frame moving with the streamwise bulk velocity $U$, and $Re_τ$ is decreased in steps of 10. As in pipe flow, the flow appears homogeneously turbulent for $Re_τ = 100$ and above. Large streamwise modulations appear around $Re_τ = 90$. These modulations travel downstream at a velocity slightly slower than the bulk velocity. They grow in amplitude and turn into increasingly disordered laminar-turbulent patterns present from $Re_τ = 90$ down to approximately 55. Below $Re_τ = 50$, a radical change appears : the individual structures forming the already irregular pattern accelerate, they dissociate, and the resulting dynamics resembles ballistic trajectories of uncorrelated localised structures. Lowering further down $Re_τ$ below 40 leads to the individual extinction of each of these structures, eventually leading to a perfectly laminar flow.

The transitional range associated with pressure-driven plane Poiseuille flow occurs grossly, according to the above visualisation, between $Re_τ = 40$ and 100. This is highly consistent quantitatively with the $100\delta_ν = 2h$ that suggests a critical value of 50. It remains to understand why there is an extended transitional range, rather than a simple critical point with laminar flow below and turbulent flow above.

1.4.3 Relevant interrogations regarding transition to turbulence in shear flows

The Moody diagram illustrates how the friction factor $C_f$ of a given flow evolves as a function of the Reynolds number $Re$ (traditionally understood as the bulk-based Reynolds number $Re_b$). $C_f$ is often
<figure><img src="image_5.png" /><figcaption>FIGURE 1.4 – Simulations of circular pipe flow for decreasing *Re* from [11].</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 1.5 – Simulations of plane Poiseuille flow for decreasing $Re_τ$ from [14].</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 1.6 – Moody diagram for pipe flow (source : www.nativdynamics.com.au)</figcaption></figure>

introduced in textbooks either as an efficiency term as compared to the inviscid equivalent, or preferably as the head loss along the considered portion of the flow, normalised by a kinetic energy term based on the bulk velocity. Although convenient for experimentalists, this last definition should still be considered when care when dealing with numerics, because it requires proper averaging (see Ref. [14] for a discussion). We prefer here the definition

$$C_f = 2 (Re_{\tau}/Re_b)^2, \qquad (1.11)$$

which has to be understood differently according to the chosen protocol. For pressure-driven flow $Re_τ$ is imposed, but $Re_b$ is not and is understood as a time average of the unsteady $Re_b$. Conversely for mass-driven flows, $Re_b$ is constant but $Re_τ$ should be defined based on the time-and-space averaged wall-shear stress. For laminar flow, both protocols are equivalent and $C_f$ can be computed analytically. For all laminar shear flows we find that $C_f \sim 1/Re$ with a coefficient specific to the flow geometry and the non-dimensionalisation only. The Moody diagram for pipe flow is depicted in figure 1.6. For high Reynolds numbers in industrial settings, $C_f$ is a constant dependent only on the type of roughness present at the walls. It is this property that justifies the constant coefficients found in all aerodynamics textbook, e.g. the $C_x$ of a car reported in car catalogues. When walls are not rough (which implies either laboratory experiments or low Reynolds number), $C_f$ depends on $Re$. Again three regimes emerge. The two limiting regimes are the laminar regime, clearly distinct at low $Re$ with the expected $O(1/Re)$ scaling, and a turbulent regime characterised by fits $O(Re^{-1/4})$. The hacked area is traditionally left with a question mark in textbooks and labelled either 'complicated', 'critical', 'uncertain' or an expression hardly clearer : 'transitional'. In the transitional range (according to visualisations) turbulence is spatially and
temporally intermittent. It is the region where mean values do not suffice to shed light on the dynamics of the flow and where the study of fluctuations becomes paramount.

Two complementary approaches will be used in this work depending on whether the focus is on the "transitional" regimes themselves, which exist in a limited range of Reynolds numbers, or on the transition towards the turbulent state, which is possible as soon as a turbulent state can be defined.

Starting from the turbulent state, we can define the "adiabatic descent procedure", inspired by the experimental methodology used e.g. by Prigent for plane Couette flow, as in figures 1.4 and 1.5. A turbulent state is achieved at a sufficient high value of $Re$, then $Re$ is decreased in very small steps (or equivalently as a ramp with $dRe/dt < 0$). By analogy with thermodynamic processes this is called an adiabatic decrease (in $Re$). It allows one to track a given turbulent state in parameter space. It is very robust when the turbulent state is unique, which is the usual case here. This process delivers the full zoology of states present in the transitional range, as well as above it if the initial value of $Re$ is high, and works well in unbounded geometries. It also yields in principle the true critical Reynolds numbers below which no turbulence is found, or at least the critical Reynolds number associated with the finite observation time chosen. Figures 1.4 and 1.5 also highlight the huge lengthscales involved in the collapse process. Such simulations are very CPU-demanding and one can be tempted (or forced) to shorten the domain at equivalent local resolution. Simulations in numerical domains smaller than such lengthscales risk either deforming the associated phenomena or simply missing them.

We can adopt a point of view starting from the well-defined laminar state instead, the "path to turbulence" procedure. There is ample mathematical, experimental and numerical evidence that, for the flows of interest here, the laminar state is always stable with respect to infinitesimal disturbances, no matter what the shape of such disturbance is. However as soon as a turbulent state is realisable for a given Reynolds number, this suggests that certain disturbances of finite amplitude are able to provoke a change of regime. Certain preferred paths towards the turbulent state exist in a reproducible way. Out of mathematical convenience, stability theory is often formulated as an initial value problem, consistently with the evolution form ($\partial_t = ...$) of the governing equations. Other approaches can be envisioned, such as the receptivity to noise [15], or to sinusoidal excitations (e.g. resolvent analysis [16]), depending mainly on their potential to match a given experimental (or numerical) protocol. They are not treated here. A large part of the work reported here deals instead with transition to turbulence as an initial value problem, based either on the full set of nonlinear equations, on their finite-dimensional truncation in numerical studies, or sometimes on a low-dimensional system of model equations when pedagogically useful. The emphasis is on the determination of physically meaningful initial conditions and the associated trajectories in a suitable state space. For such a task the theory of dynamical systems, an intrinsically nonlinear approach well adapted to complex dynamics, is of direct relevance. We will explore the current state of the art in the quest for a description of the transition process in terms of its state space description. Its limitations will also be pointed out with possible self-criticism.

## 1.5 Linear stability of canonical base flows

Symbolically the nonlinear equation system 1.1-1.2 can be rewritten in the general form :

$$ \frac{d\mathcal{U}}{dt} = \mathcal{L}\mathcal{U} + \mathcal{N}(\mathcal{U}) \qquad (1.12) $$

where $\mathcal{L}$ is a linear operator and $\mathcal{N}$ contains all nonlinear terms. The meaning of $\mathcal{U}$ does not really need to be made explicit at this stage. It is enough to know that it contains all the independent degrees of freedom of the velocity field. The pressure field $p$ can be determined uniquely (modulo an additive constant) by considering the divergence of Eq. (1.1) and does not need to be considered as an additional field.
We label here "base flow" the steady solution $\mathcal{U}_0$ of Eq. 1.1 consistent with the boundary conditions. As far as trivial solutions with a high degree of symmetry are concerned, they are usually unique. As a steady solution it must verify $\frac{d\mathcal{U}_0}{dt} = 0$. The equation for the perturbation $\mathbf{X} = \mathcal{U} - \mathcal{U}_0$ is

$$ \frac{d\mathbf{X}}{dt} = \mathbf{L}\mathbf{X} + \mathcal{N}(\mathbf{U}) - \mathcal{N}(\mathcal{U}_0). \qquad (1.13) $$

By assuming that the perturbation $\mathbf{X}$ is small enough in a given norm, the above system can be linearised in the neighbourhood of $\mathcal{U}_0$. The linearised system reads in symbolic form

$$ \frac{d\mathbf{X}}{dt} = \mathbf{L}\mathbf{X}. \qquad (1.14) $$

The issue of the stability of the base flow $\mathcal{U}_0$ with respect to infinitesimal flow disturbances is dictated by a few key properties of the linear operator $\mathbf{L}$. Writing down $\mathbf{L}$ explicitly is a tedious task. In practice, the use of a numerical discretization often leads to a large square matrix (also referred to as $\mathbf{L}$ for simplicity) whose properties can be determined computationally.

<figure><img src="image_8.png" /><figcaption>FIGURE 1.7 – Definition of different critical values of the Reynolds number $Re_E$, $Re_L$ and $Re_c$. Sketch of typical time series of the kinetic energy $\varepsilon$ for several Reynolds numbers, after Ref. [17].</figcaption></figure>

A set of three critical values of the Reynolds number must be considered, depending on the exact notion of stability considered [17]. They are defined graphically in figure 1.7.

The simplest critical value, which also corresponds to the lowest one in the hierarchy, is the energy Reynolds number $Re_E$. It is defined as the highest value of $Re$ below which all perturbations to the base flow decay monotonically to zero. By all perturbations we mean all shapes as well as all amplitudes. This is hence a nonlinear concept. The Reynolds-Orr theorem leads however to a simplified defining criterion for $Re_E$ invoking only the linear operator $\mathbf{L}$.

The usual critical value considered in linear stability theory is the critical Reynolds number $Re_L$, defined as the smallest value of $Re$ below which at least one class of perturbations grows indefinitely under the action of the linear operator. This is a linear concept only. It corresponds to the smallest value of $Re$ for which the spectrum of $\mathbf{L}$ contains at least one eigenvalue with positive real part. In shear flows, unlike e.g. Rayleigh-Bénard convection or classical Taylor-Couette flow, $Re_E$ never coincides with $Re_L$. $Re_L$ does
not necessarily match the value at which transition to turbulence is observed in experiments and numerics.

The global critical Reynolds $Re_c$, corresponds to the largest value of $Re$ below which all perturbations to the base flow eventually return to it. It is a fully nonlinear concept unlike $Re_L$ and $Re_E$. In other words it is the critical value above which nonlinearity can maintain turbulence provided the right initial condition(s) are supplied. Even if a given perturbation induces turbulence, if the turbulent state relaminarises in finite time with probability 1, there is no guarantee that $Re_c$ has been reached. $Re_c$ is, in the adiabatic limit, the value yielded by the adiabatic descent procedure.

### 1.5.1 Transient growth

It is crucial to distinguish between finite-time stability and infinite-time stability [18]. The long-time stability is fully determined by the spectrum of $L$, denoted as $Sp(L) = \{\lambda \in C, \exists X \neq 0, LX = \lambda X\}$. In finite dimensions the spectrum of $L$ is discrete and corresponds to the eigenvalues of the associated matrix. It is a classical result that if all eigenvalues in the spectrum have strictly negative real parts, then the system is asymptotically stable for long times. The decay of perturbations is also dominated on long times by the least stable eigenvalue (the one with the largest real part). However, transient growth of $||X|| = \sqrt{X} \cdot X$ is also possible even if all eigenvalues of $L$ indicate long-time stability. This arises provided $L$ and its adjoint (defined using the same scalar product as for the 2-norm) do not commute, i.e. $LL^* \neq L^*L$, in which case $L$ is called non-normal. Another equivalent measure of non-normality, for invertible operators, is the condition number of the eigenbasis, defined as $\kappa = ||L|| \cdot ||L^{-1}||$, where the matrix norm used here is the usual Frobenius norm. $\kappa = 1$ for normal operators while $\kappa > 1$ for non-normal operators [19]. A symmetric or anti-symmetric operator is for instance normal. A non-normal operator can be diagonalisable but this is not a necessary condition. If it is diagonalisable, then the eigenvectors do not form an orthogonal basis.

A numerically determined point spectrum typical of shear flows, belonging to plane Poiseuille flow at $Re = 10^4$, is displayed in figure 1.8 using dots. Any dot lying above the real axis represents an unstable mode. In the figure almost all dots lie below that line, indicating that damped modes dominate. Only one dot lies above it, and it can be established that no such unstable mode is found when $Re < 5772$ [20]. The pseudo-spectrum $\Lambda_{\epsilon} = z \in C, ||(zI - L)^{-1}|| \ge \epsilon^{-1}$ is a useful generalisation of the spectrum [21].

A popular interrogation emerging from the non-normality property is the following: given the linearised context, what is the initial condition producing the largest transient growth of perturbation kinetic energy? Let us suppose that the perturbation kinetic energy $E(t)$ is given by $\frac{1}{2}||X(t)||^2$, that the related energy gain at time $t$ is defined as $G(t) = E(t)/E(0)$, and that the optimal gain is given by

$$G_{max} = \max_{X,t} G(t). \qquad (1.15)$$

For operators that are not asymptotically stable, there is at least one eigenvalue with positive real part: if a given initial perturbation has a non-zero projection in the related eigenspace, then this perturbation will grow indefinitely with time in that direction; in that case the optimal gain is $G_{max} = \infty$. For stable normal operators, $G_{max} = 1$ since the decay of the perturbation kinetic energy is monotonic. A non-trivial situation occurs when $L$ is stable but non-normal. In such a situation there is potential for energy growth over finite times.

It was determined numerically that $G_{max}$ rises rapidly with the Reynolds number, in fact as $O(Re^2)$ [18]. This result holds apparently regardless of the flow case. It is this result, rather than the non-normal property itself, that justifies why laminar shear flows are so sensitive in practice to arbitrary perturbations such as incoming noise or small geometrical imperfections.

The so-called *optimal modes* are the flow structures achieving maximum growth in finite time. When considered in three dimensions, they correspond invariably to structured arrays of streamwise vortices. If
<figure><img src="image_9.png" /><figcaption>FIGURE 1.8 – Spectrum (dots), pseudo-spectrum for $\varepsilon = 10^{-n}$, $n = 1$ to 8 (solid lines) and numerical range (dashed line) for plane Poiseuille flow at $Re = 10^4$, from Reddy *et al.* (1993) [19].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 1.9 – Levels of $G_{max}$ in plane Poiseuille flow depending on the axial wavenumber $\alpha$, and on the Reynolds number $R$ based on laminar the centerline velocity, from Reddy *et al.* (1993) [19]. Dashed line : $G = 1$, dark area : $G = +\infty$, other lines correspond to $G = 5, 10, 20, 30, 40, 50, 60$.</figcaption></figure>
such a flow field is considered as an initial condition in the framework of an initial-value problem, then the output at the optimal time $t_{max}$ is also a structured array of streamwise streaks superimposed on the original array of streamwise vortices. The streaks are defined as spanwise variations of the streamwise velocity field, i.e. an alternation of slow and fast lanes forming a given spanwise wavelength. These opti-

<figure><img src="image_11.png" /><figcaption>FIGURE 1.10 – Sketch of the lift-up effect from F. Charru's textbook [22] showing the development in time or space of streamwise streaks from streamwise vortices.</figcaption></figure>

mal streaks, as well as the associated vortices are independent of the streamwise variable. The amplitude of the vortices themselves decreases steadily with time whereas the amplitude of the streaks reaches its maximum at $t = t_{max}$. The development of streaks from an initial distribution of streamwise vorticity is labelled *lift-up* to celebrate the fact that streamwise vorticity is able to *lift up* the streamwise velocity isocontours, although by continuity they are associated with *push-down* regions where the fluid is brought down closer to the wall by the neighbouring vortices. If the wall is at rest, upwards motion inside the vortices generates a low-speed streak, whose instability will be emphasized later. This lift-up effect, in the streamwise-independent case, is formally and mathematically similar to the rise of plumes by vortices in the advection of a passive scalar. The difference of the shear flow case with the two-dimensional convection picture is that a third velocity component arises. This allows one to generate a three-dimensional vorticity field out of an initial condition that features only streamwise vorticity. Another way to interpret the lift-up mechanism is that it demonstrates the advection of the base flow by the flow perturbation [23].

The importance of streamwise streaks in transitional as well as turbulent shear flows should not be minimised : they represent the main coherent structures found near the wall in all flows featuring boundary layers. This includes *all* the flow cases investigated here. A velocity field of plane Couette flow is computed using a visualisation routine developed by J. F. Gibson [24] in figure 1.11. The red/yellow lanes correspond to the highspeed streaks whereas the blue lanes, which appear slightly more sinuous, are the low-speed streaks. Streaks can also be found in free jets or wakes. One of the rare "formulas" of the shear flow literature is the mean spanwise wavelength of the streaks, found almost universally to be $\lambda_z^+ \approx 100$. This formula alone justifies the introduction of inner units (the non-dimensionalisation by $\delta_\nu = \sqrt{\tau_{w\nu}}$ instead of $h$). This simple scaling relation can be extracted directly, in a statistical sense, from the small-scale peak present in figure 1.2. The linear framework, more elegantly reformulated as an input-output formalism [25], is convenient to explain why streaks emerge, not only when initialised by the optimal mode, but also from a disordered vorticity field and in finite time.

The linear point of view on subcritical transition has certainly put an end to a mathematical mystery. Beyond explaining the universal emergence of streaks in high Reynolds numbers shear flows, it has also
<figure><img src="image_12.png" /><figcaption>FIGURE 1.11 – Velocity field of plane Couette flow from a computation in Ref. [24]. The red/yellow lanes correspond to the highspeed streaks whereas the blue lanes, which appear slightly more sinuous, are the low-speed streaks.</figcaption></figure>

emphasized the importance of non-normality in generic linear dynamics. This is perhaps welcome after decades of influence from quantum physics, where the emphasis is on symmetric, hence normal operators. However turbulence is a nonlinear phenomenon. A linear theory alone, although it can yield a path towards a turbulent flow, cannot explain the emergence and, more importantly, the ability of a turbulent flow state to sustain itself.

1.5.2 Nonlinear extension

For a linear system of equations, all solutions found independently can be combined linearly to form new solutions. In planar flows, the equivariance of the equations with respect to the two planar directions justifies ansätze of the form $e^{ik_x x+ik_z z}$ when solving the linearised equations. Although general solutions contain a summation of these exponential solutions, linearity allows to focus on each wavenumber separately. Such an assumption leads to the boundary conditions $u'(x + 2\pi/k_x, z + 2\pi/k_z) = u'(x, z)$ being automatically satisfied, by construction. Historically, after linear concepts have been applied to the study of the stability of fluid flows, it has become very tempting to generalise the system by restoring the neglected nonlinear terms... even if in the nonlinear context, the study of individual exponential ansätze is no longer justified! Although it is obvious that no nonlinear open flow obeys a strict spatial periodicity condition, there are some cases for which this hypothesis can be justified. On one hand, if the wavelength is sufficiently large, for instance much larger than the correlation length of the system, the natural outflow and inflows are in practice uncorrelated : in such a situation the error committed by recycling outputs into inputs is less problematic. On the other hand, if that imposed wavelength is too short, additional unphysical correlation is added into the system of equations, and we deal with a flow case that cannot be reproduced experimentally. Another case in which spatial periodicity is relevant is the numerical simulation/analysis of localised velocity fields, at least when the spatial decay is sufficiently rapid.

Not only do periodic boundary conditions make the mathematical analysis of open flows mathematically better-posed, they also imply simpler theorems. For instance, the Reynolds-Orr equation governs the time evolution of the energy $e = \frac{1}{2} \langle u' \cdot u \rangle$ of the fluctuations. It is obtained by multiplying scalarly the equation for $\frac{dX}{dt}$ by $X$ and making ample use of integration-by-parts theorems. In the presence of periodic boundary conditions in the planar variables, it can be verified that $\langle N(X), X \rangle = 0$. The
Reynolds-Orr then reads in compact form :

$$ \frac{de}{dt} = \langle LX, X \rangle. \qquad (1.16) $$

This equation, also valid in the nonlinear regime, states that the energy amplification is independent of the exact nonlinear terms, provided the energy-conserving property $\langle N(X), X \rangle = 0$ is verified. By dividing Eq. 1.16 by $|u'|^2 = |X|^2$, this conclusion can be pushed even further : the instantaneous rate of change of the perturbation kinetic energy does not depend on the instantaneous perturbation energy given by $|u'|^2$. Although this has generated controversy [26], it is now clear that this does not make the nonlinear terms less important. The norm of $X$ might not depend on $N$, but its orientation definitely does, and the term involved in the momentum equations is $X$, not $||X||$.

# 1.6 Nonlinear concepts

## 1.6.1 Self-sustaining process

Within the previous linearised framework, the main and perhaps only question addressed deals with the ability of a perturbation to escape, or not, the neighbourhood of the base flow solution. If restoring nonlinearity into the original equations is a welcome departure from linear stability theory, it does not indicate how to deal mathematically with the new nonlinear system. The true change of paradigm occurs when, instead of assessing *how far from the laminar state* the system temporarily lies, one gets interested in *how close to new alternative flow states* the system can get. This new type of interrogation brings back the question of *turbulent state* into the picture. It is mathematically consistent with the possibility that the laminar base flow, when described as a fixed point in state space, has a basin of attraction that does *not* necessarily coincide with the full state space. Now provided other dynamical regimes are possible, there must be a mechanism that maintains their dynamics *away* from the laminar regime *for all times*. Such a mechanism, called simply *self-sustaining process* [27], is necessarily nonlinear. Its identification and validation are not trivial and constitute a mechanistic way of understanding the existence of a turbulent state.

<figure><img src="image_13.png" /><figcaption>FIGURE 1.12 – Left : Reduction in domain size, taken from Jiménez & Moin [28]. Right : instantaneous wall vorticity highlighting the wavy streak.</figcaption></figure>

Crucial progress in this direction has been made by Fabian Waleffe and co-authors in a series of papers in the second half of the 1990s [9]. Its strategy can be summarised by following their 1995 paper :

- consider a flow case (they chose plane Couette flow) with a linearly stable base flow

- use direct numerical solution to simulate the turbulent regime inside a spatially periodic box,

- decrease carefully both the Reynolds number *and* the box size while keeping the dynamics non-relaminarising as in Figure 1.12,
— when the system seems on the verge of relaminarising, analyse the dynamics over a finite-time segment

The success of this strategy is to yield a new type of flow with its own dynamics, that appears at the same time prone to analysis while being strongly reminiscent of the dynamics observed in wall turbulence simulations. In particular, the velocity field found in plane Couette flow after such a reduction, shown in figure 1.13, possesses long wavy streaks flanked by shorter streamwise vortices, as found in many turbulence simulations at moderate Reynolds numbers.

<figure><img src="image_14.png" /><figcaption>FIGURE 1.13 – From a) to h) : ordered snapshots of streamwise velocity field at mid-gap in pCf, investigated by [9], illustration of one turbulent cycle consistently with the proposed self-sustaining process.</figcaption></figure>

The streamwise streaks are not steady and they are known to undergo so-called bursting events also described as "streak breakdown". It is indeed evident in figure 1.13, where the long sinuous coherent structure in figure 1.13a gets increasingly distorted with time until it almost disappears in figure 1.13f. It then gets regenerated back to the same amplitude as in a), which closes the loop. It is perhaps best understood in terms of vorticity components : a streak is characterised in the $xz$ plane of figure 1.13a as a wave of wall-normal vorticity $\omega_y$. The streak instability starting from figure 1.13b amplifies $\omega_y$ [29], up to levels where the nonlinear behaviour of the instability starts to manifest itself : the streak breaks down. Through the advection by the mean shear, $\omega_y$ gets tilted in the streamwise direction and its streamwise component $\omega_x$ increases in amplitude. Eventually, the lift-up mechanism [18], essentially the advection of the base flow by streamwise vortices, regenerates streaks and restores the original level of $\omega_y$. For the case of a spatially correlated minimal box, this sequence of events illustrates the self-sustaining process of Waleffe as sketched in figure 1.14. Variations on this mechanism can be found in the literature provided a slightly different computational domain is used, however they mainly impact the symmetries of the streak instability without modifying the sequence of events and the approximate period of the process.
As will be seen in future chapters, slightly different three-dimensional SSPs can be identified in shear flows yet only in the presence of spatial localisation.

<figure><img src="image_15.png" /><figcaption>FIGURE 1.14 – Sketch of the self-sustaining process from F. Waleffe [30].</figcaption></figure>

Dynamically speaking, the self-sustaining process is found to be weakly chaotic although the authors once hoped they would reduce to a perfectly cyclic dynamics [31]. This intuitively suggests, by analogy with simple low-dimensional chaotic attractors such as the Lorenz or the Rössler attractor -where the attractor can be covered by unstable periodic solutions-, that an ordered cyclic process is embedded in this chaotic dynamics. It should feature exactly the same ingredients as the wilder self-sustaining process (wavy streaks and rolls). This was confirmed a few years later by G. Kawahara & S. Kida who identified numerically an unstable periodic orbit (UPO) apparently embedded in the turbulent dynamics [32]. This apparent embedding was suggested by both the sequence of events visualised in videos, and by a convincing match between the statistics of the UPO of mean flow, r.m.s velocity fluctuations etc... and those from turbulent simulations. A sequence of events taken along the orbit is displayed in figure 1.15 and favourably compares with the material displayed in figures 1.13 and 1.14 Another convincing feature, although only qualitative, is the fact that the UPO seems intricately linked with the state portraits displayed in Figure 1.17.

The hunting season for other invariant solutions mimicking the turbulent dynamics was then declared open. Waleffe himself completed his description of the state space of plane shear flows by identifying solutions simpler than the UPO of Kawahara and Kida [33]. The associated numerical strategy, not so easily generalisable, assumes a forcing term able to generate wavy streaks via the subcritical instability of straight streamwise streaks. It then makes use of a numerical continuation technique to identify a solution that sustains itself without forcing. This method was then used by several groups to yield brand new solutions [34, 35]. It was at the same time understood that formerly found finite-amplitude solutions, notably the so-called Nagata solutions of plane Couette flow [36] originally identified via numerical homotopy, were from the same large family and could be found using similar techniques.

The periodic cell considered first by Hamilton *et al.* [9] was re-considered by P. Cvitanović and co-workers (a by-product of this investigation is the now freely available DNS code *channelflow* developed by J.F. Gibson). Their investigations shed a clearer light on the complexity of the state space of even that simple flow. Multiple steady states [37], heteroclinic connections [38] as well as periodic orbits [39] were identified and plotted using savant projection systems. These authors pointed out the need for symmetry reduction in order to take into account systems with advection [40]. That led ultimately to representations of the state space of pipe flow [41, 42], as visualised in an arbitrary Poincaré section in Fig. 1.19.

Although experimental verification of such concepts might seem one step more arduous (but decisive
<figure><img src="image_16.png" /><figcaption>FIGURE 1.15 – Three-dimensional visualisations of the velocity field along the turbulent-like periodic orbit found by Kawahara & Kida [32].</figcaption></figure>

for the credibility of the theory), careful work worth mentioning has been achieved by a few teams. The original theoretical idea of Waleffe to identify self-sustaining solutions by introducing a forcing term, waiting for the flow to bifurcate and then bringing back the forcing to zero was tried successfully by S. Bottin and co-workers at CEA Saclay. By using wires and beads of increasingly small size, they could confirm that streamwise vortices are indeed crucial elements of possible finite-amplitude solutions of the Navier–Stokes equations [43]. Later, excited by the recent numerical discovery [34, 35] of unstable travelling wave solutions in cylindrical pipe flow, Hof et al. [44] used cutting edge time-resolved Particle Image Velocimetry to report on furtive observations of similar travelling wave solutions from at a bulk Reynolds number of $Re=2400$. Another such travelling wave, analogous to one identified numerically by Pringle & Kerswell [45] was identified in the wake of a decaying turbulent zone by the same experimental group [46] in pipe flow. The concept of self-sustaining process itself was claimed to be observed in the wake of vortex generators in a boundary layer set-up [47]. More recently, focusing on the so-called Kolmogorov-flow driven by a sinusoidal shear profile, a detailed comparison between exact solutions computed from DNS versus those computed from experiments was made possible [48, 49].

It is important to highlight some implications of the strategy used by F. Waleffe, notably the use of smaller and smaller numerical domains in order to "confine" the dynamics. It is not strictly true that such an operation is equivalent to confinement, since the physical domain is still formally infinite. Instead it amounts to strengthening correlations. In this study, the whole flow was interpreted as a *crystalline* structure. Jimenez & Moin used a similar strategy in order to prevent the existence of velocity eddies of size larger than a specified cut-off size [28], i.e. to make turbulence less multi-scale, whereas Waleffe and subsequent authors would focus on lower $Re$ and on making turbulence less "spatiotemporal". Choosing
<figure><img src="image_17.png" /><figcaption>FIGURE 1.16 – State portrait energy input vs dissipation for plane Couette flow, taken from Ref. [32].</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 1.17 – Three-dimensional projection of the state space of plane Couette flow, including several fixed points, their unstable manifolds as well as 5 periodic orbits, from Ref. [39].</figcaption></figure>

this cut-off in inner units (the natural unit for streak width in turbulent flows) rather than in outer units allowed them to bypass the reduction in Reynolds number and focus only on the reduction in domain size. This lead to the famous concept of "Minimal Flow Unit" [28], understood as the smallest domain
size able to sustain turbulent motion, or equivalently non-laminar flow (despite the absence of any strict minimisation in the optimal sense).

1.6.2 Finite-amplitude solutions

Let us first define what is meant by "finite-amplitude solution" here in order to remove any ambiguity. The historical perspective is perhaps useful: it starts with the vague turbulent notion of coherent structure, a word used from the 1970s to label any shape recognisable in the flow that has some spatial coherence over some time interval, as opposed to the purely statistical investigations of the time [50]. F. Waleffe, aware of the contradictory nature of the concept, introduced the notion of exact coherent structure' (ECS) in order to restore the link with the Navier-Stokes equations. Once that link restored, finding ECS in practice turned into the usual numerical problem of identifying solutions to the discretised Navier-Stokes equations. The adjective "finite-amplitude" is added here in order to emphasize that the solutions are *not* obtained via the other method popular at the time, namely weakly nonlinear analysis, and thus that no assumption is made on the smallness of the amplitude. As soon as the self-sustaining process (SSP) was identified by Waleffe and co-authors, the quest for finite-amplitude solutions in the same flow set-ups became a quest for illustrations of the very same SSP using mathematically well-defined solutions.

There is some ambiguity to the fact that, especially in the numericist's jargon, any unsteady solution $X(t)$ to an initial value problem is precisely also a "solution". Behind the label "finite-amplitude solution", the notion of "simple solution" has been progressively eclipsed. "Simplicity" is not an objective property, however it still makes sense to talk about a hierarchy of solutions starting from the simplest ones.

When focusing on the temporal dynamics, the simplest solutions to the Navier-Stokes equations (or more generally any PDE system) are steady states. For all times $t, t'$ they obey the steady property

$$ \mathbf{u}(x, y, z, t) = \mathbf{u}(x, y, z, t'). \qquad (1.17) $$

The next level in the hierarchy corresponds to solutions that are steady in a given moving frame of reference: travelling wave solutions (TWs) are defined by

$$ \mathbf{u}(x, y, z, t+T) = \mathbf{u}(x+c_x T, y, z+c_z T, t), \qquad (1.18) $$

for all $t$ and $T$. The $c$'s are interpreted as phase velocities.

Periodic orbits (POs) are defined as

$$ \mathbf{u}(x, y, z, t+T) = \mathbf{u}(x, y, z, t), \qquad (1.19) $$

for all $t$ but for given values of $T$.

Next in the hierarchy, relative periodic orbits (RPOs) are defined as

$$ \mathbf{u}(x, y, z, t+T) = \mathbf{u}(x+\sigma_x, y, z+\sigma_z, t), \qquad (1.20) $$

for all $t$ but for given values of $T$, $\sigma_x$ and $\sigma_z$. Note that travelling waves are degenerate relative periodic orbits as well as degenerate steady states.

It is possible to describe such special solutions independently of the notion of space. If the system is written under the generic form $\dot{\mathbf{X}} = \mathbf{F}(\mathbf{X})$ with a propagator $\phi^t$, then steady states are simply zeros of $\mathbf{F}$, or equivalently fixed points of $\phi^t$. As for travelling waves and relative periodic orbits, the shifts in space form a vector $\sigma = (\sigma_x, \sigma_z)$ of rank two. Spatial shifts can be reinterpreted as simple application of an operator $\mathcal{T}_\sigma$ on velocity fields. That way a TW can be redefined as any solution, for any $T > 0$ to

$$ \phi^T \mathbf{u} = \mathcal{T}_\sigma \mathbf{u}, \qquad (1.21) $$
whereas an RPO is defined as a solution to the same equation yet for a specific value of T > 0.

Many methods have been used in the past to find exact solutions to the Navier–Stokes equations, especially in the context of subcritical shear flows where the laminar solution, the only one known analytically a priori, is evidently no good description of a sustained turbulent flow. Some of these methods will be described in the next chapter. The original hope is that once a finite-amplitude solution is identified, it might serve as a reduced laboratory to investigate a turbulent flow without the difficulties inherent to very unsteady signals, such as arbitrary time-integration, choice of initial condition, need for ensemble averaging, etc... in other words a laboratory where the dynamics is fully reproducible. Since in the context of turbulence most finite-amplitude solutions of interest are expected to be unstable, a pedagogic cartoon has emerged progressively to describe the dynamics of a turbulent flow : as depicted in Fig. 1.21, the dynamics of the flow is encoded in an ergodic state-space trajectory **X**(*t*). The state space is populated by pointwise finite-amplitude solutions in the shape of points and cycles, all unstable except the fixed point corresponding to the laminar base flow. A turbulent trajectory, if any, corresponds to a trajectory staying away from the basin of the laminar state forever and visting the neighbourhoods of the individual ECSs one after the other. All other trajectories eventually converge to the laminar state and are labelled *transient*, possibly *transiently turbulent* if its *lifetime* is deemed sufficiently long. The picture of transition associated with it is two-fold. To start with, only initial conditions located in the right dynamical region are likely to appear turbulent, be it only transiently. Moreover, as the Reynolds number decreases the number of unstable solutions reduces and eventually drops to zero. In such a cartoon the complexity of the state space is mainly determined by the number of unstable solutions. It is worth repeating that this view, although faithful to an original view of turbulence by Hopf [51], and no matter how pedagogically helpful, is only a cartoon. In particular it is unable to shed light on such concerns as the fractality/strangeness of the attractor and is unable to yield any information about the invariant measure of the system.

<figure><img src="image_19.png" /><figcaption>FIGURE 1.18 – Cartoon of the state space in subcritical transition to turbulence. Each cross except the laminar one denotes a locally hyperbolic finite-amplitude solution together with its stable and unstable manifolds. Courtesy of T. K. Kreilos [52].</figcaption></figure>

1.6.3 Cycle Expansions

There is a deeper reason as to why the knowledge of simple exact solutions is important : cycle
expansions. The e-book by P. Cvitanović [53] is devoted to a detailed explanation of the mathematics
underlying this very general idea. Suppose an autonomous dynamical system of the form $\dot{x} = F(x)$ in arbitrarily high dimension, associated with the flow map $\phi^t$. Suppose now also that the system above verifies the axiom A property, i.e. that in the neighbourhood of each periodic or steady solution, the system can be linearised and decomposed into stable and unstable tangent eigenspaces. This hypothesis rules out central manifolds, non-hyperbolic effects and hence rules out the proximity to a bifurcation point. For a given system, the validity of the axiom A hypothesis is notoriously difficult to predict; at best it is usually assumed as a hypothesis. Under these assumptions, a formula originally suggested by Gutzwiller [53] expresses the dynamical averages of any observable as a weighted average over the entire set of periodic orbits of the system. No information other than the location of the periodic orbits, their period and their linear stability is required. [54]. By the dynamical average of an observable $A = A(x)$, we mean all the possible moments and cumulants obtained by averaging $A(x(t))$ along an ergodic trajectory $x = x(t)$ from time $t = -\infty$ to $t = +\infty$ [55].

More mathematically, if $A^t$ is an observable evaluated at time $t$, we can always introduce the transfer operator

$$ \mathcal{L}^t(x, y) = \delta(y - \phi^t(x))e^{\beta A^t(x)}, \qquad (1.22) $$

where $x$ and $y$ are two arbitrary points of state space. $\beta$ is an auxiliary variable whose main role is to enable one to recover the desired state space average by differentiation. The kernel $\delta(y - \phi^t(x))$ is the Perron-Frobenius operator which transports densities in time (cf Chap. 6 in Ref [53]). Note that taking the trace of a linear operator containing the factor $\delta(y - \phi^t(x))$ is equivalent to considering only the points $x$ where $\delta(y - \phi^t(x)) \neq 0$, i.e. periodic orbits whose period is an integer multiple of $t$. The cycle expansion formula states that

$$ \mathrm{tr}(\mathcal{L}^t) = \sum_p T_p \sum_{r=1}^{\infty} \frac{\delta(t - rT_p)}{|\det(I - \mathcal{J}_p^r)|}, \qquad (1.23) $$

where $\mathcal{I}$ is the identity matrix, $\mathcal{J}_p$ the Jacobian at the system linearised around the periodic orbit of interest, $r$ the multiplicity of the orbit and $T_p$ its prime period.

More intuitively, the above formula is a generalisation of the trace formula $\mathrm{tr}(\log(\mathcal{L})) = \log(\det(\mathcal{L}))$ for a matrix $\mathcal{L}$. On the left-hand side of Eq. 1.23 we have the statistical information expressed as function of the right-hand side, i.e. the individual deterministic elements. Indeed the spectral theory applied to the operator $\mathcal{L} = \mathcal{L}^t$ relates its trace to a series of Policott resonances $\gamma_\alpha$ of multiplicity $m_\alpha$, i.e.

$$ \mathrm{tr}(\mathcal{L}^t) = \sum_{\alpha=0}^{\infty} m_{\alpha} e^{-\gamma_{\alpha} t}. \qquad (1.24) $$

The exponents $\gamma_\alpha$ are related to the escape rates (for the leading eigenvalue) and correlation functions (for the exponents beyond the leading one), that are all key characteristics about the asymptotics of the system as $t \to +\infty$ [56]. For instance, if averages over the state space $\Omega$ are defined by

$$ \langle e^{\beta A} \rangle = \int_{\Omega} e^{\beta A(x)} d\mathbf{x}, \qquad (1.25) $$

then the average of $A$ is given directly by

$$ \langle A \rangle = \frac{\partial}{\partial \beta} \langle e^{\beta A} \rangle |_{\beta=0}, \qquad (1.26) $$

the standard deviation of $A$ is given by the second derivative with respect to $\beta$, etc...

Owing to Eq. 1.23, such information can be evaluated using individual averages over periodic orbits, and only the periodic orbits! The inclusion of steady states is discussed in [54]. The presence of additional
symmetries is discussed in [53] and results in periodic orbits being generalised into *relative* periodic orbits.

The main limitation in Eq. 1.23 is the infinite number of periodic orbits required to *cover the attractor*. Super-exponential convergence of the main averages versus *T* is usually observed for low-dimensional systems such as the Lorenz system when only periodic orbits up to period *T* are considered [57]. Nevertheless in high dimension this still leaves one with a daunting numerical task of determining *all* periodic orbits of period less than *T* without missing any of them.

Two remarks can be made at this stage. The first one is that unstable periodic orbits are generic and are expected to be found in all dynamical systems displaying chaos. An illustration is given by the construction of the Smale's horseshoe [58]. More than an expectation, the existence of unstable periodic orbits is sometimes listed as one of the conditions defining chaos, together with the condition that they are over the attractor [59, 60]. Beyond the numerical challenge associated with it, the mere identification of a given periodic solution is hence by itself not a theoretical breakthrough. In the context of the flow inside a MFU for instance, UPOs of the kind identified by [61] are in fact expected. UPOs are a consequence of the presence of chaos and do not "explain" chaos. In the same vein, periodic orbits likely to be found embedded within turbulent flows will not *explain* turbulence either. The second one is a matter of feasibility in a finite world : it remains far from clear whether such a heavy numerical machinery is truly progress if the main goal is to compute the mean flow profile. It seems that so far, and for a number of years to come, an efficient numerical solver will be much more efficient in predicting the flow statistics.

<figure><img src="image_20.png" /><figcaption>FIGURE 1.19 – a) Projection of a Poincaré section of the state space of pipe flow. Turbulent trajectories (grey) and relative periodic orbits (black, green and red), b) and c) successive zooms. Taken from [42]. d) same as c) with the unstable manifold of one of the periodic orbits depicted in blue.</figcaption></figure>
1.6.4 Chaotic saddles

A large part of the literature on dissipative dynamical systems celebrates the notion of attractor, or
equivalently of ω-limit set : the part of the state space that will again be visited by a given trajectory, no
matter how many times it has been visited by the same trajectory already. Attractors are consistent with
i) the dissipative nature of the Navier–Stokes equations and ii) the notion of sustained dynamical regime.
A clearer link between the two notions (attractors and turbulence) emerged in the 1960s pioneered by
E.A. Lorenz, initially with thermal convection in mind, for which turbulence is indeed sustained as long as
thermal forcing is on. Direct application of this analogy to open shear flows leads to surprising results. At
large enough Reynolds number, it is certainly more straightforward to consider turbulence as sustained,
although even there some concepts become shaky due to the main advection. This is at least the case
in spatially periodic simulations of channels and pipes in which, once initiated, turbulence never dies
out. At lower Reynolds numbers, close to the onset of appearance of turbulent dynamics, the situation
is less clear. Although some initial conditions seem to lead to a sustained turbulent flow, simulations
over longer and longer observation times sometimes reveal relaminarisation. Unless transition is triggered
again, the autonomous system stays laminar forever and the turbulent episode is nothing more than
a transient. Experimentally (i.e. without periodic boundary conditions) this has a counterpart, easiest
understood in pipe flow : a turbulent patch can travel several thousands of pipe diameters and die over
the space of only a few pipe diameters [62–64]. In spatially extended systems such as plane Couette flow,
although turbulence can spread, it can also recess until the whole flow becomes laminar [65, 66]. This
relaminarisation is hence not only a numerical artefact but a genuine physical phenomenon, although
again the mechanisms at play in the presence of spatial periodicity might differ from those in free space.
Theoretical efforts to understand this relaminarisation phenomena have led to many interesting analo-
gies with other areas of theoretical or applied physics [67]. To begin with, the lifetimes corresponding
to multiple laboratory or computer experiments all pointed towards statistical distributions. Moreover
the cumulated distributions of lifetimes frequently appear exponential, indicating a *memoryless process*.
Although an analogy with radioactivity was unavoidably made, a more fruitful analogy turned out to be
classical scattering, seen as a model of idealized unimolecular fragmentation [68]. The same formalism
turned out later to be useful in other areas of hydrodynamics, namely fluid mixing [69], [70]. In the
scattering picture, trajectories bounce alternatively on an array of hard spheres in a chaotic fashion, until
they escape a given region of the associated state space. Despite the lack of dissipation, this picture of
state space in figure 1.20 is very close to that depicted schematically in Fig. 1.21. It is structured around
not one, but rather several invariant sets. In particular the rate of leakage of a delayed trajectory from the
interaction region is claimed to be analogous to the rate of decay of a metastable system. This analogy
has given a yet stronger weight to the cartoon of Fig. 1.21 : the lifetime of individual turbulent realisa-
tions corresponds to the time it takes for the corresponding trajectory to bounce through the "jungle"
of unstable solutions and eventually to escape from it. The more entangled the network the longer the
mean lifetime. A formalism of Markov chains has even been developed to account for the memoryless
statistics observed [71, 72]. The initial success of this new dynamical picture came however from i) its ori-
ginality and ii) its ability to incorporate well-chosen elements from dynamical systems theories to explain
new features of low-Reynolds number turbulence. A question that remained, however, unanswered by this
approach is about the possible divergence of turbulence's mean lifetime as the Reynolds number increases.

1.6.5 Other invariant sets

The historical focus on a set of fixed points and periodic orbits has perhaps slightly shadowed the
search for other invariant sets beyond them. Invariant manifolds are as important as fixed points in many
dissipative systems. In particular stable manifolds of objects yet-to-be-identified often act as inpenetrable
barriers to trajectories, and hence as proper state space frontiers when their codimension is one. This
will be the case for the edge manifold that will be extensively reported in the next chapter. It is as yet
unknown whether other invariant manifolds can be of interest to the description of the state space for
this class of hydrodynamic systems.
<figure><img src="image_21.png" /><figcaption>FIGURE 1.20 – Cartoon of two trajectories bouncing periodically between three repelling invariant sets. Taken from *Chaosbook.org*. Any other trajectory leaving the game after a finite number of rebounds is interpreted as a chaotic transient.</figcaption></figure>

Other orbits (understood as 'trajectories') are however of interest. The ECSs reported above correspond almost exclusively to states with both a stable and an unstable eigenspace. As a direct consequence, heteroclinic orbits leaving the neighbourhood of one given state at $t = -\infty$ and connecting with the neighbourhood of another one for $t \to +\infty$ are likely to exist. Their properties are important as they are believed to guide the dynamics of most trajectories. They have been investigated numerically in Refs [73–75].

Homoclinic trajectories are also important as a support for state-space trajectories, however their numerical identification is potentially even richer in consequences. The Smale–horseshoe, a cornerstone of dissipative dynamical systems theories [58, 76], states that the existence of a homoclinic orbit from an unstable periodic orbit automatically implies chaotic dynamics. This strong theorem has been used by [77, 78] to justify the chaoticity of plane Couette flow in the HKW cell. Interestingly, the Smale Horseshoe theorem is not expected to be valid if the periodic orbit is replaced by a steady state or a travelling wave. Given the number of unstable steady states identified in shear flows at the moment rather than UPOs, one can wonder what implication it has for the dynamics.

Another conceptually easy (but practically heavy) generalisation concerns the inclusion of unstable quasiperiodic states into the same picture. Such new solutions can emerge from Neimark–Säcker bifurcations of periodic orbits, and they can also bifurcate out of the blue. While formally the repelling role played by two-frequency solutions cannot be ignored, it remains to be shown how such information can be quantitatively useful : there is no equivalent of cycle expansions taking two-tori solutions into account. However exact toroidal solutions with more than two temporal frequencies are likely to play an increased role at slightly higher Reynolds numbers where the flow starts to display multiscale properties [79].
<figure><img src="image_22.png" /><figcaption>FIGURE 1.21 – Two homoclinic orbits of plane Couette flow starting from a time-periodic edge state (colours), together with the PDF of individual trajectories (grey). Representation in a kinetic energy/dissipation state portrait. From [77].</figcaption></figure>
2

Dynamical systems description
## 2.1 The edge of chaos

Among the concepts from dynamical systems theories developed in the preceding introductory chapter, the recently formalised concepts of edge manifold and edge state have been left aside. The present chapter represents a personal opportunity to re-define these concepts from scratch, with the relevant mathematical rigour and with a (hopefully) consistent vocabulary.

### 2.1.1 The edge as a smooth separatrix

#### 3.1.1.1 Illustration on a low-dimensional model

We consider first, mainly for pedagogic purposes, the two-dimensional autonomous model suggested by O. Dauchot and P. Manneville [17], which we label DM2D. This system is written under the general form $\dot{x} = Lx + N(x)$, where $x \in \Omega = \mathbb{R}^2$. L is a typical non-normal matrix suggested by Trefethen [80]. The components of N consist of quadratic forms that conserve energy such that $<N(X), X> = 0$ for all X in line with the Reynolds-Orr theorem. These two operators read respectively :

$$L = \begin{bmatrix} s_1 & 1 \\ 0 & s_2 \end{bmatrix}, \quad N(x_1, x_2) = \begin{bmatrix} x_1 x_2 \\ -x_1^2 \end{bmatrix}.$$

A state portrait of the DM2D model is displayed in figure 2.1 for the parameters $s_1 = -0.1875$ and $s_2 = -1$.

<figure><img src="image_23.png" /><figcaption>FIGURE 2.1 – State space of the DM2D model [17]. The laminar state $x_L$ and the turbulent state $x_T$ are represented together with their respective superstable manifolds, $W^{ss}(L)$ and $W^{ss}(T)$. The dividing line in red represents the edge manifold and coincides with the stable manifold of the edge state $x_E$, $W^s(x_E) = \Sigma$. The blue trajectories form the oriented unstable manifold of E, $W^u(x_E)$. The green diamond corresponds to the minimal seed (see Section 3.3). Taken from Ref. [81].</figcaption></figure>

The description of the state space is straightforward : two attracting fixed points $x_L$ (the "laminar state") and $x_T$ (the "turbulent state") coexist. The fact that the turbulent state has no proper "turbulent" dynamics does not matter much for the present illustration. By construction $x_L$ is linearly stable, and it can be checked analytically that $x_T$ is also linearly stable. Each attractor hence possesses its own
attraction basin, respectively $B_L$ and $B_T$. Their intersection is non-empty and we define the edge manifold as the intersection of their closures :

$$ \Sigma = \overline{B_L} \cap \overline{B_T}, \qquad (2.1) $$

displayed in red in figure 2.1. Trajectories starting at initial time from any point in $\Sigma$ stay on $\Sigma$, hence $\Sigma$ is an invariant set for the dynamics [82]. In addition, all such trajectories converge to a third fixed point labelled $x_E$. We have the identity $\Sigma = W^s(\{x_E\})$. $x_E$ is labelled the "edge state", and it is here a point of zero dimension. It should not be confused with the "edge manifold" $\Sigma$ (sometimes simply referred to as "edge") which here has dimension 1, in other words $\Sigma$ is a hypersurface of codimension one. As the governing parameter $s_1$ is decreased, the two attractors $x_E$ and $x_T$ merge together and disappear at the same time, leaving the laminar state $x_L$ alone (as fixed point) in the state space.

Within the DM2D model, the transition process from the laminar state to the turbulent one is already rich enough to be compared to the basic transition process as observed in hydrodynamical systems :

- initial conditions close inside $B_L$ can approach $x_E$ transiently but end up in $x_L$

- initial conditions close inside $B_T$ can approach $x_E$ transiently but end up in $x_T$

- the closer the initial condition to $\Sigma$, the longer the transient visit to $x_E$.

This indicates that $x_E$, although unstable with respect to two-dimensional disturbances, is stable and attracting within $\Sigma$ and unstable in any direction transverse to $\Sigma$. Numerically, $x_E$ can be hence found by a simple one-dimensional search (labelled also *shooting method* in the early Japanese literature). Such a search only requires initial conditions within a one-dimensional segment that intersects $\Sigma$ somewhere. The concept is illustrated in figure 2.10. This is an iterative process terminated once the two initial conditions one either side of $\Sigma$ are considered close enough. Note that in dissipative systems this algorithm *always* identifies an edge state, regardless of the associated dynamics.

<figure><img src="image_24.png" /><figcaption>FIGURE 2.2 – Illustration of the "classical" iterative bisection process (courtesy M. Beneitez). Grey dashed lines: lower and upper bounds chosen for an observable $a(t)$. The sets of red and blue curves bracket an edge trajectory from $t=0$ until at least $t=500$.</figcaption></figure>

### 3.1.1.2 The edge seen as a Lagrangian Coherent structure in high dimension

It is visually clear from figure 2.1 that the manifold $\Sigma$ divides the state space into two complementary parts. No initial condition located strictly within one attraction basin can end up in the other basin, which makes $\Sigma$ a barrier to transport. In addition, the hyperbolic nature of $x_E$ is responsible for repelling
all nearby trajectories, even those infinitesimally close to Σ. All these conditions are strongly reminiscent of the physical concept of hyperbolic Lagrangian Coherent Structure, where Lagrangian refers to the tracking of individual trajectories [83]. The original frame for which LCSs were developed is also motivated by hydrodynamics, namely the finite-time transport of fluid particles by known time-dependent fluid flows [84, 85]. In particular, in steady velocity fields, repelling LCSs essentially mean material barriers to transport. Other applications, for instance chemical reactions, have also been considered [86]. The class of systems under study are however very different from the current context in transition : the dynamics are always conservative, there are no attractors, the trajectories can be computed in either forward or backward time, and the dimension of the state space (which coincides with the physical space) is very low, typically two. Besides, the system is in general non-autonomous due to the time dependence of the vector field. In most applications of Lagrangian concepts, the state space associated with Lagrangian tracers (governed by an equation $\dot{x} = v$) and the physical space coincide. In the Eulerian point of view relevant to the present transition problem, both spaces differ radically in their dimension.

There are several definitions used to identify repelling LCSs. We consider here only two of the simplest quantifiers : Finite-time Lyapunov exponents (FTLEs) and Lagrangian Descriptors (LDs).

## Finite-time Lyapunov exponents

Finite-time Lyapunov exponents (FTLEs) have been popularised as a tool for the identification of LCSs in [84, 87]. Let $C_{t_0}^t$ be the (positive definite) Cauchy-Green tensor defined by

$$C_{t_0}^t = (\nabla \phi_{t_0}^t)^* \nabla \phi_{t_0}^t, \qquad (2.2)$$

where $(\cdot)^*$ represents the Hermitian transpose and the derivatives $\nabla$ are with respect to state space variables. The Cauchy-green tensor has size $n \times n$, and its entries are usually computed using second-order centered finite-differences [83]. If $\tau = t - t_0 > 0$, the average growth rate over the time interval $(t_0, t_0 + \tau)$ around an initial condition $x_0$ is given by $\sqrt{\lambda(t_0, t_0 + \tau, x_0)}$, with $\lambda$ the largest eigenvalue of the diagonalisable matrix $C_{t_0}^t$. The present system is autonomous hence $t_0 = 0$ without loss of generality. The leading finite-time Lyapunov exponent of the system at position $x_0$ is given by

$$\Lambda(x_0, \tau) = \frac{1}{\tau} \log \sqrt{\lambda(0, \tau, x_0)}. \qquad (2.3)$$

The ridges in the field of the leading FTLE at time $\tau$ can be used as a diagnostic of a hyperbolic LCS [83, 88]. Hyperbolic LCSs refer to attracting and repelling distinguished invariant manifolds [84], with only the forward time FTLEs relevant for the identification of repelling LCSs [89]. FTLEs can be computed in arbitrarily high dimension using the recent scalable algorithm based on Optimally Time Dependent (OTD) modes, which are themselves determined from a minimization principle [90]. Under generic conditions these OTD modes converge exponentially fast to the eigendirections of the Cauchy-Green tensor associated with the largest eigenvalues, i.e. the largest finite-time Lyapunov exponents [91].

## Lagrangian Descriptors

Lagrangian Descriptors (LDs) are a more recent diagnostic for Lagrangian coherence introduced in [92, 93] and its theoretical framework further developed in [94]. Like FTLEs, LDs are based on the integration of a given observable along trajectories, however they do not request any differentiation of the flow map with respect to the initial condition. The original quantity of interest is

$$M(x_0, t_0, \tau) = \int_{t_0-\tau}^{t_0+\tau} g(x(t))dt, . \qquad (2.4)$$

The observable $g$ can be defined as

$$g(x(t)) = \sum_{i=1}^{m} |f_i(x, t)|^p, \qquad (2.5)$$
where the $f_i$'s are the components of the velocity field $f$, and $p \in (0, 1]$ and $\tau \in \mathbb{R}^+$ are two parameters. $p = 0.5$ is a good choice in practice. Eq. 2.4 can be split into its forward and backwards contribution [86, 95]. For the present autonomous system, again $t_0 = 0$ without loss of generality, and these contributions write hence respectively

$$M(x_0, \tau)^+ = \int_0^\tau g(x(t))dt, M(x_0, \tau)^- = M(x_0, \tau) - M^+(x_0, 0, \tau). \quad (2.6)$$

Due to the dissipative nature of the system, numerical backwards integration is ruled out for numerical stability reasons, only $M^+$ can be considered here. Since our focus is on stable manifolds rather than unstable ones it is sufficient to focus on the computation of $M^+(x_0, t_0, \tau)$. The abrupt change also means that the derivative of $M^+$ transverse to the boundaries is discontinuous. These singular features in LD plots are often connected to the stable and unstable manifolds of saddle points. Depending on the case, it is also useful to quantify abrupt changes of $M^+$ via the $L^2$-norm of its gradient

$$B(x_0, \tau) = \left[ \sum_{i=1}^{n} \left( \frac{\partial M^+}{\partial x_{0,i}}(x_0, \tau) \right)^2 \right]^{1/2} . \quad (2.7)$$

<figure><img src="image_25.png" /><figcaption>FIGURE 2.3 – State space of the DM model represented using isolevels of a) the leading FTLE $\Lambda$ b) $M^+$ and c) $B = ||\nabla M^+||$. In each case the time horizon is $T = 60$. The edge manifold $\Sigma$ is highlighted by the ridge in a) and c).</figcaption></figure>

Other diagnostics for LCSs are routinely used, see [88] for a comparative review. Alongside the diagnostics based on a scalar field, such as FTLEs and LDs, other approaches based on transfer operators or dynamic Laplace operators seek coherent structures by formulating rigorous mathematical coherence principles (see e.g. [88, 96]). These diagnostics generally display limited scalability and are not considered here.

### 3.1.1.3 Illustration on hydrodynamic examples

The first pioneering use of the concepts of edge manifold and edge state in hydrodynamics is probably due to S. Toh and T. Itano [97]. The system they investigated is plane channel flow in a computational cell relatively close to that considered in [28]. Their pedagogic sketch, displayed in figure 2.4 shares many ressemblances with figure 2.1. The edge state identified in their study was claimed to be travelling wave solution. It was later re-interpreted, based on a longer computational, as a very long (relative) periodic orbit [98]. Eventually, it was demonstrated using better resolution and even longer edge tracking that the edge state in that system is in fact neither steady nor periodic or quasiperiodic, but apparently chaotic [99]. What was originally interpreted by Toh & Itano as the instability of the edge state was in fact part of the edge state dynamics itself. Note that in the same lab in Kyoto, G. Kawahara and S. Kida also used the same computational method in their 2001 paper on periodic orbits but in plane Couette flow : the edge state identified is the "gentle" periodic orbit visible in state portraits of the previous chapter.
<figure><img src="image_26.png" /><figcaption>FIGURE 2.4 – Sketch of state space from Toh & Itano 2001 [97]. The dividing hypersurface labelled "Basin Boundary" represents the edge manifold.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 2.5 – Three-dimensional snapshots of the unsteady velocity field for the relative periodic edge state, as identified in channel flow by Toh & Itano 2003 [98]. Light grey : isocontours of streamwise velocity, darker colours : isocontours of the vortex criterion *Q* coloured by the sign of the streamwise vorticity.</figcaption></figure>

2.1.2 Concept of edge state in chaotic saddles

In the previous example of the DM2D model, all fixed points but the edge state are attractors, and
the edge manifold appears smooth. As previously mentioned, a key property of transitional shear flows
close to their onset is the finite lifetime of the *turbulent* dynamics, observed both in experiments and in
computations [100]. This property appears incompatible with the definition of a "turbulent" attractor,
hence the only attractor left in this picture is the laminar state. The notion of attraction basin hence
needs to be generalised to that of a separatrix between *short transients* and *supertransients* [67]. The
generalisation works in practice since, for observation times shorter than the mean lifetime, the dyna-
mics is indistinguishable from the bistable dynamics detailed earlier. Similarly, most usual quantities
such as Lyapunov exponents, correlation functions etc.. can be computed over finite times only [67]. Se-
vere practical problems arise however if the observation time is longer or comparable to the mean lifetime.

There has been a strong interest in transient chaos by theoreticians, going back perhaps as early as
the foundations of complex dynamics and Fatou's dust in iterative maps in the early 20th century. Fatou's
dust is the union of the attraction basins of periodic points in simple maps such as the tent map, which
is known to have the structure of a Cantor set and whose construction is illustrated in figure 2.6.

There is a direct analogy with the edge manifold in smooth differential systems. Since the turbulent
<figure><img src="image_28.png" /><figcaption>FIGURE 2.6 – Iterative construction of the Cantor set of Fatou’s dusts of the tent map (courtesy P. Manneville).</figcaption></figure>

dynamics itself does *not* define an invariant set, the edge manifold is regarded as the only remaining
non-trivial invariant set of the system. It can hence be sought directly as a set characterised by infinite
lifetimes, i.e. the ridge in a lifetime landscape as in figure 2.7. Most related algorithms use the idea that
any segment joining two regions of the state space with different lifetimes *must* straddle the edge. Luckily
they do not differ from standard bisection algorithms [101–103].

<figure><img src="image_29.png" /><figcaption>FIGURE 2.7 – Chaotic saddle of the Henon attractor revealed by the ridge in the lifetime landscape. [103].</figcaption></figure>

Until 2005 the edge state was reported under very different names such as "chaotic saddle". The
introduction of the word *edge* (in fact, originally *edge of chaos*) dates back to the reference paper of J.
Skufca, B. Eckhardt and J. Yorke in 2006, who focus on a 9-dimensional Galerkin model of subcritical
transition. It has remained popular ever since. We would like to emphasize the distinction between the *edge* $\Sigma$ (the manifold/separatrix/barrier) and the *edge state* $E$, which is the limit state of a trajectory within $\Sigma$. $E$ is sometimes labelled the *relative attractor* on $\Sigma$. Note that, although the turbulent state is often considered unique, no principle warrants the uniqueness of the edge state for a given system. In a few cases several distinct edge states have been identified for the same set of parameters.

It did not take long until the above concept was used in a genuine Navier-Stokes context, starting with pipe flow [104–106], in situations where the mean turbulent lifetimes are either finite or unbounded. From a dynamical point of view, the edge state can host relatively simple time variations or look chaotic, however the time variations look modest compared to the associated turbulent dynamics. Experience suggests that the chaoticity of edge states is the exception rather than the rule. From a structural point of view, the perturbation velocity field associated with edge states possesses a few robust features : streamwise rolls and streaks are always present, just as for the turbulent state (see figure 2.8). Edge states are never strictly two-dimensional and display streamwise variations for finite $Re$. These elements are all consistent with the SSP picture of F. Waleffe [107]. As $Re$ goes to infinity (provided the numerical domain is unchanged), the corresponding coherent structures evolve towards a two-dimensional flow different from the laminar one, whereas the streamwise-dependent part becomes dominated by a wave part traveling at constant speed. It was later realised [108] that this is perfectly consistent with the Vortex-Wave Interaction picture suggested by P. Hall and co-workers in the 1980s to address the question of high-$Re$ shear turbulence [109]. This theory involves critical layers, whose thickness scales like $O((\alpha Re)^{-1/3})$, such that they become singularly thin in the high-$Re$ limit [110]. This apparent issue will be mentioned again in the last section of this chapter when dealing with localised edge states.

<figure><img src="image_30.png" /><figcaption>FIGURE 2.8 – Snapshot of the edge state in pipe flow with parameters $Re = 2875$ and domain length $L = 5D$, from Ref. [105]. Top : cross-sectional velocity field, bottom : three-dimensional side view (flow from left to right).</figcaption></figure>
### 2.1.3 Edge tracking techniques

#### A posteriori verification

The methods to identify numerically edge states will be referred to as "edge tracking techniques". Several methods coexist in the hydrodynamic literature with their pros and cons, although edge tracking tends to be dominated by the "classical" bisection method. The first edge states in shear flows were found accidentally even before the separatrix concept emerged. It simply turned out that, out of the multitude of finite-amplitude solutions identified numerically (no matter how), some of them have the edge state property that their stable manifold divides the state space into two separate parts. The simplest finite-amplitude solutions to date usually emerge in saddle-node bifurcations, and form pairs with one "lower-branch solution" and one "upper-branch solution" (the distinction is usually based on the perturbation kinetic energy). If one of these states has, even for a bounded parameter range, strictly one unstable eigenvalue (equivalently one unstable Floquet multiplier if the state has time periodicity), then it is very likely to be an edge state because its stable eigenspace has codimension one ($dim(E^u) = 1$). This is the case for instance with the steady lower-branch Nagata state, the first finite-amplitude solution found in plane Couette flow [111]. It was verified multiple times that i) slight perturbations to the edge state lead either to turbulent or laminar dynamics depending on the sign of the perturbation, ii) classical bisection finds this steady solution as the edge state, and iii) for low enough $Re$ this is the only such state for the corresponding numerical parameters [112]. Conversely, finite-amplitude solutions with two or more distinct instability directions ($dim(E^u) > 1$) are not eligible as edge states, however they can belong to $\Sigma$ or not.

<figure><img src="image_31.png" /><figcaption>FIGURE 2.9 – Sketch of local stability near hyperbolic points. a) $dim(E^u) > 1$ b) $dim(E^u) = 1$. Only in the case b) can the fixed point be an edge state.</figcaption></figure>

#### Bisection techniques

The core of popular edge tracking techniques is the bisection process. Geometrically speaking, a one-dimensional state space segment is first identified such that it straddles the laminar-turbulent separatrix $\Sigma$. The straddling property is identified simply by monitoring whether the flow is laminar or turbulent after some observation time. The segment is iteratively shrunk until it coincides with a state on $\Sigma$ as in figure 2.10. The process is always convergent as soon as the length of the segment can shrink to zero, irrespective of the dynamics obtained. It occurs however in practice that machine precision is reached before the bisection has led to the desired or expected dynamics. In such a case the process is simply restarted by considering another pair straddling $\Sigma$ (usually chosen from the end of the previous steps). Most searches employ this restarted technique in order to access the true asymptotic state without being fooled by finite arithmetic.

A common misunderstanding consists in claiming that edge tracking algorithms fail to converge if
<figure><img src="image_32.png" /><figcaption>FIGURE 2.10 – Sketch of edge tracking from Ref. [113]</figcaption></figure>

the end result is dynamically "complex". In addition, the seducing advantage of simple edge states is the possibility to monitor the convergence towards them. If the edge state is instead a chaotic relative attractor, it is much more delicate to assess quantitatively when the asymptotic state is 'reached' or to claim convergence.

Each bisection algorithm can be summarised as a recursive loop :

**Step 0** $k = 0$, choose a perturbation to the laminar flow $u_0$

For $k=0,1,2,...$ repeat until machine precision is reached :

**Step 1** $u_k(t=0) = \lambda_k u_0$

**Step 2** $\lambda_{k+1} = \frac{1}{2}(\min(\lambda_k, \lambda_{k-1}) + \max(\lambda_k, \lambda_{k-1}))$

**Step 3** if $u_k(t)$ ever becomes turbulent (even transiently), $\lambda_{k-1} \leftarrow \lambda_k$ and $\lambda_k \leftarrow \lambda_{k+1}$

else $\lambda_k \leftarrow \lambda_{k+1}$.

Bisection algorithms differ according to the method chosen to select the final observation time.

The "classical" bisection technique used here is based on the choice of one observable $a = a(t) = a(x(t))$ together with the choice of bounds $a_L$ and $a_T$. It is then assumed, based on physical intuition, that the dynamics has no choice but to relaminarise if $a(t) < a_L$, or no choice but to go turbulent if $a(t) > a_T$. A transient time is usually excluded at the beginning of each run before the previous property $a_L < a(t) < a_T$ is checked for. Note that a good choice of $(a_L, a_T)$ can in practice really accelerate the process by several orders of magnitude. Ideally the choices $a_L = \min\{a(x), x \in \Sigma\}$ and $a_T = \max\{a(x), x \in \Sigma\}$ are ideal, but such information is not necessarily available beforehand. If the amplitude of the edge $a_E$ is known (even when time-dependent), this can lead to further improvements
in speed.

One of the issues with classical bisection algorithm is that each step is longer than the previous one. The time needed for each iteration diverges in principle as the relative attractor is approached. This is due to the increasing transient time spent in the neighbourhood of the relative attractor before the bounds are reached. Restarting from states ever closer to $E$ can limit this divergence but not fully prevent it. The interpretation of $\Sigma$ as a Lagrangian Coherent Structure can here provide a solution : bisection now refers to the classical maximisation (or minimisation) of a given state space functional, in the hope to locate ridges in a given landscape. This functional can be for instance the leading FTLE computed over a specified time horizon $\tau$. It can also be, when using Lagrangian Descriptors, the gradient $B^+$ of $M^+$ introduced in Eq. 2.7 for a specified time horizon $\tau$.

The performances of these three different bisection algorithms have been compared in figure 2.11a for the simple two-dimensional DM2D model introduced earlier. All three algorithms converge exponentially

<figure><img src="image_33.png" /><figcaption>FIGURE 2.11 – DM2D model. Left : convergence in log scale of the three types of bisection versus the horizon time $\tau$. Right : LCS diagnostics as a function of $x_2$.</figcaption></figure>

fast towards the edge state $E$. The FTLE-based bisection outperforms the classical bisection algorithm although the improvement is not spectacular. The LD-based method, however, is much faster than the two other ones. This suggests the use of such a variant of the bisection technique for future studies. We note that the ridge structure on $\Sigma$ characterises geometrically the edge manifold only when the turbulent state is a simple fixed point. The computations above have been generalised to investigate the structure of $\Sigma$ in minimal flow units of plane Couette flow at $Re = 400$. There the landscape of both FTLEs, $M^+$ and $B$ is slightly more difficult to apprehend. In particular, turbulence for these parameters is transient : following the $\Lambda$-lemma [114], the edge manifold crosses the straddling segment multiple times. This leads to the jagged landscape depicted in figure 2.12. It is not clear whether in such cases the new bisection techniques offer any advantage compared to the original method considered by Skufca et al.

### Feedback control methods

Bisection algorithms are conceptually simple but, because they are based on multiple simulations, they can be costly and slow. Indeed one bisection requires $2^n$ individual runs until machine precision $\varepsilon_M$ is reached at $n^{th}$ iteration, i.e. $n \sim -\log(\varepsilon_M)/\log 2$, at which point the process needs to be restarted $M$ times as long as no simple edge state has been reached. The total cost is hence of $Mn \ge 100$ runs at least, with each run getting longer as accuracy improves. For instance, in [115], each bisection required a total of $Mn \approx 400$ runs, which corresponds to $\approx 10^6 - 10^7$ CPU hours. This high cost makes parametric studies infeasible in practice. Other numerical methods have been suggested as alternatives to the bisection-rootfinder combination, e.g. iterative adjoint optimisation methods [112, 116] though they
<figure><img src="image_34.png" /><figcaption>FIGURE 2.12 – Plane Couette flow computations. Left : One-dimensional landscape of state space rendered using LCS indicators (FTLEs, B) for $\tau = 300$. Right : two-dimensional cut through state space visualised using iso-levels of FTLEs. A and $\alpha$ are two coordinates centered on the laminar state and defined in Ref. [81].</figcaption></figure>

involve significant mathematical and computational complexity.

Since the goal is to pinpoint and converge towards one of the unstable regions of the state space, refined control strategies can achieve a similar task. This is done by introducing into the original system a control term that counteracts the edge instability and is able to stabilise unstable states without altering them significantly. Several theoretical control strategies have been suggested in the past in order to converge towards unstable states, some of the notable ones are Pyragas delay-based control [117] and Selective Frequency Damping [118]. The former is heavy to implement in high dimension because of memory requirements and suffers from some limitations regarding convergence. The latter is easy to implement and is convenient at least for steady states. However it turns out not to converge in shear flows, because eigenvalues with zero imaginary part are not stabilised (see e.g. [119]).

We propose a remarkably simple linear feedback control able to constrain the dynamics to the edge manifold, and to stabilise invariant solutions that are stable within it. This scheme has been applied to track unstable chimera states in systems of non-locally coupled phase oscillators [120, 121]. In its most simple form, the control scheme makes one of the system parameters $\mu$ state-dependent by imposing a linear relation between the parameter and an observable $A(t)$.

$$ \mu(t) = \mu_0 + \kappa(A_0 - A(t)) \qquad (2.8) $$

For the shear flow problems we will identify $\mu$ with Re and A with some pre-defined perturbation energy

$$ A(t) = \int_{r_i}^{r_w} (u^2 + v^2) r \, dr \, d\theta \, dz. \qquad (2.9) $$

Geometrically, the resulting dynamics can be understood by considering its representation in the $(\mu, x)$ plane. For $\kappa = 0$ we are back to the uncontrolled system, where for a fixed choice of the parameter the dynamics is constrained to a vertical line in the $(\mu, x)$ plane, determined by the choice of the parameter. The equilibria for each fixed parameter $\mu > 0$ are selected as the intersections with the corresponding vertical line. The equilibria on the upper branch of the parabola are dynamically stable, while those on the lower branch are unstable (cf figure 2.13a). By choosing the pivot point $(\mu_0, x_0)$ fixed and varying the control gain $\kappa$ we can now sweep the stable equilibrium of the controlled system along the parabola, until the straight line becomes tangential to the parabola and both equilibria disappear in a saddle-node bifurcation (cf figure 2.13b).
<figure><img src="image_35.png" /><figcaption>FIGURE 2.13 – Sketch in an (*A*, *Re*) representation of the stability of the uncontrolled (a) and controlled system using the feedback control method (b) based on the constraint (2.8) corresponding to the dashed slanted lines. The labels ‘UB’ and ‘LB’ refer to the upper and lower branch in the uncontrolled system, respectively. In fig. (b) the dashed slanted line is rotated around a pivot point at a given rate, allowing the unstable fixed point on LB to be tracked dynamically.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>FIGURE 2.14 – (a) Cross section of the stabilised TW solution with $\alpha = 1.25$ and $m_0 = 2$. Streamwise velocity perturbation in colour (from dark to white) and cross-stream components (vectors). (b) Eigenvalue spectrum for the TW solution with (blue circles) and without control (red crosses).</figcaption></figure>

This control scheme was applied to the case of pipe flow investigated in ref. [106]. Within a few hours of computation it manages to stabilise the lower-branch travelling solution which is the edge state for that system. Not only does it identify this unstable state for one parameter value, but it can trace a large part of the corresponding branch of TWs in the same computation. Although computationally very appealing, the suggested method converges exactly only when the edge state does *not* display any oscillatory instability. It would be desirable to explore further the properties of this control scheme as well as the possibilities to use unconverged results combined with another algorithm [122].

## 2.2 Edge state dynamics in terms of finite-amplitude solutions.

In the present section, we consider a finite-time segment from an edge trajectory computed in pipe flow, and analyse how the edge fluctuations are structured by the presence in state space of neighbouring finite-amplitude solutions. This analysis is based on classical recurrence analysis traditionally used in chaotic theory and signal processing [123].
## 2.2.1 Recurrences

We begin by defining the notion of spatial shift by considering a velocity field, independently of the time $t$ at which it is considered. $\sigma = (\sigma_{\theta}, \sigma_{z})$ is a pair containing the values of the shifts in the axial and azimuthal directions respectively (the generalisation to Cartesian coordinates is trivial and intuitive). Let $T_{\sigma}$ act on $\mathbf{u}$ such that

$$ (T_{\sigma}\mathbf{u})(r, \theta, z, t) = \mathbf{u}(r, \theta + \sigma_{\theta}, z + \sigma_{z}, t). \qquad (2.10) $$

Several types of finite-amplitude solutions are identified during this study. In particular, travelling wave solutions (TWs) are defined by

$$ \mathbf{u}(r, \theta, z, t + T) = (T_{\sigma}\mathbf{u})(r, \theta, z, t), \qquad (2.11) $$

with $\sigma = (c_{\theta}T, c_{z}T)$ for a specific couple $(c_{\theta}, c_{z})$ and for all $t$ and $T$. Relative periodic orbits (RPOs) are defined by the same relation (2.11) for all $t$, which however holds only for given values of $T$, $\sigma_{\theta}$ and $\sigma_{z}$.

We define now a recurrence function depending on the two time arguments $t_i$ and $t$

$$ r(t_i, t) := \frac{|\mathbf{X}(t) - T_{\sigma}\mathbf{X}(t_i)|}{|\mathbf{X}(t_i)|} \qquad (2.12) $$

Identifying all the minima in $t$ of the recurrence function yields $r_{min}$, a function of only $t_i$, defined as

$$ r_{min}(t_i) := \min_{t>t_i}\{r_i(t) : \frac{\partial r_i}{\partial t} = 0\}. \qquad (2.13) $$

The signal of $r_{min}$ is depicted in figure 2.15 for $\sigma = (0,0)$. The value of $\sigma_{\theta} = 0$ is motivated by computational simplicity. The value $\sigma_{z} = 0$ is selected arbitrarily to alleviate, in the case of travelling waves, the degeneracy between $\sigma_{z}$ and the time period. $r_{min}(t_i)$ features an alternation of relative lows and highs, with low values reaching $10^{-1}$. Note that such values are substantially lower than the corresponding $r_{min}$ found in turbulent runs for the same system, which rarely dip under 0.2–0.3.

The main idea of recurrence analysis is that local minima of $r_{min}(t_i)$, associated with a time $T$ such that $r_{min}(t_i) = r(T, t_i)$, correspond to furtive approaches to RPOs with a period close to $T$. Such RPOs would in fact correspond to exact zeros of the recurrence function $r$. Such zeros can be identified numerically within machine precision, using adaptations of classical rootfinding techniques. Rootfinders are iterative algorithms which need an initial condition, and the state corresponding to each dip of $r_{min}$ is precisely a good choice of initial condition. The classical Newton method, by far the most used rootfinder in this context, seeks for instance zeros of the functional

$$ g : \mathbf{X} \to |\phi^T(\mathbf{X}) - T_{\sigma}\mathbf{X}|^2. \qquad (2.14) $$

in a $(N+3)$-dimensional state space. $N$ is the dimension of $X$, 3 refers to the additional scalar unknowns $T$, $\sigma_{z}$ and $\sigma_{\theta}$ and $\phi^T$ is the nonlinear propagator over a time horizon $T > 0$. If convergence fails, the Newton algorithm can be complemented by a globalisation technique such as the Hookstep [124] or the double dogleg step [125]. For the case displayed in figure 2.15, a different TW solution has been found for almost every dip, and it was verified that the TWs look visually close to the initial condition given to the Newton solver. Closer inspection has revealed that one special TW is recurrent in this search (modulo azimuthal shifts), namely the TW labelled S1.

This investigation validates the sketch below (see figure 2.16) where the dynamics along the edge manifold $\Sigma$ can be described as a hyperbolic saddle, i.e. a networks of connected hyperbolic points guiding the dynamics. In the direction transverse to $\Sigma$ the dynamics is dominated by the two attractors of the system, the laminar and the turbulent state, which are outside $\Sigma$. Asymptotic convergence towards a given finite-amplitude solution belonging to $\Sigma$, i.e. towards an edge state, is possible only if, among the list of states sitting on $\Sigma$, one of them possesses a stable manifold of codimension one. This is consistent
<figure><img src="image_37.png" /><figcaption>FIGURE 2.15 – Residual $r_{min}(t_i)$ from [106]. Labels indicate the dips of $r_{min}$ likely to yield good initial conditions for a subsequent Newton search.</figcaption></figure>

with the sketches in figure 2.9.

The imposition of well-chosen discrete symmetries in the computations can help to identify more exact solutions, because it can reduce the number of unstable eigenvalues. Intuition is unfortunately a poor guide here and one relies mostly on luck. In the case of pipe flow for instance, luck was bypassed by the knowledge of one special case. A TW has been found in a former study [126] with axial wavenumber $\alpha = 1.25$. It possesses the $\mathcal{R}_2$ symmetry, i.e. it is $\pi$-periodic in the azimuthal coordinate rather than $2\pi$-periodic. When its eigenvectors are restricted to the symmetric subspace invariant under $\mathcal{R}_2$, only one eigenvalue is unstable. This TW fullfills the conditions to be an edge state according to figure 2.9a. Indeed it was identified as an asymptotic edge state using the classical bisection method. A surprise arose however : depending on the initial condition chosen for the bisection with $\mathcal{R}_2$ symmetry, either that TW is found, or another (new) solution labelled C3. This is attested to by the recurrence function $r_{min}$ shown in figure 2.17, where the convergence to a TW is clear. The convergence to a well-defined attractor such as C3 does not prevent the transient approach to other states. This is for instance the case with the TW labelled A3, identified using the recurrence function in figure 2.17, and displayed in figure 2.18 together with C3.

Edge tracking was employed in the examples above as a *tool* to reveal new exact solutions. These solutions belong by construction to $\Sigma$ and their identification allows one to construct, stone after stone, a better representation of the topology of $\Sigma$. The main advantage of this approach over e.g. homotopy, is the simplicity of the concept. Moreover, the solutions found are guaranteed to be weakly unstable (else they would not have been found) and have a dynamical role that is already understood. The edge tracking method, if complemented with a good continuation solver, yields more than edge states. Continuation allows one to track edge states (and other more unstable states) down to the saddle-node bifurcation from which they are born, and to trace the upper branch as well. This is shown in figure 2.19 as an example drawn from pipe flow. This whole diagram of new solutions [127] was started from the knowledge of one
<figure><img src="image_38.png" /><figcaption>FIGURE 2.16 – Illustration of the hyperbolic saddle structure of the edge manifold in state space, taken from Ref. [106].</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>FIGURE 2.17 – Residual $r_{min}(t_i)$ from [106]. Edge trajectory restricted to the subspace invariant by the $\mathcal{R}_2$ symmetry. $\alpha = 1.25$.</figcaption></figure>

edge state only with $\mathcal{R}_2$ symmetry. The other solutions were found by changing the index of the azimuthal
<figure><img src="image_40.png" /><figcaption>FIGURE 2.18 – Residual $r_{min}(t_i)$ from [106]. Cross-section of TW solutions A3 (left) and C3 (right). Isovalues of streamwise velocity perturbation superimposed with quivers for the transverse velocity field. $\alpha = 1.25$.</figcaption></figure>

symmetry and by applying numerical arclength continuation.

<figure><img src="image_41.png" /><figcaption>FIGURE 2.19 – Bifurcation diagram $c$ (streamwise phase velocity) vs $Re$ for travelling waves solutions in pipe flow with $\alpha = 1.25$ and different values of the azimuthal symmetry index $m$ [127].</figcaption></figure>

These findings suggest a full method in order to identify whole new families of exact solutions of the Navier-Stokes equations. It requires only a timesteper and the evidence that the dynamics is bistable, (with a laminar point as an attractor and another turbulent attractor)... and some luck too. The method can work as well if the attractor is a chaotic saddle, provided its mean lifetime is long enough with respect to the time to approach the edge state along the edge. Then if the system possesses an exact coherent state as an edge state, the bisection algorithm will asymptotically converge to it from *suitable* initial
conditions. This was the case in Refs. [106] and [128] in pipe flow, Ref. [129] in plane Couette flow, Ref. [130] in plane Poiseuille flow, Ref. [113] and Ref. [131] in ASBL (the list is non exhaustive). Convergence can be accelerated using an efficient rootfinder adapted to high dimension. Compared to the homotopy methods used beforehand, the present method offers the important advantage that the dynamic role of these solutions in state space is known by construction.

Numerical continuation of these solutions in R (or in α ) invariably reveals that the branch of the ori- ginal edge state turns around in a fold at some finite parameter value. It appears then that the edge state is a lower-branch solution and that the other solution is an upper-branch solution, the two originating from the same saddle-node bifurcation. Application of this technique in the pipe flow instance of figures 2.17, conducted for several values of the azimuthal wavenumber m , has led to many new branches of TW solutions, as in fig. 2.19. Together with the related investigations in [45], [132], [133] or [134], more than 20 different families of travelling waves have been found so far in spatially periodic pipe flow; all are exact solutions to the Navier–Stokes equations within numerical accuracy. Together they all contribute to the structure of the state space of turbulent pipe flow [126]. The new information is that the dynamics on the edge is organised around the lower-branch solutions only.

2.2.2 Chaoticity of edge trajectories

The sketch of the topology of the edge manifold in figure 2.16 is very comparable in spirit to the random walk picture suggested for the turbulent attractor/saddle in the previous chapter. It is (at least a posteriori ) no surprise for the case of a turbulent saddle, where the stable manifold of the edge state appears dense in the state space. However this is not a trivial result in the usual case where turbulence is an attractor, as evidenced by the counter-examples when the edge state is a simple relative attractor. The main implication of the dynamics sketched in figure 2.16 is the possibility that the edge dynamics be chaotic. While this has been speculated many times, the chaoticity of generic edge trajectories has never been properly demonstrated. We suggest here to focus on a given flow case in which the edge state is known to be dynamically unsteady and to test its chaoticity by using standard Lyapunov analysis, i.e. by computing its Lyapunov exponents λ$_{1}$ , λ$_{2}$ , ..., λ$_{r}$ . Lyapunov exponents are defined as the infinite time limit of FTLEs defined earlier :

$$
\lambda_i = \lim_{T \to \infty} (\lambda_{t_0}^{t_0+T})_i \qquad (2.15)
$$

independently of t$_{0}$ .

The computation of the FTLEs is achieved using the side-computation of Optimally Time-Dependent (OTD) modes [91]. The flow case chosen for this demanding computation is the minimal flow unit of pPF in a computational domain with L + x × L + z = 200 × 100 [135–137] illustrated in figure 2.20. The mass flux is kept constant and the centerline Reynolds number for the laminar flow is Re$_{cl}$ = 7 , 200 corresponding to Re$_{τ}$ = 120. This computation required at least 80 , 000 outer time units of edge tracking for the convergence of the r = 40 first exponents. The spectrum is shown in figure 2.21(left). There are at least 6 exponents that are strictly positive, which validates the former hypothesis that the edge state is chaotic. More can be deduced from this computation : one exponent dominates the spectrum whereas all other exponents (including the negative ones) form a continuum. It is interesting to compare this spectrum to the Lyapunov spectrum corresponding to the turbulent regime of channel flow, as computed in Ref. [138] in a comparable computational domain and comparable Reynolds number Re$_{τ}$ of 80 (although with a less good numerical resolution). There the exponents formed one continuum without any dominant value. The number of positive exponents was also reported to be 166, a number expected to rise for higher Re$_{τ}$ and better resolution. We now turn our attention to the dimension of the underlying attractor, defined according to the Kaplan-Yorke conjecture [139] by

$$
D_{KY} = j + \frac{1}{\lambda_{j+1}} \sum_{i=1}^{j} \lambda_{k} \qquad (2.16)
$$
<figure><img src="image_42.png" /><figcaption>FIGURE 2.20 – Edge state in minimal flow unit pPf for $Re_{\tau} = 100$. Left : 3D snapshot of isolevels of streamwise velocity perturbation. Right : space-time diagram [137] of $<u>_x$ ($y = -0.5, z, t$) (three copies concatenated together).</figcaption></figure>

where

$$ \sum_{i=1}^{j} \lambda_{k} \geq 0 $$

and

$$ \sum_{i=1}^{j+1} \lambda_{k} < 0, $$

The computation of $D_{KY}$ for the edge state at $Re_{cl} = 7,200$ suggests $D_{KY} \approx 22.5 \pm 0.1$ as can be deduced from the change of sign of $\lambda_1 + ... + \lambda_k$ in figure 2.21 (right). As with all strange attractors that number is non-integer. For the turbulent case, the corresponding Kaplan-Yorke dimension is close to 360 [138], a number likely to increase with a higher computational resolution.

### Conditions for turbulence to be transient

One interesting trend appeared during the investigation of this edge state. As previously mentioned, it is customary to define the computational domain in outer units whereas the original study of Jiménez and Moin [28] originally considered a scaling in inner units, making the box size in outer units Re-dependent. In the investigation of [137] different box sizes were tested in inner units, most of them slightly smaller than in Ref. [28], for instance $L_x^+ = 200$ instead of 300. One of the consequences was that the turbulence turned transient no matter how high the value of $Re_{\tau}$ (as far as 180). This suggests that, unlike speculated in early works, the finite lifetime is not only a trivial consequence of low Reynolds numbers. For periodic domains defined in inner units, the finiteness of the lifetimes becomes apparent when the domain size $L^+ = L \times Re_{\tau}$ itself becomes small enough once measured in outer units. This suggests that the transient character of these turbulent flows is linked, consistently with the $2h \approx 100\delta_{\nu}$ concept [9, 10], to the confinement by the domain size. In other words, turbulence in periodic domains is transient when the dynamics is spatially correlated. This is perfectly consistent with the common observation that, at fixed domain size, turbulence relaminarises only at low enough Re. Since full spatial correlation is the property that justifies the dynamical systems approach at odds with genuine spatio-temporal turbulent flows, one is left wondering how much the transiency phenomenon (at least in MFUs) is a consequence of the periodic domains only. Note a similar dilemma in the case where turbulence takes place in a triply periodic domain too [140]. All these results suggest to investigate the transiency phenomenon in the
<figure><img src="image_43.png" /><figcaption>FIGURE 2.21 – Left : Lyapunov spectrum {$\mu_i$}$_i$ for minimal flow unit pPf for $Re_{cl}$ = 7, 200. Right : sum of Lyapunov exponents $\sum_{i=1}^k \mu_i$ versus k.</figcaption></figure>

absence of finite-size effects.

## 2.3 Minimal states

### 2.3.1 Theoretical concept and limitations

The relative attractor on the edge manifold $\Sigma$ is not the only place of interest in the state space. Geometrically speaking, it is interesting to investigate how close $\Sigma$ gets to the laminar attractor. This supposes the choice of metric, usually related to the energy norm, such that the problem amounts to finding

$$ \operatorname{argmin}\{|\boldsymbol{x}|, \boldsymbol{x} \in \mathcal{B}_T\}, \qquad (2.17) $$

where $\mathcal{B}_T$ is the basin of attraction of the laminar state. In the toy problem suggested by Dauchot and Manneville in figure 2.1, the point M closest to the laminar fixed point is represented, using a green diamond, as the tangency point between $\Sigma$ and the smallest ball centered on the laminar state that touches $\Sigma$ [141]. M differs from the edge state and does not form an invariant set. The associated physical question concerns the critical perturbation energy i) below which no transition to turbulence is possible, and ii) above which transition is possible provided the right initial condition is chosen. A given perturbation is associated to this global minimum of the critical energy landscape. It is called the minimal
seed or equivalently the nonlinear optimal disturbance. It must by construction belong to $\Sigma$ which is the boundary of $B_T$, as a consequence of which it can not induce transition by itself : optimal transition requires in practice an infinitesimal perturbation of the minimal seed. A popular question in the 1990s was to identify how the energy of the minimal seed $E_c$ scales with the Reynolds number. Experimental evidence suggests that the optimal threshold should decrease with $Re$, making transition for identical initial disturbances more likely at higher $Re$. Although simple dimensionless analysis suggests a scaling $E_c = O(Re^{-2})$ independently of the flow type [142], many studies including analytical ones suggested other power-law alternatives of the kind $E_c = O(Re^{-\gamma})$ with $\gamma > 1$ [143, 144]. Solving the optimisation problem 2.17 is a direct way to determine how the critical energy or other quantities depend on $Re$. Although the definition of the norm to be minimised is straightforward, one slightly more ambiguous issue is to define a quantitative scalar indicator of the transition phenomenon itself.

A simplified sketch of the state space is shown in figure 2.22, with an emphasis on the minimal seed $M$. The laminar state is the point $O$, the edge state is represented for simplicity as a saddle point $S$, and the turbulent attractor itself does not need to be represented explicitly. The case where turbulence is not an attractor has not been fully investigated as of now. We emphasize that, as for the whole dynamics picture, the concept of minimal seed is intrinsically linked to the framework of initial value problems. Determining the minimal background disturbance likely to maintain non-laminar flow is another problem. Another possible issue, relevant to experiments where localised disturbances are considered, is to consider disturbances applied for a finite duration. This is again a different problem that would call for different mathematical approaches. We will detail below two numerical approaches to the initial value problem (2.17) and describe the transition path to turbulence starting from the neighbourhood of such initial conditions.

<figure><img src="image_44.png" /><figcaption>FIGURE 2.22 – Sketch of the state space with an underlying energy metrics, where $O$ is the laminar state, $S$ the edge state on the edge manifold $\Sigma$, and $M$ the minimal seed.</figcaption></figure>

### 2.3.2 Optimal oblique transition

During the 1990s and 2000s a large body of the numerical literature in shear flows dealt with the concept of linear optimal perturbations (LOP), those maximising the gain due to the linearised operator [18]. There is no direct theoretical link between LOPs and the minimal seed concept, a fully nonlinear concept. In particular LOPs can be computed in different specific subspaces corresponding to different symmetries, such as different wavenumbers, whereas the minimal seed is a global minimiser in state space.
However it is intuitively plausible that a transition path starting from the lowest possible amplitude would share features in common with the linear mechanisms, those commonly observed for low-amplitude (in theory infinitesimally low) disturbances. It was hence assumed for simplicity in [145] that the perturbation velocity field of the minimal seed $U_{opt}$ could be decomposed as a sum of LOPs $\bar{U}_i$, $i = 1, ..., m$, computed each for a different wavenumber :

$$U_{opt} = \sum_{i=1}^{m} A_i \bar{U}_i. \qquad (2.18)$$

The norm to be minimised is the perturbation kinetic energy

$$||U_{opt}|| = \left( \frac{1}{\mathrm{vol}(V)} \int_V |U_{opt}|^2 dV \right)^{\frac{1}{2}}, \qquad (2.19)$$

where $V$ denotes both the computational domain and $\mathrm{vol}(V)$ its three-dimensional volume. A scalar measure of transition remains to be defined. The main idea here exploits the sketch of fig. 2.22 : a perturbation in the neighbourhood of $\Sigma$ must, if it leads to transition, approach the edge state $S$, which has a well-defined energy $E_S$. Minimising the energy of $U_{opt}$ is hence equivalent to maximising the ratio $E_S/E_c$ i.e. the energy gain. This leads in the space of the $\{A_i\}$ to an minimisation of the energy needed to reach, not the turbulent attractor, but the edge state $S$. This is a much better defined optimisation problem, since the critical energy to reach the edge state is precisely a byproduct of the classical bisection process. For $m$ given, we can indeed define formally the mapping $(A_1, ..., A_m) \to E_c(A_1, ..., A_m)$. In this definition, $E_c$ is the perturbation kinetic energy of the state obtained by rescaling $\sum A_i \bar{U}_i$ using the classical bisection algorithm. A simple $m$-dimensional Newton solver is then used to make $\nabla_i E_c$ vanish, i.e. to find the explicit minima of the mapping $\{A_i\}_i \to E_c$. Note that this is a fully nonlinear approach even if LOPs are involved, and the equations used in the optimisation are the fully nonlinear ones. The finiteness of $m$ ($m=2$ or $3$ in practice) makes it only a low-price reduced approach, yet. Since LOPs are not localised in physical space, a short summation of them shares similar properties. This implies that only non-localised minimal seeds can be sought using this method, therefore the approach is restrained to minimal flow units only.

The optimal condition found in plane Couette flow for $m = 2$ is found to feature the modes $(1, 1)$ and $(1, -1)$ with equal amplitude. This result holds for all computational domains and Reynolds numbers tested. This is exactly the *oblique wave* identified earlier in Ref. [146]. This perturbation used nonlinear interactions to generate a mode $(2, 0)$ as well as $(0, 2)$. The optimal mode with $(k_x, k_z) \sim (0, 2)$ corresponds to a pair of streamwise vortices, known to be associated with a strong linear gain via the lift-up effect. This alone justifies that the oblique wave transition path is an efficient one in terms of energy gain. The present results complement this intuition by proving that this path is in fact *nonlinearly optimal*, at least for $m = 2$. The threshold energy of the oblique waves is shown in figure 2.23 as a function of $Re$. The exponent $\gamma$ resulting from the optimisation is $\gamma = 2$. This result corrects in passing the numerical evaluation $\gamma = 2.5$ from [147], the mismatch being attributed to different numerical resolutions.

The flow field resulting from the optimisation with $m = 3$ -based on the LOPs associated with wavevectors $(1, 1)$, $(1, -1)$ and $(1, 2)$ - is displayed in figure 2.24. Its structure is itself very similar to the oblique wave identified for $m = 2$, and its perturbation energy is only 1% less than $E_c$ found for $m = 2$. The corresponding transition scenario is shown in figure 2.25 as a sequence of three-dimensional plots of the wall-normal velocity field. The initial condition corresponding to $U_{opt}$ magnified in amplitude by a factor $(1 + \tilde{\epsilon})$, with $\tilde{\epsilon}$ smaller than $10^{-4}$. The different phases of the transition process unfold as follows. From $t = 0$ to $t = 10$, the oblique wave pattern gets strongly distorted by the shear. It turns rapidly from an initially upward-tilted shape to a downward-tilted one, as an illustration of the linear Orr mechanism in shear flows [18]. This is associated with the transient growth of the oblique modes. After $t = 10$, the structures elongate in the downstream direction and the disturbance is dominated from $t = 20$ on by streamwise vortices with a small but noticeable streamwise undulation. After $t = 40$, the dynamics slows down dramatically and the spatial structure of the perturbation hardly changes until $t = 220$, as an unstable steady state solution located on the edge is being transiently approached. The duration of
<figure><img src="image_45.png" /><figcaption>FIGURE 2.23 – Threshold energy $E_c$ for the oblique scenario ($m = 2$) as a function of the Reynolds number $Re$, for a computational domain of size $L_x = 4\pi$, $L_z = 2\pi$. The results obtained using low numerical resolution, similar to Ref. [147] (green line), are compared to well-resolved simulations where the number of grid points increases with $Re$ (red line). The latter is well approximated by the function $4Re^{-2}$.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>FIGURE 2.24 – (Color online) Initial condition of minimal energy leading to turbulence for $m = 3$. The state is the superposition of three linearly optimal modes with wavenumber ($k_x, k_z$) = (1, ±1) and (1, 2), where the complex amplitudes are determined by the optimisation procedure. Streamwise velocity in the mid-plane $y = 0$ (left), cross-stream velocity field in a section $x = 0$ (right).</figcaption></figure>

the transient approach increases with decreasing $\hat{\epsilon}$. Near $t \approx 240$, the low-speed streaks start to distort and the velocity field becomes rapidly unsteady. By $t \approx 270$ it has all the qualitative features of the turbulent flow reached at $Re = 400$ from other initial disturbances: elongated streamwise structures, stronger velocity fluctuations, stronger unsteadiness.

### 2.3.3 Variational technique

The former approach is difficultly exploitable with higher values of $m$. A natural alternative is to perform a more classical nonlinear adjoint optimisation, yet in the full state space of dimension $n \gg m$,
<figure><img src="image_47.png" /><figcaption>FIGURE 2.25 – Three-dimensional view of the time evolution of the wall-normal velocity, $v = \pm 4 \times 10^{-3}$, along the trajectory $T_{opt}$ initiated by the minimal perturbation M with amplitude $(1+\hat{\epsilon})$ (see text). From top to bottom and from left to right, $t = 2.5, 7.5, 17.5, 25, 90, 180, 245, 295$.</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>FIGURE 2.26 – Optimised objective function versus intial energy $\epsilon_0$. from Monokrous et al. [148]. The energy is found using classical bisection. Blue : relaminarising cases, green : transitioning cases.The minimal energy threshold corresponds to the red star.</figcaption></figure>

with $n$ of order $10^6$. Adjoint optimisation is an iterative process that requires the differentiation of a scalar Lagrangian function with respect to all the variables of the problem. It turns out that the mapping $E_c$ used in Ref. [145] is not straightforward to differentiate, therefore other indicators of the transition must be chosen. A standard choice used in all works by R.R. Kerswell, S. Cherubini and co-authors relies on the energy gain $G$, making nonlinear optimals genuine nonlinear generalisations of the linear concept of the LOPs. In the work carried out in Refs [148] and [149], the integral viscous dissipation was preferred, defined by :

$$ \mathcal{D} = 1 + \frac{1}{Re} \frac{1}{2L_x L_z} \int_V \frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i} dxdydz, \qquad (2.20) $$

where $u$ still denotes the perturbation velocity field. The objective function maximised by the constrained optimisation algorithm is the total time-averaged dissipation

$$ \mathcal{J} = \frac{1}{T} \int_{0}^{T} \mathcal{D}(t)dt. \qquad (2.21) $$

A trajectory reaching the turbulent state at finite time $T$ is expected to display a value of $\mathcal{J}$ larger than for a laminarising trajectory. The main advantage of $\mathcal{J}$ over the energy gain $G$ is that $\mathcal{J}$ contains an average contribution of positive $D(t)$ at all times, and saturates for large $T$ in a monotonous way. In contrast to $G(t)$, it neither undergoes rapid variations nor pronounced overshoots when approaching the edge. Besides, maximising $\mathcal{J}$ amounts to maximising the area between 0 and $D(t) > 0$ in the interval $[0: T]$. This implies for $\epsilon_0 > E_c$ that the algorithm seeks the trajectory reaching high dissipation in the fastest way. It was verified visually in [150] that the optimisation does not appear to depend on the choice of functional.

The optimisation algorithm relies on two steps : first an energy shell characterised by $E_0>0$ is selected, on which the adjoint optimisation is performed together with the additional constraint $E = \epsilon_0$ (for the definition of the energy used in Ref. [149] $E$ does no feature the factor $1/vol(V)$ any longer). Then $E_0$ is varied (as in the bisection process yet with less accuracy) until the maximal value of $\mathcal{J}$ displays an change of trend, indicating that turbulence is no longer reached for $\epsilon_0$ low enough, as in figure 2.26 from the actual data.
The adjoint optimisation algorithm can now be furtively described. It looks iteratively for stationary points of the Lagrange functional , where optimality is fulfilled with respect to the various design variables. The Lagrangian functional L is defined by

$$
\begin{equation}
\begin{split}
\mathcal{L}(\mathbf{X}) = \mathcal{J}(\tilde{\mathbf{u}}) - & \langle \mathbf{v}, \partial_t \mathbf{u} - \mathbf{N} \mathbf{S}(\mathbf{u}, p) \rangle \\
& - \langle q, \nabla \cdot \mathbf{u} \rangle - \lambda(\epsilon(\tilde{\mathbf{u}}) - \epsilon_0).
\end{split}
\tag{2.22}
\end{equation}
$$

(the product $\ll \cdot, \cdot \gg$ denotes here the spatial scalar product $\langle \cdot, \cdot \rangle$ together with time integration). The direct variables in the optimisation are contained in a large vector $\tilde{\mathbf{u}}$ containing the perturbation velocity field $\mathbf{u}(\mathbf{x}, t)$ and the pressure field $p(\mathbf{x}, t)$, and is considered in a suitable space of functions defined on $V \times [0 : T]$. The dynamical constraints to be satisfied include the governing equations and the energy shell constraint :

i) $\partial_t \mathbf{u} = \mathbf{N} \mathbf{S}(\mathbf{u}, p)$ at all times (the momentum equation)

ii) $\nabla \cdot \mathbf{u} = 0$ at all times (the incompressibility equation),

iii) $E = \epsilon_0$ at $t = 0$ (the energy shell constraint).

Differentiating $\mathcal{L}$ with respect to the adjoint variables $\mathbf{v}$, $q$ and $\lambda$ leads to the above constraints. Differentiating $\mathcal{L}$ leads however to a new set of equations. It includes the adjoint Navier-Stokes equations that are linear in $\mathbf{v}$ and need to be integrated backwards in time, as well as compatibility conditions to close the loop. The algorithm is convergent for a given value of $\epsilon_0$.

<figure><img src="image_49.png" /><figcaption>FIGURE 2.27 – Minimal seed for plane Couette flow at $Re = 2000$. Iso-levels of streamwise perturbation velocity. The laminar base flow profile is sketched for visual guidance.</figcaption></figure>

The minimal seed $M$ obtained in pCf for $Re = 2000$ is displayed in figure 2.27. The checkerboard structure typical of oblique waves can be recognised, again with a non-perfect symmetry. However it is now localised spatially and surrounded by laminar flow. It is indeed intuitive that, since energy is now an extensive quantity, a localised velocity field represents a better candidate than a non-localised one for an "efficient" transition process. The optimal transition path starting from a small neighbourhood of $M$ is now displayed in figure 2.28, again by using a sequence of still visualisations of the perturbation velocity field in a plane $y = -0.5$. The early times prior to $t = 80$ feature evidence for the Orr mechanism but are not analysed here. Between $t = 80$ and $t \approx 138$ a sinuous streaky structure is transiently visited, identified
with the edge state (S in the sketch of figure 2.22). A rapid instability of the streaks manifests itself in the form of intense small-scale wall-normal motions at the tips of the streaks. Their breakdown produces even smaller-scale turbulent fluctuations that form a front aligned with the streamwise direction. The front moves in the spanwise direction until the whole flow becomes and remains turbulent.

### 2.3.4 Dependence on $Re$

The dependence of the threshold energy $E_c$ versus $Re$ is displayed in figure 2.29. The emerging power-law scaling has an exponent of $\gamma \approx 2.6$, larger than the exponent $\gamma = 2$ obtained for oblique waves. It is tempting to interpret this steeper exponent as the same exponent as for oblique waves, with a $Re$-dependent correction associated with a spatial localisation that gets more pronounced as $Re$ is increased. This requires however a deeper study involving most likely larger computational domains. The present results highlight the importance of large computational domains in the presence of spatially localised dynamics. It features in passing evidence that the edge state, which is dynamically visited along the optimal transition path, displays itself spatial localisation. The localisation of edge states forms the subject of the next section.

## 2.4 Spatial localisation of edge states

### 2.4.1 Motivation

There are two main reasons to revisit all the previous theory of edge states in the absence of finite-size effects. The first one is unphysical: it is the limitation inherent to the choice of periodic computational domains. Many of the bifurcations occurring with respect to the computational box size, or the Reynolds number, appear as a by-product of the periodic boundary conditions chosen... even though the governing equations are the true ones from hydrodynamics. The second reason is simply the experimental evidence that, in a given range of Reynolds numbers, the coherent structures of the turbulent regime are indeed found to be spatially localised. For instance pipe flow features localised puffs [151] as the only possible turbulent state when $1800 \le Re \le 2500$ approximately. It makes sense to expect the liminal perturbations leading to puffs to be also spatially localised. An auxiliary motivation is that experimentalists can be tempted to excite the flow in a localised fashion; even if no localised turbulent structure sustains forever in the regime of interest, the transition from localised to delocalised features coherent structures characterised by laminar-turbulent interfaces. All these reasons have pushed researchers in subcritical transition to favour larger computational domains (without necessary changing numerical method) in an effort to reduce the impact of the periodic boundary conditions. As often, this conceptual progress has been made possible by some technical progress, notably the ever-increasing computational power, and notably by the development and spreading of parallelisation techniques (MPI). Performing edge tracking in larger domains is only a matter of computing power. The characteristics of the edge tracking algorithm(s) remain unchanged, and there is no strong reason to consider other observables for the bisection process. However, we stress that choosing an extensive quantity (e.g. perturbation kinetic energy or cross-flow energy) is a wiser choice than an intensive one, at least if domain changes are to be carried out. For the classical bisection described earlier, the choice of the bounds $a_L$ and $a_T$ can reveal practically crucial because of the more costly simulations.

To our knowledge the first localised edge states have been computed in pipe flow and shown at the Newton Institute in September 2008 by three different teams, Marburg together with Barcelona [152] at the same time as the team in Bristol [153]. The discovery was quickly extended to pCf [154–158] and later to other flows. The conclusion from all edge tracking computations so far is deceptively simple: edge states appear localised in each direction provided there is enough space for localisation to proceed, with no exception. Importantly for what follows, the choice for longer computational domains only makes spatial localisation possible, but it does not imply it from first principles: the physical reason why a
<figure><img src="image_50.png" /><figcaption>FIGURE 2.28 – Velocity field in the plane $p_y : y = -0.5$. $u_y : \text{colours}$, $(u_x, u_z) : \text{arrows}$. Iso-contours of $u_y$ normalised by their maximum, where $\max u_y = 0.0055$ (t=80), 0.0117 (t=110), 0.052 (t=138), 0.1023 (t=142), 0.2329 (t=146), 0.3333 (t=160), 0.3273 (t=176), 0.5473 (t=184).</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>FIGURE 2.29 – Energy threshold $E_c$ vs. $Re$ for $Re = 750, 1500, 2000, 3000$. The fit $125Re^{-2.7}$ is compared to the fit $4Re^{-2}$ obtained for the Oblique Wave scenario for the same computational box (data from Duguet, Brandt & Larsson (2010)).</figcaption></figure>

given flow state "prefers" to be localised requires a deeper explanation that is still pending and will be discussed in the conclusions. A few representative examples will be considered below, together with a description of the associated dynamics, in order both of increasing localisation and complex dynamics.

### 2.4.2 One-dimensional localisation

#### Striped edge state in channel flow

We begin with the localised edge state which, although not historically first, hosts the simplest dynamics. It was found and reported recently in three-dimensional computational domains of plane channel flow with extension in one direction only. The computational domain is a paralleloid, however it is not trivially aligned with the usual *x* and *z* directions of the flow (denoted *x'* and *z'* in the pedagogical sketch of figure 2.30). Instead it makes an angle *θ* with the streamwise direction. This "tilted domain" was first suggested by Tuckerman & Barkley [159] in order to reproduce numerically, at a cheaper price, the oblique laminar-turbulent patterns found in experiments of pCf. It was described initially as an extended MFU domain, the extension being along one dimension only [159]. It was later generalised without difficulty to simulations of Taylor–Couette flows and plane Poiseuille flows, both known to host oblique turbulent stripes in well-tuned parameter ranges [160–162]. The classical bisection algorithm has been run using the channelflow code [163] with the total perturbation energy as the extensive observable and a well-chosen pair of bounds. The algorithm converged easily towards a constant energy signal, which is the signature of a TW solution and is displayed in figure 2.31. This TW was converged down to accuracy $r_{min} < 10^{-12}$ using the Newton solver included in the channelflow code.

The spatial structure of the solution, shown in figure 2.32a for $\theta = 45^\circ$ and $Re = 720$ is the first property of interest. There is an active core characterised by vortices and streaks. They are approximately aligned with the physical streamwise direction although a small tilt angle is perceptible. Outside this zone the flow is streak-free and almost laminar except for the presence of a smooth secondary flow reconnecting with the laminar flow far enough from the active core. The robustness of the spatial localisation with respect to an increase in domain size has been tested in figure 2.32b by comparing different
<figure><img src="image_52.png" /><figcaption>FIGURE 2.30 – Sketch of tilted periodic domain for pPf.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>FIGURE 2.31 – Bisection Energy $Re = 720$, domain $= 40h \times 2h \times 10h$.</figcaption></figure>
<figure><img src="image_54.png" /><figcaption>FIGURE 2.32 – a) Edge state (several wavelengths concatenated b) Energy $e(z)$ for an edge state $Re = 720$ for increasing domain sizes. Localisation appears as a robust feature of edge states in large enough domains.</figcaption></figure>

converged states, computed for different values of the long domain direction $L_z$, whereas $L_x$ has been fixed to various values between 2 and 10. Here as in all other shear flows investigated, localisation turns out robust. The tails of the velocity field are exponential for all components, a result also found valid in the corresponding localised turbulent regimes [162]. There is another robust feature that distinguishes however these localised edge states from their turbulent counterparts : localisation of edge states is valid for all values of $Re$ whereas the turbulent states are only localised in a finite range of $Re$, bounded by $Re_{cl}$ approximately 2000 in the case of pPf [161].

Another recent discovery worth a deeper investigation concerns the relevance of Waleffe's SSP to the present exact states. For small computational domains it is an undisputed fact that Waleffe's SSP rules out the dynamics and structures of all exact states [33, 135]. This is justified by the fact that the SSP itself was deduced from the detailed analysis of the turbulent dynamics in MFUs. As spatial localisation is added to the picture, the interaction between streaks, rolls and waves predicted theoretically happens to be violated : as seen in figure 2.32a, the streaks in the active part do not display the streamwise undulations characteristic of the SSP. They are however not streamwise-independent either since both streak tips upstream and downstream feature some light tilting (related to the presence of a large-scale
<figure><img src="image_55.png" /><figcaption>FIGURE 2.33 – Bifurcation of TWs as function of Re for different angles.</figcaption></figure>

flow [130]. It is plausible that the feedback of the streak on the streamwise vorticity, the darkest phase in Waleffe's picture, is modified in the presence of localisation and does not even feature the streak instability phase that precedes it. As a supplementary element for discussion, we can mention that other steady localised solutions with the oblique property have been recently identified in pCf [164]. They are not edge states and emerge from modulational instabilities of MFU states. These other new solutions, unlike the edge states of pPF, do feature streamwise undulations. It is then an open question how many variants of the SSP can exist in the presence of localisation, and what the implications are as far the relevance of these solutions to transition is concerned.

Similar TWs have been converged for other values of the angle $\theta$. The range of values of $\theta$ for which TWs have been found widens with increasing $Re$. Arclength continuation has been applied to the TWs obtained for all angles, starting from initial states obtained for $Re = 720$. The structure of the bifurcation diagram for the TWs is shown in figure 2.33. It is confirmed that for a fixed value of $L_{x}$, all TWs emerge in saddle-node bifurcations at marginally low values of $Re$. The value $Re_{SN}$ of the lowest saddle-node bifurcation is less than 400, well below the value of 660 where turbulence can no longer sustain in such computational domains as well as in experiments [161, 165]. Among the many TW solutions ever identified in MFUs of pPf, these are the ones appearing at the lowest values of $Re$ We understand that, unless new solutions are found at even lower $Re$, the state space of plane Poiseuille restricts to the laminar basin for $Re < Re_{SN}$, and gets populated by more and more different unstable TW solutions as $Re$ increases beyond $Re_{SN}$. These solutions, together with a myriad of RPO states bifurcating from them [130], are understood to form the backbone of a chaotic saddle, with a real attractor forming only beyond $Re \ge 660$.
<figure><img src="image_56.png" /><figcaption>FIGURE 2.34 – Instantaneous meridian cuts of the axial perturbation vorticity field for edge states of pipe flow for $Re=2000,3000$ and $6000$ (from top to bottom). The axial dimension has been scaled by a factor of two and only $25D$ of the total domain is shown.</figcaption></figure>

### Localised edge state in pipe flow

Localised edge state in pipe flow were computed initially in Refs [152] and [153] without any imposed discrete symmetry. The bisection method in each study is strictly similar to the one used in the corresponding non-localised case. Only the streamwise length of the computational domain was increased, from e.g. $5D$ in Refs [105, 106] to $33.51D$ and $50$ to $100D$ in Refs [152] and [153], respectively. The results mimic those from the preceding plane Poiseuille case : whatever the value of $Re$ investigated, the edge state identified is localised axially. For the values of $Re$ in the transitional range, where turbulent puffs are found ($1800 < Re < 2500$ approximately), the edge state ressembles its turbulent counterpart, as shown in figure 2.34(top) for $Re = 2000$. As $Re$ is increased beyond that range, the edge state stays axially localised and starts to localise also within the cross-section. As is most evident for $Re = 6000$, ultimately the edge state is best described solely as a finite-amplitude *near-wall disturbance* without any need to be related to the turbulent puffs of the transitional range. A crucial difference with plane Poiseuille lies however in the temporal dynamics : all signals suggest chaotic dynamics. In a direct analogy with the shorter domains of pipe flow [106], an attempt to stabilise the dynamics by imposing a discrete symmetry $R_3$ was performed in Ref. [153] : the dynamics of the energetic observable stayed erratic as long as the edge state could be tracked. Fortunately, this attempt was ignored by Avila et al. who, unaware of this lack of success, tried to impose a discrete $R_2$ and observed convergence to a localised RPO state [128]. This state is quantitatively very similar to a TW, except for a temporal modulation of very small amplitude, observable at the downstream tip only [11, 166]. Note that a similar state has been also observed in a low-resolution bisection in Ref. [153] and initially interpreted as a TW. This state is displayed at the bottom in figure 2.35a using isolevels of streamwise velocity perturbation, compared to its corresponding upper-branch state and to a non-symmetric turbulent puff. Several studies of this RPO state are worth mentioning at this point. The first one is related to the origin of the localisation property : a path in parameter space was found, from the point of view of bifurcations, that relates the non-localised travelling waves of former studies to the present RPO state [167]. As far as we know no localised TW has been found in pipe flow. The second set of results concerns the emergence of complex behaviour in pipe flow [128, 168]. The scenario for the onset of chaos in pipe flow (all in the context of the imposed $R_2$ symmetry) follows a route slightly different from the usual cascade of bifurcations leading from the base flow to the turbulent state. Here the story unfolds not from the base flow, which is linearly stable, but from the RPO-branch of the edge state. Following down this branch in $Re$ leads classically to a saddle-node bifurcation at $Re = Re_{SN} \approx 1428$. After the fold, increasing $Re$ leads now to a cascade of supercritical bifurcations from the upper-branch state, also found in the non-localised cases of Refs. [169] and [170]. The type of cascade seems to depend on the details of the flow case, e.g. torus breakdown in pipe [128] versus period doubling in [169]. The chaotic attractor born at $Re \approx 1540$ is not the end of the story as in the typical Ruelle–Takens scenarios from textbooks. The chaotic attractor grows in size
<figure><img src="image_57.png" /><figcaption>FIGURE 2.35 – a) from top to bottom : turbulent state (puff), upper-branch and lower-branch RPOs for pipe flow at $Re$= (iso-levels of streamwise velocity perturbation. b) bifurcation diagram for the RPO branch.</figcaption></figure>
until it touches the stable manifold of the edge state, i.e. the edge manifold. At this point located at $Re = Re_{bc} \approx 1545 \pm 1$ a boundary crisis occurs: the attractor turns into a saddle and lifetimes become finite. Close to the boundary crisis the mean lifetimes scale as $(Re - Re_{bc})^{-1}$ and decrease with increasing $Re$, before they grow again for yet higher $Re$ [171]. The rest of the story requires even larger domains and will be detailed in the next chapter.

Summing up the pipe flow case, the edge state is again robustly localised yet with a more complex dynamics. As with smaller computational domains the edge dynamics is still of relatively low dimension, so low that a well-chosen restriction to a symmetry subspace (via the imposition of a discrete symmetry) is enough to stabilise a simple RPO state as the edge state. Again in a perfect analogy with smaller domains, a detailed bifurcation study of the RPO state (technically unfeasible as such starting from a chaotic state) reveals the dynamical origin of the chaotic behaviour of the corresponding turbulent state. Assuming that the localisation of the edge state is understood (e.g. via the amplitude modulation suggested in Ref. [167]), this justifies directly the localisation of the turbulent puffs. Besides, the transient chaos so often celebrated in pipe flow finds its origin in a boundary crisis of the upper branch state. Such a scenario starts again with the knowledge of an edge state and relies on the luck of identifying a simple tractable edge state of TW or RPO type.

### Localised edge states in ASBL

The third flow case considered here is the Asymptotic Suction Boundary Layer (ASBL) flow. As for plane Couette or plane Poiseuille flow, the simulation of ASBL in periodic domains requires the knowledge of two numerical parameters $L_x$ and $L_z$, the streamwise and spanwise wavelengths of the computational domain, respectively. An additional parameter compared to the other confined planar flows is the wall-normal extent $L_y$, although it is expected that the results should not depend on $L_y$ once it is large enough, because the base flow is uniform at a sufficiently far distance from the wall located at $y = 0$. Although this can be considered physically artificial, one-dimension spanwise extension can be enforced by choosing $L_z \gg L_x$, with the values of $L_z$ varied between $4\pi$ and $6\pi$. In a similar vein one could also choose $L_x \gg L_z$ as well as tilted domains as in pPf/pCf, although this has not been considered. The present edge tracking computations have all been performed for fixed parameters $Re = 500$ and $L_z = 40$, while the parameter $L_x$ has been varied in small steps [131, 172].

Classical edge tracking has been performed using the code SIMSON developed at KTH, with the cross-flow energy as scalar observable. As for some parameters of pPf (including those in Refs. [98, 99, 173] the energy signal on the edge displays large peak-to-peak variations typical of *bursting* dynamics. We stress that, unlike initially believed by Toh & Itano in their 2001 paper [97], the bursts are intrinsic parts of the edge state and do not belong to its unstable manifold: they do not represent an instability of the edge state but are part of its dynamics. This implies looser bounds for the bisection and eventually longer simulations. The long periods of these bursts makes Newton convergence technically impossible even when the edge state seems to be an RPO. Space-time diagrams typical of converged solutions are represented in figure 2.36. The first three frames correspond to three different asymptotic dynamics obtained for the same set of parameters from different initial conditions. It represents evidence that edge states need not be uniquely defined. The two first frames correspond to two symmetric counterparts and are a consequence of the equivariance of the system with respect to $z \leftarrow -z$ flips. For the first case, the active core where the streaks are present shifts to the right (i.e. increasing $z$) after each burst in a periodic fashion; for the second case the shifts occur towards the left (i.e. decreasing $z$). The third one displays alternate shifts towards the left and the right, also in a periodic fashion yet with a slightly shorter time period. The fourth case displayed has been found for a different streamwise wavelength $L_x$ of $4\pi$. It features very similar shifts in both direction, but the direction of the shifts appears erratic as in a random walk. The mechanics behind the shifts has been analysed in Ref. [131] and turns out to represent a generalisation of the scenario found in related MFUs [113], where low-speed and high-speed streaks periodically cross and exchange roles, therefore appearing as drifting in the spanwise direction.

It was decided to investigate intermediate values of $L_x$ between $4\pi$ and $6\pi$ to understand the dynamical origin of the chaotic edge states. Since none of the states can be converged using a rootfinder and because
<figure><img src="image_58.png" /><figcaption>FIGURE 2.36 – Space-time diagrams of cross-flow energy for edge states in wide domains of ASBL [172], Re = 500, Lz = 400, Lx = 6π except for the last frame where Lx = 4π. The branches for the Left-shifting and Left-Right-shifting states are represented with red squares and blue circles, respectively. Larger symbols represent stable (filled) and unstable (empty) states, whereas erratic behaviour is denoted with smaller symbols.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>FIGURE 2.37 – Bifurcation diagram for edge states in wide ASBL : cross-flow energy versus $L_x$ [172].</figcaption></figure>

of multistability, several new edge tracking computations had to be performed independently for each new parameter. This computationally expensive series of edge computations revealed unsuspected complexity. This is summarized in the incomplete bifurcation diagram of figure 2.37. In this diagram, coloured circles denote RPOs : blue circles represent the states shifting alternatively left and right, represent red circles the ones shifting in one direction only. All the other states are chaotic. Some of them feature clear transient yet recurrent approaches towards one type of periodic state. It was established in Ref. [172] using first-return maps that chaos emerges in various ways, including type-III intermittency [174]. Such a complex wealth of dynamical scenario makes it clear that the picture with the edge state as a uniquely defined saddle point is naive.

It is interesting to investigate the fate of supraliminal perturbations, i.e. initial perturbations that are close to the edge state and belong to the turbulent basin $B_T$. Formally speaking, this corresponds to investigating the one-dimensional unstable manifold of the edge state(s). In practice the tracking of an unstable manifold in high dimension is a delicate task (see e.g. [75, 175]. Fortunately, there are robust features common to most supraliminar perturbations. When the turbulent state is spatially localised (e.g. in the transitional range) the transition form the edge state to the turbulent state proceeds as a local phenomenon. When the turbulent state is not localised, which corresponds to all parameters beyond the transitional range, the transition proceeds in two stages [152, 153] : an initial rapid phase of local transition is followed by a slower phase of front propagation where turbulence invades laminar flow around it. This can also be seen in figures 2.28 in pCf as well as in figure 2.38 in ASBL. In a few exotic situations, the transition process can involve both localisation on the way towards the edge state and delocalisation on the way from the edge state. This is the case for instance for supraliminal noisy perturbations [153].
<figure><img src="image_60.png" /><figcaption>FIGURE 2.38 – Space-time diagram illustrating the delocalisation during the transition process initiated by the edge state in ASBL.</figcaption></figure>
### 2.4.3 Two-dimensional localisation

#### 2.4.4 plane Couette flow

The investigation of the one-dimensional localisation of edge states might seem natural for pipe flow, whose geometry has an "effective dimension" of one. Its relevance is a priori less convincing for planar flows, as attested by the dependence of the results on the parameter $L_x$. This had lead some researchers to denigrate edge states as "spurious" solutions. This is erroneous since edge states are true solutions of the Navier-Stokes equations. Nevertheless the physical relevance of formerly found solutions to experiments is, because of the of periodic boundary conditions involved in their computation, debatable. A truly robust picture of edge states exists however (in the sense that it does not depend on any input parameter). The generalisation from one-dimensional to two-dimensional localisation of edge states is theoretically straightforward. It is robustly attested in all shear flows provided the two planar dimensions (say x and z) present a sufficient extension with respect to the wall-normal direction y. The real practical challenge is the computational cost of the many individual simulations and the difficulty to reach an asymptotic state. This high cost explains why so few full computations of localised edge states have been performed to date, let alone Newton convergence and/or parametric investigations.

The first instance of such computations can be found in Ref. [154] as well as [157] for plane Couette flow. The largest dimensions considered for the computational domain in Ref. [154] are $L_x = 400$ and $L_z = 100$ in units of the half-gap $h$. These dimensions are large enough to warrant unambiguous localisation of the velocity field in $x$ and $z$ as appreciable in figure 2.39 for $Re = 375$. Halving the box dimensions also yields a state describable as localised, with an active core comparable in terms of dimensions : this suggests that the localisation is, again, a robust feature independently of the box size for large enough domains. The dimensions $L_x = 400$ and $L_z = 100$ can be compared to those of the computational box for the MFU, for instance $L_x = 5.51$ and $L_z = 2.51$ in Ref. [24], for which no localisation is detectable at least for $Re \approx 400$. No physical rule of thumb has emerged so far to distinguish between localised and non-localised cases, but evidence suggest that $L_{x,z} > 50$ is sufficient in practice. From a dynamical point of view the edge state in pCf appears chaotic as judged from global energy signals. Some efforts have been undertaken to provide one with a description of the edge trajectory in terms of finite-amplitude solutions. Although transient approaches to apparent steady states have been noted, no convergence towards an exact solution was established. Only later did J.F. Gibson and E. Brand identify a steady state solution with a localised structure similar to that of the edge state in Ref. [154], however that solution turned out neither to be an edge state nor to even lie on the edge manifold. More work is required to have a perfectly deterministic picture of the dynamics on the edge.

Moving towards higher $Re$ suggests that edge states get narrower in $z$ but also longer in $x$, see figure 2.39 for $Re = 1000$. This suggests that higher $Re$ require even longer numerical domains. At this point it is worth remembering that Vortex Wave interaction theory predicts the emergence as $Re \to \infty$ of a critical layer whose thickness scales as $(\alpha Re)^{-3/2}$ [108]. The emergence of such a singular layer is confirmed in MFUs where the streamwise wavelength $\alpha$ remains $O(1)$, however in the present computations with localisation there are several wavelengths $\alpha$ contributing. Choosing the largest wavelength as $\alpha = 2\pi/L_x \sim Re^{-1}$, the thickness of a possible critical layer is $(\alpha Re)^{-1/3} = O(1)$ and no critical layer needs to be considered [176]. Choosing $\alpha$ representative of the $O(1)$ wavelengths found in the active region makes the "local" occurrence of critical layers possible. This is the case for instance for the localised edge state found in pipe flow at $Re = 6,000$ and $10,000$ (resp. Figs 9 and 10 in Ref. [153]).

#### 2.4.5 ASBL

An analogous edge tracking computation has been considered in ASBL, seen as a model for a boundary layer flow without the technical complications due to spatial development [115]. No additional discrete symmetry was imposed. The present computation made use of a PRACE allocation of several million CPU hours exclusively dedicated to this project. Unlike the previous pCf computation, the bisection was
<figure><img src="image_61.png" /><figcaption>FIGURE 2.39 – Localised edge states in pCf for $Re =375$ and 1000 from [154]</figcaption></figure>
<figure><img src="image_62.png" /><figcaption>FIGURE 2.40 – Snapshot of the localised edge in ASBL for $Re = 500$. Isolevels of streamwise velocity perturbation (red/blue) with isolevels of $λ_2 < 0$ (grey).</figcaption></figure>

pushed as far as possible, with the hope that if a simple relative attractor exists it will not be missed. No simple state was eventually found for the parameters tested ($Re = 500$ notably) therefore the edge trajectory is claimed by default to be chaotic. The geometric parameters are $L_x = 800$ and $L_z = 100$, whereas $L_y$ is kept to 15. It has been verified that the spatial structure of the edge state (shown in figure 2.40) changes little for $Re = 500$ when these dimensions vary.

A two-dimensional space-time diagram $(x, t)$ of the cross-flow energy was constructed in a manner similar to those of figure 2.36, with the explicit aim of performing qualitative comparisons. This diagram, displayed in figure 2.41 happens to look surprisingly close to the chaotic edge state found for $L_x = 4\pi$ in figure 2.36d : the main structure is localised in $z$ over a support of approximately 20 units of $\delta_0^*$. It drifts alternately towards $z > 0$ or $z < 0$ with changes of direction over timescales of $O(10^3 \delta_0^*/U_\infty)$. Each change in the numbers of active streaks correlates well with a burst of the cross-flow energy. Such changes of direction were already found to be typical of ASBL in smaller computational domains [131, 172]. Here they similarly involve the crossing of low- and high-speed streaks preliminary to energy bursts. However these bursts and drifts were found neither in large-domain pCf nor in MFU of pCf. The qualitative conclusion from this global comparison is speculative but optimistic : edge states in large computational domains are invariably localised and are robust coherent structures, however they keep dynamically the same main features as their MFU counterparts. This is true with respect to bursts and drifts as well as for critical layers, local wavelength selection etc... More computations would be needed to confirm this speculation. If confirmed, this implies that the investigation of edge states in spatially correlated systems presents a physical interest despite their lack of robustness with respect to geometric parameters. Instead it can be used as a cheaper *laboratory* to understand the more complex yet robust dynamics of localised edge states for the same flows.

The last crucial property of edge states concerns their relevance to the transition process. It is understood that, by definition, only perturbations to the laminar flow with a very specific amplitude are attracted by the edge state. For trajectories initiated very close to the minimal seed, transient visits to the edge state were detected and even exploited to find new exact solutions [112]. This is not surprising given the saddle nature of the edge state but that property is not expected to hold beyond its linear neighbourhood. However, does the edge state play a role of any importance for generic initial conditions, in the absence of any bisection process ? This has been tested by focussing on a large domain of ASBL where different instances of noisy velocity fields (with different amplitudes) were imposed as initial conditions at $t = 0$ [115]. Some initial conditions lead to an incipient turbulent spot and some lead to laminarisation. However all trajectories, when visualised in the state portrait of figure 2.42, displayed a clear transient approach to the same state space region where the edge state is located, similar to the approach to a simple hyperbolic saddle point. A state portrait based on energetic quantities is not a rigorous projection and these results should only be seen as an indication that the edge state plays the role of an unavoidable toll station in a high-dimensional state space.
<figure><img src="image_63.png" /><figcaption>FIGURE 2.41 – Space-time diagram for the localised edge in ASBL at $R = 500$. The dynamics is qualitatively comparable to figure 2.36d.</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>FIGURE 2.42 – State portrait of ASBL during transition initiated by 20 different noisy initial conditions. [115]. Edge state : red solid line. The axis denote the r.m.s. averages of the wall-normal and spanwise velocity components, normalised by the domain size in order to account for the different computational domains used.</figcaption></figure>
## 2.5 Conclusions

This section has mainly focused on the concept of edge states and its development in the context of parallel shear flows in the last 10 years. Edge states correspond to relative attractors for the dynamics of edge trajectories, the trajectories constrained to lie on the laminar-boundary state space boundary $Σ$ [97, 177]. Edge states are not unique and do not necessarily host simple dynamics, instead they can support low-dimensional temporal chaos [106]. Chaotic edge states appear even as the generic rule rather than as an exception. The computational algorithms used to find them rely usually on the classical bisection based on energy bounds, although new variants based on Lagrangian concepts and/or feedback control can be envisioned. The principal success in the exploitation of edge states lies in their property of spatial localisation independent on the physical parameters. The localisation property emerges in practice as soon as the computational domain is large enough [154]. This makes edge states *physically robust* objects since their existence and dynamics no longer depends on numerical parameters such as the computational domain size or the numerical discretisation. The same robustness concept applies as well to minimal states, the initial perturbations of smallest energy likely to trigger actual transition [178]. In some sense, edge states and minimal states are hence the most relevant instance, if not the only one, of application of dynamical systems theory to shear flows. Other concepts such as e.g. periodic orbits have so far only been identified in MFUs. They are therefore intimately connected to periodic boundary conditions and the property of MFUs to capture only spatially correlated dynamics. There is still hope to construct a fully robust dynamical systems theory for localised turbulent states, whose simulation does not rely on a specific choice of domain size. Turbulent puffs in pipe flow appear as the best candidate for such a task and encouraging results have appeared recently [75]. At the time where the present conclusions started to emerge, there was growing evidence that more knowledge on the dynamics of localised states was required, that the study of transitional shear flows in periodic domains was reaching its limits... but also that the current numerical facilities has reached the stage of development where they could shed some light on some open questions. The next chapter is dedicated to the dynamics of localised structures in transitional shear flows. In particular it explores the recent contribution of large domain numerical simulations to the study of transitional regimes of shear flow turbulence.
3

Spatiotemporal aspects of subcritical
transition
## 3.1 Transition as a phase transition

The previous chapter was devoted to the question of the coexistence of the laminar and the turbulent regimes in state space. The present chapter attacks the question of laminar-turbulent coexistence from the point of view of their coexistence in physical space. These are two different questions that request very different angles of attack, different philosophies as well as different methodologies... and lead to different results. The state space coexistence implies, when dealing with extended systems, the possibility to achieve different regimes at different places but at a single time. These different "phases" are then separated by interfaces, whose competition rules the dynamics of the flow. Moreover, it is a well-known peculiarity of subcritical flows that, in a well chosen range of parameters, the different phases can also coexist in a sustained fashion. There is a clear analogy with the language used in equilibrium thermodynamics when dealing with different phases. The general analogy between the laminar-turbulent problem and the thermodynamics of phase change was probably proposed first by Yves Pomeau in the 1980s [179]. It rapidly developed into a school of thought around people such as P. Manneville, F. Daviaud and many others. As soon as the different phases are well defined, interpreting the transition from laminar to turbulent as a genuine phase transition becomes paramount. The use of a vocabular imported from the world of statistical physics becomes very tempting despite many caveats. For instance, should laminar-turbulent transition correspond to a continuous or a discontinuous phase transition ? And of which order ? Some of the concepts from equilibrium thermodynamics do not export very well towards the field of dissipative non-equilibrium physics. Since the laminar phase is defined without any fluctuations, a fruitful analogy must consider it as an absorbing state, whereas the turbulent state contains its own fluctuations. The concept of edge corresponds to a third regime that has the property of being unstable and not realizable. For a phase transition to make sense, a control parameter must be chosen, in general taken to be the Reynolds number $Re$ under one of its numerous incarnations. The question of which order parameter to use is more arduous. There is no free energy associated with the Navier-Stokes equations to rely on to define the type of transition in the spirit of Ehrenfest. Free energies can also be defined as limiting distributions by using the principle of large deviations but this has not yet been exploited deeply in the turbulence context, because of limitations in terms of volume of data. Instead it was chosen to use as order parameter an intensive quantity akin to a *density of turbulence*, namely a turbulent fraction that can be estimated from data (also called intermittency factor). An important question that has stirred the community in the last ten years concerns the universal and the continuous/discontinuous character of the transition using this formalism. The present chapter mentions these issues, but not only. It does not only focus on the onset of the turbulent regime but on the full transitional range, where the turbulent fraction is not necesarily small, but simply finite. Moreover, D. Barkley developed recently a low-order PDE model for pipe flow. The way he derived the model from phenomenological facts [180] is based on a constructive criticism of the simple one-component phase transition picture. This criticism invokes the need for large-scale flows to be included in the picture. The present chapter is an opportunity to investigate large-scale flows, to re-define the concept and to use it for a better understanding of morphogenesis of laminar-turbulent patterns. This chapter begins with the spatiotemporal description of the transitional regimes of pipe flow, followed by planar flows. Spatially developing boundary layer flows are treated only in the next chapter.

## 3.2 One-dimensional shear flows

We begin with the description of pipe flow, whose effective dimension is one. This makes the phenomenology geometrically simpler to describe. The comparison with laboratory experiments. is also a priori simpler. By default the forcing protocol corresponds to constant mass flux. We choose hence the bulk-based Reynolds number $Re = U_b D/\nu$ as the governing parameter. Laminar pipe flow is linearly stable for all Re, hence the laminar regime is always a possibility. A relevant bifurcation diagram for pipe flow should, for simplicity, refer only to regimes observable in the so-called thermodynamic limit of infinite domains over infinite observation times. It should feature laminar flow only when this is the only possible regime, otherwise it should indicate the turbulent regime. A sketchy diagram (taken from [181])