Time Averages and Periodic Attractors at High Rayleigh Number for Lorenz-like Models

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 ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document}, 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 nonzero 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 nonzero 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.


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 ρ 1 for the famous Lorenz equations and variants thereof.The Lorenz equations, given as follows, 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 σ > 0 is the Prandtl number, characterizing the viscosity of the fluid, and where β > 0 is a shape parameter, measuring the ratio of the length of the box to its height.Often of primary interest is the parameter ρ ≥ 0 which is the rescaled Rayleigh number, measuring the intensity of the heating.The rescaling is chosen such that a bifurcation occurs at ρ = 1, and indeed the Lorenz model accurately captures the onset of steady atmospheric convection rolls for ρ > 1.For larger values of ρ, 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 (2) H(ρ, β, σ, X 0 ) = lim sup t→∞ 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 X 0 of the solution (X, Y, Z)(t) of (1).The dependence of the Nusselt number on ρ 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, (a) (b) Figure 1.(a) Plot of the transport H(ρ, β, σ, X 0 ) for β = 8 3 , σ = 10 near the transition to chaos ρ ≈ ρ * depicting the steady convection rate H ± (ρ, β) (orange) and the rate obtained from a numerical simulation of (1) with randomly selected initial conditions X 0 (blue).(b) Plot of the transport gap ∆H(ρ) for β = 8  3 , σ = 10 and for ρ large.The vertical lines mark bifurcations that bound the region of 'large ρ' as in Fig. 5(b), in particular for ρ above the blue line the transport is dominated by stable periodic orbits.
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 ρ = 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 ρ 1 is our main motivation for this paper.The relation to stability is particularly clear in the case 0 ≤ ρ ≤ 1, where for any value of σ, β 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 ρ = 1, a pitchfork bifurcation occurs, where the aforementioned non-zero fixed points X ± = (± β(ρ − 1), ± β(ρ − 1), ρ − 1) emerge.For ρ sufficiently near 1 these non-zero fixed points seem to attract every trajectory except for X 0 belonging to the stable manifold of the origin, W s (0), so that H(ρ, β, σ, X 0 ) = 0 if X 0 ∈ W s (0), H ± (ρ, β) := β(ρ − 1) otherwise.
Here W s (0) is a surface of dimension 2, which means almost all initial conditions give a positive transport, namely that of the fixed points H ± (ρ, β).This is certainly the case for initial data in their non-trivial basins of attraction.
Upon increasing ρ 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 ρ is For 0 < λ ≤ 1, the fixed points X ± are locally stable for all ρ, cf.[23] so that the fixed point transport H ± is observed at least for X 0 belonging to a basin of attraction of positive measure.On the other hand for λ > 1, the fixed points are locally stable only for 1 < ρ < ρ * = σ(σ + β + 3)/(σ − β − 1).As ρ is increased through ρ * , the fixed points X ± lose stability via a sub-critical Andronov-Hopf bifurcation and at least for some open set containing β = 8/3, σ = 10, generic initial conditions give chaotic solutions [25].
As mentioned, it was proven in [22] that despite this complexity, for all ρ > 1 and any σ, β ∈ R >0 , X 0 ∈ R 3 one has the simple bound For the Lorenz equations this bound is actually sharp as it is realized by the steady states X ± .However, since X ± are unstable for ρ > ρ * , λ > 1, the transport that is realized by typical solutions might be much lower.Indeed, numerical experiments presented in [22] indicate that for ρ > ρ * a gap ∆H = H ± − 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 ρ 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 ρ 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 X sym = X sym (ρ, β, σ) for sufficiently large ρ and its dependence on λ > 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 ρ → ∞ 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 ∆H to leading order: We show that H sym (ρ, β, σ) := H(ρ, β, σ, X sym ) is a monotone increasing function of λ, for fixed ρ 1, whose range is a subinterval of (0, H ± ) that limits to this interval as ρ → ∞.Specifically, H(ρ, β, σ, X sym ) ∼ ρ for ρ 1, just as H ± , but with a λ-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 X ± and X sym in terms of λ.This highlights difficulty to recover from low transport when λ grows beyond the 'tipping point' at λ = 1.Moreover, for any γ ∈ (0, 1] we can choose σ(ρ), so that H sym (ρ, β, σ(ρ)) ∼ ρ γ along the periodic attractors X sym .
We employ numerical pathfollowing to corroborate the analytical results for large fixed ρ and find that for large λ 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 (λ, ρ)-plane, and find that it extends to values of ρ below 200.In Fig. 1 this region begins near ρ = 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.1.
Hamiltonian structure and the extended Melnikov theory.It is well known that the limiting system at ε = 0, is integrable with the conserved quantities 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 ≥ 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 ξ = η = 0, ζ = 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 ( 7) For the Jacobi elliptic functions, one first defines the amplitude function am(u, k) as an inverse function via and the Jacobi elliptic functions are then defined , which can be written in terms of Jacobi elliptic functions as Each (A, B) ∈ D 2 defines a pair of asymmetric periodic orbits that can be represented in elliptic functions as The region D 3 : {A = B} corresponds to the line of saddle equilibria ξ = η = 0, ζ = −B, each with a pair of homoclinic orbits L B 3,± = (ξ 3,± (τ ), η 3,± (τ ), ζ 3 (τ )) contained in D 3 and given by u = Due to the reflection symmetry of (5), without loss of generality we consider only L A,B 2 = L A,B 2,+ and L B 3 = L B 3,+ .In the following we suppress the dependence of the elliptic functions and integrals on the elliptic modulus k = k i , e.g., writing sn u for 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 ε = 0, we introduce polar coordinates with B from (6) given by (11) ζ = B cos ϕ, η = B sin ϕ, ξ = ξ which transform (4) into Notably, for ε > 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 (13) for smooth functions f i , g i .In this formulation, at ε = 0, the first two equations possess the Hamiltonian structure with A(ξ, ϕ, B) defined in (6) serving as a Hamiltonian: As noted above, the solutions at ε = 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 ε > 0. The main idea of this method is that for ε = 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 ε > 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 ε = 0 from (ξ, ϕ) to action-angle variables (I, θ).The action variable I is computed as (15) I(A, B) = and the angle θ ∈ [0, 1) increases from an 0 to 1 along the periodic orbits with constant frequency Ω(A, B) given by ( 16) Here and below we use the notation ∂R ∂P 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 (17) For small ε > 0, a trajectory starting from the initial point (I 0 , 0, B 0 ) has monotonically increasing θ from 0 to 1 with velocity ε-close to Ω, and slowly evolving (I, B) with velocity O(ε).Hence, the trajectory stays in an ε-neighborhood of the corresponding unperturbed periodic trajectory, and reaches θ = 1 after finite time T ε = 1/Ω + O(ε), thus returning to the starting plane {θ = 0}.This defines a Poincaré map (I 0 , B 0 ) 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 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 ) = 0 there exists ε 0 > 0 and an isolated fixed point of the Poincaré map (19) in an ε-neighbourhood of (I 0 , B 0 ) for any 0 < ε < ε 0 .This corresponds to a persistent periodic orbit of ( 4) and (1) for 0 < ε 1.Moreover, in case (M 1 , M 3 ) = 0 for 0 < ε 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 (20) in particular (M 1 , M 3 ) = 0 is equivalent to ( M1 , 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 ) ∈ D i there exists small ε 0 > 0 such that for 0 < ε < ε 0 the corresponding periodic orbit lies entirely in D i .Also, it is clear that upon increasing ε 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 ν ± of DM (I 0 , B 0 ) are given by ( 21) 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: 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 + ρ), where B = const and ρ 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 which is incorrect as the term ∂A ∂ρ 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 β = 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(ε) 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) = 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 < ε 1 only one saddle equilibrium (ξ, ϕ, B) = (0, π, σ) 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 ε > 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.

Symmetric periodic orbits (D 1
).In the regime D 1 the explicit solutions (8) have the following form in polar coordinates: Hence, from ( 15) and ( 16) the period and the action are ( 23) Substituting these explicit solutions into (20) using the formulas for f i , g i from ( 13) allows to compute the Melnikov integrals explicitly as Equating these expressions to zero, equivalently gives the equations 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 λ = σ+1 β+2 > 2/3 and ε sufficiently small there exists a (locally) asymptotically stable symmetric periodic orbit.It is the unique periodic orbit in an order ε-neighborhood of an explicit solution L A,B

1
, where (A, B) is the unique solution of (24) for the chosen value of λ.
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 ≤ k 1 ≤ 1.Let e 1 , e 2 be defined by (25) e and since it is clear from the definitions (7) that K ≥ E ≥ 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 ≤ 1 with k * ≈ 0.908909.Since the coefficient of B is strictly positive, one can solve the second equation for B: where d 1 = 4e 1 + βe 2 denotes the denominator of the first fraction.Substituting this into the first equation gives )K) = 0 as the remaining equation.Moving the term involving σ to the right hand side, using the definition λ = (σ + 1)/(β + 2), and after some arithmetic one obtains (27) 2λ .
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 ≤ 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 λ > 2/3.Writing the right hand side of ( 27) as 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 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 → ∞ as k 1 → 1 the right-hand side of equation ( 24) is monotonically decreasing from ∞ to 1/3 when k 1 increases from k * to 1.This range corresponds to the values λ > 2/3 on the left-hand side of (27).Also, every k 1 ∈ (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 λ > 2/3, we now consider the determinant from (22).Computing the derivatives in (22) gives where c0 , c1 , c2 are coefficients depending on σ, β, B only.Using the identity (26) to eliminate B and the identity σ = λ(β + 2) − 1 to eliminate σ, one obtains for coefficients ĉ0 , ĉ1 , ĉ2 , ĉ3 , ĉ4 depending only on λ, β.Finally, applying (27) to eliminate λ yields where This equality expresses F 1 as a sum of strictly positive terms, since for the first two terms we have 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 < ε 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 , so that by computing the derivatives in the formula (22) we find tr DM = 4 3B 3/2 ĉ1 K + ĉ0 E , where ĉ0 , ĉ1 are coefficients depending only on σ, β, B and k 1 (but not E or K).Using the first equation of ( 24) and λ > 2/3 (so that β < 2σ) we express E as and substituting this into the equation for the trace, we find (29) tr Since tr DM < 0 and det DM > 0, from (21) it follows that for ε > 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 A,B 2 (τ ) given in (9).In polar coordinates these are given as ), with period and action (30) Substituting these explicit solutions into (13) and (20) we compute the Melnikov integrals explicitly as Equating these expressions to zero, we see that we must solve 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 ε sufficiently small and all 2/3 < λ < 1 there exists a pair of saddle-type asymmetric periodic orbits.Each is unique in an order ε-neighborhood of L A,B

2
, where (A, B) corresponds to a unique solution (k 2 , B) of (31) for the chosen value of λ.
Proof.The coefficient of B in the first equation of ( 31) is always positive, since so that the unique solution in the terms of B is where d 2 denotes the denominator.Substitution into the second equation of (31) gives an equation for k 2 : (33 . Next, we prove that the right-hand side of (33) monotonically decreases from 1 to 1/3 as k 2 increases from 0 to 1, which implies precisely for every λ ∈ (2/3, 1) there exists a unique solution k 2 .To prove monotonicity, first we compute The three summands of P are each negative since for the first we have for the second K ≥ 0, 2E − K < 0, and for the third we estimate its nontrivial factor from above by the negative Regarding the determinant, by computing the derivatives in (22) we find it has the form where c0 , c1 , c2 are coefficients depending only on σ, β, B, k 2 .Again using the identities (32), σ = λ(β + 2) − 1 and (33) to eliminate B, σ, λ one finds , where F 2 is given by hence it follows that det DM < 0. This proves existence and local uniqueness for 0 < ε 1 by [31, Theorem 3.2] as mentioned above.
Turning now to the question of stability, in region D 2 one has , and hence by computing the derivatives in (22) one obtains where č1 , č0 are coefficients depending only on σ, β, B and k 2 (but not E or K).Using the identity (32) to eliminate B, the identity σ = λ(β + 2) − 1 to eliminate σ, and the identity (33) to eliminate λ we obtain (36) tr Since tr DM < 0 and det DM < 0, by (21), for any sufficiently small ε > 0, there is one eigenvalue inside the unit circle and one outside, thus these periodic orbits are of saddle type.
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, which is nonzero for ε > 0 since β > 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 τ → −∞ along the unstable direction, and lim Let us assume that a generic homoclinic orbit exists for 0 < ε 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 ν 2 = −εβ.For 0 < ε 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 ε.The latter can be computed as . 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 λ = 2/3 when ε → 0.

Implications for infinite time averages
In this section we use the angle bracket notation to denote the infinite time average: 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 X ± for a wide range of parameters (β, σ) and for all 0 < ρ < ∞: Due to the form of the Lorenz equations, certain infinite time averages must be proportional.For instance, one has In this way one has the following equalities and hence it suffices to consider the following three time averages: As stated previously, for λ > 1 the fixed points become unstable for sufficiently large ρ, 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., where T ε 1 is the period of the solution ζ 1 (τ, ε).It seems likely these values do not have a simple closed form and instead we compute the time averages via an expansion in ε.For instance, the lowest order term can be found from the explicit formulas for the unperturbed solution and period, ζ 1 (τ ) and T 1 , as follows In this way, we obtain the following expressions for the infinite time averages, expressed in terms of ρ rather than ε: Recall the expressions e 1 , e 2 , d 1 defined in ( 25), ( 26) were shown to be positive.Hence for the time averages Z and Z 2 , these expressions resolve the coefficient of the leading order term as a function of σ, β which is strictly less than one, hence less than that of the fixed point value.For Z 3 , the term d 1 e 1 is always positive, whereas one can check that the other expression inside the parentheses is positive for all λ > 2.5611..., whereas is is negative for λ less than this.Hence by fixing such λ and choosing β sufficiently large this exceeds the fixed point value.This agrees with Goluskin's result however, since the region in parameter space where Z 3 is maximized at X ± does not include large β.
As mentioned previously, the transport XY 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 and we also compute the transport obtained by the unstable, asymmetric periodic orbits Next, we cast H 1 , H 2 in terms of ρ and rescale to a finite range of transport values.This gives We first show that h ρ,1 (2/3) = h ρ,2 (2/3) = 0 understood as the limit λ 2/3, and analogously = 0 due to the following.Figure 2. Bifurcation diagrams (solid=stable, dashed=unstable) of relevant equilibria and periodic orbits in terms of λ and the rescaled transport, illustrating the hysteresis.(a) equilibria in the averaged planar system at ρ = ∞ that persist for large finite ρ as equilibria (orange) or periodic orbits (blue).
Computations are done via the elliptic integrals with Mathematica.The bullet marks λ = 2/3 at zero transport; for illustration, the thin blue line extends numerical computations to the theoretical limit at zero transport.(b) Overlay of (a) with branches of stable or unstable symmetric (red solid/dashed) and unstable asymmetric (red dashed) periodic solutions to the full system at ρ = 1000 computed with Auto [5].
0 and using (32) as well as The differences to the scaled maximum transport In particular, the stable symmetric periodic orbits L yield the same order of magnitude of transport with respect to ρ, but feature a λ dependent downshift that vanishes at λ = ∞, i.e. at X ± .In addition, from our perhaps rough error estimates we obtain a correction of order (at least) ρ −1/2 compared to the next order being ρ −1 for X ± .Numerical computations suggest that this term might in fact be of order ρ −1 , but we do not explore this further here.

Numerical computations and hysteresis loop
We present numerical results that corroborate the analytical results for ρ = ∞ and 1 ρ < ∞ 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 ≈ 0.6.The analytical prediction is that the branches terminate at λ = 2/3 in homoclinic bifurcations of the zero equilibrium and thus at zero transport.Indeed, at λ = 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 λ: For λ < 2/3 the equilibria X ± that maximize transport are stable, while for λ > 1 the symmetric periodic orbit are.Intermediate values 2/3 < λ < 1 lie in the analytically predicted region of bistability with stable equilibria X ± and stable symmetric periodic orbit L.
For large finite ρ and moderate values of λ, 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 λ = 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 ρ = 1000 this occurs at λ = λ bp (ρ) ≈ 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 λ = λ pd (ρ) ≈ 0.787 in a period doubling bifurcation.We plot the loci of λ bp , λ pd in Figure 5(b), showing that as ρ increases, λ bp , λ pd approach λ = 2/3.Further numerical simulations corroborate the hysteresis-type loop: for λ < 2/3 the maximum transport equilibria X ± appear to be global attractors, while for λ > 1 this seems to be the stable symmetric periodic orbit, as in Fig. 3(b).See also Fig. 1(b).For 2/3 < λ < 1 the situation with large finite ρ is complicated by the fact that symmetric periodic orbits are born in a homoclinic bifurcation at some λ hom ∈ (2/3, 1) and, as mentioned, are unstable until a bifurcation point λ bp ∈ (λ hom , 1).Up to the aforementioned region of stable asymmetric periodic orbits that bifurcate from λ bp , the global attractors for λ < λ bp seem to be X ± and the region of bistability with the symmetric periodic orbit is effectively λ ∈ (λ bp , 1).
Using time-varying values of λ we consistently found hysteresis as plotted in Figure 4(a): for slowly increasing λ from 0, the solution is quickly close to 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 λ.As λ 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 λ = 1 at t = 100, but with a delay, here until around t = 190.Subsequent decrease of λ causes the solution to track the stable branch of symmetric periodic orbits, cf. Figure 4(b), which decreases the observed local transport further until λ = λ bp ≈ λ pd .Upon decreasing λ below this threshold, a switch to a stable equilibrium X ± occurs, thus re-creating maximum local transport.The red curve corresponds to the extension of the branch of stable symmetric periodic orbits from Figure 2. The magenta curve is the analogue for ρ = 4000.Continuing along these branches from their lower left ends, numerically before each fold a destabilising period doubling bifurcation occurs and the solution restabilises at the fold point.
Each branch appears to terminate in a symmetric heteroclinic cycle between the pair of equilibria X ± .(b) loci of branch points of symmetric periodic orbits (blue) and period doubling points on the bifurcating branches (purple); for values of ρ above the blue curve periodic orbits appear to be stable.Gray lines mark ρ = 1000 (horizontal) and λ ≈ 2.36 at the classical Lorenz values σ = 10, β = 8/3.
While the asymptotically predicted branch of symmetric periodic orbits of Fig. 2 continues for increasing λ monotonically and unboundedly, we found that for finite ρ 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 X ± at a finite value of λ.See also Fig. 3(c).Upon increasing ρ, this turning and termination occurs at larger values of λ.Hence, this scenario is consistent with the analytical results, which concern ρ → ∞ for bounded ranges of λ.The appearance of a symmetric heteroclinic cycle between X ± in the Lorenz system has already been noticed in [23,8], albeit apparently not in the regime of large ρ.
The transport at such a heteroclinic cycle is that of the symmetric equilibria, i.e. β(1−ρ −1 ), which is indeed very closely matched at the numerical termination points.The (λ, h ρ )-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 λ; 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 ρ = ∞ 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 ρ → ∞.

Other Lorenz-like systems
The analysis for large ρ 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 with w ∈ R k , linear A, B and constant b.Upon rescaling as in (3) and w = ε −j ω, with Ξ = (ξ, η, ζ) we obtain the form where Φ ε is the right hand side in (4).For A = B = 0, i.e., in absence of ω, the difference to ( 4) is of order ε 2 .Hence, the leading order analysis is unchanged, which means that periodic orbits bifurcate/persist as for 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 ≥ 1 for a regular limit in which the right hand side of the equation for ω becomes independent of ω, and vanishes for j > 1.For j = 1 we obtain, up to terms of order ε 2 , (46) with A 1 the first row of A. For B of the form B = [B 1 |0|0|B 2 ] the equation for ω has the slow form which occurs with w ∈ R in the Lorenz-Stenflo model from [24], its magnetic variant [26], and with w ∈ R 2 in the models from [12]; for the latter we choose α = O(ε) and shift the auxiliary variables (which gives b = 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 ε −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 ω to the action-angle formulation.The additional Melnikov-integral term M 3 is then simply the integral of B 1 ξ(t) + B 2 ω 0 over the period T , with ω 0 constant.Since ξ has zero average ( φ = ξ at ε = 0 in (12)), this term becomes ω 0 /T , with the period T , so that M 3 = 0 requires ω 0 = 0.This means that the values of the other two Melnikov-integrals for (46), M1 , M2 , actually coincide with those of M 1 , M 2 from §2.The non-degeneracy condition turns into invertibility of the matrix where M2 is independent of ω, and it turns out that also M1 is: In its integrand F from ( 17), the additional term from A 1 ω 0 is constant and has a factor ∂I ∂ξ = ∂I ∂A ∂A ∂ξ = T ξ, where ξ has zero average as noted above.Hence, the matrix has lower left triangular block structure and the block 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: This system is a mode truncation of the rotating Boussinesq equations: 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 and we obtain the system of equations (50) The limiting system when = 0 coincides with (4) except trivial dynamics in the variable χ, so that the system now admits three invariants of motion (51) Using the first two invariants as for the Lorenz system, (50) can be solved at ε = 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 and (50) becomes Hence, B and χ are constant at ε = 0, and for fixed B and χ, the remaining system for (ξ, φ) possesses the same Hamiltonian structure as (5).Converting to action angle coordinates the system becomes where In this case we have a three dimensional Melnikov function given by 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 is satisfied.Since we can write Since we aim at illustration, we consider |A| ≤ B only, and then compute 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 χ 0 = 0, and M 1 , M 2 are independent of χ 0 .Hence, and the determinant is given by which is non-zero as shown in §2.

Discussion
In this paper we have revisited the dynamics of the Lorenz equation in the regime of large Rayleigh number ρ, 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 ρ > ρ * , λ > 1.However, we found that the literature concerning existence and stability theory of periodic states for large ρ was incomplete.We have therefore provided a rigorous treatment, which essentially confirms the predictions of [23].Numerical computations for large finite ρ 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 ρ we have identified a hysteresis loop in terms of the parameter λ = σ+1 β+2 , which illustrates difficulty to recover from a loss in transport once λ exceeds the 'tipping point' λ = 1.Moreover, we have computed the stability boundary of periodic attractors in the (λ, ρ)-plane and found that it extends to relatively low values of ρ below 200.For fixed ρ 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 Nu = 1 + 2 βρ H(ρ, β, σ, X 0 ), cf.[22], is bounded by 3 as ρ → ∞ due to the transport bound from [22].However, this is far from the 'ultimate' or 'classical' Nusselt number bounds of order ρ 1/2 or ρ 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 ρ regime, we have analytically determined stable structures and their transport.
(1) 1 ≥ k ≥ k1 , where k1 ≈ .9266:In this case the integrand tends to infinity sufficiently quickly as t → 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) ≤ t 0 (1) = 1/2 one has On the other hand, letting t 2 be defined as in (52), note Hence we have I + + I − > 0 as long as t 2 ≤ 3/4.But, noting that 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) ≤ t 0 ( k2 ), hence one has where k * ≈ .9089 is defined as the value of k such that K(k * ) = 2E(k * ).Note that for all k ≤ 1 one has Note that the polynomial is positive for all x > r 1 ≈ 2.34305.Note that K/E is monotonically increasing, and K/E = 2.34305 when k = k1 ≈ 0.949509.On the other hand, when k ≤ 0.94951, one has But the polynomial p 2 (x) = 1.84686x 2 − 6x + 5 is non-zero for all x, thus proving the bound.
A.3.Third elliptic integral expression.Here we show that satisfies F 2 (k) > 0 for all 0 < k < 1 in two steps, first for 0.998357 ≈ k u < k < 1 and next for 0 < k < k = 0.9984.
For the interval k u < k < 1, note that the 2nd and 3rd terms of F 2 can be rewritten as This is strictly positive for all 0 < k < 1: First, for 0 < k < 2 3 one has 3(2 − k 2 )K − 4E > 4(K − E) > 0 Second, using Mathematica, Therefore, Since (1 − k 2 )K 3 equals π 3 8 at k = 0, monotonically decreases for 0 < k < 1, and tends to zero as k → 1, there is k u with Using Mathematica, we find k u ≈ 0.998357 so that, for 1 > k > k u , 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 be the positive fourth roots.One has the identity ) 2 K 4 ).

If one defines
to be the positive fourth root, then one has x y for n ≥ 2. To see this, suppose that 2x ≤ y ≤ x < 0, n ≥ 2 and let (x, ỹ) T = R n (x, y) x ≤ x = ỹ.
Since c n < 0 for n ≥ 1, r + in (0, 1 − ) is smaller than any of its finite truncations r N := N n=0 c n k 2n .Consequently, In particular, E − r + K is positive whenever h N is.The functions h N have an expansion about k = 0 in [0, 1) Using the expansion of E(k) and K(k) in (53), we have Notice that by choosing any N ≥ 4, the leading term is always given by π 2 τ 4 k 8 , i.e. τ 4 = 3 2048 .Furthermore, we found numerically that for any choice of N we always had τ n > 0 for all n ≤ N .We therefore obtain a strategy for a proof as follows.First, considering N large but yet unspecified, split the series into three terms: The above holds for all 0 < k < 1, whereas on the interval (0, k ) one has Note that the first term on the right hand side is always positive for (54) hence by choosing N sufficiently large we can make this term positive on the entire interval (0, k ).One then obtains positivity of h N on [0, k ] if one can verify numerically that for 5 ≤ n ≤ N − 1 one has τ n > 0.

Figure 3 .Figure 4 .
Figure 3. Profiles of periodic solutions for ρ = 1000 from the red solid branch in Figure 2(b).(a) near the double homoclinic loop at the left termination point; (b) at λ ≈ 1; (c) near an apparent symmetric heteroclinic cycle λ ≈ 11.3.

Figure 5 .
Figure 5. (a) Bifurcation diagram for larger values of λ analogous to right panel of Figure2but without stability information of the periodic solutions (long dashed lines).The red curve corresponds to the extension of the branch of stable symmetric periodic orbits from Figure2.The magenta curve is the analogue for ρ = 4000.Continuing along these branches from their lower left ends, numerically before each fold a destabilising period doubling bifurcation occurs and the solution restabilises at the fold point.Each branch appears to terminate in a symmetric heteroclinic cycle between the pair of equilibria X ± .(b) loci of branch points of symmetric periodic orbits (blue) and period doubling points on the bifurcating branches (purple); for values of ρ above the blue curve periodic orbits appear to be stable.Gray lines mark ρ = 1000 (horizontal) and λ ≈ 2.36 at the classical Lorenz values σ = 10, β = 8/3.

Figure 6 .
Figure 6.Projections into the (A, B)-coordinate plane of periodic solutions for ρ = 1000 (blue) and the diagonal (orange).(a) Near the double homoclinic loop at the left termination point of the red solid branch in Figure 2(b).(b) From the red solid branch of Figure 2(b) at λ ≈ 1. (c) From the red long dashed branch of Figure 5(a) at λ = 7.2, towards the heteroclinic cycle.(d) Near the heteroclinic cycle from Figure 2(b).
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 ≥ k ≥ k1 where k1 is defined such that Note this region is deliberately chose to overlap with the region in the first case.This is easily seen, since for instance one has 2+ c k2,2 ≈ .0026> 0 T. Then it follows