On periodic solutions in the non-dissipative Lorenz model: the role of the nonlinear feedback loop

The nonlinear feedback loop and energy conservation laws

2.1.

Nonlinear terms in the 3D-NLM (and 3DLM) have been shown to result from the Jacobian term, $J(\psi ,\theta )$, in Equation (2). The nonlinear interaction of two wave modes via the Jacobian term can generate or impact a third wave mode through a downscale (or upscale) transfer process. The subsequent upscale (or downscale) transfer process associated with the third wave mode can provide feedback to the incipient wave mode(s). As illustrated in Shen (2014), XY and –XZ, respectively, represent the downscale and upscale transfer processes that form a nonlinear feedback loop. When new modes are properly included, the feedback loop can be extended. In the following subsections, we discuss the role of the 3D-NLM nonlinear feedback loop in the energy conservation and partition of available potential energy, which, in turn, helps produce oscillatory solutions.

The domain-averaged kinetic energy ($\overline{\text{KE}}$), the potential energy ($\overline{\text{PE}}$), and the available potential energy ($\overline{\text{APE}}$) are defined as follows (Thiffeault and Horton, 1996; Blender and Lucarini, 2013; Shen, 2014, 2015):

Through straightforward derivations, we obtain the following:

Equations (10) and (11) lead to:

With Equations (3)–(5), the time derivative for both Equations (13) and (14) is zero, so these two equations describe energy conservation. Both C₁ and C₂ are constants and are determined by the initial conditions. Thus, if we express Z and Y² as functions of X, they are single valued. To facilitate our discussions, the contribution to APE from an individual mode is defined as $\overline{{\text{APE}}_I}=-C_o\sigma I^2/2r$, where I = Y or Z; therefore, $\overline{\text{APE}}=\overline{{\text{APE}}_Y}+\overline{{\text{APE}}_Z}$. Note that Equations (13) and (14) are related to the two Nambu Hamiltonians, $C=-X^2/2+\sigma Z$ and $H=Y^2/2+Z^2/2-rZ$ (Nambu, 1973; Nevir and Blender, 1994; Floratos, 2012; Roupas, 2012; Blender and Lucarini, 2013).

From the initial conditions in Equation (6), we have $C_1=0$ and $C_2/C_o=-\sigma /2r$, the latter of which is –0.2 for r = 25 and –0.11 for r = 45. Figure 1 provides the time evolution of the conserved quantities: ($\overline{\text{KE}}+\overline{\text{PE}}$) and ($\overline{\text{KE}}+\overline{\text{APE}}$) in Equations (13) and (14) from the 3D-NLM. At a larger r (e.g. r = 45), a finer $△\tau$ is required to improve the numerical accuracy of simulated total energy (Fig. 1c). In this study, to facilitate discussions and unless stated otherwise, either C₁ or C₂ is assumed to be zero.

Fig. 1.

Solutions of the linear system without the nonlinear feedback loop

3.1.

By assuming no nonlinear terms, we can begin with two equations: $\text{d}X/\text{d}\tau =\sigma Y$ and $\text{d}Y/\text{d}\tau =rX$ and only one conservation law, as follows:

This linear case yields $M^2=-\sigma r$, and Equation (15) becomes ${\text{d}}^{2}X/\text{d}{\tau }^2-\sigma rX=0.$ A local analysis suggests that the origin, the only critical point in the linear system, is a saddle point. After straightforward derivations, the solution is:

Nonlinear solutions and the nonlinear restoring force for r = 0

3.2.

Here, we examine the impact of nonlinear terms using a special case with r = 0 and $(X,Y,Z)=(0,1,0)$, leading to $M^2=X^2/2$. Thus, Equation (15) and its solution are written as follows:

In Appendix B, we show that the cubic term (X³) in Equation (18a) is derived from the nonlinear feedback loop (–XZ and XY). The solution in Equations (18b) and (18c) is periodic and it has time varying frequencies. More importantly, as discussed later, when the approach is applied to the system with $r\ne 0$, we can also obtain periodic and bounded solutions, which is different from the unbounded solution of the linear system in Equation (17). Thus, we suggest that the nonlinear feedback loop (i.e. the cubic term) acts as a restoring force. Solutions for other components (Y and Z) can be obtained as follows.

As compared to the case with $r\ne 0$ in Equations (13) and (14), the energy conservation laws obtained using r = 0 are:

Since X² is defined in Equation (18b), solutions for Y and Z can be obtained as follows:

To illustrate the solution’s characteristics, Equations (22a) and (23) are solved using the following iterative method:

Fig. 2.

To verify the integral form of the solutions in Equations (24a) and (24b), the numerical solutions of the 3D-NLM using r = 0 (Equations (3)–(5)) are provided in Fig. 3. In panel (a), the blue ‘dot’ indicates the initial temporal evolution of the phase function that is calculated by performing a time integration of X using Equation (18c), where X is obtained from the 3D-NLM; and the red line indicates the phase function calculated using Equations (24a) and (24b). Both are in good agreement. The simulated trajectories in the X–Y and X–Z sections are elliptic and parabolic (Figs. 3b and 3c), respectively, consistent with the analytical relationships in Equations (21) and (19), respectively. Figure 3d provides the time evolution of oscillatory Y (red) and X (blue), consistent with the results provided in Figs. 2b and 2c, respectively. As indicated in Appendix C, the oscillatory characteristics of the solution can also be illustrated using elliptic functions.

Fig. 3.

Closed form solutions near a non-trivial critical point for $r\ne 0$

3.3.

In the previous sections, we used Equation (15) to illustrate the individual impact of the linear (heating) force and the nonlinear feedback loop using cases with M² < 0 and $M^2\ge 0$, respectively, neither of which changes the sign during the integration. Here, using the initial condition $(X,Y,Z)=(X_o,0,0)$ and $0\hspace{0.17em}<\hspace{0.17em}{X}_o$, we consider a more general case with $M^2=(X^2/2-\sigma r-X_o^2/2)$, whose sign may vary during the time integration depending on the relative magnitude of the nonlinearity and the linear (heating) force. M² has two zeros at $X=\pm X_{\text{t}}$, where $X_{\text{t}}$ is defined as $\sqrt{2\sigma r+X_o^2}$. These are called turning points. Based on an analysis using the WKB approximation (Bender and Orszag, 1978), there appears to be a growing or decaying solution for $\left|X\right|\hspace{0.17em}<\hspace{0.17em}{X}_{\text{t}}$ and an oscillatory (wave-like) solution for $\left|X\right|>X_{\text{t}}$. The former is largely impacted by the linear (heating) force while the latter is impacted by the nonlinear restoring force. Additional analyses are provided below.

To understand the stability of solutions, below, we present a local analysis by linearizing the system with respect to a non-trivial critical point. Using Equations (3)–(5), we first solve for the non-trivial critical point, which is $(X_c,Y_c,Z_c)=(X_c,0,r)$. Here, X_c can be any constant and is selected as the value of X at the turning point of Equation (15) using the justification presented below. From Equations (3) and (15), the 3D-NLM system can be written as:

The two equations, with an initial condition of (X, Y) =  $(X_0,0)$, lead to three critical points, (X_c, Y_c) at (0, 0) and $(\pm \sqrt{2\sigma r+X_0^2}\mathrm{,} 0)$. Note that the non-trivial critical points $\pm X_c$ are determined by the condition M = 0, which defines the turning points in Equation (15). Therefore, we have X_c = X_t.

Next, each of the total fields is decomposed into its basic part and perturbation, as follows:

The above system with Equations (28)–(29) yields eigenvalues of 0 and $\pm i{X}_c$, suggesting that the non-trivial critical point ($X_t,0,r$) for the linearized system is a stable node as a center. The other non-trivial critical point (−$X_t,0,r$) shares the same features. One important feature is the dependence of the solution’s period ($2\pi /X_c$) on the initial condition (X_o). Next, to illustrate its periodicity, we present a closed form solution for the nonlinear system.

Equations (13)–(14), with the initial conditions (X, Y, Z) = ($X_o,0,0$), lead to the following two equations:

Using the same procedures as those provided in Section 3.2 (for r = 0), the following closed-form solutions (for $r\ne 0)$ are obtained:

Equations (32c) and (32d) are iteratively computed in order to determine the value of $\varphi$ that is used to calculate the solutions of X, Y, and Z. Note that because $X_t^2=2\sigma r+X_o^2$, Equations (32a) and (32c) lead to:

Since $r\ge rcos(\varphi )$, X² is always positive in Equation (32c). Without a loss of generality, $X=+\sqrt{2\sigma (r-rcos(\varphi ))+X_o^2}$ is first discussed, yielding a minimum of $X_{min}=X_o$ and a maximum of $X_{max}=\sqrt{4\sigma r+X_o^2}$. In a similar manner, given an initial condition of $(X,Y,Z)=(-X_o,0,0)$, X has a maximum of $-X_o$ and a minimum of $-\sqrt{4\sigma r+X_o^2}$. Note that the average of $X_{min}^2$ and $X_{max}^2$ is equal to $2\sigma r+X_o^2$, which is equal to $X_t^2$.

Applying the initial conditions $(X,Y,Z)=(\sqrt{2\sigma r},0,0)$ that yield $X_c=X_t=2\sqrt{\sigma r}$ ($\approx 31.62$), Fig. 4 shows the closed form solutions obtained using Equation (32), as well as numerical solutions of the 3D-NLM obtained using Equations (3)–(5). The results are in good agreement. For any positive X_o, a trajectory beginning with $(X,Y,Z)=(X_o,0,0)$ is a closed curve. An oscillatory (closed) trajectory associated with a center is orbitally stable (Jordan and Smith, 2007). In Section 3.4, the time varying energy cycle in Equation (32) will be analyzed. Next, we discuss the solutions in (32) for the special case of X_o = 0.

Fig. 4.

When considering a homoclinic solution that begins with and ends at the saddle point, initial conditions of $(X,Y,Z)=(0,0,0)$ are required. At the saddle point, the time derivatives of the 3D-NLM in Equations (3)–(5) are zero. Therefore, numerically, only zero solutions for X, Y, and Z are possible. Obtaining a non-trivial solution can be attempted by initially adding a small perturbation (i.e. $(X,Y,Z)=(\mathrm{ϵ},0,0)$) and then performing numerical integrations. A similar approach can be applied in Equation (32) to obtain an approximate solution. Here, we provide an initial guess for $\varphi ,\varphi (\tau )=\tau$, over a target period in Equation (32d). We then calculate X using Equation (32c) and $\varphi$ using Equation (32d) through iterations. Figure 5 displays the solutions determined using Equation (32). The time evolution of Y is provided in Fig. 5a and suggests that Y increases slowly with time, reaches a maximum at $\tau \hspace{0.17em}\approx \hspace{0.17em}2.3$, and then decreases with time to its minimum. The X–Y plot in Fig. 5b appears reasonable and is later compared to the analytical solution. In Fig. 5a, the evolution seems ‘periodic’ but displays some irregularities that are different from those obtained using $X_o\ne 0$. This case is further analyzed using the following analytical solution, as well as the numerical solutions in Section 3.5.

Fig. 5.

From Equation (D3) of Appendix D, the ‘second’ half of the homoclinic orbital solution (e.g. the lower part of the orbit in Fig. 5b) can be obtained as follows:

When $t\to \infty$, the above solution approaches the origin. Note that a special point on the solution trajectory is $(X,Y,Z)=(2\sqrt{\sigma r},0,2r)$. Therefore, the ‘second’ half of the solution in Equation (34) describes the trajectory that begins at $(X,Y,Z)=(2\sqrt{\sigma r},0,2r)$ and that moves forward in time toward the origin, as shown in light blue in Fig. 5b. The special point provides an alternative IC for the numerical integration of the 3D-NLM required to simulate (qualitatively) the homoclinic orbit, discussed with numerical results later. When $t\to \infty$, $X(\tau )\hspace{0.17em}\approx \hspace{0.17em}4\sqrt{\sigma r}{e}^{-\sqrt{\sigma r}\tau }$ and $Y(\tau )\hspace{0.17em}\approx \hspace{0.17em}-4r{e}^{-\sqrt{\sigma r}\tau }$, yielding $X/Y=-\sqrt{\sigma /r}$. These features are the same as those for the linear (stable) solution with a heating term, as outlined in Section 3.1. As the system is invariant under $t\to -t$ and $y\to -y$ (Strogatz, 2015), the solution ($X(\tau ),-Y(\tau ),Z(\tau ))$ in backward time, which represents the ‘first’ half of the solution, displays the trajectory that begins at the origin and then moves to the point $(X,Y,Z)=(2\sqrt{\sigma r},0,2r)$ (e.g. the upper part of the orbit in Fig. 5b). The total solution that consists of the first and second half connects the stable and unstable manifolds of the saddle point and is called the homoclinic orbit. Since the trajectory ‘takes forever’ to reach the saddle point (as $t\to \infty$), it is, therefore, not a periodic solution. As a result, the ‘periodicity’ of closed-form solutions with X_o = 0, as provided in Fig. 5, is spurious, suggesting a numerical error. Unique features of the homoclinic orbit within the analytical solution, as well as the closed-form solution, will be further explored using numerical solutions below.

An energy cycle of solutions with $r\ne 0$

3.4.

Here, we analyze the time evolution of energy along the homoclinic orbit and apply it to nearby trajectories with a small energy cycle. In Equation (D3), $Y^2\ge 0$ leads to $\left|X\right|\le 2\sqrt{\sigma r}$, which yields a maximum of X ($X_{\text{m}}=2\sqrt{\sigma r}$) at Y = 0, denoted by Y_m = 0. Taking the partial derivative of Equation (D3) with respect to X indicates that Y has extrema when $X=\pm \sqrt{2\sigma r}=\pm X_{\text{t}}$. At $X=\pm X_{\text{t}}$, Z = r, and $Y=\pm r$ from Equations (13) and (14), respectively, and, thus, $\overline{{\text{APE}}_Y}$ = $\overline{{\text{APE}}_Z}$, in other words equal contributions to $\overline{\text{APE}}$ from the Y and Z modes. Furthermore, Equation (D3) suggests that Y² initially increases in association with the increase in X² but later decreases in association with an increase in X⁴. The former is consistent with the linear case in Equation (16), while the latter is consistent with the simplified nonlinear case in Equations (19) and (20). The distribution of Y as a function of X Equation ((D3)), with the aforementioned four points, is provided in Fig. 6. The energy cycle: (1) begins at point A; (2) goes through B, C, and D; and (3) returns back to A. The segment from point P to point Q is designated as Leg P–Q, where P and Q are either one of the following: $A(X,Y)=(0,0)$, $B=(X_{\text{t}},Y_{\text{t}}),C=(X_{\text{m}},Y_{\text{m}})$, or $D=(X_{\text{t}},-Y_{\text{t}})$. Point A is a saddle point as discussed in Section 3.1, $(X_{\text{t}},Y_{\text{t}})=(\sqrt{2\sigma r},r)$, and $(X_{\text{m}},Y_{\text{m}})=(2\sqrt{\sigma r},0)$. The energy cycle is referred to as a small (energy) cycle. From a perspective of potential energy partitioning, the linear force dominates in Legs A–B and D–A where $\left|\mathrm{A}\mathrm{P}{\mathrm{E}}_Y\right|\ge \left|\mathrm{A}\mathrm{P}{\mathrm{E}}_Z\right|$, while the nonlinear restoring force dominates, in Legs B–C and C–D where $\left|\mathrm{A}\mathrm{P}{\mathrm{E}}_Y\right|\le \left|\mathrm{A}\mathrm{P}{\mathrm{E}}_Z\right|$. This also illustrates the importance of the inclusion of the Z mode and thus the nonlinear terms. Here, it should be noted that the homoclinic orbit is not periodic, and that a small energy cycle may appear within the oscillatory trajectories that move close to, but do not pass through, the saddle point. In the discussion provided in the next section, the analysis discussed above is inter-compared to the numerical solutions obtained from the 3D-NLM.

Fig. 6.

Numerical solutions with $r\ne 0$

3.5.

In this section, we first discuss the numerical solutions with small cycles that were analyzed in the previous section, and then present the numerical solutions with large cycles that enclose three critical points. The former are obtained using ICs of $(X,Y,Z)=(\pm 1,0,0)$, while the latter are simulated with ICs of $(X,Y,Z)=(0,\pm 1,0)$. Since the features of solutions (e.g. homoclinic orbits) require a high level of numerical accuracy, for Fig. 7, we perform numerical integrations using a very small time step of $△\tau ={10}^{-5}$ but for a shorter period of time (τ = 4). Figure 7a displays oscillatory solutions with small cycles, each of which moves around one of the non-trivial critical points. Figure 7b displays the two oscillatory solutions with a large cycle, moving around all three critical points. The two solutions share the same orbit. The numerical result of the ‘homoclinic orbit’ is shown in Fig. 7c. Figure 7d shows the evolution of these solutions. While the evolution of energy described by the small or large (energy) cycle will be analyzed later, we first discuss the results of ‘homoclinic orbit’.

Fig. 7.

The special features of the ‘homoclinic orbit’ are illustrated using numerical simulations with ICs of $(X,Y,Z)=(2\sqrt{\sigma r},0,2r)$. Here, the ‘second’ half of the orbital solution with X > 0 and Y < 0 is a focus. As indicated by the green line in Figs. 7c and 7d, the solution Y decreases and then increases along the homoclinic orbit. Nearly zero values of Y in Fig. 7d suggest that the homoclinic trajectory approaches the saddle point and ‘remains’ close to the saddle point for a long period of time, as indicated by a comparison between the homoclinic orbit solution (in green) and the oscillatory solutions (in light blue and red). Between τ = 0 and $\tau =1.7$ in Fig. 7d, the numerical solution in green represents the ‘second’ half of the homoclinic orbit. In Fig. 7c, the corresponding trajectory in the X–Y phase space that appears in the region with positive values of X compares well with the analytical solution in Fig. 5b. Near the saddle point, (X, Y, Z) are near zero. Therefore, it is easy for numerical errors (e.g. associated with rounding or truncation) to introduce non-zero forcing to drive the flow. When this occurs, the trajectory moves into a region that has negative values of X, and moves around a non-trivial critical point, as shown in Fig. 7c. This move is inconsistent with the analytical solution. However, it is challenging to avoid this type of error since $△\tau$ is already very small, ${10}^{-5}$. Avoiding this type of error may be more complicated in a higher-dimensional space (e.g. the 5DLM), which is beyond the scope of this study.

The above discussion suggests that the homoclinic orbit separates some trajectories that‘encircle’ a non-trivial critical point from other trajectories that move around the three critical points. Therefore, the homoclinic trajectory that connects the stable and unstable manifolds of the saddle point may be viewed as a separatrix. Two different points that are initially near but not on the homoclinic path will be on periodic orbits that depart from one another, suggesting diverged trajectories. Next, we use the results of Fig. 7a,b as an example to illustrate the divergence of two trajectories starting from (ϵ, 0, 0) and (0, −ϵ, 0). Here, ϵ is a small number. Applying initial conditions $(X,Y,Z)=(1,0,0)$ and $(0,-1,0)$, two trajectories shown with red and light blue lines display ‘repeated’ divergence and convergence. (The two trajectories may give rise to ‘periodic’ divergence and convergence if the ratio of their periods is a rational number.) Therefore, while each of the oscillatory trajectories is orbitally stable, the homoclinic orbit is not. Note that this type of solution dependence on ICs is different from that of a chaotic attractor in dissipative systems (e.g. the 3DLM). Next, we analyze the energy cycle of the oscillatory solutions.

As shown in Fig. 8a using the selected ICs, the large cycle includes two ‘'small cycles’ that are anti-symmetric with respect to the saddle point A. The large cycle resembles the wing of a glasswinged butterfly. The right-hand side of the wing displays the same (or comparable) characteristics as those of the cycle in Fig. 6, while the left-hand side wing is anti-symmetric with respect to the origin (i.e. the saddle point). While the trajectories in X–Y and X–Z are shown in Figs. 8a and 8b, respectively, the distribution of $\overline{{\text{APE}}_Y}$ and $\overline{{\text{APE}}_Z}$ as a function of X is provided in Fig. 8c. Here, $\overline{{\text{APE}}_Y}$ and $\overline{{\text{APE}}_Z}$ represent averaged available potential energy from the Y and Z modes, respectively. From the perspective of the APE partition, the $\overline{{\text{APE}}_Y}$ (red) dominates in Leg A–B (D–A) while $\overline{{\text{APE}}_Z}$ dominates in Leg B–C (C–D) where nonlinearity is stronger. Alternatively, $\overline{\text{KE}}$ is predominantly converted from $\overline{{\text{APE}}_Y}$ when X² is small, and is largely converted from $\overline{{\text{APE}}_Z}$ when X² is large. Therefore, inclusion of the Z mode can lead to the oscillatory solution by enabling the partition of APE on different scales at different stages (i.e. linear and nonlinear stages). More detailed analysis of the energy cycle is provided later. Here, we briefly discuss the impact of different initial conditions.

Fig. 8.

In Fig. 8a, a solution may move ‘around’ the saddle point and the two stable critical points, and its trajectory defines the large cycle. With a different initial condition, a solution may go around one of the non-trivial stable points ($\pm \sqrt{2\sigma r+X_0^2}$, 0), and its trajectory forms a small cycle. For example, Fig. 8d displays results using three different ICs including X₀ = $\sqrt{2\sigma r}$, $\sqrt{2\sigma r}$ + 15, and $\sqrt{2\sigma r}$ + 30 as shown with black, blue, and red lines, respectively. Note that the numerical solution with the first IC was previously compared with the closed-form solution in Fig. 4. The three vertical green lines indicate the locations of the corresponding stable critical points at X_c = $\sqrt{2\sigma r+X_0^2}$. Solutions in Figs. 8a and 8d suggest that the larger the value of (initial) X_o is, the further the trajectory may move away from the saddle point. When an initial X_o is distant from zero (i.e. the saddle point), nonlinearity may become dominant as compared to the (linear) heating term. The corresponding solutions (e.g. the blue curve in Fig. 8d) become more symmetric with respect to the green line that passes through the nontrivial stable critical point.

Figure 9 provides the time evolution of $\overline{{\text{APE}}_Y}$ and $\overline{{\text{APE}}_Z}$ for the solutions in Figs. 8a–8c. The energy cycle: (1) begins near point A, at A⁺ = (0,0⁺) to be specific; (2) goes through B, C, and D; and 3) returns back to A, (i.e., A⁻ = (0,0⁻)). In Leg A–B, where the linear (heating) force dominates, the solution X grows gradually with {$\overline{\text{KE}}\uparrow ,\overline{{\text{APE}}_Y}\downarrow$, $\overline{{\text{APE}}_Z}\downarrow$} and $\left|\overline{{\text{APE}}_Y}\right|>\left|\overline{{\text{APE}}_Z}\right|$. The upper arrow (↑) and the down arrow (↓) are, respectively, used to indicate an increase and decrease. In Leg B–C (or C–D), where the nonlinear restoring force becomes dominant, the solution X increases (or decreases) rapidly with {$\overline{\text{KE}}\uparrow$, $\overline{{\text{APE}}_Y}\uparrow ,\hspace{0.17em}\overline{{\text{APE}}_Z}\downarrow$} (or {$\overline{\text{KE}}\downarrow$, $\overline{{\text{APE}}_Y}\downarrow ,\hspace{0.17em}\overline{{\text{APE}}_Z}\uparrow$}) and $\left|\overline{{\text{APE}}_Y}\right|<\left|\overline{{\text{APE}}_Z}\right|$. In Leg D–A, solution X decays slowly with {$\overline{\text{KE}}\downarrow ,\hspace{0.17em}\overline{{\text{APE}}_Y}\uparrow$, $\overline{{\text{APE}}_Z}\uparrow$} and $\left|\overline{{\text{APE}}_Y}\right|>\left|\overline{{\text{APE}}_Z}\right|$. Once the trajectory returns back to point A^–, it may experience another small cycle, going to B^–, C^–, and D^– and returning back to point A⁺. Here, $B^{-}=(-X_{\text{t}},-Y_{\text{t}}),C^{-}=(-X_{\text{m}},Y_{\text{m}})$, and $D^{-}=(-X_{\text{t}},Y_{\text{t}})$. The two ‘small cycles’ form the large cycle, and the time evolution of energy is the same for both wings.

Fig. 9.