Periodicity and Chaos of High-Order Lorenz Systems

High-order Lorenz systems with five, six, eight, nine, and eleven equations are derived by choosing different numbers of Fourier modes upon truncation. For the original Lorenz system and for the five high-order Lorenz systems, solutions are numerically computed, and periodicity diagrams are plotted on (σ,r) parameter planes with σ,r ∈ [0, 1000], and b = 8/3. Dramatic shifts of patterns are observed among periodicity diagrams of different systems, including the appearance of expansive areas of period 2 in the fifth-, eighth-, ninth-, and 11th-order systems and the disappearance of the onion-like structure beyond order 5. Bifurcation diagrams along with phase portraits offer a closer look at the two phenomena.


Introduction
In 1963, Lorenz [1963] discovered the chaotic behavior associated with the numerical solutions of the following system of three ordinary differential equations describing a simplified model of Rayleigh-Bénard convection: where X 1 , Y 1 , and Z 1 are the variables, and the dot denotes the derivative with respect to time t. There are three parameters here: namely, σ, r, and b. σ ν/κ is the Prandtl number defined as the ratio of the kinematic viscosity ν to the thermal diffusivity κ. The normalized Rayleigh number, r = Ra/Ra c , is the ratio of the Rayleigh number Ra to the critical Rayleigh number Ra c . Finally, b = 4π 2 /(π 2 + k 2 ) is the geometric parameter, where k is the horizontal wavenumber. The three equations, now called the Lorenz equations, have since spurred the development of the theory of chaos. Today, chaos theory finds its importance not only in a purely mathematical sense but also through its applications. For instance, the behavior of Earth's atmosphere is considered to exhibit chaotic [Palmer, 1993] and subharmonic behavior [Feingold et al., 2010] and due to this qualitative resemblance, Lorenz-like equations have been proposed as models for a variety of atmospheric phenomena [Vallis, 1986;Stenflo, 1996;Lorenz, 2005;Roberts et al., 2016;Koren et al., 2017]. Chaos theory is also used by some meteorologists in evaluating predictability of weather and climate models [Palmer, 1993;Shukla, 1998;Yang et al., 2006].
Despite its potential usefulness, directly taking the third-order Lorenz system as the model equations for a natural convection requires caution because much information is lost through the simplification processes involved in obtaining the ordinary differential equations. Of particular interest is the way in which Lorenz [1963] chose to truncate Fourier expansions of the streamfunction ψ and the temperature perturbation θ. Our primary motivation is to see what happens if we loosen the severity of this simplification by including additional Fourier modes upon truncation. The present study aims to investigate qualitatively how this change in truncation influences the complexity of the system, especially in the context of patterns made by possible chaotic and periodic regions in parameter spaces. The underlying assumption is that the more Fourier modes are included, the closer the system gets to the original governing equation system for Rayelgh-Bénard convection.
Choosing additional modes from Fourier series expansions of ψ and θ allows for a formulation of systems with varying numbers of equations. We hereinafter refer to these extensions of the original Lorenz system as high-order Lorenz systems. For instance, in the late 1970s, Curry [1978] jumped straight to a system of 14 equations because 14 was the largest number of equations that his machine could handle. Although a system of a very high order is an interesting subject in its own right, we are more interested, as was Shen [2014Shen [ , 2015, in the changes in behavior of the systems with respect to an increasing order. For any given order n, there is by no means a unique set of equations that can be called the Lorenz system of order n; depending on the choice of Fourier modes, one may sometimes be able to form many different systems of the same order. Roy and Musielak [2007a, 2007b, 2007c derived and compared Lorenz systems of different orders by taking into account physical viability. The fifth-and the sixth-order Lorenz systems selected in this way were further analyzed by Shen [2014Shen [ , 2015. In particular, Shen [2014Shen [ , 2015 looked for critical r values beyond which chaos would ensue. It was decided that the fifth-and the sixth-order systems possess similar critical r values, both of which are greater than the critical r value of the original Lorenz system. Critical r values, however, also depend on the choice of values of other parameters in the system such as σ and b; therefore, the behavior of high-order Lorenz systems may be better understood if pairs of parameters are considered instead of r alone. To this end, some measure of chaos and periodicity corresponding to each pair of two parameters can be plotted on a twodimensional parameter plane. Rech [2016] used the Lyapunov exponent for this purpose regarding a related system called Lorenz-Stenflo equations. We instead adopt periodicity diagrams following Dullin and Schmidt [2007], Park et al. [2015] and Park et al. [2016b].
In this paper, we consider the original Lorenz system and five high-order Lorenz systems, ranging from order 3 to order 11. We begin our discussion by introducing the equations in the next section. These equations are obtained by raising the number of modes in Fourier expansions of ψ and θ in an alternating order. The numerically computed measure of chaos for each system is plotted on the (σ, r) parameter plane in Sec. 3, followed by a more detailed look at some of the patterns that appear in the periodicity diagrams.

High-Order Lorenz Systems
Following Saltzman [1962], we write the following nondimensionalized form of the governing equations describing Rayleigh-Bénard convection in the x z plane: where ∇ 2 2 / x 2 + 2 / z 2 and ∇ 4 4 / x 4 + 4 / z 4 + 2 4 / x 2 z 2 . Let a = 1/ √ 2 be a number. Applying Fourier series expansions to the streamfunction ψ and the temperature perturbation θ yields  from which we construct systems of (N + 2M ) ordinary differential equations by selecting appropriate positive integers for M and N with the condition that 0 ≤ M − N ≤ 1. For example, upon simplification, taking M N = 1 yields the original third-order Lorenz system. The choice of N and M in (6) and (7) can be arbitrary at this stage. As mentioned in [Roy & Musielak, 2007a], there are many ways to choose different pairs of Fourier modes; however, since we are more interested in incremental changes as M and N are raised one at a time, the condition 0 ≤ M − N ≤ 1 is imposed. The order of the system obtained in this way is (N + 2M ) because there are two terms in each summand involving M and only one term in each summand involving N . Note that the same method does not produce desired ordinary differential equations without sinusoidal terms when N > M. For instance, forcing N = 2 and M = 1 in an attempt to obtain a fourth-order system does not yield an ordinary differential equation system in the same way that the other high-order systems are obtained. Usually, some terms in (6) get paired up with terms in (7) in order to eliminate all of the sinusoidal elements using trigonometric identies. In this case, however, the extra sinusoidal terms in (6) due to having N > M have no counterpart in (7) and some multiples of sines and cosines remain in the equations. The fifth-and the sixth-order systems are in agreement with Shen [2014Shen [ , 2015, and he, likewise, jumps from the third-order system to the fifth order [Shen, 2014]. The order can be raised indefinitely, but we focus on five additional systems of 5, 6, 8, 9, and 11 ordinary differential equations. In all of the following equations, parameters d n are defined by

A Lorenz system of order 5
To obtain a system of five equations out of the governing equations, we take M = 2 and N = 1 in (6) and (7). After plugging the approximations for ψ and θ back into (4) and (5) and applying trigonometric identities, we obtain

A Lorenz system of order 6
Since we raised M by 1 going from order 3 to order 5, we now raise N by 1 in turn to obtain a sixthorder system: Note that the fifth-and the sixth-order systems here are identical to those derived by Shen [2014Shen [ , 2015.

A Lorenz system of order 8
We raise M to 3 in the Fourier series while keeping N = 2 as in the case of the sixth-order system above. This gives rise to a system with two more equations as follows:

A Lorenz system of order 9
Take M = 3 and N = 3 to obtain a system of nine equations, namely,

A Lorenz system of order 11
Finally, we derive an 11th-order system by setting M = 4 and N = 3:

Some properties of high-order Lorenz systems
Following the discussion in [Strogatz, 2015], we examine whether the high-order systems listed above also possess some of the properties belonging to the original Lorenz system such as nonlinearity, volume contraction, and symmetry. The first property, nonlinearity, clearly holds for our systems. In addition to having the nonlinear terms inherited from the original Lorenz system, our high-order systems also include new nonlinear terms. Since it is thought that chaotic behavior arises from the inclusion of nonlinear terms [Shen, 2014], the important question regarding high-order systems is whether the combined effect of new and old nonlinear terms generally results in chaos. The conclusion made by Shen [2014Shen [ , 2015 is that chaotic solutions do appear in the fifth-and the sixth-order systems, but the onset of chaos may require a different set of values of parameters r, σ, and b [Shen, 2014]. This conclusion will be verified for higherorder systems in the next section.
By volume contraction, we mean the dissipative nature of our systems. Given any system Ẋ f (X), the system is dissipative if it has a negative divergence: i.e.
For instance, Lorenz [1963] had shown that the original Lorenz system satisfies this condition as That all our systems satisfy (48) is easily verified. In addition, we find that the magnitude of the divergence increases with order. For example, the divergence of the ninth-order system is given by which is less than the divergence in (49). Since each of the parameters in the above equation is a positive number, it follows that the result in (50) is less than the divergences computed for each of the systems of orders less than nine. In this sense, we can say that the greater the order of our system is, the faster the volume contraction occurs. Finally, the original Lorenz system is thought to be symmetric due to its invariance under the transformation (X, Y ) (−X, −Y ). Likewise, our high-order Lorenz systems remain invariant under

Methodology
The fourth-order Runge-Kutta method is employed to perform numerical integrations of each system with ∆t = 10 −4 . For all our numerical computations, the initial conditions are set to 0 for all variables except for X 1 which is set to 1. Figure 1 shows the time series of the variable Z 1 for all six cases on t ∈ [0, 40] with the parameter values given by Lorenz [1963]: r = 28, σ = 10, and b = 8/3. In Fig. 1, chaotic behavior is visible only in the third-order case. This is expected since the onset of chaos requires r values greater than 28 in highorder systems [Shen, 2014]. Notice that the amount of time it takes for the solution to converge to a fixed point is shorter for the 11th-order system [ Fig. 1(f)] than it takes for the fifth-order system [ Fig. 1(b)]; nevertheless, such a decrease in the time of convergence is not monotonic with respect to increasing the order. An analogous observation for chaotic solutions would be measuring how fast a slight perturbation in the initial condition gets amplified. Two solutions with slightly different initial conditions would start deviating from one another significantly if they are chaotic, and yet again numerical experiments show no monotonic relations between the order of a system and how quickly such a deviation occurs.
The aforementioned issues dictate that any comparison between Lorenz systems with different orders must be made based on their behavior after all such convergences or deviations take place; thus, from now on, we limit our observations to some time domain I such that I ⊂ [200, ∞). We assume that after the t 200 mark, the time series are representative of the behavior of solutions as t ∞. Plots in Fig. 2 are the time series of Z 1 on t ∈ [210, 212] with parameter values r 300, σ = 20, and b = 8/3. Once time series plots such as Figs. 1 and 2 are obtained, the solutions are identified as chaotic, periodic, or convergent. A solution is convergent if there exists a constant p such that lim t→∞ Z 1 (t) = p. From a numerical point of view, however, the above limit is taken with t 200. Any nonconvergent time series of Z 1 would then consist of some combination of peaks and troughs. Intuitively, a periodic solution would exhibit a repeating pattern of peaks and troughs while no such repetition of patterns would be found in chaotic solutions. These ideas are made more precise and fit for our numerical computations with the following definitions. Given the ith peak in a time series, where i = 1, 2, 3, . . . , define Z max,i to refer to the local maximum value of Z 1 corresponding to the ith peak. To account for numerical errors, we consider the ith and the jth peaks to be equivalent if |Z max,i −Z max,j | ≤ Z max,i /1000; that is, two peaks are equivalent if the difference between peak heights of the two is less than or equal to 0.1% of the first peak. If, in a time series of a nonconvergent solution, the nth peak is followed by a peak equivalent to the first peak of the series, then this n is said to be the solution's peak number. A solution with peak number greater than 16 is considered chaotic and any nonconvergent solution with peak number ≤ 16 is considered periodic. A periodic solution with peak number n is said to be of period n.
Although the above criterion for distinguishing between chaos and periodicity may seem somewhat arbitrary, the underlying assumption here is that it is not likely for 17 or more distinct peaks to reappear in exactly the same order. Strictly speaking, the peak number for any chaotic solution would necessarily be ∞ and any finite peak number would imply a periodic solution. In this sense, the higher its peak number, the closer to chaos a solution gets. Peak numbers, therefore, function as a measure of chaos and thus as an alternative to Lyapunov exponents that were used for this purpose in [Rech, 2016]. Figure 3 shows the periodicity diagrams for our Lorenz systems. We focus on parameters σ and r with σ, r ∈ [0, 1000]. Peak numbers are plotted on (σ, r) parameter planes by considering all possible integer pairs within the parameter domains. In each diagram, the numerically computed peak numbers are assigned different colors following the color scheme similar to that of Park et al. [2015] and Park et al. [2016a]. The periodicity diagrams in Fig. 3 change unpredictably from one system to next. For instance, the shapes and sizes of the period 1 regions vary widely depending on the order of the system, although it seems that overall there is an enlargement of the period 1 region associated with the increase in order. The periodicity diagram for the sixth-order system in Fig. 3(c) particularly stands out from the rest in that it exhibits a spectrum of different peak numbers along the horizontal line r ∼ 390, where σ ranges from 150 to 600. Because of such dramatic differences in periodic diagrams of systems of different orders, it is entirely possible even with a fixed value pair of (σ, r) that the solutions take on chaos, periodicity, and convergence depending on the order of the system. Such a development is also seen in Fig. 2, where the solutions are deemed chaotic in systems of orders 3, 5, 6, and 11, but periodicity is observed in systems of orders 8 and 9. As in [Dullin & Schmidt, 2007], the periodicity diagram of the original Lorenz system in Fig. 3(a) exhibits the onion-like structure resulting from the alternations between chaos and periodicity. Such a structure persists in the fifth-order case as well. One of the most notable changes going from order 5 to higher orders in Fig. 3 is the disappearance of the onion-like structure. This disappearance of white chaotic bands in systems of orders greater than 5 results in red isolated windows of period 1 in the lower left areas of Figs. 3(d)-3(f), where σ is between 50 and 150 and r ranges from 50 to 850. Disconnected from other regions of period 1, these red windows of period 1 are bounded by either chaotic regions or periodic regions of greater peak numbers.

Periodicity diagrams
Another notable feature in the periodicity diagrams is the appearance of relatively wide orange regions of period 2 in the systems of orders 5, 8, 9, and 11. Note that compared to the orange regions in Figs. 3(b) and 3(d)-3(f), the regions of period 2 in Fig. 3(a) appear merely as thin orange bands alongside the white chaotic bands that form the onionlike structure. For the fifth-order system in Fig. 3(b) in particular, the periodicity diagram exhibits both the widened orange regions and the onion-like structure consisting of chaotic bands.

Bifurcation diagrams and phase portraits
The difference between the wide orange areas and the thin orange bands is also visible in bifurcation diagrams. Each bifurcation diagram in Fig. 4 (a) consists of all local maxima of Z 1 on the parameter domains r 400 and σ ∈ [0,200]. In each bifurcation diagram, the initial onset of chaos is indicated by dense scatters at lower σ values. These concentrations of scatters correspond to the leftmost chaotic bands appearing in all six periodicity diagrams in Fig. 3. The subsequent concentrations of chaotic scatters to the right, on the other hand, appear to split into two similar but distinct blobs in Figs. 4(b), 4(d) and 4(e), followed by two nearly parallel lines signaling solutions of period 2. This separation occurs around σ = 140 for the system of order 5 [ Fig. 4(b)] and around σ = 110 for systems of orders 8 and 9 [Figs. 4(d) and 4(e)]. The 11thorder system behaves similarly to the eighth-and the ninth-order systems in that such a split is also visible in Fig. 4(f), but in this case, parallel lines also appear to the right of the initial chaotic scatter around σ = 70. The separated chaotic blobs for the 11th-order system appear briefly around σ 100 followed by again a pair of nearly parallel lines. The thin orange bands of period 2 in Fig. 3(a) seem to be, then, of fundamentally different nature from the expanded orange regions of period 2 in Figs. 3(b) and 3(d)-3(f). The quickly merging lines that follow immediately after a chaotic scatter in bifurcation diagrams correspond to the thin orange bands of period 2 in the periodicity diagram of the original Lorenz system [ Fig. 3(a)], whereas the nearly parallel lines that follow after two distinct blobs of chaotic scatter correspond to the widened orange areas of period 2 shown in Figs. 3(b) and 3(d)-3(f).
Bifurcation diagrams, however, offer little insight regarding the onion-like structure of white bands in the third-and the fifth-order cases. To visualize the effects of the onion-like structure, we  (190,400), respectively. b is fixed as 8/3. These (σ, r) pairs are chosen so that the first column represents phase portraits for systems of period 1, the second column for chaotic systems, the third column for systems of period 2, and the last column for systems of period 1.

= =
= -= turn to phase portraits. Figure 5 shows trajectories projected on the X 1 Z 1 plane. Different pairs of σ and r values are picked so that each column of Fig. 5 would consist of representative phase portraits corresponding to regions of different peak numbers that border the orange regions in periodicity diagrams. The top row of Fig. 5 corresponds to the third-order case and the middle row to the fifth-order case. Since the behavior of the systems of orders 8, 9, and 11 are similar, we only display trajectories of the eighth-order case in the bottom row of Fig. 5. Figures 5(d), 5(e), 5(l) show that if a trajectory forms a simple closed loop that has only one local maximum value of Z 1 when projected on the X 1 Z 1 plane, such a solution is of period 1. On the other hand, in Figs. 5(a), 5(h) and 5(i), it is the symmetry of the trajectories that makes these solutions to be interpreted as period 1. It is, therefore, possible for some trajectories with different structures to be associated with the same peak number. It follows that a low peak number does not necessarily require a simple closed loop trajectory.
Indeed, some transitions from period 2 to period 1 are represented by an untangling of trajectories as we can see in the transitions [5(c) 5(d)], [5(g) 5(h)] and [5(k) 5(l)]. If a trajectory for period 2 is a simple doubly tangled loop resembling Fig. 5(c) or Fig. 5(k), its transition to period 1 requires only one instance of untangling. On the contrary, a relatively complex trajectory such as Fig. 5(g) may need more than an untangling in order to transition into a simple closed loop. For example, in the fifth-order case, some loops do not get untangled entirely going from period 2 to period 1 as can be seen in the transition from Fig. 5(g) to Fig. 5(h). Despite having period 1, Fig. 5(h) is still more complicated than the simple closed loops of period 1 in other cases such as in Figs. 5(d) and 5(l). Whereas Figs. 5(c) and 5(k) are merely two different representations of the unknot, it is unclear whether the same holds for Figs. 5(g) and 5(h). In fact, when r is fixed as 400 in the fifthorder system, simple closed loops start to emerge around σ 900. This is where the (σ, r) pairs belong to an outer layer of the onion-like structure. It seems that for systems exhibiting the onion-like periodicity diagrams, the periodic trajectory loops are simpler for the pairs in the outer layers of the onion-like structure, where the difference between the two parameters is relatively large. Deeper into the onion-like structure of the fifth-order system, even the trajectories associated with low peak numbers are relatively complex [ Fig. 5(h)]. Such a development of rather complicated periodic trajectories for parameter pairs deep in the onion-like structure is also seen in the third-order case. For instance, the trajectory of the solution of the third-order system at r = 600 and σ = 400 is similar to the trajectory shown in Fig. 5(h).

Summary and Conclusion
Starting with the original Lorenz system, the present study analyzed qualitatively the behavior of high-order Lorenz systems of orders 5, 6, 8, 9, and 11. Based on the time series plots of the numerical solutions for the variable Z 1 , we concluded that the amount of time it takes for a solution to converge to an equilibrium or for solutions with different initial conditions to start diverging follows no particular pattern. Focusing on different pairs of (σ, r) parameter values, we plotted periodicity diagrams for all six systems. Owing to the remarkable diversity of the periodicity diagrams for our high-order Lorenz systems, it is possible for regimes of chaos, periodicity, and convergence to emerge in all six systems with appropriate parameter pair choices. Some characteristic changes observed in the periodicity diagrams include expanded regions of period 2 for systems of orders 5, 8, 9, and 11, and the disappearance of the onion-like structure of chaotic regions for systems of orders greater than 5. We utilized bifurcation diagrams to confirm that the widening of period 2 in the periodicity diagrams of high-order systems is a new phenomenon that is not observed in the original Lorenz system. By comparing phase portraits, we found that the trajectories of the fifth-order system are relatively complex when the parameter pairs are deeply embedded in the onion-like structure of its periodicity diagram. The same holds for the onion-like structure seen in the periodicity diagram of the third-order system. The trajectories of the periodic solutions of the third-and the fifth-order systems were contrasted with the trajectories associated with the periodic solutions of the systems of orders 8, 9, and 11, which form relatively simple representations of the unknot.
Our results show that the behavioral patterns of high-order systems in periodicity diagrams differ from the original third-order Lorenz system and 1750176-10 Int. J. Bifurcation Chaos 2017.27. Downloaded from www.worldscientific.com by 84.92.122.48 on 08/17/20. Re-use and distribution is strictly not permitted, except for Open Access articles.
--= from one another. Such behavioral changes are not monotonic with respect to the order of the system, and an increase in complexity of the equations does not necessarily lead to an increase in the system's complexity.
Although the present study considered Lorenz systems up to order 11, some of the new phenomena in high-order Lorenz systems may be studied with greater clarity if systems of even higher orders are considered. Derivations via Fourier mode truncation get quite cumbersome as the number of equations grows; thus, it may be of interest in future studies to fully generalize the Lorenz systems to systems of n equations derived from the governing equations (4) and (5), so that given any order n, corresponding Lorenz systems would be readily available. Furthermore, an expansion or a refinement of the parameter planes, or an investigation with some other pairs of parameters may also prove fruitful.