2. Numerical approach

In order to study thin-layer turbulence, we consider the idealised case of forced incompressible 3-D flow in a triply periodic box of dimensions $L \times L \times H$, where the vertical direction is thin, $H\ll L$. The flow system is identical to that described in van Kan & Alexakis (Reference van Kan and Alexakis2019). We consider a Cartesian coordinate system $(x,y,z)$ with unit vectors $(\boldsymbol {e}_x,\boldsymbol {e}_y,\boldsymbol {e}_z)$, where the thin vertical direction is chosen along $\boldsymbol {e}_z$. The flow $\boldsymbol {u}(\boldsymbol {x},t)$ is then governed by the incompressible forced Navier–Stokes equations

(2.1a)$$\begin{gather} \frac{\partial\boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u}={-}\boldsymbol{\nabla} P+\nu\nabla^{2}\boldsymbol{u}+\boldsymbol{f}, \end{gather}$$

(2.1b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

where $P(\boldsymbol {x},t)$ denotes the pressure divided by the constant density and $\nu$ represents the kinematic viscosity of the fluid. We consider a stochastic forcing $\boldsymbol {f}(\boldsymbol {x},t)$ that is vertically invariant ($\partial \boldsymbol {f}/\partial z=0$) and acts exclusively in the $(\boldsymbol {e}_x,\boldsymbol {e}_y)$ plane, i.e. in the two-dimensional, two-component (2D2C) manifold. Furthermore, the forcing is divergence-free and acts only sharply on wavenumber $k_f \equiv 2{\rm \pi} /\ell$ (specifically, exclusively the modes $(k_x,k_y)=(\pm k_f,0)$ and $(0,\pm k_f)$ are forced), where its random phase is white noise (delta-correlated) in time. This results in a fixed mean injection rate $\langle \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {f} \rangle = \epsilon$ that is solely prescribed by the forcing amplitude (Novikov Reference Novikov1965). Here, $\langle {~}\rangle$ is used to represent the ensemble average. The choice of forcing is motivated by simplicity and comparability with previous studies. In general, one may consider various 3-D forcing functions (Poujol, van Kan & Alexakis Reference Poujol, van Kan and Alexakis2020).

The input parameters of the thin-layer flow system are combined to give three dimensionless numbers. We define an injection-scale Reynolds number $Re\equiv (\epsilon \ell ^{4})^{1/3}/\nu$, the ratio between the forcing scale and the domain height $Q\equiv \ell /H$, and the ratio between the forcing scale and the width of the domain $K\equiv \ell /L$. Finally, we define a forcing time scale $\tau _f\equiv (\ell ^{2}/\epsilon )^{1/3}$ and energy scale $E_f\equiv (\epsilon \ell )^{2/3}$ that are used to non-dimensionalise the different temporal and energetic quantities reported in this work, respectively.

Equations (2.1a) and (2.1b) are solved numerically in the triply periodic domain using a pseudo-spectral code that is an adapted version of the Geophysical High-Order Suite for Turbulence (Ghost) as introduced by Mininni et al. (Reference Mininni, Rosenberg, Reddy and Pouquet2011), employing 2/3 dealiasing. In order to investigate the dependence of the condensate transition and its bistable range on $Re$ and the box size, we take the results in van Kan, Nemoto & Alexakis (Reference van Kan, Nemoto and Alexakis2019) as a starting point and extend them to smaller and larger $Re$ and $K$. For each value of $K$ and $Re$, we vary the thinness of the fluid layer $Q$ as the principal control parameter in close vicinity to the condensate transition. An overview of the full set of input parameters that are considered in this work is provided in table 1.

Table 1.

The different series of input parameters used in this work for varying box size (left columns) and varying $Re$ (right columns).

Resolutions $N_x\times N_y\times N_z$ are chosen such that we maintain the same ratio between the grid spacing $L/N_{x,y}$ and Kolmogorov length $\eta =(\nu ^{3}/\epsilon )^{1/4}$ as used in van Kan et al. (Reference van Kan, Nemoto and Alexakis2019) of $L/N_{x,y}\approx 3.2\eta$ in the horizontal directions, and we resolve finer than that in the vertical direction. For the vertical, we ensure that we keep 16 grid cells in order to maintain sufficient degrees of freedom in the thin direction.

For each unique set of input parameters ($Q,K,Re$), at least 40 independent runs are carried out by using a different random seed for the stochastic forcing in order to obtain statistics about the rare transitions from one state into the other. The combination of this need for a multitude of independent runs with the fact that the large-scale dynamics is very slow compared to the background smaller-scale 3-D dynamics renders only the rather moderate parameters that are considered here computationally accessible. The simulations in this work comprise more than $\sim$1.0 million CPU hours.

The main diagnostic quantity that we use here to probe the strength of the condensate is the 2-D large-scale energy $E_{ls}$, defined from the Fourier components $\hat {\boldsymbol {u}}(\boldsymbol {k})$ of the flow as

(2.2)\begin{equation} E_{ls}=\frac{1}{2} \sum_{\substack{\boldsymbol{k}\\|\boldsymbol{k}|\leqslant k_{max}\\k_z=0}} [|\hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot}\boldsymbol{e}_x|^{2} +|\hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot}\boldsymbol{e}_y|^{2}], \end{equation}

where the cut-off wavenumber $k_{max}=\sqrt {2}(2{\rm \pi} /L)$.

3. Bistability, bimodality and rare transitions

Bistability and bimodality are often met in dynamical systems. Here we refer to ‘bistability’ as the presence of two independent stable attractors, coexisting for the same value of parameters; whilst by ‘bimodality’ we refer to the case in which two attractors are linked by some trajectories such that, when followed, the system jumps from one attractor to the other and vice versa. In the former case, a hysteresis loop can exist when one of the parameters of the system is continuously varied. In an inherently fluctuating dynamical system like the one at hand, however, the precise extent of the bistable range in the hysteresis loop can be ambiguous, as rare sudden transitions from one hysteretic branch into the other may exist at very long time scales near both ends of the hysteresis loop. This raises the question as to how we can unambiguously define the precise extent of the bistable range in the hysteretic transition.

Earlier works that studied the bistable range in the quasi-2-D condensate transition would typically simulate up to a certain time scale that is constrained by computational limits and call a state ‘stable’ if no further transitions are observed (Yokoyama & Takaoka Reference Yokoyama and Takaoka2017; Favier et al. Reference Favier, Guervilly and Knobloch2019; van Kan & Alexakis Reference van Kan and Alexakis2019; de Wit et al. Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022). However, in view of the rare transitions, this amounts to an in principle arbitrary cut-off at a certain time scale, neglecting any possible transitions occurring at larger time scales. We refer to this as finite-time hysteresis.

However, in the case of pure hysteresis, in the strict sense, all transitions from one state into the other are prohibited at any finite time scale, such that the system is absolutely bistable. This would require the time scales of the rare transitions to diverge at a certain asymptote. The existence of such asymptotes is the principal assumption in this work. As we will show in § 4, we can define such asymptotes based on the scaling of the time scales of the rare transitions that we observe, allowing us to study how these asymptotes shift as we vary the box size and $Re$ in order to get a completely time-scale-independent and unambiguous method for comparing our results at these different parameters.

To analyse the rare transitions between the two states, we separately consider build-up events from the 3-D state into the condensate state and decay events vice versa. The build-up events are studied by initialising the simulations with a tiny perturbation onto a state of no flow and continuing until the condensate state is reached. For the decay events, we initialise the simulation with a snapshot from a condensate state at higher $Q$ and we continue the run until the condensate has decayed and the 3-D state is obtained.

Figures 2(a) and 2(b) show examples of different realisations of such rare transitions for one choice of parameters. The works of van Kan et al. (Reference van Kan, Nemoto and Alexakis2019) and de Wit et al. (Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022) have revealed that the waiting time that is spent until these sudden transitions from one state into the other commence is exponentially distributed, signifying that the transition process is memoryless. The typical mean waiting time $\tau _W$ can be obtained by defining representative thresholds in the large-scale energy and analysing the distribution of times $t_{b,d}$ after which these thresholds are crossed. Examples of the obtained empirical cumulative distributions are depicted in figure 2(c), showing that it closely follows the aforementioned exponential distribution. We can then obtain $\tau _W$ by fitting the empirical cumulative distribution function (CDF) with a (shifted) exponential as

(3.1)\begin{equation} \textrm{CDF}(t_{b,d})=1-\exp\left(-\frac{t_{b,d}-\tau_0}{\tau_W}\right). \end{equation}

This process can be repeated to obtain the waiting time scales $\tau _W$ for build-up events and decay events at different values of $Q$, varying it across the full extent of the hysteretic transition. This results in a series of typical waiting times $\tau _W$ as a function of $Q$, which can in turn be repeated for different box sizes and $Re$.

Figure 2.

(a,b) Examples of time series of different realisations for (a) build-up at $Q=1.67$ and (b) decay at $Q=1.49$ (for $1/K=9$ and $Re=192$). Horizontal dashed lines represent the threshold energy at which the build-up time $t_b$ or decay time $t_d$ is defined and the simulation is terminated. The first blue run is continued for demonstration purposes. (c) Examples of distributions of $t_{b,d}$, shown through the empirical CDF (at $1/K=8$ and $Re=192$). Dashed lines represent fits of the exponential distribution.

4. Time-scale statistics

The results for the series of transition time scales for varying box size and varying $Re$ are provided in figure 3(a,c). As the transition is approached, $\tau _W$ increases faster than exponentially (see van Kan et al. Reference van Kan, Nemoto and Alexakis2019). This implies either that $\tau _W$ increases in a non-diverging super-exponential fashion – e.g. $\tau _W \propto \exp [\exp (Q)]$, as is typical for certain transitions controlled by extreme events (Goldenfeld, Guttenberg & Gioia Reference Goldenfeld, Guttenberg and Gioia2010; Nemoto & Alexakis Reference Nemoto and Alexakis2018, Reference Nemoto and Alexakis2021; Gomé, Tuckerman & Barkley Reference Gomé, Tuckerman and Barkley2022) – or that it diverges at some critical value $Q_0$. To determine which of the two holds for the present system is beyond the scope of this work. We will thus assume that the latter case applies, although one could alternatively interpret $Q_0$ as the value at which super-exponential behaviour starts in the former case. To determine $Q_0$, we fit the transition time to a power-law divergence

(4.1)\begin{equation} \tau_W\propto\frac{1}{|Q-Q_0|^{p}}. \end{equation}

By plotting $1/\tau _W^{1/p}$ as a function of $Q$, we can then obtain $Q_0$ from a linear fit to the data (figure 3b,d). Empirically, we find that $p^{{(build\text {-}up)}}=3$ and $p^{{(decay)}}=2$ result in a satisfactory linearisation of our data. These asymptotes then predict the location in $Q$ where the transition time becomes infinite, such that we can say that, beyond the asymptote, the transition cannot occur at any finite time scale.

Figure 3.

(a,c) Waiting times $\tau _W$ for build-up (circles) and decay (triangles) events and (b,d) their power-law transformation for varying box size $1/K$ at $Re=192$ (a,b) and varying $Re$ at $1/K=8$ (c,d). Crosses on the horizontal axis in (b,d) denote the estimates for the asymptotes $Q_0$. These asymptotes are also depicted in panel (c) by the vertical lines for build-up (dashed) and decay (dashed-dotted), but are omitted in panel (a) for readability.

Comparing the results at different box sizes and different $Re$, we find first of all that the transition is observed in a similar range of $Q$ for the varying box sizes, while it clearly shifts as $Re$ is varied, in agreement with the observations in van Kan & Alexakis (Reference van Kan and Alexakis2019). More importantly, we observe that the branches of build-up and decay times move further apart as $Re$ and box size are increased: the build-up branch moves to larger $Q$, while the decay branch moves to (relatively) smaller $Q$. The branches cross for the runs with $Re\leqslant 192$ or for $K\leqslant 8$, such that, at small box size and/or small $Re$, a bimodal range of $Q$ exists for $Q_0^{(build\text {-}up)}< Q< Q_0^{(decay)}$ where the build-up and decay time scales are simultaneously finite (and, in fact, computationally accessible). Hence, in this range, the flow continually transitions back and forth between the 3-D state and the condensate state.

Conversely, for the largest box size and largest $Re$ that we consider, the branches of the build-up and decay branches never cross, as the asymptotes reside on opposite ends. This indicates a profoundly different regime, where the decay times have diverged before the build-up times become finite, such that, in the range $Q_0^{(decay)}< Q< Q_0^{(build\text {-}up)}$, both states are absolutely stable, as no transition from one state into the other can occur in any finite time. This corresponds to a regime of pure hysteresis in which the system is (absolutely) bistable. Indeed, it is this range that is arguably the most unambiguous time-scale-independent definition of the bistable range of the system. These results are summarised in figure 4. It is thus evident from our results that the bistable range grows as the box size and/or $Re$ are independently increased.

Figure 4.

Phase diagrams of the condensate for (a)  varying box size at $Re=192$ and (b) varying $Re$ at $1/K=8$. Red and blue lines denote the asymptotes of waiting times for build-up and decay events, respectively.

We argue that the strengthening of the bistability for large box sizes and $Re$ is intuitive from the increase of the condensate energy level (compared to the 3-D-state energy) as $Re$ or $1/K$ is increased, making it harder to jump from one state into the other. Our results thus indicate that bistability is not a finite-size, finite-$Re$ effect, and such states could possibly be found in the geophysical limit where both $Re$ and domain size are large.