Joran Rolland
email: ‡joran.rolland@centralelille.fr

Collapse of transitional wall turbulence captured using a rare events algorithm

Keywords: 47, 27-Cn: Transition to turbulence & 47, 27-nd: Channel flow & 05, 70-Ln: Non equilibrium and irreversible thermodynamics

This text presents one of the first successful applications of a rare events method for the study of multistability in a turbulent flow without stochastic energy injection. The trajectories of collapse of turbulence in plane Couette flow, as well as their probability and rate of occurrence are systematically computed using Adaptive Multilevel Splitting (AMS). The AMS computations are performed in a system of size Lx × Lz = 24 × 18 at Reynolds number R = 370 with an acceleration by a factor O(10) with respect to DNS and in a system of size Lx × Lz = 36 × 27 at Reynolds number R = 377 with an acceleration by a factor O(10 3 ). The AMS results are validated with a comparison to DNS in the system of size Lx × Lz = 24 × 18. Visualisations in both systems indicate that turbulence collapses because the self sustaining process of turbulence fails locally. The streamwise vortices decay first in streamwise elongated holes, leaving streamwise invariant streamwise velocity tubes that experience viscous decay. These holes then extend in the spanwise direction. The examination of more than a thousand of trajectories in the (Ec,x = u 2

x /2 d 3 x, Ec,y-z = (u 2 y /2 + u 2 z /2) d 3 x) plane in the system of size Lx × Lz = 24 × 18 confirms the faster decay of streamwise vortices and shows concentration of trajectory. This hints at an instanton phenomenology in the large size limit. The computation of turning point states, beyond which laminarisation is certain, confirms the hole formation scenario and shows that it is more pronounced in larger systems. Finally, the examination of non reactive trajectories, where a hole opens then closes, indicates that hole opening and closing are distinct processes. Both the vortices and the streaks reform concomitantly when the laminar holes close.

Introduction

Many turbulent flows of aerodynamical or geophysical interest are not homogeneous nor isotropic. They often display several possible turbulent flow configurations and rare switches between these flow configurations. For these reasons these flows are often termed multistable. Examples of multistability in turbulent flows include turbulent dynamos [START_REF] Berhanu | Magnetic field reversals in an experimental turbulent dynamo[END_REF], turbulent convection [START_REF] Podvin | Precursor for wind reversal in a square rayleigh-bénard cell[END_REF], bluff body wakes [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 1. global modes and bi-stability[END_REF], jets in the wake of a pair of cylinders [START_REF] Kim | Investigation of the flow between a pair of circular cylinders in the flopping regime[END_REF] and barotropic atmospheric-type jets [START_REF] Bouchet | Rare event algorithm links transitions in turbulent flows with activated nucleations[END_REF]. In order to understand these turbulent flows, uncovering what drive the switches is as important as explaining the mechanisms maintaining each metastable configuration. These switches are characterised by the mean first passage time before a change of configuration occurs. The switches often take place through the same chain of events, termed a transition path or a reactive trajectory, a notion originating from kinetic chemistry [START_REF] Metzner | Illustration of transition path theory on a collection of simple examples[END_REF]. It is fairly difficult to compute mean first passage times and transition paths, either experimentally or numerically for two main reasons. a) One firstly has to deal with the very large number of degrees of freedom of turbulent flows, particularly in numerical simulation. These complex flows are often difficult to simulate even by means of Large Eddy Simulations.

b) One secondly has to deal with the extremely long waiting times between each event. These waiting times are several orders of magnitude larger than the duration of a realisation of a switch and even longer than the typical eddy turnover time [START_REF] Kim | Investigation of the flow between a pair of circular cylinders in the flopping regime[END_REF]. This means that the cost of sampling more than a few events is prohibitive by classical means (see table in [START_REF] Bouchet | Rare event algorithm links transitions in turbulent flows with activated nucleations[END_REF]).

In order to propose a method to systematically study multistability, one first needs a turbulent system which has fewer effective degrees of freedom but is still complex enough to display multistability, with a moderate need for extrinsic stochastic forcing. One can thus temporarily bypass problem a). A transitional wall flow such as plane Couette flow can provide such a situation (see Fig. 1). Unlike thermal convection, for instance, a wall flow like plane Couette flow is linearly stable for all Reynolds numbers [Romanov, 1973]. Meanwhile transitional turbulence can exist at moderate Reynolds numbers, albeit transiently [START_REF] Eckhardt | Dynamical systems and the transition to turbulence in linearly stable shear flows[END_REF], Schmiegel and Eckhardt, 1997, Bottin and Chaté, 1998, Eckhardt et al., 2007]. Meanwhile, turbulence can build up from laminar flow under a finite amplitude forcing [START_REF] Wan | A dynamic-solver-consistent minimum action method: With an application to 2d navier-stokes equations[END_REF], Rolland, 2018, Liu et al., 2020], or somewhat equivalently from a finite amplitude initial condition [START_REF] Faisst | Sensitive dependence on initial conditions in transition to turbulence in pipe flow[END_REF]. This gives a first type of multistability, where the transitional turbulence of plane Couette flow can collapse down to laminar flow under its own fluctuations, and go back to turbulence if it is forced. In a way this is a purely temporal view, this image should be completed by noting that transitional wall flow turbulence tends to be localised (the type of organisation depends on the flow geometry). At moderate Reynolds numbers, wall flows where transitional turbulence extends in one dimension (such as Hagen-Poiseuille flow [START_REF] Moxey | Distinct large-scale turbulent-laminar states in transitional pipe flow[END_REF], or tilted plane Couette flow [START_REF] Gomé | Statistical transition to turbulence in plane channel flow[END_REF]) can display splitting: the turbulent puff elongates and splits in half. This can lead to an effective extension of the area occupied by turbulence. This can provide a third type of multistability events in some flow configurations. Aside from the relevance of wall flows for stable boundary layers, canonical flows such as plane Couette can thus provide a good laboratory to test rare events methods on turbulent flows. These multistable turbulent flows are not so complex that they cannot be simulated by Direct Numerical Simulations (DNS).

Even if the number of degrees of freedom is reduced, as is the case for transitional wall, one still has to deal with problem b). Transitional turbulence, like other multistable system, displays very long waiting time before a multistability event occurs [START_REF] Eckhardt | Turbulence transition in pipe flow[END_REF], Bottin and Chaté, 1998, Gomé et al., 2020]. This means that even for the study of turbulence collapse in plane Couette flow, the use of DNS to sample many events comes with a prohibitive computational cost. There are of course alternate methods to study rare events. Many of them originate from the study of kinetic chemistry and borrow much of its vocabulary. All these methods aim at performing the computation of both the reactive trajectories and the mean waiting time before they occur. These methods can be divided in two main families.

• On the one hand, one finds mostly theoretical, optimisation methods, which are applied to systems where a stochastic forcing is clearly identified. In that case, in the limit where the variance of said forcing goes to zero, all the transition paths concentrate around an instanton [Touchette, 2009] § 6. The instanton represents the most probable transition path and can be computed using action minimisation [Grafke andVanden-Eijnden, 2019, Wan andYu, 2017]. One fundamental property of said instantons is that they evolve from a first multistable state toward a saddle point of the deterministic dynamics under the action of the noise. The instanton then evolves freely from the saddle point toward the second multistable state. The mean first passage time can be estimated using the result of said action minimisation and can loosely be thought of as depending on the distance between the first metastable state and the saddle point. In gradient stochastic systems, this yields the celebrated Eyring-Kramers formula, also known as the Arrhenius law [START_REF] Hänggi | Reaction-rate theory: fifty years after kramers[END_REF], such result can be extended to non gradient systems [START_REF] Bouchet | Generalisation of the eyring-kramers transition rate formula to irreversible diffusion processes[END_REF]. While these methods give a lot of qualitative results on the structure of transition paths, they can prove tricky to implement for fluid flows [START_REF] Wan | A dynamic-solver-consistent minimum action method: With an application to 2d navier-stokes equations[END_REF]. A key property of these formula and methods is that they are formulated as large deviations [Touchette, 2009]. This means that they correspond to an asymptotic writing of the small parameter times the logarithm of probability, a probability density function or a rate of probability, which has a finite value independent of said small parameter.

• On the other hand, one finds mostly numerical cloning methods, which use the actual fluctuating dynamics of the system and push them toward realisations of the reactive trajectories. These methods compute the transition paths using N clone dynamics of the system and apply a mutation selection procedure to compute the reactive trajectories. One such method is termed Adaptive Multilevel Splitting (AMS) [Cérou andGuyader, 2007, Cérou et al., 2019]. AMS and its variants have been successfully used to compute reactive trajectories and extreme events in kinetic chemistry [START_REF] Lopes | Analysis of the adaptive multilevel splitting method on the isomerization of alanine dipeptide[END_REF] theoretical physics models [START_REF] Rolland | Computing transition rates for the 1-d stochastic ginzburg-landau-allen-cahn equation for finite-amplitude noise with a rare event algorithm[END_REF], models of transitional flows [Rolland, 2018], idealised atmospheric flows [START_REF] Bouchet | Rare event algorithm links transitions in turbulent flows with activated nucleations[END_REF], Simonnet et al., 2020]. Some variants have been applied to the study of extreme two dimensional turbulent wakes [START_REF] Lestang | Computing return times or return periods with rare event algorithms[END_REF], Lestang et al., 2020] and oceanic flow reversals [START_REF] Baars | Application of adaptive multilevel splitting to high-dimensional dynamical systems[END_REF].

In practice, these two types of methods are often used hand in hand. Once a small parameter that controls the effective noise variance has been identified, the qualitative insight from theory is used to guide numerical studies and propose manners in which results can be presented. Applying such a program to the collapse of turbulence in plane Couette flow is actually not that straightforward. Because energy is not injected in the flow through a stochastic forcing whose variance decays, the small parameter is not readily identified and cloning methods cannot be applied as such [START_REF] Lestang | Numerical study of extreme mechanical force exerted by a turbulent flow on a bluff body by direct and rare-event sampling techniques[END_REF]. Applying basic cloning rules leads to so called extinction: the flow does not separate trajectories from one another and the method does not manage to create reactive trajectory to would contain an excursion far enough from the starting metastable state. A first goal is therefore to propose a modified cloning method that can bypass this problem and can at least succeed in computing reactive trajectories faster than a DNS would. This is the purpose of anticipated AMS, which is presented and used in this text. Once the reactive trajectories are computed, if they display concentration around a typical transition path, one needs to make sense of this concentration. For this purpose, one can use both AMS and DNS to identify the right small parameter from the study of probability density functions, probability of transition, mean first passage times etc. We will keep this in mind in our study. We also note that there are examples very relevant reactive trajectories (see [START_REF] Rolland | Computing transition rates for the 1-d stochastic ginzburg-landau-allen-cahn equation for finite-amplitude noise with a rare event algorithm[END_REF] for instance) that are not necessarily of the instanton type. As a consequence, one does not always have to try to force the results in a large deviations framework.

We present the study of the collapse of transitional turbulence of plane Couette flow in the following manner. We first present Anticipated Adaptive multilevel splitting in section 2.1. We then remind the configuration of plane Couette flow in section 2.2.1, how initial conditions are generated in section 2.2.2 and what reaction coordinates are used in section 2.2.3. The DNS which are used as reference are presented in the next section ( § 2.3). The reaction coordinate used to compute reactive trajectories is presented in section 2.2.3. We then present the systematic comparison of reactive trajectories computed by AMS and DNS in a system of size L x × L z = 24 × 18 ( § 3.1). We perform the validation of the computation probability of crossing and mean first passage time in this system in section 3.2. AMS is then applied to the computation of very rare trajectories and laminar hole formation in section 4. These results are finally discussed together in the conclusion ( § 5).

Method

AMS

Before presenting the principle of anticipated AMS, let us first give a formal phase space description of the rare events we will study in this text. Let us sketch the collapse of turbulence in figure 2 (a) and consider the set A, a neighbourhood of the turbulent flow in phase space, and the set B, a neighbourhood of the laminar flow . A realisation of the dynamics which starts in A fluctuates around it, has several excursions out of C, an hypersurface closely surrounding A, and eventually crosses C and reaches B before coming back to 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 A, crosses C and reaches B before A. Precise definitions of sets A, B and hypersurface C for collapse will be given in sections 3 and 4, based on reaction coordinates defined in section 2.2.3. We then give a brief overview of the variant of AMS, termed Anticipated Adaptive Multilevel Splitting, which was used for the computation presented in this text (see [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF], Rolland and Simonnet, 2015, Bréhier et al., 2016, Rolland, 2018, Cérou et al., 2019, Lestang et al., 2020] for more details on the general methods). All variants of AMS use a reaction coordinate (or observable) φ(u) = Φ(t), which measures the position of the flow relatively to starting set A and arrival set B. All variants run N clone dynamics of the system to compute iteratively at least N -N c > 0 reactive trajectories going from a hypersurface C close to set A to the set B. The algorithm is sketched in figure 2

(b) and proceeds in the following manner

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

• The algorithm runs the stages of mutation selection. At each stage, the clones 1 ≤ i ≤ N c (with N c < N ), ordered by max t Φ i (t), are suppressed. They are branched on N c clones drawn uniformly out of the N -N c other clones at level max t Φ Nc (t), the maximal value of Φ reached by the suppressed trajectories. The branched trajectories then follow their natural dynamics until they reach either A or B, with a new realisation of the noise.

The algorithm stops after k iterations when N r ≥ N -N c + 1 clones trajectories have reached B. This yields an estimator of the probability α of reaching B before A [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF], and the corresponding mean first passage time T [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF]

α = r 1 - N c N k =α o , T = 1 α -1 (t 1 + τ ) + (t 1 + τ ) o , ( 1
)
where τ is the mean duration of reactive trajectories, t 1 is the mean duration to go from A to C and τ is the mean duration of non reactive trajectories, computed in each AMS run. The notation

• o = o l=1
• l corresponds to an average over o independent AMS runs [Cérou andGuyader, 2007, Rolland, 2018].

We will often record the velocity field, noted u last and termed the last state at the last stage, that corresponds to max t Φ Nc (t) during the last stage of the algorithm. It often gives a precise idea of the turning point in reactive trajectories. Before the flow visits the neighbourhood of that state, returning towards turbulence is more likely, beyond that point, relaminarising becomes more likely. In systems which correspond to a simple deterministic part forced by noise, that state actually corresponds to the saddle point of the deterministic part of the dynamics crossed by the instanton in the limit of the noise variance going to 0. This has been verified for models with few degrees of freedom and the one dimensional Ginzburg-Landau equation (not shown here). It can be used to educe an effective saddle between two multistable states [START_REF] Simonnet | Multistability and rare spontaneous transitions between climate and jet configurations in a barotropic model of the jovian midlatitude troposphere[END_REF]. Note that it does not a priori corresponds to an actual saddle of the Navier-Stokes equations (as computed by dichotomy or other methods [Schneider et al., 2007, Willis andKerswell, 2009]).

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

Note that in deterministic systems, we more often than not add a small background noise which helps the separation of trajectories after branching. When we branch on a trajectory, we change the realisation of this noise. This background noise will be noted as an additional stochastic forcing in the equation and its properties will be given. In AMS computations, we will use a simple "tick" a forcing exerted just at branching, white in space and time with variance σ = 10 -10 so as not to perturb the laminar flow too much. The trade-off is that we compute trajectory properties with a small error : T + δT , α + δα and τ + δτ . We will comment on the visible effects of this additional force on the trajectories and their properties.

Schematically, we can apply AMS to two types of systems. On the one hand, we find systems with a large time scale separation between some fast fluctuating degrees of freedom and slower degrees of freedom which represent the main features of the flow travelling between A and B. This is often the case for stochastically forced systems. In these systems, two slightly different initial conditions will quickly separate and the odds of creating an excursion toward B instead of A by slightly changing the noise realisation at a branching are non negligible. On the other hand, we find systems with absolutely no clear time scale separation between degrees of freedom. This is often the case of purely deterministic systems. Two slightly different initial conditions do not separate until it is too late (they both reach A). We can find situations where, no matter the structure of the small perturbation (typically at a branching), the odds of creating a further excursion toward B instead of A can be exceedingly small. This is especially the case if we perturb at the peak of an existing fluctuation. If we apply basic AMS to this second type of system, where we branch at Φ Ñc , we run the risk of a so called extinction [START_REF] Lestang | Numerical study of extreme mechanical force exerted by a turbulent flow on a bluff body by direct and rare-event sampling techniques[END_REF]. This occurs when all trajectories have the same maximum which is strictly outside the arrival set ∀1 ≤ i ≤ N, max t Φ i = Φ ext < Φ(∂B). The algorithm does not manage to proceed any further1 . In order to bypass this limitation, we can perform anticipated branching, that is to say branch the new trajectories at Φ b < Φ Nc . In that case, it may be necessary to reiterate the branching several times in order to ensure that the branched trajectories have max t Φ > Φ Ñc . In figure 2 (c), we give two examples of relations Φ b (Φ Nc ) that were tested for place Couette flow. Each one is adapted to a given situation. The converging anticipation is mostly used in this article. The point is to take advantage of higher mixing and faster separation of trajectories that take place when the flow is closer to the fully turbulent state. It converges to the line so as not to lose too much computational time rerunning trajectories when the flow is very close to turbulence collapse. The saturated anticipation takes a very different point of view. It has proved very efficient in very small domains where the collapse of turbulence is very well described by transient chaos. It uses the fact that there is much more mixing when the flow is close to the turbulent state, but that mixing gradually stops during excursions. If a trajectory is on the wrong track, there is no derailing it at large Φ. In this type of systems, the task of AMS is really about finding the right exit 12), (13) for details, both examples use a parameter ξ = 0.25).

B x A Φ 3 Φ 2 Φ 1 Observable : Φ : R d → R Φ 1 < Φ 2 < Φ 3 N = 3
point. More details on the necessity of anticipation are given in appendix A.

plane Couette flow

Set up

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

The forced incompressible Navier-Stokes equations for the field u, the departure to the laminar base flow ye x and the pressure q, read

∂u i ∂t + u j ∂u i ∂x j + y ∂u i ∂x + δ i,x u y = - ∂q ∂x i + 1 R ∂ 2 u i ∂x 2 + ∂ 2 u i ∂y 2 + ∂ 2 u i ∂z 2 + f i (x, t) , ∂ i u i = 0 , (2)
using tensorial notations. We include the term f . It is a very general forcing which can be switched on or off. When it is on, it is white in time and in y. It can be red or white in x -z: in AMS simulations, we will use a temporally localised fully white perturbation, while in perturbed DNS, we use a red forcing (see § 2.2.2, 2.3 for details). These equations are discretised in space on N x and N z de-aliased Fourier modes (so that 3 2 N x and 3 2 N z modes are used in total) in the streamwise e x and e z directions and N y Chebyshev modes in the e y . Time integration is performed using channelflow, by J. Gibson [START_REF] Gibson | Visualizing the geometry of state space in plane couette flow[END_REF]. We investigate in detail systems of two sizes. The smaller system has size L x × L z = 24 × 18 (see Fig. 1 for scale), we set N y = 27, 3 2 N x = 128 and 3 2 N z = 96. The larger system has size L x × L z = 36 × 27, we set N y = 27, 3 2 N x = 196 and 3 2 N z = 144. In both cases we will set constant dt during the time integrations, dt = 0.05 (at R = 370 and R = 377) and dt = 0.02 (at R = 600). Both these values ensure stability of the time integration and respect the CFL criterion. We do not use adaptive time step in order to have more control on trajectory reconstruction and effect of the time step on convergence.

We define the spatially averaged kinetic energy as

E c = 1 2L x L z Lx x=0 1 y=-1 Lz z=0 u 2 x + u 2 y + u 2 z 2 dxdydz . (3)
We will also distinguish the kinetic energy contained in the streamwise component, on the one hand, and the kinetic energy contained in the spanwise and wall normal components, on the other hand

E c,x = 1 2L x L z Lx x=0 1 y=-1 Lz z=0 u 2 x 2 dxdydz , E c,y-z = 1 2L x L z Lx x=0 1 y=-1 Lz z=0 u 2 y + u 2 z 2 dxdydz . (4)
The first kinetic energy E c,x roughly quantifies the energy contained in velocity streaks, while the second E c,y-z roughly quantifies the energy contained in streamwise vortices [Jiménez andMoin, 1991, Hamilton et al., 1995]. These are the two main quantities of transitional wall flow turbulence.

Initial condition generation

In this section we present the procedure used for the generation of turbulent initial conditions used to study the collapse at Reynolds number R (by mean of DNS or AMS).

a) We first create an artificial velocity field

u y = 0, u z = 0, q = 0, u x = 0.4 sin π y+1 2 cos 4π z Lz M z , with M z = max 1, Lz

5

. This corresponds to streamwise velocity tubes which are prone to streaks instability [Waleffe, 1997] and should thus lead to wall turbulence. b) On top of the streamwise velocity components we add a noise red in x and z and white in y, this yields u 0 . This red noise is such that the variance of its n x , n z Fourier mode on component i is σ i γ nx γ nz , with γ x,z = 1 for 0 ≤ |n x,z | ≤ 6, γ x,z = 6 |nx,z| for n x,z > 6. We set γ x,z = 0 if n x,z > 28 in the L x × L z = 24 × 18 system and n x,z > 36 in the L x × L z = 36 × 27 system. We use σ x = 0.05, σ y = 0.0025 and σ z = 0.015. c) This initial condition is evolved for T 0 = 500 in the L x × L z = 24 × 18 system and T 0 = 200 in the L x × L z = 36 × 27 system at R + = 600. It has been checked that this duration was long enough so that natural buffer layer turbulence forms. If the kinetic energy of this velocity field is larger than 0.03, this yields u R+,1,0 , q R+,1,0 , otherwise we go back to step b) and generate a new u 0 with a new realisation of the red noise.

d) We then generate the N i initial conditions {u n , q n } 0≤n<Ni in the following manner . We first evolve

u R+,1,n , q R+,1,n at R + for T + = 500 (L x × L z = 24 × 18) or T + = 200 (L x × L z = 36 × 27), yielding {u R+,2,n , q R+,2,n } = u R+,1,n+1 , q R+,1,n+1 .
e) We then set the Reynolds at R (where we study collapse) and we set {u

R,1,n , q R,1,n } = {u R+,2,n , q R+,2,n }. This velocity field is evolved at R during T -= 750 (L x ×L z = 24×18) or T -= 500 (L x ×L z = 36×27).
This duration is chosen so that enough mixing has occurred and each initial condition is decorrelated from the others. We obtain {u R,2,n , q R,2,n }. We then let it evolve until the kinetic energy is either within 1.25% of E 0 , in which case we have our n th initial condition u n or is below 0.03. In that case we restart at d) by setting {u R+,1,n , q R+,1,n } = {u R+,2,n , q R+,2,n }. We use E 0 = 0.055 in the system of size L x × L z = 24 × 18 and E 0 = 0.052 in the system of size

L x × L z = 36 × 27.
This approach ensures that we have N decorrelated initial conditions which verify a given constraint on kinetic energy (for instance). It is of course easily parallelised with a minimal load imbalance. If we use c threads, we generate

N i = N c or N i = N c
on each thread, each thread uses an independent u R+,1,0 .

L x × L z = 24 × 18, L x × L z = 36 × 27.

Reaction coordinates

Since the kinetic energy E c (Eq. 4) of the turbulent flow is fluctuating around a conditional average, while the kinetic energy of the laminar flow is zero (Fig. 3 (a)), a first choice to construct the reaction coordinate is to use E c (t). We can therefore propose the reaction coordinate φ E defined as

Φ E (t) = E t -E c (t) ∆E . (5)
It is natural to choose E i close to some average of the kinetic energy and ∆E E t . In order to estimate E t , one can sample E c (t) and construct the empirical probability density function conditioned on the flow experiencing no collapse of turbulence, in a system of given size and Reynolds number. One can for instance retain the part of the times series where E c (t) ≥ 0.025 (indicated by the black dashed line in figure 3 (a)). Using this, one can compute a sample mean of the kinetic energy E (Fig. 3 (b)), and then choose close enough E t and ∆E accordingly. We will state what values of E t and ∆E are chosen in AMS computations in section 3 and section 4. We also compute the conditional variance as a function of the Reynolds number for the two system sizes (Fig. 3 (c)): we will use this quantity to choose the hypersurface C, so that we typically have E t -E c (t) ≃ σ 2 on C. This corresponds to a short excursion when the kinetic energy is away from its conditional average by typically half a standard deviation. Using the reaction coordinate, we set A as all velocity fields such that Φ E ≤ 0 and B as all velocity fields such that Φ E ≥ 1.

In appendix A, we use an other reaction coordinate Φ a , based on ψ(t), the asymmetry of the streamwise velocity field with respect to the midplane y = 0

ψ(t) = 1 V Lx x=0 1 y=-1 Lz z=0 u x sgn(y) dxdydz , (6)
where sgn(y) stands for the sign of y. The choice of this reaction coordinate is motivated by the fact that we have mostly u x < 0 (coming from low speed streaks) for y > 0 and mostly u x > 0 (coming from high speed streaks) for y < 0 when velocity streaks (and thus wall turbulence) are present in the flow [START_REF] Jiménez | The minimal flow unit in near-wall turbulence[END_REF], Hamilton et al., 1995, Kawahara et al., 2003]. Note that with this definition ψ is mostly negative. We proceeded in the same manner as with the kinetic energy in order to construct a reaction coordinate: conditional averages of ψ are performed in order to calculate Φ a (t) = 1 -ψ(t) ψ .

Direct numerical simulations of collapses

We perform two types of direct numerical simulations of turbulent collapse. The first kind of simulation consists in letting the flow freely evolve from initial conditions computed following the method of section 2.2.2. These simulations are stopped when the flow has laminarised. This is deemed to have happened when the reaction coordinate has reached one. In these simulations, we only save time series, such as those of E c , E c,x , E c,yz and Φ E (an example of such time series in given in figure 3 (a)). From the last part of each DNS, we can extract the time series of these quantities in natural collapse trajectories. We will perform two types of such DNS: some free natural DNS, and some forced DNS, in order to check the effect of additional noise on the collapse of turbulence. An additional noise will necessarily be included in AMS computations, and we wish to know what is the minimal error caused on the probability of collapse, trajectory features and durations and mean first passage times by the addition of this noise.

Since we cannot add a comparable white tick in DNS, we will add a red permanent noise. In that case the forcing f (Eq. ( 2)) is exerted at all time and it is red in x and z. The forcing noise is characterised by its spatial correlation function, which is prescribed through its spectrum

f i (x, t)f j (x ′ , t ′ ) = δ ij δ(t -t ′ )C i (x -x ′ , z -z ′ ) Ci (y -y ′ ) , Ĉi (n x , n z ) = Γ i,nx,nz . (7)
In Γ i,nx,nz , i stands for the component e x , e y and e z (i = 1, 2, 3), n x stands for the Streamwise wavenumber, n z stands for the spanwise wave number. In this text we always use Ci (y-y ′ ) = δ(y-y ′ ) in numerical simulations using a finite number N y of Chebyshev modes in the wall normal direction. In DNS, we will always use Γ i,nx,nz = γ nx γ nz with the previously defined γ (see § 2.2.2). We set σ = 10 -9 . Another kind of direct numerical simulation consists in repeating the stage 0 of AMS, where we start the simulation from our initial conditions, let them cross the hypersurface C and then either laminarise or go back to φ < 0. From these, we can compute the proportion of trajectory that laminarises and thus have an unbiased estimate of the probability of collapse and trajectory durations and validate AMS estimates of these quantities.

3 System of size L x × L z = 24 × 18: reactive trajectories and validation of AMS computations

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

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

Collapse trajectories

Visualisation of the collapse trajectories

We first describe the velocity fields during the collapse, in trajectories computed by AMS. These are available in three dimensions for the three components in much greater quantity than for direct numerical simulations. We only need to record the collapse trajectories at time intervals (sampled at time intervals δt > 1) from the final set of reactive trajectories as well as non reactive trajectories obtained at the last step of AMS computations. We do not need to discard the very long fluctuations around set A contained in DNS of natural collapse, thus drastically reducing the disk space requirements for systematic studies of the collapses.

We display colour levels of the streamwise and spanwise components of the velocity field in a midplane at successive instants during a collapse trajectory in figure 4 and during one of the non reactive trajectory remaining at the last stage of AMS in figure 5. We will focus on the two fields (u x , u z ) because of the physics of wall flows and experimental constraints. On the one hand, it has been shown in models and observed in some experiments and simulations that there was a distinct behaviour of velocity components. One finds velocity streaks contained mostly in u x as well as streamwise vortices ω x that can be observed through the component u z (they are the main contribution to this velocity component) [START_REF] Hamilton | Regeneration mechanisms of near-wall turbulence structures[END_REF]. Moreover it has been observed in models, simulations and experiments that decay of turbulence was much faster in streamwise vortices (or its proxies) than in velocity streaks [Rolland, 2018, Liu et al., 2020, Gomé et al., 2020]. We will endeavour to discuss this observation in specifically computed collapse trajectories. On the other hand, most experimental observations are performed by PIV using a plane parallel to the wall, thus capturing u x and u z in a x -z plane near the midgap [START_REF] Liu | Anisotropic decay of turbulence in plane couette-poiseuille flow[END_REF], De Souza et al., 2020].

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

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

We can complement this view by observing the streamwise vorticity field ω x,last = ∂ y u z,last -∂ z u y,last leading to the larger reaction coordinate in the last suppressed state in the last stage of AMS (Fig. 6). In this system of relatively moderate size, even if one always has an impression that turbulence collapses through the formation of a laminar hole, the structure of this turning point state varies from run to run. One can observe clearly localised streamwise vortices (Fig. 6 fairly different ω x,last is that they correspond to different instants on similar trajectories that do display a hole opening. This is also caused by the fact that reactive trajectories are not yet that concentrated around a typical trajectory.

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

Comparison between trajectories computed by AMS and by DNS

We then start our comparison of the collapse trajectories computed by AMS and DNS by describing the duration of collapse trajectories, that is to say the time s elapsed between the instant where the flow crossed the hypersurface C and the instant where the reaction coordinate reaches Φ = 1. We have obtained 1138 trajectories from AMS and 189 trajectories form DNS. We firstly compute the sample mean τ and the sample variance σ of the trajectory durations, we obtain τ AMS = 413 ± 3, τ DN S = 410 ± 7 and σ AMS = 105 and σ DN S = 91. The error bars on τ are given by σ divided by the square root of the number of samples2 . We can see that both are within a few percent of one another and within the 66% confidence interval of each other. This shows that the reactive trajectories computed by AMS and DNS have almost the same duration, and thus that AMS has very little bias in the way it computes the collapse trajectories. We then compare the sampled distributions of trajectory durations. In figure 7, we display the distribution of normalised durations s = s-τ σ of collapse trajectories computed by mean of AMS and of DNS. Both are very similar and seem to be very close to a Gumbel distribution. This distribution of duration has been reported to be very close to a Gumbel in many instances of stochastically driven systems [Rolland, 2018, Rolland et al., 2016, Rolland and Simonnet, 2015]. This was often related to a simple structure of the reactive trajectory that went through a saddle of the deterministic part of the system. It is actually demonstrated that reactive trajectories durations have this distribution in one dimensional stochastic systems [START_REF] Cérou | On the length of one-dimensional reactive paths[END_REF]. To our knowledge, this is possibly the first time that such a distribution shape has been shown in turbulence without stochastic injection of energy. This suggests us that the manner in which turbulence collapses in systems of large enough size is most certainly not unrelated to reactive trajectories in stochastic systems, in that it quite possibly follows some effective dynamics agitated by background turbulence. While it is more difficult to separate the dynamics into a deterministic part and a stochastic part, this indicates that the reactive trajectories can follow a rather simple path.

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

We can compare the trajectories more precisely by computing the conditional average log(E c,y-z ) as a function of log(E c,x ). Indeed, if we define the probability of having the value log(E c,x ) during the

We can also compute the variance

σ y-z (log(E c,x )) = ρ c (log(e c,y-z ))| log(E c,x ) log(e c,y-z ) 2 d log(e c,y-z ) -log(E c,y-z ) 2 (log(E c,x )) .
(11) We display the average log(E c,y-z ) (log(E c,x )) in full line and the average plus/minus the variance log(E c,y-z ) (log(E c,x )) ± σ y-z (log(E c,x )) in dashed line for all three types of reactive trajectories in figure 8 (d). This quantitatively confirms the observation of a faster decay of the streamwise vortices and thus of E c,y-z than of the streaks, and thus of E c,x . This also confirms that up to the effect of the noise when E c,y-z is very small, using the AMS leads to the same collapse trajectories as using DNS, at a much smaller cost. We can finally note that the trajectories seem to be concentrated around this average. Note however that the constant variance of log(E c,y-z ) actually corresponds to a variance of E c,y-z which decreases exponentially as the collapse of turbulence goes on. There is a larger variance of E c,y-z among the initial conditions and thus in the amplitude and distribution of streamwise vortices in the initial conditions and stronger fluctuations in the beginning of the trajectory. On a final note, let us remind that in order to rigourously obtain instantons in stochastic systems, one has to take a limit of the variance of the added noise going to 0. In our turbulent case, this means either identifying an effective noise (which is a difficult task) or displaying large deviations in quantities like the probability density functions, the probability of collapse or the mean first passage time. Quantitatively speaking, this amounts to identifying a small parameter ǫ, proportional to the square of the suspected noise variance such that lim ǫ→0 ǫ log(α) or lim ǫ→0 ǫ log(T ) are independent on ǫ (these two limits are often the opposite of one another) [Touchette, 2009]. As it happens, such scalings have been displayed in probability density functions in simulations [Rolland, 2015], where ǫ DN S = 1

LxLz and in mean first passage times before collapse in models, where ǫ = 1 L [Rolland, 2018], when the sizes L x , L z and L of the system was increased. Thus, in the limit of large system length, one should expect to observe narrower and narrower concentration of trajectories. In order to observe more concentration of trajectories, and more generally large deviations of the probabilities and rate of probability of collapse, it may very well be necessary to observe said collapses in larger and larger domains. For this purpose increasing the streamwise size L x may be more efficient than increasing the spanwise size L z . This may be a consequence of the specific mechanism of hole formation, leading to the formation of two laminar turbulent fronts that moves on. Such situations have already been observed in theoretical physics systems [START_REF] Rolland | Computing transition rates for the 1-d stochastic ginzburg-landau-allen-cahn equation for finite-amplitude noise with a rare event algorithm[END_REF]. The rare event is the formation of the hole and of the front themselves. Once said hole is created, the probability that it will entirely open is small and decreases with size, but this decrease is slower than exponential.

We finally note that even if each AMS run has some collapse trajectories in its first iterations (approximately 3.8% of the trajectories generated in the first stage), and that it will generate other collapse trajectories by branching on these collapse trajectories, it generates even more variability than this through mutation/selection on other initial conditions. This means that each run generates distinct, independent reactive trajectories, in particular trajectories with distinct initial conditions. These trajectories are even more distinct in their later stages. The amount of variability generated can be estimated by counting the number of genealogies, that is to say the number of group of trajectories that share the same starting point. Each group will display trajectories decorrelated from those of another genealogy.

On average, we generate 24 genealogies per AMS run with 120 clones. This value is obtained by averaging over all available AMS runs the number of distinct initial conditions in computed reactive trajectories.

See [START_REF] Ragone | Rare event algorithm study of extreme warm summers and heat waves over europe[END_REF] for a view of genealogies of trajectories with a different algorithm applied to climate.

Mean first passage time and probability of collapse

We now compare the mean first passage time T and probability α of collapse computed by mean of direct numerical simulations and AMS. While the mean first passage time before collapse is the most physical quantity, the probability of collapse is always worth discussing. Looking into its statistics is always a good way of validating the AMS computations. This is all the more true because this is the one quantity for which we have mathematical results of convergence [Cérou andGuyader, 2007, Bréhier et al., 2016]. Using the methods presented in sections 2.1 and 2.3, we obtain α o = α AMS = 0.042 ± 0.003 with a variance σ α,AMS ≃ 0.016 using AMS, and α DN S = 0.037 ± 0.002 using direct numerical simulations. This indicates that the estimation of the probability of collapse is precise, we have an overlap of both 66% interval of confidence. When we suppress a fixed number (strictly larger than one) of clones at each iteration, the ideal variance of the probability of collapse is given by σ

2 id = α 2 N k Nc N -Nc + 1-ρ ρ (
where k is the average number of iterations of AMS and ρ = r is the average proportion of reactive trajectories computed at the last iteration) [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF]. This is the minimum of variance of the estimator of α. It is obtained when the ideal reaction coordinate, the committor is used [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF]. When this minimal variance is reached, AMS computations are performed with a minimal error. In practice this gives a value of the ideal variance σ AMS,id ≃ 0.007. The variance on the estimate of the probability is twice that of the ideal variance. This indicates that the quality of the computation is acceptable, but not perfect, quite possibly due to imperfections in the reaction coordinate we used.

We then compare the mean first passage times before collapse. By means of DNS, we sample the cumulated density of waiting times before collapse in simulations performed with and without additional noise, in order to estimate the mean first passage time as well as test the effect of noise on it. We can see that cumulated densities F (t) = ∞ t f (ζ) dζ, with f (ζ) the PDF of passage times, computed with and without noise decrease linearly in logarithmic scale. This is thus entirely consistent with an exponential distribution of passage times. Both cumulated density functions have very similar slopes, which is itself almost equal to the estimated mean first passage time (Fig. 9). We estimate the mean first passage time before collapse in DNS by averaging over all sampled durations, we obtain T DN S = 1.10 • 10 4 ± 8 • 10 2 without noise and T DN S,noisy = 1.3 • 10 4 ± 10 3 with additional noise. In both cases, the variance is equal to the average and the median is equal to √ 2T , further confirming that we have an exponential distribution. While the addition of noise may increase the mean first passage time, it does not increase it dramatically, so that it is acceptable to use an additive small noise in simulations and AMS to help the separation of trajectories. In turn, we estimate T AMS = 2.4•10 4 ±3•10 3 , which is clearly an overestimate. In order to understand why the estimate of the mean first passage time is not as good as what is obtained for other quantities, we examine the histograms of mean first passage times (Fig. 10 (a)) and probability of collapse (Fig. 10 (b)) computed by AMS. While the histogram of α is symmetric, that of T is skewed. We note that the estimate of T is polluted by a handful of very high values, originating from fairly small values of α, which have a much smaller effect on α AMS . On top of the improvement brought by a better reaction coordinate, the computation of the mean first passage time may very well be improved in AMS computations using more clones.

One can finally compare the relative cost of trajectory computations with AMS and DNS. Strictly speaking, the average physical duration of simulation to obtain one reactive trajectory by DNS is of order O(10 4 ), meanwhile the simulation time by means of AMS is of order O(10 3 ). Such a factor 10 of acceleration of computation already seems interesting. Moreover, one needs to remember that when we use DNS to compute reactive trajectories, the data of a large part of these O(10 4 ) time units have to be discarded (which means that more often than not, only time series can be obtained). Meanwhile a large part of the simulation time performed by AMS corresponds to directly usable data (for the purpose of studying reactive trajectories) and it makes sense to record not only the time series, but also velocity fields during the collapses. We therefore not only optimise the computation time, we also optimise the records of trajectories. This means that even in this case where the study of the system by DNS seems affordable, using AMS already leads to a substantial acceleration of computations.

4 A very rare collapse in a domain of size L x × L z = 36 × 27: spatial organisation

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

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

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

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

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

We can observe the movement of the laminar and turbulent fronts in more detail for the reactive trajectory (Fig. 14 x dx = 0.15 which respectively indicate the spanwise location of laminar-turbulent interface for the streamwise vortices and for the velocity streaks. These can be compared to earlier visualisations (Fig. 12 and Fig. 13), that justified the study of the purely spanwise location of the laminar-turbulent interface in this type of domain. In the case of the collapse trajectory, we can observe the survival of streaks where streamwise vortices have disappeared for 500 t 800 and 20 z 27 and for 900 t 1100 and 5 z 20 (Fig. 14 (a)). The concomitant advance of the laminar turbulent interfaces can be seen. Note that at z = 15, between t = 400 and t = 600, the receding laminar-turbulent interfaces for the streaks and the streamwise vortices have almost the same location. While this is partially caused by the method of detection of the interface (see Fig. 12 t = 624) and while the scenario is not at one hundred percent univocal, the situation where streamwise vortices survive where streaks have decayed is never observed in the reactive and non reactive trajectories. Similar observations can be made in the case of the non reactive trajectory (Fig. 14 (b)), while we have an apparent collapse of both streaks and vortices at t = 400 and for 15 z 20, for 500 t 800, only the laminar-turbulent interface of the vortices recedes. For 800 1200, we can see some reduction of the surface occupied by the streaks. For t 1200, we can see that the surface occupied by the vortices increases again until the laminar turbulent-interface for the streaks and the vortices is again the same and they both reinvade the domain.

Conclusion

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

A first step in this work has been to ensure the validity of the result of AMS computations. This has been done in a system of size L x × L z = 24 × 18 at Reynolds number R = 370. Computation of reactive trajectories and mean first passage time are amenable in this system, so that rare event computations could be compared to a reference. It could thus be shown that reactive trajectories (trajectories in phase space, statistics of durations) are computed precisely. Some effects of the low noise added on the flow for the separation of trajectories was identified in the very last stages of reactive trajectories, when the flow is finally laminar. The bias added to the collapse computations is small enough so that the use of this additional noise is entirely acceptable. Similarly, the probability of collapse computed using AMS was shown to be in the same 66% interval of confidence as the one computed by DNS. The computed mean first passage times are overestimated. Histograms of the computed probabilities of collapse and mean first passage time before collapse indicate that this overestimate is cause by a handful of overly large estimates of T related to underestimates of α. This is most likely an effect of the number of clones used in this computation. Other estimators of mean first passage time could be used to improve the computation [START_REF] Lestang | Numerical study of extreme mechanical force exerted by a turbulent flow on a bluff body by direct and rare-event sampling techniques[END_REF], which have their own drawbacks for trajectory computations.

A large part of this work was of course devoted to the study of the collapse trajectories properties. In the system of size L x × L z = 24 × 18, the visualisation of collapse trajectories showed that the laminarisation could be sketched in two stages. The velocity components u y and u z , and as a consequence, streamwise vortices, decayed first, so that the flow is driven out of the self sustaining process of turbulence [Waleffe, 1997]. As a consequence, streamwise velocity has less dependence on x, and then slowly decays without display streaks instability [START_REF] Marquillie | Instability of streaks in wall turbulence with adverse pressure gradient[END_REF], Kawahara et al., 2003]. Examination of the velocity fields shows that this process can happen locally, leading to a streamwise long laminar hole in u y , u z and then extends in z or more globally: streamwise vortices decay everywhere in this small domain. Viewed in the quadrant E c,x ≥ 0, E c,y-z ≥ 0, these visualisations all correspond to trajectories concentrated along the same path. While this path is more continuous than the sketch given here, it still displays a first stage where E c,y-z decays much faster than E c,x , then a second stage where E c,y-z is small (around the noise level) and E c,x decays in turn. Such concentration of trajectories are commonplace in stochastic models which display rare transitions when the noise variance goes to zero. While the separation of the dynamics of an actually turbulent flow into a deterministic part and a noise part is very difficult to make, and that the flow displays transient chaos when the domain is very small [START_REF] Eckhardt | Turbulence transition in pipe flow[END_REF], this is a first indication that the collapse of turbulence in Couette or Poiseuille flows of increasing size is best understood as a stochastic system. Given the concentration of trajectory, one may be tempted to draw a comparison to Instantons computed in stochastic systems [Bouchet et al., 2019, Grafke andVanden-Eijnden, 2019]. In order to give a meaning to such a description, one would need a low noise limit. This can possibly be the case in the large size limits as observed in models [Rolland, 2018] and in the large deviations of the kinetic energy [Rolland, 2015]. The statistics of reactive trajectories durations are another feature of the trajectories which are very close to their stochastic counterparts. Both DNS and AMS computations show that the normalised distribution of durations is very close to a Gumbel distribution. In stochastic multistable systems it is demonstrated in one dimensional systems that this is the distribution of durations [START_REF] Cérou | On the length of one-dimensional reactive paths[END_REF]. In highly dimensional systems, it is observed that this distribution is found when the multistability occurs when the system transits through a single curved saddle of its deterministic part [START_REF] Rolland | Computing transition rates for the 1-d stochastic ginzburg-landau-allen-cahn equation for finite-amplitude noise with a rare event algorithm[END_REF], Rolland, 2018]. This would place the collapse of turbulence in plane Couette flow in systems of increasing sizes in this category.

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

Another interesting point of the use of a method like AMS is the ability to identify a turning point in the dynamics of collapse [START_REF] Simonnet | Multistability and rare spontaneous transitions between climate and jet configurations in a barotropic model of the jovian midlatitude troposphere[END_REF]. This is done by recording the state leading to the largest reaction coordinate in the suppressed trajectories of the last stage of AMS. In stochastically forced systems, this state corresponds to an effective saddle between the two metastable states, and can alternatively be computed through a dichotomy procedure (see [Willis andKerswell, 2009, Schneider et al., 2007]), or even by computing the saddle point of the deterministic point of the system in systems forced by an additive noise. This turning point gives us some information on the mechanism of multistability and in our case on collapse of turbulence. They displayed localised streaks and streamwise vortices in the largest system, which backs the scenario of collapse through hole formation. In the case of collapse of turbulence in plane Couette flow, these states did not necessarily make good starting points for edge tracking. Note that edge tracking can be performed using the result of AMS computation: this has been done in the study of build up of turbulence under some additive noise and should be presented in a later text.

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

Finally, the computation of reactive trajectories in transitional wall turbulence opens the way for similar computations in high Reynolds number turbulence. Systematic computation of reactive trajectories in systems of aerodynamical interest [START_REF] Kim | Investigation of the flow between a pair of circular cylinders in the flopping regime[END_REF], or geophysical interest are more than conceivable [START_REF] Herbert | Atmospheric bistability and abrupt transitions to superrotation: wave-jet resonance and hadley cell feedbacks[END_REF]. and a handful of streamwise vortices. Their behaviour can be described by the escape away from a chaotic repeller [START_REF] Eckhardt | Dynamical systems and the transition to turbulence in linearly stable shear flows[END_REF]. They necessarily collapse globally in space, this means that the larger Φ is, the more globally quiescent the flow becomes: the mixing is less intense and there is little separation of trajectories. This is shown in figure 16 in a small system of size L x × L z = 12 × 8. The flow required some reorganisation before collapse was possible (part of the trajectory for 0 ≤ t ≤ 300), and at most two low speed and two high speed streaks are present. When the collapse takes place (between t = 350 and t = 400) u z collapses almost uniformly. Larger systems have more and more degrees of freedom, and contain more velocity streaks and vortices. An important feature is that they collapse through the formation of a laminar hole. This means that for a large portion of the reactive trajectory some area of the flow remains turbulent. This turbulent area can serve as some sort of heat bath which fuels the trajectory with perturbations and thus can greatly help the separation of trajectories. This means that in small systems, branching must indeed be performed on trajectories having greater excursions, but on their initial stages, where the flow is still agitated and thus will separate the trajectories. In that case, we use a so called saturated anticipation

Φ b,sat = ξ 1 -exp - Φ Nc ξ ( 12
)
this function is shown with ξ = 0.25 in figure 2 (c), with the "saturated" label. If we follow that strategy in the case of larger systems, it will require too many tries for branched trajectories to reach Φ Nc . A more efficient strategy is to anticipate the branching with a small shift, this can for instance be done with a function that slowly reaches its asymptote, the line of slope 1,

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

Figure 1 :

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

Figure 3 :

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

  (a)), an almost entirely quiescent flow (Fig. 6 (b)) or rather active streamwise vortices all over the domain (Fig. 6 (c)). One of the reasons we have

Figure 4 :

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

Figure 5 :

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

Figure 6 :Figure 7 :

 67 Figure 6: Streamwise vorticity in the midplane y = 0 in the system of size L x ×L z = 24×18 for last states at the last branching stage of several AMS computations. (a) Run estimating α = 0.045, T = 1.5 • 10 4 , (b) run estimating α = 0.086, T = 1.1 • 10 4 , (c) run estimating α = 0.022, T = 3.7 • 10 4 .

  reactive trajectories ρ x (log(E c,x )) = rho e (log(E c,x ), log(e c,y-z )) d log(e c,y-z ) , (8) we have the conditional probability ρ c log(e c,y-z ) log(E c,x ) = ρ e (log(e c,x ), log(E c,y-z )) ρ x (log(E c,x )) . c,y-z ) (log(E c,x )) = ρ c log(e c,y-z ) log(E c,x ) log(e c,y-z ) d log(e c,y-z ) .

Figure 8 :

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

Figure 9 :Figure 10 :

 910 Figure9: Logarithm of the cumulated distribution of waiting times in noiseless and noisy DNS for L x × L z = 24 × 18, R = 370. We add the two linear functions -t/ T for comparison.

Figure 11 :

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

  (a)) and the non reactive trajectory (Fig.14 (b)). These figures use both the

Figure 12 :Figure 13 :

 1213 Figure 12: Following a collapse trajectory in a domain of size L x × L z = 36 × 27, R = 377. The left panels show a time series of kinetic energy, with the dot indicating the instant at which the colour levels of streamwise velocity in the y = 0 plane (middle panels) and spanwise velocity (right panels) at corresponding times in the y = 0 plane are shown. 22

Figure 16 :

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

  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. (b) Sketch of the principle of AMS, showing two 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 t Φ. Trajectory one (dashed blue line) is suppressed and branched on another trajectory at level max t Φ 1 and then ran according to its natural dynamics. Trajectory 2 is then suppressed and branched on 3 at level max t Φ 2 . (c) Two examples of anticipation of branching level of reaction coordinate as a function of maximum reaction coordinate reached by suppressed trajectories Φ Nc (Φ b ) tested in anticipated branching (see § A, Eq. (

	(c)	1
				converging
		0.8	saturating
	clones	0.6
	b
		0.4
		0.2
			0
			0	0.2	0.4	0.6	0.8	1
					N	c
	Figure 2: (a)

In that case, when extinction is detected, the computation is terminated and the way AMS is used is reassessed.

It is safe enough to assume each value is decorrelated enough from the large majority of the others to assert that the central limit theorem can be applied to the sampled mean: the most general formulation only require short length correlations. There may be some degree of correlations between durations of reactive trajectories that arise from the same genealogy, but each genealogy is small enough compared to the sample size.

Acknowledgements

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

A Anticipated AMS

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

The anticipation strategy which should be followed depends on the properties of the system. In the case of the collapse of turbulence in plane Couette flow, this depends on the size of the system. The line is drawn between small systems (typically of size L x 15, and L z 10) and larger systems. Small systems which have very few degrees of freedom, display at most four or five velocity streaks (high or low velocity)