Subcritical turbulent condensate in rapidly rotating Rayleigh-B\'enard convection

TL;DR

This paper investigates the potential for subcritical transitions to large-scale vortex states in rapidly rotating Rayleigh-Bénard convection, revealing bistability and non-local energy transfer mechanisms in geostrophic turbulence.

Contribution

It demonstrates that subcritical large-scale vortex formation can occur in rotating convection even without prior inverse cascade evidence, highlighting bistability in turbulent flows.

Findings

Subcritical transition to large-scale vortices observed with finite amplitude initial conditions.

Bistability exists between vortex dipole and turbulent states.

Non-local inverse energy transfer can compete with direct cascade.

Abstract

The possibility of subcritical behaviour in the geostrophic turbulence regime of rapidly rotating thermally driven convection is explored. In this regime a non-local inverse energy transfer may compete with the more traditional and local direct cascade. We show that, even for control parameters for which no inverse cascade has previously been observed, a subcritical transition towards a large-scale vortex state can occur when the system is initialized with a vortex dipole of finite amplitude. This new example of bistability in a turbulent flow, which may not be specific to rotating convection, opens up new avenues for studying energy transfer in strongly anisotropic three-dimensional flows.

Figures

Tables

Table 1. Table 1: Summary of the parameters considered: T ​ a 𝑇 𝑎 Ta is the Taylor number, R ​ a 𝑅 𝑎 Ra is the Rayleigh number, R ​ o = R ​ a / ( P ​ r ​ T ​ a ) 𝑅 𝑜 𝑅 𝑎 𝑃 𝑟 𝑇 𝑎 Ro=\sqrt{Ra/(PrTa)} is the input Rossby number, R ​ a ~ = R ​ a ​ T ​ a − 2 / 3 ~ 𝑅 𝑎 𝑅 𝑎 𝑇 superscript 𝑎 2 3 \widetilde{Ra}=RaTa^{-2/3} is a scaled Rayleigh number and λ 𝜆 \lambda is the horizontal aspect ratio. The Prandtl number P ​ r = 1 𝑃 𝑟 1 Pr=1 in all simulations.

Case	$Ta$	$Ra$	$Ro$	$\widetilde{Ra}$	$\lambda$	$N_x\times N_y\times N_z$	$A$	Subcriticality
A1	${10}^8$	$3\times {10}^7$	$0.55$	$139$	$2$	${256}^3$	$\left[0:1000\right]$	Yes
A2	${10}^8$	$3\times {10}^7$	$0.55$	$139$	$4$	${512}^2\times 256$	$\left[0:1600\right]$	Yes
A3	${10}^8$	$3\times {10}^7$	$0.55$	$139$	$6$	${768}^2\times 256$	$800,1200$	Yes
B1	${10}^8$	$4\times {10}^7$	$0.63$	$186$	$2$	${256}^3$	$\left[1000:4000\right]$	No
B2	${10}^8$	$4\times {10}^7$	$0.63$	$186$	$4$	${512}^2\times 256$	$2000$	Yes
C1	${10}^8$	$5\times {10}^7$	$0.71$	$232$	$2$	${256}^3$	$\left[1000:4000\right]$	No
D1	${10}^8$	$4\times {10}^6$	$0.2$	$18.6$	$4$	${256}^2\times 128$	$\left[0:1000\right]$	No

Equations

$$\frac{1}{Pr}\left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p-\sqrt{Ta}{\mathbf{e}}_z\times \mathbf{u}+Ra\theta {\mathbf{e}}_z+{\nabla }^2\mathbf{u},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \theta }{\partial t}+\mathbf{u}\cdot \nabla \theta =w+{\nabla }^2\theta ,$$

$$\bm{u}=\Big{(}A\sin\left(2\pi y/\lambda\right),-A\sin\left(2\pi x/\lambda\right),0\Big{)}\quad\quad\text{and}\quad\quad\theta=0\ ,$$

$$K_0=\frac{1}{V}{\int }_V\frac{1}{2}\mathbf{u}\cdot \mathbf{u}\text{d}V=\frac{1}{2}{A}^2,$$

$$K(t)=K_0\exp (-8{\pi }^{2}Prt/{\lambda }^2).$$

$${⟨u⟩}_z(x,y)={\int }_0^{1}u(x,y,z)\text{d}z\text{and}{⟨v⟩}_z(x,y)={\int }_0^{1}v(x,y,z)\text{d}z,$$

$$u^{\prime }(x,y,z)$$

$$v^{\prime }(x,y,z)$$

$$K_{2D}=\frac{1}{2{\lambda }^2}\iint \left({⟨u⟩}_z^2+{⟨v⟩}_z^2\right)\text{d}x\text{d}y$$

$$K_{3D}=\frac{1}{2{\lambda }^2}\iiint \left(u^{\prime 2}+v^{\prime 2}+w^2\right)\text{d}x\text{d}y\text{d}z.$$

$$\frac{\partial {⟨\mathbf{u}⟩}_z}{\partial t}+{⟨\mathbf{u}⟩}_z\cdot {\nabla }_h{⟨\mathbf{u}⟩}_z$$

$$\frac{d{K}_{2D}}{dt}=\underset{\mathcal{D}}{\underbrace{\frac{Pr}{{\lambda }^2}\iint {⟨\mathbf{u}⟩}_z\cdot {\nabla }_h^2{⟨\mathbf{u}⟩}_z\text{d}x\text{d}y}}+\underset{\mathcal{F}}{\underbrace{\left(-\frac{1}{{\lambda }^2}\iint {⟨\mathbf{u}⟩}_z\cdot {⟨{\mathbf{u}}^{\prime }\cdot \nabla {\mathbf{u}}^{\prime }⟩}_z\text{d}x\text{d}y\right)}}.$$

Full text

1 Introduction

Turbulence in geophysical and astrophysical systems is a problem of major importance. Three-dimensional (3D) flows favour a forward cascade, i.e., energy flows from large to small scales, as described by the well-known Kolmogorov theory. In contrast, two-dimensional (2D) flows exhibit an inverse energy transfer, from small to large scales, that manifests itself in the appearance of large scale structures in the flow (Kraichnan, 1967; Boffetta & Ecke, 2012). However, many systems arising in geophysical and astrophysical fluid dynamics fail to be fully 3D because of the presence of strong restraints, although they are far from being 2D either. These restraints may arise from geometrical confinement (Smith et al., 1996; Celani et al., 2010; Benavides & Alexakis, 2017; Xia & Francois, 2017), rapid rotation (Smith & Waleffe, 1999; Pouquet & Marino, 2013; Campagne et al., 2014), strong stratification (Bartello, 1995; Smith & Waleffe, 2002; Oks et al., 2017) or the presence of strong magnetic fields (Alexakis, 2011; Favier et al., 2011). The detailed nature of the energy cascade in constrained 3D flows and its applications to geophysical fluid dynamics remain an open problem (Alexakis & Biferale, 2018).

Recent numerical simulations of rapidly rotating Rayleigh-Bénard convection, hereafter RRRBC, have shown that the turbulent state is susceptible to the evolution of large scale vortex (LSV) structures despite the presence of 3D fluctuations on all scales (Guervilly et al., 2014; Favier et al., 2014), in agreement with the predictions of an asymptotic description of the system valid in the limit of vanishingly small Rossby numbers (Julien et al., 2012; Rubio et al., 2014). Together these studies reveal the presence of an efficient energy transfer mechanism that extracts energy from small scale 3D fluctuations and deposits it in a box-scale barotropic (i.e. 2D) mode, bypassing the inverse energy cascade familiar from 2D dynamics. This process operates in parallel to energy transfer to large scale 2D modes by barotropic-barotropic interactions (Rubio et al., 2014) but dominates all aspects of the problem. Julien et al. (2012) and Julien et al. (2018) have suggested that the presence of the LSV introduces essential correlations among the phases of the small scale 3D fluctuations that facilitate direct energy extraction from these scales by the large scale mode, leading to a runaway that is only arrested by additional processes omitted from the simplest problem formulation. Such a runaway is characteristic of subcritical dynamics where it is triggered by finite amplitudes perturbations. The present paper is therefore devoted to a search for such subcritical behaviour in RRRBC. We do not address the question of whether geostrophic turbulence is itself linearly unstable to the generation of an LSV structure although our simulations suggest that at large enough rotation rates and large enough Rayleigh numbers this is in fact the case. We emphasize that phase correlations are missed in studies that focus on energy spectra alone, and that such correlations are inevitably absent from flows driven by a prescribed small scale force such as those studied by Chertkov et al. (2007) and Bouchet & Simonnet (2009). In these systems the only possible correlations are between the applied force and the resulting velocity field. In contrast, in the present 3D system the forcing of the large scale 2D flow can itself be dynamically affected through the action of the LSV on the small convective scales. The LSV observed in RRRBC may thus be a consequence of the proximity of the flow to 2D turbulence or due to the ability of the LSV to shape the correlations among the small scale fluctuations that appear to drive its formation.

In the present work we provide the first evidence for subcritical dynamics in turbulent RRRBC by demonstrating the coexistence of two numerically stable turbulent states at identical parameter values, one with an LSV structure and one without. Such bistability in turbulent flows is rare although it has also been found in rapidly rotating turbulence (Alexakis, 2015; Yokoyama & Takaoka, 2017), thin-layer turbulence (van Kan & Alexakis, 2018), Couette flows (Mujica & Lathrop, 2006; Zimmerman et al., 2011; Huisman et al., 2014; Xia et al., 2018) and von Kármán flows (Ravelet et al., 2004).

2 Mathematical formulation

2.1 Model and governing equations

We consider the evolution of a layer of incompressible fluid, bounded above and below by two impenetrable, fixed temperature, stress-free horizontal walls, a distance $h$ apart. The layer rotates about the $z$-axis, pointing vertically upwards, with a constant angular velocity $\bm{\Omega}=\Omega\bm{e}_{z}$ while gravity points downwards: $\bm{g}=-g\bm{e}_{z}$. The kinematic viscosity $\nu$ and thermal diffusivity $\kappa$ are assumed to be constant.

In the Boussinesq approximation, using the thermal diffusion time $h^{2}/\kappa$ as a unit of time and the depth $h$ of the layer as a unit of length, the dimensionless equations are

[TABLE]

[TABLE]

[TABLE]

where $\bm{u}\equiv\left(u,v,w\right)$ is the velocity, $p$ is the pressure and $\theta$ is the temperature fluctuation with respect to a linearly decreasing background temperature. The parameters are the Rayleigh number $Ra=\alpha g\Delta Th^{3}/(\nu\kappa)$, the Taylor number $Ta=4\Omega^{2}h^{4}/\nu^{2}$ and the Prandtl number $Pr=\nu/\kappa$. These dimensionless quantities involve $\alpha$, the coefficient of thermal expansion, and $\Delta T$, the imposed temperature difference between the two horizontal plates. For simplicity, and following earlier studies of the formation of large scales structures in this system, we take $Pr=1$. In the two horizontal directions, all variables are assumed to be periodic with the same spatial period $\lambda$ in both $x$ and $y$ directions. The boundary conditions at the upper ($z=1$) and lower ($z=0$) plates are $\partial_{z}u=\partial_{z}v=w=\theta=0$.

We solve equations (1)-(3) using the same mixed Fourier fourth-order finite difference scheme as used in Favier et al. (2014). We confirm the robustness of the present results by additionally running equivalent simulations using the fully spectral approach of Guervilly et al. (2014) and the open-source spectral-element code Nek5000 developed by Fischer et al. (2008) at the Argonne National Laboratory.

2.2 Finite amplitude initial conditions

In order to explore possible subcritical behaviour in turbulent RRRBC, we consider a particular set of initial conditions. Anticipating that the non-local upscale energy transfer eventually saturates by creating a barotropic vortex dipole at the box scale (Julien et al., 2012; Rubio et al., 2014), we consider a depth-invariant initial condition, given by

[TABLE]

and parametrised by the amplitude $A$. This initial condition corresponds to a symmetric vortex dipole at the box scale with a cyclone located at the box center $x=\lambda/2$ and $y=\lambda/2$ and an anticyclone located at the corners of the periodic domain. In addition, we add random perturbations of small amplitude $\pm 2.5\times 10^{-2}$ to the temperature field in order to initiate the Rayleigh-Bénard instability. The total kinetic energy density of this initial flow is

[TABLE]

where $V$ is the total volume. Note that $K_{0}$ is independent of the aspect ratio of the box. The purely viscous decay of this initial condition is given by

[TABLE]

2.3 Flow decomposition

We define the depth-averaged 2D horizontal flow, subsequently called the 2D mode, as

[TABLE]

where $u$ and $v$ are the velocity components in the $x$ and $y$ directions, respectively, and fast 3D fluctuations, subsequently called the 3D mode, as

[TABLE]

We can then define the (purely horizontal) kinetic energy density associated with the slow 2D mode as

[TABLE]

and the kinetic energy density associated with the fast 3D mode as

[TABLE]

3 Results

In this paper, we fix $Ta=10^{8}$ which is sufficient to sustain a non-local inverse cascade for $5\times 10^{6}\lesssim Ra\lesssim 2\times 10^{7}$ (Guervilly et al., 2014; Favier et al., 2014). For $Ra<5\times 10^{6}$ the flow is not turbulent enough while for $Ra>2\times 10^{7}$ the flow is insufficiently constrained by rotation and so not anisotropic enough, and only the traditional forward energy cascade is observed. In between these two limits, the flow is both turbulent and dynamically constrained by rotation, leading to the spontaneous emergence of a LSV from purely 3D perturbations driven by the Rayleigh-Bénard instability. In order to better understand the nature of the transition from a state with LSV to a state without, we consider here the particular case $Ra=3\times 10^{7}$, for which no systematic LSV were observed starting from random infinitesimal temperature perturbations (Guervilly et al., 2014; Favier et al., 2014), and focus on the behaviour with aspect ratio $\lambda=4$. This aspect ratio is chosen to ensure a clear scale separation between the convective eddies and the box size (see section 4 for a discussion of the effects of varying $\lambda$). For this set of control parameters, the flow is dominated by 3D small scale turbulent fluctuations (see left panels in figure 2), characterized by the Reynolds number $Re=\sqrt{\left<w^{2}\right>}\approx 535$ and the micro-Rossby number $Ro_{\omega}=\sqrt{\left<\omega_{z}^{2}\right>}/\sqrt{Ta}\approx 1.4$. Here $\omega_{z}$ is the vertical component of the vorticity and the brackets denote spatial and temporal averaging. Using this simulation as a reference, we ran a number of additional simulations, continuously varying the amplitude $A$ of the initial vortex dipole defined in Eq. (4).

Figure 1(a) shows the temporal evolution of the kinetic energy density $K_{2D}$ of the 2D flow as defined by Eq. (10). We observe a clear transition as the initial amplitude $A$ of the vortex dipole increases. For small amplitudes, typically $A\lesssim 800$, the kinetic energy of the 2D flow decreases until it reaches the equilibrium value corresponding to the reference case $A=0$ (i.e., no initial vortex). Note that this decay closely follows the purely viscous decay of the initial condition as given by Eq. (6), shown as dashed lines in figure 1(a), indicating that there is no significant energy transfer from the 3D fluctuations to the 2D flow. For larger amplitudes however, typically $A\gtrsim 800$, we observe an initial decay of $K_{2D}$ followed by an approximately linear increase until the energy eventually saturates at very long times. Note that close to the transition threshold, see case $A=875$ for example, it is not yet clear whether the vortex will grow or decay. In view of figure 1(a) which shows that the large scale vortex has not yet reached saturation, additional and longer simulations are required to determine any residual dependence of the saturated state on $A$ although we expect that all solutions involving a growing LSV will eventually saturate at the same amplitude.

These results clearly point towards subcritical behaviour in the transition between a weak subdominant 2D flow and a strong non-local inverse energy transfer efficiently feeding energy into the largest available spatial scale of the domain. Note that while the two states correspond to the same control parameters, the LSV state has a total kinetic energy density approximately 8 times that of the 3D state. The ratio between $K_{2D}$ and the total kinetic energy $K_{2D}+K_{3D}$ is shown in figure 1(b). It appears that the finite amplitude transition occurs once that $K_{2D}(t=0)\approx K_{3D}$ although it depends very subtly on the initial condition and can occur even after a long transient, as in the case $A=875$. A more accurate estimate of the critical amplitude $A$ is beyond the scope of this paper, however, since it would require running many realisations only changing the initial temperature perturbation (see section 4 below for a brief discussion concerning the nondeterministic nature of this transition). Visualizations of the saturated states for the same control parameters, but two different initial conditions, are shown in figure 2. Without the initial LSV, or when the amplitude $A$ is too small, the equilibrium state is dominated by a 3D flow at small scales while the 2D flow remains marginal. These small scale fluctuations are fully 3D, as expected from the moderate value of the Rossby number, $Ro=0.55$. In contrast, above the critical value of $A$, the LSV is continuously amplified while remaining at the box scale.

We now focus on the spectral statistics of the two different states. The kinetic energy spectra for each component of the flow, as defined by equations (7)-(9) (see also Guervilly et al. (2014) and Favier et al. (2014)), are shown in figure 3. These spectra are averaged over depth and time. With no initial vortex, the 3D flow is dominant at virtually all scales except for the smallest available wavenumber where most of the energy is contained in the 2D mode. This subdominant 2D flow is in equilibrium, meaning that there is a balance between viscous dissipation and baroclinic forcing (see below). There is no systematic growth of the 2D mode and no large scale condensate is reached. Above the critical initial amplitude, however, the LSV is able to extract energy efficiently from the small scale 3D flow leading to rapid growth of the 2D mode with horizontal wavenumber $k_{h}\leq 3$, while the 3D flow remains largely unchanged at all scales.

We examine next the energetics of the 2D depth-invariant flow. Starting from Eqs. (1)-(2), the governing equations for the purely horizontal 2D flow $\left<\bm{u}\right>_{z}$ are (see also Benavides & Alexakis (2017))

[TABLE]

where $\nabla_{h}$ is the horizontal gradient operator. The last term in Eq. (12) corresponds to the forcing of the 2D mode by the 3D fluctuations. Taking the scalar product of Eq. (12) with $\left<\bm{u}\right>_{z}$ and volume-averaging leads to the energy density equation of the 2D mode

[TABLE]

The first term on the right-hand side corresponds to the viscous dissipation rate of the 2D flow while the last term represents the 2D energy production from the 3D fluctuations. Looking at Eq. (13), it is clear that the 2D flow is in equilibrium only when the viscous dissipation $\mathcal{D}$ is balanced by the 3D forcing $\mathcal{F}$. It follows that growth of $K_{2D}$ from a given equilibrium state can only be achieved by reducing the dissipation or increasing the forcing. Figure 4(a) shows that for $A=0$, the forcing is not zero but is exactly balanced by viscous dissipation. The sum $\mathcal{D}+\mathcal{F}$ oscillates rapidly around zero, a fact consistent with the quasi-constant value of $K_{2D}$ observed for this case in figure 1(a). For $A=1600$, however, both dissipation and forcing increase in amplitude, and the sum is on average positive at least initially, corresponding to the growth of the 2D kinetic energy observed in figure 1(a). Note that the dissipation increases slowly with time while the baroclinic forcing term, while strongly fluctuating, remains quasi-constant (although it does grow very slightly). This increase in the dissipation is only observed for the 2D barotropic component; the dissipation associated with the 3D fluctuations remains largely unchanged (not shown). This is consistent with the condensation mechanism observed in purely 2D turbulence (Smith & Yakhot, 1994; Chertkov et al., 2007; Gallet & Young, 2014) whenever no large scale damping term is introduced to balance the inverse energy flux, and is a consequence of the slow growth of the dominant energy-containing scale.

At this stage, it is clear that the presence of an initially imposed vortex dipole somehow modifies the small scale convective flow, enhancing energy transfer into the 2D mode. One possible explanation for this behavior is that the 2D flow becomes significantly correlated with the 3D forcing, a correlation that can be quantified by looking at the angle $\gamma$ defined by $\cos\gamma=\left<\bm{u}\right>_{z}\cdot\left<\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}\right>_{z}/(|\left<\bm{u}\right>_{z}|\ |\left<\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}\right>_{z}|)$. However, our computations failed to reveal any significant differences in the probability density function of $\cos\gamma$ between cases with and without LSV (not shown), implying that the transfer enhancement cannot be explained by an increase in the correlation between the 3D forcing and the 2D flow. Figure 5 provides a clue. The figure reveals a clear imprint of the large scale 2D structure $\left<\bm{u}\right>_{z}$ in the vertically-averaged 3D fluctuations $\left<|\bm{u}^{\prime}|\right>_{z}$ and in the forcing term $\left<\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}\right>_{z}$: the fluctuations are locally suppressed, likely by the strong, relatively ordered large scale shear or vorticity associated with the vortex structure, implying enhanced correlations in the phases of the small scale field. This suppression of the small scale 3D fluctuations is consistent with the observed reduction in the Reynolds number (from $Re\approx 535$ for $A=0$ to $Re\approx 512$ for $A=1600$) and of the Nusselt number (from $Nu\approx 28.8$ for $A=0$ to $Nu\approx 27.6$), as already noted by Guervilly et al. (2014) for control parameters for which the LSV emerges spontaneously. It was also observed in thin-layer turbulence experiments (Xia et al., 2011) although no subcritical behaviour was reported in this study. This suppression in turn enhances interactions between two large $k_{h}$ 3D modes that transfer energy into a small $k_{h}$ 2D mode, bypassing the standard inverse cascade and enhancing the power spectrum of the forcing at these large scales, as shown in figure 5 (right panel). Note that, irrespective of the amplitude $A$, the forcing always peaks at $k_{h}\approx 30$, which is approximately the Taylor microscale characteristic of the small scale vorticity field. The 2D energy is, however, very small at these scales (figure 3) and when $A=0$ this small scale forcing is in equilibrium with the dissipation (figure 4(a)). For $A=1600$, however, the suppression of the small scale fluctuations leads directly to enhanced forcing of the 2D flow at the box scale, allowing for a runaway growth. The non-local nature of this upscale energy transfer has already been discussed in previous studies (Rubio et al., 2014; Favier et al., 2014), and we expect similar non-local energy fluxes in our subcritical state, although this remains to be fully explored in future studies.

The resulting positive feedback mechanism differs in an important respect from that present in 2D turbulence driven by a small scale prescribed stochastic force: in the present case the fluctuating force is itself dynamically modified by the LSV, and not just its rate of working as in Gallet & Young (2014). In addition, the positive feedback revealed in figure 5 provides an indication that the LSV state will not in general spontaneously jump back to the lower LSV-free state: as soon as a transition starts to take place and the LSV is observed to grow, its back-reaction on the 3D fluctuations favours energy transfer into the 2D component, leading to run-away dynamics which can only be arrested by viscous condensation at the box scale or by other large scale effects not included in our simple model.

4 Discussion

We now discuss how the subcritical behaviour identified here depends on the various control parameters. Thus far we have focused on aspect ratio $\lambda=4$, but similar simulations were performed with $\lambda=2$ (Table 1). We observe very similar results: the initial LSV is amplified only above a finite critical amplitude. For this value of the aspect ratio, however, the existence of a non-local inverse energy transfer is not systematic. When some of the simulations were repeated with different random initial temperature perturbations, some cases exhibited an inverse cascade while others did not (figure 6(a)) in a manner reminiscent of pipe flow (Darbyshire & Mullin, 1995). This is a consequence of the stochastic nature of the forcing term arising from the 3D fluctuations (figure 4) which can drastically affect the properties of the transition (Fauve et al., 2017). Indeed the transition seems much less robust for $\lambda=2$ than for $\lambda=4$, an effect we ascribe to an increase in the amplitude of the fluctuations that arises from the smaller domain size, indicating that a reasonable scale separation between the initial vortex and the small scale flow is required for a robust and reproducible transition. We confirmed this observation with simulations at $\lambda=6$, for which a case with $A=800$ decays while a case with $A=1200$ eventually grows (although with the available computational resources neither of these cases can be run until saturation owing to the presence of much longer transients, see figure 6(b)). This confirms that the observed transition is robust with respect to change in the aspect ratio of the vortex, although additional simulations (for example, decoupling the aspect ratio of the numerical domain and the initial size of the LSV) are required for a definitive conclusion about its effect on the threshold.

The results discussed above only apply for $Ta=10^{8}$ and $Ra=3\times 10^{7}$. We chose $Ta=10^{8}$ for numerical convenience, as increasing the Taylor number would be numerically prohibitive for the large aspect ratio domains and long-time dynamics considered here. Interestingly, the range of Rayleigh numbers for which a LSV is known to spontaneously emerge increases with the Taylor number (Guervilly et al., 2014; Favier et al., 2014). It is reasonable, therefore, to assume that subcritical transitions will also be observed at higher Taylor numbers, and might even be more prominent. However, this hypothesis remains to be confirmed.

We have also explored the effect of varying the Rayleigh number. First, it is known that a LSV cannot be sustained when the convection is insufficiently turbulent. For $Ta=10^{8}$, $Ra\geq 5\times 10^{6}$ is required: when we repeated our simulations for $Ra=4\times 10^{6}$, varying the amplitude $A$ of the vortex, we found no subcritical behaviour and all simulations eventually converged towards the same equilibrium dominated by 3D fluctuations. It appears that the subcritical transition is only present for sufficiently turbulent flows, well above the threshold of the linear instability. Finally, the same experiments were repeated for higher Rayleigh numbers (namely $Ra=4\times 10^{7}$ and $Ra=5\times 10^{7}$) with results that are more subtle: for $\lambda=2$, we did not observe any subcriticality, obtaining a 3D state irrespective of $A$. Surprisingly, for $\lambda=4$ and $Ra=4\times 10^{7}$, we recovered a LSV state that was not present for $\lambda=2$. This fact provides further evidence that a clear scale separation between the vortex and the small scale 3D flow favors a subcritical transition. A detailed study of the scaling of the critical amplitude $A_{c}$ with $Ta$, $Ra$ and $\lambda$ is, however, beyond the scope of this paper, although it represents an interesting line of investigation not only from a fundamental point of view, but also for possible applications to LSV structures in geophysical and astrophysical flows.

Our choice of initial conditions (4) is arbitrary and other options could be explored. In our moderate Rossby number simulations, symmetry breaking favouring cyclonic motions has been observed, in contrast with the symmetry between cyclonic and anticyclonic vortices found in the limit $Ro\to 0$ (Julien et al., 2012). Our initial condition is, however, perfectly symmetric requiring a significant transient phase to break up the large but unstable anticyclone, something that could be avoided if the calculations were to be initialized with a cyclonic structure only. In this respect, the reduced description of Julien et al. (2012) may be useful for exploring the detailed mechanism behind the subcritical transition discovered in this paper and in particular the presence of hysteresis in $Ra$, provided, of course, that the reduced equations exhibit similar behaviour.

We believe that the mechanism responsible for the observed subcritical transition is not specific to rotating convection and is likely to occur in other systems with multiple cascade scenarios, including thin-layer and magnetohydrodynamic turbulence. In particular, we emphasize that rotation is not required to observe LSV in 3D flows. Small scale anisotropy, and its possible enhancement by large scale flows, is the key and is present in all turbulent systems with multiple cascade scenarios, from thin-layer to magnetohydrodynamic turbulence. We also emphasize that the transition identified here separates two fully turbulent states, in contrast to the classical subcritical transition from laminar to turbulent shear flow in, e.g., pipe flow (Darbyshire & Mullin, 1995; Eckhardt et al., 2007). Some concepts developed in this field may nevertheless prove useful in the present context, such as the search for optimal perturbations (Kerswell, 2018).

Acknowledgments

This work was initiated during the workshop “Rotating Convection: from the Lab to the Stars” organized and supported by the Lorentz Center at the University of Leiden (http://www.lorentzcenter.nl/). This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program “Investissements d’Avenir” supervised by the Agence Nationale de la Recherche. Computations were also conducted with the support of the HPC resources of GENCI-IDRIS (Grants No. A0020407543 and A0040407543), the Rocket High Performance Computing service at Newcastle University and the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). C.G. was supported by the UK Natural Environment Research Council under grant NE/M017893/1. B.F. acknowledges funding by the European Research Council under the European Union’s Horizon 2020 research and innovation program through Grant No. 681835-FLUDYCO-ERC-2015-CoG.