1. Introduction
Vertical convection, in which a layer of fluid is confined between two vertical plates maintained at different temperatures, is relevant for industrial applications, including the control of insulation properties of double-glazed windows. Vertical convection also serves as a model system in the geophysical context to describe convectively driven flows in the Earth, the ocean and the atmosphere. Moreover, vertical convection is a fundamental hydrodynamics problem in its own right, as a prototype for studying pattern formation mechanisms within spatially extended driven dissipative nonlinear out-of-equilibrium systems. In our companion paper Zheng, Tuckerman & Schneider (Reference Zheng, Tuckerman and Schneider2024), we studied a domain in which the transverse (or spanwise) direction was taken to be of the same length as the distance between the plates (the wall-normal direction), with the vertical dimension (parallel to gravity) chosen large compared with both. Consequently, flow patterns are primarily two-dimensional (2-D), with variations predominantly in the vertical and wall-normal direction. Here, we will consider an extended three-dimensional (3-D) geometry, in which the transverse and vertical dimensions are both large compared with the inter-plate spacing and, thus, flow patterns vary in two extended directions.
We begin by briefly surveying 3-D numerical investigations of vertical convection. Chait & Korpela (Reference Chait and Korpela1989), Henry & Buffat (Reference Henry and Buffat1998) and Xin & Le Quéré (Reference Xin and Le Quéré2002) analysed the instability of 2-D nonlinear flow (transverse rolls) to 3-D perturbations in order to determine when and whether the flow could be assumed to be 2-D. In Rayleigh–Bénard convection, the stability thresholds in Rayleigh number (
$Ra$), Prandtl number (
$Pr$) and 2-D roll wavelength delimit a volume that is called the Busse balloon (Busse Reference Busse1978), named after the researcher who has been at the forefront of pattern formation research in Rayleigh–Bénard convection. Busse later also transferred his analysis to vertical convection. Using the approximation (corresponding to infinite thermal diffusivity) that the temperature retains its linear conductive profile, Nagata & Busse (Reference Nagata and Busse1983) computed a fully nonlinear 3-D solution which is probably the diamond roll state (FP2) to be described in § 3. Such 3-D solutions have sometimes been termed tertiary solutions, whereas the laminar and 2-D transverse roll solutions are called primary and secondary, respectively. Clever & Busse (Reference Clever and Busse1995) extended the computation of 3-D solutions to 
$Pr=0.71$, corresponding to convection in air, the case we study in this paper.
Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015, Reference Gao, Podvin, Sergent, Xin and Chergui2018) combined linear and weakly nonlinear theory as well as direct numerical simulations (DNS) to study the 3-D flow. Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015) studied the equilibria and periodic orbits in a computational domain of size 
$[L_x, L_y, L_z] = [1, 1, 10]$, the same domain we consider in our companion paper Zheng et al. (Reference Zheng, Tuckerman and Schneider2024). In order to study secondary instabilities in the transverse direction of the 2-D steady rolls, Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) computed their linear stability. Their analysis showed two types of instabilities, with spanwise wavelengths of about four and eight. They consequently extended the spanwise length of the domain from unity to 
$L_y=8$ to capture both instabilities. In addition, when 
$L_z=9$ the Rayleigh number thresholds of both types of 3-D instabilities are close, motivating them to decrease 
$L_z$ from 
$10$ to 
$9$ in order to study the competition between both instabilities destabilising 2-D rolls. Like a spanwise domain size of 
$L_z=10$, a domain with 
$L_z = 9$ also accommodates four co-rotating rolls in the primary instability of the base state and is large enough to allow interactions between rolls. Cimarelli & Angeli (Reference Cimarelli and Angeli2017) and Cingi, Cimarelli & Angeli (Reference Cingi, Cimarelli and Angeli2021) unsuccessfully attempted to explain the results of Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent, Xin and Chergui2018) from a bifurcation-theoretic point of view.
In this paper, we study vertical convection in air (
$Pr=0.71$) in the configuration 
$[L_x,L_y,L_z]=[1,8,9]$. Similarly to the approach described in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024), we extend previous studies by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) that were based primarily on time-stepping by using numerical continuation and stability analysis. This unravels the bifurcation-theoretic origins of complex flows and the connections between them. This approach of explaining patterns and their dynamics in terms of equilibria and periodic orbits has been applied successfully to inclined layer convection where fascinating convection patterns were previously observed in DNS and experimentally by Daniels, Plapp & Bodenschatz (Reference Daniels, Plapp and Bodenschatz2000). Through a numerical bifurcation analysis, Reetz & Schneider (Reference Reetz and Schneider2020) and Reetz, Subramanian & Schneider (Reference Reetz, Subramanian and Schneider2020) identified the invariant solutions underlying most of the patterns and constructed bifurcation diagrams connecting them. These invariant solutions capture key features and dynamics of the observed patterns and the bifurcation diagrams reveal their origin. Here, we follow the same strategy to explain flow patterns in vertical convection in a somewhat larger domain.
Using parametric continuation techniques that can follow states irrespective of their stability, we describe the discovery of three new branches of steady states, which, together with those observed by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) via time integration, brings the number of branches observed thus far to six. Several of these new states bifurcate simultaneously, at the same value of the control parameter, despite not being related by symmetry. We have shown that this otherwise non-generic phenomenon is explained by the fact that the parent branches have 
$D_4$ symmetry; see Swift (Reference Swift1985), Knobloch (Reference Knobloch1986), Chossat & Iooss (Reference Chossat and Iooss1994), Bergeon, Henry & Knobloch (Reference Bergeon, Henry and Knobloch2001) and Reetz et al. (Reference Reetz, Subramanian and Schneider2020). In our geometry, 
$D_4$ symmetry leads to simultaneous bifurcations to states that are aligned with respect to the transverse and vertical directions, and others which are diagonal with respect to them. Competition between aligned and diagonal states is also seen in two periodic orbits (observed by Gao et al. Reference Gao, Podvin, Sergent, Xin and Chergui2018), that consist of diagonal excursions from states which are more aligned. We have also discovered two new periodic orbits.
Most of the steady states and periodic orbits that we have identified are unstable. While these are not directly observed in time-dependent simulations, following unstable branches is essential for understanding the origin of stable states and for constructing a bifurcation diagram unifying the solutions to a problem. Moreover, unstable states play the role of way-stations, near which chaotic or turbulent trajectories spend much of their time. These are believed to form the core structures supporting weakly turbulent dynamics. Among the unstable periodic orbits that may influence trajectories of a fluid-dynamical system, we have discovered some whose branches terminate in global bifurcations, leading to their disappearance. Although there have been a number of computations of global bifurcations in hydrodynamic systems (Tuckerman & Barkley Reference Tuckerman and Barkley1988; Prat, Mercader & Knobloch Reference Prat, Mercader and Knobloch2002; Millour, Labrosse & Tric Reference Millour, Labrosse and Tric2003; Nore et al. Reference Nore, Tuckerman, Daube and Xin2003; Abshagen et al. Reference Abshagen, Lopez, Marques and Pfister2005; Bordja et al. Reference Bordja, Tuckerman, Witkowski, Navarro, Dwight and Bessaih2010; Bengana & Tuckerman Reference Bengana and Tuckerman2019; Reetz et al. Reference Reetz, Subramanian and Schneider2020), we are not aware of previous calculations of heteroclinic or homoclinic cycles in vertical convection.
The remainder of this article is organised as follows. In § 2 we summarise the key numerical methods used in our research which are already presented in detail in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024). The results from the bifurcation analysis are presented in § 3 for fixed points and in § 4 for periodic orbits. Concluding remarks and future research directions are outlined in § 5.
2. System and methods
We refer the reader to Reetz (Reference Reetz2019), Reetz & Schneider (Reference Reetz and Schneider2020), Reetz et al. (Reference Reetz, Subramanian and Schneider2020), and Zheng et al. (Reference Zheng, Tuckerman and Schneider2024) for detailed descriptions of the numerical methods used in the research. Here, we only summarise the key points.
2.1. The DNS of vertical convection
The vertical convection system is studied numerically by performing DNS with the ILC extension module of the Channelflow 2.0 code (Gibson et al. Reference Gibson, Reetz, Azmi, Ferraro, Kreilos, Schrobsdorff, Farano, Yesil, Schütz and Culpo2021), to solve the non-dimensionalised Oberbeck–Boussinesq equations
$$\begin{gather} \dfrac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} ={-}\boldsymbol{\nabla} p + \left(\frac{Pr}{Ra}\right)^{1/2} \nabla^2 \boldsymbol{u} + \mathcal{T}\boldsymbol e_z, \end{gather}$$
$$\begin{gather}\dfrac{\partial \mathcal{T}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla})\mathcal{T} = \left(\frac{1}{Pr\,Ra}\right)^{1/2}\nabla^2 \mathcal{T}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}= 0, \end{gather}$$
in a vertical channel, with periodic boundary conditions in 
$y$ and 
$z$, shown in figure 1. The boundary conditions in 
$x$ at the two walls are of Dirichlet type:
\begin{equation} \boldsymbol{u}(x={\pm} 0.5) = 0,\quad \mathcal{T}(x=\pm0.5)=\pm0.5. \end{equation}
Supplementary integral constraints are necessary in the periodic directions; we set the mean pressure gradient to zero in both 
$y$ and 
$z$. The laminar solution, illustrated in figure 1, is
$$\begin{gather} \boldsymbol{u}_0(x) = \frac{1}{6}\sqrt{\frac{Ra}{Pr}} \left(\frac{1}{4}x - x^3\right) \boldsymbol e_z, \end{gather}$$
$$\begin{gather}\mathcal{T}_0(x) = x, \end{gather}$$
$$\begin{gather}p_0(x) = \varPi, \end{gather}$$
with arbitrary pressure constant 
$\varPi$.
Figure 1. Schematic of the vertical convection configuration approximating 
$[L_x, L_y, L_z] = [1, 8, 9]$. The flow is bounded between two walls in 
$x$ direction at 
$x=0.5$ where the flow is heated and at 
$x=-0.5$ where the flow is cooled. The domain is periodic in 
$y$ and 
$z$ directions. Most of the visualisations that we present are taken on the 
$y$–
$z$ midplane at 
$x=0$ outlined by the dotted lines, and they are visualised in the direction of negative to positive 
$x$, as indicated by the eye and arrow. The laminar velocity and temperature are shown as the orange and green curves, respectively.
The governing equations and boundary conditions are discussed in our companion paper Zheng et al. (Reference Zheng, Tuckerman and Schneider2024). The only aspect which differs here is the domain size: instead of the narrow domain 
$[L_x, L_y, L_z] = [1, 1, 10]$ with one extended direction studied in Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013), here we study the 3-D computational domain 
$[L_x, L_y, L_z] = [1, 8, 9]$ of Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018). This domain has two extended directions and is illustrated in figure 1. This domain is spatially discretised by 
$[N_x, N_y, N_z] = [31, 96, 96]$ Chebychev–Fourier–Fourier modes.
2.2. Symmetries and computation of invariant solutions
We often refer to the symmetries of our system, the group 
$S_{VC}$, which is generated by reflection in 
$y$, combined reflection of 
$x$, 
$z$ and temperature 
$\mathcal {T}$, and translation in 
$y$ and 
$z$:
$$\begin{gather} {\rm \pi}_y[u,v,w,\mathcal{T}](x,y,z) \equiv [u,-v,w,\mathcal{T}](x, -y,z) , \end{gather}$$
$$\begin{gather}{\rm \pi}_{xz}[u,v,w,\mathcal{T}](x,y,z) \equiv [{-}u,v,-w,-\mathcal{T}]({-}x, y,-z) , \end{gather}$$
$$\begin{gather}\tau(\Delta y, \Delta z)[u,v,w,\mathcal{T}](x,y,z) \equiv [u,v,w,\mathcal{T}](x, y + \Delta y,z + \Delta z), \end{gather}$$
stated more compactly as 
$S_{VC} \equiv \langle {{\rm \pi} _y, {\rm \pi}_{xz}, \tau (\Delta y, \Delta z)} \rangle \sim [O(2)]_y \times [O(2)]_{x,z}$. The groups we use are 
$Z_n$, the cyclic group of 
$n$ elements, 
$D_n$, the cyclic group of 
$n$ elements together with a non-commuting reflection, and 
$O(2)$, the group of all rotations together with a non-commuting reflection. Here 
$[O(2)]_y$ refers to reflections and translations in 
$y$, as in (2.4a) and (2.4c), respectively, whereas 
$[O(2)]_{xz}$ refers to reflections in 
$(\mathcal {T},x,z)$ as in (2.4b) and translations in 
$z$ as in (2.4c), a convention that we use in the rest of the paper where possible. Note that the generators of a group are non-unique, as is the decomposition into direct products (indicated by 
$\times$).
We adopt the shooting-based matrix-free Newton method implemented in Channelflow 2.0 to compute invariant solutions. The only difference with respect to our description in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024) arises from the presence here of homoclinic and heteroclinic orbits. Although the Newton method can converge with one shot in most of the cases (provided that the initial guess is sufficiently close to the solution), the multi-shooting method (Sánchez & Net Reference Sánchez and Net2010; van Veen, Kawahara & Atsushi Reference van Veen, Kawahara and Atsushi2011) is required in order to converge orbits with long periods (typically 
$T > 300$ in our case) that are close to a global bifurcation point and very unstable orbits. For these periodic orbits, we employ the multi-shooting method with at most six shots.
To characterise the stability of a solution, its leading eigenvalues and eigenvectors for fixed points, or Floquet exponents and Floquet modes for periodic orbits, are determined by Arnoldi iterations. When solutions have symmetries, the resulting linear stability problem has the same symmetries, leading to multiple eigenvectors sharing the same eigenvalues. In such cases, we choose the eigenvectors appropriate to our analysis either by subtracting two nonlinear flow fields along a trajectory or branch, or by imposing symmetries.
2.3. Order parameter and flow visualisation
Once an equilibrium or time-periodic solution is converged, parametric continuation in Rayleigh number is performed to construct bifurcation diagrams. Solutions are represented via the 
$L_2$-norm of their temperature deviation 
$\theta \equiv \mathcal {T} - \mathcal {T}_0$. Branches of fixed points are represented by curves showing 
$\lvert \lvert \theta \lvert \lvert _2$ as a function of 
$Ra$; for periodic orbits, the maximum and minimum of 
$\lvert \lvert \theta \lvert \lvert _2$ along an orbit are plotted. The flows are visualized via their temperature deviation fields 
$\theta$ on the y–z plane at 
$x=0$ and on the x–z plane at 
$y=4$. The thermal energy input 
$I$ due to buoyancy and the dissipation 
$D$ due to viscosity are used to plot phase portraits.
3. Fixed points
We begin by noting that the numbering used for fixed points and for periodic orbits applies only to this paper; except for FP1, the fixed points and periodic orbits here are not the same as those in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024).
3.1. Three known fixed points: FP1–FP3
Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) observed three fixed points in the domain 
$[L_x, L_y, L_z] = [1, 8, 9]$ and presented visualisations and Fourier decompositions of them. These states have been recomputed here and their flow structures are shown in figures 2(b)–2(d). In this work, we identify the bifurcations that create and destroy these states and construct a bifurcation diagram that includes stable and unstable branches. As presented in the bifurcation diagram in figure 2(a), the laminar base flow is stable until 
$Ra=5707$, where the first fixed point, FP1, bifurcates. As in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024), FP1 is called 2-D or transverse rolls. This state contains four spanwise (
$y$)-independent co-rotating convection rolls, and is shown in figures 2(b) and 3(a). Cingi et al. (Reference Cingi, Cimarelli and Angeli2021) have reported bistability between the base flow and 2-D rolls in several Rayleigh-number ranges, but their interpretation contradicts the results obtained here and also those reported by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018). In particular, Cingi et al. (Reference Cingi, Cimarelli and Angeli2021) find the laminar flow to be bistable with 2-D rolls (FP1) over the Rayleigh number range of 
$[5708, 7000]$. We believe this reported bistability to be spurious, and to almost certainly result from the use by Cingi et al. (Reference Cingi, Cimarelli and Angeli2021) of a time-stepping code to simulate a weakly unstable state without monitoring the growth or decay of perturbations nor a complementary linear stability analysis.
Figure 2. Bifurcation diagram (a) and flow structures visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$ (b–g) of six equilibria in domain 
$[L_x, L_y, L_z] = [1, 8, 9]$. The diagram shows two supercritical pitchfork bifurcations, one from the base state to FP1 (b) and another one from FP1 to FP2 (c). Fixed point FP3 (d) bifurcates from FP2 in a subcritical pitchfork bifurcation. The unstable FP4 (e) bifurcates supercritically from FP1. The unstable FP5 branch (f) bifurcates at one end subcritically from FP2, and at the another end supercritically from FP1. FP3 and FP5 bifurcate together from FP2, whereas FP4 and FP5 bifurcate together from FP1. Two small grey rectangles surround these two simultaneous bifurcations, which are also shown in the enlarged diagrams on the right. On the lower enlarged diagram, the dashed red and brown lines are distinct, but too close to one another to be distinguished. FP6 bifurcates from FP5 in two supercritical pitchfork bifurcations and it connects FP5 at two Rayleigh numbers. In (a), solid and dashed curves signify stable and unstable states, respectively. The ranges over which FP1, FP2, FP3 and FP6 are stable are 
$[5707, 6056]$, 
$[6056, 6058.5]$, 
$[6008.5, 6140]$ and 
$[6251.4, 6257.6]$, respectively. The stars in (a) indicate the locations of the visualisations of (b–g). Fixed points FP1–FP3 are discussed in Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) whereas FP4–FP6 are newly identified in this work. Other branches of equilibria exist, which we have not followed nor shown on this diagram. Flow visualisations on the 
$x$–
$z$ plane are shown in figure 3.
Fixed point FP1 loses stability at 
$Ra = 6056$ via a circle pitchfork bifurcation that breaks the 
$y$ translation symmetry 
$\tau (\Delta y,0)$ and creates FP2, shown in figures 2(c) and 3(b). We refer to these as diamond rolls, whereas Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) called them wavy rolls. Fixed point FP2 results from the subharmonic varicose instability of FP1, which was discussed in Subramanian et al. (Reference Subramanian, Brausch, Daniels, Bodenschatz, Schneider and Pesch2016) and Reetz et al. (Reference Reetz, Subramanian and Schneider2020). FP2 undergoes subcritical pitchfork bifurcations at 
$Ra = 6058.5$, so that its stability range is only 
$[6056,6058.5]$. The time-dependent simulations of Cingi et al. (Reference Cingi, Cimarelli and Angeli2021) did not detect FP2. In contrast, Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) observed FP2 as a transient at 
$Ra=6100$ and computed its threshold via a linear stability analysis. Clever & Busse (Reference Clever and Busse1995) computed a state resembling FP2 by means of a steady-state calculation. (Their threshold of about 
$Ra\approx 6295$ can perhaps be attributed to a lack of spatial resolution available in 1995.)
The bifurcation from FP2 creates FP3, which Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) call thinning rolls. Initially unstable, FP3 is stabilised by a saddle-node bifurcation at 
$Ra = 6008.5$. At higher Rayleigh number, FP3 undergoes two additional saddle-node bifurcations at 
$Ra=6265.8$ and 
$Ra=6209.56$. As pointed out by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018), FP3 can have either of two possible diagonal orientations. Figure 2(d) shows one of the two cases: the slightly wider red portions are located along a diagonal joining the top left with the bottom right.
The symmetry (isotropy) groups of FP1–FP3 are
\begin{equation}
\left.\begin{gathered}
\text{FP1:} \ \langle{{\rm \pi}_y, \tau(\Delta y,0),{\rm \pi}_{xz},\tau(0,L_z/4)} \rangle \sim [O(2)]_y \times [D_4]_{xz};\\
\text{FP2:} \ \langle{{\rm \pi}_y, \tau(L_y/2,0), {\rm \pi}_{xz}, \tau(0, L_z/2), \tau(L_y/4, -L_z/4)} \rangle \sim D_2 \times D_4; \\
\text{FP3:} \ \langle{{\rm \pi}_y{\rm \pi}_{xz}, \tau(L_y/4,-L_z/4)} \rangle \sim D_4. \end{gathered}\right\}
\end{equation}
Note that 
$\tau (L_y/4,-L_z/4)=\tau (L_y/4,3L_z/4)$, and that the symmetry groups for FP2 and FP3 cannot be divided into those related to 
$y$ and those related to 
$x,z$. The bifurcation from 
$\textrm {FP}2\rightarrow \textrm {FP}3$ breaks the 
$D_4$ symmetry of FP2.
3.2. Three new fixed points: FP4–FP6
We have also found three new branches of fixed points, FP4–FP6. Figure 2(a) shows that there is a supercritical pitchfork bifurcation at 
$Ra=6131$, at which FP4 and FP5 bifurcate simultaneously from FP1. Both FP4 and FP5 are unstable along their entire branches. (The enlarged diagram on the bottom right of figure 2(a) contains two distinct dashed red and brown lines which are too close to be distinguished.) Since FP1 is 
$y$-independent and FP4 and FP5 are not, these are circle pitchfork bifurcations, yielding FP4 and FP5 states of any phase in 
$y$. Fixed point FP4, shown in figures 2(e) and 3(d), shares with FP3 a diagonal orientation. Fixed point FP4 also consists of rolls with a slight wavy modulation along the 
$y$ direction, but this modulation is weaker than that of FP3. FP4 plays an essential role in one of the global bifurcations that we discuss in § 4.1.2.
The FP5 branch (which we refer to occasionally as the moustache branch) is shown in figures 2(f) and 3(e). After bifurcating from FP1, the FP5 branch undergoes saddle-node bifurcations at 
$Ra=6317.5$ and 
$Ra= 6034$, towards decreasing and increasing Rayleigh number, respectively, and finally terminates at 
$Ra=6058.5$ by meeting FP2 in a subcritical pitchfork bifurcation. This is not a circle pitchfork bifurcation, since the diamond branch FP2 is also 
$y$-dependent; four possible FP5 branches emanate from FP2, related to one another by translations in 
$y$ and in 
$z$. (Fixed point FP3 is also created at 
$Ra=6058.5$, in another simultaneous bifurcation that is discussed in § 3.3.) Thus, two routes connect FP1 to FP5: one route is a single circle pitchfork bifurcation and a second route is a circle pitchfork bifurcation from FP1 to FP2 followed by an ordinary pitchfork bifurcation from FP2 to FP5. The bifurcation from FP1 to FP2 breaks 
$y$ invariance whereas that from FP2 to FP5 breaks the fourfold translation symmetry 
$\tau (L_y/4,-L_z/4)$.
The last new equilibrium, FP6, shown in figures 2(g) and 3(f), is created from FP5 at 
$Ra=6164.3$ in a supercritical pitchfork bifurcation, inheriting the instability of FP5 at the bifurcation point. Fixed point FP6 becomes stable, but only over a very short range 
$Ra \in [6251.4,6257.6]$, indicated by the slight thickening of the branch in figure 2(a). (We do not discuss or show in figure 2(a) the new branches that necessarily emanate from the stabilising bifurcation at 
$Ra=6251.4$, nor the numerous other branches created at points at which the real part of an eigenvalue crosses zero. The bifurcation at 
$Ra=6257.6$ is discussed in § 4.4.) FP6 then undergoes a saddle-node bifurcation at 
$Ra=6329$ before terminating at the FP5 branch at 
$Ra=6305.8$ in another supercritical pitchfork bifurcation.
The symmetry groups of these states are
\begin{equation}
\left.\begin{gathered}
\text{FP4:} \ \langle{{\rm \pi}_y{\rm \pi}_{xz}, \tau(L_y/4,-L_z/4)} \rangle \sim D_4; \\
\text{FP5:} \ \langle{{\rm \pi}_y,{\rm \pi}_{xz}, \tau(L_y/2, L_z/2)} \rangle \sim [Z_2]_y \times [Z_2]_{xz} \times Z_2; \\
\text{FP6:} \ \langle{{\rm \pi}_y,{\rm \pi}_{xz}\tau(L_y/2,0)} \rangle \sim [Z_2]_y\times Z_2. \end{gathered}\right\}
\end{equation}
Fixed point FP1 is homogeneous in 
$y$ and the states which branch from it, directly or indirectly, are FP2 with a 
$y$ periodicity of 
$L_y/2=4$, and FP3, FP4, FP5 and FP6 with 
$y$-periodicity 
$L_y=8$. This sets the stage for 1:2 mode interaction, as analysed in detail by Armbruster, Guckenheimer & Holmes (Reference Armbruster, Guckenheimer and Holmes1988), one of whose consequences is a robust heteroclinic cycle to be discussed in § 4.2.3.
3.3. Two simultaneous bifurcations
The two enlarged bifurcation diagrams on the right of figure 2(a) depict bifurcations at which two qualitatively different branches with different symmetries are created simultaneously. Fixed points FP3 and FP5 bifurcate simultaneously from FP2 at 
$Ra=6058.5$, and FP4 and FP5 bifurcate simultaneously from FP1 at 
$Ra= 6131$. These simultaneous bifurcations can be explained by the same 
$D_4$ scenario that is discussed in detail in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024). We repeat here the normal form corresponding to bifurcation in the presence of 
$D_4$ symmetry:
$$\begin{gather} \dot{p} = (\mu - a p^2 - b q^2) p, \end{gather}$$
$$\begin{gather}\dot{q} = (\mu - b p^2 - a q^2) q. \end{gather}$$The dynamical system (3.3) has the non-trivial solutions
$$\begin{gather} p ={\pm}\sqrt{\mu/a} \quad q = 0, \end{gather}$$
$$\begin{gather}p = 0 \quad q={\pm}\sqrt{\mu/a}, \end{gather}$$
$$\begin{gather}p ={\pm}\sqrt{\mu/(a+b)} \quad q ={\pm}\sqrt{\mu/(a+b)}, \end{gather}$$
$$\begin{gather}p ={\pm}\sqrt{\mu/(a+b)} \quad q ={\mp}\sqrt{\mu/(a+b)}, \end{gather}$$i.e. two classes of solutions, (3.4a)–(3.4b), which we call here the diagonal solutions, and (3.4c)–(3.4d), which we call here the rectangular solutions, for reasons which figure 4 makes clear. The diagonal solutions are related to one another by symmetry, as are the rectangular ones, but the diagonal solutions are not related by symmetry to the rectangular solutions.
Figure 4. (a) Eigenvector 
$e_1$ responsible for 
$\textrm {FP}2\rightarrow \textrm {FP}3$ bifurcation (obtained by subtracting FP2 at 
$Ra=6058.5$ from FP3 at 
$Ra=6056$) and (b) its 
$y$-reflected version 
${\rm \pi} _y e_1$. (c,d) Superpositions 
$(e_1 \pm {\rm \pi}_y e_1)/\sqrt {2}$. (e) Eigenvector 
$e_2$ responsible for 
$\textrm {FP}2\rightarrow \textrm {FP}5$ bifurcation (obtained by subtracting FP2 at 
$Ra=6058.5$ from FP5 at 
$Ra=6056$) and (f) its quarter-diagonal translation 
$\tau (L_y/4,-L_z/4) e_2$. (g,h) Superpositions 
$(e_2 \pm \tau (L_y/4,-L_z/4) e_2)/\sqrt {2}$. All eigenvectors are visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$. The same colour bar is used in all plots.
We begin by explaining the simultaneous bifurcation from FP2. The symmetry group 
$D_4$ of FP2 is generated by the translation operator 
$\tau (L_y/4, -L_z/4)$ together with either of the reflection operators, 
${\rm \pi} _y$ or 
${\rm \pi} _{xz}$. Fixed point FP2 is invariant under any product of these operations. In the model (3.3), FP2 corresponds to the trivial solution 
$p=q=0$ from which the other solutions bifurcate.
When FP2 loses stability at 
$Ra=6058.5$, a real eigenvalue 
$\lambda _{1,2}$ crosses the imaginary axis. This double eigenvalue has a 2-D eigenspace, spanned by any two of its linearly independent eigenvectors. Figure 4(a) shows the eigenvector 
$e_1$ of FP2 giving rise to state FP3 shown in figure 2(d), whereas figure 4(b) shows its 
$y$-reflection, 
${\rm \pi} _y e_1$. Since 
${\rm \pi} _y$ belongs to the symmetry group of FP2, 
${\rm \pi} _y e_1$ is also an eigenvector of FP2, as is any superposition of 
$e_1$ and 
${\rm \pi} _y e_1$. The diagonal solution (3.4a) represents FP3 which arises from eigenvector 
$e_1$. Solution (3.4b) represents FP3
$^\prime \equiv {\rm \pi}_y$FP3, whose diagonal is reversed and which arises from eigenvector 
${\rm \pi} _y e_1$. The amplitudes of 
$e_1$ and 
${\rm \pi} _y e_1$ are represented in the model (3.3) by variables 
$p$ and 
$q$:
\begin{equation} {\rm FP3} = {\rm FP2} + p(t) e_1 + q(t) {\rm \pi}_y e_1. \end{equation} The eigenvector 
$e_2$ of FP2 leading to state FP5 is shown in figure 4(e). Eigenvector 
$e_2$ turns out to be identical to the equal superposition of 
$e_1$ and 
${\rm \pi} _y e_1$, as shown in figure 4(c). This is a manifestation of the fact that, in the model (3.3), the rectangular solutions (3.4c) and (3.4d) contain equal amplitudes of 
$p$ and 
$q$. The shifted eigenvector 
$\tau (L_y/4,-L_z/4)e_2$ leads to 
$\textrm {FP}5^{\prime } \equiv \tau (L_y/4,-L_z/4)\textrm {FP}5$; its superposition with 
$e_2$ produces 
$e_1$. Indeed, in the model (3.3), equal superpositions of rectangular solutions of types (3.4c) and (3.4d) produce the diagonal solutions of types (3.4a) and (3.4b). (Figures 4(d) and 4(h) are also eigenvectors of FP2, identical to figures 4(f) and 4(b), respectively.) In figure 4, the eigenvectors have been approximated by subtracting FP2 from FP3 and from FP5 just beyond the bifurcation point (
$Ra=6058.5$ for FP2 and 
$Ra=6056$ for FP3 and FP5). This selects the appropriate choices out of the multitude of eigenvectors of the highly symmetric FP2.
Just as solutions (3.4c) and (3.4d) are not related to solutions (3.4a) or (3.4b) by any symmetry operation, FP5 cannot be produced by a symmetry transformation from FP3. In addition, figure 2(a) makes it clear that branches FP3 and FP5 behave differently, with a different global temperature norm and different saddle-node bifurcations.
We turn now to the simultaneous bifurcations of FP4 and FP5 from FP1 at 
$Ra = 6131$. The symmetries of FP1 are generated by reflection and translation in 
$y$ together with reflection in 
$(x,z)$ and fourfold translation in 
$z$, i.e. 
$[O(2)]_y \times [D_4]_{xz}$. We again compute the eigenvectors of FP1 responsible for these two bifurcations. Taking symmetry transformations and superpositions, we obtain the eigenvector responsible for FP5 (FP4) as the equal superposition of the eigenvector responsible for FP4 (FP5) with a symmetry-transformed version of it. Interestingly, the eigenvectors responsible for the simultaneous bifurcation from 
$\textrm {FP}1\rightarrow (\textrm {FP}4, \textrm {FP}5)$ at 
$Ra=6131$ are very similar to those responsible for the simultaneous bifurcation from 
$\textrm {FP}2\rightarrow (\textrm {FP}3, \textrm {FP}5)$ at 
$Ra=6058.5$. This can be explained as follows. The two simultaneous bifurcations occur at Rayleigh numbers which are close to each other and to 
$Ra=6056$, at which FP2 is formed via a supercritical circle pitchfork bifurcation from FP1. Fixed point FP2 inherits the spectrum of FP1, with the exception of the double eigenvalue responsible for the circle pitchfork. (Just above 
$Ra=6056$, this double eigenvalue becomes positive for FP1, whereas it splits into a zero and negative eigenvalue for FP2.) The other eigenvectors and eigenvalues of FP2 at 
$Ra=6058.5$ are close to those of FP1 at 
$Ra=6131$, including those shown in figure 4 which cause the simultaneous bifurcations. We do not show the eigenvectors of FP1 to avoid repetition.
It has been known since the mid-1980s (Swift Reference Swift1985) that 
$D_4$ symmetry leads to the simultaneous creation of non-symmetry-related branches. This has been applied to a number of situations, such as the simultaneous creation of standing and traveling waves (Knobloch Reference Knobloch1986; Borońska & Tuckerman Reference Borońska and Tuckerman2006; Reetz et al. Reference Reetz, Subramanian and Schneider2020). The application most relevant here is that of counter-rotating Taylor–Couette flow, in which spirals were first described in the classic paper by Taylor (Reference Taylor1923). The superposition of spirals of opposite helicity leads to a state called ribbons, much as the superposition of diagonal states produces the rectangular states in the current study. Exceptionally, ribbons were first predicted mathematically (Demay & Iooss Reference Demay and Iooss1984; Chossat & Iooss Reference Chossat and Iooss1994), setting off a quest to observe them experimentally, which was finally achieved by Tagg et al. (Reference Tagg, Edwards, Swinney and Marcus1989).
4. Periodic orbits
In this section, we explore four periodic orbits, PO1–PO4. Periodic orbits PO1–PO3 are created by a sequence of local bifurcations (i.e. bifurcations associated with a change in the real part of one or more eigenvalues/Floquet exponents): 
$\textrm {FP}3\rightarrow \textrm {PO}1\rightarrow \textrm {PO}2\rightarrow \textrm {PO}3$. PO1 and PO2 disappear in a global homoclinic and heteroclinic bifurcation, respectively, whereas the termination of PO3 is not discussed in this work. Periodic orbit PO4 bifurcates from and terminates on FP6 via Hopf bifurcations. The bifurcation diagram of figure 5(a) shows the six equilibria discussed in § 3 and the four periodic orbits to be discussed, whereas the periods of the limit cycles are shown in figure 5(b).
Figure 5. (a) Bifurcation diagram of fixed points (FPs) and periodic orbits (POs) and (b) periods of four periodic orbits in domain 
$[L_x, L_y, L_z] = [1, 8, 9]$. In (a), for each periodic orbit, we show two curves, the maximum and minimum of 
$\lvert \lvert \theta \lvert \lvert _2$ along an orbit. Periodic orbit PO1 appears via a Hopf bifurcation from FP3 at 
$Ra=6140$ (marked by a cyan cross) and undergoes a period-doubling bifurcation at 
$Ra=6154.7$ giving rise to PO2. Periodic orbit PO1 then undergoes a saddle-node bifurcation at 
$Ra=6157.97$ and disappears by meeting FP4 in a homoclinic bifurcation at 
$Ra_{hom}=6151.97$ at which its period diverges; see (b). Periodic orbit PO2 loses stability at 
$Ra = 6173.8$ where PO3 is created via a supercritical pitchfork bifurcation. The stability of PO2 changes multiple times along the branch for 
$6235< Ra<6255$, see details in figure 12. Periodic orbit PO2 then undergoes two closely spaced saddle-node bifurcations (at 
$Ra=6276$ and 
$6273.6$; see panel b) before terminating by meeting two symmetrically related versions of FP2 in a heteroclinic bifurcation at 
$Ra_{het}=6277.95$, at which its period diverges. Periodic orbit PO3 is continued until 
$Ra=6407.3$ (the range 
$6340< Ra<6407.3$ is not shown) and its period remains approximately constant. The apparent lack of smoothness in the curves representing PO2 and PO3 (in panel a around 
$Ra=6250$) corresponds to the overtaking of one temporal maximum or minimum of 
$\lvert \lvert \theta \lvert \lvert _2$ by another as 
$Ra$ is varied. Periodic orbit PO4 bifurcates from and terminates on FP6 at 
$Ra=6257.6$ and 
$Ra=6328.8$, and it is stable within 
$6257.6< Ra<6278$. In (a), solid and dashed curves signify stable and unstable states respectively, and the curves representing periodic orbits are slightly thicker than those of fixed points. The same colour code is used in (a,b). A schematic bifurcation diagram is shown in figure 15. Many other branches of equilibria and periodic orbits exist, which we have not followed nor shown on this diagram.
4.1. First periodic orbit (PO1)
4.1.1. Creation of PO1: Hopf bifurcation
Produced by a subcritical pitchfork bifurcation from FP2, FP3 is unstable at onset, but is stabilised by a saddle-node bifurcation at 
$Ra = 6008.5$ and then loses stability again at 
$Ra = 6140$ via a supercritical Hopf bifurcation that produces a periodic orbit PO1. PO1 inherits all of the spatial symmetries of FP3: 
$\langle {{\rm \pi} _y{\rm \pi} _{xz}, \tau (L_y/4,-L_z/4)} \rangle \sim D_4$ and, hence, no additional spatiotemporal symmetries are present. The S-shaped green curve in figure 5(a) contains the maximum (PO1
$_{max}$, above the cyan curve of FP3) and minimum (PO1
$_{min}$, below the cyan curve) values of 
$\lvert \lvert \theta \lvert \lvert _2$ over the period of each PO1 state. The period (
$T$) of PO1 increases smoothly before the saddle-node bifurcation at 
$Ra = 6157.97$. Prior to this, PO1 loses stability by undergoing a period-doubling bifurcation at 
$Ra=6154.7$ to PO2, which is discussed in § 4.2.1. The saddle-node bifurcation can be seen in both the maximum and minimum dashed green curves of figure 5(a) and leads to what we call the lower branch (because of its lower value of 
$\lvert \lvert \theta \lvert \lvert _2$).
By using the multi-shooting method with two to five shots, we have been able to continue the lower PO1 branch down in Rayleigh number to 
$Ra=6152.2041$, where the period of PO1 is very long: 
$T=955.4$ time units. We will show in § 4.1.2 that PO1 disappears via a homoclinic bifurcation, at which its period is infinite. Figures 6(a)–6(d) show snapshots of PO1 at 
$Ra = 6152.249$, on the lower branch. Among these snapshots, figures 6(a) and 6(b) capture the thinning and thickening of the rolls along the diagonal, with local waviness along the edge of the rolls. The waviness becomes weaker in figure 6(c) and finally, in figure 6(d) the edges are smoother and the roll widths almost uniform. All of the states in the cycle have a definite diagonal orientation. This implies that there exists another version of PO1 with the opposite diagonal orientation. The times at which these snapshots are taken are marked by stars in figures 6(e) and 6(f).
Figure 6. (a–d) The dynamics of PO1 (visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$) on the unstable lower branch at 
$Ra = 6152.249$ (
$Ra_{hom} = 6151.97$). Snapshot (d) converges to FP4 when used as an initial estimate for Newton solving. (e) Time series from DNS at 
$Ra =6152.249$ (
$T=900$), initialised by the unstable PO1 shown in (a) (red curve) and by FP4 with a small perturbation (black). The trajectory initialised by the unstable PO1 spends a long time near FP4 (
$250< t<800$). Both simulations converge to the stable PO1 branch (
$t>2500$) at this Rayleigh number. (f) Phase portrait illustrating the same data set as in (e). The plot shows the thermal energy input (
$I$) vs the viscous dissipation over energy input (
$D/I$). Fixed point FP4 (hollow blue circle) is located on the vertical line 
$D/I=1$, where energy dissipation and input are equal. The four red stars in (e,f) indicate the moments at which the snapshots (a–d) are taken. The same colour code is used in (e,f).
Figure 6(e) shows time series initialised with this unstable PO1 and also with a slightly perturbed FP4, at 
$Ra=6152.249$. Both of these runs eventually converge to another state: the stable upper branch of PO1, whose period 
$T = 161$ is much shorter than the period 
$T=900$ of the lower branch PO1. For 
$t< 1000$, the red curve remains close to FP4 during a large portion of the period. Figure 6(d) corresponds to the fourth star of 6(e), indicating via this projection that 6(d) is long-lived and very close to FP4. Indeed, we used figure 6(d) as the initial estimate for Newton's method to converge to FP4 at 
$Ra=6152.249$. However, figure 6(c), which only shows a transient at 
$Ra=6152.249$, resembles figure 2(e), which shows the converged FP4 at 
$Ra=6281$. We see from this that the diagonal orientation of FP4 becomes more prominent at higher Rayleigh numbers. Figure 6(f) shows a phase portrait from the same simulation as 6(e), using the thermal energy input 
$I$ and viscous dissipation 
$D$. There, FP4 is indicated as the hollow blue circle with 
$D/I=1$, showing that energy dissipation and input are equal. Near FP4, the dotted red curve looks continuous; this is due to the clustering of points near FP4.
4.1.2. Termination of PO1: homoclinic bifurcation
The close approach to FP4 implies that PO1 is close to a homoclinic cycle. We have verified that this closest approach is always to the same version of FP4 and not to another symmetry-related version. Thus, PO1 approaches a homoclinic, and not a heteroclinic cycle. A homoclinic cycle approaches a fixed point along one of its stable directions and escapes from it along one of its unstable directions. For this reason, we compute the eigenvalues and eigenvectors of FP4. Figure 7(a) shows the leading eigenvalues that we have computed at 
$Ra = 6152.249$, close to the global bifurcation point. The seven leading eigenvalues, all real, are 
$[\lambda _1, \lambda _2, \lambda _3, \lambda _4, \lambda _5, \lambda _6, \lambda _7] = [0.0212, 0.0208, 0.0026, 0, 0, -0.00017, -0.0034]$. We have set any eigenvalue whose absolute value is less than 
$10^{-7}$ to zero. Figure 7(a) shows other eigenvalues with smaller real parts as well and some of the eigenvalues are too close together to be distinguished. Certain eigenvalues of special significance are highlighted by coloured circles and their corresponding eigenvectors are shown in figures 7(c)–7(f).
Figure 7. (a) Leading eigenvalues at 
$Ra=6152.249$ of FP4: 
$[\lambda _1, \lambda _2, \lambda _3, \lambda _4, \lambda _5, \lambda _6, \lambda _7] = [0.0212, 0.0208, 0.0026, 0, 0, -0.00017, -0.0034]$. Eigenvalues 
$\lambda _2$ (escaping, red), 
$\lambda _{4,5}$ (neutral, green) and 
$\lambda _7$ (approaching, blue) are marked in colour. (b) The 
$L_2$-distance between each instantaneous flow field of PO1 and FP4 at 
$Ra=6152.249$, close to 
$Ra_{hom}=6151.97$. The evolution of PO1 (black curve) is exponential most of the time, with the escape from (red line) and approach to (blue line) FP4 governed by 
$\lambda _2$ and 
$\lambda _7$. (c–f) Four leading eigenmodes of FP4 at 
$Ra = 6152.249$, visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$: 
$e_2$, 
$e_{4}$, 
$e_{5}$ and 
$e_7$. The same colour bar is used in all plots.
The eigenvectors can be interpreted by considering FP4 and PO1, as depicted in figures 6(a)–6(d). There are two neutral directions, due to the continuous translation symmetry in the periodic directions. Eigenvalue 
$\lambda _4$ is zero and the corresponding eigenvector 
$e_4$, depicted in figure 7(d), is the neutral mode associated with 
$z$-translation (i.e. the 
$z$ derivative) of the roll-like FP4, very close to a z-translated version of figure 6(d). There must also be a marginal eigenvector corresponding to 
$y$-translation and indeed, 
$\lambda _5=0$ and we have verified numerically that 
$e_5$, depicted in figure 7(e), is the 
$y$ derivative of FP4. This is not immediately obvious, but note that for 
$z$ constant, the 
$y$ derivative of FP4 oscillates in sign and its maxima and minima are located along a diagonal. The green circle in figure 7(a) contains 
$\lambda _4$ and 
$\lambda _5$ (but also 
$\lambda _6$, whose decay rate is very small).
The other two eigenvectors shown in figures 7(c) and 7(f) are responsible for the approach to and escape from FP4. We have determined which eigenvalues are associated with approach and escape by comparing them with the observed approach and escape rates, and also by subtracting FP4 from the instantaneous flow fields and comparing the result to the eigenvectors. For the escaping dynamics of PO1 from FP4, the quantity (
$||\boldsymbol {x}(t) - \textrm {FP4}||_2$) increases exponentially at rate 
$\lambda _2=0.0208$. The corresponding eigenvector 
$e_2$ is shown in figure 7(c) and can be viewed as corresponding to widening and narrowing of the rolls. The approaching dynamics is characterised by 
$\lambda _7=-0.0034$. The corresponding eigenvector 
$e_7$, shown in figure 7(f), can be viewed as corresponding to translation in 
$y$. The portion of PO1 escaping FP4 along 
$e_2$ has been fit to the red line in figure 7(b). Although the rate of escape matches 
$\lambda _2$ closely, the approach rate only fits 
$\lambda _7$ over a short range of time (blue line). In figure 7(a), the red circle contains 
$\lambda _2$ (but also 
$\lambda _1$, which is very close to 
$\lambda _2$), and the blue circle encloses 
$\lambda _7$.
The eigendirection 
$e_7$ along which PO1 approaches FP4 is not the leading stable (least negative) one, as would be usual for a homoclinic orbit. This is because PO1 exists in the invariant symmetry-restricted subspace 
$\langle {\tau (L_y/4,-L_z/4)} \rangle$, to which 
$e_7$ also belongs. In contrast, 
$e_6$ (not shown), whose eigenvalue 
$\lambda_6$ is slightly less negative, has the opposite symmetry 
$\langle {\tau (L_y/4,L_z/4)} \rangle$. Note also that near-homoclinic orbits for which the rate of escape exceeds the rate of approach (i.e. here 
$|\lambda _2| >|\lambda _7|$) are unstable, as is already seen in the time series in figure 6(e). However, because our periodic orbits are computed using Newton's method and not time integration, we can calculate this periodic orbit despite its instability.
In addition, a homoclinic orbit bifurcating from a hyperbolic fixed point (which is the case for FP4) is structurally unstable, i.e. it exists for a single parameter value; see Kuznetsov (Reference Kuznetsov2004, Lemma 6.1) for a proof. Strictly speaking, FP4 is a relative hyperbolic fixed point, since it has zero eigenvalues along the directions of its continuous translation symmetries in 
$y$ and 
$z$, but the result applies to the evolution normal to these directions, i.e. with 
$y$ and 
$z$ phases fixed (Krupa & Melbourne Reference Krupa and Melbourne1995). Thus, the homoclinic cycle on which PO1 terminates is neither stable nor robust.
The closeness of some of the eigenvalues in figure 7(a) can be explained by the fact that the 
$y$ dependence of FP4 is extremely weak. If FP4 were entirely 
$y$-independent, like FP1, then eigenvectors would come in pairs, corresponding to a trigonometric dependence (analogous to 
$\textrm {sine}$ and 
$\textrm {cosine}$ eigenmodes) in 
$y$ with different phases, or to the choice of diagonal direction. Since the dependence in 
$y$ is weak, this is still approximately true in many cases. Eigenvalue 
$\lambda _1=0.0212$ is very close to 
$\lambda _2=0.0208$ and indeed eigenvector 
$e_1$ (not shown) resembles a 
$y$-shifted version of 
$e_2$. The near-neutral eigenvalues 
$\lambda _3=0.0026$ and 
$\lambda _6=-0.00017$ correspond to eigenvectors 
$e_3$ and 
$e_6$ (not shown), which resemble 
$e_5$ and 
$e_7$ but oriented in the opposite diagonal direction or, equivalently, reflected.
As 
$Ra$ approaches 
$Ra_{hom}$, PO1 approaches FP4 and the time spent near FP4 increases, until PO1 touches FP4 and acquires an infinite period in a homoclinic bifurcation. The period of PO1 is dominated by the time of approach to FP4, as shown in figure 7(b). This time can be estimated by the formula
\begin{equation} T \approx{-}\frac{1}{|\lambda_-|}\ln|Ra-Ra_{hom}| + c_T, \end{equation}
where 
$\lambda _-=\lambda _7=-0.0034$ is the rate of exponential approach to FP4, 
$Ra_{hom} = 6151.97$ and 
$c_T=533$ is a fitting constant. This asymptotic scaling law for 
$T$ as a function of 
$Ra$ was first derived in Gaspard (Reference Gaspard1990) and later used by various researchers including Meca et al. (Reference Meca, Mercader, Batiste and Ramírez-Piscina2004), Reetz et al. (Reference Reetz, Subramanian and Schneider2020) and Liu et al. (Reference Liu, Sharma, Julien and Knobloch2024). As shown in figure 8, we have fit the numerically computed periods of the states on the lower branch to this formula. Note that only the Rayleigh number range 
$6152.2< Ra<6152.45$ very close to 
$Ra_{hom}$ has been used for fitting and that we have extended the backwards continuation of PO1 to the lowest Rayleigh number possible (
$Ra=6152.2041$) within our numerical precision and ability.
Figure 8. Growth of the period of PO1 close to the global bifurcation point. PO1 undergoes a saddle-node bifurcation at 
$Ra = 6157.97$ where the lower branch appears. (a) The periods computed by numerical continuation and its logarithmic fit (see the text). (b) Here a logarithmic scale is used for 
$Ra-Ra_{hom}$, on which the period depends linearly. The red horizontal bar in (a,b) indicates the Rayleigh number range 
$6152.2< Ra<6152.45$ used for curve fitting.
Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) observed a periodic orbit produced by a Hopf bifurcation from a steady state; these are the solutions that we have called PO1 and FP3. Our bifurcation analysis agrees with their results. Extending their work, we have found that PO1 undergoes a saddle-node bifurcation and then terminates in a homoclinic bifurcation by meeting a new unstable fixed point, FP4.
4.2. Second periodic orbit (PO2)
4.2.1. Creation of PO2: period-doubling bifurcation and symmetry
Orbit PO2 bifurcates from PO1 in a period-doubling bifurcation at 
$Ra=6154.7$. At this value of Rayleigh number, its period (
$T=341$) is exactly twice that of PO1 (
$T=170.5$), as shown in figure 5(b). We have confirmed the threshold in two additional ways: at 
$Ra=6154.7$, the maxima and minima of 
$\lvert \lvert \theta \lvert \lvert _2$ of PO2 in the time series are extremely close in amplitude and frequency to those of PO1; and the power spectrum contains a very small component of the new frequency of PO2.
Orbit PO2 inherits all of the spatial symmetries of PO1: 
$\langle {{\rm \pi} _y{\rm \pi} _{xz}, \tau (L_y/4,-L_z/4)} \rangle \sim D_4$ and has, in addition, the spatiotemporal symmetry:
\begin{equation} (u,v,w,\theta)(x,y,z,t+T/2) = (u,v,w,\theta)(x,y-L_y/4,z,t). \end{equation}
Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) presented visualisations of PO2 in their figure 20 at 
$Ra=6250$ and noted that it satisfied (4.2). Figures 9(a)–9(d) show four snapshots of the temperature field of PO2 in which the spatiotemporal symmetry (4.2) of PO2 can clearly be seen. Figures 9(b) and 9(d) are very similar to each other and to two symmetry-related versions of the diamond-roll state, which we denote by FP2 and FP2
$^\prime \equiv \tau (L_y/4=2,0)\textrm {FP}2$. Between these instants, figures 9(a) and 9(c) show a wavy modulation of convection rolls along one of the diagonals. Since FP2 is 
$y$-reflection symmetric, there necessarily exists another version of PO2 in which the modulation occurs along the other diagonal.
Figure 9. (a–d) Snapshots of the dynamics of PO2 (visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$) at 
$Ra=6277.88$ near 
$Ra_{het} = 6277.95$. Snapshots (b,d) show states which are close to two symmetry-related versions of FP2 (figure 2c). (e) Time series from DNS at 
$Ra=6277.88$, initialised by the unstable PO2 shown in (a). The dynamics after 
$t\approx 1250$ becomes irregular and eventually terminates in chaos. (f) Phase space projection close to the global bifurcation point: shown are the PO2 at 
$Ra=6277.88$ and two symmetry-related FP2 states involved in the heteroclinic cycle. In (e,f), the four red stars indicate the moments where the snapshots (a–d) are taken and the two purple crosses mark the instants 
$t_{55}$ and 
$t_{412}$. In (f), the red arrows show the direction of the trajectory.
4.2.2. Termination of PO2: heteroclinic bifurcation and eigendirections
Figure 5(b) shows that although the period of PO2 decreases significantly with increasing 
$Ra$ until 
$Ra=6204.8$, it increases beyond that, eventually diverging. This implies that PO2 undergoes a global bifurcation. We have been able to continue PO2 until 
$(Ra,T)=(6277.88, 710.8)$ and we estimate the critical Rayleigh number for the global bifurcation to be approximately 
$Ra_{het} = 6277.95$.
The snapshots of figures 9(a)–9(d) are taken at 
$Ra=6277.88$, very close to 
$Ra_\textrm {het} = 6277.95$. PO2's alternating visits to FP2 and 
$\textrm {FP}2^\prime$ indicate that PO2 ends at a heteroclinic cycle between these two fixed points. In figure 9(e), we show a time series of PO2 at 
$Ra = 6277.88$, indicating the instants at which snapshots in figures 9(a)–9(d) are taken. It is clear that figures 9(b) and 9(d) belong to fairly long-lived plateaux. The global measurement 
$||\theta ||_2$ does not distinguish between FP2 and 
$\textrm {FP}2^\prime$, so we have plotted a phase portrait in figure 9(f), which represents each instantaneous flow field by its distance from each version of FP2. The phase portrait shows the clustering of points near FP2 and FP2
$^\prime$, confirming that PO2 is close to a heteroclinic cycle connecting these states.
The phase portrait in figure 9(f) also shows clustering of points around 
$(0.1, 0.1)$, corresponding to instants 
$t_{38}$ and 
$t_{394}$. This clustering suggests that the limit cycle might be approaching other fixed points. However, the time series in figure 9(e) does not show any other plateaux close to 
$t_{38}$ and 
$t_{394}$, and Newton's method did not converge to any new equilibria around the states from 
$t_{28}$ to 
$t_{55}$ and from 
$t_{384}$ to 
$t_{412}$. This remains true up to the highest Rayleigh number (or, equivalently, the longest period) of PO2 reached by our numerical continuation. We conclude that this heteroclinic cycle contains no other fixed points.
The dynamics along which PO2 approaches and escapes from FP2 can be described by eigenvalues and eigenvectors of FP2. Figure 10(a) shows the leading eigenvalues of FP2 at 
$Ra = 6277.88$, computed by Arnoldi iteration, and which are 
$[\lambda _{1,2}, \lambda _{3,4}, \lambda _{5,6}, \lambda _{7,8}, \lambda _{9,10}] = [0.031, 0, -0.00019, -0.00788, -0.0138]$. We previously saw that for FP4, the eigenvalues are approximately double (see figure 7a) due to the approximate symmetries of FP4. Here, FP2 has exact reflection symmetries leading to eigenvalues which are exactly double.
Figure 10. (a) Leading eigenvalues of FP2 at 
$Ra=6277.88$. The 10 leading eigenvalues are real and double: 
$[\lambda _{1,2}, \lambda _{3,4}, \lambda _{5,6}, \lambda _{7,8}, \lambda _{9,10}]=[0.031, 0, -0.00019, -0.00788, -0.0138]$. Eigenvalue 
$\lambda _{1,2}$ (escaping, red), 
$\lambda _{3,4}$ (neutral, green) and 
$\lambda _{9,10}$ (approaching, blue) are marked in colour. (b) The 
$L_2$-distance between each instantaneous flow field of PO2 and FP2 (and 
$\textrm {FP}2^\prime$) at 
$Ra=6277.88$, close to 
$Ra_{het}$. The dynamics of PO2 is exponential for most of the cycle (black and cyan curves). The approaching (blue line) and escaping (red line) dynamics of PO2 with respect to FP2 are shown and are governed by two eigenvalues of FP2. (c–f) Four leading eigenmodes of FP2 at 
$Ra=6277.88$, visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$: 
$e_1$, 
$e_3$, 
$e_4$ and 
$e_9$. The same colour bar is used in all plots.
The two neutral eigenmodes due to the continuous translation symmetries are 
$e_3$, corresponding to the 
$z$ derivative of FP2 and shown in figure 10(d) and 
$e_4$, corresponding to its 
$y$ derivative and shown in figure 10(e). The green circle in figure 10(a) contains 
$\lambda _{3,4}$, but also 
$\lambda _{5,6}$. Figure 10(c) shows the escaping eigenmode 
$e_1$, which is responsible for choosing the diagonal orientation of PO2. Looking at figures 9(a)–9(d), this is not obvious, but we have verified that subtracting FP2 from instantaneous temperature fields in the escaping phase of PO2 yields a field resembling 
$e_1$. Moreover, eigenvalues 
$\lambda _{1,2}=0.031$ capture well the escape rate from FP2, as shown in figure 10(b). Eigenmode 
$e_2$, with the same eigenvalue, is related to 
$e_1$ by reflection symmetry, as shown in figures 4(a)–4(b) for 
$Ra=6056$. The direction in which PO2 approaches FP2 is 
$e_9$, again confirmed by subtracting FP2 from the appropriate flow field in PO2, and 
$\lambda _{9,10}$ closely approximate the decay rate to FP2 shown in figure 10(b). The direction in which PO2 approaches the equilibrium is again not its leading stable eigendirection, and for the same reason as for PO1: PO2 remains within the symmetry group 
$\langle {{\rm \pi} _y{\rm \pi} _{xz}, \tau (L_y/4,-L_z/4)} \rangle$, to which 
$e_9$ belongs, but not eigenmodes 
$e_{5,6}$ and 
$e_{7,8}$ (not shown). Since 
$|\lambda _{1,2}|>|\lambda _{9,10}|$, the heteroclinic cycle is unstable, which is confirmed by the chaotic behaviour in the time series in figure 9(e) after 
$t\approx 1250$.
For FP2, since the eigenvalues are double, the eigenspace corresponding to each is 2-D; the eigenvectors that play the roles mentioned above (
$y$-translation, 
$z$-translation, escape and approach) must be selected as linear combinations of the two arbitrary eigenvectors returned by the Arnoldi method. By differentiating FP2 in 
$y$ and 
$z$ and by subtracting FP2 from the instantaneous flow fields during approaching or escaping phases, we have been able to choose the appropriate eigenvector in each case.
4.2.3. Robust heteroclinic cycle and 1:2 resonance
We now wish to show that the heteroclinic cycle that PO2 approaches is robust (also called structurally stable), i.e. that it exists over a parameter range rather than only at a single point. We have confirmed by numerical experiments that varying slightly the Rayleigh number does not affect the two transitions, also called half-cycles, 
$\textrm {FP}2\rightarrow \textrm {FP}2^\prime$ and 
$\textrm {FP}2^\prime \rightarrow \textrm {FP}2$. More rigorously, we list here the three conditions (Krupa Reference Krupa1997) that are required for a heteroclinic cycle between two fixed points, here FP2 and 
$\textrm {FP}2^\prime$, to be robust.
(i) There exist two invariant subspaces
$S$ and $S^\prime$ such that FP2 is a saddle (sink) and FP2$^\prime$ is a sink (saddle) for the flow restricted to subspace $S$ ($S^\prime$).(ii) There exist saddle-sink connections
$\textrm {FP}2\rightarrow \textrm {FP}2^\prime$ in $S$ and $\textrm {FP}2^\prime \rightarrow \textrm {FP}2$ in $S^\prime$.(iii) There exists a symmetry relation between the two fixed points.
Item (iii) is satisfied by definition: we have set FP2
$^\prime \equiv \tau (L_y/4,0)$FP2. For items (i) and (ii), we define 
$S$ and 
$S^\prime$ to be the fixed-point subspaces of two conjugate subgroups:
\begin{equation} \left.\begin{gathered} S \equiv {\rm Fix}|_{\langle{{\rm \pi}_y{\rm \pi}_{xz}\tau(L_y/2,0), \tau(L_y/4,-L_z/4)}} \rangle,\\ S^\prime \equiv {\rm Fix}|_{\langle{{\rm \pi}_y{\rm \pi}_{xz}, \tau(L_y/4,-L_z/4)} \rangle}. \end{gathered}\right\} \end{equation}
We note that 
$e_1$, depicted in figure 11(b), is an unstable eigenvector of both FP2 and FP2
$^\prime$ and belongs to subspace 
$S$ but not to 
$S^\prime$. We define 
$e_1^\prime \equiv \tau (L_y/4,0)e_1$, shown in figure 11(d), which is also an unstable eigenvector of FP2 and FP2
$^\prime$, and which belongs to subspace 
$S^\prime$ but not to 
$S$. The 
$L_2$-inner product of these two eigenmodes 
$\langle e_1, e_1^\prime \rangle$ is zero, and so they are orthogonal.
Figure 11. (a,c) Fixed points FP2 and 
$\textrm {FP}2^\prime \equiv \tau (2,0)\textrm {FP}2$ at 
$Ra=6277.8$. (b,d) Unstable eigenmodes: 
$e_1$ of FP2 and 
$e_1^\prime \equiv \tau (2,0)e_1$ of FP2
$^\prime$. The wavenumbers of the equilibria and unstable eigenmodes in the 
$y$ direction suggest a 1 : 2 mode interaction. All snapshots are visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$.
We have carried out simulations within subspaces 
$S$ and 
$S^\prime$ by numerically imposing the corresponding symmetries. When we restrict the simulation to 
$S$, the unstable eigenmode 
$e_1^\prime$ is disallowed and escape from FP2 (a saddle in subspace 
$S$) must take place along 
$e_1$. This trajectory lands on FP2
$^\prime$, which is linearly stable (a sink in subspace 
$S$). Changing the imposed subspace from 
$S$ to 
$S^\prime$, eigenvector 
$e_1$ is disallowed and escape from FP2
$^\prime$ (a saddle in subspace 
$S^{\prime }$) occurs along 
$e_1^\prime$. This trajectory lands on the stable equilibrium FP2, which is a sink in subspace 
$S^{\prime }$. Similar arguments apply to the approaches to FP2
$^\prime$ and FP2 via eigenvectors 
$e_9^\prime$ and 
$e_9$, shown in figure 10(f), respectively. Thus, we have demonstrated items (i) and (ii), proving that the heteroclinic cycle is robust. These three conditions are also discussed in Reetz & Schneider (Reference Reetz and Schneider2020), together with an example of a robust heteroclinic cycle between two symmetrically related oblique-wavy-roll equilibrium states found in inclined layer convection system.
In addition to demonstrating that the heteroclinic cycle that emerges from PO2 and FP2 is robust, we discuss its origin. We first address why FP2 has unstable and stable eigenvectors of the form 
$e_1$ and 
$e_9$. We recall that FP1 is homogeneous in 
$y$, FP2 has a 
$y$ periodicity of 
$L_y/2=4$, and FP3, FP4 and FP5 have 
$y$-periodicity 
$L_y=8$. When FP2 is created, it inherits the eigenvectors of FP1, including those which lead from FP1 to FP4 and FP5 (with 
$y$-periodicity 
$L_y=8$). The existence of such eigenvectors for FP2 is confirmed by the bifurcations from it to FP3 and FP5, which also have 
$y$-periodicity 
$L_y=8$; see, for example, figure 4. Because these 
$L_y$-periodic eigenvectors are all associated with nearby bifurcations, their corresponding eigenvalues are necessarily among the leading ones of FP2 in this range of 
$Ra$. This is the scenario of 1:2 resonance, the normal form of which was derived by Armbruster et al. (Reference Armbruster, Guckenheimer and Holmes1988):
\begin{equation} \left.\begin{gathered} \dot{z_1} = \overline{z_1} z_2 + z_1 \left(\mu_1 + e_{11} |z_1|^2 + e_{12} |z_2|^2 \right), \\ \dot{z_2} ={\pm} z_1^2 + z_2 \left(\mu_2 + e_{21} |z_1|^2 + e_{22} |z_2|^2 \right), \end{gathered}\right\} \end{equation}
where 
$z_1$ and 
$z_2$ are complex amplitudes of modes with wavenumbers 1 and 2, 
$\mu _1, \mu _2$ are control parameters and 
$e_{11}, e_{12}, e_{21}, e_{22}$ are nonlinear coefficients. These authors have demonstrated that (4.4) has a solution which is a heteroclinic orbit over a finite range of parameter values. Heteroclinic orbits of this type have been observed in full fluid dynamical configurations by, e.g., Mercader, Prat & Knobloch (Reference Mercader, Prat and Knobloch2002), Nore et al. (Reference Nore, Tuckerman, Daube and Xin2003), Bengana & Tuckerman (Reference Bengana and Tuckerman2019) and Reetz & Schneider (Reference Reetz and Schneider2020).
We now recall from § 3.3 that, in addition to the diagonally oriented eigenmode 
$e_1$ and its translation- and reflection-related versions, FP2 also has eigenmodes of type 
$e_2\equiv (e_1+{\rm \pi} _y e_1)/\sqrt {2}$, shown in figures 4(c)–4(f), which have a reflection symmetry in 
$y$ and which we have called rectangular. The diagonal eigenvector 
$e_1$ is responsible for the bifurcation to FP3, whereas the rectangular eigenvector 
$e_2$ is responsible for the bifurcation to FP5. Perturbing FP2 along 
$e_2$ can lead to a rectangular periodic orbit that retains 
$y$-reflection symmetry, which is currently under investigation.
4.2.4. Stability of PO2
The PO2 is stable over a short interval: from its onset at 
$Ra=6154.7$ until 
$Ra=6173.8$, where it becomes unstable via a pitchfork bifurcation giving rise to another periodic orbit PO3, to be discussed next in § 4.3. Just before the global bifurcation at 
$Ra=6277.95$, PO2 undergoes two saddle-node bifurcations at 
$Ra=6276$ and then at 
$Ra=6273.6$; these bifurcations do not restabilise PO2. However, Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) observed PO2 at 
$Ra=6250$ via DNS, implying that PO2 should be stable at that Rayleigh number. In order to understand this, we computed the leading Floquet exponent of PO2 over a range of 
$Ra$ surrounding 6250.
The intriguing evolution of the stability of PO2 is presented in figure 12. The leading Floquet exponent 
$\lambda _1$ is real from 
$Ra=6173.8$ to 
$Ra=6237.6$: it increases monotonically from 
$Ra=6173.8$ to 
$Ra=6225$ (not shown), and then decreases monotonically to zero over the interval 
$6226< Ra<6237.6$. The leading real exponent is then superseded by a complex conjugate pair 
$\lambda _{1,2}$ whose real part, initially negative, becomes positive over the interval 
$6239< Ra<6246.32$. At 
$Ra \approx 6246.5$, the leading exponent 
$\lambda _1$ becomes real and negative, so that there is a small interval 
$6246.5< Ra<6252$ over which PO2 is stable. It is within this very short interval that PO2 was observed by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018). In a further effort to understand the stabilisation and subsequent destabilisation of PO2 in this region, we computed the Floquet eigenmode to the left (figure 12b) and right (figure 12c) of the stable region, but we were unsuccessful in gleaning any physical insight from these. (There necessarily exist new branches bifurcating at the values at which 
$\lambda _1$ or the real part of 
$\lambda _{1,2}$ cross zero, but finding and following these new branches are beyond the scope of the current work.)
Figure 12. (a) The real part of the leading Floquet exponents of PO2 as a function of Rayleigh number. From low to high Rayleigh number, the leading Floquet exponent decreases monotonically within 
$6226< Ra<6238.75$. At 
$Ra=6238.75$, one sees the formation of a complex conjugate pair which has a positive real part for 
$6239< Ra<6246.32$. For 
$6246.5< Ra<6252$, PO2 is stable, with stability lost for 
$Ra>6252$. The apparent non-smoothness of the curve at 
$(Ra,\lambda _1)\approx (6238.7,-0.00047)$ and 
$(6250,-0.002)$ is due to the crossover of competing leading Floquet exponents. The two blue circles indicate where (b,c) are taken. (b,c) Two leading unstable Floquet eigenmodes for (b) 
$6226< Ra<6237.6$ and (c) 
$Ra>6252$, visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$. The same colour bar is used in (b,c).
4.3. Third periodic orbit (PO3): pitchfork bifurcation
As mentioned in § 4.2.4, PO2 loses stability at 
$Ra=6173.8$ via a supercritical pitchfork bifurcation which creates PO3. The visual features of PO3 resemble those of PO2 near onset, but become much less regular at higher Rayleigh numbers, for instance at 
$Ra=6407.3$, depicted in figures 13(a) and 13(d). The PO3 has spatial symmetries 
$\langle {\tau (L_y/2,-L_z/2)} \rangle \sim Z_2$, and the spatiotemporal symmetry (4.2) inherited from PO2. This spatiotemporal symmetry can be seen by comparing figures 13(b) and 13(c), for instance; the direction of drift for PO3 is from left to right. Periodic orbit PO3 loses stability at 
$Ra = 6183$. The bifurcating Floquet exponent is real, suggesting a pitchfork bifurcation leading to the creation of a pair of symmetry-related periodic orbits. However, we did not find any stable periodic orbit via DNS in the vicinity of 
$Ra=6183$, implying that such a bifurcation would be subcritical. Because PO3 is only stable for 
$6173.8 \leqslant Ra \leqslant 6183$, it is not surprising that it was not observed by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018).
Figure 13. (a–d) Snapshots of the dynamics of PO3 (visualised via the temperature field on the 
$y$–
$z$ plane at 
$x=0$) at 
$Ra=6407.3$ showing turbulent and disordered switching rolls. (e) Floquet exponent spectrum of PO3 at 
$Ra=6407.3$ showing its 51 unstable Floquet exponents. (f) Time series from DNS at 
$Ra =6407.3$, initialised by the converged unstable PO3. The temporal transition from a periodic to chaotic signal occurs at 
$t\approx 400$. The red stars indicate the moments at which the snapshots (a–d) are taken.
We continued PO3 until 
$Ra = 6407.3$, considerably into the chaotic regime (
$Ra>6300$) mentioned by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018). (The range 
$6340< Ra<6407.3$ is not included in figure 5.) Parametric continuation of PO3 for 
$Ra>6350$ was computationally challenging, probably due to the fact that the orbit is very unstable in this Rayleigh number range; see the discussion of the numerical convergence of the iterative Newton algorithm in Sánchez et al. (Reference Sánchez, Net, García-Archilla and Simo2004) and Reetz et al. (Reference Reetz, Subramanian and Schneider2020). The spectrum of PO3 at 
$Ra=6407.3$ has more than 
$50$ unstable eigendirections with a wide range of frequencies, as illustrated in figure 13(e). Moreover, integrating the converged PO3 forward in time at 
$Ra = 6407.3$, the transition from a periodic to chaotic state is triggered after fewer than two periods of the orbit; see figure 13(f). Consequently, we stopped the forward Rayleigh number continuation at 
$Ra = 6407.3$ and do not discuss how PO3 terminates.
4.4. Fourth periodic orbit (PO4): Hopf bifurcations
A new periodic orbit PO4 begins and ends on the lower branch of FP6 via two Hopf bifurcations at 
$Ra=6257.6$ and 
$Ra=6328.8$, respectively. As might be expected and as shown in figure 14, PO4 is an oscillating version of FP6. Since PO4 preserves the two reflection symmetries of FP6 
$\langle {{\rm \pi} _y,{\rm \pi} _{xz}\tau (L_y/2,0)} \rangle$, PO4 has no additional spatiotemporal symmetries. The Hopf bifurcation terminating PO4 occurs very slightly before the saddle-node bifurcation that terminates FP6 at 
$Ra=6329$. The PO4 originates from FP6 within the short range over which FP6 is stable, and at 
$Ra=6278$, PO4 is destabilised by a secondary Hopf bifurcation. Thus, PO4 is stable for 
$6257.6< Ra<6278$, as shown in figure 5(a) and in the schematic figure 15. Its period increases smoothly and monotonically throughout its range of existence, shown in figure 5(b).
Figure 14. (a) Time series from DNS of PO4 at 
$Ra =6300$ (
$T=255.7$). The red stars indicate the moments at which the snapshots (b,c) are taken. (b,c) Visualisations of PO4 at 
$Ra=6300$, via the temperature field on the 
$y$–
$z$ plane at 
$x=0$.
Figure 15. Schematic bifurcation diagram summarising the origin and stability of all of the fixed points (FPs) and periodic orbits (POs) that we identified in the computational domain 
$[L_x, L_y, L_z] = [1, 8, 9]$. PF, SN, PD, H, Het and Hom are abbreviations for pitchfork, saddle-node, period-doubling, Hopf, heteroclinic and homoclinic bifurcations. The dotted vertical lines together with the solid red lines and numbers mark the Rayleigh numbers at which bifurcations occur. Solid and dashed horizontal lines signify stable and unstable states, respectively.
Based on the bifurcation diagram in figure 5(a), the family of branches FP5, FP6 and PO4, are unusual in leaving no trace of their existence beyond the disappearance of FP6 at 
$Ra=6329$. Two FP5 branches join and terminate at 
$Ra=6317.5$; two FP6 branches, themselves created from FP5, annihilate at 
$Ra=6329$; PO4, created from FP6, disappears at 
$Ra=6328.8$. When we add to this the disappearances of periodic orbits (PO1 and PO2) via global bifurcations, we see that of the six fixed points and four periodic orbits that arise from the bifurcation of 2-D rolls at 
$Ra=5707$, only four fixed points and one periodic orbit survive past 
$Ra=6329$. Clever & Busse (Reference Clever and Busse1995) commented about simplification in another phenomenon (drifting waves) that they observed in vertical convection: ‘Of course, this is not a physically realistic scenario since there are other bifurcation points on the branch of the steady solutions 
$\ldots$ But the return from a complex structure to a more simple one with increasing control parameter is a possibility that cannot be excluded a priori’.
5. Discussion, conclusions and outlook
We have numerically investigated vertical thermal convection in the domain 
$[L_x, L_y, L_z] = [1, 8, 9]$, the configuration studied by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018), for Rayleigh number up to 
$Ra\approx 6400$. In this Rayleigh number range, the system exhibits various spatiotemporally organised flow patterns and weak turbulence. Using the computational power of parallelised numerical continuation based on matrix-free Newton methods, we have computed invariant solutions, more specifically fixed points, periodic orbits, and homoclinic and heteroclinic orbits.
We have situated all known solutions in the context of a bifurcation diagram. The diagrams shown in figures 2 and 5 are presented in schematic form in figure 15. This diagram contains the names of the states and the bifurcations between them, along with their precise thresholds, and emphasises the complexity of the bifurcation scenario. As was the case for Zheng et al. (Reference Zheng, Tuckerman and Schneider2024), all of the solution branches that we have found here are connected directly or indirectly to the laminar branch. This is not always so: our ongoing investigation has revealed branches which arise via saddle-node bifurcations and seem to be unconnected to the laminar state; see also figure 3 of Reetz et al. (Reference Reetz, Subramanian and Schneider2020).
Compared with the narrow domain 
$[1, 1, 10]$ presented in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024), the critical Rayleigh number for the primary instability of four spanwise-independent co-rotating rolls (called FP1 in both papers) in the spanwise-extended domain is only slightly lower. This is due to the slight reduction in the vertical length from 
$L_z=10$ to 
$L_z=9$ or, equivalently, from 
$\lambda =2.5$ to 
$\lambda =2.25$ in the primary roll wavelength. However, secondary and tertiary branches exist at much lower Rayleigh numbers for the 
$[1, 8, 9]$ domain than for the 
$[1,1,10]$ domain, since the larger domain accommodates a wider variety of spanwise-varying patterns.
We observe complicated bifurcation scenarios involving both spatial and temporal aspects. Spatially, parametric continuation reveals two types of branches. One set of branches consists of states which are aligned with the periodic directions 
$y$ and 
$z$: FP1 (2-D rolls), FP2 (diamond rolls), FP5 (moustache rolls) and the closely related FP6. The other set of branches consists of states which are oriented diagonally: FP3 (thinning rolls) or the similar FP4. We observed two instances of simultaneous bifurcation to branches of states with different symmetries. We were able to explain this otherwise non-generic phenomenon as the breaking of 
$D_4$ symmetry of the parent branches FP1 and FP2. (In this highly symmetric geometry, 
$D_4$ symmetry is a subgroup of the full symmetry groups of FP1 and FP2.) We confirmed this by computing and comparing the eigenvectors responsible for the simultaneous bifurcations.
Temporally, by following certain periodic orbit branches, PO1 and PO2, far from their onset via Hopf and period-doubling bifurcations, we have identified homoclinic and heteroclinic bifurcations that terminate these periodic-orbit branches. The fixed points at which these orbits spend an increasingly long time are aligned with the 
$y$ and 
$z$ axes (FP2), or nearly so (FP4), whereas the excursions are to diagonal states. Thus, these periodic orbits and global bifurcations can also be seen as a manifestation of competition between aligned and diagonal states. Although this is well understood from a mathematical group-theoretic viewpoint, there may exist some physical or phenomenological interpretation of when and why aligned or diagonal states are favoured. Another type of competition that we observe is between wavelengths: the heteroclinic orbit from FP2 can be interpreted as resulting from competition or interaction between states with wavenumbers 1 (
$y$ wavelength 
$L_y$) and 2 (
$y$ wavelength 
$L_y/2$). Indeed, the 1:2 mode interaction is a classic scenario leading to a robust heteroclinic orbit (Armbruster et al. Reference Armbruster, Guckenheimer and Holmes1988; Mercader et al. Reference Mercader, Prat and Knobloch2002; Nore et al. Reference Nore, Tuckerman, Daube and Xin2003; Bengana & Tuckerman Reference Bengana and Tuckerman2019; Reetz & Schneider Reference Reetz and Schneider2020).
The highest Rayleigh number that we have studied is 
$Ra=6407$, 
$12.3\,\%$ above the onset of convection (FP1 at 
$Ra=5707$). Even in this relatively small range of 
$Ra$, we have found a large variety of branches and bifurcation scenarios and there are certainly more to be discovered and analysed. In particular, primary bifurcations from the base state can lead to secondary states containing spanwise-independent (or 2-D) co-rotating rolls of many other wavelengths (all unstable at onset). These 2-D rolls can also undergo secondary, tertiary and global bifurcations. The increasing number of branches with Rayleigh or Reynolds number is a general feature of the Navier–Stokes and Boussinesq equations. However, branches can also disappear by the same types of local bifurcations that create them, and periodic orbits can be destroyed by global bifurcations, both of which occur in our configuration.
It has been conjectured that trajectories in chaotic and turbulent flows spend a substantial amount of time visiting unstable periodic orbits that are linked via their stable and unstable manifolds (Cvitanović & Eckhardt Reference Cvitanović and Eckhardt1991; Kawahara & Kida Reference Kawahara and Kida2001; Suri et al. Reference Suri, Kageorge, Grigoriev and Schatz2020; Crowley et al. Reference Crowley, Pughe-Sanford, Toler, Krygier, Grigoriev and Schatz2022). Computing unstable periodic orbits and understanding the bifurcations which produce and link them are thus relevant to better understanding and statistical measures of turbulent flows (Clever & Busse Reference Clever and Busse1995; Graham & Floryan Reference Graham and Floryan2021). In particular, reconstructing turbulence statistics using periodic orbits was explored by Chandler & Kerswell (Reference Chandler and Kerswell2013), using around 50 periodic orbits embedded in turbulent 2-D Kolmogorov flow; see also Cvitanović (Reference Cvitanović2013). Extending this approach to 3-D turbulent thermal convection is one of the objectives of our future research.
Acknowledgements
We thank S. Azimi, F. Reetz and O. Ashtari for fruitful discussions. We are grateful to D. Barkley, P. Beltrame, P. Gaspard and E. Knobloch for their insights on heteroclinic cycles.
Funding
This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant no. 865677).
Declaration of interests
The authors report no conflict of interest.
Open access
