A physically extended Lorenz system

The Lorenz system is a simpli ed model of Rayleigh-Bénard convection, a thermally driven uid convection between two parallel plates. Two additional physical ingredients are considered in the governing equations, namely, rotation of the model frame and the presence of a density-a ecting scalar in the uid, in order to derive a six-dimensional nonlinear ordinary di erential equation system. Since the new system is an extension of the original three-dimensional Lorenz system, the behavior of the new system is compared with that of the old system. Clear shifts of notable bifurcation points in the thermal Rayleigh parameter space are seen in association with the extension of the Lorenz system, and the range of thermal Rayleigh parameters within which chaotic, periodic, and intermittent solutions appear gets elongated under a greater in uence of the newly introduced parameters. When considered separately, the e ects of scalar and rotation manifest di erently in the numerical solutions; while an increase in the rotational parameter sharply neutralizes chaos and instability, an increase in a scalar-related parameter leads to the rise of a new type of chaotic attractor. The new six-dimensional system is found to self-synchronize, and surprisingly, the transfer of solutions to only one of the variables is needed for self-synchronization to occur. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5095466 Rayleigh-Bénard convection is a thermally driven uid convection with a variety of examples seen throughout nature. The Lorenz system is a simpli edmodel ofRayleigh-Bénard convectionwhose importance lies not only in understanding the uid convection problem but also in the context of the modern development of chaos theory. The present study extends the original Lorenz system by considering two additional physical ingredients, namely, a density-a ecting scalar and rotation of the model frame. Presented in this report are a physically extended Lorenz system and the analysis of the new system including its self-synchronization. This study contributes to the mathematical eld of nonlinear dynamics and chaos theory and is an important step toward bringing the existing Lorenz models closer to reality.


I. INTRODUCTION
Seven decades after Henri Poincaré's initial glimpse at chaos in the three-body problem in 1893, 1 Edward Lorenz discovered a strange attractor in the numerical solutions to a deceptively simple set of three ordinary di erential equations (ODEs), which was initially conceived to examine the problem of weather forecasting. 2 Its derivation is based on a model of Rayleigh-Bénard convection. 3,4 By considering only two-dimensional rolls, the governing partial di erential equations (PDEs) can be transformed into an ODE system via truncation of the Fourier series expansions of the stream function and temperature perturbation. 5 Taking only the rst mode of the Fourier series then yields the three-dimensional Lorenz system. 2 The Lorenz system belatedly triggered the explosion of interest in nonlinear dynamics and chaos theory since the 1980s 6 and led to some of the landmark results including the proof of the existence of the Lorenz attractor. 7 The governing equations from which the Lorenz system was originally derived describe the thermally driven convection of a uid. 5 Geophysical uids in nature are, however, rarely without impurities such as particulate matter in the atmosphere or salt particles dissolved in seawater. These particles can be incorporated into the model as a density-a ecting scalar in the governing equations, from which a ve-dimensional ODE system can be derived. 8 Like the three-dimensional Lorenz system it encompasses, the vedimensional system exhibits chaos through heteroclinic explosion. 9 Another important factor to consider in the context of geophysical uid convection is planetary rotation. Additional terms that take into account the Earth's rotation, or more generally, rotation of the model frame, can be added to the governing equations. 10 Via the same truncation method, a system of four ordinary di erential Chaos ARTICLE scitation.org/journal/cha equations can be obtained, 11,12 which is again an extension of the three-dimensional Lorenz system. Coincidentally, the fourdimensional system is identical to the Lorenz-Sten o equations originally designed to model acoustic-gravity waves in the atmosphere. 13 The Lorenz-Sten o system is found to exhibit interesting behavior such as closure of the chaotic parameter region within periodic regimes. 14 The Lorenz-Maas system is also an extension of the threedimensional Lorenz system applicable to ocean circulation. Its governing equations take into account both the e ects of solutes in the ocean water and rotation of the Earth. 15 This three-dimensional system happens to be identical to the conceptual model for the general atmospheric circulation proposed a decade earlier by Lorenz himself. 15,16 The derivation of the Lorenz-Maas system, however, takes a di erent approach from the aforementioned systems in that the simpli cation of the governing equations is done based on density gradients and angular momentum rather than truncation of Fourier series expansions. Furthermore, at various stages of the derivation, ocean-speci c assumptions were made to a x some of the parameters to the related variables, 15 resulting in fewer equations and parameters as part of the system than what is expected from a system with two additional ingredients; therefore, it still remains to be seen how these additional e ects can be incorporated into a system whose derivation closely follows the classical method of Lorenz. 2 The present study aims to derive rigorously and investigate a high-dimensional Lorenz system that includes the e ects of rotation of the model frame and density-a ecting scalar. The resulting system is a new six-dimensional nonlinear ODE system. It is both an extension of previously derived Lorenz systems and a generalization of the Lorenz-Maas system. The six parameters controlling this system include three additional parameters associated with rotation and density-a ecting scalar. To see the e ects of these additional parameters, the subsequent discussion centers around various parameter spaces.
The behavior of the three-dimensional Lorenz system in parameter spaces has already been extensively studied in the context of bifurcation and route to chaos. In acknowledgement of the Lorenz system's uid dynamical origin, early studies focused on nding critical Rayleigh parameters, beyond which the system experiences an onset of preturbulence 17 or chaos. 2,18 Similar discoveries are reported in the studies that explore the Prandtl number space (or σ -space). 19 Many studies are also devoted to the behavior of the original Lorenz system in the planes of two parameter pairs. [20][21][22] A detailed treatment of the bifurcation structure of the three-dimensional Lorenz system is found in a recent study. 23 In our paper, the new six-dimensional system will be examined through linear stability analysis, periodicity diagram, Lyapunov exponents, and numerically computed trajectories in various parameter spaces, encompassing the parameters both old and new.
Finally, this study will show that the new six-dimensional system self-synchronizes in the same way the three-dimensional Lorenz system does. Self-synchronization in the context of a nonlinear system generally refers to the merging of numerical solutions over time between the nonlinear system, the driver, and a slightly modi ed system (hence the pre x, "self-" in self-synchronization), the receiver, usually with some parts of the original equations replaced by solutions received from the driver. 6 It is well known that for self-synchronization of the three-dimensional Lorenz system to occur, only the solution for one of the variables from the driver needs to be passed onto the receiver. 24 Synchronization of the three-dimensional Lorenz system inspired various potential applications ranging from exchanging secret messages 25 to modeling major climate shifts. 26 While self-synchronization in high-dimensional Lorenz systems would be useful at least in enhancing certain aspects of the applications of this phenomenon, 27 it is not immediately clear whether a given high-dimensional Lorenz system self-synchronizes due to having additional variables nonlinearly interacting with each other. 28 For this reason, when a new system such as our new sixdimensional ODE system emerges, it is worthwhile to clarify whether self-synchronization occurs in the new system and if so, how much information from the driver is needed.

II. DERIVATION
Consider the following set of governing equations in a resting reference state under the Boussinesq approximation: In the governing equations, partial derivatives with respect to y are omitted because only two-dimensional rolls are considered. Equations (1)-(3) are the momentum equations in the x, y, and z directions, respectively, and Eq. (4) is the continuity equation. Equation (5) is the thermodynamic energy equation, and Eq. (6) is the scalar advection-di usion equation. In the above equations, u, v, and w are, respectively, the velocity components in the x, y, and z directions, and p denotes the departure from the reference pressure. The stream function ψ is de ned such that u = −∂ψ/∂z and w = ∂ψ/∂x. Note that T is the departure from the reference temperature and S the departure from the reference scalar concentration. In the above equations, is the angular velocity of the rotating frame, ρ 0 is the reference density, g is the gravitational acceleration, and ν m , ν T , and ν S are the kinematic viscosity, the thermal di usivity, and the scalar di usivity, respectively. The basic state temperature and scalar concentration decrease linearly with height. Let T = T b − T t and S = S b − S t , where T b , T t , S b , and S t are the temperature at the bottom boundary, the temperature at the top boundary, the scalar concentration at the bottom boundary, and the scalar concentration at the top boundary, respectively. The thermal expansion coe cient and the scalar contraction coe cient are denoted by α and β, respectively. In this study, we only consider the case when β > 0 so that the scalar negatively a ects buoyancy.

ARTICLE
scitation.org/journal/cha Similar to Saltzman, 5 de ne where H is the distance between the top and bottom boundaries. Analogously, de ne We impose boundary conditions ∂u/∂z = ψ = θ = ς = 0 at z = 0 and z = H. For nondimensionalization, the following substitutions are introduced: along with parameters de ned as follows: where R 1 and R 2 denote the Rayleigh numbers associated with temperature and scalar concentration gradients, respectively, Ta is a rotational parameter called the Taylor number, σ is the Prandtl number, and Le is the Lewis number. Then, the following nondimensional equations can be obtained: The superscripts * s are dropped in Eqs. (14)- (17) for notational convenience.
To transform the PDE system to an ODE system, the sinusoidal series expansions of ψ, v, θ , and ς in terms of ODE variables X, Y, Z, V, W, and U are truncated to the rst mode as follows: where k x and k z are the scales for wavenumbers in the x and z directions, respectively, such that k 2 = k 2 x + k 2 z . The coe cients in Eqs. (18)-(21) are given by (14)- (17) and applying the scaled parameters, yield the six-dimensional ODE system, where the dots indicate the derivatives with respect to the scaled nondimensional time τ = k 2 t.

III. SYSTEM PROPERTIES
Equations (23)-(28) contain six parameters and six variables including three new variables, V, W, and U, not in the threedimensional Lorenz system. For parameters r T , σ , and b, the canonical values originally used by Lorenz are 28, 10, and 8/3, respectively. 2 Parameter r C is the normalized Rayleigh number for scalar concentration, and parameter s is associated with rotation of the model frame. While these parameters have physical implications in the context of uid convection, here we focus instead on nonlinear dynamical aspects of the new system, viewing it as an extension of the three-dimensional Lorenz system.

A. Fixed points and stability
Setting the derivatives in Eqs. (23)- (28) to zero allows the nontrivial xed points to be expressed in terms of X 0 , which satis es the quartic equation, f (X 0 ) = 0, as follows: where A and B denote the coe cients of X 4 0 and X 2 0 , respectively, and C is the constant term in Eq. (29). The four roots of Eq. (29), lead to four nontrivial xed points X ±,∓ 0 as 6-tuples. For linear stability analysis, these xed points are then plugged into the linearized equations obtained from applying in nitesimal perturbations The real parts of the eigenvalues of matrix M are used to plot neutral stability curves in Fig. 1 corresponding to given parameter combinations. In Fig. 1(a), the solid curve for s = 0 going across the parameter space diagonally divides the r T -σ space into two based on the stability criterion; the system is unstable if the (r T , σ ) parameter pair belongs to the region on the right-hand side of the curve. In other words, if σ is too small or if r T is too large, the system loses its stability and exhibits unstable trajectories such as limit cycles or chaotic attractors. Figure 1(b) shows neutral stability curves in the r C -s parameter space for di erent values of r T . The system is unstable in the region under the curve and stabilizes when r C or s goes over the threshold. Raising r T seems to expand this unstable region in the parameter space.
Like the three-dimensional Lorenz system, the new system is symmetric in the sense that its solutions remain invariant under a mapping (X, Y, V, W) → (−X, −Y, −V, −W). Furthermore, each of the nontrivial xed points X ±,∓ 0 derived from the four roots of Eq. (29) has a symmetric counterpart with respect to X = 0 or Y = 0 if it exists in the real space. The existence of real-valued xed points in relation to f (X 0 ) of Eq. (29) is illustrated in Figs. 1(c)-1(f). The four nontrivial xed points produce at most two stability curves, one for each symmetric pair, and a neutral stability curve marks where the system becomes unstable around all, trivial and nontrivial, xed points. Within the bounds of the parameter spaces given in Fig. 1, however, the xed points X +,− It is clear from Fig. 1 that the stability of the system around the xed point is a ected by changes in parameter values. The diagonal neutral stability curves shift to the right with increasing s. Since the analogous destabilization of the nontrivial xed points via a Hopf bifurcation in the original Lorenz system is closely preceded by the heteroclinic bifurcation ushering the onset of chaos, 29 the observed shifts in the neutral stability curve raise the possibility that given a larger s, a higher critical r T is needed for the onset of chaos in the new system, which is con rmed to be the case in the subsequent analyses of numerical solutions.

B. Periodicity in the r T -σ space
The periodicity diagram is a useful tool in visualizing the numerically computed solutions' behavior on parameter spaces. 22,30 One can also plot Lyapunov exponent spectra to indirectly make inferences about the periodicity structure in parameter spaces. 20,31 For direct computation of periodicity at each parameter pair (r T , σ ) in the range r T , σ ∈ [0, 500] with resolutions r T = σ = 1, the fourth-order Runge-Kutta method with time resolution τ = 10 −4 is employed to compute numerical solutions with xed parameters given by r C = 28, b = 8/3, Le −1 = 0.1, and s = 1. We focus on the solutions for variable Z and consider any two solutions to be the same if the di erence between them is less than 5 × 10 −3 . If the solution is periodic, there will be a nite number of peaks or local maxima before the pattern repeats itself in its time series. The nite number of peaks counted in this way, say n, indicates that the solution is n-periodic. Since chaotic solutions would have in nite peaks with no patterns being repeated, the greater the number of peaks, the more likely the solution is chaotic. For example, Fig. 2(b) is the trajectory of numerical solutions from τ = 200 to  Fig. 2(a)], the solution is 3-periodic at (r T , σ ) = (320, 164), which translates to a three-looped trajectory [ Fig. 2(b)]. Where the solution is likely to be chaotic based on the periodicity diagram, such as at (r T , σ ) = (295, 164), the corresponding trajectory from τ = 200 to τ = 600 in Fig. 2(c) does not return to its starting point. Moreover, the trajectory even adds a ap to its side, with the solution in the X-Y space unpredictably alternating between two round aps that are slightly twisted in between about the Z-axis. This twisting might resemble the turning of the unstable manifold of the attractor for large r T in the original Lorenz system. 23 Like in the periodicity diagram for the three-dimensional system, the so-called "onionlike structure" 22 appears in our system, which refers to the alternating appearance of chaos and periodicity within the unstable region of the r T -σ space. The boundary marking the transition from stable (0-periodic) to unstable (n-periodic, n > 0) solutions in the periodicity diagram closely follows the neutral stability curve. Given what we have observed from Fig. 1, we expect the stable-unstable boundary in Fig. 2(a) to retreat northeastward if s is further increased. Compared to our new six-dimensional system [ Fig. 2(a)], the stable-unstable boundary for the original Lorenz system is steeper, thereby having enough room in the unstable region to accommodate more bands of the onionlike structure within the 500 × 500 parameter space window compared to what is visible in Fig. 2(a). This is in line with the behavior of the mathematically extended six-dimensional Lorenz system from an earlier study; 30 however, one notable feature in the periodicity diagram for that system 30 is having several discrete batches of regions with di erent period numbers lined up parallel to the σ -axis at r T ∼ 390. This feature is clearly not observed in the new physically extended six-dimensional system.
The apparent instability at low r T values visible as a thin white vertical band in Fig. 2 characterize the behavior of the new six-dimensional system when r T is small.

C. Effects of r C and s
By adjusting the parameters r C , Le, and s, the e ects of r C and s, loosely interpreted as the e ects of scalar gradient and rotation, respectively, can be compared against the well-known behavior of the three-dimensional Lorenz system. We focus on the two important points in the r T space, in particular. First, the three-dimensional Lorenz system is known to undergo the heteroclinic bifurcation 2,18 as the thermal Rayleigh parameter, r T , is raised beyond the critical number r T,c ∼ 24 with σ = 10 and b = 8/3. The existence of critical r T in our six-dimensional system can already be inferred from the visible shifts of the neutral stability curves in Fig. 1(a). Second, increasing r T > r T,c even further leads to alternating phases of chaos and periodicity 22 until the last bifurcation 18,33 before its return to stability takes place at r T = r T,d ∼ 313. Our goal in this subsection is to investigate how r C and s change these two points in our new system.
To get a fuller picture of how the newly introduced parameters a ect the system, we plot the Lyapunov exponent spectra and bifurcation diagrams on r T ∈ [0, 500] with a resolution of r T = 1 for di erent choices of r C and s parameter pairs as shown in Fig. 3. For each parameter choice, six Lyapunov exponents are computed based on the continuous Gram-Schmidt orthonormalization method. 34,35 For this, numerical solutions corresponding to 50 randomly selected initial conditions are computed using the fourthorder Runge-Kutta method from τ = 0 to τ = 400 with a time resolution of τ = 10 −3 . The initial data up to τ = 150 are truncated out before computing the Lyapunov exponents. The obtained exponents are then ordered so that λ 1 > λ 2 > · · · > λ 6 . The bifurcation diagrams in Figs. 3(b), 3(d), and 3(f) are plotted using the local maxima of the numerical solutions of Z, referred to as Z-peaks, obtained for the time period τ ∈ [200, 250] using the fourth-order Runge-Kutta method with a time resolution of τ = 10 −4 and the initial condition (X, Y, Z, V, W, U) = (1, 0, 0, 0, 0, 0). While the precise de nition of chaos still remains to be agreed upon, 36 it is frequently argued that if the solution is indeed chaotic, the largest Lyapunov exponent λ 1 is expected to be positive. 37  call a solution with two or more positive Lyapunov exponents hyperchaotic; 38 however, as is the case in the Lorenz-Sten o systems, 31,39 no such signs of hyperchaos in our six-dimensional system are found within the fairly extensive parameter ranges considered by this study. In bifurcation diagrams, the signature of chaos looks like a thick blob of concentrated points as opposed to points forming a nite number of, say, n clearly distinguished lines for n-periodic orbits.
With r C , s, and Le −1 set to 0 in Figs. 3(a) and 3(b), the partial system consisting of variables X, Y, and Z decouples from the rest, becoming identical to the three-dimensional Lorenz system and, therefore, the same r T,c and r T,d are expected. Indeed, the solution loses its stability in Fig. 3(a) as λ 1 soars above 0 at r T ∼ 24, the same critical r T found in the three-dimensional Lorenz system. The bifurcation diagram in Figs. 3(b) and 3(g) supports the location of this r T,c around 24. All else being equal, if we raise r C and s to 20 and 10, respectively, r T,c slightly increases from ∼24.0 to ∼24.2. If r C and s are further raised to 60 and 40, respectively, there is a greater shift in r T,c to ∼30 as shown in Figs. 3(e), 3(f), and 3(h). The changing trajectories in Fig. 4 con rm that the system transitions over to chaos as r T passes through the critical r T in each case. This increase in r T,c implies that the system's stability becomes more resistant to the destabilizing e ects of raising the thermal Rayleigh parameter under the in uence of bottom-heavy vertical distribution of scalar concentration in a rotating model frame. This does not, however, imply that the system generally becomes more stable under the in uence of rotation and density-a ecting scalar as will be shown below.
The other notable point r T,d in the three-dimensional Lorenz system marks the transition from 2-periodic orbits to symmetric 1periodic orbits, 30 shown by merging of two curves around r T ∼ 315 in Fig. 3(b). In Lyapunov spectra such as Fig. 3(a), one of the Lyapunov exponents drops below 0 at r T ∼ 315. When r C and s are set to nonzero values, this bifurcation occurs at a higher r T ; in Figs chaos is delayed with nonzero r C and s, the instability itself persists until a greater r T,d is reached. Therefore, it seems that the window of instability, the r T interval where the system undergoes chaotic, periodic, and intermittent phases, becomes elongated with increasing r C and s.
In order to see individual e ects of r C and s, Lyapunov exponent spectra are plotted in Figs. 5(a) and 5(b) for intervals r C , s ∈ [0, 500] with r C = s = 1 while keeping all other parameters constant. This setup is equivalent to cutting across the r C -s parameter space in Fig. 1(b) horizontally at s = 10 and vertically at r C = 20, respectively. For a detailed look, the three largest Lyapunov exponents have been recomputed with a higher resolution, r C = s = 0.1, as shown in Figs. 5(c) and 5(d). If s 40, the largest Lyapunov exponent dives below zero following a narrow window of intermittency [ Fig. 5(d)]. Beyond that threshold s, all Lyapunov exponents become negative, indicating stability. The Lyapunov spectra on the r C in Figs. 5(a) and 5(c) show a di erent picture. While there appears to be a dip from the initial chaos toward periodicity or stability around r C ∼ 121 in Fig. 5(a), the Lyapunov exponents do not immediately make a fully committed dive signi cantly below zero as seen in Fig. 5(c). In fact, there appears a second onset of chaos at r C ∼ 251 which lasts until r C ∼ 288 before settling down to a steady-state solution. For the chaotic regime that appears in r C ∈ [251, 288], there emerges an attractor, shown in Figs. 6(a)-6(c) and 7, which seems qualitatively distinct from the classic Lorenz attractor.
The trajectories calculated with ve initial conditions, one near 0 and four near the nontrivial xed points, are given in Figs. 6(a), 6(b), and 6(d). The trajectories for the rst onset of chaos in r C space resembles the Lorenz attractor as shown in Fig. 6(d); despite having di erent initial conditions, di erent trajectories quickly organize themselves into a nearly two-dimensional shape. In Figs. 6(a) and 6(b), on the other hand, the trajectories corresponding to the second onset of chaos wander about independently of each other for a much longer period of time than the trajectories for the Lorenz attractor, sometimes tracing an elaborate structure [for example, the black portion of the trajectory in Fig. 6(c)] before the wings of chaotic motions emerge [e.g., the blue portion of Fig. 6(c)]. Even after locking into the chaotic motions, the trajectories for the second chaotic regime maintain the distinctive disheveled appearance as shown in the chaotic trajectories in Fig. 7. It would, therefore, seem that self-synchronization, a process through which solutions with different initial conditions ultimately end up with precisely the same trajectories, is less e ective in the second chaotic regime, requiring procurement of more information from the driver compared to the original Lorenz system; nevertheless, our analysis in Subsection III D shows that this is not so. Figure 7 presents the trajectories of solutions from τ = 200 to τ = 1600 corresponding to di erent r C values under the same condition as Fig. 5(a) in the snapshot style to showcase in detail the appearance and disappearance of chaos during the "second onset" as r C is raised. Note that the intermittency of solutions inherent in our system means that there can be regime changes amongst chaotic, periodic, and steady-state solutions interspersed throughout the r C range. The color of the trajectories gets darker with τ so that for nonchaotic trajectories such as Figs. 7(a), 7(b), and 7(f), the portions of trajectories with dark colors indicate where the trajectories sink toward as τ → ∞.

D. Self-synchronization
In this subsection, we now present self-synchronization properties of the six-dimensional system. The self-synchronization in the three-dimensional Lorenz system requires one-third of the information from the driver, namely, the driver's solution to variable X. In a nonlinear system with increased complexity, it is reasonable to think that more information is needed from the driver for the receiver to replicate the driver's solution. For example, it was shown that for high-dimensional systems that exhibit hyperchaos, self-synchronization is attained by transmitting a linear combination of the variable signals. 27 For a de nite proof of self-synchronization in a given system, a Lyapunov function of the error system can be found, 40 so that the error system is globally asymptotically stable at 0. If analytically constructing an appropriate Lyapunov function is not readily available, as it seems to be the case for some high-dimensional Lorenz systems, 30 numerical experiments can be performed for demonstration. Both approaches are explored here.
Given our new six-dimensional system as the driver, the receiver system that only receives X from the driver is is a Lyapunov function so that the error e approaches 0 as τ → ∞. Therefore, the receiver's solution catches up with the driver's, regardless of the respective initial conditions. The self-synchronization of the six-dimensional system is numerically tested in Fig. 8 with parameters that produce the attractor illustrated in Fig. 7(d): r T = 28, b = 8/3, σ = 10, r C = 280, s = 10, and Le −1 = 0.1. The numerical solutions of the driver is computed using the fourth-order Runge-Kutta method with τ = 10 −4 and the initial condition (X, Y, Z, V, W, U) = (1, 0, 0, 0, 0, 0). At τ = 200, when the system comfortably exhibits chaos after some phases of transient behavior, the receiver is deployed to start out with the initial condition ( X, Y, Z, V, W, U) = (50, 50, 50, 50, 50, 50). Figure 8(a) juxtaposes the time series of Z from the driver (black) with the time series of Z from the receiver (red). The receiver's solution quickly catches up with the driver's, enough for the errors to go below 2 × 10 −2 at τ = 205 in the time series of the absolute di erence between the two in Fig. 8(b). Figure 8(c) shows how the trajectory of the receiver in red merges with the trajectory of the driver in black, which, when left evolving, produces the attractor shown in Fig. 7(d).

IV. CONCLUSION
In this study, we derived a new six-dimensional ODE system by adding the density-a ecting scalar and rotational e ects in the equations governing thermal convection and analyzed the system. Although the Lorenz systems are not precise models of real uid convection, some conceptual insights can still be gained from studying our new ODE system. For example, it may seem that the inclusion of scalar and rotational e ects stabilizes the system based on the greater critical Rayleigh parameters; however, as evidenced by Fig. 3, the increase in r T,c is also accompanied by elongated phase of chaos, periodicity, and intermittency in the r T space. For this reason, it would not be surprising if one nds partial connections between newly Chaos ARTICLE scitation.org/journal/cha observed behavior of this system and some unexplored uid phenomena observed in the Earth's atmosphere and oceans or laboratory experiments. It is also of mathematical interest to explore di erent kinds of attractors exhibited by this system. The mathematical extensions of this system to higher-dimensional systems can also be investigated, which can be a source of new discoveries as precedented by previous attempts to generalize the original Lorenz system 28,30,41 and the Lorenz-Sten o system. 31,42