Travelling waves and their bifurcations in the Lorenz-96 model

TL;DR

This paper analyzes the Lorenz-96 model's bifurcations and traveling waves, revealing supercritical Hopf bifurcations, quasi-periodic attractors, and routes to chaos through analytical and numerical methods.

Contribution

It provides explicit formulas for bifurcation analysis and demonstrates the role of Hopf-Hopf bifurcations as organizing centers in the Lorenz-96 model.

Findings

First Hopf bifurcation is always supercritical.

Traveling wave solutions emerge at bifurcation.

Routes to chaos vary with no clear pattern as dimension increases.

Abstract

In this paper we study the dynamics of the monoscale Lorenz-96 model using both analytical and numerical means. The bifurcations for positive forcing parameter $F$ are investigated. The main analytical result is the existence of Hopf or Hopf-Hopf bifurcations in any dimension $n\geq4$. Exploiting the circulant structure of the Jacobian matrix enables us to reduce the first Lyapunov coefficient to an explicit formula from which it can be determined when the Hopf bifurcation is sub- or supercritical. The first Hopf bifurcation for $F>0$ is always supercritical and the periodic orbit born at this bifurcation has the physical interpretation of a travelling wave. Furthermore, by unfolding the codimension two Hopf-Hopf bifurcation it is shown to act as an organising centre, explaining dynamics such as quasi-periodic attractors and multistability, which are observed in the original Lorenz-96 model. Finally, the region of parameter values beyond the first Hopf bifurcation value is investigated numerically and routes to chaos are described using bifurcation diagrams and Lyapunov exponents. The observed routes to chaos are various but without clear pattern as $n\rightarrow\infty$.

Tables

Table 1. Table 2: Colour coding for the Lyapunov diagram in Figure 19 .

Colour	Type of attractor
Red	Stable equilibrium
Blue	Periodic attractor for $l=2$
Green	Periodic attractor for $l=3$
Grey	Quasi-periodic attractor
Black	Chaotic attractor

Equations

$${\dot{x}}_j=x_{j-1}(x_{j+1}-x_{j-2})-x_j+F,j=1,\dots ,n,$$

$$x_{j-n}=x_{j+n}=x_j,$$

$$E=\frac{1}{2}\sum _j{x}_j^2.$$

$$x_F=(F,\dots ,F),$$

$$(c_0,c_1,\dots ,c_{n-1}),$$

$${\lambda }_j=\sum _{k=0}^{n-1}{c}_k{\rho }_j^k,{\rho }_j=\exp \left(-2\pi i\frac{j}{n}\right),j=0,\dots ,n-1,$$

$${\lambda }_j(F,n)$$

$$=-1+F\left(\exp \left(-2\pi i\frac{j}{n}\right)-\exp \left(4\pi i\frac{j}{n}\right)\right)$$

$$=-1+Ff(j,n)+Fg(j,n)i,$$

$$f(j,n)$$

$$g(j,n)$$

$$v_j=\frac{1}{\sqrt{n}}\left(\begin{array}{c}1\\ {\rho }_j\\ {\rho }_j^2\\ ⋮\\ {\rho }_j^{n-1}\end{array}\right).$$

$${\lambda }_j$$

$$v_j$$

$$F_{\min }(n)$$

$$F_{\max }(n)$$

$$\cos \frac{2\pi l_1}{n}+\cos \frac{2\pi l_2}{n}=\frac{1}{2}.$$

$${\ell }_1(l,n)=\frac{4}{n}\tan (\frac{\pi l}{n}){\sin }^2(\frac{3\pi l}{n})\frac{5\cos (\frac{2\pi l}{n})+8\cos (\frac{4\pi l}{n})-2\cos (\frac{6\pi l}{n})-8}{4\cos (\frac{2\pi l}{n})-4\cos (\frac{4\pi l}{n})+9}.$$

$$5\cos y_0+8\cos 2y_0-2\cos 3y_0-8=0,$$

$$l_1(n)=\underset{0<l<n/3}{\mathrm{argmax}}f(l,n),$$

$$\frac{n}{6}\le l_1(n)\le \frac{n}{4},$$

$${\lambda }_l=-{\omega }_{0}i\text{with}{\omega }_0=\frac{\cos \frac{\pi l}{n}}{\sin \frac{\pi l}{n}},$$

$$P(t)=\sqrt{F-F_{\operatorname{H}}}\big{(}\cos(\omega_{0}t)\operatorname{Re}(v_{l})-\sin(\omega_{0}t)\operatorname{Im}(v_{l})\big{)}.$$

$$P_j(t)=\sqrt{\frac{F-F_H}{n}}\cos \left({\omega }_{0}t-\frac{2\pi j}{n}l\right),$$

$$T=\frac{2\pi }{{\omega }_0}=2\pi \tan \frac{\pi l}{n}.$$

$$T_{\infty }=\underset{n\to \infty }{\lim }2\pi \tan \left(\frac{\pi l_1(n)}{n}\right)=2\pi \tan (\frac{1}{2}\arccos (\frac{1}{4}))\approx \mathrm{4.867.}$$

$$\underset{n\to \infty }{\lim }\frac{n}{l_1(n)}=\frac{2\pi }{\arccos (\frac{1}{4})}\approx \mathrm{4.767.}$$

$$\underset{n\to \infty }{\lim }\frac{2\pi l_1(n)}{n}=\arccos (\frac{1}{4}),$$

$${\dot{x}}_j=x_{j-1}(x_{j+1}-x_{j-2})-x_j+G(x_{j-1}-2x_j+x_{j+1})+F,j=1,\dots ,n,$$

$$(c_0,c_1,\dots ,c_{n-1})=(-1-2G,F+G,0,\dots ,0,-F,G).$$

$${\kappa }_j(F,G,n)$$

$$=-1-2G+(F+G){\rho }_j^1-F{\rho }_j^{n-2}+G{\rho }_j^{n-1}$$

$$=-1-2G+(F+G)\exp \left(-\frac{2\pi j}{n}i\right)-F\exp \left(\frac{4\pi j}{n}i\right)+G\exp \left(\frac{2\pi j}{n}i\right)$$

$$=-1-2G\left(1-\cos \frac{2\pi j}{n}\right)+Ff(j,n)+Fg(j,n)i,$$

$$G=H_l(F,n)=\frac{Ff(l,n)-1}{2(1-\cos \frac{2\pi l}{n})},$$

Full text

1 Introduction

1.1 Setting of the problem

In his 1996 paper [Lorenz06], Edward Lorenz introduced two models to study fundamental issues regarding the predictability of the atmosphere and weather forecasting. The so-called monoscale Lorenz-96 model is defined by the equations

[TABLE]

resulting in a model with circulant symmetry. For the multiscale model, which will not be discussed in this paper, the reader is referred to [Lorenz06]. The model (1) can be interpreted as a model for atmospheric waves travelling along a circle of constant latitude. Lorenz interpreted the variables $x_{j}$ as values of some meteorological quantity (e.g., temperature, pressure, or vorticity) in $n$ equal sectors of a latitude circle, where the index $j$ plays the role of longitude. The continuous parameter $F$ represents external forcing and can be used as a bifurcation parameter.

Although the Lorenz-96 model is not derived from physical principles it still has features which are commonly found in geophysical models: forcing, dissipation and energy preserving quadratic terms. Moreover, unlike the traditional Lorenz-63 model [Lorenz63] which has only one positive Lyapunov exponent, the Lorenz-96 model has multiple positive Lyapunov exponents for suitable choices of the parameters $F$ and $n$. For those reasons, and for the simplicity of the equations, the model is important and widely used nowadays and sometimes even called “the archetype of large deterministic systems displaying chaotic behavior” [Fatkullin04] or “a hallmark representative of nonlinear dynamical systems” [Frank14]. The applications of the Lorenz-96 model are broad and range from geophysical applications like data assimilation and predictability to studies in spatiotemporal chaos. Table 1 gives an overview of recent papers in which the Lorenz-96 has been used together with the values of the parameters that were used.

In contrast to its importance, only a few studies have investigated the dynamics of this model. In [Karimi10], the high-dimensional chaotic dynamics has been explored by means of the fractal dimension. A recent study on patterns of order and chaos in the multiscale model has reported the existence of regions with standing waves [Frank14]. Bifurcation diagrams in low dimensions of the Lorenz-96 model have been studied in [Orrell03], although the emphasis of their work was on methods to visualise bifurcations by means of spectral analysis, rather than exploring the dynamics itself. The previous works already revealed an extraordinarily rich structure of the dynamical behaviour of the Lorenz-96 model for specific values of $n$. However, there has been no systematic study of the dynamics of this model yet. In this paper we fill this gap by studying the dynamical nature of the Lorenz-96 model in greater detail and give analytical proofs of some basic properties for all dimensions and of the existence of Hopf and Hopf-Hopf bifurcations. These results are complemented by numerical explorations, that includes the dynamics beyond these bifurcations as well.

The Lorenz-96 model is a family of dynamical systems parameterised by the discrete parameter $n\in\mathbb{N}$ which gives the dimension of their state space. This setup is analogous to a discretised partial differential equation. In fact, in some works the Lorenz-96 model is interpreted as such [Basnarkov12, Reich15]. In [Lucarini07], a discretised quasi-geostrophic model for the atmosphere was studied. In particular, they numerically observed that the parameter value at which the first Hopf bifurcation occurs typically increases with the truncation order of their discretisation method. In pseudo-spectral discretisations of Burgers’ equation [Basto06] qualitative differences in dynamics were observed depending on whether the dimension of state space was even or odd. We may expect similar phenomena for the Lorenz-96 model. Hence, in this paper we focus in particular on the question which quantitative and qualitative features of the dynamics will persist for (almost) all $n\in\mathbb{N}$. Answers to these questions may be helpful in selecting appropriate values of $n$ and $F$ for the specific applications listed in Table 1. For example, there is a direct relation between the dimension of attractors and the statistics of extreme events [Holland12]. Although the study of this paper is unable to provide the entire picture, it offers a partial dynamical inventory using both analytical and numerical means.

1.2 Sketch of the results

The Lorenz-96 model (1) has an equilibrium solution given by $x_{F}=(F,\dots,F)$ for all $n\geq 1$ and $F\in\mathbb{R}$. Clearly, for $F=0$ this equilibrium is stable. Numerical simulations show that for $F=1.2$ the dynamics of the model is periodic for all $n\geq 4$. This suggests that for $0<F<1.2$ a supercritical Hopf bifurcation occurs at which the equilibrium $x_{F}$ loses its stability and gives birth to a periodic attractor. Figure 1 shows that the period of the periodic attractor at $F=1.2$ is an oscillating function of the dimension $n$. Observe that the oscillations decay with $n$ and that the period seems to converge to a value of approximately 4.5. The spatiotemporal properties of these periodic orbits are further explored in Figure 2 by means of so-called Hovmöller diagrams [Hovmoller49]. In these diagrams the value of the variables $x_{j}(t)$ is plotted as a function of time $t$ and “longitude” $j$. Clearly, the waves are travelling in the direction of decreasing $j$ and their wave number increases with $n$. The numerical results clearly indicate that spatiotemporal properties of travelling waves in the Lorenz-96 model depend on the dimension $n$.

A natural question is then which properties stabilize in the limit $n\to\infty$. In this paper we will present the following results which provide a (partial) answer to this question:

• •

For all $n\geq 4$ we prove that the trivial equilibrium $x_{F}=(F,\dots,F)$ exhibits several Hopf or Hopf-Hopf bifurcations for the parameter value $F_{\operatorname{H}}\in\mathbb{R}$. In case of a Hopf bifurcation, we also prove whether the bifurcation is sub- or supercritical.

• •

For $F>0$ we prove that the first Hopf bifurcation takes place for $F\in(\tfrac{8}{9},1.19)$ and is always supercritical. The periodic attractor born at this bifurcation has the physical interpretation of a travelling wave. At the Hopf bifurcation the period of this wave is an oscillating function of $n$ which tends to $T_{\infty}\approx 4.86$ as $n\to\infty$ and the wave number increases linearly with $n$.

• •

We encounter further bifurcations of the stable orbit beyond the value $F_{\operatorname{H}}$ for which the first Hopf bifurcation takes place. Eventually, this leads to chaotic behaviour. The diagram in Figure 3 shows the bifurcations following only the stable orbit for various dimensions $n$. Also, the parameter value where chaos sets in is indicated. A clear pattern for all $n$ can not be discerned from the diagram, though a pattern is observed for dimensions $n\leq 100$ where $n$ is a multiple of $5$.

• •

To unfold the codimension two Hopf-Hopf bifurcation we add an extra parameter to the original model (1) via a Laplace-like diffusion term in such a way that the original model is easily retrieved. The thus obtained two-parameter system clarifies the role of the Hopf-Hopf bifurcation as organising centre and so it sheds more light on the original model, especially for $n=12$ in which case the Hopf-Hopf bifurcation is the first bifurcation for $F>0$.

From our results we can conclude that the spatiotemporal properties of the Lorenz-96 model depend on $n$. This again shows the importance of selecting an appropriate value of $n$ in specific applications. The presence of a Hopf or Hopf-Hopf bifurcation persists for all $n\geq 4$ and the $F_{\operatorname{H}}$-value of the first of these bifurcations converges to $\tfrac{8}{9}$ as $n\to\infty$. However, since the resulting waves have different wave numbers the subsequent bifurcation patterns vary with $n$. The case $n=5m$ seems to follow a more regular pattern and will be discussed in Section 3.2. Furthermore, the linear increase of the wave number with $n$ indicates that the Lorenz-96 model cannot be interpreted as a discretised partial differential equation.

1.3 Outline

This paper has been organised as follows. Section 2 presents the analytical results of the research. It starts with some general properties of the Lorenz-96 model, followed by the main theorems of this paper, which concern the Hopf and Hopf-Hopf bifurcations. Subsequently, we show that the periodic orbits, resulting from these bifurcations, can be interpreted as travelling waves. After this, we describe the unfolding of the Hopf-Hopf bifurcation. Next, Section 3 is concerned with a numerical analysis of the dynamics after the first bifurcation for positive $F$. The proofs of all analytical results follow in Section 4.

2 Analytical results

In this section we give an overview of the basic properties of the Lorenz-96 model and state the main theorems of this paper.

2.1 General properties

Trapping Region

The Lorenz-96 model is constructed such that the quadratic part does not affect the total energy of the system,

[TABLE]

This property implies the existence of a trapping region in all dimensions [Lorenz84a]. We will prove this in Proposition 7 in the more general setting of a two-parameter unfolding of the Lorenz-96 model.

Equilibrium

The previous result on the existence of a global attractor shows that the attractors are found in a neighbourhood of the origin. It is easy to see that system (1) has in any dimension the trivial equilibrium

[TABLE]

where we take $F\in\mathbb{R}$.

Due to the symmetry in our system we are able to determine the eigenvalues of $x_{F}$ explicitly. In any dimension, the Jacobian matrix at this equilibrium is a circulant matrix, which means that each row is a right cyclic shift of the row above it. Let us first look at the case $n<4$, which differs from the general case $n\geq 4$. If $n=1$ or $3$, then the Jacobian matrix is equal to minus the identity matrix and so the eigenvalues are all equal to $-1$. In case $n=2$, the first row of the Jacobian matrix is given by $(-1-F,F)$ and, hence, its eigenvalues are $\lambda_{0}=-1$ and $\lambda_{1}=-1-2F$. Since all these eigenvalues never become complex, one can conclude that we do not have any Hopf bifurcation of the trivial equilibrium for dimensions $n<4$.

Next, consider the case $n\geq 4$. Denote the first row of the Jacobian matrix at $x_{F}$ by

[TABLE]

where $c_{0}=-1$, $c_{1}=F$, $c_{n-2}=-F$ and $c_{k}=0$ for $k\neq 0,1,n-2$. It follows that the eigenvalues are given by [Gray06]

[TABLE]

which can be expressed in terms of $F,j$ and $n$:

[TABLE]

where $f$ and $g$ are defined as

[TABLE]

The eigenvector $v_{j}$ corresponding to $\lambda_{j}$ can also be expressed in terms of $\rho_{j}$:

[TABLE]

Observe that the eigenvalue $\lambda_{0}$ equals $-1$ and its corresponding eigenvector has all entries equal to 1. Due to the fact that $\rho_{n-j}=\bar{\rho}_{j}$, all the other eigenvalues and eigenvectors form complex conjugate pairs since

[TABLE]

except when $n$ is even, in which case the eigenvalue for $j=\tfrac{n}{2}$ is real and equals $\lambda_{n/2}=-1-2F$. Note as well that when $n$ is a multiple of 3 the eigenvalues for $j=\tfrac{n}{3},\tfrac{2n}{3}$ both equal $-1$. We call the pair $\{\lambda_{j},\lambda_{n-j}\}$ the $j$-th eigenvalue pair. In the next subsection we will see that each complex eigenvalue pair has a particular parameter value $F$ for which it crosses the imaginary axis and thus causes a Hopf bifurcation.

2.2 Hopf Bifurcations

The main part of this section is devoted to bifurcations of the trivial equilibrium (2) for positive values of $F$. To find bifurcations in the Lorenz-96 model, we take $F$ as the bifurcation parameter and vary it along the real line. In this way we discover the occurrence of several Hopf and Hopf-Hopf bifurcations, which we summarize here in two theorems. These are preceded by a lemma which proves that we have the desired eigenvalue crossing that is needed for both cases. The proofs of these results are postponed to Section 4.

Before we formulate our results on the Hopf-Hopf and Hopf bifurcation in Theorem 2 and 4, respectively, let us first state the preliminary result:

Lemma 1** (Eigenvalue crossing).**

Let $n\geq 4$ and $l\in\mathbb{N}$ such that $0<l<\tfrac{n}{2},l\neq\tfrac{n}{3}$, then the following holds:

The $l$-th eigenvalue pair $\{\lambda_{l},\lambda_{n-l}\}(F,n)$ of the trivial equilibrium $x_{F}$ of system (1) crosses the imaginary axis transversally at the parameter value $F_{\operatorname{H}}(l,n):=1/f(l,n)$ and thus the equilibrium changes stability. 2. 2.

$F_{\operatorname{H}}(l,n)$* lies in the domain $\left(F_{\min}(n),-\tfrac{1}{2}\right)\cup\left[\tfrac{8}{9},F_{\max}(n)\right)$ with*

[TABLE]

where $r$ is the quotient of $n$ after division by 3.

Due to the shape of the function $f(j,n)$, at most two eigenvalue pairs can have zero real part simultaneously for a particular value of $F$ (see Figure 24). This indicates that the crossing of eigenvalue pairs, described in Lemma 1, can lead to Hopf bifurcations and Hopf-Hopf bifurcations only.

Hopf-Hopf bifurcations

Let us first describe the Hopf-Hopf case: suppose that we have two distinct eigenvalue pairs with $l_{1}$ and $l_{2}$ which both cross the imaginary axis at the same parameter value $F_{\operatorname{HH}}:=F_{\operatorname{H}}(l_{1},n)=F_{\operatorname{H}}(l_{2},n)$. In that case we have a Hopf-Hopf bifurcation:

Theorem 2** (Hopf-Hopf Bifurcation).**

Let $l_{1},l_{2}$ and $n$ satisfy the assumptions in Lemma 1 with $l_{1}\neq l_{2}$. Then the trivial equilibrium $x_{F}$ exhibits a Hopf-Hopf bifurcation at $F_{\operatorname{HH}}$ if and only if $l_{1}$ and $l_{2}$ satisfy

[TABLE]

From equation (7) we can deduce two infinite sequences of dimensions for which a Hopf-Hopf bifurcation takes place:

Corollary 3**.**

Let $m\in\mathbb{N}$, then a Hopf-Hopf bifurcation occurs if we select $l_{1},l_{2}$ and $n$ according to one of the following criteria:

$n=10m$* and $l_{1}=m,l_{2}=3m$, which corresponds to $F_{\operatorname{HH}}=2$;* 2. 2.

$n=12m$* and $l_{1}=2m,l_{2}=3m$, which corresponds to $F_{\operatorname{HH}}=1$.*

Remark*.*

Equation (7) gives a necessary and sufficient condition for the occurrence of a Hopf-Hopf bifurcation. However, the explicit values of $l_{1},l_{2}$ and $n$, given in Corollary 3, possibly do not provide all occasions where a Hopf-Hopf bifurcation occurs.

Hopf bifurcations

On the other hand, if equation (7) is not satisfied then we have only one eigenvalue pair crossing the imaginary axis, which implies a Hopf bifurcation:

Theorem 4** (Hopf Bifurcation).**

Let $l$ and $n$ be as in Lemma 1. If the $l$-th eigenpair is the only one crossing the imaginary axis at the corresponding parameter value $F_{\operatorname{H}}(l,n)$, then the equilibrium $x_{F}$ exhibits a Hopf bifurcation at $F_{\operatorname{H}}$. The first Lyapunov coefficient for this bifurcation is given by

[TABLE]

Fix $y_{0}\in(0,\pi)$ such that

[TABLE]

then $\ell_{1}(l,n)$ is

• •

positive if $l$ and $n$ satisfy $0<\tfrac{l}{n}<\tfrac{y_{0}}{2\pi}\approx 0.08825746$, which corresponds to a subcritical bifurcation;

• •

negative if $\tfrac{l}{n}\in\left(\tfrac{y_{0}}{2\pi},\tfrac{1}{2}\right)\setminus\{\tfrac{1}{3}\}$ holds, which corresponds to a supercritical bifurcation.

The number of possible Hopf bifurcations for a given dimension $n$ is exactly equal to the number of conjugate eigenvalue pairs which satisfy Lemma 1. Using equation (6), we can count the number of such eigenvalue pairs by the number of eigenvalues with $0<j<\tfrac{n}{2}$, which gives the number $\lceil n/2-1\rceil$ (we need the ceiling-function here if $n$ is odd). However, as described above, if $n$ is a multiple of 3, then the eigenvalue pair with $j=\tfrac{n}{3}$ is not complex, so in this case the number of such eigenvalue pairs equals $\lceil n/2-2\rceil$. For the actual number of Hopf bifurcations, these numbers should be reduced by the number of Hopf-Hopf bifurcations.

We now restrict our attention to the parameter range $F\geq 0$. For $F=0$ the equilibrium $x_{F}$ is clearly stable. We are interested in the smallest value of $F>0$ at which the equilibrium bifurcates and becomes unstable. From the previous results we know that this must be either a Hopf or a Hopf-Hopf bifurcation.

Proposition 5**.**

Let $n\geq 4$ be fixed. For $F>0$ the first Hopf or Hopf-Hopf bifurcation occurs for the eigenpair with index

[TABLE]

which satisfies the bounds

[TABLE]

except for $n=7$, in which case we have to take $l_{1}=1$.

In particular, if the first bifurcation is a Hopf bifurcation, then this bifurcation is supercritical.

Remark*.*

For $n=12$ the first bifurcation is not a Hopf bifurcation, but a Hopf-Hopf bifurcation. In this case we have $l_{1}(n)=\{2,3\}$ as shown in Corollary 3; see Section 2.4.2 for a more detailed exposition of the dynamics in this case.

2.3 Travelling waves

A fluid is said to be hydrodynamically unstable when small perturbations of the flow can grow spontaneously, drawing energy from the mean flow. At a Hopf bifurcation an equilibrium loses its stability and gives birth to a periodic orbit. In the context of a fluid this can be interpreted as a steady flow becoming unstable to an oscillatory perturbation, such as a travelling or standing wave. Hopf bifurcations are found in many geophysical models as the first bifurcation which destabilizes a steady flow [Lucarini07, Broer11, Dijkstra05, Sterk10]. In this section we show that periodic orbits of the Lorenz-96 model born at a Hopf bifurcation, as described in the previous section, can be interpreted as travelling waves.

Under the conditions of Theorem 4 a Hopf bifurcation associated with the $l$-th eigenpair occurs at $F_{\operatorname{H}}=1/f(l,n)$. In this case equation (2.1) gives

[TABLE]

where we take $\omega_{0}$ by convention to be the absolute value of the imaginary part at the bifurcation value. Observe that

[TABLE]

is a periodic orbit of the Lorenz-96 model which is linearized around the trivial equilibrium $x_{F}$. For $F-F_{\operatorname{H}}>0$ sufficiently small the function $P(t)$ is a good approximation of the periodic orbit of the nonlinear system, which allows us to determine the physical properties of the wave.

Using equation (5) gives the $j$-th component of $P(t)$ as

[TABLE]

which is indeed the expression for a travelling wave. In this expression the integer $l$ is the wave number and $\tfrac{2\pi j}{n}$ plays the role of discrete longitude. The temporal frequency of the wave is given by $\omega_{0}$ which implies that its period is given by

[TABLE]

The level curves of $P_{j}(t)$ are given by the lines $\omega_{0}t-\tfrac{2\pi j}{n}l=\text{constant}$, which are decreasing in the $(j,t)$-plane. This implies that waves travel in the direction of decreasing $j$, which is indeed observed in Figure 2.

Furthermore, it is easy to see that a small $l$ implies a small period and a large wave length (i.e. a long and fast wave), whereas a larger $l$ gives a larger period and a smaller wave length (i.e. a short and slow wave). Proposition 5 shows that the wave number for the first bifurcation for $F\geq 0$, $l_{1}(n)$, increases linearly with $n$. In the limit $n\rightarrow\infty$ the following relations hold for the wave number $l_{1}(n)$ and the period:

Proposition 6**.**

In the limit $n\to\infty$, the period of the periodic attractor born at the first Hopf bifurcation is given by

[TABLE]

Similarly, the quotient of $n$ with the wave number $l_{1}(n)$ satisfies

[TABLE]

Proof.

Both results follow immediately from the fact that

[TABLE]

by the maximum of $\tilde{f}(y)=\cos y-\cos 2y$, combined with equation (8). ∎

Proposition 6 explains the features that can be observed in Figures 1 and 2. The period of the waves tends to a finite limit as $n\to\infty$ whereas the spatial wave number is unbounded.

Figure 4 compares the period of the periodic attractor for parameter values near and at the bifurcation value. The black curve is taken from Figure 1 and shows the period of the periodic attractor, computed numerically at the value $F=1.2$ for all $n\in[4,100]$. The red curve is the period of the periodic attractor computed via the theoretical formula (8) for $l=l_{1}(n)$ with $n\in[4,100]$, i.e. exactly at the bifurcation value $F_{\operatorname{H}}(l_{1},n)$. The difference is caused by the fact that the value $F=1.2$ used for the numerical computation is in most cases more than $0.2$ apart from the exact value $F_{\operatorname{H}}(l_{1},n)$.

2.4 Unfolding for the Hopf-Hopf bifurcation

In many systems one finds bifurcations of codimension less than or equal to the dimension $p$ of the parameter space that are subordinate to certain bifurcation points of codimension $p+1$. To study the qualitative dynamics around such codimension $p+1$ points, one can embed the system in a family of systems parameterised by $p+1$ parameters [Wiggins03]. In this way, the codimension $p+1$ bifurcations act as organising centres of the bifurcation diagram.

Corollary 3 shows that for $n=12$ the trivial equilibrium $x_{F}$ loses stability through a codimension two Hopf-Hopf bifurcation at $F_{\operatorname{HH}}=1$. In order to unfold this codimension two bifurcation completely an extra parameter is needed. In the following we propose an unfolding for general $n$ and take the codimension two Hopf-Hopf bifurcation as an organising centre for this family of systems. Thereafter we will describe the case $n=12$ in more detail and show the role of the Hopf-Hopf bifurcation as organising centre. In section 3.4 we will show by numerical computations how it dominates the dynamics in its neighbourhood and even influences the phase space for larger parts of the parameter space. A similar approach is taken in [BroerSaleh05, BroerSaleh07, Broer02], showing the existence and influence of several codimension three, respectively two, points that act as an organising centre.

We choose to add a Laplace-like diffusion term multiplied by a new parameter $G$ to our original equations (1a) to obtain the two-parameter system for general dimension $n$ with equations

[TABLE]

and keep the boundary condition (1b) as it is. The parameter value $G=0$ returns the original Lorenz-96 model.

2.4.1 General dimensions

Trapping region

We start the exploration of the two-parameter Lorenz-96 system with the following observation:

Proposition 7** (Trapping Region).**

The two-parameter system (9) has a trapping region for any dimension $n\in\mathbb{N}$ and for all $F\in\mathbb{R}$ and $G>-\tfrac{1}{4}$.

This result means that we can expect an attractor to exist in the region with $G>-\tfrac{1}{4}$, but in the half-plane below this line an attractor does not necessarily exist. Therefore, the results for the two-parameter system are relevant only for parameter values $G>-\tfrac{1}{4}$.

Equilibrium and Hopf bifurcations

The two-parameter system (9) has the same trivial equilibrium (2) for all $F,G\in\mathbb{R}$. The Jacobian matrix at this equilibrium is a circulant matrix with first row equal to

[TABLE]

We can express its eigenvalues $\kappa_{j}$ with $j=0,\dots,n-1$ easily in terms of $F$ and $G$ [Gray06]:

[TABLE]

where $f$ and $g$ are given by equation (4). The corresponding eigenvector $v_{j}$ is again given by (5). Note that equation (6) holds for $\kappa$ instead of $\lambda$ with similar restrictions, that is, if $j=\tfrac{n}{2}$, we have $\kappa_{n/2}=-1-2F-4G$; and if $j=\tfrac{n}{3},\tfrac{2n}{3}$, then we have $\kappa_{j}=-1-3G$. This shows that the case $l=\tfrac{n}{3}$ can not give a Hopf bifurcation and that the additional restriction $0<l<\tfrac{n}{2}$ provides all possible complex eigenvalue pairs. Note also that at $F=0$ all eigenvalues are real, so no Hopf bifurcation is then possible.

The following lemma demonstrates that the two-parameter system (9) can exhibit as many different Hopf bifurcations as system (1).

Lemma 8**.**

Let $n\geq 4$ and $l\in\mathbb{N}$ such that $0<l<\tfrac{n}{2},l\neq\tfrac{n}{3}$, then the trivial equilibrium $x_{F}$ of system (9) exhibits a Hopf bifurcation on the linear bifurcation curves

[TABLE]

where $F\in\mathbb{R}\backslash\{0\}$.

The Hopf bifurcation points are now turned into straight lines in the $(F,G)$-plane. Along these curves (11) it is possible to determine the type of the bifurcation by computing the first Lyapunov coefficient explicitly in a similar manner as done in the proof of Theorem 4 (see Section 4.3), but we will not repeat the procedure here.

Hopf-Hopf bifurcations

It is obvious that the intersections of the Hopf-lines cause Hopf-Hopf bifurcations. One can find all Hopf-Hopf bifurcation points of the trivial equilibrium by equating two Hopf-lines from formula (11) with different $l$. Under the assumption that all nondegeneracy conditions are satisfied, the truncated normal form for a Hopf-Hopf bifurcation reads [Kuznetsov04]

[TABLE]

after a suitable choice of phase variables. (Note that we use indices 2 and 3 for these phase variables which is most convenient for our discussion of the case $n=12$ in Section 2.4.2.) Here, $\mu_{j}$ is defined as $\mu_{j}:=\operatorname{Re}\kappa_{j}$, $\sigma=\pm 1$ and $\vartheta,\delta,\Theta,\Delta$ are other normal form coefficients. The sign $\sigma$ mainly determines the type of behaviour near the Hopf-Hopf point [Kuznetsov04]. In any case, two Neimark-Sacker bifurcation (NS) curves emanate from the Hopf-Hopf point. The directions of these NS-curves depend on $\vartheta$ and $\delta$ and can be computed up to first order via the real part of the eigenvalues at the Hopf-Hopf point [Kuznetsov04]. Note that the type of dynamics we have around the Hopf-Hopf bifurcation point does not depend on the choice of unfolding (9), since the normal form coefficients should be evaluated at the bifurcation point.

2.4.2 Unfolding for $n=12$

In this section we describe the interesting situation of dimension $n=12$ more explicitly, using the results we obtained for general dimensions. Therefore, consider system (9) and let $n=12$. The eigenvalues of the Jacobian matrix belonging to the trivial equilibrium $x_{F}$ are given by equation (2.4.1).

A Hopf bifurcation for the $l$-th eigenvalue pair (with $0<l<6,l\neq 4$) occurs along the Hopf-lines in equation (11). So, we obtain explicitly the following Hopf bifurcation curves as a function of $F\in\mathbb{R}\setminus\{0\}$:

[TABLE]

The curves for $l=2,3$ intersect each other at $(F,G)=(1,0)$, which is the Hopf-Hopf bifurcation point we discovered in the original system. It is easy to see that for $G=0$ this is the first bifurcation of $x_{F}$ one encounters by increasing the parameter $F$, since the only Hopf bifurcation with positive $F_{\operatorname{H}}$-value has $F_{\operatorname{H}}(1,12)\approx 2.732051>F_{\operatorname{HH}}$. Observe that there is no resonance present at the Hopf-Hopf point (since there $\omega_{2}=\operatorname{Im}(\kappa_{2}(1,H_{2}(1,12),12))=\sqrt{3}$ and $\omega_{3}=\operatorname{Im}(\kappa_{3}(1,H_{3}(1,12),12))=1$). Furthermore, in our particular situation, numerical computation of the normal form coefficients using MatCont [Dhooge11] yields the following values

[TABLE]

showing that the bifurcation is indeed nondegenerate and that the normal form (12) is valid.

From the value of $\sigma$ it follows that the dynamics of system (9) is of “type I in the simple case” as described by [Kuznetsov04]. This means that the two NS-curves are the only bifurcation curves that emanate from the Hopf-Hopf point and in between these curves there exists a region in the $(F,G)$-plane where two stable periodic orbits coexist with an unstable 2-torus. These NS-curves correspond to the limit cycles with $l=2,3$ and are approximated by, respectively,

[TABLE]

Removing higher order terms and solving for $G$ gives the following linear curves in $F$:

[TABLE]

These lines are tangent to the real NS-curve at the Hopf-Hopf point.

Figure 5 displays the local bifurcation diagram with the Hopf-lines and approximated NS-curves for $l=2,3$ together with the phase portraits for each region. In the next section, we verify these results numerically.

3 Dynamics beyond the Hopf and Hopf-Hopf bifurcation

In Section 2 we have proven that the trivial equilibrium $x_{F}$ of the Lorenz-96 model (1) loses stability through either a supercritical Hopf or a Hopf-Hopf bifurcation for all dimensions $n\geq 4$. At these bifurcations a periodic attractor is born which has the physical interpretation of a travelling wave (see Section 2.3).

In this section we explore the dynamics of the Lorenz-96 model numerically for dimensions up to $n=100$ and beyond the first Hopf bifurcation. Our emphasis is on the bifurcations through which the periodic attractor loses stability. A natural question is to what extent these bifurcations depend on the dimension $n$. Moreover, for a few dimensions we comment on routes to chaos using tools such as continuation, integration, Poincaré sections and Lyapunov exponents. The numerical analysis was carried out using mainly the original Lorenz-96 model (1). The two-parameter model, however, turns out to be useful to explain the features observed in the one-parameter model. Whenever the two-parameter system (9) was used, this is stated explicitly. Otherwise, $G$ is assumed to be equal to [math]. (Recall that we retrieve the original model from the two-parameter system by setting $G=0$.)

An overview of the bifurcations for several dimensions is presented in the diagram in Figure 3. As can be seen from this diagram there are various routes to chaos. We will discuss a few of them below. Moreover, for all dimensions shown, except for $n=4$, chaos sets in for $F\in(3,7)$.

3.1 Dimension $n=4$

In the four dimensional Lorenz-96 model there is only one Hopf bifurcation at $F_{\operatorname{H}}(1,4)=1$. Continuing the periodic attractor originating from this bifurcation in $F$ and plotting its period against $F$ gives the diagram in Figure 6. The original periodic orbit disappears through a fold bifurcation at $F_{SN}\approx 11.83823$. Chaos is observed for parameter values $F\geq 11.84$. Figure 7 compares the periodic attractor for $F=11.83$ with the chaotic attractor for $F=11.9$, while Figure 8 shows time series of the first variable for both parameter values. Observe that the dynamics for $F=11.9$ alternates between approximate periodic behaviour and chaotic behaviour. This is the classical type 1 intermittency scenario as described in [Pomeau80, Eckmann81]. Note that for intermittency we not only need an attractor that has disappeared, but we also need the global dynamics to be such that a typical evolution recurrently visits the part of state space where the attractor disappeared. In our case, this global dynamics might be provided by a heteroclinic structure, as we will show below.

At $F\approx 8.540498$ an additional limit cycle appears through a fold bifurcation of limit cycles, which is stable for only a short interval. This bifurcation is followed by more fold bifurcations, which accumulate at $F\approx 11.77$, as can be seen from Figure 6. This phenomenon suggests a homoclinic or heteroclinic structure [Kuznetsov04]. Similar behaviour has been observed in other atmospheric models [Veen03]. Analysis of the system for this parameter value indicates a heteroclinic structure. At $F\approx 8.898979$, namely, two pairs of four equilibria appear through a fold bifurcation. By numerical continuation we found the following coordinates for these equilibria at $F=12.081216$ (the importance of this value will become clear in a moment):

[TABLE]

while the other six equilibria can be obtained by a cyclic shift of the entries. Both types of equilibria are hyperbolic saddles with three, resp. two, stable eigenvalues. However, only at $F\approx 12.081216$ (which is in the chaotic region) we have numerically detected a heteroclinic cycle between the equilibria $x^{1}$, using MatCont. A continuation of these connections in the $(F,G)$-plane for the two-parameter system does not yield any other value $F$ for which a heteroclinic cycle exist at $G=0$. The heteroclinic cycle for $(F,G)\approx(12.081216,0)$ is shown in Figure 9. Notice the similarity between the right panel and the attractor in Figure 7.

3.2 Dimension $n=5$

For $n=5$, the first bifurcation after the Hopf bifurcation at $F_{\operatorname{H}}(1,5)\approx 0.8944272$ is a period-doubling bifurcation (PD) which occurs at $F_{\operatorname{PD},1}\approx 3.937853$. The next three PDs occur for the parameter values $F_{\operatorname{PD},2}\approx 4.981884$, $F_{\operatorname{PD},3}\approx 6.371496$, $F_{\operatorname{PD},4}\approx 6.640968$, respectively.

The bifurcation diagrams in Figure 10 suggest that a cascade of period-doubling bifurcations takes place. After the cascade, a chaotic attractor is detected at $F=6.72$, see Figure 11. The Poincaré section of this attractor appears to have the structure of a fattened curve. This suggests that the attractor is of Hénon-like type, which means that it is the closure of an unstable manifold of an unstable periodic point of the Poincaré map. We have numerically detected an unstable periodic orbit at $F=6.72$ which corresponds to an unstable period-3 point for the Poincaré return map to the section $\Sigma=\{x_{1}=5\}$. The unstable manifold of this period-3 point was computed with standard numerical techniques which are described in [Simo90]. Figure 12 shows a magnification of the unstable manifold along with the attractor of the Poincaré map. The two plots are in very good agreement with each other. Therefore, we conjecture the attractor in Figure 11 to be the closure of the unstable manifold of a saddle periodic orbit.

The PDs persist for each $n=5m$ with $m=1,\dots,10$. In all these cases, the bifurcation values of the first Hopf bifurcation and the first period-doubling are exactly the same as in the case of $n=5$ (see Figure 13). From $n=55$ on the pattern deviates, because the parameter value of the first Hopf bifurcation changes: now a Neimark-Sacker bifurcation (NS) is the first bifurcation after the Hopf bifurcation, but the torus originating from this bifurcation disappears for slightly larger $F$ and we seem to have again the PD-pattern with a PD at $F_{\operatorname{PD}}=3.937853$ as before.

We conjecture that this phenomenon can be explained by the wave number of the periodic attractor after the Hopf bifurcation. It turns out that for $n=5m$ with $m=1,\dots,10$ the wave number of this attractor is exactly equal to $m$. Hence, $n/l\in\mathbb{N}$ and therefore the attractor has repeating coordinates with $x_{j+5}=x_{j}$ for all $1\leq j\leq n$ and indices modulo $n$. This implies that the travelling wave is contained in an invariant subspace of $\mathbb{R}^{n}$ defined by

[TABLE]

where the index $j$ has to be taken modulo $n$. For any $n=5m$, the dynamics restricted to $W^{5}$ is governed by the Lorenz-96 model for $n=5$. We will further explore this property in forthcoming work [Kekem17b]. However, for $n\geq 55$ this phenomenon breaks down, since the wave number $l$ of the periodic attractor no longer satisfies the relation $n/l\in\mathbb{N}$. Nevertheless, it can happen that this periodic attractor becomes unstable and that again the periodic attractor with wave number $m$ takes up stability. However, this is not guaranteed to happen in general, especially for high dimensions, since the quotient $n/l_{1}(n)$ converges to a non-integer number for $n\rightarrow\infty$ as shown in Proposition 6. Therefore, for increasing $n$ there is an increasing number of periodic attractors whose wave numbers $l$ satisfy $\tfrac{n}{l_{1}(n)}<\tfrac{n}{l}<5$ and so it will become more rare to find a stable periodic attractor which inherits its dynamics from the case $n=5$.

3.3 Dimension $n=6$

For $n=6$ the first bifurcation after the Hopf bifurcation at $F_{\operatorname{H}}(1,6)=1$ is an NS, which occurs at $F_{\operatorname{NS}}\approx 5.456661$. At this bifurcation the periodic attractor loses stability and gives birth to a quasi-periodic attractor in the form of a 2-dimensional torus, see Figure 14. The Lyapunov diagram in Figure 15 clearly shows alternating intervals of periodic behaviour and quasi-periodic behaviour. This phenomenon can be clarified by the two-parameter system. In the $(F,G)$-plane this alternation organizes itself in the form of the well-known Arnol’d resonance tongues, which emanate from the NS-curve [Kuznetsov04]. For a better visualization of these tongues the affine transformation $(F,G)=(U+6V+1,0.35\,V-0.25)$ has been used to obtain the right panel of Figure 15.

3.4 Dimension $n=12$

Part of the dynamics for this dimension is already explained in Section 2.4.2. Here we present the results of our numerical exploration which support the analytical results very well.

Recall that the first bifurcation of the trivial equilibrium for $G=0$ is a Hopf-Hopf bifurcation, which is rather exceptional for the Lorenz-96 model. This codimension two point acts as an organising centre, as explained in Section 2.4. Two codimension one NS-curves originate from this bifurcation point, each corresponding to one of the wave numbers $l=2$ or $3$. The local bifurcation diagram obtained by MatCont is presented in Figure 16 and should be compared with the analytically computed bifurcation diagram in Figure 5.

In the region enclosed by both NS-curves multistability occurs, due to the coexistence of the periodic attractors for both $l=2$ and $l=3$. Both attractors are plotted for the same parameter values $(F,G)=(1.5,0)$ in Figure 17. Together with their Hovmöller diagrams in Figure 18, this shows that both waves are of a different nature. Multistability is also reflected by the Lyapunov diagrams in Figure 19. The left (resp. right) panel is obtained by fixing the value of the parameter $F$ and increasing (resp. decreasing) the parameter $G$. Along each vertical line in the parameter plane we have used the last point on the attractor detected in the previous step as an initial condition for the next one. In both diagrams we have used a grid of size 1000 by 1000. See Table 2 for an explanation of the colouring for each region. Figure 19 clearly shows that there is a region in the parameter plane where two different periodic attractors coexist. Also note that the bifurcation curves of Figure 16 are clearly visible in these diagrams. Lastly, Figure 19 shows that the Hopf-Hopf bifurcation influences a large portion of the parameter space as well as the phase space.

3.5 Dimension $n=36$

For $n=36$, we observe again coexistence of attractors, like in the case $n=12$. The Hopf-Hopf bifurcation that induces this phenomenon occurs at the intersection of the Hopf-lines for wave numbers $l=7$ and $l=8$ where $(F,G)=(0.919586,0.0144084)$, i.e. close to the $F$-axis. Note that these wave numbers correspond to the first two Hopf bifurcations of the trivial equilibrium for $F>0$ and $G=0$. The Hopf-Hopf bifurcation is of the same type as for $n=12$ (see Section 2.4.2), meaning that only two NS-curves arise from the codimension two point. The local bifurcation diagram in Figure 20 shows these two curves together with their corresponding Hopf-lines. The blue NS-curve (corresponding to $l=7$) intersects the line $G=0$ at $F\approx 0.9092541$, so we can observe multistability in the one-parameter model (1) for $F$ larger than this value. Again, the Hopf-Hopf bifurcation point acts as an organising centre.

In Figures 21 and 22 the Lyapunov diagrams are shown for $l=7$ and $l=8$, respectively, with $G=0$ fixed. For wave number $l=7$, the first bifurcation after the Hopf bifurcation at $F_{\operatorname{H}}(7,36)\approx 0.902474$ is the mentioned NS at $F\approx 0.9092541$, which is followed by another NS at $F\approx 4.389139$. The resulting quasi-periodic attractor then bifurcates to a 3-torus (see below). For $l=8$, a stable periodic attractor originates from a Hopf bifurcation at $F_{\operatorname{H}}(8,36)\approx 0.898198$. This attractor exhibits a PD at $F\approx 3.155456$ and becomes unstable via a subcritical NS at $F\approx 3.162597$, which can be seen from the right panel of Figure 22. The only stable attractor after $F\approx 3.162597$ is the one with wave number $l=7$. This is reflected in the Lyapunov diagrams of Figure 22, where the Lyapunov exponents take up the values for $l=7$ right after the subcritical NS at $F\approx 3.162597$ (compare with Figure 21). These observations show that the region of multistability is bounded for $G=0$.

The Lyapunov diagram in the right panel of Figure 21 suggests that for $G=0$ a 3-torus exists in a small interval of $F$-values before chaotic attractors are observed. Figure 23 shows a 3-torus attractor for $(n,F,G)=(36,4.45,0)$ together with a Poincaré section defined by $\Sigma=\{x_{1}=2\}$. This type of behaviour has also been observed for $n=24$ (not shown). Newhouse, Ruelle and Takens [Newhouse78] proved that small perturbations of a quasi-periodic flow on the 3-torus can lead to strange Axiom A attractors. Concrete routes of the NRT-scenario were reported in [Broer08, Broer08A] in the setting of a model map for the Hopf-saddle-node bifurcation in diffeomorphisms. Some techniques to study bifurcations of 3-tori in continuous-time dynamical systems are described in [Kamiyama15]. Unravelling the bifurcations of 3-tori and the associated routes to chaos in the Lorenz-96 model is left for future research.

4 Proofs of the analytical results

In this section we prove the statements in Section 2, in the same order.

4.1 Proof of Lemma 1

Proof of statement 1.

Let $n\geq 4$ and $l\in\mathbb{N}$. To investigate the $l$-th eigenvalue pair of $x_{F}$ we need to consider $0<l<\tfrac{n}{2}$ only, by equation (6). Let the $l$-th eigenvalue be written as $\lambda_{l}(F,n):=\mu(F)+\omega(F)i$, with real and imaginary parts as in (2.1). Then $\lambda_{l}(F,n)$ has zero real part if and only if

[TABLE]

Since the value of $f$ is already fixed by choosing $l$ and $n$, this is achieved only if $F$ equals

[TABLE]

where we need the additional condition that $l\neq\tfrac{n}{3}$, since $f(n/3,n)=0$ (see Figure 24).

The eigenvalues cross the imaginary axis with nonzero speed by the fact that

[TABLE]

due to the constraints on $l$. Moreover, let $\omega_{0}=|\omega(F_{\operatorname{H}})|:=-F_{\operatorname{H}}g(l,n)$ denote the absolute value of the imaginary part at the Hopf bifurcation. Then, by the restrictions on $l$, $\omega_{0}$ is nonzero as well. ∎

For later purposes, we note that at the eigenvalue crossing the dependence of the eigenvalues on $F$ can be replaced by the dependence on $l$, by substituting $F=F_{\operatorname{H}}(l,n)$. This gives

[TABLE]

for all $j=0,\ldots,n-1$. One can easily see that for $j=l$, $\lambda_{l}$ is purely imaginary at $F_{\operatorname{H}}$, i.e. $\lambda_{l}(l,n)=-i\omega_{0}$, where $\omega_{0}$ is also expressed in $l$ and $n$:

[TABLE]

for $0<l<\tfrac{n}{2},l\neq\tfrac{n}{3}$.

Proof of statement 2.

Let $\tilde{f}$ be defined by $\tilde{f}(y):=\cos y-\cos 2y$ such that $\tilde{f}(\tfrac{2\pi l}{n})=f(l,n)$, i.e. $\tilde{f}$ is equal to the function $f$ with the discrete points $\tfrac{2\pi j}{n}$ replaced by the continuous variable $y$, see Figure 24. By the definition of $F_{\operatorname{H}}$ in (17), its positive, respectively negative, values with the smallest absolute value occur at the maximum, respectively minimum, of the function $\tilde{f}$.

The extreme values of $\tilde{f}$ are obtained by

[TABLE]

and, hence, $y=k\pi$ for $k\in\mathbb{Z}$ or $\cos y=\tfrac{1}{4}$. The cases $y=k\pi$ give the global minimum of $\tilde{f}$ if $k$ is odd: $\tilde{f}(\pi)=-2$ (if $k$ is even, we have a local minimum). Therefore, the upper bound on the negative values of $F_{\operatorname{H}}$ is equal to $\tfrac{1}{\tilde{f}(\pi)}=-\tfrac{1}{2}$. This value can be obtained from $f(l,n)$ only if we take $l=\tfrac{n}{2}$, which is, however, excluded by the assumptions on $l$ and, hence, will never be attained. The solution $y_{\text{top}}=\cos^{-1}\left(\tfrac{1}{4}\right)$ gives the maximum of $\tilde{f}$ as $\tilde{f}(y_{\text{top}})=\tfrac{9}{8}$. Therefore the lowest possible positive value for which a Hopf bifurcation can occur is $\tfrac{1}{\tilde{f}(y_{\text{top}})}=\tfrac{8}{9}$.

Upper bound of $F_{\operatorname{H}}$

The positive values with the largest absolute value are obtained as follows. Since there are only finitely many $l$ satisfying $0<l<\tfrac{n}{2}$, we know that these values of $F_{\operatorname{H}}$ should be bounded for any $n$. The largest value for $F_{\operatorname{H}}$ is obtained when $f(l,n)$ is the smallest.

Claim**.**

The smallest value of $f(l,n)$ for every $n\geq 4$ is obtained at $l=1$ except when $n=7$, in which case we have to take $l=2$.

Proof of the claim.

For $n=4,5$ or $6$, this is trivial since $l=1$ is the only integer satisfying $0<l<\tfrac{n}{3}$. In the case $n=7$, we have two integers that satisfy $0<l<\tfrac{n}{3}$ and it is easily checked that $f(1,7)>f(2,7)$ gives the desired exception. So, all we need to show in order to verify our claim is that

[TABLE]

holds for any $n\geq 8$, where $l^{\prime}$ is the largest integer for which $l^{\prime}<\tfrac{n}{3}$. Since the function $f$ becomes negative for $l>\tfrac{n}{3}$, this $l^{\prime}$ will be the largest integer for which $f$ is positive and becomes close to 0, see Figure 24. This gives rise to the following three cases:

If $n=0\bmod 3$, then $l^{\prime}=\tfrac{n}{3}-1$. Then equation (20) can be simplified to

[TABLE]

Since we have to consider $n\geq 9$ here, the first sine-term is always positive, while the second one is always negative (by the fact that its entry is negative and bigger than $-\tfrac{\pi}{6}$). Hence, $f(1,n)-f(l^{\prime},n)<0$ holds for these particular values of $n$. 2. 2.

If $n=1\bmod 3$, then $l^{\prime}=\tfrac{n-1}{3}$. Then equation (20) reduces to

[TABLE]

In this case, we have to take $n\geq 10$, for which the first sine-term and the cosine-term are both always positive. The second sine-term, $\sin\tfrac{10-n}{6n}\pi$ is exactly equal to 0 if $n=10$ (compare this with criterion 1 in Corollary 3) and strictly less than 0 for $n>10$. This gives the desired inequality. 3. 3.

If $n=2\bmod 3$, then $l^{\prime}=\tfrac{n-2}{3}$ and we have that

[TABLE]

Observe that the part in brackets is monotonically decreasing, since the first sine-term decreases as $n$ increases, while both components $\cos\tfrac{\pi}{n}$ and $\sin\tfrac{\pi}{3}(1-\tfrac{2}{n})$ increase. Moreover, for $n=8$ — the least possible $n$ in this case — we find that

[TABLE]

which implies that also in this case equation (20) holds for any $n\geq 8$.

Now, we have established equation (20) for any $n\geq 4$, except $n=7$, and for these $n$ we can conclude that $l=1$ is the right choice to get the lowest value of $f(l,n)$. ∎

Continuation of proof of statement 2. We can conclude that the upper bound (actually, a maximum) on $F_{\operatorname{H}}$ is given by

[TABLE]

Lower bound of $F_{\operatorname{H}}$

In the case of negative $F_{\operatorname{H}}$-values, we do not have an eigenvalue crossing for $n=4$ and 6, so $F_{\min}=-\tfrac{1}{2}$ (this gives the empty set). For all other $n$ we need the integer $l>\tfrac{n}{3}$ which is the closest to $\tfrac{n}{3}$ on the right to get the largest value of $f$, see Figure 24. If we write $n=3r+s$, where $r,s\in\mathbb{N}$ with $s=n\bmod 3$, then we see that

[TABLE]

Therefore, we have to take $l=r+1$ to obtain the lowest integer satisfying $l>\tfrac{n}{3}$. Hence, the lower bound on $F_{\operatorname{H}}$ is given by

[TABLE]

which is again an actual minimum for given $n$. ∎

4.2 Proofs for the Hopf-Hopf bifurcation

The occurrence of Hopf-Hopf bifurcations is stated in Theorem 2 and Corollary 3; we will prove both here.

Proof of Theorem 2.

Throughout the proof we assume that $l_{1}<l_{2}$, without loss of generality.

Suppose that at a certain parameter value $F$ a Hopf-Hopf bifurcation occurs, that means, the $l_{1}$-th and $l_{2}$-th eigenvalue pairs both have real part equal to 0. So, we need to have that $F_{\operatorname{H}}(l_{1},n)=F_{\operatorname{H}}(l_{2},n)$, or, equivalently,

[TABLE]

where

[TABLE]

(The second condition follows from the fact that if $l_{1}$ and $l_{2}$ give the same value for $f$, then for the continuous function $\tilde{f}$ we need $y_{1}=\tfrac{2\pi l_{1}}{n}$ to be left and $y_{2}=\tfrac{2\pi l_{2}}{n}$ right of the top $y_{\text{top}}=\cos^{-1}\left(\tfrac{1}{4}\right)$ in the domain of consideration, $(0,\pi)$. So, $y_{1}$ and $y_{2}$ have to satisfy $0<y_{1}<y_{\text{top}}<y_{2}<\tfrac{2\pi}{3}$ (see Figure 24 for the picture). This is equivalent to equation (22).)

Since $f(l,n)$ can be written as

[TABLE]

the substitution $x=\cos\tfrac{2\pi l}{n}$ gives the function

[TABLE]

Condition (21) then becomes

[TABLE]

By condition (22) on $l_{1}$ and $l_{2}$, $\cos\tfrac{2\pi l_{1}}{n}$ is on the left and $\cos\tfrac{2\pi l_{2}}{n}$ is on the right of the maximum $x=\tfrac{1}{4}$ of $h$. Since the function $h(x)$ is symmetric around the maximum, $l_{1},l_{2}$ and $n$ should satisfy

[TABLE]

This provides the condition for a Hopf-Hopf bifurcation to occur.

Conversely, suppose that equation (7) holds for $l_{1}$ and $l_{2}$ satisfying $0<l_{1},l_{2}<\tfrac{n}{2}$, $l_{1},l_{2}\neq\tfrac{n}{3}$ and $l_{1}\neq l_{2}$. The existence of a Hopf-Hopf bifurcation can be derived from Lemma 1 by showing that the bifurcation parameters $F_{\operatorname{H}}$ corresponding to each of these eigenvalue pairs coincide, i.e. $F_{\operatorname{H}}(l_{1},n)=F_{\operatorname{H}}(l_{2},n)$.

Let us denote $y_{1}=\tfrac{2\pi l_{1}}{n}$ and $y_{2}=\tfrac{2\pi l_{2}}{n}$, then equation (7) becomes:

[TABLE]

From this equation, we obtain $\cos y_{2}=\tfrac{1}{2}-\cos y_{1}$ and, with a little trigonometry, its double angle reads

[TABLE]

Now, observe that the following holds:

[TABLE]

Hence, it holds that $f(l_{1},n)=\tilde{f}(y_{1})=\tilde{f}(y_{2})=f(l_{2},n)$ and therefore $F_{\operatorname{H}}(l_{1},n)=F_{\operatorname{H}}(l_{2},n)$ as desired. ∎

Proof of Corollary 3.

From equation (7) we can determine explicit combinations of $l_{1},l_{2}$ and $n$ for which a Hopf-Hopf bifurcation will occur. To begin with the easiest one, criterion 2: choose $l_{2}/n$ such that $\cos\tfrac{2\pi l_{2}}{n}=0$, i.e. $l_{2}/n=1/4$. This implies that $l_{1}/n$ has to be equal to $1/6$ to satisfy equation (7). Since all numbers have to be integers, we should take $n=12m$ with $m\in\mathbb{N}$ and hence, $l_{1}=2m$ and $l_{2}=3m$.

Criterion 1 is obtained by observing that

[TABLE]

so that we have the relations $2l_{1}/n=1/5$ and $2l_{2}/n=3/5$. These are satisfied by taking multiples of $m\in\mathbb{N}$ as follows: $n=10m$ and $l_{1}=~{}m,l_{2}=3m$. ∎

4.3 Proof of Theorem 4

To prove the occurrence of a Hopf bifurcation we need to show that (a) there is an eigenvalue pair crossing the imaginary axis; and (b) the first Lyapunov coefficient $\ell_{1}(l,n)$ is nonzero at $F_{\operatorname{H}}$ [Kuznetsov04]. Lemma 1 shows that the transversality condition (a) holds. Below, we prove the nondegeneracy condition (b), while we also clarify the condition for $\ell_{1}(l,n)$ to be positive or negative. In the following, let $l$ and $n$ be as in Lemma 1 and assume that there is no $l_{2}\neq l$ which satisfies both Lemma 1 and equation (7) (with $l_{1}=l$).

The first Lyapunov coefficient $\ell_{1}$ corresponding to a Hopf bifurcation of the equilibrium $x_{F}$ for the $l$-th eigenvalue pair is given by the following invariant expression [Kuznetsov04]:

[TABLE]

where $A$ is the Jacobian matrix and $B$ and $C$ are multilinear functions obtained via the Taylor expansion of the nonlinear part of (1a). The vectors $q$ and $p$ are complex eigenvectors of $A$ and $A^{\top}$, respectively, and have to be taken such that $q$ is associated with the eigenvalue which crosses the imaginary axis at $F_{\operatorname{H}}(l,n)$ and which has positive imaginary part $\omega_{0}$, while $p$ is its adjoint eigenvector. In other words, $q$ and $p$ have to be eigenvectors corresponding to $\lambda_{n-l}$ and $\bar{\lambda}_{n-l}$, respectively. We will specify them later on. Furthermore, note that the inner product on $\mathbb{C}^{n}$ is defined such that it is antilinear in the first component, i.e.

[TABLE]

In this section we will simplify formula (4.3) to an analytic expression which depends on the variables $l$ and $n$ only and whose sign is easily determined. We do this by taking advantage of the fact that the Jacobian matrix at $x_{F}$ is circulant. The first step is to simplify the expression as much as possible in a general setting. Next, we plug in all the known terms specific to our system. In the last part we determine the condition to have either a positive or a negative Lyapunov coefficient, proving the remaining part of Theorem 4.

4.3.1 Simplifying the expression

First of all, note that by a change of coordinates, $y_{j}=x_{j}-F$ (which translates the equilibrium $x_{F}$ to the origin), we can write the Lorenz-96 model (1) in the following form:

[TABLE]

where $A$ is the $n\times n$ Jacobian matrix at the origin and $B:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a bilinear map whose $k$-th component is given by

[TABLE]

Since the cubic terms are absent in this model, we can immediately simplify (4.3) to

[TABLE]

We split this equation into two components as follows:

[TABLE]

where

[TABLE]

In the Lorenz-96 case the matrix $A$ is circulant for the trivial equilibrium $x_{F}$ and therefore unitarily equivalent with a diagonal matrix $D$ [Gray06]:

[TABLE]

Here, the diagonal entries of $D$ are the eigenvalues of $A$ and the columns of $U$ are the eigenvectors of $A$ in the same order. Moreover, since the matrix $A$ is real, we have that $A^{\top}=UD^{*}U^{*}$ which means that

[TABLE]

Since $q$ and $p$ are eigenvectors corresponding to the $(n-l)$-th eigenvalue of $A$ and $A^{\top}$, respectively, this means that at the eigenvalue crossing we have

[TABLE]

so that we can take

[TABLE]

the normalized eigenvector (5) of $A$ corresponding to the $(n-l)$-th eigenvalue. This fact already shows that the only values we need to compute formula (25) are $F_{\operatorname{H}},\rho_{j}$, and $\omega_{0}$, which are all determined by the choice of $l$ and $n$.

Elimination of inverse matrices

The next step is to remove the inverse matrices in formula (25). The fact that the eigenvectors of $A$ form a unitary matrix implies that we can express any $x\in\mathbb{C}^{n}$ in terms of the eigenvectors of $A$ using a standard Fourier decomposition

[TABLE]

This makes it easy to determine how $A$ and its inverses act on any vector $x$:

[TABLE]

where we use the relation $A^{-1}=UD^{-1}U^{*}$ in the second line. In the following we will implement these relations in both equations (26) and (27). Note that up to this step, we only used the property that $A$ is normal.

First component $\ell_{1,a}$

By the bilinearity of the operator $B$, the linearity of the inner product in the second component and the expression for the inverse of $A$, the first part of the first Lyapunov coefficient (26) can be written as

[TABLE]

In the Lorenz-96 case, the inner product terms in the last line become

[TABLE]

by equation (24).

We can fill in the explicit relations from equations (5) and (28), i.e. $v_{j}^{k}=\rho_{j}^{k}/\sqrt{n}$ and $p_{k}=q_{k}=v_{n-l}^{k}$, to obtain

[TABLE]

Note that the sums at the end of both equations are conjugate to each other. For $j\neq 0\bmod n$, we can use a result from finite geometric series and the fact that $\rho_{j}^{n}=1$ to compute

[TABLE]

For $j=0$, we have that $\rho_{0}=1$ and so the sum becomes

[TABLE]

In brief, we found that

[TABLE]

Thus, equations (30) are equal to 0 for all $j$ but 0. Hence, only the term for $j=0$ is left in the summation for $\ell_{1,a}$, which then reduces to

[TABLE]

In the computation of formula (25) we only need the real part of the last expression. It mainly consists of powers of $\rho_{l}$, so Euler’s formula yields

[TABLE]

Second component $\ell_{1,b}$

The second part can be simplified similarly. Using the bilinearity of the operator $B$, the linearity of the inner product in the second component and the expression for the inverse matrix, the second part of the first Lyapunov coefficient (27) is given by

[TABLE]

In the case of the Lorenz-96 model, each of the inner product parts can be written as

[TABLE]

by a small adjustment of equations (29).

As before, we can replace $p,q$ and $v$ by powers of $\rho$ times a constant:

[TABLE]

Note that the sums in both equations are again conjugate to each other. The summand of the sum in the second equation can be written as $\rho_{j}^{k}\rho_{l}^{2k}=\rho_{j+2l}^{k}$. Formula (31), with $j$ replaced by $j+2l$, then shows that

[TABLE]

The sum in the first equation of (33) gives exactly the same result, by conjugacy. Therefore, in both cases only the terms with $j=n-2l$ are nonzero:

[TABLE]

These results reduce $\ell_{1,b}$ to

[TABLE]

Note that — as with the first component — the summation and indices $j$ disappeared.

Again, we need the real part of (4.3.1) only. This is a bit more complicated than in the case of $\ell_{1,a}$. First of all, we can reduce the numerator to

[TABLE]

We can take the real part of the complex denominator easily via {\rm Re}\big{(}\tfrac{1}{a+bi}\big{)}=\tfrac{a}{a^{2}+b^{2}}. Thus we find, using the expressions (18) and (19),

[TABLE]

Finally, by substituting these intermediate results in equation (4.3.1), the real part of the second component $\ell_{1,b}$ results as

[TABLE]

4.3.2 Sign of the first Lyapunov coefficient

Observe that both main components ${\rm Re}\ \ell_{1,a}$ and ${\rm Re}\ \ell_{1,b}$ only depend on $l$ and $n$. So, combining equations (4.3.1) and (35) gives an expression of the first Lyapunov coefficient (25) merely in terms of $l$ and $n$:

[TABLE]

where $l$ should be taken such that $0<l<\tfrac{n}{2},l\neq\tfrac{n}{3}$. It is easy to see already that $\ell_{1}=0$ if we would choose $l=0$ or $\tfrac{n}{3}$ and, moreover, that for fixed $n$ we have $\lim_{l\rightarrow n/2}\ell_{1}(l,n)=-\infty$, by the tangent function. Figure 25 shows these properties in the (continuous) graph of formula (4.3.2).

As it turns out, equation (4.3.2) is a useful and easy way to compute the first Lyapunov coefficient. In order to conclude whether the bifurcation is sub- or supercritical, we need to show which combinations of $l$ and $n$ yield a positive or negative value of $\ell_{1}(l,n)$.

Firstly, observe that the factors in front of the quotient in (4.3.2) are always positive, either by the square or by the fact that the tangent function is positive for $0<\tfrac{l}{n}<\tfrac{1}{2}$. Secondly, it is easy to check that the denominator of the big quotient of (4.3.2) is positive on the entire domain. It remains to determine where the numerator of the quotient is positive or negative. Let

[TABLE]

be the numerator defined as a continuous function in $y\in[0,\pi]$, where we replaced $\tfrac{2\pi l}{n}$ by the variable $y$, see Figure 26. Its derivative satisfies

[TABLE]

It is easy to see that $L^{\prime}(y)$ has a zero if $\sin y=0$ or if the second order polynomial $24x^{2}-32x-11$, obtained by the substitution $x=\cos y$, has a zero. In the first case, $\sin y=0$, we have zeroes at $y=0$ and $y=\pi$, which give a global and a local maximum of $L(y)$. The second case provides us with exactly two zeroes of the polynomial, namely $x_{\pm}=\tfrac{2}{3}\pm\sqrt{\tfrac{65}{72}}$. However, only $x_{-}$ is a true solution of $L^{\prime}(y)$, since the value $x_{+}$ lies outside the range of the cosine function. Hence, y_{\min}=\arccos\Big{(}\tfrac{2}{3}-\sqrt{\tfrac{65}{72}}\Big{)} gives the global minimum of $L(y)$ (see Figure 26). These three solutions are the only zeroes of $L^{\prime}(y)$ on the domain $[0,\pi]$. It follows that the sign of $L^{\prime}(y)$ on each of the intervals $(0,y_{\min})$ and $(y_{\min},\pi)$ does not change.

In order to determine the zeroes of $L(y)$, let us consider first the interval $(0,y_{\min})$. Observe that the derivative $L^{\prime}$ is negative on this interval (because $L^{\prime}\left(\tfrac{\pi}{2}\right)=-11$, for example). Since $L(0)=3$, this means that $L$ can have at most one zero on the interval $(0,y_{\min})$. Likewise, $L^{\prime}(y)>0$ for all $y\in(y_{\min},\pi)$ (by the fact that $L^{\prime}\left(\tfrac{2\pi}{3}\right)=\tfrac{11}{2}\sqrt{3}$). Since $L(\pi)=-3$, this implies that $L$ has no zero on $(y_{\min},\pi)$, while it should have at least one zero on the entire interval $(0,\pi)$. Consequently, there exists exactly one $y_{0}\in(0,y_{\min})$ such that $L(y_{0})=0$, namely $y_{0}\approx 0.5545380$ and, moreover, we have that

[TABLE]

This can also be seen from Figure 26.

To conclude: the simplified expression (4.3.2) for the first Lyapunov coefficient $\ell_{1}$ is positive for $\tfrac{l}{n}<\tfrac{y_{0}}{2\pi}$ and negative for $\tfrac{l}{n}>\tfrac{y_{0}}{2\pi}$, where $\tfrac{y_{0}}{2\pi}\approx 0.08825746$. Therefore, the first Lyapunov coefficient itself is positive (hence, a subcritical Hopf bifurcation) for any $\tfrac{l}{n}\in(0,\tfrac{y_{0}}{2\pi})$ and negative (hence, a supercritical Hopf bifurcation) for $\tfrac{l}{n}\in\left(\tfrac{y_{0}}{2\pi},\tfrac{1}{2}\right)\setminus\{\tfrac{1}{3}\}$ (see Figure 25). ∎

4.4 Proof of Proposition 5

Lemma 1 implies that the trivial equilibrium undergoes a Hopf bifurcation at the parameter value $F_{\operatorname{H}}(l,n)=1/f(l,n)$. The first Hopf bifurcation for $F>0$ takes place for the integer $l\in(0,\tfrac{n}{3})$ that minimizes the value of $F_{\operatorname{H}}(l,n)$, which is equivalent to maximizing $f(l,n)$. (If $l\in(\tfrac{n}{3},\tfrac{n}{2})$, then we obtain a negative value $F_{\operatorname{H}}$.)

For all $n\geq 4$ except $n=7$ there exists at least one integer $l\in[\tfrac{n}{6},\tfrac{n}{4}]$. Indeed, for $n=4,5$ and $6$ this follows by simply taking $l=1$, and for $n=8,9,10$ and $11$ it follows by taking $l=2$. For $n\geq 12$ the width of the interval is larger than 1. We now claim that this implies that

[TABLE]

as well. To that end we use $\tilde{f}(y)=\cos y-\cos 2y$. Note that $y\in[\tfrac{\pi}{3},\tfrac{\pi}{2}]$ implies that $\tilde{f}(y)\geq 1$ and $y\in(0,\tfrac{\pi}{3})\cup(\tfrac{\pi}{2},\tfrac{2\pi}{3})$ implies that $0<\tilde{f}(y)<1$. Moreover, $l\in[\tfrac{n}{6},\tfrac{n}{4}]$ implies that $\frac{2\pi l}{n}\in[\tfrac{\pi}{3},\tfrac{\pi}{2}]$. Therefore, $f(l,n)$ is maximized for some integer $l\in[\tfrac{n}{6},\tfrac{n}{4}]$.

In case $n=7$, we can easily compute the smallest value $F_{\operatorname{H}}(l,n)$ for which a Hopf bifurcation occurs. We have shown in the proof of Lemma 1 (see Section 4.1) that this is the case for $l=1$.

Finally, assume that the bifurcation is a Hopf bifurcation, i.e., only one eigenvalue pair crosses the imaginary axis. Since $l_{1}(n)/n\in[\tfrac{1}{7},\tfrac{1}{4}]$ it follows immediately from Theorem 4 that the first Lyapunov coefficient is negative, which means that the bifurcation is supercritical. ∎

4.5 Proofs for the two-parameter system

For the two-parameter system (9) we have stated the existence of trapping regions in Proposition 7 and the occurrence of Hopf bifurcations in Lemma 8. These two results are proven below.

Proof of Proposition 7.

Recall that $E=\tfrac{1}{2}\sum_{j}x_{j}^{2}$ and let $V^{2}=2E$, which means that $V=\|x\|$, using the Euclidean norm. Fix $n\in\mathbb{N}$ and $F\in\mathbb{R}$. Take the time derivative of $E$ to obtain

[TABLE]

using the standard inner product on $\mathbb{R}^{n}$ with $A$ the matrix equal to the Jacobian at the origin. This matrix is circulant and symmetric with first row equal to

[TABLE]

Therefore, all eigenvalues of $A$ are real and given by [Gray06]:

[TABLE]

To estimate the first term in the right-hand-side of equation (4.5), we use that the Rayleigh quotient satisfies

[TABLE]

Therefore, we have to find the largest eigenvalue $\lambda^{A}_{\max}$ of $A$. It is easy to check that for negative values of $G$ the largest eigenvalue occurs for $j=\tfrac{n}{2}$ (or the nearest integer). For non-negative values of $G$, the largest eigenvalue is always equal to $\lambda^{A}_{0}=-1$.

The second term of equation (4.5) can be estimated by the Cauchy-Schwarz inequality. Equation (4.5) then becomes

[TABLE]

which can be simplified to

[TABLE]

Since $V$ is defined as the norm of $x$, $\lambda^{A}_{\max}$ needs to be negative in order to obtain a solution for $V$ such that $\frac{\textrm{d}V}{\textrm{d}t}<0$. This is certainly the case when $G\geq 0$. For negative $G$ and even $n$ the largest eigenvalue is obtained by $\lambda^{A}_{\max}=\lambda^{A}_{n/2}=-1-4G$, which is less than 0 if $G>-\tfrac{1}{4}$ for any $n$. For odd values of $n$ we may allow slightly smaller values of $G$, which makes $G=-\tfrac{1}{4}$ a suitable upper bound for the range of $G$ without trapping region.

It follows that $\tfrac{\textrm{d}V}{\textrm{d}t}<0$, whenever $G>-\tfrac{1}{4}$ and whenever $V$ is larger than the radius

[TABLE]

Therefore, under these conditions, each sphere with radius $r>R(n)$ is a trapping region for dimension $n$. In particular, if $G=0$ then $R=\sqrt{n}|F|$. ∎

Proof of Lemma 8.

Let $n$ and $l$ be as given. We choose $F$ as the bifurcation parameter. In order to have a Hopf bifurcation for the $l$-th eigenvalue pair $\{\kappa_{l},\kappa_{n-l}\}(F,G,n)$, the real part $\mu_{l}:=\operatorname{Re}\kappa_{l}$ has to vanish. This occurs if $F$ equals

[TABLE]

For these parameter values we have a purely imaginary eigenvalue pair with the absolute value of the imaginary part given by

[TABLE]

Note that $\omega_{0}=0$ if $l=\tfrac{n}{2}$ or if $G=-\tfrac{1}{2}(1-\cos\tfrac{2\pi l}{n})^{-1}$ for some $l$. The last condition is equivalent to $F_{\operatorname{H}}=0$ by formula (39) and implies even that $\kappa_{l}=0$. Therefore, this parameter value needs to be excluded.

Furthermore, the eigenvalue pair crosses the imaginary axis with nonzero speed, by the fact that for our restriction of $l$ we have

[TABLE]

where the derivative is with respect to the bifurcation parameter $F$. (If we would have taken $G$ as bifurcation parameter, then $\mu_{l}^{\prime}$ will be nonzero as well.)

Equation (39) thus gives us for general $n$ and for each allowed $l$ a whole line of Hopf bifurcations, which is linear in $G$. Rewritten in terms of $F$ gives the linear curves (11). ∎

5 Conclusions and outlook

The main goal of the current study was to investigate the dynamical features of the Lorenz-96 model and subsequently prove their existence by analytical means. In this partial inventory, we have proven the existence of Hopf and Hopf-Hopf bifurcations for the trivial equilibrium in all dimensions $n\geq 4$. It is also shown that if the first bifurcation for $F>0$ is a Hopf bifurcation, then it is supercritical by providing an exact formula for the first Lyapunov coefficient, which holds for all possible $n$ and $l$. Our analytical results coincide with MatCont’s numerical estimates of the parameter value $F_{\operatorname{H}}$ and (up to a scaling) the value of the first Lyapunov coefficient. Furthermore, the periodic orbits born at the first Hopf bifurcation have the physical interpretation of travelling waves.

We have proved a necessary and sufficient condition for the occurrence of a Hopf-Hopf bifurcation in the original Lorenz-96 model. To unfold these Hopf-Hopf bifurcations we introduced an extra parameter $G$ via a Laplace-like diffusion term. For this particular unfolding the Hopf bifurcations of the trivial equilibrium are given by straight lines in the $(F,G)$-plane and their intersections give Hopf-Hopf bifurcations. In the special case of dimension $n=12$ a Hopf-Hopf bifurcation lies on the line $G=0$ and is in fact the first bifurcation through which the trivial equilibrium bifurcates for $F>0$. We have shown that this codimension two bifurcation point acts as an organising centre. For dimension $n=36$ we find a Hopf-Hopf bifurcation close to the line $G=0$. In both cases normal form analysis shows that two periodic attractors coexist in a region of the $(F,G)$-plane between two Neimark-Sacker bifurcation curves. This region actually intersects the line $G=0$, which means that multistability also occurs in the original Lorenz-96 model. We expect that this mechanism leading to multistability can be found in other dimensions as well. Therefore, by adding a new parameter to the Lorenz-96 model the dynamics observed in the original model can be explained better.

Finally, we have numerically investigated the dynamics of the model for dimensions up to $n=100$ and parameter values $F$ beyond the first Hopf bifurcation value. In dimension $n=4$, a periodic attractor disappears through a fold bifurcation. After this fold bifurcation intermittency is detected, which is possibly explained through a nearby heteroclinic cycle between four equilibria. For dimensions up to $n=100$ that are multiples of $5$, a pattern is discovered with a persisting period-doubling bifurcation. However, we do not expect this pattern to persist for high dimensions. For general $n$, the routes to chaos are numerous and can comprise intermittent transitions, period-doubling cascades and possibly Newhouse-Ruelle-Takens scenarios.

Based on our results we can draw three major conclusions. Firstly, contrary to the persistence of Hopf and Hopf-Hopf bifurcations for any dimension $n\geq 4$, no clear pattern on bifurcations of periodic attractors is found. Secondly, the dependence of the dynamics on $n$ shows the importance of choosing appropriate values of the parameters in specific applications of the Lorenz-96 model such as those listed in Table 1. Lastly, the observation that the wave number of the travelling waves increases with $n$ indicates that the Lorenz-96 cannot be interpreted as a discretised PDE model.

Despite the lack of a clear bifurcation pattern for all dimensions $n$, the Lorenz-96 remains an interesting model to study for its rich dynamics. There are several open questions that have not been addressed in the present paper. For example, how can chaotic attractors with more than one positive Lyapunov exponent arise? Or, how typical is the multistability property in higher dimensions and how does it influence the dynamics? What dynamics can be expected for $F<0$? What is the role of symmetries in the Lorenz-96 model? We are therefore continuing to investigate this model on its dynamical behaviour using both analytical and numerical methods.

Acknowledgements

We would like to thank the two anonymous reviewers for their useful comments and suggestions that have helped to improve this paper.