2.1 Doubly diffusive convection

Following Beaume et al. (Reference Beaume, Bergeon and Knobloch2013a ) we consider doubly diffusive convection in a binary fluid placed within a three-dimensional box of square cross-section in the horizontal and large extent in the vertical. We use the coordinate $x$ for the vertical direction and $(y,z)$ for the horizontal coordinates. Convection is produced by competing but balanced horizontal gradients of temperature and concentration that result from the imposition of appropriate Dirichlet boundary conditions at $z=0,l$ leading to the conduction state $\boldsymbol{u}^{\ast }=0$ , $T^{\ast }=T_{r}+z\unicode[STIX]{x0394}T/l$ , $C^{\ast }=C_{r}+z\unicode[STIX]{x0394}C/l$ , where $\boldsymbol{u}^{\ast }$ is the velocity, $T^{\ast }$ is the temperature, $C^{\ast }$ is the concentration of the heavier component of the binary mixture (hereafter shortened to concentration), $T_{r}$ and $C_{r}$ are reference temperature and concentration and $\unicode[STIX]{x0394}T>0$ and $\unicode[STIX]{x0394}C>0$ are the imposed temperature and concentration differences across the box. We non-dimensionalize the governing equations using $l$ for lengths, $\unicode[STIX]{x0394}T$ for the temperature, $\unicode[STIX]{x0394}C$ for the concentration and $l^{2}/\unicode[STIX]{x1D705}$ for the time, where $\unicode[STIX]{x1D705}$ is the thermal diffusivity. In the absence of cross-diffusion, the Boussinesq equations describing the system are

Here $\boldsymbol{u}=(u,v,w)$ is the dimensionless velocity field in $(x,y,z)$ coordinates, $p$ is the dimensionless pressure, $T\equiv (T^{\ast }-T_{r})/\unicode[STIX]{x0394}T$ is the dimensionless temperature and $C\equiv (C^{\ast }-C_{r})/\unicode[STIX]{x0394}C$ is the dimensionless concentration. The above equations involve three dimensionless numbers: the Rayleigh number $Ra\equiv g|\unicode[STIX]{x1D70C}_{T}|\unicode[STIX]{x0394}Tl^{3}/\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}$ , the Prandtl number $Pr\equiv \unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ and the (inverse) Lewis number $\unicode[STIX]{x1D70F}\equiv D/\unicode[STIX]{x1D705}$ , where $g$ is the gravitational acceleration, $\unicode[STIX]{x1D70C}_{T}/\unicode[STIX]{x1D70C}_{0}<0$ is the thermal expansion coefficient in the Boussinesq approximation evaluated at $T_{r}$ , $\unicode[STIX]{x1D70C}_{0}$ is the density of the fluid at $T_{r}$ and $\unicode[STIX]{x1D708}$ and $D$ are, respectively, the kinematic viscosity and the molecular diffusivity of the heavier fluid component, both assumed to be constant. In writing this set of equations, we have assumed that the buoyancy ratio $N\equiv \unicode[STIX]{x1D70C}_{C}\unicode[STIX]{x0394}C/\unicode[STIX]{x1D70C}_{T}\unicode[STIX]{x0394}T=-1$ , where $\unicode[STIX]{x1D70C}_{C}/\unicode[STIX]{x1D70C}_{0}>0$ is the solutal expansion coefficient evaluated at $T_{r}$ , $C_{r}$ . With this choice of $N$ the system (2.1)–(2.4) possesses the stationary conduction solution alluded to above.

Figure 1.

Sketch of the physical domain from an angle similar to that used in subsequent figures. For convenience the system has been rotated anticlockwise by $90^{\circ }$ around the $y$ axis, i.e. clockwise on the page. As a result, the acceleration due to gravity points towards the left. Shaded boundaries are maintained at fixed temperature and concentration with the darker wall maintained at higher values of $T$ and $C$ than the lighter one. Only a portion of the domain in $x$ is shown.

The physical domain is a closed container of square horizontal cross-section and aspect ratio $L\gg 1$ in the vertical direction, as shown in figure 1. We prescribe no-slip boundary conditions for the velocity everywhere and no-flux boundary conditions for the temperature and concentration on all walls except $z=0,1$ :

The system (2.1)–(2.4), together with the boundary conditions (2.5)–(2.8), admits the conduction solution $(u,v,w,T,C)=(0,0,0,z,z)$ for all values of the parameters $Ra$ , $Pr$ and $\unicode[STIX]{x1D70F}$ . In terms of the perturbation quantities $(u,v,w,\unicode[STIX]{x1D6E9},\unicode[STIX]{x1D6F4})$ with $\unicode[STIX]{x1D6E9}\equiv T-z$ and $\unicode[STIX]{x1D6F4}\equiv C-z$ , the system is equivariant with respect to the dihedral group $D_{2}$ generated by the two reflections

The group structure implies that the equations are also equivariant with respect to $S_{c}=S_{y}\circ S_{\unicode[STIX]{x1D6E5}}$ :

These symmetries play a key role in the properties of the solutions described below. In particular, the symmetry $S_{\unicode[STIX]{x1D6E5}}$ plays a role that is analogous to the role played by spatial reflection symmetry in the Swift–Hohenberg equation.

When $N\neq -1$ the properties of the system change dramatically. In particular, there is no motionless conduction state and the base state now consists of a single large scale roll that spans the whole domain $0\leqslant x\leqslant L$ . This circulation is responsible for changing the primary bifurcation to an imperfect bifurcation for values of $N$ close to $-1$ . Additional effects can be identified for $N$ farther away from $-1$ . In this paper we focus on the case $N=-1$ only.

We solve the system (2.1)–(2.4) numerically using a splitting scheme for time evolution and a finely tuned numerical continuation algorithm (Beaume Reference Beaume2017). Throughout this paper, the parameters are chosen following Xin, Le Quéré & Tuckerman (Reference Xin, Le Quéré and Tuckerman1998), Bergeon & Knobloch (Reference Bergeon and Knobloch2002): $Pr=1$ , $\unicode[STIX]{x1D70F}=1/11$ and $L=19.8536$ , corresponding to eight wavelengths of the primary instability of the two-dimensional problem. The Rayleigh number $Ra$ is kept as a control parameter.

2.2 State of the art

Spatially localized solutions of doubly diffusive convection in a vertically extended domain have been studied for over 10 years. Their origin was first investigated in a study of near-onset behaviour in a two-dimensional vertically periodic slot of relatively small spatial period (Bergeon & Knobloch Reference Bergeon and Knobloch2008a ). These authors found a number of spatially periodic states yielding a complex bifurcation diagram and identified states composed of a single roll occupying a part of the domain. This study was followed up by a study describing the evolution of spatially localized doubly diffusive convectons in the same configuration but with a large spatial period in the vertical direction (Bergeon & Knobloch Reference Bergeon and Knobloch2008b ). The authors showed that the primary bifurcation from the conduction state is subcritical. The resulting spatially periodic solution is $S_{\unicode[STIX]{x1D6E5}}$ -symmetric and transforms rapidly into an array of corotating rolls as one proceeds away from the primary bifurcation, before turning towards larger values of $Ra$ at a fold and generating an array of large-amplitude corotating rolls. During this process the solution remains $S_{\unicode[STIX]{x1D6E5}}$ -symmetric. Subsequent branches of spatially periodic states with different wavenumbers but also $S_{\unicode[STIX]{x1D6E5}}$ symmetry behave in the same way but may acquire stability at large amplitude. As expected from earlier studies of the Swift–Hohenberg equation, the primary subcritical branch of spatially periodic states becomes unstable almost immediately to spatially modulated states – still $S_{\unicode[STIX]{x1D6E5}}$ -symmetric – and these evolve rapidly with decreasing $Ra$ into strongly localized states exhibiting homoclinic snaking (Woods & Champneys Reference Woods and Champneys1999; Burke & Knobloch Reference Burke and Knobloch2006). As in the Swift–Hohenberg equation there are two branches of localized states, one consisting of an odd number of corotating rolls all rotating anticlockwise (clockwise on the page, see figure 1) about the $y$ axis and the other of an even number of corotating rolls also rotating anticlockwise. In the following we refer to these states as $L^{+}$ and $L^{-}$ , respectively, and emphasize that each is $S_{\unicode[STIX]{x1D6E5}}$ -symmetric with respect to the centre of the domain. In contrast to the Swift–Hohenberg equation, however, the doubly diffusive problem is not gradient-like and therefore all asymmetric solutions on the rungs of the snakes-and-ladders structure of the bifurcation diagram correspond to states that drift in the $x$ direction. States of this type were not computed in Bergeon & Knobloch (Reference Bergeon and Knobloch2008a ,Reference Bergeon and Knobloch b ) as they are expected to be unstable, and are likewise not computed here. Such drifting states were computed in related systems, however, confirming the basic bifurcation scenario (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2013; Lo Jacono, Bergeon & Knobloch Reference Lo Jacono, Bergeon and Knobloch2017).

The above properties of the system persist in two-dimensional domains with end walls at $x=0,L$ . However, such walls destroy the translation invariance of the system in the $x$ direction and prevent the formation of spatially periodic solutions (Beaume et al. Reference Beaume, Bergeon and Knobloch2013a ). As a result localized states bifurcate directly from the conduction state (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2009) rather than appearing in a secondary bifurcation from a periodic state. This behaviour persists in three-dimensional domains with no-slip walls (Sezai & Mohamad Reference Sezai and Mohamad2000; Bergeon & Knobloch Reference Bergeon and Knobloch2002). In this case all the primary branches have an additional symmetry, the reflection symmetry $S_{y}$ . This symmetry exerts a substantial influence on the properties of the resulting solutions since it can be broken in secondary bifurcations, a fact responsible for the formation of twisted localized states (Beaume et al. Reference Beaume, Bergeon and Knobloch2013a ), as discussed further below.

The bifurcation diagram for steady doubly diffusive convectons in a three-dimensional domain of square cross-section and aspect ratio $L=19.8536$ is shown in figure 2. The solutions are represented in terms of the total kinetic energy:

The figure shows that the convectons bifurcate directly from the conduction state in a transcritical bifurcation at $Ra\approx 850.86$ . This bifurcation gives rise to two branches of convectons, $L^{+}$ and $L^{-}$ , both of which snake within the interval $703\lesssim Ra\lesssim 807$ . Both types of convectons are fully symmetric, i.e. invariant under $S_{y}$ , $S_{\unicode[STIX]{x1D6E5}}$ and, by extension, $S_{c}$ . Convectons from the branch $L^{+}$ possess an odd number of corotating, anticlockwise rolls centred in the domain while those from the branch $L^{-}$ possess an even number of corotating, anticlockwise rolls also centred (see figure 3). These states differ from those found in a two-dimensional periodic domain in minor ways only: (i) because the boundary conditions at $x=0,L$ are no longer periodic, no spatially periodic state is present and the convectons therefore cannot and do not connect to any spatially periodic state, and (ii) owing to the no-slip boundary condition in $y$ , the convecton structure is fully three-dimensional. Despite these differences, the snaking in the two-dimensional and three-dimensional systems is similar: in the direction of increasing energy, going from a left fold to the next right fold (hereafter called a positive segment, alluding to the fact that $Ra$ increases), the corotating rolls strengthen, while going from a right fold to the next left fold (hereafter called a negative segment: $Ra$ decreases) results in the nucleation of one new roll on either side of the existing roll structure.

Figure 2.

Bifurcation diagram representing the total kinetic energy $E$ as a function of the Rayleigh number $Ra$ . The solution $E=0$ corresponds to the conduction state and convectons ( $E>0$ ) emerge through a sequence of bifurcations from the conduction state at $Ra\approx 850.86$ . The primary branches $L^{\pm }$ consist of corotating rolls with axes parallel to the $y$ axis. Solutions on the secondary branches $L_{1}^{+}$ , $L_{2}^{+}$ and $L_{1}^{-}$ , $L_{2}^{-}$ bifurcate from $L^{\pm }$ and break the symmetry $S_{y}$ . Such states are referred to as twisted.

Figure 3.

Convectons with the lowest kinetic energy along the branch $L^{+}$ at $Ra=710.4464$ (a) and $L^{-}$ at $Ra=703.9686$ (b) represented in terms of the isocontours $u=\pm 0.3$ . Positive values are shown in red. In this and subsequent figures $x$ points to the right (gravitational acceleration is to the left) and $z$ points upward.

Secondary bifurcations occur along the negative segments of these snaking branches. These symmetry-breaking bifurcations all break the $S_{y}$ symmetry together with one other reflection and yield roll structures we refer to as twisted. In these states the axes of adjacent rolls are rotated in opposite directions around the vertical axis. One such pitchfork bifurcation occurs along the first negative segment of $L^{+}$ (counting from low to high energy) and yields the branch $L_{1}^{+}$ of twisted convectons, while two bifurcations occur in very close succession along the first negative segment of $L^{-}$ . One of these yields the branch $L_{1}^{-}$ shown in figure 2 but we have not been able to continue the states created in the other. The states on the $L_{1}^{+}$ branch near its birth consist of a single $S_{c}$ -symmetric roll, i.e. a roll whose axis is rotated about the vertical, while the $L_{1}^{-}$ branch consists of two-roll $S_{\unicode[STIX]{x1D6E5}}$ -symmetric states rotated in opposite directions about this axis. These states are shown in figure 4(a,b).

Figure 4.

Convectons with the lowest kinetic energy along the branch $L_{1}^{+}$ at $Ra=817.9884$ (a), $L_{1}^{-}$ at $Ra=818.0300$ (b), $L_{2c}^{+}$ at $Ra=820.0688$ (c), $L_{2\unicode[STIX]{x1D6E5}}^{+}$ at $Ra=820.1873$ (d), $L_{2c}^{-}$ at $Ra=823.9576$ (e) and $L_{2\unicode[STIX]{x1D6E5}}^{-}$ at $Ra=823.4842$ (f) represented in the same way as in figure 3.

Successive secondary bifurcations yielding both $S_{c}$ - and $S_{\unicode[STIX]{x1D6E5}}$ -symmetric twisted convectons can be found along each subsequent negative segment of each primary branch. We refer to the $S_{c}$ -symmetric states bifurcating from the second negative segment of the $L^{+}$ branch as $L_{2c}^{+}$ and to the $S_{\unicode[STIX]{x1D6E5}}$ -symmetric states bifurcating from this segment as $L_{2\unicode[STIX]{x1D6E5}}^{+}$ . These branches are indistinguishable in the bifurcation diagram in figure 2 and are therefore labelled $L_{2}^{+}$ for simplicity. Both consist of solutions composed of one broad untwisted roll inherited from the $L^{+}$ branch, surrounded by two narrower rolls of smaller amplitude and rotated about the $x$ axis, as shown in figure 4(c,d). This figure also shows solutions $L_{2c}^{-}$ and $L_{2\unicode[STIX]{x1D6E5}}^{-}$ produced in successive bifurcations along the second negative segment of $L^{-}$ (figure 4 e,f). These secondary branches produce secondary snaking in $745\lesssim Ra\lesssim 819$ but the mechanism differs from that responsible for the primary snaking. Going from a left fold to the next right fold in the direction of increasing energy (positive segment), the dominant effect is a rotation of the roll axes about the $x$ axis. Past the right fold, each structure gradually nucleates a pair of new rolls, one on either side, with rotation that is always opposite to that of the roll adjacent to it.

Having recapitulated existing results we now present three examples of the evolution of the system initialized with the three-dimensional states just described. In figure 5(a) we present the temporal evolution of a quasi-two-dimensional localized $L^{-}$ convecton, followed in figure 5(b) by that of an initial condition taken from the same branch but that is domain-filling and in figure 5(c) by that of a secondary twisted convecton from the $L_{1}^{+}$ branch. The time evolution is represented in terms of the convection amplitude quantified through the quantity

This measure of convection strength is appropriate since $u$ is the dominant velocity component. The interesting and remarkable fact is that all our initial conditions ultimately collapse to the conduction state, whether they are very localized, domain-filling or twisted. This observation suggests that all the non-trivial states in the subcritical regime of this problem are unstable. The main goal of this paper is to explain this remarkable fact and explore its consequences.

Figure 5.

Space–time evolution of three subcritical states: (a)  $L^{-}$ convecton obtained at $Ra\approx 747.14104$ and time integrated at $Ra=747.1410$ , (b) $L^{-}$ convecton obtained at $Ra\approx 756.04067$ and time integrated at $Ra=756.0407$ and (c) $L_{1}^{+}$ convecton obtained at $Ra\approx 775.76716$ and time integrated at $Ra=775.7672$ .

3.1 Stability of the primary convectons

Our eigenvalue computations show that all unstable modes of the convecton branches $L^{\pm }$ are associated with monotonic growth: $\unicode[STIX]{x1D706}_{i}=0$ . Figure 8 shows the resulting growth rates $\unicode[STIX]{x1D706}_{r}$ along both branches as a function of the arclength $s$ . The $L^{+}$ branch is subcritical and is thus characterized by one unstable amplitude mode in addition to the unstable eigendirection inherited from the primary pitchfork (see table 1). The branch becomes three times unstable at $s\approx 100$ ( $Ra\approx 750$ ) owing to the onset of the twist instability responsible for the secondary branch $L_{1}^{+}$ . The twist instability is triggered by the eigenmode shown in figure 9. The eigenmode represents a perturbation roll that is rotated around the $x$ axis by $90^{\circ }$ as compared to that in the base state, resulting in the progressive rotation of the convection roll around the $x$ axis as the instability develops. At $s\approx 139$ , the branch undergoes a saddle-node bifurcation where the mode shown in figure 10 becomes stable. This mode has a similar profile to that of the base state and is thus an amplitude mode: any instability arising from this mode either strengthens or weakens the base convection roll. This eigenmode remains stable along the whole subsequent snaking branch as it acts only on the central roll. The next two saddle nodes are caused by the destabilization and restabilization of a different eigenmode, shown in figure 11(b). The figure shows that this mode is responsible for the addition of one convection roll on either side of the localized structure, a process that takes place just beyond this fold ( $s\approx 236$ ). This eigenmode becomes stabilizing at the next saddle node ( $s\approx 337$ ) and remains so thereafter since it always determines the stability of the same two rolls. Consequently, each subsequent pair of saddle nodes is generated by a new and different eigenmode. This is so for the phase mode (figure 11 c) as well, and is in contrast to the Swift–Hohenberg equation where the same two amplitude and phase eigenvalues repeatedly pass through zero as one follows the snaking branches (Burke & Knobloch Reference Burke and Knobloch2006; Kao, Beaume & Knobloch Reference Kao, Beaume and Knobloch2014; Knobloch Reference Knobloch2015). The saddle node at $s\approx 236$ is followed shortly thereafter by a symmetry-breaking bifurcation triggered by a phase mode, breaking the $S_{\unicode[STIX]{x1D6E5}}$ and $S_{c}$ symmetries but preserving the $S_{y}$ symmetry. This bifurcation creates a corotating roll on one side of the structure and a counter-rotating roll on the other side. Since $\unicode[STIX]{x1D70F}<1$ , the growth of the latter is suppressed by nonlinear terms (Thangam, Zebib & Chen Reference Thangam, Zebib and Chen1982) resulting in the creation of a branch of asymmetric rung states. Owing to the non-gradient nature of the present system this asymmetric state will necessarily drift (in the $x$ direction) until such time as it comes into contact with one of the end walls and comes to rest (see below). The next change in stability is located farther away on the negative segment at $s\approx 300$ ( $Ra\approx 740$ , see table 1) where two successive bifurcations take place in rapid succession. These bifurcations are of the same nature: they are both produced by eigenmodes breaking the symmetry $S_{y}$ and one other symmetry, and result in rotation of the axes of the end rolls either in the same or opposite direction as shown in figure 12. Since these modes act on the end rolls, their regions of influence are separated by the central roll where the eigenmode vanishes. As a result, there is little dynamical difference between them and the associated bifurcations occur at similar values of $s$ , equivalently $Ra$ . These twist modes destabilize the solution yet further. However, near the next left saddle node (at $s\approx 335$ ) the amplitude and phase mode stabilize again, albeit at slightly different values of $Ra$ . This scenario then repeats during each and every oscillation of the $L^{+}$ branch. Overall, the $L^{+}$ solution gains two unstable eigendirections during each such oscillation, as summarized in table 1.

Figure 11.

(a) $L^{+}$ convecton at $Ra\approx 802.1785$ ( $s\approx 239$ ) represented using the same colour scheme as in figure 3. (b) Amplitude and (c) phase eigenmodes for the state in (a) represented in the same way as in figure 9(b). The former is responsible for the second and third saddle nodes while the latter is responsible for the emergence of drifting rung states via parity-breaking bifurcations.

Figure 12.

(a) $L^{+}$ convecton at $Ra\approx 738.1785$ ( $s\approx 303$ ) represented using the same colour scheme as in figure 3. (b) $S_{\unicode[STIX]{x1D6E5}}$ -preserving twist eigenmode for the state in (a). (c) $S_{c}$ -preserving twist eigenmode for the state in (a). Both eigenmodes are represented in the same way as in figure 9(b).

The branch $L^{-}$ emerges supercritically from the transcritical bifurcation that produces $L^{+}$ but is initially once unstable, owing to the pitchfork bifurcation that precedes it. The branch undergoes a saddle node at low amplitude as shown in figure 7 where the solution gains a second unstable eigendirection and thereafter undergoes a similar sequence of bifurcations as $L^{+}$ , as summarized in table 1. The eigenvalue count is initialized at the small-amplitude saddle node. At $s\approx 12$ ( $Ra\approx 838$ ), the solution gains a third unstable direction with the corresponding eigenmode shown in figure 13. The figure shows that this mode strengthens one of the end rolls comprising the solution at the expense of the other. This instability likely corresponds to the termination of the branch created in the primary pitchfork bifurcation (Bergeon & Knobloch Reference Bergeon and Knobloch2002). Proceeding up the negative branch segment, one next encounters a pair of near-simultaneous bifurcations at $s\approx 108$ ( $Ra\approx 742$ ) responsible for generating the twisted states labelled $L_{1}^{-}$ following similar notation for $L_{1}^{+}$ (figure 14). After this original set of bifurcations, the $L^{-}$ branch follows a similar behaviour to $L^{+}$ , gaining two unstable eigendirections during each back and forth oscillation of the branch in parameter space: the phase mode stabilizes just below the left fold followed by the amplitude mode before destabilizing in reverse order at the right fold. A pair of twist instabilities occurs along each negative segment of the branch leading to the net generation of two unstable eigendirections per back and forth oscillation.

Figure 13.

(a) $L^{-}$ convecton at $Ra\approx 837.8913$ ( $s\approx 12$ ) represented using the same colour scheme as in figure 3 but for isocontours $u=\pm 0.1$ . (b) Eigenmode for the state in (a) represented in the same way as in figure 9(b).

Figure 14.

(a) $L^{-}$ convecton at $Ra\approx 741.8996$ ( $s\approx 108$ ) represented using the same colour scheme as in figure 3. (b) $S_{\unicode[STIX]{x1D6E5}}$ -preserving twist eigenmode for the state in (a). (c) $S_{c}$ -preserving twist eigenmode for the state in (a). Both eigenmodes are represented in the same way as in figure 9(b).

It will have been observed that both $L^{+}$ and $L^{-}$ possess a weakly unstable eigenvalue, with growth rate around $\unicode[STIX]{x1D706}_{r}\approx 10^{-4}$ for a one-roll (strongly localized) state and $\unicode[STIX]{x1D706}_{r}\approx 10^{-3}$ for a six-roll (domain-filling) state. The corresponding eigenfunction for a $L^{+}$ solution consisting of three convection rolls is shown in figure 15(b) and is of odd parity. Here and in the following we use this term to refer to eigenfunctions that are odd with respect to $S_{\unicode[STIX]{x1D6E5}}$ , the symmetry of $L^{\pm }$ , but invariant under $S_{y}$ . In a vertically unbounded system bifurcations triggered by the destabilization of an odd parity mode are associated with the onset of drift and are sometimes called drift bifurcations.

Figure 15.

(a) $L^{+}$ convecton at $Ra\approx 805.7820$ ( $s\approx 440$ ) represented using the same colour scheme as in figure 3. (b) Odd parity drift eigenmode for the state in (a) represented in the same way as in figure 9(b).

To unveil the consequences of the presence of this mode, we selected a solution that allows us to suppress all other unstable eigendirections by imposing appropriate symmetry conditions on a time simulation. For this purpose we selected a solution on the second positive segment of $L^{-}$ (counting the very small-amplitude segment shown in figure 7) at $Ra\approx 757.1410$ consisting of a pair of convection rolls. This state is shown in figure 16(a) and possesses three unstable eigendirections, two of which are associated with the twist instability and are even in $x$ (see figure 8) and one with odd parity in $x$ . Since the twist instability breaks the symmetry $S_{y}$ while the odd parity mode does not, we isolate the dynamics of the latter by imposing the symmetry $S_{y}$ in our time simulation. The effect of the odd parity mode is represented in figure 17 using the quantity $A_{conv}$ defined in expression (2.13). The figure shows that drift occurs on a time scale of the order of $10^{4}$ time units, a time scale that compares well with the magnitude of the odd parity mode eigenvalue. In this simulation, the convecton drifts downward until it meets the bottom wall. By symmetry, an upward drifting convecton can also be obtained. When the convecton meets the wall, it stops drifting and settles into a steady state in which the roll closest to the wall is stronger than that farther away. This wall-attached convecton has been converged using a Newton method (see figure 16 b) and is twice unstable, each unstable eigendirection being responsible for the twist of one of the rolls. We expect that these states also snake, much like similar states in binary fluid convection in a horizontal layer (Mercader et al. Reference Mercader, Batiste, Alonso and Knobloch2011).

Figure 16.

Initial condition (a) and final state (b) obtained through direct numerical simulation with imposed $S_{y}$ symmetry at $Ra=757.141$ ( $s\approx 200$ ). The states are represented by the isocontours $|u|=1.1$ .

Figure 17.

Space–time evolution at $Ra=757.141$ of the quantity $A_{conv}$ defined in expression (2.13) and initialized with the state in figure 16(a).

We can understand these results as follows. In an unbounded, translation-invariant system a growing odd parity perturbation would render an even state asymmetric. Since in such systems asymmetric states generically drift, the appearance of such an unstable mode is associated with a parity-breaking or drift bifurcation. In the present case, however, the walls at $x=0,L$ prevent translation, implying that the presence of a growing odd mode must lead to a stationary state. A simple model equation, studied by Knobloch, Hettel & Dangelmayr (Reference Knobloch, Hettel and Dangelmayr1995) and Dangelmayr, Hettel & Knobloch (Reference Dangelmayr, Hettel and Knobloch1997), captures the effect of broken translation invariance on the parity-breaking bifurcation. The model takes the form

where $\unicode[STIX]{x1D707}$ is the bifurcation parameter and $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D702}$ , $\unicode[STIX]{x1D716}$ are parameters proportional to the amplitude of the spatial inhomogeneity, assumed small and of period $2\unicode[STIX]{x03C0}$ . When these terms are absent equation (3.2) reduces to the normal form for the parity-breaking bifurcation, assumed here to occur at $\unicode[STIX]{x1D707}=0$ and to be supercritical, with $c$ the drift speed and $\unicode[STIX]{x1D719}$ the spatial phase or displacement. When $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D702}$ , $\unicode[STIX]{x1D716}$ are non-zero $c$ is no longer a drift speed but must be interpreted as the degree of asymmetry of the state. There are then generically two fixed points $(c,\unicode[STIX]{x1D719})=(0,0)$ and $(0,\unicode[STIX]{x03C0})$ , corresponding to a symmetric state at the ‘wall’ ( $\unicode[STIX]{x1D719}=0$ ) and a symmetric state at the centre of the domain ( $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}$ ), respectively. Other, asymmetric fixed points (with $c\neq 0$ ) may be present as well (Knobloch et al. Reference Knobloch, Hettel and Dangelmayr1995; Dangelmayr et al. Reference Dangelmayr, Hettel and Knobloch1997) but we focus here on the symmetric states. These have opposite stability properties implying that if $(0,0)$ is stable then $(0,\unicode[STIX]{x03C0})$ is unstable and vice versa. We identify the symmetry-breaking perturbation of a symmetric convecton at the centre as the perturbation that initially generates translation, but that is arrested by pinning arising from the wall either at $x=0$ or $x=L$ that brings this state to rest. In our case the final state that results is not symmetric but this is a consequence of the non-periodic boundary conditions employed in the simulation.

3.2 Stability of the twisted convectons

The stability properties along the lowest branches of twisted convectons, $L_{1}^{\pm }$ , are summarized in figure 18. The $S_{c}$ -symmetric twisted convectons on $L_{1}^{+}$ bifurcate towards higher Rayleigh numbers and are thus initially three times unstable: in addition to the unstable eigendirection associated with the twist instability and responsible for the formation of the branch, the $L_{1}^{+}$ convectons are also unstable with respect to an amplitude mode (responsible for the primary transcritical bifurcation) and a phase mode inherited from the odd parity drift mode of $L^{+}$ . Very close to the bifurcation point to $L_{1}^{+}$ the amplitude eigenvalue decreases while the twist eigenvalue increases resulting in an eigenvalue collision at $s\approx 4.7$ ( $Ra\approx 756.1$ ) forming a complex conjugate pair. Such a collision is possible because both modes have the same symmetry and is shown in the left panel of figure 19; the imaginary part of the eigenvalues along the $L_{1}^{+}$ branch is shown in figure 20. The resulting $(1,1)$ unstable $L_{1}^{+}$ convecton continues to larger values of $Ra$ with the unstable eigenmodes at $Ra\approx 781.9884$ shown in figure 21.

Figure 18.

Linear stability results for $L_{1}^{+}$ (a) and $L_{1}^{-}$ (b). The results are presented in terms of the quantity $e^{\unicode[STIX]{x1D706}_{r}}$ and shown as a function of the arclength $s$ along the branch. Folds correspond to saddle-node bifurcations ( $\unicode[STIX]{x1D706}_{r}=0$ ) and are denoted by vertical lines. Eigenvalue collisions are indicated by a red dot.

Figure 19.

Enlargement of the $s\in [1,6]$ interval (a) and of the $s\in [143.4,144]$ interval (b) for $L_{1}^{+}$ (figure 18 a). One eigenvalue collision occurs at $s\approx 4.7$ while two successive collisions are observed at $s\approx 143.7$ . After each of these collisions, the associated eigenvalues become complex conjugate eigenvalues.

Figure 20.

Imaginary part of the eigenvalues from the linear stability analysis of $L_{1}^{+}$ (a) and $L_{1}^{-}$ (b) as a function of the arclength $s$ along the branch. Saddle nodes are denoted by vertical lines.

Figure 21.

(a) $L_{1}^{+}$ convecton at $Ra\approx 781.9884$ ( $s\approx 31$ ) represented using the same colour scheme as in figure 3. All unstable eigenmodes are represented as in figure 9(b): (b, c) show the real and imaginary parts of the complex conjugate modes, while (d) shows the odd parity drift mode.

The odd parity drift mode of the twisted convecton shown in figure 21(a) is generally associated with slow dynamics while the observable fast dynamics is generated by unstable oscillatory modes (figure 21 b,c) that rotate the convection roll about the $x$ axis first in one direction and then back to its original state. To illustrate this behaviour, we initialized a simulation using an $L_{1}^{+}$ convecton at $Ra\approx 801.9884$ placed in the centre of a larger domain in order to reduce the influence of the walls. This convecton is three times unstable. The simulation is run at $Ra=801.9884$ and we let the instability grow from numerical error. Selected snapshots from the simulation are shown in figure 22. The snapshots show that the amplitude of the roll first grows to a maximum, its axis then rotates about the $x$ axis before its amplitude and rotation angle decrease again and the process repeats. This process can be seen in the transition from panel (b) at $t=75$ to (c) at $t=85$ where the roll straightens and shrinks, before regrowing as in (d) at $t=95$ , and then rotating again as in (e) at $t=99$ before its ultimate decay to the conduction state.

Figure 22.

(a) $L_{1}^{+}$ convecton at $Ra\approx 801.9884$ ( $s\approx 51$ ) used as initial condition for a time simulation in a domain of aspect ratio $2L$ . Snapshots of the simulation taken at $t=76$ (b), $t=85$ (c), $t=95$ (d) and $t=96$ (e) span a quasi-period of the oscillations observed. All states are represented using the same colour scheme as in figure 3.

In the vicinity of the right fold, the $L_{1}^{+}$ solution gains two more unstable eigendirections in rapid succession: at $Ra\approx 818.29402$ , a second mode with a real growth rate and same symmetry as the odd parity drift mode becomes unstable and its eigenvalue collides with that of the odd parity eigenmode to form a complex pair of eigenvalues at $Ra\approx 818.29409$ . At this stage and right before the saddle node, the solution stability is thus $(0,2)$ . The saddle-node bifurcation corresponding to the right fold occurs at $Ra\approx 818.29457$ and beyond this point the solution is $(1,2)$ unstable. By $Ra\approx 818.29446$ the eigenvalues that collided right before the saddle node recombine on the positive real axis and then split, producing a state that is $(3,1)$ unstable. This behaviour is shown in figure 23. The solution remains $(3,1)$ unstable until the leftmost fold where the stability again changes in a non-trivial way. This fold occurs at $Ra\approx 745.35871$ where the convecton gains an additional unstable eigendirection. At $Ra\approx 745.35876$ , above the fold, a further eigenmode becomes unstable making the solution $(5,1)$ unstable. The two new unstable modes collide with the two unstable modes inherited from the dynamics close to $s=143.7$ , as shown in the right panel of figure 19, forming two pairs of unstable complex modes and making the solution $(1,3)$ unstable. This process repeats all along the $L_{1}^{+}$ branch and the solution gains two pairs of unstable complex conjugate eigenvalues after each snaking oscillation.

Figure 23.

Enlargement of the $s\in [67.152,67.16]$ interval of the top panel of figure 18. The instability threshold $e^{\unicode[STIX]{x1D706}_{r}}=1$ is indicated by a thin horizontal line and the saddle node corresponds to the change of stability indicated by the black circle.

The stability of the $S_{\unicode[STIX]{x1D6E5}}$ -symmetric $L_{1}^{-}$ solutions follows the same pattern: after a pair of initial collisions (two collisions occur, see figure 18), the solution is $(1,2)$ unstable along most of the positive segment of the branch. Like $L_{1}^{+}$ this state acquires an additional unstable mode close to the rightmost saddle node and the associated eigenvalue collides with that of the (unstable) drift eigenmode shortly thereafter to form an additional unstable complex conjugate pair of eigenvalues. Approaching the first fold, the $L_{1}^{-}$ solution is $(0,3)$ unstable. It acquires a new unstable eigendirection at the saddle node and shortly thereafter the most recent complex pair of eigenvalues collides on the positive real axis and splits. The state thus becomes $(3,2)$ unstable. Overall, the passage through this saddle node (and its vicinity) has added two unstable eigenvalues to the $L_{1}^{-}$ state and it remains $(3,2)$ unstable along most of the subsequent negative segment. The branch then gains two more unstable eigendirections which collide with the previous two around the left fold, where the solution becomes $(1,4)$ unstable and the process then repeats.

From the previous discussion, it follows that all the secondary localized states are unstable and become more so as the structure broadens. The most unstable eigenmodes along the positive segments of the $L_{1}^{+}$ and $L_{1}^{-}$ branches correspond to oscillatory instabilities arising after eigenvalue collisions (see figure 18) with growth rate in the vicinity of $e^{\unicode[STIX]{x1D706}_{r}}=1.2$ . As shown in figure 20, the magnitude of the corresponding imaginary part of the complex eigenvalues reaches a maximum at each right fold, and then decreases towards the next left fold. Importantly, as we go up the snake, the number of unstable complex conjugate eigenvalues increases while their real parts remain comparable, and so do their imaginary parts. The presence of these eigenvalues is typical of non-gradient systems (Burke & Dawes Reference Burke and Dawes2012) and implies the emergence of non-trivial oscillatory dynamics in the vicinity of the corresponding fixed points. In contrast, the leading eigenvalues are real along most of the negative segments of the snake. Figure 24 summarizes all the unstable eigenmodes associated with an $L_{1}^{+}$ convecton computed at $Ra\approx 780.7143$ along the first negative segment of the branch as it begins snaking (see figure 24 a). The oscillatory instability corresponding to panels (b) and (c) has complex growth rate $0.1485\pm 0.5400i$ and, if triggered, is responsible for oscillatory tilt of the central roll about the vertical. This oscillatory instability is inherited from the eigenvalue collision right after the bifurcation of the $L_{1}^{+}$ state from the $L^{+}$ state but does not dominate along most of the negative branch segment. The dominant modes are instead the eigenmodes shown in (d) and (e), with real growth rates $0.2861$ and $0.2860$ , respectively. The former leads to solutions with no remaining reflection symmetry and is akin to the eigenmode represented in figure 11(c): this mode creates a bifurcation to drifting rung states close to each fold. The latter is the amplitude mode that changes stability at the fold. Lastly, the odd parity drift mode in panel (f) is once again close to marginal (its growth rate is $7.715\times 10^{-4}$ ). Thus, the growth of the rolls is observed first and is then followed by their rotation about the $x$ axis. To understand the stability changes along the secondary snaking branches, we represent in figure 25 all the unstable eigenmodes of a $L_{1}^{+}$ convecton at $Ra\approx 779.7672$ corresponding to the next positive branch segment. All three rolls are now of approximately the same strength (compare figure 25 a to figure 24 a) and the convecton is now 7 times unstable as opposed to being 5 times unstable as in figure 24(a). The mode responsible for the oscillations of the central roll is still unstable with growth rate $0.1379+0.5579i$ . Two new unstable eigenmodes have appeared and undergone collisions with those shown in figure 24(d,e). These modes are shown next to those they are coupled to: the eigenmodes in figure 25(d,e) have growth rate $0.1605\pm 0.4970i$ and are responsible for in-phase oscillations of the axes of the outer rolls of the convecton with respect to the $y$ direction, while those in figure 25(f,g) have growth rate $0.1571\pm 0.4959i$ and generate out-of-phase oscillations on the part of the outer rolls. Lastly, the odd parity drift mode is still present with the same small growth rate of $1.067\times 10^{-3}$ (figure 25 h).

Figure 24.

(a) $L_{1}^{+}$ convecton at $Ra\approx 780.7143$ ( $s\approx 105$ ) on the first negative segment represented using the same colour scheme as in figure 3. All unstable eigenmodes are represented as in figure 9(b): panels (b) and (c) show the real and imaginary parts of the complex conjugate modes, while (d) and (e) show the new steady eigenmodes and (f) the odd parity drift mode.

Figure 25.

(a) $L_{1}^{+}$ convecton at $Ra\approx 779.7672$ ( $s\approx 174$ ) on the second positive segment represented using the same colour scheme as in figure 3. All unstable eigenmodes are represented as in figure 9(b): panels (b) and (c) show the real and imaginary parts of the complex conjugate modes associated with the first eigenvalue collision (matching figure 24 b,c), while panels (d) and (e), and (f) and (g), show the real and imaginary parts of the newly formed complex conjugate modes; (h) is the odd parity drift mode.

The process by which the stability properties evolve during the secondary snaking scenario repeats after each back and forth oscillation and each time is responsible for the addition of two pairs of unstable complex conjugate eigenvalues to the convecton spectrum. The evolution described here repeats for all subsequent branches of secondary states (not shown) with one main difference – the number of unstable eigenvalues inherited from the primary convecton branch: $L_{2\unicode[STIX]{x1D6E5}}^{+}$ and $L_{2c}^{+}$ convectons are $(6,0)$ unstable before a double eigenvalue collision making them $(2,2)$ unstable while $L_{2\unicode[STIX]{x1D6E5}}^{-}$ and $L_{2c}^{-}$ convectons are originally $(7,0)$ unstable and turn into a $(3,2)$ unstable state. Moreover, since the new rolls are smaller than those on the primary convecton branches, it follows that the secondary snaking proceeds further than the primary snaking for given length of structure, implying that the twisted states can be more unstable than the primary states.

3.3 Depinning

The main differences between the two-dimensional system studied by Bergeon & Knobloch (Reference Bergeon and Knobloch2008b ) and the three-dimensional system under investigation here can be understood in terms of the additional $S_{y}$ symmetry. This symmetry gives rise to twist instabilities that are absent in the two-dimensional problem and these instabilities are in turn responsible for the emergence of secondary convectons. Moreover, it also impacts the dynamics outside the snaking region. In a typical depinning scenario, the convecton length grows in time via successive nucleation of additional rolls on either side of the structure. These nucleation events are interspersed with intervals of stasis during which the structure is almost steady. The time spent in the latter state grows with decreasing distance from the edge of the pinning region and does so as the inverse square root of the distance (Burke & Knobloch Reference Burke and Knobloch2006) but this relationship has not been tested in fluid systems.

Figure 26.

Space–time dynamics at $Ra=808$ initialized with an $L^{+}$ convecton at the lower right fold corresponding to $Ra\approx 805.2252$ with no imposed symmetry (a) and with $S_{y}$ symmetry imposed (b), and at $Ra=809$ with no imposed symmetry (c). The solution is represented in the same way as in figure 17.

Figure 26 shows the dynamics at $Ra=808$ of a convecton just outside the pinning region. The simulation was initialized using a convecton at the lower right saddle node on the $L^{+}$ branch ( $Ra\approx 805.2252$ ) and run with and without imposing the symmetry $S_{y}$ . The figure highlights the fundamental differences between $S_{y}$ -symmetric dynamics and unconstrained direct numerical simulation (DNS). The symmetric dynamics behaves like the simulations reported by Bergeon & Knobloch (Reference Bergeon and Knobloch2008b ): the front depins and rolls are sequentially added to both sides of the convecton until the roll pattern fills the domain. This is not the case in the unconstrained DNS where the pattern decays to the conduction state when $t$ reaches $t\approx 60$ , as also found by Bergeon & Knobloch (Reference Bergeon and Knobloch2008b ). The decay time in fact depends quite sensitively on the state used to initialize the simulation and does not vary monotonically with increasing $Ra$ . For example, at $Ra=807$ , it takes $40$ units of time to decay and no nucleation takes place (not shown but see figure 28 a for a similar situation). At $Ra=808$ , it takes $60$ units of time to decay but there are ‘some’ nucleation events, see figure 26(a). At $Ra=809$ (figure 26 c), it takes approximately $50$ units of time: the nucleation event observed at $t=35$ for $Ra=808$ now occurs at $t=30$ and the nucleated rolls survive until $t\approx 48$ (a similar time to $Ra=808$ ) but the central roll reappears at $t=40$ only to vanish at $t=45$ , whereas it only appeared at $t=50$ for $Ra=808$ . Another fundamental difference concerns the behaviour of the central roll: this roll remains essentially unchanged in the symmetric simulation while it tilts and then decays between $t=30$ and $t=35$ in the DNS before temporarily regrowing. Despite this novel behaviour nucleation of new rolls still takes place and does so on a similar time scale as in the symmetric simulation. However, the DNS state at $t=40$ takes the form of two rolls separated by a roll-wide gap as opposed to a three-roll pattern in the symmetric simulation. The reappearance in the DNS of the central roll at $t\approx 50$ coincides with the decay of the two previously nucleated rolls. In each case, the decay of a roll is preceded by rotation of its axis suggesting that the decay is due to increased dissipation arising from the resulting closer approach to the adiabatic walls. Both simulations reveal a characteristic overshoot when a roll is nucleated, a property indicative of a non-gradient system and observed by Bergeon & Knobloch (Reference Bergeon and Knobloch2008b ).

Figure 27.

Time $\unicode[STIX]{x1D6FF}$ at which the first nucleation event occurs as a function of the Rayleigh number $Ra$ for direct numerical simulations (‘ $\times$ ’ signs) and for simulations with the symmetry $S_{y}$ imposed (‘ $+$ ’ signs). The straight line corresponds to the relation $\unicode[STIX]{x1D6FF}=65(Ra-805.22523)^{-0.5}$ expected on the basis of general theory.

To compare the dynamics of our DNS to established depinning dynamics, we report in figure 27 the time at which the first nucleation event occurs in simulations initialized by the same solution as in figure 26. The figure shows a relatively good agreement between the inverse square root law found in other systems (Knobloch Reference Knobloch2015) and the simulations with the symmetry $S_{y}$ imposed. This agreement indicates that the two-dimensional system of Bergeon & Knobloch (Reference Bergeon and Knobloch2008b ) behaves in a similar fashion to the Swift–Hohenberg equation (Burke & Knobloch Reference Burke and Knobloch2006). On the other hand, the DNS provides substantially different results. Sufficiently far away from the fold, we recover the inverse square root law. However, for $Ra<807.6$ , this law breaks down and no depinning event is observed. The reason for this failure is that with no imposed symmetry, all the computed convecton states are unstable with respect to the twist instability and in the vicinity of a fold the time scale involved in the depinning process is so long that this instability kicks in before any nucleation can occur. A precursor to this behaviour can be seen in figure 26(a) where the central roll dies shortly after undergoing a rotation about the $x$ axis, and does so at the same time as the nucleation of a pair of side rolls is taking place. For comparison, at $Ra=807.6$ (not shown), the central roll decays on a similar time scale but the nucleation of the side rolls takes much longer, up to $45$ units of time. If the Rayleigh number is decreased below $807.6$ , the central roll decays even further and there is no seed left for nucleation (see figure 28 a for $Ra=806$ ). Further away from the snaking interval, nucleation occurs much more quickly than twist and several nucleation events can be observed before any of the rolls decay, as exemplified in figure 28(b) for $Ra=834$ .

Figure 28.

Space–time dynamics at (a) $Ra=806$ and (b) $Ra=834$ initialized with an $L^{+}$ convecton at the lower right fold corresponding to $Ra\approx 805.2252$ with no imposed symmetry. The solution is represented in the same way as in figure 17.

It remains to consider what happens when the Rayleigh number $Ra$ moves past $Ra\approx 850.78$ , the location of the primary pitchfork bifurcation. For $Ra\lesssim 850.78$ all solutions, as far as we can tell, decay to the stable conduction state. This is no longer possible once $Ra\gtrsim 850.78$ since the primary state is now unstable. In this regime (figure 29) we find that the system exhibits relaxation oscillations resembling bursts. The figure shows that the flow frequently collapses towards the conduction state but does so in a spatially incoherent manner. As a result the rolls never completely disappear, and since the conduction state is now weakly unstable, they eventually grow back. As revealed in the figure both the collapse and the regrowth proceed by front propagation with well-defined speeds whereby the conduction state invades the roll state and vice versa. Moreover, because the collapse is spatially incoherent this behaviour results in irregular oscillations of relaxation type. Evidently this behaviour is a consequence of the competition between roll decay as a result of the twist instability and the linear instability of the conduction state acting as a regrowth mechanism, mediated via back-and-forth front propagation that maintains spatial asymmetry at all times, thereby guaranteeing spatially incoherent collapse and regrowth. The length of time the system spends near the conduction state depends on the initial condition prior to the initiation of the collapse event. Figure 29 shows episodes when the domain is almost devoid of convection ( $t\approx 40$ and $t\approx 175$ ) as well as others when it is populated by a single roll ( $t\approx 60$ and $t\approx 95$ ). Indeed figure 29(a) reveals the presence of a kind of ‘quantization’ reflecting the presence of $0,1,2,\ldots$ rolls in the domain with abrupt and irregular transitions between these states. Despite these types of regularity the final state represents spatio-temporal chaos since rolls may regrow at random locations owing to the instability of the background state. Figure 29(b) shows the appearance at $t\approx 190$ of an isolated roll near $x=18$ .

Figure 29.

Relaxation oscillations observed at $Ra=850.80$ , in the vicinity of the primary instability ( $Ra\approx 850.78$ ), shown through the temporal evolution of the kinetic energy (a) and through the corresponding space–time plot (b). The origin of time is chosen to fit the beginning of the visualized time frame.

The situation resembles the behaviour of binary fluid convection in a horizontal layer just prior to the appearance of stable stationary convectons (Batiste et al. Reference Batiste, Knobloch, Alonso and Mercader2006) where it is a focusing instability that leads to intermittent dynamics from which the stationary convectons emerge.