Time averages and periodic attractors at high Rayleigh number for Lorenz-like models

TL;DR

This paper analyzes periodic solutions and time averages in Lorenz-like models at high Rayleigh numbers, proving existence and stability of symmetric periodic orbits, and exploring transport properties and hysteresis phenomena.

Contribution

It provides a rigorous proof of symmetric periodic orbit existence and stability in Lorenz '63 equations at high Rayleigh numbers, expanding time averages using elliptic integrals.

Findings

Existence and stability of symmetric periodic orbits confirmed.

Transport can be arbitrarily small on periodic attractors.

Hysteresis loop observed between periodic attractors and equilibria.

Abstract

Revisiting the Lorenz '63 equations in the regime of large of Rayleigh number, we study the occurrence of periodic solutions and quantify corresponding time averages of selected quantities. Perturbing from the integrable limit of infinite $\rho$, we provide a full proof of existence and stability of symmetric periodic orbits, which confirms previous partial results. Based on this, we expand time averages in terms of elliptic integrals with focus on the much studied average `transport', which is the mode reduced excess heat transport of the convection problem that gave rise to the Lorenz equations. We find a hysteresis loop between the periodic attractors and the non-zero equilibria of the Lorenz equations. These have been proven to maximize transport, and we show that the transport takes arbitrarily small values in the family of periodic attractors. In particular, when the non-zero equilibria are unstable, we quantify the difference between maximal and typically realized values of transport. We illustrate these results by numerical simulations and show how they transfer to various extended Lorenz models.

Equations

$$\begin{array}[]{l}X^{\prime}=\sigma(Y-X)\\ Y^{\prime}=\rho X-Y-XZ\\ Z^{\prime}=-\beta Z+XY,\end{array}$$

$$H(\rho ,\beta ,\sigma ,{\text{X}}_0)=\underset{t\to \infty }{\mathrm{limsup}}\frac{1}{t}{\int }_0^{t}X(t)Y(t)dt,$$

$$H(\rho,\beta,\sigma,\textbf{X}_{0})=\left\{\begin{array}[]{cl}0&\text{if }\textbf{X}_{0}\in W^{s}(\textbf{0}),\\ H_{\pm}(\rho,\beta):=\beta(\rho-1)&\text{otherwise. }\end{array}\right.$$

$$\lambda =\frac{\sigma +1}{\beta +2}.$$

$$H(\rho ,\beta ,\sigma ,{\text{X}}_0)\le H_{\pm }(\rho ,\beta ).$$

$$\epsilon ={\rho }^{-\frac{1}{2}},X={\epsilon }^{-1}\xi ,Y={\epsilon }^{-2}{\sigma }^{-1}\eta ,Z={\epsilon }^{-2}({\sigma }^{-1}\zeta +1),t=\epsilon \tau ,$$

$$\begin{array}[]{l}\dot{\xi}=\eta-\varepsilon\sigma\xi\\ \dot{\eta}=-\xi\zeta-\varepsilon\eta\\ \dot{\zeta}=\xi\eta-\varepsilon\beta(\zeta+\sigma).\end{array}$$

$$\begin{array}[]{l}\dot{\xi}=\eta\\ \dot{\eta}=-\xi\zeta\\ \dot{\zeta}=\xi\eta,\end{array}$$

$$A=\frac{1}{2}{\xi }^2-\zeta ,B=({\eta }^2+{\zeta }^2{)}^{1/2},$$

$$K(k)={\int }_0^1\frac{1}{\sqrt{(1-t^2)(1-k^2{t}^2)}}dt\text{ , }E(k)={\int }_0^1\sqrt{\frac{1-k^2{t}^2}{1-t^2}}dt.$$

$$u={\int }_0^{\varphi }\frac{d\theta }{\sqrt{1-k^2{\sin }^2\theta }}\iff \mathsf{am}(u,k)=\varphi$$

$$\operatorname{sn}(u,k)=\sin\big{[}\mathsf{am}(u,k)\big{]}\hskip 14.22636pt\text{ , }\hskip 14.22636pt\operatorname{cn}(u,k)=\cos\big{[}\mathsf{am}(u,k)\big{]}\hskip 14.22636pt\text{ , }\hskip 14.22636pt\operatorname{dn}(u,k)=\Big{(}1-k^{2}\sin^{2}\big{[}\mathsf{am}(u,k)\big{]}\Big{)}^{1/2}.$$

$$\begin{array}{rl}u& =\sqrt{B}\tau \text{ , }{k}_1^2=\frac{A+B}{2B}\\ {\xi }_1(\tau )& =2k_1\sqrt{B}\mathrm{cn}(u,k_1)\\ {\eta }_1(\tau )& =-2k_{1}B\mathrm{dn}(u,k_1)\mathrm{sn}(u,k_1)\\ {\zeta }_1(\tau )& =B(1-2k_1^2{\mathrm{sn}}^2(u,k_1)).\end{array}$$

$$\begin{array}{rl}u& =\sqrt{B}{k}_2^{-1}\tau \text{ , }{k}_2^2=2{\left(1+\frac{A}{B}\right)}^{-1}\\ {\xi }_{2,\pm }(\tau )& =\pm 2\sqrt{B}{k}_2^{-1}\mathrm{dn}(u,k_2)\\ {\eta }_{2,\pm }(\tau )& =\mp 2B\mathrm{sn}(u,k_2)\mathrm{cn}(u,k_2)\\ {\zeta }_2(\tau )& =B(1-2{\mathrm{sn}}^2(u,k_2)).\end{array}$$

$$\begin{split}\displaystyle u&=\sqrt{B}\tau\\ \xi_{3,\pm}(\tau)&=\pm 2\sqrt{B}\operatorname{sech}(u)\\ \displaystyle\eta_{3,\pm}(\tau)&=\mp 2B\tanh(u)\operatorname{sech}(u)\\ \displaystyle\zeta_{3}(\tau)&=B\big{(}1-2\tanh^{2}(u)\big{)}\end{split}$$

$$\zeta =B\cos \phi ,\eta =B\sin \phi ,\xi =\xi$$

$$\begin{array}[]{l}\dot{\xi}=B\sin\varphi-\varepsilon\sigma\xi\\ \displaystyle\dot{\varphi}=-\xi+\varepsilon\frac{\sin\varphi}{B}((\beta-1)B\cos\varphi+\beta\sigma)\\ \dot{B}=-\varepsilon(B+(\beta-1)B\cos^{2}\varphi+\beta\sigma\cos\varphi).\end{array}$$

$$\begin{array}[]{l}\dot{\xi}=f_{1}(\xi,\varphi,B)+\varepsilon g_{1}(\xi,\varphi,B)\\ \dot{\varphi}=f_{2}(\xi,\varphi,B)+\varepsilon g_{2}(\xi,\varphi,B)\\ \dot{B}=\varepsilon g_{3}(\xi,\varphi,B).\end{array}$$

$$A(\xi ,\phi ,B)=\frac{1}{2}{\xi }^2-B\cos \phi ,f_1=\frac{\partial A}{\partial \phi },f_2=-\frac{\partial A}{\partial \xi }.$$

$$I(A,B)={\oint }_{L_i^{A,B}}\phi d\xi =-{\oint }_{L_i^{A,B}}\xi d\phi =\underset{0}{\overset{T_i}{\int }}{\xi }^{2}d\tau ,$$

$$\frac{1}{\Omega(A,B)}=\frac{\partial I}{\partial A}\Bigr{|}_{B}.$$

$$\begin{array}[]{l}\dot{I}=\varepsilon F(I,\theta,B)\\ \dot{\theta}=\Omega+\varepsilon G(I,\theta,B)\\ \dot{B}=\varepsilon\tilde{g}_{3}(I,\theta,B),\end{array}$$

$$\begin{array}[]{l}\displaystyle F(I,\theta,B)=\frac{\partial I}{\partial\xi}\tilde{g}_{1}+\frac{\partial I}{\partial\varphi}\tilde{g}_{2}+\frac{\partial I}{\partial B}\tilde{g}_{3}\\[8.61108pt] \displaystyle G(I,\theta,B)=\frac{\partial\theta}{\partial\xi}\tilde{g}_{1}+\frac{\partial\theta}{\partial\varphi}\tilde{g}_{2}+\frac{\partial\theta}{\partial B}\tilde{g}_{3}\\[8.61108pt] \tilde{g}_{i}(I,\theta,B)=g_{i}(\xi(I,\theta,B),\varphi(I,\theta,B),B),\;\;i=1,2,3.\end{array}$$

$$(I_1,B_1)=(I_0,B_0)+\epsilon (M_1,M_3)+O({\epsilon }^2),$$

$$DM(I_0,B_0)=\frac{\partial (M_1,M_3)}{\partial (I,B)}(I_0,B_0).$$

$$\begin{array}[]{l}\displaystyle M_{1}=\frac{1}{\Omega}\tilde{M}_{1}+\frac{\partial I}{\partial B}\Bigr{|}_{A}M_{3},\;\;\;\tilde{M}_{1}=\int\limits_{0}^{T_{i}}\left[f_{1}\tilde{g}_{2}-f_{2}\tilde{g}_{1}+\frac{\partial A}{\partial B}\Bigr{|}_{\xi,\varphi}\tilde{g}_{3}\right](L_{i}^{A,B}(\tau))d\tau\\ \displaystyle M_{3}=\int\limits_{0}^{T_{i}}\tilde{g}_{3}(L_{i}^{A,B}(\tau))d\tau;\end{array}$$

$$\nu_{\pm}=1+\frac{\varepsilon}{2}\big{(}\operatorname{tr}DM\pm\sqrt{(\operatorname{tr}DM)^{2}-4\det DM}\big{)}+O(\varepsilon^{2}).$$

$$\begin{split}\operatorname{tr}DM&=\left(\frac{\partial I}{\partial k_{i}}\Bigr{|}_{B}\right)^{-1}\Big{[}\frac{1}{\Omega}\frac{\partial\tilde{M}_{1}}{\partial k_{i}}\Bigr{|}_{B}+\frac{\partial I}{\partial B}\Bigr{|}_{A}\frac{\partial M_{3}}{\partial k_{i}}\Bigr{|}_{B}\Big{]}+\frac{\partial M_{3}}{\partial B}\Bigr{|}_{k_{i}}-\frac{\partial M_{3}}{\partial k_{i}}\Bigr{|}_{B}\frac{\partial I}{\partial B}\Bigr{|}_{k_{i}}\left(\frac{\partial I}{\partial k_{i}}\Bigr{|}_{B}\right)^{-1}\\ \det DM&=\frac{1}{\Omega}\left(\frac{\partial I}{\partial k_{i}}\Bigr{|}_{B}\right)^{-1}\Big{[}\frac{\partial\tilde{M}_{1}}{\partial k_{i}}\Bigr{|}_{B}\frac{\partial M_{3}}{\partial B}\Bigr{|}_{k_{i}}-\frac{\partial\tilde{M}_{1}}{\partial B}\Bigr{|}_{k_{i}}\frac{\partial M_{3}}{\partial k_{i}}\Bigr{|}_{B}\Big{]}.\end{split}$$

$$A(\xi ,\theta ,\rho )=\frac{1}{2}{\xi }^2-B\cos \theta ,$$

$$\begin{array}{rl}\xi & =2k_1\sqrt{B}\mathsf{cn}(\sqrt{B}\tau )\\ \sin \phi & =-2k_1\mathsf{dn}(\sqrt{B}\tau )\mathsf{sn}(\sqrt{B}\tau )\\ \cos \phi & =1-2k_1^2{\mathsf{sn}}^2(\sqrt{B}\tau ).\end{array}$$

Full text

Acknowledgments: This paper is a contribution to the project M7 of the Collaborative Research Centre TRR 181 "Energy Transfers in Atmosphere and Ocean" funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 274762653. IO is also supported by the grant of the Russian Science Foundation 19-11-00280. A large part of this research was conducted while JR, BL and IO were affiliated with the University of Bremen. We thank Camilla Nobili for helpful discussions on scaling bounds and Dmitry Turaev for discussions on homoclinic bifurcations.

1. Introduction

In physical processes, infinite time-averaged quantities are often of more interest than particular solutions. Their dependence on parameters is of fundamental theoretical interest and also practical importance. The derivation of bounds for such averages by algebraic optimization has received increasing attention in recent years [6, 9, 15, 14]. Bounds derived by estimates need not be sharp, but also sharp bounds may be misleading, if they are realized by dynamically unstable states [28, 29]. Inspired by [22, 9, 15], in this paper we study such a situation in the regime of large Rayleigh number $\rho\gg 1$ for the famous Lorenz equations and variants thereof. The Lorenz equations, given as follows,

[TABLE]

arose first in the context of atmospheric convection. In the original derivation [11], Lorenz considered a fluid in a periodic box, heated from below and cooled from above, and obtained (1) from a PDE model by retaining only the lowest order Fourier modes. The parameters in (1) stem from the PDE, where $\sigma>0$ is the Prandtl number, characterizing the viscosity of the fluid, and where $\beta>0$ is a shape parameter, measuring the ratio of the length of the box to its height. Often of primary interest is the parameter $\rho\geq 0$ which is the rescaled Rayleigh number, measuring the intensity of the heating. The rescaling is chosen such that a bifurcation occurs at $\rho=1$, and indeed the Lorenz model accurately captures the onset of steady atmospheric convection rolls for $\rho>1$. For larger values of $\rho$, the Lorenz equations do not accurately capture the full PDE model, but the relation between the hierarchy of higher order mode truncations with the convection PDE model continue to be of interest, e.g. [7, 18, 15, 14].

Although physically unrealistic for large regions of parameter space, (1) is frequently used as a benchmark and testbed for nonlinear dynamics, in particular in the famous chaotic regime, but also in the context of time averages [15]. Of particular interest is the average

[TABLE]

which we refer to as transport, following [22]. The quantity $H$ is the mode truncated form of excess heat transport, which defines the Nusselt number by a scaling factor and constant shift. $H$ is well-defined due to the dissipation at infinity in (1) and in particular depends on the initial condition $\textbf{X}_{0}$ of the solution $(X,Y,Z)(t)$ of (1). The dependence of the Nusselt number on $\rho$ is of major physical interest, but is difficult to determine or bound analytically and numerically [28, 29]. Examining $H$ and its parameter dependence provides a tractable, non-trivial case study which can provide insight into the analysis of the Nusselt number for the full PDE. As such, an optimal bound for $H$ has been a longstanding question, which was settled by Souza and Doering [22]. They proved that transport is maximal in the non-trivial equilibria of (1), which emerge in the bifurcation at $\rho=1$. However, since these equilibria are unstable in large parameter regimes, the question remains what values the transport takes in attractors.

The inclusion of stability in the study of transport and the resulting scaling for $\rho\gg 1$ is our main motivation for this paper. The relation to stability is particularly clear in the case $0\leq\rho\leq 1$, where for any value of $\sigma,\beta$ the system admits a Lyapunov function and the origin is the global attractor (see for instance [23]). Hence, for any initial condition we obtain $H=0$. For higher Rayleigh numbers, the functional form of transport becomes more complicated. At $\rho=1$, a pitchfork bifurcation occurs, where the aforementioned non-zero fixed points $\textbf{X}_{\pm}=(\pm\sqrt{\beta(\rho-1)},\pm\sqrt{\beta(\rho-1)},\rho-1)$ emerge. For $\rho$ sufficiently near $1$ these non-zero fixed points seem to attract every trajectory except for $\textbf{X}_{0}$ belonging to the stable manifold of the origin, $W^{s}(\textbf{0})$, so that

[TABLE]

Here $W^{s}(\textbf{0})$ is a surface of dimension 2, which means almost all initial conditions give a positive transport, namely that of the fixed points $H_{\pm}(\rho,\beta)$. This is certainly the case for initial data in their non-trivial basins of attraction.

Upon increasing $\rho$ further, additional periodic orbits emerge and further complicate the function $H$. It has been noticed in [23] that a decisive parameter for (1) at higher $\rho$ is

[TABLE]

For $0<\lambda\leq 1$, the fixed points $\textbf{X}_{\pm}$ are locally stable for all $\rho$, cf. [23] so that the fixed point transport $H_{\pm}$ is observed at least for $\textbf{X}_{0}$ belonging to a basin of attraction of positive measure. On the other hand for $\lambda>1$, the fixed points are locally stable only for $1<\rho<\rho^{*}=\sigma(\sigma+\beta+3)/(\sigma-\beta-1)$. As $\rho$ is increased through $\rho^{*}$, the fixed points $\textbf{X}_{\pm}$ lose stability via a sub-critical Andronov-Hopf bifurcation and at least for some open set containing $\beta=8/3$, $\sigma=10$, generic initial conditions give chaotic solutions [25].

As mentioned, it was proven in [22] that despite this complexity, for all $\rho>1$ and any $\sigma,\beta\in\mathbb{R}_{>0}$, $\textbf{X}_{0}\in\mathbb{R}^{3}$ one has the simple bound

[TABLE]

For the Lorenz equations this bound is actually sharp as it is realized by the steady states $\textbf{X}_{\pm}$. However, since $\textbf{X}_{\pm}$ are unstable for $\rho>\rho^{*},\lambda>1$, the transport that is realized by typical solutions might be much lower. Indeed, numerical experiments presented in [22] indicate that for $\rho>\rho^{*}$ a gap $\Delta H=H_{\pm}-H>0$ occurs, cf. Fig. 1. To the best of our knowledge, quantitative results for the size of the observed transport gap for large Rayleigh number that we provide in this paper have not yet been established previously.

It has been observed already in [19] that for sufficiently large $\rho$ the chaotic attractor collapses and a periodic attractor occurs, but it appears that the implications for transport and other time averages have not been studied. Moreover, we found the arguments given in [19] and also [10] to be incomplete. Sparrow [23] devotes a chapter to the regime of large $\rho$ and obtains various results based on a formal application of the method of averaging. In this paper we rigorously confirm several of these results and provide a complete proof of the existence and stability of symmetric periodic attractors $\textbf{X}_{\rm sym}=\textbf{X}_{\rm sym}(\rho,\beta,\sigma)$ for sufficiently large $\rho$ and its dependence on $\lambda>2/3$. We also prove the existence and instability of a pair of asymmetric periodic orbits and obtain some results on homoclinic orbits. Our analysis is based on the observation that the well known limit system as $\rho\to\infty$ possesses a Hamiltonian structure that allows one to apply the extended Melnikov theory of Wiggins and Holmes [31].

For the periodic orbits it becomes tractable to compute the transport analytically and we quantify the gap $\Delta H$ to leading order: We show that $H_{\rm sym}(\rho,\beta,\sigma):=H(\rho,\beta,\sigma,\textbf{X}_{\rm sym})$ is a monotone increasing function of $\lambda$, for fixed $\rho\gg 1$, whose range is a subinterval of $(0,H_{\pm})$ that limits to this interval as $\rho\to\infty$. Specifically, $H(\rho,\beta,\sigma,\textbf{X}_{\rm sym})\sim\rho$ for $\rho\gg 1$, just as $H_{\pm}$, but with a $\lambda$-dependent downshift that can bring it arbitrarily close to zero, and we provide the leading order term of the downshift in terms of elliptic integrals. The resulting bifurcation diagram contains a hysteresis loop between $\textbf{X}_{\pm}$ and $\textbf{X}_{\rm sym}$ in terms of $\lambda$. This highlights difficulty to recover from low transport when $\lambda$ grows beyond the ‘tipping point’ at $\lambda=1$. Moreover, for any $\gamma\in(0,1]$ we can choose $\sigma(\rho)$, so that $H_{\rm sym}(\rho,\beta,\sigma(\rho))\sim\rho^{\gamma}$ along the periodic attractors $\textbf{X}_{\rm sym}$.

We employ numerical pathfollowing to corroborate the analytical results for large fixed $\rho$ and find that for large $\lambda$ the symmetric periodic orbits terminiate in a symmetric heteroclinic cycle, akin to cycles found in a different regime in [23]. We also compute the stabillity boundary of the symmetric orbits in the $(\lambda,\rho)$-plane, and find that it extends to values of $\rho$ below $200$. In Fig. 1 this region begins near $\rho=313$. It is well known that beyond this boundary various period doubling bifurcations occur [19].

Finally, we turn to variants and extensions of the Lorenz equations and identify regimes in which our analytical results system remain valid and explicitly illustrate this for the Lorenz-Stenflo system.

2. Periodic orbits at large Rayleigh number

It is well known that (1) possesses a semi-Lyapunov function at infinity, cf. [22], and thus a bounded trapping region exist. However, $\textbf{X}_{\pm}$ grow unboundedly with $\rho$, and as pointed out by Sparrow [23, Chapter 7], see also [19], a suitable scaling of the variables with respect to $\rho$ is given by

[TABLE]

which yields the equivalent system

[TABLE]

Notably, the symmetry $(X,Y,Z)\mapsto(-X,-Y,Z)$ of (1) turns into $(\xi,\eta,\zeta)\mapsto(-\xi,-\eta,\zeta)$ in (4).

2.1. Hamiltonian structure and the extended Melnikov theory

It is well known that the limiting system at $\varepsilon=0$,

[TABLE]

is integrable with the conserved quantities

[TABLE]

and we recap known results, essentially as provided in [23], in preparation of the existence and stability proofs.

The character of a solution is completely determined by its location in the $(A,B)$-half-plane with $B\geq 0$. In particular, $B=0$ consists of a line of equilibria. When $B>0$ there are two domains with distinct behavior, $D_{1}=\{0<|A|<B\}$, $D_{2}=\{A>B\}$, and two boundaries $D_{3}=\{A=B\}$, $D_{4}=\{A=-B\}$, since the region $A<-B$ has no real solutions. $D_{4}$ has the simplest solutions, consisting only of the line of equilibria $\xi=\eta=0$, $\zeta=B$. In other regions, the solutions are given in terms of Jacobi elliptic functions and complete elliptic integrals. Following [2], for a given elliptic modulus $0<k<1$, the complete elliptic integrals of the first and second kind are defined by

[TABLE]

For the Jacobi elliptic functions, one first defines the amplitude function $\mathsf{am}(u,k)$ as an inverse function via

[TABLE]

and the Jacobi elliptic functions are then defined

[TABLE]

Each $(A,B)\in D_{1}$ defines a symmetric periodic orbit $L_{1}^{A,B}=(\xi_{1}(\tau),\eta_{1}(\tau),\zeta_{1}(\tau))$ with period T_{1}=4KB^{\scalebox{0.5}[1.0]{-}\frac{1}{2}}, which can be written in terms of Jacobi elliptic functions as

[TABLE]

Each $(A,B)\in D_{2}$ defines a pair of asymmetric periodic orbits $L_{2,\pm}^{A,B}=(\xi_{2,\pm}(\tau),\eta_{2,\pm}(\tau),\zeta_{2}(\tau))$ with period \displaystyle T_{2}=4K(k_{2})k_{2}B^{\scalebox{0.5}[1.0]{-}\frac{1}{2}} that can be represented in elliptic functions as

[TABLE]

The region $D_{3}:\{A=B\}$ corresponds to the line of saddle equilibria $\xi=\eta=0,\;\;\zeta=-B$, each with a pair of homoclinic orbits $L_{3,\pm}^{B}=(\xi_{3,\pm}(\tau),\eta_{3,\pm}(\tau),\zeta_{3}(\tau))$ contained in $D_{3}$ and given by

[TABLE]

Due to the reflection symmetry of (5), without loss of generality we consider only $L_{2}^{A,B}=L_{2,+}^{A,B}$ and $L_{3}^{B}=L_{3,+}^{B}$.

In the following we suppress the dependence of the elliptic functions and integrals on the elliptic modulus $k=k_{i}$, e.g., writing $\operatorname{sn}{u}$ for $\operatorname{sn}(u,k_{i})$ and $K$ for $K(k_{i})$.

Next we deviate from the approach of Sparrow, who proceeds with a formal use of the method of averaging, and instead follow that of Li and Zhang [10] with corrections. In order to exploit the Hamiltonian structure of (5) that is available when $\varepsilon=0$, we introduce polar coordinates with $B$ from (6) given by

[TABLE]

which transform (4) into

[TABLE]

Notably, for $\varepsilon>0$ the radial variable $B$ is no longer a conserved quantity, and must be included as a dynamical variable.

We will use that (12) has the form

[TABLE]

for smooth functions $f_{i}$, $g_{i}$. In this formulation, at $\varepsilon=0$, the first two equations possess the Hamiltonian structure with $A(\xi,\varphi,B)$ defined in (6) serving as a Hamiltonian:

[TABLE]

As noted above, the solutions at $\varepsilon=0$ are given by families of periodic orbits, saddle equilibria and their homoclinic orbits. Therefore, as already noticed in [10], we can apply the extended Melnikov perturbation theory of Wiggins and Holmes [31, 30] in order to identify which periodic and homoclinic orbits of (12) persist for small $\varepsilon>0$. The main idea of this method is that for $\varepsilon=0$ the phase space is represented as a one-parametric family of two-dimensional manifolds (parametrized by $B$ in our case), and on each manifold the system is Hamiltonian. Then, generically, in the phase space there exist two-parameter families of periodic orbits, parametrized by $A$ (the Hamiltonian) and $B$, and one-parameter families of homoclinic orbits. Upon perturbing by $\varepsilon>0$, this structure is destroyed, and, generically, one can expect that only isolated periodic orbits and isolated homoclinic orbits will exist. The existence, local uniqueness and the topological type of these objects can be established with the help of Melnikov integrals.

For the periodic orbits $L^{A,B}_{i}$, the analysis simplifies when changing canonical Hamiltonian variables for $\varepsilon=0$ from $(\xi,\varphi)$ to action-angle variables $(I,\theta)$. The action variable $I$ is computed as

[TABLE]

and the angle $\theta\in[0,1)$ increases from an [math] to $1$ along the periodic orbits with constant frequency $\Omega(A,B)$ given by

[TABLE]

Here and below we use the notation \displaystyle\frac{\partial R}{\partial P}\Bigr{|}_{Q} from thermodynamics and elsewhere to emphasize that we differentiate the quantity $R$ only with respect to its explicit dependence on $P$, neglecting implicit relationships between $P$ and $Q$.

In these new variables (12) turns into

[TABLE]

where

[TABLE]

For small $\varepsilon>0$, a trajectory starting from the initial point $(I_{0},0,B_{0})$ has monotonically increasing $\theta$ from [math] to $1$ with velocity $\varepsilon$-close to $\Omega$, and slowly evolving $(I,B)$ with velocity $O(\varepsilon)$. Hence, the trajectory stays in an $\varepsilon$-neighborhood of the corresponding unperturbed periodic trajectory, and reaches $\theta=1$ after finite time $T_{\varepsilon}=1/\Omega+O(\varepsilon)$, thus returning to the starting plane $\{\theta=0\}$. This defines a Poincaré map $(I_{0},B_{0})\to(I_{1},B_{1})$,

[TABLE]

whose fixed points correspond to periodic orbits, and where $M_{1},M_{3}$ will be Melnikov integral terms. We denote the linearization matrix at $(I_{0},B_{0})$ as

[TABLE]

Theorem 3.2 of [31] states that for any $(I_{0},B_{0})$ for which $M_{1}=M_{3}=0$, $\det DM(I_{0},B_{0})\neq 0$ there exists $\varepsilon_{0}>0$ and an isolated fixed point of the Poincaré map (19) in an $\varepsilon$-neighbourhood of $(I_{0},B_{0})$ for any $0<\varepsilon<\varepsilon_{0}$. This corresponds to a persistent periodic orbit of (4) and (1) for $0<\varepsilon\ll 1$. Moreover, in case $(M_{1},M_{3})\neq 0$ for $0<\varepsilon\ll 1$ there is no periodic orbit in a neighbourhood of $(I_{0},B_{0})$. The expressions for $M_{1}$ and $M_{3}$ are given in [31], and in our notation read

[TABLE]

in particular $(M_{1},M_{3})=0$ is equivalent to $(\tilde{M}_{1},M_{3})=0$. We emphasize that we apply [31, Theorem 3.2] separately for two topologically different families of periodic orbits, lying respectively in domains $D_{1}$ and $D_{2}$. For every fixed point of the Poincaré map $(I_{0},B_{0})\in D_{i}$ there exists small $\varepsilon_{0}>0$ such that for $0<\varepsilon<\varepsilon_{0}$ the corresponding periodic orbit lies entirely in $D_{i}$. Also, it is clear that upon increasing $\varepsilon$ further, the trajectory can touch the boundary of both domains – $D_{3}$, and cross it. This scenario is confirmed below by numerical experiments in section 4, see Fig. 6.

The stability type of such a fixed point, and thus the periodic orbit, is determined by the eigenvalues $\nu_{\pm}$ of $DM(I_{0},B_{0})$ are given by

[TABLE]

Although $(M_{1},M_{3})$ depend on $(I,B)$, it is more convenient to consider them as functions of $k_{i},B$, where $k_{i}$ are the moduli defined in (8), (9). The following simplified formulas for the trace and determinant can be obtained by using the fact that the Melnikov functions are zero at the appropriate value of $(I_{0},B_{0})$, and by changing variables:

[TABLE]

Remark 2.1**.**

We briefly comment on the previous existence and stability studies in the literature. The authors of [10] also follow the method from [31] with the difference that in formula (11) they define the radius as $(B+\rho)$, where $B=\operatorname{const}$ and $\rho$ being the third dynamical variable in system (12) (not the Rayleigh number, as in our notation). However, in the new coordinates they then write formula (14) as

[TABLE]

which is incorrect as the term $\displaystyle\frac{\partial A}{\partial\rho}$ is missing from the Melnikov integrals (20). This is why the results of [10] differ from [23], [19] and our results.

In [19] the Lorenz system of another form is considered. It can be obtained from system (1) via a coordinate transformation and setting $\beta=1$, thus the parameter space is reduced. Moreover, there is a gap in the argumentation. Namely, for the unperturbed periodic orbit $(x_{0}(t),y_{0}(t),z_{0}(t))$ a perturbation $(x_{1}(t),y_{1}(t),z_{1}(t))=O(\varepsilon)$ is considered. The author solves the system of differential equations for $(x_{1}(t),y_{1}(t),z_{1}(t))$ under an assumption $z_{0}(t)\neq 0$ (see [19, formula (6)]), but on the symmetric orbit function $z_{0}(t)$ obviously vanishes two times. Thus, while the results of the computation seem correct, they are not sufficiently justified.

In the book [23] system (4) is analyzed via formally averaging over the unperturbed periodic orbits. Isolated equilibrium points of the averaged system correspond then to isolated periodic orbits of the original system. Formulas (27) and (33), that determine the existence and uniqueness of periodic orbits in domains $D_{1}$ and $D_{2}$, were also obtained there, however without a rigorous proof of monotonicity. Also, analogues of (29) and (36) were obtained in [23], giving the sign of the trace of the linearization matrix. However, the determinant was not computed, and thus the stability of the symmetric periodic orbit and instability of the pair of non-symmetric periodic orbits could not be determined.

Remark 2.2**.**

The question of persistence also arises for the homoclinic orbits (10) of (12) in region $D_{3}$. From the line of the corresponding saddle equilibria, for $0<\varepsilon\ll 1$ only one saddle equilibrium $(\xi,\varphi,B)=(0,\pi,\sigma)$ remains. Hence, it is possible that its stable and unstable manifolds will form homoclinic orbits, and periodic orbits may appear in bifurcations from the homoclinic orbits at $\varepsilon>0$. In [30] such a phenomenon was studied, however, the results are not valid in the claimed generality and do not apply in our case as discussed in §2.4 below.

2.2. Symmetric periodic orbits ($D_{1}$)

In the regime $D_{1}$ the explicit solutions (8) have the following form in polar coordinates:

[TABLE]

Hence, from (15) and (16) the period and the action are

[TABLE]

Substituting these explicit solutions into (20) using the formulas for $f_{i},g_{i}$ from (13) allows to compute the Melnikov integrals explicitly as

[TABLE]

Equating these expressions to zero, equivalently gives the equations

[TABLE]

for $A,B$ as necessary conditions for the existence of persistent periodic orbits. Indeed, the same equations were obtained by C. Sparrow and P. Swinnerton-Dyer (see [23, Appendix K]) by formally using the method of averaging. In fact, part of the following analysis is an extension of that in [23] regarding (24).

We are now ready to formulate and prove our main result for symmetric periodic orbits.

Proposition 2.3**.**

For every $\lambda=\frac{\sigma+1}{\beta+2}>2/3$ and $\varepsilon$ sufficiently small there exists a (locally) asymptotically stable symmetric periodic orbit. It is the unique periodic orbit in an order $\varepsilon$-neighborhood of an explicit solution $L_{1}^{A,B}$, where $(A,B)$ is the unique solution of (24) for the chosen value of $\lambda$.

Proof.

We start by showing that (24) possesses a unique solution in the terms of $k_{1}$ and $B$. In order to solve (24), we first show that the coefficient of $B$ in the second equation is strictly positive for $0\leq k_{1}\leq 1$. Let $e_{1},e_{2}$ be defined by

[TABLE]

For $e_{1}$ write

[TABLE]

and since it is clear from the definitions (7) that $K\geq E\geq 0$, this quantity is strictly positive except when $k_{1}=0$, in which case it is zero. The function $e_{2}$ is also strictly positive, although it is more work to prove it. Rather than interrupt the Melnikov analysis, we include this proof in Appendix A. Thus, in order to have positive solutions for $B$, one should have $2E-K<0$. As noticed in [23, Appendix K], this is possible precisely when $k_{*}<k_{1}\leq 1$ with $k_{*}\approx 0.908909$.

Since the coefficient of $B$ is strictly positive, one can solve the second equation for $B$:

[TABLE]

where $d_{1}=4e_{1}+\beta e_{2}$ denotes the denominator of the first fraction. Substituting this into the first equation gives

[TABLE]

as the remaining equation. Moving the term involving $\sigma$ to the right hand side, using the definition $\lambda=(\sigma+1)/(\beta+2)$, and after some arithmetic one obtains

[TABLE]

Next, we prove the claim of [23] that the right-hand side in equation (27) is a monotonically decreasing function of $k_{1}$ in the domain $k_{*}<k_{1}\leq 1$ and the right hand side limits to $1/3$ at $k_{1}=1$. This immediately implies the existence and uniqueness of the solution of equations (24) for $\lambda>2/3$. Writing the right hand side of (27) as

[TABLE]

we claim that the first factor and both summands in the parentheses are non-negative and monotonically decreasing. Indeed, we compute for the prefactor and the first summand that

[TABLE]

cf. [2]. The second summand is decreasing due the monotonicity of $E$ and $K$, and is positive since $K>2E$. Together with $E(1)=0$ and $K\to\infty$ as $k_{1}\to 1$ the right-hand side of equation (24) is monotonically decreasing from $\infty$ to $1/3$ when $k_{1}$ increases from $k_{*}$ to $1$. This range corresponds to the values $\lambda>2/3$ on the left-hand side of (27). Also, every $k_{1}\in(k_{*},1)$ gives a unique positive value for $B$ via formula (26).

Having found the unique solution to the condition $M_{1}=M_{3}=0$ for each $\lambda>2/3$, we now consider the determinant from (22). Computing the derivatives in (22) gives

[TABLE]

where $\tilde{c}_{0},\tilde{c}_{1},\tilde{c}_{2}$ are coefficients depending on $\sigma,\beta,B$ only. Using the identity (26) to eliminate $B$ and the identity $\sigma=\lambda(\beta+2)-1$ to eliminate $\sigma$, one obtains

[TABLE]

for coefficients $\hat{c}_{0},\hat{c}_{1},\hat{c}_{2},\hat{c}_{3},\hat{c}_{4}$ depending only on $\lambda,\beta$. Finally, applying (27) to eliminate $\lambda$ yields

[TABLE]

where

[TABLE]

This equality expresses $F_{1}$ as a sum of strictly positive terms, since for the first two terms we have

[TABLE]

and we prove positivity of the last term in Appendix A.2. This proves positivity of the determinant and thus existence and local uniqueness for $0<\varepsilon\ll 1$ by [31, Theorem 3.2] as mentioned above.

Having proven the existence of a persistent periodic orbit, we can now prove its stability. Since we already have $\det DM>0$ we study the trace from (22). In region $D_{1}$ one has \displaystyle\frac{\partial I}{\partial B}\big{|}_{A}=\frac{4}{\sqrt{B}}(2E-K), so that by computing the derivatives in the formula (22) we find

[TABLE]

where $\hat{c}_{0},\hat{c}_{1}$ are coefficients depending only on $\sigma$, $\beta$, $B$ and $k_{1}$ (but not $E$ or $K$). Using the first equation of (24) and $\lambda>2/3$ (so that $\beta<2\sigma$) we express $E$ as

[TABLE]

and substituting this into the equation for the trace, we find

[TABLE]

Since $\operatorname{tr}DM<0$ and $\det DM>0$, from (21) it follows that for $\varepsilon>0$ sufficiently small the eigenvalues lie inside the unit circle so that the periodic orbit is asymptotically stable. ∎

2.3. Asymmetric periodic orbits ($D_{2}$)

In the regime $D_{2}$ we consider the explicit solutions $L_{2}^{A,B}(\tau)$ given in (9). In polar coordinates these are given as

[TABLE]

with period and action

[TABLE]

Substituting these explicit solutions into (13) and (20) we compute the Melnikov integrals explicitly as

[TABLE]

Equating these expressions to zero, we see that we must solve

[TABLE]

for $A,B$. As above, the same conditions were obtained in [23, Appendix K] by formal averaging. We can now state and prove our existence and instability result for asymmetric periodic orbits:

Proposition 2.4**.**

For every $\varepsilon$ sufficiently small and all $2/3<\lambda<1$ there exists a pair of saddle-type asymmetric periodic orbits. Each is unique in an order $\varepsilon$-neighborhood of $L_{2}^{A,B}$, where $(A,B)$ corresponds to a unique solution $(k_{2},B)$ of (31) for the chosen value of $\lambda$.

Proof.

The coefficient of $B$ in the first equation of (31) is always positive, since

[TABLE]

so that the unique solution in the terms of $B$ is

[TABLE]

where $d_{2}$ denotes the denominator. Substitution into the second equation of (31) gives an equation for $k_{2}$:

[TABLE]

Next, we prove that the right-hand side of (33) monotonically decreases from $1$ to $1/3$ as $k_{2}$ increases from [math] to $1$, which implies precisely for every $\lambda\in(2/3,1)$ there exists a unique solution $k_{2}$. To prove monotonicity, first we compute

[TABLE]

where

[TABLE]

The three summands of $P$ are each negative since for the first we have

[TABLE]

for the second $K\geq 0$, $2E-K<0$, and for the third we estimate its nontrivial factor from above by the negative $-4K^{2}+2E^{2}$ using the quadratic estimate $2KE\leq K^{2}+E^{2}$.

Regarding the determinant, by computing the derivatives in (22) we find it has the form

[TABLE]

where $\tilde{c}_{0},\tilde{c}_{1},\tilde{c}_{2}$ are coefficients depending only on $\sigma,\beta,B,k_{2}$. Again using the identities (32), $\sigma=\lambda(\beta+2)-1$ and (33) to eliminate $B,\sigma,\lambda$ one finds

[TABLE]

where $F_{2}$ is given by

[TABLE]

In Appendix A.3 we prove that $F_{2}(k_{2})>0$ for $0<k_{2}<1$, hence it follows that $\det DM<0$. This proves existence and local uniqueness for $0<\varepsilon\ll 1$ by [31, Theorem 3.2] as mentioned above.

Turning now to the question of stability, in region $D_{2}$ one has \displaystyle\frac{\partial I}{\partial B}\big{|}_{A}=4(2E-K(2-k_{2}^{2}))/(k_{2}\sqrt{B}), and hence by computing the derivatives in (22) one obtains

[TABLE]

where $\check{c}_{1},\check{c}_{0}$ are coefficients depending only on $\sigma$, $\beta$, $B$ and $k_{2}$ (but not $E$ or $K$). Using the identity (32) to eliminate $B$, the identity $\sigma=\lambda(\beta+2)-1$ to eliminate $\sigma$, and the identity (33) to eliminate $\lambda$ we obtain

[TABLE]

Since $\operatorname{tr}DM<0$ and $\det DM<0$, by (21), for any sufficiently small $\varepsilon>0$, there is one eigenvalue inside the unit circle and one outside, thus these periodic orbits are of saddle type. ∎

2.4. Homoclinic orbits ($D_{3}$)

In $D_{3}$, the homoclinic solutions from (10) are given for every $B>0$ in polar coordinates as

[TABLE]

For $\varepsilon>0$ the slow evolution of the quantities $A$ and $B$ can be written as (cf. [23]):

[TABLE]

where zero index denotes the unperturbed solutions from (10).

As already mentioned in Remark 2.2, in the present case the only equilibrium point which survives the perturbation $0<\varepsilon\ll 1$ is $(\xi,\eta,\zeta)=(0,0,-\sigma)=:\Xi_{\sigma}$, corresponding to the origin in the original Lorenz system (1), and we have $B=A=\sigma$ at that point. For $\varepsilon>0$ the equilibrium $\Xi_{\sigma}$ has two stable eigenvalues $\displaystyle\nu_{1}=-\sqrt{\sigma}+O(\varepsilon)$ and $\nu_{2}=-\varepsilon\beta$ and the unstable one $\displaystyle\nu_{3}=\sqrt{\sigma}+O(\varepsilon)$. The leading stable direction corresponding to the eigenvalue $\nu_{2}$ is the invariant line $\{\xi=0,\;\eta=0\}$ of (4), which is $\{A=B\}$ in (37).

We note that the second equation of system (37) determines the dynamics of variable $B$, and is in particular equivalent to the last equation of system (13). Integration over the unperturbed homoclinic orbit gives, to leading order,

[TABLE]

which is nonzero for $\varepsilon>0$ since $\beta>0$. This is inconsistent with the persistence result for perturbed homoclinic orbits in [30], where it is claimed that this integral always vanishes. Hence, [30] is not applicable in the present situation (and is not valid in the claimed generality).

Towards a correct prediction, first note that a homoclinic orbit converges to the equilibrium point as $\tau\to-\infty$ along the unstable direction, and $\lim\limits_{\tau\to-\infty}A(\tau)=\lim\limits_{\tau\to-\infty}B(\tau)=\sigma$. Let us assume that a generic homoclinic orbit exists for $0<\varepsilon\ll 1$, which approaches the equilibrium along the leading direction. In the present case this is the line $\{A=B\}$, along which it slowly converges to the equilibrium with rate $\nu_{2}=-\varepsilon\beta$. For $0<\varepsilon\ll 1$ this slow part dominates the Melnikov integral, which means a connection of the unstable manifold along $\{A=B\}$ requires that the jumps of $B^{2}$ in (38), and of $A^{2}$ coincide to leading order in $\varepsilon$. The latter can be computed as

[TABLE]

which equals (38) if and only if $\sigma=\frac{1}{3}(1+2\beta)$, that is, $\lambda=2/3$. This is consistent with the discussion of periodic orbits in Proposition 2.4 and, although we do not give a full proof, we thus expect there exists a homoclinic bifurcation curve that tends to $\lambda=2/3$ when $\varepsilon\to 0$.

3. Implications for infinite time averages

In this section we use the angle bracket notation to denote the infinite time average:

[TABLE]

In [9], Goluskin illustrates an application of semi-definite programming to dynamical systems by obtaining bounds on time averages for the Lorenz equations. He proves that the time averages of the following monomials are maximized by the value obtained at the fixed points $\textbf{X}_{\pm}$ for a wide range of parameters $(\beta,\sigma)$ and for all $0<\rho<\infty$:

[TABLE]

Due to the form of the Lorenz equations, certain infinite time averages must be proportional. For instance, one has

[TABLE]

In this way one has the following equalities

[TABLE]

and hence it suffices to consider the following three time averages:

[TABLE]

As stated previously, for $\lambda>1$ the fixed points become unstable for sufficiently large $\rho$, hence while these sharp upper bounds are indeed valid, they do not represent the values that most trajectories obtain. However, with the results of the previous section in hand we can provide complementary results which give values for these time averages which are observed for all trajectories within the nontrivial basin of attraction of the symmetric periodic orbit. For such initial conditions the infinite time averages above are given by the average over the periodic orbit, i.e.,

[TABLE]

where $T_{1}^{\varepsilon}$ is the period of the solution $\zeta_{1}(\tau,\varepsilon)$. It seems likely these values do not have a simple closed form and instead we compute the time averages via an expansion in $\varepsilon$. For instance, the lowest order term can be found from the explicit formulas for the unperturbed solution and period, $\zeta_{1}(\tau)$ and $T_{1}$, as follows

[TABLE]

In this way, we obtain the following expressions for the infinite time averages, expressed in terms of $\rho$ rather than $\varepsilon$:

[TABLE]

Recall the expressions $e_{1},e_{2},d_{1}$ defined in (25), (26) were shown to be positive. Hence for the time averages $\langle Z\rangle$ and $\langle Z^{2}\rangle$, these expressions resolve the coefficient of the leading order term as a function of $\sigma,\beta$ which is strictly less than one, hence less than that of the fixed point value. For $\langle Z^{3}\rangle$, the term $d_{1}e_{1}$ is always positive, whereas one can check that the other expression inside the parentheses is positive for all $\lambda>2.5611...$, whereas is is negative for $\lambda$ less than this. Hence by fixing such $\lambda$ and choosing $\beta$ sufficiently large this exceeds the fixed point value. This agrees with Goluskin’s result however, since the region in parameter space where $\langle Z^{3}\rangle$ is maximized at $\textbf{X}_{\pm}$ does not include large $\beta$.

As mentioned previously, the transport $\langle XY\rangle$ is of particular interest, since this is the truncated version of the Nusselt number from fluid dynamics and hence the most well studied such average. Towards understanding the organization of solution branches and the associated transport, we denote the transport of the stable symmetric periodic orbit by

[TABLE]

and we also compute the transport obtained by the unstable, asymmetric periodic orbits

[TABLE]

Next, we cast $H_{1},H_{2}$ in terms of $\rho$ and rescale to a finite range of transport values. This gives

[TABLE]

We first show that $h_{\rho,1}(2/3)=h_{\rho,2}(2/3)=0$ understood as the limit $\lambda\searrow 2/3$, and analogously $h_{\rho,1}(1)=h_{\rho,2}(\infty)=\beta$. Indeed, $R_{1}(2/3)=R_{2}(2/3)=1$, $R_{1}(\infty)=R_{2}(1)=0$ due to the following.

• •

$R_{1}(2/3)=1$: $\lambda\to 2/3$ gives $k_{1},k_{2}\to 1$ so that $E(1)=0$, $K(1)=\infty$ and (26), (32) imply $B/\sigma\to 1$;

• •

$R_{2}(2/3)=1$: the previous also implies $(K(2-k_{2}^{2})-2E)/(Kk_{2}^{2})\to 1$;

• •

$R_{1}(\infty)=0$: $\lambda\to\infty$ means $k_{1}\to k_{*}$, $B$ is bounded and $E/K\to 1/2$;

• •

$R_{2}(1)=0$: $\lambda\to 1$ gives $k_{2}\to 0$ and using (32) as well as $E=K=\pi/2$ at $k_{2}=0$ implies

[TABLE]

The differences to the scaled maximum transport $H_{\pm}/\rho=\beta(1-\rho^{-1})$ are the positive quantities

[TABLE]

In particular, the stable symmetric periodic orbits $L$ yield the same order of magnitude of transport with respect to $\rho$, but feature a $\lambda$ dependent downshift that vanishes at $\lambda=\infty$, i.e. at $\textbf{X}_{\pm}$. In addition, from our perhaps rough error estimates we obtain a correction of order (at least) $\rho^{-1/2}$ compared to the next order being $\rho^{-1}$ for $\textbf{X}_{\pm}$. Numerical computations suggest that this term might in fact be of order $\rho^{-1}$, but we do not explore this further here.

4. Numerical computations and hysteresis loop

We present numerical results that corroborate the analytical results for $\rho=\infty$ and $1\ll\rho<\infty$ of the previous section, and that highlight the occurrence of a hysteresis loop. In Figure 2(a) we plot the numerical evaluation of (41) for symmetric periodic orbits and (42) for asymmetric ones. However, the elliptic integral routines of the current version of the software Mathematica for $B,E,K$ have failed to numerically converge for transport below $\approx 0.6$. The analytical prediction is that the branches terminate at $\lambda=2/3$ in homoclinic bifurcations of the zero equilibrium and thus at zero transport. Indeed, at $\lambda=2/3$ the intersection of the level sets of the conserved quantities $A,B$ forms a symmetric pair of homoclinic loops, cf. Fig. 3(a), which is the limit of the branch of symmetric periodic orbits, and each branch of asymmetric periodic orbits limits on one of the homoclinic loops.

The arrangement of branches in Figure 2(a) together with the stability properties suggest a hysteresis loop of equilibria and periodic orbits in terms of $\lambda$: For $\lambda<2/3$ the equilibria $\textbf{X}_{\pm}$ that maximize transport are stable, while for $\lambda>1$ the symmetric periodic orbit are. Intermediate values $2/3<\lambda<1$ lie in the analytically predicted region of bistability with stable equilibria $\textbf{X}_{\pm}$ and stable symmetric periodic orbit $L$.

For large finite $\rho$ and moderate values of $\lambda$, numerical pathfollowing computations using Auto corroborate that branches of symmetric and asymmetric periodic orbits persist as predicted. See Figure 2(b). Towards zero transport the branches of symmetric and asymmetric periodic orbits appear to terminate in homoclinic bifurcations to the zero equilibrium near $\lambda=0.688$. See also Fig. 3(a,b). The asymmetric periodic orbits are unstable as predicted, but the symmetric periodic orbits lose stability at low transport. For $\rho=1000$ this occurs at $\lambda=\lambda_{\rm bp}(\rho)\approx 0.79$ in a supercritical pitchfork bifurcation. A branch of stable periodic orbits bifurcates, which are asymmetric in a different sense, but these lose stability at $\lambda=\lambda_{\rm pd}(\rho)\approx 0.787$ in a period doubling bifurcation. We plot the loci of $\lambda_{\rm bp},\lambda_{\rm pd}$ in Figure 5(b), showing that as $\rho$ increases, $\lambda_{\rm bp},\lambda_{\rm pd}$ approach $\lambda=2/3$.

Further numerical simulations corroborate the hysteresis-type loop: for $\lambda<2/3$ the maximum transport equilibria $\textbf{X}_{\pm}$ appear to be global attractors, while for $\lambda>1$ this seems to be the stable symmetric periodic orbit, as in Fig. 3(b). See also Fig. 1(b). For $2/3<\lambda<1$ the situation with large finite $\rho$ is complicated by the fact that symmetric periodic orbits are born in a homoclinic bifurcation at some $\lambda_{\rm hom}\in(2/3,1)$ and, as mentioned, are unstable until a bifurcation point $\lambda_{\rm bp}\in(\lambda_{\rm hom},1)$. Up to the aforementioned region of stable asymmetric periodic orbits that bifurcate from $\lambda_{\rm bp}$, the global attractors for $\lambda<\lambda_{\rm bp}$ seem to be $\textbf{X}_{\pm}$ and the region of bistability with the symmetric periodic orbit is effectively $\lambda\in(\lambda_{\rm bp},1)$.

Using time-varying values of $\lambda$ we consistently found hysteresis as plotted in Figure 4(a): for slowly increasing $\lambda$ from [math], the solution is quickly close to $\textbf{X}_{+}$ so that maximum local transport is realized, i.e., transport computed over a time interval of finite length, which can be chosen longer for slower change of $\lambda$. As $\lambda$ increases beyond $1$, the solution eventually approaches the stable symmetric periodic orbit, so that the realized local transport is smaller than the theoretical maximum. Analogous to delayed bifurcations, this transition to the periodic orbit does not occur immediately after crossing $\lambda=1$ at $t=100$, but with a delay, here until around $t=190$. Subsequent decrease of $\lambda$ causes the solution to track the stable branch of symmetric periodic orbits, cf. Figure 4(b), which decreases the observed local transport further until $\lambda=\lambda_{\rm bp}\approx\lambda_{\rm pd}$. Upon decreasing $\lambda$ below this threshold, a switch to a stable equilibrium $\textbf{X}_{\pm}$ occurs, thus re-creating maximum local transport.

While the asymptotically predicted branch of symmetric periodic orbits of Fig. 2 continues for increasing $\lambda$ monotonically and unboundedly, we found that for finite $\rho$ this is not the case. As plotted in Figure 5, the branch of stable symmetric periodic orbits turns around, oscillates, and appears to terminate in a symmetric heteroclinic bifurcation of $\textbf{X}_{\pm}$ at a finite value of $\lambda$. See also Fig. 3(c). Upon increasing $\rho$, this turning and termination occurs at larger values of $\lambda$. Hence, this scenario is consistent with the analytical results, which concern $\rho\to\infty$ for bounded ranges of $\lambda$. The appearance of a symmetric heteroclinic cycle between $\textbf{X}_{\pm}$ in the Lorenz system has already been noticed in [23, 8], albeit apparently not in the regime of large $\rho$.

The transport at such a heteroclinic cycle is that of the symmetric equilibria, i.e. $\beta(1-\rho^{-1})$, which is indeed very closely matched at the numerical termination points. The $(\lambda,h_{\rho})$-loci of the termination points lie near the curve of symmetric periodic orbits (blue solid), which therefore appear to predict the loci of the heteroclinic cycles. The oscillating stability along the branch creates multi-stable regions in $\lambda$; we note that generic unfoldings of the type of heteroclinic cycle with leading oscillating dynamics yield chaotic attractors [1].

We plot the projection of a solution near the symmetric heteroclinic cycle into the $(A,B)$-plane in Figure 6. This corroborates the conjecture by Sparrow in [23] that orbits bifurcate from $\rho=\infty$ which cross through the diagonal $A=B$. We find that also the solutions near the double homoclinic loop with small transport cross the diagonal. In contrast, the solutions for moderate transport remain in $D_{1}$ as predicted by the limit $\rho\to\infty$.

5. Other Lorenz-like systems

The analysis for large $\rho$ carries over to other models related to the Lorenz equations (1). For the general context of extensions, we refer to [4, 23, 18, 14] and the references therein. For illustration purposes, let us consider linear additions to (1) in the form

[TABLE]

with $w\in\mathbb{R}^{k}$, linear $\textbf{A},\textbf{B}$ and constant b. Upon rescaling as in (3) and $w=\varepsilon^{-j}\omega$, with $\mathbf{\Xi}=(\xi,\eta,\zeta)^{\intercal}$ we obtain the form

[TABLE]

where $\mathbf{\Phi}_{\varepsilon}$ is the right hand side in (4). For $\textbf{A}=\textbf{B}=0$, i.e., in absence of $\omega$, the difference to (4) is of order $\varepsilon^{2}$. Hence, the leading order analysis is unchanged, which means that periodic orbits bifurcate/persist as for $\textbf{b}=0$, although their symmetry properties may be broken. In particular, this applies to the Lorenz models with offsets from [27, 17] for which one can also show that the transport is maximized in an equilibrium [16].

Non-zero B generally requires $j\geq 1$ for a regular limit in which the right hand side of the equation for $\omega$ becomes independent of $\omega$, and vanishes for $j>1$. For $j=1$ we obtain, up to terms of order $\varepsilon^{2}$,

[TABLE]

with $\textbf{A}_{1}$ the first row of A. For B of the form $\textbf{B}=[B_{1}|0|0|B_{2}]$ the equation for $\omega$ has the slow form

[TABLE]

which occurs with $w\in\mathbb{R}$ in the Lorenz-Stenflo model from [24], its magnetic variant [26], and with $w\in\mathbb{R}^{2}$ in the models from [12]; for the latter we choose $\alpha=O(\varepsilon)$ and shift the auxiliary variables (which gives $\textbf{b}\neq 0$) to obtain the form (47). An extension of Lorenz-Stenflo with nonlinear additional equations is considered in [13], but still fits into the present framework when, e.g., scaling the variables in addition to Lorenz-Stenflo with $j=3$ and choosing Lewis number of order $\varepsilon^{-2}$. Other extensions of the Lorenz model with two nonlinear auxiliary equations are studied in [3, 21, 7], which also fit into the present framework when suitably scaling the auxiliary modes and parameters. However, in many cases the situation is more complicated, for instance for the three-dimensional extension in [20, 7].

We next show that for the case (47) the results of the previous sections also carry over; the following analysis is more explicit in §5.1 for the model from [24]. In the case (47) the Melnikov analysis of §2 can be simply extended by adding the slow equation for $\omega$ to the action-angle formulation. The additional Melnikov-integral term $\textbf{M}_{3}$ is then simply the integral of $B_{1}\xi(t)+B_{2}\omega_{0}$ over the period $T$, with $\omega_{0}$ constant. Since $\xi$ has zero average ($\dot{\varphi}=\xi$ at $\varepsilon=0$ in (12)), this term becomes $\omega_{0}/T$, with the period $T$, so that $\textbf{M}_{3}=0$ requires $\omega_{0}=0$. This means that the values of the other two Melnikov-integrals for (46), $\tilde{M}_{1},\tilde{M}_{2}$, actually coincide with those of $M_{1},M_{2}$ from §2. The non-degeneracy condition turns into invertibility of the matrix

[TABLE]

where $\tilde{M}_{2}$ is independent of $\omega$, and it turns out that also $\tilde{M}_{1}$ is: In its integrand $F$ from (17), the additional term from $\textbf{A}_{1}\omega_{0}$ is constant and has a factor $\frac{\partial I}{\partial\xi}=\frac{\partial I}{\partial A}\frac{\partial A}{\partial\xi}=T\xi$, where $\xi$ has zero average as noted above. Hence, the matrix has lower left triangular block structure and the block $\partial_{\omega}\textbf{M}_{3}=TB_{2}$ is invertible, if $B_{2}$ is. In that case the non-degeneracy condition is therefore the same as for $M_{1},M_{2}$ from the original Lorenz system.

5.1. Lorenz-Stenflo

The Lorenz-Stenflo system is given as follows:

[TABLE]

This system is a mode truncation of the rotating Boussinesq equations:

[TABLE]

where one is considering convection in a fluid in a rotating frame, and a term representing the Coriolis force has been added. One obtains (48) by making the analogous reduction to a system of ODE’s for the Fourier coefficients, but one must include an additional Fourier coefficient $V(t)$ in the expansion of the velocity, which couples to the $X$-mode via the Coriolis force. The parameter $s$ measures the speed of the rotation.

Since $X$ and $V$ both represent velocity variables, we expect they have the same scaling in $r$, hence we scale

[TABLE]

and we obtain the system of equations

[TABLE]

The limiting system when $\epsilon=0$ coincides with (4) except trivial dynamics in the variable $\chi$, so that the system now admits three invariants of motion

[TABLE]

Using the first two invariants as for the Lorenz system, (50) can be solved at $\varepsilon=0$, where the solutions have a different form depending on the choice of $(A,B)$ as described in §2. Analogous to §2, we change coordinates via

[TABLE]

and (50) becomes

[TABLE]

Hence, $B$ and $\chi$ are constant at $\varepsilon=0$, and for fixed $B$ and $\chi$, the remaining system for $(\xi,\phi)$ possesses the same Hamiltonian structure as (5). Converting to action angle coordinates the system becomes

[TABLE]

where

[TABLE]

In this case we have a three dimensional Melnikov function given by

[TABLE]

In order to find the persistent periodic orbits, we need to find the zeros of this vector valued Melnikov function such that the non-degeneracy condition

[TABLE]

is satisfied. Since we can write

[TABLE]

it suffices to find $A,B$ such that $\tilde{M}_{1}=M_{3}=M_{4}=0$. Explicitly the Melnikov integrals are given via

[TABLE]

Since we aim at illustration, we consider $|A|\leq B$ only, and then compute

[TABLE]

As noticed a priori for such an extension of the Lorenz system, the first two Melnikov functions are the same as for (4), the third vanishes if and only if $\chi_{0}=0$, and $M_{1}$, $M_{2}$ are independent of $\chi_{0}$. Hence,

[TABLE]

and the determinant is given by

[TABLE]

which is non-zero as shown in §2.

6. Discussion

In this paper we have revisited the dynamics of the Lorenz equation in the regime of large Rayleigh number $\rho$, which is known to feature periodic attractors rather than the famous chaotic dynamics [19, 23, 3]. Our main motivation was to study properties of transport of attractors in a parameter regime where states that maximize transport are dynamically unstable. For the Lorenz equations it was proven in [22] that maximal transport is realized by the non-zero fixed points, which are unstable for $\rho>\rho^{*},\lambda>1$. However, we found that the literature concerning existence and stability theory of periodic states for large $\rho$ was incomplete. We have therefore provided a rigorous treatment, which essentially confirms the predictions of [23]. Numerical computations for large finite $\rho$ based on continuation methods and direct simulations have further corroborated these findings. In addition, we have quantified the transport of the periodic attractors and thus the gap of transport compared with the maximum possible. In particular, the transport of the periodic attractors can be arbitrarily small in a parameter range of bistability, where the states that maximize transport are also stable. Indeed, for fixed $\rho$ we have identified a hysteresis loop in terms of the parameter $\lambda=\frac{\sigma+1}{\beta+2}$, which illustrates difficulty to recover from a loss in transport once $\lambda$ exceeds the ‘tipping point’ $\lambda=1$. Moreover, we have computed the stability boundary of periodic attractors in the $(\lambda,\rho)$-plane and found that it extends to relatively low values of $\rho$ below $200$. For fixed $\rho$ we also found a relation to well-known period-doubling bifurcations and symmetric heteroclinic cycles, which produce further regions of bi- and multi-stability of local attractors.

The Lorenz equations are the crudest mode truncation of the physical model, and there are numerous extensions. For several such generalisations, we have found that our results apply in suitable parameter regimes, in particular for the Lorenz-Stenflo system [24]. Although our results have no immediate implications in the context of atmospheric convection, we believe they provide a relevant case study for the relation of theoretical bounds and dynamically realized transport. The approach by perturbing selected solutions from the infinite Rayleigh number limit by exploiting structural properties would be interesting to explore for higher mode truncations and even the viscous Boussinesq equations. Indeed, recent numerical investigations for meaningful bounds in the Boussinesq equation are based on specific solutions and consider stability properties [28, 29]. We remark that the mode reduced Nusselt number $\mathrm{Nu}=1+\frac{2}{\beta\rho}H(\rho,\beta,\sigma,\textbf{X}_{0})$, cf. [22], is bounded by $3$ as $\rho\to\infty$ due to the transport bound from [22]. However, this is far from the ’ultimate’ or ‘classical’ Nusselt number bounds of order $\rho^{1/2}$ or $\rho^{1/3}$ for the PDE model [28, 29].

The present paper makes a step towards completely settling the question of transport for the Lorenz model. The set of parameter values for which the transport has not been analytically determined is now reduced to a compact set for which the dynamics are chaotic. In the large $\rho$ regime, we have analytically determined stable structures and their transport. Although we have found numerical evidence for further stable invariant structures, it numerically appears (but remains to be proven) that for fixed $\lambda>1$ and sufficiently large $\rho$ the symmetric periodic orbits are the only attractors. In the chaotic regime for intermediate Rayleigh numbers the transport is also reduced compared to the non-zero steady states. However, despite the numerous analytical results for the Lorenz attractor, it seems difficult to quantify the transport in that case. It would also be interesting to explore the possible emergence of discrete Lorenz attractors in extended Lorenz systems such as Palmer’s [17], which is close to a periodic forcing of the Lorenz in a suitable parameter regime.

Appendix A Positivity of elliptic integral expressions

A.1. First elliptic integral expression

Here we prove that

[TABLE]

for any $0\leq k\leq 1$. First, note from the explicit formulas

[TABLE]

For $0\leq k^{2}\leq 3/4$, the integrand is pointwise positive a.e., hence the integral is positive. Indeed, for $0\leq k^{2}\leq 1/2$ one has

[TABLE]

whereas for $1/2<k^{2}\leq 3/4$, one has

[TABLE]

On the other hand, for $3/4<k^{2}\leq 1$ the integrand is no longer pointwise positive. Instead, note that

[TABLE]

and since for $k^{2}\in[3/4,1]$ it follows that

[TABLE]

Hence the integrand is a strictly monotonically increasing function of $t$ for $t\in[0,1]$. The minimum of the integrand for $t\in[0,1]$ is thus achieved at $t=0$, with minimum equal to $3-4k^{2}$. On the other hand note that the integrand tends toward positive infinity as $t\to 1$. Thus for any $p\geq 0$ we can define $t_{p}=t_{p}(k)$ to be the point such that the integrand is equal to $p$, ie

[TABLE]

Note that $0\leq t_{p_{1}}<t_{p_{2}}<1$ for all $0\leq p_{1}<p_{2}<\infty$ follows trivially, and furthermore that (52) is equivalent to

[TABLE]

with

[TABLE]

In particular, the integrand is less than zero for $t<t_{0}$ and greater than zero for $t>t_{0}$, where

[TABLE]

which is easily seen to be a monotonically increasing function of $k$ for $k\in[\sqrt{3}/2,1]$. One can therefore split the domain of integration via

[TABLE]

and the desired positivity will then follow from lower bounds on $I_{-}$ and $I_{+}$. We consider two cases:

1. (1)

$1\geq k\geq\hat{k}_{1}$, where $\hat{k}_{1}\approx.9266$: In this case the integrand tends to infinity sufficiently quickly as $t\to 1$ that we can use the most naive bounds on $I_{-}$. As mentioned the integrand achieves its minimum at $t=0$, and since $0<t_{0}(k)\leq t_{0}(1)=1/2$ one has

[TABLE]

On the other hand, letting $t_{2}$ be defined as in (52), note

[TABLE]

Hence we have $I_{+}+I_{-}>0$ as long as $t_{2}\leq 3/4$. But, noting that

[TABLE]

for $k=1$, it follows from monotonicity of the integrand and the continuity of $t_{p}(k)$ with respect to $k$ that $t_{2}<3/4$ for all $1\geq k\geq\hat{k}_{1}$ where $\hat{k}_{1}$ is defined such that

[TABLE] 2. (2)

$\hat{k}_{2}=.927\geq k\geq\sqrt{\frac{3}{4}}$: Note this region is deliberately chose to overlap with the region in the first case. This is easily seen, since for instance one has

[TABLE]

However in this region we have tighter bounds on $I_{-}$. Again due the monotonicity the integrand achieves its minimum at $t=0$, and this time $0<t_{0}(k)\leq t_{0}(\hat{k}_{2})$, hence one has

[TABLE]

Again we have

[TABLE]

hence we have $I_{+}+I_{-}>0$ at least as long as $t_{2}\leq 1-\frac{.185}{2}=.9075$. But it is easily seen that

[TABLE]

for all $.927\geq k\geq\sqrt{\frac{3}{4}}$, hence the result follows.

A.2. Second elliptic integral expression

Here we prove

[TABLE]

for $1\geq k\geq k_{*}$, where $k_{*}\approx.9089$ is defined as the value of $k$ such that $K(k_{*})=2E(k_{*})$. Note that for all $k\leq 1$ one has

[TABLE]

Note that the polynomial

[TABLE]

is positive for all $x>r_{1}\approx 2.34305$. Note that $K/E$ is monotonically increasing, and $K/E=2.34305$ when $k=\hat{k}_{1}\approx 0.949509$. On the other hand, when $k\leq 0.94951$, one has

[TABLE]

But the polynomial

[TABLE]

is non-zero for all $x$, thus proving the bound.

A.3. Third elliptic integral expression

Here we show that

[TABLE]

satisfies $F_{2}(k)>0$ for all $0<k<1$ in two steps, first for $0.998357\approx k_{\mathrm{u}}<k<1$ and next for $0<k<k_{\ell}=0.9984$.

For the interval $k_{u}<k<1$, note that the 2nd and 3rd terms of $F_{2}$ can be rewritten as

[TABLE]

This is strictly positive for all $0<k<1$: First, for $0<k<\sqrt{\frac{2}{3}}$ one has

[TABLE]

Second, using Mathematica, $\frac{K(\sqrt{2/3})}{E(\sqrt{2/3})}\approx 1.6088<8/3$, and since $K/E$ is monotonically increasing one has

[TABLE]

Therefore,

[TABLE]

Since $(1-k^{2})K^{3}$ equals $\frac{\pi^{3}}{8}$ at $k=0$, monotonically decreases for $0<k<1$, and tends to zero as $k\to 1$, there is $k_{u}$ with

[TABLE]

Using Mathematica, we find $k_{u}\approx 0.998357$ so that, for $1>k>k_{u}$,

[TABLE]

which means $F_{2}>0$ on this interval.

On the other hand, for $k$ in any neighborhood bounded away from $1$, we can use uniformly convergent series expansions to prove $F_{2}$ is positive. First, let

[TABLE]

be the positive fourth roots. One has the identity

[TABLE]

If one defines

[TABLE]

to be the positive fourth root, then one has

[TABLE]

Hence, positivity of $F_{2}$ is reduced to showing the positivity of the much simpler expression

[TABLE]

Since we consider $k\in[0,1-\epsilon]$, we can use the series expansions for $E$ and $K$:

[TABLE]

Since $r_{+}$ is real analytic on $|k|<1-\epsilon$ for any $\epsilon>0$ we can expand

[TABLE]

One easily finds that $c_{0}=1$, $c_{1}=-\frac{1}{2}$, $c_{2}=-\frac{1}{16}$ and $c_{3}=-\frac{1}{32}$. Furthermore, the coefficients satisfy the recurrence relation

[TABLE]

To see this, let h(x)=\Big{(}(1-x)(1-\frac{x}{2})^{2}\Big{)}^{1/4}=\sum_{n=0}^{\infty}c_{n}x^{n}, for which one has

[TABLE]

Hence,

[TABLE]

from which the recursion relation follows. Next, we claim that $2c_{n+1}\leq c_{n}\leq c_{n+1}<0$ for all $n\geq 2$. This follows inductively from the fact that $2c_{3}\leq c_{2}\leq c_{3}<0$ and that the set $2x\leq y\leq x<0$ is invariant under the maps that generate the recursion

[TABLE]

for $n\geq 2$. To see this, suppose that $2x\leq y\leq x<0$, $n\geq 2$ and let $(\tilde{x},\tilde{y})^{T}=\textbf{R}_{n}(x,y)^{T}$. Then it follows

[TABLE]

Since $c_{n}<0$ for $n\geq 1$, $r_{+}$ in $(0,1-\epsilon)$ is smaller than any of its finite truncations $r_{N}:=\sum_{n=0}^{N}c_{n}k^{2n}$. Consequently,

[TABLE]

In particular, $E-r_{+}K$ is positive whenever $h_{N}$ is. The functions $h_{N}$ have an expansion about $k=0$ in $[0,1)$

[TABLE]

Using the expansion of $E(k)$ and $K(k)$ in (53), we have

[TABLE]

Notice that by choosing any $N\geq 4$, the leading term is always given by $\frac{\pi}{2}\tau_{4}k^{8}$, i.e. $\tau_{4}=\frac{3}{2048}$. Furthermore, we found numerically that for any choice of $N$ we always had $\tau_{n}>0$ for all $n\leq N$. We therefore obtain a strategy for a proof as follows. First, considering $N$ large but yet unspecified, split the series into three terms:

[TABLE]

Since $c_{n}<0$ for $n\geq 1$, we have $\tau_{n}>\frac{2n}{1-2n}(P_{2n})^{2}>-(P_{2n})^{2}$ for $n\geq N$. One also has monotonicity $P_{2(n+1)}<P_{2n}$, hence

[TABLE]

The above holds for all $0<k<1$, whereas on the interval $(0,k_{\ell})$ one has

[TABLE]

Note that the first term on the right hand side is always positive for

[TABLE]

hence by choosing $N$ sufficiently large we can make this term positive on the entire interval $(0,k_{\ell})$. One then obtains positivity of $h_{N}$ on $[0,k_{\ell}]$ if one can verify numerically that for $5\leq n\leq N-1$ one has $\tau_{n}>0$.

By way of example, in the case where $k_{\ell}=\frac{1}{2}$ we find that choosing $N=9$ gives \big{(}\frac{3\pi}{4096\cdot K(\frac{1}{2})}\big{)}^{\frac{1}{8}}\approx 0.517, so this covers the entire interval $(0,\frac{1}{2})$. Using Mathematica, we find that

[TABLE]

On the other hand for the interval $k_{\ell}=0.9984$ we used Mathematica to find that for $N=2360$ the interval of positivity in (54) includes all of $(0,0.9984)$, and that $\tau_{n}>0$ for $5\leq n\leq N-1$.