1. Introduction
A fluid subjected to a horizontal temperature gradient, often called natural or vertical convection, is encountered in a wide range of geophysical (Hart Reference Hart1971), meteorological and engineering (Kaushika & Sumathy Reference Kaushika and Sumathy2003; Arici, Karabay & Kan Reference Arici, Karabay and Kan2015) applications. Scientific research on natural convection with its many variants has a long history. Motivated by the crucial application of thermal insulation, Batchelor (Reference Batchelor1954) sought to determine the width that maximized the insulating properties of an air-filled cavity within a wall or window, i.e. the double-glazing problem; this solution was later amended by Poots (Reference Poots1958) and Gershuni, Zhukhovitskii & Tarunin (Reference Gershuni, Zhukhovitskii and Tarunin1966). Elder (Reference Elder1965) observed, both experimentally and theoretically, the oblique convection rolls that form in a tall enclosure that we will study here. These rolls, in particular their onset, were further investigated by, e.g., Eckert & Carlson (Reference Eckert and Carlson1961), Vest & Arpaci (Reference Vest and Arpaci1969), Gershuni & Zhukhovitskii (Reference Gershuni and Zhukhovitskii1972), Korpela, Gözüm & Baxi (Reference Korpela, Gözüm and Baxi1973) and Mizushima & Gotoh (Reference Mizushima and Gotoh1976). After nonlinear numerical simulations became feasible, a number of researchers studied the evolution of the number of rolls in an air-filled cavity with height-to-width aspect ratio eight to twenty by means of time integration (Roux et al. Reference Roux, Grondin, Bontoux and de Vahl Davis1980; Le Quéré Reference Le Quéré1990; Wakitani Reference Wakitani1998) or by Newton–Krylov methods (Mizushima & Tanaka Reference Mizushima and Tanaka2002; Salinger et al. Reference Salinger, Lehoucq, Pawlowski and Shadid2002; Gelfgat Reference Gelfgat2004; Xin & Le Quéré Reference Xin and Le Quéré2006). Among the phenomena that they observed were hysteresis, multiplicity of solutions, and a non-monotonic evolution in the number of rolls with Rayleigh number. Quasiperiodicity and chaos in a cavity whose height is twice its width have been studied by Oteski, Duguet & Pastur (Reference Oteski, Duguet and Pastur2014) and Oteski et al. (Reference Oteski, Duguet, Pastur and Le Quéré2015).
Attesting to its importance, natural convection has been chosen as a computational benchmark for evaluating the accuracy and efficiency of fluid dynamics codes. The benchmark competition on an air-filled square cavity (de Vahl Davis & Jones Reference de Vahl Davis and Jones1983) attracted thirty-seven contributions, comparing results from codes described in, for example, Gilly, Bontoux & Roux (Reference Gilly, Bontoux and Roux1981), Winters (Reference Winters1982, Reference Winters1987), Le Quéré & Alziary de Roquefort (Reference Le Quéré and Alziary de Roquefort1985) and Le Quéré (Reference Le Quéré1991). An entire conference and journal volume were devoted to the benchmark problem of an air-filled height-to-width-ratio of eight (Christon, Gresho & Sutton Reference Christon, Gresho and Sutton2002).
Continuing our survey of archetypal configurations in vertical convection, low-Prandtl-number liquid metals are used in the process of growing semiconductor crystals; the goal is to avoid transition to oscillatory flow that engenders imperfections. A shallow cavity of height-to-width-ratio $1\,{:}\,4$ filled with a low-Prandtl-number liquid metal was the topic of yet another benchmark (Roux Reference Roux1990). Bifurcation analyses of this configuration have been carried out by, e.g., Winters (Reference Winters1988), Pulicani et al. (Reference Pulicani, Del Arco, Randriamampianina, Bontoux and Peyret1990), Henry & Buffat (Reference Henry and Buffat1998) and Gelfgat, Bar-Yoseph & Yarin (Reference Gelfgat, Bar-Yoseph and Yarin1999). We will continue our survey of the literature in Zheng, Tuckerman & Schneider (Reference Zheng, Tuckerman and Schneider2024), where we will focus on three-dimensional patterns.
Vertical convection is a special case of inclined layer convection, in which the container is tilted against gravity so that both buoyancy and shear forces drive the flow (Poots Reference Poots1958; Fujimura & Kelly Reference Fujimura and Kelly1993; Daniels, Plapp & Bodenschatz Reference Daniels, Plapp and Bodenschatz2000; Subramanian et al. Reference Subramanian, Brausch, Daniels, Bodenschatz, Schneider and Pesch2016; Grayer et al. Reference Grayer, Yalim, Welfert and Lopez2020). Extrapolating in inclination angle from the well-understood buoyancy-driven Rayleigh–Bénard case to shear-dominated vertical convection may give insights into transition in pure shear flows such as plane Couette flow (Nagata & Busse Reference Nagata and Busse1983). This was one of the motivations for the extensive study of inclined layer convection by Reetz & Schneider (Reference Reetz and Schneider2020) and Reetz, Subramanian & Schneider (Reference Reetz, Subramanian and Schneider2020), whose spirit and methods are carried over to the present study.
Rayleigh–Bénard convection, in which the layer is horizontal, is probably the most studied case of inclined layer convection, but it is exceptional in several important respects: the Prandtl number $Pr$ (ratio of kinematic viscosity to thermal diffusivity) plays no role at threshold, and the primary instability is always steady. In contrast, Korpela et al. (Reference Korpela, Gözüm and Baxi1973) showed that the primary instability in vertical convection is steady for $Pr<12.7$ and oscillatory for $Pr>12.7$, a value that was refined to $Pr=12.45$ by Fujimura & Kelly (Reference Fujimura and Kelly1993).
Rayleigh–Bénard convection is also exceptional in that its basic state is motionless, so that lateral boundaries can be assumed to affect only the regions immediately adjacent to them. The interior of a finite domain can therefore be approximated as homogeneous in the horizontal directions parallel to the rigid boundaries, so periodic boundary conditions can be used in these directions. In contrast, in vertical convection, the basic state includes a velocity that is vertical and hence normal to boundaries situated at the top and bottom. Such boundaries can have a substantial influence on the basic solution in the bulk if the aspect ratio is low or the Rayleigh number is high. This undermines the approximation of vertical homogeneity, without which theoretical or numerical treatment becomes much more difficult. Batchelor (Reference Batchelor1954), Eckert & Carlson (Reference Eckert and Carlson1961), Gill (Reference Gill1966), Vest & Arpaci (Reference Vest and Arpaci1969), Mizushima & Gotoh (Reference Mizushima and Gotoh1976) and Bergholz (Reference Bergholz1978) distinguished two regimes for the basic solution in the bulk of a finite cavity: the conductive regime, in which the temperature depends only on the distance from the walls, and the double boundary layer regime, in which the temperature also has a vertical gradient resulting from the flow meeting the upper and lower boundaries. The researchers cited above showed that even in the boundary layer regime, a cavity of finite height can be approximated by a vertically homogeneous problem if modified boundary conditions are imposed on the temperature at the two vertical bounding plates, either a finite vertical gradient or else horizontal isoflux conditions (Kimura & Bejan Reference Kimura and Bejan1984; Le Quéré Reference Le Quéré2022). The configuration that we study here, with aspect ratio 10 and Rayleigh numbers lower than 14 000, falls safely into the conductive regime. This means that our simulations using periodic vertical boundary conditions and bounding plates each of constant temperature resemble the flow in the interior regions of cavities of finite height. For a full treatment of the differences between periodic domains and those with boundaries (free-slip or rigid), see Hirschberg & Knobloch (Reference Hirschberg and Knobloch1997).
Our investigation is based on a series of studies by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015, Reference Gao, Podvin, Sergent, Xin and Chergui2018). These authors used direct numerical simulations (DNS) combined with linear and weakly nonlinear approaches to study vertical convection in air ($Pr=0.71$). By systematically increasing the Rayleigh number, Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) surveyed the regimes in a three-dimensional domain whose periodic vertical height was ten times that of the other two. They observed that the flow transitioned from the conductive state to steady rolls, then to oscillatory flow, and finally to a chaotic state. After acquiring four identical stacked co-rotating rolls, the flow continued to have a vertical periodicity of a quarter of the domain length over a fairly large Rayleigh number range. By subsequently confining the domain to this height to suppress large-scale instabilities, Gao et al. (Reference Gao, Podvin, Sergent and Xin2015) observed a period-doubling cascade leading to chaotic dynamics as the Rayleigh number was increased. However, a quantitative numerical bifurcation analysis corresponding to these studies has not been performed, thus the bifurcation-theoretic origins of the observed complex flow patterns remain to be fully explored. This motivates the present study.
We consider the domain $[L_x,L_y,L_z]=[1,1,10]$, where $x$, $y$ and $z$ represent the direction between the two bounding plates, the transverse direction, and the direction of gravity, respectively, as shown in figure 1. Here, only one of the three spatial directions is extended, thus the resulting flow structures, while fascinating and surprisingly complex, predominantly vary only in the vertical direction. Although the domain has only one spatially extended direction, weakly two-dimensional patterns have also been observed. Note, however, that all computations of Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015, Reference Gao, Podvin, Sergent, Xin and Chergui2018) as well as those presented here are fully three-dimensional. The domain $[L_x,L_y,L_z]=[1,8,9]$ studied by Gao et al. (Reference Gao, Podvin, Sergent, Xin and Chergui2018) has two extended directions and correspondingly fully two-dimensional patterns and behaviour. A bifurcation-theoretic analysis of these will be the subject of our companion paper Zheng et al. (Reference Zheng, Tuckerman and Schneider2024).
Above onset, as the Rayleigh number is increased, a sequence of convective patterns emerges from the conductive state. At each bifurcation point, symmetries of the previous state are in general sequentially broken, leading to patterns of increasing complexity. Those sequentially broken symmetries include continuous or $n$-fold translation symmetry, reflection symmetry, centro-symmetry and so on. The transition to $n$-fold translation symmetry in an effectively one-dimensional and reflection-symmetric domain generically leads to $D_n$ symmetry. The phenomenon of competition between steady branches with different wavenumbers is the essence of the Eckhaus instability (Eckhaus Reference Eckhaus1965; Tuckerman & Barkley Reference Tuckerman and Barkley1990), especially in the idealized case of long domains. For the particular finite vertical aspect ratio 10 in our convection problem, four co-rotating rolls are favoured, competing with three rolls at increasing $Ra$. As it happens, symmetry groups $D_4$ and $D_3$ have very particular properties; it is this combination of group theory, topology and fluid mechanics that shapes the resulting bifurcation diagram. The competition between three and four rolls is also manifested by branches of time-dependent states in which the number of rolls alternates periodically or chaotically.
More generally, Dangelmayr (Reference Dangelmayr1986) carried out a comprehensive investigation using weakly nonlinear model equations of the scenarios resulting from the competition between periodic patterns with different wavenumbers. Crawford, Knobloch & Riecke (Reference Crawford, Knobloch and Riecke1990) applied similar equations to a Faraday wave experiment. Among the features of these scenarios are steady pure-mode (single wavenumber) and mixed-mode (multiple wavenumber) branches, as well as travelling and standing waves. One of the main topics of our investigation is the numerical simulation and qualitative interpretation of the mixed-mode branches in our hydrodynamic system.
We begin by reproducing the equilibria and periodic orbits computed by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013). By following the branches to which these solutions belong, we discover new solution branches and identify the bifurcations giving rise to all of them, from onset to the chaotic regime. The remainder of the paper is structured as follows. In § 2, we discuss the governing equations, numerical aspects, symmetries, and the measurements and visualizations to be presented. Section 3 presents the steady solutions or equilibria, together with a detailed interpretation of the observed bifurcation scenarios using $D_4$ and $D_3$ symmetry (Gambaudo Reference Gambaudo1985; Swift Reference Swift1985; Knobloch Reference Knobloch1986; Golubitsky, Stewart & Schaeffer Reference Golubitsky, Stewart and Schaeffer1988; Dawes Reference Dawes2005). Time-periodic solutions are presented in § 4. Finally, we conclude with a summary of key results and a discussion in § 5.
2. Vertical convection system and numerical aspects
2.1. Governing equations
We used the inclined layer convection (ILC) extension module of the MPI-parallel pseudo-spectral code Channelflow 2.0 (Reetz Reference Reetz2019; Gibson et al. Reference Gibson, Reetz, Azimi, Ferraro, Kreilos, Schrobsdorff, Farano, Yesil, Schütz and Culpo2021) to carry out DNS of the non-dimensionalized Oberbeck–Boussinesq equations
in a vertical channel as depicted in figure 1. In (2.1), $\boldsymbol {u} = [u, v, w](x,y,z,t)$ and $\mathcal {T}=\mathcal {T}(x,y,z,t)$ stand for total velocity and temperature, respectively. The constant buoyancy term has been omitted from (2.1a); correspondingly, the pressure $p=p(x,y,z,t)$ is relative to the hydrostatic pressure. Bold symbols denote vector quantities, and $\boldsymbol {e}_z$ is the vertical unit vector. The equations have been non-dimensionalized with respect to the temperature difference $\Delta \vartheta$ and the distance $W$ between the walls, and the free-fall time unit $(W/g\alpha \,\Delta \vartheta )^{1/2}$, where $\alpha$ is the thermal expansion coefficient, and $g$ is the gravitational acceleration. Two independent dimensionless parameters appear: Rayleigh number $Ra = g \alpha \,\Delta \vartheta \, W^3/(\nu \kappa )$ and Prandtl number $Pr = \nu /\kappa$, where $\nu$ is the kinematic viscosity, and $\kappa$ is the thermal diffusivity.
Periodic boundary conditions are imposed in the $y$ and $z$ directions with spatial periods $L_y$ and $L_z$, respectively. The walls are no-slip and have prescribed temperatures
A supplementary integral constraint on either the pressure gradient or the mean flux must be set in the periodic directions. In order to match the simulations of Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015, Reference Gao, Podvin, Sergent, Xin and Chergui2018), we impose a mean pressure gradient of zero in $y$ and in $z$. Equations (2.1) together with the boundary conditions admit the conductive solution sketched in figure 1(a):
with arbitrary pressure constant $\varPi$.
2.2. Numerical methods
Channelflow-ILC adopts Chebychev–Fourier–Fourier (in $x$, $y$ and $z$) expansions for representing flow fields in space, and a finite differencing method for time integration (see detailed description in Appendix A of Reetz & Schneider Reference Reetz and Schneider2020). We have simulated the three-dimensional computational domain studied in Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013). This narrow domain $[L_x, L_y, L_z] = [1, 1, 10]$ is discretized by $[N_x, N_y, N_z] = [31, 32, 96]$ collocation points, resulting in a state space dimension $N=4 \times N_x \times N_y \times N_z \times (\frac {2}{3})^2$ of the order of $2\times 10^5$. The factor four stems from three components of velocity field and one in temperature field, and $(\frac {2}{3})^2$ is due to dealiasing in two Fourier directions (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang2006). Although our resolution is slightly less than that reported in Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013), we find it to be sufficient, since the ratio of the $L_2$-norm of the last resolved mode to the first mode of the velocity and temperature fields is less than $10^{-6}$ in the $y$ and $z$ directions, and less than $10^{-9}$ in the $x$ direction, a criterion also employed by Gibson & Schneider (Reference Gibson and Schneider2016).
As an extension to the studies based on DNS observations (Gao et al. Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015, Reference Gao, Podvin, Sergent, Xin and Chergui2018), our objective is to construct the invariant solutions such as equilibria and periodic orbits underlying the complex spatio-temporal flow dynamics. For identifying linearly stable states, time-marching (DNS) appropriate initial conditions gives access to these solutions, which is how Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015, Reference Gao, Podvin, Sergent, Xin and Chergui2018) proceeded. However, the root-finding technique is required for constructing unstable states. Invariant solutions are state vectors $\boldsymbol {x}^{*}(t)$ satisfying
where $\sigma$ is a symmetry operator, and $\mathcal {F}^T$ is the time-evolution operator integrating (2.1a)–(2.1c) from an initial state $\boldsymbol {x^*}$ over a finite time period $T$ (where $T$ is the period of a periodic solution, which is arbitrary for a steady solution). The shooting-based Newton–Raphson method in Channelflow-ILC uses a matrix-free Krylov method in which successive Krylov vectors are generated by time-marching initial conditions (Kelley Reference Kelley2003; Sánchez et al. Reference Sánchez, Net, García-Archilla and Simo2004). It is usually combined with a hookstep trust-region optimization based on the Krylov vectors, leading to a greatly increased radius of convergence (Viswanath Reference Viswanath2007, Reference Viswanath2009). Convergence is considered to be reached once the norm of the residual of (2.4) is sufficiently close to machine precision (of the order of $10^{-12}$). The converged solutions are subsequently continued parametrically along a range of Rayleigh numbers to form bifurcation diagrams (Sánchez et al. Reference Sánchez, Net, García-Archilla and Simo2004; Dijkstra et al. Reference Dijkstra2014) so as to understand their bifurcation structure.
The stability of each converged state is evaluated by using the Arnoldi algorithm (Arnoldi Reference Arnoldi1951; Antoulas Reference Antoulas2005) to determine its leading eigenvalues and eigenvectors for fixed points, or Floquet exponents and Floquet modes for periodic orbits. In a highly symmetric problem like this one, most eigenvalues are multiple, since symmetry operations applied to non-symmetric eigenvectors can yield other eigenvectors. For multiple eigenvalues, the Arnoldi algorithm returns an arbitrary set of linearly independent eigenvectors. We take linear combinations of these to construct those eigenvectors within the respective eigenspaces that are appropriate for our purposes.
2.3. Symmetries of the system
The vertical convection system is equivariant under $y$-reflection (2.5a), combined $x$- and $z$-reflection (2.5b), and translation in $y$ and $z$ (2.5c):
Since $y$ and $z$ are periodic directions, the centre of reflection $(y_0,z_0)$ is arbitrary, so reflections $y\rightarrow -y$ and $z\rightarrow -z$ in (2.5a) and (2.5b) should be more generally written as $y_0+y\rightarrow y_0-y$ and $z_0+z\rightarrow z_0-z$ for some $y_0$ and $z_0$. We will write these merely as ${\rm \pi} _y$ and ${\rm \pi} _{xz}$, while for visualizations we will choose whatever axis of reflection seems most appropriate for $y_0$ and $z_0$, usually $L_y/2$ and $L_z/2$.
The symmetry transformations (2.5) form the equivariance group of the system, which consists of all products of the generators $S_{VC} \equiv \langle {{\rm \pi} _y, {\rm \pi}_{xz}, \tau (\Delta y, \Delta z)} \rangle \sim [O(2)]_y \times [O(2)]_{xz}$. (Although symmetry groups cannot always be associated with only $y$ or $(x,z)$, we will do so occasionally when this is convenient and possible.) The groups that arise in this study are $Z_n$, the cyclic group of $n$ elements, $D_n$, the cyclic group of $n$ elements together with a non-commuting reflection, and $O(2)$, the group of all rotations (or equivalently translations in our periodic domain) together with a non-commuting reflection. We note that $D_1 =Z_2$ and $D_2=Z_2\times Z_2$. Aside from the conductive solution, which is invariant under the full group $S_{VC}$, other solutions may be invariant only under proper (smaller) subgroups of $S_{VC}$. Trajectories that begin in an invariant subspace remain so under exact arithmetic, but may depart due to instability. At times in this study, we have imposed reflection symmetries or periodicity over an interval shorter than $L_y$ or $L_z$ in order to restrict the dynamics to the desired invariant subspace or to expedite numerical continuation.
2.4. Numerical measurements and visualizations
We define the deviation from the conductive solution $\theta \equiv \mathcal {T} - \mathcal {T}_0$, which we will usually refer to merely as the temperature, and employ its $L_2$-norm
as an observable for plotting the bifurcation diagrams. For fixed points, a single curve representing $\| \theta \|_2$ as a function of the Rayleigh number is plotted. For periodic orbits, the maximum and minimum of $\| \theta \|_2$ along an orbit are plotted, resulting in two different curves representing one solution. Multiple solutions related by symmetry, in particular those resulting from pitchfork bifurcations, share the same value of $\| \theta \|_2$. In order to distinguish between symmetry-related flow fields, we use a local measurement $\theta _{local}$ based on the temperature at a single point. Here and in Zheng et al. (Reference Zheng, Tuckerman and Schneider2024), the bifurcation diagrams contain apparent intersections of curves indicating solution branches that are not related to bifurcations but result from projecting the high-dimensional flow fields onto a one-dimensional scalar quantity. Apparent intersections that are not labelled as bifurcations are of this spurious type.
In addition, we also calculate the thermal energy input ($I$) due to buoyancy forces, and the dissipation ($D$) due to viscosity, both averaged over the domain, for phase portrait visualizations. We refer readers to Reetz & Schneider (Reference Reetz and Schneider2020) for more details. In order to visualize instantaneous flow fields or eigenvectors, we plot their temperature fields $\theta$ on the $y$–$z$ plane on the midplane at $x=0$, and/or on the $x$–$z$ plane at $y=0.5$.
3. Equilibria
Our goal is to understand the formation and instabilities of convection rolls in the computational domain $[L_x, L_y, L_z] = [1, 1, 10]$, the domain studied by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013). Figure 2 displays the equilibria that we have studied. Many more unstable branches undoubtedly exist that are not shown in this figure, since a new branch is formed whenever the real part of an eigenvalue traverses zero. Since some of the states that we discuss can also exist in domains $[1, 1, 2.5]$ and $[1, 0, 10]$, we will also mention their existence and stability ranges in these smaller domains.
3.1. Two primary and one secondary circle pitchfork bifurcations
The conductive base flow becomes linearly unstable at $Ra=5826$, close to the threshold $Ra=5800$ reported by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013), where it bifurcates to a two-dimensional state containing four co-rotating transverse convection rolls. Each roll has height (or wavelength) $\Delta z=L_z/4=10/4=2.5$, and we will use both decimal and fractional notation as seems appropriate. The critical wavelength and Rayleigh number for $Pr=0.71$ computed by Vest & Arpaci (Reference Vest and Arpaci1969) are $2.37$ and $5595$, respectively; since our wavelength is constrained by our imposed vertical periodicity to be a divisor of 10, the threshold in $Ra$ is necessarily higher.
The four-roll state, called FP1 in figure 2(a), is illustrated in figures 2(d) and 1(c), which show the temperature field $\theta$. We recall that we have defined $\theta$ to be the deviation from the conductive solution, which we show in figure 1(b); the full temperature field is shown in figure 1(d). Examination of figure 1 along with the corresponding velocity fields shows that the motion of the deviation fields is clockwise (figure 1c), but that when added to the base flow (figure 1b), the full motion in each roll (figure 1d) is anticlockwise: colder fluid on the left ($x=-0.5$) crosses the cavity towards the right and then rises, while warmer fluid on the right ($x=0.5$) crosses towards the left and then descends. This instability is driven by the shear in the vertical velocity, in contrast to the buoyancy-driven rolls that occur in Rayleigh–Bénard convection. Here, FP1 has reflection and translation symmetries $S_{FP_{1}} \equiv \langle {{\rm \pi} _y,{\rm \pi} _{xz}, \tau (\Delta y,2.5)}\rangle \sim [O(2)]_y \times [D_4]_{xz}$, where the translation symmetry in $L_z/4=2.5$ results from its four vertically stacked identical rolls in figure 2(d).
We have found another fixed point, FP2, containing three identical rolls, which is shown in figure 2( f). Fixed point FP2 bifurcates from the unstable conductive base flow at $Ra=6868.7$ and remains unstable over its entire range of existence. It is invariant under reflection and translation symmetries $S_{FP_{2}} \equiv \langle {{\rm \pi} _y,{\rm \pi} _{xz}, \tau (\Delta y,10/3)}\rangle \sim [O(2)]_y \times [D_3]_{xz}$. Fixed point FP1 is stable until $Ra=10\,166$, when it bifurcates to a state containing four three-dimensional (3-D) steady rolls, which we have called FP3. This state, observed by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) and shown in figure 2(e), is in turn stable until $Ra=11\,261$. Fixed point FP3 is invariant under $S_{FP_{3}} \equiv \langle {{\rm \pi} _y, {\rm \pi}_{xz}\tau (0.5,0), \tau (0,2.5)}\rangle \sim [D_1]_y \times [D_4]_{xz}$; symmetry $\tau (\Delta y,0)$ is broken at the circle pitchfork bifurcation point at $Ra = 10\,166$.
3.2. Fixed point FP4: connector between FP1 and FP2 states
Figures 2(a–c) show another equilibrium, which we have called FP4, bifurcating from FP1 at $Ra= 13\,383.9$, and intersecting FP2 at $Ra=11\,283$. Two sets of solutions, figures 2j,k, appear from the same FP1 state via simultaneous subcritical pitchfork bifurcations. We will call these half-branches; the reason for this and for their simultaneous bifurcation will become clear below.
We will need to consider another translation-symmetry-related version of FP1, shifted by a half-roll ($\Delta z=\pm 1.25$) from FP1, which we will call FP1$^\prime \equiv \tau (0, 1.25)$ FP1. Because the global quantity $\|\theta \|_2$ cannot distinguish between symmetry-related states, we represent FP4 in figure 3(a) by the local and normalized quantity
To emphasize the variation of $\theta _{local}$ as FP4 is traversed, for visualization, we suppress the variation along the FP1 and FP2 branches.
The two endpoints of the FP4 branch, related by a half-roll shift of $\Delta z = 1.25$, are shown in figures 3( j) and 3(k). In the bifurcations from FP1 to FP4, the fourfold translational symmetry in $z$ is lost, but $(x,z)$ reflection symmetry is retained, leading to $S_{FP_4} \equiv \langle {{\rm \pi} _y, {\rm \pi}_{xz}, \tau (\Delta y, 0)} \rangle \sim [O(2)]_y \times [Z_2]_{xz}$. We have chosen the spatial phase such that the centres of symmetry of figures 3( j) and 3(k) are located at $z$ values that are multiples of $10/8=1.25$. During the numerical continuation of the FP4 branch, the phase in $z$ has been fixed by imposing two reflection symmetries.
3.2.1. Tour of FP4: two methods for eliminating one roll
We begin our tour of FP4 from figure 3( j), which displays one end of the FP4 branch, or equivalently, FP1. Going from figure 3( j) to figure 3(i), the between-roll boundary at $z=2.5$ becomes weaker. In contrast, at $z=7.5$ the roll boundary is strengthened, while the far edges of the two surrounding rolls are weakened. By figure 3(h), the two rolls formerly surrounding $z=2.5$ have merged into a single large roll. For this reason, we call ( j,i,h,g,f) in figure 3(a) the roll-merging half-branch. Starting from FP1$^\prime$ in figure 3(k), i.e. the opposite endpoint of the FP4 branch, the roll centred around $z=7.5$ weakens in figure 3(l) and has almost disappeared by the saddle–node bifurcation of figure 3(m). We call (k,l,m,n,f) in figure 3(a), the roll-disappearing half-branch.
At $Ra=11\,283$, figure 3( f) has three equally spaced rolls and belongs to branch FP2. This is why we choose to call this state the dividing point of branch FP4 into two half-branches. The meeting between FP2 and the two half-branches is a transcritical bifurcation that will be the topic of § 3.2.4. Both FP4 half-branches lead from four rolls to three rolls, but in different ways. In the pathway from figure 3( j) to figure 3( f), the space between two rolls blurs, and the two rolls merge. In the pathway from figure 3(k) to figure 3( f), one roll weakens and disappears. These two types of transitions can occur at any of the four roll centres and roll boundaries. Thus eight half-branches bifurcate simultaneously from any FP1 state: four roll-merging half-branches like figures 3( j)–3( f), and four roll-disappearing half-branches like figures 3(k)–3( f). These eight branches connect an FP1 state with its half-roll-shifted state FP1$^\prime$.
This scenario is schematized in figure 4. Each line of longitude (meridian) on the globe-like figure represents a branch connecting FP1 (top square) and FP1$^\prime$ (bottom square), like that shown in figure 3. The roll-merging half-branches are coloured in red and emerge from the corners of a square; the roll-disappearing half-branches in blue emerge from the sides of a square. The fact that four of each emerge at each of the squares corresponds to the fact that each of FP1 and FP1$^\prime$ contains four rolls and four inter-roll spaces that can undergo roll disappearance or roll merging. Each half-branch of one colour emanating from FP1 meets a half-branch of the opposite colour emanating from FP1$^\prime$ at the equator, which contains transcritical bifurcation points of different phases in $z$.
3.2.2. Eigenvectors of FP1 and FP1$^\prime$
Figures 3(b) and 3(c) show the unstable eigenmodes $e_3$ of FP1 and $e_4^\prime$ of FP1$^\prime$ at $Ra=13\,384$ that are responsible for the two simultaneous subcritical pitchfork bifurcations that create the two half-branches of FP4. We call these $e_3$ and $e_4^\prime$ because the FP1 (and FP1$^\prime$) branch at $Ra=13\,384$ has two larger positive eigenvalues resulting from the circle pitchfork bifurcation to FP3. Eigenvectors $e_3$ and $e_4^\prime$ have the same eigenvalue, $\lambda _{3,4}$. The Arnoldi method computes two linearly independent eigenmodes of the double eigenvalue $\lambda _{3,4}$; for $e_3$ we select the linear combination of these that most resembles the difference between FP4 and FP1 at the bifurcation point of figure 3( j). The eigenvectors of FP1$^\prime$ are related to those of FP1 by a translation $\tau (0,1.25)$. We select as $e_4^\prime$ the analogous combination of eigenmodes of FP1$^\prime$ that most resembles the difference between FP4 and FP1$^\prime$ at the bifurcation point of figure 3(k). Eigenmodes $e_3$ and $e_4^\prime$ differ qualitatively: $e_3$ has two narrow intense temperature extrema surrounded by wide diffuse patches of the opposite sign, while $e_4^\prime$ has two wide diffuse patches surrounded by narrow extrema. Each eigenmode has two centres of ${\rm \pi} _{xz}$ symmetry, at $z=2.5$ and $7.5$.
Examining the eigenvectors helps us to understand the progression from the fourfold translation-symmetric FP1 (and FP1$^\prime$) to FP4. The eigenvectors describe defects that, when added to FP1 (or FP1$^\prime$), lead to roll merging or roll disappearance. The red (hot) and blue (cold) diffuse patches of $e_3$ are in opposition to those of FP1 at the boundary between two rolls at $z=2.5$; compare figures 3(b) and 3( j). This implies that the between-roll boundary at $z=2.5$ becomes weaker along this half-branch. In contrast, at $z=7.5$, FP1 and $e_3$ have temperatures of the same sign, so this roll boundary is strengthened. Turning now to the pathway (k,l,m,n,f), this roll-disappearing half-branch is associated with eigenvector $e_4^\prime$ in figure 3(c). Eigenvector $e_4^\prime$ is very weak at $z=2.5$ and at $z=7.5$, around which rolls of FP1$^\prime$ are centred. However, the temperature of $e_4^\prime$ surrounding $z=2.5$ is such as to reinforce the roll at $z=2.5$ of FP1$^\prime$, whereas $e_4^\prime$ and FP1$^\prime$ display opposite temperatures surrounding $z=7.5$. The roll of FP1$^\prime$ at $z=7.5$ consequently disappears along this half-branch.
Eigenmodes of the Jacobian matrix describe the temporal dynamics near a fixed point $\bar {u}$, but we have used them above to describe the tangent along a branch (or half-branch) near a bifurcation. We now explain the justification for this. For a dynamical system ${\rm d}u/{\rm d}t=f(u,Ra)$, a curve of fixed points $\bar {u}(Ra)$ is defined via $0=f(\bar {u}(Ra),Ra)$. Differentiating in $Ra$ yields
where $[{D}f]_{\bar {u}}$ is the Jacobian evaluated at $\bar {u}$. Near a bifurcation, the Jacobian has an eigenvalue $\lambda _{bif}$ near zero so that multiplication by the inverse Jacobian projects onto the bifurcating eigenvector $e_{bif}$:
where $\langle \cdot,\cdot \rangle$ is an inner product, and $(\lambda _j,e_j)$ are the eigenpairs of $[Df]_{\bar {u}}$, with $|1/\lambda _{bif}| \gg |1/\lambda _j|$ for the other eigenvalues of $[Df]_{\bar {u}}$. This leads to the expression
for the evolution of a branch near a bifurcation.
3.2.3. Normal form of $D_4$ symmetry
The simultaneous occurrence of two pitchfork bifurcations described above is precisely the scenario seen in pattern formation on a square domain, which, like FP1 (and FP1$^\prime$), has the symmetry group $D_4$, generated by ${\rm \pi} _{xz}$ and $\tau (0,2.5)$. In the square, rolls can be oriented horizontally or vertically, and these are equivalent because they are related by a rotation by ${\rm \pi} /2$. The eigenvectors associated with vertical and horizontal rolls can also be combined to form diagonal eigenvectors. The nonlinear equations that are equivariant (compatible) with $D_4$ symmetry predict the existence of branches of diagonal states (Swift Reference Swift1985; Bergeon, Henry & Knobloch Reference Bergeon, Henry and Knobloch2001) that originate from eigenvectors that are equal superpositions of vertical and horizontal eigenmodes, as will be discussed below. The diagonal roll branches bifurcate simultaneously with the horizontal and vertical roll branches, but the nonlinear diagonal states are not related to the horizontal or vertical states by symmetry operations and are therefore not equivalent. Both types of branches have a reflection symmetry – vertical or horizontal in one case, and diagonal in the other case – so that their symmetry groups are $Z_2$.
This scenario for pattern formation on a square domain also exists for the FP1 branch, with the four co-rotating rolls playing the role of the four sides of a square, and the four inter-roll intervals playing the role of the corners. Instead of considering the FP4 branch with its endpoints FP1 and FP1$^\prime$ as we did in figure 3, we now consider a single phase of FP1 and its two bifurcations to roll-merging and roll-disappearing half-branches corresponding to its eigenvectors $e_3$ and $e_4$. Four bifurcating branches resembling figures 3( j)–3( f) result from eigenvector $e_3$ along with shifted and reflected versions, and four branches resembling figures 3(k)–3( f) result from $e_4$ along with shifted and reflected versions. The bifurcations occur at the same value of the Rayleigh number, but the branches are not equivalent, as seen in figures 2(a) and 3(a), for example, by the different locations of the saddle–node bifurcations emanating from FP1 and FP1$^\prime$.
We now give a more quantitative explanation of this scenario. Consider the dynamical system governing $(p,q)\in \mathcal {R}^2$:
with $\mu,a,b$ all real parameters. The bifurcation parameter is $\mu$, and $a,b$ are nonlinear coefficients that saturate the instability. System (3.5) is a projection of a larger system near a bifurcation onto the bifurcating eigenmodes. A normal form is the smallest system, in terms of both number of variables and polynomial order, that is able to reproduce the behaviour of the larger system near the bifurcation. The form of system (3.5) is dictated by the requirements that it be equivariant under (consistent with) change in sign of $p$ or $q$, and interchange of $p$ and $q$, which defines the group $D_4$.
The Jacobian of (3.5) is
Evaluated at the trivial solution $(p,q)=(0,0)$, this becomes $\mu {\boldsymbol {I}}$, i.e. a double eigenvalue $\mu$. The non-trivial fixed points of system (3.5) are
Thus (3.5) has eight non-trivial solutions, two each of types (3.7a), (3.7b), (3.7c) and (3.7d). Although solutions (3.7a) and (3.7b) are related to one another by the symmetry $(p,q) \rightarrow (-q,p)$, as are solutions (3.7c) and (3.7d), solutions (3.7c) and (3.7d) are not related to solutions (3.7a) and (3.7b) by interchanging $p$ and $q$, or by changing their signs.
The scenario by which FP1 gives rise to FP4 is analogous to system (3.5) and (3.7), with FP1 playing the role of the trivial solution $p=q=0$. The assumption of normal form theory is that FP4 solutions can be approximated as superpositions of the base state FP1 and its eigenvectors $e_3$ and $\tau (0,2.5)\,e_3$ at the bifurcation:
with $p(t)$ and $q(t)$ governed by the amplitude equations or normal form (3.5). The quantity $p$ measures the amplitude of eigenvector $e_3$ in figure 3(b), which gives rise to the half-branch in which two rolls merge at $z=2.5$. The phase-shifted $\tau (0,2.5)\,e_3$, whose amplitude is measured by $q$, gives rise to a different half-branch in which roll merging occurs at $z=5.0$. Figures 5(a,b) show these eigenvectors, while figure 5(c) shows their normalized sum. Further shifts, $\tau (0,5)\,e_3 = -e_3$ and $\tau (0,7.5)\,e_3 = -\tau (0,2.5)\,e_3$, correspond to $-p$ and $-q$, respectively.
Turning now to the four roll-disappearing half-branches bifurcating from FP1, these are produced by eigenvectors $\tau (0, 2.5 n)\,e_4$ for $n=0,1,2,3$. The normalized sum of the roll-merging eigenvectors $(e_3 + \tau (0, 2.5)\,e_3)/\sqrt {2}$ turns out to be the roll-disappearing eigenvector $e_4$, analogously to the fact that the sum of a $p$ solution and a $q$ solution yields a $p+q$ solution. (Because these are eigenvectors, their amplitudes have no importance.) The sum of two neighbouring roll-disappearing eigenvectors (not shown) is a roll-merging eigenvector, analogously to the combinations $(p+q)+(p-q)\propto p$ and $(p+q)-(p-q)\propto q$. This confirms the correspondence between the normal form (3.5) and our hydrodynamic system with its four-roll branch FP1, its connector branch FP4, and its eigenvectors $e_3$ and $e_4$.
3.2.4. Transcritical bifurcation between FP4 and $D_3$-symmetric FP2
Figure 3(a) shows the intersection between FP4 and the three-roll branch FP2 at figure 3( f) in a transcritical bifurcation (TC). From the point of view of FP2, the FP4 roll-merging half-branch ( j,i,h,g,f) can be called a roll-splitting half-branch when traversed in the opposite order ( f,g,h,i,j). Similarly, the FP4 roll-disappearing half-branch (k,l,m,n,f) can be called a roll-creation half-branch when traversed in the order ( f,n,m,l,k). Because FP2 has threefold translation symmetry $\tau (0,10/3)$, any of the three rolls in figure 3( f) can be the site of a roll-splitting or a roll-creation event, so six FP4 half-branches, three of each type, emanate from FP2 at TC. These join pairwise at FP2: for example, in figure 3, the upper half-branch with roll splitting at $z=2.5$ (figures 3h,i) meets the lower half-branch in which roll creation occurs at $z=7.5$ (figures 3m,l). (We assume that the saddle–node bifurcations have no effect on this scenario.)
The bifurcation from FP2 to FP4 breaks threefold translation symmetry but retains reflection symmetry ${\rm \pi} _{xz}$. This can be seen in figure 3(h), for example, where the roll centred at $z=2.5$ is stretched, while the other two rolls remain of the same size and related to one another by reflection in $z=2.5$. For figure 3(m), the reasoning is the same, but is applied to the inter-roll space at $z=7.5$.
We now turn to the eigenmodes of FP2 and FP4 close to TC. To the right of figure 3( f), figure 3(d) displays the eigenmode $e_5$ of FP2 at $Ra\gtrsim 11\,283$, leading to FP4. As previously, the name $e_5$ is used because FP2 has four eigenvalues with larger real parts. By a slight abuse of notation, we use $-e_5$ to denote the direction in which FP2 is approached from FP4 for $Ra\lesssim 11\,283$, and visualize it in figure 3(e). We obtain $\pm e_5$ by subtracting the temperature fields of FP2 from FP4 states to the right and left of point ( f) in figure 3(a), as well as from the Arnoldi method. Using (3.4) again, we can superpose FP2 in figure 3( f) with eigenvector $e_5$ in figure 3(d) to yield FP4 in figures 3(g) and 3(h), since $e_5$ opposes the roll in FP4 centred at $z=2.5$, leading to an expanded roll and eventually to roll splitting. Similarly, we can superpose figure 3( f) with eigenvector $-e_5$ in figure 3(e) to yield FP4 in figures 3(n) and 3(m). Since $-e_5$ opposes the rolls on either side of $z=7.5$ in FP4, this inter-roll space expands, eventually making room for roll creation.
As FP2 has threefold translation symmetry, $e_5$ of FP2 can be shifted by $\pm 10/3$ in $z$, yielding a triple of eigenvectors only two of which are linearly independent, since $\tau (0,10/3)\,e_5 + \tau (0,-10/3)\,e_5 = -e_5$. These share the same eigenvalue $\lambda _{5,6}$, depicted in figure 6(a). Along branch FP4, these eigenvectors are modified, so that they are no longer related by $\tau (0,\pm 10/3)$ and hence have different eigenvalues, shown as $\lambda _5$ and $\lambda _6$ in figure 6(b). Eigenvectors $e_5$ and $e_6$ of FP4 at $Ra=11\,292.2$, shown in figures 6(d) and 6(e), are symmetric and anti-symmetric, respectively, with respect to $xz$-reflection symmetry about $z=2.5$ and $z=7.5$.
3.2.5. $D_3$ symmetry
We now consider the equations governing bifurcation in the presence of $D_3$ symmetry. In this case, we will consider not the normal form, but a related system, i.e. the universal unfolding of the degenerate case of the normal form, because these are the equations that best describe our results. See Gambaudo (Reference Gambaudo1985), Golubitsky et al. (Reference Golubitsky, Stewart and Schaeffer1988) and Dawes (Reference Dawes2005) for details. These equations are
with $a, b$ real parameters, $\mu$ the bifurcation parameter, and $(p,q)$ the amplitudes of eigenmodes of the FP2 branch. As stated in § 3.2.3, the correspondence of (3.9) with our high-dimensional hydrodynamic system consists of approximating FP4 by a superposition of FP2 with its eigenvectors $e_5$ and $\tau (0, 10/3)\,e_5$ at the bifurcation point, whose amplitudes are represented here by $p$ and $q$:
Let us begin by studying steady solutions with $q=0$:
These are shown in figure 7(a). The trivial solution $p=0$ corresponds to the FP2 branch. Two sets of solutions corresponding to FP4 are created at $\mu =b^2/(4a)$; this is a saddle–node bifurcation. The set of solutions closer to zero intersects the trivial FP2-type branch $p=0$ in a transcritical bifurcation at $\mu =0$. These two bifurcations are marked as SN and TC in the parabola in figure 7(a). The saddle–node and transcritical bifurcations are also seen in the hydrodynamic case and are labelled by (g,f) in the bifurcation diagram of figure 3.
The Jacobian of (3.9) is
Evaluated at $(p,q)=(0,0)$ corresponding to FP2, this becomes $-\mu {\boldsymbol {I}}$, i.e. one double eigenvalue $-\mu$. Evaluated at the $q=0$, $\mu - bp+ ap^2=0$ solution corresponding to FP4, we obtain
Since $ap^2$ can be neglected near TC, the FP4-type solutions (3.11b) take the form $p\approx \mu /b$. (In accordance with the previous nomenclature, these are two half-branches, one for $\mu >0$, and the other for $\mu <0$.) The Jacobian becomes
Thus the FP4 states emanating from TC each have two eigenvalues of opposite signs, which are approximately $\mu$ (in the $p$ direction connecting FP2 and FP4) and $-3\mu$ (in the $q$ direction, perpendicular to $p$). This is precisely the behaviour seen in figures 6(a) and 6(b). Indeed, if we define $-\mu$ to be the eigenvalue of FP2 in figure 6(a), then we find that the eigenvalues of FP4 in figure 6(b) are approximately $\mu$ and $-3\mu$.
Dropping now the requirement that $q=0$, two more solutions to (3.9) of FP4-type exist, related to the $q=0$ solution by rotation by $\pm 2{\rm \pi} /3$ in the $p\unicode{x2013}q$ plane:
Each can be assigned an amplitude $\sqrt {p^2+q^2} \approx \mu /b$ and a phase $\tan ^{-1}(q/p)$. The phase can in turn be mapped to a vertical location in $[0,L_z)$ of a defect in one of the three rolls or inter-roll spaces. Such a defect is a precursor to roll splitting or roll creation, respectively, as one leaves TC along one of the half-branches. The eigenvalues of the other two FP4 solutions are again $\mu$ and $-3\mu$. The eigenvector associated with $\mu$ resembles the defect itself; i.e. it corresponds to a change in its amplitude. The eigenvector associated with $-3\mu$ corresponds to a change in phase, i.e. a tendency for the defect to translate in $z$. Like all eigenvectors, this tendency is local to TC, and nonlinear trajectories deviate from the phase-translation path before any phase change is actually achieved.
A schematic illustration of these solutions and the transcritical bifurcation is shown in figure 7(b). Roll-splitting and roll-creation half-branches are shown in red and blue, respectively. Three of each type of half-branch exist on each cone. The thick red and blue half-branches comprise the branch followed in figure 3 from FP1 to FP1$^\prime$, along which roll merging and then roll creation occur. Another two half-branches comprise a branch from FP1 to FP1$^\prime$ along which roll disappearance and then roll splitting occur. The other four branches have starting or ending points that are shifted by $\Delta z\pm 10/3$ from FP1 or FP1$^\prime$. These six branches are all included in figure 7(b) in order to give a full picture of the transcritical bifurcation; their depiction is not necessary for the understanding of the single path of figure 3.
3.3. Wider perspectives
In the previous subsections, we have extensively discussed the $D_4$ and $D_3$ bifurcation scenarios separately. We now take a wider perspective, discussing various aspects of the interaction between the $D_4$ and $D_3$ bifurcations.
The double-cone visualization of figure 7 may seem incompatible with the globe-like visualization of figure 4. In fact, each figure is local, and the two visualizations have only two branches in common. Each branch belongs simultaneously to a sphere and to a double cone. Four of the branches traversing TC through the double cone are not present in figure 4; two more spheres would be required to contain them. Similarly, more double cones would be required to contain all the meridians of the globe. A total of $2 \times 3 \times 8 = 48\text { half-branches} = 24\text { branches}$ are necessary to close the system, i.e. to include all branches that emanate from all bifurcations encountered by branches created at FP1. (Other phase changes in $z$ lead to a continuous infinity of branches.) It is not possible (or we have not been able) to depict the entire scenario in a single diagram. However, we again emphasize that figure 3 can be understood without recourse to this large number of other symmetry-related branches.
We have depicted FP4 as connecting two versions of FP1 related by a shift $\Delta z=1.25$, while passing continuously through a transcritical bifurcation at FP2. However, there exists another construction of this scenario: the transcritical bifurcation could be broken apart in such a way that rather than traversing FP2 smoothly, FP4 enters the transcritical bifurcation at FP2 but then exits at FP2$^\prime \equiv \tau (1.25,0)$ FP2. In figure 3, the second row would contain not a repeated version of figure 3( f), but instead a shifted version of it. Figures 3(n,m,l,k) would then also be shifted, with the result that figure 3(k) would be identical to figure 3( j) instead of being a shifted version of it. In the schematic figure 7, the left cone and right cone would each be reflected in such a way as to separate the two vertices and to join the two bases. In the schematic figure 4, the equator would be cut open, while the north pole and a rotated south pole would be joined. In summary, the FP4 branches start and end in states with a shift $\Delta z=1.25$ between them, but this can happen either by connecting FP1 and FP1$^\prime$ while passing continuously through FP2, as discussed in §§ 3.2.1–3.2.5, or alternatively, by connecting FP2 and FP2$^\prime$ while passing continuously through FP1. Yet another alternative point of view is to double the FP4 branch, passing continuously through FP1, FP2, FP1$^\prime$, FP2$^{\prime }$, FP1 without any phase jumps.
Dangelmayr (Reference Dangelmayr1986) determined the normal form for the occurrence of bifurcations to periodic patterns of two different wavenumbers, and its unfolding. The equations for the $D_3$–$D_4$ mode interaction predict a number of the features that we observe for our branch FP4 connecting the four-roll and three-roll branches FP1 and FP2, such as the existence of two qualitatively different connecting branches (like our roll-disappearing and roll-merging branches), a nearby saddle–node bifurcation on one of them (like that of figure 3g), and a Hopf bifurcation giving rise to temporally periodic solutions (such as the PO1 branch to be discussed in § 4.2). Some of the features of our scenario are not present in the normal form, in particular the subcriticality of the bifurcations from FP1 to FP4 and the possibly related two additional saddle–node bifurcations of FP4. Another feature that is not mentioned in Dangelmayr (Reference Dangelmayr1986) is the involvement of the two $L_z/8$ phase-shifted versions of FP1. The normal form predicts the stable and unstable eigenvalues of the solution branches, such as those discussed in the next subsection.
3.4. Stability
We start by discussing the stability of the FP4 state, which changes multiple times along its branch. Bifurcating subcritically from FP1, the upper (roll-merging) branch of FP4 is initially unstable and inherits the four unstable eigendirections of FP1 at $Ra=13\,383.9$ (two of which are the $y$-dependent modes that give rise to FP3). Three of these four positive eigenvalues are successively stabilized along the upper branch by undergoing tertiary bifurcations to quaternary states that are not discussed in this work. The last positive eigenvalue is stabilized after undergoing a saddle–node bifurcation at $Ra=8255$. This stability is short-lived, however, ending when FP4 undergoes a Hopf bifurcation at $Ra=9980$ to a periodic orbit to be discussed in the next section. The FP4 branch undergoes two more saddle–node bifurcations, at $Ra=11\,437$ and then at $Ra=10\,020$, before it finally terminates on FP1$^\prime$, a translated version of the FP1 branch that also has four unstable eigendirections. The connector branch FP4 is thus stable over the interval $8255< Ra<9980$, along with the four-roll branch FP1. However, it is not surprising that the FP4 branch was overlooked by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013), since it is unstable over most of its range of existence, and its bifurcation occurs from a point at which FP1 is unstable.
As mentioned in § 3.1, in domain $[L_x, L_y, L_z] = [1, 1, 10]$, branch FP1 is stable at onset, while FP2 is unstable. Both FP1 and FP2 exist in domain $[1, 0, 10]$ as well. After the bifurcation to FP1, the conductive state acquires two unstable $y$-independent eigenmodes. These are inherited by FP2 at onset, and so FP2 is also unstable when computed in domain $[1, 0, 10]$. Concerning the stability of FP1, since the bifurcation to FP3 breaks $y$-translation symmetry, it does not occur in domain $[1,0,10]$, and FP1 remains stable until the bifurcation to FP4 (in § 3.2) at $Ra=13\,383.9$. Regarding FP4, its range of stability does not change if computed in domain $[1,0,10]$, since the unstable part of the FP4 branch always has at least one unstable $y$-independent eigenmode. Due to their fourfold symmetry, FP1 and FP3 can exist in domain $[1, 1, 2.5]$, in which their existence and stability ranges are the same as in $[1, 1, 10]$. These ranges are summarized and compared in table 1 below. Since every zero-crossing of an eigenvalue or of its real part is accompanied by a bifurcation, there necessarily exist many more branches that we have not computed, for example along the FP4 branch.
4. Periodic orbits
Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) report that the flow becomes oscillatory at $Ra=11\,270$ through a Hopf bifurcation from the three-dimensional steady rolls (FP3), and that this is followed by a period-doubling bifurcation at approximately $12\,100 < Ra < 12\,200$. In our simulations of the domain $[1, 1, 10]$ in the range $9980 < Ra < 11\,270$, the flow can be either steady or time-periodic, depending on the initial conditions. This suggests that in addition to the stable steady solution FP3, a limit cycle also exists in this configuration and Rayleigh number range. We have performed simulations at multiple Rayleigh numbers far from the onset of convection. The time-periodic states have been numerically identified and converged to periodic orbits via the standard Newton shooting approach. These converged time-periodic solutions were subsequently extended in Rayleigh number by parametric continuation. The connections between the periodic orbits and the previously discussed fixed points are shown in the bifurcation diagram in figure 8(a). (As stated previously, other unstable equilibria and periodic orbits exist that we did not investigate or include in figure 8.) This section is devoted to explaining this figure.
4.1. Period-doubling cascade: PO3–PO6
Periodic orbit PO3 arises from FP3 in a supercritical Hopf bifurcation at $Ra=11\,261$, at which all of the spatial symmetries of FP3 are preserved. This is depicted in the upper left inset of figure 8(a), where the two red branches bifurcating from the FP3 branch correspond to the maximum and minimum of $\| \theta \|_2$ along PO3. Orbit PO3 was observed by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015), and the threshold that they reported is $Ra=11\,270$. Comparing the vorticity isosurfaces of FP3 at $Ra = 11\,000$ and PO3 at $Ra = 11\,500$, Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) show in their figures 10 and 13 that PO3 conserves most of the spatial structure of FP3, with the addition of strands connecting the rolls that appear and disappear periodically.
Orbit PO3 undergoes a period-doubling bifurcation leading to PO4 at $Ra = 12\,066$, with further period-doubling bifurcations leading to PO5 and PO6 at $Ra=12\,257$ and 12 306 (very similar to the thresholds $Ra=12\,068$, 12 258 and 12 306 found by Gao et al. Reference Gao, Podvin, Sergent and Xin2015). The temperature norms $\| \theta \|_2$ of PO3–PO6 are all close, as is typical for period-doubling cascades. The periods of these limit cycles are shown in figure 8(b), where period-doubling bifurcation points are indicated by stars. We were able to continue all of PO3–PO6 until at least $Ra=13\,300$. The dynamics of PO4–PO6 are very similar to those of PO3, so we do not show visualizations of them. Orbits PO3–PO6 inherit the fourfold symmetry $[D_4]_{xz}$ of FP3, hence their spatial and temporal variations all take place within a single roll. These states and the transitions between them are the only phenomena that we report that are not related to competition between three and four rolls.
We refer readers to Gao et al. (Reference Gao, Podvin, Sergent and Xin2015) for their measurements of the convergence to the Feigenbaum number characterizing the accumulation of period-doubling bifurcations until chaos. Indeed, Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) observed that the flow in $[L_x, L_y, L_z] = [1, 1, 10]$ becomes irregular following subharmonic oscillations at $Ra=12\,200$. Gao et al. (Reference Gao, Podvin, Sergent and Xin2015) estimated that a chaotic regime is reached at approximately $Ra = 12\,320$ (in $[L_x, L_y, L_z] = [1, 1, 2.5]$), and our DNS results confirm this. Gao et al. (Reference Gao, Podvin, Sergent and Xin2015) also discovered and reported five other periodic windows in their table II, each corresponding to another period-doubling cascade leading to chaotic behaviour. Although we have converged and continued these periodic windows, we omit them from the bifurcation diagram to avoid its becoming even more crowded.
Domain size has a major effect on the stability of PO3–PO6. When computed in a domain of the size of one wavelength $[1, 1, 2.5]$, PO3 is stable until it is succeeded by PO4, but in the domain $[1, 1, 10]$, it becomes unstable at $Ra\approx 11\,700$ by undergoing a large-wavelength instability which breaks the fourfold symmetry. Since this is prior to the period-doubling bifurcation, PO4, PO5, PO6 and the subsequent states resulting from the period-doubling cascade are also unstable in $[1, 1, 10]$. The stability properties that we have chosen to indicate in figure 8(a) for PO3–PO6 are those for the domain $[1, 1, 2.5]$. We have summarized the range of existence and stability of PO3–PO6 in both three-dimensional domains in table 1. Because FP3 is a three-dimensional state, PO3–PO6 do not exist in the two-dimensional domain $[1, 0, 10]$.
4.2. Wavelength-changing periodic orbits: PO1 and PO2
For $Ra>12\,200$ in domain $[1, 1, 10]$, Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) observed that numerical simulations started from a random initial condition settled down to three rolls, while four rolls were found intermittently. They also observed hysteresis in a simulation at $Ra=12\,000$, initiated from an instantaneous three-roll flow field at $Ra = 12\,200$, which finally settled down to a periodic orbit with three rolls. We have carried out DNS at nearly the same parameter values, illustrated in figure 9(a), which show long intervals of two types of time-periodic behaviour, which we call PO1 and PO2. The simulation is started from a random initial condition. The initial chaotic behaviour for $t \lesssim 250$ is succeeded by the weakly unstable PO1 ($300 \lesssim t \lesssim 650$). Afterwards, PO2 is visited ($1000 \lesssim t \lesssim 1600$) before a transition via subcritical period-doubling to a fully chaotic state ($t \gtrsim 1800$), which then persists without relaminarization. A phase portrait of these DNS is shown in figure 9(b), where instantaneous flow fields are represented by grey dots in the $D/I\unicode{x2013}I$ plane, where $D$ is the dissipation due to viscosity, and $I$ is the thermal energy input due to buoyancy. It can be seen that the two periodic orbits, with relatively low input/dissipation, are surrounded by the fully chaotic dynamics in this projection.
Continuing PO1 backwards in Rayleigh number, we find that PO1 is created via a supercritical reflection-symmetry-breaking Hopf bifurcation from FP4 at $Ra=9980$, where FP4 is stable, so that PO1 is stable at onset. The complex conjugate neutral eigenvector pair of FP4 is anti-$xz$-reflection symmetric, so PO1 has the spatial symmetry group $S_{PO_1} \equiv \langle {{\rm \pi} _y, \tau (\Delta y,0)} \rangle \sim [O(2)]_y$, and the spatio-temporal symmetry
where here, $z_0 \approx 4$ and $T/2 \approx 14$; compare figures 9(d) and 9(e), as well as figure 10(a).
Orbit PO1 arises from FP4 at a Rayleigh number at which it has three unequal rolls; its temporal dynamics consists of periodic lengthening and near-fragmentation of the longest roll. Figure 10 shows temperature profiles in $z$ for PO1 (and PO2) at fixed $x$ and $y$ locations, and at two instants separated in time by approximately a half-period. In these profiles, two successive zero-crossings ($\theta _{local}(z)=0$) correspond to one roll. Figure 10(a) shows a close, but unsuccessful approach to roll creation in PO1 at $z\approx 3$ ($t=12$) and then again at $z\approx 5$ ($t=26$).
Orbit PO1 loses stability in a supercritical pitchfork bifurcation at $Ra=12\,013$ (emphasized in the upper right inset of figure 8a), leading to the creation of PO2, in which $y$-translation symmetry is broken, but $y$-reflection symmetry is retained: $S_{PO_2} \equiv \langle {{\rm \pi} _y} \rangle \sim [Z_2]_y$. Unlike PO1, in PO2 the central roll succeeds in splitting, as can be seen in figures 11(b) and 11(h), as well as by the zero-crossings in figure 10(b), but only temporarily (figures 11(b,e,h,k). Orbit PO2 has the spatio-temporal symmetry
where $z_0 \approx 4$ and $T/2 \approx 17$. (The wavy structure of PO2 leads to the requirement that a half-domain shift in $y$ be included in (4.2). The combination of $\tau (0.5,0)$ and ${\rm \pi} _{xz}$ relates the red upward-facing ‘tongues’ centred at $y=0.5$ to the blue downward-facing ‘tongues’ centred at $y=0$.)
Orbit PO2 undergoes a saddle–node bifurcation at $Ra=12\,832$, and we followed this PO2 branch until $Ra=12\,804$, shortly after the saddle–node bifurcation. Also, PO2 undergoes a secondary Hopf bifurcation (also called a Neimark–Sacker or torus bifurcation) at $Ra=12\,082$, at which a complex conjugate pair of eigenvalues crosses the real axis, which generally leads to a torus that we will not discuss in the present work. Thus PO2 is stable only over a very small range of Rayleigh number ($12\,013 < Ra < 12\,082$), as can be seen in the upper right inset of figure 8(a). It is therefore not surprising that PO2 was not observed by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013). The period of PO1 remains nearly constant over a wide range of $Ra$, while that of PO2 increases with $Ra$ and then stays almost constant.
Orbits PO1 and PO2 capture the oscillatory dynamics of convection rolls. Orbit PO1 has three non-uniform rolls of fluctuating size and intensity (figures 9c–e). Although some rolls stretch and become quite weak, they never actually split anywhere along the branch. For PO2, the number of rolls varies between three (figures 11(a,c,d,f) along with figures 11(g,i,j,l)) and four (figures 11(b,e) along with figures 11(h,k)). We suggest that the intermittency and hysteresis observed by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) are a manifestation of visits to the coexisting unstable periodic orbits PO1 and PO2.
Since PO1 is $y$-independent, it can also exist in domain $[L_x, L_y, L_z] = [1, 0, 10]$. While the existence ranges of PO1 in domains $[1, 0, 10]$ and $[1, 1, 10]$ are the same, their stabilities differ: in $[1, 0, 10]$, the bifurcation to $y$-dependent PO2 does not occur, so PO1 remains linearly stable at least until $Ra=14\,000$, the upper limit at which we stopped the continuation. Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013) also observed this oscillating three-roll flow in the two-dimensional domain $[1, 0, 10]$, for $13\,500< Ra<15\,300$. Orbits PO1 and PO2 cannot exist in domain $[1,1,2.5]$, and PO2 cannot exist in domain $[1,0,10]$. The existence and stability intervals that we have computed for all of these flows are stated and compared with those of Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015) in table 1.
5. Discussion and conclusions
Vertical convection supports a large variety of flow patterns and thus can serve as a paradigm for nonlinear pattern formation in driven dissipative out-of-equilibrium systems. In this work, we have investigated thermal convection between two vertical plates held at different temperatures via both numerical simulation and continuation. We have computed the stable and unstable invariant solutions of the fully nonlinear three-dimensional Oberbeck–Boussinesq equations leading to the spatio-temporal complex convection patterns observed in experiments and simulations, far beyond the onset of convection.
We have also discovered previously unknown fixed points and periodic orbits. For these and the previously known solutions, we have identified the bifurcations responsible for their generation and termination, as well as their stabilization and destabilization. Summarizing the bifurcations and regimes in the computational domains $[L_x, L_y, L_z] = [1, 1, 10]$, $[1, 1, 2.5]$ and $[1, 0, 10]$, we compare in table 1 the results from Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013, Reference Gao, Podvin, Sergent and Xin2015) with those that we have obtained by computing unstable states and their bifurcations. Good quantitative agreement is achieved for those states that are observed in both studies. We note that all of the solution branches that we have found are connected by one or more bifurcations to the laminar branch.
Despite the complexity of the bifurcation diagrams shown in figures 2 and 8, almost all of the dynamics of this system results from one simple physical phenomenon: the competition between three and four rolls. (The exception is the steady three-dimensional state FP3 and the subsequent period-doubling sequence PO3–PO6, whose dynamics involves a single roll.) This competition takes different forms for the steady and time-dependent states.
For steady states, the most basic scenario in pattern formation is essentially one-dimensional and consists of a sequence of primary bifurcations from a featureless state to branches with different numbers of regularly spaced rolls, cells or waves. Of these, only the first to bifurcate is stable at onset; primary branches often undergo secondary instabilities to mixed-mode branches that transfer stability between the different branches (e.g. Knobloch & Guckenheimer Reference Knobloch and Guckenheimer1983). In an idealized context, this is the Eckhaus instability (e.g. Eckhaus Reference Eckhaus1965; Tuckerman & Barkley Reference Tuckerman and Barkley1990). A mathematical framework that covers wavelength competition is that of mode interaction (Dangelmayr Reference Dangelmayr1986).
Our hydrodynamic configuration involving the four-roll branch FP1, the three-roll branch FP2, and the connector branch FP4 presents a complicated version of this basic scenario. The FP1 branch with four equally spaced rolls has $D_4$ symmetry and hence gives rise to two sets of mixed-mode branches FP4. These two sets are associated with two qualitatively different paths for passing from four to three rolls: the merging of two rolls, and the disappearance of a roll. (This is kinematically possible because these are co-rotating rolls, rather than the counter-rotating rolls of the standard Rayleigh–Bénard configuration.) Indeed, dual sets of branches are a typical feature of bifurcation in the presence of $D_4$ symmetry (Swift Reference Swift1985; Knobloch Reference Knobloch1986). The $D_4$ bifurcation scenario is present in many other situations and will be encountered again in our companion paper Zheng et al. (Reference Zheng, Tuckerman and Schneider2024), in the more classic two-dimensional context of competition between straight and diagonal orientations (Demay & Iooss Reference Demay and Iooss1984; Tagg et al. Reference Tagg, Edwards, Swinney and Marcus1989; Chossat & Iooss Reference Chossat and Iooss1994; Bengana & Tuckerman Reference Bengana and Tuckerman2019; Reetz et al. Reference Reetz, Subramanian and Schneider2020).
A new feature of the $D_4$ scenario seen here is that each FP4 half-branch of disappearing-roll type meets and merges smoothly with an FP4 half-branch of merging-roll type. To the best of our knowledge, this phenomenon has not been observed previously. The FP4 branches that merge do not emanate from the same four-roll branch, but from two branches, FP1 and FP1$^\prime$, that are phase-shifted by a half-roll with respect to one another. The simultaneous existence of branches FP1 and FP1$^\prime$ is in turn due to the fact that FP1 branches of all phases are created by a circle pitchfork bifurcation. At this meeting point, the two types of half-branches also meet the FP2 branch, which contains states with three equal rolls, in a transcritical bifurcation. The phase jump can instead be assigned to the transcritical bifurcation, i.e. between three-roll branches FP2 and FP2$^\prime$, rather than to the pitchfork bifurcations from the four-roll branches. A last alternative is to follow the FP4 branch twice, via FP1, FP2, FP1$^\prime$, FP2$^\prime$, FP1 without any phase jumps. The $D_3$ symmetry of the FP2 states governs the details of the transcritical bifurcation. The physical phenomena of roll merging and roll disappearance, roll creation and roll splitting, provide visual illustrations of the group-theoretic $D_4$ and $D_3$ scenarios, and of the $D_3$–$D_4$ mode interaction equations in Dangelmayr (Reference Dangelmayr1986) and Crawford et al. (Reference Crawford, Knobloch and Riecke1990).
Turning now to time-dependent solutions, periodic orbits PO1 and PO2 could be considered to be temporal versions of the variation with Rayleigh number along the connector branch FP4, from which PO1 bifurcates. Figure 9 shows that PO1 contains temporal phases in which one of its three rolls widens or weakens, resembling the precursors to four rolls seen in figure 3 as $Ra$ is varied along the FP4 branch. Although PO1 does not succeed in creating a fourth roll, figure 11 shows that in PO2, these events culminate in the periodic formation and destruction of a fourth small and fragile roll. This may or may not be related to the breaking of $y$-translation symmetry from PO1 to PO2; perhaps roll formation or destruction is facilitated when rolls become wavy. The competition between three and four rolls continues to dominate for Rayleigh numbers above 12 082 when PO2 is destabilized, since the dynamics beyond this point consists of chaotic three-roll flow with infrequent and irregular bursts of four rolls, as illustrated by Gao et al. (Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013). In the present study in domain $[L_x, L_y, L_z] = [1, 1, 10]$, we report no stable fixed points for $Ra>11\,261$ and no stable periodic orbits for $Ra>12\,082$. Nevertheless, we have been able to numerically continue these unstable solutions into the chaotic regime, far from the parameter regime in which they are stable.
Although DNS can give access to complex flow dynamics at specific control parameters, numerical continuation organizes these solutions and determines their bifurcation-theoretic origin, by situating them in the context of a bifurcation diagram. Our work bridges the gap between purely DNS-based observations and numerical bifurcation analysis, leading to a better description and understanding of complex convective flows.
Acknowledgements
We thank S. Azimi, F. Reetz and O. Ashtari for fruitful exchanges. We thank P. Le Quéré for sharing his encyclopedic knowledge of vertical convection. We are grateful to D. Barkley, J. Dawes, E. Knobloch and A. Rucklidge for their insights on symmetry. We are grateful to the anonymous referee who added the alternative construction of the connector branch, and to all of the referees for their very valuable comments.
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.