Interactions between two adjacent convection rolls in turbulent Rayleigh-Benard convection

TL;DR

This study investigates the interactions and stability of two adjacent convection rolls in turbulent Rayleigh-Benard convection, revealing how turbulent thermal diffusivity influences their dynamics and stability, including the effects of temperature differences and spontaneous state switching.

Contribution

The paper introduces a stochastic differential equation model predicting stable states and interactions of neighboring convection rolls, linking turbulent thermal diffusivity to their stability and dynamics.

Findings

Stable co- and counter-rotating states identified

Temperature differences influence the stability of convection roll states

Spontaneous switching between states observed and explained

Abstract

Rayleigh-B{\'e}nard convection experiments were done with two adjacent cubic cells with a partial wall in between to force the generation of two interacting convection rolls. Observed stable states include both counter-rotating and co-rotating states. The stability of each of these states and their dynamics were modeled by stochastic ordinary differential equations of motion in terms of the orientation, amplitude, and mean temperature of each convection roll. The form of the interaction terms is predicted based on an effective turbulent diffusion of temperature between the adjacent rolls. Predictions are made for stable fixed points of the co- and counter-rotating states. This suggests that the same turbulent thermal diffusivity that describes macroscopically averaged heat transport also controls the interactions between neighboring convection rolls. The surprising stability of co-rotating states is due to the temperature difference between the neighboring rolls becoming large enough that the heat flux between the rolls stabilizes the temperature profile of aligned co-rotating states. This temperature difference can be driven by heating the plates of the two cells to different mean temperatures. This shifts the orientations of the rolls of counter-rotating states in opposite directions, and for large temperature differences only co-rotating states are stable Spontaneous switching between co-rotating and counter-rotating states is also observed. Switching to counter-rotating states occurs mainly due to cessation (a significant weakening of a convection roll), which reduces damping on changes in orientation, allowing the orientation to change rapidly due to diffusive fluctuations. Switching to co-rotating states is mainly driven by smaller diffusive fluctuations, which have a positive feedback that destabilizes the counter-rotating state.

Tables

Table 1. Table 1: Parameters that are relevant to the model of turbulent thermal diffusion for the interaction between neighboring rolls. Parameters referred to as ‘known’ have known values from previous experiments with a single LSC. The ‘measured’ parameters are those that describe the LSCs. The ‘dependent’ parameters are useful in results expressed in terms of differences between neighboring rolls and cells. The parameter κ t subscript 𝜅 𝑡 \kappa_{t} is unknown, and is the single new fit parameter.

	definition	type
$L$	cell height	controlled
$A$	fraction of interface open	controlled
$\kappa$	thermal diffusivity	controlled
$T_{m,i}$	mean temp. of plates	controlled
$\Delta T_m$	temp. diff. of adjacent plates	dependent ($\Delta T_{m,i}$)
${\theta }_i$	LSC orientation	measured
${\delta }_i$	LSC amplitude	measured
$T_{0,i}$	LSC mean temperature	measured
$\Delta T_0$	temp. diff. between LSCs	dependent ($\Delta T_{0,i}$)
${\delta }_0$	stable LSC amplitude	known Ji and Brown (2020b)
${\tau }_{\delta }$	damping timescale of ${\delta }_i$	known Ji and Brown (2020b)
${\tau }_{\theta }$	damping timescale of ${\theta }_i$	known Ji and Brown (2020b)
$D_{\delta }$	turbulent diffusivity for ${\delta }_i$	known Ji and Brown (2020b)
$D_{\dot{\theta }}$	turbulent diffusivity for ${\theta }_i$	known Ji and Brown (2020b)
$D_T$	turbulent diffusivity for $T_{0,i}$	not tested
Nu	Nusselt number	known Funfschiling et al. (2005)
${\kappa }_t$	turbulent thermal diffusivity	fit

Table 2. Table 2: Fractional changes of various parameters towards their new steady state values before switching events to counter- and co-rotating states. Switching to counter-rotating states is mostly driven by cessations ( d c subscript 𝑑 𝑐 d_{c} ) during the event, in which a reduction in one of the δ i subscript 𝛿 𝑖 \delta_{i} reduces the damping on motion in θ i subscript 𝜃 𝑖 \theta_{i} which can be driven to large changes by random fluctuations. Switching to co-rotating states is driven by a mix of diffusive fluctuations in θ 1 − θ 2 subscript 𝜃 1 subscript 𝜃 2 \theta_{1}-\theta_{2} , T 0 , 2 − T 0 , 1 subscript 𝑇 0 2 subscript 𝑇 0 1 T_{0,2}-T_{0,1} , and δ i subscript 𝛿 𝑖 \delta_{i} which have a positive feedback that destabilizes the counter-rotating state.

parameter	change in	to counter-rotating	to co-rotating
$d_{\theta }$	$\left|{\theta }_1-{\theta }_2\right|$	0.12	0.37
$d_T$	$\left|T_{0,2}-T_{0,1}\right|$	0.01	0.52
$d_{\delta }$	$\left({\delta }_1+{\delta }_2\right)/2$	-0.23	0.53
$d_c$	${\overline{\delta }}_i-{\delta }_i$	0.36	0.31

Equations

$$T_i=T_{0,i}+{\delta }_{i}cos(\theta -{\theta }_i)$$

$${\dot{\delta }}_i=\frac{{\delta }_i}{{\tau }_{\delta }}-\frac{{\delta }_i^{3/2}}{{\tau }_{\delta }\sqrt{{\delta }_0}}+f_{\delta }(t)+{\dot{\delta }}_{\kappa ,i}.$$

$${\ddot{\theta }}_i=-\frac{{\dot{\theta }}_i{\delta }_i}{{\tau }_{\dot{\theta }}{\delta }_0}-\nabla V_g({\theta }_i)+f_{\dot{\theta }}(t)+{\ddot{\theta }}_{\kappa ,i}.$$

$$\dot{\delta}_{\kappa,1}=\frac{A\kappa_{t}}{\pi L^{2}}\bigg{[}\sqrt{2}(T_{0,2}-T_{0,1})\cos\theta_{1}\\ +\delta_{2}\sin(\theta_{1}-\theta_{2})\sin(2\theta_{1})/2\\ +[\delta_{2}\cos(\theta_{1}-\theta_{2})-\delta_{1}][\pi/4+\cos(2\theta_{1})/2]\bigg{]}\\ -\frac{\sqrt{2}\dot{T}_{0,1}}{\pi}\cos\theta_{1}\ .$$

$$\ddot{\theta}_{\kappa,1}=-\frac{A\kappa_{t}}{\pi L^{2}\tau_{\dot{\theta}}\delta_{0}}\bigg{[}\delta_{2}\sin(\theta_{1}-\theta_{2})[\pi/4-\cos(2\theta_{1})/2]\\ +[\delta_{2}\cos(\theta_{1}-\theta_{2})-\delta_{1}]\sin(2\theta_{1})/2+\sqrt{2}(T_{0,2}-T_{0,1})\sin\theta_{1}\bigg{]}\\ +\frac{\sqrt{2}\dot{T}_{0,1}}{\pi\delta_{1}}\sin\theta_{1}$$

$$\begin{gathered} {\dot{T}}_{0,1}=f_T(t)+\frac{4\text{mbox}Nu\kappa }{L^2}(T_{m,1}-T_{0,1})+ \\ \frac{A{\kappa }_t}{L^2}\left[T_{0,2}-T_{0,1}+0.88({\delta }_2\cos {\theta }_2-{\delta }_1\cos {\theta }_1)\right]. \end{gathered}$$

$${\ddot{\theta }}_{1,counter}\approx -0.79\frac{A{\kappa }_t}{\pi L^2{\tau }_{\dot{\theta }}}({\theta }_1-{\theta }_2).$$

$$\ddot{\theta}_{1,co}\approx\frac{A\kappa_{t}\bar{\delta}_{co}}{\pi L^{2}\delta_{0}\tau_{\dot{\theta}}}\bigg{[}+0.28(\theta_{1}-\theta_{2}-\pi)\\ -\left(\frac{1.41\Delta\bar{T}_{0,co}}{\bar{\delta}_{co}}-2\right)(\theta_{1}-\theta_{p,1})\ .\bigg{]}$$

$${\ddot{\theta }}_i=a\sin ({\theta }_i-{\theta }_{p,i})+b\sin ({\theta }_1-{\theta }_2)$$

$$p({\theta }_i)=\exp \left(-\frac{V({\theta }_i){\bar{\delta }}^2}{D_{\dot{\theta }}{\tau }_{\dot{\theta }}{\delta }_0^2}\right)$$

$$\Delta\bar{T}_{0,co}=\pm\frac{1.76}{1+2\mbox{Nu}\kappa/A\kappa_{t}}\bar{\delta}_{co}\$$

$${\bar{\delta }}_{co}\approx {\delta }_0+\frac{(2A{\kappa }_t{\tau }_{\delta }/\pi L^2)(1.41|\Delta {\bar{T}}_{0,co}|-2.57{\delta }_0)}{1+2.57A{\kappa }_t{\tau }_{\delta }/\pi L^2}.$$

$$\Delta {\theta }_p=\pm \frac{A{\kappa }_t\Delta {\bar{T}}_0}{\pi L^2\bar{\delta }{\omega }_r^2{\tau }_{\dot{\theta }}}.$$

$$\Delta {\bar{T}}_{0,co}(\Delta T_m)=\frac{\pm 1.76{\bar{\delta }}_{co}+(2\text{mbox}Nu\kappa /A{\kappa }_t)\Delta T_m}{1+2\text{mbox}Nu\kappa /A{\kappa }_t}.$$

$$\Delta {\bar{T}}_0=\frac{-0.88({\bar{\delta }}_2\cos {\theta }_{p,2}-{\bar{\delta }}_1\cos {\theta }_{p,1})+(2\text{mbox}Nu\kappa /A{\kappa }_t)\Delta T_m}{1+2\text{mbox}Nu\kappa /A{\kappa }_t}.$$

$$\Delta \bar{\delta }=2{\tau }_{\delta }{\left(\frac{{\dot{\delta }}_{\kappa ,2}-{\dot{\delta }}_{\kappa ,1}}{2}\right)}_{counter}=\frac{2{\tau }_{\delta }A{\kappa }_t}{\pi L^2}(\bar{\delta }\Delta {\theta }_p-\Delta {\bar{T}}_0).$$

$${\bar{\delta }}_{counter}={\delta }_0+2{\tau }_{\delta }{\overline{\dot{\delta }}}_{counter}={\delta }_0+\frac{3.6A{\kappa }_t{\tau }_{\delta }}{\pi L^2}\Delta T_m^2$$

$$S=\frac{10A{\kappa }_t{\delta }_0[\sqrt{2}|\Delta {\bar{T}}_{0,co}|/{\bar{\delta }}_{co}-2]}{{\pi }^2{L}^2{\tau }_{\dot{\theta }}{\omega }_r^2{\bar{\delta }}_{co}}.$$

$${\dot{T}}_1(\theta )=\frac{A{\kappa }_t}{L^2}\left[T_2(\theta ,{\theta }_2)-T_1(\theta ,{\theta }_1)\right]+\frac{4\text{mbox}Nu\kappa }{H^2}[T_{m,1}-T_{0,1}].$$

$${\dot{T}}_1(\theta )={\dot{T}}_{0,1}+{\delta }_1{\dot{\theta }}_1\sin (\theta -{\theta }_1)+{\dot{\delta }}_1\cos (\theta -{\theta }_1).$$

$${\dot{\theta }}_{\kappa ,1}=\frac{1}{\pi {\delta }_1}{\int }_0^{2\pi }[{\dot{T}}_1(\theta )-{\dot{T}}_{0,1}]\sin (\theta -{\theta }_1)d\theta$$

$${\dot{\delta }}_{\kappa ,1}=\frac{1}{\pi }{\int }_0^{2\pi }[{\dot{T}}_1(\theta )-{\dot{T}}_{0,1}]\cos (\theta -{\theta }_1)d\theta .$$

$$\dot{\theta}_{\kappa,1}=-\frac{A\kappa_{t}}{\pi L^{2}\delta_{1}}\bigg{[}\delta_{2}\sin(\theta_{1}-\theta_{2})[\pi/4-\cos(2\theta_{1})/2]\\ +[\delta_{2}\cos(\theta_{1}-\theta_{2})-\delta_{1}]\sin(2\theta_{1})/2+\sqrt{2}(T_{0,2}-T_{0,1})\sin\theta_{1}\bigg{]}\\ +\frac{\sqrt{2}\dot{T}_{0,1}}{\pi\delta_{1}}\sin\theta_{1}$$

$$\dot{\delta}_{\kappa,1}=\frac{A\kappa_{t}}{\pi L^{2}}\bigg{[}\sqrt{2}(T_{0,2}-T_{0,1})\cos\theta_{1}\\ +\delta_{2}\sin(\theta_{1}-\theta_{2})\sin(2\theta_{1})/2\\ +[\delta_{2}\cos(\theta_{1}-\theta_{2})-\delta_{1}][\pi/4+\cos(2\theta_{1})/2]\bigg{]}\\ -\frac{\sqrt{2}\dot{T}_{0,1}}{\pi}\cos\theta_{1}\ .$$

$${\ddot{\theta }}_{\kappa ,1}=\frac{{\dot{\theta }}_{\kappa ,1}{\delta }_1}{{\tau }_{\dot{\theta }}{\delta }_0}.$$

$$\begin{gathered} {\dot{T}}_{0,1}=f_T(t)+\frac{4\text{mbox}Nu\kappa }{H^2}(T_{m,1}-T_{0,1})+\frac{A{\kappa }_t}{2L^2}\times \\ {\int }_{-\pi /4}^{\pi /4}\frac{T_{0,2}-T_{0,1}+{\delta }_2\cos (\theta -{\theta }_2)-{\delta }_1\cos (\theta -{\theta }_1)}{{\cos }^2\theta }d\theta . \end{gathered}$$

$$\begin{gathered} {\dot{T}}_{0,1}=f_T(t)+\frac{4\text{mbox}Nu\kappa }{H^2}(T_{m,1}-T_{0,1})+ \\ \frac{A{\kappa }_t}{L^2}\left[T_{0,2}-T_{0,1}+0.88({\delta }_2\cos {\theta }_2-{\delta }_1\cos {\theta }_1)\right]. \end{gathered}$$

$$\ddot{\theta}_{1,co}\approx\frac{A\kappa_{t}\bar{\delta}_{co}}{\pi L^{2}\delta_{0}\tau_{\dot{\theta}}}\bigg{[}+0.28(\theta_{1}-\theta_{2}-\pi)\\ -\left(\frac{1.41\Delta\bar{T}_{0,co}}{\bar{\delta}_{co}}-2\right)(\theta_{1}-\theta_{p,1})\ .\bigg{]}$$

Full text

I Introduction

While turbulent flows are often thought of as irregular and erratic, large-scale coherent flow structures are commonplace in turbulence. An example is convection rolls driven by buoyancy in natural convection. Such structures and their dynamics can play a significant role in heat and mass transport. A particular challenge that is the focus of this manuscript is to develop a model for how these large-scale flow structures interact with each other, for example to result in neighboring convection rolls that are counter-rotating or co-rotating.

We investigate this in the model system of turbulent Rayleigh-Bénard convection. In Rayleigh-Bénard convection, a fluid is heated from below and cooled from above to generate buoyancy-driven flow Ahlers et al. (2009); Lohse and Xia (2010). This system exhibits convection rolls which are robust large-scale coherent structures that retain a similar organized flow structure over a long time. For example, in containers of aspect ratio near 1, a large-scale circulation (LSC) forms. This LSC consists of localized blobs of coherent fluid known as plumes. The plumes collectively form a single convection roll in a vertical plane that can be identified by averaging over the flow field or timescales longer than the circulation period Krishnamurti and Howard (1981). This LSC spontaneously breaks the symmetry of symmetric containers, but turbulent fluctuations cause the LSC orientation $\theta_{0}$ in the horizontal plane to meander spontaneously and erratically, and allow it to sample different orientations to recover the symmetry over long times Brown and Ahlers (2006a). While the LSC exists nearly all of the time, on rare occasions these fluctuations lead to spontaneous cessations followed by reformation of the LSC Brown and Ahlers (2006a); Xi and Xia (2007). The LSC exhibits oscillation modes Heslot et al. (1987); Sano et al. (1989); Castaing et al. (1989); Ciliberto et al. (1996); Takeshita et al. (1996); Cioni et al. (1997); Qiu et al. (2000); Qiu and Tong (2001); Niemela et al. (2001); Qiu and Tong (2002); Qiu et al. (2004); Sun et al. (2005); Tsuji et al. (2005), which in circular cylindrical containers consists of twisting and sloshing Funfschilling and Ahlers (2004); Xi et al. (2009); Zhou et al. (2009); Brown and Ahlers (2009), and at some aspect ratios a jump-rope-like mode Vogt et al. (2018); Horn et al. (2022).

The qualitative behavior of the LSC depends on the geometry of the cell, which is necessary to account for if we are to understand the interactions between neighboring rolls in some geometry. In containers with rectangular horizontal cross-sections, the preferred alignment of the LSC orientation $\theta_{0}$ is along the longest diagonals, and the LSC orientation can spontaneously switch between adjacent corners Liu and Ecke (2009); Song et al. (2014); Bai et al. (2016); Foroozani et al. (2017); Giannakis et al. (2018); Vasiliev et al. (2016, 2019). A regular oscillation of $\theta_{0}$ can occur between nearest-neighbor diagonals in a non-square rectangular cross-section Song et al. (2014). In a cubic container, the oscillation structure of the LSC corresponds to an advected oscillation mode with one oscillation per LSC turnover period Ji et al. (2020), and an oscillation in the shape of the temperature profile that does not occur in a circular cross-section cell Ji and Brown (2020a).

In principle, the Navier-Stokes equations describe fluid flow, but they are impractical to solve for such complex turbulent flows, so low-dimensional models are desired. It has long been recognized that the states and dynamics of a single large-scale coherent structure are similar to those of low-dimensional dynamical systems models Lorenz (1963) and stochastic ordinary differential equations Brown and Ahlers (2007a); de la Torre and Burguete (2007); Thual et al. (2014); Rigas et al. (2015). Low-dimensional models can potentially describe parameters of stable states, dynamics such as oscillation modes, and behavior of transitions between states. While low-dimensional models lack detail of smaller structures, they have the advantage of being much simpler to solve and understand the behavior because they involve simpler equations – such as ordinary differential equations instead of partial differential equations.

There are several low-dimensional models for single roll LSCs. Early models tried to characterize flow reversals in simplified flow in a two-dimensional plane Sreenivasan et al. (2002); Benzi (2005); Fontenele Araujo et al. (2005) However, these two-dimensional models could not characterize more complex three-dimensional dynamics with motion in $\theta_{0}$ such as reorientations and oscillation modes. Some models are obtained by transforming high-dimensional fluid velocity field data (usually from direct numerical simulation) and reducing it to lower-dimensional models consisting of a few highest-energy Fourier modes or eigenmodes. These models have been able to characterize the detailed shape of the LSC and its dynamics, including flow reversals in two dimensions Chandra and Verma (2011); Podvin and Sergent (2015) spontaneous corner-switching and oscillations in cubic cells Giannakis et al. (2018); Vasiliev et al. (2019), and the twisting, sloshing, and jump-rope oscillation modes of the LSC Horn et al. (2022). Because these models are obtained from high-dimensional data, they are descriptive in higher detail, but the models are not formulated in such a way to predict flows when detailed data is not already available. We desire a low-dimensional modeling technique that can be predictive and generalizable to other systems with more limited or no experimental input required.

We build off the low-dimensional model of Brown & Ahlers, where model terms are derived from approximations of the Navier-Stokes equations Brown and Ahlers (2008a). The model consists of a pair of stochastic differential equations for the LSC, in terms of the orientation $\theta_{0}$ and amplitude $\delta$. The model of Brown & Ahlers and its extensions have successfully described most of the known dynamical modes of the LSC in including the meandering, cessations, and twisting and sloshing oscillation modes described above Brown and Ahlers (2008a, b, 2009); Zhong et al. (2017); Sterl et al. (2016), with the exception of the jump rope mode. The combination of twisting and sloshing oscillation modes Funfschilling and Ahlers (2004); Xi et al. (2009); Zhou et al. (2009); Brown and Ahlers (2009) can alternatively be described in this model as a single advected oscillation mode, with two oscillation periods per LSC turnover period in a circular cross-section Brown and Ahlers (2009), or with one oscillation period per LSC turnover period in a cubic cross-section which is excited by a potential $V_{g}$ due to the shape of the container Ji et al. (2020). This potential $V_{g}$ can be predicted as a function of container cross-section geometry without experimental input Brown and Ahlers (2008b); Song et al. (2014); Ji and Brown (2020b). This same potential $V_{g}$ also explains the preferred orientation along diagonals of a rectangular cross-section container, the oscillations between diagonals Song et al. (2014), and the stochastic switching between diagonals Liu and Ecke (2009); Bai et al. (2016); Foroozani et al. (2017); Giannakis et al. (2018); Vasiliev et al. (2016, 2019). This model requires experimental measurements of two diffusivity parameters that characterize the strength of turbulent fluctuations in $\theta_{0}$ and the LSC strength, but these are relatively simple parameters that can be obtained from short term measurements of a single state. Predictions of oscillation frequencies, average rates of stochastic switching between states, widths of probability distributions, and state boundaries are typically accurate within a factor of 3 Brown and Ahlers (2007a, 2008a); Bai et al. (2016); Ji and Brown (2020b), but can be more accurate when more fit parameters are used Assaf et al. (2011).

While systems of aspect ratio close to one tend to have a single convection roll, horizontally extended convection systems tend to consist of multiple convection rolls, each of aspect ratio on the order of one, arranged side-by-side, typically counter-rotating relative to their neighbors Bailon-Cuba et al. (2010); van der Pel et al. (2012). Counter-rotating behavior is claimed to be prevalent in nature, for example, in textbook pictures of convection rolls in the atmosphere, in the convection layer of stars, or planetary cores. However, a few simulations have been able to produce co-rotating states in horizontally adjacent rolls in non-turbulent convection, with lateral heating Pérez-García et al. (2004), or with an inclination relative to gravity of 0.01 rad Mercader et al. (2019). Vertically extended systems often have counter-rotating rolls stacked on top of each other Zwirner et al. (2020). A turbulent convection experiment with two fluids, one on top of the other, found a convection roll in each fluid with two stable states, one where the rolls are counter-rotating, and one where the rolls are co-rotating Xie and Xia (2013), with rare stochastic switching between the two states. A variation on stacked co-rotating rolls was observed in simulations of a non-turbulent vertical channel where vertical flow in opposite directions occurred along opposite side walls. In this case, the co-rotating rolls are forced by the opposite vertical flows on either side Gao et al. (2013).

Two-dimensional theories using linear stability analysis in the non-turbulent regime have been able to predict the existence of stable counter- and co-rotating states. In one case, a lateral heating was shown to provide a forcing to produce stable co-rotating states in horizontally neighboring rolls Pérez-García et al. (2004). Another model considered vertically stacked rolls, and showed stable counter-rotating states where the coupling was dominated by viscous forces, as well as stable co-rotating states where the coupling was mainly through vertical buoyancy forces Petry and Busse (2003). Since these earlier models were focused on two-dimensions Petry and Busse (2003); Pérez-García et al. (2004), they can identify co-rotating and counter-rotating states, but they are missing the orientation of rolls, which is important because stability in three dimensions requires stability in the orientation as well as the flow strength, and many dynamics of rolls involve changes in orientation Funfschilling and Ahlers (2004); Brown and Ahlers (2006a, b); Vogt et al. (2018); Akashi et al. (2022). These previous models also focused on linear stability analysis of laminar flow, but many natural flows are turbulent, well beyond the range where linear stability analysis rigorously applies. It has yet to be determined whether three-dimensional systems of interacting turbulent convection rolls can be captured by low-dimensional models, which is the primary goal of this work.

We hypothesize that the low-dimensional model of Brown & Ahlers Brown and Ahlers (2008a, b); Bai et al. (2016); Ji and Brown (2020b); Ji et al. (2020) can be extended to systems of multiple convection rolls using an existing set of ordinary differential equations of motion for each roll, and adding interaction terms to the equations of motion for each neighboring roll. We seek to test whether a low-dimensional model can predict or describe the preferred states (e.g. counter-rotating and/or co-rotating), their parameter values (e.g. orientation, flow strength), and their dynamics (e.g. how does switching between states occur). We also seek to understand the physical origin of the interaction terms, in particular due to turbulent thermal diffusion across the interface between the rolls so that a general predictive model can be made. This is in the spirit of low-dimensional models with heat transport in one dimension that have already proven effective, for example, to explain an oscillation in the shape of the temperature profile with an orientation-dependent vertical heat transport Ji and Brown (2020a), and to explain thermal waves with a height-dependent vertical heat transport Urban et al. (2023). In principle, if we can write equations for convection rolls with interaction terms for a neighbor, results from two convection rolls can be extended to systems of more rolls assuming there are no significant interaction terms that involve three or more rolls.

The remainder of this manuscript is organized as follows. Section II explains the experimental apparatus, and methods used to characterize the LSC. Section III summarizes the existing low-dimensional model of Brown & Ahlers for a single LSC to use as a starting point. Section IV describes initial observations of counter-rotating and co-rotating states of two neighboring rolls. Section V presents a physical derivation of a model for neighboring roll interactions based on effective turbulent thermal diffusion. This model is tested in Sec. VI for the symmetric case where the control of both cells is nominally identical. Section VII extends this model to the cases with asymmetric forcing on the two cells due to a difference in mean temperature of the plates, to better understand co-rotating states, and to test some of the model terms. Section VIII presents observations of spontaneous switching between co-rotating and counter-rotating states, and explains how these can be understood with the model.

II Methods

II.1 Experiment setup

The experimental apparatus is the same one used in Ji and Brown (2020b), with the relevant details and modifications presented here. The apparatus consists of two adjacent cells as illustrated in Fig. 1. Each cell is nearly cubic with $H=20.32$ cm and horizontal lengths $L=20.02\pm 0.03$ cm. Since two rolls may be preferred at an aspect ratio 2 without a wall in between van der Pel et al. (2012), the apparatus is modified from Ji and Brown (2020b) with a partial opening in the middle wall to enforce a flow with two convection rolls that can interact with each other.

To control the temperature difference $\Delta T$ between the top and bottom of each cell, water is circulated through top and bottom plates. Each top and bottom plate is controlled by its own temperature-controlled water bath so that $\Delta T$ could be controlled in each cell more precisely, and a difference $\Delta T_{m}$ could be imposed between the mean plate temperatures of the two cells when desired. The plates are aluminum, with double-spiral water-cooling channels as in Brown et al. (2005), except that each plate has its own double-spiral, and the inlet and outlet of each plate were adjacent to minimize the spatial temperature variation within the plates. The baths pump water in a pattern symmetric around the middle wall to minimize temperature asymmetries between the two cells. Each plate has 5 thermistors to record $\Delta T$ and $\Delta T_{m}$, with one thermistor at the center and four on the diagonals, halfway between the center and each corner of the plate. The standard deviation of temperatures within each plate is $0.005\Delta T$ Ji and Brown (2020b). The top and bottom plates are parallel within 0.06*∘*.

The sidewalls of the cubic cells are plexiglas to thermally insulate the cell from the surroundings. The outer sidewalls have a thickness of 0.55 cm, and are further insulated from the room by foam insulation. The middle wall is shared in between the two flow cells. The middle wall has a thickness of $0.91$ cm to thermally insulate the cells from each other, and extends 2.2 cm into the gap between the top and bottom plates to thermally insulate the plates of different cells from each other. The middle wall has gaps in it to allow flow in the two cells to interact with each other. Three sections of height $H_{1}=13.6\pm 0.1$ cm by $L_{1}=4.9\pm 0.1$ cm are cut out of the middle wall, as shown in Fig. 2. A quarter circle with radius $0.35\pm 0.05$ cm was cut out of each corner of the cell to allow gas bubbles to escape the cell during a degassing process to prepare the working fluid for experiments. This results in a fraction $A=0.50$ of the middle wall open to allow interaction between the two cells. All of the cutouts are symmetric around the mid-height of the cell and symmetric from left to right.

Two vertical grooves were cut in the middle wall as shown in Fig. 2 to place thermistors and run their wiring out the top of the cell. The remaining area of the grooves was filled with epoxy to keep the walls as flat as possible, as explained in Ji and Brown (2020b). Detailed internal dimensions are given in Ji and Brown (2020b), where it is shown that their effect on the symmetry of the flow is negligible compared to the effects of the cubic shape, as well as compared to effects of non-uniformity of the heating and cooling plate temperatures.

The apparatus was further insulated from the room as in Brown et al. (2005) by surrounding it with 5 cm thick closed-cell foam, which itself was surrounded on the sides by a copper shield with water circulating through a pipe welded to the shield. The circulating water temperature was set to match the mean temperature $\langle T_{0}\rangle$ averaged over the two cells, with a standard deviation of $0.006^{\circ}$ K. The shield was surrounded by another layer of 2.5 cm thick open-cell foam.

To measure the LSC, thermistors were mounted in the sidewalls as in Brown and Ahlers (2006b). There are three rows of thermistors: at heights $+H/4$, [math], and $-H/4$ relative to the mid-height of the cell, as shown in Fig. 2. In each row, eight thermistors are equally spaced in the angle $\theta$ around the mid-plane, as illustrated in Fig. 1. The coordinate $\theta$ in the horizontal plane is measured relative to the axis going through the center of both cells, as illustrated in Fig. 1. The LSC orientation $\theta_{1}$ when corresponding to the flow direction near the bottom plate when viewed from above is measured counter-clockwise in cell 1 and $\theta_{2}$ is measured clockwise in cell 2, so that counter-rotating states correspond to $\theta_{1}=\theta_{2}$ when the orientation vectors align head-to-head, and co-rotating states correspond to $\theta_{1}=\theta_{2}+\pi$ rad when they align head-to-tail.

The cell was leveled in the direction perpendicular to the axis going through the center of both cells, with an uncertainty of $0.03^{\circ}$. For some experiments, the cell was intentionally tilted by an angle $\beta$ relative to the level cell along the axis going through the center of both cells to introduce a forcing from buoyancy that breaks the symmetry of the two cells.

The working fluid was degassed and deionized water with mean temperature of 23.0*∘C, for a Prandtl number $Pr=\nu/\kappa=6.41$ where $\nu=9.36\times 10^{-7}$ m2/s is the kinematic viscosity and $\kappa=1.46\times 10^{-7}$ m2/s is the thermal diffusivity. The Rayleigh number is given by $Ra=g\alpha\Delta TH^{3}/\kappa\nu$ where $g$ is the acceleration of gravity, and $\alpha=0.000238$ K-1* is the thermal expansion coefficient.

We calibrated thermistors at five mean temperatures from 21*∘* C to 25*∘* C. The calibrations were run with a small $\Delta T=0.04$ K to enhance mixing. The uncertainty on temperature measurements is 1.9 mK for the sidewall thermistors and 1.2 mK for the plate thermistors, adding in quadrature the contributions from the standard deviation in thermistor temperatures during calibration runs, and the root-mean-square differences between the mean temperatures and the calibration fit.

II.2 Obtaining the LSC parameters

The LSC orientation $\theta_{i}$, amplitude $\delta_{i}$, and mean temperature $T_{0,i}$ are the main parameters used to characterize the LSC. They were obtained using the same methods as Brown and Ahlers (2006a), but are now labeled with indices $i=1,2$ to differentiate the two cells. As the LSC moves hot fluid from near the bottom plate up one side and cold fluid from near the top plate down the other side, a temperature difference is detected along a horizontal direction at the mid-height of the container. We fit the thermistor temperatures at the middle row in each cell to the function

[TABLE]

to get the orientation $\theta_{i}$, temperature amplitude $\delta_{i}$, and spatial mean temperature $T_{0,i}$ of the LSC roll. To obtain a time series, these fits are done at every measured time step, which is typically 7 s.

Due to the frequent failure of thermistors in the interior of the cell, we only report data from the sidewall thermistors, leaving only six thermistors in each cell at the middle row. For example, at $\Delta T=19$ K, where the time-averaged amplitude $\bar{\delta}_{i}$ is 0.42 K in a counter-rotating state, this fit results in average uncertainties on instantaneous measurements of 0.18 rad on $\theta_{i}$ and 0.06 K on $\delta_{i}$ from fitting Eq. 1. Since the uncertainties on individual thermistor measurements are only 1.9 mK, and deviations from the sinusoidal profile due to oscillation of the structure are only about 6% of the amplitude Ji and Brown (2020a), the largest contribution to these random uncertainties is large turbulent fluctuations around the mean temperature profile. The small amplitude of $\Delta T=0.04$ K during calibrations is expected to produce an LSC with mean amplitude $\bar{\delta}_{i}=2.9$ mK based on an extrapolation of measured values of $\bar{\delta}_{i}$ Ji and Brown (2020b), which introduces a systematic error of 2.9 mK on measurements of $\delta_{i}$. Using Eq. 1, this error propagates to a systematic error on $\theta_{i}$ of 0.007 rad at $\Delta T=19$ K where $\bar{\delta}_{i}=0.42$ K, or an error on $\theta_{i}$ of 0.024 rad when $\bar{\delta}_{i}=0.12$ K at $\Delta T=4$K.

III Review of model for a single LSC, i.e. without neighboring roll interactions

In this section we summarize the model of Brown & Ahlers Brown and Ahlers (2008a), which we use as a baseline of comparison to the behavior of a single-roll LSC. The model consists of a pair of stochastic ordinary differential equations, using the empirically known, robust LSC structure as an approximate solution to the Navier-Stokes equations to obtain equations of motion for parameters that describe the LSC dynamics. The effects of fast, small-scale turbulent fluctuations are separated from the slower, large-scale coherent motion when obtaining this approximate solution, then added back in as a stochastic term in the low-dimensional model. The flow strength in the direction of the LSC is represented by the temperature amplitude $\delta_{i}$, which is proportional to the mean flow speed in the LSC. The equation of motion for $\delta_{i}$ is

[TABLE]

The first forcing term on the right side of the equation corresponds to buoyancy, which strengthens the LSC. The second term is a non-linear damping approximated for a boundary-layer dominated flow, which weakens the LSC. $\delta_{0}$ is the stable fixed point value of $\delta_{i}$ where buoyancy and damping balance each other, and $\tau_{\delta}$ is a damping timescale for changes in the strength of the LSC. $f_{\delta}(t)$ is a stochastic forcing term representing the effect of small-scale turbulent fluctuations and is modeled as Gaussian white noise with diffusivity $D_{\delta}$. Thus, $\delta_{i}$ tends to exhibit strong fluctuations around the stable fixed point value of $\delta_{0}$. $\dot{\delta}_{\kappa}$ is a placeholder for forcing terms due to the interaction between neighboring rolls, to be derived in Sec. V.

The equation of motion for the LSC orientation $\theta_{0}$ is

[TABLE]

The first term on the right side of the equation is a damping term which comes from the advective term of the Navier-Stokes equations, where $\tau_{\dot{\theta}}$ is a damping time scale for changes of orientation of the LSC. $f_{\dot{\theta}}$ is another stochastic forcing term with diffusivity $D_{\dot{\theta}}$. $V_{g}$ is a potential which represents the pressure of the sidewalls acting on the LSC, and is a function of the geometry of the cell, so $-\nabla V_{g}(\theta_{i})\equiv-\partial V_{g}/\partial\theta_{i}$ is the forcing due to this geometric potential. $\ddot{\theta}_{\kappa}$ is another placeholder for forcing terms due to the interaction between neighboring rolls, to be derived in Sec. V.

Since the terms of Eqs. 2 and 3 were derived from the Navier-Stokes equations, functional predictions for $\delta_{0}$, $\tau_{\delta}$ and $\tau_{\dot{\theta}}$ exist Brown and Ahlers (2008a); Ji and Brown (2020b). The diffusivities $D_{\delta}$ and $D_{\dot{\theta}}$ have so far been measured from data Brown and Ahlers (2008a); Ji and Brown (2020b).

The geometric potential $V_{g}(\theta_{i})$ can be expressed as a function of the $\theta$-dependent diameter of the cell, and thus can be calculated for convex cell geometries. $V_{g}(\theta_{i})$ is predicted to be inversely proportional to the diameter across the cell as a function of LSC orientation $\theta_{0}$, and thus the lowest potentials are aligned with the longest diagonals Brown and Ahlers (2008b). In the case of a cubic cell, this corresponds to four potential minima, one for each corner of the cell, which correspond to preferred orientations of the LSC. This was tested in the same apparatus for a single LSC in which there was an insulating middle wall which totally isolated the roll from the neighboring cell Ji and Brown (2020b).

To characterize interaction terms, we will look for differences from these baseline results from Eqs. 2 and 3.

IV Observations of preferred states

Before developing complicated models, we categorize the observed states of behavior to identify what solutions a good model should have. An example time series of the LSC orientation $\theta_{i}$, amplitude $\delta_{i}$, and mean temperature $T_{0,i}$ lasting 7 days is shown in Fig. 3 at Ra $=2.7\times 10^{9}$ ($\Delta T=19$ K). To show the long time series, only 1 in every 10 data points is shown in this figure. The system is seen to spontaneously switch between two distinct types of states. Before 160,000 s and after 370,000 s, there are fluctuations around a mean state with $\theta_{1}\approx\theta_{2}$. This corresponds to a counter-rotating state with the orientation vectors at $\theta_{i}$ aligned head-to-head, with flow in the same direction adjacent to the interface between the rolls, as illustrated at the top of Fig. 3a. Counter-rotating states are usually aligned along a diagonal, which is the preferred state for a single roll in a cubic cell, and in these examples is at $\pi/4$ rad before 160,000 s, and $7\pi/4$ rad after 370,000s. From 160,000-370,000 s, there are fluctuations around a mean state in which $\theta_{1}$ and $\theta_{2}$ instead differ by $\pi$ rad, with stable values of $\theta_{1}\approx 0$ and $\theta_{2}\approx\pi$ rad. This corresponds to a co-rotating state where both rolls are rotating in the same direction, with counter-flow adjacent to each other at the interface between the cells, as illustrated at the top of Fig. 3a.

The LSC amplitude $\delta_{i}$ also varies with the flow state, as shown in Fig. 3b. For comparison, the mean amplitude $\delta_{0}$ for a single, non-interacting roll is drawn as a dashed line Ji and Brown (2020b). In the counter-rotating state, the mean amplitude is 9% smaller than $\delta_{0}$, while in the co-rotating state, the mean amplitude is 38% smaller than $\delta_{0}$. This indicates the interaction between neighboring rolls reduces the temperature amplitude $\delta$, especially in the co-rotating state. While strong turbulent fluctuations in both states hide any detailed patterns that might appear within a state, the widths of the distributions of $\theta_{i}$ and $\delta_{i}$ around their mean values in a state are both narrower in the counter-rotating state, suggesting the stabilizing forces are stronger in the counter-rotating state. In the case of $\delta_{i}$, the non-linear damping term in Eq. 2 produces less of a stabilizing force when $\bar{\delta}_{i}$ is smaller, which may account for the smaller width of $\delta_{i}$ in the co-rotating state.

The LSC mean temperature $T_{0,i}$ also varies with the flow state, as shown in Fig. 3c. $T_{0,i}$ acquires a large systematic offset between the two cells in the co-rotating state, such that the cell whose colder side is at the interface ($\theta_{i}=\pi$ rad) acquires a higher mean temperature $T_{0,i}$. The fact that $T_{0,i}$ has any systematic change is notable, since for a single LSC, no patterns in $T_{0}$ were found or reported due to cessations Brown and Ahlers (2008a) or changes in preferred orientation Ji and Brown (2020b), except a small periodic modulation due to the jump-rope oscillation mode of the LSC Horn et al. (2022).

To identify the possible states over a wide parameter space, we carried out experiments at different values of control parameters, including a temperature difference $\Delta T$ between the top and bottom plates in the range $0.53$ K $\leq\Delta T\leq 24.4$ K ($7.6\times 10^{7}\leq$ Ra $\leq 3.5\times 10^{9}$), a difference $\Delta T_{m}$ between the mean plate temperatures of each cell as large as $\Delta T$, and tilt angles $\beta$ ranging from $0\leq|\beta|\leq 10^{\circ}$. Experiments typically lasted about 1 day, and switching between counter-rotating and co-rotating states occurred on average about once every other day. The dynamics of this switching will be discussed more thoroughly in Sec. VIII.

To characterize the preferred orientations of the states, a state diagram is made by plotting the values of the preferred orientation $\theta_{p,1}$ of cell 1 against the preferred orientation $\theta_{p,2}$ of cell 2 in Fig. 4. Preferred orientation were obtained from fits of the peaks of probability distributions of $\theta_{i}$. Experiments with switching between clear preferred orientations were divided up into separate datasets based on the preferred orientation to identify the preferred stable states. Counter-rotating states are defined by $\theta_{1,p}\approx\theta_{2,p}$ (dashed line in Fig. 4) such that they both flow in the same direction at the interface between the rolls. In the most symmetrically-driven cases, where $|\Delta T_{m}|/\delta_{0}<0.1$ and $|\beta|<0.03^{\circ}$ (solid circles in Fig. 4), usually both cells are aligned near a diagonal (illustrated as the red Xs in Fig. 4) which we refer to more specifically as counter-diagonal states. This matches the preferred diagonal state of a single cell Liu and Ecke (2009); Bai et al. (2016); Foroozani et al. (2017); Giannakis et al. (2018); Vasiliev et al. (2016, 2019); Ji and Brown (2020b). The alignment of the two rolls in a counter-rotating state with $\theta_{1,p}\approx\theta_{2,p}$ suggests there is a stable forcing in $\theta_{1}-\theta_{2}$.

However, in a few nominally symmetric cases we found counter-rotating states with preferred orientations roughly halfway between a diagonal and on-axis alignment. In some of these cases, there were switching events between the two corresponding symmetric preferred orientations on either side of $\theta=0$ or $\pi$ rad. This alignment suggests there is a component of force that aligns the flow of the two cells in opposite directions toward the axis going through $\theta=0$ and $\theta=\pi$ rad, which combines with the geometric forcing $-\nabla V_{g}(\theta_{i})$ to produce intermediate preferred orientations. Since this forcing likely comes from an unintended asymmetry of our setup, these counter-axis states will not be the focus of the work, although they are discussed in Sec. VII.7.

With a slight asymmetric forcing added to the basic counter-diagonal state with either $0.1<|\Delta T_{m}|/\delta_{0}<1.7$ or $0.03^{\circ}<|\beta|<2^{\circ}$, we find asymmetric offsets from the basic counter-rotating case where $\theta_{1,p}$ and $\theta_{2,p}$ shift in opposite directions away from the corner (open circles in Fig. 4), so we refer to these as counter-offset states.

Stable co-rotating states are always found to be aligned with one of $\theta_{1}$ and $\theta_{2}$ aligned nearly with 0, and the other nearly with $\pi$ rad (open squares in Fig. 4). Co-rotating states are clustered around these orientations regardless of the value of $\Delta T_{m}$ and $\beta$. While we are able to find co-rotating states some fraction of the time for any parameter values in our experiment, when a large enough asymmetry is introduced, either $|\Delta T_{m}|/\delta_{0}\stackrel{{\scriptstyle>}}{{{}_{\sim}}}1.7$ or $|\beta|\stackrel{{\scriptstyle>}}{{{}_{\sim}}}2^{\circ}$, we find only co-rotating states and no more counter-offset or other counter-rotating states. Since these preferred orientations $\theta_{p,i}=0$ or $\pi$ rad are along the axis between the center of both cells, this indicates the interaction between neighboring rolls adds a significant forcing stable in $\theta_{i}-\theta_{p,i}$, with $\theta_{p,i}=0$ or $\pi$ rad that can be dominant over the geometric potential $V_{g}(\theta_{i})$ to align both rolls along this axis in co-rotating states.

V Model for neighboring roll interaction: effective turbulent thermal diffusion

V.1 Forcing terms $\dot{\delta}_{\kappa,i}$ and $\ddot{\theta}_{\kappa,i}$ due to neighboring roll interactions

The model for the neighboring roll interaction terms is derived in detail in the appendix, and we summarize the physical mechanisms and resulting equations here. We start with the assumption that Eqs. 2 and 3 describe a baseline model for each LSC, and that we only need to derive terms for additional forcing due to neighboring-roll interactions. The forcing on the temperature profile along the interface between neighboring convection rolls is assumed to come from the turbulent diffusion of heat across the interface with effective turbulent thermal diffusivity $\kappa_{t}$, which is assumed to be uniform and constant. The turbulent thermal diffusion represents the enhancement of heat transport by eddies relative to the thermal diffusivity $\kappa$ due to advection, analogous to a turbulent viscosity or eddy viscosity. This formulation is mathematically analogous to boundary layer approximations with thermal diffusion in an interfacial mixing layer. We include the unknown mixing layer thickness in the value of $\kappa_{t}$, which will be a fit parameter, so we can use the known cell size $L$ as the lengthscale in the thermal diffusion equation. Because of the middle insulating wall blocking half of the interface, the heat transport acts over the exposed fractional area $A=0.50$ of the interface between the two cells. The model terms are calculated by applying this turbulent thermal diffusion in a $\theta$-dependent heat equation for heat flux in the direction perpendicular to the interface between neighboring cells. The resulting forcing terms (as calculated in appendix A) are:

[TABLE]

and

[TABLE]

For brevity, we write all equations in this section for the forcing on cell 1 only, as the equations for cell 2 are identical other than an exchange of the subscripts 1 and 2 in each equation. These forcings $\dot{\delta}_{\kappa,1}$ and $\ddot{\theta}_{\kappa,i}$ can be inserted directly to the existing stochastic equations of motion for a single LSC (Eqs. 2 and 3, respectively).

V.2 Equation of motion for $T_{0,i}$

The model for a single LSC did not include an equation of motion for the mean temperature $T_{0,i}$ because it was observed to have trivial behavior. With neighboring convection cells, we observe a systematic difference between $T_{0,2}$ and $T_{0,1}$ (Fig. 3c), and this difference may drive dynamics in $\theta_{i}$ and $\delta_{i}$ according to Eqs. 5 and 4. Thus, we need an equation of motion for $T_{0,i}$.

In an equation of motion for $T_{0,i}$, we should also expect to have some diffusive fluctuations driven by turbulence, analogous to Eqs. 2 and 3, so we include a fluctuation term $f_{T}(t)$ with diffusivity $D_{T}$. The deterministic part of the equation for $T_{0,1}$ is assumed to be due to the same net heat flux between the cells that contributes to the forcings $\dot{\delta}_{\kappa,i}$ and $\ddot{\theta}_{\kappa,i}$. An additional vertical heat transport from each of the top and bottom plates to each LSC is calculated using the standard boundary layer approximation in terms of the Nusselt number Nu. This results in (as calculated in Appendix B):

[TABLE]

While this equation could be plugged into Eqs. 4 and 5, we find it more insightful to leave Eqs. 4 and 5 as is for a more direct physical interpretation of their dependence on $T_{0,2}-T_{0,1}$ at the stable fixed points where $\dot{T}_{0,i}=0$.

VI Testing the model for symmetric cases ($\Delta T_{m}=0$)

VI.1 Forcing terms on $\ddot{\theta}_{i}$ in $\theta_{1}-\theta_{2}$ and $\theta_{i}-\theta_{p}$

In Sec. IV, the observed counter- and co-rotating states suggest forcing terms on $\ddot{\theta}_{i}$ that are stable in $\theta_{1}-\theta_{2}$ and $\theta_{i}-\theta_{i,p}$, respectively. In this section, we test the predictions for these forcing terms from Eq. 5 for counter- and co-rotating states, starting with the simpler special case of $\Delta T_{m}=0$.

VI.1.1 Stable fixed points and linearized forcing

Qualitatively, Eq. 5 has stable fixed points that are linearly stable for both the observed counter-rotating states ($\theta_{1}=\theta_{2})$, and co-rotating states ($\theta_{1}=0$ and $\theta_{2}=\pi$ rad or $\theta_{1}=\pi$ rad and $\theta_{2}=0$), as shown in Appendix C and D, respectively, confirming that these states are predicted by the model.

In the counter-rotating state, a linear expansion of Eq. 5 around the stable fixed point at a preferred orientation of the cubic cell (e.g. $\theta_{1}=\pi/4$ rad) and assuming $\bar{\delta}_{i}=\delta_{0}$ is

[TABLE]

In the co-rotating states, assuming constant $\Delta\bar{T}_{0,co}\equiv\bar{T}_{0,2}-\bar{T}_{0,1}$ and $\bar{\delta}_{co}\equiv\bar{\delta}_{1}=\bar{\delta}_{2}$, a linear approximation around $\theta_{p,1}=0$ or $\pi$ rad for co-rotating states of Eq. 5 is

[TABLE]

VI.1.2 Multiple linear regression

To test the functional form of the net forcing around stable fixed points in the Eqs. 7, and 8 added to the stochastic equation of motion (Eq. 3), we carried out a multiple linear regression of the equation

[TABLE]

with stiffnesses $a$ and $b$. We insert the sine functions as the simplest way to account for periodicity of the coordinate $\theta$ while remaining linear in the lowest order expansion around the stable fixed points.

To be able to carry out an accurate linear regression, we first precisely determine the peak locations $\theta_{p,i}$ by fitting the probability distribution $p(\theta_{i})$ for each state to a Gaussian function, shown in Fig. 5 for the same data as Fig. 3. To calculate the forcing $\ddot{\theta}_{i}$, we calculate the 2nd derivative from discrete data as $\ddot{\theta}_{i}=[\theta_{i}(t+\Delta t)-2\theta_{i}(t)+\theta_{i}(t-1)]/(\Delta t)^{2}$. We apply a multiple linear regression of Eq. 9 to measured values of $\ddot{\theta}_{i}$ to obtain the coefficients $a$ and $b$, using only data with $\delta>0.5\bar{\delta}$ to avoid bias during cessations where the acceleration can be much larger than typical due to weaker damping in Eq. 3 Brown and Ahlers (2008a).

From the linear regression of Eq. 9, we generally find $a$ to be negative in the counter-rotating and co-rotating states, confirming that the system is stable against displacements in $\theta_{i}-\theta_{p,i}$, as predicted due to the geometric potential $V_{g}$ Ji and Brown (2020b) and Eq. 8, respectively. We find $b$ negative for cell 1 and positive for cell 2 in counter-rotating states, confirming they are stable in $\sin(\theta_{1}-\theta_{2})$ as predicted in Eq. 7. We also find $b$ positive for co-rotating states in cell 1 and negative in cell 2, which – due to the flip in the sign of the sine function with a phase shift of $\pi$ rad – means that co-rotating states are also stable in $\sin(\theta_{1}-\theta_{2})$, corresponding to a forcing to align the orientation vectors head-to-tail in the co-rotating state, but opposite the sign predicted in Eq. 8.

If we calculate the forcing in terms of $\dot{\theta}_{i}(\theta_{1}-\theta_{2})$ instead of $\ddot{\theta}_{i}(\theta_{1}-\theta_{2})$, similar qualitative results are found, but the strongest forcing is found when $\dot{\theta}_{i}$ is calculated with about a 20 s delay after the time that $\theta_{1}-\theta_{2}$ is recorded. This is comparable to the damping timescale $\tau_{\dot{\theta}}$ that corresponds to the ratio of these two terms in Eq. 3. This confirms that a forcing in terms of $\ddot{\theta}_{i}$ better describes the effects of a neighboring roll interaction, justifying the conversion of $\dot{\theta}_{\kappa,i}$ into $\ddot{\theta}_{\kappa,i}$ in Eq. 25.

VI.1.3 Functional form of forcing $\ddot{\theta}_{i}(\theta_{i}-\theta_{p,i})$

In co-rotating states and counter-diagonal states, we typically find the stiffness $a$ much larger than $b$ from the regression, which indicates that the forcing in $a\sin(\theta_{i}-\theta_{p,i})$ usually determines the overall stability of the state. For example, for the data in Fig. 3, $a=10b$ in the counter-rotating state, and $a=2.6b$ in the co-rotating state 111For counter-offset states with $\Delta T_{m}/\delta_{0}\stackrel{{\scriptstyle>}}{{\sim}}0.1$, $\theta_{p,1}$ and $\theta_{p,2}$ shift away from the corner, so that the forcing is not necessarily centered on the peak of $p(\theta_{i})$. In these cases the apparent value of $b$ drops significantly even for small $\Delta T_{m}$. We attempted to add an extra constant offset term to be fit in the linear regression to fix this. While the values of $b$ do not systematically drop with small $\Delta T_{m}$ with a constant offset in the linear regression, the extra regression parameter causes uncertainties to increase to be comparable to fit parameters. We do not report results for $b$ with $\Delta T_{m}\neq 0$, and caution about this limitation of the algorithm for states where $\theta_{p,1}$ is offset from $\theta_{p,2}$.. The dominance of the forcing $a\sin(\theta_{i}-\theta_{p,i})$ provides an opportunity to analyze the functional form of the forcing terms separately. Since $a\sin(\theta_{i}-\theta_{p,i})$ is much larger, then the probability distribution $p(\theta_{i})$ shown in Fig. 5 is determined mostly by this term, so we can first obtain an approximate forcing in $\theta_{i}-\theta_{p,i}$ while disregarding the forcing in $\theta_{1}-\theta_{2}$. The forcing in $\theta_{i}-\theta_{p,i}$ can be related to a probability distribution of $\theta_{i}-\theta_{p,i}$ in the overdamped limit of the stochastic Eq. 3 (i.e. $\ddot{\theta}_{i}$ is small compared to the damping and forcing terms) and if $\delta$ is nearly a constant by a Fokker-Planck equation Ji and Brown (2020b)

[TABLE]

where $V(\theta_{i})=V_{g}(\theta_{i})-\int\ddot{\theta}_{\kappa,i}d\theta_{i}$ is the net potential from geometric and neighboring roll interaction forces. A linear forcing in $\theta_{i}-\theta_{p,i}$ corresponds to a quadratic potential $V(\theta_{i})$, and a Gaussian probability distribution. Each probability distribution $p(\theta_{i})$ is fit to a Gaussian function, as shown in Fig. 5. Errors are shown in Fig. 5 calculated assuming Poisson statistics. Errors on $\theta_{i}$ of 0.18 rad from the fit of the temperature profile (Eq. 1) divided by the square root of the number of counts are too small to see in Fig. 5, but remain significant in the fitting of the data. Fitting data up to 2.7 standard deviations away from the peak yields reduced $\chi^{2}$ of 1.4,1.3, 1.2, and 0.9 for the four fits shown in Fig. 5, confirming the data are consistent with Gaussian probability distributions in this range. The Gaussian shape of $p(\theta_{i})$ confirms the forcing in $\theta_{i}-\theta_{p,i}$ is linear around the stable fixed point, confirming the assumption made in Eq. 9, and matching the prediction of Eq. 8 for co-rotating states, as well as the predicted linear forcing for counter-diagonal states due to the geometric potential $V_{g}(\theta_{i})$, which has been previously observed to be linear for a single LSC Ji and Brown (2020b).

VI.1.4 Functional form of forcing $\ddot{\theta}_{i}(\theta_{1}-\theta_{2})$

Since the multiple linear regression in Sec. VI.1.2 showed that the forcing in $\theta_{1}-\theta_{2}$ is smaller than the forcing in $\theta_{i}-\theta_{p,i}$, then a probability distribution of $\theta_{1}-\theta_{2}$ will not reveal the forcing in $\theta_{1}-\theta_{2}$ because it will be overwhelmed by the stronger forcing in $\theta_{i}-\theta_{p,i}$. Instead we calculate a corrected forcing $\ddot{\theta}_{i,c}=\ddot{\theta}_{i}(\theta_{1}-\theta_{2})-a\sin(\theta_{i}-\theta_{p,i})$ as a function of $\theta_{1}-\theta_{2}$. This method allows us to correct for the linear dependence on $\theta_{i}-\theta_{p,i}$ that was confirmed by Fig. 5, so that only $b\sin(\theta_{1}-\theta_{2})$ is expected to remain, based on Eq. 9. To do this, we calculate the azimuthal acceleration rate from discrete data as the 2nd order difference $\ddot{\theta}_{i}(t)=[\theta_{i}(t+1)-2\theta_{i}(t)+\theta_{i}(t-1)]/(\Delta t)^{2}$. The forcing we report with a measurement timestep of $\Delta t=7$ s is consistently about 20% smaller than runs with a shorter timestep of 2 s, due to the approximation of a second derivative with a non-zero timestep. For each data point, we subtract $a\sin(\theta_{i}-\theta_{p,i})$ using the value of $a$ from the regression (Sec. VI.1.2). We exclude data with $\delta<0.5\bar{\delta}$ to avoid biasing data by cessations where the acceleration can be much larger than typical values Brown and Ahlers (2008a). We then bin values of $\theta_{i,c}$ over small ranges of $\theta_{1}-\theta_{2}$ to calculate forcing as a function of $\theta_{1}-\theta_{2}$ and average $\langle\ddot{\theta}_{i,c}(\theta_{1}-\theta_{2})\rangle$ in each bin to reduce the contribution of the large stochastic fluctuations in Eq. 3. The error on the average in each bin is reported as the standard deviation of the mean, assuming the data points are independent.

This corrected forcing $\langle\ddot{\theta}_{i,c}(\theta_{1}-\theta_{2})\rangle$ is shown in Fig. 6 for the same dataset as in Fig. 3. Data from the counter-rotating state are shown in panel a, and data from the co-rotating state are shown in panel b. In both states, cell 1 has a negative slope in $\theta_{1}-\theta_{2}$, and cell 2 has a positive slope in $\theta_{1}-\theta_{2}$, confirming they are stable at their intercepts. The intercepts are near $\theta_{1}-\theta_{2}=0$ for the counter-rotating state, and $\theta_{1}-\theta_{2}=\pi$ rad for the co-rotating state, corresponding to their stable relative orientations (Fig. 4).

The corrected forcing is fit to $\langle\ddot{\theta}_{i,c}(\theta_{1}-\theta_{2})\rangle=b\sin(\theta_{1}-\theta_{2})$ where the stiffness $b$ is a fit parameter. Only bins with at least 10 data points are included in the fit. The reduced $\chi^{2}$ values are 1.2 and 1.2 for the counter-rotating data sets, 1.3 and 2.4 for the co-rotating data, indicating the sine function is consistent with the data for the counter-rotating states, but not necessarily for the co-rotating states. The magnitudes of the values of $b$ are consistent with the values obtained when using the linear regression analysis within a couple of standard deviations of the mean, which is on average 16% of the mean. This self-consistency in the magnitudes obtained from both methods confirms the earlier assumptions that the net forcing on the LSC orientation can be represented by Eq. 9, consistent with the functional forms predicted by $V_{g}$ and Eqs. 7 for counter-rotating states, and the terms are independent enough to use the conditional average in Fig. 6. The magnitude of $b$ for counter-rotating states will be used to obtain $\kappa_{t}$ and compared to other model terms in Sec. VI.3. However, the negative sign of the stiffness $b$ for co-rotating states is inconsistent with the prediction of Eq. 8, which predicted $b>0$ (destabilizing) for co-rotating states. Since $a>b$, then $b$ is not responsible for the stability of co-rotating states, but this remains a minor disagreement with the model Eq. 5.

VI.1.5 Comparison of forcing term $\ddot{\theta}_{i}(\theta_{i}-\theta_{p,i})$ to the geometric forcing for the counter-diagonal state

The forcing $-\nabla V_{g}$ due to the cell geometry in Eq. 3 has been observed and predicted to be linear and stable in an expansion of $\theta_{i}-\theta_{p,i}$ around the potential minima along the diagonals Ji and Brown (2020b), which is found to be the preferred orientation for most counter-rotating states. Thus, it is predicted to be responsible for the measured value of the stiffness $a$ for counter-rotating states. To determine how much of the measured value of $a$ comes from the geometric potential, we calculated the forcing $\langle\ddot{\theta}_{i}(\theta_{i})\rangle$ by averaging values of $\ddot{\theta}_{i}$ in bins of $\theta_{i}$ as in Fig. 6 in Sec. VI.1.4 from an experiment with a single LSC in the same apparatus, but with a insulating middle wall that completely blocks the two cells off from each other (i.e. $A=0$) at Ra $=2.7\times 10^{9}$ Ji and Brown (2020b). For a single LSC, we find the stiffness $a=-1.02\pm 0.11$ mrad/s2, and for neighboring counter-rotating rolls, we find $a=-1.11\pm 0.17$ mrad/s2. The consistency of these values suggest that the stiffness $a$ comes entirely from the geometric forcing $-\nabla V_{g}$ for counter-rotating states, and there is no significant neighboring roll interaction contribution to $a$ for counter-rotating states, in agreement with the prediction of Eq. 7.

For co-rotating states, the stiffness $a$ is expected to come from a competition between the second and third terms of Eq. 5, and be reduced by the geometric forcing $-\nabla V_{g}$ around its unstable fixed point. Since the source of the stability of co-rotating states is a more fundamental question and involves more parameters, we put this off until Sec. VII.5 after $\kappa_{t}$ is fit and more terms of the model are tested.

VI.2 Diffusive fluctuations of $T_{0}$

While the stochastic equations of motion for the LSC orientation $\theta_{0}$ (Eq. 3) and temperature amplitude $\delta$ (Eq. 2) have been tested previously for single convection rolls Brown and Ahlers (2008a); Ji and Brown (2020b), the equation of motion for the mean temperature (Eq. 6) is a entirely new. To test whether the stochastic term representing turbulent fluctuations is diffusive, we measure the mean-square displacement $(dT_{0,i})^{2}$ of $T_{0,i}$ over different time intervals $dt$. The mean-square displacement $(dT_{0,i})^{2}$ is averaged over the two cells and different starting times, and plotted as a function of time interval $dt$ in Fig. 7 for a counter-rotating state at Ra $=2.7\times 10^{9}$. The apparent linearity in the limit of short time intervals is an indicator of diffusive fluctuations, as we assumed in the model Eq. 6. A linear function $\langle(dT_{0,i})^{2}\rangle=D_{T}dt$ is fit to the data for $dt<35$ s to obtain a diffusivity $D_{T}$. The time scale of 35 s is approximately the damping timescale $\tau_{\delta}$, which is likely relevant because $\delta_{i}$ appears in Eq. 6, so fluctuations in $\delta_{i}$ may affect changes in $T_{0,i}$ at larger $dt$.

The inset of Fig. 7 shows the mean-square displacement $(dT_{0,i})^{2}$ over a larger range of $dt$. At large $dt$, the mean-square displacement reaches a plateau $(dT_{0,i})^{2}\equiv\tau_{T}D_{T}$ which defines a damping timescale $\tau_{T}$ as the time of intersection of the two limiting scaling laws.

VI.3 The effective turbulent thermal diffusivity $\kappa_{t}$

Since we introduced the effective turbulent thermal diffusivity $\kappa_{t}$ as an unknown fit parameter, we do not have first-principles predictions of model terms. On the other hand, since nearly all of the new deterministic forcing terms in Eqs. 4, 5, and 6 are proportional to $\kappa_{t}$, we can check the consistency of values of $\kappa_{t}$ required for different forcing terms to fit with data.

VI.3.1 $\kappa_{t}$ from $b$

As one measure of $\kappa_{t}$, we compare the measured stabilizing force in $b\sin(\theta_{1}-\theta_{2})$ from Eq. 9 to the predicted linearized forcing in Eq. 7 to obtain $\kappa_{t}/\kappa=4bL^{2}\tau_{\dot{\theta}}/A\kappa$. We use measurements of $b$ from Sec. VI.1.4 for counter-diagonal states with small $\Delta T_{m}/\bar{\delta}<0.1$, and measurements of $\tau_{\dot{\theta}}$ from the crossover time of limiting scalings of of the mean-square displacement of $\dot{\theta}_{i}$ (analogous to $\tau_{T}$) obtained in previous measurements of a single LSC Ji and Brown (2020b). Counter-axis states, which are found at $\Delta T\leq 2$ K are not included in the analysis because they usually have multiple peaks in the probability distribution of $\theta_{i}$, making the Gaussian fitting and linear regression algorithm to obtain $b$ unrepresentative of the forcing. Errors are propagated from errors on $b$ obtained from the linear regression. These values of $\kappa_{t}/\kappa$ and shown as a function of Ra in Fig. 8. Unfortunately, the uncertainty on $b$ makes the uncertainties too large to draw strong conclusions about the trend over a small range of Ra from this data alone.

VI.3.2 $\kappa_{t}$ from the damping timescale $\tau_{T}$

For another measure of $\kappa_{t}$, we use the damping timescale $\tau_{T}$ for fluctuations in $T_{0}$ seen in Fig. 7. The coefficient $A\kappa_{t}/L^{2}$ of the second term of Eq. 6 corresponds to an inverse damping timescale $1/\tau_{T}$ in the limit where $4\mbox{Nu}\kappa/A\kappa_{t}\ll 1$ (confirmed in Sec. VII), and assuming the $\delta_{i}$ terms do not contribute significantly in the counter-rotating state (the $\delta_{i}$ terms cancel each other out at the stable fixed point). This relation leads to $\kappa_{t}/\kappa=L^{2}/A\kappa\tau_{T}$, using the measured values of $\tau_{T}$ from the intersection of the two limiting scaling laws in Fig. 7. The corresponding values of $\kappa_{t}/\kappa$ vs. Ra are plotted in Fig. 8 for counter-rotating states. Errors are propagated from fits of $(\Delta T_{0,i})^{2}$, which were calculated as half the difference between two cells, plus the error from fits of $D_{T}$, plus half the difference over a fit range from 0 to $0.6\tau_{\delta}$ or $0.8\tau_{\delta}$. This resulted in an average error of 12%. A power law fit to the data obtained from $\tau_{T}$ yields $\kappa_{t}/\kappa=(0.8\pm 0.3)Ra^{0.41\pm 0.02}$ with a reduced $\chi^{2}=1.5$. These measurements of $\kappa_{t}/\kappa$ based on the damping timescale $\tau_{T}$ are consistent with the data for $\kappa_{t}/\kappa$ based on measurements of the stiffness $b$, confirming that the same value of $\kappa_{t}$ drives both the $T_{0,i}$-term in Eq. 6 and the $\theta_{1}-\theta_{2}$-term in Eq. 5.

In principle we could plot the forcing for other terms of Eqs. 4, 5, and 6. However, other measured factors depend on a small difference between two terms. For example, Eq. 4 has a difference in terms depending on $\delta_{2}$ and $T_{0,2}-T_{0,1}$, so that uncertainties on the data and model can heavily bias scaling laws from these small differences. The values of $\kappa_{t}$ from measurements of $\tau_{T}$ remain our most precise method of measuring $\kappa_{t}$. Thus, we will use this value of $\kappa_{t}$ obtained from measurements of $\tau_{T}$ for predictions of $\kappa_{t}$ in all other model terms.

VI.4 Stable fixed points of temperature amplitude $\bar{\delta}$ and mean temperature difference $\Delta\bar{T}_{0}$

VI.4.1 Predictions for $\bar{\delta}_{co}$ and $\Delta\bar{T}_{0,co}$ in co-rotating states

For co-rotating states where $\theta_{1}=0$ and $\theta_{2}=\pi$ ra for $\theta_{1}=\pi$ rad and $\theta_{2}=0$, the symmetry of Eqs. 4 and 6 still requires $\delta_{1}=\delta_{2}$ at stable fixed points for co-rotating states, but they do not require $T_{0,1}=T_{0,2}$. For the symmetric driving case of $T_{m,2}=T_{m,1}$, in the Boussinesq limit where the mean temperature of the bulk equals the mean temperature of the plates, requires $(T_{0,2}+T_{0,1})/2=T_{m,i}$. Defining the fixed point value of $\delta_{i}$ in the co-rotating state as $\bar{\delta}_{co}$, Eq. 6 and its corresponding equation for $\dot{T}_{2}$ have stable fixed points ($\dot{T}_{0,i}=0$) when

[TABLE]

where the $+$ sign corresponds to the co-rotating state orientation with $\theta_{1}=0$ and $\theta_{2}=\pi$ rad, and the $-$ sign corresponds to swapped orientations. To obtain the the stable fixed point $\bar{\delta}_{co}$ in co-rotating states, we evaluate Eq. 4 at $\theta_{1}=0$ and $\theta_{2}=\pi$ rad when $\Delta\bar{T}_{0,co}>0$ or $\theta_{1}=\pi$ rad and $\theta_{2}=0$ when $\Delta\bar{T}_{0,co}<0$. This simplifies Eq. 4 to $(\bar{\dot{\delta}}_{\kappa,i})_{co}=(A\kappa_{t}/\pi L^{2})\left(1.41|\Delta\bar{T}_{0,co}|-2.57\delta_{i}\right)$. The stable fixed point value of $\bar{\delta}_{co}$ is obtained by using this for $\dot{\delta}_{\kappa,i}$ in Eq. 2, setting $\dot{\delta}_{i}=0$, and linearly expanding around $\delta_{i}=\delta_{0}$ to obtain

[TABLE]

Note that a simultaneous evaluation of Eq. 11 and Eq. 12 would lead to a subtraction of two comparable numbers when $2\mbox{Nu}\kappa/A\kappa_{t}$ is small, which can result in predictions of differing signs and arbitrarily small magnitudes when including uncertainties from assumptions in the model, so we cannot make accurate predictions of $\bar{\delta}_{co}$. Instead, we will use the measured relationships between $\Delta\bar{T}_{0,co}$ and $\bar{\delta}_{co}$ to text for self-consistency of Eqs. 11 and 12.

VI.4.2 Measurements of $\bar{\delta}$ and $\Delta\bar{T}_{0}$

To characterize the change in the mean LSC amplitude $\bar{\delta}$ due to the interaction between neighboring rolls, we plot values normalized by $\delta_{0}$ in Fig. 9 as a function of Ra. The average $\bar{\delta}$ is calculated as the average of $\delta_{i}$ over both cells and over time while the system remained in one state. The normalization values of $\delta_{0}$ are obtained from a fit of a power law to $\delta_{0}(Ra)$ from previous experiments with a single LSC in a cubic cell Ji and Brown (2020b).

Counter-rotating states consistently have a mean amplitude $\bar{\delta}$ close to that of the single-cell, with $\bar{\delta}_{counter}/\delta_{0}=0.94$ with a standard deviation of 0.03, only slightly lower than the predicted value of 1 from Eq. 4 when $\theta_{1}=\theta_{2}$. This applies to both counter-diagonal and counter-axis states, so it is independent of the preferred orientation of the counter-rotating state. The slight decrease in $\bar{\delta}$ from $\delta_{0}$ may be partly due to fluctuations in $\theta_{1}-\theta_{2}$ away from the mean of zero, which can make the third term of Eq. 4 negative. Co-rotating states have a smaller mean amplitude with a consistent ratio $\bar{\delta}_{co}/\delta_{0}=0.57$ with a standard deviation of 0.04.

We measured $\Delta\bar{T}_{0}\equiv\bar{T}_{0,2}-\bar{T}_{0,1}$ as an average over time while the system remained in one state. Data are normalized by $\bar{\delta}$ for a sense of scale. Counter-rotating states are scattered around zero, but the values are much larger than uncertainties on temperature measurements, which will be explained due to small differences in $\Delta T_{m}$ in Sec. VII. Co-rotating states have a consistent average $|\Delta\bar{T}_{0,co}|/\bar{\delta}_{co}=2.22\pm 0.13$ (the $\pm$ refers to one standard deviation), 30% larger than predicted due to effective turbulent thermal diffusion from Eq. 11 assuming $2B/A\ll 1$. This is well within the typical errors of this modeling approach of a factor of 3 Brown and Ahlers (2008a).

Although the measured ratio of $|\Delta\bar{T}_{0,co}|/\bar{\delta}_{co}$ is close to the prediction, the value $|\Delta\bar{T}_{0,co}|/\bar{\delta}_{co}>2$ has some significance. It indicates that the cold cell is colder than the hot cell along the entire interface, as well as being colder than the mean of the plates, so turbulent thermal diffusion is expected to provide a net heat flux into the cold cell at this value of $|\Delta\bar{T}_{0,co}|/\bar{\delta}_{co}$, regardless of assumptions about how the shape of the temperature profile affects the model terms. Thus, this magnitude $|\Delta\bar{T}_{0,co}|/\bar{\delta}_{co}>2$ is a violation of the assumption in Eq. 20 that $T_{0,i}$ is driven by turbulent thermal diffusion only. An additional heat transport mechanism is required, which is likely advection of fluid in a coherent flow from one cell to the other. For example, a coherent advective flow parallel to the plates has been observed in other co-rotating states Petry and Busse (2003); Pérez-García et al. (2004); Xie and Xia (2013); Gao et al. (2013); Mercader et al. (2019). This coherent advection can transport heat from the hot plate of the cold cell directly to the hot cell and from the cold plate of the hot cell directly to the cold cell to balance the diffusive heat transport between cells.

While we could not accurately predict the value of $\bar{\delta}_{co}$ from Eq. 12, we can provide a self-consistency test of the model by using the measured ratio $|\Delta\bar{T}_{0,co}|/\bar{\delta}_{co}$ as empirical input into Eq. 12. Since the product of $\kappa_{t}$ and $\tau_{\delta}$ has a mild dependence with $Ra^{-0.17}$, we estimate a value of $\bar{\delta}_{co}$ at the middle of the range of Ra. Using empirical input into Eq. 12 of $|\Delta\bar{T}_{0,co}|/\bar{\delta}_{co}=2.22\pm 0.13$ from Fig. 10, $\bar{\delta}_{co}=0.57\delta_{0}$ from Fig. 9, $\kappa_{t}/\kappa=3200\pm 400$ from the average and standard deviation of five measurements at Ra $=5.5\times 10^{8}$ in Fig. VI.3, and $\tau_{\delta}=81$ s from a fit of measurements of a single LSC in a cubic cell Ji and Brown (2020b), results in $\bar{\delta}_{co}=(0.83\pm 0.07)\delta_{0}$. This underestimates the reduction in $\bar{\delta}$ in the co-rotating state by 60% with input of $\Delta\bar{T}_{0,co}$. The variation of this prediction at the extremes of the measured range of Ra results in a difference of $0.05\delta_{0}$.

VII Testing the model with a difference in mean temperature of the plates of the two cells $\Delta T_{m}$

In this section we extend the predictions of Eqs. 4, 5, and 6 and test them in cases where the difference in the mean temperatures of the plates of the two cells $\Delta T_{m}\neq 0$. Controlling $\Delta T_{m}$ allows testing several terms of these equations which depend on $\Delta\bar{T}_{0}$, explains how some asymmetries in counter-rotating states come about due to a small unintentional $\Delta T_{m}$ in experiments, and helps identify the conditions for stable counter- or co-rotating states.

Figure 11 shows measurements of the stable fixed points $\theta_{p,i}$, $\bar{\delta}\equiv(\delta_{1}+\delta_{2})/2$, and $\Delta\bar{T}_{0}$ for different values of $\Delta T_{m}$ at Ra=$5.5\times 10^{8}$ ($\Delta T=3.8$ K). A smaller value of Ra was chosen for this set of experiments because switching between states is more frequent at this Ra Bai et al. (2016), making the sampling of different states easier. A counter-rotating state is found for small $\Delta T_{m}$, with $\theta_{p}=\pi/4$ rad in both cells at $\Delta T_{m}=0.04$ K. This offset from zero suggests a slight asymmetry in the cells resulting in a horizontal temperature difference that is 1% of $\Delta T$. For small $\Delta T_{m}$, $\theta_{p}$ varies with $\Delta T_{m}$ in opposite directions in the two cells, corresponding to counter-offset states. $\bar{\delta}$ and the magnitude of $\Delta\bar{T}_{0}$ also increase with $|\Delta T_{m}|$ for these counter-offset states. At larger $|\Delta T_{m}|$, only co-rotating states are found.

The different states in Fig. 11 can be understood qualitatively from Eqs. 5 and 6. When $\Delta T_{m}>0$ and thus $\Delta\bar{T}_{0}>0$ from Eq. 6, the interface between cells is hotter than cell 1, so there will be a forcing on $\theta_{1}$ towards the interface with the hotter cell according to Eq. 5. $\theta_{2}$ is forced in the opposite direction as the interface is colder than cell 2. In the counter-rotating state, this results in the opposite shifts in $\theta_{p,i}$ with increasing $|\Delta T_{m}|$. When the forcing in Eq. 5 becomes large enough as $|\Delta T_{m}|$ and $|\Delta\bar{T}_{0}|$ increase, it pushes $\theta_{p}$ across an unstable fixed point of Eq. 5 at $\theta_{p}=\pi/2$ rad so that only a co-rotating state is stable. In agreement with this prediction, the measured values of $\theta_{p}$ for counter-offset states in Fig. 11a approach but do not cross $\pi/2$ rad, and beyond that value of $\Delta T_{m}$ only co-rotating states are found.

In principle we could solve the coupled Eqs. 4, 5, and 6, to make direct predictions of how $\theta_{p}$ and $\Delta\bar{T}_{0}$ depend directly on $\Delta T_{m}$. However, the predictions end up being due to small differences between terms, which are sensitive to uncertainties in parameter values. Since the parameters have uncertainties as large as a factor of 3, resulting model predictions can have much larger uncertainties and even an uncertain sign. Instead, we evaluate stable fixed point values for one parameter at a time, using other parameters as input, to test the direct dependences of the model and check for self-consistency.

VII.1 Offset of the preferred orientations $\Delta\theta_{p}$ in counter-offset states

In counter-rotating states when $|\Delta\bar{T}_{0}|>0$, the third term of Eq. 5 gives a forcing on $\ddot{\theta}_{1}$ of $-\sqrt{2}A\kappa_{t}\Delta\bar{T}_{0}\sin\theta_{1}/\pi L^{2}\bar{\delta}\tau_{\dot{\theta}}$, which pushes the orientation $\theta_{p}$ of the LSC in the colder cell towards alignment with the hotter cell, and pushes $\theta_{p}$ in the hotter cell away from the colder cell, so that the $\theta_{p}$ is pushed equally in opposite directions in the two cells. We define the shift in preferred orientations $\Delta\theta_{p}\equiv(\theta_{p,2}-\theta_{p,1})/2$. Treating this as a small perturbation away from a counter-diagonal state where $\sin\theta_{1}\approx(\theta_{1}-\theta_{p})/\sqrt{2}$, $\Delta\theta_{p}$ is obtained by balancing this forcing from Eq. 5 with the overdamped forcing obtained in a single cubic cell $-\nabla V_{g}=-\omega_{r}^{2}(\theta_{1}-\theta_{p})$, where $\omega_{r}$ is the natural frequency of oscillation in the geometric potential around the corner Ji and Brown (2020b). At lowest order, this leads to a shift in $\theta_{p}$ of

[TABLE]

The $\pm$ sign corresponds to the sign of $\sin\theta_{p}$.

To test the prediction of $\Delta\theta_{p}$ from Eq. 13, we use as input $\tau_{\dot{\theta}}=22$ s , and $\omega_{r}^{2}=(2.1\pm 0.4)\times 10^{-4}$ s*-2* from the linear regression of Sec. VI.1.2, both from data at Ra=$5.4\times 10^{8}$ for a single LSC in a cubic cell Ji and Brown (2020b) 222This value of $\omega_{r}^{2}$ is about half of the value reported in Ji and Brown (2020b), which was obtained less directly through measurements of the mean-square displacement of $\dot{\theta}_{0}$ divided by the variance of the probability distribution of $\theta_{0}$. We also use as input the measured value $\bar{\delta}=0.12$ K, a linear fit of $\Delta\bar{T}_{0}=(0.57\pm 0.05)\Delta T_{m}$ from Fig. 11c, and $\kappa_{t}/\kappa=3200\pm 400$ at Ra $=5.5\times 10^{8}$ from Fig. 8. This yields the prediction $\Delta\theta_{p}=(1.9\pm 0.2)\Delta T_{m}$ rad/K. A linear fit to $\theta_{p}$ in Fig. 11a yields $\Delta\theta_{p}=(3.7\pm 0.4)\Delta T_{m}$ rad/K, which is within a factor of 2 of the prediction of Eq. 13. This discrepancy also means Eq. 13 overpredicts the maximum $\Delta\bar{T}_{0}$ where counter-offset states are stable by about the same factor, by extrapolating $\Delta\theta_{p}$ to where it crosses the unstable fixed point of Eq. 5 at $\theta_{i}=\pi/2$ rad.

VII.2 Increase of steady-state temperature difference $\Delta\bar{T}_{0,co}$ with $\Delta T_{m}$ in co-rotating states

In the co-rotating state, Eq. 6 has a stable fixed point solution

[TABLE]

The $\pm$ corresponds to the sign of $\Delta T_{m}$. Values of Nu are known from fits of the Grossmann-Lohse scaling model Ahlers et al. (2009). Specifically, at $Ra=5.5\times 10^{8}$, Nu=55 in a cylindrical cell at the same Ra Funfschiling et al. (2005) (values of Nu in cubic and cylindrical cells are found to agree within 2% in this range Kaczorowski and Xia (2013)). This corresponds to $\mbox{Nu}\kappa/\kappa_{t}=0.017\pm 0.002$. Using this value, and the fit slope $|\partial\bar{\delta}_{co}/\partial\Delta T_{m}|=(0.042\pm 0.006)\Delta T_{m}$ for $\Delta T_{m}<-0.02$ K, yields a prediction of $|\partial\Delta\bar{T}_{0,co}/\partial\Delta T_{m}|=0.14\pm 0.02$, about a factor of 2 larger than the measured slope $|\partial\Delta\bar{T}_{0,co}/\partial\Delta T_{m}|=(0.062\pm 0.014)$. However, we note that the trends start to change for values of $\Delta T_{m}>0.2$, such that the discrepancy becomes a factor of 5 in the range $0.2<\Delta T_{m}<1$.

VII.3 Increase of $\Delta\bar{T}_{0}$ with $\Delta T_{m}$ in counter-rotating states

While ideal symmetric counter-rotating states correspond to $\Delta T_{m}=0$, in practice the LSC is very sensitive to asymmetries Krishnamurti and Howard (1981); Brown and Ahlers (2006b). Even in nominally symmetric counter-rotating states, $\Delta T_{m}$ is not exactly zero. In our experiment, $\Delta T_{m}$ is not directly controlled. Rather, we control the temperature of the water baths pumping water through the top and bottom plates, but the plate temperatures are also affected by heat losses in the piping from the baths to the plates, and the finite conductivity of the aluminum top and bottom plates in our experiments allows them to be coupled to the temperature profile of the LSC. In practice, this coupling results in a typical change in $\Delta T_{m}$ between counter-rotating and co-rotating states of $0.002\Delta T$. That small asymmetry in driving has a large consequence, with a $|\Delta\bar{T}_{0}|/\bar{\delta}$ typically varying by 0.1 in cases where we intended to produce symmetric counter-rotating states, which can lead to significant terms in the model Eqs. 4, 5, and 6 in counter-rotating states. In particular, a general steady-state solution of Eq. 6 for both counter- and co-rotating states that includes variations in $\theta_{p,i}$ and $\bar{\delta}_{i}$ with $\Delta T_{m}$ is

[TABLE]

To test the $\theta_{p,i}$-dependent terms of Eq. 15, we plot measured values of $\Delta\bar{T}_{0}/\bar{\delta}$ as a function of measured values of $[\bar{\delta}_{2}\cos\theta_{p,2}-\bar{\delta}_{1}\cos\theta_{1}]/\bar{\delta}$ in Fig. 12. The normalization of both axes by $\bar{\delta}$ allows us to include all of our counter-rotating data at different Ra on the same scale to better interpret the magnitude of the asymmetry. Co-rotating state data with $\Delta T_{m}/\bar{\delta}<1.3$ are shown to be tightly clustered because $\delta_{2}\cos\theta_{p,2}-\delta_{1}\cos\theta_{1}=2\bar{\delta}_{co}$ consistently in co-rotating states with small $\Delta T_{m}$. For counter-rotating data, a linear function plus a constant is fit to the data in Fig. 12. The fit yields a slope of $-0.88\pm 0.05$ plus a constant offset of $0.20\pm 0.02$. The offset must be the result of some undetermined asymmetry in the nominally symmetric system. The predicted slope in Fig. 12 is $-0.93$, obtained by rearranging Eq. 15 to isolate $\Delta\bar{T}_{0}$, $\Delta T_{m}=\Delta\bar{T}_{0}/0.57$ from a fit of Fig. 11c, and the measured value $\mbox{Nu}\kappa/\kappa_{t}=0.017\pm 0.002$. The predicted slope is within a standard deviation of the fit slope of $-0.88\pm 0.05$ in Fig. 12, so this confirms the validity of Eq. 15, and thus the detailed functional form of the third term of Eq. 6. Since Eqs. 4 and 5 for $\dot{\delta}_{i,\kappa}$ and $\dot{\theta}_{i,\kappa}$ depend directly on $T_{0}$, the small uncontrolled differences in $\Delta T_{m}$ and resulting changes in $\Delta\bar{T}_{0}$ likely account for much of the scatter in plots such as Fig. 4, and corresponding asymmetries in Fig. 3.

VII.4 Trends of $\bar{\delta}_{i}$ with $\Delta T_{m}$

VII.4.1 Co-rotating states

In the co-rotating state, $\bar{\delta}$ increases linearly with $\Delta T_{m}$ in Fig. 11b. This change is predicted from the linear expansion in Eq. 12 since $\Delta\bar{T}_{0,co}$ also increases linearly with $\Delta T_{m}$. Using the fit $\partial\Delta\bar{T}_{0,co}/\partial\Delta T_{m}=(0.061\pm 0.014)$ in the co-rotating state for $|\Delta T_{m}|<0.2$ K from Fig. 11c, $\tau_{\delta}=81$ s from a cubic cell with a single LSC at the same Ra $=5.5\times 10^{8}$ Ji and Brown (2020b), and $\kappa_{t}/\kappa=3200\pm 400$ from the fit in Fig. 8 yields a prediction $|\partial\bar{\delta}_{co}/\partial\Delta T_{m}|=(0.019\pm 0.004)$, a factor of 2 smaller than the measured slope $|\partial\bar{\delta}_{co}/\partial\Delta T_{m}|=0.042\pm 0.006$ for $\Delta T_{m}<-0.2$ K in Fig. 11b. There is a reduction in the slope of both $\Delta\bar{T}_{0,co}$ and $\bar{\delta}_{co}$ with $\Delta T_{m}$ for 0.2 K $<\Delta T_{m}<1$ K. This change in slope is not predicted by the model. In the range 0.2 K $<\Delta T_{m}<1$ K, using the measured slope $\partial\Delta\bar{T}_{0,co}/\partial\Delta T_{m}=(0.019\pm 0.007)$ the predicted slope $|\partial\bar{\delta}_{co}/\partial\Delta T_{m}|=0.006\pm 0.001$ is 5 times smaller than the measured slope $0.031\pm 0.002$.

VII.4.2 Asymmetry $\delta_{2}-\delta_{1}$ in Counter-rotating states

When a small $\Delta T_{m}$ produces counter-offset states, there is a shift in $\bar{\delta}$ with increasing $\Delta T_{m}$ in opposite directions in the two cells. This shift $\Delta\bar{\delta}\equiv(\bar{\delta}_{2}-\bar{\delta}_{1})/2$ is shown in Fig. 13 for the counter-offset data from Fig. 11. Such a difference is predicted from Eq. 4. $\Delta\bar{\delta}$ is obtained by inserting Eq. 4 with a linear expansion in $\Delta\theta_{p}$ around a diagonal (e.g. $\theta_{p}=\pi/4$ rad) into a linear expansion around the stable fixed point of $\delta=\delta_{0}$ in Eq. 2, resulting in

[TABLE]

Using $\tau_{\delta}=81$ s from measurements in a single cubic cell Ji and Brown (2020b), the fits of $\Delta\bar{T}_{0}=0.57\Delta T_{m}$ and $\Delta\theta_{p}=(3.7\pm 0.4)\Delta T_{m}$ rad/K from Fig. 11, and measurements $\kappa_{t}/\kappa=3200\pm 400$ and $\bar{\delta}=0.12$ K, yields the prediction $\Delta\bar{\delta}=(-0.039\pm 0.006)\Delta T_{m}$. The data in Fig. 13 is fit by a linear function plus a constant, yielding $\Delta\bar{\delta}=(-0.046\pm 0.006)\Delta T_{m}+0.005$ K. The offset of 0.005 K is indicative of an asymmetry of the setup, comparable to the uncertainty on thermistor measurements. The fit slope is consistent with the predicted value and within 20%, which confirms the validity of the detailed form of the first and second terms of Eq. 4.

VII.4.3 Second order expansion of $\bar{\delta}(\Delta T_{m})$

Figure 11b showed an increase in $\bar{\delta}$ with $|\Delta T_{m}|$ in counter-offset states. Since it increases with both positive and negative $\Delta T_{m}$, we fit a quadratic function, obtaining a curvature $\partial^{2}\bar{\delta}/\partial\Delta T_{m}^{2}=0.26\pm 0.11$/K. This quadratic trend can be predicted from Eq. 4 using the prediction and observation from Fig. 11 that both $\Delta\theta_{p}$ and $\Delta\bar{T}_{0}$ are linear in $\Delta T_{m}$, so that $\cos\theta_{p,i}=\sqrt{2}(1\pm\Delta\theta_{p})$. Averaging over cells 1 and 2, the first-order expansions in $\Delta\bar{T}_{0}$ or $\theta_{p}$ of the first two terms of Eq. 4 have opposite signs for the two different cells at a diagonal (e.g. $\theta_{i}=\pi/4$ rad), and cancel out in the average over the two cells (these are the terms that were responsible for $\Delta\bar{\delta}$). The expansion of the first and third terms of Eq. 4 up to second-order in $\Delta\theta_{p}$ and $\Delta\bar{T}_{0}$ yields $\overline{\dot{\delta}}_{counter}=(A\kappa_{t}/\pi L^{2})(2\Delta\theta_{p}\Delta\bar{T}_{0}-\pi\delta_{0}\Delta\theta_{p}^{2}/8)$. Plugging this into a first-order expansion of Eq. 2, using the fit $\Delta\bar{T}_{0}=(0.57\pm 0.05)\Delta T_{m}$ from Fig. 11c, the fit $\Delta\theta_{p}=(3.7\pm 0.4)\Delta T_{m}$ rad/K from Fig. 11a, and using the measured value of $\bar{\delta}$ at $\Delta T_{m}=0$ to approximate $\delta_{0}=0.12$ K yields the prediction

[TABLE]

Using $\tau_{\delta}=81$ s Ji and Brown (2020b) and $\kappa_{t}/\kappa=3200\pm 400$, the predicted curvature in $\Delta T_{m}^{2}$ is $0.54\pm 0.18$, consistent with the fit curvature $0.26\pm 0.11$/K. This confirms that Eq. 4 is consistent with measurements even to 2nd order.

The self-consistency of all of the model predictions of the dependencies of $\Delta\theta_{p}$, $\Delta\bar{T}_{0}$, $\bar{\delta}$, and $\Delta\bar{\delta}$ on $\Delta T_{m}$ was confirmed within a factor of 3 in Sec. VII.1, VII.2, VII.3, and VII.4. This confirms the validity of the dependence of Eqs. 4, 5, and 6 on $\Delta T_{m}$.

VII.5 Does the model quantitatively predict the conditions for stable co-rotating states?

In this section, we check whether the model can quantitatively predict the conditions in which co-rotating states are stable. To be stable at the alignment $\theta_{p}=0$ or $\pi$ rad, there has to be a forcing to overcome the geometric forcing $-\nabla V_{g}$, which has unstable fixed points at $\theta_{p}=0$ and $\pi$ rad Ji and Brown (2020b). A first-order expansion in $\theta_{i}-\theta_{p}$ at fixed $\delta_{i}$ of the predicted unstable geometric forcing is $-\nabla V_{g}\approx\pi\omega_{r}^{2}\bar{\delta}_{co}^{2}(\theta_{1}-\theta_{p})/10\delta_{0}^{2}$ Ji and Brown (2020b). This expansion for the geometric potential was untested due to limited data near the unstable fixed points of the potential $V_{g}$ in previous experiments Ji and Brown (2020b). We compare this predicted unstable forcing to the stabilizing forcing from the neighboring roll interaction, using a linear expansion in $\theta_{i}-\theta_{p}$ of Eq. 5 for a fixed $\delta_{i}$ (this corresponds to the second term of Eq. 8). We calculate the magnitude of the ratio $S$ of stabilizing forcing for co-rotating states to the destabilizing forcing from the expansion of the geometric forcing$-\nabla V_{g}$ to obtain the forcing ratio

[TABLE]

To predict values of $S$ for different datasets, we use $\kappa_{t}/\kappa=0.8Ra^{0.41}$ from the fit in Fig. 8, $\tau_{\dot{\theta}}=2.7\times 10^{4}Ra^{-0.35}$ from a fit of data with a single LSC in a cubic cell data Ji and Brown (2020b), $\omega_{r}^{2}=6.3\times 10^{-13}Ra^{0.98}$ from linear regression of data with a single LSC in a cubic cell Ji and Brown (2020b), $\bar{\delta}/\delta_{0}=0.57$ for co-rotating states from Fig. 9, and the measured $\Delta\bar{T}_{0,co}/\bar{\delta}_{co}$ for each dataset. The ratio of stabilizing to destabilizing forces is plotted Fig. 14a for each co-rotating state, including data at all measured Ra and $\Delta T_{m}$ for completeness. Since counter-diagonal states are predicted to be stable to both forcings, the stabilizing ratio is not relevant for them. Co-rotating states are found to be concentrated in Fig. 14 with $S>1.13$. This concentration of many points just above this measured minimum of $1.13$ suggests a critical value slightly below this is required for co-rotating states. This threshold is near the predicted threshold of 1, well within typical uncertainties of this model. This confirms that the model successfully predicts the conditions for co-rotating states to be stable due to turbulent thermal diffusion, when it overcomes the geometric forcing that destabilizes co-rotating state orientations. In co-rotating states, the stabilizing factor $b\sin(\theta_{2}-\theta_{1})$ is typically small (on average, $b=0.36a$ in co-rotating states), thus while an understanding of the source of the stabilizing $b$ in co-rotating states might adjust the threshold value of $S$, the stability of co-rotating states can be explained independent of the unexpectedly stabilizing values of $b$.

VII.6 Unpredicted stable state at large $\Delta T_{m}$

In Fig. 14, the ratio of stabilizing to destabilizing forcings for co-rotating states $S$ is seen to drop towards the threshold of stability at $S=1$ for large $\Delta T_{m}$. This linearly decreasing $S$ and loss of stability is predicted because the stabilizing forcing (which scales as $\Delta\bar{T}_{0}/\delta_{0}$) is growing linearly with $\Delta T_{m}$, but not as fast as the linearly increasing destabilizing term (which scales as $\bar{\delta}_{co}^{2}/\delta_{0}^{2}$) (see trends in Fig. 11). A linear fit of $S$ for a series of points at $\Delta T=3.8$ K for $\Delta T_{m}>0.14$ K is shown to intercept the predicted threshold of stability at $\Delta T_{m}/\bar{\delta}_{co}=10$ in Fig. 14. At this point, both the co-rotating state and counter-rotating state are predicted to be unstable, so some new state should appear.

To see if the co-rotating state destabilizes at larger $\Delta T_{m}$, we plot $\bar{\delta}$ and $\Delta\bar{T}_{0}/\bar{\delta}$ as a function of $\Delta T_{m}$ for fixed $\Delta T=3.8$ K (Ra $=5.5\times 10^{8})$ in Fig. 15, extending the range of experiments from Fig. 11. For $\Delta T_{m}<0.9$ K, the data is the same as Fig. 11, with the predicted instability threshold of the co-rotating state at $\Delta T_{m}=0.9$ K. When $\Delta T_{m}>0.9$ K, new scaling behaviors appear where $\bar{\delta}$ starts to decrease with $\Delta T_{m}$, and $\Delta\bar{T}_{0}/\bar{\delta}$ starts to increase significantly above its typical value of $2.22\pm 0.13$ for co-rotating states. Values of $\theta_{p,i}$ are not shown because they remain constant as $\Delta T_{m}$ increases. While the values of $\theta_{p,i}$ appear as if this is still a co-rotating state, the decrease in $\bar{\delta}$ and large increase in $\Delta\bar{T}_{0}/\bar{\delta}$ are inconsistent with Eqs. 4 and 6. In particular, values of $\Delta\bar{T}_{0}/\bar{\delta}\gg 2$ mean the temperature on the two sides of the interface between neighboring rolls is very different, and no longer dominated by a balance of turbulent thermal diffusion across the interface. A likely candidate for this behavior is more coherent advection of heat between the cells, which can transport heat more coherently between the cells to balance the larger turbulent thermal diffusion for $\Delta\bar{T}_{0}/\bar{\delta}\gg 2$. However, we cannot make specific predictions without more detailed knowledge of such advective flow fields. The rapid increase of $\Delta\bar{T}_{0}/\bar{\delta}$ means $S$ increases again for $\Delta T_{m}>0.9$ K and does not drop below the threshold of instability for co-rotating orientations, so the stability of this high-$\Delta T_{m}$ state is still consistent with Eq. 5 for $\dot{\theta}_{\kappa,i}$, and perhaps only Eq. 6 for $\dot{T}_{\kappa,i}$ requires modification.

VII.7 Stability of counter-axis states

Counter-axis states – in which the preferred orientation is somewhere between a corner and axis – were found for $\Delta T<2$ K instead of counter-diagonal states. However, counter-axis states are not predicted to be stable from Eq. 5.

What might cause counter-axis states? Equations of motion in $\ddot{\theta}_{i}$ could in principle produce stable counter-axis states if there is a force pushing both cells to be stable in $\pm\sin\theta_{i}$. For example, if there was an asymmetry in the temperature profile of the plates such that the both plates are hotter near the interface. It is notable that our counter axis states all have an unexpectedly large $\Delta\bar{T}_{0}>\Delta T_{m}$, which is inconsistent with Eq. 15 and thus Eq. 6 for turbulent thermal diffusion. Thus, the large $\Delta\bar{T}_{0}$ must come from some other mechanism for these counter-axis states. This unexpectedly large $\Delta\bar{T}_{0}$ might be an indication of asymmetric heating within the plates, although the forcing that would push both cells toward a counter-axis state would not be best represented by the parameter $\Delta\bar{T}_{0}$. We found counter-axis states to be more likely when the flow direction in some of the cooling baths controlling the plate temperatures was switched (individual plate thermistor mean temperatures typically changed by $0.001\Delta T$), suggesting the plate temperature profile can produce such an asymmetry. We also found counter-axis states to be more prominent in an early version of the apparatus with a different middle wall covering only the top half of the interface between the two rolls. This suggests that there was some stabilizing force in $\pm\sin\theta_{i}$ with the same sign in both cells, presumably due to an interaction of the LSCs with the middle wall Ji (2019). This asymmetry could have resulted from an interaction between a corner-roll and the middle wall, as corner-rolls are more prominent at the top only on the side of the up-flow (not down-flow), so the interaction could be different for a top-half middle wall depending on whether the flow in the middle is upward or downward. Whatever the source, it seems likely that the observed counter-axis states are due to some unpredicted asymmetry of our setup, and the forcing can be surprisingly significant in systems with mild-seeming asymmetries.

VII.8 Tilt

Tilting the cell relative to gravity by an angle $\beta$ in the vertical plane going between the center of both cells results in a component of buoyant forcing aligned with $\theta=0$ in one cell and $\theta=\pi$ rad in the other cell Chillà et al. (2004); Brown and Ahlers (2008b); Ji et al. (2020). Since this is similar to the effect of a temperature difference $\Delta T_{m}$, it should be no surprise that a tilt of about $\beta=2$ degrees generally results in co-rotating states, with a mix of co-rotating and counter-offset states at smaller $|\beta|$ (data is included in Fig. 4), analogous to Fig. 11a. A detailed analysis of tilting a single cubic cell was presented previously Ji et al. (2020), and the same forcing terms are expected to apply here, in addition to the forcing in Eqs. 4 and 5. Some quantitative data on the tilt-dependence is presented for the earlier version of our experiment with the middle wall only in the top half of the cell, which also included an extra forcing term that is responsible for the counter-axis states Ji (2019). These observations serve as a confirmation that forcing terms from different physical mechanisms can be added in the low-dimensional model.

VIII Dynamics of switching between co-rotating and counter-rotating

We observed numerous stochastic switching events between co-rotating and counter-rotating states such as in Fig. 3. To obtain statistics to characterize these events, we searched through experiments with a total run time of 148 days with different Ra, $\Delta T_{m}$, and tilt angles $\beta$. To intentionally produce more switching events, we also occasionally enforced a co-rotating state in one direction or the other by biasing $|\Delta T_{m}|/\delta_{0}\stackrel{{\scriptstyle>}}{{{}_{\sim}}}1.7$ using Fig. 11 as a guide.

To systematically identify switching events between counter-rotating and co-rotating states, we define a transition based on the time-dependent difference in LSC orientations $|\theta_{1}-\theta_{2}|$. We first smoothed values of $\theta_{i}$ over a duration of $1.5\tau_{\delta}$. We define a transition from counter-rotating to co-rotating as starting when $|\theta_{1}-\theta_{2}|$ last exceeds $\pi/2$ rad, and ending when $|\theta_{1}-\theta_{2}|$ first exceeds $3\pi/4$ rad, without returning below $\pi/4$ rad in between. These first and last crossing times at different thresholds are used to avoid counting jitter around these threshold values. We set the threshold $\pi/2$ rad further from the ideal counter-diagonal state value to include counter-axis states where $|\theta_{1}-\theta_{2}|$ is larger. Likewise, a change from co-rotating to counter-rotating is defined to start when $|\theta_{1}-\theta_{2}|$ last drops below $3\pi/4$ rad, and end when $|\theta_{1}-\theta_{2}|$ first drops below $\pi/2$ rad, without returning above $3\pi/4$ rad in between. To make sure the change in $|\theta_{1}-\theta_{2}|$ is not a false event due to a temporary fluctuation in one variable, we also calculated short-time averages of parameters $|\theta_{1}-\theta_{2}|$, $(\delta_{1}+\delta_{2})/2$, and $|T_{0,2}-T_{0,1}|$ in the steady states before and after the event – in the final analysis we used the range of $14\tau_{\delta}$ to $56\tau_{\delta}$ before and after the event, excluding data where $\cos(\theta_{2}-\theta_{1})$ has the wrong sign for the expected state, and excluding data with either $\delta_{i}<0.3\delta_{0}$. We only considered a switching event to have occurred if the parameters $|\theta_{1}-\theta_{2}|$, $(\delta_{1}+\delta_{2})/2$, and $|T_{0,2}-T_{0,1}|$ averaged over the specified time range changed from before to after the event by more than 60% toward the expected counter- or co-rotating state: the expected differences going from counter- to co-rotating states are $|\theta_{1}-\theta_{2}|$ increased by $\pi$ rad (when values of $\theta_{1}-\theta_{2}$ are reduced to the range $-\pi$ rad to $\pi$ rad), $|T_{0,2}-T_{0,1}|$ increased by $2.22\bar{\delta}_{co}$ based on Fig. 10, and $(\delta_{1}+\delta_{2})/2$ decreased by $0.37\delta_{0}$ based on Fig. 9.

We found 89 switching events between counter- and co-rotating states, corresponding to an average frequency of 0.60 per day, and only 0.48 per day if we do not count events that occurred shortly after we forced a switch to a co-rotating state by applying a large $|\Delta T_{m}|$. We only observed a few direct switching events from one counter-diagonal state to another corner. We never observed a direct switch between two co-rotating orientations. Even in cases where a co-rotating state was stable, and we changed $\Delta T_{m}$ to drive the co-rotating state in the opposite direction, the system always switched to a counter rotating state and resided there for at least 1000s before switching to the new preferred co-rotating state. Switching between two orientations of counter-axis states is relatively frequent, and statistics of switching between counter-axis orientations are reported elsewhere Ji (2019). Mechanisms of switching between counter-axis states will not be addressed here because they are driven by some uncharacterized asymmetry.

VIII.1 Different driving mechanisms for switching events

In previous work, changes in LSC orientation between stable fixed point states were often found to be due to dramatic reorientations in $\theta_{0}$, driven by extreme fluctuations in $\theta_{0}$ coming from the stochastic term $f_{\dot{\theta}}(t)$ Brown and Ahlers (2008a); Xi and Xia (2008); Bai et al. (2016); Ji (2019). Other events are driven by cessations – a reduction in $\delta$ to near zero driven by stochastic fluctuations to overcome a potential barrier due to Eq. 2. Cessations reduce the damping in Eq. 3 so that the stochastic term $f_{\dot{\theta}}(t)$ allows $\theta_{0}$ to change more rapidly Brown and Ahlers (2007a, 2008a); Xi and Xia (2007); Xie and Xia (2013). For either mechanism the behavior could be quantified statistically as a rare event due to diffusive fluctuations crossing a large potential barrier in Eq. 2 or 3, respectively Brown and Ahlers (2008a); Bai et al. (2016); Ji and Brown (2020b). However, with only 89 switching events and a large number of variables, we do not have enough statistics to meaningfully test quantitative predictions of barrier crossing rates as a function of any parameter.

Instead, to identify the relevant mechanisms for switching between counter- and co-rotating states, we calculate correlations between changes in LSC parameters before an event, with the expected change towards the new stable value after switching. These fractional changes $d_{\theta}$, $d_{\delta}$, and $d_{T}$ are calculated as the average fractional change of each parameter ($|\theta_{1}-\theta_{2}|$, $(\delta_{1}+\delta_{2})/2$, and $|T_{0,2}-T_{0,1}|$, respectively) in the time period of $2\tau_{\delta}$ before the start of a switching event, relative to the expected change towards the stable state value after switching. These fractions are averaged separately for events switching to either a counter- or co-rotating state, and summarized in Table 2. We also include a row for a fractional change toward cessation driving $d_{c}$, calculated as the fractional change in $\delta$ toward zero ($\bar{\delta}_{i}-\delta_{i}$) relative to the stable value of $\bar{\delta}$ before switching . We only count the maximum $d_{c}$ among the two cells, since a cessation of one of them is enough to cause a change in state.

Table 2 suggests that switching from the co-rotating to counter-rotating state is most strongly correlated with cessation, with an average $d_{c}=0.36$, indicating a significant drop in $\delta$ from its previous stable fixed point in the co-rotating state, which weakens the damping of motion in $\theta_{i}$ (Eq. 3) Brown and Ahlers (2008a). This strong correlation is especially notable because the new stable fixed point in the counter-rotating state has a larger $\delta_{i}$, meaning the cessation initially moves $\delta_{i}$ further from its new stable fixed point value at the start of the transition between states. All but one of the switching events to counter-rotating states has a reduction in one of the $\langle\delta_{i}\rangle$ either within the $2\tau_{\delta}$ before the event, or during the event, The single event with negative $d_{c}$ is telling, because it occurred in the first 4000 s of one of the counter-offset states in Fig. 11a, which we did not find to have a longer stable lifetime in a counter-rotating state, suggesting the potential barrier that needed to be crossed was relatively low in that case.

An example of a switching to a counter-rotating state Ra $=2.7\times 10^{9}$ is shown in Fig. 16. To represent a typical state, we chose the event with the median value of $d_{c}$ among switching events to counter-rotating states. In this example, there is a large drop in $\delta_{1}$ in cell 1 that allows a rapid change of $\theta_{1}$ due to smaller damping on $\theta_{1}$ from Eq. 3 Brown and Ahlers (2008a). Once $\theta_{1}$ reaches near a new stable fixed point orientation for a counter-rotating state, the other parameters $\theta_{2}$, $T_{0,1}$, $T_{0,2}$, $\delta_{1}$ and $\delta_{2}$ follow to their new stable fixed point values within a few turnover times, or a few $\tau_{\delta}$.

Generally, switching from co-rotating to counter-rotating states due to cessation (reduced $\delta_{i}$) can be understood to be due to reduced damping in Eq. 3, which allows $f_{\dot{\theta}}$ to drive large changes in $\theta_{i}$ in the same cell to cross the potential barrier that stabilizes co-rotating states in Eq. 5. We generally find during switches to counter-rotating states that the drop in $\delta_{i}$ causes $\theta_{i}$ in the same cell to change further than the neighboring cell, consistent with this explanation. Note that the reduction in $\delta_{i}$ only allows the motion of $\theta_{i}$ to be dominated by the stochastic term if the damping factor in Eq. 25 includes the factor of $\delta_{i}/\delta_{0}$ to cancel out the $1/\delta_{i}$ in Eq. 23, so this observation is our best justification for using that factor $\delta_{i}/\delta_{0}$ in Eq. 25. In principle, Eq. 5 also allows for the co-rotating state to become unstable due to large fluctuations that reduce $|T_{0,2}-T_{0,1}|$. We see no evidence that this is a major driving factor in switching from co-rotating to counter-rotating states ($d_{T}=0.01$) without the aid of a cessation, perhaps because fluctuations in $T_{0}$ tend to be relatively small (see Fig. 3).

In the other case of switching from a counter-rotating state to a co-rotating state, Table 2 shows that cessation ($d_{c}$) has the weakest correlations with these switching events, and is difficult to distinguish from a driving towards $\delta_{co}$, which also involves a reduction in $\delta$ in this case. The correlations with an increase in $|T_{0,2}-T_{0,1}|$, decrease in $(\delta_{1}+\delta_{2})/2$ and increase in $|\theta_{1}-\theta_{2}|$ are all stronger than the driving due to cessation, and much stronger than when switching to a counter-rotating state.

An example of a switching event from a counter-rotating state to a co-rotating state led by random diffusive motion in $T_{0,2}-T_{0,1}$ at Ra $=2.7\times 10^{9}$ is shown in Fig. 17. To represent a typical switching to a co-rotating state, we chose the event with nearly median values in each of $d_{\theta}$, $d_{T}$, and $d_{\delta}$. In this example, $T_{0,1}$ and $T_{0,2}$ start slowly drifting apart about 66,200 s, well before the change in $\theta_{i}$ and $\delta_{i}$ to their new stable values in a co-rotating state. Generally, as $|T_{0,2}-T_{0,1}|$ increases, it produces an increasing forcing on $\theta_{i}$ in the third term of Eq. 5 to push $\theta_{1}$ and $\theta_{2}$ apart. Additionally, fluctuations in both $\theta_{i}$ and $\delta_{i}$ away from the counter-rotating state fixed point drive an increase in $|T_{0,2}-T_{0,1}|$ in Eq. 6, providing a positive feedback that destabilizes a counter-rotating state to drive a transition to a co-rotating state. This positive feedback can be seen in Fig. 17 as the changes in $\theta_{i}$, $\delta_{i}$, and $T_{0,i}$ all start to increase more rapidly at about the same times. In Eq. 5, when $\sin\theta\approx\sqrt{2}/2$ near the stable counter-rotating state, the third term that destabilizes the counter-rotating state becomes dominant over the first term that stabilizes the counter-rotating state when $T_{0,2}-T_{0,1}>\pi\delta/4$ (alternatively, the threshold of instability in Eq. 5 can also be identified as $|\theta_{1}-\theta_{2}|>\pi/2$ as seen in Fig. 11). In the example in Fig. 17, both of these thresholds are crossed roughly in between the first 2 timing lines, which is about when the rate of change in $\theta_{i}$ starts to increase. Once the threshold of instability is crossed, the values of $\theta_{i}$, $\delta_{i}$, and $T_{0,2}-T_{0,1}$ quickly move towards their stable fixed point values for co-rotating states. Since the positive feedback between increasing $|T_{0,2}-T_{0,1}|$, decreasing $(\delta_{1}+\delta_{2})/2$, and increasing $|\theta_{1}-\theta_{2}|$ is what destabilizes counter-rotating states, this explains why $d_{T}$, $d_{\delta}$, and $d_{\theta}$ are all correlated with switching to the co-rotating state in Table 2.

The behavior described above for switching from counter-rotating to co-rotating states is qualitatively similar in most events. Notably, during switching from counter- to co-rotating states, we generally observe that $T_{0,1}$ and $T_{0,2}$ are driven in opposite directions at nearly equal rates, as seen for example in Fig. 17, and as expected from Eq. 6. The cell that starts slightly hotter also generally continues to get hotter, as seen for example in Fig. 17c. This determines which orientation is easier to switch to, as smaller fluctuations in $\theta_{i}$ and $T_{0,2}-T_{0,1}$ in this direction are required to destabilize the counter-offset state in Eq. 6.

If switching to counter-rotating states is due to cessation and other large drops in $\delta_{i}$, we can predict that switching to counter-rotating would occur in some fraction of cessations, which depends mainly on the mean value of $\delta_{i}$. Any forcing that leads to a larger $\bar{\delta}_{i}$ is expected to decrease the likelihood of events Brown and Ahlers (2008a); Bai et al. (2016), such as increasing $\Delta T_{m}$ or tilting the cell. On the other hand, a prediction of switching to co-rotating states would instead depend on crossing barriers in a three-dimensional space of $|T_{0,2}-T_{0,1}|$, $(\delta_{1}+\delta_{2})/2$, and $|\theta_{1}-\theta_{2}|$. For either type of switching, due to the rarity of events, we do not have enough statistics to test such predictions.

IX Discussion

IX.1 Rejected model: modified drag at the interface between two rolls

We originally considered a modification of viscous drag at the interface to model the interactions between neighboring rolls. Since the effects of viscous drag at the interface between neighboring rolls could be modeled as a turbulent viscous diffusion of a vertical velocity $U\cos(\theta_{i}-\theta_{0})$ at the interface, some of its effects could be in part mathematically analogous to turbulent thermal diffusion of temperature $T_{i}=T_{0,i}+\delta\cos(\theta_{i}-\theta_{0})$.

One important difference between the velocity and temperature is the offset temperature $T_{0,2}-T_{0,1}$ between the neighboring cells, as there is no analogous offset for velocity. The difference between the mean temperature of neighboring rolls $T_{0,2}-T_{0,1}$ is important because it is the term that stabilizes the observed co-rotating state in Eq. 5. Without a mathematically analogous offset, viscous drag would only tend to align neighboring rolls toward a counter-rotating state. $T_{0,2}-T_{0,1}$ is also the intermediary parameter for effects of the difference in mean plate temperature of the two cells $\Delta T_{m}$. Thus, we rule out the hypothesis that the drag force is the most important interaction between neighboring convection rolls in our experiments, at least when it comes to determining the temperature profile and stable states of the flow structures.

On the other hand, viscous drag may play a smaller role in modifying parameter values. Since the Prandtl number of the water is $Pr=\nu/\kappa=6.41$, viscous diffusion should in principle be stronger than thermal diffusion by that factor. Since we could represent both the temperature difference and velocity profiles as sine functions in the coordinate $\theta$ Brown and Ahlers (2008a), the diffusion equations and effective thermal diffusivity $\kappa_{t}$ measured may actually contain contributions from viscous diffusion. Larger coefficients depending on $Pr$ might be expected for terms without $T_{0,2}-T_{0,1}$. However, since the threshold for stability of co-rotating states is just barely met by the third term of Eq. 5 overcoming the first and second terms, a significant destabilizing drag added to other terms would overwhelm the stabilizing term which depends on $T_{0,2}-T_{0,1}$, resulting in co-rotating states being unstable, in disagreement with observations. Thus, any contribution of drag is probably small compared to the stabilizing term in our experiments. It remains to be seen if systems with larger $Pr$, or those driven by forced convection rather than buoyancy, are dominated by turbulent viscous diffusion rather than turbulent thermal diffusion.

In Eq. 4 for $\dot{\delta}_{\kappa,i}$, a modified drag term would be smaller in a counter-rotating state relative to a single LSC due to the flow in the same direction at the interface. On the other hand, this modified drag may also change the time a fluid parcel spends near the top and bottom plates to lose or gain heat, respectively, and change the rate of turbulent mixing. These effects may cancel each other out in the temperature profile, since there is little change in $\bar{\delta}$ [we observed a 6% reduction of $\bar{\delta}$ in the counter-rotating state compared to a single roll (Fig. 9)]. The fact that there is not a larger change in $\bar{\delta}$ in the counter-rotating state suggests that any change in drag due to neighboring rolls has a small effect on $\bar{\delta}$.

IX.2 Discussion of previous work on co-rotating states

Previous observations of co-rotating rolls in systems with horizontal forcing due to lateral heating Pérez-García et al. (2004) or tilt Mercader et al. (2019) are analogous to a one-dimensional version of our model with forcing due to $\Delta T_{m}$. We found that motion in $\theta_{i}$ is important for determining the stability of states in three dimensions. In particular, the stability of counter- and co-rotating states is determined by Eq. 5, which has terms that both stabilize and destabilize flow in terms of $\theta_{i}$, and stabilize in $\theta_{1}-\theta_{2}$ (Sec. VI.1). In addition, in any non-circular cross-section, the forcing $-\nabla V_{g}(\theta_{i})$ due to the geometry of the cell can potentially add to the stabilizing or destabilizing terms (Sec. VII.5).

The behaviors of stacked co-rotating rolls Xie and Xia (2013) are not explained by our model, as the forcing terms in the vertical direction would be different, and vertically stacked co-rotating states likely have more coherent advection between the two cells.

IX.3 Preferred states in other systems

Our observations of stable co-rotating states leads to the question of why are co-rotating states not reported more often in other convective flows? Equation 5 predicts stable co-rotating states when the stabilizing forcing in $\theta_{1}$ from the 3rd terms of Eq. 5 are larger than the 1st and 2nd terms that stabilize counter-rotating states and any destabilizing forcing $-\nabla V_{g}(\theta_{i})$ due to the geometry of the cell (Sec. VII.5). For co-rotating states to be stable with $\theta_{2}=\theta_{1}+\pi$ rad and assuming $\bar{\delta}_{2}=\bar{\delta}_{1}$, this requires $\Delta\bar{T}_{0,co}/\bar{\delta}_{co}>1.62$, not including contributions from $V_{g}$ which can increase this threshold by an amount depending on the geometry. Equation 6 gave a prediction of $\Delta\bar{T}_{0,co}/\bar{\delta}_{co}=1.76$, which is only 9% greater than the threshold of stability for co-rotating states, and we measured $T_{0,co}/\bar{\delta}_{co}=2.22$, which is 40% larger than this threshold. We do not yet know how $\Delta\bar{T}_{0,co}/\bar{\delta}_{co}$ varies in other systems. Since the numerical coefficients in Eq. 5 are expected to change with the shape of the interface and the shape of the convection patterns, it seems plausible that a version of Eqs. 5 and 6 for another system might not produce a large enough $\Delta\bar{T}_{0,co}/\bar{\delta}_{co}$ to be above the threshold required for stability of co-rotating states. These variables, might be part of the reason that co-rotating states are not often reported. In other words, while this idealized laboratory system can reproduce many qualitative and quantitative features of natural convection rolls, the large number of model terms and the sensitivity of behavior to parameter values means that behavior may not be representative of all natural convection systems.

Equation 5 by itself does not predict any preferred orientation for counter-rotating states. Thus, counter-rotating rolls can in principle be aligned at any orientation. Specific systems likely have a preference for certain orientations based on their geometric potential $V_{g}$ (i.e. the longest diagonals of a closed cell), or other forces. The counter-rotating predicted by Eq. 5 corresponds to $\theta_{1}=\theta_{2}$ such that the LSC orientation vectors line up head-to-head, and do not in general have to be parallel. For example, a hexagonal convection pattern with 3 counter-rotating rolls at 120 degree angles relative to each other meeting at a vertex with adjacent upward-flowing regions meet this criteria for counter-rotating states – in this case it is the tessellation of a hexagonal pattern that determines the preferred orientations of counter-rotating states Bailon-Cuba et al. (2010).

X Conclusions

We observed both counter-rotating and co-rotating states to be stable for adjacent convection rolls in cubic Raleigh Bénard convection cells with $7.6\times 10^{7}<$ Ra $<3.5\times 10^{9}$ and Prandtl number 6.41 (Figs. 3, 4). Notably co-rotating states were observed in nominally symmetric systems which are usually argued to support only counter-rotating states. The stability of the observed counter-rotating states and co-rotating states can be explained by a low-dimensional model consisting of three stochastic ordinary differential equations for each roll. The interactions between neighboring rolls are predicted from the scaling of heat transport by turbulent thermal diffusion, and added to existing model equations Brown and Ahlers (2008a, b); Ji and Brown (2020b) for the LSC amplitude $\delta_{i}$ (Eqs. 2, 4), orientation $\theta_{i}$ (Eqs. 3, 5), and a new equation for mean temperature $T_{0,i}$ (Eq. 6). Despite there only being one new mechanism in the model, projecting turbulent thermal diffusion for two neighboring rolls with just simple sinusoidal temperature profiles into the parameters $\theta_{i}$, $\delta_{i}$ and $T_{0,i}$ resulted in eleven new model terms.

To make quantitative predictions, the model required a parameter as experimental input: a turbulent thermal diffusivity $\kappa_{t}$ in the horizontal direction across the interface between neighboring rolls. Once $\kappa_{t}$ was measured based on a fit of one model term, it could be used to predict the magnitudes of all eight model terms in which $\kappa_{t}$ affects steady-state solutions for $5.5\times 10^{8}<$ Ra $<3.5\times 10^{9}$. Predictions were not tested in the lower measured range of Ra due to an unknown additional forcing from some asymmetry that was unquantified. One term of the thermal diffusivity model depends on the Nusselt number for vertical heat transport, and the remaining two predicted terms of the model affect transients only, which were not tested. Specifically, the turbulent thermal diffusion model predicts the stabilizing forcing in $\sin(\theta_{1}-\theta_{2})$ for counter-rotating states and the stabilizing forcing in $\theta_{i}$ for parallel co-rotating states (Figs. 5, 6), given the value of $\kappa_{t}$ obtained from an independent measurement of correlation time of $T_{0,i}$ (Figs. 7, 8). The co-rotating states are stable if the difference $\Delta\bar{T}_{0,co}$ in the mean temperature of the two rolls is large enough that the heat flux between the rolls stabilizes the temperature profile of aligned co-rotating states. It was confirmed that the model correctly predicts the value $\Delta\bar{T}_{0,co}$ is large enough to produce a stabilizing forcing (Figs. 9, 10), and that we find stable co-rotating states in our cubic cell if and only if this forcing is large enough to overcome the forcing due to the cell geometry (Fig. 14). This large $\Delta\bar{T}_{0,co}$ could be driven by a difference $\Delta T_{m}$ between the mean temperatures of the plates of the two cells (Fig. 11), tilt of the cell (Fig. 4), or even by stochastic fluctuations from a stable counter-rotating state in a nominally symmetric setup where both states are stable (Fig. 3). These predictions are all consistent within a factor of 3 of measurements when using the same value of $\kappa_{t}$ for each term, a typical uncertainty of this modeling approach Brown and Ahlers (2008a, b); Ji and Brown (2020b). The model can predict the decrease in the LSC amplitude $\bar{\delta}_{co}$ in co-rotating states if the measured $\Delta\bar{T}_{0,co}$ is used as an input (Fig. 9), however, the model is not accurate enough to predict the sign or lower bound of this change without experimental input due to the predictions depending on differences of terms with comparable magnitude and having large uncertainties.

When an asymmetry between neighboring cells is introduced through an imposed temperature difference $\Delta T_{m}$ between the mean plate temperatures of the two cells, the model for effective turbulent thermal diffusion explains an observed increase in $\Delta\bar{T}_{0}$ with increasing $\Delta T_{m}$, which shifts the preferred orientations of counter-rotating rolls in opposite directions by $\Delta\theta_{p}$, until the shift is large enough that the state becomes unstable and switches to a co-rotating state (Fig. 11). The model self-consistently predicted the trends $\Delta\theta_{p}$, $\bar{\delta}$, and $\Delta\bar{T}_{0}$ in $\Delta T_{m}$ for both counter- and co-rotating states, and a difference $\Delta\bar{\delta}$ between the two cells in counter-rotating states, when using the other trends as input (Figs. 11, 12, 13) all within a factor of 3. However, the sensitivity of the model to differences in estimated terms makes it unable to predict these trends without any experimental input.

We also observed spontaneous switching events of the LSC between counter-rotating and co-rotating states. Events switching from co-rotating to counter-rotating states were almost all driven by cessations of the LSC – or at least a significant reduction in $\delta$ (Fig. 16) – which reduces the damping term in Eq. 3, so that the orientation of the LSC can be driven easily by random turbulent fluctuations Brown and Ahlers (2008a); Bai et al. (2016). We found that switching from counter-rotating states to co-rotating states is generally driven by a combination of fluctuations that cause increases in $|T_{0,2}-T_{0,1}|$, decreases in $\delta_{i}$, and increases in $|\theta_{1}-\theta_{2}|$, which result in a positive feedback in Eqns. 5 and 6 that destabilizes the counter-rotating state (Fig. 17).

The observations above tested each of the model terms that affect stable fixed point values (including the first 3 terms of Eq. 5, the first 3 terms of Eq. 4, and all 3 deterministic terms of Eq. 6), confirming the validity of the turbulent thermal diffusion model for the interaction between neighboring rolls. This level of detail suggests that the hypothesized mechanism of turbulent thermal diffusion is dominant, and the $\kappa_{t}$ that controls the interactions between neighboring convection rolls is the same turbulent thermal diffusivity that describes macroscopically averaged heat transport. It remains to be confirmed whether the $\kappa_{t}$ measured here based on the interaction between convection rolls has the same magnitude as an effective thermal diffusivity for macroscopically averaged heat transport, or if the value of $\kappa_{t}$ can be predicted from a general scaling model. It also remains to be confirmed if the prediction for the term of Eq. 6 corresponding to vertical heat transport is determined by the Nusselt number, as our data had very limited resolution to test this (Sec. VII.2). If these can be confirmed, then it would suggest that low-dimensional models can include both coherent flow structures and transport in the same model, using equations such as Eqs. 2, 3, 4, 5, and 6 for flow structure interactions (including interactions with boundaries through predictable terms such as $V_{g}$), and relating them to equations such as Eq. 19 for heat transport.

While most of the features of the interaction between neighboring LSCs could be explained by an effective turbulent thermal diffusion between the cells, some stabilizing features of the co-rotating state could not be explained. This includes the stabilizing forcing on the orientation of co-rotating states in $\theta_{1}-\theta_{2}$, and the value of $\Delta\bar{T}_{0,co}/\bar{\delta}_{co}>2$, which increases to much higher values at large $\Delta T_{m}$ (Fig. 15). These observations are not consistent with turbulent thermal diffusion. It is likely that these discrepancies are due to an additional coherent advection of heat from one cell to the other, which was not included in the model.

The interaction between two neighboring rolls should be straightforwardly generalizable to systems of more rolls with the addition of similar interaction terms to Eqs. 4, 5, and 6 for each neighboring roll, assuming that there are no higher-order interaction terms that involve 3 or more rolls. In other words, systems with multiple rolls – such as the cases with large lateral extend with a series of adjacent rolls – would then have a set of equations of motion in $\theta_{i}$, $\delta_{i}$, and $T_{0,i}$ for each roll (the number of rolls would have to be known), and each equation of motion would have a set of interaction terms for each neighbor. When considering flows in different geometries, the geometry-dependence is expected to be mostly addressed by the separate term $V_{g}(\theta_{i})$ characterizing the interaction with the boundary, which was predicted in earlier papers Brown and Ahlers (2008b); Ji and Brown (2020b), while the neighboring-roll interaction terms are expected to be separate from boundary geometry effects. The model terms were derived here based on a sinusoidal temperature profile, which could in principle have higher order Fourier terms in other flows, adding even more terms to the model, although the sinusoidal temperature profile has been found to be a good approximation (within a few percent) for flows in circular cross-sections of different aspect ratios as well Brown and Ahlers (2007b); Vogt et al. (2018). The bounds of the integrals when evaluating terms of Eqs. 21, 22, and 26 depend on the length of the interaction region, but would only affect the coefficients of each model term, and not the form of the terms. These coefficients could be re-derived for an interaction region of different length following the procedure in the Appendix.

XI Acknowledgements

Preliminary experimental investigations and model analysis were done by D. Ji, which was essential to develop our understanding of the system. Experimental investigations for published data, analysis, and writing of the manuscript was done by E. Brown. We thank the University of California, Santa Barbara machine shop and K. Faysal for helping with construction of the experimental apparatus. This work was supported by Grant CBET-1255541 of the U.S. National Science Foundation.

Appendix: model derivations

XI.1 Forcing terms $\dot{\delta}_{\kappa,i}$ and $\ddot{\theta}_{\kappa,i}$ due to neighboring roll interactions

We start with the assumption that Eqs. 2 and 3 describe a baseline model for each LSC, and that we only need to derive terms for additional forcing due to neighboring-roll interactions. The forcing on the temperature profile along the interface between neighboring convection rolls is assumed to come from the turbulent diffusion of heat in across the interface with effective turbulent thermal diffusivity $\kappa_{t}$, which is assumed to be uniform and constant. The turbulent thermal diffusion represents the enhancement of heat transport by eddies relative to the thermal diffusivity $\kappa$ due to advection, analogous to a turbulent viscosity or eddy viscosity. This thermal diffusion is approximated from a heat equation for heat flux only in the direction perpendicular to the interface between neighboring cells, as a function of the coordinate $\theta$ along the interface, as heat transport in other directions should already be accounted for in the existing Eqs. 2 and 3. We estimate the thermal diffusion term as $\kappa_{t}\nabla^{2}T\approx\kappa_{t}[T_{1}(\theta,\theta_{1})-T_{2}(\theta,\theta_{2})]/L^{2}$, where $T_{1}(\theta,\theta_{1})-T_{2}(\theta,\theta_{2})$ is the temperature difference between the neighboring cells over an interfacial mixing layer of unknown thickness. We include the unknown mixing layer thickness in the value of $\kappa_{t}$, which will be a fit parameter, so we can use the known cell size $L$ as the lengthscale in the thermal diffusion equation. This formulation is mathematically analogous to boundary layer approximations with thermal diffusion in an interfacial mixing layer of thickness $L\kappa/\kappa_{t}$ between the two cells, although it seems unlikely that such a thin mixing layer exists in this case. Because of the middle insulating wall blocking half of the interface, the heat transport acts over the exposed fractional area $A=0.50$ of the interface between the two cells. An additional vertical heat transport from each of the top and bottom plates is calculated using the standard boundary layer approximation that the heat transport from each plate at temperature $T_{m,i}\pm\Delta T/2$ to the bulk at temperature $T_{0,i}$ is diffusive across thermal boundary layers of thickness $H/2\mbox{Nu}$, where the Nusselt number Nu is the dimensionless turbulent heat transport in the vertical direction. In the vertical heat transport term we ignore the variation in bulk temperature with $\theta$, which is a small correction since $\delta\ll\Delta T$. The resulting simplification of the heat equation for cell 1 is

[TABLE]

For brevity, we write all equations in this section for the forcing on cell 1 only, as the equations for cell 2 are identical other than an exchange of the subscripts 1 and 2 in each equation.

Equation 19 is a function of the coordinate $\theta$ along the interface between cells. To instead put it in terms of parameters $\theta_{i}$, $\delta_{i}$, and $T_{0,i}$ of the LSC, we need to transform the equation into that parameter space. To do this, we make use of the time derivative of the temperature profile of Eq. 1:

[TABLE]

In this equation, the lowest order terms of a Fourier series of $\dot{T}_{1}(\theta)$ relate directly to the forcing terms $\dot{\theta}_{1}$ and $\dot{\delta}_{1}$. We assume for simplicity that these lowest order Fourier series terms are a good approximation of the LSC – i.e. the only parameters needed to describe the LSC are $\theta_{i}$, $\delta_{i}$, and $T_{0,i}$ from the sinusoidal temperature profile of Eq. 1, and ignore the small higher order terms in the Fourier series Brown and Ahlers (2007b); Vogt et al. (2018); Ji and Brown (2020a). Thus, we perform a Fourier transform of Eq. 20 in terms of sines and cosines to obtain expressions for forcing terms on the LSC parameters

[TABLE]

and

[TABLE]

Equations 21 and 22 are evaluated by inserting Eq. 19 for $\dot{T}_{1}(\theta)$, which uses as input Eq. 1 for $T_{1}(\theta,\theta_{1})$ and $T_{2}(\theta,\theta_{2})$. Expanding the term $\cos(\theta-\theta_{2})=\cos(\theta-\theta_{1})\cos(\theta_{1}-\theta_{2})-\sin(\theta-\theta_{1})\sin(\theta_{1}-\theta_{2})$ from Eq. 19 makes the integrands of Eqs. 21 and 22 functions only of $\theta-\theta_{1}$. The neighboring roll interaction terms along the interface only need to be integrated over the range $-\pi/4<\theta<\pi/4$ rad, while the second term of Eq. 19 due to vertical heat transport drops out when integrated over $\theta$ from 0 to $2\pi$ rad. Integrating Eqs. 21 and 22 over $\theta$ results in:

[TABLE]

and

[TABLE]

This forcing $\dot{\delta}_{\kappa,1}$ can be inserted directly to the existing equation of motion (Eq. 2). Since Eq. 3 is in terms of $\ddot{\theta}_{1}$, but we derived a forcing in $\dot{\theta}_{1}$ (Eq. 23), we assume the forcing in $\dot{\theta}_{1}$ can be treated as the forcing in the overdamped limit, corresponding to an offset in the $\dot{\theta}_{1}$-term of Eq. 3, so that

[TABLE]

XI.2 Equation of motion for $T_{0,1}$

In an equation of motion for $T_{0,i}$, we should also expect to have some diffusive fluctuations driven by turbulence, analogous to Eqs. 2 and 3, so we include a fluctuation term $f_{T}(t)$ with diffusivity $D_{T}$. The deterministic part of the equation for $T_{0,1}$ is assumed to be due to the net heat flux between the cells and from the top and bottom plates to the cells, calculated by integrating both sides of Eq. 19 over $d\theta/\cos^{2}\theta$. The denominator of $\cos^{2}\theta$ accounts for the relative interface area at each infinitesimal $d\theta$ to satisfy conservation of energy. The term of Eq. 19 with $A$ due to the neighboring roll interaction only contributes in the interaction region $|\theta|<\pi/4$ rad so is only integrated over that range. This results in:

[TABLE]

Using the same substitution as before for $\cos(\theta-\theta_{2})$, and integrating, results in

[TABLE]

For simplicity of presenting equations for our cubic cells, we present equations in the main paper for aspect ratio 1 ($H=L$).

XI.3 Stable fixed points for counter-rotating states for $\Delta T_{m}=0$

In a counter-rotating state where $\theta_{1}=\theta_{2}$, the equations for the two cells are symmetric in cell 1 and cell 2, so we expect a fixed point when $\delta_{1}=\delta_{2}$ and $T_{0,1}=T_{0,2}=T_{m,1}=T_{m,2}$ in the Boussinesq approximation. Inspection of Eqs. 23, 24 and 27 shows that each term is zero in this case, so $\dot{\theta}_{i}=0$, $\dot{\delta}_{i}=0$ and $\dot{T}_{i}=0$, corresponding to a fixed point for counter-rotating states. Since Eq. 4 is a forcing on $\delta_{i}$ in addition to Eq. 2, the zero forcing means that $\delta_{i}$ has the same stable fixed point $\delta_{0}$ as Eq. 2.

Is this fixed point stable? The first term of Eq. 23 has a restoring force proportional to $-\sin(\theta_{1}-\theta_{2})$, where the proportionality remains negative for all possible values of $\pi/4-\cos(2\theta_{1})/2$, which is stable when $\theta_{1}=\theta_{2}$. A Taylor expansion of the second term of Eq. 23 in $\theta_{1}-\theta_{2}$ has second-order terms, but no first order terms, so is linearly stable in $\theta_{1}-\theta_{2}$. The third term of Eq. 23 is zero in the symmetric case when $T_{0,1}=T_{0,2}$, and the fourth term is zero at a fixed point where $\dot{T}_{0,i}=0$. Thus, Eq. 23 is linearly stable in $\theta_{1}-\theta_{2}$, and counter-rotating states are a stable fixed point solution of this equation. The first term of Eq. 23 which is responsible for this stability represents the heat transfer from effective turbulent diffusion from the LSC of roll 2 projected onto an orientation perpendicular to roll 1 to drive changes in orientation.

XI.4 Stable fixed points for co-rotating states for $\Delta T_{m}=0$

Observed stable fixed points for co-rotating states correspond to $\theta_{1}=0$ and $\theta_{2}=\pi$ rad or $\theta_{1}=\pi$ rad and $\theta_{2}=0$, $\delta_{1}=\delta_{2}$, and $\dot{T}_{0,i}=0$ in steady state. All four terms of Eq. 23 are zero in in this case, confirming these are fixed points of Eq. 23 as well. There are forcings on $\dot{\delta}_{i}$ in Eq. 4 and $\dot{T}_{0,i}$ in Eq. 6 that shift their stable fixed point values, which are addressed in Sec. VI.4. Assuming constant $\Delta\bar{T}_{0,co}\equiv\bar{T}_{0,2}-\bar{T}_{0,1}$ and $\bar{\delta}_{co}\equiv\bar{\delta}_{1}=\bar{\delta}_{2}$ in the co-rotating state, a linear approximation around $\theta_{p,1}=0$ or $\pi$ rad for co-rotating states of Eq. 23 plugged into Eq. 25 is

[TABLE]

While the first term of Eq. 23 is what stabilized the counter-rotating state, for the co-rotating state where $\theta_{1}=\theta_{2}+\pi$ rad, this forcing instead is predicted to be destabilizing in the first term of Eq. 28. The second term of Eq. 23 results in a restoring force proportional to $+\sin(2\theta_{1})$ which also destabilizes co-rotating states where $\theta_{1}=0$ or $\pi$ rad. The third term of Eq. 23 results in a restoring force proportional to $\Delta\bar{T}_{0,co}\sin\theta_{1}$, which is stabilizing for $\theta_{1}=0$ and $\Delta\bar{T}_{0,co}>0$, or $\theta_{1}=\pi$ rad and $\Delta\bar{T}_{0,co}<0$. These combinations correspond to same combinations as observed co-rotating states. This stabilizing term represents the heat transfer by effective turbulent thermal diffusion from the difference in mean temperatures $\Delta\bar{T}_{0,co}$, which provides a uniform heat flux across the interface, which when projected onto the LSC temperature profile, tends to align the LSC of cell 1 with $\theta_{1}=0$ if cell 2 is hotter ($\Delta\bar{T}_{0,co}>0$) or $\theta_{1}=\pi$ rad if cell 2 is colder ($\Delta\bar{T}_{0,co}<0$). A stable co-rotating state requires the third term of Eq. 23 to be larger than the second term (which are combined in the second term of Eq. 8), and the first term of Eq. 23, as well as the forcing from the geometric potential $-\nabla V_{g}(\theta_{i})$ which is destabilizing at $\theta_{i}=0$ or $\pi$ rad Ji and Brown (2020b). Whether $\Delta\bar{T}_{0,co}/\bar{\delta}_{co}$ is large enough to make the third term of Eq. 23 dominant over all other terms so the model has stable fixed points for co-rotating states is tested in Sec. VII.5.