3.1. Linear stability analysis

We first consider small-amplitude stationary normal mode perturbations to the conduction state

(3.1)\begin{equation} (u,w,p,\varTheta,\varPhi)^{\rm T} = \epsilon(\left(U_1(x),W_1(x),P_1(x),\varTheta_1(x),\varPhi_1(x)\right)^{\rm T}{\rm e}^{{\rm i}kz} + {\rm c.c.}) + O(\epsilon^2), \end{equation}

where c.c. denotes the complex conjugate of the preceding term, $\epsilon \ll 1$, and $k$ is the vertical wavenumber and $\lambda$ is the temporal growth rate of the perturbation. Inserting expansion (3.1) into system (2.17)–(2.20) and linearising the resulting system yields the eigenvalue problem

(3.2)\begin{equation} \mathcal{L}(Ra)\boldsymbol{\varPsi}_{\boldsymbol{1}} = \boldsymbol{0}, \end{equation}

for $Ra$ and $\boldsymbol {\varPsi }_1$ where

(3.3)\begin{equation} \boldsymbol{\varPsi}_{\boldsymbol{1}} = \left(U_1,W_1,P_1,\varTheta_1,\varPhi_1\right)^{\rm T}{\rm e}^{{\rm i}kz} + {\rm c.c.}, \end{equation}

and

(3.4)\begin{equation} \mathcal{L}(Ra) = \begin{pmatrix} {\nabla}^2 & 0 & -\partial_x & 0 & 0 \\ 0 & {\nabla}^2 & -\partial_z & Ra & -Ra \\ \partial_x & \partial_z & 0 & 0 & 0 \\ -1 & 0 & 0 & {\nabla}^2 & 0 \\ -1 & 0 & 0 & 0 & \dfrac{1}{Le}{\nabla}^2 \end{pmatrix}.\end{equation}

The complex functions $U_1$, $W_1$, $P_1$, $\varTheta _1$ and $\varPhi _1$ satisfy Dirichlet boundary conditions on the sidewalls for the velocity, temperature and concentration perturbations

(3.5)\begin{equation} U_1(x) = W_1(x) = \varTheta_1(x) = \varPhi_1(x) = 0\quad \text{on } x = 0,1, \end{equation}

and the projection of the Navier–Stokes equation onto the boundary for the pressure perturbation

(3.6)\begin{equation} 0 ={-}\frac{\partial P_1}{\partial x} + \frac{\partial^2U_1}{\partial^2x} \quad \text{on } x = 0,1.\end{equation}

Solutions to (3.2)–(3.6) are independent of $Pr$ and satisfy $\varPhi _1 = Le\,\varTheta _1$. Consequently, the only parameter dependence in the linear problem comes from the buoyancy term in the momentum equation, which takes the form $Ra(1-Le)\varTheta _1$. The accordingly simplified version of (3.2) is then solved using a spectral eigenvalue solver based on a Chebyshev–Legendre collocation method for a range of $k$ to determine the marginal stability curve in figure 2. This curve reveals a minimum at $k_c \approx 2.5318$ and $Ra_c|Le-1| \approx 6509$, which corresponds to the primary instability of the conduction state. The absolute value here comes from the fact that the system resulting from left wall heating and that resulting from right wall heating are equivalent. Contour plots presenting the fields for the velocity components, streamfunction and convective variables of this eigenmode for $Le=11$ are shown in figure 3. The conduction state is thus first unstable to a spatially periodic state constituted of counter-rotating convection rolls that, when $Le > 1$, slant downwards from the hotter wall, filling the domain and extending to the cold wall.

Figure 2.

Marginal stability curve for the onset of doubly diffusive convection. The conduction state is stable to modes with wavenumber $k$ below the curve, and unstable to them above. The minimum of this curve is $Ra_c|Le-1| \approx 6509$ with wavenumber $k_c \approx 2.5318$ and corresponds to the location of the primary bifurcation.

Figure 3.

Contour plots of a single wavelength of the real critical eigenvector $\boldsymbol {\varPsi }_1$ for $Le = 11$ ($k_c \approx 2.53$). The profiles show the perturbations in (a) horizontal velocity, $u$, (b) vertical velocity, $w$, (c) velocity streamfunction, $\psi$, where $u = -{\psi }_z$ and $w = {\psi }_x$, and in the convective variables (d) $\varTheta$ and (e) $\varPhi$. Black (grey, dotted) lines indicate positive (negative, zero) values and are separated by 20 % of the maximum absolute value.

The conduction state can also undergo Hopf bifurcations, where the growth rate is purely imaginary. However, for the parameter values tested, these bifurcations occurred for Rayleigh numbers that are orders of magnitude larger than that of the primary stationary bifurcation and are therefore out of the scope of the present work.

3.2. Weakly nonlinear analysis

To investigate the weakly nonlinear regime around this primary bifurcation, we set $Ra = Ra_c + \epsilon ^2 r$ with $r = {O}(1)$ and $\epsilon \ll 1$ and assume that the system evolves on a slow temporal scale $T_1 = \epsilon ^2 t$. We also introduce a long spatial scale, $Z = \epsilon z$, to allow small-amplitude states with long spatial modulations. The small-aspect-ratio domains considered in the numerical computations in § 4 do not permit these long-scale modulations, so terms involving derivatives with respect to $Z$ may be ignored in the subsequent analysis with no effect on the main result of this section: the criticality of the primary bifurcation. However, this long spatial variable has been included here to broaden the scope of the analysis and will be considered in future work. We emphasise that each of the state variables of our system ($u$, $w$, $p$, $\varTheta$ and $\varPhi$) depend upon the independent variables: $x$, $z$, $Z$ and $T_1$. Using this multiple-scale approach, the partial derivatives become

(3.7a,b)\begin{equation} \dfrac{\partial }{\partial t} \mapsto \epsilon^2 \dfrac{\partial}{\partial T_1} \quad \text{and} \quad \dfrac{\partial }{\partial z} \mapsto \dfrac{\partial}{\partial z} + \epsilon \dfrac{\partial}{\partial Z}. \end{equation}

Introducing the notation $\boldsymbol {\varPsi } = (u,w,p,\varTheta,\varPhi )^{\rm T}$, we can express each of the variables as a perturbation expansion in $\epsilon$ about the conduction state $\boldsymbol {\varPsi }_0 = (0,0,0,0,0)^{\rm T}$

(3.8)\begin{equation} \boldsymbol{\varPsi} = \boldsymbol{\varPsi}_0 + \epsilon \boldsymbol{\varPsi}_1 + \epsilon^2 \boldsymbol{\varPsi}_2 + \dots, \end{equation}

where $\boldsymbol {\varPsi }_j = (u_j,w_j,p_j,\theta _j,\phi _j)^{\rm T}$ is the correction at $O(\epsilon ^j)$ for $j = 1,2,\dots$.

The corrections are periodic in $z$ and also satisfy homogeneous boundary conditions at each order in $\epsilon$

(3.9)\begin{equation} u_j = w_j = \theta_j = \phi_j = 0 \quad \text{on } x = 0,1, \quad j = 1,2,\dots, \end{equation}

as well as the pressure boundary condition

(3.10)\begin{equation} \frac{\partial^2 u_j}{\partial x^2} - \frac{\partial p_j}{\partial x} = 0 \quad \text{on } x = 0,1, \quad j = 1,2,\dots. \end{equation}

The expansion (3.8) is substituted into the full system (2.17)–(2.20) and the perturbations are solved numerically order by order in $\epsilon$ using an extension of the aforementioned collocation method. By further extracting the parameter dependence of the perturbations at each order, we obtain a Ginzburg–Landau equation that can be applied for all parameter values and will indicate the criticality of the primary bifurcation.

We proceed by detailing this formulation, which should be applied to the cases $Le > 1$ and $Le < 1$ separately, owing to the parameter combination $Ra(1-Le)$ changing sign between them. However, the results for $Le <1$ can be related to those for $Le>1$ by using an alternative non-dimensionalisation to (2.2a–f) involving the solutal diffusivity, $D$, instead of thermal diffusivity, $\kappa$, which results in a set of equations like (2.3)–(2.6), except with $T$ and $C$ exchanged and the Lewis, Prandtl and Rayleigh numbers replaced by the inverse Lewis number, Schmidt number, $Sc = LePr$, and solutal Rayleigh number, $Ra_S = -RaNLe$, respectively.

The conduction state solves the system at leading order. At ${O}(\epsilon )$, the correction is given by the solution to linear system (3.2)

(3.11)\begin{equation} \boldsymbol{\varPsi}_1 = A_1(Z,T_1) \left(U_1(x),W_1(x),P_1(x),\varTheta_1(x),Le \varTheta_1(x)\right)^{\rm T}{\rm e}^{{\rm i}k_cz} + {\rm c.c.},\end{equation}

where $k_c$ is the critical wavenumber. No phase constraint is applied at this point, but the amplitude of the eigenfunction is fixed using

(3.12)\begin{equation} \langle U_1,U_1\rangle + \langle W_1,W_1\rangle+ \langle P_1,P_1\rangle + \langle \varTheta_1,\varTheta_1\rangle = 1, \end{equation}

with the inner product

(3.13)\begin{equation} \langle\, f,g\rangle = \frac{1}{\lambda_c}\int_0^{\lambda_c}\int_0^1 \bar{f}^{\rm T}\, g\,{{\rm d}\kern0.06em x}\,{\rm d}z, \end{equation}

where $\lambda _c = 2{\rm \pi} /k_c$ is the wavelength of the critical eigenvector, the overbar denotes complex conjugation and the superscript T denotes the transposition operation when $f$ is a vector. Due to the lack of available explicit expressions for the solutions to this perturbation problem, these inner products need to be computed numerically, which we achieved by using the Clenshaw–Curtis quadrature on the collocation nodes used in § 3.1. The amplitude $A_1$ evolves over both long spatial and temporal scales according to an amplitude equation that will be determined at higher order.

At ${O}(\epsilon ^2)$, the linear operator $\mathcal {L}$ acts on the second-order terms and is forced by both the nonlinear terms between the ${O}(\epsilon )$ corrections and terms proportional to the slow spatial derivative of the $O(\epsilon )$ correction $A_{1Z}$

(3.14)\begin{equation} \mathcal{L}(Ra_c)\boldsymbol{\varPsi}_2 =\underbrace{\begin{pmatrix} \displaystyle\frac{1}{Pr}f_{10} |A_1|^2 + \left(A_{1Z}f_{11}{\rm e}^{{\rm i}k_cz} + {\rm c.c.}\right) + \displaystyle\frac{1}{Pr}\left(f_{12} A_1^2{\rm e}^{2{\rm i}k_cz} + {\rm c.c.}\right)\\ \displaystyle\frac{1}{Pr}f_{20} |A_1|^2 + \left(A_{1Z}f_{21}{\rm e}^{{\rm i}k_cz} + {\rm c.c.}\right)+\displaystyle\frac{1}{Pr}\left(f_{22} A_1^2{\rm e}^{2{\rm i}k_cz} + {\rm c.c.}\right)\\ \hphantom{f_{20} |A_1|^2 +}\left(A_{1Z}f_{31}{\rm e}^{{\rm i}k_cz} + {\rm c.c.}\right)\hphantom{+\left(f_{22} A_1^2{\rm e}^{2{\rm i}k_cz} + {\rm c.c.}\right)}\\ f_{40} |A_1|^2 + \left(A_{1Z}f_{41}{\rm e}^{{\rm i}k_cz} + {\rm c.c.}\right)+\left(f_{42} A_1^2{\rm e}^{2{\rm i}k_cz} + {\rm c.c.}\right)\\ Lef_{40} |A_1|^2 + \left(A_{1Z}f_{41}{\rm e}^{{\rm i}k_cz} + {\rm c.c.}\right)+ Le\left(f_{42} A_1^2{\rm e}^{2{\rm i}k_cz} + {\rm c.c.}\right) \end{pmatrix}}_{\mathcal{N}_2}, \end{equation}

where the functions $f_{ij}(x)$ for $i = 1,2,3,4$ and $j = 0,1,2$ are independent of $Pr$ and $Le$ and are given in table 1.

Table 1.

Functions $f_{ij}$ ($i = 1,2,3,4$, $j = 0,1,2$) in the nonlinear term $\mathcal {N}_2$ at ${O}(\epsilon ^2)$ in (3.14). The overbar denotes complex conjugation.

To ensure the existence of a unique solution at this order, we derive a solvability condition using the Fredholm alternative theorem. This involves the adjoint operator to $\mathcal {L}$, $\mathcal {L}^{{\dagger} }$, defined through the relationship

(3.15)\begin{equation} \langle\, f,\mathcal{L}g\rangle = \langle \mathcal{L}^{{{\dagger}}}f,g\rangle, \end{equation}

which holds for all vector functions $f$ and $g$. Integrating the left-hand side by parts, we find that the adjoint operator takes the form

(3.16)\begin{equation} \mathcal{L}^{{{\dagger}}} = \begin{pmatrix} {\nabla}^2 & 0 & -\partial_x & -1 & -1 \\ 0 & {\nabla}^2 & -\partial_z & 0 & 0 \\ \partial_x & \partial_z & 0 & 0 & 0 \\ 0 & Ra_c & 0 & {\nabla}^2 & 0\\ 0 & -Ra_c & 0 & 0 & \dfrac{1}{Le}{\nabla}^2 \end{pmatrix}, \end{equation}

together with the adjoint boundary conditions

(3.17)$$\begin{gather} u^{{{\dagger}}} = 0,\quad w^{{{\dagger}}} = 0,\quad \theta^{{{\dagger}}} = 0,\quad \phi^{{{\dagger}}} = 0 \quad \text{on}\ x = 0,1, \end{gather}$$

(3.18)$$\begin{gather}\frac{\partial^2 u^{{{\dagger}}}}{\partial x^2} -\frac{\partial p^{{{\dagger}}}}{\partial x} = 0\quad \text{on}\ x = 0,1, \end{gather}$$

and periodicity in the vertical direction.

The Fredholm alternative allows us to pose the adjoint problem

(3.19)\begin{equation} \mathcal{L}^{{{\dagger}}} \boldsymbol{\varPsi}^{{{\dagger}}} = \boldsymbol{0}, \end{equation}

whose solution is unique up to a vertical translation and a multiplicative constant. This solution may be written in the form

(3.20)\begin{equation} \boldsymbol{\varPsi}^{{{\dagger}}} = \left(U^{{{\dagger}}}(x), W^{{{\dagger}}}(x),P^{{{\dagger}}}(x), \frac{1}{1-Le}\varTheta^{{{\dagger}}}(x),-\frac{Le}{1-Le}\varTheta^{{{\dagger}}}(x) \right)^{\rm T}{\rm e}^{{\rm i}k_cz} + {\rm c.c.}, \end{equation}

where the parameter dependence of the components has been extracted. The amplitude and phase are fixed by imposing the conditions

(3.21)\begin{equation} \langle U^{{{\dagger}}},U^{{{\dagger}}}\rangle + \langle W^{{{\dagger}}},W^{{{\dagger}}}\rangle+ \langle P^{{{\dagger}}},P^{{{\dagger}}}\rangle + \langle \varTheta^{{{\dagger}}},\varTheta^{{{\dagger}}}\rangle = 1, \end{equation}

and

(3.22)\begin{equation} {\rm Im}(\langle U^{{{\dagger}}},U_1\rangle) = 0, \end{equation}

where Im represents the imaginary part, respectively.

Using this adjoint solution, we then apply the ${O}(\epsilon ^2)$ solvability condition

(3.23)\begin{equation} \langle \varPsi^{{{\dagger}}},\mathcal{N}_2\rangle = 0. \end{equation}

Owing to the vertical wavenumber dependence of terms in $\mathcal {N}_2$, the only non-trivial contributions come from those proportional to $A_{1Z}$ and their complex conjugates and (3.23) then reduces to

(3.24)\begin{equation} -2{\rm i}k_c\langle U^{{{\dagger}}},U_1\rangle -2{\rm i}k_c\langle W^{{{\dagger}}},W_1\rangle + \langle W^{{{\dagger}}},P_1\rangle - \langle P^{{{\dagger}}},W_1\rangle -2{\rm i}k_c \langle \varTheta^{{{\dagger}}},\varTheta_1\rangle= 0, \end{equation}

which may be further simplified to

(3.25)\begin{equation} \left\langle \boldsymbol{\varPsi}^{\dagger},\dfrac{\partial \mathcal{L}\boldsymbol{\varPsi}_1}{\partial k_c}\right\rangle = 0. \end{equation}

Thus, this solvability condition is automatically satisfied as the primary bifurcation occurs at a quadratic minimum of the marginal stability curve (see figure 2).

The ${O}(\epsilon ^2)$ system (3.14) can be solved to find that the second-order correction to the conduction state is

(3.26)\begin{equation} \boldsymbol{\varPsi}_2 = |A_1|^2\boldsymbol{\varPsi}_2^0 + ((A_2\boldsymbol{\varPsi}_1 +A_{1Z}\boldsymbol{\varPsi}_2^1){\rm e}^{{\rm i}k_cz} + {\rm c.c.}) + (A_1^2\boldsymbol{\varPsi}_2^2{\rm e}^{2{\rm i}k_cz} + {\rm c.c.}),\end{equation}

where $\boldsymbol{\varPsi}_2 = (u_2,w_2,p_2,\theta _2,\phi _2)^{\rm T}$ and the functions $\boldsymbol{\varPsi}_2^i$ for $i=0,1,2$ have the following parameter dependence:

(3.27)\begin{gather} \boldsymbol{\varPsi}_2^0 = \left(0,\frac{1}{Pr}\tilde{w}_2 + (1+Le)\tilde{w}_3,\frac{1}{Pr}\tilde{p}_2,\tilde{\theta}_3, Le^2\tilde{\theta}_3\right)^{\rm T},\end{gather}

(3.28)\begin{gather} \boldsymbol{\varPsi}_2^1 = \left(\tilde{u}_7,\tilde{w}_7,\tilde{p}_7,\tilde{\theta}_7,Le\,\tilde{\theta}_7\right)^{\rm T}, \end{gather}

(3.29)\begin{gather} \boldsymbol{\varPsi}_2^2 = \left(\frac{1}{Pr}\left(\tilde{u}_4,\tilde{w}_4,\tilde{p}_4,\tilde{\theta}_4,Le\tilde{\theta}_4\right)+(1+Le)(\tilde{u}_5,\tilde{w}_5,\tilde{p}_5,0,0) \right.\nonumber\\ \left.\qquad \vphantom{\frac{1}{Pr}}+ (0,0,0,\tilde{\theta}_5,Le^2\tilde{\theta}_5)+Le(0,0,0,\tilde{\theta}_6, \tilde{\theta}_6)\right)^{\rm T}. \end{gather}

The newly introduced functions $\tilde {u}_i$, $\tilde {w}_i$, $\tilde {p}_i$ and $\tilde {\theta }_i$ for $i = 2,\ldots,7$ are independent of $Le$ and $Pr$ and satisfy (A1)–(A5) in § A.1 of the Appendix A.

Continuing to ${O}(\epsilon ^3)$, both the deviation away from the critical Rayleigh number and the slow time dependence of the solution appear in the right-hand side of the resulting system, in addition to nonlinear terms between first and second-order corrections and terms with slow spatial derivatives. The system to solve at third order is

(3.30)\begin{equation} \mathcal{L}(Ra_c)\boldsymbol{\varPsi}_3 = \underbrace{\begin{pmatrix}\dfrac{1}{Pr}\left(\dfrac{\partial u_1}{\partial_{T_1}} + J(\boldsymbol{u},u) \right) -2\dfrac{\partial^2u_2}{\partial z\partial Z} - \dfrac{\partial^2u_1}{\partial Z^2}\\ \dfrac{1}{Pr}\left(\dfrac{\partial w_1}{\partial_{T_1}} + J(\boldsymbol{u},w)\right) -r\left(\theta_1 - \phi_1\right) -2\dfrac{\partial^2w_2}{\partial z\partial Z} - \dfrac{\partial^2w_1}{\partial Z^2} +\dfrac{\partial p_2}{\partial Z}\\ -\dfrac{\partial w_2}{\partial Z}\\ \left(\dfrac{\partial \theta_1}{\partial_{T_1}} +J(\boldsymbol{u},\theta) \right) -2\dfrac{\partial^2\theta_2}{\partial z\partial Z} - \dfrac{\partial^2\theta_1}{\partial Z^2} \\ \left(\dfrac{\partial \phi_1}{\partial_{T_1}} +J(\boldsymbol{u},\phi)\right) -\dfrac{2}{Le}\dfrac{\partial^2\phi_2}{\partial z\partial Z} - \dfrac{1}{Le}\dfrac{\partial^2 \phi_1}{\partial Z^2} \end{pmatrix}}_{\mathcal{N}_3},\end{equation}

where the advective terms are

(3.31)\begin{equation} J(\boldsymbol{u},f) = \boldsymbol{u}_1\boldsymbol{\cdot} \boldsymbol{\nabla} f_2 + \boldsymbol{u}_2\boldsymbol{\cdot} \boldsymbol{\nabla} f_1 + w_1\partial_Zf_1, \end{equation}

where $f_1$ and $f_2$, respectively, refer to the first- and second-order corrections of the variables $f= u,w,\theta$ and $\phi$.

The solvability condition at this order

(3.32)\begin{equation} \langle \boldsymbol{\varPsi}^{{{\dagger}}},\mathcal{N}_3\rangle = 0, \end{equation}

is no longer trivially satisfied because some nonlinear terms contained in $\mathcal {N}_3$ have $e^{ik_cz}$ dependence arising from terms proportional to $A_1$, $|A_1|^2A_1$, $A_{1ZZ}$, $A_{2Z}$ and their complex conjugates. However, the contributions to (3.32) from terms proportional to $A_{2Z}$, cancel for the same reason that the solvability condition at $O(\epsilon ^2)$ was satisfied. Consequently, $A_2$ remains arbitrary at this order. Collecting the remaining terms from (3.32), we obtain the Ginzburg–Landau equation, that holds for both $Le > 1$ and $Le < 1$

(3.33)\begin{equation} \alpha A_{1T_1} = \gamma r A_1 + \beta |A_1|^2A_1 + \delta A_{1ZZ}, \end{equation}

where table 2 indicates which terms of $\mathcal {N}_3$ contribute to each term above. This equation is equivariant under the $O(2)$ symmetry so we may choose the phase of the ${O}(\epsilon )$ correction so that these coefficients are real. The coefficient $\delta$ is independent of the physical parameters $Pr$ and $Le$, while $\alpha, \beta$ and $\gamma$ satisfy the relations

(3.34)$$\begin{gather} \alpha = \frac{1}{Pr}\alpha_1 + (1+Le)\alpha_2, \end{gather}$$

(3.35)$$\begin{gather}\beta = \frac{1}{Pr^2}\beta_1 + \frac{1+Le}{Pr}\beta_2 + (1+Le^2)\beta_3+Le\beta_4, \end{gather}$$

(3.36)$$\begin{gather}\gamma = (1-Le)\gamma_1, \end{gather}$$

where full expressions used to obtain $\alpha _i$, $\beta _i$, $\gamma _1$ and $\delta$ are provided in (A6)–(A9) and evaluated in table 4 in § A.2. By dividing (3.33) through by $\alpha$, (3.33) is more conveniently written as

(3.37)\begin{equation} A_{1T_1}= a_1rA_1 + a_2 |A_1|^2A_1 + a_3A_{1ZZ} , \end{equation}

where $a_1 = \gamma /\alpha$, $a_2 = \beta /\alpha$ and $a_3 = \delta / \alpha$.

Table 2.

Terms from $\mathcal {N}_3$ (see (3.30)) contributing to the Ginzburg–Landau equation (3.33). The column in which these terms are placed informs on the term to which they contribute. Here, $f_1$ and $f_2$, respectively refer to first- and second-order corrections of the variables $f= u,w,\theta$ and $\phi$.

The solutions to the Ginzburg–Landau equation (3.37) are good approximations of the small-amplitude solutions of the full doubly diffusive system (2.17)–(2.20). Of particular interest here are the two steady solutions that are invariant with respect to the long spatial scale $Z$. The first of these solutions is the trivial solution

(3.38)\begin{equation} A_1 = 0. \end{equation}

This solution is valid for all $r$ and corresponds to the conduction state (2.16a–c). The second important solution is

(3.39)\begin{equation} A_1 = \left(-\frac{a_1 r}{a_2}\right)^{1/2}{\rm e}^{{\rm i}\chi}, \end{equation}

where $\chi$ is an arbitrary phase. This solution relates to states of small-amplitude spatially periodic convection that can be found near the primary bifurcation. These fluid states can then be approximated by

(3.40)\begin{equation} \left(u,w,p,\varTheta,\varPhi\right)^{\rm T} \approx \sqrt{-\frac{a_1 (Ra-Ra_c)}{a_2}} \left(U_1(x),W_1(x),P_1(x),\varTheta_1(x),Le \varTheta_1(x)\right)^{\rm T}{\rm e}^{{\rm i}k_cz} + {\rm c.c.}, \end{equation}

where the phase $\chi$ has been absorbed into $z$ via a vertical translation. These states only exist at small amplitude for Rayleigh numbers that satisfy

(3.41)\begin{equation} \frac{a_1}{a_2}\left(Ra_c-Ra\right) > 0. \end{equation}

Consequently, the sign of the ratio $a_1/a_2$ determines the criticality of the primary bifurcation and the initial direction of branching.

Using the numerical values in table 4, we computed $a_i$ over a range of parameter values. The coefficients $a_1$ and $a_3$ are positive for all $Pr$ provided $Le \neq 1$, whereas the sign of $a_2$ changes as these parameters are varied. As a result, there exists a boundary in parameter space that separates regions where the primary bifurcation is subcritical ($a_2 >0$) from those where it is supercritical ($a_2 < 0$). This boundary is shown in figure 4 and implies that, for any value of the Lewis number, there exists a critical value of the Prandtl number, $Pr_c(Le)$, expressed in terms of the physical parameters in (A15), above which the bifurcation is subcritical. This critical value tends to $0.376$ for small Lewis numbers while it approaches the asymptotic relation $Pr_c \sim 0.376/Le$ as the Lewis number tends to infinity. We further note that the parameter values for physical doubly diffusive systems from Schmitt (Reference Schmitt1983) all lie within the region where the primary bifurcation is subcritical. While we are unaware of further fluid systems lying within the supercritical region of parameter space, we expect that they exist since some of the physical systems identified in figure 4, including humidity/heat and stellar interiors (marked (6) and (8), respectively), have parameter values that are within an order of magnitude of the sub/supercritical boundary.

Figure 4.

Boundary $a_2 = 0$ in $(Le,Pr)$ parameter space separating the region where the primary bifurcation from the conduction state is subcritical (above) from that where it is supercritical (below). The conduction state is linearly stable for all $Ra$ at $Le = 1$ and this point is indicated by the open circle. The grey regions indicate parameter values from Schmitt (Reference Schmitt1983) for the physical doubly diffusive systems: (1) salt/sugar, (2) magmas, (3) oxide semiconductors, (4) heat/salt $0\,^{\circ }$C, (5) heat/salt $30\,^{\circ }$C, (6) humidity/heat, (7) liquid metals and (8) stellar interiors.

We can gain physical insight into the criticality of the primary bifurcation by examining the contributions that each of the nonlinear terms from (2.3)–(2.6) make to $a_2$, using a similar approach to the one Requilé et al. (Reference Requilé, Hirata, Ouarzazi and Barletta2020) applied to plane Poiseuille and plane Couette flows with viscous dissipation. The expression of the coefficient $\beta$ (3.35) and the corresponding numerical values provided in table 4 in the Appendix A, show that the inertial term $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}$ (contributing to $\beta _1$ and $\beta _2$) provides a negative contribution to $a_2$, whereas thermal $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } T$ and solutal $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } C$ advective terms (mostly contributing to $\beta _3$ and $\beta _4$) provide a positive contribution to $a_2$. The latter statement is further justified in § A.3 in the Appendix A. It is therefore solutal and thermal advection in the system that drives the subcriticality of the primary bifurcation, while inertial effects drive the supercriticality. Thus, reducing the Prandtl number reduces the subcriticality of the bifurcation since the effects of inertia are strengthened.

The final term in the Ginzburg–Landau equation (3.37), $a_3A_{1ZZ}$, allows small-amplitude solutions of the doubly diffusive system to exhibit long-scale amplitude modulation. These solutions include phase-winding states that describe patterns whose wavenumbers are close to the critical wavenumber $k_c$ (Cross et al. Reference Cross, Daniels, Hohenberg and Siggia1983), and spatially modulated states that can develop into localised states away from the primary bifurcation (Bergeon & Knobloch Reference Bergeon and Knobloch2008a). These are out of the scope of the present work, but we will consider the effect of the term $a_3A_{1ZZ}$ on spatially localised states in future work.

4.1. Bifurcation structure

The results can be summarised by dividing parameter space according to the qualitative nature of the bifurcation diagram. Figure 5(a) indicates the four main regimes found. Region (1) describes the moderate and large $Pr$ behaviour for all $Le$. In this region, the primary bifurcation is strongly subcritical and the primary branch has a single saddle node, as shown in figure 5(b). Parameter values within this region have received the most attention in previous studies focusing on subcritical pattern formation (e.g. see Xin et al. Reference Xin, Le Quéré and Tuckerman1998; Bergeon & Knobloch Reference Bergeon and Knobloch2008b). Region (2) occupies a small region of parameter space above the boundary $Pr = Pr_c$, where the primary bifurcation is weakly subcritical, and separates the typical subcritical behaviour in region (1) from the supercritical behaviour in regions (3) or (4). The steady convection state branches typically have three saddle nodes in region (2), as exemplified in figure 5(c). Regions (3) and (4) identify the two qualitatively different types of bifurcation diagrams observable when the primary bifurcation is supercritical. In both cases, the primary branch has two saddle nodes, with the first lying in the supercritical region $Ra > Ra_c$. The difference between the regions is the location of the second saddle node: in region (3), it is found for $Ra < Ra_c$ (see figure 5d), whereas, in region (4), it is found in $Ra > Ra_c$ (see figure 5e). Consequently, a large-amplitude convection state may coexist with the stable conduction state when the primary bifurcation is supercritical, but, for sufficiently small $Pr$, steady convection states are found entirely within the supercritical region, where the conduction state is unstable. There may exist a fifth region, where the primary branch increases monotonically in both Rayleigh number and in amplitude, but we have not identified it in this study.

Figure 5.

(a) Enlargement of a subset of the parameter space shown in figure 4 showing four regions where the bifurcation diagrams exhibit qualitatively different behaviours. The thick line separates subcritical from supercritical branching, while the additional region boundaries are identified with either dotted or dashed lines. (b–e) Representative bifurcation diagrams for parameter values within each of the four regions. The stability of the branch segments is also indicated using thick lines for stable solutions, thin lines for solutions unstable to amplitude perturbations and dashed lines for solutions unstable to drift. The location of bifurcations depend upon the specific parameter values used, so those used for each sketch have been marked in panel (a).

We now determine the structure of the primary branch as $Pr$ decreases for a fixed value of $Le$. To achieve this, we follow the locations of its three saddle nodes with respect to $Ra$ and $Pr$. In doing so, we observed two different scenarios according to whether the pair of saddle nodes are created on the lower or upper part of the primary branch. These are exemplified in figure 6 for $Le = 11$ (representative of $5 \leqslant Le \lesssim 15$) and $Le = 20$ (representative of $19 \lesssim Le < 100$). Since the transition between the two scenarios occupies a small region of parameter space within region (2) for $15 \lesssim Le \lesssim 19$, we did not investigate it any further.

Figure 6.

Location of the three saddle-node bifurcations of the primary branch in $(Ra,Pr)$ parameter space for (a) $Le = 11$ and (b) $Le = 20$. The dashed line marks the critical Rayleigh number at which the primary bifurcation is found, $Ra_c$, and the cross marks the codimension two point $(Ra_c,Pr_c)$ explained in the text. The insets provide enlargements of the area around $Ra_c$ in each case. Arrows mark the bifurcation diagrams shown in figure 7.

To help interpret the plots in figure 6, figure 7 demonstrates the evolution of the bifurcation diagrams as $Pr$ decreases for $Le = 11$ (a–f) and $Le = 20$ (g–l). The structure of the bifurcation diagrams in region (1), for high $Pr$, remain similar, as shown in figures 7(a) and 7(g). From the primary bifurcation, the primary branch extends towards lower Rayleigh numbers and proceeds to turn around at a saddle node, hereafter referred to as $SN_1$, before heading towards large amplitude convection states at large $Ra$. Figure 6 suggests that, as $Pr\rightarrow \infty$, the location of $SN_1$ tends to a constant Rayleigh number, dependent upon $Le$. This figure also shows that $SN_1$ occurs at larger $Ra$ as the Prandtl number is decreased and the primary bifurcation becomes less subcritical.

Figure 7.

Bifurcation diagrams showing the primary branch of steady convection and the stability of the related states across the four regions, indicated in the top left corner. Thick lines indicate stable solutions, thin lines indicate solutions unstable to amplitude perturbations and dashed lines indicate solutions unstable to drift. Saddle nodes are marked by symbols: $SN_1$ (filled circle), $SN_2$ (asterisk) and $SN_3$ (triangle). The open circle corresponds to the destabilising drift bifurcation. The parameter values, also indicated by the arrows in figure 6, are: $Le = 11$, and (a) $Pr = 1$, (b) $Pr = 0.1$, (c) $Pr = 0.042$, (d) $Pr = 0.032$, (e) $Pr = 0.01$, (f) $Pr = 0.005$; as well as $Le = 20$, and (g) $Pr = 1$, (h) $Pr = 0.1$, (i) $Pr = 0.023$, (j) $Pr = 0.02$, (k) $Pr = 0.01$ and (l) $Pr = 0.005$. For $Le = 11$ (respectively $Le = 20$), $Pr_c \approx 0.031$, $Pr_{cusp} \approx 0.033$ (respectively $Pr_c \approx 0.018$, $Pr_{cusp}\approx 0.023$).

Upon decreasing the Prandtl number, the primary branch undergoes a cusp bifurcation at $Pr \approx Pr_{cusp}(Le)> Pr_c(Le)$, while still subcritical, denoting the beginning of region (2). The cusp produces two additional saddle nodes along the primary branch: $SN_2$ and $SN_3$. The exact process by which this is achieved depends on the Lewis number. For $Le \lesssim 15$, the cusp bifurcation occurs at smaller amplitude than $SN_1$ and the saddle nodes are labelled $SN_3$, $SN_2$, $SN_1$ as the branch is followed in the direction of increasing energy (see, for e.g. figure 7(d) for $Le = 11$ and $Pr = 0.032$, near the cusp parameter value: $Pr_{cusp}\approx 0.033$). In contrast, for $Le \gtrsim 19$, the cusp bifurcation occurs at higher amplitude than $SN_1$ and saddle nodes are labelled $SN_1$, $SN_2$, $SN_3$, as shown in figure 7(i) for $Le = 20$, $Pr = 0.023 \approx Pr_{cusp}$.

Continuing to reduce $Pr$ across region (2) (from $Pr_{cusp}$ to $Pr_c$), the Rayleigh number associated with $SN_2$ increases so that it reaches the supercritical region before $Pr = Pr_c$. During this transition, the saddle node with smallest amplitude ($SN_3$ for $Le \lesssim 15$; $SN_1$ for $Le \gtrsim 19$) moves to larger Rayleigh numbers but with decreasing amplitude until it collides with the primary bifurcation at ${Pr = Pr_c}$ and ${Ra = Ra_c}$, where the primary bifurcation changes from subcritical to supercritical. This process is highlighted in the insets of figure 6 and results in the primary branch possessing only two saddle nodes in the supercritical regime ($Pr < Pr_c$).

The locations of the remaining two saddle nodes go toward larger $Ra$ as $Pr$ decreases and are found in the supercritical region ($Ra > Ra_c$) in region (4), as shown in figure 5. It is therefore clear that multiple steady convection states can exist for the same parameter values near the onset of convection, regardless of the criticality of the primary bifurcation. This result extends earlier observations on the number of saddle-node bifurcations occurring along the primary branch in related systems (Tsitverblit & Kit Reference Tsitverblit and Kit1993).

More insight into these results can be obtained by representing, as in figure 8, the location of the saddle nodes for various Lewis numbers as a function of the reduced Prandtl number $Pr/Pr_c$ and combined parameter $Ra|Le-1|$. These reduced parameters allow us to identify the location where the criticality of the primary bifurcation changes as the single coordinate point: $Pr/Pr_c = 1$, $Ra|Le-1|\approx 6509$.

Figure 8.

Saddle-node locations for $Le = 5$ (green), $Le=11$ (purple), $Le=20$ (blue) and $Le=50$ (red). The black dashed line marks the location of the primary bifurcation and the black cross marks the codimension two point where the criticality of the primary bifurcation changes, at $Pr=Pr_c$. The black dotted line represents the relationship (4.2). Saddle nodes are marked by circles for $SN_1$, asterisks for $SN_2$ and triangles for $SN_3$.

Figure 8 shows that, for $Pr < Pr_c$ and the chosen values of the Lewis number, the location of the first supercritical saddle node $SN_2$ can be approximated by

(4.2)\begin{equation} Ra_{SN_2} \approx \frac{6460}{|Le-1|}\left(\frac{Pr}{Pr_c}\right)^{{-}0.24}. \end{equation}

For $Pr < 10^{-2}$ (not shown), the location of saddle node $SN_2$ deviates from the relation above, indicating a potentially different asymptotic regime. These results also illustrate the large $Pr$ behaviour of the subcritical saddle node $SN_1$: $Ra_{SN_1}|Le - 1|$ tends to a constant as the Prandtl number tends to infinity. This constant increases with $Le$ and saturates for large values of the Lewis number. These results echo those obtained in doubly diffusive convection in a two-dimensional vertical porous enclosure, where Mamou, Vasseur & Bilgen (Reference Mamou, Vasseur and Bilgen1998) used a parallel flow approximation to demonstrate that the Rayleigh number at which the subcritical saddle node occurs is proportional to $1/(1-Le)$ for large enough Lewis numbers.

4.2. Solution profiles

Despite the different scenarios obtained at different values of the Prandtl number (see figure 5), the steady convection states undergo similar structural changes along their branch, as evidenced in figure 9 for $Le = 11$ and $Pr = 1$, 0.032, 0.01 and $0.005$. The streamfunction profiles are similar near the primary bifurcation regardless of the value of the Prandtl number (see second column of figure 9), which is in agreement with the linear stability results from figure 3(c). Moving along the branches in the direction of increasing energy, the first change that we observe is the strengthening of the anticlockwise roll, where fluid near the hotter wall moves upwards. This occurs in both the subcritical and the supercritical regimes, as can be seen in the third column of figure 9. Continuing the branches to the large-amplitude saddle node and beyond, the amplitude of the weaker roll decreases, leaving room for the stronger roll to straighten. At large enough amplitude, an anticlockwise roll occupies the domain, irrespective of the value of $Pr$. Its amplitude grows as the upper branch is followed to larger values of $Ra$, where the Prandtl number starts to impact the flow: the roll occupies a smaller area at lower values of the Prandtl number, as seen within the final column of figure 9. This resembles the fly-wheel convection, with nearly circular streamlines, seen in low Prandtl Rayleigh–Bénard convection as studied by Clever & Busse (Reference Clever and Busse1981).

Figure 9.

Streamfunctions of the steady states on the primary branch when ${Le = 11}$ for different values of the Prandtl number: top row ${Pr = 1}$, second row ${Pr = 0.032}$, third row ${Pr = 0.01}$ and bottom row ${Pr = 0.005}$. The left column shows the respective bifurcation diagrams and indicates with a cross the solutions that have been represented in the subsequent panels. Black (grey, dotted) contours indicate positive (negative, zero) values of the streamfunction. Contour intervals: first column, top two rows – $10^{-4}$; first column, third row – $10^{-5}$; first column, bottom row – $2\times 10^{-5}$; second column, top two rows – $0.02$; second column, bottom two rows – $0.01$; third column – $0.02$; fourth and fifth columns – $0.05$.

To characterise these observations in more detail, figure 10 reports the horizontal velocity profiles observed on the upper branch for $Pr = 1$, $0.1$, $0.032$ and $0.005$. The decrease in roll size is apparent when $Pr$ is decreased. This is particularly evident for $Pr = 0.005$, where the horizontal velocity remains small except within the range $0.6 \lesssim z \lesssim 1.9$, in such a way that the roll only occupies about half of the domain's extent. Figure 10(a) additionally shows the transition to these states from the large rolls observed at $O(1)$ Prandtl numbers. For $Pr = 1$, the maximum horizontal velocity is achieved far from the centre of the roll, at $z \approx 0.44,\,2.04$, producing a region of strong shear between the rolls and gentle quasi-linear velocity variations inside the rolls. As $Pr$ is lowered, these maxima move towards the centre of the roll by initially becoming less pronounced and creating flatter extrema (see figure 10c), followed by the emergence of peaks at $z = 1$ and $z\approx 1.5$. The maximum horizontal velocity does not change significantly within this range of Prandtl number values in such a way that the low $Pr$ rolls represent narrow regions of strong shear surrounded by low-amplitude flow.

Figure 10.

Horizontal velocity and streamfunction of solutions from the upper segment of the primary branch at $Ra = 700$ for $Pr = 1$, 0.1, 0.032, $0.005$ and $Le = 11$ represented via (a) the midline horizontal velocity ($u(x=0.5,z)$) and streamfunction contours plots for (b) $Pr = 1$, (c) $Pr = 0.1$, (d) $Pr = 0.032$ and (e) $Pr = 0.005$ with contour intervals $0.1$.

4.3. Stability of the nonlinear states

The stability of states on the primary branch is controlled by two eigenmodes: an amplitude mode that preserves the $S_{\varDelta }$ symmetry of the system and a drift mode that breaks the $S_{\varDelta }$ symmetry. The translation mode, associated with vertical translations due to the symmetry $T_{\delta }$, remains marginal along the branch and none of the other eigenmodes become destabilising over the range of parameters considered.

Close to the onset of convection, the amplitude mode is initially destabilising when the bifurcation is subcritical ($Pr > Pr_c$), whereas it is stabilising when the bifurcation is supercritical ($Pr < Pr_c$). This mode subsequently changes stability at successive saddle nodes. In particular, it becomes stabilising at saddle nodes $SN_1$ and $SN_3$, where the branch turns towards higher $Ra$, but becomes destabilising at $SN_2$, where the branch turns towards lower $Ra$. As a result, the upper branches of steady convection states are always stable to amplitude perturbations for all $Le$ and $Pr$.

The drift mode is stabilising near the primary bifurcation at $Ra = Ra_c$ for all $Pr$, but becomes destabilising at a drift-pitchfork bifurcation further along the branch at $Ra = Ra_d$, whose location depends upon both $Le$ and $Pr$, as can be seen in figure 7. The marginal mode is identical to the translation mode at this bifurcation and its destabilisation leads to a pair of branches of travelling wave solutions, as shown in figure 11(a) for $Pr = 0.1$ and $Le = 11$. Close to their onset, these states take the form of a single large-amplitude convection roll (see figure 11c) that slowly drifts either upwards or downwards. As these branches are followed beyond the drift bifurcation, an asymmetric streaming flow strengthens while the convection roll weakens and moves toward the wall where the streaming flow is the weakest. This transition is shown from figure 11(b) at $Ra = 630$ to figure 11(d) at $Ra= 700$. At the same time, the drift speed increases at a rate approximately proportional to $\sqrt {Ra-Ra_d}$, as shown in figure 11(e). This result extends the findings obtained for $Le = 1.2$, $Pr = 1$ by Xin et al. (Reference Xin, Le Quéré and Tuckerman1998) to a wider range of parameter values.

Figure 11.

Drift bifurcation and downward-travelling waves for $Pr = 0.1$, $Le =11$, for which $Ra_d \approx 638$. (a) Bifurcation diagram showing the kinetic energy $E$ as a function of the Rayleigh number $Ra$ for steady states and travelling waves. Thick lines indicate stable solutions, thin lines indicate solutions unstable to amplitude perturbations and dashed lines indicate solutions unstable to drift. The drift bifurcation is shown by the open circle. (b) Stable convection state at $Ra = 630$ shown by contours of its streamfunction with intervals $0.1$ (first red cross in (a)). Further panels show similar representations of stable travelling waves at: (c) $Ra = 645$ and (d) $Ra = 700$. (e) Squared drift speed along the stable branch as a function of the Rayleigh number. The dotted line shows the fitting law: $v_d \approx 0.12\sqrt {Ra - 640}$.

The stability of the travelling waves is determined by the location of the drift bifurcation: these states are initially stable when the bifurcation occurs on the upper branch of steady convection states, whereas they are unstable when the bifurcation occurs along the lower branch. Both cases can be achieved for a given $Le$ when $Pr$ is varied, as figure 7 illustrates for selected values of the Prandtl number with $Le = 11$ and $Le = 20$. For large values of the Prandtl number, the drift bifurcation occurs on the upper branch at large Rayleigh numbers. This location moves closer to the saddle node with decreasing Prandtl numbers so that the two coincide at $Pr = Pr^*$ and $Ra = Ra^*$. For $Le = 11$, we found that $Pr^*\approx 0.042$ and $Ra^* \approx 614.9$ (see figure 7(c) for a bifurcation diagram at similar values of the parameters). For smaller values of the Prandtl number, the drift bifurcation occurs along the lower branch of convection states and at a value of the Rayleigh number that increases as $Pr$ is decreased. For all the parameter values tested, this bifurcation was found to occur at larger amplitude than saddle node $SN_2$ and, consequently, the small-amplitude steady convection states remain stable to drift.

4.4. Dynamical attractors

The temporal dynamics of the system changes as the drift bifurcation passes below the subcritical saddle node since all the steady convection states from the upper branch and the travelling wave states are destabilised in the process. Many initial conditions will consequently decay towards the conduction state at low $Pr$ and $Ra$. This decay is not possible when the conduction state is unstable for $Ra > Ra_c$, where we find that the dynamics converge on time-dependent states.

To understand how this behaviour arises, we unfold the saddle-node-pitchfork normal form near the codimension two point $(Ra^*,Pr^*)$ where the drift bifurcation and saddle node coincide. This unfolding takes the form (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983)

(4.3)$$\begin{gather} \dot{x} ={-}\mu_1 x + b_1 xz, \end{gather}$$

(4.4)$$\begin{gather}\dot{z} = \mu_2 - x^2 - z^2 + b_2 z^3, \end{gather}$$

where $x$ represents the extent to which the state drifts, $z$ represents the amplitude of the convection states, $b_1>0$, $b_2<0$, and $\mu _1$ and $\mu _2$ are two unfolding parameters that are introduced to respectively represent the deviations ${Pr-Pr^*}$ and ${Ra-Ra^*}$.

When $\mu _1 = \mu _2 = 0$, the trivial state, $(x,z) = (0,0)$, undergoes a codimension two bifurcation. One of five phase portraits is observed in the vicinity of this bifurcation and these are shown in figure 12(a). In addition to the steady states previously discussed, the unfolding reveals the presence of periodic orbits. Relating the unfolding back to doubly diffusive convection, these correspond to relative periodic orbits consisting of drifting states that originate either from a travelling wave undergoing a Hopf bifurcation or from a global bifurcation where two steady convection states connect heteroclinically.

Figure 12.

(a) Unfolding near the codimension two saddle-node-pitchfork bifurcation at $\mu _1=\mu _2=0$ given by system (4.3), (4.4), after Guckenheimer & Holmes (Reference Guckenheimer and Holmes1983). The different phase portraits are classified in five different regions labelled using Roman numerals and accompanied with a sketch of the corresponding phase space. In each of these sketches, the fixed points on the vertical line represent steady convection states. The vertical (respectively horizontal) direction is the eigendirection related to the amplitude (respectively drift) mode. (b) Analogy with the doubly diffusive convection problem is made by replacing $\mu _1$ by $Pr-Pr^*$ and $\mu _2$ by $Ra-Ra^*$ and regions of the $(Ra,Pr)$ parameter space are shown as a function of the observed temporal behaviour for $Le = 11$. (c) Magnification of panel (b) near $(Ra^*,Pr^*)$. Arrows indicate the values of $Pr$ used to produce the bifurcation diagrams in figure 13. In the panels, the bifurcations are represented by: black, red and blue solid lines (saddle nodes), blue dotted lines (drift bifurcation), red dotted lines (Hopf bifurcation), red dot-dashed lines (heteroclinic connection) and in (b,c), the vertical dashed lines (primary stationary bifurcation of the conduction state).

Although the normal form (4.3), (4.4) only formally represents the dynamics of the full system close to the codimension two point, each of the regions shown in figure 12 continues to be observed an appreciable distance away from this point. Figures 12(b) and 12(c) illustrate the extent of the corresponding regions in the doubly diffusive system when $Le = 11$ and we anticipate that similar results will hold for other values of the Lewis number. In this figure, the regions have been subdivided according to the types of stable attracting states that they display. The subdivisions occur owing to the instability of the conduction state at $Ra_c$ and the creation of a pair of saddle nodes at $(Ra_{cusp},Pr_{cusp})$, which enrich the previous unfolding. The resulting subregions, together with their associated attracting states, are summarised in table 3 and on the bifurcation diagrams in figure 13. As $Pr$ varies, the system admits one of seven qualitatively distinct bifurcation diagrams. Six of these are presented in figure 13, which also indicate the range of kinetic energies over each relative periodic orbit attained via time stepping. The seventh type of bifurcation diagram, where the primary branch lies entirely within the supercritical regime, is not shown but possesses similar features to that seen for $Pr = 0.02$ in figure 13(f) including stable small-amplitude steady convection states and relative periodic orbits.

Table 3.

Stability of the observed doubly diffusive states within each region of the parameter space from figure 12. The naming convention used is as follows: $O$, conduction state; $SOC_s$, small-amplitude stationary overturning convection; $SOC_l$, large-amplitude stationary overturning convection; $TW$, travelling wave; and $PO$, relative periodic orbit. The regions ${\rm Oa},\ldots, {\rm Vb}$ refer to the regions introduced in figure 12.

Figure 13.

Bifurcation diagrams showing the primary branch and other stable attracting states for $Le = 11$, and (a) $Pr =1$, (b) $Pr = 0.1$, (c) $Pr = 0.043$, (d) $Pr = 0.04$, (e) $Pr = 0.032$ and (f) $Pr = 0.02$. The solid circles mark the saddle nodes and open circles indicate where the drift bifurcation occurs. Thick (thin) lines represent states stable (unstable) to the amplitude mode whilst solid (dashed) lines show those stable (unstable) to the drift mode. Thick blue lines indicate the minimal and maximal energies achieved in the stable limit cycle, which starts in a Hopf bifurcation in (c) and in a heteroclinic bifurcation in (d–f). The unstable branches of travelling waves are not shown.

The three most relevant stable attracting states close to the primary bifurcation at high $Pr$ ($Pr>Pr^*$ here) are: the conduction state ($O$), the large-amplitude steady convection states ($SOC_l$) and the travelling wave states ($TW$). Below the onset of convection (region Oa), all initial conditions decay towards the first of these. In region Ia, above subcritical onset but before the drift instability, initial conditions converge towards $SOC_l$, as evidenced by the energy–time and drift speed–time plot in figure 14(a). Increasing $Ra$ beyond the drift instability into region IIa, $SOC_l$ is now unstable and the flow converges towards $TW$. Figure 14(b) shows that the former state may still be observed in the temporal dynamics as the initial condition first rapidly changes amplitude to approach $SOC_l$ before it builds vertical drift and converges to $TW$.

Figure 14.

Energy–time (black) and drift speed–time (red) plots illustrating regions I–V in figure 12 with $Le = 11$. In each case, the initial state was the large-amplitude convection state at $Ra = 700$ for $Pr = 0.1$ that was perturbed in the direction of its unstable drift eigenmode. States approached during the trajectory are labelled as follows: (a) region Ia, convergence to $SOC_l$ when $Pr = 0.1$ and $Ra = 630$; (b) region IIb, convergence to $TW$ when $Pr = 0.1$ and $Ra = 660$; (c) region IIIb, convergence to $PO$ when $Pr = 0.032$ and $Ra = 660$, (d) region IVa, convergence to $O$ when $Pr = 0.032$ and $Ra = 630$; (e) region IVc, convergence to $SOC_s$ when $Pr = 0.02$ and $Ra = 700$; and (f) region Va, convergence to $O$ when $Pr = 0.032$ and $Ra = 620$.

The stable branch of travelling waves destabilises in a supercritical Hopf bifurcation that leads to a stable relative periodic orbit, as shown in figure 13 for $Pr = 0.043$. Figures 15(a)–15(e) depict such an orbit shortly after the bifurcation at $Ra = 650$ and $Pr = 0.043$, where we see that the states exhibit small oscillations about a drifting state. The Hopf bifurcation moves towards lower Rayleigh numbers as $Pr$ approaches $Pr^*$ from above, which reduces the extent over which stable $TW$ are found. This continues until $Pr = Pr^*$, when stable $TW$ cease to exist and the relative periodic orbit bifurcates directly from the codimension two bifurcation at the saddle node.

Figure 15.

Temporal evolution of downward-travelling states across one cycle of two relative periodic orbits at: (a–e) $Ra = 650$ with $Pr = 0.043$ and $Le = 11$ and (f–j) $Ra = 660$ with $Pr = 0.032$ and $Le = 11$ (as in figure 14c). (a,f) Anticlockwise trajectory of the periodic orbit in drift speed–energy phase space. Blue dots in ( f) mark the conduction and steady convection states. (b,g) Energy–time (top) and drift speed–time (bottom) plots. (c–e) Streamfunctions of states along the orbit in (a) at (c) $t = 0$, (d) $t = 8$ and (e) $t=20$ with contour intervals $0.1$. (h–j) Streamfunctions of states along the orbit in (f) at (h) $t = 26$, (i) $t = 52$ and (j) $t=68$, with contour intervals $0.05$. The streamfunctions have been translated vertically for better visual representation.

Upon further decreasing of the Prandtl number, so that the drift bifurcation occurs on the lower branch of steady convection, the system admits neither stable $SOC_l$ nor stable $TW$. Instead, the bifurcation diagrams are similar to that shown for $Pr = 0.04$ in figure 13(d), where a branch of unstable $TW$ extends from the drift bifurcation towards higher Rayleigh numbers and stable relative periodic orbits exist after a global bifurcation, where the stable manifold of $SOC_l$ connects heteroclinically with the unstable manifold of the convection state on the lower branch and vice versa.

The lack of stability of the nonlinear states before the heteroclinic connection lead all initial conditions to decay down to the conduction state in regions IVa and Va. Figures 14(d) and 14(f) illustrate this tendency for $Pr = 0.032$ when $Ra = 630$ and $Ra = 620$, respectively. In both cases, the amplitude of the initially imposed roll rapidly decreases to approach that of $SOC_l$ rolls. Afterwards, the drift speed of the state increases, as $SOC_l$ is unstable to drift, and reaches a maximum around $t \approx 50$. The drift speed subsequently decays down to zero, due to the instability of $TW$, and the time-dependent state converges on the conduction state, which is the only stable attractor in these regions.

Beyond the heteroclinic connection, initial conditions tend to converge towards the relative periodic orbit, as they invariably do in region IIIb, where the conduction state is unstable. Figure 14(c) illustrates this convergence starting from a large-amplitude roll with $Pr = 0.032$ and $Ra = 660$ perturbed in the direction of its unstable drift eigenmode. A single cycle of this orbit is shown in further detail in figures 15(f)–15(j). This relative periodic orbit cycles between the three states: $SOC_l$, $TW$ and a steady small-amplitude convection state, in the following manner. The first stage of the orbit, from $15\lesssim t \lesssim 40$, resembles the temporal behaviour seen in region IIa (figure 14b), where the solution remains close to $SOC_l$ in profile (figure 15h) while the drift speed slowly increases in magnitude. Following this, between $t\approx 40$ and $t\approx 54$, the drift speed and kinetic energy rapidly increase as the profile of the state exhibits properties of the travelling wave ($TW$) solution (figure 15i). Between $t\approx 54$ and $t \approx 68$, both the drift speed and kinetic energy decrease as the state approaches a small-amplitude, non-drifting convection state with inclined rolls (figure 15j). The final stage of this orbit is the transition from the small-amplitude back to large-amplitude steady convection, which is indicated by the monotonic increase in kinetic energy while maintaining $v_d \approx 0$ for $t \gtrsim 70$ and $t \lesssim 15$ in figure 15(g).

The heteroclinic connection leading to these orbits moves towards higher Rayleigh numbers as $Pr$ decreases and coincides with $SN_2$ for $Pr \lesssim 0.032$ (see figure 13e,f). This suggests that a saddle-node infinite period (SNIPER) bifurcation explains the origin of the relative periodic orbits at low Prandtl and high Rayleigh numbers. However, by considering various properties of the relative periodic orbits for $Pr = 0.032$ and $Le = 11$ as $Ra$ approaches $Ra_{SN_2}$ from above (figure 16), we additionally find that a gluing bifurcation occurs in the vicinity of the SNIPER bifurcation.

Figure 16.

Relative periodic orbits for $Le = 11$, $Pr = 0.032$ where $Ra_{SN_2} \approx 650.82$. (a) Trajectories in $(v_d,E)$ phase space for $Ra = 650.85$ (blue), $Ra=660$ (red) and $Ra=700$ (black). For $Ra = 700$ and $Ra = 660$, a pair of relative periodic orbits associated with either negative or positive drift velocity are shown, while for $Ra = 650.85$, a single periodic orbit with alternating negative and positive drift velocities is shown. (b) Period $t_P$ of orbits for selected $Ra>Ra_{SN_2}$. The red dashed line shows that approximately $t_P \propto (Ra-Ra_{SN_2})^{-0.56}$. (c) Energy–time plots for $Ra = 700$ (top), $Ra = 660$ (middle) and $Ra = 650.85$ (bottom).

At large Rayleigh numbers, a pair of relative periodic orbits with states drifting either upwards or downwards are related by the reflection symmetry. The maximal energy and drift speed attained along these orbits decrease with decreasing Rayleigh number, and the trajectories approach the stable and unstable manifolds of $SOC_l$, as seen in figure 16(a). This leads to the two relative periodic orbits connecting in a gluing bifurcation around $Ra \approx 652$ so that the trajectories become a single periodic orbit where states alternately drift in opposite directions. This is reminiscent of the pulsating waves seen in nonlinear magnetoconvection (Matthews et al. Reference Matthews, Proctor, Rucklidge and Weiss1993).

The resulting single periodic orbit persists until $Ra_{SN_2}$, where it terminates in the SNIPER bifurcation. This is evidenced by the period of a single loop of the orbit scaling like $t_P \propto (Ra-Ra_{SN_2})^{-0.56}$ as $SN_2$ is approached, which is close to the expected $t_P \approx |Ra-Ra_{SN_2}|^{-0.5}$ scaling. The energy–time plots in figure 16(c) illustrate that the predominant increase in duration occurs near the small-amplitude steady convection state as the orbit approaches the steady state at $SN_2$ in phase space. We also find that the time spent near $SOC_l$ increases, whilst the time where the state has large drift speed remains small, implying that the global bifurcation is due to the collision of the periodic orbit with the stable manifold of $SOC_l$.

The final attracting state that the flow may converge to is $SOC_s$, as figure 14(e) illustrates for $Pr = 0.02$ and $Ra = 700$. This is possible for $Pr < Pr_{cusp}$ in the supercritical regions Ob, IVc and Vb, where it is the only stable attracting state, and in the subcritical region IVb, where convergence towards the stable conduction state is also possible.