2. Problem formulation

We consider two-dimensional Rayleigh–Bénard convection governed by the dimensionless equations

(2.1a)$$\begin{gather} \nabla^{2}\psi_t + \{ \psi, \nabla^{2}\psi\} = {Ra}Pr\theta_x + Pr\nabla^4 \psi, \end{gather}$$

(2.1b)$$\begin{gather}\theta_t + \{ \psi, \theta\} = \psi_x + \nabla^2 \theta, \end{gather}$$

where $\psi (x,y,t)$ is the streamfunction and $\theta (x,y,t)$ is the field of temperature fluctuations from the heat-conducting temperature profile $T = 1 - y + \theta$. Here, subscripts denote partial derivatives, and $\{f, g\} = f_x g_y - g_x f_y$ is the usual Poisson bracket. Equations (2.1) have been non-dimensionalised using $d$, $d^2/\kappa$ and ${\rm \Delta} T$ as the relevant scales for length, time and temperature, respectively, where $d$ is the height of the fluid layer, $\kappa$ is the thermal diffusivity and ${\rm \Delta} T$ is the temperature difference imposed across the layer. The two dimensionless parameters in the system (2.1) are the Prandtl and Rayleigh numbers:

(2.2a,b)\begin{equation} Pr = \frac{\nu}{\kappa},\quad Ra = \frac{\alpha {\rm \Delta} T g d^3}{\nu \kappa}, \end{equation}

where $\nu$ is the kinematic viscosity, $\alpha$ is the thermal expansion coefficient and $g$ is the gravitational acceleration.

The dimensionless spatial domain of our problem is $(x, y) \in [0, \varGamma ] \times [0, 1]$, where $\varGamma$ is the width-to-height aspect ratio. The aspect ratio in this study is fixed at $\varGamma = 2$. At the lower ($y = 0$) and upper ($y = 1$) boundaries the temperature satisfies isothermal conditions while the velocity field satisfies no-penetration and stress-free boundary conditions, i.e.

(2.3a,b)\begin{equation} \theta |_{y=0,1} = 0, \quad \psi |_{y=0,1} = \psi_{yy} |_{y=0,1} = 0,\end{equation}

and periodic boundary conditions are imposed in the $x$ direction.

Given the above boundary conditions, it is convenient to decompose the streamfunction into basis functions with Fourier modes in the $x$ direction and sine modes in the $y$ direction:

(2.4)\begin{align} \psi(x,y,t) &= \sum_{k_x\in\mathbb{Z}} \sum_{k_y\in\mathbb{Z}_{{>}0}} \hat{\psi}_{k_x,k_y}(t)\, \mathrm{e}^{2{\rm \pi} \,\mathrm{i} \,k_x x/\varGamma} \sin{({\rm \pi} k_y y)} \nonumber\\ &= \sum_{k_x,k_y}\hat{\psi}_{k_x,k_y}(t) \phi_{k_x,k_y}(x,y), \end{align}

where $\hat { \psi }_{k_x,k_y}$ is the amplitude of the ($k_x, k_y$) mode of $\psi$, and we have the constraint $\hat { \psi }_{k_x,k_y} = \hat { \psi }^*_{-k_x,k_y}$ (with the asterisk denoting complex conjugation). We decompose $\theta$ in the same way.

Rayleigh (Reference Rayleigh1916) has shown that the static conduction state ($\psi = \theta = 0$) bifurcates supercritically to a steady pattern of counter-rotating convection rolls vertically spanning the layer $(k_y = 1)$ when $Ra > Ra_c$ with

(2.5)\begin{equation} Ra_c = \min_{k_x}{\left(\frac{{\rm \pi}^4}{4\varGamma^4} \frac{(4k_x^2+\varGamma^2)^3}{k_x^2} \right)}. \end{equation}

The minimum is taken over integer wavenumbers $k_x$, and occurs at $k_x = 1$ provided that $\varGamma < 2^{4/3}\sqrt {1+2^{2/3}} \approx 4.05$. Consequently, for our set-up with $\varGamma =2$, the modes $\hat \psi _{1,1}$ and $\hat \theta _{1,1}$ from (2.4) become excited through a supercritical pitchfork bifurcation at $Ra_c = 8{\rm \pi} ^4$. In the remainder of the article, we refer to the resulting steady pattern as the steady convection roll state (SCRS). An instance of the SCRS with the parameters $(Pr,Ra) = (1,10^5)$ is shown in figure 1. A further decomposition which will aid our analysis is to partition the modes into odd modes indicated by $O$ with $k_x + k_y \in 2\mathbb {Z} + 1$, and even modes indicated by $E$ with $k_x + k_y \in 2\mathbb {Z}$ similarly to previous studies (Chandra & Verma Reference Chandra and Verma2011; Verma, Ambhire & Pandey Reference Verma, Ambhire and Pandey2015). Consistent with this notion of odd and even modes, we decompose $\psi$ as

(2.6a)$$\begin{gather} \psi^O(x,y,t) = \tfrac{1}{2}[\psi(x,y,t) + \psi(x+\varGamma/2,1-y,t)], \end{gather}$$

(2.6b)$$\begin{gather}\psi^E(x,y,t) = \tfrac{1}{2}[\psi(x,y,t) - \psi(x+\varGamma/2,1-y,t)], \end{gather}$$

and similarly for $\theta ^O$ and $\theta ^E$. Thus, $\psi ^O$ and $\theta ^O$ consist only of odd modes as defined above, while $\psi ^E$ and $\theta ^E$ consist only of even modes. Applying the decomposition (2.6) to (2.1), we find that

(2.7a)$$\begin{gather} \nabla^{2}\psi^O_t + \{ \psi^E, \nabla^{2}\psi^O\} + \{ \psi^O, \nabla^{2}\psi^E\} = {Ra}Pr\theta^O_x + Pr\nabla^4 \psi^O, \end{gather}$$

(2.7b)$$\begin{gather}\theta^O_t + \{ \psi^E, \theta^O\} + \{ \psi^O, \theta^E\}= \psi^O_x + \nabla^2 \theta^O \end{gather}$$

and

(2.8a)$$\begin{gather} \nabla^{2}\psi^E_t + \{ \psi^E, \nabla^{2}\psi^E\} - {Ra}Pr\theta^E_x - Pr\nabla^4 \psi^E ={-} \{ \psi^O, \nabla^{2}\psi^O\}, \end{gather}$$

(2.8b)$$\begin{gather}\theta^E_t +\{ \psi^E, \theta^E\} - \psi^E_x - \nabla^2 \theta^E ={-}\{ \psi^O, \theta^O\}. \end{gather}$$

We observe that $( \psi ^O, \theta ^O )$ satisfy the linear partial differential equations (2.7), with coefficients that depend on $( \psi ^E, \theta ^E )$, while $( \psi ^E, \theta ^E )$ satisfy the autonomous nonlinear partial differential equations (2.8), with forcing terms that depend on $( \psi ^O, \theta ^O )$.

Figure 1.

The SCRS calculated with the parameters $(Pr,Ra) = (1,10^5)$ using the numerical scheme outlined in Appendix A.

At the onset of steady convection, the excited modes $\hat \psi _{1,1}$ and $\hat \theta _{1,1}$ are even modes, and indeed the SCRS has the property that $\psi ^O \equiv \theta ^O \equiv 0$, which remains true as we increase $Ra$. In this case, the system of (2.7) and (2.8) reduces to

(2.9a)$$\begin{gather} \nabla^{2}\psi^E_t + \{ \psi^E, \nabla^{2}\psi^E\} = {Ra}Pr\theta^E_x + Pr\nabla^4 \psi^E, \end{gather}$$

(2.9b)$$\begin{gather}\theta^E_t + \{ \psi^E, \theta^E\} = \psi^E_x + \nabla^2 \theta^E, \end{gather}$$

where the even modes are decoupled from the odd modes. The SCRS has a reflection symmetry about the vertical line which separates the counter-rotating convection rolls. Without loss of generality, we choose this vertical line to be at $x=0$. Thus, the steady-state solution $(\psi ^E, \theta ^E)$ of the system (2.9) is invariant under the symmetry

(2.10a,b)\begin{equation} \psi^E(x,y) \mapsto -\psi^E({-}x,y), \quad \theta^E(x,y) \mapsto \theta^E({-}x,y). \end{equation}

Since the zonal modes $\hat \psi _{0,k_y}$ are $x$-independent, (2.10a,b) implies that $\hat \psi ^E_{0,k_y}\equiv 0$. In other words, the SCRS has no zonal flow, and zonal modes can emerge only once the SCRS has become unstable.

3. Linear stability analysis

In this section, we present the linear stability analysis of the SCRS. For given values of $Pr$ and $Ra$, we first solve numerically the steady-state problem of (2.9), following the procedure described in Appendix A. We then exploit the partial decoupling of odd and even modes to consider odd and even perturbations separately. We analyse odd perturbations by setting

(3.1)\begin{align} \psi(x,y,t) &= \psi^E(x,y) + \delta\mathrm{e}^{\sigma^o t} \psi^o(x,y) \nonumber\\ &= \psi^E(x,y) + \delta\mathrm{e}^{\sigma^o t}\sum_{k_x+k_y\in2\mathbb{Z}+1} \hat{\psi}^o_{k_x,k_y}\phi_{k_x,k_y}(x,y), \end{align}

(3.2)\begin{align} \theta(x,y,t) &= \theta^E(x,y) + \delta\mathrm{e}^{\sigma^o t}\theta^o(x,y) \nonumber\\ &=\theta^E(x,y) + \delta\mathrm{e}^{\sigma^o t}\sum_{k_x+k_y\in2\mathbb{Z}+1} \hat{\theta}^o_{k_x,k_y}\phi_{k_x,k_y}(x,y), \end{align}

where $(\psi ^E,\theta ^E)$ is the SCRS, $\delta \ll 1$ and $\psi ^o$, $\theta ^o$ consist only of odd modes. By substituting (3.1) into (2.1) and linearising with respect to $\delta$, we obtain an eigenvalue problem for the odd growth rates $\sigma ^o$ (see Appendix B for details). Then, we solve this eigenvalue problem numerically and determine that the SCRS is unstable with respect to odd perturbations if $\max {[\mathcal {R}(\sigma ^o)]} > 0$. An analogous approach is used to determine the corresponding even growth rates $\sigma ^e$ and thus the stability of the SCRS with respect to even perturbations.

Figure 2 displays the stability of the SCRS in the $(Pr, Ra)$ plane, highlighting the stable ($S$) regime in white and the unstable ($US$) regime in blue. The pink region denotes $Ra < Ra_c$ where the static conduction state is stable. Although the SCRS does become unstable to even perturbations within the region $US$, we always observe that $\max {[\mathcal {R}(\sigma ^o)]}>\max {[\mathcal {R}(\sigma ^e)]}$ and the dominant unstable perturbation therefore consists of odd modes. Moreover, for $Pr\lesssim 7\times 10^5$ the most unstable eigenfunction has the property that $\hat {\psi }^o_{k_x,k_y}$ and $\mathrm {i}\hat {\theta }^o_{k_x,k_y}$ are purely real, so that $\psi ^o$ and $\theta ^o$ break the reflection symmetry (2.10a,b). Zonal modes $\hat \psi _{0,k_y}$ with $k_y\in 2\mathbb {Z}+1$ are present in the most unstable perturbation until we reach very large Prandtl numbers (see § 5). In particular, we find that the largest scale zonal mode $\hat \psi _{0,1}$ is present in the most unstable eigenfunction at the stability boundary for $Pr\lesssim 7\times 10^5$, and is the dominant mode for $Pr\lesssim 4.79$ and $8.58\lesssim Pr\lesssim 40$.

Figure 2.

The $(Pr, Ra)$ plane highlighting when the SCRS is linearly stable (white, $S$) and unstable (blue, $US$). Only the static conduction state is stable in the red shaded region below $Ra = Ra_c$ (dashed line). When $Pr \ll 1$, the SCRS becomes unstable at $0< Ra-Ra_c \ll 1$. The highlighted points in red are at $p_1 = (5.48 \times 10^{-2},1.16 \times 10^3)$, $p_2 = (1.81 \times 10^{-2},4.12 \times 10^3)$, $p_3 = (0.2175,2.59 \times 10^5)$, $q_1 = (0.2,3.32 \times 10^5)$, $q_2 = (4.97,3.949 \times 10^6)$, $q_3 = (8.58,1.27 \times 10^6)$, $p_4 = (6.16,5.43 \times 10^4)$, $p_5 = (9.53, 3.24 \times 10^4)$ and $q_4 = (7 \times 10^5,2.54 \times 10^7)$. Visualisations of the SCRS at all the points $p_i$ and $q_i$ are given as supplementary figures available at https://doi.org/10.1017/jfm.2022.185.

The points $p_i$ annotated in figure 2 highlight turning points in the stability boundary. The behaviour of $\sigma ^o$ near each such point is plotted in figure 3. Figure 3(c) shows that, close to $Pr=0.2$ and with increasing Rayleigh number, a real eigenvalue crosses the origin twice, and the SCRS first loses and then regains stability. The steady state unexpectedly regains stability as the Rayleigh number increases at point $p_1$ with $Pr\approx 0.055$. Similar behaviour has been observed by Paul et al. (Reference Paul, Verma, Wahi, Reddy and Kumar2012), with $Pr = 6.8$ and $\varGamma = 2\sqrt {2}$. At point $p_2$, the opposite occurs with the SRCS first regaining and then losing stability as $Ra$ increases. At points $p_4$ and $p_5$ we instead find complex eigenvalues crossing the imaginary axis twice as $Ra$ increases with fixed $Pr$ or vice versa.

Figure 3.

The largest odd growth rate $\sigma ^o$, plotted with one dimensionless parameter varied and the other fixed. The plots show the behaviour close to each of the points $p_i$ highlighted in figure 2. In (a–c) $\sigma ^o$ is real, while in (d,e) the motion of $\sigma ^o$ in the complex plane is tracked, with the arrows indicating the direction of increasing $Ra$ or $Pr$.

At each of the points labelled $q_i$, the stability boundary is not smooth, and there is a qualitative change in the eigenfunction to which the SCRS becomes unstable. Figure 4(a) shows how, with $Pr=0.2$, the SCRS regains stability with a real eigenvalue crossing zero at $Ra \approx 3.32859 \times 10^5$. With further increase in $Ra$, the eigenvalue splits into a complex conjugate pair, which then recrosses the imaginary axis as the SCRS loses stability in a Hopf bifurcation. The cusp at $q_1$ occurs when these three events happen simultaneously. These observations have been confirmed with direct numerical simulation (DNS) of the system (2.7) and (2.8), using the pseudospectral scheme described by Winchester et al. (Reference Winchester, Dallas and Howell2021). Figure 5 shows numerical results for the largest scale odd mode $\hat {\psi }_{0,1}$ obtained with the right-hand side of (2.8) set to $0$, so that we have persistent exponential growth or decay in the odd modes. We observe that, with increasing $Ra$, the odd perturbations regain stability, then become oscillatory and finally become unstable again.

Figure 4.

The largest odd growth rate(s) $\sigma ^o$, plotted with one dimensionless parameter varied and the other fixed. The plots show the behaviour close to each of the points $q_i$ highlighted in figure 2. The motion of $\sigma ^o$ in the complex plane is tracked, with the arrows indicating the direction of increasing $Ra$ or $Pr$.

Figure 5.

Time series of the $\hat \psi _{0,1}$ mode from the linearised DNS with the right-hand side of (2.8) set to zero. Here we fix $Pr = 0.2$ ($q_1$) and have $Ra = 3.32 \times 10^5$ (a), $Ra = 3.32859 \times 10^5$ (b), $Ra = 3.33 \times 10^5$ (c) and $Ra = 3.35 \times 10^5$ (d). The initial conditions are the SCRS with a small perturbation.

Figure 4(c) shows that, with $Pr=8.58$, the SCRS regains stability through a Hopf bifurcation at $Ra = 1.28 \times 10^6$ but loses stability shortly after at $Ra = 1.29 \times 10^6$ with much larger oscillation frequency. The corner $q_3$ occurs when these two pairs of eigenvalues cross the imaginary axis simultaneously. Again, this observation has been verified with DNS close to $q_3$. In figure 6 we fix $Pr = 8.57$ and see that the odd mode $\hat \psi _{0,1}$ grows at $Ra = 1.2 \times 10^6$, decays at $Ra = 1.27 \times 10^6$ and grows again at $Ra = 1.3 \times 10^6$, but now with a much greater oscillation frequency. As can been seen in figure 4(b), similar behaviour occurs at corner $q_2$. The final corner point $q_4$, which lies in the large-$Pr$ region, is considered in more detail in § 5. How the dominant eigenfunction changes at each of the points $q_i$ is shown in the supplementary material.

Figure 6.

Time series of the $\hat \psi _{0,1}$ mode close to $q_3$ from the linearised DNS with the right-hand side of (2.8) set to zero. We fix $Pr = 8.57$ and the initial conditions are the SCRS with a small perturbation. Longer time series clearly showing the exponential growth and decay are in the supplementary material.

Figure 2 hints that the Rayleigh number at which the SCRS loses stability approaches $Ra_c$ as $Pr \to 0$. In figure 7 we plot

(3.3)\begin{equation} \delta Ra:= Ra - Ra_c \end{equation}

versus $Pr$ and indeed see that the stability boundary follows a clear power law with $\delta Ra \propto Pr^2$. When this boundary is crossed, at small $Pr$ the SCRS becomes unstable through a pitchfork bifurcation as a new steady state containing both even and odd modes emerges.

Figure 7.

The $(Pr, \delta Ra)$ plane for small $Pr$ highlighting when the SCRS is linearly stable (white, $S$) and unstable (blue, $US$). On the vertical axis we plot $\delta Ra = Ra - Ra_c$, where $Ra_c = 8{\rm \pi} ^4$. The red dashed line highlights the power law $\delta Ra \propto Pr^2$, where the prefactor has been calculated by fitting.

This prediction has been verified with fully nonlinear DNS solving (2.1) using the pseudospectral scheme described in Winchester et al. (Reference Winchester, Dallas and Howell2021). In figure 8(a) we plot a bifurcation diagram showing the amplitudes of the even mode $\hat \psi _{1,1}$ and the odd mode $\hat \psi _{0,1}$ versus $\delta Ra$ with $Pr=0.01$. The SCRS emerges at $\delta Ra = 0$ and becomes unstable at $\delta Ra = \delta Ra_c \approx 6.135$ as the new steady state emerges. In figure 8(b) we plot the time series of $\lvert \hat \psi _{1,1} \rvert$ and $\lvert \hat \psi _{0,1} \rvert$ at $Pr=0.01$ and $\delta Ra = 10.7$. At this point the mixed steady state has become unstable and the mode $\lvert \hat \psi _{1,1} \rvert$ undergoes nonlinear relaxation oscillations whilst $\lvert \hat \psi _{0,1} \rvert$ exhibits bursts that quickly decay away. At $\delta Ra = 2.22 \times 10^3$, the relaxation oscillations and the bursts persist, but are separated by a much slower buildup phase, as shown in figure 8(c).

Figure 8.

(a) Bifurcation diagram plotting $\lvert \hat \psi _{1,1} \rvert$ (blue) and $\lvert \hat \psi _{0,1} \rvert$ (red) against $\delta Ra$ with $Pr = 10^{-2}$. Corresponding time series of $\lvert \hat \psi _{1,1} \rvert$ and $\lvert \hat \psi _{0,1} \rvert$ when (b) $\delta Ra = 10.7$ and (c) $\delta Ra = 2.2 \times 10^3$.

In summary, the dynamics at low Prandtl number with $\delta Ra >\delta Ra _c$ can be partitioned into three phases, as follows.

1. (i) In the first phase we are close to the static conduction state with $\lvert \hat \psi _{1,1} \rvert,\lvert \hat \psi _{0,1} \rvert \approx 0$. However, when $\delta Ra > 0$, the static conduction state is unstable, so that $\hat \psi _{1,1}$ and other even modes grow exponentially.

2. (ii) The even modes grow until we are in the SCRS. Since $\delta Ra > \delta Ra_c$, the SCRS is unstable to odd perturbations, so that $\hat \psi _{0,1}$ and other odd modes grow exponentially.

3. (iii) Once $\hat \psi _{1,1}$ and $\hat \psi _{0,1}$ are of comparable amplitude, the system quickly collapses back to the static conduction state and the cycle continues.

This dynamics differs markedly from the chaotic behaviour found close to $Ra_c$ in small-$Pr$ and zero-$Pr$ convection in three dimensions (Thual Reference Thual1992; Pal et al. Reference Pal, Wahi, Paul, Verma, Kumar and Mishra2009). In the following section, we carry out a weakly nonlinear analysis to help us explain the observed power law (see figure 7) as well as the dynamics of the nonlinear relaxation oscillations and of the bursts (see figure 8).

4. Weakly nonlinear analysis at small Prandtl number

We perform a weakly nonlinear analysis with the Prandtl number being our small parameter, $\epsilon = Pr \ll 1$. The analysis is carried out for arbitrary aspect ratio $\varGamma <2^{4/3}\sqrt {1+2^{2/3}}$, and the other parameters and variables are scaled as follows:

(4.1a)$$\begin{gather} Ra = Ra_c+r\epsilon^2= {\rm \pi}^4\frac{(4+\varGamma^2)^3}{4\varGamma^4}+r\epsilon^2, \end{gather}$$

(4.1b)$$\begin{gather}t= \frac{\tau_{_O}}{\epsilon}= \frac{\tau_{_E}}{\epsilon^3} , \end{gather}$$

(4.1c)$$\begin{gather}\psi^E = \epsilon\psi_1^E(\tau_E, x,y) + \epsilon^2 \psi_2^E(\tau_E, x,y) + {O}(\epsilon^3), \end{gather}$$

(4.1d)$$\begin{gather}\theta^E = \epsilon\theta_1^E(\tau_E, x,y) + \epsilon^2 \theta^E(\tau_E, x,y) + {O}(\epsilon^3), \end{gather}$$

(4.1e)$$\begin{gather}\psi^O = \epsilon^2 \psi_1^O(\tau_O, x,y) + \epsilon^3 \psi_2^O(\tau_O, x,y) + {O}(\epsilon^4), \end{gather}$$

(4.1f)$$\begin{gather}\theta^O = \epsilon^2 \theta_1^O(\tau_O, x,y) + \epsilon^3 \theta_2^O(\tau_O, x,y) + {O}(\epsilon^4). \end{gather}$$

We have introduced two time scales $\tau _{_O}$ and $\tau _{_E}$, due to the distinct time-scale separation between the odd and even modes, observed for example in figure 8(c). Equation (4.1a) is inspired by the power law $\delta Ra \propto Pr^2$. The scalings of the dependent variables are constructed such that weak nonlinearity and coupling between even and odd modes enter at the same order, as we will see below.

First we examine the evolution of the even modes. At ${O}(\epsilon )$ and ${O}(\epsilon ^2)$ we find

(4.2a)$$\begin{gather} \psi^E_1 = \frac{4+\varGamma^2}{\varGamma}{\rm \pi} A \sin{(2{\rm \pi} x/\varGamma)} \sin{({\rm \pi} y)}, \end{gather}$$

(4.2b)$$\begin{gather}\theta^E_1 = 2 A \cos{(2{\rm \pi} x/\varGamma)} \sin{({\rm \pi} y)} \end{gather}$$

and

(4.3a)$$\begin{gather} \psi^E_2 = 0, \end{gather}$$

(4.3b)$$\begin{gather}\theta^E_2 ={-}{\rm \pi} A^2 \frac{4+\varGamma^2}{2\varGamma^2}\sin{(2 {\rm \pi}y)}, \end{gather}$$

which describe the asymptotic behaviour of the SCRS as $\delta Ra \to 0$.

At ${O}(\epsilon ^3)$, we have

(4.4a)$$\begin{gather} \nabla^{2}\psi^E_{1\tau_{_E}} + \{ \psi_1^O, \nabla^{2}\psi^O_1\} = Ra_c\theta^E_{3x} + \nabla^4 \psi^E_3 + r\theta^E_{1x}, \end{gather}$$

(4.4b)$$\begin{gather}\{ \psi_1^E,\theta^E_2\}= \psi^E_{3x} + \nabla^2 \theta^E_3, \end{gather}$$

the solvability condition for which gives us

(4.5)\begin{equation} \frac{(4+\varGamma^2)^2}{2\varGamma^3}{\rm \pi}^3\, \frac{\mathrm{d}A}{\mathrm{d}\tau_{E}} - r\frac{2{\rm \pi}}{\varGamma} A + \frac{(4+\varGamma^2)^4}{4\varGamma^7}{\rm \pi}^7 A^3 = \frac{2\mathrm{i}}{\varGamma}\langle \{ \psi_1^O, \nabla^{2}\psi^O_1\}, \phi_{1,1} \rangle. \end{equation}

The right-hand side of (4.5) is the inner product of the nonlinear forcing term with the $(1,1)$ basis function, and captures the net effect of the odd modes on the amplitude $A$ of the dominant even mode. If this term is set to zero, then (4.5) reduces to the standard Landau equation governing the amplitude $A$ of weakly nonlinear perturbations to the static conduction state (Fowler Reference Fowler1997).

To evaluate the right-hand side of (4.5), we now turn to the odd modes. At lowest order we obtain the problem

(4.6a)$$\begin{gather} \nabla^{2}\psi^O_{1\tau_{_O}} + \{ \psi^O_1, \nabla^{2}\psi^E_1\} + \{ \psi^E_1, \nabla^{2}\psi^O_1\} = Ra_c\theta^O_{1x} + \nabla^4 \psi^O_1, \end{gather}$$

(4.6b)$$\begin{gather}0= \psi^O_{1x} + \nabla^2 \theta^O_1, \end{gather}$$

which can be reduced to

(4.7)\begin{align} &\nabla^{2}\psi^O_{1\tau_{_O}} + \frac{{\rm \pi}(4+\varGamma^2)}{\varGamma}A\left\{ \sin\left(\frac{2{\rm \pi} x}{\varGamma}\right)\sin({\rm \pi} y), \nabla^2\psi_1^O+\frac{{\rm \pi}^2(4+\varGamma^2)}{\varGamma^2}\psi_1^O \right\} \nonumber\\ &\quad = \nabla^4 \psi^O_1 - Ra_c\nabla^{{-}2}\psi^O_{1xx}. \end{align}

Using the Fourier decomposition (2.4), i.e.

(4.8)\begin{equation} \psi^O_1 = \sum_{k_x+k_y \in 2\mathbb{Z}+1} \hat{\psi}_{k_x,k_y}\phi_{k_x,k_y}, \end{equation}

we can express (4.7) as a linear dynamical system of the form

(4.9)\begin{equation} \frac{\mathrm{d}\boldsymbol{\psi}^O_1}{\mathrm{d}\tau_{O}} = (A\boldsymbol{\mathsf{M}}+\boldsymbol{\mathsf{D}})\boldsymbol{\psi}^O_1. \end{equation}

Here $\boldsymbol{\mathsf{M}}$ and $\boldsymbol{\mathsf{D}}$ are known constant matrices, calculated as described in Appendix C, and $\boldsymbol {\psi }^O_1=(\hat {\psi }_{k_x,k_y})_{k_x+k_y\in 2\mathbb {Z}+1}$ is the vector of odd modes.

We recall that $A$ evolves on the much slower time scale $\tau _E=\epsilon ^2\tau _O$. As justified in Appendix C, to leading order in $\epsilon$, it follows that $\boldsymbol {\psi }^O_1$ is proportional to the most unstable eigenvector of the problem (4.9). We can therefore write

(4.10)\begin{equation} \boldsymbol{\psi}^O_1 (\tau_O;A)\sim B(\tau_O)\boldsymbol{v}(A), \end{equation}

where

(4.11)\begin{equation} (A\boldsymbol{\mathsf{M}}+\boldsymbol{\mathsf{D}})\boldsymbol{v}(A)=\sigma(A)\boldsymbol{v}(A), \end{equation}

with $\sigma (A)$ the eigenvalue of (4.11) with the largest real part, and $\boldsymbol {v}$ normalised such that $\|\boldsymbol {v}\|=1$. At leading order in $\epsilon$, the net amplitude $B$ of the odd modes thus evolves according to

(4.12)\begin{equation} \frac{\mathrm{d}B}{\mathrm{d}\tau_{O}} = \sigma(A)B. \end{equation}

With $\psi ^O_1$ decomposed as in (4.10), we can express the right-hand side of (4.5) in the form

(4.13)\begin{equation} \frac{2\mathrm{i}}{\varGamma}\langle \{ \psi_1^O, \nabla^{2}\psi^O_1\}, \phi_{1,1} \rangle={-}B^2\lambda(A), \end{equation}

where

(4.14)\begin{equation} \lambda(A)={-}\frac{2\mathrm{i}}{\varGamma} \langle\{v(A),\nabla^2v(A)\},\phi_{1,1}\rangle. \end{equation}

For a given value of $\varGamma$, the functions $\sigma (A)$ and $\lambda (A)$ can both be computed once and for all from (4.11) and (4.14). Equations (4.5) and (4.12) then provide a closed two-dimensional autonomous system for the amplitudes $A$ and $B$ of the even and odd modes, respectively.

By using the scalings

(4.15) \begin{equation} \left.\begin{array}{c@{}} A = \left( \dfrac{4\varGamma^7}{(4+\varGamma^2)^4{\rm \pi}^7} \right)^{1/3}\hat{A} = a\hat{A}, \quad B = a \dfrac{(4+\varGamma^2)^2}{2\varGamma^3} {\rm \pi}^3 \hat{B} = b\hat{B}, \\ r =\dfrac{\varGamma}{2{\rm \pi} a} \hat{r}, \quad \tau_E = b \widehat{\tau_E}, \quad \lambda = \dfrac{\hat{\lambda}}{b^2}, \quad \sigma = \dfrac{\hat{\sigma}}{b}, \end{array}\right\} \end{equation}

we normalise the system (4.5) and (4.12) and are left with the two ordinary differential equations (after dropping hats)

(4.16a)$$\begin{gather} \dot{A} = rA - A^3 - B^2\lambda(A), \end{gather}$$

(4.16b)$$\begin{gather}\epsilon^2 \dot{B} = \sigma(A)B, \end{gather}$$

where the time derivative is taken with respect to the slow time scale $\tau _E$. In figure 9 we show the functions $\sigma (A)$ and $\lambda (A)$ (normalised according to (4.15)) in the instance when $\varGamma = 2$. The highlighted point $A_*$ where $\sigma (A_*) = 0$ takes the value $A_* \approx 83.1$.

Figure 9.

The functions (a) $\sigma (A)$ and (b) $\lambda (A)$ as they appear in (4.16) with $\varGamma =2$. Highlighted is the point $A_*$ at which $\sigma (A_*) = 0$.

Since $\sigma (A)$ is an even function and $\lambda (A)$ is an odd function, we need only consider solutions of the phase plane problem (4.16) in the quadrant $A,B\geq 0$. The critical points are the following.

1. (i) Point $(A,B) = (0,0)$ corresponds to the pure conduction state, which loses stability through the primary pitchfork bifurcation as $r$ increases through zero. This bifurcation excites the even modes in the system.

2. (ii) Point $(\sqrt {r},0)$ exists for $r > 0$ and represents the SCRS, with no odd modes. This state loses stability through a secondary pitchfork bifurcation at $r =r_c= A_*^{2} \approx 6.91 \times 10^3$ (when $\varGamma =2$), at which point the SCRS becomes unstable to odd perturbations.

3. (iii) Point $(A_*, \sqrt {(rA_* - A_*^3)/\lambda (A_*)})$ exists for $r>r_c$, when a mixed steady state including odd and even modes emerges. Finally, this steady state loses stability in a Hopf bifurcation at $r = r_* = A_*^2 + (2A_*^2\lambda (A_*))/(\lambda (A_*)-A_*\lambda '(A_*))$. When $\varGamma = 2$, we calculate $r_* \approx 1.53\times 10^4$.

We can use the value of $r_c$ from (ii) to infer the prefactor in the power law found in § 4 for the stability boundary at small $Pr$. Rescaling back to our original variables using (4.15), we find that $r_c=\delta Ra _c/Pr^2 \approx 6.35 \times 10^4$ which indeed gives excellent agreement with the fit obtained in figure 7.

To see the bifurcations itemised above in practice, we solve (4.16) numerically. Figure 10(a) is a bifurcation diagram showing how the non-trivial steady states identified in (ii) and (iii) above emerge at $r = 0$ and $r = r_c=A_*^2$, respectively. The Hopf bifurcation occurs at $r=r_*$, beyond which point the system undergoes the nonlinear relaxation oscillations and bursts displayed in figures 10(b) and 10(c) for $r = 2 \times 10^4$ and $r = 3 \times 10^4$, respectively.

Figure 10.

(a) Bifurcation diagram plotting $A/A_*$ (blue) and $\epsilon B$ (red) against $r$ with $\varGamma =2$ and $\epsilon = 10^{-2}$. Annotated are $r=r_c= A_*^2$, the point of the secondary pitchfork bifurcation, and $r= r_*$, where the Hopf bifurcation occurs. Time series of $A/A_*$ and $\epsilon B$ when (b) $r = 2 \times 10^4$ and (c) $r = 3 \times 10^4$.

There is great qualitative agreement between the behaviours of our vastly simplified two-dimensional system and of the full system, displayed in figures 10 and 8, respectively. Although the system (4.16) was formally derived in the asymptotic limit as $Pr\rightarrow 0$, similar relaxation oscillations are found in many parts of the $(Pr, Ra)$ plane (Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014). The resemblance of these oscillations to the behaviour of simple predator–prey population models and the Lotka–Volterra equations has often been studied (Leboeuf, Charlton & Carreras Reference Leboeuf, Charlton and Carreras1993; Malkov, Diamond & Rosenbluth Reference Malkov, Diamond and Rosenbluth2001; Garcia et al. Reference Garcia, Bian, Paulsen, Benkadda and Rypdal2003; Decristoforo, Theodorsen & Garcia Reference Decristoforo, Theodorsen and Garcia2020). In this analogy, the place of the predator population is taken by a quantity undergoing relaxation oscillations, and the place of the prey population is taken by a bursting quantity. In contrast, our model (4.16) has been derived systematically from the governing equations without the need for any ad hoc closure assumptions.

5. Stability at large Prandtl number

Figure 11 displays the $(Pr, Ra)$ plane for $10^5< Pr<10^6$, highlighting where the SCRS is stable (white, $S$) or unstable (blue, $US$). At $q_4 = (7 \times 10^5,2.54 \times 10^7)$ there is a corner in the stability boundary and two complex conjugate pairs of eigenvalues cross the imaginary axis simultaneously, as was discussed in § 3.

Figure 11.

The $(Pr, Ra)$ parameter space for large $Pr$ highlighting when the SCRS is linearly stable (white, $S$) and unstable (blue, $US$). The point $q_4$ is where the stability boundaries corresponding to the symmetry breaking (red, $E_{SB}$) and period doubling (green, $E_{PD}$) eigenfunctions cross. The insets show the odd eigenfunctions, $\psi ^O$, associated with $E_{SB}$ and $E_{PD}$ where red, green and blue indicate $\psi ^O > 0$, $\psi ^O \approx 0$ and $\psi ^O < 0$, respectively. The parameters are $(Ra,Pr) = (1.56 \times 10^7, 10^5)$ and $(Ra,Pr) = (2.54 \times 10^7, 7 \times 10^5)$.

The stability boundaries for the two most unstable eigenfunctions therefore cross at $q_4$, as indicated by the dashed curves in the inset in figure 11. As noted in § 3, the symmetry breaking eigenfunction (labelled $E_{SB}$) that is excited as we cross the red dashed curve has the property that $\hat \psi ^o_{k_x,k_y}$ and $\mathrm {i}\hat \theta ^o_{k_x,k_y}$ are purely real and therefore breaks the reflection symmetry between counter-rotating convection rolls in the SCRS. It follows that the horizontal velocity is an even function of $x$ which allows for instantaneous zonal flow in the solution.

The period doubling eigenfunction (labelled $E_{PD}$) corresponding to the green dashed curve has the complementary symmetry that $\hat \psi ^o_{k_x,k_y}$ and $\mathrm {i}\hat \theta ^o_{k_x,k_y}$ are purely imaginary, and the zonal modes are therefore all zero: $\hat \psi ^o_{0,k_y}=0$ for all $k_y$. The horizontal velocity is now an odd function of $x$ and, although the $E_{PD}$ eigenfunction introduces odd modes to the solution, it does not break the reflection symmetry (2.10a,b) in the SCRS, but instead causes a period doubling in $x$. Consequently, zonal modes at the instability switch off as we pass the point $q_4$. The qualitative change in the structure of the dominant eigenfunction is demonstrated by the insets in figure 11, and how these eigenfunctions evolve in time is shown in the supplementary material. To measure the proportion of energy held in zonal modes in the eigenfunction $E_{SB}$, we define

(5.1)\begin{equation} R_{{zonal}} = \frac{\displaystyle \sum_{k_x} |K_{k_x,0} \hat \psi_{k_x,0}|^2}{\displaystyle \sum_{k_x,k_y} |K_{k_x,k_y} \hat \psi_{k_x,k_y}|^2}, \end{equation}

where $K_{k_x,k_y} = \sqrt {({2{\rm \pi} k_x}/{\varGamma })^2 + ({\rm \pi} k_{y})^2}$ is the non-dimensional wavenumber. As shown in figure 12, $R_{{zonal}}$ decreases as we increase $Pr$ along the stability boundary. Nevertheless, $E_{SB}$ still has a small but non-zero zonal component at the point $q_4$ where the eigenfunction $E_{PD}$ takes over and zonal modes disappear completely. The instance of $E_{SB}$ in figure 11 is at the large Prandtl number $Pr = 10^5$, and hence the zonal modes do not carry much energy. Supplementary movie 6 shows $E_{SB}$ with $Pr = 100$, in which case the presence of zonal modes is clear.

Figure 12.

The proportion of energy held in zonal modes, $R_{{zonal}}$ defined in (5.1), in the eigenfunction $E_{SB}$ at the stability threshold plotted versus $Pr$. The red dashed line highlights the point $q_4$ beyond which the eigenfunction $E_{PD}$, with no zonal modes, becomes the most unstable.

As we increase the Prandtl number past $q_4$, the green dashed line in figure 11 appears to approach a limiting value $Ra = 2.54 \times 10^7$, at which the SCRS is unstable for all $Pr$. To examine this statement, we consider the $Pr \to \infty$ limit in which the governing equations (2.1) reduce to

(5.2a)$$\begin{gather} \nabla^4 \psi ={-}Ra\theta_x, \end{gather}$$

(5.2b)$$\begin{gather}\theta_t + \{ \psi, \theta\} = \psi_x + \nabla^2 \theta. \end{gather}$$

We analyse the linear stability of the SCRS for the reduced system (5.2) as in § 3, and we find that it becomes unstable at $Ra = 2.54 \times 10^7$, consistent with the results presented in figure 11. It is clear from (5.2a) that $\hat \psi _{0,k_y} = 0$ in the $Pr \to \infty$ limit, again in agreement with the period doubling eigenfunction identified above. We conclude that $Pr > 7 \times 10^5$ is a large enough Prandtl number such that the system exhibits the asymptotic behaviour of identically zero zonal modes, at least linearly at the SCRS instability.

We note that the disparity between the thermal and viscous time scales makes it very difficult to reproduce the behaviour close to $q_4$ using DNS. The linear stability analysis predicts that $E_{PD}$ oscillates with a frequency of order 10 with respect to the thermal time scale, and it is therefore necessary to integrate the underlying equations over at least $Pr/10\approx 7\times 10^4$ viscous time units to observe just one complete oscillation. At the required spatial resolution, it proved unfeasible to compute enough cycles to reliably observe the exponential growth or decay associated with the Hopf bifurcation.