Rare transitions to thin-layer turbulent condensates

TL;DR

This study investigates rare transitions between 2-D condensate and 3-D flow states in thin-layer turbulence, revealing exponential distributions of transition times and modeling the dynamics with a stochastic Langevin equation near the critical layer height.

Contribution

It provides a detailed statistical analysis of bistable transitions in thin-layer turbulence and introduces a stochastic model to describe the large-scale energy dynamics near the critical threshold.

Findings

Transition times follow exponential distributions with mean values increasing rapidly near the threshold.

The large-scale kinetic energy dynamics can be modeled by a stochastic Langevin equation.

Effective potential analysis reveals two minima corresponding to 2-D and 3-D states.

Abstract

Turbulent flows in a thin layer can develop an inverse energy cascade leading to spectral condensation of energy when the layer height is smaller than a certain threshold. These spectral condensates take the form of large-scale vortices in physical space. Recently, evidence for bistability was found in this system close to the critical height: depending on the initial conditions, the flow is either in a condensate state with most of the energy in the two-dimensional (2-D) large-scale modes, or in a three-dimensional (3-D) flow state with most of the energy in the small-scale modes. This bistable regime is characterised by the statistical properties of random and rare transitions between these two locally stable states. Here, we examine these statistical properties in thin-layer turbulent flows, where the energy is injected by either stochastic or deterministic forcing. To this end, by using a large number of direct numerical simulations (DNS), we measure the decay time $\tau_d$ of the 2-D condensate to 3-D flow state and the build-up time $\tau_b$ of the 2-D condensate. We show that both of these times $\tau_d,\tau_b$ follow an exponential distribution with mean values increasing faster than exponentially as the layer height approaches the threshold. We further show that the dynamics of large-scale kinetic energy may be modeled by a stochastic Langevin equation. From time-series analysis of DNS data, we determine the effective potential that shows two minima corresponding to the 2-D and 3-D states when the layer height is close to the threshold.

Figures

Equations

$${\mathbf{f}}_d(\mathbf{x},t)=ϵ\sum _{\mathbf{k}\in K_F}\frac{{\hat{\mathbf{v}}}_{2D}(\mathbf{k})e^{i\mathbf{k}\cdot \mathbf{x}}}{{\sum }_{{\mathbf{k}}^{\prime }\in K_F}|{\hat{\mathbf{v}}}_{2D}({\mathbf{k}}^{\prime }){|}^2}+i(ϵk_f^2{)}^{1/3}\sum _{{\mathbf{k}}^{\prime }\in K_F}{\Omega }_{\mathbf{k}}{\hat{\mathbf{v}}}_{2D}(\mathbf{k})e^{i\mathbf{k}\cdot \mathbf{x}},$$

$$V_{ls}^2=\sum _{\genfrac{}{}{0ex}{}{\mathbf{k},k_z=0}{|\mathbf{k}|<k_{max}}}\left[|{\hat{v}}_x(\mathbf{k}){|}^2+|{\hat{v}}_y(\mathbf{k}){|}^2\right],$$

$$\Delta E(t)=-\frac{\partial U(E)}{\partial E}\Delta t+\sqrt{2B(E)}\Delta W(t),$$

$${\tau }_b=\frac{2}{a+b}\left[E_{\min }-E_{\max }+\frac{b}{a}{E}_{\max }\left(-1+(E_{\max }/E_{\min }{)}^{a/b}\right)\right],$$

$$p(E,t+\Delta t|E^{\prime },t)-\delta (E-E^{\prime })=-\frac{\partial }{\partial E}A(E)\delta (E-E^{\prime })\Delta t+\frac{{\partial }^2}{\partial E^2}B(E)\delta (E-E^{\prime })\Delta t+O(\Delta t^2)$$

$$\int dE(E-E^{\prime })p(E,t+\Delta t|E^{\prime },t)=$$

$$\int dE{\left(E-E^{\prime }\right)}^{2}p(E,t+dt|E^{\prime },t)=$$

Full text

1 Introduction

Turbulence is ubiquitous in the universe, from stars to tea cups. Many astrophysical and geophysical turbulent flows, such as planetary oceans and atmospheres, are subject to geometrical constraints, e.g. thinness in one spatial direction (Pedlosky, 2013). Such constraints significantly change the properties of the flow, which therefore deviate from those of classical three-dimensional (3-D) homogeneous and isotropic turbulence. Fully 3-D turbulence is characterised by a forward cascade of energy from large to small scales (Frisch, 1995), while in two dimensions (2D), an inverse energy cascade from small to large scales occurs due to additional inviscid invariants such as enstrophy (Boffetta & Ecke, 2012). Turbulence in thin layers combines properties of both cases, as large-scale dynamics are constrained to be 2-D, while small-scale dynamics are not (Smith et al., 1996; Celani et al., 2010; Benavides & Alexakis, 2017; Musacchio & Boffetta, 2017). As a consequence, in thin-layer turbulence, energy may cascade both to small and large scales depending on the layer height $H$: Above a critical height $H_{{}_{3D}}$ there is no inverse cascade, while below this critical height an inverse cascade develops, whose amplitude (measured by the inverse energy flux) increases continuously. Furthermore, a second critical height $H_{{}_{2D}}<H_{{}_{3D}}$ was observed where the flow become exactly two-dimensional and no forward cascade was observed. Similar transitions towards an inverse cascade occur in rotating turbulence (Smith & Waleffe, 1999; Deusebio et al., 2014), stratified turbulence (Métais et al., 1996; Marino et al., 2013, 2014; Sozza et al., 2015) and magnetohydrodynamic systems (Alexakis, 2011; Seshasayanan & Alexakis, 2016; Seshasayanan et al., 2014), among others. (see the review articles (Alexakis & Biferale, 2018; Pouquet et al., 2018)).

In a finite domain, an inverse cascade saturates at late times, forming a condensate in which the energy is concentrated in the largest scales. This condensation has been extensively studied in 2-D turbulence (see Hossain et al., 1983; Smith & Yakhot, 1993, 1994; Chertkov et al., 2007). In quasi-2D systems it has been observed in rapidly rotating convection (Rubio et al., 2014; Favier et al., 2019), rotating turbulence (Alexakis, 2015; Yokoyama & Takaoka, 2017; Seshasayanan & Alexakis, 2018) and thin-layer turbulence (Xia et al., 2011; van Kan & Alexakis, 2019). In many of these cases, the amplitude of the condensate state (measured by the energy in the large scales) has been shown to vary discontinuously with the system parameters. Furthermore, close to the transition, bistability has been observed: the system was either attracted or not to the condensate state depending on the initial conditions. In particular in thin-layer flows, for values of $H$ close to $H_{{}_{3D}}$, the system was attracted to either a 2-D condensate state (where most of the energy is concentrated in two counter-rotating, large-scale, 2-D vortices) or a 3-D flow state (where energy is mostly contained in 3-D small-scale fluctuations) (van Kan & Alexakis, 2019; Musacchio & Boffetta, 2019). The bistability in this system was accompanied by sudden ‘jumps’ between these two states. These transitions occur randomly with the waiting times that are, presumably, stochastic, following a statistical distribution that characterises the bistable regime.

In this paper, we present the first analysis of the statistical properties of thin-layer turbulence close to the critical height. We use a very large number of direct numerical simulations (DNS) and calculate the probability distribution functions (PDFs) of the transition times: the decay time $\tau_{d}$ from a 2-D condensate state to a 3-D flow state and the build-up time $\tau_{b}$ from a 3-D flow state to a 2-D condensate state. We examine their dependence to $H$ and attempt to model the transitions in terms of a particle in a one-dimensional potential using a Langevin equation.

2 Setup and results from direct numerical simulations

In this study, we consider forced incompressible 3-D flow in a triply periodic domain of dimensions $L\times L\times H$ with $H\ll L$. The setup is identical with the one studied in (van Kan & Alexakis, 2019). The thin direction is referred to as the vertical ‘$z$’ direction and the remaining two as the horizontal ‘$x,y$’ directions. The flow obeys the incompressible Navier-Stokes equation {subeqnarray} ∂_t v + v⋅∇v =& - ∇P + ν∇^2 v + f ,

∇⋅v = 0,

where $\mathbf{v}(\mathbf{x},t)$ is the velocity field, $P(\mathbf{x},t)$ is physical pressure divided by constant density, $\nu$ is kinematic viscosity and $\mathbf{f}$ is the external body force injecting energy into the flow. In this work, we use two different forcing functions: stochastic $\mathbf{f}_{s}$ and deterministic $\mathbf{f}_{d}$. Both forcing functions depend only on $x$ and $y$ and have only $x$ and $y$ components, i.e., are two-dimensional-two-component (2D2C) fields. In both cases, the force is divergence-free and only acts on a shell of wavenumbers $|\mathbf{k}|=k_{f}=2\pi/\ell$. The stochastic force is delta-correlated in time, which leads to a fixed mean injection rate $\langle\mathbf{v\cdot\mathbf{f}_{s}\rangle=}\epsilon$.

The deterministic force $\mathbf{f}_{d}$ is written in terms of the Fourier transform of the velocity field, $\mathbf{\hat{v}}(\mathbf{k})$, as

[TABLE]

where the sums are taken over all modes in the set $K_{F}$ that consists of all wavenumbers with $k_{z}=0$ and $k_{f}\leq|{\bf k}|<k_{f}+1$ and only in-plane components $\hat{\mathbf{v}}_{2D}=(\hat{v}_{x},\hat{v}_{y},0)$ are forced. There are two terms in this forcing. The first term is responsible for injecting energy at a fixed rate $\epsilon$, matching the mean injection rate of $\mathbf{f}_{s}$. The second term is not injecting energy, but is responsible for de-correlating the phases of the forced modes at a timescale given by $(\epsilon k_{f}^{2})^{-1/3}$. This is achieved by the linear term $i(\epsilon k_{f}^{2})^{1/3}\Omega_{\mathbf{k}}\mathbf{\hat{v}(k)}$ where $\Omega_{\mathbf{k}}$ are time-independent random numbers that are uniformly distributed over $[-1,1]$. They are kept fixed throughout the simulation time and are the same for all simulations. Since $\Omega_{\mathbf{k}}$ are time-independent, the forcing is indeed deterministic, i.e. is fully determined by the velocity field at any time. Without the second term, the flow at late times is attracted to a degenerate state where the velocity and the forcing are strongly correlated so that two strong vortices with opposite signs, of the size of the forcing scale, are formed. This forcing (1) has previously been used successfully by (Benavides & Alexakis, 2017).

For both forcing functions, the system is characterised by three non-dimensional parameters: the injection scale Reynolds number $\Rey=(\epsilon\ell^{4})^{1/3}/\nu$, the ratio between forcing scale and box height $Q=\ell/H$, and the ratio between forcing scale and the horizontal domain size $K=\ell/L$. In all simulations, we focus on the horizontal large-scale kinetic energy, defined as

[TABLE]

where $k_{max}=\sqrt{2}\frac{2\pi}{L}$. The simulations performed for this work use an adapted version of the Geophysical High-Order Suite for Turbulence (GHOST) which uses pseudo-spectral methods including 2/3 de-aliasing to solve for the flow in the triply periodic domain, (see Mininni et al., 2011). For all experiments, we fix $K=1/8$ and $\Rey=203$ at a resolution of $128\times 128\times 16$, varying $Q$ over the interval [{\color[rgb]{0,0,1}1.4},2.0]. We choose this low resolution and $\Rey$ because very long-duration runs are needed for this study. In addition, we also made runs at a resolution $256\times 256\times 16$ at $\Rey=406$, which qualitatively showed the same dynamics, even though reliable statistical analysis was not done in this case due to longer CPU times required to integrate (2).

Depending on $Q$ and the initial conditions, the system is attracted either to the 2-D condensate state or the 3-D flow state. Two snapshots of the vorticity for these two states are displayed in figure 1 where in the left panel the flow is in the 3-D state and has a small value of $V_{ls}$ while in the right panel the flow is in the 2-D condensate state and has large value of $V_{ls}$. Figure 2 shows a sketch of these states in the $Q-V_{ls}$ plane. The upper branch corresponds to the 2-D condensate state while the lower branch corresponds to the 3-D flow state. The dash-dotted lines indicate regions where one of the two states is unstable and the flow can jump from that state to the other. As detailed below, we perform decay and build-up experiments for each value of $Q$. In the build-up experiments, initial conditions corresponding to the 3-D flow state are used (see figure 2). In these simulations, we observe the build-up of 2-D condensates after a certain simulation time, which we denote by $\tau_{b}$. In the decay experiments, initial conditions corresponding to the 2-D condensate state are used. The system is evolved until the 2-D condensate decays, where the decay time is denoted by $\tau_{d}$. When decay or build-up events occur, the integration is interrupted and the next independent experiment is initiated. For the stochastic forcing, runs are started from a fixed initial condition but with different random number sequences, while for deterministic forcing, the initial conditions are altered by a small random perturbation that is different for every run. Figure 3 illustrates typical time-series for build-up and decay experiments. We note that similar procedures have been used to determine a critical Reynolds number for turbulent-laminar transitions in pipe flows (see Avila et al., 2011).

3 Build-up and decay times

We measure the statistical properties of the decay and build-up times $\tau_{d},\tau_{b}$ by simulating the system for more than ten million ($10^{7}$) eddy turnover times, using a total of 2.5 million CPU hours. This corresponds to between 40 and over 1000 independent runs per value of $Q$. Figure 4 shows PDFs of $\tau_{d}$ and $\tau_{b}$ for the stochastic forcing $\mathbf{f}_{s}$ for different values of $Q$. All PDFs have an exponential tail, whose slope (in log-linear scale) increases as the transition is approached. The PDFs for the deterministic forcing $\mathbf{f}_{d}$ show qualitatively the same results. The PDFs remain almost unchanged when half of the data points in the sample are removed, which indicates that they are sufficiently well sampled. An exponential PDF of waiting times implies that the waiting mechanism can be modeled using a memoryless process, in which long-time correlations are absent (Billingsley, 2008). We explore this possibility in Section 4.

The resulting mean build-up and decay times, $\langle\tau_{b}\rangle,\langle\tau_{d}\rangle$, are shown in figure 5 for stochastic (left panel) and deterministic (right panel) forcing in log-linear coordinates. The uncertainty in the transition times due to the finite sample size is estimated based on exponential PDFs and shown by the error bars. The ascending branch (left, in blue) represents the mean decay time $\langle\tau_{d}\rangle$ and the descending branch (right, in red) represents the build-up time $\langle\tau_{b}\rangle$. For each case, both transition times increase drastically when a certain value of $Q$ is approached. This increase is faster than exponential (super-exponential) and could be either diverging at some critical values $Q_{b}$ and $Q_{d}$ or staying finite (similar to what is observed for the decay and split time of turbulent puffs in pipe flows (Hof et al., 2006, 2008; Avila et al., 2011)).

In the former case, where the two times diverge, for $Q<Q_{b}$ the 2-D condensate will eventually decay to the 3-D flow state, while for $Q>Q_{d}$ the 3-D state will transition to a condensate. Thus $Q_{d}$ corresponds to the smallest value of $Q$ for which a condensate can be observed at late times, while $Q_{b}$ is the largest value of $Q$ at which 3D-flow state can be observed. Both $Q_{b}$ and $Q_{d}$ are close to the value $Q_{3D}\propto 1/H_{{}_{3D}}$ that gives the onset of the inverse cascade as observed in (Benavides & Alexakis, 2017; van Kan & Alexakis, 2019). We note however that in (Benavides & Alexakis, 2017), where the presence of a hypo-viscosity and the large box limit was considered, no bistability was observed, while in (van Kan & Alexakis, 2019) due to the limited statistics in that work the distinction between $Q_{b}$ and $Q_{d}$ could not be made. In the latter case of transition times staying finite despite the super-exponential increase, a double-exponential function might be used to fit this super-exponential increase, supported by an argument using extreme value statistics (Fisher & Tippett, 1928; Gumbel, 1935; Goldenfeld et al., 2010; Goldenfeld & Shih, 2017). This can be justified if we assume that the transition to the condensate state is triggered when the maximum value of the small-scale vorticity exceeds a certain threshold value. Testing this assumption requires knowledge of the entire flow field at times near rare transition events. We believe that it would be important to clarify this question in a future study.

Four different scenarios may be envisaged for the $Q$-dependence of the transition times, illustrated in figure 6. In the first scenario, both $\langle\tau_{b}\rangle$ and $\langle\tau_{d}\rangle$ diverge at $Q_{b},Q_{d}$ with $Q_{d}<Q_{b}$. In the range $(Q_{d},Q_{b})$, where both transition times diverge, either the 3-D flow or 2-D condensate state is selected, depending on initial conditions. In the second scenario, $\langle\tau_{b}\rangle$ and $\langle\tau_{d}\rangle$ diverge at the same point, i.e., $Q_{d}=Q_{b}$. This is an analogue of standard equilibrium phase transitions with a single power-law singularity. In the third scenario, a crossover is observed, i.e., $Q_{b}<Q_{d}$, where random transitions between the two states are possible within a finite range of $Q$ between $Q_{b}$ and $Q_{d}$. Finally, in the fourth scenario, a crossover is again observed, but without any divergence of $\langle\tau_{b}\rangle$ and $\langle\tau_{d}\rangle$ for finite $Q$. This scenario is compatible with the double-exponential fitting function explained above, where transitions between the two states are possible for all $Q$ (although they are extremely rare). From the data in figure 5, we can see that for the stochastic forcing, the two branches are intersecting around $Q=Q_{s}\approx 1.55$ at a value of around $10^{5}$ eddy turnover times. Therefore we can exclude cases 1 and 2 for $\mathbf{f}_{s}$, while all four scenarios are possible for $\mathbf{f}_{d}$. In addition, it may be proven (see Gallet & Doering, 2015) and has been confirmed numerically (van Kan & Alexakis, 2019), that beyond a second critical value $Q_{2D}$, the flow two-dimensionalises and the condensate is stable to 3-D perturbations in the long-time limit. Therefore, condensate decay time diverges at least for $Q>Q_{2D}$ and we can exclude a double-exponential behavior extending to all $Q$ for $\langle\tau_{d}\rangle$. Note that based on our available data, we cannot exclude that condensate build-up time never becomes infinite. In a future study, rare event algorithms (see Section 5) may help elucidating these questions, including the problem of how rare transitions are triggered.

An alternative way of reducing the computational cost of the problem could be to consider not the statistics of transition events which become extremely rare, but rather the statistics of bursting events, which can be seen in figure 3 to be significantly more frequent than transition events in the bistable regime. While we have not studied them in very much detail, we have found that waiting times for bursting events also follow exponential PDFs, whose mean increases towards larger $Q$. Bursting events may prove more fruitful than transition events at higher Reynolds numbers and we believe they are a promising object of study for future work.

4 Effective Markovian modelling

The exponential PDFs of the transition times (figure 4) indicate that these times are stochastically determined by a mechanism that is memoryless, in other words, which leads to vanishing long-time correlations for these times. Since the transitions are quantitatively characterized by a single macroscopic variable, the horizontal large-scale kinetic energy $V_{ls}^{2}$ (figure 3), this observation implies that the dynamics of $V_{ls}^{2}$ could be described by a Markov (memory-less) process, such as an inertia-less particle moving in a double-well potential in the presence of white noise. Within this effective description, the transitions are characterized as rare jumps of the particle between the two wells of the potential.

Motivated by this observation, in this section, we discuss to what extent the dynamics of $E\equiv V_{ls}^{2}/(\epsilon\ell)^{2/3}$ can be described by one-dimensional Markov process. Assuming the continuity of the trajectory of $E(t)$, a general form of this process is written as the following Langevin equation (Gardiner, 1986): denoting by $\Delta E(t)$ the increment of $E(t)$ during a small time interval $\Delta t$, the equation is written as

[TABLE]

where $\Delta W(t)$ is a Gaussian noise satisfying $\langle\Delta W(t)\rangle=0$ and $\langle\Delta W(t)\Delta W(t^{\prime})\rangle=\Delta t\ \delta_{t,t^{\prime}}$ and we use the Itô rule for the multiplication between $\Delta W(t)$ and $\sqrt{2B(E)}$, i.e., $\langle\sqrt{2B(E)}\Delta W(t)\rangle=0$. The potential $U(E)$ and the $E$-dependent diffusion constant $B(E)$ characterise the dynamics and are to be determined. Note that the noise term models the smaller-scale turbulent motions. Eq. (3) indicates that the average and the variance of $\Delta E(t)$ conditional on $E$ are equal to $-\left(\partial U(E)/\partial E\right)\Delta t$ and $2B(E)\Delta t$ respectively. From the time-series data of the DNS, we thus evaluate these statistical properties and estimate $U(E)$ and $B(E)$ (see Appendix A for more detail). The results are shown in figures 7 and 8. When $Q$ is well below the transition, $U(E)$ has a single minimum at $E\simeq 0.1$. Positive slopes at larger values of $E$ indicate a mean drift towards the 3-D flow state. As $Q$ approaches the transition, $U(E)$ becomes flatter, meaning that this drift vanishes, and a second minimum appears at $E\simeq 7$ for $Q\simeq 1.5$. This indicates that the system is in a bistable regime. The second minimum becomes more pronounced as $Q$ is further increased. Eventually for $Q\gtrsim 1.65$ the first minimum disappears and the system is left with a single minimum corresponding to the condensate state. The time scale for the system to jump from one minimum to the other can be estimated using a first-passage time formula in terms of $U(E)$ and $B(E)$ [see Section 5.2.7 “First Passage Times for Homogeneous Processes” in (Gardiner, 1986)]. For example, we consider a transition time for build-up events at $Q=1.563$, where the corresponding $U(x)$ and $B(x)$ are shown in the right panel of figure 7 and in the left panel of figure 8, respectively. The build-up event is then characterized as the jump from the potential (local) minimum around $E_{\rm min}=0.3$ to the potential maximum around $E_{\rm max}=1.8$. In this range of $E$, $B(E)$ can be well-described as a linear function $B(E)=bE$ (with $b\simeq 0.001$). We also approximate $U(E)$ as a linear function $U(E)=aE+c$ (with $a\simeq 0.002$) for simplicity111A generalization for more realistic potentials is straightforward (Gardiner, 1986). Here we choose this approximation because the obtained formula (4) becomes simple.. Using the first-passage time formula222We set a reflecting boundary at $E=E_{\rm min}$ and an absorbing boundary at $E=E_{\rm max}$. Substituting the explicit expressions of $B(E)$ and $U(E)$ into the mean first passage time formula in p.139 of (Gardiner, 1986), we get (4). Note that, for more detailed analysis, we need to model $B(E)$ and $U(E)$ outside $E_{\rm min}<E<E_{\rm max}$, and define a reflecting boundary at $E=0$ and an absorbing boundary at another local minimum of $U(E)$ corresponding to the 2-D stable state. This will possibly increase the estimation of $\tau_{b}$ compared with (4). , we get the time scale $\tau_{b}$ as

[TABLE]

which estimates $\tau_{b}=2\times 10^{4}$. This shows a rough agreement with the data in the left panel of figure 5 in the previous section.

The interpretation of our results using a Langevin equation with a double-well potential is the simplest model to describe a discontinuous phase transition. The presence of a double-well potential implies that the system can always jump from one state to the other if one waits long enough. However we cannot conclude that $\langle\tau_{d}\rangle$ and $\langle\tau_{b}\rangle$ do not diverge for finite $Q$ (Scenario 4 in figure 6) from this observation. First of all, the range of $Q$ that we have studied is highly limited. Second, an oversimplification has been made when we assumed that the system could be described by a single variable $E$, while the system in reality evolves in an extremely high-dimensional space. Indeed, disagreements between the 1D model (3) and DNS results can be detected when we look at the distribution of the energy increment $\Delta E$ over a small time interval $\Delta t$. According to the Langevin equation (3), this distribution needs to be Gaussian, but this is not exactly true for the DNS results: using the time-series data of DNS, we evaluate the transition probability $p(E,t+\Delta t|E^{\prime},t)$ from the state $E^{\prime}$ at time $t$ to $E$ at time $t+\Delta t$. In figure 8 we plot this probability with $E^{\prime}=0.5$ and $\Delta t=0.2$ for stochastic forcing with $Q=1.5625$, showing deviations from a Gaussian distribution.

5 Conclusions

In this work we have studied the statistical properties of thin-layer flows close to the transition between a 3-D flow state and the formation of a 2-D condensate. Such transitions have recently been discovered in a variety of systems (Seshasayanan & Alexakis, 2018; Favier et al., 2019; van Kan & Alexakis, 2019) and this work is the first attempt to systematically study their statistics. We have measured the probabilities of the transition times between the two states, where the mean transition times were shown to increase by three orders of magnitude by a relatively small change (10%) of the control parameter $Q$. We point out the qualitative similarity between figure 5 of (Avila et al., 2011) for the turbulent-laminar transitions in pipe flows and figure 5 of our work. Although the physical situations are different, in both cases we observe super-exponentially growing time scales of two competing processes. Our results could neither exclude nor confirm whether this sharp increase has a double-exponential scaling form $\exp(\beta\exp(\alpha x))$, supported by an argument based on extreme value statistics (Goldenfeld et al., 2010). This leaves a possibility of the divergence of mean transition times, i.e., transitions from one state to the other could become impossible for a certain range of $Q$.

Our results show that the system can be modeled to some extent as an inertia-less particle trapped in a one-dimensional potential in the presence of stochastic noise. The model revealed that close to the transition the potential displays two minima implying the existence of a bistable state. Discrepancies in the noise statistics of the DNS and of the stochastic model were observed, that were attributed to the multi-dimensional nature of the real problem.

Several simplifying assumptions have been made in this work in order to render the problem more tractable. For example, the domain was triply periodic and also the forcing was 2-D. These simplifications could limit the applicability of these results to laboratory or natural flows. Further investigations with more realistic boundary conditions and forcing are thus necessary. More importantly, the present investigation was limited to a single, relatively small, value of $\Rey$ and $K$ (the scale separation between the forcing length scale and the horizontal domain scale). Examining the observed behaviour at larger values of $\Rey$ and $K$ is crucial for validating the robustness of our results and for their applications to natural flows. Another important limitation is the low vertical resolution of only 16 points imposed by the long simulation times. It is known that the minimum resolution required for correct DNS of turbulence is around 256, for a detailed study of resolution issues in turbulence see (Donzis et al., 2008). The impact of the low vertical resolution used here is that the small-scale vertical motion will be more strongly affected by viscosity. This likely affects the statistics of the rare transitions, the latter being induced by turbulent fluctuations. We use the term ”3-D flow state” instead of ”3-D turbulent state” to highlight this. Unfortunately, the rareness of decay and build-up events close to the bistable regime limits the range of parameters and the resolution that can be examined with DNS.

However in recent years, alternative methods have been developed that could address this problem. In particular, our understanding of the problem could benefit from studies using rare event sampling algorithms, such as a method calculating the instanton based on Freidlin-Wentzell theory (Chernykh & Stepanov, 2001; Heymann & Vanden-Eijnden, 2008; Grafke et al., 2015b, a; Grigorio et al., 2017), splitting methods that copy rare event realizations to efficiently accumulate statistics (Allen et al., 2005; Giardinà et al., 2006; Cérou & Guyader, 2007; Tailleur & Kurchan, 2007; Teo et al., 2016; Nemoto et al., 2016; Lestang et al., 2018; Bouchet et al., 2019) and also a recently proposed method that relies on feedback control of Reynolds number (Nemoto & Alexakis, 2018). Such studies can help to overcome the difficulty caused by the extremely long computational time required to accurately describe the rare transition events close to the onset. Studies at larger $\Rey$ or scale separations could therefore become tractable using these methods. In addition to more efficient simulation, statistics at higher Reynolds numbers may be accessible if one chooses not to study transitions but bursting events out of one of the attractors. Bursts are more frequent than transitions in the bistable regime and their statistics therefore are more easily accessible.

Furthermore, given the large experimental literature (see Xia et al., 2011; Xia & Francois, 2017) on the transition between 3-D turbulence and condensation in thin layers, it would be exciting and very important to study the observed bistability experimentally. The biggest advantage of an experiment compared to DNS would be that much longer observation times are possible as well as higher $\Rey$. The same remarks apply to rotating turbulence.

Acknowledgements

The authors gratefully acknowledge the comments of three anonymous referees, which have helped improve the clarity of the paper. This work was granted access to the HPC resources of MesoPSL financed by the Région Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by Agence Nationale pour la Recherche and HPC resources of GENCI-TGCC-IRENE & CINES Occigen (Project A0030506421 & A0050506421) where the present numerical simulations have been performed. This work has also been supported by the Agence nationale de la recherche (ANR DYSTURB project No. ANR-17-CE30-0004). We thank Freddy Bouchet for helpful discussions.

Appendix A Determining $U(E)$ and $B(E)$

In this appendix we express $U(E)$ and $B(E)$ in terms of transition probabilities to obtain expressions which allow determining them from the DNS time series. For convenience, let $A(E)\equiv-\frac{\mathrm{d}U}{\mathrm{d}E}$, the drift velocity. For a small time increment $\Delta t$, we denote by $p(E,t+\Delta t|E^{\prime},t)$ the transition probability from $E(t)=E^{\prime}$ at time $t$ to $E(t+\Delta t)=E$ at time $t+\Delta t$. It obeys, (see Gardiner, 1986),

[TABLE]

By multiplying both sides by $E-E^{\prime}$ or $\left(E-E^{\prime}\right)^{2}$ and integrating over $E$, we get

[TABLE]

The left-hand sides (and thus $A$ and $B$) are measurable from a time-series $E_{DNS}(t)$ of large-scale energy by computing the mean of $E_{DNS}(t+\Delta t)-E_{DNS}(t)$ and $(E_{DNS}(t+\Delta t)-E_{DNS}(t))^{2}$, over all times $t$ for which $E_{DNS}(t)=E^{\prime}$, for all values of $E^{\prime}$. Finally the potential $U(E)$ is obtained from $A(E)$ by simple integration.