Bifurcation analysis of two-dimensional Rayleigh–Bénard convection using deflation

Steady states of Rayleigh–Bénard convection are computed using a recent numerical technique called deflation [20]. This algorithm is based on Newton’s method and can compute multiple solutions to a nonlinear system of equations $F(\phi,\lambda)=0$, where $\phi$ is the solution and $\lambda\in\mathbb{R}$ is a bifurcation parameter. First, fix $\lambda\in\mathbb{R}$ and assume a solution $\phi_{1}$ to $F(\cdot,\lambda)$ has been previously discovered by Newton’s method. One can then construct a deflated problem

such that $G(\phi,\lambda)$ does not converge to zero as $\phi\to\phi_{1}$. The term $\mathcal{M}(\phi,\phi_{1}):=\|\phi-\phi_{1}\|^{-2}+1$ is called a deflation operator and ensures that Newton’s method applied again to $G$ will not converge to a previously computed solution. By adding one to the expression of the deflation operator, we impose that the new problem $G$ behaves similarly to $F$ away from the root $\phi_{1}$. This process can be repeated to obtain a set of solutions $S(\lambda)$ to a nonlinear problem $F(\phi,\lambda)=0$ at the bifurcation parameter $\lambda$. Deflated continuation is the combination of this idea with continuation. The previously discovered solutions are continued to a new bifurcation parameter $\lambda-\Delta\lambda$: each $\phi\in S(\lambda)$ is used as an initial guess for the deflation procedure in order to construct a new set of solutions $S(\lambda-\Delta\lambda)$ at $\lambda-\Delta\lambda$.

We apply this algorithm repeatedly to discover multiple solutions to the steady state Rayleigh–Bénard problem:

where $\phi=({\bf{u}},p,T)$ and $Ra$ is the bifurcation parameter. This equation is discretised with Taylor–Hood finite elements for the velocity and pressure (piecewise biquadratic and piecewise bilinear respectively), together with piecewise bilinear polynomials for the temperature, using the Firedrake finite element library [38]. The unit square domain is represented by a mesh with $50\times 50$ square cells to preserve the symmetries of the problem and avoid spurious symmetry breaking solutions. A coarse discretization is needed due to the complexity of the problem and number of distinct solutions discovered. However, we did not observe significant differences in selected experiments where we refined the mesh.

Deflated continuation requires two different computational tasks. The first is the continuation of a previously discovered branch by the continuation step, which usually results in the convergence of the Newton solver in an average of three iterations. The convergence criterion is defined to be when the Euclidean $l^{2}$-norm of the discretised PDE residual is below $10^{-8}$. The second task is the deflation step to discover new branches. Most of the time, the solver does not converge, which means that it does not discover any new branch within the maximum number of Newton iterations allowed, which is set to one hundred. Typically, when the deflation step is successful, it requires around $80$ iterations to find a new solution. Since the deflation step is employed for each of the steady states in the different branches, this yields at least hundreds of thousands of Newton iterations in total, and a few weeks of calculation on a desktop with 32 cores. The linear systems resulting from Newton’s method are solved using the sparse direct solver MUMPS [39].

A crucial question in the deflation technique is the choice of the deflation operator $\mathcal{M}$ because it can affect the convergence of Newton’s method. We use the following operator,

where $\|\cdot\|_{2}$ denotes the $L^{2}$-norm. Note that the deflation operator $\mathcal{M}$ does not depend on the pressure $p$ because we wish to deflate all the solutions of the form $(u_{1},p_{1}+c,T_{1})$, where $c\in\mathbb{R}$, to avoid discovering solutions trivially related to known ones.

We perform bifurcation analysis of Rayleigh–Bénard convection and compute solutions to the nonlinear system of Eq. 8 for $0\leq Ra\leq 10^{5}$ using deflated continuation for $Ra=10^{5},10^{5}-\Delta R_{a},\ldots,0$. Here, $\Delta Ra$ denotes the continuation step size in the bifurcation parameter, which is chosen as $\Delta Ra=100$. In a previous work on the computation of solutions to the 3D nonlinear Schrödinger equation [26], we found that having good initial guesses (e.g. solutions of the linearised problem) for the initial deflation steps facilitates the convergence of Newton’s method and leads to more complex and interesting states. However, contrary to [26] we do not have the analytic expressions of the linear solutions to the Rayleigh–Bénard convection due to the no-slip boundary conditions. Therefore, we numerically solve an eigenvalue problem (see Section III.2) to obtain the first ten unstable eigenmodes $({\bf{u}}_{1},T_{1}),\ldots,({\bf{u}}_{10},T_{10})$ linearised around the conducting steady state Eq. 4 (see Fig. 2). Then, the sums of the conducting state with the respective normalized perturbations: $({\bf{u}}_{1},1-z+T_{1}),\ldots,({\bf{u}}_{10},1-z+T_{10})$ are used as initial guesses for deflated continuation.

The stability analysis of a given steady state $({\bf{u}}_{0},T_{0})$ is performed by considering the perturbation ansatz

where the velocity perturbation ${\bf{v}}\ll 1$ and the temperature perturbation $\theta\ll 1$ have an eigenvalue $\lambda$, whose real part is the growth rate and imaginary part is the frequency of the corresponding eigenvector. When at least one of the eigenvalues has a positive real part $\mathcal{R}(\lambda)>0$, we can have two types of instabilities: to a real eigenmode, occurring when the imaginary part of the eigenvalue $\mathcal{I}(\lambda)=0$, and to a complex eigenmode, occurring when $\mathcal{I}(\lambda)\neq 0$. If all the eigenvalues have negative real parts, then the steady state is stable. Inserting Eq. 10 into Eq. 1 and considering the linearised system of equations, we obtain the following generalised eigenvalue problem at leading order

where $I$ is the relevant identity operator and $A$ is the linear operator defined by

The eigenvalue problem described in Eq. 11 is solved with a Krylov–Schur method [40] using the SLEPc library [41].

To obtain the critical value of the Rayleigh number $Ra_{c}$ at which the conducting state Eq. 4 becomes unstable, we modify the eigenvalue problem in Eq. 11. Inserting Eq. 4 into Eq. 11 we obtain the following generalized eigenvalue problem

A steady state becomes unstable when $\lambda=0$ and this corresponds to the critical Rayleigh number $Ra_{c}$ which is the smallest $Ra$ satisfying Eq. 12 for $\lambda=0$. So, to obtain the critical value of the Rayleigh number for different base states we have reformulated the problem to a generalised eigenvalue problem for $Ra_{c}$ as follows

In the following sections this linear analysis will be used to study the stability of the conducting state Eq. 4 and of the nonlinear steady states that we obtain from deflated continuation.

In this section, we analyse the stability of the conducting steady state defined by Eq. 4. We found that the first instability of this state arises at $Ra_{c}^{(1)}:=Ra_{c}\approx 2586$. It is well known [42] that the critical Rayleigh number is approximately equal to $Ra_{c}^{*}\approx 1707.762$ for a domain with periodic side boundaries and no-slip top and bottom boundaries. This difference in the value of $Ra_{c}$ is due to the effect of the side walls on the instability. This effect should decrease when the length of the domain, $L$, becomes large such that the behaviour predicted for an unbounded domain is recovered. To demonstrate this, we systematically increase the aspect ratio

Figure 1 (left) clearly shows that for $\Gamma\gg 1$ the $Ra_{c}^{(1)}$ converges to $Ra_{c}^{*}$.

In addition, for aspect ratios $\Gamma<1$ we observe $Ra_{c}\propto\Gamma^{-3.88}$. The critical Rayleigh number for Rayleigh–Bénard convection with stress-free boundary conditions is $Ra_{c}^{*}(k)=(\pi^{2}+k^{2})^{3}/k^{2}$ [43], where $k$ is the horizontal wavenumber. For $k\gg 1$ the critical Rayleigh number scales like $Ra_{c}^{*}(k)\propto k^{4}\propto\Gamma^{-4}$ using Eq. 14. This scaling relation is very close to the observed one, which implies that the exponent is not affected much by the no-slip boundary conditions. The exponent we find is in agreement with [44].

The case on which we focus for the rest of the paper is the square cell ($\Gamma=1$), which is denoted with a red dot in Fig. 1 (left). We first study the linear stability of the conducting state Eq. 4. In the range $0\leq Ra\leq 10^{5}$ we observe 10 supercritical stationary bifurcations that arise from the conducting state. Figure 1 (right) shows the growth rates of the unstable eigenmodes as a function of $Ra$. These bifurcations occur when $\mathcal{R}(\lambda)=0$. The critical values of the Rayleigh number at the onsets are listed in Table 1. Fig. 1 (right) demonstrates that the growth rates of the eigenmodes vary with Rayleigh number such that the curves can cross. The eigenmode which bifurcated from the conducting state at $Ra_{c}^{(2)}=6746$ ends up being the most unstable as $Ra\to 10^{5}$.

The velocity magnitude and temperature fields of the eigenmodes from these primary bifurcations of the flow are shown in Fig. 2 with the corresponding critical Rayleigh numbers indicated on the top of each flow pattern. In Table 1 we characterise the flow patterns by the number of rolls $(m_{x},m_{z})$, where $m_{x}$ and $m_{z}$ are the number of rolls in the temperature field in the horizontal and vertical direction, respectively. Note that due to the aspect ratio $\Gamma=1$ of the domain the flow pattern of the first bifurcation at $Ra_{c}^{(1)}=2586$ is just a single convection roll (see Fig. 2) according to the classical result [33]. In the following sections, we first analyse the branches of solutions that are created from some of these primary bifurcations before moving to disconnected branches.

The combination of deflation and the initialisation strategy presented in Section III.1 identified $43$ branches for the two-dimensional steady-state Rayleigh–Bénard convection at $Ra=10^{5}$, leading to $129$ branches when reflections are included. Here, we analyse the evolution and the linear stability of some of the branches arising from the first four excited states (see Fig. 2) as well as two disconnected branches. Moreover, we discard the mirror symmetric solutions with respect to the $x=1/2$ and $z=1/2$ axes (see symmetries in Section II). Our diagnostics for the bifurcation diagrams, representing the evolution of the steady states as a function of $Ra$, are the kinetic energy $\|{\bf{u}}\|_{2}^{2}$, the potential energy $\|T\|_{2}^{2}$, and the Nusselt number $Nu$, which are defined as

The Nusselt number characterises the efficiency of convective heat transfer and is given by the ratio of the total heat transfer (i.e. both advective and diffusive) to the conductive heat transfer $Nu\equiv HL/\kappa d=\langle{wT}\rangle-\kappa\partial_{z}\langle{T}\rangle=\langle{(\nabla T)^{2}}\rangle$, where $H$ is the heat flux and $\langle{\cdot}\rangle$ stands for the area average. A Nusselt number of one represents heat transfer by pure conduction. Here, the Nusselt number is essentially defined as the dissipation rate of the temperature variance (see second equality). This is obtained by integrating the equation for potential energy over the entire area [45].

The bifurcation diagrams are presented in Fig. 3, where the numbers denote the different branches of steady solutions, indicating the values of kinetic energy, potential energy and Nusselt number that correspond to each branch on each diagram.

The first steady solutions that we obtain from deflated continuation as $Ra$ is varied are four branches that arise from the first four bifurcations of the conducting steady state: branch (1) at $Ra_{c}^{(1)}=2586$, branch (2) at $Ra_{c}^{(2)}=6746$, branch (3) at $Ra_{c}^{(3)}=19655$ and branch (4) at $Ra_{c}^{(4)}=23346$. Moreover a secondary branch (10) bifurcates at $Ra\approx 13550$ from branch (2), which originates from the second instability of the conducting state. A similar behaviour is observed for the fourth instability, where branch (6) bifurcates at $Ra\approx 27300$ from branch (4). We also focus on two disconnected branches (11) and (8) obtained directly by deflated continuation. In Fig. 3(a), branch (11) is close to branch (12). However, note that the bifurcation diagram using the potential energy (see Fig. 3(b)) shows that branch (11) does not bifurcate from branch (12). On the other hand, branch (8) bifurcates from branch (11) at $Ra\approx 37000$ (see Fig. 3). It is interesting to point out that in the range $0\leq Ra\leq 10^{5}$ branch (1) has the highest kinetic energy while branch (6) has the highest potential energy (see Fig. 3(a) and 3(b)). The steady states that are most effective in convecting heat transfer for $0\leq Ra\lesssim 3\times 10^{4}$ are in branch (1) from the first instability of the conducting state while for $3\times 10^{4}\lesssim Ra\leq 9.5\times 10^{4}$ are in the branch (2), and in branch (11) for $Ra\gtrsim 9.5\times 10^{4}$ (see Fig. 3(c)).

In this section we analyse the states that emanate from the conducting steady state and their evolution as a function of the bifurcation parameter $Ra$. Our analysis involves the detailed evolution of the velocity magnitude and the temperature on the branches of the bifurcation diagrams using the kinetic and potential energy, respectively. In addition, we present the largest growth rates and corresponding frequencies from the linear stability analysis we have performed on the steady states.

Figure 4 shows results of the aforementioned analysis for branch (1), which arises at $Ra_{c}^{(1)}=2586$ from the first eigenmode (see Fig. 2). This bifurcation breaks the $x$-reflection symmetry but the branch has the Boussinesq symmetry. As $Ra$ increases, the magnitude of the velocity field evolves from a circular to a bent convection roll at $Ra\approx 4\times 10^{4}$ (see Fig. 4(a)). On the other hand, the $z$-shaped interface of the temperature field is sheared (see Fig. 4(b)), enhancing the heat transfer in the convection cell (see also branch (1) in Fig. 3(c)). Fig. 4(c) demonstrates that the solutions on this branch remain stable ($\mathcal{R}(\lambda)<0$) with an almost constant growth rate over the whole range of Rayleigh numbers we consider. The frequency corresponding to the largest growth rate is zero from the stationary bifurcation at $Ra_{c}^{(1)}=2586$ to the $Ra=10^{5}$ (see Fig. 4(d)), while the next eigenvalue has a negative real part and a large imaginary part.

The second eigenmode introduced in Fig. 2 gives birth to branch (2) at $Ra_{c}^{(2)}=6746$, which is depicted in Fig. 5(a)-(b). The bifurcation is transcritical and preserves the symmetries of the eigenmode in the velocity and temperature profiles. Then, the largest growth rate in Fig. 5(c) indicates that the solutions on this branch are stable over the interval $Ra\in[13500,85500]$ and unstable outside. A bifurcation occurs at $Ra\approx 13500$ when $\mathcal{R}(\lambda)=\mathcal{I}(\lambda)=0$, which is analysed later and illustrated in Fig. 6. The bifurcation from branch (2) to branch (10) is accompanied by a loss of Boussinesq symmetry. Another bifurcation is observed at $Ra\approx 85500$ and has been obtained by deflated continuation.

Figure 6(a) demonstrates that the two symmetries of the second eigenmode (see Section II) are rapidly broken as $Ra$ increases after the bifurcation from branch (2) to the secondary branch (10). This symmetry breaking leads to a primary large scale circulation spanning the domain and a secondary vortex with smaller amplitude in one of the corners depending on which symmetric solution one is referring to. As $Ra\to 10^{5}$, the difference in the flow pattern of the velocity magnitude between branch (1) and (10) is the secondary vortex in the corner of the convection cell. The symmetry breaking is also obvious in the temperature field (see Fig. 6(b)) with the potential energy of the flow increasing and then decreasing slightly in contrast to the increase of the kinetic energy and the convective heat transfer (see also Fig. 3(c)). The largest growth rate shows that the states are unstable to a real eigenmode (see Fig. 6(d)). The subsequent growth rates shown in Fig. 6(c) are negative and almost constant in this range of Rayleigh numbers.

Branch (3) originates from the eigenmode of the conducting state at $Ra_{c}^{(3)}=19655$, which is depicted in Fig. 7. This branch preserves the centro-symmetry of the 3rd eigenmode and has a turning point located at $Ra\approx 42130$. Interestingly, the lower part of the branch, coloured in red in Fig. 7(a,b) bifurcates from branch (5) around $Ra\approx 27000$. The latter branch emanates from the fourth eigenmode and is illustrated later in Fig. 10. The largest growth rates of the upper and lower parts of branch (3) are displayed in Fig. 7(c,d), together with their associated frequencies. We find that the steady states in the branch are unstable with a zero frequency throughout the range of Rayleigh numbers for which the branch exists. Moreover, we observe that the eigenmode coloured in purple in Fig. 7(c) traverses zero at $Ra\approx 23800$, indicating the presence of a secondary bifurcation, which we will analyse in the following paragraph. We point out that this branch has not been originally discovered by deflated continuation (our method does not guarantee to find all the steady states of a problem) but obtained by continuing the third linear state with finer steps in the Rayleigh number. This illustrates that deflated continuation may be complemented by standard bifurcation analysis techniques, such as arclength continuation and branch-switching algorithms.

We now focus on branch (12), which bifurcates from branch (3) at $Ra\approx 23800$, and display the evolution of the steady states in the branch in Fig. 8 along with its stability analysis. This branch breaks the centro-symmetry of branch (3) and the pattern of the velocity magnitude becomes more complex as $Ra$ increases. This complex pattern is one reason for the meandering of the potential energy in this range of Rayleigh numbers and enhances the convective heat transfer (see Fig. 3(c)) by mixing the temperature field. Fig. 8(c) shows that the steady states of branch (12) are unstable to two eigenmodes until $Ra\approx 39500$, where they coalesce. At this Rayleigh number, we see that the two most unstable eigenmodes coalesce into an unstable complex conjugate pair of eigenvalues (see Fig. 8(d)). All the subsequent eigenmodes are found to be stable with either $\mathcal{I}(\lambda)=0$ or $\mathcal{I}(\lambda)\neq 0$.

Fig. 10 shows bifurcation diagrams and the linear stability results of branch (4) which originates from the fourth eigenmode at $Ra_{c}^{(4)}=23346$ through a transcritical bifurcation. The flow pattern of the velocity magnitude has a symmetric form of an array of four vortices, which is conserved over the whole range of the Rayleigh numbers we consider as Fig. 10(a) suggests. Similarly, the symmetric pattern of the temperature field remains unaffected in the regime $23346\leq Ra\leq 10^{5}$ (see Fig. 10(b)). We have verified numerically that the symmetries are conserved within this range of $Ra$. Fig. 10(c) shows the four largest growth rates as a function of $Ra$. At Rayleigh numbers close to $Ra_{c}^{(4)}=23346$ we see that three of the eigenmodes are real and unstable, while the fourth one is real and stable. As $Ra$ increases, the two most unstable eigenmodes (coloured in red and green, respectively, in Fig. 10(c)) coalesce at $Ra\approx 56000$ into a complex unstable eigenmode whose growth rate and oscillation frequency increases as $Ra\to 10^{5}$. Similarly, the other two real unstable modes coalesce at $Ra\approx 40000$ into a complex conjugate pair of eigenmodes which are unstable for $40000<Ra\leq 10^{5}$. At $Ra\approx 52000$ and $Ra\approx 73500$, we observe in Fig. 10(c) that the growth rate coloured in blue vanishes with a nonzero frequency, leading to Hopf bifurcations where periodic solutions to the time-dependent problem (1) arise.

We close our discussion on the branches arising from the conducting steady state with branch (5), presented in Fig. 10. This branch originates from the fourth eigenmode of the conducting state at $Ra\approx 23346$ through the same transcritical bifurcation as branch (4). The steady states in branch (5) conserve the Boussinesq and $x$ symmetries, with a velocity field consisting of four vortices, located symmetrically at each corner of the domain. The velocity fields of branches (4) and (5), illustrated in Fig. 11, are different, indicating that branch (5) is not a symmetric version of branch (4). The stability analysis reveals that the steady states in this branch are very sensitive to perturbations and unstable to a real eigenmode throughout the range of $Ra$ considered in this study. Moreover, the largest growth rate at $Ra=10^{5}$ is equal to $\mathcal{R}(\lambda)\approx 188$, which is approximately three times larger than the largest growth rate of the steady state in branch (4) at the same Rayleigh number (see Fig. 10(c) and Fig. 10(c)). One of the eigenmodes vanishes at $Ra\approx 27000$ giving rise to branch (3), depicted in Fig. 7. Finally, the real part of the eigenvalue depicted in blue in Fig. 10(c-d) crosses zero at $Ra\approx 57500$, which shows the existence of a Hopf bifurcation.

We now analyze the branches bifurcating from branches (4) and (5). The fourth largest growth rate of branch (4), coloured in blue in Fig. 10(c), crosses zero at $Ra\approx 38150$, leading to the bifurcating branch (6) depicted in Fig. 13. We observe that the Boussinesq symmetry of branch (4) is broken by the bifurcation, while the symmetry with respect to the $x$ axis is preserved. Additionally, the linear stability analysis indicates that this branch is unstable with a zero frequency except in the interval $Ra\in[72500,93000]$, where the largest growth rate is associated with a positive frequency, as shown by the eigenvalues coloured in blue in Fig. 13(a) and (c). Moreover, the third largest growth rate (coloured in orange) crosses zero around $Ra\approx 63500$, indicating the presence of a bifurcation which has not been discovered by deflation.

Finally, a branch, denoted by branch (7), bifurcates from branch (5) at $Ra\approx 46300$. The steady states in this branch are illustrated in Fig. 13. We observe on the velocity magnitude and temperature field that the symmetry of branch (5) with respect to the $x$ axis is maintained while the Boussinesq symmetry is broken. Branches (6) and (7) might seem nearly symmetric but they arise at different Rayleigh numbers and have different eigenvalues. The linear stability analysis of branch (7) shows that the steady states are unstable to a real eigenmode, with a largest growth rate around $\mathcal{R}(\lambda)\approx 150$ throughout the range of $Ra$ considered in this study. In Fig. 13(c) and (d), we observe that the third eigenvalue (plotted in green) crosses zero at $Ra\approx 92000$, resulting in a tertiary branch that breaks the remaining symmetry. The branch has been discovered by deflation but we chose to not report it in the study.

We now focus on the disconnected branch (11) (see Fig. 3), which comprises an upper and a lower branch, depicted in blue and red in Fig. 14, respectively. At the upper (blue) branch the flow pattern of the velocity magnitude is a three vortex state with one large titled vortex spanning the domain and two smaller vortices located at the corners. As $Ra\to 10^{5}$ we notice that the state looks similar to the branches (1) and (10). The difference is that branch (10) involves only one smaller vortex at the corner of the domain and branch (1) has no vortices at the corners, while branch (11) has two smaller vortices located at the corners. At $Ra\approx 52000$ the branch takes the form of an S-shaped curve with hysteresis, similar to the cusp bifurcation [46].

As the saddle-node bifurcation is approached from below, the growth rate of the most unstable eigenmode decreases toward zero at the critical point (see Fig. 14(c)). Above the bifurcation point $\mathcal{R}(\lambda)$ increases and seems to saturate as $Ra\to 10^{5}$. The frequency corresponding to the growth rate is zero in the whole range of $Ra$ numbers suggesting that the steady state is unstable to a real eigenmode in this regime. The potential energy of the upper (blue) branch decreases with $Ra$ in contrast to the kinetic energy. As $Ra$ increases the flow pattern of the temperature becomes more efficient at transferring heat due to the large scale circulation in the flow (see Fig. 3(c)). The flow pattern of the velocity of the upper (blue) branch, shown for $Ra=10^{5}$ (see Fig. 14(a)), exhibits a large scale vortex and two smaller vortices at the corners of the domain. This is reminiscent of the flow patterns from studies on the dynamics of random-in-time flow reversals [47, 29, 30, 31]. These reversals in the flow direction of the large scale circulation at irregular time intervals have been observed in domains with free-slip or no-slip boundaries and in a narrow band of Rayleigh numbers in the range $10^{7}\leq Ra\leq 10^{9}$ depending on the Prandtl number [29, 30, 31]. In this regime the system fluctuates in time between essentially the state we show for $Ra=10^{5}$ (see Fig. 14(a)) and its mirror symmetric state. These large changes in the dynamics of the flow are caused by bifurcations over a turbulent background [27]. Such bifurcations have been observed in many other types of flows [48, 49, 50, 51, 52] and understanding their dynamics is an open question of great interest.

At the lower (red) branch illustrated in Fig. 14(a) we observe the central vortex to be oriented more vertically, giving space to the smaller vortices at the corners to extend along the $z$ axis and almost reach the height of the domain. This flow pattern is an extension of the two-vortex pattern observed in branch (5) (see Fig. 6) and it has similar features to the three-vortex pattern found near the linear limit of branch (3) (see Fig. 7). The patterns of the kinetic and potential energy of this branch are almost invariant within the range $30000\lesssim Ra\leq 10^{5}$. The flow pattern of the temperature clearly manages to effectively push the hot fluid toward the cold plate at the top and vice versa. This efficiency in convective heat transfer is clearly depicted in Fig. 3(c), where this part of branch (11) has the highest Nusselt number for $Ra>95000$. Analysing the stability of the lower (red) branch, we find that the steady solutions become increasingly more stable as the Rayleigh number increases, with $\mathcal{R}(\lambda)\to 0$ at $Ra\approx 10^{5}$, leading to a Hopf bifurcation (see Fig. 14(e)). On the other hand, the frequency corresponding to the growth rate satisfies $\mathcal{I}(\lambda)=0$ at the threshold of the disconnected branch and as $Ra$ increases the instability becomes progressively oscillatory (see Fig. 14(f)). The subsequent growth rates presented in Fig. 14(e) remain negative as $Ra$ increases suggesting that this part of branch (11) is another preferred steady state of the system particularly for $Ra>10^{5}$.

Branch (8) illustrated in Fig. 15 bifurcates from the lower part of the disconnected branch (11) at $Ra\approx 37000$ (see also Fig. 3). In this case, the pattern of the velocity transitions into a main bent vortex with a vertically elongated vortex on its left and a smaller vortex at the bottom right corner (see Fig. 15(a)). From Fig. 3 we can see that even though the kinetic energy of branch (8) is equal to that of the upper part of branch (11) for $Ra>50000$, the potential energy is much smaller and apparently the smallest among the branches we have considered. The pattern of the temperature is essentially a hot plume that rises effectively all the way to the cold plate (see Fig. 15(b)). This flow pattern becomes increasingly more efficient at transferring heat across the domain as $Ra$ increases (see Fig. 3(c)). However, it is not as efficient in convecting heat as the pattern of branch (11) (see lower branch in 14(b)) because two hot plumes are obviously more effective than one. The steady state of branch (8) is unstable to a real eigenmode over the range $37000\lesssim Ra\leq 10^{5}$ except for a small range around $Ra=50000$ where it becomes unstable to a complex eigenmode. This can be seen from Fig. 15(c) and 15(d), where the two largest growth rates (coloured in red and blue) swap with each other in this small range of Rayleigh numbers. Finally, the real part of the eigenvalue depicted in blue in Fig. 15(c-d) crosses zero at $Ra\approx 62000$, giving birth to a Hopf bifurcation.

The second disconnected branch that we find in this range of Rayleigh numbers is branch (9) and arises at $Ra\approx 70000$ (see Fig. 16). From Fig. 16(a), we see that all the symmetries of the flow are broken in this branch and the flow pattern consists of a number of vortices with different orientations. These flow structures are not persistent over the range of parameters that we consider, and hence the heat is not transferred effectively across the domain (see also Fig. 3(c)). Here, Fig. 16 shows that the potential energy of branch (9) increases along with the kinetic energy unlike in most of the other branches we considered and from Fig. 3(b) we see that it is significant in comparison to many other branches. Moreover, from the linear stability analysis of the steady states we find that both the upper and lower parts of this disconnected branch are unstable to a real eigenmode, i.e. $\mathcal{R}(\lambda)>0$ and $\mathcal{I}(\lambda)=0$, with large growth rates. This suggests that this branch is very unlikely to be observed in a laboratory experiment, and highlights the effectiveness of deflated continuation in identifying such unstable disconnected branches and steady states.