# Collapse of transitional wall turbulence captured using a rare events
algorithm

Joran Rolland 111<EMAIL_ADDRESS>222Laboratoire de Physique à
l’ENS de Lyon, UMR CNRS 5672, Univ. Lyon, Univ. Claude Bernard Lyon 1, France.
333 UMR CNRS 9014 - LMFL - Laboratoire de Mécanique des Fluides de Lille -
Kampé de Fériet, Univ. Lille, Centrale Lille, ENSMA, ONERA, France.

###### Abstract

This text presents one of the first successful applications of a rare events
sampling method for the study of multistability in a turbulent flow without
stochastic energy injection. The trajectories of collapse of turbulence in
plane Couette flow, and their probability and rate of occurrence are
systematically computed using Adaptive Multilevel Splitting (AMS). The AMS
computations are performed in a system of size $L_{x}\times L_{z}=24\times 18$
at Reynolds number $R=370$ with an acceleration by a factor $\mathcal{O}(10)$
with respect to DNS and in a system of size $L_{x}\times L_{z}=36\times 27$ at
Reynolds number $R=377$ with an acceleration by a factor
$\mathcal{O}(10^{3})$. The AMS results are validated with a comparison to DNS
in the smaller system. Visualisations indicate that turbulence collapses
because the self sustaining process of turbulence fails locally. The
streamwise vortices decay first in streamwise elongated holes, leaving
streamwise invariant streamwise velocity tubes that experience viscous decay.
These holes then extend in the spanwise direction. The examination of more
than a thousand trajectories in the $(E_{k,x}=\int u_{x}^{2}/2\,{\rm
d}^{3}\mathbf{x},E_{k,y-z}=\int(u_{y}^{2}/2+u_{z}^{2}/2)\,{\rm
d}^{3}\mathbf{x})$ plane in the smaller system confirms the faster decay of
streamwise vortices and shows concentration of trajectories. This hints at an
instanton phenomenology in the large size limit. The computation of turning
point states, beyond which laminarisation is certain, confirms the hole
formation scenario and shows that it is more pronounced in larger systems.
Finally, the examination of non-reactive trajectories indicates that both the
vortices and the streaks reform concomitantly when the laminar holes close.

##  1 Introduction

Many turbulent flows of aerodynamical or geophysical interest are not
homogeneous nor isotropic. They often display several possible turbulent flow
configurations and rare switches between these flow configurations. For these
reasons these flows are often termed multistable. Examples of multistability
in turbulent flows include turbulent dynamos (Berhanu et al.,, 2007),
turbulent convection (Podvin and Sergent,, 2017), bluff body wakes
(Grandemange et al.,, 2013), jets in the wake of a pair of cylinders (Kim and
Durbin,, 1988) and barotropic atmospheric-type jets (Bouchet et al., 2019a, ;
Simonnet et al.,, 2021). In order to understand these turbulent flows,
uncovering what drives the switches is as important as explaining the
mechanisms maintaining each metastable configuration. These switches are
characterised by the mean first passage time before a change of configuration
occurs. The switches often take place through the same chain of events, termed
a _transition path_ or a _reactive trajectory_ , a notion originating from
kinetic chemistry (Metzner et al.,, 2006). It is fairly difficult to compute
mean first passage times and transition paths, either experimentally or
numerically for two main reasons.

  1. a)

One firstly has to deal with the very large number of degrees of freedom of
turbulent flows, particularly in numerical simulation. These complex flows are
often difficult to simulate even by means of Large Eddy Simulations.

  2. b)

One secondly has to deal with the extremely long waiting times between each
event. These waiting times are several orders of magnitude larger than the
duration of a realisation of a switch and even longer than the typical eddy
turnover time (Kim and Durbin,, 1988). This means that the cost of sampling
more than a few events is prohibitive by classical means (see table in
(Bouchet et al., 2019a, )).

In order to propose a method to systematically study multistability, one first
needs a turbulent system which has fewer effective degrees of freedom but is
still complex enough to display multistability, with a moderate need for
extrinsic stochastic forcing. One can thus temporarily bypass problem a). A
transitional wall flow such as plane Couette flow can provide such a situation
(see Fig. 1). Unlike thermal convection, for instance, a wall flow like plane
Couette flow is linearly stable for all Reynolds numbers (Romanov,, 1973).
Meanwhile transitional turbulence can exist at moderate Reynolds numbers,
albeit transiently (Eckhardt et al.,, 2008; Schmiegel and Eckhardt,, 1997;
Bottin and Chaté,, 1998; Eckhardt et al.,, 2007). Transitional wall turbulence
takes the form of velocity streaks: long wavy streamwise tubes flanked by
short streamwise vortices. The streaks and vortices sustain one another in a
cycle termed the Self-Sustaining Process (SSP) of wall turbulence (Waleffe,,
1997; Hamilton et al.,, 1995). The streamwise vortices extract energy from the
mean flow to regenerate the streaks, a process termed lift-up. Meanwhile the
streaks regenerate the streamwise vortices through a Kelvin-Helmholtz
instability and vorticity tilting by the base flow. How does wall turbulence
collapse, that is to say whether the SSP fails because vortices disappear
first or streaks disappear first or both collapse at the same time, is a
question that could only be address very recently (Gomé et al.,, 2020; Liu et
al.,, 2020; Rolland,, 2018). Conversely, turbulence can build up from laminar
flow under a finite amplitude forcing (Wan and Yu,, 2017; Rolland,, 2018; Liu
et al.,, 2020), or somewhat equivalently from a finite amplitude initial
condition (Faisst and Eckhardt,, 2003). This gives a first type of
multistability, where the transitional turbulence of plane Couette flow can
collapse down to laminar flow under its own fluctuations, and go back to
turbulence if it is forced. In a way this is a purely temporal view. This
image should be completed by noting that transitional wall turbulence tends to
be localised (the type of organisation depends on the flow geometry). At
moderate Reynolds numbers, wall flows where transitional turbulence extends in
one dimension (such as Hagen–Poiseuille flow (Moxey and Barkley,, 2010),
tilted plane Couette (Shi et al.,, 2013) flow or tilted plane Poiseuille flow
(Gomé et al.,, 2020)) can display splitting: the turbulent puff elongates and
splits in half. This can lead to an effective extension of the area occupied
by turbulence. This can provide a third type of multistability events in some
flow configurations. When all these processes are taken into account, one can
find laminar-turbulent coexistence in large flow domains containing several
puffs. In models, it can be shown that the laminar-turbulent coexistence is
also transient, and that its lifetime grows exponentially with the system size
(Rolland,, 2018). The laminar-turbulent coexistence and the possibility of
reinvasion of laminar holes through puff splitting are actually preeminent
ingredients in the collapse scenario and size scaling of lifetimes in the
model. Transitional wall flows are relevant for high Reynolds number wall
turbulence, since the processes described here still take place very close to
the wall (up to thirty wall units Pope, (2001)), in a layer termed the buffer
layer, or the viscous layer when a slightly thicker layer is taken into
account. Moreover, canonical flows such as plane Couette can thus provide a
good laboratory to test rare events methods on turbulent flows. These
multistable turbulent flows are not so complex that they cannot be simulated
by Direct Numerical Simulations (DNS).

Even if the number of degrees of freedom is reduced, as is the case for
transitional wall flows, one still has to deal with problem b). Transitional
turbulence, like other multistable systems, displays very long waiting times
before a multistability event occurs (Eckhardt et al.,, 2007; Bottin and
Chaté,, 1998; Shi et al.,, 2013; Gomé et al.,, 2020). This means that even for
the study of turbulence collapse in plane Couette flow, the use of DNS to
sample many events comes with a prohibitive computational cost. There are of
course alternative methods to study rare events. Many of them originate from
the study of kinetic chemistry and borrow much of its vocabulary. All these
methods aim at performing the computation of the reactive trajectories, the
probability that a reactive trajectory occurs, and the mean waiting time
before they occur. These methods can be divided in three main families
(Bouchet et al., 2019b, ), two of them (action minimisation and cloning
methods) will be invoked in this text.

  * •

Firstly, one finds mostly theoretical, optimisation methods, which are applied
to systems where a stochastic forcing is clearly identified. In that case, in
the limit where the variance of said forcing goes to zero, all the transition
paths concentrate around an _instanton_ (Touchette,, 2009) § 6. The instanton
represents the most probable transition path and can be computed using action
minimisation (Grafke and Vanden-Eijnden,, 2019; Wan and Yu,, 2017). One
fundamental property of said instantons is that they evolve from a first
multistable state toward a saddle point of the deterministic dynamics under
the action of the noise. The instanton then evolves freely from the saddle
point toward the second multistable state. The mean first passage time can be
estimated using the result of said action minimisation and can loosely be
thought of as depending on the distance between the first metastable state and
the saddle point. In gradient stochastic systems, this yields the celebrated
Eyring–Kramers formula, also known as the Arrhenius law (Hänggi et al.,,
1990). Such results can be extended to non gradient systems (Bouchet and
Reygner,, 2016). While these methods give a lot of qualitative results on the
structure of transition paths, they can prove tricky to implement for fluid
flows (Wan and Yu,, 2017). A key property of these formulae and methods is
that they are formulated as _large deviations_ (Touchette,, 2009). This means
that probabilities $\alpha$, probability density functions $\rho$, rates of
transitions $1/T$ are considered in the limit of a vanishing small parameter
$\epsilon\rightarrow 0$. In that case, we have that $\lim_{\epsilon\rightarrow
0}-\epsilon\log(\alpha)=I_{\alpha}$, $\lim_{\epsilon\rightarrow
0}-\epsilon\log(\rho)=I_{\rho}$, $\lim_{\epsilon\rightarrow
0}\epsilon\log(T)=I_{T}$, with $I_{\alpha}$, $I_{\rho}$ and $I_{T}$
independent on $\epsilon$. In other words, the leading dependence is
exponential, and noted as $\alpha\underset{\epsilon\rightarrow
0}{\asymp}e^{-\frac{I_{\alpha}}{\epsilon}}$.

  * •

Secondly, one finds mostly numerical cloning methods, which use the actual
fluctuating dynamics of the system and push them toward realisations of the
reactive trajectories. These methods compute the transition paths using $N$
clone dynamics of the system and apply a mutation selection procedure to
compute the reactive trajectories. One such method is termed _Adaptive
Multilevel Splitting_ (AMS) (Cérou and Guyader,, 2007; Cérou et al., 2019b, ).
AMS and its variants have been successfully used to compute reactive
trajectories and extreme events in kinetic chemistry (Lopes and Lelièvre,,
2019) theoretical physics models (Rolland et al.,, 2016), models of
transitional flows (Rolland,, 2018), idealised atmospheric flows (Bouchet et
al., 2019a, ; Simonnet et al.,, 2021). Some variants have been applied to the
study of extreme two dimensional turbulent wakes (Lestang et al.,, 2018, 2020)
and oceanic flow reversals (Baars et al.,, 2021). Thirdly, one can use
importance sampling methods (l’Ecuyer et al.,, 2009; Hartmann et al.,, 2019).
These methods modify the dynamics so that the events of interest can then be
sampled according to a new probability distribution function (see (Lestang et
al.,, 2018), Fig. 8 (Ragone and Bouchet,, 2019) Fig. 3 or (Ragone and
Bouchet,, 2020), Fig 1 (a,b)). Under this new distribution, the event of
interest is much more probable and is therefore sampled much more precisely.
Conversion factors are then used to rescale the estimated mean first passage
times and probabilities. Importance sampling is often performed in stochastic
systems, where the modifications are applied on the noise. As a consequence,
this often makes importance sampling hard to apply to deterministic systems.
The key question is then what modification to apply and how to compute the
rescaling factors. When answering these questions, one can note that the
boundaries between optimisation, cloning and importance sampling methods are
porous. For instance, results of action minimisation can be used to design
importance sampling methods, when studying extreme events on a given time
interval (Grafke and Vanden-Eijnden,, 2019; Ebener et al.,, 2019). Similarly,
some cloning methods lead to a situation comparable to importance sampling
(Ragone and Bouchet,, 2019, 2020), where the rescaling of probabilities is
based on Large Deviations for time averaged variables.

Action minimisation and cloning methods are often used hand in hand. Once a
small parameter that controls the effective noise variance has been
identified, the qualitative insight from theory is used to guide numerical
studies and propose manners in which results can be presented, as it has been
done in models of pipe flow (Rolland,, 2018). Applying such a program to the
collapse of turbulence in plane Couette flow is actually not that
straightforward. Because energy is not injected in the flow through a
stochastic forcing whose variance decays, the small parameter is not readily
identified. For similar reasons, cloning methods cannot be applied as such
(Lestang et al.,, 2020). Applying basic cloning rules leads to so called
extinction: the flow does not separate trajectories from one another and the
method does not manage to create a reactive trajectory which contains an
excursion far enough from the starting metastable state. A first goal is
therefore to propose a modified cloning method that can bypass this problem
and can at least succeed in computing reactive trajectories faster than a DNS
would. This is the purpose of anticipated AMS, which is presented and used in
this text. Once the reactive trajectories are computed, if they display
concentration around a typical transition path, one needs to make sense of
this concentration. For this purpose, one can use both AMS and DNS to identify
the right small parameter from the study of probability density functions,
probabilities of transition, mean first passage times _etc_. We will keep this
in mind in our study. We also note that there are examples of very relevant
reactive trajectories (see (Rolland et al.,, 2016) for instance) that are not
necessarily of the instanton type. As a consequence, one does not always have
to try to force the results into a large deviations framework.

We present the study of the collapse of transitional turbulence of plane
Couette flow in the following manner. We describe plane Couette flow in
section 2.1. We then present Anticipated Adaptive Multilevel Splitting in
section 2.2. The generation of initial conditions is presented in section 2.3.
The reaction coordinate used to compute reactive trajectories is presented in
section 2.4. The DNS which are used as reference are presented in the next
section (§ 2.5). We then present the systematic comparison of reactive
trajectories computed by AMS and DNS in a system of size $L_{x}\times
L_{z}=24\times 18$ (§ 3.1). We perform the validation of the computation of
the probability of crossing and mean first passage time in this system in
section 3.2. AMS is then applied to the computation of very rare trajectories
and laminar hole formation in section 4. These results are finally discussed
together in the conclusion (§ 5).

Figure 1: Overview of transition to turbulence: the bottom axis indicates the
Reynolds number and the text on the top indicating typical states. The left
panel displays a sketch of plane Couette flow. The middle panel displays
colour levels of the kinetic energy density in a $y=0$ plane at $R=370$,
showing banded laminar turbulent coexistence. The green box indicates a domain
of size $L_{x}\times L_{z}=36\times 27$, the blue box indicates a domain of
size $L_{x}\times L_{z}=24\times 18$. The right panel displays colour levels
of the kinetic energy in a $y=0$ plane at $R=500$.

##  2 Flow configuration and numerical procedures

###  2.1 Plane Couette flow

We will perform the study of collapse in plane Couette flow, the flow between
two parallel walls located at $y=h$ and $y=-h$, respectively moving at
velocities $U\textbf{e}_{x}$ and $-U\textbf{e}_{x}$ (Fig. 1, left). We term
$\textbf{e}_{x}$ the streamwise direction, $\textbf{e}_{y}$ the wall normal
direction and $\textbf{e}_{z}$ the spanwise direction. Lengths are
nondimensionalised by $h$, velocities are nondimensionalised by $U$ and times
by $h/U$. The first and foremost control parameter is the Reynolds number
$R=hU/\nu$, with $\nu$ the kinematic viscosity. The nondimensional streamwise
and spanwise sizes $L_{x}$ and $L_{z}$ are two other control parameters of the
system. The full velocity field is written
$\textbf{v}=y\textbf{e}_{x}+\textbf{u}$, where $y\textbf{e}_{x}$ is the
laminar baseflow.

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

$\frac{\partial u_{l}}{\partial t}+u_{m}\frac{\partial u_{l}}{\partial
x_{m}}+y\frac{\partial u_{l}}{\partial x}+\delta_{l,x}u_{y}=-\frac{\partial
q}{\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)+f_{l}(x,t)\,,\,\frac{\partial u_{m}}{\partial x_{m}}=0\,,$ (1)

using tensorial notations. We include the term $\mathbf{f}$. It is a very
general forcing which can be switched on or off. When it is on, it is white in
time and in $y$. It can be red or white in $x-z$: in AMS simulations, we will
use a temporally localised fully white perturbation (see § 2.2.1), while in
perturbed DNS, we use a red forcing (see § 2.3, 2.5 for details). These
equations are discretised in space on $N_{x}$ and $N_{z}$ de-aliased Fourier
modes (so that $\frac{3}{2}N_{x}$ and $\frac{3}{2}N_{z}$ modes are used in
total) in the streamwise $\textbf{e}_{x}$ and $\textbf{e}_{z}$ directions and
$N_{y}$ Chebyshev modes in the $\textbf{e}_{y}$ direction. Time integration is
performed using channelflow, by J. Gibson (Gibson et al.,, 2008). We
investigate in detail systems of two sizes. The smaller system has size
$L_{x}\times L_{z}=24\times 18$ (see Fig. 1 for scale), we set $N_{y}=27$,
$\frac{3}{2}N_{x}=128$ and $\frac{3}{2}N_{z}=96$. The larger system has size
$L_{x}\times L_{z}=36\times 27$, we set $N_{y}=27$, $\frac{3}{2}N_{x}=196$ and
$\frac{3}{2}N_{z}=144$. In both cases we will set a constant timestep $\Delta
t$ during the time integrations, $\Delta t=0.05$ (at $R=370$ and $R=377$) and
$\Delta t=0.02$ (at $R=600$). Both these values ensure stability of the time
integration and respect the CFL criterion. We do not use adaptive time steps
in order to have more control on trajectory reconstruction and effect of the
time step on convergence.

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\,.$ (2)

We will also distinguish the kinetic energy contained in the streamwise
component $E_{k,x}$, on the one hand, and the kinetic energy contained in the
spanwise and wall normal components $E_{k,y-z}$, on the other hand

$\displaystyle
E_{k,x}=\frac{1}{2L_{x}L_{z}}\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=0}^{L_{z}}\frac{u_{x}^{2}}{2}\,{\rm
d}x{\rm d}y{\rm d}z\,,$ $\displaystyle
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\,.$ (3)

The first kinetic energy $E_{k,x}$ roughly quantifies the energy contained in
velocity streaks, while the second $E_{k,y-z}$ roughly quantifies the energy
contained in streamwise vortices (Jiménez and Moin,, 1991; Hamilton et al.,,
1995). These are the two main flow structures of transitional wall flow
turbulence. They regenerate one another in the cyclic self sustaining process
of turbulence (Waleffe,, 1997).

###  2.2 Adaptive Multilevel Splitting

####  2.2.1 Principle of the algorithm

Before presenting the principle of anticipated AMS, let us first give a formal
phase space description of the rare events we will study in this text. Let us
sketch the collapse of turbulence in figure 2 (a) and consider the set
$\mathcal{A}$, a neighbourhood of the turbulent flow in phase space, and the
set $\mathcal{B}$, a neighbourhood of the laminar flow. A realisation of the
dynamics which starts in $\mathcal{A}$ fluctuates around it, has several
excursions out of $\mathcal{C}$, a 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}$. Precise definitions of sets $\mathcal{A}$, $\mathcal{B}$
and hypersurface $\mathcal{C}$ for collapse will be given in sections 3 and 4,
based on reaction coordinates defined in section 2.4.

We then give a brief overview of the variant of AMS, termed _Anticipated
Adaptive Multilevel Splitting_ , which was used for the computation presented
in this text (see (Cérou and Guyader,, 2007; Rolland and Simonnet,, 2015;
Bréhier et al.,, 2016; Rolland,, 2018; Cérou et al., 2019b, ; Lestang et al.,,
2020) for more details on the general methods). All variants of AMS use a
reaction coordinate (or observable) $\phi(\mathbf{u})$, a real-valued function
of the velocity field. The reaction coordinate gives a relative distance in
phase space between $\mathbf{u}$ and the sets $\mathcal{A}$ and $\mathcal{B}$.
For practical reasons, we will see the reaction coordinate as a function of
time $\Phi(t)=\phi(\mathbf{u}(t))$ on the trajectories. The reaction
coordinate is often rescaled such that $\phi(\partial\mathcal{A})=0$, on the
boundary of set $\mathcal{A}$, $\phi(\partial\mathcal{B})=1$, on the boundary
of set $\mathcal{B}$, and grows monotonously inbetween. All variants run $N$
clone dynamics of the system to compute iteratively at least $N-N_{c}>0$
reactive trajectories going from a hypersurface $\mathcal{C}$, close to set
$\mathcal{A}$, towards the set $\mathcal{B}$. At each iteration, $N_{c}<N$
clones are replaced. The algorithm is sketched in figure 2 (b) and proceeds in
the following manner:

  * •

There is a first stage of natural dynamics, where each clone dynamics starts
inside set $\mathcal{A}$. As much as possible, these initial conditions should
be distributed according to the natural flow, restricted to $\mathcal{A}$ (see
§ 2.3 for an example of procedure). We let all the initial conditions evolve
according to their natural dynamics until they cross $\mathcal{C}$ and we stop
them when they reach either $\mathcal{A}$ or $\mathcal{B}$. We set the number
of iterations $\kappa=0$.

  * •

In a second stage, the algorithm iterates the mutation-selections as long as
there are less than $N-N_{c}$ clones that transit to $\mathcal{B}$.

At each iteration, we set $\kappa=\kappa+1$. Then all the clones are ordered
with index $i$, $1\leq i\leq N$, according to the maximum of the reaction
coordinate on the trajectory $\max_{t}\Phi_{i}(t)$. The clones $1\leq i\leq
N_{c}$, that have the shortest excursion out of $\mathcal{A}$ are removed from
the set of clones. In order to keep a constant number of clones, $N_{c}$ new
clones are created by branching on $N_{c}$ clones, labeled $i^{\prime}$, drawn
uniformly out of the $N-N_{c}$ non-removed clones at level $\Phi_{b}$, a
function (specified later in the text) of $\Phi_{N_{c}}$ (the maximal value of
reaction coordinate reached by the removed trajectories). The branched
trajectories first share the dynamics of the clone on which they are branched,
_i.e._ we set
$\\{\mathbf{u}(t),q(t)\\}_{i}=\\{\mathbf{u}(t),q(t)\\}_{i^{\prime}}$ from
$t=0$ to the time where clone $i^{\prime}$ first reach $\Phi_{b}$. Then, the
branched clones follow their natural dynamics until they reach either
$\mathcal{A}$ or $\mathcal{B}$, with a new realisation of the noise. One may
have to repeat this branching operation several times to make sure that the
maximum of $\Phi$ on the branched trajectory is strictly larger than
$\Phi_{N_{c}}$.

The algorithm stops when $Nr\geq N-N_{c}+1$ clones trajectories have reached
$\mathcal{B}$. We usually perform $o>1$ independent AMS runs. In each of the
runs labeled by the number $a$, $1\leq a\leq o$, we obtain the total number of
iterations $\kappa_{a}$, and $r_{a}$ the proportion of clones that transit to
$\mathcal{B}$. This yields an estimator of the probability $\alpha$ of
reaching $\mathcal{B}$ before $\mathcal{A}$ (Cérou and Guyader,, 2007), and
the corresponding mean first passage time $T$ (Cérou et al.,, 2011)

$\alpha_{a}=r_{a}\left(1-\frac{N_{c}}{N}\right)^{\kappa_{a}}\,,\,\hat{\alpha}=\left\langle\alpha_{a}\right\rangle_{o}\,,\,\hat{T}=\left\langle\left(\frac{1}{\alpha_{a}}-1\right)(t_{1,a}+\tilde{\tau}_{a})+(t_{1,a}+\tau_{a})\right\rangle_{o}\,,$
(4)

where $\tau_{a}$ is the mean duration of reactive trajectories, $t_{1,a}$ is
the mean duration to go from $\mathcal{A}$ to $\mathcal{C}$, and
$\tilde{\tau}_{a}$ is the mean duration of non-reactive trajectories, computed
in each AMS run $1\leq a\leq o$. The notation
$\langle\cdot\rangle_{o}=\frac{1}{o}\sum_{a=1}^{o}\cdot_{a}$ corresponds to an
average over the $o$ independent AMS runs (Cérou and Guyader,, 2007; Rolland,,
2018).

We will often record the velocity field, noted $\mathbf{u}_{\rm last}$ and
termed the last state at the last stage, that corresponds to
$\max_{t}\Phi_{N_{c}}(t)$ during the last stage of the algorithm. It often
gives a precise idea of the turning point in reactive trajectories. Before the
flow visits the neighbourhood of that state, returning towards turbulence is
more likely, beyond that point, relaminarising becomes more likely. In systems
which correspond to a simple deterministic part forced by noise, that state
actually corresponds to the saddle point of the deterministic part of the
dynamics crossed by the instanton in the limit of the noise variance going to
$0$. This has been verified for models with few degrees of freedom and the one
dimensional Ginzburg–Landau equation (not shown here). It can be used to educe
an effective saddle between two multistable states (Simonnet et al.,, 2021).
Dichotomy procedures have been started from states seen during turbulence
collapse (De Lozar et al.,, 2012). However, we have _a priori_ no certainty
that the field $\mathbf{u}_{\rm last}$ corresponds to an actual saddle of the
Navier–Stokes equations (as computed by dichotomy or other methods (Schneider
et al.,, 2007; Willis and Kerswell,, 2009)).

The algorithm is naturally parallelised over the $N_{c}$ suppressed clones.
This will be done for the calculation presented in this text. We usually
choose the number of threads $p$ such that $N_{c}/p$ is an integer larger than
or equal to two. Since the trajectories have a random duration, we cannot have
a perfect load balancing in this parallelisation. Note however that as
$N_{c}/p$ increases, it has been observed that the differences in trajectory
durations average out and that we can reach a reasonable load balancing
between threads.

Note that in AMS computations in deterministic systems, we add a small
background noise which helps the separation of trajectories after branching.
In practice we will switch on the forcing $\mathbf{f}$ (see Eq. (1)), white in
time and space $\langle
f_{l}(\mathbf{x},t)f_{m}(\mathbf{x}^{\prime},t^{\prime})\rangle=\frac{1}{\beta_{f}}\delta_{lm}\delta(\mathbf{x}-\mathbf{x}^{\prime})\delta(t-t^{\prime})$,
with inverse variance $\beta_{f}=10^{10}$, at branching timestep and switch it
off afterwards. The realisation of this noise is decorrelated from one
branching to another. The variance is small enough so as not to perturb the
laminar flow too much. The trade-off is that we compute trajectory properties
with a small error : $T+\delta T$, $\alpha+\delta\alpha$ and
$\tau+\delta\tau$. We will comment on the visible effects of this additional
force on the trajectories and their properties in section 3.1.2 and 3.2.

####  2.2.2 Anticipated AMS and branching level $\Phi_{b}$

Schematically, we can apply AMS to two types of systems. On the one hand, we
find systems with a large time scale separation between some fast fluctuating
degrees of freedom and slower degrees of freedom which represent the main
features of the flow travelling between $\mathcal{A}$ and $\mathcal{B}$. This
is often the case for stochastically forced systems. In these systems, two
slightly different initial conditions will quickly separate and the odds of
creating an excursion toward $\mathcal{B}$ instead of $\mathcal{A}$ by
slightly changing the noise realisation at a branching are non-negligible. On
the other hand, we find systems with absolutely no clear time scale separation
between degrees of freedom. This is often the case of purely deterministic
systems. Two slightly different initial conditions do not separate until it is
too late (they both reach $\mathcal{A}$). We can find situations where, no
matter the structure of the small perturbation (typically at a branching), the
odds of creating a further excursion toward $\mathcal{B}$ instead of
$\mathcal{A}$ can be exceedingly small. This is especially the case if we
perturb at the peak of an existing fluctuation. If we apply basic AMS to this
second type of system, where we branch at $\Phi_{N_{c}}$, we run the risk of a
so called extinction (Lestang et al.,, 2020). This occurs when all
trajectories have the same maximum of reaction coordinate but none of them
reach the arrival set, so that $\forall 1\leq i\leq
N,\max_{t}\Phi_{i}=\Phi_{\rm ext}<\phi(\partial\mathcal{B})$, where
$\partial\mathcal{B}$ is the boundary of set $\mathcal{B}$. The algorithm does
not manage to proceed any further444In that case, when extinction is detected,
the computation is terminated and the way AMS is used is reassessed.. In order
to bypass this limitation, we can perform anticipated branching, that is to
say branch the new trajectories at $\Phi_{b}<\Phi_{N_{c}}$. In that case, it
may be necessary to reiterate the branching several times in order to ensure
that the branched trajectories have $\max_{t}\Phi>\Phi_{N_{c}}$: all the
branched trajectories have to go further than the maximum level of reaction
coordinate reached by removed trajectories. In figure 2 (c), we give two
examples of relations $\Phi_{b}(\Phi_{N_{c}})$ that were tested for plane
Couette flow. Each one is adapted to a given situation. The converging
anticipation is mostly used in this article. The point is to take advantage of
higher mixing and faster separation of trajectories that take place when the
flow is closer to the fully turbulent state. For this matter, one has first
has a small branching level $\Phi_{b}\ll\Phi_{N_{c}}$ for $\Phi_{N_{c}}<0.5$.
We then almost branch at the maximal possible level
$\Phi_{b}\lesssim\Phi_{N_{c}}$ for $\Phi_{N_{c}}\geq 0.5$, so as not to lose
too much computational time rerunning trajectories when the flow is very close
to turbulence collapse. The saturated anticipation takes a very different
point of view. It has proved very efficient in very small domains where the
collapse of turbulence is very well described by transient chaos. This
modification of AMS uses the fact that there is much more mixing when the flow
is close to the turbulent state, but that mixing gradually stops during
excursions. If a trajectory is on the wrong track, there is no derailing it at
large $\Phi$. In this type of systems, the task of AMS is really about finding
the right exit point. More details on the necessity of anticipation are given
in appendix A.

(a) (b)

Figure 2: (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 two 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)). (c) Two
examples of anticipation of branching level of reaction coordinate as a
function of maximum reaction coordinate reached by suppressed trajectories
$\Phi_{N_{c}}(\Phi_{b})$ tested in anticipated branching (see § A, Eq. (12),
(13) for details, both examples use a parameter $\xi=0.25$).

####  2.2.3 On the goodness of AMS computations

We use AMS (or any other rare event simulation method) on top of a numerical
discretisation of the Navier–Stokes equations to compute (by order of
priority):

  1. a)

a large number of reactive trajectories, in our case from turbulent to laminar
flow,

  2. b)

provide an estimate of the probability that these reactive trajectories occur,

  3. c)

estimate the mean first passage time before a reactive trajectory occurs.

We wish these computations to be precise, with an emphasis firstly on reactive
trajectories, secondly on probabilities and ideally on mean first passage
times. Like every other numerical procedure, the use of AMS can lead to errors
on estimates which are deemed to large because of an inadapted reaction
coordinate or an insufficient number of clones. This is similar to the effect
of time or space steps too large that lead to excessive errors in a numerical
discretisation. Similarly, an unadapted scheme can lead to a numerical
convergence with time and space steps which is too slow.

A successful application of AMS first means that a) the resulting reactive
trajectories should faithfully represent the actual reactive trajectories that
would be computed in very long DNS (in term of path chosen, of duration of the
paths _etc_.). This can be checked by comparing the paths computed by AMS to
some reference, for instance a DNS. The comparison is usually done in a space
with few dimensions (two or three), where we display the most probable paths
as computed by AMS and by the reference and show that they go through the same
stages. If several clearly distinct types of paths are possible, one should
check that in the set of trajectories computed by AMS, there is the correct
proportion of trajectories going through each path (see for instance (Rolland
and Simonnet,, 2015), Fig. 8). In order to validate quantities like the
average duration of trajectories $\langle\tau\rangle_{o}$ (averaged over all
the sampled durations over a series of independent AMS runs), one compares the
sample mean to the sample mean of the reference. Since these means are the sum
of independent random variables, drawn from identical distributions with a
finite variance, one often invokes the Central Limit Theorem to state that the
actual average has a 66% chance of being within $\pm$ the variance of the
distribution
$\sigma_{\tau}=\sqrt{\frac{1}{o}\sum_{a=1}^{o}\tau_{a}^{2}-\langle\tau\rangle_{o}^{2}}$
divided by the square root of the number of samples (or 98% chances of being
$\pm$ twice the variance). The interval
$\left(\langle\tau\rangle_{o}-\frac{\sigma_{\tau}}{\sqrt{o}};\langle\tau\rangle_{o}+\frac{\sigma_{\tau}}{\sqrt{o}}\right)$
is termed the confidence interval. The averaged quantities should be close in
that their respective confidence interval should overlap: in that case, we
cannot assert that the estimation is biased. The result of AMS computations
should first pass this first test in order to be validated. Indeed, it may be
that incorrect trajectories are selected, usually because of a very poor
reaction coordinate (Rolland and Simonnet,, 2015; Bréhier et al.,, 2016). If
one can avoid this phenomenon, one can ensure that the more clones $N$ are
used in AMS computations and the more independent runs $o$ are performed, the
more precise the result is going to be. Furthermore the correctly computed
trajectories can be used to improve the computation of more sensitive
quantities like the probability $\alpha$ and the mean first passage time $T$

A successful application of AMS secondly means that b) the estimate of the
probability that the trajectory occurs is precise. One can perform an estimate
of $\alpha$ by averaging the result over AMS runs. The AMS runs are
independent. Moreover, one can ensures that the distribution of $\alpha_{a}$
computed in AMS runs has a finite variance
$\sigma_{\alpha}=\sqrt{\frac{1}{o}\sum_{a=1}^{o}\alpha_{a}^{2}-\langle\alpha\rangle_{o}}$.
One can provide an estimate with a 66% interval of confidence
$\langle\alpha\rangle_{o}\pm\frac{\sigma_{\alpha}}{\sqrt{o}}$ for the
probability, using the sample variance and the number of samples, and compare
this to a reference. Of course this estimate is tainted by the finite number
of clones $N$, the finite number of AMS realisations $o$, and possibly a poor
reaction coordinate. It is often observed that $\langle\alpha\rangle_{o}$
underestimates $\alpha$ with a probability close to one. However, it is
demonstrated that for most version of AMS, this quantity is estimated without
bias (Bréhier et al.,, 2016) or a controlled bias decaying with $N$ (Cérou and
Guyader,, 2007). It is conjectured that this discrepancy is caused by an
effect called the _apparent bias_ in multilevel splitting (Glasserman et al.,,
1998). A correct estimation would then require an infinite number of
realisations. Note that this also occurs when Importance Sampling is performed
(Devetsikiotis and Townsend,, 1993). The quality of the estimate of $\alpha$
can further be tested using the sample variance of $\hat{\alpha}$, which
should not be too large compared to the ideal variance. The ideal variance is
obtained numerically if one uses the ideal reaction coordinate, termed the
committor (Cérou et al., 2019a, ). If a problem is suspected, because the
confidence intervals of $\alpha$ estimated by AMS and the reference absolutely
do not overlap and/or the sample variance of $\alpha$ is too large, one can
compute the histograms of $\alpha$. When the apparent bias phenomenon occurs,
heavy powerlaw tails appear, usually toward large $\alpha$: a few very large
overestimate compensate the large majority of underestimates. If the apparent
bias phenomenon is avoided, one can ensure that the estimate of $\alpha$ is
unbiased: that this to say that there is no additive irreducible error on the
$\alpha$ computed after each AMS run (dependent on $N$ or not) and that the
average over AMS realisations will converge toward the reference when the
number of realisations $o$ goes to infinity. When errors occur, the number of
clones in each realisation can be increased. One can also rely on the
information brought by correctly computed reactive trajectories to construct a
better reaction coordinate.

The most successful applications of AMS finally means that c) the estimated
mean first passage time is precise. Again, one can provide an estimate of $T$
and an interval of confidence. Note however that unlike the properties of
trajectories or the probability of collapse, one cannot demonstrate that the
estimator of equation (4) is unbiased. In practice the relative error on the
estimate of $T$ using a small number of clones can be quite larger that the
error on the estimate of each of the separate terms involved in equation (4).
In order to reduce said bias, it has been observed that increasing the number
of clones and improving the reaction coordinate will improve the estimate of
$T$. Biases in the estimate of $T$ generally lead to overestimates. These
biases arise because the estimate of $T$ is not direct and comes from the
product of several other random variables. There are version of AMS that lead
to a more direct estimate of $T$. However, this rewriting comes with
additional constraints, such as fixed durations for trajectories (Lestang et
al.,, 2018).

All things considered, we can assert that AMS computations give realiable
results when the reactive trajectories and the probability of crossing are
computed precisely, with clear accelerations of computations with respect to
DNS. When these two quantities are correctly estimated, we can use the AMS
results to discuss the physics of the multistability of the problem we
investigate. This also means that we can reuse the information on reactive
trajectories from these computations in order to improve the reaction
coordinate. Since the probability of crossing follows the same exponential
scalings as the PDFs or the mean first passage times, it can be used as a
proxy to investigate the Large Deviations or PDF tails if we are not satisfied
with the estimate of the mean first passage time.

###  2.3 Initial condition generation

In this section we present the procedure used for the generation of turbulent
initial conditions used to study the collapse of turbulence at Reynolds number
$R$ (by means of DNS or AMS). This procedure uses mixing at a higher Reynolds
number $R_{+}=600$ to naturally decorrelate turbulent initial conditions. It
is easily parallelised with a minimal load imbalance. If we use $p$ threads,
we generate $N_{j}=\left\lceil\frac{N}{p}\right\rceil$ or
$N_{j}=\left\lfloor\frac{N}{p}\right\rfloor$ initial conditions on each thread
$j$. Each thread uses distinct seeds for random number generation.

  1. a)

On each thread $j$, $1\leq j\leq p$, we first create an artificial velocity
field $u_{y}=0$, $u_{z}=0$, $q=0$,
$u_{x}=0.4\sin\left(\pi\frac{y+1}{2}\right)\cos\left(4\pi\frac{z}{L_{z}}M_{z}\right)$,
with spanwise wavenumber
$M_{z}=\max\left(1,\left\lfloor\frac{L_{z}}{5}\right\rfloor\right)$. This
corresponds to streamwise velocity tubes which are prone to the streaks
instability, a key process in the cyclic SSP (Waleffe,, 1997) and should thus
lead to wall turbulence.

  2. b)

On top of the velocity components, we add a noise red in $x$ and $z$ and white
in $y$. This yields an initial velocity field $\mathbf{u}_{0}$. This red noise
is such that the variance of its Fourier mode for streamwise wavenumber
$n_{x}$ and spanwise wavenumber $n_{z}$ on component $l$ is
$\sigma_{l}\gamma_{n_{x}}\gamma_{n_{z}}$, with shape factors
$\gamma_{n_{m}}=1$ for $0\leq|n_{m}|\leq 6$,
$\gamma_{n_{m}}=\frac{6}{|n_{m}|}$ for $n_{m}>6$, with $m=x$ or $z$. We set
$\gamma_{n_{m}}=0$ if $n_{m}>28$ in the $L_{x}\times L_{z}=24\times 18$ system
and if $n_{m}>36$ in the $L_{x}\times L_{z}=36\times 27$ system, with $m=x$ or
$z$. We use the variances $\sigma_{x}=0.05$, $\sigma_{y}=0.0025$ and
$\sigma_{z}=0.015$.

  3. c)

This initial condition is evolved for $T_{0}=500$ in the $L_{x}\times
L_{z}=24\times 18$ system and $T_{0}=200$ in the $L_{x}\times L_{z}=36\times
27$ system at mixing Reynolds number $R_{+}=600$. It has been checked that
this duration was long enough so that natural buffer layer turbulence forms
(Pope,, 2001). If the kinetic energy of this velocity field is larger than
$0.03$, this yields the fields $\\{\mathbf{u}_{R_{+},1,0},q_{R_{+},1,0}\\}$,
otherwise we go back to step b) and generate a new $\mathbf{u}_{0}$ with a new
realisation of the red noise.

  4. d)

We then generate the $N_{j}$ initial conditions
$\\{\mathbf{u}_{n},q_{n}\\}_{0\leq n<N_{j}}$ in the following manner. We first
evolve the fields $\\{\mathbf{u}_{R_{+},1,n},q_{R_{+},1,n}\\}$ at $R_{+}$ for
$T_{+}=500$ ($L_{x}\times L_{z}=24\times 18$) or $T_{+}=200$ ($L_{x}\times
L_{z}=36\times 27$), yielding the fields
$\\{\mathbf{u}_{R_{+},2,n},q_{R_{+},2,n}\\}=\\{\mathbf{u}_{R_{+},1,n+1},q_{R_{+},1,n+1}\\}$,
which will be used to generate the initial condition at $R$ and generate a
subsequent decorrelated field at $R_{+}$.

  5. e)

We then set the Reynolds at $R$ (where we study collapse) and we set
$\\{\mathbf{u}_{R,1,n},q_{R,1,n}\\}=\\{\mathbf{u}_{R_{+},2,n},q_{R_{+},2,n}\\}$.
These velocity and pressure fields are evolved at $R$ during $T_{-}=750$
($L_{x}\times L_{z}=24\times 18$) or $T_{-}=500$ ($L_{x}\times L_{z}=36\times
27$). This duration is chosen so that enough mixing has occurred and each
initial condition is decorrelated from the others. We obtain
$\\{\mathbf{u}_{R,2,n},q_{R,2,n}\\}$. We then let it evolve until the kinetic
energy is either within $1.25\%$ of $E_{0}$, in which case we have our $n^{\rm
th}$ initial condition $\\{\mathbf{u}_{n},q_{n}\\}$ or is below $0.03$. In
that case we restart at d) by setting
$\\{\mathbf{u}_{R_{+},1,n},q_{R_{+},1,n}\\}=\\{\mathbf{u}_{R_{+},2,n},q_{R_{+},2,n}\\}$.
We use $E_{0}=0.055$ in the system of size $L_{x}\times L_{z}=24\times 18$ and
$E_{0}=0.052$ in the system of size $L_{x}\times L_{z}=36\times 27$.

This approach ensures that we have $N$ decorrelated initial conditions which
verify a given constraint on kinetic energy (for instance).

###  2.4 Reaction coordinates

Since the kinetic energy $E_{k}$ (Eq. 3) of the turbulent flow is fluctuating
around a conditional average, while the kinetic energy of the laminar flow is
zero (Fig. 3 (a)), a first choice to construct the reaction coordinate is to
use $E_{k}(t)$. The simplest reaction coordinate based on $E_{k}$ is affine.
We can therefore propose the reaction coordinate $\Phi_{E}$ defined directly
as a function of time by

$\Phi_{E}(t)=\frac{E_{\rm turb}-E_{k}(t)}{\Delta E}\,.$ (5)

In this affine function, we use the shift to the kinetic energy $E_{\rm turb}$
and the normalisation of the reaction coordinate $\Delta E$. These two
quantities are chosen so that $\Phi_{E}\simeq 0$ when the flow is turbulent
and $\Phi_{E}\simeq 1$ when the flow is laminar. It is natural to choose
$E_{\rm turb}$ close to some average of the kinetic energy and $\Delta
E\lesssim E_{\rm turb}$. In order to estimate $E_{\rm turb}$, one can sample
$E_{k}(t)$ and construct the empirical probability density function
conditioned on the flow experiencing no collapse of turbulence, in a system of
given size and Reynolds number. One can for instance retain the part of the
time series where $E_{k}(t)\geq 0.025$ (indicated by the black dashed line in
figure 3 (a)). Using this, one can compute a sample mean of the kinetic energy
$E$ (Fig. 3 (b)), and then choose close enough $E_{\rm turb}$ and $\Delta E$
accordingly. We will state what values of $E_{\rm turb}$ and $\Delta E$ are
chosen in AMS computations in section 3 and section 4. We also compute the
conditional variance $\sigma$ as a function of the Reynolds number for the two
system sizes (Fig. 3 (c)): we will use this quantity to choose the
hypersurface $\mathcal{C}$, so that we typically have $E_{\rm
turb}-E_{k}(t)\simeq\frac{\sigma}{2}$ on $\mathcal{C}$. This corresponds to a
short excursion where the kinetic energy is away from its conditional average
by typically half a standard deviation. Using the reaction coordinate, we set
$\mathcal{A}$ as all velocity fields such that $\phi_{E}\leq 0$ and
$\mathcal{B}$ as all velocity fields such that $\Phi_{E}\geq 1$.

One can note that the variance of the kinetic energy decreases with Reynolds
number in the range of interest (Fig. 3 (c)). The larger variance for $R\leq
380$ are caused by fluctuations of kinetic energy toward low values at lower
Reynolds numbers. The distributions of kinetic energy are asymmetric for
$330\leq R\leq 380$ with a relative skewness that can be below $-0.5$. As the
Reynolds number is increased, the excursions of $E_{k}$ toward low values
become less and less probable and the PDF of kinetic energy become narrower
and more symmetric (see (Rolland,, 2015) Fig. 9 (c), in a case where the
domain can accommodate laminar-turbulent bands). In domain large enough to
contain laminar-turbulent bands (Fig. 1, centre), one can observe a narrow
peak of the variance in the range $378\lesssim R\lesssim 382$ (at our
numerical resolution), because temporary closing of laminar holes can be
observed on timescales smaller than $\mathcal{O}\left(10^{4}\right)$. At
higher Reynolds numbers, the variance of kinetic energy increases again
because of turbulent fluctuations.

Figure 3: (a) Example of a time series of kinetic energy in a domain of size
$L_{x}\times L_{z}=24\times 18$ at Reynolds number $R=370$ The black dashed
lines indicate were the sampling is stopped for the construction of the
empirical probability density function of kinetic energy conditioned to not
collapse happening. (b) Conditional average of the kinetic energy as a
function of the Reynolds number for domains of size $L_{x}\times
L_{z}=24\times 18$, $L_{x}\times L_{z}=36\times 27$. (c) Conditional variance
of the kinetic energy as a function of the Reynolds number for domains of size
$L_{x}\times L_{z}=24\times 18$, $L_{x}\times L_{z}=36\times 27$.

In appendix A, we use an other reaction coordinate $\Phi_{\rm asym}$, based on
$\psi(t)$, the asymmetry of the streamwise velocity field with respect to the
midplane $y=0$

$\psi(t)=\frac{1}{V}\int_{x=0}^{L_{x}}\int_{y=-1}^{1}\int_{z=0}^{L_{z}}u_{x}{\rm
sgn}(y)\,{\rm d}x{\rm d}y{\rm d}z\,,$ (6)

where ${\rm sgn}(y)$ stands for the sign of $y$. The choice of this reaction
coordinate is motivated by the fact that we have mostly $u_{x}<0$ (coming from
low speed streaks) for $y>0$ and mostly $u_{x}>0$ (coming from high speed
streaks) for $y<0$ when velocity streaks (and thus wall turbulence) are
present in the flow (Jiménez and Moin,, 1991; Hamilton et al.,, 1995; Kawahara
et al.,, 2003). Note that with this definition $\psi$ is mostly negative. We
proceeded in the same manner as with the kinetic energy in order to construct
a reaction coordinate: conditional averages of $\psi$ are performed in order
to calculate $\Phi_{\rm asym}(t)=1-\frac{\psi(t)}{\langle\psi\rangle}$.

###  2.5 Direct numerical simulations of collapses

We perform two types of direct numerical simulations of turbulence collapse.
The first kind of simulation consists in letting the flow freely evolve from
initial conditions computed following the method of section 2.3. These
simulations are stopped when the flow has laminarised. This is deemed to have
happened when the reaction coordinate has reached one. In these simulations,
we only save time series, such as those of $E_{k}$, $E_{k,x}$, $E_{k,y-z}$ and
$\Phi_{E}$ (an example of such time series in given in figure 3 (a)). From the
last part of each DNS, we can extract the time series of these quantities in
natural collapse trajectories. We will perform two types of such DNS: some
free natural DNS, and some forced DNS, in order to check the effect of
additional noise on the collapse of turbulence. An additional noise will
necessarily be included in AMS computations, and we wish to know what is the
minimal error caused on the probability of collapse, trajectory features and
durations and mean first passage times by the addition of this noise. Since we
cannot add a comparable temporally localised perturbation in DNS, we will add
a noise red in space and white in time. In that case the forcing f (Eq. (1))
is exerted at all time and it is red in $x$ and $z$. The forcing noise is
characterised by its spatial correlation function

$\langle
f_{l}(\textbf{x},t)f_{m}(\textbf{x}^{\prime},t^{\prime})\rangle=\delta_{lm}\delta(t-t^{\prime})C_{l}(x-x^{\prime},z-z^{\prime})\tilde{C}_{l}(y-y^{\prime})\,,\,\hat{C}_{l}(n_{x},n_{z})=\Gamma_{l,n_{x},n_{z}}\,,$
(7)

with $C_{l}$, the correlation function of the noise for variables $x$ and $z$,
along component $l$, $\tilde{C}_{l}$, the correlation function of the noise
for variable $y$ along component $l$. The correlation function $C_{l}$ is
prescribed using its Fourier transform $\hat{C}_{l}$, with amplitude
$\Gamma_{l,n_{x},n_{z}}$, where $n_{x}$ stands for the Streamwise wavenumber,
$n_{z}$ stands for the spanwise wave number. In this text we always use
$\tilde{C}_{l}(y-y^{\prime})=\delta(y-y^{\prime})$ in numerical simulations
using a finite number $N_{y}$ of Chebyshev modes in the wall normal direction.
In DNS, we will always set the amplitude
$\Gamma_{i,n_{x},n_{z}}=\sigma_{l}\gamma_{n_{x}}\gamma_{n_{z}}$ using the
previously defined shape factors $\gamma_{n_{m}}$ for $m=x$ or $z$ (see §
2.3). We set the variance $\sigma_{l}=10^{-9}$ for all components.

Another kind of direct numerical simulation consists in repeating the first
stage of AMS, where we start the simulation from our initial conditions, let
them cross the hypersurface $\mathcal{C}$ and then either laminarise or go
back to turbulence as detected by $\Phi_{E}<0$. From these, we can compute the
proportion of trajectory that laminarises and thus have an unbiased estimate
of the probability of collapse and trajectory durations and validate AMS
estimates of these quantities.

##  3 System of size $L_{x}\times L_{z}=24\times 18$: reactive trajectories
and validation of AMS computations

In this section, we will compare the properties of the reactive trajectories
computed by means of AMS and by means of DNS as well as discuss the trajectory
properties in a system of size $L_{x}\times L_{z}=24\times 18$ at Reynolds
number $R=370$. We choose rectangular shaped domains: a shape comparable to
what is feasible in experiments, albeit at smaller size (see Fig. 1) (tilted
periodical domains being a numerical experiment). The values of the domain
size and Reynolds number imply that the study of turbulence collapse is
affordable by means of direct numerical simulations. Results of DNS and AMS
computations will be compared to assert what degree of trust can be placed in
general outputs of AMS computations. For estimates (of the probability of
collapse $\alpha$, the average duration of trajectories $\tau$, and the mean
first passage time $T$ and paths followed by trajectories) coming from both
types of computations, we will provide intervals of confidence and check
whether they overlap or if biases are present in the results of AMS
computations. From these we will be able to deduce what should be the error
bars _a minima_ that should be placed on the results of anticipated AMS
computations. This will be useful when no DNS are available for comparison,
for instance in larger domains at larger Reynolds numbers, in section 4.

In these computations, we will use $E_{\rm turb}=0.05$, $\Delta E=0.048$ for
the reaction coordinate $\Phi_{E}$ (Eq. (5)), in conjunction with the
conditional average of the kinetic energy (Fig. 3 (b)). The properties of the
initial conditions in $\mathcal{A}$ are given in section 2.3. The hypersurface
$\mathcal{C}$ will correspond to the set of velocity fields such that
$\Phi_{E}=0.06$, which roughly corresponds to half a standard deviation from
the conditional average (Fig. 3 (c)). We set the parameter $\xi=0.2$ with the
converging $\Phi_{b}$ (Eq. (13), § A, Fig. 2 (c)) for the anticipation of
branching. This corresponds to a good trade-off between the need for mixing
and separation of trajectories, and the minimisation of the number of retries
when branching trajectories. The AMS computations use $N=120$ clones and
suppress $N_{c}=32$ clones at each iteration. They are ran on 16 threads,
which leads to reasonable load balancing.

###  3.1 Collapse trajectories

####  3.1.1 Visualisation of the collapse trajectories

We first describe the velocity fields during the collapse, in trajectories
computed by AMS. With the use of AMS, we could produce far more time and space
resolved three components velocity fields than was previously shown when using
direct numerical simulations. Reactive and non-reactive trajectories are
simply reconstructed and saved at time intervals $\delta t>1$ from the
checkpoints used in AMS computations.

We display colour levels of the streamwise and spanwise components of the
velocity field in a midplane at successive instants during a collapse
trajectory in figure 4 and during one of the non-reactive trajectory remaining
at the last stage of AMS in figure 5. We will focus on the two fields
$(u_{x},u_{z})$ because of the physics of wall flows and experimental
constraints:

  * •

Indeed, on the one hand, it has been shown in models and observed in some
experiments and simulations that there was a distinct behaviour of velocity
components. One finds velocity streaks contained mostly in $u_{x}$ as well as
streamwise vortices, visible in the streamwise vorticity field
$\omega_{x}=\partial_{y}u_{z}-\partial_{z}u_{y}$, but that can be observed
through the component $u_{z}$ (they are the main contribution to this velocity
component) (Hamilton et al.,, 1995). Moreover it has been very recently
observed in models, simulations and experiments that decay of turbulence was
much faster in streamwise vortices (or its proxies) than in velocity streaks
(Rolland,, 2018; Liu et al.,, 2020; Gomé et al.,, 2020). We will endeavour to
discuss this observation in specifically computed collapse trajectories.

  * •

On the other hand, most experimental observations are performed by PIV using a
plane parallel to the wall, thus capturing $u_{x}$ and $u_{z}$ in a $x-z$
plane near the midgap (Liu et al.,, 2020; De Souza et al.,, 2020).

We first visualise the velocity fields as turbulence collapses. The reactive
trajectory starts from uniform buffer layer turbulence in the whole domain
(Fig. 4, $t=50$). The kinetic energy is of course close to the conditional
average. As time goes on ($t=174$), the kinetic energy decreases, and
observation indicates that the spanwise component of the velocity field is
much less intense than in the initial condition. Moreover, the largest values
of $u_{z}$ are spatially localised (intense small spatial scales for
$9\lesssim z\lesssim 14$, less intense, largest spatial scales for $0\lesssim
z\lesssim 9$). The streamwise velocity field has not yet decayed in amplitude,
however it is streamwise invariant where $u_{z}$ is almost $0$. As the flow
laminarises ($t=200$), the spanwise velocity field further decays while the
streamwise velocity field becomes even more streamwise invariant. The
amplitude of $u_{x}$ remains comparable to what was found in the initial
condition $|u_{x}|\lesssim 0.6$. This finally leads to a situation where
$u_{z}$ is negligible and only streamwise velocity tubes are left in the flow
($t=300$). These tubes then undergo viscous decay.

Figure 4: Time series of the kinetic energy (left panels) with a dot
indicating the point in time of each line ($t=50$, $t=174$, $t=200$, $t=400$),
alongside colour levels of the streamwise velocity (central panels) and
spanwise velocity (right panels), at corresponding times, in the midplane
$y=0$ during a collapse trajectory in a domain of size $L_{x}\times
L_{z}=24\times 18$ computed by AMS.

We can comparably examine a non-reactive trajectory, that is to say a
realisation of the dynamics that undergoes a large enough excursion of kinetic
energy to be retained at the last stage of AMS computations, but that still
goes back to a fully turbulent flow (Fig. 5). As with all other trajectories,
this one starts with turbulence in the whole domain (Fig. 5, $t=2$). We again
observe that $u_{z}$ decays faster than $u_{x}$ ($t=74$). In that case the
spanwise localisation of the decay is not disputable. Along with the decay of
$u_{z}$, the streamwise velocity component becomes streamwise invariant
($t=200$). The kinetic energy fluctuates for some time around a plateau
($150\lesssim t\lesssim 350$), then increases again. The values taken by
$u_{z}$ are getting more intense over an area of the domain that increases
($t=400$). Note that the streamwise velocity field had decayed locally during
the plateau (for $0\lesssim z\lesssim 7$). We observe an asymmetry between the
decay and the reinvasion processes. The component $u_{z}$ retracts to smaller
areas than $u_{x}$ during the decay: the active-quiescent interface (typically
planes $z=\text{cst}$) is not the same for $u_{x}$ and $u_{z}$. A strict
definition of these interfaces will be proposed in section 4 and figure 14.
Then, both $u_{z}$ and a streamwise dependent $u_{x}$ then fully reinvades the
domain ($t=600$), in that case the active-quiescent interface is the same for
both components. We thus observe a concomitant restart of both components of
the self sustaining process of turbulence.

Figure 5: Time series of the kinetic energy (left panels) with a dot
indicating the point in time of each line ($t=2$, $t=74$, $t=200$, $t=400$,
$t=600$), alongside colour levels of the streamwise velocity (central panels)
and spanwise velocity (right panels) in the midplane $y=0$ during a non-
reactive trajectory (hole opening then closing), at corresponding times, in a
domain of size $L_{x}\times L_{z}=24\times 18$ computed by AMS.

We can complement this view by observing the streamwise vorticity field
$\omega_{x{\rm,last}}=\partial_{y}u_{z{\rm,last}}-\partial_{z}u_{y{\rm,last}}$
leading to the largest value of the reaction coordinate in the last suppressed
state in the last stage of AMS (Fig. 6). In this system of relatively moderate
size, even if one always has an impression that turbulence collapses through
the formation of a laminar hole, the structure of this turning point state
varies from one run to another. One can observe clearly localised streamwise
vortices (Fig. 6 (a)), an almost entirely quiescent flow (Fig. 6 (b)) or
rather active streamwise vortices all over the domain (Fig. 6 (c)). One of the
reasons we have fairly different $\omega_{x{\rm,last}}$ is that they
correspond to different instants on similar trajectories that do display a
hole opening. This is also caused by the fact that the size of the domain is
not large enough for the reactive trajectories to be concentrated around a
typical trajectory.

In any case, we can have a first view of the path followed during the collapse
of turbulence. Of all the manner for the cyclic self sustaining process of
wall turbulence to fail (streaks first, vortices first, or concomitant decay
of streaks and vorticies), the observed scenario is that of a first decay of
streamwise vortices followed by a decay of the velocity streaks. No more
energy is extracted from the base flow so that the velocity streaks, turned
into tubes, slowly decay. They do not undergo any new streaks instability that
would refuel the streamwise vortices (Marquillie et al.,, 2011; Jiménez and
Moin,, 1991). We will examine all the reactive trajectory to show that far
from being some picked visualisations, this description is statistically
significant. We observed spanwise localisation of this failure of the self
sustaining process of turbulence. From the observations and monitoring of
turbulent fraction, one can note that there is more variability in this
respect. We will examine a larger system in section 4 to check whether this
observation is disputable or not.

Figure 6: Streamwise vorticity in the midplane $y=0$ in the system of size
$L_{x}\times L_{z}=24\times 18$ for last states at the last branching stage of
several AMS computations. (a) Run estimating $\hat{\alpha}=0.045$,
$\hat{T}=1.5\cdot 10^{4}$, (b) run estimating $\hat{\alpha}=0.086$,
$\hat{T}=1.1\cdot 10^{4}$, (c) run estimating $\hat{\alpha}=0.022$,
$\hat{T}=3.7\cdot 10^{4}$.

Figure 7: Normalised distribution of duration of collapse trajectories (in
logarithmic scale) sampled by AMS ($1590$ samples), by DNS ($189$ samples) and
compared to a normalised Gumbel distribution, in a system of size $L_{x}\times
L_{z}=24\times 18$ at Reynolds number $R=370$.

####  3.1.2 Comparison between trajectories computed by AMS and by DNS

We then start our comparison of the collapse trajectories computed by AMS and
DNS by describing the duration of collapse trajectories, that is to say the
time $s$ elapsed between the instant where the flow crossed the hypersurface
$\mathcal{C}$ and the instant where the reaction coordinate reaches $\Phi=1$.
We have obtained 1590 trajectories from AMS and $189$ trajectories form DNS.
We firstly compute the sample mean $\langle\tau\rangle_{o}$ and the sample
variance $\sigma_{\tau}$ over AMS realisations of the trajectory durations, we
obtain $\tau_{AMS}=413\pm 3$, $\tau_{DNS}=410\pm 7$ and
$\sigma_{\tau,AMS}=105$ and $\sigma_{\tau,DNS}=91$. As stated in section
2.2.3, the error bars on $\tau$ are given by $\sigma_{\tau}$ divided by the
square root of the number of samples555It is safe enough to assume each value
is decorrelated enough from the large majority of the others to assert that
the central limit theorem can be applied to the sampled mean: the most general
formulation only require short length correlations. There may be some degree
of correlations between durations of reactive trajectories that arise from the
same genealogy, but each genealogy is small enough compared to the sample
size.. We can see that both are within a few percent of one another and within
the $66\%$ confidence interval of each other. This shows that the reactive
trajectories computed by AMS and DNS have almost the same duration, and thus
that AMS has very little bias in the way it computes the collapse
trajectories. We then compare the sampled distributions of trajectory
durations. In figure 7, we display the distribution of normalised durations
$\tilde{s}=\frac{s-\tau}{\sigma_{\tau}}$ of collapse trajectories computed by
means of AMS and of DNS. Both are very similar and seem to be very close to a
Gumbel distribution. This distribution of duration has been reported to be
very close to a Gumbel in many instances of stochastically driven systems
(Rolland,, 2018; Rolland et al.,, 2016; Rolland and Simonnet,, 2015). This was
often related to a simple structure of the reactive trajectory that went
through a saddle of the deterministic part of the system. It is actually
demonstrated that reactive trajectories durations have this distribution in
one dimensional stochastic systems (Cérou et al.,, 2013). To our knowledge,
this is possibly the first time that such a distribution shape has been shown
in turbulence without stochastic injection of energy. This suggests us that
the manner in which turbulence collapses in systems of large enough size can
be compared to reactive trajectories in stochastic systems, in that it quite
possibly follows some effective dynamics agitated by background turbulence.
While it is more difficult to separate the dynamics into a deterministic part
and a stochastic part, this indicates that the reactive trajectories can
follow a rather simple path.

In order to investigate to what extent the dynamics could follow a stochastic
system-type transition path (see (Metzner et al.,, 2006)) when turbulence
collapses, we can then examine the trajectories themselves, computed by mean
of AMS and DNS. Since we obtained mostly time series from direct numerical
simulations, we will compare the path followed in the $(E_{k,x},E_{k,y-z})$
plane by the reactive trajectories (Fig 8). We produce a visualisation in a
small dimension space comparable to figure 6 of (Bouchet et al., 2019a, ).
Therefore, in order to examine the most probable path in this plane and how
far trajectories depart from these, we construct the probability density
functions $\rho_{e}$ in the $(E_{k,x},E_{k,y-z})$ plane using solely the
collapse trajectories computed in free Direct Numerical Simulations (Fig. 8
(a)), in noisy direct numerical simulations (Fig. 8 (b)) and AMS computations
(Fig. 8 (c)). We first note that the trajectories computed in all three cases
have very similar beginnings: we observe in all three cases a decrease of
$E_{k,y-z}$ much faster than that of $E_{k,x}$. While $E_{k,x}$ is divided by
$5$ (from approximately $0.05=\exp(-3)$ to $0.01\simeq\exp(-4.5)$),
$E_{k,y-z}$ is divided by $400$ (from approximately $0.0025\simeq\exp(-6)$ to
$5\cdot 10^{-6}\simeq\exp(-12)$). This difference in decay rate quantifies
what is observed in the velocity field: the streamwise vortices (whose
amplitudes mostly contribute to $u_{y}$ and $u_{z}$) decay before the velocity
streaks (whose amplitudes mostly contribute to $u_{x}$). We note that as
$E_{k,y-z}$ becomes smaller, the three PDF deviate in shape. In the case of
free DNS (Fig. 8 (a)), both $E_{k,x}$ and $E_{k,y-z}$ keep decreasing, albeit
at a smaller rate for $E_{k,y-z}$. However, in situations where some noise in
added to the flow (noisy DNS (Fig. 8 (b)) or AMS computations (Fig. 8 (c))),
one way or the other, the kinetic energy contained in the wall normal and
spanwise components $E_{k,y-z}$ reaches a neighbourhood of a minimal value
while $E_{k,x}$ keeps decreasing (the now streamwise invariant streaks keep
decaying in amplitude). This minimum value is the consequence of the forcing
exerted on the flow: $u_{y}$ and $u_{z}$ respond almost linearly to this
forcing and fluctuate around these small values. In the case of AMS, this is
the partially numerical effect of what is left in the spectrum of the white
perturbations placed at branching. The value around which the components
fluctuate is a function of the noise variance and spectrum shape. It has been
checked in an AMS study of the build up of turbulence that the amount of
energy given to the flow by this noise is small enough so that the probability
turbulence restarts from this forcing is negligible. The differences of
amplitude along the paths in theses PDF are (among other things) the
consequence of the number of bins used, the number of sampled trajectories and
the time spent by trajectories near low values of kinetic energy.

We can compare the trajectories more precisely by computing the conditional
average $\langle\log(E_{k,y-z})\rangle$ as a function of $\log(E_{k,x})$.
Indeed, if we define the probability $\rho_{x}$ of having the value
$\log(E_{k,x})$ during the reactive trajectories

$\rho_{x}(\log(E_{k,x}))=\int\rho_{e}(\log(E_{k,x}),\log(e_{k,y-z}))\,{\rm
d}\log(e_{k,y-z})\,,$ (8)

we have the conditional probability $\rho_{c}$ of observing
$\log(E_{k,y-z})=\log(e_{k,y-z})$ if one has $\log(E_{k,x})$

$\rho_{c}\left(\log(e_{k,y-z})\big{|}\log(E_{k,x})\right)=\frac{\rho_{e}(\log(E_{k,x}),\log(e_{k,y-z}))}{\rho_{x}(\log(E_{k,x}))}\,.$
(9)

and the conditional average is

$\langle\log(E_{k,y-z})\rangle(\log(E_{k,x}))=\int\rho_{c}\left(\log(e_{k,y-z})\big{|}\log(E_{k,x})\right)\log(e_{k,y-z})\,{\rm
d}\log(e_{k,y-z})\,.$ (10)

We can also compute the conditional variance

$\displaystyle\sigma_{y-z}(\log(E_{k,x}))=$
$\displaystyle\sqrt{\int\left(\rho_{c}(\log(e_{k,y-z}))|\log(E_{k,x})\log(e_{k,y-z})^{2}\right)\,{\rm
d}\log(e_{k,y-z})-\langle\log(E_{k,y-z})\rangle^{2}(\log(E_{k,x}))}\,.$ (11)

We display the average $\langle\log(E_{k,y-z})\rangle(\log(E_{k,x}))$ in full
line and the average plus/minus the variance
$\langle\log(E_{k,y-z})\rangle(\log(E_{k,x}))\pm\sigma_{y-z}(\log(E_{k,x}))$
in dashed line for all three types of reactive trajectories in figure 8 (d).
This quantitatively confirms the observation of a faster decay of the
streamwise vortices, and thus of $E_{k,y-z}$, than of the streaks, and thus of
$E_{k,x}$. This also confirms that using the AMS leads to the same collapse
trajectories as using DNS, at a much smaller cost, up to the effect of the
noise added at branching on $E_{k,y-z}$ when it is very small. We can finally
note that the trajectories seem to be concentrated around this average. Note
however that the constant variance of $\log(E_{k,y-z})$ actually corresponds
to a variance of $E_{k,y-z}$ which decreases exponentially as the collapse of
turbulence goes on. There is a larger variance of $E_{k,y-z}$ among the
initial conditions and thus in the amplitude and distribution of streamwise
vortices in the initial conditions and stronger fluctuations in the beginning
of the trajectory. On a final note, we state that in order to rigorously
obtain instantons in stochastic systems, one has to take a limit of the
variance of the added noise going to $0$. In the case of our turbulent flows,
this means either identifying an effective noise (which is a difficult task)
or displaying large deviations in quantities like the probability density
functions, the probability of collapse or the mean first passage time.
Quantitatively speaking, this amounts to identifying a small parameter
$\epsilon$, proportional to the square of the suspected noise variance such
that $\lim_{\epsilon\rightarrow 0}\epsilon\log(\alpha)$ or
$\lim_{\epsilon\rightarrow 0}\epsilon\log(T)$ are independent on $\epsilon$
(these two limits are often the opposite of one another) (Touchette,, 2009).
As it happens, such scalings have been displayed in probability density
functions of kinetic energy in simulations when the sizes $L_{x}$, $L_{z}$
were increased so that domains contained up to 6 wavelengths of the laminar
turbulent bands, with $\epsilon_{DNS}=\frac{1}{L_{x}L_{z}}$ (Rolland,, 2015).
In that case, the rate functions of the PDF of kinetic energy showed the tails
for $E_{k}<\left\langle E_{k}\right\rangle$, describing the fluctuations of
kinetic energy toward low values and eventually toward laminar flows. This
type of scaling was also shown in mean first passage times before collapse of
laminar-turbulent coexistence in models of pipe flow where the length $L$ of
the system was increased so that it contained up to ten puffs, with
$\epsilon=\frac{1}{L}$, (Rolland,, 2018). One thus observes large deviations
in the large system size limit, and in that limit one should expect to observe
narrower and narrower concentration of trajectories, and an increase of the
mean first passage time before collapse which is exponential with the length
of the system. Note that models indicate that nothing in laminar-turbulent
coexistence prevents the collapse of turbulence: if the Reynolds number is
increased slightly above the threshold of transition, this only requires the
concomitant collapse of all the puffs. This condition is at the origin of the
exponential dependence of the lifetime of turbulence on size. It may very well
be that the very same thing happens in actual pipe flows or Couette flows,
where one would require the rare, but not impossible, concomitant collapse of
all the puffs or of the bands of the flow. This synchronisation could be at
the origin of an exponential scaling in size. Note that in that respect, there
is a similar phenomenology in our rectangular boxes and larger domains: in
both case, the collapse of turbulence on a streamwise long corridor is
required, and in both cases, that collapse can be countered by contamination
from the sides of the hole. This would indicate that laminar turbulent
coexistent is strictly permanent only in domains of infinite size.

Given these considerations, we can conclude that in order to observe more
concentration of trajectories, and more generally large deviations of the
probabilities and rate of probability of collapse, it may very well be
necessary to observe said collapses in larger and larger domains. For this
purpose increasing the streamwise size $L_{x}$ may be more efficient than
increasing the spanwise size $L_{z}$. This may be a consequence of the
specific mechanism of hole formation, leading to the formation of two active-
quiescent fronts that move on. Such situations have already been observed in
theoretical physics systems (Rolland et al.,, 2016). The rare event is the
formation of the hole and of the fronts themselves. Once said hole is created,
the probability that it will entirely open is small and decreases with size,
but this decrease is slower than exponential.

(a) (b) (c)

Figure 8: Probability density functions in the $(E_{k,x},E_{k,y-z})$ plane,
conditioned on being in a collapse trajectory, for (a) Trajectories from free
DNS, (b) trajectories from noisy DNS, (c) 1138 trajectories from AMS. (d)
Conditional Average trajectories and their variances.

We finally note that even if each AMS run has some collapse trajectories in
its first iterations (approximately $3.8\%$ of the trajectories generated in
the first stage), and that it will generate other collapse trajectories by
branching on these collapse trajectories, it generates even more variability
than this through mutation/selection on other initial conditions, with the
help of anticipated branching. This means that each run generates distinct,
independent reactive trajectories, in particular trajectories with distinct
initial conditions. These trajectories are even more distinct in their later
stages. The amount of variability generated can be estimated by counting the
number of genealogies, that is to say the number of group of trajectories that
share the same starting point. Each group will display trajectories
decorrelated from those of another genealogy. On average, we generate 24
genealogies per AMS run with $120$ clones. This value is obtained by averaging
over all available AMS runs the number of distinct initial conditions in
computed reactive trajectories. See (Ragone and Bouchet,, 2020) for a view of
genealogies of trajectories with a different algorithm applied to climate.

###  3.2 Mean first passage time and probability of collapse

We now compare the mean first passage time before collapse $T$ and probability
$\alpha$ of collapse computed by means of direct numerical simulations and
AMS. While the mean first passage time before collapse is the most physical
quantity, the probability of collapse nevertheless contains a lot of
information on the precision of the computations and the rarity of the event.
Looking into its statistics is always a good way of validating the AMS
computations. This is all the more true because this is the one quantity for
which we have mathematical results of convergence (Cérou and Guyader,, 2007;
Bréhier et al.,, 2016; Cérou et al., 2019a, ). Using the methods presented in
sections 2.2 and 2.5, we obtain the averages with the 66% confidence intervals
$\langle\alpha\rangle_{o}=\alpha_{AMS}=0.042\pm 0.003$ with a variance
$\sigma_{\alpha,AMS}\simeq 0.016$ using AMS, and $\alpha_{DNS}=0.037\pm 0.002$
using direct numerical simulations. This indicates that the estimation of the
probability of collapse is precise, we have an overlap of both 66% interval of
confidence (as presented in section 2.2.3). When we suppress a fixed number
(strictly larger than one) of clones at each iteration, the ideal variance of
the probability of collapse is given by
$\sigma_{id}^{2}=\frac{\alpha^{2}}{N}\left(\langle\kappa\rangle_{o}\frac{N_{c}}{N-N_{c}}+\frac{1-\langle
r\rangle_{o}}{\langle r\rangle_{o}}\right)$ (where $\langle\kappa\rangle_{o}$
is the average number of iterations of AMS and $\langle r\rangle_{o}$ is the
average proportion of reactive trajectories computed at the last iteration)
(Cérou and Guyader,, 2007). This is the minimum of variance of the estimator
of $\alpha$. It is obtained when the ideal reaction coordinate, the committor
is used (Bréhier et al.,, 2016; Cérou et al., 2019a, ). When this minimal
variance is reached, AMS computations are performed with a minimal error. In
practice this gives a value of the ideal variance $\sigma_{AMS,id}\simeq
0.007$. The variance on the estimate of the probability is twice the ideal
variance. This indicates that the quality of the computation is acceptable,
but not perfect, quite possibly due to imperfections in the reaction
coordinate we used.

We then compare the mean first passage times before collapse. By means of DNS,
we sample the cumulated density of waiting times before collapse in
simulations performed with and without additional noise, in order to estimate
the mean first passage time as well as test the effect of noise on it. We can
see that cumulated densities $F(t)=\int_{t}^{\infty}f(\zeta)\,{\rm d}\zeta$,
with $f(\zeta)$ the PDF of passage times, computed with and without noise
decrease linearly in logarithmic scale (Fig. 9). This is thus entirely
consistent with an exponential distribution of passage times. Both cumulated
density functions have very similar slopes, which is itself almost equal to
the estimated mean first passage time. We estimate the mean first passage time
before collapse in DNS by averaging over all sampled durations, we obtain
$T_{DNS}=1.10\cdot 10^{4}\pm 8\cdot 10^{2}$ without noise and
$T_{DNS,noisy}=1.3\cdot 10^{4}\pm 10^{3}$ with additional noise, where the
error bars correspond to the 66% confidence interval. In both cases, the
variance is equal to the average and the median is equal to $\sqrt{2}T$,
further confirming that we have an exponential distribution. While the
addition of noise may increase the mean first passage time, it does not
increase it dramatically, so that it is acceptable to use an additive small
noise in simulations and AMS to help the separation of trajectories. In turn,
we estimate $T_{AMS}=2.4\cdot 10^{4}\pm 3\cdot 10^{3}$. Again, we used the 66%
interval of confidence as error bars. We note that $T_{AMS}$ over estimates
the mean first passage time. The confidence intervals do not overlap, which
indicates that this overestimate is most likely an effect of a bias in the
computation of $T_{AMS}$. In order to understand why the estimate of the mean
first passage time is not as good as what is obtained for other quantities, we
examine the histograms of mean first passage times (Fig. 10 (a)) and
probability of collapse (Fig. 10 (b)) computed by AMS. While the histogram of
$\alpha$ is symmetric, that of $T$ is skewed. We note that the estimate of $T$
is polluted by a handful of very high values, originating from fairly small
values of $\alpha$, which have a much smaller effect on the average over
realisations $\alpha_{AMS}$. On top of the improvement brought by a better
reaction coordinate, the computation of the mean first passage time may very
well be improved in AMS computations using more clones.

Figure 9: Logarithm of the cumulated distribution of waiting times in
noiseless and noisy DNS for $L_{x}\times L_{z}=24\times 18$, $R=370$. We add
the two linear functions $-t/\langle T\rangle$ for comparison.

Figure 10: Histograms of $T$ (a) and $\alpha$ (b) sampled by AMS in the system
of size $L_{x}\times L_{z}=24\times 18$, at Reynolds number $R=370$.

One can finally compare the relative cost of trajectory computations with AMS
and DNS. Strictly speaking, the average physical duration of simulation to
obtain one reactive trajectory by DNS is of order $\mathcal{O}(10^{4})$,
meanwhile the simulation time by means of AMS is of order
$\mathcal{O}(10^{3})$. Such a factor $10$ of acceleration of computation
already seems interesting. This means that even in this case where the study
of the system by DNS seems affordable, using AMS already leads to a
substantial acceleration of computations. Of course, with the use of
checkpoints to store the trajectories during AMS computations, the reactive
trajectories and non-reactive trajectories can be recorded without unnecessary
use of RAM and disk-space. Indeed, in many cases, multistability (or more
precisely collapse of turbulence) is studied using long records of time
series, where most of the information has to be eventually discarded. Using
such records to compute the velocity fields during collapse would be very
inefficient. Velocity fields presenting collapse are scarcely presented in DNS
or experimental studies. This does not have to be a consequence of the
numerical procedure, as one can argue that only records spanning the last one
thousand time units or so of the velocity have to be kept. This requires
careful memory management and has very seldom been explicitly presented in the
literature to the author’s knowledge.

##  4 A very rare collapse in a domain of size $L_{x}\times L_{z}=36\times
27$: spatial organisation

The observation of collapse of turbulence in the domain of size $L_{x}\times
L_{z}=24\times 18$ indicated that, locally, collapse occurred by a more rapid
failure of streamwise vortices than velocity streaks. In that domain, however,
it was not entirely clear whether this collapse was always a local process,
leading to a streamwise laminar hole that then extended in the spanwise
direction, or whether this could happen in a synchronous manner everywhere.
The question is open as both behaviours have been observed in a finite number
of decays forced by decrease of the Reynolds number to very low value
(Manneville,, 2011; De Souza et al.,, 2020). We can therefore ask the question
of whether we would always see both types of collapse when the domain size is
increased and the simulation is closer to what can be observed in experiments.
For this reason we will consider the collapse of turbulence in a domain of
size $L_{x}\times L_{z}=36\times 27$ at Reynolds number $R=377$ (see Fig. 1
for an illustration of the domain size). This will also be the occasion to
test anticipated AMS in a situation where the mean first passage time before
collapse is much larger than in section 3.2. In that case direct numerical
simulations are not affordable and there will be no collapse trajectory
computed at stage $0$ of AMS. This will indicate how expensive an AMS
computation is in that situation and whether anticipated AMS manages to create
collapse trajectories and ensures that there is some variability between these
trajectories.

For the AMS computations, we use $E_{\rm turb}=0.05$, $\Delta E=0.047$ in
$\Phi_{E}$ (Eq. (5)). We set the hypersurface $\mathcal{C}$ as the set of
velocity fields such that $\Phi=0.07$. The initial conditions within
$\mathcal{A}$ are prepared as stated in section 2.3. We use the parameter
$\xi=0.15$, with the converging $\Phi_{b}$ (Eq. (13), § A) for the
anticipation of branching. In each computation, we use $60$ clones and
suppressed 16 clones at each iterations. On 8 threads, this gives reasonable
load balancing for the simulations. Due to the rarity of collapse of
turbulence for these control parameters, there is no collapse trajectory in
the very first stage of AMS computation: the method selects them naturally. A
downside is that there is less variability between the trajectories computed
in the same AMS run, in their first 100 time units. There is of course
variability from one AMS run to the other. We note that the larger the domain
is, the more decorrelated fluctuations are adding up in the flow and the less
need there is for anticipation of branching. With these parameters, we
estimate $\alpha=8\cdot 10^{-5}\pm 7\cdot 10^{-5}$ and $T=8\cdot 10^{7}\pm
7\cdot 10^{7}$, where the error bars are given by the 66% interval of
confidence. We have to simulate the flow for approximately $5000$ time units
per trajectory obtained. Even if one takes into account the fact that mean
first passage times are somewhat overestimated with a finite number of clones,
this still represents an acceleration of computations by a factor of more than
one thousand.

We display the last state at the last branching stage $\mathbf{u}_{\rm
last},\omega_{x{\rm,last}}$ from one of the AMS computations in figure 11.
This state (which is not steady in any way) displays a laminar hole localised
in $z$ which occupies the whole streamwise length of the domain. For the
streamwise component of the flow (Fig. 11 a,e), this hole is rather narrow and
is situated for $17\lesssim z\lesssim 24$. It is flanked by streamwise
invariant streaks for $8\lesssim z\lesssim 17$ and $24\lesssim 27$, $0\lesssim
3$ (due to periodic boundary conditions in $z$). In the wall normal component
(Fig. 11 (b)), the spanwise component (Fig. 11 (c)) and subsequently the
streamwise vorticity (Fig. 11 (d)), the laminar hole is actually much wider.
These components are away from zero only for $3\lesssim z\lesssim 8$, that is
to say where the streamwise component displays fluctuations. This is the only
area left in the flow where the self sustaining process of turbulence is still
active. This shift in the active-quiescent interface between the streaks and
streamwise vortices flow components highlights again the scenario of collapse
where the vortices disappear before the streaks. Unlike in the domain of size
$L_{x}\times L_{z}=24\times 18$, all the last state at the last branching
stage that have been computed display such an opening of hole. This has also
been observed at other Reynolds numbers ($R=351$, $357$ and $370$, not shown
here). The observation of this state confirms that for this regime of
parameters, the scenario of collapse of turbulence is far more univocal. As we
will see in the observation of collapse and partial collapse trajectories,
there is again a local failure of streamwise vortices that disappear in a
streamwise long hole leaving streamwise invariant streaks that decay
viscously. These holes then extend in the streamwise direction.

Figure 11: Colour levels of the streamwise (a), wall normal (b) and spanwise
(c) components of the velocity field, in the $y=0$,
$\vec{e}_{x}\times\vec{e}_{z}$ plane for the last stage at the last branching
step in an AMS computation of collapse trajectories at $R=377$ in a domain of
size $L_{x}\times L_{z}=36\times$, along with (d) the streamwise vorticity.
(e) Colour levels of streamwise averaged streamwise velocity as a function of
$y$ and $z$ for the same state.

We display a typical collapse trajectory obtained using AMS in figure 12. We
follow the decay of the kinetic energy (left panels) along with the streamwise
velocity and spanwise velocity in the midplane (central and right panels). The
trajectory starts from typical uniform wall turbulence at $t=2$, where the
streamwise velocity is organised in low speed and high speed streaks, and the
spanwise velocity varies at much shorter lengthscales, as would be expected
from its participation in spanwise vortices. After some slight decay in
amplitude of the spanwise velocity component from $t=2$ to $t=292$, a laminar
hole forms in $u_{z}$ from $t=292$ to $t=440$ as seen in the decay of kinetic
energy. In that case, this hole is located in $14\lesssim z\lesssim 25$ for
all $x$. As to $u_{x}$, the velocity streaks are still present at that
location, albeit almost streamwise invariant. At $t=624$, this hole has
extended in the spanwise direction. The velocity streaks are now streamwise
invariant for $z\gtrsim 13$ and decay in amplitude in this part of the flow.
Meanwhile the spanwise component of the flow has fallen to zero at the center
of the laminar hole. As time moves forward, at $t=876$, $u_{x}$ totally falls
to zero for $27\gtrsim z\gtrsim 19$, one can notice that the interface between
the active and quiescent domain has moved in the spanwise direction from
$z\simeq 16$ to $z\simeq 19$. Observation indicates that since the active-
quiescent interface is the same for $u_{z}$ and $u_{x}$ in that case. After
this event, one nevertheless can note that at $t=824$ the amplitude of the
spanwise velocity field has decreased in the turbulent domain indicating a
failure of the streamwise vortices in the whole $3\lesssim z\lesssim 19$ area.
The remaining streamwise vortices then completely decay, leaving only totally
streamwise invariant streaks localised in $z$ at $t=1076$, that then viscously
decay.

Figure 12: Following a collapse trajectory in a domain of size $L_{x}\times
L_{z}=36\times 27$, $R=377$. The left panels show a time series of kinetic
energy, with the dot indicating the instant at which the colour levels of
streamwise velocity in the $y=0$ plane (middle panels) and spanwise velocity
(right panels) at corresponding times in the $y=0$ plane are shown. _A movie
of this collapse trajectory is provided in the supplementary material_.

We display one of the $n<N_{c}$ non-reactive trajectories left at the last
stage of an AMS computation in figure 13. This trajectory corresponds to an
excursion toward laminar flow that evolves back to a fully turbulent flow. We
follow the excursion toward low values then increase of the kinetic energy
(left panels) along with the streamwise velocity and spanwise velocity in the
midplane (central and right panels). This trajectory starts like the one of
figure 12, by the opening of a laminar hole, first by decay of $u_{z}$ (for
all $x$, for $20\lesssim z\lesssim 27$ at $t=352$), leaving almost streamwise
invariant streaks at $t=432$. Part of this similarity may be explained by the
fact that the two trajectories share a common beginning, being branched from
the same genealogy. However they quickly separate, while displaying common
features of laminar hole opening, and of simplification of the streak preceded
by failure of the streamwise vortices. At $t=806$, there is a laminar hole in
$u_{x}$ for all $x$ and $17\lesssim z\lesssim 22$ and the vortices have almost
decayed. At $t=1076$ the streaks have further receded, however, the self
sustaining process of turbulence has restarted in the middle of the turbulent
domain, as seen by the increase of amplitude of $u_{z}$ in the $2\lesssim
z\lesssim 10$. From this point on, turbulence reinvades the domain. We can
note that at $t=1370$, the turbulent hole is closing, as seen in the surge of
kinetic energy. However the opening and the closing are different processes:
now the active-quiescent interfaces are almost at the same position in $z$ for
both $u_{x}$ and $u_{z}$. In the end, turbulence reinvades the whole domain.

We can observe the movement of the active-quiescent interfaces in more detail
for the reactive trajectory (Fig. 14 (a)) and the non-reactive trajectory
(Fig. 14 (b)). These figures use both the contours of
$\sqrt{\int_{x=0}^{36}u_{z}^{2}\,{\rm d}x-\left(\int_{x=0}^{36}u_{z}\,{\rm
d}x\right)^{2}}=0.03$ and $\sqrt{\int_{x=0}^{36}u_{x}^{2}\,{\rm d}x}=0.15$
which respectively indicate the spanwise location of active-quiescent
interface for the streamwise vortices and for the velocity streaks. These can
be compared to earlier visualisations (Fig. 12 and Fig. 13), that justified
the study of the purely spanwise location of the active quiescent interface in
this type of domain. In the case of the collapse trajectory, we can observe
the survival of streaks where streamwise vortices have disappeared for
$500\lesssim t\lesssim 800$ and $20\lesssim z\lesssim 27$ and for $900\lesssim
t\lesssim 1100$ and $5\lesssim z\lesssim 20$ (Fig. 14 (a)). While this is
partially caused by the method of detection of the interface (see Fig. 12
$t=624$) and while the scenario is not at one hundred percent univocal, the
situation where streamwise vortices survive where streaks have decayed is
never observed in the reactive and non-reactive trajectories. Similar
observations can be made in the case of the non-reactive trajectory (Fig. 14
(b)). We have an apparent collapse of both streaks and vortices at $t=400$ and
for $15\lesssim z\lesssim 20$. For $500\lesssim t\lesssim 800$, only the
active-quiescent interface of the vortices recedes. For $800\lesssim 1200$, we
can see some reduction of the surface occupied by the streaks. For $t\gtrsim
1200$, we can see that the surface occupied by the vortices increases again
until the active-quiescent interface for the streaks and the vortices is again
the same and they both reinvade the domain.

Figure 13: Following a non-reactive trajectory (hole opening then closing) in
a domain of size $L_{x}\times L_{z}=36\times 27$, $R=377$. The left panels
show a time series of kinetic energy, with the dot indicating the instant at
which the colour levels of streamwise velocity in the $y=0$ plane (middle
panels) and spanwise velocity (right panels) at corresponding times in the
$y=0$ plane are shown. _A movie of this non-reactive trajectory is provided in
the supplementary material_.

Figure 14: Contours of $\sqrt{\int_{x=0}^{36}u_{z}^{2}\,{\rm
d}x-\left(\int_{x=0}^{36}u_{z}\,{\rm d}x\right)^{2}}=0.03$ and
$\sqrt{\int_{x=0}^{36}u_{x}^{2}\,{\rm d}x}=0.15$ as a function of $z$ and $t$
for (a) the collapse trajectory of figure 12 and (b) the non-reactive
trajectory of figure 13.

##  5 Conclusion

In this work we have presented one of the first applications of a rare event
sampling method to study multistability in a turbulent flow not forced
stochastically. Namely, Adaptive Multilevel Splitting was used with
anticipated branching to compute and study turbulence collapse trajectories in
transitional plane Couette flow. A large number of collapse trajectories could
thus be computed using a dramatically reduced amount of computational time. In
some of the cases considered here, simulations can be more than one thousand
times faster. Calculations may be even more accelerated for rarer events as
the computational cost of AMS typically evolves like the logarithm of the
computational time of DNS. Moreover, the collapse trajectories can be
straightforwardly recorded in full details using the checkpoints created
during AMS computations. The reactive trajectories in phase space, statistics
of reactive trajectories durations, velocity fields during reactive and non-
reactive trajectories as well as estimates of the probability of turbulence
collapse and mean first passage time before turbulence collapse (with their
respective error bars) are then available at low cost using this approach.

A first step in this work has been to ensure the validity of the result of AMS
computations. This has been done in a system of size $L_{x}\times
L_{z}=24\times 18$ at Reynolds number $R=370$. Computation of reactive
trajectories and mean first passage time are amenable in this system, so that
rare event computations could be compared to a reference. It could thus be
shown that reactive trajectories (trajectories in phase space, statistics of
durations) are computed precisely. Some effects of the low noise added on the
flow for the separation of trajectories were identified in the very last
stages of reactive trajectories, when the flow is finally laminar. The bias
added to the collapse computations is small enough so that the use of this
additional noise is entirely acceptable. Similarly, the probability of
collapse computed using AMS was shown to be in the same 66% interval of
confidence as the one computed by DNS. The computed mean first passage times
are overestimated. Histograms of the computed probabilities of collapse and
mean first passage time before collapse indicate that this overestimate is
caused by a handful of overly large estimates of $T$ related to underestimates
of $\alpha$. This is most likely an effect of the number of clones used in
this computation. Other estimators of mean first passage time could be used to
improve the computation (Lestang et al.,, 2020), which have their own
drawbacks for trajectory computations, such as constraints on trajectory
durations.

A large part of this work was of course devoted to the study of the collapse
trajectories properties. In the system of size $L_{x}\times L_{z}=24\times
18$, the visualisation of collapse trajectories showed that the laminarisation
could be sketched in two stages. The velocity components $u_{y}$ and $u_{z}$,
and as a consequence, streamwise vortices, decayed first, so that the flow is
driven out of the self sustaining process of turbulence (Waleffe,, 1997). As a
consequence, streamwise velocity has less dependence on $x$, and then slowly
decays without displaying streaks instability (Marquillie et al.,, 2011;
Kawahara et al.,, 2003). Examination of the velocity fields shows that this
process can happen locally, leading to a streamwise long laminar hole in
$u_{y}$, $u_{z}$ and then extends in $z$ : streamwise vortices decay
everywhere in this small domain. Viewed in the quadrant $E_{k,x}\geq
0,E_{k,y-z}\geq 0$, these visualisations all correspond to trajectories
concentrated along the same path. While this path is more continuous than the
sketch given here, it still displays a first stage where $E_{k,y-z}$ decays
much faster than $E_{k,x}$, then a second stage where $E_{k,y-z}$ is small
(around the noise level) and $E_{k,x}$ decays in turn. Such concentration of
trajectories are commonplace in stochastic models which display rare
transitions when the noise variance goes to zero. While the separation of the
dynamics of an actually turbulent flow into a deterministic part and a noise
part is very difficult to make, and that the flow displays transient chaos
when the domain is very small (Eckhardt et al.,, 2007), this is a first
indication that the collapse of turbulence in Couette or Poiseuille flows of
increasing size is best understood as a stochastic system. Given the
concentration of trajectories, one may be tempted to draw a comparison to
Instantons computed in stochastic systems (Bouchet et al., 2019a, ; Grafke and
Vanden-Eijnden,, 2019). In order to give a meaning to such a description, one
would need a low noise limit. This can possibly be the case in the large size
limits as observed in models (Rolland,, 2018) and in the large deviations of
the kinetic energy (Rolland,, 2015). The statistics of reactive trajectories
durations are another feature of the trajectories which are very close to
their stochastic counterparts. Both DNS and AMS computations show that the
normalised distribution of durations is very close to a Gumbel distribution.
In stochastic multistable systems it is demonstrated in one dimensional
systems that this is the distribution of durations (Cérou et al.,, 2013). In
extended systems, it is observed that this distribution is found when the
multistability occurs when the system transits through a single curved saddle
of its deterministic part (Rolland et al.,, 2016; Rolland,, 2018). This would
place the collapse of turbulence in plane Couette flow in systems of
increasing sizes in this category.

The second part of the study of trajectories concerned the question of the
localisation of the laminar holes. For this matter the size of the studied
system has been increased, from $L_{x}\times L_{z}=24\times 18$ to $36\times
27$. Increasing the domain size led to the computation of a single type of
reactive trajectories. They all display the local opening of a spanwise
localised streamwise elongated hole, first in streamwise vorticity, then in
streamwise velocity as this hole extends in the spanwise direction.
Laminar–turbulent front movements have been presented by (Duguet et al.,,
2011). There has been some observations of similar hole opening, as well as
global collapse in numerical simulations and in experiments (Manneville,,
2011; De Souza et al.,, 2020). However, there has been so far no general way
to ascertain the weight of each scenario in turbulence collapse for given
values of control parameters. Using a method like AMS and sampling a large
number of trajectories at a comparably low cost can therefore help give some
statistical weight to each scenario of collapse. Moreover, this suppressed the
need to perform a quench or a more continuous decrease of the Reynolds number
in order to observe said collapses. This reopens the way for the study of
collapse of non isolated turbulent puffs and spots, which had been considered
experimentally (Bottin and Chaté,, 1998) and more recently in models. The
dynamics of couette flow, where turbulence can extend in two direction, are
even richer than in the case of the collapse of laminar-turbulence coexistence
in models of the canonical pipe flow puff (Rolland,, 2018).

Another interesting point of the use of a method like AMS is the ability to
identify a turning point in the dynamics of collapse (Simonnet et al.,, 2021).
This is done by recording the state leading to the largest reaction coordinate
in the suppressed trajectories of the last stage of AMS. In stochastically
forced systems, this state corresponds to an effective saddle between the two
metastable states, and can alternatively be computed through a dichotomy
procedure (see (Willis and Kerswell,, 2009; Schneider et al.,, 2007)), or even
by computing the saddle point of the deterministic point of the system in
systems forced by an additive noise. This turning point gives us some
information on the mechanism of multistability and in our case on collapse of
turbulence. It displayed localised streaks and streamwise vortices in the
largest system, which backs the scenario of collapse through hole formation.
In the case of collapse of turbulence in plane Couette flow, these states did
not necessarily make the best starting points for edge tracking. Edge tracking
has already been performed with the results of collapse of turbulence in pipe
flow (De Lozar et al.,, 2012). Note that edge tracking can be performed using
the result of AMS computation: this has been done in the study of build up of
turbulence under some additive noise and should be presented in a later text.

While the use of AMS to compute the reactive trajectories has proven
successful, there is still some room for improvement, in particular for the
estimation of the mean first passage time before collapse. With a given set
up, one can always bring improvement on the estimate by increasing the number
of clones used to compute the trajectories (and thus their numerical cost).
This has been repeatedly measured in stochastically driven systems (Rolland
and Simonnet,, 2015) and demonstrated for such systems asymptotically (Bréhier
et al.,, 2016; Cérou et al., 2019a, ). However, one could wish to perform
these improvements at little additional numerical cost. Said studies have also
shown that much improvement can be brought by improving the reaction
coordinate used in computations (Rolland and Simonnet,, 2015; Bréhier et al.,,
2016; Rolland et al.,, 2016; Cérou et al., 2019a, ). This has been so far
performed by hand with trial and error by integrating physical properties of
the system in the reaction coordinate in order to mimic the committor
function: the probability of reach the arrival state before the departure
state from any given state (Onsager,, 1938; Bréhier et al.,, 2016; Rolland and
Simonnet,, 2015; Cérou et al., 2019a, ). Indeed it has been demonstrated that
it was the ideal reaction coordinate: the one that lead to the best estimates
(Cérou et al., 2019a, ). However, this committor function is in a way the
answer to the question we ask when we study multistability. Systematic methods
to approximate the committor and improve the reaction coordinate are being
validated in few degrees of freedom systems. They could be directly applied to
transitional turbulence. Indeed, the requirement for these methods to function
are met: trajectories, reactive or not, are reasonably well estimated. These
approaches may provide another estimate of the mean first passage time. This
could improve the estimates of $T$ provided in this text and avoid the
accumulation of error that occur when equation (4) is used.

Finally, the computation of reactive trajectories in transitional wall
turbulence opens the way for similar computations in high Reynolds number
turbulence. Systematic computation of reactive trajectories in systems of
aerodynamical interest (Kim and Durbin,, 1988), or geophysical interest are
more than conceivable (Herbert et al.,, 2020).

## Acknowledgements

The author thanks F. Bouchet, C.-E. Bréhier, T. Lestang and E. Simonnet as
well as T. Liu and B. Semin for interesting discussions. The author thanks the
help of the PSMN platform of ENS de Lyon, where most of the computation of
this work have been performed. The author finally thanks the hospitality of
Institut PPrime and ENSMA and their computational resources, where this work
has been started, as well as CICADA (now AZZURA, centre de calcul interactif
de l’Université de Nice). The author also thanks Felipe Alves Portela for help
proofreading the manuscript. The author also acknowledges the comments of the
anonymous referees that helped improve the manuscript.

##  Appendix A Anticipated AMS

The various anticipated branching strategies proposed in section 2.2 help
increasing the variability amongst trajectories. More importantly this
approach helps avoiding so called extinction. Indeed, when using AMS, we wish
to avoid a situation where no trajectories are able to go further than
$\Phi_{\rm ext}<\Phi(\partial\mathcal{B})$. In that case, no matter how many
additional steps are performed, we have at stage $k$
$\Phi_{\max,k}=\max_{i}\max_{t}(\Phi_{i}(t))=\Phi_{\rm ext}$ and reactive
trajectories cannot be computed. We illustrate this situation in figure 15
where we display the time series of the reaction coordinate $\Phi$ for all
trajectories successively computed by AMS in a run where no anticipation is
performed. The reaction coordinate used here is based on the asymmetry of the
flow (see Eq. 6). The figure was obtained in a sequential run of AMS where all
the time series of the reaction coordinate for all the trajectories simulated
were successively saved. We draw a distinction between the trajectories
computed in the initial step of the algorithm (times before the black vertical
line) and those computed during the mutation selection process (times after
the black line). During the initial free run stage of AMS (“stage 0”), the
maximum of reaction coordinate over all trajectories is
$\Phi_{\max,0}=0.248\pm 0.001$. Note that this maximum of $\Phi$ on a
trajectory corresponds to a turning point where $\Phi(t)$ will decrease with
an extremely high probability even if slightly perturbed. In the first steps
of the classical AMS computation, the levels of reaction coordinate at which
trajectories are branched are rather below $\Phi_{\max,0}$, so that
variability among trajectories is initially created and the maximum of $\Phi$
eventually reaches $\Phi_{\max,k}=0.430\pm 0.001$ after $k$ stages. Extinction
occurs because branching are performed at higher and higher levels of $\Phi$
which are closer and closer to the value $\Phi_{\max,k}$. Eventually all
trajectories reach that value, all branching is performed at $\Phi_{\max,k}$
and even with small perturbation, all trajectories subsequently decay in
$\Phi$. Closer examination shows that they all differ from one another after
some time. The issue here is not that trajectories do not separate from one
another, but that they do so during a decay phase of $\Phi$ and thus cannot
lead to a progress of the reaction coordinate. This has also been seen in the
study of turbulent wakes (Lestang et al.,, 2020). Anticipation of branching
consists in branching trajectories before these turning points so that they
all have the time to separate and generate new trajectories going further and
further in $\Phi$.

The anticipation strategy which should be followed depends on the properties
of the system. In the case of the collapse of turbulence in plane Couette
flow, this depends on the size of the system. The line is drawn between small
systems (typically of size $L_{x}\lesssim 15$, and $L_{z}\lesssim 10$) and
larger systems. Small systems which have very few degrees of freedom, display
at most four or five velocity streaks (high or low velocity) and a handful of
streamwise vortices. Their behaviour can be described by the escape away from
a chaotic repeller (Eckhardt et al.,, 2008). They necessarily collapse
globally in space, this means that the larger $\Phi$ is, the more globally
quiescent the flow becomes: the mixing is less intense and there is little
separation of trajectories. This is shown in figure 16 in a small system of
size $L_{x}\times L_{z}=12\times 8$. The flow required some reorganisation
before collapse was possible (part of the trajectory for $0\leq t\leq 300$),
and at most two low speed and two high speed streaks are present. When the
collapse takes place (between $t=350$ and $t=400$) $u_{z}$ collapses almost
uniformly. Larger systems have more and more degrees of freedom, and contain
more velocity streaks and vortices. An important feature is that they collapse
through the formation of a laminar hole. This means that for a large portion
of the reactive trajectory some area of the flow remains turbulent. This
turbulent area can serve as some sort of heat bath which fuels the trajectory
with perturbations and thus can greatly help the separation of trajectories.
This means that in small systems, branching must indeed be performed on
trajectories having greater excursions, but on their initial stages, where the
flow is still agitated and thus will separate the trajectories. In that case,
we use a so called saturated anticipation

$\Phi_{b,sat}=\xi\left(1-\exp\left(-\frac{\Phi_{N_{c}}}{\xi}\right)\right)$
(12)

this function is shown with $\xi=0.25$ in figure 2 (c), with the “saturated”
label. If we follow that strategy in the case of larger systems, it will
require too many tries for branched trajectories to reach $\Phi_{N_{c}}$. A
more efficient strategy is to anticipate the branching with a small shift,
this can for instance be done with a function that slowly reaches its
asymptote, the line of slope $1$,

$\Phi_{b,conv}=\frac{\Phi_{N_{c}}\left(1+\tanh\left(\frac{\Phi_{N_{c}}-\xi}{\xi}\right)\right)}{2}\,,$
(13)

This function is displayed in figure 2 (c), with the “converging” label.

On a final note, we state that the strong necessity for anticipated branching
may come from the imperfection of the reaction coordinate. Indeed, we have so
far used a reaction coordinate constructed with a formula under simple
physical arguments. For instance, in the case of the kinetic energy, we
considered that $E$ decreases as the flow goes toward the laminar state. In
the case of large systems, it may very well be that the committor, the
function which associates the probability of reaching the final state to (in
our case) an instantaneous velocity field, may not be so far from an affine
function of $E$. In the case of very small systems, which display a behaviour
close to temporal chaos, it is likely that the committor has very large
variation on isosurfaces of constant $E$, in particular when one moves
slightly away from turbulence. This is caused by the complex structure of the
laminar-turbulent boundary in small scale systems (Schmiegel and Eckhardt,,
1997). Both situations are relevant to collapse of turbulence in shear flows
in general. Indeed, pipe flows will most of the time display a localised
turbulent puff which keeps a small number of degrees of freedom as the length
of the pipe is increased. Meanwhile spatially extended channel flows such as
Couette or Poiseuille flow will display localised turbulence which will have
more and more degrees of freedom as the domain size is increased. The option
of anticipated branching should remain on the table as there is no certainty
that methods to estimate better reaction coordinate could work without prior
runs of AMS.

Figure 15: Time series of $\Phi$ during an AMS calculation with the $24\times
16$ system leading to extinction (the laminar state corresponds to $\phi\geq
1$). The black line indicates the end of the stage 0 of AMS and the red line
the maximal $\phi=0.248\pm 0.001$ reached during this stage (“contained in the
initial conditions”). The green line indicate the maximum $\phi=0.430\pm
0.001$ reached during the whole AMS calculation.

Figure 16: Collapse of turbulence in a very small system of size $L_{x}\times
L_{z}=12\times 8$. Left column: time series of kinetic energy with dot
indicating the instant in the simulation. middle column, colour levels of
streamwise velocity in the horizontal midplane. Right column: colour levels of
spanwise velocity in the horizontal midplane.

## References

  * Baars et al., (2021) Baars, S., Castellana, D., Wubs, F., and Dijkstra, H. (2021).  Application of adaptive multilevel splitting to high-dimensional dynamical systems.  Journal of Computational Physics, 424:109876.
  * Berhanu et al., (2007) Berhanu, M., Monchaux, R., Fauve, S., Mordant, N., Pétrélis, F., Chiffaudel, A., Daviaud, F., Dubrulle, B., Marié, L., Ravelet, F., Bourgoin, M., Odier, P., Pinton, J.-F., and Volk, R. (2007).  Magnetic field reversals in an experimental turbulent dynamo.  Eur. Phys. Lett., 77(5):59001.
  * Bottin and Chaté, (1998) Bottin, S. and Chaté, H. (1998).  Statistical analysis of the transition to turbulence in plane couette flow.  The European Physical Journal B-Condensed Matter and Complex Systems, 6(1):143–155.
  * Bouchet and Reygner, (2016) Bouchet, F. and Reygner, J. (2016).  Generalisation of the eyring–kramers transition rate formula to irreversible diffusion processes.  In Annales Henri Poincaré, volume 17, pages 3499–3532. Springer.
  * (5) Bouchet, F., Rolland, J., and Simonnet, E. (2019a).  Rare event algorithm links transitions in turbulent flows with activated nucleations.  Physical review letters, 122(7):074502.
  * (6) Bouchet, F., Rolland, J., and Wouters, J. (2019b).  Rare event sampling methods.  Chaos: An interdisciplinary journal of nonlinear science, 29(8):2.
  * Bréhier et al., (2016) Bréhier, C.-E., Gazeau, M., Goudenège, L., Lelièvre, T., Rousset, M., et al. (2016).  Unbiasedness of some generalized adaptive multilevel splitting algorithms.  The Annals of Applied Probability, 26(6):3559–3601.
  * (8) Cérou, F., Delyon, B., Guyader, A., and Rousset, M. (2019a).  On the asymptotic normality of adaptive multilevel splitting.  SIAM/ASA Journal on Uncertainty Quantification, 7(1):1–30.
  * Cérou and Guyader, (2007) Cérou, F. and Guyader, A. (2007).  Adaptive multilevel splitting for rare event analysis.  Stochastic Analysis and Applications, 25(2):417–443.
  * Cérou et al., (2013) Cérou, F., Guyader, A., Lelievre, T., and Malrieu, F. (2013).  On the length of one-dimensional reactive paths.  ALEA, 10(1):359–389.
  * Cérou et al., (2011) Cérou, F., Guyader, A., Lelievre, T., and Pommier, D. (2011).  A multiple replica approach to simulate reactive trajectories.  The Journal of chemical physics, 134(5):054108.
  * (12) Cérou, F., Guyader, A., and Rousset, M. (2019b).  Adaptive multilevel splitting: Historical perspective and recent results.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(4):043108.
  * De Lozar et al., (2012) De Lozar, A., Mellibovsky, F., Avila, M., and Hof, B. (2012).  Edge state in pipe flow experiments.  Physical review letters, 108(21):214502.
  * De Souza et al., (2020) De Souza, D., Bergier, T., and Monchaux, R. (2020).  Transient states in plane couette flow.  Journal of Fluid Mechanics, 903.
  * Devetsikiotis and Townsend, (1993) Devetsikiotis, M. and Townsend, J. K. (1993).  Statistical optimization of dynamic importance sampling parameters for efficient simulation of communication networks.  IEEE/ACM transactions on networking, 1(3):293–305.
  * Duguet et al., (2011) Duguet, Y., Le Maitre, O., and Schlatter, P. (2011).  Stochastic and deterministic motion of a laminar-turbulent front in a spanwisely extended couette flow.  Physical Review E, 84(6):066315.
  * Ebener et al., (2019) Ebener, L., Margazoglou, G., Friedrich, J., Biferale, L., and Grauer, R. (2019).  Instanton based importance sampling for rare events in stochastic pdes.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(6):063102.
  * Eckhardt et al., (2008) Eckhardt, B., Faisst, H., Schmiegel, A., and Schneider, T. M. (2008).  Dynamical systems and the transition to turbulence in linearly stable shear flows.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 366(1868):1297–1315.
  * Eckhardt et al., (2007) Eckhardt, B., Schneider, T. M., Hof, B., and Westerweel, J. (2007).  Turbulence transition in pipe flow.  Annu. Rev. Fluid Mech., 39:447–468.
  * Faisst and Eckhardt, (2003) Faisst, H. and Eckhardt, B. (2003).  Sensitive dependence on initial conditions in transition to turbulence in pipe flow.  arXiv preprint physics/0312078.
  * Gibson et al., (2008) Gibson, J., Halcrow, J., and Cvitanović, P. (2008).  Visualizing the geometry of state space in plane couette flow.  J. Fluid Mech., 611:107–130.
  * Glasserman et al., (1998) Glasserman, P., Heidelberger, P., Shahabuddin, P., and Zajic, T. (1998).  A large deviations perspective on the efficiency of multilevel splitting.  IEEE Transactions on Automatic Control, 43(12):1666–1679.
  * Gomé et al., (2020) Gomé, S., Tuckerman, L. S., and Barkley, D. (2020).  Statistical transition to turbulence in plane channel flow.  Physical Review Fluids.
  * Grafke and Vanden-Eijnden, (2019) Grafke, T. and Vanden-Eijnden, E. (2019).  Numerical computation of rare events via large deviation theory.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(6):063118.
  * Grandemange et al., (2013) Grandemange, M., Gohlke, M., and Cadot, O. (2013).  Turbulent wake past a three-dimensional blunt body. part 1. global modes and bi-stability.  Journal of Fluid Mechanics, 722:51–84.
  * Hamilton et al., (1995) Hamilton, J. M., Kim, J., and Waleffe, F. (1995).  Regeneration mechanisms of near-wall turbulence structures.  Journal of Fluid Mechanics, 287(1):317–348.
  * Hänggi et al., (1990) Hänggi, P., Talkner, P., and Borkovec, M. (1990).  Reaction-rate theory: fifty years after kramers.  Reviews of modern physics, 62(2):251.
  * Hartmann et al., (2019) Hartmann, C., Kebiri, O., Neureither, L., and Richter, L. (2019).  Variational approach to rare event simulation using least-squares regression.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(6):063107.
  * Herbert et al., (2020) Herbert, C., Caballero, R., and Bouchet, F. (2020).  Atmospheric bistability and abrupt transitions to superrotation: wave–jet resonance and hadley cell feedbacks.  Journal of the Atmospheric Sciences, 77(1):31–49.
  * Jiménez and Moin, (1991) Jiménez, J. and Moin, P. (1991).  The minimal flow unit in near-wall turbulence.  Journal of Fluid Mechanics, 225:213–240.
  * Kawahara et al., (2003) Kawahara, G., Jiménez, J., Uhlmann, M., and Pinelli, A. (2003).  Linear instability of a corrugated vortex sheet-a model for streak instability.  Journal of fluid mechanics, 483:315.
  * Kim and Durbin, (1988) Kim, H.-J. and Durbin, P. (1988).  Investigation of the flow between a pair of circular cylinders in the flopping regime.  Journal of Fluid Mechanics, 196:431–448.
  * l’Ecuyer et al., (2009) l’Ecuyer, P., Mandjes, M., Tuffin, B., et al. (2009).  Importance Sampling in Rare Event Simulation. Wiley Online Library.
  * Lestang et al., (2020) Lestang, T., Bouchet, F., and Lévêque, E. (2020).  Numerical study of extreme mechanical force exerted by a turbulent flow on a bluff body by direct and rare-event sampling techniques.  Journal of Fluid Mechanics, 895.
  * Lestang et al., (2018) Lestang, T., Ragone, F., Bréhier, C.-E., Herbert, C., and Bouchet, F. (2018).  Computing return times or return periods with rare event algorithms.  Journal of Statistical Mechanics: Theory and Experiment, 2018(4):043213.
  * Liu et al., (2020) Liu, T., Semin, B., Klotz, L., Godoy-Diana, R., Wesfreid, J. E., and Mullin, T. (2020).  Anisotropic decay of turbulence in plane couette-poiseuille flow.  arXiv preprint arXiv:2008.08851.
  * Lopes and Lelièvre, (2019) Lopes, L. J. and Lelièvre, T. (2019).  Analysis of the adaptive multilevel splitting method on the isomerization of alanine dipeptide.  Journal of computational chemistry, 40(11):1198–1208.
  * Manneville, (2011) Manneville, P. (2011).  On the decay of turbulence in plane couette flow.  Fluid Dynamics Research, 43(6):065501.
  * Marquillie et al., (2011) Marquillie, M., Ehrenstein, U., Laval, J.-P., et al. (2011).  Instability of streaks in wall turbulence with adverse pressure gradient.  Journal of Fluid Mechanics, 681(205-240):30.
  * Metzner et al., (2006) Metzner, P., Schütte, C., and Vanden-Eijnden, E. (2006).  Illustration of transition path theory on a collection of simple examples.  The Journal of chemical physics, 125(8):084110.
  * Moxey and Barkley, (2010) Moxey, D. and Barkley, D. (2010).  Distinct large-scale turbulent-laminar states in transitional pipe flow.  Proceedings of the National Academy of Sciences, 107(18):8091–8096.
  * Onsager, (1938) Onsager, L. (1938).  Initial recombination of ions.  Physical Review, 54(8):554.
  * Podvin and Sergent, (2017) Podvin, B. and Sergent, A. (2017).  Precursor for wind reversal in a square rayleigh-bénard cell.  Physical Review E, 95(1):013112.
  * Pope, (2001) Pope, S. B. (2001).  Turbulent flows.  IOP Publishing.
  * Ragone and Bouchet, (2019) Ragone, F. and Bouchet, F. (2019).  Computation of extreme values of time averaged observables in climate models with large deviation techniques.  Journal of Statistical Physics, pages 1–29.
  * Ragone and Bouchet, (2020) Ragone, F. and Bouchet, F. (2020).  Rare event algorithm study of extreme warm summers and heat waves over europe.  arXiv preprint arXiv:2009.02519.
  * Rolland, (2015) Rolland, J. (2015).  Mechanical and statistical study of the laminar hole formation in transitional plane couette flow.  The European Physical Journal B, 88(3):66.
  * Rolland, (2018) Rolland, J. (2018).  Extremely rare collapse and build-up of turbulence in stochastic models of transitional wall flows.  Physical Review E, 97(2):023109.
  * Rolland et al., (2016) Rolland, J., Bouchet, F., and Simonnet, E. (2016).  Computing transition rates for the 1-d stochastic ginzburg–landau–allen–cahn equation for finite-amplitude noise with a rare event algorithm.  Journal of Statistical Physics, 162(2):277–311.
  * Rolland and Simonnet, (2015) Rolland, J. and Simonnet, E. (2015).  Statistical behaviour of adaptive multilevel splitting algorithms in simple models.  Journal of Computational Physics, 283:541–558.
  * Romanov, (1973) Romanov, V. A. (1973).  Stability of plane-parallel couette flow.  Functional analysis and its applications, 7(2):137–146.
  * Schmiegel and Eckhardt, (1997) Schmiegel, A. and Eckhardt, B. (1997).  Fractal stability border in plane couette flow.  Physical review letters, 79(26):5250.
  * Schneider et al., (2007) Schneider, T. M., Eckhardt, B., and Yorke, J. A. (2007).  Turbulence transition and the edge of chaos in pipe flow.  Physical review letters, 99(3):034502.
  * Shi et al., (2013) Shi, L., Avila, M., and Hof, B. (2013).  Scale invariance at the onset of turbulence in couette flow.  Physical review letters, 110(20):204502.
  * Simonnet, (2016) Simonnet, E. (2016).  Combinatorial analysis of the adaptive last particle method.  Statistics and Computing, 26(1-2):211–230.
  * Simonnet et al., (2021) Simonnet, E., Rolland, J., and Bouchet, F. (2021).  Multistability and rare spontaneous transitions between climate and jet configurations in a barotropic model of the jovian mid-latitude troposphere.  Journal of the Atmospheric Sciences, 78:1889–1911.
  * Touchette, (2009) Touchette, H. (2009).  The large deviation approach to statistical mechanics.  Physics Reports, 478(1-3):1–69.
  * Waleffe, (1997) Waleffe, F. (1997).  On a self-sustaining process in shear flows.  Physics of Fluids, 9(4):883–900.
  * Wan and Yu, (2017) Wan, X. and Yu, H. (2017).  A dynamic-solver–consistent minimum action method: With an application to 2d navier–stokes equations.  Journal of Computational Physics, 331:209–226.