Finite Cascades of Pitchfork Bifurcations and Multistability in Generalized Lorenz-96 Models

Abstract

In this paper, we study a family of dynamical systems with circulant symmetry, which are obtained from the Lorenz-96 model by modifying its nonlinear terms. For each member of this family, the dimension n can be arbitrarily chosen and a forcing parameter F acts as a bifurcation parameter. The primary focus in this paper is on the occurrence of finite cascades of pitchfork bifurcations, where the length of such a cascade depends on the divisibility properties of the dimension n. A particularly intriguing aspect of this phenomenon is that the parameter values F of the pitchfork bifurcations seem to satisfy the Feigenbaum scaling law. Further bifurcations can lead to the coexistence of periodic or chaotic attractors. We also describe scenarios in which the number of coexisting attractors can be reduced through collisions with an equilibrium.

1. Introduction

The Lorenz-96 model, which was constructed by E.N. Lorenz [1], is a frequently used toy model in studies that are related to predictability and weather forecasting. The so-called monoscale version of the model is given by the equations

$${\dot{x}}_j=x_{j-1}(x_{j+1}-x_{j-2})-x_j+F,j\in \mathbb{Z},$$

(1)

where the indices of the variables $x_j$ are taken modulo a fixed integer $n\ge 1$:

$$x_{j+n}=x_j\mathrm{for}\mathrm{all}j\in \mathbb{Z}.$$

(2)

Equation (2) can be interpreted as a periodic boundary condition. In fact, Lorenz [1] interpreted the variables $x_j$ as values of some atmospheric quantity in n equispaced sectors of a latitude circle, where the index j plays the role of longitude. The dimension $n\in \mathbb{N}$ and forcing parameter $F\in \mathbb{R}$ are free parameters.

Low-dimensional atmospheric models that were used in earlier predictability studies, such as the Lorenz-63 and Lorenz-84 models, were derived as Galerkin projections of partial differential equations describing the laws of physics [2,3,4,5]. In contrast, the Lorenz-96 model was not constructed as a physically realistic model, but rather as a model that is easy to use in numerical experiments. Nevertheless, the model has physically relevant components, such as advection terms, damping terms, and external forcing. The feature that the dimension n of the model can be chosen arbitrarily large allows for much richer dynamics in comparison to the aforementioned Lorenz-63 and Lorenz-84 models. The latter property has made the Lorenz-96 model an attractive model for studies on forecasting [6,7,8,9], predictability [10,11,12], high-dimensional chaos [13,14,15,16], and data assimilation schemes [17,18,19,20].

The broad interest in the Lorenz-96 model has inspired several authors to introduce and study modifications of the model. Kerin and Engler [21] identified the desired properties for the nonlinear advection terms and provided a classification of all advection terms that are quadratic, energy-preserving, and equivariant with respect to circulant permutations of coordinates, and localized up to some degree. Vissio and Lucarini [22] supplement the Lorenz-96 model with temperature-like variables. This addition allows for the existence of an energy cycle, in which conversion between kinetic and potential energy is possible.

In this paper, we introduce our own variant of the Lorenz-96 model. Instead of deriving our modification from physical considerations, we stay closer to the structure of the nonlinear terms in Equation (1). Specifically, we modify the nonlinear terms by simply changing indices: for a given triple $(\alpha ,\beta ,\gamma )\in {\mathbb{Z}}^3$, we consider the system

$${\dot{x}}_j=x_{j+\alpha }(x_{j+\beta }-x_{j+\gamma })-x_j+F,j\in \mathbb{Z},$$

(3)

which is again subject to the condition in Equation (2). Note that the boundary condition in Equation (2) implies that the numbers $\alpha$, $\beta$, and $\gamma$ can always be taken modulo the dimension n, whenever this is necessary. In the remainder of this paper, these systems will be identified by the symbol ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$. Note that ${\mathcal{L}}_{-1,1,-2}\left(n\right)$ is just the original Lorenz-96 model of Equation (1). The dynamics of the latter model has been studied in detail in the papers [23,24,25]. A particular phenomenon that was discovered was the succession of pitchfork bifurcations, which, in combination with further bifurcations, leads to the coexistence of attractors. The main purpose of this paper is to determine what extent this scenario occurs in the systems ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$. Rather than providing an exhaustive analysis of all possible cases, we will highlight the differences and similarities for three concrete choices of $(\alpha ,\beta ,\gamma )$.

This paper is structured, as follows. In Section 2, we first discuss some general properties of the system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$, such as the boundedness of orbits and stability properties. In Section 3, we first derive sufficient conditions under which a pitchfork bifurcation occurs. Next, we numerically investigate whether the first pitchfork bifurcation is followed by a cascade of such bifurcations. Section 4 shows how pitchfork bifurcations in conjunction with other bifurcations lead to the coexistence of periodic or chaotic attractors. Section 5 concludes the paper with a discussion of the open questions that arise from our results.

2. General Properties

In this section, we study the properties that all members of the family ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ have in common, such as boundedness of orbits and the stability of equilibria.

2.1. Boundedness of Orbits

In this section, we determine sufficient conditions under which orbits of the system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ remain bounded for all time. This is a necessary first step, since, unlike the generalized Lorenz-96 models that were introduced in [21], the quadratic advection terms in our model do not necessarily preserve the total energy. Although unbounded orbits are not of physical relevance, they can still appear in physically relevant models, such as those that are derived from the shallow water equations, see [26,27].

The discussion in this section closely follows the arguments that are presented in [26]. Let

$${\mathbb{S}}^n=\{x=(x_0,\cdots ,x_{n-1})\in {\mathbb{R}}^n:x_0^2+\cdots +x_{n-1}^2=1\}$$

denote the unit sphere in ${\mathbb{R}}^n$ and define the following quantities:

$$\begin{array}{cc}\hfill A_n& :=\underset{x\in {\mathbb{S}}^n}{max}\sum _{j=0}^{n-1}{x}_j{x}_{j+\alpha }(x_{j+\beta }-x_{j+\gamma }),\hfill \\ \hfill C_n& :=\underset{x\in {\mathbb{S}}^n}{max}F\sum _{j=0}^{n-1}{x}_j,\hfill \end{array}$$

where we take the indices of the variables $x_j$ modulo n. Note that $A_n,C_n\ge 0$, since both of the quantities are obtained by maximizing a polynomial of odd degree. If we define $R^2={\sum }_{i=0}^{n-1}{x}_i^2$, then it follows that

$$R\frac{dR}{dt}=\frac{1}{2}\frac{d}{dt}\left(R^2\right)=\sum _{j=0}^{n-1}{x}_j{\dot{x}}_j=\sum _{j=0}^{n-1}{x}_j{x}_{j+\alpha }(x_{j+\beta }-x_{j+\gamma })-\sum _{j=0}^{n-1}{x}_j^2+\sum _{j=0}^{n-1}F{x}_j,$$

which implies that

$$\frac{dR}{dt}\le A_n{R}^2-R+C_n.$$

Now, assume that $1-4A_n{C}_n>0$, or, equivalently,

$$A_n{C}_n<\frac{1}{4}.$$

(4)

If $A_n=0$, then this condition trivially holds and $dR/dt<0$ whenever $R>C_n$, which means that the orbits of the system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ are bounded for all $t>0$. If $A_n>0$, then $dR/dt<0$ whenever $R_{-}<R<R_{+}$, where

$$R_{\pm }=\frac{1\pm \sqrt{1-4A_n{C}_n}}{2A_n}.$$

Orbits for which $R\left(t_1\right)<R_{+}$ for some $t_1\ge 0$ satisfy $R\left(t\right)<R_{+}$ for all $t\ge t_1$, which gives a sufficient condition for the boundedness of orbits. Orbits for which $R\left(t_1\right)\ge R_{+}$ for some $t_1\ge 0$ are potentially unbounded. Hence, we are mainly interested in the region $R<R_{+}$.

The next question is for which the values of n and F the condition in Equation (4) holds. For the original Lorenz-96 model ${\mathcal{L}}_{-1,1,-2}\left(n\right)$, it can be easily verified that $A_n=0$ for all $n\in \mathbb{N}$, which, in particular, implies that all of the orbits are bounded. Another example is given by the system ${\mathcal{L}}_{-1,-2,1}\left(n\right)$. More generally, if the system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ has the property that $A_n=0$ for all $n\in \mathbb{N}$, then so does the system ${\mathcal{L}}_{\alpha ,\gamma ,\beta }\left(n\right)$. However, for some systems, it may occur that $A_n>0$ for some $n\in \mathbb{N}$. A concrete example is the system ${\mathcal{L}}_{-1,0,-2}\left(n\right)$ for which $A_n=0$ when $n=4$, but $A_n>0$ when $n=5$. In such cases, it is useful to know how $A_n$ and $C_n$ vary with n in order to check the condition in Equation (4).

The equality case of the Cauchy–Schwarz inequality implies that $C_n=\sqrt{n}\left|F\right|$. However, an exact expression for the quantity $A_n$ is not so easy to derive analytically: for example, when applying the method of Lagrange multipliers, a system of quadratic equations needs to be solved as a result of the cubic terms in the expression for $A_n$. Therefore, we perform a numerical experiment in order to determine the dependence of $A_n$ on n. For a fixed choice of $(\alpha ,\beta ,\gamma )$ and n, the quantity $A_n$ is approximated by taking the largest value of the sum

$$\sum _{j=0}^{n-1}{x}_j{x}_{j+\alpha }(x_{j+\beta }-x_{j+\gamma })$$

over ${10}^5$ randomly chosen points $(x_0,\cdots ,x_{n-1})\in {\mathbb{S}}^n$. These maxima are plotted as a function of n in Figure 1; increasing the number of points in the maximization procedure only produces very minor differences. The figure clearly suggests, for several choices of the triple $(\alpha ,\beta ,\gamma )\in {\mathbb{Z}}^3$, that $A_n=O\left(n^p\right)$ with $p\approx -1$, which implies that $A_n{C}_n=O\left(\right|F|n^{p+1/2})$. In these cases, the condition presented in Equation (4) will be satisfied for sufficiently large n and sufficiently small $\left|F\right|$. Moreover, a larger n allows for a larger value of $\left|F\right|$.

Note that Figure 1 only shows the results for small values of the triple $(\alpha ,\beta ,\gamma )\in {\mathbb{Z}}^3$. For the larger range $-50\le \alpha ,\beta ,\gamma \le 50$, the decay of $A_n$ with n is the same as in Figure 1 (not shown). In the remainder of this paper, we will restrict ourselves to small values of $(\alpha ,\beta ,\gamma )$. This is motivated by the observation that the Lorenz-96 model and its generalizations shown in Equation (3) have a similar structure as finite-difference discretisations of partial differential equations in which the interaction between the nonlinear terms are often local in nature; also see the discussion in [21].

2.2. Stability of Equilibria

All of the systems ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ have an equilibrium solution that is given by $x_F=(F,F,\cdots ,F)$. The eigenvalues of the Jacobian matrix at this equilibrium can be expressed explicitly in terms of $(\alpha ,\beta ,\gamma )$ and n, as the next result shows.

Lemma 1.

The eigenvalues of the Jacobian matrix of ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ at $x_F$ are given by

$$\begin{array}{c}\hfill {\lambda }_j=-1+F\left[\eta (j,n;\beta ,\gamma )+i\xi (j,n;\beta ,\gamma )\right],j=0,1,\cdots ,n-1,\end{array}$$

(5)

where

$$\begin{array}{cc}\hfill & \eta (j,n;\beta ,\gamma )=cos\left(\frac{2\pi j\beta }{n}\right)-cos\left(\frac{2\pi j\gamma }{n}\right),\hfill \\ \hfill & \xi (j,n;\beta ,\gamma )=sin\left(\frac{2\pi j\gamma }{n}\right)-sin\left(\frac{2\pi j\beta }{n}\right).\hfill \end{array}$$

Moreover, the eigenvector $v_j$ corresponding to the eigenvalue ${\lambda }_j$ is given by

$$v_j=\frac{1}{\sqrt{n}}{(1,{\rho }_j,\cdots ,{\rho }_j^{n-1})}^{\top },$$

where ${\rho }_j=e^{-2\pi ij/n}$ are the n-th roots of unity.

Proof. 

Without a loss of generality we may assume that $0\le \alpha ,\beta ,\gamma \le n-1$. Note that the Jacobian matrix of System (3) evaluated at $x_F$ is circulant, which means that each row is a cyclic right shift of the row above. If we denote the first row by $(c_0,c_1,\cdots ,c_{n-1})$, then it follows from [28] that the eigenvalues and corresponding eigenvectors are given by

$${\lambda }_j=\sum _{k=0}^{n-1}{c}_k{\rho }_j^k,v_j=\frac{1}{\sqrt{n}}{(1,{\rho }_j,\cdots ,{\rho }_j^{n-1})}^{\top },j=0,1,\cdots ,n-1.$$

Because we have that $c_0=-1$, $c_{\beta }=F$, $c_{\gamma }=-F$, and $c_k=0$ for all $k\ne 0,\beta ,\gamma$, it follows that

$$\begin{array}{cc}\hfill {\lambda }_j& =-1+F{\rho }_j^{\beta }-F{\rho }_j^{\gamma }\hfill \\ \hfill & =-1+Fexp\left(\frac{-2\pi ij\beta }{n}\right)-Fexp\left(\frac{-2\pi ij\gamma }{n}\right)\hfill \\ \hfill & =-1+F\left(cos\left(\frac{-2\pi j\beta }{n}\right)+isin\left(\frac{-2\pi j\beta }{n}\right)\right)-F\left(cos\left(\frac{-2\pi j\gamma }{n}\right)+isin\left(\frac{-2\pi j\gamma }{n}\right)\right)\hfill \\ \hfill & =-1+F\left(cos\left(\frac{2\pi j\beta }{n}\right)-cos\left(\frac{2\pi j\gamma }{n}\right)\right)+iF\left(sin\left(\frac{2\pi j\gamma }{n}\right)-sin\left(\frac{2\pi j\beta }{n}\right)\right).\hfill \end{array}$$

This completes the proof. □

Note that the parameter $\alpha$ does not influence the stability of the equilibrium $x_F$. Hence, in the remainder of the paper, we will mainly focus on the case $\alpha =-1$. The next result follows directly from Lemma 1.

Proposition 1.

For a fixed integer $n\ge 1$, we have the following:

• If $\eta (j,n;\beta ,\gamma )=0$ for some $j=0,1,\cdots ,n-1$, then $\Re \left({\lambda }_j\right)=-1$ and, hence, the eigenvalue ${\lambda }_j$ cannot cross the imaginary axis upon varying F. In particular, if $\beta =\pm \gamma$, then $\eta (j,n;\beta ,\gamma )=0$ for all $j=0,\cdots ,n-1$, which implies that the equilibrium $x_F$ is stable for all $F\in \mathbb{R}$.
• If $\eta (j,n;\beta ,\gamma )\ne 0$ for some $j=0,1,\cdots ,n-1$, then the eigenvalue ${\lambda }_j$ will cross the imaginary axis at the parameter value $F=1/\eta (j,n;\beta ,\gamma )$.
• If $\eta (j,n;\beta ,\gamma )\le 0$ for all $j=0,\cdots ,n-1$, then bifurcations of the equilibrium $x_F$ can only occur for $F<0$. In particular, this holds when $\gamma =0$.
• If $\eta (j,n;\beta ,\gamma )\ge 0$ for all $j=0,\cdots ,n-1$, then bifurcations of the equilibrium $x_F$ can only occur for $F>0$. In particular, this holds when $\beta =0$.

Lemma 1 also implies that the equilibrium $x_F$ is stable for $\left|F\right|$ sufficiently small. The shape of the graphs of the functions $\eta$ and $\xi$ strongly determines the nature of the first bifurcation of $x_F$.

2.3. Degenerate Cases

For specific choices of $(\alpha ,\beta ,\gamma )$, the bifurcations of the system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ can be degenerate. For example, consider the system ${\mathcal{L}}_{-1,0,2}\left(n\right)$. If $n=4k$, where $k\in \mathbb{N}$, Lemma 1 implies that ${\lambda }_k={\lambda }_{3k}=0$ for $F=\frac{1}{2}$, which means that the equilibrium $x_F$ becomes unstable by two real eigenvalues crossing zero. The occurrence of a Bogdanov–Takens bifurcation can be ruled out, since the Jacobian matrix at $x_F$ is diagonalizable [29]. In fact, the next result shows that two lines of equilibria appear at $F=\frac{1}{2}$.

Proposition 2.

For system ${\mathcal{L}}_{-1,0,-2}\left(n\right)$ with $n=4k$ and $F=\frac{1}{2}$, the following sets are equilibrium solutions:

$$\begin{array}{cc}\hfill L_1& =\left\{(\frac{1}{2},t,\frac{1}{2},1-t,\frac{1}{2},t,\frac{1}{2},1-t,\cdots ,\frac{1}{2},t,\frac{1}{2},1-t):t\in \mathbb{R}\right\},\hfill \\ \hfill L_2& =\left\{(t,\frac{1}{2},1-t,\frac{1}{2},t,\frac{1}{2},1-t,\frac{1}{2},\cdots ,t,\frac{1}{2},1-t,\frac{1}{2}):t\in \mathbb{R}\right\}.\hfill \end{array}$$

Proof. 

For $F=\frac{1}{2}$, the right hand side of System (3) reads as

$$f_j=x_{j-1}(x_j-x_{j-2})-x_j+\frac{1}{2},j=0,\cdots ,n-1.$$

A vector $(x_0,\cdots ,x_{n-1})\in {\mathbb{R}}^n$ belongs to $L_1$ if and only if one the following cases is satisfied:

• $(x_{j-2},x_{j-1},x_j)=(1-t,\frac{1}{2},t)$;
• $(x_{j-2},x_{j-1},x_j)=(t,\frac{1}{2},1-t)$;
• $(x_{j-2},x_{j-1},x_j)=(\frac{1}{2},t,\frac{1}{2})$; and,
• $(x_{j-2},x_{j-1},x_j)=(\frac{1}{2},1-t,\frac{1}{2})$.

Straightforward computations show that, in each of theses cases, $f_j=0$ for all $j=0,\cdots ,n-1$. This proves that elements of $L_1$ are indeed equilibria of System (3). The proof for $L_2$ is similar. □

The previous result implies that the equilibrium $x_F$ undergoes a degenerate bifurcation at $F=\frac{1}{2}$. For $F>\frac{1}{2}$, the equilibrium $x_F$ is unstable and numerical experiments show that orbits can become unbounded when their initial point is not contained in a compact invariant set, such as an equilibrium or periodic orbit.

The occurrence of multiple zero eigenvalues is not limited to the specific system ${\mathcal{L}}_{-1,0,-2}\left(n\right)$. More generally, the Jacobian matrix of the system ${\mathcal{L}}_{\alpha ,0,\gamma }\left(n\right)$ at $x_F$ has precisely $\left|\gamma \right|$ eigenvalues that are equal to zero when n is a multiple of $2\left|\gamma \right|$ and $F=\frac{1}{2}$. Indeed, if $\beta =0$ and $\gamma \ne 0$, then Lemma 1 implies that the eigenvalues are given by

$${\lambda }_j=-1+F\left(1-cos\left(\frac{2\pi j\gamma }{n}\right)\right)+iFsin\left(\frac{2\pi j\gamma }{n}\right).$$

Solving the equation ${\lambda }_j=0$ for $F=\frac{1}{2}$ gives

$$j=\frac{n}{2\left|\gamma \right|}(2k+1)\mathrm{where}k\in \mathbb{Z}.$$

The restriction $0\le j\le n-1$ leads to the restriction $0\le 2k+1<2\left|\gamma \right|$. Hence, if $2\left|\gamma \right|$ divides n, then there precisely exist $\left|\gamma \right|$ integers $0\le j\le n-1$ for which ${\lambda }_j=0$.

In the remainder of this paper, we will only consider systems ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ for which degeneracies, as described above, do not occur.

3. Finite Cascades of Pitchfork Bifurcations

In this section, we first discuss under which conditions on the triple $(\alpha ,\beta ,\gamma )$ the equilibrium $x_F=(F,\cdots ,F)$ loses stability through a pitchfork bifurcation. Next, we discuss the resulting bifurcation scenario for specific choices of $(\alpha ,\beta ,\gamma )$. Rather than providing an exhaustive analysis of all possible cases, we will restrict the discussion to three concrete examples. In particular, we discuss to what extent these examples differ from the original Lorenz-96 system.

3.1. Conditions for a First Pitchfork Bifurcation

For $x_F$ to loose stability through a pitchfork bifurcation, it is necessary that one of the eigenvalues equals zero, whereas all other eigenvalues have a negative real part. A sufficient condition on $(\alpha ,\beta ,\gamma )$ for which this holds is given in the next result.

Proposition 3.

Assume that, for $\beta ,\gamma \in \mathbb{Z}$, the function

$$\tilde{\eta }\left(x\right)=cos\left(\beta x\right)-cos\left(\gamma x\right)$$

attains a unique global maximum or minimum on the interval $[0,2\pi ]$ at $x=\pi$. Subsequently, it follows that:

• If n is odd, then zero eigenvalues of the equilibrium $x_F$ must occur at least in pairs; and,
• If n is even, then one eigenvalue equals zero for $F=-1/\tilde{\eta }\left(\pi \right)$, whereas the other eigenvalues have a negative real part.

Proof. 

Recall, from Lemma 1, that the eigenvalues ${\lambda }_j$, where $j=0,\cdots ,n-1$, of $x_F$ satisfy $\Re \left({\lambda }_j\right)=-1+\tilde{\eta }(2\pi j/n)$. It is straightforward to check that $\tilde{\eta }(2\pi -x)=\tilde{\eta }\left(x\right)$. This implies that $\Re \left({\lambda }_j\right)=\Re \left({\lambda }_{n-j}\right)$, which means that eigenvalues cross the imaginary axis in pairs.

When n is even, it follows that ${\lambda }_{n/2}=0$ for $F=-1/\tilde{\eta }\left(\pi \right)$. Because $\tilde{\eta }$ is assumed to have a unique global maximum or minimum at $x=\pi$, it immediately follows that all of the other eigenvalues have a negative real part at $F=-1/\tilde{\eta }\left(\pi \right)$. □

It is straightforward to verify that the Lorenz-96 model satisfies the conditions of Proposition 3. Other choices of $(\beta ,\gamma )$ given by $(-1,2)$, $(1,0)$, and $(0,-1)$ which will be discussed in more detail below. In addition, observe that

$$\tilde{\eta }\left(\pi \right)=\left\{\begin{array}{cc}-2\hfill & \mathrm{if}\gamma \mathrm{is}\mathrm{odd}\mathrm{and}\beta \mathrm{is}\mathrm{even},\hfill \\ \phantom{-}0\hfill & \mathrm{if}\gamma \mathrm{and}\beta \mathrm{have}\mathrm{the}\mathrm{same}\mathrm{parity},\hfill \\ \phantom{-}2\hfill & \mathrm{if}\gamma \mathrm{is}\mathrm{even}\mathrm{and}\beta \mathrm{is}\mathrm{odd}.\hfill \end{array}\right.$$

Therefore, if the equilibrium $x_F$ loses stability through a zero eigenvalue crossing, then this will necessarily occur at $F=1/2$ or $F=-1/2$.

A zero eigenvalue crossing is a signature of a saddle-node bifurcation, a transcritical bifurcation, or a pitchfork bifurcation. The saddle-node bifurcation can be ruled out, since the equilibrium $x_F$ continues to exist after the eigenvalue crossing has taken place. However, further analysis is needed in order to distinguish between a transcritical and a pitchfork bifurcation. A typical approach is to compute a normal form while using a center manifold reduction [29], but these computations are rather long. For the Lorenz-96 model, such computations are provided in [30], but below we shall adopt a more elementary and quicker approach. The next result shows, for certain members of the family ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$, a supercritical pitchfork bifurcation of the equilibrium $x_F$ takes place when $n=2k$ and $F=\pm 1/2$.

Proposition 4.

If $n=2k$, with $k\in \mathbb{N}$, then it follows that, for the systems ${\mathcal{L}}_{-1,1,-2}\left(n\right)$, ${\mathcal{L}}_{-1,-1,-2}\left(n\right)$, and ${\mathcal{L}}_{-1,-1,0}\left(n\right)$, the equilibrium $x_F$ loses stability through a supercritical pitchfork bifurcation at $F=-1/2$, while, in the system ${\mathcal{L}}_{-1,0,-1}\left(n\right)$, a supercritical pitchfork bifurcation occurs at $F=1/2$.

Proof. 

As an ansatz, we assume that System (3) has an equilibrium of the form $x_{P,1}=(a,b,a,b,\cdots ,a,b)$, in which case by circulant symmetry $x_{P,2}=(b,a,b,a,\cdots ,b,a)$ is also an equilibrium. For the system ${\mathcal{L}}_{-1,1,-2}\left(n\right)$, it then follows that

$$\begin{array}{cc}\hfill b(b-a)-a+F& =0,\hfill \\ \hfill a(a-b)-b+F& =0.\hfill \end{array}$$

(6)

Of course, $a=b=F$ is a solution, but this would lead to the already known equilibrium $x_F=(F,\cdots ,F)$. Solving a from the first equation gives $a=(b^2+F)/(b+1)$, so that the second equation yields a cubic equation for b:

$$-2b^3+(2F-2)b^2+(F-1)b+F(F+1)=0.$$

If $r_1$, $r_2$, and $r_3$ denote the roots of the latter equation, then Vieta’s formulas give

$$\begin{array}{cc}\hfill r_1+r_2+r_3& =F-1,\hfill \\ \hfill r_1{r}_2+r_2{r}_3+r_3{r}_1& =(1-F)/2,\hfill \\ \hfill r_1{r}_2{r}_3& =F(F+1)/2.\hfill \end{array}$$

Without loss of generality, we can take $r_3=F$, in which case we find

$$\begin{array}{cc}\hfill r_1+r_2& =-1,\hfill \\ \hfill r_1{r}_2& =(F+1)/2.\hfill \end{array}$$

Solving $r_1$ and $r_2$ from these equations is straightforward and it gives the following expressions for a and b:

$$a=\frac{-1+\sqrt{-1-2F}}{2}\mathrm{and}b=\frac{-1-\sqrt{-1-2F}}{2}.$$

This means that, for $F<-1/2$, two new equilibria appear that coalesce with $x_F$ at $F=-1/2$. Because the two new branches of equilibria extend in the direction, in which $x_F$ becomes unstable, we conclude that a supercritical pitchfork bifurcation takes place.

For the systems ${\mathcal{L}}_{-1,-1,-2}\left(n\right)$ and ${\mathcal{L}}_{-1,-1,0}\left(n\right)$, we obtain the same result and, therefore, the computations are omitted. For the system ${\mathcal{L}}_{-1,0,-1}\left(n\right)$, we obtain the equations

$$\begin{array}{cc}\hfill b(a-b)-a+F& =0,\hfill \\ \hfill a(b-a)-b+F& =0,\hfill \end{array}$$

which reduces to Equation (6) by substituting $(-a,-b,-F)$ for $(a,b,F)$. Hence, we obtain the solutions

$$a=\frac{1+\sqrt{-1+2F}}{2}\mathrm{and}b=\frac{1-\sqrt{-1+2F}}{2}.$$

This means that, for $F>1/2$, two new equilibria appear, which coalesce with $x_F$ at $F=1/2$. Again, because the two new branches of equilibria extend in the direction in which $x_F$ becomes unstable, we conclude that a supercritical pitchfork bifurcation takes place. □

3.2. Beyond the First Pitchfork Bifurcation

In [25], it was numerically shown for the Lorenz-96 model ${\mathcal{L}}_{-1,1,-2}\left(n\right)$ that the first supercritical pitchfork bifurcation is, in fact, followed by a finite cascade of subsequent supercritical pitchfork bifurcations. More precisely, for dimension $n=2^p$, the equilibrium $x_F$ and the subsequent branches that emanate from it undergo a finite cascade of p pitchfork bifurcations at the parameter values $F_{P,k}$ for $k=1,\cdots ,p$. The numerically computed parameter values that are listed in

Table 1. Numerically computed parameter values $F_{P,k}$ for which a pitchfork bifurcation takes place in the system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ with dimension $n=2^k$ for $k=1,\cdots ,9$. Note that ${\mathcal{L}}_{-1,1,-2}\left(n\right)$ is the original Lorenz-96 model. In addition, the ratio $r_k=(F_{P,k-1}-F_{P,k-2})/(F_{P,k}-F_{P,k-1})$ is listed, except for the system ${\mathcal{L}}_{-1,-1,0}\left(n\right)$, in which case the second pitchfork bifurcation is subcritical and the cascade terminates. See the main text for explanations.

	${\mathcal{L}}_{-1,1,-2}\left(\mathit{n}\right)$		${\mathcal{L}}_{-1,-1,-2}\left(\mathit{n}\right)$		${\mathcal{L}}_{-1,-1,0}\left(\mathit{n}\right)$	${\mathcal{L}}_{-1,0,-1}\left(\mathit{n}\right)$
$\mathit{k}$	${\mathit{F}}_{\mathit{P},\mathit{k}}$	${\mathit{r}}_{\mathit{k}}$	${\mathit{F}}_{\mathit{P},\mathit{k}}$	${\mathit{r}}_{\mathit{k}}$	${\mathit{F}}_{\mathit{P},\mathit{k}}$	${\mathit{F}}_{\mathit{P},\mathit{k}}$	${\mathit{r}}_{\mathit{k}}$
1	$-0.5000000$		$-0.5000000$		$-0.5000000$	$0.5000000$
2	$-3.0000000$		$-1.0000000$		$-1.0000000$	$0.6000000$
3	$-6.6000000$	$0.694$	$-1.0908967$	$5.501$		$0.6178351$	$5.607$
4	$-8.0107123$	$2.552$	$-1.1113355$	$4.447$		$0.6214900$	$4.880$
5	$-8.4360408$	$3.317$	$-1.1156935$	$4.690$		$0.6222663$	$4.708$
6	$-8.5275625$	$4.647$	$-1.1166297$	$4.655$		$0.6224320$	$4.685$
7	$-8.5474569$	$4.600$	$-1.1168299$	$4.676$		$0.6224675$	$4.668$
8	$-8.5517234$	$4.663$	$-1.1168729$	$4.656$		$0.6224751$	$4.671$
9	$-8.5526377$	$4.663$	$-1.1168821$	$4.674$		$0.6224767$	$4.750$

suggest that the following scaling law is satisfied:

$$\underset{k\to \infty }{lim}\frac{F_{P,k-1}-F_{P,k-2}}{F_{P,k}-F_{P,k-1}}={\delta }_F,$$

(7)

where ${\delta }_F\approx 4.669$ is Feigenbaum’s constant.

In addition to this scaling law, the pitchfork cascade satisfies another property, which has an analogy with the classical period doubling cascade of periodic points in iterated maps. We say that a vector $x=(x_0,\cdots ,x_{n-1})\in {\mathbb{R}}^n$ is periodic with period p if $x_{i+p}=x_i$ for all $i=0,\cdots ,n-1-p$. Numerical computations show that, after each pitchfork bifurcation, the period of newly born equilibria is being doubled. This observation forms the motivation for the ansatz, which was used in the proof of Proposition 4.

The following result implies that, if a pitchfork bifurcation occurs in system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ for dimension n, then this pitchfork bifurcation will also occur for dimensions that are multiples of n.

Proposition 5.

Assume that $x=(x_0,\cdots ,x_{n-1})\in {\mathbb{R}}^n$ is an equilibrium of system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$.

• For $k\in \mathbb{N}$, the point $x^{\prime }=(x,x,\cdots ,x)\in {\mathbb{R}}^{kn}$, where the coordinates of x are repeated k times, is an equilibrium of system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(kn\right)$.
• Denote the Jacobian matrices evaluated at these equilibria with $J_x\in {\mathbb{R}}^{n\times n}$ and $J_{x^{\prime }}\in {\mathbb{R}}^{kn\times kn}$, respectively. Subsequently, the spectrum of $J_x$ is contained in the spectrum of $J_{x^{\prime }}$.

In particular, if the equilibrium x undergoes a bifurcation at some parameter value $F=F_0$, then the equilibrium $x^{\prime }$ will undergo the same bifurcation at the same parameter value.

Proof. 

Note that the components of $x^{\prime }$ are related to those of x by

$$x_{j+in}^{\prime }=x_j\mathrm{for}\mathrm{all}i=0,\cdots ,k-1,j=0,\cdots ,n-1.$$

Denote, by $f_j$ and $f_j^{\prime }$, the right hand sides of System (3) for dimensions n and $kn$, respectively. Subsequently, by using the relation between the components of x and $x^{\prime }$, it follows that

$$\begin{array}{cc}\hfill f_{j+in}^{\prime }& =x_{j+in+\alpha }^{\prime }(x_{j+in+\beta }^{\prime }-x_{j+in+\gamma }^{\prime })-x_{j+in}^{\prime }+F\hfill \\ & =x_{j+\alpha }(x_{j+\beta }-x_{j+\gamma })-x_j+F\hfill \\ & =f_j.\hfill \end{array}$$

In particular, $f_j=0$ for all $j=0,\cdots ,n-1$ implies that $f_j^{\prime }=0$ for all $j=0,\cdots ,kn-1$. This proves statement 1.

If $v\in {\mathbb{R}}^n$ is any vector, then the j-th components of the righthand side of (3) evaluated at $x+v$ and $x-v$ are given by

$$\begin{array}{cc}\hfill g_j& :=(x_{j+\alpha }+v_{j+\alpha })(x_{j+\beta }+v_{j+\beta }-x_{j+\gamma }-v_{j+\gamma })-x_j-v_j+F,\hfill \\ \hfill h_j& :=(x_{j+\alpha }-v_{j+\alpha })(x_{j+\beta }-v_{j+\beta }-x_{j+\gamma }+v_{j+\gamma })-x_j+v_j+F,\hfill \end{array}$$

respectively. Because System (3) only has uadratic nonlinearities, it follows that the j-th component of the Jacobian matrix evaluated at x multiplied by v is given by

$$\frac{g_j-h_j}{2}=v_{j+\alpha }(x_{j+\beta }-x_{j+\gamma })+x_{j+\alpha }(v_{j+\beta }-v_{j+\gamma })-v_j.$$

In particular, if v is an eigenvector of the Jacobian matrix that corresponds to an eigenvalue $\lambda$, then it follows that

$$v_{j+\alpha }(x_{j+\beta }-x_{j+\gamma })+x_{j+\alpha }(v_{j+\beta }-v_{j+\gamma })-v_j=\lambda v_j\mathrm{for}\mathrm{all}j=0,\cdots ,n-1.$$

A similar reasoning as for statement 1 then shows that $v^{\prime }=(v,v,\cdots ,v)\in {\mathbb{R}}^{kn}$ is an eigenvector of the Jacobian matrix of (3) with dimension $kn$ that corresponds to the eigenvalue $\lambda$. This completes the proof.

However, note that Proposition 5 does not guarantee that the equilibrium $x^{\prime }$ in system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(kn\right)$ is stable whenever the equilibrium x in system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ is stable. Indeed, it may happen that the equilibrium $x^{\prime }$ undergoes bifurcations that do not occur for the equilibrium x. Indeed, becuase the matrix $J_{x^{\prime }}$ has more eigenvalues than the matrix $J_x$, there are more possibilities for the equilibrium $x^{\prime }$ to bifurcate. A concrete example for which this happens is system ${\mathcal{L}}_{-1,-1,-2}\left(n\right)$, which will be discussed below.

In the next sections, we will discuss to what extent a cascade of pitchfork bifurcation is observed in system ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ for three particular choices of $(\alpha ,\beta ,\gamma )$. The analysis will mainly rely on numerical computations that are performed by the continuation software AUTO-07p [31].

3.3. The System ${\mathcal{L}}_{-1,-1,-2}\left(n\right)$

In dimension $n=2^p$, a cascade of p pitchfork bifurcations occurs; see Table 1 for the parameter values up to $p=9$. Note that the numerical results suggest that the parameter values of the pitchfork bifurcations satisfy the Feigenbaum scaling of Equation (7). This bifurcation scenario is the same as for the Lorenz-96 model, albeit that the parameter values of the pitchfork bifurcations are different.

However, apart from this quantitative difference, there is also a qualitative difference. In the Lorenz-96 model with $n=8$ the four equilibria created at the second pitchfork bifurcation lose stability through a Hopf bifurcation before they bifurcate again through a third pitchfork bifurcation, and this leads to 8 unstable equilibria after the pitchfork cascade. For $(\alpha ,\beta ,\gamma )=(-1,-1,-2)$ and $n=8$ a Hopf bifurcation only takes place after the third pitchfork bifurcation; see Figure 2. For $n=16$, the scenario is again the same as for the Lorenz-96 model: now, a Hopf bifurcation occurs again after the second pitchfork bifurcation, as in the case $n=4$; see Figure 2.

3.4. The System ${\mathcal{L}}_{-1,1,0}\left(n\right)$

This case is very different from the Lorenz-96 model: for dimension $n=4$, the second pitchfork bifurcation is subcritical. This means that the two stable equilibria that were born at the first pitchfork bifurcation for $F_{P,1}=-1/2$ become unstable at $F_{P,2}=-1$ and four new unstable equilibria come into existence for $F>F_{P,2}$, as opposed to stable equilibria being created for $F<F_{P,2}$; see Figure 3. By Proposition 5, this bifurcation scenario will carry over to all dimensions n that are a multiples of four.

3.5. The System ${\mathcal{L}}_{-1,0,-1}\left(n\right)$

In dimension $n=2^p$, a succession of p pitchfork bifurcations occurs; see Table 1 for the parameter values up to $p=9$. Again, the numerical results suggest that the parameter values of the pitchfork bifurcations satisfy the Feigenbaum scaling of Equation (7). However, from a qualitative point of view, this case is different from the Lorenz-96 model. Numerical computations show that, up to $p=9$, the equilibria that were created in the pitchfork bifurcations do not lose stability before the last pitchfork bifurcation has occurred. Hence, there is a coexistence of $2^p$ stable equilibria. A Hopf bifurcation occurs only after the last pitchfork bifurcation, which leads to the coexistence of $2^p$ stable periodic orbits. The latter may bifurcate into chaotic attractors. Section 4 discusses some particular routes to chaos.

4. Multistability in the System ${\mathcal{L}}_{-1,0,-1}\left(n\right)$

The system ${\mathcal{L}}_{-1,0,-1}\left(n\right)$ is of particular interest. For odd dimensions n, it follows from Lemma 1 that the equilibrium $x_F=(F,\cdots ,F)$ loses stability through a complex conjugate pair of eigenvalues that cross the imaginary axis. The periodic orbits that are born through this Hopf bifurcation typically bifurcate further through period doublings or Neĭmark–Sacker bifurcations. Because these scenarios are very common, they will not be discussed further in the present paper. Instead, we will focus on dimension $n=2^p$ for which a cascade of p pitchfork bifurcations is followed by a Hopf bifurcation. In this section, we discuss the fate of the $2^p$ coexisting stable periodic orbits that arise in this way. To that end, we will numerically compute Lyapunov exponents as a function of the parameter F to detect the relevant bifurcations. In particular, we will show that the number of coexisting attractors may become smaller than $2^p$ due to collisions.

4.1. Lyapunov Exponents for Coexisting Attractors

A convenient way to detect bifurcations of attractors is by computing Lyapunov exponents as a function of the continuation parameter. In systems that exhibit multistability, i.e., the coexistence of attractors, care must be taken in selecting the initial condition. However, in System (3), coexisting attractors are related by circulant shifts of coordinates. We will now explain that the Lyapunov exponents for such attractors are equal.

Denote, with C, the $n\times n$ matrix that has the following circulant shift property:

$$C{(x_0,x_1,\cdots ,x_{n-1})}^{\top }=(x_1,\cdots ,x_{n-1},x_0).$$

If $f:{\mathbb{R}}^n\to {\mathbb{R}}^n$ denotes the right-hand side of Equation (3), then $Cf\left(x\right)=f\left(Cx\right)$ for all $x\in {\mathbb{R}}^n$. In particular, if $x\left(t\right)$ is a solution of (3), then it is $Cx\left(t\right)$. If two vectors $u,v\in {\mathbb{R}}^n$ satisfy the relation $u=Cv$, then the fact that f only has quadratic nonlinearities implies that

$$CDf\left(x\right)v=\frac{1}{2}C\left[f(x+v)-f(x-v)\right]=\frac{1}{2}\left[f(Cx+u)-f(Cx-u)\right]=Df\left(Cx\right)u.$$

In particular, if $v\left(t\right)$ satisfies the variational equation $\dot{v}=Df\left(x\left(t\right)\right)v$, then $u\left(t\right)=Cv\left(t\right)$ satisfies the variational equation $\dot{u}=Df\left(Cx\left(t\right)\right)u$. The usual algorithm for computing Lyapunov exponents consists of integrating the variational equations and applying the Gram–Schmidt procedure to the variational vectors [32,33]; see [34,35] for alternative methods. It is clear that the matrix C is unitary and, thus, preserves the Euclidean inner product. Hence, the Lyapunov exponents that are computed for the orbits $x\left(t\right)$ and $Cx\left(t\right)$ must be the same. By induction, the Lyapunov exponents will be the same for the orbits $C^{k}x\left(t\right)$ for all $k=0,\cdots ,n-1$.

For a fixed parameter value F, we first perform a transient integration in order to obtain an initial condition $x\left(0\right)$ on an attractor. For the attractor, we then compute the Lyapunov exponents while using the algorithm described in [32,33]. As explained above, the attractors with initial condition $C^{k}x\left(0\right)$ have the same Lyapunov exponents for all $k=0,\cdots ,n-1$. During the computation, we also minimize the $x_0$-coordinate along the attractors:

$${\mu }_k=min\left\{{\left[C^{k}x\left(t\right)\right]}_0:0\le t\le T\right\},$$

(8)

where T is the total integration time is used in our computations. If ${\mu }_j\ne {\mu }_k$ for $j\ne k$, then the attractors with initial points $C^{j}x\left(0\right)$ and $C^{k}x\left(0\right)$ are different. In case that ${\mu }_j={\mu }_k$, we have a strong indication that the attractors with initial points $C^{j}x\left(0\right)$ and $C^{k}x\left(0\right)$ are equal. Hence, the quantities ${\mu }_k$ help to detect collisions of attractors. To the final point $x\left(T\right)$ on the attractor, a small random perturbation is added, the parameter F is increased by a small step size, and the computations are repeated. The two largest Lyapunov exponents ${\lambda }_1\ge {\lambda }_2$ and the quantities ${\mu }_0,\cdots ,{\mu }_{n-1}$ are then plotted as a function of F. In the resulting diagram, bifurcations of attractors manifest themselves as parameter values for which a Lyapunov exponent becomes zero, and collisions of attractors manifest themselves as parameter values for which at least two of the quantities ${\mu }_k$ attain the same value.

4.2. Dimension $n=4$

Figure 4 shows the Lyapunov diagram for dimension $n=4$. After the second pitchfork bifurcation at $F=0.6$, four stable equilibria coexist. At $F\approx 0.64546$, these equilibria become unstable through a supercritical Hopf bifurcation, which leads to the coexistence of four stable periodic orbits; for $F=0.672$, these orbits are shown in Figure 5 (left panel). At $F\approx 0.6740$, one pair of periodic orbits collides with an equilibrium, while another pair of periodic orbits collides with a different equilibrium. This leads to the creation of four homoclinic orbits that form two curves in a “figure eight” configuration, see Figure 5 (middle panel). After the collision, each of the “figure eight” curves split into two coexisting stable periodic orbits, see Figure 5 (right panel). This shows how the number coexisting attractors can change due to the collision of an attractor with another invariant object.

A similar scenario, as described above, can also increase the number of coexisting attractors. At $F\approx 0.7080$, two stable periodic orbits collide with an equilibrium. This gives rise to four homoclinic orbits forming two “figure eight” curves. After the homoclinic bifurcation, there are again four coexisting stable periodic orbits, see Figure 6.

4.3. Dimension $n=6$

Figure 7 shows the Lyapunov diagram for dimension $n=6$. Despite the dimension being larger, the bifurcation scenario is somewhat simpler in this case, since periodic orbits do not collide with equilibria, as in the case $n=4$. There is only one pitchfork bifurcation at $F=0.6$ after which two stable equilibria coexist. Each of these equilibria becomes unstable through a supercritical Hopf bifurcation at $F\approx 0.61611$, which leads to the coexistence of two stable periodic orbits. Each of these bifurcates through a period doubling cascade, which leads to the coexistence of two chaotic attractors. Figure 8 shows periodic orbits after one period doubling ($F=0.6495$), after two period doublings ($F=0.6539$), and the chaotic attractors after the cascade ($F=0.6545$). Chaotic attractors that arise after a period doubling cascade are typically observed to have a Hénon-like structure, which means that the attractors are formed by the closure of the unstable manifold of saddle periodic orbit [23,36,37].

4.4. Dimension $n=8$

Figure 9 shows the Lyapunov diagram for dimension $n=8$. After the last pitchfork bifurcation at $F\approx 0.61784$, eight stable equilibria coexist. Each of these equilibria becomes unstable at $F\approx 0.62336$ through a supercritical Hopf bifurcation, which leads to the coexistence of eight stable periodic orbits. These periodic orbits undergo a period doubling bifurcation at $F\approx 0.6255$ and a period halving bifurcation at $F\approx 0.6263$. At $F\approx 0.6268$, the periodic orbits disappear through a saddle-node bifurcation. For $F>0.6268$, eight different stable periodic orbits coexist. The latter have more windings than the formerly existing periodic orbits, which suggests that the new periodic orbits already have bifurcated through a succession of period doubling bifurcations. Finally, a period doubling cascade takes place, which leads to the formation of eight coexisting chaotic attractors; for $F=0.62698$, these are plotted in Figure 10.

5. Discussion

In this paper, we have considered a family of generalized Lorenz-96 models, of which each member is parameterized by a triple $(\alpha ,\beta ,\gamma )\in {\mathbb{Z}}^3$. In particular, we have shown for two concrete members in this family that a finite cascade of pitchfork bifurcations takes place, just as in the original Lorenz-96 model. In particular, in the system ${\mathcal{L}}_{-1,0,-1}\left(n\right)$, the equilibria remain stable up to the last pitchfork bifurcation. This is a qualitative difference with the Lorenz-96 model, in which equilibria already undergo a Hopf bifurcation after the second pitchfork bifurcation. The Hopf bifurcations of the stable equilibria lead to the coexistence of periodic attractors. Further bifurcations of the latter can result in the coexistence of chaotic attractors. The number of attractors that coexist can change due to collisions.

The results that are presented in this paper give rise to several open questions that warrant further research. The first question is whether the pitchfork cascade in the models ${\mathcal{L}}_{-1,1,-2}\left(n\right)$, ${\mathcal{L}}_{-1,-1,-2}\left(n\right)$, and ${\mathcal{L}}_{-1,0,-1}\left(n\right)$ continues indefinitely when the dimension n increases. Numerical computations show the occurrence of pitchfork bifurcations up to dimension $n=512$, but this does not guarantee that the pitchfork bifurcations will persist for $n>512$. Indeed, in the model ${\mathcal{L}}_{-1,1,0}\left(n\right)$, the second pitchfork bifurcation becomes subcritical and the cascade comes to a halt.

How to check whether the pitchfork bifurcations persist for dimensions $n>512$? An attempt at explicitly computing all of the equilibria and their bifurcations would not be useful for at least two reasons. The first reason is that an analytical computation is simply not feasible due to the fact that the equilibria born through pitchfork bifurcations will have complicated expressions involving radicals; already after two pitchfork bifurcations, analytical expressions become unwieldy. The second reason is that an explicit computation would only provide information regarding one particular model, whereas our numerical computations suggest that properties of the pitchfork cascade are common to at least three members of the family ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$.

Another open question is whether there exist dynamical systems with a finite-dimensional state space, in which an infinite cascade of pitchfork bifurcations occurs. Most likely, such systems cannot have only polynomial nonlinearities. Indeed, as a consequence of Bezout’s theorem, generically one expects only finitely many equilibria in such systems. Exceptions are degenerate cases, such as system ${\mathcal{L}}_{-1,0,-2}\left(n\right)$, which has infinitely many equilibria for certain parameter values. However, for systems with transcendental nonlinearities, the situation might be different. These open questions will have to be addressed in future research.

It is an intriguing numerical observation that the parameter values of the pitchfork bifurcations that are described above satisfy the Feigenbaum scaling in Equation (7). The latter scaling was originally discovered in period doubling cascades of periodic points in unimodal maps [38,39]. A key feature is its universal nature, which means that the scaling holds for an entire class of unimodal maps possesing a quadratic maximum, rather than just one particular map. Lanford [40] and Eckmann and Wittwer [41] provided computer assisted proofs of the Feigenbaum conjectures. A computer-free proof was given by Lyubich [42]. Briggs [43,44] carried out high-precision computations of the Feigenbaum constants.

The fact that the Feigenbaum scaling for the pitchfork bifurcations occurs in at least three different members of the family ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ suggests that it might be a universal property within a suitably chosen class of systems. However, it is not clear a priori how this could be analytically proven. Proof strategies, as in [40,41,42], are not applicable in the context of System (3) because of several conceptual differences. Firstly, (3) is not a map but a differential equation. Secondly, the pitchfork cascade in System (3) is finite and the dimension needs to increase in order to observe a longer cascade. If one is interested in universality, then the first step is to identify an appropriate class of systems for which a universal result can be conjectured. Perhaps pitchfork cascades can occur in systems that are different from Equation (3) provided a sufficient amount of symmetry is present. However, to the best of our knowledge, examples of such systems are unknown to us.

Finally, the results that are presented in this paper also illustrate the important point that bifurcation scenarios in the models ${\mathcal{L}}_{\alpha ,\beta ,\gamma }\left(n\right)$ strongly depend on divisibility properties of the dimension n. The relevance of this phenomenon extends beyond the context of this paper. As an example, for discretisations of Burgers’ equation, it has been observed that, for odd degrees of freedom, the dynamics are confined to an invariant subspace, whereas, for even degrees of freedom, this is not the case. In turn, this has an effect on the bifurcation diagrams that one observes in this model [45]. Similar phenomena can be expected in the discretisations of partial differential equations, in which periodic boundary conditions lead to circulant symmetries.