<figure><img src="image_1.png" /></figure>

HAL
open science

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

Joran Rolland

<figure><img src="image_2.png" /></figure>

## To cite this version:

Joran Rolland. Does rare, noise-induced, bypass transition in plane Couette flow bypass instantons?.
2024. hal-04386611

HAL Id: hal-04386611

https://hal.science/hal-04386611v1

Preprint submitted on 10 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Does rare, noise-induced, bypass transition in plane Couette flow bypass instantons ?

Joran Rolland *†

January 10, 2024

## Abstract

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

## 1 Introduction

Many physical systems, after they underwent a subcritical bifurcation or a first order transition (or a similar transition in more complex systems), display two or more metastable states. This situation means that the physical system can switch from one metastable state to another if it undergoes some statistically steady forcing, usually by crossing a boundary separating the basin of attraction of one state from another. Conversely, this physical system will relax toward either states if it is prepared with an initial condition

<footnote>*Univ. Lille, CNRS, ONERA, Arts et Métiers Institute of Technology, Centrale Lille, UMR 9014 - LMFL - Laboratoire de Mécanique des Fluides de Lille - Kampé de Fériet, F-59000 Lille, France.</footnote>

<footnote>†joran.rolland@centralelille.fr</footnote>
on the corresponding side of the separatrix. Traditional bistable systems are found in kinetic chemistry [Onsager, 1938, Hänggi et al., 1990, Van Erp et al., 2003, Lopes and Lelièvre, 2019], where the two states are the reactants and the products or different configurations or conformations [Lopes and Lelièvre, 2019] of a complex molecule, while the force driving transitions is the thermal noise. Much insight on bistability has been gained from formal noise driven systems [Metzner et al., 2006, Cérou et al., 2011, Bouchet and Touchette, 2012, Rolland and Simonnet, 2015, Rolland et al., 2016, Lucente et al., 2022]. More recently, a lot of effort has also been directed toward the study of bistable climate systems [Lucarini and Bódai, 2017, Lucarini and Bódai, 2019, Baars et al., 2021, Herbert et al., 2020] and bistable turbulent flows [Podvin and Sergent, 2017] of aerodynamical [Grandemange et al., 2013, Kim and Durbin, 1988] or geophysical interest [Berhanu et al., 2007, Bouchet et al., 2019]. In these systems, the metastable states are usually large scale coherent circulations, like jets, or vortices, while the force driving transition is intrinsic to the system and related to the fluctuating small scales of turbulence.

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

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

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

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

This was the case of the most recent works, concerning the collapse of turbulence in plane Couette flow [Rolland, 2022, Gomé et al., 2022]. Plane Couette flow (PCF) is the flow between two parallel moving walls (Fig. 1, left). In PCF, like in other wall flows such as plane Poiseuille flow [Wan, 2013, Wan et al., 2015, Wan and Yu, 2017], Couette-Poiseuille flow [Liu et al., 2021], Hagen-Poiseuille pipe flow [Willis and Kerswell, 2009] or boundary layers [Spangler and Wells Jr, 1968, Rigas et al., 2021], wall turbulence can coexist at Reynolds number (or positions in the case of boundary layers) for which the laminar baseflow is linearly stable (see [Romanov, 1973] for the case of PCF). As a consequence, transitional turbulence can collapse under its own fluctuations, while the flow can go from laminar to turbulent under a forcing [Spangler and Wells Jr, 1968, Rigas et al., 2021] or if it is given an initial condition “close enough” to turbulence [Schmiegel and Eckhardt, 1997]. The study of the collapse of turbulence using rare event methods gave the occasion to investigate the physical mechanism of laminar hole opening as well as the technical solutions to study rare multistability events in a three dimensional flow not stochastically forced [Rolland, 2022, Gomé et al., 2022]. Conversely, the study of the forcing of laminar flow toward turbulence, under a noise or another, would give the occasion to investigate the crossing of a very complex separatrix<fnref n="1" />.

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

In wall flows, the complexity of this boundary is such that the saddle on it is not the point of the separatrix which has the smallest amplitude of departure to the laminar baseflow. Experimental and numerical studies of initial conditions, generated using various specific rules, of minimal energy that evolve into turbulence have highlighted that. This leaves the question of whether the most probable forced path would visit this saddle at very small but finite forcing noise variance. Alternate investigations have thus been aimed at identifying a path forced and/or triggered by a finite amplitude initial condition that would drive the flow from laminar to turbulent. Linear studies of optimal initial conditions of infinitesimal amplitude leading to maximal amplification of energy identified lift-up: the extraction of energy from

<footnote n="1">¹We use the term separatrix for lack of a better word. strictly speaking, an unforced wall flow starting from any initial condition eventually reaches the laminar state: what changes is the relaxation time before it does so. This duration dramatically increases as an hypersurface is crossed [Schmiegel and Eckhardt, 1997].</footnote>
the laminar base flow from streamwise vortices. The lift-up mechanism was thus identified in transient
growth study that highlight the fact that energy growth can happen while all linear modes are stable
is the linearised Navier–Stokes equations lead to a non normal operator [Farrell and Ioannou, 1993a,
Schmid and Henningson, 1994, Trefethen et al., 1993]. Such an optimal (in shape) initial condition whose
amplitude places it on the laminar side of the boundary would remain on the laminar side of said boundary.
However, for initial conditions starting on the turbulent side of the boundary, lift-up is an important
part of the self-sustaining process [Waleffe, 1997]. Refinement of optimisation of initial condition energy
amplification technics can be used to compute the optimal finite amplitude initial condition leading to
maximal non-linear energy growth [Cherubini et al., 2011, Rabin et al., 2012, Monokrousos et al., 2011].
This initial condition is termed a minimal seed. Such computations are still performed on the laminar side
of the separatrix, however a slight modification of the seed, in the right direction, can generate an initial
condition which has very little energy but can nevertheless travel to turbulence extremely efficiently.

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

This motivates our present study. As to the rare event simulation method, we use a method al-
ternate to action minimisation termed Adaptive Multilevel Splitting (AMS) [Cérou and Guyader, 2007,
Cérou et al., 2011], to compute the rare transitions from laminar baseflow to turbulence in plane Cou-
ette flow, in domains of increasing sizes at transitional Reynolds numbers. By decreasing the variance
of the forcing, proportional to the energy injection rate, we will investigate the possible concentration
of reactive trajectories and whether this trajectories are structured by an instanton. Using several noise
spectra, we will test whether forcing directly the most responsive structure is a sine qua non condition
to travel to turbulence, while increasing the system size will highlight possible spatial localisation of the
growing turbulence. We thus present this study in the following manner. We first present the forced
Navier–Stokes equations (§ 2.1), then the version of Adaptive Multilevel Splitting that is used to compute
the trajectories (§ 2.2), including the details of the set up of AMS computations and additional diag-
nostics used to study the trajectories such as bisections. We describe the computed paths in section 3,
first visually in relevant examples (§ 3.1), then statistically, on the whole collection of computed path
(§ 3.2). We finally study the effect of the energy injection rate, on the transition rate (§ 4.1) and on the
trajectories duration (§ 4.2). We will discuss these results in view of the literature in the conclusion (§ 5).
<figure><img src="image_4.png" /><figcaption>Figure 1: A picture of transition to turbulence in plane Couette flow indicating the type of flow obtained as the Reynolds number is increased. Plane Couette flow is sketched on the left: we indicate the two walls moving at velocities ±U, the domain size $L_x, L_z$ and h as well as the coordinate system. For Reynolds numbers larger than a $R_g \simeq 325$, laminar-turbulent coexistence can be sustained, as indicated by the colour levels of kinetic energy $\frac{1}{2}u^2$ in the horizontal midplane $y = 0$ obtained by means of direct numerical simulations in a domain of size $L_x \times L_z = 110 \times 32$ at $R = 370$ (centre). The coloured frames indicate the domain sizes considered here (cyan: $L_x \times L_z = 6 \times 4$, blue $L_x \times L_z = 18 \times 12$, green $L_x \times L_z = 36 \times 24$). As the Reynolds number is further increased, wall turbulence can invade the whole domain, as illustrated by the colour levels of kinetic energy $\frac{1}{2}u^2$ in the horizontal midplane $y = 0$ obtained by means of direct numerical simulations in a domain of size $L_x \times L_z = 110 \times 32$ at $R = 370$ (right). The visualisations are performed using numerical results originally performed for [Rolland, 2015], and the whole figure was originally created for [Rolland, 2022].</figcaption></figure>

# 2 Method

## 2.1 Forced Navier-Stokes equations

### 2.1.1 Governing equations

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

The forced incompressible Navier-Stokes equations for the field $u$, the departure to the laminar base flow $y e_x$, and $p$ the dimensionless pressure, read

$$ \frac{\partial u_l}{\partial t} + u_j \frac{\partial u_l}{\partial x_j} + y \frac{\partial u_l}{\partial x} + \delta_{l,x} u_y = - \frac{\partial p}{\partial x_l} + \frac{1}{R} \left( \frac{\partial^2 u_l}{\partial x^2} + \frac{\partial^2 u_l}{\partial y^2} + \frac{\partial^2 u_l}{\partial z^2} \right) + \sqrt{\frac{2}{\beta}} f_l(x, t), \quad \frac{\partial u_l}{\partial x_l} = 0. \quad (2) $$

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

* $\frac{3}{2}N_x = 32$, $\frac{3}{2}N_z = 24$, $N_y = 27$ and $\Delta t = 0.05$, that we termed standard resolution,

* $\frac{3}{2}N_x = 64$, $\frac{3}{2}N_z = 48$, $N_y = 47$, and $\Delta t = 0.0125$, that we termed high resolution,

* $\frac{3}{2}N_x = 64$, $\frac{3}{2}N_z = 48$, $N_y = 27$ and $\Delta t = 0.05$,

* $\frac{3}{2}N_x = 128$, $\frac{3}{2}N_z = 96$, $N_y = 35$ and $\Delta t = 0.025$,

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

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

$$ \langle f_l(\mathbf{x}, t)f_l(\mathbf{x}', t') \rangle = \delta(t-t')C_l(x-x', z-z')\delta(y-y'), \quad \hat{C}_l(n_x, n_z) = \Gamma_{l,n_x,n_z}. \quad (3) $$

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

### 2.1.2 Diagnostics of the flow

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

$$ E_k = \frac{1}{2L_x L_z} \int_{x=0}^{L_x} \int_{y=-1}^{1} \int_{z=0}^{L_z} \frac{u_x^2 + u_y^2 + u_z^2}{2} dxdydz. \quad (4) $$
<figure><img src="image_5.png" /><figcaption>Figure 2: Colour levels in the wavenumber $2\pi n_x/L_x$, $2\pi n_z/L_z$ plane of shape factors of divergence free forcing used for AMS computations in domains of size $L_x \times L_z = 6 \times 4$ at resolution $N_x = 64, N_y = 48, N_z = 48$ at Reynolds number $R = 500$: (a) logarithm of the product (($2\pi n_z/L_z)\gamma_{n_x}\gamma_{n_z}$)$^2$ for forcing applied on $u_x$, (b) logarithm of the product (($2\pi n_x/L_x)\gamma_{n_x}\gamma_{n_z}$)$^2$ for forcing applied on $u_z$ (the details of $\gamma_x$ and $\gamma_z$ are given in appendix A). See appendix A for full detail on how this generates the vector potential used to exert a divergence free forcing on the flow. From the velocity field obtained in a DNS of a turbulent flow in a domain of $L_x \times L_z = 6 \times 4$ at resolution $N_x = 64, N_y = 48, N_z = 48$ at Reynolds number $R = 500$, (c) logarithm of the modulus of $x-z$ spectrum of $u_x$ in the midplane $y=0$, (d) logarithm of the modulus of $x-z$ spectrum of $u_z$ in the midplane $y=0$.</figcaption></figure>
The kinetic energy is used to monitor the global state of the flow. When the flow is turbulent in the whole domain, we find $E_k = O(0.1)$. Meanwhile, when the flow is globally laminar, $E_k$ is much smaller and close to zero. We will also distinguish $E_{k,l}$ the kinetic energy contained in a given component $l = x, y, z$, on the one hand, and $E_{k,y-z}$ the kinetic energy contained in the spanwise and wall normal components, on the other hand

$$E_{k,l} = \frac{1}{2L_x L_z} \int_{x=0}^{L_x} \int_{y=-1}^{1} \int_{z=0}^{L_z} \frac{u_l^2}{2} dxdydz, \quad 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} dxdydz. \quad (5)$$

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

$$e(\tilde{n}_z) = \frac{1}{2L_x \delta z} \int_{x=0}^{L_x} \int_{y=-1}^{1} \int_{z=\tilde{n}_z \delta z}^{(\tilde{n}_z+1)\delta z} (u_y^2 + u_z^2) dxdydz, \quad (6)$$

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

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

$$\frac{\partial \langle E_k \rangle}{\partial t} = - \underbrace{\left\langle \frac{1}{2L_x L_z} \int u_x u_y \, dxdydz \right\rangle}_{\text{extraction}} - \underbrace{\left\langle \frac{1}{R2L_x L_z} \int \frac{\partial u_l}{\partial x_j} \frac{\partial u_l}{\partial x_j} \, dxdydz \right\rangle}_{\text{dissipation}} + \epsilon, \quad (7)$$

where we have the energy injection rate $\epsilon$

$$\epsilon = \frac{1}{\beta} \left( N_y \sum_{n_x = -\frac{N_x}{2} + 1}^{\frac{N_x}{2}} \sum_{n_z = -\frac{N_z}{2} + 1}^{\frac{N_z}{2}} \sum_{l=1}^{3} \Gamma_{l,n_x,n_z} \right), \quad (8)$$

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

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

$$ 
\begin{aligned}
\frac{dE_{tot}}{dt} &= \frac{1}{R}(I - D) + \epsilon, \quad E_{tot} = \int_{x=0}^{L_x} \int_{y=-1}^{1} \int_{z=0}^{L_z} \frac{(u_x + y)^2 + u_y^2 + u_z^2}{2} dxdydz \\
I &= 1 + \frac{1}{2L_xL_z} \int_{x=0}^{L_x} \int_{z=0}^{L_z} \left( \frac{\partial u_x}{\partial y}(y = +1) + \frac{\partial u_x}{\partial y}(y = -1) \right) dxdz, \\
D &= \frac{1}{2L_xL_z} \int_{x=0}^{L_x} \int_{y=-1}^{1} \int_{z=0}^{L_z} \|\nabla \times (\mathbf{u} + y\mathbf{e}_x)\|_2^2 dxdydz,
\end{aligned}
\quad (9) $$

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

## 2.2 Adaptive Multilevel Splitting

### 2.2.1 The algorithm

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

We then present the generic formulation for two of the variants of Adaptive Multilevel Splitting (AMS) which were used for the computation presented in this text. They were initially proposed by Cérou &
<figure><img src="image_6.png" /><figcaption>Figure 3: (a) Sketch of two bistable states **A** and **B** and the hypersurface **C** closely surrounding **A**. two realisations of the dynamics are sketched: a single excursion in blue and a first passage trajectory in black and red. The red part of the first passage trajectory is the reactive trajectory (figure originally made for [Rolland and Simonnet, 2015]). (b) Sketch of the principle of AMS, showing to iterations of the algorithm, with *N* = 3 clones, indicating the starting state **A** and its neighbourhood, the arrival state and its neighbourhood **B**, three trajectories are ordered by their max<sub>*t*</sub> Φ<sub>*</sub>. Trajectory one (dashed blue line) is suppressed and branched on another trajectory at level max<sub>*t*</sub> Φ<sub>1</sub> and then ran according to its natural dynamics. Trajectory 2 is then suppressed and branched on 3 at level max<sub>*t*</sub> Φ<sub>2</sub> (figure originally made for [Simonnet, 2016]).</figcaption></figure>

Guyader [Cérou and Guyader, 2007]. Both variants of AMS use a reaction coordinate (or observable) $\phi(\mathbf{u}(t)) = \Phi(t)$. The reaction coordinate $\phi$ is a mapping from $\mathbb{R}^D$, where the discretised velocity field $\mathbf{u}$ lives, to $[0, 1]$. The reaction coordinate can be viewed as a function of time $\Phi(t)$ on trajectories. The reaction coordinate measures the position of the flow relatively to starting set $\mathcal{A}$, where is it 0 and arrival set $\mathcal{B}$ where it is 1. We will define our reaction coordinate quantitatively in section 2.2.2. All variants of AMS run $N$ clone dynamics of the system to compute iteratively from $N - N_c + 1$ up to $N$ reactive trajectories going from set $\mathcal{A}$ to the hypersurface $\mathcal{C}$, closely surrounding set $\mathcal{A}$, then from $\mathcal{C}$ to the set $\mathcal{B}$. Here $1 \le N_c < N$ is the number of trajectories that will be removed at each iteration. We will label the velocity field and pressure field of each trajectory by $1 \le i \le N$ as in $(\mathbf{u}_i \ p_i)(t)$. The reaction coordinate as a function of time for the corresponding trajectory is then denoted by $\Phi_i(t) = \phi(\mathbf{u}_i(t))$. The algorithm is sketched in figure 3 (b) and proceeds in the following manner

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

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

The algorithm stops after iteration number $\kappa \in \mathbb{N}$, the first iteration where at least $N - N_c + 1$ clones trajectories have reached $\mathcal{B}$. We denote by $r \in [1 - \frac{N_c-1}{N}, 1]$ the proportion of clones that have reached $\mathcal{B}$ at the last stage. Note that $\kappa$ and $r$ are a priori random numbers, different from one AMS run to the other, whose value depend on the independent noise realisations drawn during simulation of trajectories.
<figure><img src="image_7.png" /><figcaption>Figure 4: Average kinetic energy conditioned on no build up having occurred and corresponding variance as a function of the energy injection rate for several shape of the forcing spectrum and/or several resolutions when (a) noise on components is applied, (b) when divergence free noise is applied.</figcaption></figure>

This yields an estimator of the probability $\alpha$ of reaching $\mathcal{B}$ before $\mathcal{A}$ starting from the distribution of first exit points on $\mathcal{C}$, and the corresponding mean first passage time $T$

$$ \alpha = r \left(1 - \frac{N_c}{N}\right)^{\kappa}, \quad \hat{\alpha} = \langle \alpha \rangle_o, \quad \hat{T} = \left\langle \left(\frac{1}{\alpha} - 1\right) \tilde{\tau} + \tau \right\rangle_o, \qquad (10) $$

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

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

### 2.2.2 Reaction coordinate

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

$$ \Phi(t) = \frac{E_k(t) - E_A}{E_B}. \qquad (11) $$

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

$$C = \{ \mathbf{u} \setminus \phi(\mathbf{u}) = \frac{\sigma}{E_B} \}, \qquad (12)$$

where $\sigma$ is the standard deviation of $E_k(t)$ conditioned on no transition occurring, sampled from the same dataset as $E_A$. The rationale is that the flow leaves the neighbourhood of the laminar state when $E_k \ge E_A + \sigma$.

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

In order to propose the interpolated values, we adjust the linear dependence of $E_A$ and $\sigma$ on $\beta$ in logarithmic scale for each set of parameters setting the forcing spectrum. We fit the constants $d_E$, $o_E$, $d_S$ and $o_S$ from data such that

$$\log(E_A) = \log(\beta)d_E + o_E, \quad \log(\sigma) = \log(\beta)d_S + o_S. \qquad (13)$$

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

$$d_E = \frac{M \sum_{a=1}^{M} (\log(E_a) \log(\beta_a)) - (\sum_{a=1}^{M} \log(E_a)) (\sum_{a=1}^{M} \log(\beta_a))}{M \sum_{a=1}^{M} (\log(\beta_a))^2 - (\sum_{a=1}^{M} \log(\beta_a))^2}, \quad (14)$$

$$o_E = \frac{\left(\sum_{a=1}^{M} \log(\beta_a)\right) \left(\sum_{a=1}^{M} \log(\beta_a) \log(E_a)\right) - \left(\sum_{a=1}^{M} \log(E_a)\right) \left(\sum_{a=1}^{M} \log(\beta_a)\right)^2}{\left(\sum_{a=1}^{M} \log(\beta_a)\right)^2 - M \sum_{a=1}^{M} \log(\beta_a)^2}, \quad (15)$$

$$d_S = \frac{M \sum_{a=1}^{M} (\log(\sigma_a) \log(\beta_a)) - \left(\sum_{a=1}^{M} \log(\sigma_a)\right) \left(\sum_{a=1}^{M} \log(\beta_a)\right)}{M \sum_{a=1}^{M} (\log(\beta_a))^2 - \left(\sum_{a=1}^{M} \log(\beta_a)\right)^2}, \quad (16)$$

$$o_S = \frac{\left(\sum_{a=1}^{M} \log(\beta_a)\right) \left(\sum_{a=1}^{M} \log(\beta_a) \log(\sigma_a)\right) - \left(\sum_{a=1}^{M} \log(\sigma_a)\right) \left(\sum_{a=1}^{M} \log(\beta_a)\right)^2}{\left(\sum_{a=1}^{M} \log(\beta_a)\right)^2 - M \sum_{a=1}^{M} \log(\beta_a)^2}. \quad (17)$$

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

We give the value used to define the sets $A$, $C$ and $B$ for each set of simulations in appendix B. This comprises of specific cases where we will give $E_A$, $E_B$ and $\sigma$ and generic cases where we give $d_E$, $o_E$, $d_S$ and $o_S$, which mostly concern simulations performed in the largest systems.

### 2.2.3 Direct Numerical simulations of transitions

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

### 2.2.4 Methodology for bisections

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

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

$$ \begin{pmatrix} \mathbf{u}_{m,n} \\ p_{m,n} \end{pmatrix} (t_{m,0}) = \eta_n \begin{pmatrix} \mathbf{u}_{\text{lamin},m} \\ p_{\text{lamin},m} \end{pmatrix} + (1 - \eta_n) \begin{pmatrix} \mathbf{u}_{\text{turb},m} \\ p_{\text{turb},m} \end{pmatrix}. \quad (18) $$

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

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

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

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

$$ \begin{align} \begin{pmatrix} \mathbf{u}_{m,-} \\ p_{m,-} \end{pmatrix} (t = t_{m,0}) &= \eta_{\text{lamin}} \begin{pmatrix} \mathbf{u}_{\text{lamin},m} \\ p_{\text{lamin},m} \end{pmatrix} + (1 - \eta_{\text{lamin}}) \begin{pmatrix} \mathbf{u}_{\text{turb},m} \\ p_{\text{turb},m} \end{pmatrix}, \notag \\ \begin{pmatrix} \mathbf{u}_{m,+} \\ p_{m,+} \end{pmatrix} (t = t_{m,0}) &= \eta_{\text{turb}} \begin{pmatrix} \mathbf{u}_{\text{lamin},m} \\ p_{\text{lamin},m} \end{pmatrix} + (1 - \eta_{\text{turb}}) \begin{pmatrix} \mathbf{u}_{\text{turb},m} \\ p_{\text{turb},m} \end{pmatrix}. \tag{19} \end{align} $$

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

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

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

# 3 Paths toward turbulence

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

## 3.1 Velocity fields in selected examples of transition paths

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

We will observe the transition paths on four examples:

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

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

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

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

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

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

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

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

3.1.2 Last state at the last stage

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

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

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

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

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

In the aforementioned figure, we include the streamwise velocity in the $y = 0$ plane of this state (Fig. 9 (a,c,e,g)) and the streamwise vorticity in the $y = 0$ plane of this state (Fig. 9 (b,d,f,h)), obtained from the AMS computations that yielded the reactive trajectories presented in section 3.1.1. This diagnostic indicates that at this turning point of AMS computations, streamwise elongated streamwise velocity tubes are formed and display turbulent like intensity with $\max_{\text{space}} |u_x| \ge 0.5$, although they present little trace of streamwise modulation in the smaller, MFU type, system (Fig. 9 (a,c)). The corresponding streamwise vorticity (Fig. 9 (b,d)) thus has a maximum intensity of $|\omega_x| \lesssim 0.5$, which is a fraction of what can be found when the flow is actually turbulent. Again, we note that there is no striking qualitative difference in those states obtained two distinct types of forcing. Streamwise velocity tubes are also formed in larger systems (Fig. 9 (e,g)). However, we note that in a restricted range of $z$ (taking into account periodic boundary conditions, $0 \lesssim z \lesssim 3$ and $10 \lesssim z \lesssim 12$ for $L_x \times L_y = 18 \times 12$ Fig. 9 (e), $15 \lesssim z \lesssim 24$ for $L_x \times L_y = 36 \times 24$ Fig. 9 (g)), wall turbulence like streamwise modulation of $u_x$ is present. We note that
<figure><img src="image_8.png" /><figcaption>Figure 5: View of a reactive trajectory for a system of size $L_x \times L_z = 6 \times 4$ at Reynolds number $R = 500$ computed at high resolution under a divergence free forcing at $\beta = 5.5 \cdot 10^4$, along with the result of relaxations and bisections started from this trajectory. (a) Colour levels of the streamwise velocity in the plane $y=0$ when the flow has reached the turbulent state. (b) Colour levels of the streamwise vorticity in the plane $y=0$ at $t = 440$ (first time step relaxing to turbulence). (c) Colour levels of the streamwise velocity field in the plane $y=0$ at $t = 440$. (d) Streamwise average of the streamwise velocity field in the $y-z$ plane at a point on the unstable limit cycle reached after from the bisections. (e) Colour levels of the streamwise vorticity field in the plane $y=0$ at a point on the unstable limit cycle reached after from the bisections. (f) Time series of the kinetic energy during a reactive trajectory (blue), during two relaxation from states taken on the reactive trajectory, at $t = 439$ (last relaxation to laminar state, red), $t = 440$ (first relaxation to turbulent state, orange) and on the bisections started from these two successive velocity fields (green).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: View of a reactive trajectory, relaxations and bisections started from said trajectory, for a system of size $L_x \times L_z = 6 \times 4$ at Reynolds number $R = 500$ under a forcing on component $x$ at $\beta = 1.1 \cdot 10^5$. (a) Isovalues $u_x = -0.3$ (blue) and $u_x = 0.3$ (red) of streamwise velocity in the whole domain when the flow has reached the turbulent state. (b) Colour levels of the streamwise average of the streamwise velocity field in the vertical plane $z-y$ for the state relaxing to laminar flow used in the bisection. (c) Colour levels of the streamwise velocity field in the $y=0$ plane when bisections have reached the neighbourhood of the unstable limit cycle. (d) Colour levels of the streamwise vorticity field in the $y=0$ plane when bisections have reached the neighbourhood of the unstable limit cycle. (e) Time series of the kinetic energy during a reactive trajectory (blue), during two relaxation from states taken on the reactive trajectory, at $t = 1048$ (last relaxation to laminar state, red), $t = 1050$ (first relaxation to turbulent state, orange) and on the dichotomies started from these two successive velocity fields (green).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7: View of a reactive trajectory for a system of size $L_x \times L_z = 18 \times 12$ at Reynolds number $R = 500$ under a divergence free forcing at $\beta = 6.25 \cdot 10^4$, along with the relaxations and bisections started from this trajectory. (a) Colour levels of the streamwise velocity in the plane $y = 0$ when the flow has reached the turbulent state. (b) Right: colour levels of streamwise velocity in the plane $y = 0$ at $t = 50$ in the early stage of the reactive trajectory, left: colour levels of the streamwise average of the streamwise velocity in a $y-z$ plane at the same time. (c) Colour levels of the streamwise vorticity in the plane $y = 0$ at time $t = 182$ (last time step relaxing to laminar flow) from which bisections are started. (d) Colour levels of the streamwise velocity in the plane $y = 0$ at time $t = 182$. (e) Colour levels of the streamwise velocity in the plane $y = 0$ when the bisections have reached a saddle. (f) Colour levels of the streamwise vorticity in the plane $y = 0$ when the bisections have reached a saddle. (g) Time series of the kinetic energy during the reactive trajectory (blue), during two relaxations from states taken on the reactive trajectory, at $t = 182$ (last relaxation to laminar state, red), $t = 183$ (first relaxation to turbulent state, orange) and on the bisections started from these two successive velocity fields (green).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: View of a reactive trajectory for a system of size $L_x \times L_z = 36 \times 24$ at Reynolds number $R = 500$ under a divergence free forcing at $\beta = 3 \cdot 10^4$, and the bisections started from the trajectory. (a) Colour levels of streamwise velocity in the plane $y = 0$ when the flow has reached the turbulent state. (b) Right: colour levels of streamwise velocity in the plane $y = 0$ at $t = 50$ in the early stage of the reactive trajectory, left: colour levels of the streamwise average of the streamwise velocity in a $y-z$ plane at the same time. (c) Colour levels of the streamwise vorticity in the plane $y = 0$ at time $t = 217$ at which bisections are started. (d) Colour levels of the streamwise velocity in the plane $y = 0$ at the time $t = 217$ at which bisections are started. (e) Colour levels of the streamwise velocity in the plane $y = 0$ when the bisections have reached a saddle. (f) Colour levels of the streamwise vorticity in the plane $y = 0$ when the bisections have reached a saddle. (g) Time series of the kinetic energy during a reactive trajectory (blue), during two relaxation from states taken on the reactive trajectory, at $t = 217$ (last relaxation to laminar state, red), $t = 218$ (first relaxation to turbulent state, orange) and on the dichotomies started from these two successive velocity fields (green).</figcaption></figure>
streamwise vorticity of turbulent like intensity $\omega_x = O(1)$ and with turbulent like scale is now only found
in a similar range of $z$ (taking into account periodic boundary conditions, $10 \lesssim z \lesssim 12$ and to a lesser
extent $0 \lesssim z \lesssim 3$ for $L_x \times L_z = 18 \times 12$ Fig. 9 (f), $18 \lesssim z \lesssim 24$ for $L_x \times L_z = 36 \times 24$ Fig. 9 (h)). Vorticity
is more intense in the larger system. This further indicates that while the reactive trajectory corresponds
to a response to the forcing through the development of streamwise velocity tube, the flow only reach
a turning point when locally, the tubes become sufficiently intense and display intense enough shear
between one another to trigger the formation of streamwise modulation and corresponding streamwise
vortices. Once this has happened, this streamwise modulation and said vortices will contaminate nearby
regions in $z$ until they extend to the whole flow.

3.1.3 Bisections

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

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

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

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

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

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

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

$$ \psi(t) = L_x \left( \frac{\arg(\hat{u}_z(n_x = 1, y = 0, z = 2, t))}{2\pi} + \frac{1}{2} \right). \qquad (20) $$

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

As the system size is increased to larger values $12 \times 8$ (not shown) and $L_x \times L_z = 18 \times 12$, several changes occur in this saddle. The saddle becomes time invariant, both for kinetic energy (Fig. 7 (g)) and velocity and vorticity fields (Fig. 7 (e,f)): it is thus an actual saddle point. It still presents streamwise modulated streamwise velocity tubes that occupy the whole domain, however, these velocity tubes are less intense with $\max_{\text{space}} |u_x| = 0.44 \pm 0.01$. This leads to a steady value of kinetic energy of $E_k = 0.2$, a third of the time average found in the MFU type system. Additionally, the streamwise modulation in the adjacent streamwise velocity tubes is of varicose type. The spatial organisation and average kinetic energy of this saddle is consistent with saddles computed at $R = 400$ in domains of comparable size [Schneider et al., 2008]. The spatial variation of the streamwise vorticity is less sharp than in the MFU type systems and its amplitude is still moderate $|\omega_x| = O(0.1)$ (Fig. 7 (f)). As the system size is further increased to $L_x \times L_z = 36 \times 24$, yet another change appears in the saddle which is reached after bisections: the state becomes localised in $z$ (Fig. 8 (e)) and modulated in $x$ to a smaller extent. the most intense
<figure><img src="image_13.png" /><figcaption>Figure 10: (a) Rescaled phase of the first *x* Fourier mode of $u_z$ computed at $z = 2$ and $y = 0$ as a function of time computed for the limit cycle found in MFU type domains (Fig. 5 (d,e), Fig. 6 (c,d)). (b) Rescaled phase of the third and fourth *x* Fourier mode of $u_z$ computed at $z = 12$ and $y = -0.75$ as a function of time computed for the travelling wave found in the domain of size $L_x \times L_z = 36 \times 24$ (Fig. 8 (e,f)).</figcaption></figure>

streamwise velocity tube has $u_x \ge 0$ in our computation, with $\max_{\text{space}} |u_x| = 0.78 \pm 0.01$, it is flanked by pairs of negative then positive streamwise velocity tube of comparable, though decreasing intensity. The streamwise velocity then decays rapidly with $z$. Moreover the modulation of streamwise velocity is localised in $x$, leading to a streamwise localisation of streamwise vorticity (Fig. 8 (f)). The streamwise propagation of this modulation makes this saddle a travelling wave. This travelling wave is consistent with those observed at smaller Reynolds numbers $R = 400$, in domains slightly larger in the spanwise direction and larger in the streamwise direction [Schneider et al., 2010].

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

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

## 3.2 Statistical description of trajectories

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

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

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

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

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

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

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

For one degree of freedom systems, in the limit of vanishing noise variance, it has been demonstrated that this is the distribution of duration [Cérou et al., 2013]. It has been observed that the sampled distribution is close to a Gumbel in more complex stochastic systems [Rolland et al., 2016, Rolland, 2018] when we have instanton phenomenology (with clear deviation when we do not, for instance when the
<figure><img src="image_14.png" /><figcaption>Figure 11: (a) Plot of a reactive trajectory obtained by means of AMS in the MFU type domain at $\beta = 600000$ (blue line), relaxations toward laminar (red line) and turbulent flow (orange line) started from the reactive trajectory as well as the result of the bisections (green line) in the plane formed by the logarithm of kinetic energy on components $\log(E_{k,x})$ and $\log(E_{k,y-z})$. An inset zooms on the area where the dynamics winds around the unstable limit cycle. (b) Plot of the same reactive trajectory, relaxation and result of bisections, in the plane $I, D$ where $I$ is the total energy injection and $D$ is the total dissipation (Eq. (9)), using the same colour code. The dashed line $D=I$ is added. An inset zooms on the area where the dynamics winds around the unstable limit cycle. (c) Average plus and minus standard deviation of the logarithm of the kinetic energy in $y-z$ components $\log(E_{k,y-z})$ conditioned on logarithm of the kinetic energy in $x$ component $\log(E_{k,x})$ obtained from reactive trajectories computed by means of AMS with a divergence free forcing and increasing $\beta$ as well as with forcing on the $x$ component. The path in the $(\log(E_{k,x}), \log(E_{k,y-z}))$ plane of a relaxation from the saddle toward turbulence is added in orange, with a dot indicating its starting point. (d) Three dimensional plot of 151 reactive trajectories obtained by means of AMS in the domain of size $L_x \times L_z = 18 \times 12$ at $\beta = 80000$ (blue) as well as the relaxation from the fixed point (black) in the space formed by the logarithm of kinetic energy contained in components $\log(E_{k,x})$, $\log(E_{k,y-z})$ and turbulent fraction $F$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 12: Distribution of duration of reactive trajectories rescaled by their sample average $\tau$ and sample standard deviation $\sigma$. The trajectories are computed in a system of size $L_x \times L_z = 6 \times 4$ at Reynolds number $R = 500$ under forcing on components or divergence free forcing for several energy injection rates. A normalised Gumbel distribution (Eq. (21)) is added to the plots. (a) in linear scale. (b) in Logarithmic scale. 5320 trajectories were sampled for $\beta = 50000$ with forcing on $x$ component. 3220 trajectories were sampled for $\beta = 200000$ with forcing on $y$ and $z$ components. 1271 trajectories were sampled for $\beta = 110000$ with forcing on $x$ component. 1600 trajectories were sampled for $\beta = 200000$, 1399 were sampled for $\beta = 300000$ and 1055 were sampled for $\beta = 400000$ with divergence free forcing.</figcaption></figure>

saddle becomes degenerate [Rolland et al., 2016]). This distribution is also observed in systems not forced by noise [Rolland, 2022, Gomé et al., 2022]. In our case, we observe a very good match between the sampled distribution and a Gumbel, with small differences for $t \le \tau$ (which was also observed for instanton phenomenology in extended gradient systems [Rolland et al., 2016]).

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

# 4 Effect of energy injection rate

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

## 4.1 Probability of crossing and mean waiting time

In this section, we present the probability of crossing from a neighbourhood of the laminar flow to turbulent flow under a divergence free forcing as well as the corresponding mean first passage time. Both of them will be presented as a function of $\beta$ with a fixed spectrum shape. The mean first passage times are computed by means of AMS, following the procedure described in section 2.2.1. The probability of crossing is computed by means of AMS for all $\beta$ of interest as well as by means of DNS (following the procedure of § 2.2.3) up to $\beta = 3 \cdot 10^5$. These computations were performed in the MFU type system of size $L_x \times L_z = 6 \times 4$ at Reynolds number $R = 500$. The estimates are performed by averaging over the values of several realisations, while the error bars are constructed from the 66% interval of confidence obtained from the sample standard deviations. Up to $\beta = 2.8 \cdot 10^5$, $N = 200$ clones are used in AMS
<figure><img src="image_16.png" /><figcaption>Figure 13: From computations in system of size $L_x \times L_z = 6 \times 4$ at Reynolds number $R = 500$ under divergence free forcing, by means of AMS and up to $\beta = 3 \cdot 10^5$ by means of DNS. (a) : Probability to cross from the defined neighbourhood of the laminar baseflow to the turbulent flow as a function of $\beta$, an affine fit is performed for $\beta \ge 2.7 \cdot 10^5$ and is added to the plot. (b) : Mean first passage time before transition to the turbulent flow as a function of $\beta$, an affine fit is performed for $\beta \ge 2.7 \cdot 10^5$ and is added to the plot. (c) : average duration of reactive trajectories as a function of $\beta$. (d) : normalised standard deviation of the estimate of $\alpha$ as a function of $\beta$, the black dashed line indicates the ideal value 1.</figcaption></figure>

computations, $N = 100$ clones are used for $\beta = 2.9 \cdot 10^5$ and $\beta = 3 \cdot 10^5$ and $N = 50$ clones are used for computations at $3.1 \cdot 10^5 \le \beta \le 4 \cdot 10^5$. In these simulations, one clone was suppressed at each iteration. We performed further computations at $\beta = 5 \cdot 10^5$ using $N = 800$ clones and removing $N_c = 160$ clones per iteration and at $\beta = 6 \cdot 10^5$ using $N = 400$ clones and removing $N_c = 80$ clones per iteration. At least 10 realisations of AMS computations per control parameters are performed for $\beta < 4 \cdot 10^5$. This number goes up to more than 30 for $2.9 \cdot 10^5 \le \beta \le 3.3 \cdot 10^5$, while it is of the order of several units for $\beta \ge 4 \cdot 10^5$.

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

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

## 4.2 Duration of trajectories

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

### 4.2.1 Forcing on components

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

From all this we are drawn to address two matters

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

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

### 4.2.2 Divergence free forcing

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

## 5.1 Summary

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

## 5.2 Mechanical description of the paths

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

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

## 5.3 Statistical description of the paths

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

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

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

The comparison between the reactive trajectory and the nearest fixed point on the crossed separatrix is very often the focus of attention. This is because, as stated earlier, in systems which do not have a degenerate saddle [Bouchet and Touchette, 2012], in the limit of energy injection and/or forcing noise variance going to zero, the reactive trajectories are structured by instantons, which cross the separatrix right at a saddle point (see the sketch in figure 15 (a)). In said regime, one finds a structuration of the bistability: an exponential dependence of the mean first passage time and crossing probability in the noise variance, a concentration of trajectories in phase space and a specific distribution of reactive trajectories duration, to mention the most important features. Surprisingly enough, all these properties were found in our reactive trajectories.
This raises the question of why do we have this apparent contradiction. There can be several causes
for such a situation:

i) The paths sampled by means of AMS could be incorrect. This can happen when the reaction coordinate is not adapted. One then observes the apparent bias phenomenon: most sample reactive trajectories do not follow the most likely path [Bréhier et al., 2016, Lucente et al., 2022]. However, that type of behaviour would show a clear distinction between the trajectory durations computed by means of AMS and those computed by means of DNS as well as lead to a very large relative sample variance of the probability of crossing, which is not the case here. Moreover, when this phenomenon manifests itself, the reactive trajectories tend to select a path that still visit the neighbourhood of a saddle [Rolland and Simonnet, 2015].

ii) One can observe correctly computed reactive trajectories which are outside of the instanton regime if the relevant saddle is very flat [Rolland et al., 2016]. In that situation, one does observes an exponential growth of mean first passage time and an exponential decay of the crossing probability, which is our case. Meanwhile, one observes a Gaussian distribution of trajectory duration and a linear increase of their average duration and variance, which is not our case. Moreover, this regime still visits the neighbourhood of a saddle, which is not our case.

iii) When the saddle in the system can become degenerate, one can observe sub-instantons, that do not follow the expected path [Bouchet and Touchette, 2012]. However, in that case, the probability of crossing and the mean first passage time have a very distinct scaling in exponential of the square root of the noise variance, which is not our case.

iv) Finally, a peculiar situation has been reported following the computation of the instanton going from a non linear traveling wave to the laminar flow in two dimensional plane Poiseuille flow [Wan et al., 2015]. In that case, sketched in figure 15 (b), the instanton reaches the neighbourhood of the separatrix near a flow configuration termed a mediator state, travels along the separatrix toward the unstable fixed point then relaxes toward its destination. The peculiarity of that situation is that once the flow has reached the mediator point, a small numerical error in instanton calculation, or a small amount of noise in a stochastic system, can push the dynamics on the other side of the separatrix : the dynamics then travels near deterministically toward the final state B. Let us term $\mathcal{I}_{A \to B}$ the full action of the instanton going from A to B, $\mathcal{I}_{A \to m}$, the action of a trajectory following the instanton from A to the mediator state, $\mathcal{I}_{m \to s}$ the action of a trajectory following the instanton from the mediator state to the saddle and $\mathcal{I}_{m \to B}$ the action of a trajectory starting on the mediator and being slightly perturbed to cross the separatrix and then travel deterministically toward B. The action of deterministic paths is always zero: this means that in the limit of noise variance going to zero, the action of any path starting on the B side of the separatrix to B that has any chance of being observed is zero, because it will follow a deterministic path. This means that $\mathcal{I}_{A \to B} = \mathcal{I}_{A \to m} + \mathcal{I}_{m \to s} \simeq \mathcal{I}_{A \to m}$: the action of the instanton is that of traveling toward the mediator point then to the saddle and then relaxing. Moreover the part of the path along the separatrix has a negligible action $\mathcal{I}_{m \to s} \ll \mathcal{I}_{A \to m}$ because it is a near deterministic path in the stable manifold of the saddle. Similarly, the action of a path traveling toward the mediator along the instanton then receiving a small amount of noise to cross the separatrix there (as sketched in blue sketched in figure 15 (b)) is $\mathcal{I}_{A \sim B} = \mathcal{I}_{A \to m} + \mathcal{I}_{m \sim B} \simeq \mathcal{I}_{A \to m}$. The action $\mathcal{I}_{A \to B}$ and $\mathcal{I}_{A \sim B}$ are very similar and the corresponding probabilities and mean first passage time respectively go to zero and to infinity exponentially as very similar rates. At finite though small noise variance, the probability of observing a trajectory going along the instanton all the way becomes small compared to observing a trajectory going toward the mediator then crossing toward B there, because the instanton requires the trajectory not to be perturbed once the mediator has been reached. Therefore, computations by means of numerical simulation of the stochastic process (like DNS or AMS) would give the trajectory crossing near the mediator. Moreover, trajectories would still concentrate up to the mediator. However, the estimate of the probability and the mean first passage time would go like $\exp(\pm\beta\mathcal{I}_{A \to m})$ and thus be exponential and very close to the one predicted by the Freidlin Wentzel principle of large deviations.
<figure><img src="image_18.png" /><figcaption>Figure 15: (a) Sketch of the two metastable states **A** and **B** in phase space, separated by a separatrix (green curve), on which we have the transition state: the most accessible saddle point (red arrows indicating stable and unstable directions). In the limit of vanishing energy injection rate $\epsilon \to 0$, the reactive trajectory should follow an instanton (in cyan): on this first sketch, we represent an instanton that nears the separatrix only as it reaches the saddle along mostly unstable directions. For small but finite energy injection rate $\epsilon$, we would expect that trajectories concentrate around the instanton (in yellow). Instead, we systematically observe trajectories that concentrate around a path that crosses the separatrix away from the instanton and the most accessible fixed point (in blue). For these two types of paths: the full line indicates the average trajectory, while the dashed lines indicate where the zone that contains the majority of paths. (b) Modification of sketch (a) in the case where the instanton approaches the separatrix near a mediator state, then goes along the separatrix approaches the nearest fixed point and finally travels towards **B**, with a scenario similar to the one described by [Wan et al., 2015].</figcaption></figure>

In our case, the exponential dependence of $\alpha$ and $T$, the crossing of the separatrix away from the saddle and the concentration of paths are consistent with this scenario. We note however that the computed path do not concentrate on top of one another in planes computed using kinetic energy. This could be an effect of the finite noise variance in these quadratic spatially averaged quantities. Indeed, the mechanical scenario is not changed from one value of $\beta$ to another.

As far as we can deduce from the data presented here, scenario iv is the most similar to build up in plane Couette flow. If it is the case, this would mean that we would observe an actual instanton all along the trajectory, not just until the mediator or turning point only, for exceedingly small values of the noise. Moreover, if we follow this line of reasoning, we cannot rule out the possibility that there could be a mediator even closer to the laminar state, comparable to a minimal seed, that would only be reached at minuscule values of noise variance. It could be even beyond the reach of rare events simulations methods. The case of build up in plane Couette flow thus helped us identifying a different regime of convergence toward instantons, already pointed out in two dimensional plane Poiseuille flow [Wan et al., 2015] and can be linked to observations made in other systems where crossing away from the saddle and very flat dependence of the action on the crossing point in model of oceanic flows [Börner et al., 2023]. This behaviour could thus manifest itself in more complex bistable geophysical flows.

## Acknowledgment

The author thanks Institut PPrime, ENSMA and their computational facilities, where this work was started, for its hospitality. The author thanks the help of the computational resources of the PSMN platform of ENS de Lyon, as well as the hospitality of the laboratoire de Physique de l'ENS de Lyon, where part of the computation of this work have been performed. Additional computational resources
were granted by Mésocentre de Calcul Scientifique Intensif of the University of Lille and this work was granted access to the HPC resources of IDRIS under the allocation 2022-A0122A01741 made by GENCI. The author also thanks the support from CICADA (centre de calcul interactifs) (now mésocentre Azzurra) of University Nice Sophia-Antipolis where some early computations have been performed. The author finally thanks F. Bouchet, R. Börner, J.-C. Vassilicos and J.-P. Laval for inspiring discussions.

A Structure of the forcing noise

In section 2.1, we have introduced a forcing f in the Navier–Stokes questions for the departure to laminar plane Couette flow (Eq. (2)). We had stated that this forcing was prescribed using the spectrum of its correlation function (Eq. (3)), governed by the shape factor $\Gamma_{l,n_x,n_z}$ on component $l$ and in plane wavenumbers $n_x$ and $n_z$. We had mentioned the two families of forcing used to generate the build up trajectories, illustrated with a visualisation of $\Gamma_{l,n_x,n_z}$ (Fig. 2 (a,b)). In this appendix, we will give a quantitative description of the procedure followed to generate these two forcings:

* Forcing on components. In that case, we prescribe the shape factor of the Fourier transform of the correlation function: we rewrite $\Gamma_{l,n_x,n_z}$ using tensors $\sigma_l \gamma_{n_x} \gamma_{n_z}$ on Fourier component $n_x, n_z$. The tensor $\sigma_l$ scales the correlation function on each component: it can also be used to select a forcing on one component rather than the other, using $\sigma_l = 0$. Meanwhile, the tensors $\gamma$ prescribe the shape of the spectrum of the correlation function. They all have in common the property that $\gamma_{n_x} = 0$ if $n_x > \lfloor \frac{N_x}{2} \rfloor$, $\gamma_{n_z} = 0$ if $n_z > \lfloor \frac{N_z}{2} \rfloor$, $\gamma_{n_x} = 0$ if $n_x < \lfloor \frac{L_x}{L_c} \rfloor$, $\gamma_{n_z} = 0$ if $n_z < \lfloor \frac{L_z}{L_c} \rfloor$. This first restriction is used to avoid forcing aliased modes. We always used $L_c = 1.0$ for the largest scale of cut-off of the forcing.This second restriction ensures that the low wavelength of transitional wall turbulence are not directly forced. Using a cut-off wavenumber $N_r$, we set $\gamma_{n_x} = 1$ for $\lfloor \frac{L_x}{L_c} \rfloor \le n_x \le N_r$ and $\gamma_{n_z} = 1$ for $\lfloor \frac{L_z}{L_c} \rfloor \le n_z \le N_r$: we have a white spectrum at small wavenumbers. We then set $\gamma_{n_x} = \frac{N_x}{n_x}$ for $N_r \le n_x \le \lfloor \frac{N_x}{2} \rfloor$ and $\gamma_{n_z} = \frac{N_z}{n_z}$ for $N_r \le n_z \le \lfloor \frac{N_z}{2} \rfloor$: we have a red spectrum at larger wavenumbers. The forcing on component $l$ at wavenumber $n_x, n_z$ is then a Gaussian random number of variance $\frac{1}{\beta}$ times $(\sigma_l \gamma_{n_x} \gamma_{n_z})^2$. This means that $\Gamma_{l,n_x,n_z} = (\sigma_l \gamma_{n_x} \gamma_{n_z})^2$. We considered such a forcing on the three components (termed “iso”, with $\sigma_x = \sigma_y = \sigma_z$) or only on the streamwise component (termed “x” with $\sigma_x = 1$, $\sigma_y = \sigma_z = 0$). We only used it in MFU type systems of size $L_x \times L_z = 6 \times 4$ for the standard resolution. We have used several values of $N_r = 4$, 6 and 8 (which are specified when the results are presented.

* Divergence free forcing obtained from a curl. For this matter, we use a potential vector $A_{n_x,n_y}$, such that forcing on component $x$ at wavenumber $n_x$, $n_z$ is $-\frac{1}{\sqrt{\beta}} \frac{2\pi n_x}{L_z} A_{n_x,n_z}$, there is no forcing on component $y$ and the forcing on component $z$ is $\frac{1}{\beta} \frac{2\pi n_x}{L_x} A_{n_x,n_z}$, with $k_x = \frac{2\pi n_x}{L_x}$ and $k_z = \frac{2\pi n_z}{L_z}$. We then have $A_{n_x,n_z} = \gamma_{n_x} \gamma_{n_z} \xi$, with $\xi$ a Gaussian random number of variance 1. These factors are non zero in a set range of wavenumber: we set that $\gamma_{n_x} = 0$ if $n_x > \lfloor \frac{N_x}{2} \rfloor$, $\gamma_{n_z} = 0$ if $n_z > \lfloor \frac{N_z}{2} \rfloor$, $\gamma_{n_x} = 0$ if $n_x < \lfloor \frac{L_x}{L_c} \rfloor$, $\gamma_{n_z} = 0$ if $n_z < \lfloor \frac{L_z}{L_c} \rfloor$. We still set $A_{n_x,n_z} = \gamma_x \gamma_z$. This time we have $\gamma_x = \frac{1}{n_x}$ for $\lfloor \frac{L_x}{L_c} \rfloor \le n_x \le N_r$ and $\gamma_z = \frac{1}{n_z}$ for $\lfloor \frac{L_z}{L_c} \rfloor \le n_z \le N_r$. We then set $\gamma_{n_x} = \frac{N_x}{n_x}$ for $N_r \le n_x \le \lfloor \frac{N_x}{2} \rfloor$ and $\gamma_{n_z} = \frac{N_z}{n_z}$ for $N_r \le n_z \le \lfloor \frac{N_z}{2} \rfloor$. The Fourier transform of the correlation function $\Gamma_{x,n_x,n_z} = \frac{4\pi^2 n_x^2}{L_z^2} A_{n_x,n_z}^2$ is displayed in the $k_x, k_z$ plane in figure 2 (a), $\Gamma_{y,n_x,n_z} = 0$ and $\Gamma_{z} = \frac{4\pi^2 n_z^2}{L_x^2} A_{n_x,n_z}^2$ is displayed in the $k_x, k_z$ in figure 2 (b).

In the MFU type domain of size $L_x \times L_z = 6 \times 4$, for the standard resolution, we used the parameters $L_c = 1.0$ and $N_r = 4$. In that domain, for the finer resolution, we used the parameters $L_c = 1.0$ and $N_r = 4$. In the domain of intermediate size $L_x \times L_z = 18 \times 12$, we used the parameters $L_c = 1.0$ and $N_r = 12$. In the largest domain, of size $L_x \times L_z = 36 \times 24$, we used the parameters $L_c = 1.0$ and $N_r = 24$.
<table><caption>Table 1: Table of parameters $d_E$, $o_E$, $d_S$, $o_S$ used in systematic AMS simulations to define the reaction coordinate for the study of build up at Reynolds number $R$ in a system of size $L_x \times L_z$ using numerical resolution $N_x \times N_y \times N_z$ in space and time step $\Delta t$ in space, and a forcing spectrum controlled by $\beta$, $L_c$, $N_r$.</caption><thead><tr><th>L<sub>x</sub> × L<sub>z</sub></th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub></th><th>Δt</th><th>β<sub>min</sub></th><th>β<sub>max</sub></th><th>L<sub>c</sub></th><th>N<sub>r</sub></th><th>d<sub>E</sub></th><th>o<sub>E</sub></th><th>d<sub>S</sub></th><th>o<sub>S</sub></th><th>type</th></tr></thead><tbody><tr><td>6 × 4</td><td>32 × 27 × 24</td><td>0.05</td><td>1.5 × 10<sup>5</sup></td><td>6 × 10<sup>5</sup></td><td>1.0</td><td>4</td><td>-1.05</td><td>7.1</td><td>-1.11</td><td>6.4</td><td>div(f) = 0</td></tr><tr><td>6 × 4</td><td>32 × 27 × 24</td><td>0.05</td><td>1.0 × 10<sup>5</sup></td><td>1.1 × 10<sup>5</sup></td><td>1.0</td><td>4</td><td>-1.0</td><td>4.52</td><td>-1.0</td><td>1.44</td><td>x</td></tr><tr><td>6 × 4</td><td>32 × 27 × 24</td><td>0.05</td><td>1.0 × 10<sup>5</sup></td><td>4.0 × 10<sup>5</sup></td><td>1.0</td><td>4</td><td>-1.02</td><td>6.33</td><td>-1.03</td><td>4.46</td><td>iso</td></tr><tr><td>6 × 4</td><td>32 × 27 × 24</td><td>0.05</td><td>1.0 × 10<sup>5</sup></td><td>5.0 × 10<sup>5</sup></td><td>1.0</td><td>6</td><td>-1.02</td><td>6.64</td><td>-0.98</td><td>3.58</td><td>iso</td></tr><tr><td>6 × 4</td><td>32 × 27 × 24</td><td>0.05</td><td>1.0 × 10<sup>5</sup></td><td>5.5 × 10<sup>5</sup></td><td>1.0</td><td>8</td><td>-1.05</td><td>7.48</td><td>-1.12</td><td>6.30</td><td>iso</td></tr><tr><td>6 × 4</td><td>64 × 27 × 48</td><td>0.05</td><td>1.0 × 10<sup>5</sup></td><td>4.0 × 10<sup>5</sup></td><td>1.0</td><td>4</td><td>-1.01</td><td>6.47</td><td>-1.02</td><td>4.30</td><td>iso</td></tr><tr><td>6 × 4</td><td>128 × 35 × 96</td><td>0.025</td><td>1.0 × 10<sup>5</sup></td><td>3.0 × 10<sup>5</sup></td><td>1.0</td><td>4</td><td>-1.02</td><td>6.58</td><td>-1.04</td><td>4.63</td><td>iso</td></tr></tbody></table>

<table><caption>Table 2: Table of parameters $E_A$, $E_B$ and $\sigma$ used in non systematic AMS simulations to define the reaction coordinate for the study of build up at Reynolds number $R$ in a system of size $L_x \times L_z$ using numerical resolution $N_x \times N_y \times N_z$ in space and time step $\Delta t$ in space, and a spectrum of the divergence free forcing controlled by $\beta$, $L_c$ and $N_r$.</caption><thead><tr><th>R</th><th>L<sub>x</sub> × L<sub>z</sub></th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub></th><th>Δt</th><th>β</th><th>L<sub>c</sub></th><th>N<sub>r</sub></th><th>E<sub>A</sub></th><th>E<sub>B</sub></th><th>σ</th></tr></thead><tbody><tr><td>500</td><td>6 × 4</td><td>64 × 47 × 48</td><td>0.0125</td><td>5.5 × 10<sup>5</sup></td><td>1.0</td><td>4</td><td>0.0096</td><td>0.06</td><td>0.0021</td></tr><tr><td>500</td><td>18 × 12</td><td>96 × 27 × 72</td><td>0.05</td><td>62500</td><td>1.0</td><td>12</td><td>0.0117</td><td>0.06</td><td>0.0034</td></tr><tr><td>500</td><td>18 × 12</td><td>96 × 27 × 72</td><td>0.05</td><td>80000</td><td>1.0</td><td>12</td><td>0.0081</td><td>0.06</td><td>0.0023</td></tr><tr><td>500</td><td>36 × 24</td><td>192 × 27 × 144</td><td>0.05</td><td>30000</td><td>1.0</td><td>24</td><td>0.0123</td><td>0.0477</td><td>0.0028</td></tr></tbody></table>

## B Parameter used in definition of reaction coordinate in each set of simulations

In section 2.2.2, we presented the set up of the reaction coordinates (Eq. (11)), using three free parameters $E_A$ to define the laminar state $A$, $\sigma$ to define the hypersurface $C$ that surrounds it and $E_B$ to define the turbulent state (Fig. 3 (a)). Most AMS computations used a systematic definition of these three parameters with $E_B = 0.06$, and the parameters controlling $E_A$ and $\sigma$ reported in table 1

Some AMS computations did not use a systematic definition for $E_A$, $\sigma$ and $E_B$, because the value of $\beta$ was too low for the system to be in the regime where $E$ is linear in $\epsilon$. For those cases, the control parameters used are reported in table 2.

## C Convergence of AMS computations

In this appendix, in complement to the precision tests of section 4.1 (Fig. 13 (d)), we will check that the number of clones used in AMS computation was sufficient. For this matter, we present convergence results with the number of clones for the two types of forcings in domains of size $L_x \times L_z = 6 \times 4$ at Reynolds number $R = 500$ using standard resolution $N_x = 32$, $N_y = 27$, $N_z = 24$ and $\Delta t = 0.05$. Isotropic forcing was added using $\beta = 4 \cdot 10^5$ and $\beta = 5 \cdot 10^5$ with $N_r = 6$, and $L_c = 1.0$. Divergence free forcing was added using $\beta = 2.5 \cdot 10^5$ and $3 \cdot 10^5$ using $L_c = 1$ and $N_r = 4$. In all four cases, the crossing probabilities are not that small, so that the reactive trajectories and their properties can be computed by DNS with a high enough degree of precision.

We compare the probability $\alpha$ in logarithmic scale computed by means of AMS and DNS in figure 16 (a). For all four examples considered, we observe the overlap of the DNS and AMS 66% intervals of confidence, obtained using $\pm$ the sample standard deviation, as soon as $N \ge 20$. This indicates that no error are performed in these AMS computation if the number of clones is sufficient. For $N = 5$ and
<figure><img src="image_19.png" /><figcaption>Figure 16: In a domain of size $L_x \times L_z = 6 \times 4$ at Reynolds number $R = 500$ for four different forcings (see § C for details). (a) Crossing $\alpha$ probability as a function of the number of clones, computed by AMS or DNS. The error bars correspond to the 66% confidence interval of AMS computations, the dotted line corresponds to the 66% interval of confidence of DNS computations. (b) normalised standard deviation (Eq. (22)) of the crossing probability $\alpha$ as a function of the number of clones, compared to 1. (c) Average duration $\tau$ of reactive trajectories as a function of the number of clones, computed by means of AMS and DNS for the same cases and the same line codes as (a).</figcaption></figure>

$N = 10$, we observe underestimates of $\alpha$ which are expected. Due to the increased rarity of the event for $\beta = 3 \cdot 10^5$ with divergent free forcing, the expensive computations means that the dataset is smaller, which explains the wider confidence intervals for DNS estimates.

We then present the normalised standard deviation of the estimate of $\alpha$ in the case where $N_c = 1$ and in the case where $N_c > 1$

$$ \sigma_{norm} = \frac{\sqrt{\langle \alpha^2 \rangle_o - \langle \alpha \rangle_o^2} \sqrt{N}}{\alpha \sqrt{|\log(\alpha)|}}, \quad \sigma_{norm,N_c} = \frac{\sqrt{\langle \alpha^2 \rangle_o - \langle \alpha \rangle_o^2} \sqrt{N}}{\alpha \sqrt{\left(\langle \kappa \rangle_o \frac{N_c}{N-N_c} + \frac{1-\langle r \rangle_o}{\langle r \rangle_o}\right)}}. \quad (22) $$

where $\langle \kappa \rangle_o$ is the average number of iterations over AMS runs and $\langle r \rangle_o$ is the average proportion of clones reaching $\mathcal{B}$ over AMS runs. It has been demonstrated that this quantity is larger than or equal to one, and it is one when the committor is used as reaction coordinate [Cérou et al., 2019]. In parallel, it has also been shown that the closer a reaction coordinate is to the committor, the most precise the AMS estimates are when a finite number of runs are used. As a consequence, this quantity is used as an estimate of the quality of the computation because the better the reaction coordinate, the smaller is $\sigma$. A very large $\sigma$ would be that the estimates of $\alpha$ and $T$ are incorrect, or even that the paths are incorrectly selected [Rolland and Simonnet, 2015, Lucente et al., 2022]. In that case what is typically observed in the selection of shorter though less probable paths. In the two test cases considered here (Fig. 16 (c)), we can see that $\sigma$ is slightly larger than one. This ensures correct estimates, but also tells us that there is room for improvement of the reaction coordinate.

We finally test the computation of the average duration of reactive trajectories. The result of the AMS computations of $\tau$ as a function of $N$ in our four cases is presented in figure 16 (c). We add the result of DNS computation and the 66% interval of confidence to this figure. We note the overlap of the intervals of confidence in all four cases as soon as $N \ge 10$, indicating that no discernable bias can be observed in the estimates of $\tau$ by means of AMS. If those estimates are averaged over enough repetitions, they are precise enough in those cases. In conclusion to this appendix, we note that we do not discern any irreducible error in the computations performed by means of AMS, even if the number of clones used is relatively modest (at least 20). In those case, a larger number of repetitions of the computations is required to compute said estimates. This confirm the observation performed in section 4.1 for several values of $\beta$ (and the corresponding $\alpha$) with only one clone number per case.
References

[Baars et al., 2021] Baars, S., Castellana, D., Wubs, F. W., and Dijkstra, H. A. (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. c., Chiffaudel, A., Daviaud, F. c., Dubrulle, B., Marié, L., Ravelet, F., Bourgoin, M., Odier, P., Pinton, J.-F. c., and Volk, R. (2007). Magnetic field reversals in an experimental turbulent dynamo. *Eur. Phys. Lett.*, 77(5):59001.

[Börner et al., 2023] Börner, R., Deeley, R., Römer, R., Grafke, T., Lucarini, V., and Feudel, U. (2023). Saddle avoidance of noise-induced transitions in multiscale systems. *arXiv preprint arXiv:2311.10231*.

[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.

[Bouchet et al., 2019] Bouchet, F., Rolland, J., and Simonnet, E. (2019). Rare event algorithm links transitions in turbulent flows with activated nucleations. *Physical review letters*, 122(7):074502.

[Bouchet and Touchette, 2012] Bouchet, F. and Touchette, H. (2012). Non-classical large deviations for a noisy system with non-isolated attractors. *Journal of Statistical Mechanics: Theory and Experiment*, 2012(05):P05028.

[Bréhier et al., 2016] Bréhier, C.-E., Gazeau, M., Goudenège, L., Lelièvre, T., and Rousset, M. (2016). Unbiasedness of some generalized adaptive multilevel splitting algorithms. *The Annals of Applied Probability*, 26(6):3559–3601.

[Brzeźniak et al., 2015] Brzeźniak, Z., Cerrai, S., and Freidlin, M. (2015). Quasipotential and exit time for 2d stochastic navier-stokes equations driven by space time white noise. *Probability Theory and Related Fields*, 162(3-4):739–793.

[Cérou et al., 2019] Cérou, F., Delyon, B., Guyader, A., and Rousset, M. (2019). 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., Lelièvre, 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., Lelièvre, T., and Pommier, D. (2011). A multiple replica approach to simulate reactive trajectories. *The Journal of chemical physics*, 134(5):054108.

[Chantry and Schneider, 2014] Chantry, M. and Schneider, T. M. (2014). Studying edge geometry in transiently turbulent shear flows. *Journal of fluid mechanics*, 747:506–517.

[Cherubini et al., 2011] Cherubini, S., De Palma, P., Robinet, J.-C., and Bottaro, A. (2011). The minimal seed of turbulent transition in the boundary layer. *Journal of Fluid Mechanics*, 689:221.

[Ducimetière et al., 2022] Ducimetière, Y.-M., Boujo, E., and Gallaire, F. (2022). Weakly nonlinear evolution of stochastically driven non-normal systems. *Journal of Fluid Mechanics*, 951:R3.

[Farrell and Ioannou, 1993a] Farrell, B. F. and Ioannou, P. J. (1993a). Optimal excitation of three-dimensional perturbations in viscous constant shear flow. *Physics of Fluids A: Fluid Dynamics*, 5(6):1390–1400.
[Farrell and Ioannou, 1993b] Farrell, B. F. and Ioannou, P. J. (1993b). Stochastic forcing of the linearized navier–stokes equations. *Physics of Fluids A: Fluid Dynamics*, 5(11):2600–2609.

[Freidlin and Wentzell, 1998] Freidlin, M. I. and Wentzell, A. D. (1998). Random perturbations. In *Random perturbations of dynamical systems*, pages 15–43. Springer.

[Gibson et al., 2008] Gibson, J. F., Halcrow, J., and Cvitanović, P. (2008). Visualizing the geometry of state space in plane couette flow. *J. Fluid Mech.*, 611:107–130.

[Gomé et al., 2022] Gomé, S., Tuckerman, L. S., and Barkley, D. (2022). Extreme events in transitional turbulence. *Philosophical Transactions of the Royal Society A*, 380(2226):20210036.

[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.

[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 and Kida, 2001] Kawahara, G. and Kida, S. (2001). Periodic motion embedded in plane couette turbulence: regeneration cycle and burst. *Journal of Fluid Mechanics*, 449:291–300.

[Kim and Durbin, 1988] Kim, H.-J. and Durbin, P. A. (1988). Investigation of the flow between a pair of circular cylinders in the flopping regime. *Journal of Fluid Mechanics*, 196:431–448.

[Kreilos and Eckhardt, 2012] Kreilos, T. and Eckhardt, B. (2012). Periodic orbits near onset of chaos in plane couette flow. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 22(4).

[Liu et al., 2021] Liu, T., Semin, B., Klotz, L., Godoy-Diana, R., Wesfreid, J. E., and Mullin, T. (2021). Decay of streaks and rolls in plane couette–poiseuille flow. *Journal of Fluid Mechanics*, 915.

[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.

[Lucarini and Bódai, 2017] Lucarini, V. and Bódai, T. (2017). Edge states in the climate system: exploring global instabilities and critical transitions. *Nonlinearity*, 30(7):R32.

[Lucarini and Bódai, 2019] Lucarini, V. and Bódai, T. (2019). Transitions across melancholia states in a climate model: Reconciling the deterministic and stochastic points of view. *Physical review letters*, 122(15):158701.

[Lucente et al., 2022] Lucente, D., Rolland, J., Herbert, C., and Bouchet, F. (2022). Coupling rare event algorithms with data-based learned committor functions using the analogue markov chain. *Journal of Statistical Mechanics: Theory and Experiment*, 2022(8):083201.
[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.

[Monokrousos et al., 2011] Monokrousos, A., Bottaro, A., Brandt, L., Di Vita, A., and Henningson, D. S. (2011). Nonequilibrium thermodynamics and the optimal path to turbulence in shear flows. *Physical review letters*, 106(13):134502.

[Nusse and Yorke, 1989] Nusse, H. E. and Yorke, J. A. (1989). A procedure for finding numerical trajectories on chaotic saddles. *Physica D: Nonlinear Phenomena*, 36(1-2):137-156.

[Onsager, 1938] Onsager, L. (1938). Initial recombination of ions. *Physical Review*, 54(8):554.

[Onsager and Machlup, 1953] Onsager, L. and Machlup, S. (1953). Fluctuations and irreversible processes. *Physical Review*, 91(6):1505.

[Philip and Manneville, 2011] Philip, J. and Manneville, P. (2011). From temporal to spatiotemporal dynamics in transitional plane couette flow. *Physical Review E*, 83(3):036308.

[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.

[Rabin et al., 2012] Rabin, S. M. E., Caulfield, C.-C. P., and Kerswell, R. R. (2012). Triggering turbulence efficiently in plane couette flow. *Journal of Fluid Mechanics*, 712:244.

[Rigas et al., 2021] Rigas, G., Sipp, D., and Colonius, T. (2021). Nonlinear input/output analysis: application to boundary layer transition. *Journal of Fluid Mechanics*, 911.

[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, 2022] Rolland, J. (2022). Collapse of transitional wall turbulence captured using a rare events algorithm. *Journal of Fluid Mechanics*, 931.

[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.

[Schmid and Henningson, 1994] Schmid, P. J. and Henningson, D. S. (1994). Optimal energy density growth in hagen-poiseuille flow. *Journal of Fluid Mechanics*, 277:197-225.

[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.

[Schneider et al., 2008] Schneider, T. M., Gibson, J. F., Lagha, M., De Lillo, F., and Eckhardt, B. (2008). Laminar-turbulent boundary in plane couette flow. *Physical Review E*, 78(3):037301.

[Schneider et al., 2010] Schneider, T. M., Marinc, D., and Eckhardt, B. (2010). Localized edge states nucleate turbulence in extended plane couette cells. *Journal of Fluid Mechanics*, 646:441-451.
[Simonnet, 2016] Simonnet, E. (2016). Combinatorial analysis of the adaptive last particle method. *Statistics and Computing*, 26(1):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.

[Spangler and Wells Jr, 1968] Spangler, J. G. and Wells Jr, C. S. (1968). Effects of freestream disturbances on boundary-layer transition. *AIAA Journal*, 6(3):543–545.

[Toh and Itano, 2003] Toh, S. and Itano, T. (2003). A periodic-like solution in channel flow. *Journal of Fluid Mechanics*, 481:67–76.

[Touchette, 2009] Touchette, H. (2009). The large deviation approach to statistical mechanics. *Physics Reports*, 478(1-3):1–69.

[Trefethen et al., 1993] Trefethen, L. N., Trefethen, A. E., Reddy, S. C., and Driscoll, T. A. (1993). Hydrodynamic stability without eigenvalues. *Science*, 261(5121):578–584.

[Van Erp et al., 2003] Van Erp, T. S., Moroni, D., and Bolhuis, P. G. (2003). A novel path sampling method for the calculation of rate constants. *The Journal of chemical physics*, 118(17):7762–7774.

[Waleffe, 1997] Waleffe, F. (1997). On a self-sustaining process in shear flows. *Physics of Fluids*, 9(4):883–900.

[Wan, 2013] Wan, X. (2013). A minimum action method for small random perturbations of two-dimensional parallel shear flows. *Journal of Computational Physics*, 235:497–514.

[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.

[Wan et al., 2015] Wan, X., Yu, H., and Weinan, E. (2015). Model the nonlinear instability of wall-bounded shear flows as a rare event: a study on two-dimensional poiseuille flow. *Nonlinearity*, 28(5):1409.

[Willis and Kerswell, 2009] Willis, A. P. and Kerswell, R. R. (2009). Turbulent dynamics of pipe flow captured in a reduced model: puff relaminarization and localized 'edge' states. *Journal of Fluid Mechanics*, 619:213–233.