# Does rare, noise-induced, bypass transition in plane Couette flow bypass
instantons ?

Joran Rolland  Univ. Lille, CNRS, ONERA, Arts et Métiers Institute of
Technology, Centrale Lille, UMR 9014 - LMFL - Laboratoire de Mécanique des
Fluides de Lille - Kampé de Fériet, F-59000 Lille, France.
<EMAIL_ADDRESS>
###### Abstract

This text presents the study of rare noise induced transitions from stable
laminar flow to transitional turbulence in plane Couette flow, that we will
term _build up_. We wish to study forced paths that go all the way from
laminar to turbulent flow and to focus the investigation on whether these
paths share the properties of noise induced transitions in simpler systems.
The forcing noise has a red spectrum without any component in the natural,
large scale, linear receptivity range of the flow. As we decreased the forcing
energy injection rate, the transitions became rare. The rare paths from
laminar to turbulent flow are computed using Adaptive Multilevel Splitting
(AMS), a rare event simulation method, and are validated against Direct
Numerical Simulations (DNS) at moderately small energy injection rate. On the
computed trajectories, the flow manages to non-linearly redistribute energy
from the small forced scales to the unforced large scales so that the reactive
trajectories display forced streamwise velocity tubes at the natural scale of
velocity streaks. As the trajectory proceeds, these tubes gradually grow in
amplitude until they cross the separatrix between laminar and turbulent flow.
Streamwise vortices manifest themselves only after velocity tubes have reached
near turbulent amplitude, displaying a two stage process reminiscent of the
“backward” path from turbulence to laminar flow. We checked that these were
not time reversed turbulence collapse paths. As the domain size is increased
from a Minimal Flow Unit (MFU) type flow at $L_{x}\times L_{z}=6\times 4$ (in
half gap units) to a large domain $L_{x}\times L_{z}=36\times 24$, spatial
localisation then extension of the generated coherent streaks and vortices in
the spanwise direction is observed in the reactive paths. The paths
systematically computed in MFU display many of the characteristics of
instantons, that often structure noise induced transitions: such as
concentration of trajectories, exponentially increasing waiting times before
transition, and Gumbel distribution of trajectory durations. However,
bisections started from successive states on the reactive trajectories
indicate that for all sizes and energy injection rates investigated, the
trajectory lack two key ingredients of instantons. Firstly, they do not visit
the neighbourhood of the nearest saddle point and do not display the natural
relaxation path from that saddle to transitional wall turbulence. This
discrepancy is observed for all system sizes. Secondly, the reactive paths do
not concentrate more and more around the same trajectory as energy injection
rate is decreased, but instead gradually move in phase space. They might
reconnect with instantons at very small energy injection rate and exceedingly
long waiting times. They would explain why classical instanton calculations
have proved to be tremendously difficult in wall flows.

##  1 Introduction

Many physical systems, after they underwent a subcritical bifurcation or a
first order transition (or a similar transition in more complex systems),
display two or more metastable states. This situation means that the physical
system can switch from one metastable state to another if it undergoes some
statistically steady forcing, usually by crossing a boundary separating the
basin of attraction of one state from another. Conversely, this physical
system will relax toward either states if it is prepared with an initial
condition on the corresponding side of the separatrix. Traditional bistable
systems are found in kinetic chemistry [Onsager, 1938, Hänggi et al., 1990,
Van Erp et al., 2003, Lopes and Lelièvre, 2019], where the two states are the
reactants and the products or different configurations or conformations [Lopes
and Lelièvre, 2019] of a complex molecule, while the force driving transitions
is the thermal noise. Much insight on bistability has been gained from formal
noise driven systems [Metzner et al., 2006, Cérou et al., 2011, Bouchet and
Touchette, 2012, Rolland and Simonnet, 2015, Rolland et al., 2016, Lucente et
al., 2022]. More recently, a lot of effort has also been directed toward the
study of bistable climate systems [Lucarini and Bódai, 2017, Lucarini and
Bódai, 2019, Baars et al., 2021, Herbert et al., 2020] and bistable turbulent
flows [Podvin and Sergent, 2017] of aerodynamical [Grandemange et al., 2013,
Kim and Durbin, 1988] or geophysical interest [Berhanu et al., 2007, Bouchet
et al., 2019]. In these systems, the metastable states are usually large scale
coherent circulations, like jets, or vortices, while the force driving
transition is intrinsic to the system and related to the fluctuating small
scales of turbulence.

A central object of study for metastable systems is the succession of states
visited by the system when it transits from one bistable state to another. In
kinetic chemistry, this part of the dynamics is termed a reactive trajectory
[Hänggi et al., 1990] or a transition path [Van Erp et al., 2003]. In kinetic
chemistry, and more generally, in a wide class of stochastic systems, as the
thermal noise is decreased and the transition becomes rare, the transition
paths tend to concentrate around a specific path termed the instanton. The
instanton generally consists in a fluctuation path under noise toward a
specific transition state which is a saddle of the deterministic part of the
dynamics, followed by a deterministic relaxation path [Touchette, 2009] (Fig.
15 (a)). This path is important because it highlights the physical mechanisms
of the transition, while the properties of the transition state control the
parametric dependence of the transition rate [Hänggi et al., 1990, Touchette,
2009, Bouchet and Reygner, 2016]. While the instanton always visits the
transition state, complex trajectories can be found in non-gradient systems
[Wan et al., 2015]. The rarity and the concentration of trajectories is
quantitatively studied in the case where noise variance is a constant times
$\epsilon$, a vanishing parameter. One then finds that the transition rate, or
its inverse the mean first passage time before a transition occurs $T$, varies
exponentially with $\epsilon$, as

$\lim_{\epsilon\rightarrow
0}\epsilon\log(T)=\mathcal{I}(\cancel{\epsilon})\,,\,T\underset{\epsilon\rightarrow
0}{\asymp}\exp\left(\frac{\mathcal{I}}{\epsilon}\right)\,,$ (1)

where $\mathcal{I}$ is termed a rate function: it is finite (the limits are
defined), not everywhere zero and is independent on $\epsilon$. This rate
function is directly controlled by the properties of the escaped multistable
state and the saddle. In kinetic chemistry, for relatively low temperatures,
this relation was identified as Arrhenius law. It was then formalised in the
context of stochastic gradient systems as Eyring–Kramers theory [Hänggi et
al., 1990]. Such a formula can also be obtained in non gradient systems
[Bouchet and Reygner, 2016]. In general, the asymptotic study of exponential
dependences in probabilities or rate of probabilities (Eq. (1)) is termed
Large Deviations, and encompasses the case of rare reactive trajectories. In
the case of rare paths, this Large Deviation principle is named after Freidlin
and Wentzell [Freidlin and Wentzell, 1998]. The difficulty of computing the
transition path directly from the dynamics, either because of its complexity
or because direct simulation times are too long for molecular dynamics, has
led to the development and application of rare event simulation methods for
kinetic chemistry [Van Erp et al., 2003, Lopes and Lelièvre, 2019]. Another
approach for these computations consists in exploiting the fact that the
instanton minimises the so-called Onsager–Machlup action [Onsager and Machlup,
1953] in optimisation procedures [Grafke and Vanden-Eijnden, 2019, Wan, 2013,
Wan et al., 2015, Wan and Yu, 2017, Börner et al., 2023]. Note that the noise
covariance matrix is involved in the action, so that the type of forcing can
influence the fluctuation path and the transition rate. The transition state,
however, is a property of the deterministic part of the dynamics and thus has
a degree of universality.

The success of rare event simulation methods in chemistry has encouraged their
use in fluid dynamics. For that purpose, they have been applied to spatially
extended systems of gradual complexity, starting from one dimensional
stochastically forced systems [Rolland, 2018, Baars et al., 2021], moving to
two dimensional stochastically forced turbulence [Bouchet et al., 2019] and
recently reaching three dimensional deterministic transitional turbulence
[Rolland, 2022, Gomé et al., 2022]. The progress has been milestoned in this
manner because of the size of the leap between stochastic differential
equation and fully developed turbulent flows. This means a drastic increase of
the number of degrees of freedom (DoF) involved, a complete change of the
force driving the metastability, and sometimes a strong increase in complexity
of the transition state between the two metastable states. While the extended
systems studied were firstly chosen for their physical interest, the order of
studies was also motivated by the need to tackle this changes separately from
one another at first.

This was the case of the most recent works, concerning the collapse of
turbulence in plane Couette flow [Rolland, 2022, Gomé et al., 2022]. Plane
Couette flow (PCF) is the flow between two parallel moving walls (Fig. 1,
left). In PCF, like in other wall flows such as plane Poiseuille flow [Wan,
2013, Wan et al., 2015, Wan and Yu, 2017], Couette–Poiseuille flow [Liu et
al., 2021], Hagen–Poiseuille pipe flow [Willis and Kerswell, 2009] or boundary
layers [Spangler and Wells Jr, 1968, Rigas et al., 2021], wall turbulence can
coexist at Reynolds number (or positions in the case of boundary layers) for
which the laminar baseflow is linearly stable (see [Romanov, 1973] for the
case of PCF). As a consequence, transitional turbulence can collapse under its
own fluctuations, while the flow can go from laminar to turbulent under a
forcing [Spangler and Wells Jr, 1968, Rigas et al., 2021] or if it is given an
initial condition “close enough” to turbulence [Schmiegel and Eckhardt, 1997].
The study of the collapse of turbulence using rare event methods gave the
occasion to investigate the physical mechanism of laminar hole opening as well
as the technical solutions to study rare multistability events in a three
dimensional flow not stochastically forced [Rolland, 2022, Gomé et al., 2022].
Conversely, the study of the forcing of laminar flow toward turbulence, under
a noise or another, would give the occasion to investigate the crossing of a
very complex separatrix111We use the term separatrix for lack of a better
word. strictly speaking, an unforced wall flow starting from any initial
condition eventually reaches the laminar state: what changes is the relaxation
time before it does so. This duration dramatically increases as an
hypersurface is crossed [Schmiegel and Eckhardt, 1997]..

As such, this has brought the question of the separatrix between laminar and
turbulent flow at the forefront of the study of transition. The first
descriptions by Waleffe and collaborators of the mechanisms maintaining
transitional wall turbulence was the self sustaining process (SSP) of wall
turbulence [Hamilton et al., 1995, Waleffe, 1997]. This process represented
wall turbulence as a stable limit cycle, which was therefore separated from
the laminar baseflow in phase space by an actual separatrix with a saddle on
it. The SSP describes how streamwise velocity streaks and streamwise vortices
regenerate one another. It can take place provided the Reynolds number is
large enough and the domain size is larger than the so called Minimal Flow
Unit (MFU) [Jiménez and Moin, 1991]. Refinements of that picture indicated
that this process was not periodic nor an attractor (thus leaving the room for
collapse of turbulence) but nevertheless confirmed the idea that there was a
separatrix, on which saddles could be found. The structure of the separatrix,
which is a key matter for the forced path from laminar to turbulent flow, was
often studied using methods of dichotomy or bisections [Nusse and Yorke,
1989]. This has led to the study of edge states, saddles on that separatrix,
in wall flows [Toh and Itano, 2003, Schneider et al., 2007, Schneider et al.,
2008, Schneider et al., 2010, Willis and Kerswell, 2009] as well as in climate
models [Lucarini and Bódai, 2017]. Note that additional unstable limit cycles
and fixed points can be found on the turbulent side of the separatrix [Gibson
et al., 2008]. In wall flows, these unstable fixed points and limit cycles are
often reminiscent of the turbulent flow, in that they display velocity streaks
and streamwise vortices. However, they display less small scale fluctuations
and are often very spatially regular. This regularity can be characterised
using the symmetry groups of the flow [Gibson et al., 2008, Schneider et al.,
2010]. Note finally that, as the separatrix and the unstable states on it are
a property of the flow in the configuration of interest and not of initial
conditions or forcing type, they also represent an universal feature of
transition to turbulence.

In wall flows, the complexity of this boundary is such that the saddle on it
is not the point of the separatrix which is has the smallest amplitude of
departure to the laminar baseflow. Experimental and numerical studies of
initial conditions, generated using various specific rules, of minimal energy
that evolve into turbulence have highlighted that. This leaves the question of
whether the most probable forced path would visit this saddle at very small
but finite forcing noise variance. Alternate investigations have thus been
aimed at identifying a path forced and/or triggered by a finite amplitude
initial condition that would drive the flow from laminar to turbulent. Linear
studies of optimal initial conditions of infinitesimal amplitude leading to
maximal amplification of energy identified lift-up: the extraction of energy
from the laminar base flow from streamwise vortices. The lift-up mechanism was
thus identified in transient growth study that highlight the fact that energy
growth can happen while all linear modes are stable is the linearised
Navier–Stokes equations lead to a non normal operator [Farrell and Ioannou,
1993a, Schmid and Henningson, 1994, Trefethen et al., 1993]. Such an optimal
(in shape) initial condition whose amplitude places it on the laminar side of
the boundary would remain on the laminar side of said boundary. However, for
initial conditions starting on the turbulent side of the boundary, lift-up is
an important part of the self-sustaining process [Waleffe, 1997]. Refinement
of optimisation of initial condition energy amplification technics can be used
to compute the optimal finite amplitude initial condition leading to maximal
non-linear energy growth [Cherubini et al., 2011, Rabin et al., 2012,
Monokrousos et al., 2011]. This initial condition is termed a minimal seed.
Such computations are still performed on the laminar side of the separatrix,
however a slight modification of the seed, in the right direction, can
generate an initial condition which has very little energy but can
nevertheless travel to turbulence extremely efficiently.

The question of whether a random forcing of a wall flow can generate a minimal
seed, a flow organisation whose shape and amplitude is precisely defined,
remains open. Indeed, the study of forced wall flows gives little weight to
effects like the Orr mechanism, involved in the growth of minimal seeds. In
parallel with studies of optimal initial conditions leading to transient
growth, the study of the response of the linearised Navier–Stokes equations to
a forcing with a wide range of spatial spectra, which are white in time,
identified a reacting flow structure [Farrell and Ioannou, 1993b]. The weakly
non linear behaviour of flows under such stochastic forcing can also be
performed [Ducimetière et al., 2022]. The inclusion of the strong non-
linearities may very well strongly change this picture because they introduce
the aforementioned complex separatrix. The experimental study of the response
of boundary layers to harmonic forcing of different spectrum has a long
history and shows a strong dependence on the forcing frequency and amplitude
of the point of appearance of wall turbulence [Spangler and Wells Jr, 1968].
Numerical equivalent of such studies of the response to an harmonic forcing
have only been recently undertaken [Rigas et al., 2021]. Meanwhile, to refocus
on our subject of interest, studying the response to a stochastic forcing of
vanishing variance requires a rare event simulation method. This has been
performed in the case of a model of Hagen–Poiseuille pipe flow [Rolland,
2018]. In that case, we have termed the forcing of the flow model from laminar
to turbulent _build up_. However the separatrix in that model is nowhere near
as complex as the one found in actual Hagen–Poiseuille pipe flow. A direct
instanton computation through action minimisation in plane Poiseuille flow can
be extremely complex, and only instantons linking the laminar baseflow to non
linear Tolmien–Schlichting type waves have been performed so far [Wan, 2013,
Wan et al., 2015].

This motivates our present study. As to the rare event simulation method, we
use a method alternate to action minimisation termed Adaptive Multilevel
Splitting (AMS) [Cérou and Guyader, 2007, Cérou et al., 2011], to compute the
rare transitions from laminar baseflow to turbulence in plane Couette flow, in
domains of increasing sizes at transitional Reynolds numbers. By decreasing
the variance of the forcing, proportional to the energy injection rate, we
will investigate the possible concentration of reactive trajectories and
whether this trajectories are structured by an instanton. Using several noise
spectra, we will test whether forcing directly the most responsive structure
is a _sine qua non_ condition to travel to turbulence, while increasing the
system size will highlight possible spatial localisation of the growing
turbulence. We thus present this study in the following manner. We first
present the forced Navier–Stokes equations (§ 2.1), then the version of
Adaptive Multilevel Splitting that is used to compute the trajectories (§
2.2), including the details of the set up of AMS computations and additional
diagnostics used to study the trajectories such as bisections. We describe the
computed paths in section 3, first visually in relevant examples (§ 3.1), then
statistically, on the whole collection of computed path (§ 3.2). We finally
study the effect of the energy injection rate, on the transition rate (§ 4.1)
and on the trajectories duration (§ 4.2). We will discuss these results in
view of the literature in the conclusion (§ 5).

##  2 Method

###  2.1 Forced Navier–Stokes equations

(6,1.5)(6,1.4) (10.5,1.5)(10.5,1.4) (6,2.85)(5.9,2.85)
0$110h$$32h$$\vec{e}_{z}$$\vec{e}_{x}$$\vec{e}_{z}$$\vec{e}_{x}$$R=R_{\rm
g}\simeq 325$$R\simeq R_{\rm t}\simeq 400$Spatiotemporal intermittency(Banded
structure)Uniform Buffer layer turbulence Figure 1: A picture of transition to
turbulence in plane Couette flow indicating the type of flow obtained as the
Reynolds number is increased. Plane Couette flow is sketched on the left: we
indicate the two walls moving at velocities $\pm U$, the domain size $L_{x}$,
$L_{z}$ and $h$ as well as the coordinate system. For Reynolds numbers larger
than a $R_{\rm g}\simeq 325$, laminar-turbulent coexistence can be sustained,
as indicated by the colour levels of kinetic energy
$\frac{1}{2}\mathbf{u}^{2}$ in the horizontal midplane $y=0$ obtained by means
of direct numerical simulations in a domain of size $L_{x}\times
L_{z}=110\times 32$ at $R=370$ (centre). The coloured frames indicate the
domain sizes considered here (cyan: $L_{x}\times L_{z}=6\times 4$, blue
$L_{x}\times L_{z}=18\times 12$, green $L_{x}\times L_{z}=36\times 24$). As
the Reynolds number is further increased, wall turbulence can invade the whole
domain, as illustrated by the colour levels of kinetic energy
$\frac{1}{2}\mathbf{u}^{2}$ in the horizontal midplane $y=0$ obtained by means
of direct numerical simulations in a domain of size $L_{x}\times
L_{z}=110\times 32$ at $R=370$ (right). The visualisations are performed using
numerical results originally performed for [Rolland, 2015], and the whole
figure was originally created for [Rolland, 2022].

####  2.1.1 Governing equations

We consider plane Couette flow, the flow between two parallel walls separated
by a distance $2h$ moving at velocities $\pm U$ (sketched in figure 1, left).
We use our standard notations [Rolland, 2022] so that $\mathbf{e_{x}}$ is the
streamwise direction, $\mathbf{e_{y}}$ is the wall normal direction and
$\mathbf{e_{z}}$ is the spanwise direction. Velocities are made dimensionless
with $U$, lengths are made dimensionless with $h$ and times are made
dimensionless with $h/U$. The main dimensionless control parameter is the
Reynolds number $R=\frac{hU}{\nu}$ where $\nu$ is the kinematic viscosity. The
flow is also controlled by the dimensionless domain sizes $L_{x}$ and $L_{z}$.

The forced incompressible Navier–Stokes equations for the field u, the
departure to the laminar base flow $y\textbf{e}_{x}$, and $p$ the
dimensionless pressure, read

$\frac{\partial u_{l}}{\partial t}+u_{j}\frac{\partial u_{l}}{\partial
x_{j}}+y\frac{\partial u_{l}}{\partial x}+\delta_{l,x}u_{y}=-\frac{\partial
p}{\partial x_{l}}+\frac{1}{R}\left(\frac{\partial^{2}u_{l}}{\partial
x^{2}}+\frac{\partial^{2}u_{l}}{\partial
y^{2}}+\frac{\partial^{2}u_{l}}{\partial
z^{2}}\right)+\sqrt{\frac{2}{\beta}}f_{l}(\mathbf{x},t)\,,\,\frac{\partial
u_{l}}{\partial x_{l}}=0\,.$ (2)

Using tensorial notations, we denote each component of $\mathbf{u}$ by
$u_{l}$, where subscript $l$ and $j$ stands for the component
$\textbf{e}_{x}$, $\textbf{e}_{y}$ and $\textbf{e}_{z}$ (respectively
$l,j=1,2,3$). The forcing $\mathbf{f}$ has zero ensemble average and has in
plane spatial correlations. we will see that its energy injection rate is
controlled by $\beta$. Before giving further details on $\mathbf{f}$, we need
to present the numerical discretisation that we use. The Navier–Stokes
equations are discretised in space using $N_{x}$ and $N_{z}$ dealiased Fourier
modes in the $x$ and $z$ direction (so that $\frac{3}{2}N_{x}$ and
$\frac{3}{2}N_{z}$ modes are used in total) and $N_{y}$ Chebyshev Polynomials
in the $y$ direction. This discretisation and the time integration are
performed using a code based on channelflow, by J. Gibson [Gibson et al.,
2008]. We will use a fixed time step $\Delta t$ for more control on our
numerical procedure. For each system size, Reynolds number and resolution, the
value of $\Delta t$ is chosen so that the CFL criterion is respected. We will
present results of simulations performed in domains of size $L_{x}\times
L_{z}=36\times 24$, where we used $\frac{3}{2}Nx=192$, $\frac{3}{2}N_{z}=144$,
$N_{y}=27$ and $\Delta t=0.05$, $L_{x}\times L_{z}=18\times 12$, where we used
$\frac{3}{2}Nx=96$, $\frac{3}{2}N_{z}=72$, $N_{y}=27$ and $\Delta t=0.05$.
Finally, in MFU type domains of size $L_{x}\times L_{z}=6\times 4$, we
considered four resolutions:

  * •

$\frac{3}{2}N_{x}=32$, $\frac{3}{2}N_{z}=24$, $N_{y}=27$ and $\Delta t=0.05$,
that we termed standard resolution,

  * •

$\frac{3}{2}N_{x}=64$, $\frac{3}{2}N_{z}=48$, $N_{y}=47$, and $\Delta
t=0.0125$, that we termed high resolution,

  * •

$\frac{3}{2}N_{x}=64$, $\frac{3}{2}N_{z}=48$, $N_{y}=27$ and $\Delta t=0.05$,

  * •

$\frac{3}{2}N_{x}=128$, $\frac{3}{2}N_{z}=96$, $N_{y}=35$ and $\Delta
t=0.025$,

to test numerical convergence with spatial resolution of the simulations. We
will comment on the effect of resolution throughout the text. We present the
domain sizes in figure 1, centre, using concentric frames.

The forcing $\mathbf{f}$ is quantitatively prescribed by the correlations in
space and time of each of its component $f_{l}$

$\langle
f_{l}(\textbf{x},t)f_{l}(\textbf{x}^{\prime},t^{\prime})\rangle=\delta(t-t^{\prime})C_{l}(x-x^{\prime},z-z^{\prime})\delta(y-y^{\prime})\,,\,\hat{C}_{l}(n_{x},n_{z})=\Gamma_{l,n_{x},n_{z}}\,.$
(3)

The force $\mathbf{f}$ is white (as indicated by the delta correlations) in
time and in the wall normal direction. Our numerical simulations use a finite
number $N_{y}$ of Chebyshev modes in the wall normal direction. This means
that with a red spectrum in $x$ and $z$, the force $\mathbf{f}$ injects a
finite amount of energy. The correlation function
$C_{l}(x-x^{\prime},z-z^{\prime})$ in the streamwise and spanwise direction on
component $l$ is defined using its Fourier transform. Said transform is
expressed using a shape factor $\Gamma_{l,n_{x},n_{z}}$, where $n_{x}$ stands
for the streamwise wavenumber, $n_{z}$ stands for the spanwise wavenumber. We
have followed two different options as to the details of the forcing,
presented in details in appendix A. One type of noise is forcing selected
components (either a single components or two or all three components). In
that case we have that different components are decorrelated $\langle
f_{l}f_{j}\rangle=0$ if $l\neq j$. Another type of forcing is divergence free
and forces components $\mathbf{e}_{x}$ and $\mathbf{e}_{z}$: it is the curl of
a vector potential along $\mathbf{e}_{y}$. In that case components
$\mathbf{e}_{x}$ and $\mathbf{e}_{z}$ are forced and a correlation can exists
between these two components of the forcing. Numerically, we chose to add the
forcing after the predictor-corrector stage of time integration, in line with
the formulation of an Itô Process. We checked _a posteriori_ that
incompressibility was enforced up to the same numerical precision as in
unforced simulations. Note that in studies of action minimisation in plane
Poiseuille flow, both an unconstrained noise white in time and space [Wan,
2013] and a noise defined as a Wiener process on the space of divergent free
vector fields [Wan et al., 2015, Wan and Yu, 2017], still white in time, were
considered. Linear [Farrell and Ioannou, 1993b] and weakly non linear
[Ducimetière et al., 2022] studies of flows used noise with little constraint
on spectrum shape and component correlations. In the large majority of cases,
we will not directly force the spatial scales expected in the flow of order
$\mathcal{O}(1)$ or larger: we force below a corresponding cut-off scale
$L_{c}$. The value of this scale will be specified in each forcing case. What
we do is that we let the flow non linearly select its scales from energy
injected solely at moderately small scales (some comparison with forcing all
large scales have been performed). Both types of forcing typically give a flat
spectrum for the forcing at intermediate wave numbers and after a wavenumber
cut-off $N_{r}$, we will use a decaying spectrum for the forcing at larger
wavenumbers. We display an example of shape factor $\Gamma_{x,n_{x},n_{z}}$ in
the divergence free case in figure 2 (a) and $\Gamma_{z,n_{x},n_{z}}$ in
figure 2 (b), we we chose $L_{c}=1.0$. As a comparison, the spectra of $u_{x}$
and $u_{z}$ computed in the midplane $y=0$ for transitional turbulence in a
MFU type domain at the largest resolution are displayed in figure 2 (c) and
figure 2 (d) respectively. We can notice that these two components have their
intense modes in the wavenumber range $\left|\frac{2\pi
n_{x}}{L_{x}}\right|\leq 3$, $\left|\frac{2\pi n_{z}}{L_{z}}\right|\leq 6$,
where the forcing is absent. The forcing starts to act in a wavenumber range
where the natural Fourier modes of the velocity field are at least one order
of magnitude smaller than their maximum.

Figure 2:  Colour levels in the wavenumber $2\pi n_{x}/L_{x},2\pi n_{z}/L_{z}$
plane of shape factors of divergence free forcing used for AMS computations in
domains of size $L_{x}\times L_{z}=6\times 4$ at resolution $N_{x}=64$,
$N_{y}=48$, $N_{z}=48$ at Reynolds number $R=500$: (a) logarithm of the
product $((2\pi n_{z}/L_{z})\gamma_{n_{x}}\gamma_{n_{z}})^{2}$ for forcing
applied on $u_{x}$, (b) logarithm of the product $((2\pi
n_{x}/L_{x})\gamma_{n_{x}}\gamma_{n_{z}})^{2}$ for forcing applied on $u_{z}$
(the details of $\gamma_{x}$ and $\gamma_{z}$ are given in appendix A). See
appendix A for full detail on how this generates the vector potential used to
exert a divergence free forcing on the flow. From the velocity field obtained
in a DNS of a turbulent flow in a domain of $L_{x}\times L_{z}=6\times 4$ at
resolution $N_{x}=64$, $N_{y}=48$, $N_{z}=48$ at Reynolds number $R=500$, (c)
logarithm of the modulus of $x-z$ spectrum of $u_{x}$ in the midplane $y=0$,
(d) logarithm of the modulus of $x-z$ spectrum of $u_{z}$ in the midplane
$y=0$.

####  2.1.2 Diagnostics of the flow

We define some quantities of interest that are useful to follow the
transition. We define the spatially averaged kinetic energy as

$E_{k}=\frac{1}{2L_{x}L_{z}}\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=0}^{L_{z}}\frac{u_{x}^{2}+u_{y}^{2}+u_{z}^{2}}{2}\,{\rm
d}x{\rm d}y{\rm d}z\,.$ (4)

The kinetic energy is used to monitor the global state of the flow. When the
flow is turbulent in the whole domain, we find $E_{k}=\mathcal{O}(0.1)$.
Meanwhile, when the flow is globally laminar, $E_{k}$ is much smaller and
close to zero. We will also distinguish $E_{k,l}$ the kinetic energy contained
in a given component $l=x,y,z$, on the one hand, and $E_{k,y-z}$ the kinetic
energy contained in the spanwise and wall normal components, on the other hand

$E_{k,l}=\frac{1}{2L_{x}L_{z}}\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=0}^{L_{z}}\frac{u_{l}^{2}}{2}\,{\rm
d}x{\rm d}y{\rm
d}z\,,\,E_{k,y-z}=\frac{1}{2L_{x}L_{z}}\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=0}^{L_{z}}\frac{u_{y}^{2}+u_{z}^{2}}{2}\,{\rm
d}x{\rm d}y{\rm d}z\,.$ (5)

The kinetic energy in the streamwise component $E_{k,x}$ will be the most
used. The energy $E_{k,x}$ roughly quantifies the amount of energy contained
in the velocity streaks, while the energy $E_{k,y-z}$ roughly quantifies the
amount of energy contained in the streamwise vortices, which are the two main
coherent structures in transitional wall flows. We finally used a turbulent
fraction $F$ to observe the possible spatial localisation of turbulence. This
indicator is computed in the following manner: we spatially average the
kinetic energy contained in components $u_{y}$ and $u_{z}$ in $\tilde{N_{z}}$
streamwise bands of width $\delta z=L_{z}/\tilde{N_{z}}=2.0$, to define

$e(\tilde{n}_{z})=\frac{1}{2L_{x}\delta
z}\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=\tilde{n}_{z}\delta
z}^{(\tilde{n}_{z}+1)\delta z}(u_{y}^{2}+u_{z}^{2})\,{\rm d}x{\rm d}y{\rm
d}z\,,$ (6)

with $0\leq\tilde{n}_{z}<\tilde{N}_{z}$. We then compute the proportion of
bands where $e(\tilde{n}_{z})>0.0005$, which gives our turbulent fraction $F$.
This definition will be adapted to the spatial localisation we will observe on
the trajectories.

We now introduce a parameter alternative to $\beta$ which is possibly better
suited to follow the transitions: the energy injection rate. For this purpose,
we multiply equation (2) by $u_{l}$, sum over the three components and
simplify using incompressibility. We use Itô’s Lemma to manage this non linear
change of variables. We perform a spatial average and an ensemble average with
respect to the forcing noise $\mathbf{f}$. This yields the average energy
budget

$\frac{\partial\langle E_{k}\rangle}{\partial
t}=\underbrace{-\left\langle\frac{1}{2L_{x}L_{z}}\int u_{x}u_{y}\,{\rm d}x{\rm
d}y{\rm
d}z\right\rangle}_{\text{extraction}}\underbrace{-\left\langle\frac{1}{R2L_{x}L_{z}}\int\frac{\partial
u_{l}}{\partial x_{j}}\frac{\partial u_{l}}{\partial x_{j}}\,{\rm d}x{\rm
d}y{\rm d}z\right\rangle}_{\text{dissipation}}+\epsilon\,,$ (7)

where we have the energy injection rate $\epsilon$

$\epsilon=\frac{1}{\beta}\left(N_{y}\sum_{n_{x}=-\frac{N_{x}}{2}+1}^{\frac{N_{x}}{2}}\sum_{n_{z}=-\frac{N_{z}}{2}+1}^{\frac{N_{z}}{2}}\sum_{l=1}^{3}\Gamma_{l,n_{x},n_{z}}\right)\,,$
(8)

in the simulation. In this energy budget (Eq. (7)), the extraction is
performed by lift-up, in the streamwise direction, thus involving $u_{x}$ and
$u_{y}$ and the wall normal gradient of the laminar baseflow equal to $1$.
Dissipation is performed by viscosity. In the case of an unforced turbulent
plane Couette flow, these two terms compensate one another after time average.
For small energy injection rate $\epsilon$, it is generally expected that the
flow remains close to laminar and that equation (2) can be treated linearly
[Farrell and Ioannou, 1993b] or could be processed using a weakly non linear
approach [Ducimetière et al., 2022]. However, this behaviour is conditioned on
no transition happening and is thus not really represented by the actual
ensemble average in equation (7). We will see in the rest of this text that
even in that case, transitions can occur. This means that starting from the
initial condition $\langle E_{k}\rangle(t=0)=0$ even with a small $\epsilon$,
the solution of equation (7) would display a very slow but steady increase of
$\langle E_{k}\rangle(t>0)$ until the flow becomes turbulent. Let us discuss
how transition would be represented in this energy budget. We would expect
that $\epsilon$ would act against dissipation and could maintain some energy
extraction (that would not be sustainable without energy injection, regardless
of the spectrum shape). It is possible that a large part of the energy
transfers would not be visible in this equation as they would be erased by the
spatial averages. We would also expect that the larger $\epsilon$ is, the more
efficiently transition is forced, similarly to the dependence on amplitude of
perturbations to boundary layers [Spangler and Wells Jr, 1968, Rigas et al.,
2021]. The energy injection rate is controlled by two parameters whose effect
on the build-up can be presented independently. Indeed, from equation (8), we
deduce that $\beta$ and the spectrum shape are the two independent control
parameters of this energy injection rate $\epsilon$. Thus, on the one hand,
decreasing $\beta$ increases $\epsilon$, and on the other hand, increasing the
sum of the $\Gamma_{l,n_{x},n_{z}}$ also increases $\epsilon$: in both cases,
we force the flow more intensely. In our definitions (§ A), these shape
factors $\Gamma$ will be normalised and will decrease at large wavenumbers.
Thus the sum will be larger if that decrease starts at larger wave numbers and
if that decrease is slow. Note that we cannot rule out that there could be
effect of the precise shape of the spectrum. These complex effects would be
hidden in the way energy is transferred from the forced scales to the most
receptive scales. There can also be effects of a direct forcing of scales
versus energy transfer. As a final remark, equation (8) reminds that in order
to have a finite or even vanishing energy injection rate in cases where the
flow is forced with a noise white in space (see [Wan, 2013]), one would need
to take the limit of $\beta\rightarrow 0$ first and
$N_{x},N_{y},N_{z}\rightarrow\infty$ next (or use specific relation between
those parameters). With our spectrum shape, we keep the energy injection rate
finite, similarly to what is done in study of linear response [Farrell and
Ioannou, 1993b] and non linear response [Ducimetière et al., 2022] to
stochastic forcing. This physical argument leads to the same mathematical
condition for well posedness of the Freidlin–Wentzell Large Deviation
principle in two dimensional Navier–Stokes equations as noted by [Wan et al.,
2015] and demonstrated by [Brzeźniak et al., 2015].

Note that a similar energy budget is also commonly used in the study of
transition to turbulence [Kawahara and Kida, 2001, Gibson et al., 2008,
Kreilos and Eckhardt, 2012]. This time, it is obtained for the total spatially
averaged kinetic energy $E_{\rm tot}$. For this matter, we apply Itô’s Lemma
to the Navier–Stokes equations for the total flow
$\mathbf{u}+y\mathbf{e}_{x}$, followed by spatial and ensemble average gives
us

$\displaystyle\frac{dE_{\rm tot}}{dt}=\frac{1}{R}(I-D)+\epsilon\,,\,E_{\rm
tot}=\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=0}^{L_{z}}\frac{(u_{x}+y)^{2}+u_{y}^{2}+u_{z}^{2}}{2}\,{\rm
d}x{\rm d}y{\rm d}z$ $\displaystyle
I=1+\frac{1}{2L_{x}L_{z}}\int_{x=0}^{L_{x}}\int_{z=0}^{L_{z}}\left(\frac{\partial
u_{x}}{\partial y}(y=+1)+\frac{\partial u_{x}}{\partial y}(y=-1)\right)\,{\rm
d}x{\rm d}z\,,$ $\displaystyle
D=\frac{1}{2L_{x}L_{z}}\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=0}^{L_{z}}\left\lVert\nabla\times(\mathbf{u}+y\mathbf{e}_{x})\right\rVert_{2}^{2}\,{\rm
d}x{\rm d}y{\rm d}z\,,$ (9)

where $I$ is the total energy injection at the wall and $D$ is the total
dissipation, written using total vorticity because of incompressibility,
periodic and no slip boundary conditions. As such, both are rescaled by their
value in the laminar flow where $\mathbf{u}=\mathbf{0}$ and thus $I=D=1$. This
pair of variable is convenient to observe transition in wall flows. The line
$I=D$ has a special importance because all invariant states lie on it (since
$\frac{dE_{\rm tot}}{dt}=0$ on them), while limit cycles lie over it (since
the time average over a period of $\frac{dE_{\rm tot}}{dt}$ is $0$ on them).
Meanwhile turbulent flow is close to $I=D=3$ [Kawahara and Kida, 2001, Gibson
et al., 2008].

###  2.2 Adaptive Multilevel Splitting

####  2.2.1 The algorithm

(a) (b)

Figure 3: (a) Sketch of two bistable states $\mathcal{A}$ and $\mathcal{B}$
and the hypersurface $\mathcal{C}$ closely surrounding $\mathcal{A}$. two
realisations of the dynamics are sketched: a single excursion in blue and a
first passage trajectory in black and red. The red part of the first passage
trajectory is the reactive trajectory (figure originally made for [Rolland and
Simonnet, 2015]). (b) Sketch of the principle of AMS, showing to iterations of
the algorithm, with $N=3$ clones, indicating the starting state $\mathcal{A}$
and its neighbourhood, the arrival state and its neighbourhood $\mathcal{B}$,
three trajectories are ordered by their $\max_{t}\Phi$. Trajectory one (dashed
blue line) is suppressed and branched on another trajectory at level
$\max_{t}\Phi_{1}$ and then ran according to its natural dynamics. Trajectory
$2$ is then suppressed and branched on 3 at level $\max_{t}\Phi_{2}$ (figure
originally made for [Simonnet, 2016]).

In order to present the numerical method used to compute rare trajectories,
let us first give a formal phase space description of the rare events we will
study in this text. We sketch the build up trajectories in figure 3 (a) and
consider $\mathcal{A}$ a neighbourhood of the laminar flow in phase space, an
hypersurface $\mathcal{C}$ that closely surrounds $\mathcal{A}$ and
$\mathcal{B}$ a neighbourhood of the turbulent flow in phase space. A
realisation of the dynamics which starts in $\mathcal{A}$ fluctuates around it
(for instance has several excursions out of $\mathcal{C}$, an hypersurface
closely surrounding $\mathcal{A}$) and eventually crosses $\mathcal{C}$ and
reaches $\mathcal{B}$ before coming back to $\mathcal{A}$ is termed a first
passage. Its average duration is termed the mean first passage time $T$. The
last stage of the dynamics is termed a reactive trajectory: this is the part
of the dynamics that starts in $\mathcal{A}$, crosses $\mathcal{C}$ and
reaches $\mathcal{B}$ before $\mathcal{A}$. Quantitative definitions of sets
$\mathcal{A}$, $\mathcal{B}$ and hypersurface $\mathcal{C}$ for our case will
be given in section 2.2.2.

We then present the generic formulation for two of the variants of Adaptive
Multilevel Splitting (AMS) which were used for the computation presented in
this text They were initially proposed by Cérou & Guyader [Cérou and Guyader,
2007]. Both variants of AMS use a reaction coordinate (or observable)
$\phi(\mathbf{u}(t))=\Phi(t)$. The reaction coordinate $\phi$ is a mapping
from $\mathbb{R}^{D}$, where the discretised velocity field $\mathbf{u}$
lives, to $[0,1]$. The reaction coordinate can be viewed as a function of time
$\Phi(t)$ on trajectories. The reaction coordinate measures the position of
the flow relatively to starting set $\mathcal{A}$, where is it $0$ and arrival
set $\mathcal{B}$ where it is $1$. We will define our reaction coordinate
quantitatively in section 2.2.2. All variants of AMS run $N$ clone dynamics of
the system to compute iteratively from $N-N_{c}+1$ up to $N$ reactive
trajectories going from set $\mathcal{A}$ to the hypersurface $\mathcal{C}$,
closely surrounding set $\mathcal{A}$, then from $\mathcal{C}$ to the set
$\mathcal{B}$. Here $1\leq N_{c}<N$ is the number of trajectories that will be
removed at each iteration. We will label the velocity field and pressure field
of each trajectory by $1\leq i\leq N$ as in
$\left(\begin{matrix}\mathbf{u}_{i}\\\ p_{i}\end{matrix}\right)(t)$. The
reaction coordinate as a function of time for the corresponding trajectory is
then denoted by $\Phi_{i}(t)=\phi(\mathbf{u}_{i}(t))$. The algorithm is
sketched in figure 3 (b) and proceeds in the following manner

  * •

There is a first stage of natural dynamics, where each clone dynamics starts
in $\mathcal{A}$ from $\mathbf{u}(t=0)=\mathbf{0}$, $p(t=0)=0$, is ran with an
independent realisation of the noise $\mathbf{f}$ and is let to exit and
reenter $\mathcal{A}$ until it crosses the hypersurface $\mathcal{C}$. Once
this has happened, this dynamics is run until the first time it reaches either
$\mathcal{A}$ or $\mathcal{B}$.

  * •

The algorithm then runs the stages of mutation selection. At each stage, the
clones $1\leq i\leq N$, are ordered by $\max_{t}\Phi_{i}(t)$: the maximum of
the reaction coordinate achieved on the trajectory. The $N_{c}$ clones $1\leq
i^{\prime}\leq N_{c}$ having the smallest values of
$\max_{t}\Phi_{i^{\prime}}(t)$ are removed. In order to keep a constant number
of clones, $N_{c}$ new clones are generated by branching. We draw $N_{c}$
clones $i^{\prime\prime}$ uniformly from the clones $N_{c}+1\leq i\leq N$ such
that $\max_{t}\Phi_{i}(t)\geq\max_{t}\Phi_{N_{c}}(t)$. We denote by
$t_{i^{\prime\prime}}$ the first time such that clone $i^{\prime\prime}$
reaches the reaction coordinate value $\max_{t}\Phi_{i^{\prime}}(t)$:
$t_{i^{\prime\prime}}=\min\left\\{t^{\prime}\geq
0\backslash\Phi_{i^{\prime\prime}}(t^{\prime})>\max_{t}\Phi_{i^{\prime}}(t)\right\\}$.
The new clone $i^{\prime}$ is such that
$\mathbf{u}_{i^{\prime}}(t)=\mathbf{u}_{i^{\prime\prime}}(t)$,
$p_{i^{\prime}}(t)=p_{i^{\prime\prime}}(t)$, for $0\leq t\leq
t_{i^{\prime\prime}}$. For $t>t_{i^{\prime\prime}}$ the clone $i^{\prime}$
then follows its natural dynamics under a new, independent realisation of the
forcing $\mathbf{f}$ until it reaches either $\mathcal{A}$ or $\mathcal{B}$.

The algorithm stops after iteration number $\kappa\in\mathbb{N}$, the first
iteration where at least $N-N_{c}+1$ clones trajectories have reached
$\mathcal{B}$. We denote by $r\in\left[1-\frac{N_{c}-1}{N},1\right]$ the
proportion of clones that have reached $\mathcal{B}$ at the last stage. Note
that $\kappa$ and $r$ are _a priori_ random numbers, different from one AMS
run to the other, whose value depend on the independent noise realisations
drawn during simulation of trajectories. This yields an estimator of the
probability $\alpha$ of reaching $\mathcal{B}$ before $\mathcal{A}$ starting
from the distribution of first exit points on $\mathcal{C}$, and the
corresponding mean first passage time $T$

$\alpha=r\left(1-\frac{N_{c}}{N}\right)^{\kappa}\,,\,\hat{\alpha}=\left\langle\alpha\right\rangle_{o}\,,\,\hat{T}=\left\langle\left(\frac{1}{\alpha}-1\right)\tilde{\tau}+\tau\right\rangle_{o}\,,$
(10)

where $\tau$ is the mean duration of reactive trajectories and $\tilde{\tau}$
is the mean duration of non reactive trajectories, computed in each AMS run
and $\langle\cdot\rangle_{o}$ corresponds to an arithmetic average over $o$
AMS runs [Cérou and Guyader, 2007, Rolland, 2018]. In this text, we both
performed calculations with $N_{c}=1$ (systematically leading to $r=1$) and
$N_{c}>1$. The algorithm is naturally parallelised over the $N_{c}$ suppressed
clones: the speed up using $N_{p}$ cores has been measured (not shown here)
and requires $N_{c}/N_{p}$ to be strictly larger than one to remain close to
linear.

As it was done in other studies [Simonnet et al., 2021, Rolland, 2022], we
save a state at the last step $\kappa$ of the AMS computation. This is the
state where the trajectory realising the largest $\Phi$ among the last non
passing reactive trajectory suppressed is branched. We use the same
terminology as in [Rolland, 2022], we term it the last stage at the last step
and denote it by $\mathbf{u}_{\rm last}$. Using the notations introduced
above, this state is computed using the clone $\mathbf{u}_{i^{\prime\prime}}$
that replaces $\mathbf{u}_{i^{\prime}=N_{c}}$. We then have $\mathbf{u}_{\rm
last}=\mathbf{u}_{i^{\prime\prime}}(t_{i^{\prime\prime}}$). In complex
spatially extended systems, it has been verified that this state is very
similar to the saddle crossed by the trajectory [Simonnet et al., 2021]. In
gradient systems, (with few DoF or with spatial extension [Rolland et al.,
2016]), it can be verified that this gives a finer and finer approximation of
the actual, perfectly defined saddle as the noise variance is decreased (not
shown here). In unforced system, this indicates a turning point for the
trajectory [Rolland, 2022].

####  2.2.2 Reaction coordinate

Figure 4: Average kinetic energy conditioned on no build up having occurred
and corresponding variance as a function of the energy injection rate for
several shape of the forcing spectrum and/or several resolutions when (a)
noise on components is applied, (b) when divergence free noise is applied.

Because the kinetic energy $E_{k}$ is almost $0$ for the laminar state, takes
its largest possible values when the flow is turbulent and grows in between,
we used a reaction coordinate affine in spatially averaged kinetic energy

$\Phi(t)=\frac{E_{k}(t)-E_{\mathcal{A}}}{E_{\mathcal{B}}}\,.$ (11)

This reaction coordinate is parameterised by two specific values of kinetic
energy $E_{\mathcal{A}}$ and $E_{\mathcal{B}}$ which depend indirectly on the
control parameters used in the AMS calculations: the Reynolds number, the
energy injection rate and the shape of the noise spectrum. The value
$E_{\mathcal{A}}$ should be such that $E_{k}(t)\simeq E_{\mathcal{A}}$ and
$\Phi(t)\simeq 0$ when the flow fluctuates near the laminar state under
forcing. The exact value of $E_{\mathcal{A}}$ is computed using time series of
kinetic energy sampled at the control parameters of interest:
$E_{\mathcal{A}}$ is defined as the average of $E_{k}(t)$ conditioned on no
transition occurring. The other parameter $E_{\mathcal{B}}$ is such that the
reaction coordinate is one when the flow has reached the neighbourhood of the
turbulent state. Equation (11) then implies that
$E_{\mathcal{B}}+E_{\mathcal{A}}\simeq E_{k,\text{turb}}$. Since
$E_{\mathcal{A}}$ will often be relatively small, we thus chose
$E_{\mathcal{B}}=0.06$ regardless of control parameters, aside of a few
specified cases where we used $E_{\mathcal{B}}=0.06-E_{\mathcal{A}}$.

We furthermore define the hypersurface $\mathcal{C}$ as

$\mathcal{C}=\left\\{\mathbf{u}\backslash\phi(\mathbf{u})=\frac{\sigma}{E_{\mathcal{B}}}\right\\}\,,$
(12)

where $\sigma$ is the standard deviation of $E_{k}(t)$ conditioned on no
transition occurring, sampled from the same dataset as $E_{\mathcal{A}}$. The
rationale is that the flow leaves the neighbourhood of the laminar state when
$E_{k}\geq E_{\mathcal{A}}+\sigma$.

While some AMS computations used _ad hoc_ parameters sampled on a case by case
basis, for most values of the control parameters $R$ and $\epsilon$, We did
not resample $E_{k}(t)$ to compute $E_{\mathcal{A}}$ and $\sigma$. Instead, we
used the parametric dependence of $E_{\mathcal{A}}$ and $\sigma$ on $\beta$ to
interpolate values. Indeed, we can note that provided $\beta$ is large enough,
or conversely $\epsilon$ small enough (the two being related, as seen in (Eq.
(8))) for build up events to be rare, we find that $E_{\mathcal{A}}$ and
$\sigma$ are proportional to $\frac{1}{\beta}$. Moreover, for the two distinct
families of forcing, we find that $E_{\mathcal{A}}$ and $\sigma$ are
proportional to $\epsilon$ in good approximation (Fig. 4 (a) for forcing on
component and Fig. 4 (b) for divergence free forcing). This parametric
dependence is in line with upper bounds of kinetic energy found for the
linearised flow [Farrell and Ioannou, 1993b]. Note however that the
coefficient of this linear dependence depend on the main class of forcing. The
proportionality coefficient in the divergence free case is larger than the one
is the case of forcing on components by a factor more than 10. We also note
that within a class of forcing we do not get exactly the same proportionality
coefficient between $E_{\mathcal{A}}$ and $\epsilon$. While $E_{\mathcal{A}}$
and $\sigma$ grow with the energy injection rate, this simplified dependence
is an approximation. Owing to the different spectrum shape and the variety of
process bringing energy from the forcing to the receptive scale (direct
forcing and transfer across scales), the actual dependence is actually more
complex, and the relation between $E_{\mathcal{A}}$ and $\epsilon$ is not
universal. We finally note that changing the numerical resolution does not
lead to a change in the way $E_{\mathcal{A}}$ and $\sigma$ depend on
$\epsilon$. Note that in more resolved simulations, we have pushed the cut-off
of the spectrum to larger scales, so that the forcing spectrum is different.
Because of that and because the forcing is white in the wall normal direction,
a higher resolved simulation receives more energy for the same value of
$\beta$, as expected from equation (8). We have also tested this relation with
a larger cut-off size of $L_{c}=6$ (not shown here), which leads to a forcing
that can directly affect the most energetic scales. The same linear dependence
of $E_{\mathcal{A}}$ on $\frac{1}{\epsilon}$ is found, but in that case the
values of $E$ are slightly above those computed when $L_{c}=1$ in the forcing,
possibly because large scales now do not only receive energy from non linear
redistribution.

In order to propose the interpolated values, we adjust the linear dependence
of $E_{\mathcal{A}}$ and $\sigma$ on $\beta$ in logarithmic scale for each set
of parameters setting the forcing spectrum. We fit the constants $d_{E}$,
$o_{E}$, $d_{S}$ and $o_{S}$ from data such that

$\log(E_{\mathcal{A}})=\log(\beta)d_{E}+o_{E}\,,\,\log(\sigma)=\log(\beta)d_{S}+o_{S}\,.$
(13)

For this matter, for $M$ values of $\beta_{a}\in\left[4\cdot
10^{5},10^{10}\right]$, $1\leq a\leq M$, we sample $E_{a}$ and $\sigma_{a}$ in
time series where the flow does not transit to turbulence (this is extremely
rare except maybe for $\beta=4\cdot 10^{5}$). Linear regression gives

$\displaystyle
d_{E}=\frac{M\sum_{a=1}^{M}(\log(E_{a})\log(\beta_{a}))-\left(\sum_{a=1}^{M}\log(E_{a})\right)\left(\sum_{a=1}^{M}\log(\beta_{a})\right)}{M\sum_{a=1}^{M}(\log(\beta_{a}))^{2}-\left(\sum_{a=1}^{M}\log(\beta_{a})\right)^{2}}\,,$
(14) $\displaystyle
o_{E}=\frac{\left(\sum_{a=1}^{M}\log(\beta_{a})\right)\left(\sum_{a=1}^{M}\log(\beta_{a})\log(E_{a})\right)-\left(\sum_{a=1}^{M}\log(E_{a})\right)\left(\sum_{a=1^{M}}\log(\beta_{a})^{2}\right)}{\left(\sum_{a=1}^{M}\log(\beta_{a})\right)^{2}-M\sum_{a=1}^{M}\log(\beta_{a})^{2}}\,,$
(15) $\displaystyle
d_{S}=\frac{M\sum_{a=1}^{M}(\log(\sigma_{a})\log(\beta_{a}))-\left(\sum_{a=1}^{M}\log(\sigma_{a})\right)\left(\sum_{a=1}^{M}\log(\beta_{a})\right)}{M\sum_{a=1}^{M}(\log(\beta_{a}))^{2}-\left(\sum_{a=1}^{M}\log(\beta_{a})\right)^{2}}\,,$
(16) $\displaystyle
o_{S}=\frac{\left(\sum_{a=1}^{M}\log(\beta_{a})\right)\left(\sum_{a=1}^{M}\log(\beta_{a})\log(\sigma_{a})\right)-\left(\sum_{a=1}^{M}\log(\sigma_{a})\right)\left(\sum_{a=1}^{M}\log(\beta_{a})^{2}\right)}{\left(\sum_{a=1}^{M}\log(\beta_{a})\right)^{2}-M\sum_{a=1}^{M}\log(\beta_{a})^{2}}\,.$
(17)

We always have $d_{E}\simeq d_{S}\simeq-1$, owing to the proportionality
relation between $E_{\mathcal{A}}$, $\sigma$ and $\frac{1}{\beta}$ (as
reported in table 1, § B). Using equation (13), we can interpolate/extrapolate
$\log(E_{\mathcal{A}})$ and $\log(\sigma)$ for most values of $\beta$, we then
exponentiate.

We give the value used to define the sets $\mathcal{A}$, $\mathcal{C}$ and
$\mathcal{B}$ for each set of simulations in appendix B. This comprises of
specific cases where we will give $E_{\mathcal{A}}$, $E_{\mathcal{B}}$ and
$\sigma$ and generic cases where we give $d_{E}$, $o_{E}$, $d_{S}$ and
$o_{S}$, which mostly concern simulations performed in the largest systems.

####  2.2.3 Direct Numerical simulations of transitions

Along with the AMS computations, we perform Direct Numerical Simulations (DNS)
of the build up, following a procedure we systematically use [Rolland and
Simonnet, 2015, Lucente et al., 2022, Rolland, 2022]. These DNS consist in
systematically repeating the first stage of AMS computations, without the
additional mutation selection. Using the same definition of sets
$\mathcal{A}$, $\mathcal{B}$ and hypersurface $\mathcal{C}$, we start these
simulations at $\mathbf{u}=\mathbf{0}$, $p=0$ in the whole domain, let them
evolve until they reach $\mathcal{C}$ and then stop them when they reach
either $\mathcal{A}$ or $\mathcal{B}$. The proportion of trajectories that
reach $\mathcal{B}$ before $\mathcal{A}$ provides an unbiased estimate of the
crossing probability. Meanwhile the average of the duration of trajectories
that reach $\mathcal{B}$ before $\mathcal{A}$ provides an unbiased estimate of
the average duration of reactive trajectories. These quantities are used to
validate AMS computation and detect possible problems up up to a given degree
of rarity of the event. The estimates rapidly become much more costly than AMS
computations as $\beta$ is increased. For very rare events other tests
(mentioned in § C) are used to probe the validity of computations.

####  2.2.4 Methodology for bisections

We expect that in the later stages of reactive trajectories, the flow crosses
the separatrix between the laminar and the turbulent basins of attractions. We
will first test this by performing relaxation simulations using states taken
on a reactive trajectory. In practice this consists in using states
$\left(\begin{matrix}\mathbf{u}(t_{l}=l\delta t)\\\ p(t_{l}=l\delta
t)\end{matrix}\right)$ sampled every $\delta t>\Delta t$ along a reactive
trajectory and using it as an initial condition for an unforced simulation
($\beta=0$ in equation (2)) and letting it naturally evolve, either toward the
laminar flow $\mathbf{u}\rightarrow\mathbf{0}$, or toward turbulence. We can
thus first detect a $t_{s}$ such that
$\left(\begin{matrix}\mathbf{u}(t_{s})\\\ p(t_{s})\end{matrix}\right)$ relaxes
toward the laminar flow and $\left(\begin{matrix}\mathbf{u}(t_{s}+\delta t)\\\
p(t_{s}+\delta t)\end{matrix}\right)$ relaxes toward the turbulent flow.

Starting from these two fields, successive shootings of a bisection procedure
are then performed [Toh and Itano, 2003]. During shooting $m\geq 1$, we
initialise $\eta_{0}=\frac{1}{2}$, $\eta_{\text{lam}}=1$ and
$\eta_{\text{turb}}=0$. For $0\leq n\leq 36$, we run unforced simulations for
$t\geq t_{m,0}$ started from

$\left(\begin{matrix}\mathbf{u}_{m,n}\\\
p_{m,n}\end{matrix}\right)(t_{m,0})=\eta_{n}\left(\begin{matrix}\mathbf{u}_{\text{lam},m}\\\
p_{\text{lam},m}\end{matrix}\right)+(1-\eta_{n})\left(\begin{matrix}\mathbf{u}_{\text{turb},m}\\\
p_{\text{turb},m}\end{matrix}\right)\,.$ (18)

During the first shooting $m=1$ we use
$\left(\begin{matrix}\mathbf{u}_{\text{lam},1}\\\
p_{\text{lam},1}\end{matrix}\right)=\left(\begin{matrix}\mathbf{u}(t_{s})\\\
p(t_{s})\end{matrix}\right)$ and
$\left(\begin{matrix}\mathbf{u}_{\text{turb},1}\\\
p_{\text{turb},1}\end{matrix}\right)=\left(\begin{matrix}\mathbf{u}(t_{s}+\delta
t)\\\ p(t_{s}+\delta t)\end{matrix}\right)$. We set
$t_{1,0}=t_{s}+\frac{\delta t}{2}$ for presentation of the time series. These
simulations are run until

  1. i)

The flow laminarises, which is deemed to happen if the kinetic energy is below
the threshold $E_{k}(t)\leq 0.001$. In that case we set
$\eta_{n+1}=\eta_{n}-\frac{1}{2^{n+1}}$ and $\eta_{\text{lam}}=\eta_{n}$.

or

  2. ii)

The flow first reaches turbulence, which is deemed to happen if the kinetic
energy is above the threshold $E_{k}(t)\geq 0.07$. In that case we set
$\eta_{n+1}=\eta_{n}+\frac{1}{2^{n+1}}$ and $\eta_{\text{turb}}=\eta_{n}$.

When the shooting is over at $n=36$, we prepare the initial conditions for the
next shooting. We set

$\displaystyle\left(\begin{matrix}\mathbf{u}_{m,-}\\\
p_{m,-}\end{matrix}\right)(t=t_{m,0})=\eta_{\text{lam}}\left(\begin{matrix}\mathbf{u}_{\text{lam},m}\\\
p_{\text{lam},m}\end{matrix}\right)+(1-\eta_{\text{lam}})\left(\begin{matrix}\mathbf{u}_{\text{turb},m}\\\
p_{\text{turb},m}\end{matrix}\right)\,,$
$\displaystyle\left(\begin{matrix}\mathbf{u}_{m,+}\\\
p_{m,+}\end{matrix}\right)(t=t_{m,0})=\eta_{\text{turb}}\left(\begin{matrix}\mathbf{u}_{\text{lam},m}\\\
p_{\text{lam},m}\end{matrix}\right)+(1-\eta_{\text{turb}})\left(\begin{matrix}\mathbf{u}_{\text{turb},m}\\\
p_{\text{turb},m}\end{matrix}\right)\,.$ (19)

We evolve them in parallel until a separation time $t_{m,\text{sep}}$. Setting
$E_{k,m,-}(t)$ and $E_{k,m,+}(t)$ the spatial averaged kinetic energy of these
two fields at $t$, this separation time is the first time multiple of $\delta
t$ such that
$\frac{\left|E_{k,m,-}(t_{m,\text{sep}})-E_{k,m,+}(t_{m,\text{sep}})\right|}{E_{k,m,-}(t_{m,\text{sep}})}>10^{-2}$.
We then define $\left(\begin{matrix}\mathbf{u}_{\text{turb},m+1}\\\
p_{\text{turb},m+1}\end{matrix}\right)=\left(\begin{matrix}\mathbf{u}_{m,-}\\\
p_{m,-}\end{matrix}\right)(t=t_{m,\text{sep}})$ and
$\left(\begin{matrix}\mathbf{u}_{\text{turb},m+1}\\\
p_{\text{turb},m+1}\end{matrix}\right)=\left(\begin{matrix}\mathbf{u}_{m,+}\\\
p_{m,+}\end{matrix}\right)(t=t_{m,\text{sep}})$. We also set
$t_{m+1,0}=t_{m,\text{sep}}$. We run shootings until we have reached an
unstable fixed point and spent a long enough time near it, or we have reached
an unstable limit cycle and enough periods have elapsed, for a precise
computation of its properties.

For the presentation of results, we will concatenate the successive times
series of $E_{k,+}$ and $E_{k,-}$ as the result of the successive shootings.
We will also present velocity fields $\mathbf{u}_{\text{lam},m}$ or
$\mathbf{u}_{\text{turb},m}$ for $m$ large enough such that they are close
enough to either an unstable fixed point or an unstable limit cycle. In
practice we used $m=8$ for all domain sizes except for the largest at
$L_{x}\times L_{z}=36\times 24$ where we used $m=16$.

On these bisections, we use the same spatial and temporal resolution as in the
AMS and DNS computations generating $\left(\begin{matrix}\mathbf{u}(t_{s})\\\
p(t_{s})\end{matrix}\right)$ and $\left(\begin{matrix}\mathbf{u}(t_{s}+\delta
t)\\\ p(t_{s}+\delta t)\end{matrix}\right)$. We used $\delta t=1$ for all
cases, except in the MFU type system at higher resolution where we used
$\delta t=2$.

##  3 Paths toward turbulence

In that section we describe the paths followed by the flow from the
neighbourhood of the laminar state to wall turbulence under forcing, computed
by means of AMS. In the first subsection 3.1, We will first describe full
velocity fields on selected paths (§ 3.1.1), in the last state at the last
stage (§ 3.1.2) and after bisections (§ 3.1.3). In the second subsection 3.2,
we will place these observations in the statistical description of all the
paths we have computed.

###  3.1 Velocity fields in selected examples of transition paths

We will first visualise selected paths to examine the physical mechanism of
crossing from laminar to turbulent flow under forcing. The visualisations will
be complemented by examination of the last state at the last stages as well as
bisections started from velocity fields taken from the reactive trajectories.
We will investigate whether the trajectories follow the classical instanton
phenomenology, whether they cross a separatrix and in that case, whether they
cross it near an unstable fixed point.

####  3.1.1 Description of fields on paths

We will observe the transition paths on four examples:

  * •

In a domain of size $L_{x}\times L_{x}=6\times 4$ at $R=500$ under divergence
free forcing at high resolution (Fig. 5) using $N=20$, $N_{c}=1$,
$\beta=5.5\cdot 10^{4}$, leading to $\alpha=0.35$ and $T=1.2\cdot 10^{3}$.

  * •

In a domain of size $L_{x}\times L_{x}=6\times 4$ at $R=500$ under a forcing
on component $x$ at our standard resolution (Fig. 6), using $\beta=1.1\cdot
10^{5}$, $N=50$, $N_{c}=1$, leading to $\alpha=2.0\cdot 10^{-3}\pm 4.0\cdot
10^{-4}$ and $T=5.3\cdot 10^{5}\pm 1.8\cdot 10^{5}$.

  * •

In a domain of size $L_{x}\times L_{z}=18\times 12$ (Fig. 7) at $R=500$ and
under a divergence free forcing at standard resolution using $N=20$,
$N_{c}=1$, $\beta=6.25\cdot 10^{4}$, leading to $\alpha=0.32$ and $T=6\cdot
10^{3}$.

  * •

In a domain of size $L_{x}\times L_{z}=36\times 24$ (Fig. 8) at $R=500$ and
under a divergence free forcing at standard resolution using $\beta=3\cdot
10^{4}$, $N=20$ and $N_{c}=4$ leading to $\alpha=0.24$ at $T=3.9\cdot 10^{3}$.

Note that computations at $R=400$ in the MFU type domain and in a domain of
moderate size $L_{x}\times L_{z}=12\times 8$ at $R=500$ at standard resolution
were also performed, yielding similar results (not shown here). Note also that
regimes of parameters leading to much smaller probability of build up $\alpha$
and much larger waiting times before build up will be presented later in this
text (in this section and in § 3.2 and § 4). The same phenomenology was
observed in those cases. The parameter defining the reaction coordinate for
these computations are presented in appendix B.

Our visualisations will be centred on the time series of kinetic energy during
reactive trajectories and relaxation/bisection trajectories constructed from
the reactive trajectories (Fig. 5 (f), Fig. 6 (e), Fig. 7 (g), Fig. 8 (g)). In
all four cases, the time series of $E_{k}(t)$ indicate that after crossing
$\mathcal{C}$, the kinetic energy first grows gradually. It does so by
receiving energy in streamwise velocity tubes that show little trace of
streamwise modulation and have a typical spacing slightly shorter than that of
velocity streaks (Fig. 5 (c), Fig. 6 (b), Fig. 7 (b)), Fig. 8 (b)). The energy
contained in the flow increases rapidly after $t\simeq 400$ in the MFU type
system at higher resolution (Fig. 5 (f)), after $t\simeq 1000$ in the MFU type
system forced along $\mathbf{e}_{x}$ (Fig. 6 (e)), after $t\simeq 150$ in the
system of size $L_{x}\times L_{z}=18\times 12$ (Fig. 7 (g)) and after $t\simeq
200$ in the system of size $L_{x}\times L_{z}=36\times 24$ (Fig. 8 (g)). At
that time, streamwise modulation of streamwise velocity, with a spanwise
localisation in these larger domains, can be can be observed (Fig. 7 (d), Fig.
8 (d)). Concomitantly, streamwise vorticity with a near turbulent amplitude is
also observed where $u_{x}$ is modulated (Fig. 5 (b), Fig. 7 (c), Fig. 8 (c)).
The flow then has received enough energy to start the SSP on its own, even
locally and finally evolve toward transitional turbulence globally. We note
that this two stage process is visible in the $\log(E_{k,x}),\log(E_{k,y-z})$
plane for an example of reactive trajectory computed in a MFU type system
(Fig. 11 (a)) at standard resolution for divergence free forcing using
$\beta=6\cdot 10^{5}$ (leading to $\alpha=2.2\cdot 10^{-25}\pm 0.8\cdot
10^{-25}$): we first observe an increase of $E_{k,x}$, corresponding to the
growth in intensity of velocity streaks, at rather low $E_{k,y-z}$ (negligible
streamwise vortices). The next stage of rapid growth of kinetic energy then
corresponds to the concomitant growth of $E_{k,x}$ and $E_{k,y-z}$. The same
reactive trajectory leaves a less clear trace in the $I,D$ plane (Fig. 11
(b)), where both $I$ and $D$ grow in a correlated manner from $1$. We note
however that $D$ fluctuates more than $I$ and that we have $D>I$ on the
reactive trajectory. Being quadratic, $D$ is much more sensitive to
fluctuations of velocity caused by forcing noise, which can average out in the
linear $I$ (Eq. (9)). We find that $D-I\simeq 0.4\pm 0.1$ for $\beta=6\cdot
10^{6}$, corresponding to $(D-I)/R\simeq 8\cdot 10^{-4}$. Since we have
$\epsilon\geq 10^{-3}$ for our control parameters, the growth of total energy
is indeed driven by energy injection by the forcing, but part of that injected
energy is spent in dissipation of the smaller scales generated by the noise
forcing. We further checked this assertion by computing the average and
standard deviation of $I-D$ on reactive trajectories in the range $3\cdot
10^{5}\leq\beta\leq 6\cdot 10^{5}$, both are affine when displayed in log-log
scale as a function of $\beta$ with a $-1$ slope, indicating that they scale
like $\frac{1}{\beta}$.

In the MFU type domain, the comparison of trajectories obtained at higher
resolution with divergent free forcing (Fig. 5), at standard resolution with
divergent free forcing (no detailed example shown here) and with forcing on
the streamwise component at standard resolution (Fig. 6), indicate that the
scenario observed on the trajectory is the same in all three types of case: we
find the same succession of events, the time series of kinetic energy have the
same shape. This indicates that the mechanism generating build up during the
computed reactive trajectories is independent of the forcing family and that
they can be obtained at standard resolution. Note that reactive trajectories
are random realisations of a given process. The trajectories differ from a
computation to the other in duration, in spatial organisation of the
fluctuations, because distinct realisations of the forcing noise are drawn.
This means that we cannot perform a strict one to one comparison of reactive
trajectories are different resolutions at this point. In order to evaluate
whether the computation is converged, we will perform additional comparisons
in the next sections, between distributions of values, averaged values or
paths (§ 3.2), saddle organisation (§ 3.1.3) as well as the scaling laws
$\alpha$, $T$ and $\tau$ follow (§ 4.2).

In extended domains, at the beginning of the rapid increase of kinetic energy
(which will be identified as the crossing of separatrix between laminar and
turbulent flow in § 3.1.3), one can notice that streamwise modulation of
streamwise velocity can be observed in a narrow band in $z$ ($5\lesssim
z\lesssim 10$ for $L_{x}\times L_{z}=18\times 12$, Fig. 7 (d), $13\lesssim
z\lesssim 18$ for $L_{x}\times L_{z}=36\times 24$ Fig. 8 (d)). Streamwise
vorticity of intensity comparable to what is found when the flow is turbulent
appears in the same band ($L_{x}\times L_{z}=18\times 12$, Fig. 7 (c),
$L_{x}\times L_{z}=36\times 24$ Fig. 8 (c)), while its amplitude is almost ten
times smaller outside of it. The streamwise velocity tubes observed outside
said band correspond to a response to the forcing, but they tend to be
narrower, less intense and showing less streamwise modulation than the streaks
showing some streamwise modulation. The subsequent growth of kinetic energy
corresponds to the extension of the modulation in the $z$ direction, until the
whole flow reaches the turbulent state, visualised in figure (Fig. 5 (a)),
(Fig. 7 (a)), and (Fig. 8 (a)), with colour levels of the streamwise velocity
in the midplane $y=0$ and in (Fig. 6 (a)) with isosurfaces of streamwise
velocity in the simulation box. The final turbulent states comprises of the
expected velocity streaks, as visible in the visualisations, and streamwise
vortices.

(f)(a)(b)(c)(c)(d)(e)

Figure 5: View of a reactive trajectory for a system of size $L_{x}\times
L_{z}=6\times 4$ at Reynolds number $R=500$ computed at high resolution under
a divergence free forcing at $\beta=5.5\cdot 10^{4}$, along with the result of
relaxations and bisections started from this trajectory. (a) Colour levels of
the streamwise velocity in the plane $y=0$ when the flow has reached the
turbulent state. (b) Colour levels of the streamwise vorticity in the plane
$y=0$ at $t=440$ (first time step relaxing to turbulence). (c) Colour levels
of the streamwise velocity field in the plane $y=0$ at $t=440$. (d) Streamwise
average of the streamwise velocity field in the $y-z$ plane at a point on the
unstable limit cycle reached after from the bisections. (e) Colour levels of
the streamwise vorticity field in the plane $y=0$ at a point on the unstable
limit cycle reached after from the bisections. (f) Time series of the kinetic
energy during a reactive trajectory (blue), during two relaxation from states
taken on the reactive trajectory, at $t=439$ (last relaxation to laminar
state, red), $t=440$ (first relaxation to turbulent state, orange) and on the
bisections started from these two successive velocity fields (green).

(e)(a)(b)(c)(d)

Figure 6:  View of a reactive trajectory, relaxations and bisections started
from said trajectory, for a system of size $L_{x}\times L_{z}=6\times 4$ at
Reynolds number $R=500$ under a forcing on component $x$ at $\beta=1.1\cdot
10^{5}$. (a) Isovalues $u_{x}=-0.3$ (blue) and $u_{x}=0.3$ (red) of streamwise
velocity in the whole domain when the flow has reached the turbulent state.
(b) Colour levels of the streamwise average of the streamwise velocity field
in the vertical plane $z-y$ for the state relaxing to laminar flow used in the
bisection. (c) Colour levels of the streamwise velocity field in the $y=0$
plane when bisections have reached the neighbourhood of the unstable limit
cycle. (d) Colour levels of the streamwise vorticity field in the $y=0$ plane
when bisections have reached the neighbourhood of the unstable limit cycle.
(e) Time series of the kinetic energy during a reactive trajectory (blue),
during two relaxation from states taken on the reactive trajectory, at
$t=1048$ (last relaxation to laminar state, red), $t=1050$ (first relaxation
to turbulent state, orange) and on the dichotomies started from these two
successive velocity fields (green).

(g)(a)(b)(c)(d)(e)(f)

Figure 7: View of a reactive trajectory for a system of size $L_{x}\times
L_{z}=18\times 12$ at Reynolds number $R=500$ under a divergence free forcing
at $\beta=6.25\cdot 10^{4}$, along with the relaxations and bisections started
from this trajectory. (a) Colour levels of the streamwise velocity in the
plane $y=0$ when the flow has reached the turbulent state. (b) Right: colour
levels of streamwise velocity in the plane $y=0$ at $t=50$ in the early stage
of the reactive trajectory, left: colour levels of the streamwise average of
the streamwise velocity in a $y-z$ plane at the same time. (c) Colour levels
of the streamwise vorticity in the plane $y=0$ at time $t=182$ (last time step
relaxing to laminar flow) from which bisections are started. (d) Colour levels
of the streamwise velocity in the plane $y=0$ at time $t=182$. (e) Colour
levels of the streamwise velocity in the plane $y=0$ when the bisections have
reached a saddle. (f) Colour levels of the streamwise vorticity in the plane
$y=0$ when the bisections have reached a saddle. (g) Time series of the
kinetic energy during the reactive trajectory (blue), during two relaxations
from states taken on the reactive trajectory, at $t=182$ (last relaxation to
laminar state, red), $t=183$ (first relaxation to turbulent state, orange) and
on the bisections started from these two successive velocity fields (green).

(a)(b)(d)(c)(e)(f)(g)

Figure 8:  View of a reactive trajectory for a system of size $L_{x}\times
L_{z}=36\times 24$ at Reynolds number $R=500$ under a divergence free forcing
at $\beta=3\cdot 10^{4}$, and the bisections started from the trajectory. (a)
Colour levels of streamwise velocity in the plane $y=0$ when the flow has
reached the turbulent state. (b) Right: colour levels of streamwise velocity
in the plane $y=0$ at $t=50$ in the early stage of the reactive trajectory,
left: colour levels of the streamwise average of the streamwise velocity in a
$y-z$ plane at the same time. (c) Colour levels of the streamwise vorticity in
the plane $y=0$ at time $t=217$ at which bisections are started. (d) Colour
levels of the streamwise velocity in the plane $y=0$ at the time $t=217$ at
which bisections are started. (e) Colour levels of the streamwise velocity in
the plane $y=0$ when the bisections have reached a saddle. (f) Colour levels
of the streamwise vorticity in the plane $y=0$ when the bisections have
reached a saddle. (g) Time series of the kinetic energy during a reactive
trajectory (blue), during two relaxation from states taken on the reactive
trajectory, at $t=217$ (last relaxation to laminar state, red), $t=218$ (first
relaxation to turbulent state, orange) and on the dichotomies started from
these two successive velocity fields (green).

####  3.1.2 Last state at the last stage

We add visualisations of the last state at the last stage to our
visualisations of the reactive trajectories, following the definition of this
state in section 2.2. Those states were obtained from AMS computations using
the following parameters

  * •

In a domain of size $L_{x}\times L_{x}=6\times 4$ at $R=500$ under divergence
free forcing at high resolution (Fig. 9 (a,b)) using $N=400$, $N_{c}=80$,
$\beta=6\cdot 10^{5}$, leading to $\alpha=2.2\cdot 10^{-25}\pm 0.8\cdot
10^{-25}$ and $T=6\cdot 10^{27}\pm 2\cdot 10^{27}$.

  * •

In a domain of size $L_{x}\times L_{x}=6\times 4$ at $R=500$ under a forcing
on component $x$ at our standard resolution (Fig. 9 (c,d)), using
$\beta=1.1\cdot 10^{5}$, $N=50$, $N_{c}=1$, leading to $\alpha=2.0\cdot
10^{-3}\pm 4.0\cdot 10^{-4}$ and $T=5.3\cdot 10^{5}\pm 1.8\cdot 10^{5}$.

  * •

In a domain of size $L_{x}\times L_{z}=18\times 12$ (Fig. 9 (e,f)) at $R=500$
and under a divergence free forcing at standard resolution using $N=80$,
$N_{c}=16$, $\beta=8\cdot 10^{4}$, leading to $\alpha=7.7\cdot 10^{-6}$ and
$T=1.8\cdot 10^{7}$.

  * •

In a domain of size $L_{x}\times L_{z}=36\times 24$ (Fig. 9 (g,h)) at $R=500$
and under a divergence free forcing at standard resolution using $\beta=3\cdot
10^{4}$, $N=20$ and $N_{c}=4$ leading to $\alpha=0.24$ at $T=3.9\cdot 10^{3}$.

In the aforementioned figure, we include the streamwise velocity in the $y=0$
plane of this state (Fig. 9 (a,c,e,g)) and the streamwise vorticity in the
$y=0$ plane of this state (Fig. 9 (b,d,f,h)), obtained form the AMS
computations that yielded the reactive trajectories presented in section
3.1.1. This diagnostic indicates that at this turning point of AMS
computations, streamwise elongated streamwise velocity tubes are formed and
display turbulent like intensity with $\max_{\text{space}}|u_{x}|\geq 0.5$,
although they present little trace of streamwise modulation in the smaller,
MFU type, system (Fig. 9 (a,c)). The corresponding streamwise vorticity (Fig.
9 (b,d)) thus has a maximum intensity of $|\omega_{x}|\lesssim 0.5$, which is
a fraction of what can be found when the flow is actually turbulent. Again, we
note that there is no striking qualitative difference in those states obtained
two distinct types of forcing. Streamwise velocity tubes are also formed in
larger systems (Fig. 9 (e,g)). However, we note that in a restricted range of
$z$ (taking into account periodic boundary conditions, $0\lesssim z\lesssim 3$
and $10\lesssim z\lesssim 12$ for $L_{x}\times L_{z}=18\times 12$ Fig. 9 (e),
$15\lesssim z\lesssim 24$ for $L_{x}\times L_{z}=36\times 24$ Fig. 9 (g)),
wall turbulence like streamwise modulation of $u_{x}$ is present. We note that
streamwise vorticity of turbulent like intensity $\omega_{x}=\mathcal{O}(1)$
and with turbulent like scale is now only found in a similar range of $z$
(taking into account periodic boundary conditions, $10\lesssim z\lesssim 12$
and to a lesser extent $0\lesssim z\lesssim 3$ for $L_{x}\times L_{z}=18\times
12$ Fig. 9 (f), $18\lesssim z\lesssim 24$ for $L_{x}\times L_{z}=36\times 24$
Fig. 9 (h)). Vorticity is more intense in the larger system. This further
indicates that while the reactive trajectory corresponds to a response to the
forcing through the development of streamwise velocity tube, the flow only
reach a turning point when locally, the tubes become sufficiently intense and
display intense enough shear between one another to trigger the formation of
streamwise modulation and corresponding streamwise vortices. Once this has
happened, this streamwise modulation and said vortices will contaminate nearby
regions in $z$ until they extend to the whole flow.

Figure 9:  From the AMS computation in the system of size $L_{x}\times
L_{z}=6\times 4$ at $R=500$ and $\beta=5.5\cdot 10^{5}$ at high resolution,
(a) colour levels of the streamwise velocity in the plane $y=0$ on the last
state at the last stage for the AMS computation, (b) colour levels of the
streamwise vorticity in the plane $y=0$ on the last state at the last stage
for this AMS computation. From the AMS computation in the system of size
$L_{x}\times L_{z}=6\times 4$ at $R=500$ and $\beta=6\cdot 10^{5}$ at sandard
resolution, (c) colour levels of the streamwise velocity in the plane $y=0$ on
the last state at the last stage for the AMS computation, (d) colour levels of
the streamwise vorticity in the plane $y=0$ on the last state at the last
stage for this AMS computation. From the AMS computation in the system of size
$L_{x}\times L_{z}=18\times 12$ at $R=500$, and $\beta=8\cdot 10^{4}$, (e)
colour levels of the streamwise velocity in the plane $y=0$ on the last state
at the last stage for the AMS computation, (f) colour levels of the streamwise
vorticity in the plane $y=0$ on the last state at the last stage for this AMS
computation. From the AMS computation in the system of size $L_{x}\times
L_{z}=36\times 24$ at $R=500$, at $\beta=3\cdot 10^{4}$, (g) Colour levels of
the streamwise velocity in the plane $y=0$ on the last state at the last stage
for this AMS computation, (h) Colour levels of the streamwise vorticity in the
plane $y=0$ on the last state at the last stage for this AMS computation.

####  3.1.3 Bisections

Following the procedure presented in section 2.2.4, we perform relaxations and
bisections using reactive trajectories computed by means of AMS. We will
examine the type of unstable attractor (unstable fixed point, unstable limit
cycle, unstable travelling wave) which is reached and describe its main
features. The result of these bisections will be compared to reactive
trajectories and the last state at the last stage examined in the former
subsection. We will check whether the system displays a standard instanton
phenomenology, which would be characterised by greater and greater
correspondence between $\mathbf{u}_{\rm last}$ and the result of bisections:
this would mean the trajectories would cross the separatrix very near the
saddle. The time series of these relaxations and the result of following
bisections are added to the plots:

  * •

of figure 5 (f) for a MFU type system at high resolution with divergence free
noise with $\beta=5.5\cdot 10^{4}$ and a sampling time of $\delta t=1$ on the
trajectory, starting bisections at $t=439$,

  * •

of figure 6 (e) for a MFU system at standard resolution with noise on
component $\mathbf{e}_{x}$ with $\beta=1.1\cdot 10^{5}$, and a sampling time
of $\delta t=2$ on the trajectory, starting bisections at $t=1048$,

  * •

of figure 7 (g) for a system of size $L_{x}\times L_{z}=18\times 12$ with
divergence free noise with $\beta=6.25\cdot 10^{4}$ and a sampling time of
$\delta t=1$ on the trajectory, starting bisections at $t=182$,

  * •

of figure 8 (g) for a system of size $L_{x}\times L_{z}=36\times 24$ with
divergence free noise with $\beta=3\cdot 10^{4}$ and a sampling time of
$\delta t=1$ on the trajectory, starting bisections at $t=217$.

The corresponding values of $\alpha$ and $T$ are given in section 3.1.1. The
first two relaxations as well as the result of bisections is also added to the
plots in the $\log(E_{k,x}),\log(E_{k,y-z})$ plane (Fig. 11 (a)) and in the
$(I,D)$ plane (Fig. 11 (b)). Relaxation toward the laminar flow naturally
corresponds to $\log(E_{k,x}),\log(E_{k,y-z})\rightarrow-\infty$, and
$I,D\rightarrow 1$, while relaxation toward turbulence leads to time
fluctuations, naturally found near $I\simeq D\simeq 3$. Note that during
relaminarisation, we have $I>D$, energy injection is larger than dissipation
leading to $E_{\rm tot}$ is actually growing (Eq. (9)). Indeed, this energy is
dominated by the streamwise flow which happen to be less and less depleted as
the flow goes from turbulent to laminar. As the size of the system is
increased (from $6\times 4$ to $36\times 24$), longer and longer transients
have to be passed until said saddle is reached. During these transients, the
smallest scales visible on snapshots of the reactive trajectories are smoothed
out, in particular in the wall normal and spanwise components on the velocity
field.

For the smallest, MFU type system, we reach an unstable limit cycle that
displays small oscillations of kinetic energy (Fig. 5 (f), Fig. 6 (e)). This
cycle can also be viewed in the $\log(E_{k,x})$, $\log(E_{k,y-z})$ plane (Fig.
11 (a)), this indicates that these two components have a phase shift
($E_{k,x}$ is almost in phase opposition with $E_{k,z}$ and in phase
quadrature with $E_{k,y}$). This limit cycle consists in oscillations of the
streamwise modulation of a pair of streamwise velocity tubes (Fig. 5 (d), Fig.
6 (c)). As expected the low speed tubes are located in the upper part of the
channel $y>0$ and the high speed tubes is located in the lower part of the
channel $y<0$. The velocity in the tubes is rather intense, leading to
$\max_{\text{space}}|u_{x}|=0.75\pm 0.05$. The modulation in the two velocity
tubes are best described as in phase quadrature with one another. This
streamwise modulation leads to a corresponding streamwise vorticity which has
a rather low magnitude $|\omega_{x}|=\mathcal{O}(0.1)$. Note that the spatial
variation of streamwise vorticity is sharp, with an oblique $\omega_{x}=0$
stripe striking each vortex in the $y=0$ plane (Fig. 5 (e), Fig. 6 (d)). This
cycle is found both at lower and higher resolution and regardless of the
forcing type (Fig. 5 (f), Fig. 6 (e), Fig. 11 (a)). This indicates that the
reactive trajectories cross the separatrix in the same region for both
families of forcing, and that the crossing location of the separatrix can be
viewed as converged with spatial resolution. The relative amplitude of
oscillations of the streamwise component is small: at high resolution the time
average of kinetic energy is $\langle E_{k,x}\rangle_{t}=0.0620$, while the
rescaled root mean square is $\sqrt{2}\sqrt{\langle
E_{k,x}^{2}\rangle_{t}-\langle E_{k,x}\rangle_{t}^{2}}=0.0023$. The relative
amplitude of oscillation of the wall normal and spanwise components is much
larger, we find that the time averages of kinetic energy on these components
are $\langle E_{k,y}\rangle_{t}=3.55\cdot 10^{-5}$ and $\langle
E_{k,z}\rangle_{t}=1.44\cdot 10^{-4}$, while the rescaled root mean square are
$\sqrt{2}\sqrt{\langle E_{k,y}^{2}\rangle_{t}-\langle
E_{k,y}\rangle_{t}^{2}}=1.22\cdot 10^{-5}$ and $\sqrt{2}\sqrt{\langle
E_{k,z}^{2}\rangle_{t}-\langle E_{k,z}\rangle_{t}^{2}}=9.14\cdot 10^{-5}$. The
increase of resolution leads to minute changes to this limit cycle: we find a
relative increase of 2% on $\langle E_{k,x}\rangle_{t}$ and 6% on $\langle
E_{k,y}\rangle_{t}$ and $\langle E_{k,z}\rangle_{t}$, with smaller changes on
the amplitude of the cycles, as measured by much smaller change in the RMS
kinetic energy. The visualisation of $\omega_{x}$ at both resolution (Fig. 5
(e), Fig. 6 (d)) confirms that the difference is minute. We also observe the
cycle in the total energy injection $I$, total dissipation $D$ plane (see Eq.
(9)) in figure 11 (b). Note that in line with the observation of the
disappearance of very small scale fluctuations, dissipation rapidly decrease
from higher fluctuating value to a value close to $I$ during the relaxation.
Bisections lead to the limit cycle which lies slightly under to the $I=D$
line. This is an effect of discretisation errors in the evaluation of $I$ and
$D$ combined with the very small amplitude of the limit cycle. As resolution
is increased, this cycle moves toward the line. Finally, we measured the
period of oscillations of kinetic energy by determining the frequency of the
largest amplitude peak discrete Fourier transforms of the concatenated times
series of $E_{k}(t)$ obtained during successive shootings. This gives
$\mathcal{T}=59.8\pm 0.6$ at low resolution (using a time series of duration
$4537$) and $\mathcal{T}=58.2\pm 0.6$ at higher resolution (using a time
series of duration $5249$). This again confirms the convergence of the
numerical simulation. Observation of the velocity fields indicate that every
half period of kinetic energy, the streamwise modulation of velocity fields
nears zero, then grows with a $L_{x}/2$ phase shift. This means that the
actual time period of the velocity fields is twice that of the kinetic energy.
We observe this in figure 10 (a) using the rescaled phase of the first Fourier
mode of $u_{z}$ computed in the midplane at $z=2$ defined as

$\psi(t)=L_{x}\left(\frac{\arg(\hat{u}_{z}(n_{x}=1,y=0,z=2,t))}{2\pi}+\frac{1}{2}\right)\,.$
(20)

This rescaling matches with a $x$ position. It is incremented by $L_{x}$ each
time $\arg(\hat{u}_{z})$ goes from $2\pi$ to $0$ and decremented by $-L_{x}$
each time $\arg(\hat{u}_{z})$ goes from $0$ to $2\pi$. The succession of
plateau separated by a rapid jump every period of the kinetic energy
quantifies the fact that we periodically see the pattern in $\omega_{x}$
disappear and reappear at a distance $L_{x}/2$. This edge state in the form of
an unstable limit cycle has not been specifically studied in the literature to
our knowledge, possibly because of the smallness of the domain. However, a
similar spatial organisation of the flow, similar periods and average
amplitudes have been observed in the $(I,D)$ plane for unstable limit cycles
in comparably small domains at smaller Reynolds numbers [Kreilos and Eckhardt,
2012]. Note however that the amplitude of the observed cycles was larger than
what we display here by an order of magnitude. A very similar state has been
observed at $R=400$ (not shown here).

As the system size is increased to larger values $12\times 8$ (not shown) and
$L_{x}\times L_{z}=18\times 12$, several changes occur in this saddle. The
saddle becomes time invariant, both for kinetic energy (Fig. 7 (g)) and
velocity and vorticity fields (Fig. 7 (e,f)): it is thus an actual saddle
point. It still presents streamwise modulated streamwise velocity tubes that
occupy the whole domain, however, these velocity tubes are less intense with
$\max_{\text{space}}|u_{x}|=0.44\pm 0.01$. This leads to a steady value of
kinetic energy of $E_{k}=0.2$, a third of the time average found in the MFU
type system. Additionally, the streamwise modulation in the adjacent
streamwise velocity tubes is of varicose type. The spatial organisation and
average kinetic energy of this saddle is consistent with saddles computed at
$R=400$ in domains of comparable size [Schneider et al., 2008]. The spatial
variation of the streamwise vorticity is less sharp than in the MFU type
systems and its amplitude is still moderate $|\omega_{x}|=\mathcal{O}(0.1)$
(Fig. 7 (f)). As the system size is further increased to $L_{x}\times
L_{z}=36\times 24$, yet another change appears in the saddle which is reached
after bisections: the state becomes localised in $z$ (Fig. 8 (e)) and
modulated in $x$ to a smaller extent. the most intense streamwise velocity
tube has $u_{x}\geq 0$ in our computation, with
$\max_{\text{space}}|u_{x}|=0.78\pm 0.01$, it is flanked by pairs of negative
then positive streamwise velocity tube of comparable, though decreasing
intensity. The streamwise velocity then decays rapidly with $z$. Moreover the
modulation of streamwise velocity is localised in $x$, leading to a streamwise
localisation of streamwise vorticity (Fig. 8 (f)). The streamwise propagation
of this modulation makes this saddle a travelling wave. This travelling wave
is consistent with those observed at smaller Reynolds numbers $R=400$, in
domains slightly larger in the spanwise direction and larger in the streamwise
direction [Schneider et al., 2010].

Figure 10:  (a) Rescaled phase of the first $x$ Fourier mode of $u_{z}$
computed at $z=2$ and $y=0$ as a function of time computed for the limit cycle
found in MFU type domains (Fig. 5 (d,e), Fig. 6 (c,d)). (b) Rescaled phase of
the third and fourth $x$ Fourier mode of $u_{z}$ computed at $z=12$ and
$y=-0.75$ as a function of time computed for the travelling wave found in the
domain of size $L_{x}\times L_{z}=36\times 24$ (Fig. 8 (e,f)).

When the saddle state occupies the whole domain (at any given time for
$L_{x}\times L_{z}=6\times 4$ and any time for $L_{x}\times L_{z}=12\times 8$
and $L_{x}\times L_{z}=18\times 12$), it has spatial symmetries: it is
invariant under the operation $y\rightarrow-y$, $z\rightarrow-z$
$x\rightarrow-x+\delta x$ with $\delta x\lesssim 1$. The the domain of size
$L_{x}\times L_{z}=6\times 4$, the time periodicity, the time periodicity can
be viewed as a time symmetry. Moreover, advancing in time by half a period and
applying the afore mentioned spatial transformation leads to a space time
symmetry. Moreover, increasing the system size to $L_{x}\times L_{z}=12\times
8$ and to $L_{x}\times L_{z}=18\times 12$, means that the saddle has
respectively two and three wavelengths of the same pattern in both $x$ and $z$
direction and thus respectively introduces the symmetry groups
$\mathbb{Z}/2\mathbb{Z}\times\mathbb{Z}/2\mathbb{Z}$ and
$\mathbb{Z}/3\mathbb{Z}\times\mathbb{Z}/3\mathbb{Z}$. The symmetry and the
parity of velocity tubes numbers is broken in the larger domain $L_{x}\times
L_{z}=36\times 24$, where the state is localised. This broken symmetry is
accompanied by a breaking of the steadiness in $x$ of the saddle. Our
trajectory led to a saddle which is a wave travelling in the $-\mathbf{e}_{x}$
direction. We note that this system can also sustain a travelling wave
obtained under the operation $y\rightarrow-y$,
$\mathbf{u}\rightarrow-\mathbf{u}$ [Schneider et al., 2010]. This conditions
the direction of propagation as the celerity of the wave is related to
$\int_{y,z}u_{x}\,{\rm d}y{\rm d}z\neq 0$ because of the spanwise
localisation. More precisely, we find a carrier wave of wave number $n_{x}=4$
which is almost steady is space with an almost constant amplitude, with a
travelling modulation at lower wavenumbers. We can observe this using equation
(20, but this time at $y=-0.75$ and $z=12$) giving the rescaled phase of
Fourier modes $n_{x}=4$ (part of the carrier) and $n_{x}=3$ (modulation) in
figure 10 (b), this time at $z=12$ and $y=-0.75$ where $u_{z}$ is the most
intense and most concentrated. The rescaled phase of mode $n_{x}=4$ is near
constant and confirms us the spatial steadiness of the modulation. Meanwhile
the monotonous decrease of the phase of mode $n_{x}=3$ confirms the travel of
the pattern toward smaller $x$. The length $\Delta\psi=-24$ by which it has
decreased in the time window is consistent with the observed movement of the
pattern in the periodic domain. Note that this propagation occurs at constant
spatially averaged kinetic energy (Fig. 8 (g)).

We further our comparison of the respective properties of the reactive
trajectories: the corresponding last state at the last stage on the one hand
and the saddles which are the closest to the crossing point of the separatrix
on the other hand. In particular, we perform this comparison in view of the
properties of last states and saddles in many of observed noise induced
transitions. This comparison highlights a striking difference. In most known
stochastically forced systems, where the saddle is clearly identified, both
types of states are very similar: in those cases, all the trajectories cross
the separatrix in the neighbourhood of the saddle, in particular when the
noise variance is smaller and trajectories concentrate in phase space. In the
case of larger noise variance, trajectories tend to spread around the saddle,
instead of all avoiding it. In our case, for all system sizes, the reactive
trajectories cross the separatrix with velocity fields (Fig. 5 (b,c), Fig. 6
(b), Fig. 7 (c,d), Fig. 8 (c,d)) which remain fairly similar in scale and in
velocity and vorticity amplitude. These amplitudes and scales are all clearly
distinct from the nearest saddle on the separatrix (Fig. 5 (d,e), Fig. 6
(c,d), Fig. 7 (e,f), Fig. 8 (e,f)), whose stable manifold they cross. For all
system sizes one common difference is that the amplitude of streamwise
vortices is stronger ($\max_{\text{space}}|\omega_{x}|\lesssim 1$) on the
trajectory as it crosses the separatrix as it is in the saddle. The spatial
characteristic scale of said vortices is also smaller on the trajectory than
it is on the saddle. The same can be said of the modulation of the streamwise
velocity tubes. Other size dependent distinctions can be found. For the MFU
type system, the trajectory cross the separatrix at a kinetic energy
$E_{k}\simeq 0.02$ lower than that of the unstable limit cycle where
$E_{k}\simeq 0.06$. This value increases very little with $\beta$ in the
tested range. In the system of size $L_{x}\times L_{z}=12\times 8$ most
reactive trajectories start localising in $z$ (not shown here). All
trajectories localise in $z$ in the system of size $L_{x}\times L_{z}=18\times
12$, while the saddle is spread over the whole domain in those two systems. In
the largest system of size $L_{x}\times L_{z}=36\times 24$, streamwise
vorticity is more localised in $z$ than it is the case for the traveling wave.
We will see in the next sections (§ 3.2, 4.1) that this happens while all
trajectories have the properties of concentrated trajectories for small enough
forcing noise variance.

###  3.2 Statistical description of trajectories

We complement our description of trajectories with the discussion of the
statistical properties of large ensembles of trajectories obtained by means of
AMS. We will use diagnostics proposed in earlier studies of multistability and
rare events. We will confirm that the cases presented above for the discussion
of spatial organisation of the flow on paths are representative and we will
further compare our reactive trajectories to typical rare noise induced
transitions.

We first present the average path and standard deviation around it in the
plane formed by $\log(E_{k,y-z})$ and $\log(E_{k,x})$ (Fig. 11 (c)), in a
system of size $L_{x}\times L_{z}=6\times 4$ at $R=500$. The figure is
generated using the same procedure as in [Rolland, 2022] § 3.1.2: we compute
the average and standard deviation of $\log(E_{k,y-z})$, conditioned on the
value of $\log(E_{k,x})$ on reactive trajectories. The tube formed by the
average plus/minus the standard deviation thus indicates us where $66\%$ of
paths go (using the reasonable assumption of a Gaussian spread around the
mean). We compare the type of paths followed in the case of forcing on
component $x$ (with $554$ paths) to paths obtained under divergence free
forcing at increasing $\beta$ ($1600$ paths for $\beta=200000$, $687$ paths
for $\beta=300000$, $738$ paths for $\beta=330000$, $718$ paths for
$\beta=400000$, $740$ paths for $\beta=5000000$ and 755 paths for
$\beta=600000$). These paths have an estimated probability of occurrence that
range from $0.59\pm 0.01$ to $2.2\cdot 10^{-25}\pm 0.8\cdot 10^{-25}$ (which
can be considered as rare). Each set of paths starts from a slightly different
neighbourhood. Indeed, since $\mathcal{A}$ and $\mathcal{C}$ are defined using
the average and standard deviation of kinetic energy conditioned on no
transition occurring (§ 2.2.2), the smaller the energy injection rate is, the
smaller are $E_{k-x}$ and $E_{k,y-z}$ at the beginning of the trajectories
(Fig. 4). However, each group of paths (for decreased $\beta$ or changed
forcing type) displays the same type of evolution, with clear convexity of the
curve $\langle\log(E_{k,y-z})\rangle(\log(E_{k,x}))$. This confirms in a
statistical manner the observation that the two types of forcing eventually
lead to flows that follow the same type of paths as they build up toward
turbulence, in that they cross the separatrix at similar locations, with
similar flow patterns. Schematically, we first observe a growth of
$\log(E_{k,x})$ then a growth of $\log(E_{k,y-z})$. This indicates that the
velocity streaks would first receive energy from the forcing, mostly
indirectly, and are in an externally sustained process. This part of the
trajectory is similar to what is observed using deterministic targeted forcing
[Rigas et al., 2021]. They do so until they are intense enough to naturally
start the streaks instability and to go from an interplay of streaks and
vortices sustained by energy injection to the self sustaining process of
turbulence. In all cases, the trajectories are more and more concentrated
around the corresponding average path as $\beta$ is increased and $\epsilon$
is decreased. For the lowest values of $\epsilon$, we sampled these
trajectories from several independent AMS runs to ensure that this narrower
variance was not an effect of a smaller variability on the first part of the
path within an AMS run.

We perform a similar observation, this time using the curves of all the
reactive trajectories computed by means of AMS along with a relaxation from
the fixed point computed using bisections started from a reactive trajectory
in a system of size $L_{x}\times L_{z}=18\times 12$ at $\beta=80000$. The
paths are represented in the three dimensional phase space formed by
$\log(E_{k,x})$, $\log(E_{k,y-z})$ and the turbulent fraction $F$ (Fig. 11
(d)). This last quantity was added to discuss the effect of spatial
localisation on the paths. With a probability of observing the transition at
this energy injection rate at $\alpha=\mathcal{O}\left(10^{-5}\right)$, we
observe a concentration of paths in phase space (this can even be observed at
$\beta=6.25\cdot 10^{4}$ where $\alpha=\mathcal{O}\left(10^{-1}\right)$). In
that case we first have the growth of $E_{k,x}$ (from
$\log(E_{k,x})\lesssim-4.5$ to $\log(E_{k,x})\gtrsim-4$) followed by that of
$E_{k,y-z}$ already observed in smaller systems, while $F=0$. This corresponds
to the formation of a streamwise band where velocity streaks are intense
enough, modulated and in interplay with spanwise localised streamwise
vortices. Once this band is intense enough ($\log(E_{k,x}\gtrsim-4$,
$\log(E_{k,y-z})\gtrsim-10$) the turbulent fraction $F$ become non-zero and
turbulence expands in the $z$ direction. This is indicated by the correlated
increase of $\log(E_{k,x})$, $\log(E_{k,y-z})$ and $F$.

We note two peculiar facts. Firstly, while the paths are more and more
concentrated around their average in the MFU sized system, these averages are
distinct from one value of $\epsilon$ to another and the bands given by
plus/minus the standard deviation do not overlap. We do not observe the
typical narrower and narrower concentration around an instanton. We note
however that these bands move in phase space as $\epsilon$ is decreased, which
means that said instanton could be reached asymptotically. Secondly, we added
to the graph of figure 11 (c) the path followed in the
$\log(E_{k,y-z}),\log(E_{k,x})$ plane by the relaxation from the unstable
limit cycle which is the nearest to the transition path to turbulence, as
computed in section 3.1.3. At that size, the relaxation starts from values of
$E_{k,x}$ close to what is found in turbulence and displays a burst of
$E_{k,y-z}$ along with an overshoot of $E_{k,x}$: the path followed in that
plane is distinct from reactive trajectories.

This can be viewed as a corollary of the first observation. We similarly added
the curve followed in the three dimensional space $\log(E_{k,x})$,
$\log(E_{k,y-z})$, $F$ by the relaxation from the fixed point in the system of
size $L_{x}\times L_{z}=18\times 12$ (Fig. 11 (d)). This path can be viewed as
the relaxation path that would be followed by an instanton in this system. We
note that none of the group of transitions paths follows that relaxation path.
The relaxation starts from a low value of $E_{k,y-z}$
($\log(E_{k,y-z})\simeq-9.7$) and moderate value of $E_{k,x}$
($\log(E_{k,x})\simeq-3.9$). The zero value of $F$ is a consequence of the
value of the threshold used to compute the turbulent fraction, as the saddle
is delocalised in the whole domain. As the flow relaxes and follows the
unstable direction, the departure to the edge state grows globally in the
domain, not locally to the point at $\log(E_{k,x})\simeq-3.6$,
$E_{k,y-z}=-8.3$, where values of $e$ (see Eq. (6)) are large enough
everywhere to be above the detection threshold and the turbulent fraction
jumps to $1$. The kinetic energy contained in both streaks and vortices then
keep growing globally until the flow reaches transitional wall turbulence.
Following both relaxations thus confirms what had been visualised on the field
taken from examples.

(d)

Figure 11:  (a) Plot of a reactive trajectory obtained by means of AMS in the
MFU type domain at $\beta=600000$ (blue line), relaxations toward laminar (red
line) and turbulent flow (orange line) started from the reactive trajectory as
well as the result of the bisections (green line) in the plane formed by the
logarithm of kinetic energy on components $\log(E_{k_{,}x})$ and
$\log(E_{k,y-z})$. An inset zooms on the area where the dynamics winds around
the unstable limit cycle. (b) Plot of the same reactive trajectory, relaxation
and result of bisections, in the plane $I,D$ where $I$ is the total energy
injection and $D$ is the total dissipation (Eq. (9)), using the same colour
code. The dashed line $D=I$ is added. An inset zooms on the area where the
dynamics winds around the unstable limit cycle. (c) Average plus and minus
standard deviation of the logarithm of the kinetic energy in $y-z$ components
$\log(E_{k,y-z})$ conditioned on logarithm of the kinetic energy in $x$
component $\log(E_{k,x})$ obtained from reactive trajectories computed by
means of AMS with a divergence free forcing and increasing $\beta$ as well as
with forcing on the $x$ component. The path in the
$(\log(E_{k,x}),\log(E_{k,y-z}))$ plane of a relaxation from the saddle toward
turbulence is added in orange, with a dot indicating its starting point. (d)
Three dimensional plot of 151 reactive trajectories obtained by means of AMS
in the domain of size $L_{x}\times L_{z}=18\times 12$ at $\beta=80000$ (blue)
as well as the relaxation from the fixed point (black) in the space formed by
the logarithm of kinetic energy contained in components $\log(E_{k,x})$,
$\log(E_{k,y-z})$ and turbulent fraction $F$.

We then consider the distribution of durations of reactive trajectories, in a
system of size $L_{x}\times L_{z}=6\times 4$ at $R=500$. They are sampled
using AMS with forcing on components and divergence free forcing. We display
the distribution of durations rescaled by the sample average duration $\tau$
and the sample standard deviation $\sigma$ in figure 12 (a) in linear scale
and in figure 12 (b) in logarithmic scale. In both cases, the distributions
are compared to a normalised Gumbel distribution defined using the $\gamma$
Euler constant

$\mu(s)=\frac{\pi}{\sqrt{6}}\exp\left(-\left(\frac{\pi}{\sqrt{6}}s+\gamma+\exp\left(\frac{\pi}{\sqrt{6}}s+\gamma\right)\right)\right)\,,\,\gamma=0.577\ldots\,.$
(21)

For one degree of freedom systems, in the limit of vanishing noise variance,
it has been demonstrated that this is the distribution of duration [Cérou et
al., 2013]. It has been observed that the sampled distribution is close to a
Gumbel in more complex stochastic systems [Rolland et al., 2016, Rolland,
2018] when we have instanton phenomenology (with clear deviation when we do
not, for instance when the saddle becomes degenerate [Rolland et al., 2016]).
This distribution is also observed in systems not forced by noise [Rolland,
2022, Gomé et al., 2022]. In our case, we observe a very good match between
the sampled distribution and a Gumbel, with small differences for $t\leq\tau$
(which was also observed for instanton phenomenology in extended gradient
systems [Rolland et al., 2016]).

Along with the concentration of trajectories and the shape of the distribution
of duration of trajectories, we saw that the computed reactive trajectories
still displayed some properties of classical noise induced transitions. We
will test this comparison further in the next section with the dependence of
transition on the energy injection rate.

Figure 12: Distribution of duration of reactive trajectories rescaled by their
sample average $\tau$ and sample standard deviation $\sigma$. The trajectories
are computed in a system of size $L_{x}\times L_{z}=6\times 4$ at Reynolds
number $R=500$ under forcing on components or divergence free forcing for
several energy injection rates. A normalised Gumbel distribution (Eq. (21)) is
added to the plots. (a) in linear scale. (b) in Logarithmic scale. 5320
trajectories were sampled for $\beta=50000$ with forcing on $x$ component.
3220 trajectories were sampled for $\beta=200000$ with forcing on $y$ and $z$
components. 1271 trajectories were sampled for $\beta=110000$ with forcing on
$x$ component. 1600 trajectories were sampled for $\beta=200000$, 1399 were
sampled for $\beta=300000$ and 1055 were sampled for $\beta=400000$ with
divergence free forcing.

##  4 Effect of energy injection rate

In this section, we analyse the results obtained from a series of AMS runs
performed for increasing $\beta$ and decreasing $\epsilon$ in the MFU type
system using divergence free forcing and forcing on components. We will follow
quantities such as the probability of observing a transition $\alpha$, the
mean first passage time before observing a transition $T$ as well as the
average duration of reactive trajectories $\tau$ as a function of $\beta$ and
$\epsilon$. We will investigate whether we find exponential scaling for the
first two, in our ongoing comparison of the transition with typical behaviour.

###  4.1 Probability of crossing and mean waiting time

In this section, we present the probability of crossing from a neighbourhood
of the laminar flow to turbulent flow under a divergence free forcing as well
as the corresponding mean first passage time. Both of them will be presented
as a function of $\beta$ with a fixed spectrum shape. The mean first passage
times are computed by means of AMS, following the procedure described in
section 2.2.1. The probability of crossing is computed by means of AMS for all
$\beta$ of interest as well as by means of DNS (following the procedure of §
2.2.3) up to $\beta=3\cdot 10^{5}$. These computations were performed in the
MFU type system of size $L_{x}\times L_{z}=6\times 4$ at Reynolds number
$R=500$. The estimates are performed by averaging over the values of several
realisations, while the error bars are constructed from the $66$% interval of
confidence obtained from the sample standard deviations. Up to $\beta=2.8\cdot
10^{5}$, $N=200$ clones are used in AMS computations, $N=100$ clones are used
for $\beta=2.9\cdot 10^{5}$ and $\beta=3\cdot 10^{5}$ and $N=50$ clones are
used for computations at $3.1\cdot 10^{5}\leq\beta\leq 4\cdot 10^{5}$. In
these simulations, one clone was suppressed at each iteration. We performed
further computations at $\beta=5\cdot 10^{5}$ using $N=800$ clones and
removing $N_{c}=160$ clones per iteration and at $\beta=6\cdot 10^{5}$ using
$N=400$ clones and removing $N_{c}=80$ clones per iteration. At least 10
realisations of AMS computations per control parameters are performed for
$\beta<4\cdot 10^{5}$. This number goes up to more than 30 for $2.9\cdot
10^{5}\leq\beta\leq 3.3\cdot 10^{5}$, while it is of the order of several
units for $\beta\geq 4\cdot 10^{5}$.

We present the probability to observe a build up computed by means of AMS and
DNS in figure 13 (a). The corresponding mean first passage time before
crossing to turbulence, computed solely by means of AMS, is displayed in
figure 13 (b). Both are presented in logarithmic scale. These two quantities
indicate that the noise induced transition becomes rapidly very rare as
$\beta$ is increased: the estimate of $\log(\alpha)$ decreases by more than
five decades, while that of $\log(T)$ increases by five decades in the range
of $\beta$ considered here. This corresponds to variations by more than twenty
orders of magnitude in linear scale. Most of this rapid variation takes place
between $\beta=2.7\cdot 10^{5}$ and $\beta=6\cdot 10^{5}$. In that range both
$\log(\alpha)$ and $\log(T)$ are very well approximated by affine functions.
The fit are added to figure 13 (a,b) and respectively give the slopes
$(-1.53\pm 0.01)\cdot 10^{-4}$ and $(1.53\pm 0.01)\cdot 10^{-4}$. This
confirms that the probability $\alpha$ decays exponentially and that the mean
first passage time $T$ increases exponentially and that their rate are
opposite of one another. This scaling of $T$ and $\alpha$ in $\beta$ is the
same as the one observed in most systems where rare bistability is present.

Figure 13: From computations in system of size $L_{x}\times L_{z}=6\times 4$
at Reynolds number $R=500$ under divergence free forcing, by means of AMS and
up to $\beta=3\cdot 10^{5}$ by means of DNS. (a) : Probability to cross from
the defined neighbourhood of the laminar baseflow to the turbulent flow as a
function of $\beta$, an affine fit is performed for $\beta\geq 2.7\cdot
10^{5}$ and is added to the plot. (b) : Mean first passage time before
transition to the turbulent flow as a function of $\beta$, an affine fit is
performed for $\beta\geq 2.7\cdot 10^{5}$ and is added to the plot. (c) :
average duration of reactive trajectories as a function of $\beta$. (d) :
normalised standard deviation of the estimate of $\alpha$ as a function of
$\beta$, the black dashed line indicates the ideal value $1$.

We perform two tests to validate the computation of $\alpha$. Up to
$\beta=3\cdot 10^{5}$, computations by means of DNS are conceivable since
$\alpha\geq 10^{-3}$. The results are added to figure 13 (a) and show an
overlap between the AMS and DNS intervals of confidence. This indicates that
the values computed by means of AMS can be trusted in that range. This
indicates that the exponential decrease of $\alpha$ is valid in that range. We
complement this test by the computation of $\sigma_{\rm norm}$, the normalised
variance of the estimate of $\alpha$ (see Eq. (22), § C), which is presented
in figure 13 (d). In the whole range of interest, we can see that this
quantity is around one for the whole range of $\beta$ considered here. This
confirms that the estimates performed by means of AMS are not tainted by an
uncontrolled error, and that the exponential decrease of $\alpha$ is valid in
the whole range of interest. Further tests of validity of computations using
this indicator are presented in appendix C.

###  4.2 Duration of trajectories

After examining whether the probability of crossing and the mean first passage
time before crossing displayed typical exponential behaviour with $\beta$, we
will test another statistical property of reactive trajectories: their average
duration. We will investigate whether we find the typical slow growth with
$\beta$, and we will separate the case of the divergence free forcing and the
case of the forcing on components, where we will further test the effect of
the energy injection rate $\epsilon$.

Figure 14: (a) Logarithm of the duration of reactive trajectories as a
function of $\beta$ for various forcing shapes (spectrum, resolution, forced
components) computed by means of AMS in the MFU size system. (b) Logarithm of
these duration of reactive trajectories as a function of $1/\epsilon$ for
various forcing shapes (spectrum, resolution, forced components).

####  4.2.1 Forcing on components

The average duration $\tau$ of build up trajectories is computed as a function
of $\beta$ for several shape of the spectrum of the forcing on components
($\sigma_{y,z}=0$ or not, several values of $N_{r}$ at several resolutions).
We first observe the logarithm of $\tau$ as a function of $\beta$, for all
these cases in figure 14, (a). We note that, as is commonly reported, $\tau$
increases with $\beta$ at fixed spectrum shape. Moreover, the growth of the
logarithm of $\tau$ is faster than linear with the parameter $\beta$
controlling the energy injection rate independently of the structure of the
forcing. We can observe this plot differently and consider $\tau$ at fixed
$\beta$ as the spectrum shape is changed. In particular, we consider the
effect of a decreased energy injection rate as less modes are forced (as
indicated in Eq. 8). Among forcing on all three components, we note that as we
decrease $N_{r}$ from $8$ to $4$, thus reducing the range of low wavenumbers
where the spectrum is white and consequently reducing $\epsilon$, $\tau$ is
larger and larger. At fixed $N_{r}=4$, as we go from forcing on three
components to forcing on the one streamwise component, $\tau$ is even larger.
From these two observations, we can conclude that our estimates of $\tau$ all
increase as the energy injection rate $\epsilon$ is decreased: indeed,
equation (8) indicates that both increasing $\beta$ and decreasing $N_{r}$
reduces the energy injection rate. We now test this dependence even further
and ask whether $\tau$ is solely a function of the energy injection rate in
this forcing family. For this matter, we display $\tau$ as a function of
$1/\epsilon$ (Fig. 14, (b)). The curves of
$\log(\tau)\left(\frac{1}{\epsilon}\right)$ are regrouped when we use this
rescaling. Moreover, four of the six cases tested here, which have the same
numerical resolution, collapse together. As such, we have a behaviour similar
to what had been observed on the averaged kinetic energy conditioned on no
transition to occur (§ 2.2.2, Fig. 4). The matter is not trivial: we expect
$\tau$ to generally grow with $\frac{1}{\epsilon}$: more energy injected per
unit time should mean faster travel toward turbulence. However, there is also
no reason why the shape of the spectrum of the forcing noise should have no
subtle effect, we have not examined the energy transfer, direct and indirect,
in detail. From these AMS computations, we also estimated the probability of
crossing and the mean first passage time before crossing (not shown here).
Similarly to what has been observed in the former section, we have
respectively a rapid decrease of $\alpha$ and a rapid increase of $T$ as
$\beta$ is increased. Moreover, we also note that the effect of the spectrum
shape is similar to what has been observed with $\tau$ : at fixed $\beta$,
injecting less energy leads to a decrease of $\alpha$ and an increase of $T$.
Observing those quantities as a function of $1/\epsilon$ regroups the plots
together, but as tightly as for $\tau$.

From all this we are drawn to address two matters

  1. i)

Firstly, since we find the same dependence of average kinetic energy, the same
type of path (in term of stages) and the same type of dependence of their
duration on $\epsilon$, regardless of the forcing structure within this
forcing family, we can further back the assertion that the computed reactive
trajectories are selected by the non linear dynamics and not by the forcing.
Behind these tests, we aim to investigate the question of universality in
these paths: to what extent are they a product of the Navier–Stokes equations
alone and the energy injection rate as a sole parameter, and to what extent do
they depend on the forcing spectrum.

  2. ii)

Secondly, the growth of $\tau$ with the parameter controlling the noise
variance (or energy injection rate) is apparently exponential: its logarithm
is convex. As such it is quite faster then what had been previously observed
in simpler systems [Rolland and Simonnet, 2015, Rolland et al., 2016]. Indeed,
it has been demonstrated in one degree of freedom systems having a simple
saddle between both metastable states [Cérou et al., 2013] and verified
numerically in such systems having more DoF [Rolland and Simonnet, 2015,
Rolland et al., 2016] that the duration of reactive trajectories grew
logarithmically with the variance of the noise forcing. When the saddle
between both states become extremely flat, the reactive trajectories undergoes
a random walk on the plateau between both states. In that case, a faster
linear increase of the duration of reactive trajectories is measured [Rolland
et al., 2016]. However, this leads to a Gaussian distribution of trajectory
durations, unlike the Gumbel observed here (Fig. 12). An exponential increase
can be seen in simple models if there is a local (weakly) metastable minimum
standing between the two main metastable states. In that case, the exponential
increase is caused by the waiting time before escape of the secondary
metastable state [Rolland and Simonnet, 2015]. In the case of build up of
turbulence, there is no known intermediate metastable state. However the
complexity of phase space could play a similar role.

####  4.2.2 Divergence free forcing

We present the average duration of reactive trajectories obtained under
divergence free forcing as a function of $\beta$ in figure 13 (c). In contrast
with the behaviour of reactive trajectories obtained with forcing on
components, the average duration of reactive trajectories grows rather slowly.
Indeed, we can observe an increase of $\tau$ with $\beta$ up to $\beta\lesssim
260000$, followed by a plateau for $260000\lesssim\beta\lesssim 3300000$.
Above that value of $\beta\gtrsim 330000$, the average duration of reactive
trajectories estimated by means of AMS is smaller. As such, it is not
straightforward to precisely fit a known functional dependence of
$\tau(\beta)$ on that data, even if the slow logarithmic growth would be the
closest to what we estimated. There are apparent distinctions with that
behaviour, such as the local maximum. Such maxima are not commonly reported in
simpler noise driven transitions. An effect of the large error bars cannot be
ruled out, especially since no change of tendency is observed in
$\log(\alpha)$, $\log(T)$ nor in the trajectories.

##  5 Conclusion

###  5.1 Summary

In this text, we studied the trajectories going from laminar flow to
transitional turbulence in plane Couette flow under a stochastic forcing, in
systems of increasing sizes. The energy injection rate was decreased so that
the mean waiting time before observing said trajectory became too large to be
estimated by means of DNS computations. A rare events simulation method,
Adaptive Multilevel Splitting, was therefore used for the systematic
computations of trajectories along with the probability that they occur and
the mean first passage time before they occur. In order to determine whether
trajectories crossed a separatrix, and in that case where they crossed the
separatrix, several diagnostics were used. Firstly, the last state at the last
branching stage was systematically computed. Moreover, dynamics without
forcing, that relax either to the laminar or to the turbulent flow, were
performed from the successive states on the trajectory. Finally, using the
last state relaxing to the laminar flow and the first state relaxing to the
turbulent flow, bisections were performed until they converged to an unstable
fixed point or an unstable limit cycle. These computations were performed in
domains of increasing sizes, from a MFU type flow of size $L_{x}\times
L_{z}=6\times 4$ to a system of size $L_{x}\times L_{z}=36\times 24$, to
investigate effects of localisation. In the MFU type domains, the effect of
decreasing energy injection rate was investigated to study how these
transitions became rare. The numerical computations were validated: doubling
of modes in each direction indicated that numerical convergence with respect
to space and time discretisation was obtained, while comparison of AMS
computations with DNS and tests of the variance of the estimator indicated
that reactive trajectories were correctly computed. The flow displayed on the
paths was described, in view of the literature on transition to turbulence.
The avoidance of the saddles by the computed reactive trajectories was
particularly noted. Finally, the rarity of transitions was highlighted through
a parametric study in energy injection rate in the MFU type domain.

###  5.2 Mechanical description of the paths

Observation of the trajectories indicated that no matter the domain size and
the structure of forcing, the flow responded to the added noise by generating
streamwise invariant streamwise velocity tubes, along with weak streamwise
vorticity. On the trajectory, these tubes gradually accumulate energy obtained
from the forcing, in that way the non trivial flow is sustained by the
statistically steady energy injection. The tubes finally become intense enough
for the self sustaining process of turbulence to start, after the trajectory
has crossed the separatrix between laminar and turbulent flow. It is notable
that plane Couette flow was efficiently able to generate these tubes, even if
the wavenumber forcing range was distinct from the natural wavenumber range of
transitional turbulence. In that respect, these trajectories are distinct from
the one where the wavenumber of velocity streaks and streamwise vortices are
directly non-linearly forced [Rigas et al., 2021]. This indicates that the
flow is able to redistribute the energy from one range of wavenumbers to
another, even in the regime where the energy injection rate is small.
Moreover, this hints that the resulting coherent structures (in the form of
streamwise velocity tubes) are naturally selected by plane Couette flow, not
imposed by the forcing. In that respect, this would mean that we observed a
non linear receptivity of the flow in the first part of the trajectory.
Because that receptivity was independent from the several types of forcing we
applied, it could indicate some universality in the response of the flow to
the forcing. That search for universal features in the transition, such as a
mediator flow state arising solely from the Navier–Stokes equations in plane
Couette flow, and the corresponding mechanical processes, is then in line with
the computation of the properties of the separatrix such as the saddle on it
or the point on it which is the closest to the laminar flow, as measured by
kinetic energy. The build up is distinct from these last two flow states
because minimal seeds are precisely defined non linear flows and would not
naturally emerge from finite variance noise and, as we will further discuss,
edge states are avoided by the computed build up trajectories. What we observe
here is closer to (non)-linear response, with the distinction that we have not
optimised the forcing, and we are not observing the optimal response to any
regular enough forcing.

Similarly to paths to turbulence identified using other approaches (forcing,
bisections, optimisations), the reactive trajectories computed here display a
spatial localisation of the triggered velocity streaks as the domain size is
increased. This localisation is visible in the last state at the last stage.
The main difference is that in those flow states the actual self sustaining
process of turbulence is triggered in a narrow region of the domain.

###  5.3 Statistical description of the paths

After regrouping our observations on the spatial organisation of the coherent
flow structures encountered during the build up, we now discuss the build up
in a statistical manner. This means considering trajectories as a group,
histograms and averages, their parametric dependence, in view of what is known
on bistability in other noise driven system.

On several aspects, there is a similarity between the trajectories of build up
of turbulence under forcing and that of collapse of turbulence. The paths
followed in the plane $\log(E_{k,x})-\log(E_{k,y-z})$ both show a strong
convexity. The succession of events : first the generation of streamwise
velocity tubes followed by the appearance of streamwise vorticity (with
spanwise localisation when domain size is increased) can qualitatively seem as
a rewind of those observed in the collapse of turbulence [Rolland, 2022, Liu
et al., 2021, Gomé et al., 2022]. Elements of this scenario had been noted in
earlier studies without being at the centre of attention yet [Philip and
Manneville, 2011, Chantry and Schneider, 2014]. However, the build-up and
collapse trajectories are _by no mean_ time reversed paths of one another.
Firstly the aforementioned similarity is not exact. Secondly, they take place
under clearly distinct conditions: the build up is forced, while collapses are
purely deterministic. A corollary is that build up paths cross a separatrix,
with the help of the forcing, while the collapse paths do not: they avoid it.
This induces further differences, such as the duration of each phase: the
generation of the streamwise vorticity tubes shows a polynomial time
dependence, while the decay of said tube is exponential in time [Rolland,
2022]. The distinction between the forward and backward types of paths is also
different from what is found in noise driven non-gradient system. In those
system, there is a clearly distinction between the backward and forward path.
However, they both comprise of a fluctuation path toward the same saddle,
followed by a relaxation path. Meanwhile, the build up has a fluctuation path
toward the separatrix then relaxes, while the collapse can be viewed as a pure
relaxation path.

The relaxations that preceded bisections indicate that the reactive
trajectories do cross the separatrix between laminar and turbulent flow. They
do so with a kinetic energy which is several orders of magnitude larger than
that of minimal seeds, typically of the order of
$\mathcal{O}\left(10^{-6}\right)$. It is doubtful that this difference could
be solely explained by the fact that the typical slow power law decrease of
such energies starts at larger Reynolds numbers. The ensuing bisection
indicates that the separatrix is crossed in the stable manifold of the edge
state, which also has a relatively large kinetic energy. For all sizes and
Reynolds numbers considered, this edge is fairly distinct from the flow
configuration observed on the last state at the last stage and in the two flow
states on the reactive trajectories just before and just after the crossing of
the separatrix. The edge has larger spatial scales and less intense, more
coherent structures in wall normal and spanwise components than what is
observed in the reactive trajectories. The edges are less localised than the
trajectories in larger domains. In the MFU type domain, The unstable limit
cycle has more intense velocity streaks than what is seen on trajectories.