3.1. Linear stability and the linear onset branch

Although a strong magnetic field can completely suppress convection, the motionless conducting state becomes linearly unstable with sufficient thermal driving in an infinitely wide and long domain, with the critical Rayleigh number ${{\textit{Ra}}}_c \sim {{\textit{Ha}}}^2$ (Chandrasekhar Reference Chandrasekhar1961) (shown by the open black star in figure 1 a). However, including insulating sidewalls results in a different linear instability giving rise to wall modes beneath the onset of bulk convection (Houchens et al. Reference Houchens, Witkowski and Walker2002; Busse Reference Busse2008; Bhattacharya et al. Reference Bhattacharya, Boeck, Krasnov and Schumacher2024), even in large-aspect ratio domains (also with finite sidewall conductivity). The critical Rayleigh number for a sidewall in a semi-infinite domain with free-slip top and bottom boundaries is ${{\textit{Ra}}}_c \sim {{\textit{Ha}}}^{3/2}$ for large ${\textit{Ha}}$ (Busse Reference Busse2008) (shown by the filled black star in figure 1 a). Since no stability theory has been performed for the domain and BCs we consider, or any other fully confined domain, we first perform a linear stability analysis for the conducting state. The linear stability threshold is obtained by perturbing the conducting state with random noise of magnitude $O(10^{-10})$ and scanning ${\textit{Ra}}$ to find an approximate threshold. Two equilibria are converged on each side of the threshold with $|\partial _t E_k| = O(10^{-14})$ where $E_k = (1/2)\int _\varOmega (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}) \,\mathrm{d} \boldsymbol{x}$ , ensuring agreement between the fastest growing/slowest decaying modes. This modal structure (similar to figure 1 b for the linear onset branch (LB) at ${{\textit{Ra}}}=8\times 10^5$ ) is used as the spatial eigenmode to bound the linear stability threshold ${{\textit{Ra}}}_{c,L}$ by bisection. For ${{\textit{Ha}}}=500$ we have $6.89 \times 10^5 \lt {{\textit{Ra}}}_{c,L} \lt 6.90 \times 10^5$ and at ${{\textit{Ha}}}=1000$ we have $ 1.83 \times 10^6 \lt {{\textit{Ra}}}_{c,L} \lt 1.84 \times 10^6$ . Locally, based on these two points, ${{\textit{Ra}}}_{c,L}\sim {{\textit{Ha}}}^{1.412\pm 0.005}$ , where the error corresponds to the uncertainty in ${{\textit{Ra}}}_{c,L}$ at each ${\textit{Ha}}$ . Although this is a local measurement of the scaling at finite ${\textit{Ha}}$ , it is close to the asymptotic result for a single sidewall obtained by Busse (Reference Busse2008) where ${{\textit{Ra}}}_{c}\sim {{\textit{Ha}}}^{3/2}$ at leading order. The stable equilibria born from the linear instability (LB in the bifurcation diagram in figure 1 a) exhibit a discrete fourfold rotational symmetry in the vertical velocity field (figure 1 b (LB)), and is thus invariant under the rotation $r: rw(x,y,z) = \mathcal{R}_{\pi /2} w(x,y,z),$ where $\mathcal{R}_\theta$ represents a rotation $\theta$ about $\boldsymbol{e}_z$ . The solution additionally has reflection symmetries in the $x$ - $y$ plane (about the $x$ and $y$ axes, and about the two diagonals). Thus, these symmetries may be described by a group isomorphic to the dihedral group $D_4$ , generated by $r$ , the identity $\mathrm{id}$ and the reflection $s$ about the $x=y$ diagonal, which we will denote $D_4^+ = \{\mathrm{id}, r, r^2, r^3, s, sr, sr^2, sr^3\}$ . These solutions have thin convective wall modes with a two-layer structure near the sidewalls, with convection almost completely suppressed in the bulk. Each sidewall features two counter-rotating rolls, which lead to the spatially periodic pattern of up/down flows shown in the vertical velocity isosurfaces in figure 1(b) (LB). A secondary solution can be obtained by changing the direction of the rolls and flipping the vertical coordinate $ \textit{z}\kern-5.5pt{\scriptstyle -}\,:\,\textit{z}\kern-5.5pt{\scriptstyle -} w(x,y,z)= -w(x,y,-z)$ , and thus the bifurcation is a supercritical pitchfork, which has broken the $\mathbb{Z}_2$ Boussinesq symmetry of the $D_4 \times \mathbb{Z}_2 = D_4^+ \times \{\mathrm{id},\textit{z}\kern-5.5pt{\scriptstyle -}\}$ conducting state. Increasing ${\textit{Ra}}$ results in protrusions growing from the wall modes into the bulk of the flow, similar to other observed solutions (Liu et al. Reference Liu, Krasnov and Schumacher2018; McCormack et al. Reference McCormack, Teimurazov, Shishkina and Linkmann2023; Xu et al. Reference Xu, Horn and Aurnou2023), and the equilibria eventually lose stability in the range $ 3 \times 10^6 \lt {{\textit{Ra}}} \lt 4\times 10^6$ . Since solutions on this branch are only weakly unstable up to ${{\textit{Ra}}}\approx 10^7$ , we have estimated the continuation of the unstable branch in figure 1(a).

3.2. Subcritical branch

Having determined the threshold for a linear instability, we now assess the nonlinear stability of the conducting state and search for subcritical solutions. We note that the solutions on the LB have a different spatial symmetry than other wall-mode solutions in the same domain (McCormack et al. Reference McCormack, Teimurazov, Shishkina and Linkmann2023). Thus, to try to construct a different branch of solutions, we initialise from a solution obtained far from the threshold of the linear instability ( ${{\textit{Ra}}}=5\times 10^7$ ) at the given ${\textit{Ha}}$ , which displays chaotic dynamics, and then slowly reduce ${\textit{Ra}}$ until we encounter equilibria. The obtained solutions again feature wall-mode structures but indeed have a different symmetry in the vertical velocity field (the same as in McCormack et al. (Reference McCormack, Teimurazov, Shishkina and Linkmann2023)), now invariant under the rotation $\textit{z}\kern-5.5pt{\scriptstyle -} r: (\textit{z}\kern-5.5pt{\scriptstyle -} r) w(x,y,z) = -\mathcal{R}_{\pi /2} w(x,y,-z)$ which features a reflection in sign and vertical coordinate (figure 1 b for the subcritical branch (SB)). This state still exhibits reflection symmetries about the diagonals in the $x{-}y$ plane, but also along the $x$ and $y$ axes, this time accompanied by a reflection in sign and vertical coordinate. Thus, the symmetry group is again isomorphic to $D_4$ but now with the different generators, $r_- = (\textit{z}\kern-5.5pt{\scriptstyle -} r), \mathrm{id}$ , and $s$ , which we will denote $D_4^- = \{\mathrm{id}, r_-, r_-^2, r_-^3, s, sr_-, sr_-^2, sr_-^3\}$ . Importantly, we have found that this solution branch extends beneath the linear stability threshold ${{\textit{Ra}}}_{c,L}$ (figure 1 a) to as low as ${{\textit{Ra}}}=6\times 10^5$ $(\approx 0.87 {{\textit{Ra}}}_c)$ for ${{\textit{Ha}}}=500$ , and ${{\textit{Ra}}}=1.7\times 10^6$ $(\approx 0.93 {{\textit{Ra}}}_c)$ for ${{\textit{Ha}}}=1000$ , and is thus subcritical. These simulations required lengths of approximately 800 or 700 free-fall times (6.5 or 3.4 diffusion times), respectively, to assess stability. Solutions at lower ${\textit{Ra}}$ evolve extremely slowly, and we have not been able to determine whether they are stable. The simplest explanation from the given data is that this SB stems from a saddle-node bifurcation, although we have been unable to determine where the unstable branch from this bifurcation reconnects with the other branches due to computational cost. We note that tests conducted using an under-resolved $60^2\times 80$ grid did not obtain the correct qualitative transition scenario exhibited by the properly resolved simulations we present here, and thus resolving the thin boundary layers on all walls is imperative to assessing the stability of the solutions. Importantly, the SB at both magnetic field strengths extends quite close to the conducting state, suggesting that experimentally small disturbances with the right modal structure could initiate the transition subcritically. Although the magnitude of the subcriticality ( ${{\textit{Ra}}}/{{\textit{Ra}}}_{c,L} \approx 0.9$ in both cases) is quite small, the clear difference in the symmetries of the solutions is promising and should be exploited to measure this phenomenon in experiments. Additionally, the difference between the lowest ${\textit{Ra}}$ subcritical solution found and $Ra_{c,L}$ grows with increased magnetic field strength by a factor of approximately 1.5 from ${{\textit{Ha}}}=500$ to ${{\textit{Ha}}}=1000$ , suggesting that the gap in ${\textit{Ra}}$ between the point of linear instability and the point of the expected saddle-node bifurcation grows with increased magnetic field strength as approximately ${{\textit{Ha}}}^{3/2}$ . Thus, the subcriticality becomes increasingly important at higher ${\textit{Ha}}$ .

3.3. Mixed symmetry branch

A third branch of equilibrium solutions denoted the mixed symmetry branch (MB) has also been identified, first at ${{\textit{Ha}}}=1000$ by trying a number of different initial conditions at parameter values where stable equilibria exist. This branch has then been continued to ${{\textit{Ha}}}=500$ (figure 1 a). The solutions have a symmetry that is mixed between the two other solutions, featuring one/two roll combinations on adjacent walls with matched symmetry on opposing walls (figure 1 b (MB)). The vertical velocity field is invariant under the $180^{\circ }$ rotation $rr_-: (rr_-)w(x,y,z) = -\mathcal{R}_{\pi } w(x,y,-z)$ , a reflection in $y$ , and a reflection in $x$ accompanied by a reflection in sign and vertical coordinate. Thus, the symmetry group is isomorphic to the Klein 4-group $K_4 = \{\mathrm{id}, rr_-, sr^3, sr_-\}$ . These stable solutions are interesting as they suggest that narrow/wide rolls characteristic of the LB/SB, respectively, may coexist. This is relevant to larger-aspect-ratio domains where recent simulations from Wu et al. (Reference Wu, Chen and Ni2025) show near-onset solutions that have various combinations of wide/narrow rolls. This suggests that the results are unlikely to be specific to this geometry or aspect ratio.

3.4. Amplitude equations from symmetry

A model amplitude system which captures the main features of the bifurcation diagram (figure 1 a) may be deduced from symmetries of the system. Associated with each of the three equilibrium solutions found, a secondary solution can be obtained by reversing the direction of the rolls and flipping in the vertical direction $z\mapsto -z$ . As ${\textit{Ra}}$ is increased, the LB first bifurcates from the conducting state through a supercritical pitchfork bifurcation. We then suppose that the MB and SB subsequently bifurcate through subcritical pitchfork bifurcations as ${\textit{Ra}}$ is increased further. Close to onset, when the reduced Rayleigh number $R \sim ({{\textit{Ra}}} - {{\textit{Ra}}}_{c,L})/{{\textit{Ra}}}_{c,L} \ll 1$ , we assume that only these three modes are relevant to the dynamics of the system, which can thus be described by three amplitudes $A_n(t)$ corresponding to the LB, SB and MB ( $n={1,2,3}$ ) equilibria, respectively:

(3.1) \begin{equation} \boldsymbol{\varPsi }(x,y,z,t) = \sum _{n=1}^3 \boldsymbol{\varPsi }_n(x,y,z,t) + \mathrm{s.m.} = \sum _{n=1}^3 A_n(t) \boldsymbol{\psi }_n(x,y,z) + \mathrm{s.m.}, \end{equation}

where $\boldsymbol{\varPsi }$ is the state vector of the system and $\mathrm{s.m.}$ are stable modes which will be neglected. Since for each solution $\boldsymbol{\varPsi}_{\!n}(x,y,z,t)$ there exists a distinct antisymmetric solution $\textit{z}\kern-5.5pt{\scriptstyle -} \boldsymbol{\varPsi }_n(x,y,z,t) = -\boldsymbol{\varPsi }_n(x,y,-z,t)$ , the equation that governs the evolution of the amplitudes must respect this symmetry. Thus, writing a general nonlinear evolution equation for the amplitudes as $\partial _t{\boldsymbol{A}} = \mathcal{N}(\boldsymbol{A})$ , where $\boldsymbol{A} = [A_1, A_2, A_3]^\textrm T$ , we note that $\mathcal{N}$ must be equivariant with respect to the reflection $\varpi \boldsymbol{A} \equiv -\boldsymbol{A}$ , meaning the commutative property $(\varpi \mathcal{N})(\boldsymbol{A}) = (\mathcal{N}\varpi )(\boldsymbol{A)}$ must hold (Crawford & Knobloch Reference Crawford and Knobloch1991). This implies that $\mathcal{N}$ is odd $-\mathcal{N}(\boldsymbol{A}) = \mathcal{N}(-\boldsymbol{A})$ , and thus the evolution equation must have the form $\partial _t{\boldsymbol{A}} = \boldsymbol{\mathcal{F}}(A_k^2,R)\boldsymbol{A}$ , for $k=\{1,2,3\}$ . We note that these equations are also equivariant with respect to the additional symmetries of the different states. The information about these symmetries is contained in the eigenfunctions $\boldsymbol{\psi }_n$ and do not change the amplitudes themselves. We then Taylor expand $\boldsymbol{\mathcal{F}}$ , for $R\ll 1$ , and truncate the system to the lowest order that qualitatively captures the desired dynamics, obtaining

(3.2a) \begin{align} \partial _t{A}_1 &= \big[R - b_1A_2^2 - c_1A_3^2 - d_1A_1^2\big]A_1, &(\textrm{LB}) \end{align}

(3.2b) \begin{align} \partial _t{A}_2 &= \big[(R-a_2) - b_2A_1^2 - c_2A_3^2 + d_2A_2^2 - e_2A_2^4\big]A_2, &(\textrm {SB}) \end{align}

(3.2c) \begin{align} \partial _t{A}_3 &= \big[(R-a_3) - b_3A_1^2 - c_3A_2^2 + d_3A_3^2 - e_3A_3^4\big]A_3, &(\textrm{MB}) \end{align}

where $a_n, b_n, c_n, d_n, e_n\gt 0$ are undetermined real coefficients. The positions of the subcritical pitchfork bifurcations of the SB and MB branches are set by the coefficients $a_2$ and $a_3$ . The saddle-node bifurcations occur at $R_{sn,n} = a_n - d_n^2/4e_n$ at amplitudes $A_n = \pm \sqrt {d_n/2e_n}$ for $n=\{2,3\}$ . The coupling coefficients $b_i, c_i$ determine the bifurcation points, stability and domains of existence of the mixed-mode solutions. For illustration, we have set $a_2=0.5, a_3=0.25, b_1 = 1.5, e_2=b_2=c_n=2, d_2=3, d_3 = 3.5, b_3=6,e_3=8$ . The bifurcation diagram for this choice of parameters is shown in figure 2, exhibiting good qualitative agreement with the bifurcation diagram of the full system in figure 1(a). Thus, although additional complexity may exist in the full system, the behaviour we observe in the DNS can be rationalised by a simplified system, derived from assumptions about the origin of the SB and MB branches, and a consideration of the symmetries of the problem.

Figure 2. Bifurcation diagram of the amplitude (3.2), showing the norm of the amplitudes $\|\boldsymbol{A}\|_2$ as a function of the reduced Rayleigh number $R$ for the various stable/unstable equilibria denoted by solid/dotted lines. Single mode solutions corresponding to the LB, MB and SB states are shown in green, blue and orange, respectively. Mixed mode solutions are shown in grey. Markers show the bifurcation points.

A natural question that arises is whether this behaviour is also expected in a cylindrical geometry, often used in experiments, where the system has a continuous rotational symmetry i.e. $O(2)\times \mathbb{Z}_2$ instead of the discrete $D_4\times \mathbb{Z}_2$ as before. We describe the competition between $N$ unstable modes of azimuthal wavenumber $n$ in cylindrical coordinates by expanding the state vector as

(3.3) \begin{equation} \boldsymbol{\varPsi }(r,\theta ,z,t) = \textrm{Re}\left [\sum _{n=1}^N a_n(t) \textrm{e}^{\textit{in}\theta }\boldsymbol{\psi }_n(r,z)\right ] + \mathrm{s.m.}, \end{equation}

where $\textrm{Re}[\boldsymbol{\cdot }]$ takes the real part of the now complex eigenfunction decomposition. Again, the equation evolving the amplitudes $a_n$ must be equivariant with respect to the symmetries of the system, and thus must commute with rotations $\theta \rightarrow \theta + \theta _0$ , for $\theta_0$ in $\mathbb{R}$ , and the reflection. Thus, $\partial _t{a}_n = \mathcal{F}(|a_k|^2,R)a_n$ , for $k\in \{1,\ldots ,N\}$ . Since the conducting state now bifurcates through a circle pitchfork bifurcation, meaning the solutions are neutrally stable with respect to azimuthal rotations, we write the complex amplitudes in polar form $a_n = A_n \textrm{e}^{\mathrm{i}\varphi _n}$ , noting that due to this neutrality, $\partial _t{\varphi }_n=0$ . Again, Taylor expanding $\mathcal{N}$ , for $R\ll 1$ , we arrive at a system for the real amplitudes $A_n$ which is equivalent to (3.2) for $N$ equations. This means that the continuous symmetry of a cylinder does not preclude the possibility of having the same behaviour we have observed in the cube. Whether the same behaviour occurs as in the cylinder depends on the sign of the coefficients as described before, which must be determined to make qualitative statements about the behaviour of the system. However, we observe that the continuous symmetry of the cylinder does not affect the final structure of the equations we obtain compared with the cube. Furthermore, any finite subgroup is isomorphic to a dihedral $D_n$ or cyclic $\mathbb{Z}_n$ group and thus, the solutions bifurcating from the conducting state through these circle pitchforks will generically have $D_n$ or $\mathbb{Z}_n$ symmetries. This suggests that similar symmetry breaking phenomena observed in the cube, may plausibly occur in the cylinder.

Figure 3. Phase portrait at ${{\textit{Ha}}}=500$ of the (a) invariant 2-torus at ${{\textit{Ra}}}=4\times 10^6$ and (c) chaotic solution at ${{\textit{Ra}}}=10^8$ constructed using time-delay embedding with a time-delay $\tau \approx 4$ . The colour map corresponds to $\|w(t+3\tau )\|_2$ . The unstable equilibrium point (LB) is shown in black. (b,d) Corresponding power spectral density of $\|w(t)\|_2$ , respectively, as a function of the frequency $f$ . Fundamental frequencies are labelled $f_i$ . Zoomed out spectra are shown in the insets.

4.1. Consequences of the subcritical transition

At all values of ${\textit{Ha}}$ studied, the route to chaos originates through a sequence of bifurcations on the SB. Often in supercritical transitions, the equilibrium solution born from the first linear instability in the system imprints the turbulent flow, organising large-scale features (e.g. Taylor vortices in turbulent TCF (Taylor Reference Taylor1923; Grossmann, Lohse & Sun Reference Grossmann, Lohse and Sun2016) or convection rolls in turbulent RBC (Chandrasekhar Reference Chandrasekhar1961; Busse Reference Busse1978; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009)). Although the subcriticality is reasonably small in this system at the considered ${\textit{Ha}}$ , it has a significant effect on the flow in the turbulent regime. At ${{\textit{Ra}}}=10^9$ ( ${\approx} 1000\,{{\textit{Ra}}}_{c,L}$ ), the flow field exhibits clear multiscale behaviour with small-scale vortices being observed in an instantaneous snapshot of the flow in figure 4(a). However, unlike the supercritical cases, this turbulent flow has been formed through bifurcations on the SB and thus the turbulent flow retains properties of the subcritical equilibria and not the linear instability. This is clearly observed in the mean flow of the solution which retains the same $D_4^-$ symmetry as the subcritical equilibria, and not the $D_4^+$ symmetry of the solutions born from the linear instability. This means that understanding the subcritical nature of the transition is imperative to understanding properties of turbulent MC, such as the globally organising flow structures.