Natural convection in a vertical channel. Part 3.
Bifurcations of many (additional) unstable periodic
orbits and their dynamical relevance

Zheng Zheng, Laurette S Tuckerman, Tobias M Schneider

► To cite this version:

Zheng Zheng, Laurette S Tuckerman, Tobias M Schneider. Natural convection in a vertical channel. Part 3. Bifurcations of many (additional) unstable periodic orbits and their dynamical relevance. Journal of Fluid Mechanics, 2025, 1022, pp.A42. 10.1017/jfm.2025.10800. hal-05349853

HAL Id: hal-05349853

https://hal.science/hal-05349853v1

Submitted on 5 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Natural convection in a vertical channel. Part 3.
## Bifurcations of many (additional) unstable periodic orbits and their dynamical relevance

Zheng Zheng¹ ID, Laurette S. Tuckerman² ID and Tobias M. Schneider¹ ID

¹Emergent Complexity in Physical Systems Laboratory (ECPS), École Polytechnique Fédérale de Lausanne, Lausanne CH 1015, Switzerland

²Physique et Mécanique des Milieux Hétérogènes (PMMH), CNRS, ESPCI Paris, PSL University, Sorbonne Université, Université de Paris, 75005 Paris, France

**Corresponding author:** Zheng Zheng, zheng.zheng@epfl.ch

(Received 7 April 2025; revised 30 September 2025; accepted 3 October 2025)

Vertical thermal convection exhibits weak turbulence and spatio-temporally chaotic behaviour. For this configuration, we report seven new equilibria and 26 new periodic orbits. These orbits, together with four previously studied in Zheng *et al*. (J. Fluid Mech., 2024b, vol. 1000, p. A29) bring the number of periodic-orbit branches computed so far to 30, all solutions to the fully nonlinear three-dimensional Navier-Stokes equations. These new and unstable invariant solutions capture intricate spatio-temporal flow patterns including straight, oblique, wavy, skewed and distorted convection rolls, as well as bursts and defects. These interesting and important fluid mechanical processes in a small flow unit are shown to also appear locally and instantaneously in a chaotic simulation in a large domain. Most of the solution branches show rich spatial and/or spatio-temporal symmetries. The bifurcation-theoretic organisation of these solutions is discussed; the bifurcation scenarios include Hopf, pitchfork, saddle-node, period-doubling, period-halving, global homoclinic and heteroclinic bifurcations, as well as isolas. Furthermore, these orbits are shown to be able to reconstruct statistically the core part of the attractor, so that these results may contribute to a quantitative description of transitional fluid turbulence using periodic orbit theory.

**Key words:** convection in cavities, bifurcation, pattern formation

## 1. Introduction

Thermal convection plays a crucial role in a wide range of natural processes, including in atmospheric and oceanic circulation, in mantle convection, as well as in engineering

© The Author(s), 2025. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
applications such as cooling and heat transfer (Kaushika & Sumathy 2003; Arici, Karabay & Kan 2015). A better fundamental understanding of convection is thus essential for tackling a diverse array of academic and industrial challenges (Bodenschatz, Pesch & Ahlers 2000; Ahlers, Grossmann & Lohse 2009; Lohse & Shishkina 2024). One particular case is vertical convection, the motion of a fluid bounded by two vertical walls maintained at different temperatures, driven by both buoyancy and shear. Like the well-known and well-studied Rayleigh-Bénard convection, vertical convection is an ideal system for studying pattern formation. Recent experimental, numerical and theoretical studies have continued to refine our understanding of the stability, transition, turbulence and heat transport mechanisms which are governed by the deterministic flow equations; these make vertical convection a subject of ongoing interest in fluid mechanics research. We refer readers to the introductions of the two previous papers in this series (Zheng, Tuckerman & Schneider 2024a, 2024b) for a more complete literature review of this field.

Using the methodology described in Zheng *et al.* (2024a, 2024b), we consider vertical convection as a (very) high-dimensional nonlinear dynamical system and employ the dynamical-systems-based approach to investigate the flow. This approach has been used to study transition to turbulence in various shear-dominated flows; see reviews in Kawahara *et al.* (2012), Graham & Floryan (2021) and references therein. (We sometimes use the word turbulence to refer to spatio-temporally chaotic dynamics instead of to a fully developed turbulent flow with an energy cascade across multiple spatial scales.) Building on ideas by Smale, Ruelle, Bowen and Sinai, turbulence is sometimes viewed as a chaotic walk through a forest of non-chaotic invariant solutions (particularly equilibria and periodic orbits) (Lanford 1982). While equilibria may reproduce characteristic features of the flow, they are time-independent and so the information contained within such solutions is limited. However, unstable periodic orbits are much more dynamically important and are believed to be transiently visited by a weakly turbulent flow, and to form the skeleton and building blocks of the chaotic dynamics of transitional turbulence (Cvitanović 1991; Kawahara & Kida 2001).

In our vertical convection system, the control parameters are the Rayleigh number (Ra), which is proportional to the temperature difference between the two walls, and the Prandtl number (Pr), which is the ratio between kinematic viscosity and thermal diffusivity. In addition to individual invariant solutions that are identified at fixed control parameters of the system, a bifurcation analysis via parametric continuations in one of the control parameters (Ra in this work) may reveal the bifurcation-theoretic origins of solutions and connections between them. A prominent example is in plane Couette flow; Reetz, Kreilos & Schneider (2019) constructed the first equilibrium solution underlying the self-organised oblique turbulent-laminar stripe pattern, and suggested that it emerges from the well-studied Nagata equilibrium (Nagata 1990). Many other examples in convective systems can be found in Borońska & Tuckerman (2010a, 2010b), Reetz & Schneider (2020a), Reetz, Subramanian & Schneider (2020b) and Zheng *et al.* (2024a, 2024b).

The present work follows two previous numerical investigations (Gao *et al.* 2018; Zheng *et al.* 2024b). Gao *et al.* (2018) have surveyed the flow regimes in a three-dimensional computational domain of size [$L_x, L_y, L_z$] = [1, 8, 9], depicted in figure 1, by systematically increasing the Rayleigh number from the onset of convection at Ra = 5707 to Ra ≈ 6300 (with Pr = 0.71 corresponding to air). Zheng *et al.* (2024b) used the same domain size, constructed invariant solutions and extended the upper limit to Ra ≈ 6400; a sequence of bifurcations was determined, and six equilibria and four time-periodic solutions were analysed in detail. In addition to these known solutions, we present here 33 new unstable invariant solutions, including seven equilibria and 26 periodic
<figure><img src="image_3.png" /><figcaption>Figure 1. Vertical convection cell with size $[L_x, L_y, L_z] = [1, 8, 9]$. The flow is bounded between two fixed walls at $x = \pm 0.5$ at which the flow is heated and cooled, respectively. We visualise the flow on the $y-z$ plane at $x = 0$ (dotted), from left to right as indicated by the eye and arrow. The laminar velocity $u_0(x) = \sqrt{Ra/Pr}(x/4 - x^3)/6 e_z$ and temperature $T_0(x) = x$ of this system are traced as an orange curve and a green line, respectively. Gravity, denoted by $g$, is in the vertical direction $e_z$.</figcaption></figure>

orbits, and we have extended the Rayleigh-number range to $Ra \approx 6650$. Even though the increase in $Ra$ in each paper may seem insignificant and negligible compared with fully turbulent convection, the emerging complexity of the bifurcation problem is indeed already overwhelming.

The new solutions that we will discuss are mainly found by the standard recurrent flow analysis which uses time-dependent simulations to locate states at which nearly periodic solutions or near recurrences are detected, and uses them as initial conditions for Newton solving. By construction, then, the new solutions are embedded in the trajectories followed by the flow. The new periodic orbits not only capture important dynamics of the transitional flow, but may also act as a basis to predict the statistical properties of the dynamics (Hopf 1948); these two points will be discussed in detail in § 5. For recent analysis of statistical descriptions based on periodic orbits in two-dimensional Kolmogorov flow, see for instance Chandler & Kerswell (2013), Cvitanović (2013), Lucas & Kerswell (2015) and Page et al. (2024). For alternative techniques for constructing initial estimates of unstable periodic orbits, see Page & Kerswell (2020), Redfern, Lazer & Lucas (2024) and Engel et al. (2025).

The rest of the manuscript is structured as follows. The numerical methods are summarised in § 2. We discuss in § 3 the new equilibria and in § 4 the new periodic orbits, with a focus on their bifurcation scenarios. Section 5 explores the dynamical relevance of the identified orbits and discusses a statistical description based on unstable orbits. Our conclusion is in § 6.

## 2. System, computation of invariant solutions and symmetries

As the numerical methods used in the current research are the same as those described in our preceding papers, we refer readers to Zheng *et al.* (2024a, 2024b) as well as to the introduction of Zheng (2025) for detailed descriptions of the governing equations, laminar base solutions, boundary conditions, symmetries of the system, computation (including parametric continuation and linear stability analysis) of invariant solutions, as well as numerical visualisations. Here, we will succinctly summarise the key ingredients.
The Oberbeck-Boussinesq equations that govern our vertical convection system are

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \left(\frac{\mathrm{Pr}}{\mathrm{Ra}}\right)^{1/2} \nabla^2 \mathbf{u} + \mathcal{T}\mathbf{e}_z, \quad (2.1a) $$

$$ \frac{\partial \mathcal{T}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathcal{T} = \left(\frac{1}{\mathrm{Pr} \mathrm{Ra}}\right)^{1/2} \nabla^2 \mathcal{T}, \quad (2.1b) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (2.1c) $$

where $\mathbf{u} = [u, v, w](x, y, z, t)$ and $\mathcal{T} = \mathcal{T}(x, y, z, t)$ are total velocity and temperature, respectively, and $p = p(x, y, z, t)$ is the pressure. We impose periodic boundary conditions (in $y$ and $z$), Dirichlet boundary conditions (in $x$) and a zero mean pressure gradient (in $y$ and $z$) integral constraint. We solve (2.1) with these boundary conditions by direct numerical simulation (DNS), using the ILC extension module of the Channelflow 2.0 code (Gibson et al. 2019). As in Zheng et al. (2024b), we use a computational domain of size $[L_x, L_y, L_z] = [1, 8, 9]$ and discretise it by $[N_x, N_y, N_z] = [31, 96, 96]$ Chebychev-Fourier-Fourier modes, see figure 1.

The symmetries of the system are

$$ \pi_y[u, v, w, \mathcal{T}](x, y, z) \equiv [u, -v, w, \mathcal{T}](x, -y, z), \quad (2.2a) $$

$$ \pi_{xz}[u, v, w, \mathcal{T}](x, y, z) \equiv [-u, v, -w, -\mathcal{T}](x, y, -z), \quad (2.2b) $$

$$ \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). \quad (2.2c) $$

Equations (2.2a)–(2.2c) define reflection in $y$, combined reflection of $x, z$ and temperature $\mathcal{T}$ and translation in $y$ and $z$, respectively. These symmetry operations generate the group $S \equiv \langle \pi_y, \pi_{xz}, \tau(\Delta y, \Delta z) \rangle \sim [O(2)]_y \times [O(2)]_{x,z}$, where $[O(2)]_y$ refers to reflections and translations in $y$, as in (2.2a) and (2.2c), respectively, while $[O(2)]_{xz}$ refers to reflections in $(\mathcal{T}, x, z)$ as in (2.2b) and translations in $z$ as in (2.2c). The laminar base flow has symmetry $[O(2)]_y \times [O(2)]_{x,z}$, but symmetries are broken at each bifurcation. We will discuss the symmetries of each solution in §3 and §4. In addition, symmetries are also used as a tool to find invariant solutions in this work, because if solutions have such symmetries, they are usually less unstable in the constrained symmetry subspace, and thus less difficult to find and converge.

Invariant solutions (equilibria and periodic orbits) are flow fields $\mathbf{x}^* \equiv [u, v, w, \mathcal{T}](x, y, z)$ satisfying

$$ \mathcal{G}(\mathbf{x}^*) = \sigma \mathcal{F}^T(\mathbf{x}^*) - \mathbf{x}^* = 0, \quad (2.3) $$

where $\sigma$ is a symmetry operator and $\mathcal{F}^T$ is the time-evolution operator integrating (2.1) from an initial state $\mathbf{x}^*$ over a finite time period $T$. The period $T$ is arbitrary for a steady solution, and is the period of a time-periodic solution. The periodic orbits that we will discuss in §4 are of three types. Periodic orbits are solutions which recur exactly after a period ($\sigma$ is the identity in (2.3)). Relative periodic orbits are orbits whose shortest recurrence occurs for a non-trivial symmetry operation, e.g. $\sigma \equiv \tau(\Delta y, \Delta z)$. Pre-periodic orbits are relative periodic orbits which recur exactly after some finite number of periods, see e.g. Cvitanovic et al. (2005) and Budanur & Cvitanović (2017). (That is, pre-periodic orbits are relative periodic orbits in which $\sigma$ does not contain translations by irrational multiples of $L_y$ or $L_z$.) In most of the later figures, we use $T$ to denote periods for periodic orbits, relative periods for relative periodic orbits and pre-periods for pre-periodic orbits. Whenever we discuss a period $T$ of an orbit in the text, we specify to which type of period it refers; in most cases, it is the shortest period for which (2.3) holds for some $\sigma$.

Invariant solutions are computed by the shooting-based Newton method, with initial guesses generated by a systematic recurrent flow analysis. The success rate for converging
to an invariant solution from one of these initial guesses is roughly 40%. These solutions are then parametrically continued in Rayleigh number to construct bifurcation diagrams. We define the temperature deviation $\theta \equiv T - T_0$ from the conductive state $T_0$ shown in figure 1. We use its $L_2$-norm $||\theta||_2$ (i.e. the square root of the integral of $\theta^2$ over the domain) to plot bifurcation diagrams throughout §3 and §4. The linear stability of converged states is evaluated by the Arnoldi algorithm. All of the 33 new solution branches are linearly unstable; they will be shown as solid curves in bifurcation diagrams. The thermal energy input (I) due to buoyancy and the dissipation (D) due to viscosity, both averaged over the domain, are used for phase portrait visualisations; we refer readers to §2.3 of Reetz & Schneider (2020a) for formulas. The flows (and eigenvectors) are visualised via their temperature fields $\theta$ on the $y-z$ plane at $x=0$, see figure 1. In order to avoid overcrowding the figures, we have omitted colour bars in all snapshots while insuring that all snapshots in a single figure share the same colour bar.

### 3. Unstable equilibria

In addition to the six equilibrium solutions (FP1–FP6) presented in Zheng et al. (2024b), we discuss here seven new unstable equilibria: FP7–FP13. These equilibria are relevant for the discussion on Hopf and global bifurcations of periodic orbits in §4. A bifurcation diagram including all of the new steady states as well as FP1, FP2 and FP4 (FP3, FP5 and FP6 are not plotted to avoid clutter) is depicted in figure 2(a). All of the equilibria are continued forward in Rayleigh number until at least $Ra = 6750$. Note that many other branches of equilibria (and periodic orbits) exist, which we have not found, followed or shown on these diagrams.

#### 3.1. Equilibria FP7–FP8

We begin our survey by briefly discussing FP1 and FP2. Equilibrium FP1 (two-dimensional rolls) is shown in figure 2(b) and contains four straight convection rolls of wavelength $\lambda_{FP1} = 9/4 = 2.25$ whose axes are oriented in the $y$-direction. Equilibrium FP2 bifurcates from FP1 at $Ra = 6056$. Equilibrium FP2, shown in figure 2(c), is called wavy rolls in Gao et al. (2018) and diamond rolls in Zheng et al. (2024b). (The list of generators for FP2 in (3.1) omits $\tau(0, L_z/2)$, contained in (3.1) of Zheng et al. (2024b), because it can be produced by the other generators and is hence redundant.) Equilibrium FP7, shown in figure 2(e), bifurcates from FP2 at $Ra = 6279.5$ in a supercritical pitchfork bifurcation, in which the $\pi_y$ reflection and fourfold translation (along both diagonals) symmetries are broken. Equilibrium FP8, shown in figure 2(f), bifurcates from FP7 in a supercritical pitchfork bifurcation at $Ra = 6282.9$, in which the $\pi_y\tau(0, L_z/2)$ symmetry is broken. Equilibrium FP8 gives rise to PO23 and PO24 in two Hopf bifurcations, see §4.7. The symmetry groups of FP1, FP2, FP7 and FP8 are

$$ 
\begin{aligned}
\text{FP1: } & \langle \pi_y, \tau(\Delta y, 0), \pi_{xz}, \tau(0, L_z/4) \rangle \sim [O(2)]_y \times [D_4]_{xz}; \\
\text{FP2: } & \langle \tau(L_y/2, 0), \pi_{xz}, \pi_y, \tau(L_y/4, -L_z/4) \rangle \sim D_2 \times D_4; \\
\text{FP7: } & \langle \tau(L_y/2, 0), \pi_y\pi_{xz}, \pi_y\tau(0, L_z/2) \rangle \sim D_2 \times Z_2; \\
\text{FP8: } & \langle \tau(L_y/2, 0), \pi_y\pi_{xz} \rangle \sim D_2.
\end{aligned} 
\quad (3.1) $$

#### 3.2. Equilibria FP9–FP12

Equilibrium FP9, shown in figure 2(g), bifurcates from the homogeneous unstable base flow at $Ra = 5941$ (not shown in figure 2a) in a supercritical pitchfork bifurcation. Equilibrium FP9 has four pairs of oblique but straight convection rolls of wavelength $\lambda_{FP9} = 2L_z/\sqrt{L_z^2/16 + L_y^2} \approx 2.166$ each, in the direction perpendicular to the rolls.
<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Bifurcation diagram of equilibria and (b–k) flow structures visualised via the midplane temperature field. In (a), all branches shown are unstable, with the exception of FP1 for Ra < 6056 and of FP2 for 6056 < Ra < 6058.5; on the right are two enlarged diagrams, zooming in on the FP2 → FP7 → FP8 and FP9 → FP10 → FP11 bifurcations. The solutions (b) FP1, (c) FP2 and (d) FP4 have been presented in Zheng *et al.* (2024b) and are shown with thinner curves in (a). The solution (e) FP7 bifurcates from FP2 at Ra = 6279.5; (f) FP8 bifurcates from FP7 at Ra = 6282.9. The solution (g) FP9 bifurcates from the unstable base state at Ra = 5941; (h) FP10 bifurcates from FP9 at Ra = 6360; (i) FP11 bifurcates from FP10 at Ra = 6369.2 and undergoes a saddle-node bifurcation at Ra = 6213.5; (j) FP12 bifurcates from FP9 at Ra = 6184. The solution (k) FP13 undergoes a saddle-node bifurcation at Ra = 6449 and both upper and lower branches exist at least until Ra = 6800.</figcaption></figure>
The oblique angle with respect to the y-direction is $γ = \arctan(0.25L_z/L_y) \approx 15.7°$. Because FP9–FP12 all share this oblique orientation, we introduce tilted coordinates

$$ \begin{pmatrix} y' \\ z' \end{pmatrix} = \begin{pmatrix} \cos \gamma & \sin \gamma \\ -\sin \gamma & \cos \gamma \end{pmatrix} \begin{pmatrix} y \\ z \end{pmatrix}, \qquad (3.2) $$

which are drawn in figure 2(g). In the tilted coordinates, we consider a virtual computational domain having length $L'_z = 4λ_{FP9} \approx 8.664$ in $z'$ and $L'_y \approx 15$ in $y'$. (The length $L'_y \approx 15$ is three times the wavelength corresponding to the prominent structure along $y'$ in FP10 and FP11, and four times the wavelength corresponding to the wavy structure in FP12. Introducing this length will be convenient for the description of symmetry groups.) In this tilted domain, FP9 has $O(2)$ symmetry in $y'$ and $D_4$ symmetry in $xz'$; see (3.3).

Equilibrium FP10, shown in figure 2(h), bifurcates from FP9 at $Ra = 6360$ in a supercritical pitchfork bifurcation. In this bifurcation, the $O(2)$ symmetry of FP9 along $y'$ is broken and succeeded by a discrete (twofold) translation. In $z'$, the fourfold translation must now be combined with a discrete (fourfold) translation in $y'$ in order to remain a symmetry of the flow. Finally, the two independent reflection symmetries are replaced by a single combined reflection $\pi_{y'}\pi_{xz'}$.

Equilibrium FP11, shown in figure 2(i), bifurcates from FP10 at $Ra = 6369.2$ in a subcritical pitchfork bifurcation. Equilibrium FP11 then undergoes a saddle-node bifurcation at $Ra = 6213.5$ and continues to exist at least until $Ra = 7000$. In going from FP10 to FP11, the spatial periodicity along $y'$ changes from $L'_y/2$ to $L'_y$, while other symmetries are retained. Equilibrium FP11 gives rise to PO14 in a Hopf bifurcation, see § 4.5.1 in which FP9 (and FP12 below) will also be relevant.

Equilibrium FP12, shown in figure 2(j), also bifurcates from FP9, in a supercritical pitchfork bifurcation at $Ra = 6184$. In $y'$, the $O(2)$ symmetry of FP9 is succeeded by a fourfold translation, the fourfold translation in $z'$ is retained, and the two reflection symmetries are replaced by the single combined reflection $\pi_{y'}\pi_{xz'}$.

The symmetry groups of FP9–FP12 are

$$ \begin{aligned} \text{FP9:} & \quad \langle \pi_{y'}, \tau(\Delta y', 0), \pi_{xz'}, \tau(0, L'_z/4) \rangle \sim [O(2)]_{y'} \times [D_4]_{xz'}; \\ \text{FP10:} & \quad \langle \pi_{y'}\pi_{xz'}, \tau(L'_y/2, 0), \tau(L'_y/4, L'_z/4) \rangle \sim [D_2]_{y'} \times [Z_4]_{xz'}; \\ \text{FP11:} & \quad \langle \pi_{y'}\pi_{xz'}, \tau(L'_y/4, L'_z/4) \rangle \sim [Z_2]_{y'} \times [Z_4]_{xz'}; \\ \text{FP12:} & \quad \langle \pi_{y'}\pi_{xz'}, \tau(L'_y/4, 0), \tau(0, L'_z/4) \rangle \sim [D_4]_{y'} \times [Z_4]_{xz'}. \end{aligned} \qquad (3.3) $$

### 3.3. Equilibrium FP13

Equilibrium FP13 is shown in figure 2(k) and exists beyond $Ra = 6800$, where we stopped the continuation; its bifurcation-theoretic origin remains unclear. In the Rayleigh-number range that we consider, FP13 undergoes one saddle-node bifurcation at $Ra = 6449$. Equilibrium FP13 is relevant for PO23 in § 4.7.1, and its symmetry group is

$$ \text{FP13: } \langle \pi_y \pi_{xz} \rangle \sim Z_2. \qquad (3.4) $$

## 4. Unstable periodic orbits

In this section, we discuss 26 newly identified unstable periodic orbits, RPO5–PPO30. Our naming convention is that the letters (e.g. R and P) describe the type of orbit, while the numbers (e.g. 5 and 30) identify different states based on the sequential order in which the orbits were found. Figure 3 includes all of the periodic orbits PO2–PPO30 (PO1–PO4 are discussed in Zheng et al. 2024b) only to give an impression of the complexity of the full bifurcation diagram. To explain the various branches and scenarios, separated bifurcation
<figure><img src="image_5.png" /><figcaption>Figure 3. Temperature norms (a) and periods (b) of periodic orbits. Orbits PO2–PO4 are discussed in detail in Zheng *et al.* (2024b). In (a), for each orbit, we show two curves, the maximum and minimum of ||*θ*||₂ along an orbit. All of RPO5–PPO30 are linearly unstable. The upper limit of (b) is set to *T* = 700, even though some orbits are continued to higher period. The bifurcation scenarios include Hopf, pitchfork, saddle-node, period-doubling, period-halving and global homoclinic/heteroclinic bifurcations and isolas. For more clarity, bifurcation diagrams for selected sets of orbits will be shown in figures 4, 5, 7, 8, 11, 15 and 20. The apparent lack of smoothness in some ||*θ*||₂ curves corresponds to the overtaking of one temporal maximum or minimum of ||*θ*||₂ by another as *Ra* is varied.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4. Temperature norms (a) and periods (b) of RPO13, RPO15, RPO26 and RPO28. Branch RPO13 bifurcates from and terminates on RPO18 (which is shown more completely in figure 5) in two period-doubling bifurcations. The bifurcation points are indicated by stars on the right plot. Branches RPO15, RPO26 and RPO28 begin and terminate at saddle-node bifurcations and form isolas.</figcaption></figure>

diagrams for smaller groups of periodic orbits will be shown in figures 4, 5, 7, 8, 11, 15 and 20. Given the complexity of all of the bifurcation diagrams, we recommend reading each diagram by first focusing on the plot of the temporal period and then comparing it with that of the temperature norm. The reason is that the quantity $||\theta||_2$ might sometimes be close for multiple orbits and for one orbit along the branch (due to saddle-node bifurcations), but the periods of the orbits are more distinct and thus lead to a better understanding of the bifurcation scenarios. In this work, we focus on Rayleigh numbers up to $Ra \approx 6650$, thus ~16.5 % above the onset of convection, even though some orbits can be continued to much higher Rayleigh numbers. We do not discuss if and how their branches end there.

The bifurcation scenarios explored include Hopf, pitchfork, saddle-node, period-doubling, period-halving and global homoclinic/heteroclinic bifurcations and isolas. Given the large number of orbits that we will discuss, this section is organised in terms of the symmetries of the orbits. The eight subsections below will discuss orbits identified in the following symmetry subspaces: fourfold translation along the domain diagonal with a non-commuting reflection: $\langle \pi_y \pi_{xz}, \tau(L_y/4, L_z/4) \rangle$ in §4.1; twofold translation with a commuting reflection: $\langle \pi_y \pi_{xz}, \tau(L_y/2, L_z/2) \rangle$ in §4.2; fourfold translation: $\langle \tau(L_y/4, L_z/4) \rangle$ in §4.3; twofold translation: $\langle \tau(L_y/2, L_z/2) \rangle$ in §4.4; threefold translation with a non-commuting reflection: $\langle \pi_y \pi_{xz}, \tau(L_y/3, L_z/3) \rangle$ in §4.5; fivefold translation with a non-commuting reflection: $\langle \pi_y \pi_{xz}, \tau(L_y/5, L_z/5) \rangle$ in §4.6; and single reflection: $\langle \pi_y \pi_{xz} \rangle$ in §4.7. Orbits without any spatial symmetry will be presented in §4.8. Table 1 provides a summary of these solutions in terms of their symmetries, bifurcation scenarios, section of coverage and figures in which they are shown.

## 4.1. Symmetry subspace: reflection with fourfold translation

Seven orbits identified in the symmetry subspace $\langle \pi_y \pi_{xz}, \tau(L_y/4, L_z/4) \rangle \sim D_4$ will be discussed in this subsection. Due to this imposed symmetry constraint, the dynamics of each of these seven orbits has a diagonal orientation and consists of diagonal excursions from more aligned states. We only show snapshots of RPO18 (figure 6) for illustration.

### 4.1.1. Orbit RPO13: period-doubling bifurcations

Orbit RPO13 was found at $Ra = 6285$. Forward and backward continuation in Rayleigh number reveals that RPO13 bifurcates from and ends at RPO18 in two period-doubling
<table><caption>Table 1. Summary of spatial symmetries and bifurcation scenarios of 30 periodic orbits found in domain [L<sub>x</sub>, L<sub>y</sub>, L<sub>z</sub>] = [1, 8, 9], with PO1–PO4 discussed in Zheng et al. (2024b). Abbreviations PF, SN, PD, PH, H and GB stand for pitchfork, saddle-node, period-doubling, period-halving, Hopf and global bifurcations.</caption><thead><tr><th>Spatial symmetry</th><th>Periodic orbit</th><th>Bifurcations</th><th>Discussed in</th><th>Shown in figure(s)</th></tr></thead><tbody><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>PO1</td><td>H, SN, GB</td><td>Zheng et al. (2024b)</td><td>5, 6, 7, 8, 15</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>PO2</td><td>PD, SN, GB</td><td>Zheng et al. (2024b)</td><td>5, 9, 10, 11, 15</td></tr><tr><td>⟨τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>PO3</td><td>PF</td><td>Zheng et al. (2024b)</td><td>5, 13, 15</td></tr><tr><td>⟨π<sub>y</sub>, π<sub>xz</sub>τ(L<sub>y</sub>/2, 0)⟩</td><td>PO4</td><td>H</td><td>Zheng et al. (2024b)</td><td>5, 14, 15</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO13</td><td>PD</td><td>§ 4.1.1</td><td>4</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO15</td><td>SN, isola</td><td>§ 4.1.2</td><td>4</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO17</td><td>SN, isola</td><td>§ 4.1.2</td><td>5</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO26</td><td>SN, isola</td><td>§ 4.1.2</td><td>4</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO28</td><td>SN, isola</td><td>§ 4.1.2</td><td>4</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO18</td><td>SN, GB</td><td>§ 4.1.3</td><td>5, 6</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO19</td><td>PD, SN</td><td>§ 4.1.4</td><td>7</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>PO6</td><td>SN, GB</td><td>§ 4.2.1</td><td>8, 9, 10</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>PPO7</td><td>SN, isola</td><td>§ 4.2.2</td><td>11, 12</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>RPO10</td><td>SN, GB</td><td>§ 4.2.3</td><td>8, 13</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>RPO25</td><td>PF, SN, PH</td><td>§ 4.2.4</td><td>7</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>RPO27</td><td>PF, SN</td><td>§ 4.2.5</td><td>5</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>RPO29</td><td>SN, GB</td><td>§ 4.2.6</td><td>8, 14</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>PPO30</td><td>PD</td><td>§ 4.2.7</td><td>11</td></tr><tr><td>⟨τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO12</td><td>SN, GB</td><td>§ 4.3.1</td><td>15, 16</td></tr><tr><td>⟨τ(L<sub>y</sub>/4, L<sub>z</sub>/4)⟩</td><td>RPO20</td><td>SN, GB</td><td>§ 4.3.2</td><td>15, 16</td></tr><tr><td>⟨τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>RPO5</td><td>SN, GB</td><td>§ 4.4.1</td><td>11, 17</td></tr><tr><td>⟨τ(L<sub>y</sub>/2, L<sub>z</sub>/2)⟩</td><td>RPO8</td><td>SN</td><td>§ 4.4.2</td><td>11, 17</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/3, L<sub>z</sub>/3)⟩</td><td>PO14</td><td>H, GB</td><td>§ 4.5.1</td><td>8, 18</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/5, L<sub>z</sub>/5)⟩</td><td>PPO16</td><td>GB</td><td>§ 4.6.1</td><td>8, 19</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/5, L<sub>z</sub>/5)⟩</td><td>RPO21</td><td>SN</td><td>§ 4.6.2</td><td>20</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>, τ(L<sub>y</sub>/5, L<sub>z</sub>/5)⟩</td><td>RPO22</td><td>SN</td><td>§ 4.6.2</td><td>20</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>⟩</td><td>PO23</td><td>H, SN, GB</td><td>§ 4.7.1</td><td>20, 21</td></tr><tr><td>⟨π<sub>y</sub>π<sub>xz</sub>⟩</td><td>PO24</td><td>H</td><td>§ 4.7.2</td><td>20</td></tr><tr><td>No spatial symmetry</td><td>PO9</td><td>SN</td><td>§ 4.8.1</td><td>11, 22</td></tr><tr><td>No spatial symmetry</td><td>RPO11</td><td>SN</td><td>§ 4.8.2</td><td>15, 23</td></tr></tbody></table>

bifurcations, at Ra = 6280.01 and 6285.12, as indicated by the stars in figure 4. (Orbit RPO18 will be discussed in § 4.1.3.)

### 4.1.2. Orbits RPO15, RPO17, RPO26 and RPO28: saddle-node bifurcations and isolas

We identify four isolas in this symmetry subspace. As the name implies, they do not bifurcate from any other states but only undergo several saddle-node bifurcations to turn back in Ra, as evidenced by figure 4 for RPO15, RPO26, RPO28, and by figure 5 for RPO17. One more isola with slightly different symmetry will be discussed in § 4.2.2.

### 4.1.3. Orbit RPO18: saddle-node and global bifurcations

We found RPO18 at Ra = 6350 and continued it forwards up to Ra = 6686. Continuing backwards, RPO18 undergoes several saddle-node bifurcations and finally terminates in a global bifurcation by meeting FP2. The dynamics of RPO18 at Ra = 6277.958, the lowest Rayleigh number that we have reached for this state, is shown in figure 6(a–e). Orbit RPO18 resembles PO2, as presented in Zheng et al. (2024b): it has a clear oblique
<figure><img src="image_7.png" /><figcaption>Figure 5. Temperature norms (a) and periods (b) of RPO17, RPO18 and RPO27. The RPO17 branch forms an isola. Orbit RPO18 bifurcates from FP2 in a global homoclinic bifurcation at Ra = 6277.96 and continues to exist up to at least Ra = 6686. Orbit RPO27 is generated from RPO18 in a pitchfork bifurcation at Ra = 6279.7 and continues to exist up to at least Ra = 6650.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Dynamics of RPO18 at Ra = 6277.958 (close to the global bifurcation point) with relative period T = 476.31. (a–e) Snapshots of the midplane temperature field. (f) Time series from DNS. The five red stars indicate the moments at which the snapshots (a–e) are taken.</figcaption></figure>

orientation. The global bifurcation is evidenced by the plateau in the time series in figure 6(f) and the divergence of its period in figure 5.

We have determined the eigenvector along which RPO18 approaches and escapes from FP2 by computing the leading eigenvalues of FP2 at Ra = 6277.958 within the symmetry subspace $\langle \pi_y \pi_{xz} , \tau(L_y/4, L_z/4) \rangle$, as was done in Zheng *et al.* (2024b). We find that RPO18 escapes from FP2 along eigendirection $e_1$ associated with eigenvalue $\lambda_1 = 0.031285$ and approaches FP2 via eigendirection $e_2$ associated with eigenvalue $\lambda_2 = -0.0138$. These eigenvectors turn out to be the same or symmetrically related versions of those which are responsible for PO2 escaping from and approaching FP2. We do not show these to avoid repetition and refer readers to figures 9 and 10 of Zheng *et al.* (2024b) for
<figure><img src="image_9.png" /><figcaption>Figure 7. Temperature norms (a) and periods (b) of PO2, RPO19 and RPO25. Branch RPO19 bifurcates from and ends on PO2 (discussed in Zheng et al. 2024b) at Ra = 6252 and Ra = 6274 in two period-doubling bifurcations. Branch RPO25 bifurcates from RPO19 at Ra = 6260.5 in a pitchfork bifurcation, undergoes saddle-node bifurcations and terminates in a period-halving bifurcation (marked by PH) on another branch that is not shown or studied in this paper.</figcaption></figure>

details. Interestingly, the global bifurcations of PO2 and RPO18 occur at almost the same Rayleigh number, and we will see in §§ 4.2.1 and 4.2.6 that FP2 plays a role in other global heteroclinic and homoclinic bifurcations.

### 4.1.4. Orbit RPO19: period-doubling and saddle-node bifurcations

As shown in figure 7, RPO19 bifurcates from and terminates on PO2 in two period-doubling bifurcations. As discussed in § 4.2 of Zheng et al. (2024b), PO2 has a spatio-temporal symmetry and contains two pre-periodic orbits, each of whose period is half that of PO2. The periods of PO2 that we show in figure 7 (and figure 3b) correspond to full periods or twice the pre-periods. Orbit RPO19 also undergoes two saddle-node bifurcations at Ra ≈ 6254.6, which are difficult to see in the figure.

### 4.2. Symmetry subspace: reflection with twofold translation

Seven orbits identified in the symmetry subspace $\langle \pi_y \pi_{xz}, \tau(L_y/2, L_z/2) \rangle \sim D_2$ will be discussed in this subsection. Similarly to § 4.1, this imposed symmetry leads the dynamics of most of the orbits to acquire a diagonal orientation. We will show snapshots of PO6 (diagonal), PPO7 (non-diagonal) and RPO10 (diagonal), for illustration.

#### 4.2.1. Orbit PO6: saddle-node and global bifurcations

We first observed PO6 at Ra = 6280. Orbit PO6, which has the spatio-temporal symmetry

$$ (u, v, w, \theta)(x, y, z, t + T/2) = \pi_y(u, v, w, \theta)(x, y + L_y/2, z, t), \quad (4.1) $$

should be understood as a pre-periodic orbit and be properly called PPO6. Exceptionally, we use PO6 here for consistency with PO2 in Zheng et al. (2024b) whose dynamics closely resembles PO6. The period of PO6 shown in figure 3(b) and figure 8 is twice its pre-period. Continuing PO6 backwards, it approaches a heteroclinic cycle linking two symmetrically related versions of FP2. From Ra = 6280 to Ra = 6218.6, the period of PO6 increases monotonically to T = 1069.1 (we have restricted the range in figure 3(b) to T ≤ 700 to better show the periods of other orbits). At the global bifurcation point at slightly lower Rayleigh number, the period should diverge; compare with, for example, the slopes of PO2, PO14, PPO16, RPO18 and RPO29. We have been unable to continue PO6 further due to limits on numerical precision.
<figure><img src="image_10.png" /><figcaption>Figure 8. Temperature norms (a) and periods (b) of PO6, RPO10, PO14, PPO16 and RPO29. Branch PO6 approaches a heteroclinic cycle linking two symmetrically related versions of FP2 in a global bifurcation at $Ra \approx 6218.6$, at which its period diverges. At higher Rayleigh numbers, PO6 undergoes saddle-node bifurcations and continues to exist at least until $Ra = 6615$. Branch RPO10 possibly bifurcates from FP4 in a global bifurcation at $Ra \approx 6298.7$ and continues to exist at least until $Ra = 6650$. Branch PO14 bifurcates from FP11 in a Hopf bifurcation and terminates in a global bifurcation by meeting FP9. Branch PPO16 is created from FP9 in a global bifurcation at $Ra \approx 6240.6$ and continues to exist until at least $Ra = 6656.5$. Branch RPO29 bifurcates from FP4 at $Ra \approx 6274.14$ and terminates on FP2 at $Ra \approx 6402$ in two global bifurcations.</figcaption></figure>

The dynamics of PO6 at $Ra = 6218.6$, close to the global bifurcation point, is shown in figure 9(a-e). Orbit PO6 resembles RPO18 above. The dynamics of each half-period of PO6 has fourfold translation symmetry along one of the diagonals, either $\langle \tau(L_y/4, L_z/4) \rangle$ or $\langle \tau(L_y/4, -L_z/4) \rangle$; compare figures 9(b) and 9(e). The only symmetry possessed instantaneously by all members of PO6 is $\langle \tau(L_y/2, L_z/2) \rangle$. The heteroclinic bifurcation is evidenced by the plateaus in the time series in figure 9(f) as well as by the phase space projection in figure 9(g). The eigenvectors along which PO6 approaches and escapes from FP2 (and FP2' $\equiv \pi_y \tau(L_y/2, 0)$) are shown in figure 10, together with the phase of FP2 (and FP2').

To show that PO6 approaches a robust heteroclinic cycle as is the case for PO2, we identify two subspaces within the symmetry space of FP2

$$ 
\begin{aligned}
S &\equiv \mathrm{Fix}|_{\langle \pi_y \pi_{xz}, \tau(L_y/4, -L_z/4) \rangle}, \\
S' &\equiv \mathrm{Fix}|_{\langle \pi_y \pi_{xz}, \tau(L_y/4, L_z/4) \rangle}.
\end{aligned} 
\quad (4.2) 
$$

For the flow restricted to subspace $S$ ($S'$), FP2 is a saddle (sink), FP2' is a sink (saddle) and there exists a saddle-sink connection FP2→FP2' (FP2'→FP2). More details of the conditions required for a robust heteroclinic cycle can be found in Krupa (1997), Reetz & Schneider (2020a) and § 4.2.3 of Zheng et al. (2024b). Like PO2, this robust cycle results from a 1:2 mode interaction (Armbruster, Guckenheimer & Holmes 1988). As the arguments on robustness and 1:2 interaction are fundamentally the same as for PO2, in order to avoid repetition, we refer readers to § 4.2.3 of our previous paper (Zheng et al. 2024b).

We have continued PO6 forward in Rayleigh number up to $Ra = 6618.12$ where it has period $T = 712.76$. Along the branch, PO6 undergoes a sequence of saddle-node bifurcations in the range $6450 \lesssim Ra \lesssim 6460$ and then at $Ra \approx 6505$, as shown in figure 8. Beyond $Ra = 6550$, the period of PO6 increases monotonically and seems to diverge. Integrating PO6 at $Ra = 6618.12$ forward in time, we observe that its trajectory tends to visit two symmetrically related instances of FP2, with the time spent near FP2 substantially shorter than at $Ra = 6218.6$. We suggest that PO6 may bifurcate again from FP2 in another heteroclinic bifurcation slightly beyond $Ra = 6618.12$. However, since the
<figure><img src="image_11.png" /></figure>

Figure 9. (a–e) Snapshots of the dynamics of PO6 at Ra = 6218.6. Snapshots (a) and (d) show states which are close to two symmetry-related versions of FP2. (f) Time series of PO6 at Ra = 6218.6 (with period T = 1069.1). (g) Phase space projection at Ra = 6218.6 close to the global bifurcation point. The curve shows PO6 and triangles show two symmetry-related FP2 states involved in the heteroclinic cycle. In (f) and (g), the five red stars indicate the moments at which the snapshots (a–e) are taken. In (g), the red arrows show the direction of the trajectory.

<figure><img src="image_12.png" /></figure>

Figure 10. Equilibria and eigenmodes at Ra = 6218.6. Case (a) FP2, (b) its unstable eigenmode $e_1$ and (c) its stable eigenmode $e_2$. Case (d) $FP2' \equiv \pi_y \tau(4, 0)FP2$, (e) its unstable eigenmode $e_1' \equiv \pi_y \tau(4, 0)e_1$ and (f) its stable eigenmode $e_2' \equiv \pi_y \tau(4, 0)e_2$. The wavenumbers of the equilibria and eigenmodes in the y-direction suggest a 1:2 mode interaction.

continuation becomes computationally difficult for higher Rayleigh numbers, we have not been able to confirm this. We do not show snapshots of PO6 at Ra = 6618.12, as they do not differ substantially from those in figure 9(a–e).

### 4.2.2. Orbit PPO7: saddle-node bifurcations and isola

Like RPO15, RPO17, RPO26 and RPO28, which are discussed in § 4.1.2, branch PPO7 also forms an isola, with two saddle-node bifurcations at Ra = 6280.1 and 6417.2, as
<figure><img src="image_13.png" /><figcaption>Figure 11. Temperature norms (a) and periods (b) of RPO5, PPO7, RPO8, PO9 and PPO30. Branches RPO5, RPO8 and PO9 undergo saddle-node bifurcations and are continued until $Ra = 6635$ for one of their endpoints. At the other endpoints, RPO5 possibly bifurcates from FP4 in a global bifurcation, and the termination of RPO8 and PO9 are unclear. Branch PPO7 undergoes saddle-node bifurcations and forms an isola. Branch PPO30 bifurcates from and terminates on PPO7 in two period-doubling bifurcations.</figcaption></figure>

shown in figure 11. Orbit PPO7 gives rise to PPO30 via two period-doubling bifurcations that will be discussed later in § 4.2.7.

The dynamics of PPO7 at $Ra = 6280.38$ (the period-doubling bifurcation point) is shown in figure 12(a)–(j). The simulation starts from a state close to the moustache rolls (FP5 in Zheng et al. 2024b); the roll at the domain centre (as well as corners due to $\tau(L_y/2, L_z/2)$ symmetry) then becomes more intense, distorted and ramified (figure 12(b–e)). At $t = 115$, the central roll pinches off and merges with neighbouring rolls at $t = 130$. After a smooth transition towards nearly straight rolls at $t = 180$, the trajectory returns to the distorted-roll state at $t = 267$.

Orbit PPO7 has the spatio-temporal symmetry

$$ (u, v, w, \theta)(x, y, z, t + T) = \pi_y(u, v, w, \theta)(x, y + L_y/4, z - L_z/4, t), \quad (4.3) $$

where $T = 226$ is the pre-period of PPO7 at $Ra = 6280.38$. After a pre-period, the state at $t = 267$, figure 12(i), is a reflected and translated version of the state at $t = 41$, figure 12(b); after integrating over a second pre-period, the states at $t = 493$ and $t = 41$ are related by $\sigma \equiv \tau(\pm L_y/2, 0)$ or $\tau(0, \pm L_z/2)$. Finally, after integrating during four pre-periods, the initial state matches the final state, i.e. PPO7 is a periodic orbit.

A remarkable feature of PPO7 in figure 3(a) is that its minimum temperature norm along the branch ($||\theta||_2 \approx 0.01$) is almost the lowest among all orbits. The faint figure 12(f)–(g) corresponds to the moments of lowest temperature norm. These are the moments of a cutting-joining-like dynamics of convection rolls, very similar to the longitudinal burst pattern observed by Daniels, Plapp & Bodenschatz (2000) experimentally and by Reetz et al. (2020b) numerically, at slightly different control parameters. The cutting-joining dynamics induces defects, disordered structures in rolls and roll bursting; these contribute to the transition to turbulence. To the best of our knowledge, PPO7 may be the first (pre-) periodic orbit that captures these aspects of the dynamics.

### 4.2.3. Orbit RPO10: saddle-node and global bifurcations

We first found RPO10 at $Ra = 6400$. Forward continuation in $Ra$ shows that it exists until at least $Ra = 6675$, where we stopped the continuation. Continuing RPO10 backwards in $Ra$, it undergoes a sequence of saddle-node bifurcations after which its period increases
<figure><img src="image_14.png" /><figcaption>Figure 12. Dynamics of PPO7 at Ra = 6280.38 with pre-period T = 226. (a–j) Snapshots of the midplane temperature field. (k) Time series from DNS. The ten red stars indicate the moments at which the snapshots (a)–(j) are taken.</figcaption></figure>

monotonically, as evidenced by figure 8. We have been able to continue RPO10 until Ra = 6298.686 with relative period T = 623.35, shortly after the last saddle-node bifurcation at (Ra, T) = (6298.39, 622.97).

Integrating RPO10 at Ra = 6298.686 in time, we observe a long plateau (250 < t < 460) in the time series shown in figure 13(a). The dynamics of RPO10 at this Ra is shown in figure 13(b)–(d). The states corresponding to the location of the plateau are very similar to FP4, in terms of both flow structure and temperature norm. This suggests that RPO10 disappears in a global homoclinic bifurcation by meeting FP4, although the slope of the last computed portion of T (Ra) in figure 8 is not as close to vertical as the corresponding slopes for PO6, PO14, PPO16 and RPO29.

### 4.2.4. Orbit RPO25: pitchfork, saddle-node and period-halving bifurcations

As shown in figure 7, RPO25 bifurcates from RPO19 at Ra = 6260.5 in a pitchfork bifurcation, in which the translation symmetry $τ(L_y/4, L_z/4)$ is broken to $τ(L_y/2, L_z/2)$ and the reflection symmetry $π_yπ_{xz}$ is retained. Orbit RPO25 then undergoes several saddle-node bifurcations and finally terminates in a period-halving bifurcation at Ra = 6432.34 (indicated in figure 7). We do not show or discuss the period-halved orbit in this work.
<figure><img src="image_15.png" /><figcaption>Figure 13. Dynamics of RPO10 at Ra = 6298.686 with relative period T = 623.35. (a) Time series of RPO10; the three red stars indicate the moments at which the snapshots (b)–(d) of the midplane temperature field are taken.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14. (a) Periods and (b) time series of RPO29. (Branch RPO29 also appears as part of figure 8.) The inset in (a) shows a sequence of saddle-node bifurcations before the global bifurcation at Ra ≈ 6274.14. (b) Time series from the last continuation point (longest period) at Ra = 6274.144 and Ra = 6402.012. Branch RPO29 approaches FP2 and FP4 in two different global homoclinic bifurcations at its two endpoints.</figcaption></figure>

4.2.5. **Orbit RPO27: pitchfork and saddle-node bifurcations**

As shown in figure 5, RPO27 bifurcates from RPO18 at Ra = 6279.7 in a pitchfork bifurcation at which τ(L_y/4, L_z/4) symmetry is broken to τ(L_y/2, L_z/2) (with π_y, π_z retained). It undergoes a sequence of saddle-node bifurcations, particularly between Ra = 6420 and 6500, and we have continued it until Ra = 6650.

4.2.6. *Orbit RPO29: saddle-node and global bifurcations*

We first observed RPO29 at Ra = 6300. Its bifurcation diagram is shown in figure 8; we
show its period again in figure 14(a). Backward continuation reveals the interesting feature
of at least 13 saddle-node bifurcations (and probably even more if it could be continued
further towards higher periods) before the global bifurcation at which the period diverges.
We have been able to continue RPO29 down to Ra = 6274.144, where it has relative period
T = 563.38. Orbit RPO29 undergoes successive approaches to FP4 (see time series in
figure 14b), as is the case for RPO10 in § 4.2.3.

Multiple regularly spaced saddle-node bifurcations, as in figure 14(a), are seen in at least
two situations. One is the homoclinic snaking of branches of spatially localised equilibria
(Batiste et al. 2006; Burke & Knobloch 2006) and periodic orbits (Al Saadi et al. 2024).
At each saddle-node bifurcation, a new pair of rolls or structures is added to the solution.
For an increase in periods, as in figure 14(a), we could conjecture an increasing number of
Z. Zheng, L.S. Tuckerman and T.M. Schneider

temporal maxima. However, in our case, we have checked that the states along the RPO29 branch do not show an increasing number of rolls or structures, nor an increasing number of peaks after each saddle-node bifurcation.

Another situation in which multiple regularly spaced saddle-node bifurcations of periodic orbits occur is the Shil'nikov bifurcation. This is the approach to a homoclinic orbit that is formed at the collision of a limit cycle with a saddle focus which has one real eigenvalue and one complex eigenpair (Glendinning & Sparrow 1984). In our case, RPO29 and FP4 are in the symmetry subspace $\langle \pi_y, \pi_{xz}, \tau(L_y/2, L_z/2) \rangle$ and hence we consider the leading eigenvalues of FP4 in this symmetry subspace. At $Ra = 6274.144$, these are $[\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_{5,6}] = [0.0384, 0.0364, 0.0096, 0.0019, -0.0172 \pm 0.1167i]$. We have determined that the eigenvalues controlling the escape from and approach to FP4 are those closest to zero, $\lambda_4$ and $\lambda_{5,6}$. Glendinning & Sparrow (1984) show that multiple saddle-node bifurcations occur all the way to the global bifurcation if the rate of escape exceeds the rate of approach, but in our case, the rate of approach exceeds the rate of escape by a factor of $| -0.0172/0.0019| \approx 9$. This implies that the oscillations in figure 14(a) will cease when the period becomes large enough, to be replaced by a monotonic increase in period. Unfortunately, we have not been able to continue RPO29 further towards higher periods and thus cannot confirm the potential monotonic increase.

Continuing RPO29 forward from $Ra = 6300$, many saddle-node bifurcations are also seen and the branch also terminates in a global homoclinic bifurcation, this time by meeting FP2. We have been able to continue RPO29 until $(Ra, T) = (6402.012, 477.93)$. Even though we believe that this is still far from the actual global bifurcation point as the period is not yet very long, a close approach to FP2 is evidenced by a clear plateau (whose corresponding norm is very close to that of FP2) in the time series in figure 14(b) and by inspection of flow fields (not shown).

### 4.2.7. Orbit PPO30: period-doubling bifurcations

Orbit PPO30 bifurcates from and terminates on PPO7 in two period-doubling bifurcations, at $Ra = 6280.38$ and 6417.25, indicated in figure 11. For PPO30

$$ \begin{aligned} (u, v, w, \theta)(x, y, z, t+T) &= (u, v, w, \theta)(x, y \pm L_y/2, z, t), \\ &= (u, v, w, \theta)(x, y, z \pm L_z/2, t), \end{aligned} \quad (4.4) $$

where $T$ is the pre-period of PPO7. After two pre-periods, the initial state matches the final one. The quantities $||\theta||_2$ of PPO7 and PPO30 are almost indistinguishable as can be seen in figure 11 and as is usual for period-doubling bifurcations.

### 4.3. Symmetry subspace: fourfold translation

In this subsection, we discuss RPO12 and RPO20, two time-periodic solutions identified in the symmetry subspace $\langle \tau(L_y/4, L_z/4) \rangle \sim Z_4$. These branches are shown in the bifurcation diagram of figure 15.

#### 4.3.1. Orbit RPO12: saddle-node and global bifurcations

Orbit RPO12 undergoes saddle-node bifurcations at $Ra = 6293.9$, in the range $6643 \lesssim Ra \lesssim 6646$ and at $Ra = 6675.8$. The lower branch (in terms of period) of RPO12 has been continued until $(Ra, T) = (6680, 111.4)$ where we stopped the continuation. The upper branch is continued until $Ra = 6654.865$ where the relative period $T = 301.9$ is the highest that we were able to attain numerically.

Integrating RPO12 at $Ra = 6654.865$, we observe that the dynamics slows down slightly near a state that is close to FP4; see figure 16(a) and its inset. We subsequently used
<figure><img src="image_17.png" /><figcaption>Figure 15. Temperature norms (a) and period (b) of RPO11, RPO12 and RPO20. Branch RPO11 undergoes saddle-node bifurcations and both the lower and upper branches are continued beyond Ra = 6680; its bifurcation structure remains unclear. The lower RPO12 branch exists beyond Ra = 6680, while the upper branch seems to terminate in a global bifurcation by meeting FP4, close to Ra = 6655. Branch RPO20 bifurcates from FP4 in a global bifurcation at Ra ≈ 6561 at which its period seems to diverge; its termination is unclear.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 16. (a) Time series of RPO12 with relative period $T = 301.9$ at Ra = 6654.865. A snapshot of the midplane temperature field at instant $t = 170$ is shown in the inset and is close to FP4. (b) Time series of RPO20 with relative period $T = 343$ at Ra = 6561.2. The five red stars indicate the moments at which the snapshots (c)–(g) are taken. Snapshot (e) is close to FP4.</figcaption></figure>

this state ($t = 170$) as the initial guess for Newton's method to converge to FP4 at Ra = 6654.865. Even though the plateau in figure 16(a) is not obvious and the period of RPO12 in figure 15 is not yet very long, our observations suggest that RPO12 terminates in a global homoclinic bifurcation by meeting FP4. This scenario is very similar to that of RPO20 to be discussed just after in § 4.3.2, for which we will show more snapshots.
### 4.3.2. Orbit RPO20: saddle-node and global bifurcations

We have continued RPO20 to $Ra = 6561.2$ with relative period $T = 343.6$ and to $Ra = 6660.2$ with $T = 245.5$. Close to $Ra = 6561.2$, its period seems to diverge, slightly more so than that of RPO12. Figure 16(b) shows the time series from a simulation of RPO20 at $Ra = 6561.2$, with the corresponding snapshots of the temperature field shown in figure 16(c)–(g). The plateau-like behaviour ($150 \lesssim t \lesssim 250$) in the time series as well as the close resemblance between figure 16(e) and FP4 suggest that RPO20 bifurcates from FP4 in a global homoclinic bifurcation at a nearby Rayleigh number.

By looking at the snapshots in figure 16(c)–(g), we notice that the dynamics is not very different from previous cases, as the reflection symmetry $\langle \pi_y \pi_{xz} \rangle$ is only weakly broken by the global bifurcation from FP4. Since it is clear that the RPO20 (and RPO12 in § 4.3.1) with the longest period that we have succeeded in computing is still far from the actual homoclinic cycle, we do not discuss or show the eigendirections of FP4 along which RPO20 (and RPO12) may approach and escape from FP4. The termination of the other end of the RPO20 branch is unclear and not discussed.

### 4.4. Symmetry subspace: twofold translation

In this subsection, we discuss two relative periodic orbits (RPO5 and RPO8) identified in the symmetry subspace $\langle \tau(L_y/2, L_z/2) \rangle \sim Z_2$. Their branches are contained in the bifurcation diagram of figure 11, while their time series and snapshots are shown in figure 17.

#### 4.4.1. Orbit RPO5: saddle-node and global bifurcations

In the $Ra$ range we study, RPO5 undergoes a sequence of saddle-node bifurcations. The lower branch (in period) continues to exist until at least $Ra = 6635$. For the upper branch, the seemingly diverging period at $Ra \approx 6510.4$ suggests that RPO5 might disappear in a global bifurcation. We have been able to continue RPO5 until $Ra = 6510.4$ with relative period $T = 400.5$. Integrating RPO5 at $Ra = 6510.4$ in time, the dynamics slows down slightly close to FP4 ($t \approx 300$), as shown by the time series in figure 17(f) and the snapshot in figure 17(d). We expect that the time spent near FP4 would increase if we were able to continue RPO5 further.

#### 4.4.2. Orbit RPO8: saddle-node bifurcations

As shown in figure 11, RPO8 undergoes a sequence of saddle-node bifurcations and is continued until $Ra = 6636.26$ (lower branch, where relative period $T = 246.97$) and $Ra = 6416.88$ (upper branch, where $T = 373.92$). With the available information, we have not been able to determine the origin of RPO8. Figure 17(h)–(l) shows five snapshots of RPO8 at $Ra = 6388.46$. Like PPO7, rolls in RPO8 tend to distort and to develop defects, and the variation of $||\theta||_2$ along the orbit is large; compare for instance figures 17(h) and 17(j). Note also that figure 17(k) is similar to a less symmetric version of FP8 shown in figure 2(f).

### 4.5. Symmetry subspace: reflection with threefold translation

Only one orbit is identified in the symmetry subspace $\langle \pi_y \pi_{xz}, \tau(L_y/3, L_z/3) \rangle \sim D_3$.

#### 4.5.1. Orbit PO14: Hopf and global bifurcations

As shown in figure 8, PO14 bifurcates from FP11 at $Ra = 6289.6$ in a symmetry-preserving Hopf bifurcation. (For consistency with symmetry groups of other orbits, we do not introduce tilted coordinates for PO14 as we did for FP9–FP12 in § 3.2. It can be verified that FP9, FP11 and FP12 all have the symmetry $\langle \pi_y \pi_{xz}, \tau(L_y/3, L_z/3) \rangle$ in the $y-z$ coordinate.) The period of PO14 close to the bifurcation point closely matches
<figure><img src="image_19.png" /><figcaption>Figure 17. Dynamics of RPO5 with relative period $T = 400.5$ at $Ra = 6510.4$ and of RPO8 with $T = 375.4$ at $Ra = 6388.46$. (a–e, h–l) Snapshots of the midplane temperature field. Snapshot (d) of RPO5 is similar to FP4 (figure 2d) and converges to FP4 when used as an initial guess for Newton’s method. (f–g) Time series, initialised by the states shown in (a) and (h).</figcaption></figure>

the value predicted by the bifurcating imaginary eigenvalue pair of FP11. Equilibria FP9, FP10 and FP11 are called secondary, tertiary and quaternary states, respectively. Orbit PO14 can thus be called a quinary state. Forward continuation of PO14 suggests that it terminates in a global homoclinic bifurcation by meeting FP9, close to $Ra = 6313$ at which its period diverges.

Figure 18(a)–(d) shows four snapshots of PO14 at $Ra = 6313$, the highest Rayleigh number that we have reached. The corresponding time series and phase space projection are shown in figure 18(e)–(f), with special instants indicated by red stars. The long plateau (250 $\lesssim t \lesssim 550$) in the time series and the clustering of points close to FP9 in the phase space projection suggest that PO14 approaches and slows down near FP9. Interestingly, the time series also shows another short plateau near $t \approx 690$, whose corresponding state, shown in figure 18(d), resembles FP12 shown in figure 2(j). However, the non-negligible difference of norms in figure 18(e)–(f) between FP12 and the state shown in figure 18(d) suggests that PO14 does not visit FP12 closely. It might be possible that, if we were able to continue PO14 further with longer periods, FP12 would be visited more closely by PO14. But based on the available data, we conclude that PO14 ends by meeting FP9 in a homoclinic cycle.
<figure><img src="image_20.png" /><figcaption>Figure 18. Dynamics of PO14 with period $T = 775.62$ at $Ra = 6313$ (close to the global bifurcation point). (a–d) Snapshots of the midplane temperature field. Snapshots (c) and (d) converge to FP9 and FP12 when used as initial guesses. (e) Time series from DNS. (f) Phase space projection: shown are PO14 (curve with dots) as well as FP9 and FP12 (triangles). In (e) and (f), the four red stars indicate the moments at which the snapshots (a)–(d) are taken. (g) The $L_2$-distance between each instantaneous flow field of PO14 and FP9 (and FP12). The dynamics of PO14 is exponential for most of the cycle (blue curve). The approaching (black dashed line) and escaping (red dashed line) dynamics of PO14 with respect to FP9 are shown and are governed by two eigenvalues, $\lambda_1$ and $\lambda_2$, of FP9. (h–i) Two eigenmodes $e_1$ and $e_2$ of FP9, visualised via the midplane temperature field.</figcaption></figure>

In order to analyse this homoclinic cycle, we have computed the spectrum of FP9 at $Ra = 6313$. (Rather than referring back to figure 2(g), the reader can look at figure 18(c), which closely resembles FP9. For a detailed explanation of the correspondence between global bifurcations and the eigenvalues of the equilibria that are approached by the trajectories, see Zheng et al. 2024b.) Restricting the computation to the symmetry subspace $(\pi_y, \pi_{xz}, \tau(L_y/3, L_z/3))$ gives three leading eigenvalues, all real: $[\lambda_1, \lambda_2, \lambda_3] = [0.0162, -0.0077, -0.012]$. Since we imposed FP9’s symmetries in the Arnoldi iterations,
<figure><img src="image_21.png" /><figcaption>Figure 19. Dynamics of PPO16 with pre-period $T = 1007.05$ at $Ra = 6240.6429$ (close to the global bifurcation point). (a–e) Snapshots of the midplane temperature field. Snapshot (b) converges to FP9 when used as an initial guess. (f) Time series from DNS. The five red stars indicate the moments at which the snapshots (a–e) are taken. (g) The $L_2$-distance between each instantaneous flow field of PPO16 and FP9. The dynamics of PPO16 is exponential for most of the cycle (blue curve). The approaching (black dashed line) and escaping (red dashed line) dynamics of PPO16 with respect to FP9 are shown to be governed by two eigenvalues, $\lambda_1$ and $\lambda_2$, of FP9. (h–i) Two eigenmodes $e_1$ and $e_2$ of FP9, visualised via the midplane temperature field.</figcaption></figure>

the two neutral eigenvalues corresponding to translations in the periodic directions $y$ and $z$ are not present. We have determined that the eigendirections along which PO14 leaves and approaches FP9 are $e_1$ and $e_2$ associated with $\lambda_1$ and $\lambda_2$, respectively. These eigendirections are confirmed by subtracting FP9 from the instantaneous flow fields of PO14, and comparing the resulting fields with eigenmodes obtained by Arnoldi iterations, as well as by the close matches for the exponential growth and decay rate between PO14 and FP9, see figure 18(g). Since $\lambda_1 > |\lambda_2|$, the homoclinic orbit is unstable, which is confirmed by the chaotic behaviour in the time series after time integrating about two periods (not shown in figure 18e).

The two relevant eigendirections can be interpreted and analysed by comparing PO14 and FP9. Eigenmode $e_1$, shown in figure 18(h), breaks the $O(2)$ symmetry of FP9 along its straight and homogeneous rolls by introducing alternating red and blue patches in this tilted direction; these patches lead to wavy-roll structures. It is not difficult to imagine that superposing FP9 (figure 18(c)) and $e_1$ gives approximately figure 18(d). Eigenmode $e_2$, shown in figure 18(i), consists of a grid of red and blue rhombi, while figure 18(b) consists of rolls with bulges and constrictions. The colours of the rhombi are opposite to those of the bulges and the same as those of the constrictions. Thus, superposing $e_2$ on figure 18(b) reduces both distortions, restoring the broken symmetries of FP9.

## 4.6. Symmetry subspace: reflection with fivefold translation

Three orbits identified in the symmetry subspace $\langle \pi_y, \pi_{xz}, \tau(L_y/5, L_z/5) \rangle \sim D_5$ are discussed in this subsection. The dynamics of these three orbits appears to be similar and we only show snapshots of PPO16 (figure 19(a–e)) for illustration.
<figure><img src="image_22.png" /></figure>

Figure 20. Temperature norms (a) and periods (b) of RPO21, RPO22, PO23 and PO24. On the left, the minima of $||\theta||_2$ of RPO21 and RPO22 are too close to be distinguished; the lack of smoothness in the maxima of $||\theta||_2$ of RPO21 corresponds to the overtaking of one temporal maximum or minimum of $||\theta||_2$ by another as Ra is varied. The creation and termination of RPO21 and RPO22 are not discussed. Both PO23 and PO24 bifurcate from FP8 in two Hopf bifurcations; PO23 possibly terminates in a global bifurcation at $Ra \approx 6589.5$ by meeting FP13, and PO24 exists until at least $Ra = 6667$.

### 4.6.1. Orbit PPO16: global bifurcation

Orbit PPO16 is a pre-periodic orbit; its spatial phase shifts by $\langle \tau(L_y/10, L_z/10) \rangle$ after a pre-period; compare figures 19(a) and 19(e). After ten such pre-periods, the final state matches the initial state. The branch of PPO16 states is included in the bifurcation diagram of figure 8. We have continued PPO16 towards increasing $Ra$ to $Ra = 6656.54$. Towards decreasing $Ra$, the period of PPO16 increases monotonically and eventually diverges, suggesting a global bifurcation. Figure 19(a–e) shows five snapshots of PPO16 at $Ra = 6240.6429$ with pre-period $T = 1007.05$, the lowest Rayleigh number we have reached. The corresponding time series in figure 19(f) indicates that PPO16 slows down significantly between $150 \lesssim t \lesssim 800$ and spends a long time near an oblique-roll state (figure 19b). This oblique-roll state is subsequently converged via Newton’s method to FP9; figures 19(b) and 2(g) are related by $\pi_y$. Similarly to PO14 described in § 4.5.1, PPO16 also bifurcates from FP9 in a global homoclinic bifurcation.

We computed the eigenvectors and eigenvalues of FP9 at $Ra = 6240.6429$ in the symmetry subspace $\langle \pi_y\pi_{xz}, \tau(L_y/5, L_z/5) \rangle$. The Arnoldi iterations return five leading eigenvalues, all real: $[\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5] = [0.01651, -0.00631, -0.057, -0.0628, -0.07664]$. Clearly, PPO16 escapes from FP9 along $e_1$, associated with $\lambda_1$ and shown in figure 19(h), the only unstable eigendirection in this subspace. The direction along which PPO16 approaches FP9 is $e_2$, associated with the second eigenvalue $\lambda_2$ and shown in figure 19(i). These two eigendirections are subsequently confirmed by subtracting FP9 from PPO16, as well as by the exponential decay and growth rates shown in figure 19(g). Similarly to the scenario for PO14, here $e_1$ breaks the O(2) symmetry of FP9 in the tilted direction and $e_2$ restores them.

### 4.6.2. Orbits RPO21 and RPO22: saddle-node bifurcations

Two relative periodic orbits, RPO21 and RPO22, are shown in the bifurcation diagram of figure 20; both undergo a sequence of saddle-node bifurcations, and the termination and/or creation of both orbits remain unclear. Branch RPO21 exists over the short range $6291.23 < Ra < 6448.75$; its two endpoints based on our continuation are quite close together, at $(Ra, T) = (6352.99, 402.83)$ and $(Ra, T) = (6354.95, 376.35)$. Integrating
<figure><img src="image_23.png" /><figcaption>Figure 21. Dynamics of PO23 at Ra = 6589.47 with period T = 404.6. (a–e) Snapshots of the midplane temperature field. (f) Time series from DNS. The five red stars indicate the moments at which the snapshots (a–e) are taken.</figcaption></figure>

RPO21 in time at these two values of Ra does not show remarkable behaviour or a close approach to an equilibrium that would indicate a global bifurcation. Branch RPO22 originates in a saddle-node bifurcation at Ra = 6259; the upper and lower branches which emanate from it can be continued at least until Ra = 6667.

4.7. Symmetry subspace: reflection

In this subsection, we discuss two periodic orbits, PO23 and PO24, that are in the symmetry subspace $\langle \pi_y, \pi_{xz} \rangle \sim Z_2$. As shown in figure 20, both bifurcate from FP8 at Ra = 6367.9 and Ra = 6321, in two Hopf bifurcations which preserve the reflection symmetry of FP8 and break its translation symmetry $\langle \tau(L_y/2, 0) \rangle$.

4.7.1. *Orbit PO23: Hopf, saddle-node and global bifurcations*

After its creation, PO23 undergoes a sequence of saddle-node bifurcations. Figure 21(a–e) shows five snapshots of PO23 at Ra = 6589.47. The initial phase (25 ≲ t ≲ 175) of PO23 resembles FP13, compare figures 21(a) and 2(k); we used the state in figure 21(a) as an initial guess for Newton’s method to converge to FP13. We have been able to continue PO23 until Ra = 6589.47, where its period is T = 404.6. Figure 20 shows that its period seems to diverge; we believe that PO23 terminates on FP13 in a global bifurcation point at a nearby Ra.

4.7.2. *Orbit PO24: Hopf bifurcation*

After bifurcating from FP8 at Ra = 6321, PO24 is continued until Ra = 6667 where we stopped the continuation. It is clear that PO24 oscillates around FP8 and the oscillation amplitude is smaller than that of PO23. We do not show snapshots of PO24.

4.8. No spatial symmetry

In this subsection, we discuss two orbits without any spatial symmetries.
<figure><img src="image_24.png" /><figcaption>Figure 22. Dynamics of PO9 with period $T = 311.18$ at $Ra = 6413.11$. (a) Time series from DNS. The four red stars indicate the moments at which the snapshots (b)–(e) are taken.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 23. Dynamics of RPO11 with relative period $T = 209.26$ at $Ra = 6500$. (a) Time series from DNS. The four red stars indicate the moments at which the snapshots (b)–(e) are taken.</figcaption></figure>

### 4.8.1. Orbit PO9: saddle-node bifurcations

Orbit PO9 is shown in the bifurcation diagram of figure 11. Orbit PO9 undergoes a sequence of saddle-node bifurcations. We have continued the lower branch (in period) of PO9 until ($Ra, T$) = (6631.4, 188.7) and the upper branch until ($Ra, T$) = (6475.16, 314.47); the origin of PO9 remains unclear. Snapshots and time series of PO9 at $Ra = 6413.11$ are shown in figure 22. Comparing snapshots (b) and (e), we notice that PO9 has the spatio-temporal symmetry

$$ (u, v, w, \theta)(x, y, z, t + T/2) = \pi_y \pi_{xz}(u, v, w, \theta)(x, y + 0.11 L_y, z + 0.03 L_z, t), \quad (4.5) $$

where $T$ is the period of PO9. Since (4.5) contains the combined reflection operator $\pi_y \pi_{xz}$, after two such pre-periods, the discrete translations in $L_y$ and $L_z$ cancel out and the final state is identical to the initial state – an actual periodic orbit which we converged in our study.

### 4.8.2. Orbit RPO11: saddle-node bifurcations

Branch RPO11 is contained in the bifurcation diagram of figure 15. Branch RPO11 undergoes a sequence of saddle-node bifurcations. We stopped the continuation at ($Ra, T$) = (6680, 212) for the lower branch and at ($Ra, T$) = (6680, 230.6) for the upper branch; its origin remains unknown. Figure 23(b)–(e) shows four snapshots of RPO11 at ($Ra, T$) = (6500, 209.26). Like several of the other periodic orbits, RPO11 contains a state (figure 23c) with four relatively straight but deformed rolls, which breaks along a diagonal fault line (figure 23d), leading to disordered states (figure 23e,b), which then reform into approximate rolls (figure 23c).

## 5. Dynamical relevance of unstable periodic orbits and statistical descriptions

In order to illustrate the possible dynamical relevance of the periodic orbits that we have computed and discussed in § 4, we show in this section some phase space projections
<figure><img src="image_26.png" /><figcaption>Figure 24. Phase space projection at Ra = 6300. The plot shows the projection onto the thermal energy input (I) and the viscous dissipation over energy input (D/I) of 34 periodic orbits and of instantaneous flow fields, separated by $\Delta t = 1$, of the chaotic dynamics during a DNS of length $2 \times 10^5$ time units. The inset shows the chaotic dynamics only (the dots appear slightly denser due to the inset's smaller size). The subscripts n in POXn indicate different orbits on the same solution branch related by saddle-node bifurcations.</figcaption></figure>

and statistical evidence supporting the close visit of orbits by the chaotic dynamics. Here, we focus on the dynamics and statistics at a fixed Rayleigh number Ra = 6300 which is approximately 10% above the onset of convection. At Ra = 6300, we identified 19 different periodic-orbit branches. Because of the many saddle-node bifurcations causing branches to zigzag back and forth in Ra, several different solutions can be located at the same value of Ra on the same branch. In particular, 34 periodic orbits exist at Ra = 6300. We use the notation POXn to refer to the nth solution on branch POX.

## 5.1. Phase space projections

The projections of a long DNS ($2 \times 10^5$ time units, initiated with random small-amplitude noise) and of 34 orbits are shown in figure 24. This projection employs two global quantities, the thermal energy input (I) and viscous dissipation (D), as observables; instantaneous flow fields are represented by dots (for DNS) and closed loops (for orbits) in the $D/I - I$ plane. It can be seen that most of the orbits lie on the core part ($0.05 \lesssim I \lesssim 0.1$) of the attractor in the current projection. In addition, three orbits shown with thicker curves (PPO7$_{1,2}$ and PPO30) are able to capture some of the very low input events ($I \approx 0.01$) of the DNS trajectory. (Recall in § 4.2.7 that PPO30 bifurcates from the PPO7 branch in period-doubling bifurcations.) To verify that the close match is not merely an artefact of the projection, we examine the flow fields corresponding to these instants. The comparison is shown in figure 25. In figure 25(a), we select a short portion of the DNS trajectory (550 time units, belonging to the long DNS trajectory shown in figure 24) closely following PPO7$_{1,2}$ in the projection, and we show three temperature fields along the trajectory (marked by crosses in figure 25(a) in figure 25(b)–(d). The corresponding flow fields of PPO7$_2$ (which are closest to those of the DNS trajectory both in terms of projection and $L_2$-norm) are shown in figure 25(e)–(g). Comparing figures 25(b,e),
<figure><img src="image_27.png" /><figcaption>Figure 25. Chaotic dynamics and PPO7 at Ra = 6300. (a) Projection as in figure 24 but with only a short portion of DNS and two orbits. The inset zooms in on the slow dynamics close to D = I and I ≈ 0.063. The crosses indicate instants at which the snapshots (b–g) are taken and the triangles indicate the beginning and end of the selected DNS trajectory. (b–d) Temperature fields corresponding to three instants of the chaotic dynamics shadowing PPO7₂. (e–g) Temperature fields corresponding to three instants of PPO7₂.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 26. Chaotic dynamics and RPO12 at Ra = 6300. (a) Projection as in figure 24 but with only a short portion of DNS and two orbits. The crosses indicate instants at which the snapshots (b)–(g) are taken and the triangles indicate the beginning and end of the DNS trajectory. (b–d) Temperature fields corresponding to three instants of the chaotic dynamics shadowing RPO12₁. (e–g) Temperature fields corresponding to three instants of RPO12₁.</figcaption></figure>

(c,f) and (d,g) convincingly suggests that PPO7₂ is very closely visited by the chaotic dynamics. (Note that we have shifted the flow fields in y and/or z to facilitate the visual comparison; the optimal shifts are determined by minimising the L₂ difference between two flow fields.)

The same analysis can be carried out with almost all other orbits and we will illustrate another example on RPO12 which lies in the core part of the attractor. Figure 26(a) shows
<figure><img src="image_29.png" /><figcaption>Figure 27. Probability density functions (PDFs) of I and ||θ||₂ at Ra = 6300. Shown are the data from DNS and predicted statistics based on 34 periodic orbits. A total of 80 bins is used for each PDF.</figcaption></figure>

a short DNS trajectory (830 time units) approaching RPO12₁ and then RPO12₂. The almost indistinguishable flow fields between figure 26(b–d) and figure 26(e–g) suggest that RPO12 plays a basic role in the spatio-temporal dynamics of the system. Four orbits (RPO21₁,₂ and RPO22₁,₂) reaching very high input range (I ≈ 0.12) are not visited at all by the flow, based on the projection in figure 24. The reason is that, as stated in § 4.6, these four orbits are identified in the symmetry subspace ⟨π_y,π_xz, τ(L_y/5, L_z/5)⟩. The fivefold translation symmetry along the domain diagonal greatly constrains the possible dynamics and these four orbits are very unstable in the full phase space.

## 5.2. Reconstruction of flow statistics from periodic orbits

In § 5.1, we suggest that many unstable periodic orbits that we have identified are embedded in the trajectories followed by the flow, are closely visited by the chaotic dynamics, and cover the core part of the chaotic attractor. A logical next step is to quantify and to understand how the statistical properties of orbits are related to the statistics of the flow. For this reason, we plot in figure 27 two PDFs in terms of the quantities I and ||θ||₂, reconstructed from 34 orbits at Ra = 6300, together with the PDF of a long DNS (2 × 10⁵ time units). To reconstruct a PDF from periodic orbits, we use the formula

$$ \Gamma_{prediction}^{N} := \frac{\sum_{i=1}^{N} w_{i} \Gamma_{i}}{\sum_{i=1}^{N} w_{i}}, \qquad (5.1) $$

where Γ is the reconstructed PDF, N = 34 the number of orbits, Γᵢ the PDF of the i'th orbit and wᵢ the weight of the i'th orbit. Here, we consider all orbits to be equally important in contributing to the dynamics; in other words, each orbit has the same weight wᵢ := 1. For other heuristic choices of weights based on the stability and period of orbits, see § 5 of Chandler & Kerswell (2013).

Based on figure 27, it can be seen that the peaks of the PDF from DNS, or equivalently the core part of the attractor, are correctly captured by the prediction from periodic orbits. This is consistent with the previous observation on figure 24. For both large I and ||θ||₂ values (rightmost part of each PDF), the predictions show again that some orbits are not visited by the flow, as discussed in § 5.1. The PDFs in this region are dominated by four (very unstable) orbits (RPO21₁,₂ and RPO22₁,₂), and are 1–2 order(s) of magnitude smaller than that in the core regions. For 0.005 ≲ I ≲ 0.04 and 0.005 ≲ ||θ||₂ ≲ 0.025, we
<figure><img src="image_30.png" /><figcaption>Figure 28. (a) one snapshot of temperature field from DNS in a large spatial domain $[L_x, L_y, L_z] = [1, 80, 90]$ at $Ra = 6300$ and at a fixed time. Smaller boxes of size $[L_x, L_y, L_z] = [1, 8, 9]$ surround patterns that are (approximately) captured by invariant solutions, shown in the eight snapshots in (b), that have been studied in this work or in Zheng *et al.* (2024b). The grid used for the large domain computation has the same density of points as that for the smaller domain.</figcaption></figure>

observe non-negligible fluctuations in PDFs reconstructed from periodic orbits, while that of the DNS is relatively smooth. Looking back again at the projection in figure 24, there are indeed very few identified orbits in the region $0.005 \lesssim I \lesssim 0.04$, suggesting that we are probably missing some important orbits covering these parts of the attractor.

Even though the predictions of the statistics (with equal weights) already show reasonably good results, we will leave a deeper study on higher-order statistics, different weighting protocols and, eventually, periodic orbit theory (Cvitanović & Eckhardt 1991) to a future occasion. But as the results in this section suggest, to tackle a quantitative description of transitional turbulence via periodic orbits, a concerted effort is still required to identify a sufficient number of dynamically relevant periodic orbits embedded in the chaotic attractor.

### 5.3. Relevance to large domain chaos

To conclude this section, we show in figure 28(a) an instantaneous temperature field from a chaotic simulation in a large spatial domain of size $[L_x, L_y, L_z] = [1, 80, 90]$ at
$Ra = 6300$. Despite the much greater freedom allowed in the large domain, the traces of invariant solutions that we studied in the small domain of size $[L_x, L_y, L_z] = [1, 8, 9]$ are still present. The small numbered boxes contain portions of the pattern which resemble the equilibria and periodic orbits studied in this work and in Zheng *et al.* (2024b). The eight snapshots in figure 28(b) are labelled by the number of each box and by the corresponding solution names. In particular, the different phases of the roll-bursting dynamics (shown in figure 12 for PPO7) are seen. (Here, we determined these locations by eye, but a more systematic and quantitative approach could also be used.) This correspondence, while qualitative, suggests that unstable solutions in small domains are relevant to the spatio-temporal dynamics in a large domain and provide a promising framework for understanding the bifurcation-theoretic origins of various complex states.

## 6. Discussion and conclusions

In this work, we have discussed 26 newly identified periodic-orbit branches in the vertical thermal convection system at fixed Prandtl number $Pr = 0.71$ in a fixed-size domain $[L_x, L_y, L_z] = [1, 8, 9]$ and in the Rayleigh-number range $6225 \lesssim Ra \lesssim 6650$. These new branches, together with four previously studied in Zheng *et al.* (2024b), bring the total to 30 periodic-orbit branches. To the best of our knowledge, this is the largest number of solutions found thus far in three-dimensional Navier-Stokes systems and there certainly exist many more solution branches that we have not followed.

The bifurcations that these branches undergo and their spatial symmetries are summarised in table 1; eight different symmetry groups have been identified. Several orbits (RPO8, PO9, RPO11, RPO21 and RPO22) are of unknown bifurcation-theoretic origins in the Rayleigh-number range we have studied. Five isolas are found, which might be connected to other branches if other parameters were varied (e.g. inclination angle, Prandtl number or domain size); we do not explore this.

As mentioned throughout § 4 and in table 1, almost all of the branches undergo multiple saddle-node bifurcations. Because of this, one branch often contains several portions in the same Rayleigh-number range. Taking this into account, we find that 34 periodic orbits exist at $Ra = 6300$, 45 at $Ra = 6400$, and 23 at $Ra = 6500$, all with different periods, as well as a different dynamics and thus different statistics. In addition, there also exist ‘ghost states’ that emerge from saddle-node bifurcations (Zheng *et al.* 2025), which resemble the solutions at the saddle-node bifurcation. Like nearby unstable states, ghosts influence the trajectory of the chaotic dynamics near a saddle-node bifurcation and are relevant for the spatio-temporal patterns observed in weakly turbulent flows.

Local (Hopf) and global (heteroclinic or homoclinic) bifurcations are two standard scenarios by which periodic orbits are created or destroyed. We have found three orbits which are born from equilibria in Hopf bifurcations, and eight orbits that disappear or appear via global bifurcations. The orbit period diverges at global bifurcation points, which makes them challenging to compute, but we have been able to continue heteroclinic and homoclinic orbits, as we did in Zheng *et al.* (2024b).

In addition to various (complicated) bifurcation scenarios of periodic orbits, in § 5.1 we have presented evidence via phase space projections that suggest close approaches of the chaotic trajectory to certain periodic orbits, highlighting the dynamical relevance of unstable orbits to a chaotic flow. These time-periodic solutions can be said to locally guide the flow trajectory and their statistics. Using the PDFs of these orbits to reconstruct that of the flow, we find that the core part of the chaotic attractor is well captured by the prediction and that the overall agreement is satisfactory. Finally, we have located versions of our solutions in chaotic simulations in a much larger box. Together, the results in this work
Z. Zheng, L.S. Tuckerman and T.M. Schneider

emphasise the power of the nonlinear dynamical systems approach for shedding light on the convection patterns observed in high-dimensional spatio-temporally chaotic systems, and for quantitatively describing them.

**Acknowledgements.** We thank O. Ashtari for helpful discussions. We acknowledge E. Knobloch, D. Barkley and C. Beaume for their valuable comments on the paper. We are grateful to the Scientific IT and Application Support (SCITAS) team of EPFL for providing computational resources and high performance computing (HPC) expertise. We thank the referees for contributing their expertise.

**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.

## REFERENCES

AHLERS, G., GROSSMANN, S. & LOHSE, D. 2009 Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. *Rev. Mod. Phys.* **81**, 503–537.

AL SAADI, F., KNOBLOCH, E., NELSON, M., H., UECKER 2024 Time-dependent localized patterns in a predator-prey model. *Chaos Interdisciplinary J. Nonlinear Sci.* **34** (4), 043143.

ARICI, M., KARABAY, H. & KAN, M. 2015 Flow and heat transfer in double, triple and quadruple pane windows. *Energ. Build.* **86**, 394–402.

ARMBRUSTER, D., GUCKENHEIMER, J. & HOLMES, P. 1988 Heteroclinic cycles and modulated travelling waves in systems with O(2) symmetry. *Phys. D* **29** (3), 257–282.

BATISTE, O., KNOBLOCH, E., ALONSO, A. & MERCADER, I. 2006 Spatially localized binary-fluid convection. *J. Fluid Mech.* **560**, 149–158.

BODENSCHATZ, E., PESCH, W. & AHLERS, G. 2000 Recent developments in Rayleigh-Bénard convection. *Annu. Rev. Fluid Mech.* **32** (1), 709–778.

BOROŃSKA, K. & TUCKERMAN, L.S. 2010a Extreme multiplicity in cylindrical Rayleigh-Bénard convection. I. Time dependence and oscillations. *Phys. Rev. E* **81**, 036320.

BOROŃSKA, K. & TUCKERMAN, L.S. 2010b Extreme multiplicity in cylindrical Rayleigh-Bénard convection. II. Bifurcation diagram and symmetry classification. *Phys. Rev. E* **81**, 036321.

BUDANUR, N.B. & CVITANOVIĆ, P. 2017 Unstable manifolds of relative periodic orbits in the symmetry-reduced state space of the Kuramoto-Sivashinsky system. *J. Stat. Phys.* **167** (3–4), 636–655.

BURKE, J. & KNOBLOCH, E. 2006 Localized states in the generalized Swift-Hohenberg equation. *Phys. Rev. E* **73**, 056211.

CHANDLER, G.J. & KERSWELL, R.R. 2013 Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow. *J. Fluid Mech.* **722**, 554–595.

CVITANOVIĆ, P. 1991 Periodic orbits as the skeleton of classical and quantum chaos. *Phys. D* **51**, 138–151.

CVITANOVIĆ, P. 2013 Recurrent flows: the clockwork behind turbulence. *J. Fluid Mech.* **726**, 1–4.

CVITANOVIC, P., ARTUSO, R., MAINIERI, R., TANNER, G., VATTAY, G., WHELAN, N. & WIRZBA, A. 2005 *Chaos: Classical and Quantum. ChaosBook. Org.* Niels Bohr Institute. vol. **69**, p. 25.

CVITANOVIĆ, P. & ECKHARDT, B. 1991 Periodic orbit expansions for classical smooth flows. *J. Phys. A: Math. Gen.* **24**, L237–L241.

DANIELS, K., PLAPP, B. & BODENSCHATZ, E. 2000 Pattern formation in inclined layer convection. *Phys. Rev. Lett.* **84**, 5320–5323.

ENGEL, M., ASHTARI, O., SCHNEIDER, T.M. & LINKMANN, M. 2025 Search for unstable relative periodic orbits in channel flow using symmetry-reduced dynamic mode decomposition. *J. Fluid Mech.* **1013**, A45.

GAO, Z., PODVIN, B., SERGENT, A., XIN, S. & CHERGUI, J. 2018 Three-dimensional instabilities of natural convection between two differentially heated vertical plates: linear and nonlinear complementary approaches. *Phys. Rev. E* **97**, 053107.

GIBSON, J.F. 2019 Channelflow 2.0. Available at https://www.channelflow.ch.

GLENDINNING, P. & SPARROW, C. 1984 Local and global behavior near homoclinic orbits. *J. Stat. Phys.* **35**, 645–696.

GRAHAM, M.D. & FLORYAN, D. 2021 Exact coherent states and the nonlinear dynamics of wall-bounded turbulent flows. *Annu. Rev. Fluid Mech.* **53**, 227–253.

HOPF, E. 1948 A mathematical example displaying features of turbulence. *Commun. Pur. Appl. Math.* **1** (4), 303–322.

KAUSHIKA, N.D. & SUMATHY, K. 2003 Solar transparent insulation materials: a review. *Renew. Sust. Energy Rev.* **7**, 317–351.
KAWAHARA, G. & KIDA, S. 2001 Periodic motion embedded in plane Couette turbulence: regeneration cycle and burst. *J. Fluid Mech.* **449**, 291–300.

KAWAHARA, G., UHLMANN, M. & VAN VEEN, L. 2012 The significance of simple invariant solutions in turbulent flows. *Annu. Rev. Fluid Mech.* **44**, 203–225.

KRUPA, M. 1997 Robust heteroclinic cycles. *J. Nonlinear Sci.* **7** (2), 129–176.

LANFORD, O.E. 1982 The strange attractor theory of turbulence. *Annu. Rev. Fluid Mech.* **14**, 347–364.

LOHSE, D. & SHISHKINA, O. 2024 Ultimate Rayleigh–Bénard turbulence. *Rev. Mod. Phys.* **96**, 035001.

LUCAS, D. & KERSWELL, R.R. 2015 Recurrent flow analysis in spatiotemporally chaotic 2-dimensional Kolmogorov flow. *Phys. Fluids* **27** (4), 045106.

NAGATA, M. 1990 Three-dimensional finite-amplitude solutions in plane Couette flow: bifurcation from infinity. *J. Fluid Mech.* **217**, 519–527.

PAGE, J. & KERSWELL, R.R. 2020 Searching turbulence for periodic orbits with dynamic mode decomposition. *J. Fluid Mech.* **886**, A28.

PAGE, J., NORGAARD, P., BRENNER, M.P. & KERSWELL, R.R. 2024 Recurrent flow patterns as a basis for two-dimensional turbulence: predicting statistics from structures. *Proc. Natl. Acad. Sci. USA* **121** (23), 1–10.

REDFERN, E.M., LAZER, A.L. & LUCAS, D. 2024 Dynamically relevant recurrent flows obtained via a nonlinear recurrence function from two-dimensional turbulence. *Phys. Rev. Fluids* **9**, 124401.

REETZ, F., KREILOS, T. & SCHNEIDER, T.M. 2019 Exact invariant solution reveals the origin of self-organized oblique turbulent-laminar stripes. *Nat. Commun.* **10** (1), 1–6.

REETZ, F. & SCHNEIDER, T.M. 2020a Invariant states in inclined layer convection. Part 1. Temporal transitions along dynamical connections between invariant states. *J. Fluid Mech.* **898**, A22.

REETZ, F., SUBRAMANIAN, P. & SCHNEIDER, T.M. 2020b Invariant states in inclined layer convection. Part 2. Bifurcations and connections between branches of invariant states. *J. Fluid Mech.* **898**, A23.

ZHENG, Z. 2025 Pattern formation in transitional turbulent thermal convection: invariant solutions and their bifurcation structures. PhD thesis, Ecole Polytechnique Fédérale de Lausanne.

ZHENG, Z., BECK, P., YANG, T., ASHTARI, O., PARKER, J.P. & SCHNEIDER, T.M. 2025 Ghost states underlying spatial and temporal patterns: how nonexistent invariant solutions control nonlinear dynamics. *Phys. Rev. E* **112**, 024212.

ZHENG, Z., TUCKERMAN, L.S. & SCHNEIDER, T.M. 2024a Natural convection in a vertical channel. Part 1. Wavenumber interaction and Eckhaus instability in a narrow domain. *J. Fluid Mech.* **1000**, A28.

ZHENG, Z., TUCKERMAN, L.S. & SCHNEIDER, T.M. 2024b Natural convection in a vertical channel. Part 2. Oblique solutions and global bifurcations in a spanwise-extended domain. *J. Fluid Mech.* **1000**, A29.