1. Introduction
Thermal convection and the associated large-scale circulation (LSC) play an instrumental role in applications diverse as atmospheric and oceanic flows (Salmon Reference Salmon1998; Zhong, Funfschilling & Ahlers Reference Zhong, Funfschilling and Ahlers2009), mantle and liquid-core convection (Whitehead Reference Whitehead1972; Zhang & Libchaber Reference Zhang and Libchaber2000; Zhong & Zhang Reference Zhong and Zhang2005; Whitehead & Behn Reference Whitehead and Behn2015; Huang et al. Reference Huang, Zhong, Zhang and Mertz2018) and solar magneto-hydrodynamics (de Wit et al. Reference de Wit2020). In these settings, it is known that the LSC can spontaneously reverse direction (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009), manifesting, for example, as a sudden change in wind direction (van Doorn et al. Reference van Doorn, Dhruva, Sreenivasan and Cassella2000) or potentially a reversal of the Earth's magnetic dipole (Glatzmaier et al. Reference Glatzmaier, Coe, Hongre and Roberts1999).
Large-scale circulation reversals have been observed in controlled laboratory experiments (Creveling et al. Reference Creveling, Paz, Baladi and Schoenhals1975; Gorman, Widmann & Robbins Reference Gorman, Widmann and Robbins1984, Reference Gorman, Widmann and Robbins1986; Castaing et al. Reference Castaing, Gunaratne, Heslot, Kadanoff, Libchaber, Thomae, Wu, Zaleski and Zanetti1989; Brown & Ahlers Reference Brown and Ahlers2007; Xi & Xia Reference Xi and Xia2007; Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Song, Villermaux & Tong Reference Song, Villermaux and Tong2011; Wang et al. Reference Wang, Lai, Song and Tong2018; Chen et al. Reference Chen, Huang, Xia and Xi2019) and numerical simulations (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Xu, Chen & Xi Reference Xu, Chen and Xi2021), where a progressive increase of the Rayleigh number 
$Ra$ triggers a sequence of transitions. Depending on underlying conditions, this sequence can include: (i) the onset of fluid motion giving rise to steady circulation; (ii) the destabilization of this circulatory state giving rise to chaotic reversals of the LSC; (iii) a return to order at high 
${Ra}$ in which LSC reversals recur periodically despite small-scale turbulence.
Despite much progress, LSC reversals remain poorly understood. Current theory can be broadly categorized as application of the Lorenz equations or phenomenological models. The Lorenz equations, originally derived in the context of planar upper and lower boundaries with unbounded horizontal periodicity (Lorenz Reference Lorenz1963), captures many of the transitions listed above. However, when applied to the finite geometries accessible to experiments, the Lorenz system only describes the spatially averaged flow (Welander Reference Welander1967; Gorman et al. Reference Gorman, Widmann and Robbins1986; Tritton Reference Tritton1988; Widmann, Gorman & Robbins Reference Widmann, Gorman and Robbins1989; Ehrhard & Müller Reference Ehrhard and Müller1990; Singer, Wang & Bau Reference Singer, Wang and Bau1991), resulting in substantial quantitative differences with experiments (Gorman et al. Reference Gorman, Widmann and Robbins1986). More recent phenomenological models account for additional physical effects, such as detached thermal plumes or corner rolls, by supplementing fundamental conservation laws with nonlinear or stochastic terms (Araujo, Grossmann & Lohse Reference Araujo, Grossmann and Lohse2005; Brown & Ahlers Reference Brown and Ahlers2007; Ni, Huang & Xia Reference Ni, Huang and Xia2015). While these models lend great physical insight, the connection to first principles may not be self-evident due to the ad hoc nature of the extra terms. In rectangular and cylindrical domains, it has been suggested that corner vortices and the associated turbulent fluctuations can perturb the LSC structure, causing it to switch between two bistable states (Brown & Ahlers Reference Brown and Ahlers2007; Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010). Switching may occur irregularly, partly due to intermittent heat accumulation and release (Wang et al. Reference Wang, Lai, Song and Tong2018).
In this article, we discuss one physical scenario in which a first-principled and precise understanding of LSC reversals can be gained. The scenario is thermal convection in a narrow annulus, closely related to a so-called closed-loop thermosyphon (Welander Reference Welander1967; Gorman et al. Reference Gorman, Widmann and Robbins1986; Tritton Reference Tritton1988; Widmann et al. Reference Widmann, Gorman and Robbins1989; Ehrhard & Müller Reference Ehrhard and Müller1990; Singer et al. Reference Singer, Wang and Bau1991; Basu, Bhattacharyya & Das Reference Basu, Bhattacharyya and Das2014). The annular geometry reinforces the dominant circular structure of natural convection, while eliminating corner-induced effects that tend to introduce greater complexity and tend to be geometric specific. The confined nature of the annular flow is more amenable to low-dimensional characterization, while also exhibiting sufficiently complex behaviour, such as chaotic and periodic LSC reversals, to be useful for understanding convection. The results of our study complement those conducted in rectangular geometries, where domain corners impact convective behaviour (Brown & Ahlers Reference Brown and Ahlers2007; Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Ni et al. Reference Ni, Huang and Xia2015; Chen et al. Reference Chen, Huang, Xia and Xi2019).
Rather than beginning with the Lorenz equations and making adjustments to suit the annulus, we derive a low-order system directly from a Fourier–Laurent truncation of the governing Navier–Stokes Boussinesq (NSB) equations. The resulting system resembles the Lorenz equations, but differs in a few important ways. Notably, the Laurent expansion accounts for spatial dependence of the flow (see also Yorke, Yorke & Mallet-Paret Reference Yorke, Yorke and Mallet-Paret1987), permitting exact enforcement of boundary conditions without the need for empirically estimated friction or heat-transfer coefficients. Comparison with direct numerical simulations (DNS) shows this system captures the entire sequence of transitions, including regimes of chaotic and periodic LSC reversals.
Importantly, we show equivalence between this low-order system and an externally driven mechanical pendulum. Two length scales naturally emerge from this correspondence: the fulcrum 
$y_1$ of the pendulum and the point 
$y_0$ towards which the fluid centre of mass (CoM) is externally driven, both given by explicit formulas. Knowledge of these quantities yields a simple formula for the frequency of periodic LSC reversals in the high-
${Ra}$ regime, shown to be in excellent agreement with DNS. Further, this correspondence reveals a transparent mechanism for reversals, namely the inertial overshoot of the fluid CoM as equivalent to a forced pendulum. The clean characterization afforded by annular convection may offer a new point of approach for understanding convection in other geometries.
2. Convective states revealed by direct numerical simulations
Figure 1(a) depicts the problem setup in which a two-dimensional annular fluid domain is heated from below. Thermal exchange occurs along the outer boundary with an imposed temperature that decreases linearly with height, while the inner boundary remains adiabatic. Dimensionless temperature 
$T$, velocity 
${\boldsymbol {u}}$ and pressure 
$p$ fields are governed by the incompressible NSB equations
$$\begin{gather} \frac{\partial {\boldsymbol{u}}}{ \partial t } + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} ={-}\boldsymbol{\nabla} p + {Pr}\,\nabla^2 {\boldsymbol{u}} + {Ra}\, {Pr} \, T {\boldsymbol{e_y}}, \end{gather}$$
$$\begin{gather}\frac{\partial T}{ \partial t } + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \nabla^2 T, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}} = 0 , \end{gather}$$
holding in the dimensionless annulus, 
$r_0< r<1/2$, where 
$r_0$ is the radius of the inner boundary shown in figure 1(a). Both the inner and outer rings are no-slip boundaries. Parameters include the Rayleigh number 
${Ra} = \beta _T \Delta T h^3 g/(\nu \kappa )$ and Prandtl number 
${Pr} = \nu /\kappa$ (Appendix A), where 
$h{}$ is the dimensional outer boundary diameter and 
$\Delta T$ is the temperature difference between the bottom and top points of the outer boundary. Other physical parameters include 
$\beta _T$, 
$g$, 
$\nu$, 
$\kappa$, which are the thermal expansion coefficient, acceleration due to gravity, kinematic viscosity and thermal diffusivity, respectively. When 
${Ra}$ is sufficiently high, the destabilizing action of buoyancy can give rise to natural convection.
Direct numerical simulations of natural convection in an annulus. (a) Schematic of an annular fluid domain heated from below. (b) At low 
$Ra$ (
$3.9\times 10^5$), the conductive state is stable, resulting in a raised CoM (green dot). Any initial angular momentum quickly dissipates as shown in the plot of 
$L(t)$ below. (c) At higher 
$Ra$ (
$3.1\times 10^6$) the system transitions to steady circulation with offset CoM and non-zero 
$L$. (d) At yet higher 
$Ra$ (
$5\times 10^7$), the LSC can spontaneously reverse direction. The fluid CoM wanders erratically (green trajectory) and 
$L(t)$ reverses chaotically. (e) At the highest 
$Ra$ (
$1.6\times 10^9$), the LSC reversals recur periodically, even though the small-scale flow is turbulent. ( f) The temperature power spectrum of (e) peaks at frequency 
$f^*$, corresponding to the LSC reversal frequency. At higher frequency, the decay rate is consistent with a 
$-1.4$ power law. See supplementary movies (b)–(e) available at https://doi.org/10.1017/jfm.2024.584. In all cases, 
${Pr} = 4$ and 
$r_0 = 0.4$.
To quantify different convective states, we will examine the spatially averaged fluid angular momentum 
$L(t)$ and the fluid CoM coordinates 
$(X(t), Y(t))$, given by
$$\begin{gather} L(t) = \frac{1}{\text{A}_0}\int_{\varOmega} r u \, {\rm d}A , \end{gather}$$
$$\begin{gather}X(t) ={-} \frac{1}{\text{A}_0} \int_{\varOmega} x T \, {\rm d}A , \end{gather}$$
$$\begin{gather}Y(t) ={-} \frac{1}{\text{A}_0} \int_{\varOmega} y T \, {\rm d}A . \end{gather}$$
Here, 
$\text {A}_0 = {\rm \pi}(1-4r_0^2)/4$ is the area of the annulus 
${\varOmega }$ and 
${\rm d}A = r\,{\rm d}r \,{\rm d}\theta$ is the area element. The fluid CoM coordinates above are expressed in terms of the temperature field owing to the fact that fluid density varies as the negative of 
$T$. We note that 
$L>0$ corresponds to a counter-clockwise rotating flow.
The range of convective states are revealed by DNS of the NSB system, as shown in figure 1. Simulations are based on a Chebyshev–Fourier pseudo-spectral discretization (Appendix B) of (2.1) to (2.3) in streamfunction–vorticity form with implicit–explicit time stepping (Peyret Reference Peyret2002; Huang, Shelley & Stein Reference Huang, Shelley and Stein2021; Huang & Zhang Reference Huang and Zhang2023). At low 
${Ra}$, figure 1(b) shows the existence of a stable conductive state with no fluid motion and with raised CoM (green dot). In this regime, perturbations to the conductive state decay rapidly, as seen in the plot below showing 
$L(t) \to 0$. Increasing 
${Ra}$ eventually destabilizes the system, leading to the state shown in figure 1(c), where the fluid circulates either clockwise (CW) or counterclockwise (CCW) at a constant rate. In this state, the fluid CoM is fixed and offset from centre. By further increasing 
${Ra}$, this steady circulating state also destabilizes; the direction of circulation now alternates over time and the flow reverses chaotically, as shown in the time series of 
$L$ in figure 1(d). The fluid CoM wanders erratically in this regime. Interestingly, large-scale chaos disappears when 
${Ra}$ becomes sufficiently high, and figure 1(e) reveals an oscillating state with periodic LSC reversals. Here, the oscillatory CoM trajectory resembles pendulum motion. Although the reversals are periodic, the small-scale flow is turbulent and resolved by the DNS. The turbulent fluctuations are characterized by the frequency power spectrum of the temperature field, shown in figure 1( f) to follow the turbulent Bolgiano–Obukhov power law of natural convection (Wu et al. Reference Wu, Kadanoff, Libchaber and Sano1990; Lohse & Xia Reference Lohse and Xia2010).
3. Low-order model for LSC reversals
All of these states can be recovered by a low-dimensional system that arises directly from the governing NSB equations. As further detailed in Appendix C, the brief derivation is as follows. In polar coordinates, 
${\boldsymbol {u}} = u(r,\theta,t) {\boldsymbol {e_\theta }} + v(r,\theta,t) {\boldsymbol {e_r}}$ and 
$T = T(r,\theta,t)$, consider a Fourier expansion in 
$\theta$ and a Laurent expansion in 
$r$, and truncate each to a desired order while enforcing all boundary conditions (BCs). The choice of Laurent expansion is guided by the form of the conductive-state solution and recovers this basic state with no approximation made. Inserting the truncated variables into (2.1) to (2.3) and projecting onto the Fourier–Laurent basis yields a finite-dimensional system.
Truncating the combined Fourier–Laurent expansion at the lowest order capable of satisfying all BCs and casting in terms of the physically relevant variables 
$L(t), X(t), Y(t)$ produces the dynamical system
$$\begin{gather} \dot{L} ={-} {Ra} \,{Pr} \, X - \alpha \,{Pr} \, L , \end{gather}$$
$$\begin{gather}\dot{X} ={-}k L (Y - y_1) - \beta X , \end{gather}$$
$$\begin{gather}\dot{Y} = k L X - \beta (Y - y_0) . \end{gather}$$
This system governs the evolution of the fluid angular momentum 
$L(t)$ and fluid CoM coordinates 
$(X(t), Y(t))$ defined in (2.4) to (2.6). The coefficients 
$\alpha, \beta, k, y_0, y_1 >0$ are purely geometric in that they depend on 
$r_0$ only. Formulas for these coefficients are given in Appendix C. Although no assumption is made on the width of the annulus, the low-order truncation is most accurate for a relatively narrow annulus. We therefore set 
$r_0 = 0.4$ in all numerical examples. Higher-order truncations could be used for a wider annulus.
The above system exhibits the same quadratic nonlinearity as the Lorenz equations. The parameter structure, however, arises directly from the annular geometry and differs from that of Lorenz. Differences, therefore, exist in the parameter regimes accessible by each system and especially in the threshold values separating different states.
While the comparison with the Lorenz equations can be illuminating, even more physical insight can be gained by recognizing how (3.1) to (3.3) relate to a mechanical pendulum. In particular, if one artificially sets 
$\beta = 0$, then (3.1) to (3.3) are identical to those of a pendulum with fulcrum 
$y_1$, length 
$l = \sqrt {X^2 + (Y-y_1)^2}$, effective gravitational constant 
$g_{{eff}} = k l^2 \, {Ra} \, {Pr}$ and damping rate 
$\alpha {Pr}$. The system is simply written in terms of the pendulum CoM and angular momentum rather than the more familiar angular displacement. In this equivalence, the pendulum CoM corresponds exactly to the fluid CoM 
$(X,Y)$ as defined in (2.5) and (2.6), and the pendulum angular momentum corresponds exactly to the fluid angular momentum 
$L$ as defined in (2.4). Note that the effective gravitational constant 
$g_{{eff}}$ of the pendulum system is unrelated to the actual gravitational constant 
$g$ of the thermal convection system. Also note that if 
$\beta \ne 0$, the two additional driving terms present in (3.2) and (3.3) can cause the pendulum length 
$l = \sqrt {X^2 + (Y-y_1)^2}$ to vary dynamically, opening the possibility of chaotic dynamics.
The driving terms involving 
$\beta$ arise from the interaction of boundary heating and buoyancy. As seen in (3.1) to (3.3), these terms drive the CoM towards the point 
$(0,y_0)$, which Appendix A shows is the CoM of the conductive state. Thus, 
$(L,X,Y)=(0,0,y_0)$ corresponds to the conductive-state solution that is given explicitly by (A4) and depicted in figure 1(b). Stability analysis discussed in § 4 shows that this state is stable up to a critical Rayleigh number.
The most important parameters in the pendulum correspondence of (3.1) to (3.3) are 
$y_0$ and 
$y_1$, corresponding to the height of the conductive-state CoM and the pendulum fulcrum, respectively. Appendices A and C give explicit formulas for these two length scales in terms of the geometric parameter 
$r_0$. When the fluid CoM lies at the pendulum fulcrum, 
$(X,Y) = (0,y_1)$, the restoring torque in (3.1) to (3.3) vanishes but the driving terms that push 
$(X,Y)$ towards 
$(0,y_0)$ do not vanish. The state 
$(L,X,Y) = (0,0,y_1)$ is therefore not an equilibrium of the system due to the continual injection of thermal energy from the boundary. Appendix C further shows that 
$y_0 > y_1 > 0$ for any choice of 
$r_0$, implying that the thermal injection always acts to raise the CoM above the fulcrum and, hence, tends to destabilize the system.
4. Bifurcations and comparison with DNS
How well do these simple ordinary differential equations (ODE) describe the dynamics of convection? Figure 2 shows trajectories of 
$(L,X,Y)$ computed by fully resolved DNS (top) vs those computed by the ODE model (bottom) for the same Rayleigh numbers as figure 1(c–e). Figure 2(a–c) shows that the trajectories from DNS and the ODE model are remarkably similar across the range of 
${Ra}$, exhibiting (a) convergence to a stable circulating state, (b) chaotic dynamics near a strange attractor and (c) periodic orbits at the highest 
${Ra}$. The trajectories in figure 2(b,c) indicate reversals of the LSC, as can be seen by the sign change of 
$L$. The LSC reversals are chaotic in figure 2(b) and periodic in figure 2(c).
Trajectories of ODE system (3.1)–(3.3) in comparison with fully resolved DNS. The trajectories of 
$(L,X,Y)$ are remarkably similar across the range of Rayleigh numbers, showing (a) convergence to a stable circulating state for 
${Ra} = 3.1\times 10^6$, (b) strange-attractor dynamics for 
${Ra} = 5\times 10^7$ and (c) periodic dynamics for 
${Ra} = 1.1\times 10^9$. In all cases, 
${Pr} = 4$ and 
$r_0 = 0.4$.
The bifurcation diagram in figure 3 shows that a pitchfork bifurcation occurs at a critical value 
${Ra}_1^*$. At this value, the conductive state loses stability, and, simultaneously, the bistable circulating states appear (CW and CCW circulation). At a second critical value, 
${Ra}_2^*$, these circulating states lose stability through a Hopf bifurcation. Immediately past 
${Ra}_2^*$, the dynamics is fractal-like and chaotic, characteristic of a strange attractor. These observations are further supported by measurements of the fractal dimension 
$D_2$ (Ott Reference Ott2002) and Lyapunov exponent 
$\lambda$ shown in the inset. At much higher 
${Ra}$, order reemerges and trajectories resemble the arc-like path of a pendulum.
Bifurcation diagram shows a pitchfork bifurcation at 
${Ra}_1^*$ and a Hopf bifurcation at 
${Ra}_2^*$. Trajectories corresponding to figure 2(a–c) are marked on the diagram; for all trajectories, 
${Pr} = 4$ and 
$r_0 = 0.4$. Inset: the fractal dimension 
$D_2$ and Lyapunov exponent 
$\lambda$ distinguish chaotic states from orderly ones.
The ODE model yields exact formulas for both critical values (Appendix D)
\begin{equation} {Ra}_1^* = \frac{\alpha\beta}{k\Delta y}, \quad {Ra}_2^* = \frac{\alpha^2 \, {Pr}}{k \Delta y} \left( \frac{\alpha \,{Pr} + 4\beta}{\alpha \,{Pr} - 2\beta}\right), \end{equation}
where 
$\Delta y = y_0-y_1>0$ is the distance between the conductive-state CoM and the pendulum fulcrum. Briefly, the value 
${Ra}_1^*$ is found through linear stability analysis of the conductive state 
$(L,X,Y) = (0,0,y_0)$. As 
${Ra}$ crosses 
${Ra}_1^*$, the conductive state loses stability and the circulating states appear. Immediately past 
${Ra}_1^*$, the Jacobian of each circulating state possesses three real, negative eigenvalues. As 
${Ra}$ increases further, two eigenvalues become complex, 
$z_{2,3} = \sigma \pm {\rm i} \omega$, with 
$\sigma <0$ initially. As 
${Ra}$ crosses 
${Ra}_2^*$, 
$\sigma$ becomes positive and thus the circulating states lose stability, giving way to the strange attractor seen in figure 2(b). This analysis is similar to that conducted for the Lorenz equations (Welander Reference Welander1967; Creveling et al. Reference Creveling, Paz, Baladi and Schoenhals1975; Gorman et al. Reference Gorman, Widmann and Robbins1986; Ehrhard & Müller Reference Ehrhard and Müller1990), the main difference being that modelling choices made early on (e.g. accounting for the flow's spatial dependence) enable greater accuracy in predicting the bifurcation parameters than obtained previously (Gorman et al. Reference Gorman, Widmann and Robbins1986).
The phase diagram in figure 4 gives a bird’s eye view of the different convective states. In the figure, coloured dots correspond to fully resolved DNS, showing regions of a stable conductive state (blue), bistable circulating states (green) and LSC reversals, both chaotic (orange) and periodic (red). The boundaries between these regions are well predicted by the formulas for 
${Ra}_1^*$ and 
${Ra}_2^*$ given in (4.1). In particular, 
${Ra}_1^*$ is independent of the Prandtl number, giving the vertical green line. The predicted value of 
${Ra}_1^*$ agrees with DNS in all cases to within the grid resolution of figure 4. Meanwhile, the orange curve shows the 
${Pr}$ dependence of 
${Ra}_2^*$. For the important case of water, 
$4 < {Pr} < 8$, this threshold is also predicted to within the grid resolution. While some discrepancy is visible for other values of 
${Pr}$, the curve captures the qualitative shape of the boundary, in particular the horizontal asymptote 
${Pr}^* = 2 \beta /\alpha$ obtained by setting the denominator of 
${Ra}_2^*$ equal to zero. We note that, even though it is a very different geometry, thermal convection in a rectangular domain yields a phase diagram with the same states and with similar orders of magnitude for the thresholds (Araujo et al. Reference Araujo, Grossmann and Lohse2005).
Phase diagram of different convective states. Coloured dots are from DNS, where blue indicates a stable conductive state, green indicates bistable circulating states, orange indicates chaotic LSC reversals and red indicates periodic LSC reversals. Formulas for 
${Ra}_1^*$ and 
${Ra}_2^*$ from the ODE model predict the boundaries between the regions well.
5. Periodic LSC reversals at high 
${Ra}$
As 
${Ra}$ increases well beyond 
${Ra}_2^*$, large-scale chaos subsides and gives way to the nearly periodic trajectories seen in figure 2(c). While the bifurcations discussed in the previous section have been qualitatively described by related models, the periodic regime has received less attention and, so far, has resisted clean characterization. It is precisely this regime where the novel mechanical-pendulum correspondence becomes most valuable.
As seen in the figure 3 inset, the return to order at high 
${Ra}$ is indicated by the fractal dimension dropping to one and the Lyapunov exponent dropping to zero at the same Rayleigh number, roughly 
${Ra} = 10^9$. At this value, a stable limit cycle emerges in the ODE system, giving CoM orbits that resemble pendulum motion. Figure 5(a) shows four such orbits of the fluid CoM 
$(X(t),Y(t))$ for Rayleigh numbers in the range 
$1/4 \times 10^{10} < {Ra} < 16 \times 10^{10}$. At the lowest 
${Ra}$, the corresponding pendulum length 
$l(t) = \sqrt {X^2 + (Y-y_1)^2}$ varies somewhat over the period. At higher 
${Ra}$, though, the orbit tightens and 
$l$ remains nearly constant throughout. Recall that, even though the large-scale dynamics of the fluid angular momentum 
$L(t)$ and CoM 
$(X(t),Y(t))$ is regular in this regime, the DNS shows that turbulent fluctuations inhabit the small scales (see figure 1e).
At very high 
${Ra}$, order reemerges and the large-scale dynamics becomes periodic. (a) The fluid CoM follows the swinging motion of a pendulum about fulcrum 
$(0,y_1)$. (b) The frequency of LSC reversals is well predicted by (5.1) for high 
${Ra}$. In all cases, 
${Pr} = 4$ and 
$r_0 = 0.4$.
Each swing of the pendulum seen in figure 5(a) corresponds to a sign change of the fluid angular momentum 
$L$ and, therefore, a reversal of the LSC. This simple observation offers a way to predict the dominant frequency 
$f^*$ of LSC reversals. In the pendulum correspondence of (3.1) to (3.3), the gravitational constant is 
$g_{{eff}} = k l^2 \, {Ra} \, {Pr}$, giving a small-amplitude frequency of 
$\sqrt {k l \, {Ra} \, {Pr}}/(2{\rm \pi} )$. The amplitudes seen in figure 5(a), however, are not small, implying that the frequency depends on both the pendulum length 
$l$ and the maximum swing angle 
$\phi _{{max}}$.
As detailed in the Appendix E, both of these quantities can be estimated through an energy balance with 
$E_{{eff}} = \frac {1}{2} k L^2 + {Ra}\, {Pr} (Y-y_1)$ representing the sum of kinetic and potential energy of the mechanical-pendulum system. This quantity is an effective energy of the pendulum system and does not directly represent the actual energy of the thermal convection system. Although 
$E_{{eff}}$ is not necessarily conserved over the dynamics, it does satisfy a precise energy law with energy injection (from boundary heating) and dissipation (from fluid viscosity). In the case of a limit cycle, the total energy injected must balance that dissipated, and the period-averaged 
$E_{{eff}}$ is conserved. With a few additional approximations made, this principle allows one to solve for pendulum length and maximum swing angle that arise in the case of a limit cycle. The resulting values of 
$l$ and 
$\phi _{{max}}$ depend on geometry and 
${Pr}$, but not on 
${Ra}$ (see Appendix E). With these values known, the (dimensionless) frequency of high-
${Ra}$ LSC reversals is given by
\begin{equation} f^* = \frac{\sqrt{kl \, {Ra} \, {Pr}}}{ 4 K(\sin^2 ({\phi_{{max}}}/{2})) }, \end{equation}
where 
$K$ is the complete elliptic integral of the first kind.
Figure 5(b) shows a comparison between this simple formula and the reversal frequency measured in the fully resolved DNS. In the DNS, the reversal frequency is measured as the peak location in the temperature power spectrum (see figure 1f). Figure 5(b) shows that (5.1) predicts the reversal frequency measured in DNS remarkably well over the largest decade of 
${Ra}$ run (roughly 
${Ra} = 2\times 10^8$ to 
$2\times 10^9$). At higher 
${Ra}$, DNS becomes computationally prohibitive but numerical solution of the ODE model is feasible, and the corresponding measurements of 
$f^*$ also agree with (5.1). The close agreement between DNS, the ODE model, and (5.1) suggests the primary mechanism for LSC reversals has been properly accounted for by the mechanical-pendulum correspondence.
In dimensional terms, (5.1) gives a reversal frequency of 
$F^* = c N^*$, where 
$N^* = \sqrt {\beta _T \Delta Tg/h}$ is the Brunt–Väisälä frequency, i.e. the inverse of the free-fall time scale, and 
$c$ is a constant that depends on geometry and 
${Pr}$, but not 
${Ra}$ (Appendix E). We note that 
$c=0.04$ for the case of figure 5, indicating roughly 25 free-fall time scales per reversal event. Other geometries may yield different values of this ratio.
6. Discussion
The low-order system given by (3.1) to (3.3) arises from systematic analysis of the governing NSB equations and has been shown to accurately describe a range of convective states in the annular domain. In contrast with related Lorenz-type models (Welander Reference Welander1967; Gorman et al. Reference Gorman, Widmann and Robbins1986; Tritton Reference Tritton1988; Widmann et al. Reference Widmann, Gorman and Robbins1989; Ehrhard & Müller Reference Ehrhard and Müller1990; Singer et al. Reference Singer, Wang and Bau1991; Araujo et al. Reference Araujo, Grossmann and Lohse2005), the Laurent expansion underlying (3.1) to (3.3) accounts for spatial dependence of the flow and temperature fields, precluding the need for empirically estimated friction or heat-transfer coefficients. This modelling choice enables greater accuracy in predicting parameter bifurcations, as demonstrated by direct comparison with fully resolved DNS. Many of these related models have been used as the foundation for control (Singer et al. Reference Singer, Wang and Bau1991; Wang, Singer & Bau Reference Wang, Singer and Bau1992), data assimilation (Harris et al. Reference Harris, Ridouane, Ridouane, Hitt and Danforth2012; Chen & Majda Reference Chen and Majda2018) and machine learning (Chen Reference Chen2020). The accuracy and conceptually transparency afforded by (3.1) to (3.3) could further such endeavours.
Importantly, this low-order system reveals a previously unrecognized pendulum structure underlying natural convection. In particular, (3.1) to (3.3) correspond to a damped pendulum with CoM driven upwards towards the conductive-state CoM. In addition to its physical elegance, this equivalence enables accurate predictions for the frequency of regular LSC reversals observed at high 
${Ra}$. Furthermore, it provides a transparent mechanism for the reversals. Just like a mechanical pendulum, inertia causes the fluid CoM to overshoot equilibrium. The CoM eventual reaches a zenith, at which point the restoring torque reverses the system's angular momentum, thereby creating a LSC reversal. The driving terms in (3.1) to (3.3) are necessary to counteract the damping from viscous dissipation; it is the interplay between these two terms that selects the effective pendulum length, the swing angle and thus the frequency of reversals.
This low-order model is closely related to the Lorenz equations; indeed, a change of variables can map (3.1) to (3.3) to a system with the same variable structure but a different parameter structure as Lorenz. A corollary of this fact is that the Lorenz system too could be usefully regarded as an externally driven mechanical pendulum. This observation may lend new insights into the study of the Lorenz equations. We note that, although periodic solutions of the Lorenz equations have been explored mathematically (Robbins Reference Robbins1979; Sparrow Reference Sparrow2012), the physical connection to LSC reversals was not made.
A benefit of the annular domain is that it accentuates the dominant circulating pattern of natural convection, while suppressing other, geometric-specific effects. Related to this fact, we have found that the annulus yields a simple law for the 
${Ra}$-dependence of the periodic LSC reversal frequency 
$f^*$, described by an explicit formula, (5.1), and corroborated by comparison with DNS. Previous theoretical analysis in a rectangular geometry suggests the power law 
$f^* \sim {Ra}^{0.44}$ (Araujo et al. Reference Araujo, Grossmann and Lohse2005), and laboratory experiments with cryogenic helium gas in a cylindrical container suggest 
$f^* \sim {Ra}^{0.71}$ (Araujo et al. Reference Araujo, Grossmann and Lohse2005). In the present case of an annulus, the modelling prediction and DNS are in agreement, both unambiguously showing a scaling of 
$f^* \sim {Ra}^{0.5}$. Therefore, annular convection may be considered an ideal ‘ground state’ that yields a precisely determined scaling law for the frequency of reversals. Perhaps future studies could build upon this model to determine how the scaling law is modified by various geometric effects.
Curiously, the scaling law 
$f^* \sim {Ra}^{0.5}$ is seen in experimental measurements of thermal convection in a disk (Song et al. Reference Song, Villermaux and Tong2011), but for the frequency of oscillations in the strength of the LSC. The period of this oscillation is much shorter than the average reversal time seen in the experiments (Wang et al. Reference Wang, Lai, Song and Tong2018), suggesting differences in the LSC reversals that occur in disk convection. In the annular domain, the inner boundary at 
$r = r_0$ serves as a confinement that regulates the flow. The recent study of Li, Chen & Xi (Reference Li, Chen and Xi2024) demonstrates this concept experimentally, as the inclusion of a central obstruction in Rayleigh–Bénard convection substantially modifies the flow structures and enhances heat transfer.
Here, we have focused on the lowest-order system capable of satisfying the BCs on the annulus, but the truncation procedure can in principle be carried out to any order. At extremely high 
$Ra$, turbulent effects (Lohse & Xia Reference Lohse and Xia2010) are associated with fine-scale structures, which could potentially be captured by retaining higher-order terms in the Fourier–Laurent expansion, either directly or through stochastic parameterization. Moreover, extension of the model into three dimensions could account for azimuthal rotations of the LSC plane, which experiments have shown take a stochastic character (Brown, Nikolaenko & Ahlers Reference Brown, Nikolaenko and Ahlers2005). Finally, we hope to couple the model to slowly moving boundaries to examine phase-change processes, such as melting or dissolution, that couple to the action of natural convection (Huang, Moore & Ristroph Reference Huang, Moore and Ristroph2015; Moore Reference Moore2017; Huang et al. Reference Huang, Tong, Shelley and Ristroph2020; Huang & Moore Reference Huang and Moore2022; Weady et al. Reference Weady, Tong, Zidovska and Ristroph2022).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.584.
Funding
The authors thank J. Zhang for useful discussions. J. M. H. acknowledges support from the National Natural Science Foundation of China (12272237, 92252204). N.J.M. acknowledges support from the National Science Foundation (DMS-2012560).
Declaration of interests
The authors report no conflict of interest.
Author contributions
N.J.M. and J.M.H. contributed equally to this work.
Appendix A. Formulation and background
To obtain the dimensionless (2.1) to (2.3), space has been rescaled on 
$h$ (the diameter of the outer annulus boundary), time on 
$h^2 / \kappa$ (the thermal diffusive time scale), velocity on 
$\kappa / h$ and density variations on 
$\Delta \rho = \rho _0 \beta _T \Delta T$. The prescribed temperature on the outer boundary, 
$r=1/2$, decreases linearly with height, while the inner boundary, 
$r=r_0$, remains adiabatic. The velocity field, expressed as 
${\boldsymbol {u}} = u {\boldsymbol {e_\theta }} + v {\boldsymbol {e_r}}$ in polar coordinates, satisfies no-slip conditions on both boundaries. The BCs are thus
$$\begin{gather} u = v = 0 \quad \text{at } r=r_0 \text{ and } r=1/2, \end{gather}$$
$$\begin{gather}\frac{\partial T}{ \partial r } = 0 \quad \text{at } r=r_0, \end{gather}$$
$$\begin{gather}T = \frac{1-\sin \theta}{2} \quad \text{at } r=1/2 . \end{gather}$$ Equations (2.1) to (2.3) support a conductive-state solution with no fluid motion 
$(u,v) = (0,0)$. The corresponding temperature field that satisfies BCs (A2)–(A3) is given by
\begin{equation} T = \frac{1}{2} - \frac{r_0}{1+4r_0^2} \left(\frac{r}{r_0}+\frac{r_0}{r}\right)\sin{\theta}. \end{equation} The fluid angular moment 
$L$ and fluid CoM coordinates (see (2.4) to (2.6)) associated with the conductive-state solution are
\begin{equation} L = 0, \quad X = 0, \quad Y = y_0 = \frac{ 1+12 r_0^2 }{16 (1+4r_0^2)}. \end{equation}
That is, 
$y_0$ represents the CoM of the conductive state. As seen above, 
$y_0 > 0$ for any value of 
$r_0$, indicating that the conductive-state CoM always lies above the annulus center.
Appendix B. Numerical methods
Equations (2.1) to (2.3) can be written in the streamfunction–vorticity form
$$\begin{gather} \frac{\partial \omega}{ \partial t } + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega = {Pr} \nabla^2 \omega + {Pr}\, {Ra}\, \left(\frac{\partial T}{\partial r}\cos\theta-\frac{1}{r}\frac{\partial T}{\partial \theta}\sin\theta\right), \end{gather}$$
$$\begin{gather}\frac{\partial T}{ \partial t } + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \nabla^2 T, \end{gather}$$
$$\begin{gather}-\nabla^2 \psi = \omega,\quad {\boldsymbol{u}} = \nabla_\perp \psi. \end{gather}$$
Where the vorticity is 
$\omega = r^{-1} [\partial _r(r u)-\partial _\theta v]$, and the streamfunction 
$\psi$ recovers the flow velocity through 
${\boldsymbol {u}} = \nabla _\perp \psi = r^{-1}\psi _\theta {\boldsymbol {e_r}} - \psi _r{\boldsymbol {e_\theta }}$.
Discretizing time with the second-order Adam–Bashforth backward differentiation method, (B1) to (B3) become
$$\begin{gather} \nabla^2 \omega^{(n)} -\sigma_1\omega^{(n)} = f^{(n)}, \end{gather}$$
$$\begin{gather}\nabla^2 T^{(n)} -\sigma_2 T^{(n)} = g^{(n)}, \end{gather}$$
$$\begin{gather}-\nabla^2 \psi^{(n)} = \omega^{(n)}, \end{gather}$$
at time step 
$t = n\Delta t$. Here
$$\begin{align} \nabla^2 &= \frac{\partial^2}{\partial r^2} + \frac{1}{r}\frac{\partial }{ \partial r } + \frac{1}{r^2} \frac{\partial^2}{\partial\theta^2}, \end{align}$$
$$\begin{align}\sigma_1 &= \frac{3}{2\,{Pr}\,\Delta t},\quad \sigma_2 = \frac{3}{2\Delta t}, \end{align}$$
\begin{align} f^{(n)} &= {Pr}^{{-}1}[ 2 ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega)^{(n-1)} - ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega)^{(n-2)}] \nonumber\\ &\quad - (2\, {Pr}\, \Delta t)^{{-}1} (4\omega^{(n-1)}-\omega^{(n-2)}) \nonumber\\ &\quad -{Ra}(T_r\cos\theta-r^{{-}1}T_\theta\sin\theta)^{(n)}, \end{align}
\begin{align} g^{(n)} &= [ 2 ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} T)^{(n-1)} - ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} T)^{(n-2)}]\nonumber\\ &\quad - (2\Delta t)^{{-}1} (4T^{(n-1)}-T^{(n-2)}) . \end{align}Explicit and nonlinear terms in 
$f^{(n)}$ and 
$g^{(n)}$ are computed pseudo-spectrally with an efficient anti-aliasing filter (Hou & Li Reference Hou and Li2007). Equations (B4) to (B6) are then solved by the Fourier–Chebyshev method detailed in Peyret (Reference Peyret2002), Huang & Zhang (Reference Huang and Zhang2023) and Huang & Moore (Reference Huang and Moore2023). There are typically 1024 Fourier modes and 128 Chebyshev nodes in each DNS of this article. The time step is 
$\Delta t = 5\times 10^{-4} \,{Ra}^{-1/2}$, considering the flow velocity 
$|{\boldsymbol {u}}|\sim \sqrt {{Ra}}$. These parameters are tested to yield resolved and accurate solutions. For the ODE model, we use the ode45 package of MATLAB.
Figure 6 shows the convergence test of the DNS scheme. In the spatial convergence test figure 6(a), the time step 
$\Delta t = 10^{-4}$ is fixed and 
$N_r = N_\theta = N$. A high-resolution solution with 
$N = 1024$ is computed first and the convergence towards this solution is tested by letting 
$N$ progressively increase. Figure 6(a) shows the results of this test and demonstrates spectral convergence. That is, the error decreases exponentially with 
$N$ until a limiting error of roughly 
$10^{-9}$ is reached. In the temporal convergence test shown in figure 6(b), 
$N_r = N_\theta = 100$ is fixed and 
$\Delta t$ is decreased by half during each test and then the error between each refinement is compared. Shown in figure 6(b), the refinement error decays as 
$ {O}(\Delta t^2)$, demonstrating a second-order temporal convergence.
Convergence of the numerical solver. (a) Spatial convergence test shows the error decays spectrally. (b) Temporal convergence test demonstrates a second-order convergence in time stepping.
Appendix C. Derivation of the ODE model
In polar coordinates, the 
$\theta$ components of (2.1) to (2.3) are given by
\begin{align} &u_t + v u_r + \frac{1}{r} u u_\theta + \frac{1}{r} u v ={-}\frac{1}{r} p_\theta + {Ra}\, {Pr} \, T \cos \theta\nonumber\\ &\quad +{Pr} \left( u_{rr} + \frac{1}{r} u_r + \frac{1}{r^2} u_{\theta \theta} - \frac{1}{r^2} u + \frac{2}{r^2} v_\theta \right) , \end{align}
$$\begin{align} & T_t + \frac{u}{r} T_\theta + v T_r = \frac{1}{r} \frac{\partial }{ \partial r } ( r T_r ) + \frac{1}{r^2} T_{\theta \theta} , \end{align}$$
$$\begin{align} & v_r + \frac{1}{r} v + \frac{1}{r} u_\theta = 0 . \end{align}$$
Multiplying (C1) by 
$r^2$, integrating over the fluid domain, applying incompressibility (C3) and the no-slip condition (A1) and using the CoM definition (2.5)–(2.6) gives
\begin{equation} \dot{L} ={-}{Ra}\, {Pr} \, X + \left.\frac{{Pr}}{\text{A}_0} \int_0^{2{\rm \pi}} (r^2 u_r )\right|_{r_0}^{r_1} \, {\rm d}\theta. \end{equation}This evolution equation for the fluid angular momentum is exact within the NSB framework.
The temperature and velocity fields 
$T(r,\theta,t)$, 
$u(r,\theta,t)$, 
$v(r,\theta,t)$ are each periodic in 
$\theta$ and so can be written as Fourier series
$$\begin{gather} T(r,\theta,t) = a_0(r,t) + \sum_{n=1}^{\infty} a_n(r,t) \cos n\theta + b_n(r,t) \sin n\theta , \end{gather}$$
$$\begin{gather}u(r,\theta,t) = \sum_{n={-}\infty}^{\infty} \hat{u}_n(r,t) \,{\rm e}^{{\rm i} n \theta}, \end{gather}$$
$$\begin{gather}v(r,\theta,t) = \sum_{n={-}\infty}^{\infty} \hat{v}_n(r,t) \,{\rm e}^{{\rm i} n \theta}. \end{gather}$$The temperature BCs (A2)–(A3) imply
$$\begin{gather} \partial_{r} a_n = \partial_{r} b_n = 0 \quad \text{at } r=r_0, \end{gather}$$
$$\begin{gather}a_0 = 1/2,\quad b_1 ={-}1/2, \quad \text{all others vanish at } r=1/2. \end{gather}$$The no-slip BC (A1) and incompressibility (C3) yield conditions
$$\begin{gather} \hat{u}_n(r,t)=\hat{v}_n(r,t) = 0 \quad \text{at } r=r_0 \text{ and } r=1/2, \end{gather}$$
$$\begin{gather}{\rm i} n \hat{u}_n + \hat{v}_n + r \partial_{r} \hat{v}_n = 0 \quad \text{for } r\in(r_0,1/2) . \end{gather}$$Given the constraints (C10) and (2.3), the lowest-order truncation of (C5) and (C6) possible is
\begin{gather} u(r,t) = \hat{u}_0(r,t),\quad v(r,t) = 0, \end{gather}
\begin{gather} T(r,\theta,t) = a_0(r,t) + a_1(r,t) \cos \theta + b_1(r,t) \sin \theta , \end{gather}
where 
$a_0, a_1, b_1, \hat {u}_0$ satisfy the BCs (C8) to (C10).
Inserting (C12) and (C13) into (C2), multiplying by 
$r^2$ and projecting onto respective Fourier mode gives
$$\begin{gather} \dot{a}_0 = r^{{-}1} \partial_{r} ( r \partial_{r} a_0 ) , \end{gather}$$
$$\begin{gather}r^2 \dot{a}_1 ={-} r \hat{u}_0 b_1 - a_1 + r \partial_{r} ( r \partial_{r} a_1 ) , \end{gather}$$
$$\begin{gather}r^2 \dot{b}_1 ={+} r \hat{u}_0 a_1 - b_1 + r \partial_{r} ( r \partial_{r} b_1) . \end{gather}$$
Boundary conditions (C8)–(C9) then imply 
$\lim _{t\to \infty } a_0(r,t) = 1/2$ irrespective of the initial condition. We therefore set 
$a_0 = 1/2$, since variations from this value simply reflect the transient dynamics that is decoupled from the rest of the system.
From (2.4) to (2.6), 
$L,\ X,\ Y$ can be evaluated as
$$\begin{gather} L(t) = \frac{2{\rm \pi}}{\text{A}_0}\int_{r_0}^{1/2} r^2 \hat{u}_0(r,t) \, {\rm d}r , \end{gather}$$
$$\begin{gather}X(t) ={-}\frac{\rm \pi}{\text{A}_0} \int_{r_0}^{1/2} r^2 a_1(r,t) \, {\rm d}r , \end{gather}$$
$$\begin{gather}Y(t) ={-}\frac{\rm \pi}{\text{A}_0} \int_{r_0}^{1/2} r^2 b_1(r,t) \, {\rm d}r . \end{gather}$$Differentiating equations (C18) and (C19) with respect to time and inserting (C15) and (C16) gives
$$\begin{gather} \dot{X} = \frac{\rm \pi}{\text{A}_0} \int_{r_0}^{1/2} r \hat{u}_0(r,t) b_1(r,t) \, {\rm d}r - \left.\frac{\rm \pi}{\text{A}_0} \left(r^2 \frac{\partial a_1}{ \partial r } - r a_1 \right) \right|_{r_0}^{r_1}, \end{gather}$$
$$\begin{gather}\dot{Y} ={-} \frac{\rm \pi}{\text{A}_0} \int_{r_0}^{1/2} r \hat{u}_0(r,t) a_1(r,t) \, {\rm d}r - \left.\frac{\rm \pi}{\text{A}_0} \left(r^2 \frac{\partial b_1}{ \partial r } - r b_1 \right) \right|_{r_0}^{r_1}. \end{gather}$$ Given that the conductive-state solution (A4) takes the form of a Laurent polynomial, we consider a Laurent expansion of the variables 
$\hat {u}_0(r,t)$, 
$a_1(r,t)$, 
$b_1(r,t)$. Truncating each series to the lowest order capable of satisfying BCs (C8)–(C10) gives
$$\begin{gather} \hat{u}_0(r,t) = C(t) (r-r_0) (1-2r) r^{{-}1}, \end{gather}$$
$$\begin{gather}a_1(r,t) = \tfrac{1}{2} A(t) (2r-1) (1 - 2 r_0^2 r^{{-}1} ), \end{gather}$$
$$\begin{gather}b_1(r,t) ={-}\tfrac{1}{2} + \tfrac{1}{2} B(t) (2r-1) (1 - 2 r_0^2 r^{{-}1} ). \end{gather}$$
where 
$(A,B,C)$ are time-dependent coefficients. We note that setting 
$A(t) = C(t) = 0$, 
$B(t) = -(4r_0^2+1)^{-1}$ recovers the conductive-state solution, (A4), exactly.
Inserting (C22) to (C24) into (C17), (C18) and (C19) gives linear relationships between 
$(L,\, X,\, Y)$ and 
$(A,\, B,\, C)$:
$$\begin{gather} L(t) = \frac{(1-2 r_0)^2}{12} C(t) , \end{gather}$$
$$\begin{gather}X(t) = \frac{(1-2r_0)^2 (1 + 6r_0 + 16r_0^2)}{48(1+2 r_0)} A(t), \end{gather}$$
$$\begin{gather}Y(t) = \frac{1 + 2r_0 + 4 r_0^2}{12(1+2r_0)} + \frac{(1-2r_0)^2 (1 + 6r_0 + 16r_0^2)}{48(1+2 r_0)} B(t) . \end{gather}$$ Evaluating the right-hand sides of (C4), (C20) and (C21) using (C22) to (C27) gives the dynamical system (3.1) to (3.3) that is the main focus of this article. The coefficients 
$\alpha,\ \beta,\ y_0,\ y_1$, and 
$k$ can be found by analytical integration (accelerated by symbolic programming). Each of these coefficients depends on 
$r_0$ only as given by
\begin{gather} \alpha = \frac{48}{(1-2r_0)^2}, \quad \beta = \frac{48 (1+4r_0^2)}{(1-2r_0)^2 (1+6r_0+16r_0^2)}, \quad y_0 = \frac{1+12r_0^2}{16 (1+4r_0^2)}, \end{gather}
$$\begin{gather} k = 24 \frac{(1-2r_0) (1 -6r_0 -4r_0^2 -88r_0^3 + 32r_0^4) -96 r_0^3 \ln{(2r_0)}} {(1-2r_0)^5 (1+6r_0+16r_0^2)}, \end{gather}$$
$$\begin{gather}y_1 = \frac{(1-4r_0^2) (1 -8r_0 -224r_0^3 -80r_0^4) - 192r_0^3(1+2r_0+4r_0^2)\ln{(2r_0)}}{24 (1-4r_0^2) (1 -6r_0 -4r_0^2 -88r_0^3 + 32r_0^4) -2304 (1+2r_0)r_0^3 \ln{(2r_0)}}. \end{gather}$$Appendix D. Stability analysis of the ODE system
The fixed points of (3.1) to (3.3) are obtained by setting the right-hand sides equal to zero. There can be up to three fixed points, given by the following.
(i) The conductive state
(D1a–c)\begin{equation} L =0, \quad X = 0, \quad Y = y_0. \end{equation}(ii) The circulating states
(D2a–c)where\begin{equation} L ={\pm} L_1, \quad X ={\mp}\frac{\alpha}{{Ra}} L_1, \quad Y = y_1 + \frac{\alpha\beta}{k {Ra}}, \end{equation}(D3)\begin{equation} L_1 ={\pm} \frac{\beta}{k} \sqrt{ \frac{k \,{Ra}}{\alpha \beta} \Delta y -1 }. \end{equation}
We note that the circulating states only exist if 
${Ra} \geq \alpha \beta / (k \Delta y) = {Ra}^*_1$.
The Jacobian of (3.1) to (3.3) is
\begin{equation} J(L,X,Y) = \left[ \begin{array}{@{}ccc@{}} -\alpha\,{Pr} & -{Ra}\,{Pr} & 0 \\ -k(Y-y_1) & -\beta & -kL \\ kX & kL & -\beta \end{array} \right]. \end{equation}Evaluating the Jacobian at fixed point (D1), one can show that all eigenvalues are negative if 
${Ra} < {Ra}^*_1$, thereby confirming the conductive state is stable when 
${Ra} < {Ra}^*_1$. Above 
${Ra}^*_1$, (D1) becomes unstable and the two circulating states given by (D2) appear, indicating a pitchfork bifurcation. The circulating states are stable provided that 
${Ra}^*_1 < {Ra} < {Ra}^*_2$, with 
${Ra}^*_2$ defined in (4.1). Above 
${Ra}^*_2$, the real part of the complex eigenvalues become positive, rendering the circulating states unstable.
Appendix E. High-${Ra}$ LSC reversal frequency
At very high 
${Ra}$, order reemerges in the system and the large-scale dynamics becomes nearly periodic. The dominant frequency 
$f^*$ of the LSC reversals can be obtained by the pendulum equivalence of (3.1) to (3.3). In this correspondence, the gravitational constant is 
$g_{{eff}} = k l^2 \, {Ra} \, {Pr}$, which gives a small-amplitude frequency of 
$\sqrt {k l \, {Ra} \, {Pr}}/(2{\rm \pi} )$. The oscillation amplitude, however, is not small, which implies that 
$f^*$ depends on both the pendulum length 
$l$ and the maximum swing angle 
$\phi _{{max}}$, both of which can be estimated through an energy balance.
Multiplying (3.2) by 
$X$, (3.3) by 
$Y$, and adding gives the exact relation
\begin{equation} \frac{{\rm d} }{{\rm d} t}{l^2} ={-}2 \beta l^2 + 2\beta \Delta y (Y-y_1). \end{equation}
Assuming periodicity implies that the time average of 
${\rm d}l^2/{\rm d}t$ vanishes, giving
\begin{equation} \langle l^2 \rangle = \Delta y \langle Y - y_1 \rangle , \end{equation}
where 
$\langle {\cdot } \rangle$ indicates a time average.
Consider the effective energy of the pendulum system
\begin{equation} E_{{eff}} = \tfrac{1}{2} k L^2 + {Ra} \,{Pr} (Y-y_1) . \end{equation}Here, the terms on the right-hand side represent the kinetic and potential energy of the mechanical pendulum, respectively. Taking a time derivative and using (3.1) and (3.3), gives the energy law
\begin{equation} \dot{E}_{eff} ={-}\alpha \,{Pr} \, k L^2 + \beta \,{Ra} \,{Pr} (y_0-Y) . \end{equation}
The first term above represents energy dissipation due to damping while the second term represents positive energy injected into the system from the external driving. The assumption of periodicity implies 
$\langle \dot {E}_{{eff}} \rangle = 0$, which gives
\begin{equation} k \alpha \langle L^2 \rangle = {Ra}\, \beta \langle y_0 - Y \rangle . \end{equation}Meanwhile, directly averaging (E3) gives
\begin{equation} \langle E_{{eff}} \rangle = \tfrac{1}{2} k \langle L^2 \rangle + {Ra} \,{Pr} \langle Y-y_1 \rangle . \end{equation} At the bottom of the swing, 
$Y = y_1 - l$, 
$L = L_{max}$, the energy is
\begin{equation} E_{{bot}} = \tfrac{1}{2} k L_{{max}}^2 - {Ra} \,{Pr} \, l . \end{equation}
Due to periodicity, 
$\langle L^2 \rangle = L_{{max}}^2/m$, where 
$m$ is a constant. Although energy is not strictly conserved, it is conserved on average for periodic dynamics. We therefore make the assumption of nearly constant energy, 
$E_{{bot}} = \langle E \rangle$ in order to solve for 
$l$ and 
$\phi _{{max}}$. Setting (E6) equal to (E7) and using (E2) and (E5) gives
\begin{equation} l = \Delta y \left( \frac{(m-1) \beta} {(m-1)\beta + 2 \alpha {Pr}} \right) . \end{equation}
At the apex, 
$\phi = \phi _{{max}}$ and 
$\dot {X}=\dot {Y}=0$. Equations (3.2) and (3.3) then simplify to 
$X^2 = (y_0-Y)(Y-y_1)$. As 
$X = l \sin \phi$ and 
$Y = y_1 - l \cos \phi$, we have 
$l = -\Delta y \cos \phi _{{max}}$, thus providing the value of 
$\phi _{{max}}$ once 
$l$ is given by (E8). Based on figure 1(e), we choose the value 
$m=2.5$ as the midpoint of a sinusoidal (
$m=2$) and a triangular (
$m=3$) waveform.
With the values of 
$l$ and 
$\phi _{{max}}$ determined, the frequency of (large-amplitude) pendulum oscillation, and therefore the frequency of regular LSC reversals, is given by (5.1) in the text. We note that the parameters of figure 4, 
$Pr = 4$ and 
$r_0 = 0.4$, give values 
$l = 0.005$ and 
$\phi _{{max}} = 1.62$.
Converting (5.1) to a dimensional frequency gives
\begin{equation} F^* =\frac{\kappa}{h^2}f^* = c N^*, \end{equation}
where 
$N^* = \sqrt {\beta _T \Delta Tg/h}$ is the Brunt–Väisälä frequency and
\begin{equation} c = \frac{\sqrt{kl}}{ 4 K(\sin^2 ({\phi_{{max}}}/{2})) } \end{equation}
is a constant that depends only on 
$r_0$ and 
${Pr}$. In particular, 
$c$ is independent of the Rayleigh number.
Open access
