2.2.1. Principle of the algorithm

Before presenting the principle of anticipated AMS, let us first give a formal phase space description of the rare events we will study in this text. Let us sketch the collapse of turbulence in figure 2(a) and consider the set $\mathcal {A}$, a neighbourhood of the turbulent flow in phase space, and the set $\mathcal {B}$, a neighbourhood of the laminar flow. A realisation of the dynamics which starts in $\mathcal {A}$ fluctuates around it, has several excursions out of $\mathcal {C}$, a hypersurface closely surrounding $\mathcal {A}$, and eventually crosses $\mathcal {C}$ and reaches $\mathcal {B}$ before coming back to $\mathcal {A}$, is termed a first passage. Its average duration is termed the mean first passage time $T$. The last stage of the dynamics is termed a reactive trajectory: this is the part of the dynamics that starts in $\mathcal {A}$, crosses $\mathcal {C}$ and reaches $\mathcal {B}$ before $\mathcal {A}$. Precise definitions of sets $\mathcal {A}$, $\mathcal {B}$ and hypersurface $\mathcal {C}$ for collapse will be given in §§ 3 and 4, based on reaction coordinates defined in § 2.4.

Figure 2. (a) Sketch of two bistable states $\mathcal {A}$ and $\mathcal {B}$ and the hypersurface $\mathcal {C}$ closely surrounding $\mathcal {A}$. Two realisations of the dynamics are sketched: a single excursion in blue and a first passage trajectory in black and red. The red part of the first passage trajectory is the reactive trajectory (figure originally made for Rolland & Simonnet Reference Rolland and Simonnet2015). (b) Sketch of the principle of AMS, showing two iterations of the algorithm, with $N=3$ clones, indicating the starting state $\mathcal {A}$ and its neighbourhood, the arrival state and its neighbourhood $\mathcal {B}$, three trajectories are ordered by their $\max _t\varPhi$. Trajectory one (dashed blue line) is suppressed and branched on another trajectory at level $\max _t\varPhi _1$ and then ran according to its natural dynamics. Trajectory $2$ is then suppressed and branched on 3 at level $\max _t\varPhi _2$ (figure originally made for Simonnet Reference Simonnet2016). (c) Two examples of anticipation of branching level of reaction coordinate as a function of maximum reaction coordinate reached by suppressed trajectories $\varPhi _{N_c}(\varPhi _b)$ tested in anticipated branching (see Appendix A, (A1), (A2) for details, both examples use a parameter $\xi =0.25$).

We then give a brief overview of the variant of AMS, termed anticipated AMS, which was used for the computation presented in this text (see Cérou & Guyader (Reference Cérou and Guyader2007), Rolland & Simonnet (Reference Rolland and Simonnet2015), Bréhier et al. (Reference Bréhier, Gazeau, Goudenège, Lelièvre and Rousset2016), Rolland (Reference Rolland2018), Cérou et al. (Reference Cérou, Guyader and Rousset2019b) and Lestang et al. (Reference Lestang, Bouchet and Lévêque2020) for more details on the general methods). All variants of AMS use a reaction coordinate (or observable) $\phi (\boldsymbol {u})$, a real-valued function of the velocity field. The reaction coordinate gives a relative distance in phase space between $\boldsymbol {u}$ and the sets $\mathcal {A}$ and $\mathcal {B}$. For practical reasons, we will see the reaction coordinate as a function of time $\varPhi (t)=\phi (\boldsymbol {u}(t))$ on the trajectories. The reaction coordinate is often rescaled such that $\phi (\partial \mathcal {A})=0$, on the boundary of set $\mathcal {A}$, $\phi (\partial \mathcal {B})=1$, on the boundary of set $\mathcal {B}$ and grows monotonically in between. All variants run $N$ clone dynamics of the system to compute iteratively at least $N-N_c>0$ reactive trajectories going from a hypersurface $\mathcal {C}$, close to set $\mathcal {A}$, towards the set $\mathcal {B}$. At each iteration, $N_c< N$ clones are replaced. The algorithm is sketched in figure 2(b) and proceeds in the following manner:

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

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

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

The algorithm stops when $Nr\ge N-N_c+1$ clone trajectories have reached $\mathcal {B}$. We usually perform $o>1$ independent AMS runs. In each of the runs labelled by the number $a$, $1\le a\le o$, we obtain the total number of iterations $\kappa _a$, with $r_a$ the proportion of clones that transit to $\mathcal {B}$. This yields an estimator of the probability $\alpha$ of reaching $\mathcal {B}$ before $\mathcal {A}$ (Cérou & Guyader Reference Cérou and Guyader2007), and the corresponding mean first passage time $T$ (Cérou et al. Reference Cérou, Guyader, Lelievre and Pommier2011)

(2.4a–c)\begin{align} \alpha_a=r_a\left(1-\frac{N_c}{N} \right)^{\kappa_a},\quad \hat{\alpha}=\left\langle \alpha_a\right\rangle_o,\quad \hat{T}=\left\langle\left(\frac{1}{\alpha_a}-1\right)(t_{1,a}+\tilde{\tau}_a)+ (t_{1,a}+\tau_a)\right\rangle_o , \end{align}

where $\tau _a$ is the mean duration of reactive trajectories, $t_{1,a}$ is the mean duration to go from $\mathcal {A}$ to $\mathcal {C}$ and $\tilde {\tau }_a$ is the mean duration of non-reactive trajectories, computed in each AMS run $1\le a \le o$. The notation $\langle \cdot \rangle _o=({1}/{o})\sum _{a=1}^o\cdot _a$ corresponds to an average over the $o$ independent AMS runs (Cérou & Guyader Reference Cérou and Guyader2007; Rolland Reference Rolland2018).

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

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

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

2.2.2. Anticipated AMS and branching level $\varPhi _b$

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

2.2.3. On the goodness of AMS computations

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

1. (a) a large number of reactive trajectories, in our case from turbulent to laminar flow;

2. (b) provide an estimate of the probability that these reactive trajectories occur;

3. (c) estimate the mean first passage time before a reactive trajectory occurs.

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

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

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

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

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

3.1.1. Visualisation of the collapse trajectories

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

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

1. (i) Indeed, on the one hand, it has been shown in models and observed in some experiments and simulations that there was a distinct behaviour of velocity components. One finds velocity streaks contained mostly in $u_x$ as well as streamwise vortices, visible in the streamwise vorticity field $\omega _x=\partial _yu_{z}-\partial _zu_{y}$, but that can be observed through the component $u_z$ (they are the main contribution to this velocity component) (Hamilton et al. Reference Hamilton, Kim and Waleffe1995). Moreover, it has been very recently observed in models, simulations and experiments that decay of turbulence is much faster in streamwise vortices (or its proxies) than in velocity streaks (Rolland Reference Rolland2018; Gomé et al. Reference Gomé, Tuckerman and Barkley2020; Liu et al. Reference Liu, Semin, Klotz, Godoy-Diana, Wesfreid and Mullin2021). We will endeavour to discuss this observation in specifically computed collapse trajectories.

2. (ii) On the other hand, most experimental observations are performed by particle image velocimetry (PIV) using a plane parallel to the wall, thus capturing $u_x$ and $u_z$ in an ($x$,$z$) plane near the midgap (De Souza, Bergier & Monchaux Reference De Souza, Bergier and Monchaux2020; Liu et al. Reference Liu, Semin, Klotz, Godoy-Diana, Wesfreid and Mullin2021).

Figure 4. Time series of the kinetic energy (a,d,g,j) 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 (b,e,h,k) and spanwise velocity (c,f,i,l), at corresponding times, in the midplane $y=0$ during a collapse trajectory in a domain of size $L_x\times L_z=24\times 18$ computed by AMS.

Figure 5. Time series of the kinetic energy (a,d,g,j,m) 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 (b,e,h,k,n) and spanwise velocity (c,f,i,l,o) in the midplane $y=0$ during a non-reactive trajectory (hole opening then closing), at corresponding times, in a domain of size $L_x\times L_z=24\times 18$ computed by AMS.

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

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 (figure 5). As with all other trajectories, this one starts with turbulence in the whole domain (figure 5, $t=2$). We again observe that $u_z$ decays faster than $u_x$ ($t=74$). In that case, the spanwise localisation of the decay is not disputable. Along with the decay of $u_z$, the streamwise velocity component becomes streamwise invariant ($t=200$). The kinetic energy fluctuates for some time around a plateau ($150\lesssim t\lesssim 350$), then increases again. The values taken by $u_z$ are getting more intense over an area of the domain that increases ($t=400$). Note that the streamwise velocity field had decayed locally during the plateau (for $0\lesssim z\lesssim 7$). We observe an asymmetry between the decay and the reinvasion processes. The component $u_z$ retracts to smaller areas than $u_x$ during the decay: the active–quiescent interface (typically planes $z=\text {cst}$) is not the same for $u_x$ and $u_z$. A strict definition of these interfaces will be proposed in § 4 and figure 14. Then, both $u_z$ and a streamwise dependent $u_x$ fully reinvade the domain ($t=600$), in that case the active–quiescent interface is the same for both components. We thus observe a concomitant restart of both components of the self-sustaining process of turbulence.

We can complement this view by observing the streamwise vorticity field $\omega _{x,last}=\partial _yu_{z,last}-\partial _zu_{y,last}$ leading to the largest value of the reaction coordinate in the last suppressed state in the last stage of AMS (figure 6). In this system of relatively moderate size, even if one always has an impression that turbulence collapses through the formation of a laminar hole, the structure of this turning point state varies from one run to another. One can observe clearly localised streamwise vortices (figure 6a), an almost entirely quiescent flow (figure 6b) or rather active streamwise vortices all over the domain (figure 6c). One of the reasons we have fairly different $\omega _{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 the size of the domain is not large enough for the reactive trajectories to be concentrated around a typical trajectory.

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

In any case, we can have a first view of the path followed during the collapse of turbulence. In terms of the manner in which the cyclic self-sustaining process of wall turbulence fails (streaks first, vortices first or concomitant decay of streaks and vortices), the observed scenario is that of the first decay of streamwise vortices followed by decay of the velocity streaks. No more energy is extracted from the base flow so that the velocity streaks, turned into tubes, slowly decay. They do not undergo any new streak instability that would refuel the streamwise vortices (Jiménez & Moin Reference Jiménez and Moin1991; Marquillie, Ehrenstein & Laval Reference Marquillie, Ehrenstein and Laval2011). We will examine all the reactive trajectories 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 § 4 to check whether this observation is disputable or not.

3.1.2. Comparison between trajectories computed by AMS and by DNS

We then start our comparison of the collapse trajectories computed by AMS and DNS by describing the duration of collapse trajectories, that is to say the time $s$ elapsed between the instant where the flow crossed the hypersurface $\mathcal {C}$ and the instant where the reaction coordinate reaches $\varPhi =1$. We have obtained 1590 trajectories from AMS and $189$ trajectories form DNS. We firstly compute the sample mean $\langle \tau \rangle _o$ and the sample variance $\sigma _\tau$ over AMS realisations of the trajectory durations, we obtain $\tau _{AMS}=413\pm 3$, $\tau _{DNS}=410\pm 7$ and $\sigma _{\tau,AMS}=105$ and $\sigma _{\tau,DNS}=91$. As stated in § 2.2.3, the error bars on $\tau$ are given by $\sigma _\tau$ divided by the square root of the number of samples. (It is safe enough to assume each value is decorrelated enough from a 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 with the sample size.) We can see that both are within a few per cent of one another and within the $66\,\%$ confidence interval of each other. This shows that the reactive trajectories computed by AMS and DNS have almost the same duration, and thus that AMS has very little bias in the way it computes the collapse trajectories. We then compare the sampled distributions of trajectory durations. In figure 7, we display the distribution of normalised durations $\tilde {s}={(s-\tau )}/{\sigma _\tau }$ of collapse trajectories computed by means of AMS and of DNS. Both are very similar and seem to be very close to a Gumbel distribution. This distribution of duration has been reported to be very close to a Gumbel in many instances of stochastically driven systems (Rolland & Simonnet Reference Rolland and Simonnet2015; Rolland et al. Reference Rolland, Bouchet and Simonnet2016; Rolland Reference Rolland2018). 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 trajectory durations have this distribution in one-dimensional stochastic systems (Cérou et al. Reference Cérou, Guyader, Lelievre and Malrieu2013). To our knowledge, this is possibly the first time that such a distribution shape has been shown in turbulence without stochastic injection of energy. This suggests us that the manner in which turbulence collapses in systems of large enough size can be compared with 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.

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

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

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

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

(3.1)\begin{equation} \rho_x(\log(E_{k,x}))=\int \rho_e(\log(E_{k,x}),\log(e_{k,y-z}))\,\textrm{d}\log(e_{k,y-z}), \end{equation}

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

(3.2)\begin{equation} \rho_c\left(\log(e_{k,y-z})|\log(E_{k,x})\right)=\frac{\rho_e(\log(E_{k,x}), \log(e_{k,y-z}))}{\rho_x(\log(E_{k,x}))}, \end{equation}

and the conditional average is

(3.3)\begin{equation} \langle \log(E_{k,y-z})\rangle (\log(E_{k,x}))=\int \rho_c\left(\log(e_{k,y-z})|\log(E_{k,x})\right)\log(e_{k,y-z})\,\textrm{d}\log(e_{k,y-z}). \end{equation}

We can also compute the conditional variance

(3.4)\begin{align} & \sigma_{y-z}(\log(E_{k,x})) \nonumber\\ &\quad =\sqrt{\int \left(\rho_c(\log(e_{k,y-z})|\log(E_{k,x}))\log(e_{k,y-z})^2\right)\,\textrm{d} \log(e_{k,y-z})-\langle \log(E_{k,y-z})\rangle^2 (\log(E_{k,x}))}. \end{align}

We display the average $\langle \log (E_{k,y-z})\rangle (\log (E_{k,x}))$ in full line and the average plus/minus the variance $\langle \log (E_{k,y-z})\rangle (\log (E_{k,x}))\pm \sigma _{y-z}(\log (E_{k,x}))$ by the dashed line for all three types of reactive trajectories in figure 8(d). This quantitatively confirms the observation of a faster decay of the streamwise vortices, and thus of $E_{k,y-z}$, than of the streaks, and thus of $E_{k,x}$. This also confirms that using the AMS leads to the same collapse trajectories as using DNS, at a much smaller cost, up to the effect of the noise added at branching on $E_{k,y-z}$ when it is very small. We can finally note that the trajectories seem to be concentrated around this average. Note, however, that the constant variance of $\log (E_{k,y-z})$ actually corresponds to a variance of $E_{k,y-z}$ which decreases exponentially as the collapse of turbulence goes on. There is a larger variance of $E_{k,y-z}$ among the initial conditions and thus in the amplitude and distribution of streamwise vortices in the initial conditions and stronger fluctuations in the beginning of the trajectory. On a final note, we state that, in order to rigorously obtain instantons in stochastic systems, one has to take a limit of the variance of the added noise going to $0$. In the case of our turbulent flows, this means either identifying an effective noise (which is a difficult task) or displaying large deviations in quantities like the probability density functions, the probability of collapse or the mean first passage time. Quantitatively speaking, this amounts to identifying a small parameter $\epsilon$, proportional to the square of the suspected noise variance such that $\lim _{\epsilon \rightarrow 0}\epsilon \log (\alpha )$ or $\lim _{\epsilon \rightarrow 0}\epsilon \log (T)$ are independent on $\epsilon$ (these two limits often have opposite signs) (Touchette Reference Touchette2009). As it happens, such scalings have been displayed in probability density functions of kinetic energy in simulations of plane Couette flow when the sizes $L_x$, $L_z$ were increased so that domains contained up to 6 wavelengths of the laminar turbulent bands, with $\epsilon _{DNS}={1}/{L_xL_z}$ (Rolland Reference Rolland2015). In that case, the rate functions of the PDF of kinetic energy showed the tails for $E_k<\langle E_k\rangle$, describing the fluctuations of kinetic energy toward low values and eventually toward laminar flow. This type of scaling was also shown in mean first passage times before collapse of laminar–turbulent coexistence in models of pipe flow where the length $L$ of the system was increased so that it contained up to ten puffs, with $\epsilon ={1}/{L}$, (Rolland Reference Rolland2018). In that case, collapse trajectories started from laminar–turbulent coexistence. This situation can have a highly fluctuating dynamics, especially at Reynolds numbers slightly above the transition threshold, with events like single puff decay or splitting. However, the collapse trajectories themselves followed a rather systematic scenario (collapse puff by puff for Reynolds numbers slightly above the transition threshold and synchronous collapse at higher Reynolds numbers) and showed typical paths in phase space constituted by global quantities.

One thus observes large deviations in the large system size limit, and in that limit one should expect to observe an increasingly narrower concentration of trajectories, and an increase of the mean first passage time before collapse which is exponential with the length of the system. Note that models indicate that nothing in laminar–turbulent coexistence prevents the collapse of turbulence: if the Reynolds number is increased slightly above the threshold of transition, this only requires the concomitant collapse of all the puffs. This condition is at the origin of the exponential dependence of the lifetime of turbulence on size. It may very well be that the very same thing happens in actual pipe flows or Couette flows, where one would require the rare, but not impossible, concomitant collapse of all the puffs or of the bands of the flow. This synchronisation could be at the origin of an exponential scaling in size. Note that, in that respect, there is a similar phenomenology in our rectangular boxes and larger domains: in both cases, the collapse of turbulence on a streamwise long corridor is required, and in both cases, that collapse can be countered by contamination from the sides of the hole. This would indicate that laminar turbulent coexistent is strictly permanent only in domains of infinite size.

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

We remark that the dynamics leading to collapse of turbulence can be different from the typical laminar–turbulent coexistence. In plane Couette flow, the laminar–turbulent coexistence can be fairly complex, in particular at the lowest Reynolds number where it is permitted: in flows extended in two dimensions, one can observe growth or decay of turbulent spots or stubs of bands. These can be separated in two, drift etc. (see Chantry, Tuckerman & Barkley Reference Chantry, Tuckerman and Barkley2017). Visually, this can give fairly distinct initial conditions for the collapse of turbulence, which nevertheless have comparable statistical properties. In these cases, the most probable route to laminar flow will most likely remain to suppress spots or band stubs one by one at lower Reynolds numbers, and collapse them all in a synchronous manner at higher Reynolds numbers. Time series of global quantities, such as turbulent fraction, kinetic energy of different components etc. would then have a very similar structure.

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 it will generate other collapse trajectories by branching on these collapse trajectories, it generates even more variability than this through mutation/selection on other initial conditions, with the help of anticipated branching. This means that each run generates distinct, independent reactive trajectories, in particular trajectories with distinct initial conditions. These trajectories are even more distinct in their later stages. The amount of variability generated can be estimated by counting the number of genealogies, that is to say the number of group of trajectories that share the same starting point. Each group will display trajectories decorrelated from those of another genealogy. On average, we generate 24 genealogies per AMS run with $120$ clones. This value is obtained by averaging over all available AMS runs the number of distinct initial conditions in computed reactive trajectories. See Ragone & Bouchet (Reference Ragone and Bouchet2021) for a view of genealogies of trajectories with a different algorithm applied to climate.