THE FOLD-FLIP BIFURCATION

A. Yu. V. Kuznetsov, H. G. E. Meijer, L. van Veen

► To cite this version:

A. Yu. V. Kuznetsov, H. G. E. Meijer, L. van Veen. THE FOLD-FLIP BIFURCATION. International journal of bifurcation and chaos in applied sciences and engineering, 2004, 14 (7), pp.2253-2282. 10.1142/S0218127404010576 . hal-01352127

HAL Id: hal-01352127

https://hal.science/hal-01352127

Submitted on 5 Aug 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# THE FOLD-FLIP BIFURCATION

YU. A. KUZNETSOV and H. G. E. MEIJER

Mathematical Institute, University of Utrecht,
P.O. Box 80.010, 3508 TA Utrecht, the Netherlands

L. VAN VEEN

Department of Mechanical Engineering,
Graduate School of Engineering, Kyoto University,
Yoshida-Honmachi Sakyo-ku, Kyoto 606-8501, Japan

The fold-flip bifurcation occurs if a map has a fixed point with multipliers +1 and -1 simultaneously. In this paper the normal form of this singularity is calculated explicitly. Both local and global bifurcations of the unfolding are analyzed by exploring a close relationship between the derived normal form and the truncated amplitude system for the fold-Hopf bifurcation of ODEs. Two examples are presented, the generalized Hénon map and an extension of the Lorenz-84 model. In the latter example the first-, second- and third-order derivatives of the Poincaré map are computed using variational equations to find the normal form coefficients.

Keywords: Bifurcations of fixed points; normal forms; center manifold.

## 1. Introduction

Consider a family of discrete-time dynamical systems generated by iteration of

$$x \mapsto F(x, \alpha), \quad (1)$$

where map $F: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$ is sufficiently smooth.

It is well known (see e.g. [Arnold, 1983] for the general theory and [Kuznetsov, 1998] for computational formulas) that in generic one-parameter families (1) only the following three bifurcations of fixed points happen:

1. **Fold:** The fixed point has a simple eigenvalue $\lambda_1=1$ and no other eigenvalues on the unit circle, while the restriction of (1) to a one-dimensional center manifold at the critical parameter value has the form

$$\xi \mapsto \xi + \frac{1}{2} a \xi^2 + O(\xi^3), \quad (2)$$

where $a \neq 0$. When the parameter crosses the critical value, two fixed points coalesce

and disappear. This bifurcation is often called the saddle-node bifurcation. If $Av = F_xv$ and $B(u, v) = F_{xx}[u, v]$ are evaluated at the critical fixed point, then

$$a = \langle q^*, B(q, q) \rangle, \quad (3)$$

where $Aq = q$, $A^Tq^* = q^*$, and $\langle q^*, q \rangle = 1$. Here and in what follows $\langle u, v \rangle = \bar{u}^T v$ is the standard scalar product in $\mathbb{C}^n$ (or $\mathbb{R}^n$).

2. **Flip:** The fixed point has a simple eigenvalue $\lambda_1 = -1$ and no other eigenvalues on the unit circle, while the restriction of (1) to a one-dimensional center manifold at the critical parameter value can be transformed to the normal form

$$\xi \mapsto -\xi + \frac{1}{6} b \xi^3 + O(\xi^4), \quad (4)$$

where $b \neq 0$. When the parameter crosses the critical value, a cycle of period 2 bifurcates from the fixed point. This phenomenon is often called the period-doubling bifurcation. If $C(u, v, w) =$
$F_{xxx}[u, v, w]$ is evaluated at the critical fixed point, then

$$b = \langle p^*, C(p, p, p) + 3B(p, (I_n - A)^{-1}B(p, p)) \rangle, \quad (5)$$

where $I_n$ is the unit $n \times n$ matrix, $Ap = -p$, $A^T p^* = -p^*$, and $\langle p^*, p \rangle = 1$.

2. Neimark-Sacker: The fixed point has simple critical eigenvalues $\lambda_{1,2} = e^{\pm i\theta_0}$ and no other eigenvalues on the unit circle. Assume that

$$e^{iq\theta_0} - 1 \neq 0, \\ q = 1, 2, 3, 4 \text{ (no strong resonances).}$$

Then, the restriction of (1) to a two-dimensional center manifold at the critical parameter value can be transformed to the normal form

$$\eta \mapsto \eta e^{i\theta_0} \left( 1 + \frac{1}{2} d|\eta|^2 \right) + O(|\eta|^4), \quad (6)$$

where $\eta$ is a complex variable and $d$ is a complex number. Further assume that

$$c = \mathrm{Re} d \neq 0.$$

Under the above assumptions, a unique closed invariant curve around the fixed point appears when the parameter crosses the critical value. One has the following expression for $d$:

$$d = e^{-i\theta_0} \langle v^*, C(v, v, \bar{v}) + 2B(v, (I_n - A)^{-1}B(v, \bar{v})) + B(\bar{v}, (e^{2i\theta_0}I_n - A)^{-1}B(v, \bar{v})) \rangle, \quad (7)$$

where $Av = e^{i\theta_0}v$, $A^T v^* = e^{-i\theta_0}v^*$, and $\langle v^*, v \rangle = 1$.

In generic two-parameter families, the above listed codim 1 bifurcations occur when we cross the corresponding bifurcation curve (defined by the single condition on the eigenvalues) at a typical point. Moreover, one can choose the eigenvectors so that the normal form coefficients $a, b$ and $c$ will be smooth functions along these curves.<fnref n="1" /> While tracing a bifurcation curve in a generic two-parameter family (1), one may encounter a double-degenerate singularity if either (i) extra eigenvalues approach the unit circle, thus changing the dimension of the center manifold; or (ii) one of the nonequalities on the normal form coefficients mentioned above becomes an equality. Therefore, the following eleven doubly degenerate (or codim 2) bifurcation points

can be met in generic two-parameter families of maps:

$$\begin{align*}
D_1: \lambda_1 &= 1, \quad a = 0 \text{ (cusp)}; \\
D_2: \lambda_1 &= -1, \quad b = 0 \text{ (generalized flip)}; \\
D_3: \lambda_{1,2} &= e^{\pm i\theta_0}, \quad c = 0 \text{ (Chenciner point)}; \\
D_4: \lambda_1 &= \lambda_2 = 1 \text{ (1:1 resonance)}; \\
D_5: \lambda_1 &= \lambda_2 = -1 \text{ (1:2 resonance)}; \\
D_6: \lambda_{1,2} &= e^{\pm i\theta_0}, \quad \theta_0 = 2\pi/3 \text{ (1:3 resonance)}; \\
D_7: \lambda_{1,2} &= e^{\pm i\theta_0}, \quad \theta_0 = \pi/2 \text{ (1:4 resonance)}; \\
D_8: \lambda_1 &= 1, \quad \lambda_2 = -1; \\
D_9: \lambda_1 &= 1, \quad \lambda_{2,3} = e^{\pm i\theta_0}; \\
D_{10}: \lambda_1 &= -1, \quad \lambda_{2,3} = e^{\pm i\theta_0}; \\
D_{11}: \lambda_{1,2} &= e^{\pm i\theta_0}, \quad \lambda_{3,4} = e^{\pm i\theta_1}.
\end{align*}$$

The cases $D_1-D_7$ are well understood and presented in many textbooks (see, e.g. [Arnold, 1983; Arrowsmith & Place, 1990; Kuznetsov, 1998]). Notice that these are the only possibilities if the map (1) is a Poincaré map associated with a limit cycle in a three-dimensional autonomous system of ODEs: Then there are just two eigenvalues $\lambda_{1,2}$ and their product is always positive due to the Liouville Theorem, i.e. $\lambda_1\lambda_2 > 0$. This is why cases $D_8-D_{11}$ received much less attention.

This paper is devoted to a detailed analysis of case $D_8$, that we call the fold-flip bifurcation. This bifurcation has been first treated by [Gheiner, 1994]. In the present paper we clarify, correct, and extend that analysis. In Sec. 2, we derive the parameter-dependent normal form for a generic fold-flip bifurcation and express explicitly the critical normal form coefficients in two- and n-dimensional systems. In Sec. 3 we complete the analysis of local bifurcations of the normal form and study global bifurcations approximating the map by time-shifts along orbits of an auxiliary planar ODE, present (as complete as possible) bifurcation diagrams of the truncated normal form, including computer-generated phase portraits, and then discuss effects of the truncation. Section 4 is devoted to the normal form analysis of the fold-flip bifurcation in two examples: A generalized Hénon map and the Poincaré map for a four-dimensional extension of the Lorenz-84 model. The generalized Hénon map appears in the study of

<footnote n="1">¹The coefficient $\theta_0$ also has this property.</footnote>
bifurcations of diffeomorphisms with codim 2 homoclinic and heteroclinic tangencies, while the second model describes the atmospheric circulation and, to our best knowledge, is the first autonomous ODE system appearing in applications that demonstrates the fold-flip bifurcation. For this model, we compute numerically the partial derivatives of the Poincaré map (see Appendix) and then use them to evaluate the critical normal form coefficients and, thus,

demonstrate the nondegeneracy of the fold-flip in this case.

## 2. Normal Form for the Fold-Flip Bifurcation

### 2.1. Planar normal form

**Proposition 2.1.1** (Critical normal form). Suppose a smooth map $F_0: \mathbb{R}^2 \to \mathbb{R}^2$ has the form

$$ \begin{pmatrix} \xi_1 \\ \xi_2 \end{pmatrix} \mapsto \begin{pmatrix} \xi_1 + \sum_{i+j=2,3} \frac{1}{i!j!} g_{ij} \xi_1^i \xi_2^j \\ -\xi_2 + \sum_{i+j=2,3} \frac{1}{i!j!} h_{ij} \xi_1^i \xi_2^j \end{pmatrix} + O(\|\xi\|^4) \quad (8) $$

and $h_{11} \neq 0$. Then $F_0$ is smoothly equivalent near the origin to a map

$$ \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \mapsto \begin{pmatrix} x_1 + \frac{1}{2}a(0)x_1^2 + \frac{1}{2}b(0)x_2^2 + \frac{1}{6}c(0)x_1^3 + \frac{1}{2}d(0)x_1x_2^2 \\ -x_2 + x_1x_2 \end{pmatrix} + O(\|x\|^4), \quad (9) $$

where

$$ \begin{aligned} a(0) &= \frac{g_{20}}{h_{11}}, & b(0) &= g_{02}h_{11}, \\ c(0) &= \frac{1}{h_{11}^2} \left( g_{30} + \frac{3}{2}g_{11}h_{20} \right), & & \end{aligned} \quad (10) $$

$$ \begin{aligned} d(0) &= \frac{3g_{02}(h_{02}h_{20}+2h_{21}-2g_{11}h_{20})-g_{20}(3h_{02}^2+2h_{03})}{6h_{11}} \\ & \quad -g_{11}^2+g_{12}+\frac{1}{2}g_{11}h_{02}-h_{02}^2-\frac{2}{3}h_{03}. & & \end{aligned} \quad (11) $$

*Proof*

**Step 1** (Quadratic terms). Applying to (8) a polynomial coordinate transformation

$$ \left\{ \begin{aligned} \xi_1 &= x_1 + \frac{1}{2}G_{20}x_1^2 + G_{11}x_1x_2 + \frac{1}{2}G_{02}x_2^2, \\ \xi_2 &= x_2 + \frac{1}{2}H_{20}x_1^2 + H_{11}x_1x_2 + \frac{1}{2}H_{02}x_2^2, \end{aligned} \right. \quad (12) $$

we obtain:

$$ x_1 \mapsto x_1 + \frac{1}{2}g_{20}x_1^2 + (g_{11} + 2G_{11})x_1x_2 + \frac{1}{2}g_{02}x_2^2 + \dots, $$

$$ x_2 \mapsto -x_2 + \frac{1}{2}(h_{20} - 2H_{20})x_1^2 + h_{11}x_1x_2 + \frac{1}{2}(h_{02} - 2H_{02})x_2^2 + \dots, $$

where dots stand for higher-order terms. By setting

$$ G_{11} = -\frac{1}{2}g_{11}, \quad H_{20} = \frac{1}{2}h_{20}, \quad H_{02} = \frac{1}{2}h_{02}, \quad (13) $$

we eliminate as much quadratic terms as possible. The remaining quadratic terms are called resonant.

**Step 2** (Cubic terms). Assume now that Step 1 is already done, so that (8) has only resonant quadratic and all cubic terms. Consider a polynomial transformation:

$$ \left\{ \begin{aligned} \xi_1 &= x_1 + \frac{1}{6}G_{30}x_1^3 + \frac{1}{2}G_{21}x_1^2x_2 + \frac{1}{2}G_{12}x_1x_2^2 \\ & \quad + \frac{1}{6}G_{03}x_2^3, \\ \xi_2 &= x_2 + \frac{1}{6}H_{30}x_1^3 + \frac{1}{2}H_{21}x_1^2x_2 + \frac{1}{2}H_{12}x_1x_2^2 \\ & \quad + \frac{1}{6}H_{03}x_2^3, \end{aligned} \right. \quad (14) $$
Obviously, it does not change the quadratic terms. After this transformation, we get

$$x_1 \mapsto x_1 + \frac{1}{2}g_{20}x_1^2 + \frac{1}{2}g_{02}x_2^2 \\ + \frac{1}{6}g_{30}x_1^3 + \frac{1}{2}(g_{21} + 2G_{21})x_1^2x_2 \\ + \frac{1}{2}g_{12}x_1x_2^2 + \frac{1}{6}(g_{03} + 2G_{03})x_2^3 + \dots$$

and

$$x_2 \mapsto -x_2 + h_{11}x_1x_2 + \frac{1}{6}(h_{30} - 2H_{30})x_1^3 \\ + \frac{1}{2}h_{21}x_1^2x_2 + \frac{1}{2}(h_{12} - 2H_{12})x_1x_2^2 \\ + \frac{1}{6}h_{03}x_2^3 + \dots$$

By setting

$$G_{21} = -\frac{1}{2}g_{21}, \quad G_{03} = -\frac{1}{2}g_{03}, \\ H_{30} = \frac{1}{2}h_{30}, \quad H_{12} = \frac{1}{2}h_{12},$$

we eliminate four cubic terms. The remaining cubic terms are also called resonant. They are not altered by (14).

**Step 3** (More cubic terms). The coefficients $H_{11}$, $G_{20}$, and $G_{02}$ of (12) do not affect quadratic terms of (8) but alter its cubic terms. Taking into account (13) while computing the cubic terms of the transformed map, we obtain

$$x_1 \mapsto x_1 + \frac{1}{2}g_{20}x_1^2 + \frac{1}{2}g_{02}x_2^2 + \frac{1}{6}\left(g_{30} + \frac{3}{2}g_{11}h_{20}\right)x_1^3 \\ + \frac{1}{2}\left(2g_{02}H_{11} - g_{02}G_{20} + (g_{20} + 2h_{11})G_{02} + \frac{1}{2}g_{11}h_{02} + g_{12} - g_{11}^2\right)x_1x_2^2 + \dots$$

and

$$x_2 \mapsto -x_2 + h_{11}x_1x_2 \\ + \frac{1}{2}\left(g_{20}H_{11} + h_{11}G_{20} - g_{11}h_{20} + \frac{1}{2}h_{02}h_{20} + h_{21}\right)x_1^2x_2 \\ + \frac{1}{6}\left(3g_{02}H_{11} + 3G_{02}h_{11} + h_{03} + \frac{3}{2}h_{02}^2\right)x_2^3 + \dots,$$

where only the resonant cubic terms are displayed. Thus, we can try to eliminate three altered terms by selecting $H_{11}$, $G_{20}$ and $G_{02}$. This implies solving the following linear system:

$$\begin{pmatrix} 2g_{02} & -g_{02} & 2h_{11} + g_{20} \\ g_{20} & h_{11} & 0 \\ 3g_{02} & 0 & 3h_{11} \end{pmatrix} \begin{pmatrix} H_{11} \\ G_{20} \\ G_{02} \end{pmatrix} = \begin{pmatrix} -\frac{1}{2}g_{11}h_{02} - g_{12} + g_{11}^2 \\ g_{11}h_{20} - \frac{1}{2}h_{02}h_{20} - h_{21} \\ -h_{03} - \frac{3}{2}h_{02}^2 \end{pmatrix}.$$

Its matrix has zero determinant. However, using the nondegeneracy condition $h_{11} \neq 0$, we can eliminate the resonant cubic terms in the second component of the normal form. Thus, we set

$$H_{11} = 0 \qquad (15)$$

and from the above linear system obtain:

$$G_{02} = -\frac{1}{h_{11}}\left(\frac{1}{3}h_{03} + \frac{1}{2}h_{02}^2\right), \\ G_{20} = \frac{1}{h_{11}}\left(g_{11}h_{20} - h_{21} - \frac{1}{2}h_{02}h_{20}\right). \qquad (16)$$

**Step 4** (Final transformation). Transform now the original map (8) using (12) with the coefficients (13) defined in Step 1, and (15), (16) defined in Step 3. This results in a map with resonant quadratic terms, nonresonant cubic terms, and only two remaining resonant cubic terms in the first component. Transformation (14) from Step 2 allows then to eliminate all nonresonant cubic terms, while keeping unchanged all remaining quadratic and cubic resonant terms. Finally, make the linear scaling

$$x_1 \mapsto \frac{x_1}{h_{11}}$$

to put the coefficient in front of $x_1x_2$ in the second component equal to one. As a result, we obtain the expressions (10) and (11) for the critical normal form coefficients. ■

**Remark 2.1.2.** We have obtained the same critical normal form (9) as Gheiner [1994]. Moreover, we have derived explicit expressions for the critical normal form coefficients.
**Proposition 2.1.3** (Parameter-dependent normal form). Consider a smooth two-parameter family of smooth planar maps

$$ \xi \mapsto F(\xi, \alpha), \quad \xi \in \mathbb{R}^2, \quad \alpha \in \mathbb{R}^2, $$

where $F: \mathbb{R}^2 \times \mathbb{R}^2 \to \mathbb{R}^2$ is smooth and such that

1. $F_0: \mathbb{R}^2 \to \mathbb{R}^2, \xi \mapsto F_0(\xi) = F(\xi, 0)$ satisfies Proposition 2.1.1;

$$ \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \mapsto \left( \mu_1 + (1+\mu_2)x_1 + \frac{1}{2}a(\mu)x_1^2 + \frac{1}{2}b(\mu)x_2^2 + \frac{1}{6}c(\mu)x_1^3 + \frac{1}{2}d(\mu)x_1x_2^2 \right) + O(\|x\|^4), \quad (18) $$

where all coefficients are smooth functions of $\mu$ and their values at $\mu_1 = \mu_2 = 0$ are given by (10) and (11).

**Remark 2.1.4.** Map (17) can be substituted by the map

$$ \begin{pmatrix} \xi \\ \alpha \end{pmatrix} \mapsto \begin{pmatrix} F(\xi, \alpha) - \xi \\ \det(F_{\xi}(\xi, \alpha) + I_2) \\ \det(F_{\xi}(\xi, \alpha) - I_2) \end{pmatrix}. $$

*Proof of Proposition 2.1.3.* Expand $F$ in $\xi$ at $\xi = 0$ for any small $\|\alpha\|$:

$$ \xi \mapsto F(\xi, \alpha) = \gamma(\alpha) + A(\alpha)\xi + R(\xi, \alpha), $$

where $\gamma(0) = 0$ and $R(\xi, \alpha) = O(\|\xi\|^2)$. The first assumption implies the existence of two eigenvectors, $q(\alpha)$ and $p(\alpha)$ in $\mathbb{R}^2$, such that

$$ A(\alpha)q(\alpha) = \lambda_1(\alpha)q(\alpha), \quad A(\alpha)p(\alpha) = \lambda_2(\alpha)p(\alpha), $$

where $\lambda_1(0) = 1$ and $\lambda_2(0) = -1$. Note that, due to the simplicity of the eigenvalues $\pm 1$ of $A(0)$, $\lambda_1, \lambda_2$ depend smoothly on $\alpha$, and $q, p$ can also be assumed to be smooth functions of $\alpha$. Any $\xi \in \mathbb{R}^2$ can now be represented for all small $\|\alpha\|$ as

$$ \xi = \eta_1 q(\alpha) + \eta_2 p(\alpha), $$

where $\eta = (\eta_1, \eta_2)^T \in \mathbb{R}^2$. One can compute the components of $\eta$ explicitly:

$$ \eta_1 = \langle q^*(\alpha), \xi \rangle, \quad \eta_2 = \langle p^*(\alpha), \xi \rangle, $$

where

$$ A^T(\alpha)q^*(\alpha) = \lambda_1(\alpha)q^*(\alpha), $$

$$ A^T(\alpha)p^*(\alpha) = \lambda_2(\alpha)p^*(\alpha) $$

2. The map $T: \mathbb{R}^2 \times \mathbb{R}^2 \to \mathbb{R}^2 \times \mathbb{R} \times \mathbb{R}$, defined by

$$ \begin{pmatrix} \xi \\ \alpha \end{pmatrix} \mapsto T(\xi, \alpha) = \begin{pmatrix} F(\xi, \alpha) - \xi \\ \det F_{\xi}(\xi, \alpha) + 1 \\ \mathrm{Tr} F_{\xi}(\xi, \alpha) \end{pmatrix} \quad (17) $$

is regular at $(\xi, \alpha) = (0, 0)$.

Then $F$ is smoothly equivalent near the origin to a family

and $\langle q^*(\alpha), q(\alpha) \rangle = \langle p^*(\alpha), p(\alpha) \rangle = 1$. Since $\langle q^*(\alpha), p(\alpha) \rangle = \langle p^*(\alpha), q(\alpha) \rangle = 0$, the map $F$ in the $\eta$-coordinates takes the form

$$ \begin{pmatrix} \eta_1 \\ \eta_2 \end{pmatrix} \mapsto \begin{pmatrix} \sigma_1(\alpha) + \lambda_1(\alpha)\eta_1 + S_1(\eta, \alpha) \\ \sigma_2(\alpha) + \lambda_2(\alpha)\eta_2 + S_2(\eta, \alpha) \end{pmatrix}, \quad (19) $$

where

$$ \begin{pmatrix} \sigma_1(\alpha) \\ \sigma_2(\alpha) \end{pmatrix} = \begin{pmatrix} \langle q^*(\alpha), \gamma(\alpha) \rangle \\ \langle p^*(\alpha), \gamma(\alpha) \rangle \end{pmatrix}, $$

$$ \begin{pmatrix} S_1(\eta, \alpha) \\ S_2(\eta, \alpha) \end{pmatrix} = \begin{pmatrix} \langle q^*(\alpha), R(\eta_1 q(\alpha) + \eta_2 p(\alpha), \alpha) \rangle \\ \langle p^*(\alpha), R(\eta_1 q(\alpha) + \eta_2 p(\alpha), \alpha) \rangle \end{pmatrix}. $$

Expanding $S_{1,2}(\eta, \alpha)$ further, we can write (19) as

$$ \begin{pmatrix} \eta_1 \\ \eta_2 \end{pmatrix} \mapsto \begin{pmatrix} \sigma_1(\alpha) + \lambda_1(\alpha)\eta_1 + \sum_{i+j=2,3} \frac{1}{i!j!} g_{ij}(\alpha) \eta_1^i \eta_2^j \\ \sigma_2(\alpha) + \lambda_2(\alpha)\eta_2 + \sum_{i+j=2,3} \frac{1}{i!j!} h_{ij}(\alpha) \eta_1^i \eta_2^j \end{pmatrix} + O(\|\eta\|^4). \quad (20) $$

Now we wish to put (20) in the form of (18) up to and including cubic terms by means of a smoothly dependent parameter coordinate transformation. Consider the following change of variables:

$$ \left\{ \begin{aligned} \eta_1 &= x_1 + \varepsilon_0(\alpha) + \varepsilon_1(\alpha)x_2 + \frac{1}{2}G_{20}(\alpha)x_1^2 \\ &\quad + G_{11}(\alpha)x_1x_2 + \frac{1}{2}G_{02}(\alpha)x_2^2 \\ &\quad + \frac{1}{2}G_{21}(\alpha)x_1^2x_2 + \frac{1}{6}G_{03}(\alpha)x_2^3, \\ \eta_2 &= x_2 + \delta_0(\alpha) + \delta_1(\alpha)x_1 + \frac{1}{2}H_{20}(\alpha)x_1^2 \\ &\quad + \frac{1}{2}H_{02}(\alpha)x_2^2 + \frac{1}{6}H_{30}(\alpha)x_2^3 + \frac{1}{2}H_{12}(\alpha)x_1x_2^2, \end{aligned} \right. \quad (21) $$
where all coefficients are yet unknown smooth functions of $\alpha$ such that $\varepsilon_i(0) = \delta_i(0) = 0$ for $i = 0, 1$. Obviously, for $\alpha = 0$ (21) reduces to the transformation introduced in Step 4 of the

proof of Proposition 2.1.1 just before the final scaling.

Require now that the Taylor expansion of (20) in the x-coordinates up to and including cubic terms takes the form

$$ \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \mapsto \left( \mu_1(\alpha) + (1 + \mu_2(\alpha))x_1 + \frac{1}{2}A(\alpha)x_1^2 + \frac{1}{2}B(\alpha)x_2^2 + \frac{1}{6}C(\alpha)x_1^3 + \frac{1}{2}D(\alpha)x_1x_2^2 \right), -x_2 + E(\alpha)x_1x_2 $$

where, $\mu_1(0) = \mu_2(0) = 0$. After all substitutions, this requirement translates into a system of algebraic equations:

$$ Q(\varepsilon_0, \varepsilon_1, \delta_0, \delta_1, \mu_1, \mu_2, G_{20}, G_{11}, G_{02}, G_{21}, G_{03}, H_{20}, H_{02}, H_{30}, H_{12}, A, B, C, D, E) = 0, $$

where $Q: \mathbb{R}^{20} \to \mathbb{R}^{20}$ results from equating the corresponding Taylor coefficients. For the Jacobi matrix $J = DQ$ of this system evaluated at

$$ \varepsilon_0 = \varepsilon_1 = \delta_0 = \delta_1 = \mu_1 = \mu_2 = 0 $$

we have $\det(J) = -3072h_{11}^3 \neq 0$. Therefore, the Implicit Function Theorem guarantees the local existence and smoothness of the coefficients of transformation (21) as functions of $\alpha$.

Moreover, one can show that

$$ \left\{ \begin{aligned} \mu_1 &= A_1\alpha_1 + A_2\alpha_2 + O(\|\alpha\|^2), \\ \mu_2 &= \frac{1}{h_{11}(0)}\left(\begin{aligned}[t] &[(g_{11}(0)B_1 + 2A_3)h_{11}(0) \\ &-(h_{02}(0)B_1 + 2B_3)g_{20}]\alpha_1 \\ &+[(g_{11}(0)B_2 + 2A_4)h_{11}(0) \\ &-(h_{02}(0)B_2 + 2B_4)g_{20}]\alpha_2 \end{aligned}\right) + O(\|\alpha\|^2), \end{aligned} \right. \quad (22) $$

where $A_i$ and $B_i$ are defined by the following expansions of the coefficients of (20):

$$ \begin{aligned} \sigma_1(\alpha) &= A_1\alpha_1 + A_2\alpha_2 + O(\|\alpha\|^2), \\ \lambda_1(\alpha) &= 1 + A_3\alpha_1 + A_4\alpha_2 + O(\|\alpha\|^2), \\ \sigma_2(\alpha) &= B_1\alpha_1 + B_2\alpha_2 + O(\|\alpha\|^2), \\ \lambda_2(\alpha) &= -1 + B_3\alpha_1 + B_4\alpha_2 + O(\|\alpha\|^2). \end{aligned} \quad (23) $$

The functions $\mu_1$ and $\mu_2$ can be taken as the unfolding parameters if the determinant of the Jacobi matrix

$$ \Delta = \det \left( \frac{\partial \mu}{\partial(\alpha_1, \alpha_2)} \right) \Big|_{\alpha=0} \neq 0. $$

From (22) we have

$$ \Delta = \frac{1}{2h_{11}} \left[ \begin{aligned} & (A_1 B_2 - A_2 B_1)(g_{11} h_{11} - g_{20} h_{02}) \\ & + 2(A_1 A_4 - A_2 A_3)h_{11} \\ & + 2(A_2 B_3 - A_1 B_4)g_{20} \end{aligned} \right]_{\alpha=0}. $$

On the other hand, taking into account (23), we obtain by direct computations:

$$ \det \left( \frac{\partial T}{\partial(\eta, \alpha)} \right) \Big|_{\eta=\alpha=0} = 4h_{11}\Delta, $$

where $T$ is the map (17) written in the $(\eta, \alpha)$-coordinates, i.e. for the map (20). Thus, if $h_{11}(0) \neq 0$, the regularity of (17) at the origin is equivalent to $\Delta \neq 0$.

The scaling

$$ x_1 \mapsto \frac{x_1}{E(\mu)}, \quad \mu_1 \mapsto \frac{\mu_1}{E(\mu)}, $$

where $E(\mu) = h_{11}(0) + O(\|\alpha\|)$, gives finally (18). Obviously, the critical coefficients are the same as in Proposition 2.1.1. ■

## 2.2. Center manifold reduction

We now apply a special version of the center manifold reduction combined with the normalization to our model (see [Kuznetsov, 1999; Beyn et al., 2002] and references therein for the case of ordinary differential equations). Consider a two-parameter family of maps $F: \mathbb{R}^n \times \mathbb{R}^2 \to \mathbb{R}^n$,

$$ x \mapsto F(x, \alpha), \quad (24) $$

having at $\alpha = 0$ a fixed point $x = 0$ with two simple $\lambda_{1,2} = \pm 1$ and no other critical eigenvalues. We can write:

$$ F(x, \alpha) = Ax + A_1\alpha + \frac{1}{2}B(x, x) + B_1(x, \alpha) \\ + \frac{1}{6}C(x, x, x) + \frac{1}{2}C_1(x, x, \alpha) + \dots, $$
where only relevant homogeneous linear, quadratic and cubic terms are displayed. The multilin-ear forms B and C are defined in Sec. 1, while $B_1(q, \beta) = F_{xx\alpha}[q, \beta]$ and $C_1(q, p, \beta) = F_{xx\alpha}[q, p, \beta]$ are evaluated at $(x, \alpha) = (0, 0)$. The matrix A has the eigenvalues $\lambda_{1,2} = \pm 1$. Introduce the associated eigenvectors $q, q^*, p, p^*$, such that

$$Aq = q, \quad A^T q^* = q^*, \quad \langle q^*, q \rangle = 1,$$

$$Ap = -p, \quad A^T p^* = -p^*, \quad \langle p^*, p \rangle = 1.$$

We wish to restrict (24) to its two-dimensional center manifold

$$x = H(w, \beta), \quad H: \mathbb{R}^2 \times \mathbb{R}^2 \to \mathbb{R}^n, \quad (25)$$

depending on two parameters $(\beta_1, \beta_2)$. Expand H as

$$H(w, \beta) = \sum_{j_1+j_2+k_1+k_2 \ge 1} \frac{1}{j_1!j_2!k_1!k_2!} h_{j_1j_2k_1k_2} w_1^{j_1} w_2^{j_2} \beta_1^{k_1} \beta_2^{k_2}.$$

Using the freedom of choosing the w-coordinates on the center manifold, we can assume that the restric-tion of (24) to (25),

$$w \mapsto G(w, \beta), \quad G: \mathbb{R}^2 \times \mathbb{R}^2 \to \mathbb{R}^2, \quad (26)$$

has been put into the form

$$G(w, \beta) = \left( \begin{array}{l} \beta_1 + (1 + \beta_2)w_1 + \frac{1}{2}a_1(\beta)w_1^2 + \frac{1}{2}b_1(\beta)w_2^2 + \frac{1}{6}c_1(\beta)w_1^3 + \frac{1}{2}c_2(\beta)w_1w_2^2 \\ -w_2 + e_1(\beta)w_1w_2 + \frac{1}{2}c_3(\beta)w_1^2w_2 + \frac{1}{6}c_4(\beta)w_2^3 \end{array} \right) + O(\|w\|^4)$$

Furthermore, we assume a dependence

$$\alpha = K(\beta) = K_1\beta + O(\|\beta\|^2),$$

where $K: \mathbb{R}^2 \to \mathbb{R}^2$ and $K_1$ is a 2 × 2-matrix.

Remark 2.2.1. Note that we have included all terms which become resonant at $\alpha = 0$ in the second com-ponent of the normal form (26). In the proof of Proposition 2.1.1 we used primary and secondary normalization. By the primary normalization we ap-plied a quadratic change of coordinates to remove nonresonant quadratic terms, and a cubic change for nonresonant cubic terms. The secondary nor-malization used quadratic terms to remove some resonant cubic terms, if an additional nondegener-acy condition is satisfied. Sometimes this is called a *hypernormalization*. The coefficient $d(0)$ is changed by this secondary normalization. A *priori* we do not know if the nondegeneracy condition $h_{11} \neq 0$ is satisfied. Therefore we compute all cubic coef-ficients in the second component of G, and then use formulas (10) and (11) to obtain the critical coefficients.

The condition of invariance of the center man-ifold is the homological equation

$$F(H(w, \beta), K(\beta)) = H(G(w, \beta), \beta). \quad (27)$$

Collecting quadratic terms in Eq. (27) gives for

$$\beta = 0$$

$$\begin{align} (A - I_n)h_{2000} &= a_1(0)q - B(q, q), \tag{28} \\ (A + I_n)h_{1100} &= e_1(0)p - B(p, q), \nonumber \\ (A - I_n)h_{0200} &= b_1(0)q - B(p, p), \nonumber \end{align}$$

while the cubic coefficients are obtained from

$$\begin{align*} (A - I_n)h_{3000} &= c_1(0)q + 3a_1(0)h_{2000} - 3B(q, h_{2000}) \\ &\quad - C(q, q), \\ (A - I_n)h_{1200} &= c_2(0)q + b_1(0)h_{2000} - 2e_1(0)h_{0200} \\ &\quad - B(q, h_{0200}) - 2B(p, h_{1100}) \\ &\quad - C(p, q), \\ (A + I_n)h_{2100} &= c_3(0)p + (2e(0) - a(0))h_{1100} \\ &\quad - 2B(q, h_{1100}) - B(p, h_{2000}) \\ &\quad - C(q, q, p), \\ (A + I_n)h_{0300} &= c_4(0)p - 3b_1(0)h_{1100} \\ &\quad - 3B(p, h_{0200}) - C(p, p, p). \end{align*}$$

Since 1 and −1 are simple eigenvalues of A, the Fredholm solvability condition yields the critical quadratic coefficients

$$\begin{align*} a_1(0) &= \langle q^*, B(q, q) \rangle, \quad e_1(0) = \langle p^*, B(p, q) \rangle, \\ b_1(0) &= \langle q^*, B(p, p) \rangle, \end{align*}$$
and using the bordering technique described in [Kuznetsov, 1999] we find

$$h_{2000} = (A - I)^{\text{INV}} (\langle q^*, B(q, q) \rangle q - B(q, q)),$$

$$h_{1100} = (A + I)^{\text{INV}} (\langle p^*, B(q, p) \rangle p - B(q, p)),$$

$$h_{0200} = (A - I)^{\text{INV}} (\langle q^*, B(p, p) \rangle q - B(p, p)).$$

Here $h = L^{\text{INV}}R$ denotes the unique solution of the nonsingular bordered linear system:

$$\begin{pmatrix} L & v \\ v^* & 0 \end{pmatrix} \begin{pmatrix} h \\ s \end{pmatrix} = \begin{pmatrix} R \\ 0 \end{pmatrix}.$$

It is assumed that $L$ has a one-dimensional null-space, so that $Lv = L^T v^* = 0$, $\langle v^*, v \rangle = 1$.

It follows from (28) that $\langle q^*, h_{2000} \rangle = \langle q^*, h_{0200} \rangle = \langle p^*, h_{1100} \rangle = 0$. Therefore we can simply write the expressions for the critical cubic coefficients

$$c_1(0) = \langle q^*, C(q, q, q) + 3B(q, h_{2000}) \rangle,$$

$$c_2(0) = \langle q^*, C(q, p, p) + B(q, h_{0200}) + 2B(p, h_{1100}) \rangle,$$

$$c_3(0) = \langle p^*, C(q, q, p) + B(p, h_{2000}) + 2B(q, h_{1100}) \rangle,$$

$$c_4(0) = \langle p^*, C(p, p, p) + 3B(p, h_{0200}) \rangle.$$

**Remark 2.2.2.** Note that we do calculate $e_1(0)$ in the center manifold reduction. If $e_1(0) \neq 0$, we can scale and using the previous remark obtain

$$a(0) = \frac{a_1(0)}{e_1(0)}, \quad b(0) = b_1(0)e_1(0), \quad c(0) = \frac{c_1(0)}{e_1^2(0)},$$

$$d(0) = c_2(0) + \frac{1}{e_1(0)} \left( b_1(0)c_3(0) - \frac{1}{3}(2e_1(0) + a_1(0))c_4(0) \right).$$

Then $a(0)$, $b(0)$, $c(0)$ and $d(0)$ are the coefficients for (9).

For $\beta \neq 0$ we obtain from Eq. (27) the coefficients $h_{j_1j_2k_1k_2}$ and $(K_1)_{ij}$. This shows us how the center manifold and $\alpha$ depend on $\beta$ in linear approximation. The w-independent terms give

$$(A - I_n)[h_{0010} \ h_{0001}]\beta = ([q \ 0] - A_1K_1)\beta.$$

Here $[h_{0010} \ h_{0001}]$ should be interpreted as a $2 \times n$ matrix and so on. We take the matrix product with $q^*$ and find $q^*A_1K_1 = (1, 0)$. We have some freedom and can scale $K_1$ and take as a solution

$$K_1 = \frac{1}{\gamma_1^2 + \gamma_2^2} \begin{pmatrix} \gamma_1 & -\gamma_2 \\ \gamma_2 & \gamma_1 \end{pmatrix}, \\ \text{with } (\gamma_1, \gamma_2) = q^*A_1. \qquad (29)$$

Then one can also solve for $h_{0010}$, $h_{0001}$. The equations for other coefficients are

$$(A - I_n)[h_{1010} \ h_{1001}] = [2h_{2000} \ q] - 2[B(q, h_{0010}) \ B(q, h_{0001})] - B_1(q, K_1),$$

$$(A + I_n)[h_{0110} \ h_{0101}] = [-h_{1100} \ 0] - 2[B(p, h_{0010}) \ B(p, h_{0001})] - B_1(p, K_1),$$

$$(A - I_n)[h_{2010} \ h_{2001}] = a(0)[h_{1010} \ h_{1001}] + 2[a_1(0) \ h_{2000} \ h_{2000}] - 2[B(h_{0010}, h_{2000})B(h_{0001}, h_{2000})] - B_1(h_{2000}, K_1) + C_1(q, q, K_1),$$

$$(A + I_n)[h_{1110} \ h_{1101}] = e_1(0)[h_{0110} \ h_{0101}] + [e_1(0)h_{1100} - h_{1100}] - 2[B(h_{0010}, h_{1100})B(h_{0001}, h_{1100})] - B_1(h_{1100}, K_1) - 2C_1(q, p, K_1),$$

$$(A - I_n)[h_{0210} \ h_{0201}] = b_1(0)[h_{1010} \ h_{1001}] + 2[b_1(0)h_{2000} \ 0] - 2[B(h_{0010}, h_{0200})B(h_{0001}, h_{0200})] - B_1(h_{0200}, K_1) - C_1(p, p, K_1).$$

### 3. Analysis of the Normal Form

Introduce the truncated normal form

$$\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \mapsto N(x, \mu) = \begin{pmatrix} \mu_1 + (1 + \mu_2)x_1 + \frac{1}{2}a(\mu)x_1^2 + \frac{1}{2}b(\mu)x_2^2 + \frac{1}{6}c(\mu)x_1^3 + \frac{1}{2}d(\mu)x_1x_2^2 \\ -x_2 + x_1x_2 \end{pmatrix} \quad (30)$$
**Remark 3.0.3.** Note that (30) is invariant under the reflection in the $x_1$-axis:

$$x \mapsto Rx, \quad R = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \qquad (31)$$

for which $R^2 = I_2$. The phase portraits below will reflect this $\mathbb{Z}_2$-symmetry.

Denote the critical values of the normal form coefficients by

$$a_0 = a(0), \quad b_0 = b(0), \quad c_0 = c(0), \quad d_0 = d(0).$$

In this section, we study local and global bifurcations of the truncated normal form (30), present its bifurcation diagrams, and then briefly discuss relationships between (30) and (9).

### 3.1. Local codim 1 bifurcations

**Proposition 3.1.1.** The family of maps (30) has the following local codim 1 bifurcations in a sufficiently small neighborhood of $(x, \mu) = (0, 0)$.

$$A(x, \mu) = N_x(x, \mu) = \begin{pmatrix} 1 + \mu_2 + a(\mu)x_1 + \frac{1}{2}c(\mu)x_1^2 + \frac{1}{2}d(\mu)x_2^2 & b(\mu)x_2 + d(\mu)x_1x_2 \\ x_2 & -1 + x_1 \end{pmatrix}.$$

**Fold bifurcation.** The map (30) has a fixed point $x$ with multiplier 1 if

$$\begin{cases} N(x, \mu) = x, \\ \det(A(x, \mu) - I_2) = 0. \end{cases}$$

Using the Implicit Function Theorem, we see that this algebraic system has a unique solution curve $t_{\text{fold}}$ near the origin that is given in statement 1. The critical (adjoint) eigenvectors are $q = q^* = (1, 0)^T$, while

$$B(p, q) = \begin{pmatrix} a_0 p_1 q_1 + b_0 p_2 q_2 + O(\mu_2) \\ p_1 q_2 + p_2 q_1 \end{pmatrix}.$$

With (3) from Sec. 1, we obtain

$$a_{\text{fold}} = a_0 + O(\mu_2),$$

and if $a_0 \neq 0$ we have a nondegenerate (quadratic) fold when $\mu_2 \to 0$.

**Flip bifurcation.** Another look at the Jacobi matrix above shows that along the curve $t_{\text{flip}}$ defined

1. There is a curve

$$t_{\text{fold}} : (x_1, x_2, \mu_1) = \left(-\frac{\mu_2}{a_0} + O(\mu_2^2), 0, \frac{\mu_2^2}{2a_0} + O(\mu_2^3)\right),$$

on which a nondegenerate fold bifurcation occurs if $a_0 \neq 0$.

2. There is a curve $t_{\text{flip}} : (x_1, x_2, \mu_1) = (0, 0, 0)$ on which a nondegenerate flip bifurcation occurs if $b_0 \neq 0$.

3. If $b_0 > 0$ and $\mu_1 < 0$, there is a curve

$$t_{\text{NS}} : (x_1, x_2, \mu_2) = \left(0, \sqrt{-\frac{2\mu_1}{b_0}} + O(\mu_1^{3/2}), \frac{(d_0 + 2b_0)\mu_1}{b_0} + O(\mu_1^2)\right),$$

on which a nondegenerate Neimark-Sacker bifurcation of the second iterate of (30) occurs, provided

$$a_0^2 b_0 + 3a_0 b_0 + a_0 d_0 - b_0 c_0 \neq 0. \qquad (32)$$

*Proof.* The Jacobi matrix of (30) is

in the statement 2, the truncated normal form (30) has a fixed point with multiplier $-1$, i.e. this curve satisfies the algebraic system

$$\begin{cases} N(x, \mu) = x, \\ \det(A(x, \mu) + I_2) = 0. \end{cases}$$

Now we have $p = p^* = (0, 1)^T$ as (adjoint) eigenvector. Clearly $C(p, p, p) = 0$ and we compute the flip coefficient (5) as

$$b_{\text{flip}} = \langle p^*, 3B(p, (I_2 - A)^{-1}B(p, p)) \rangle = -\frac{3b_0}{\mu_2} + o(1),$$

where $o(1)$ is bounded as $\mu_2 \to 0$. Then, the flip bifurcation is nondegenerate if $b_0 \neq 0$.

**Neimark-Sacker bifurcation.** Considering the second iterate of (30) we solve for its fixed point with the determinant of the Jacobi matrix equal to one, i.e.

$$\begin{cases} N(N(x, \mu), \mu) = x, \\ \det A(N(x, \mu), \mu) \det A(x, \mu) - 1 = 0. \end{cases}$$
We find the following exact solution to this system:

$$x_1 = 0, \quad b(\mu)x_2^2 + 2\mu_1 = 0, \quad \mu_2 = \frac{d(\mu) + 2b(\mu)}{b(\mu)}\mu_1, \tag{33}$$

which implies the expansion for $t_{NS}$ in the statement 3. Evaluating the Jacobian matrix of the second iterate of (30) on (33), we find

$$A = (N(N(x, \mu), \mu))_x \\ = \begin{pmatrix} 1 + 6\mu_1 + 4\mu_1^2 & 2b(\mu)\sqrt{-\frac{2\mu_1}{b(\mu)}}(1 + \mu_1) \\ -2\sqrt{-\frac{2\mu_1}{b(\mu)}}(1 + \mu_1) & 1 + 2\mu_1 \end{pmatrix}.$$

For small $\mu_1 < 0$, it has complex eigenvalues $e^{\pm i\theta_0} = 1 + 4\mu_1 + 2\mu_1^2 \pm \sqrt{\mu_1(2 + \mu_1)(1 + \mu_1)^2}$. Therefore, we find that in the case of $b_0 > 0$, $\mu_1 < 0$ there is a Neimark-Sacker bifurcation.

We wish to know the sign of the first Lyapunov coefficient $c_{NS}$ along (33). We take

$$q = \left( \frac{b(\mu)}{2} \sqrt{-\frac{2\mu_1}{b(\mu)}} \left( 1 + \sqrt{\frac{2+\mu_1}{\mu_1}} \right), 1 \right)^T \text{ and}$$

$$q^* = \left( -\frac{1}{2} \sqrt{-\frac{2\mu_1}{b(\mu)}} \left( 1 - \sqrt{\frac{2+\mu_1}{\mu_1}} \right), 1 \right)^T,$$

such that $Aq = e^{i\theta_0}q$ and $A^Tq^* = e^{-i\theta_0}q^*$, with $e^{i\theta_0} = 1 + 4\mu_1 + 2\mu_1^2 - \sqrt{\mu_1(2 + \mu_1)(1 + \mu_1)^2}$. We should still scale $q^*$, since $\langle q^*, q \rangle \neq 1$. Next we compute the first Lyapunov coefficient $c_{NS}$ on $t_{NS}$, using (7) from Sec. 1, where the multilinear forms $B$ and $C$ correspond to the second iterate of (30). This gives:

$$c_{NS} = \frac{1}{2}(-a_0^2 b_0 - 3a_0 b_0 - a_0 d_0 + b_0 c_0) + o(\mu_1)$$

with

$$d_1 = \frac{1}{6} \left( c_0 - \frac{3}{2} a_0^2 \right), \quad d_2 = \frac{1}{2} \left( d_0 + \frac{1}{2} b_0 (2 - a_0) \right),$$

$$d_3 = \frac{1}{4}(a_0 - 2), \quad d_4 = \frac{1}{4}b_0.$$

as $\mu_1 \to 0$. Therefore, the Neimark-Sacker bifurcation of the period-2 cycle of (30) is nondegenerate near the origin, if (32) holds. ■

**Remark 3.1.2.** [Gheiner, 1994] obtained similar results of the fold and flip bifurcations, but some calculations are puzzling. The Neimark-Sacker bifurcation was proven to exist, but no attempt had been made to analyze it. As we have seen, this bifurcation is nondegenerate under the condition (32) on the critical normal form coefficients. Choosing, for numerical example, $b_0 = c_0 = d_0 = 1$ we find $a_0 \neq -2 \pm \sqrt{5}$ as a condition for $c_{NS}$ to be nonzero.

## 3.2. Global bifurcations

As we shall see, the map (30) has two saddle fixed points, which can possess a heteroclinic structure. To study this global bifurcation phenomenon, we derive a vector field, such that the unit shift along its orbits approximates (30). Bifurcations of this vector field are easy to analyze, since it is similar to an amplitude system for the fold-Hopf bifurcation (see [Chow et al., 1994; Kuznetsov, 1998]).

**Proposition 3.2.1** (Approximating vector field). In a small neighborhood of $(x, \mu) = (0, 0)$, the truncated normal form (30) satisfies

$$RN(x, \mu) = \varphi^1(x, \mu) + O(\|\mu\|^2) + O(\|x\|^2\|\mu\|) \\ + O(\|x\|^4). \tag{34}$$

Here $R$ is the matrix defined by (31), $\varphi^t$ is the flow generated by the system

$$\dot{x} = X(x, \mu), \quad x \in \mathbb{R}^2, \quad \mu \in \mathbb{R}^2,$$

where the vector field $X$ is given by

$$X(x, \mu) = \begin{pmatrix} \mu_1 + \left(-\frac{1}{2}a_0\mu_1 + \mu_2\right)x_1 + \frac{1}{2}a_0x_1^2 + \frac{1}{2}b_0x_2^2 + d_1x_1^3 + d_2x_1x_2^2 \\ \frac{1}{2}\mu_1x_2 - x_1x_2 + d_3x_1x_2^2 + d_4x_2^3 \end{pmatrix}, \tag{35}$$

*Proof.* We construct $\varphi^t$ as the first two-dimensional component of the flow

$$\xi \mapsto \phi^t(\xi) = \begin{pmatrix} \varphi^t(x, \mu) \\ \mu \end{pmatrix}, \quad \xi = \begin{pmatrix} x \\ \mu \end{pmatrix} \in \mathbb{R}^4,$$
generated by a four-dimensional system with the
parameters considered as constant variables:

$$
\dot{\xi} = Y(\xi). \tag{36}
$$

Here

$$
Y(\xi) = J\xi + Y_2(\xi) + Y_3(\xi) + \dots,
$$

$$
J = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad Y_k(\xi) = \begin{pmatrix} X_k(\xi) \\ 0 \end{pmatrix},
$$

where each $X_k$ is an order-$k$ homogeneous function
from $\mathbb{R}^4$ to $\mathbb{R}^2$ with unknown coefficients. Define

$$
M(\xi) = \begin{pmatrix} N(x, \mu) \\ \mu \end{pmatrix}
$$

and introduce the 4 × 4 block-diagonal matrix

$$
S = \begin{pmatrix} R & 0 \\ 0 & I_2 \end{pmatrix},
$$

where *R* is given in (31). We look for a vector field
*Y* such that *S* *M*(*ξ*) = *φ*¹(*ξ*) + *O*(||*ξ*||⁴) (cf. [Takens,
1974]).

To find the vector field $Y$ explicitly, perform
three Picard iterations for (36) as described in
[Kuznetsov, 1998, Chap. 9]. We start with setting
$\phi_1^t(\xi) = e^{Jt}\xi$. Then, clearly, the linear part of
$SM(\xi)$ coincides with $\phi_1^t(\xi)$.

Since we know how the result $\phi_2^t$ of the second
Picard iteration should look like, we set some coef-
ficients of $Y_2$ equal to zero immediately:

$$
Y_2 = \begin{pmatrix}
A_{10}\mu_1 x_1 + A_{01}\mu_2 x_1 + \frac{1}{2}A_{20}x_1^2 + \frac{1}{2}A_{02}x_2^2 \\
B_{11}x_1 x_2 + B_{10}\mu_1 x_2 \\
0 \\
0
\end{pmatrix}.
$$

Then

$$
\phi_2^t(\xi) = e^{Jt}\xi + \int_0^t e^{J(t-\tau)} Y_2(\phi_1^\tau(\xi)) d\tau \\
= \begin{pmatrix} x_1 + t\mu_1 \\ x_2 \\ \mu_1 \\ \mu_2 \end{pmatrix} + \begin{pmatrix} \left(A_{10}t + \frac{1}{2}A_{20}t^2\right)\mu_1 x_1 + A_{01}\mu_2 x_1 t + \frac{1}{2}A_{20}x_1^2 t + \frac{1}{2}A_{02}x_2^2 t \\ B_{11}x_1 x_2 t + \left(\frac{1}{2}B_{11}t^2 + B_{10}t\right)\mu_1 x_2 \\ 0 \\ 0 \end{pmatrix} + O(\|\mu\|^2).
$$

Comparing quadratic terms in $SM(\xi)$ and $\phi_2^1(\xi)$,
we find the coefficients of $Y_2$

$$
A_{10} = -\frac{1}{2}a_0, A_{20} = a_0, A_{01} = 1,
$$

$$
A_{02} = b_0, B_{10} = \frac{1}{2}, B_{11} = -1.
$$

Passing on to the cubic part we remark that we
are only interested in cubic terms in $x$. Therefore,

we put

$$
Y_3 = \begin{pmatrix}
\sum_{i+j=3} \frac{1}{i!j!} A_{ij} x_1^i x_2^j \\
\sum_{i+j=3} \frac{1}{i!j!} B_{ij} x_1^i x_2^j \\
0 \\
0
\end{pmatrix}
$$

and get

$$
\phi_3^t(\xi) = e^{Jt}\xi + \int_0^t e^{J(t-\tau)} [Y_2(\phi_2^\tau(\xi)) + Y_3(\phi_2^\tau(\xi))] d\tau \\
= \begin{pmatrix} x_1 + t\mu_1 + \frac{1}{2}a_0(t^2-t)x_1\mu_1 + x_1\mu_2t + \frac{1}{2}a_0x_1^2t + \frac{1}{2}b_0x_2^2t \\ x_2 - x_1x_2t + \frac{1}{2}(t-t^2)\mu_1x_2 \\ \mu_1 \\ \mu_2 \end{pmatrix}
$$
$$ 
\begin{aligned}
& + \begin{pmatrix}
\left(\frac{1}{6}A_{30}t + \frac{1}{4}A_{20}^2t^2\right)x_1^3 + \frac{1}{2}tA_{21}x_1^2x_2 + \left(\frac{1}{2}A_{12}t + \frac{1}{4}t^2A_{02}(A_{20} + 2B_{11})\right)x_1x_2^2 + \frac{1}{6}tA_{03}x_2^3 \\
& + \frac{1}{6}tB_{30}x_1^3 + \left(\frac{1}{2}tB_{21} + \frac{1}{4}t^2B_{11}(2B_{11} + A_{20})\right)x_1^2x_2 + \frac{1}{2}tB_{12}x_1x_2^2 + \left(\frac{1}{6}tB_{03} + \frac{1}{4}t^2A_{02}B_{11}\right)x_2^3 \\
& 0 \\
& 0
\end{pmatrix} \\
& + O(\|\mu\|^2) + O(\|x\|^2\|\mu\|).
\end{aligned}
$$

Comparing cubic terms in $SM(\xi)$ and $\phi_3^1(\xi)$, we find the coefficients of $Y_3$:

$$ 
\begin{aligned}
A_{30} &= c_0 - \frac{3}{2}a_0^2, \quad A_{21} = 0, \\
A_{12} &= d_0 + \frac{1}{2}b_0(2 - a_0), \quad A_{03} = 0, \\
B_{30} &= 0, \quad B_{21} = \frac{1}{2}(a_0 - 2), \quad B_{12} = 0, \quad B_{03} = \frac{3}{2}b_0.
\end{aligned}
$$

This gives the vector field (35) from the Proposition. ■

To explore relationships between the map (30) and the vector field (35), consider first local bifurcations of the vector field $X$. One can check that there are two curves, $t_{\text{fold}} : (x_1, x_2, \mu_1) = (-\mu_2/a_0 + O(\mu_2^2), 0, \mu_2^2/(2a_0) + O(\mu_2^3))$ and $t_{\text{flip}} : (x_1, x_2, \mu_1) = (0, 0, 0)$, on which equilibria of (35) have a zero eigenvalue. These are the same expansions as we computed for the map (30). The center manifold reduction shows that a fold (saddle-node) bifurcation occurs on the first curve, while a pitchfork bifurcation occurs on the second. Next, we computed a Hopf bifurcation curve for (35). We get indeed the same expression

$$ t_{\text{NS}} : (x_1, x_2, \mu_2) = \left(0, \sqrt{-\frac{2\mu_1}{b_0}} + O(\mu_1^{3/2}), \frac{(2b_0 + d_0)\mu_1}{b_0} + O(\mu_1^2)\right) $$

as in Proposition 3.1.1.

Next, we can classify the critical phase portraits of the vector field (35). In the usual polar coordinates the vector field at $\mu = 0$ becomes

$$ \begin{pmatrix} \dot{r} \\ \dot{\theta} \end{pmatrix} = \begin{pmatrix} r^2 \left( \frac{1}{2} a_0 \cos^2 \theta + \left( \frac{1}{2} b_0 - 1 \right) \sin^2 \theta \right) + O(r^3) \\ -r \sin \theta \left( \left( 1 + \frac{1}{2} a_0 \right) \cos^2 \theta + \frac{1}{2} b_0 \sin^2 \theta \right) + O(r^2) \end{pmatrix}. $$

We see that there are invariant lines in the critical normal form if $\dot{\theta} = 0$. This equation is satisfied if $\theta = 0, \pi$, which is expected due to the invariance of

the vector field under the map (31). Another possibility is that

$$ \tan^2 \theta = -\frac{2 + a_0}{b_0}. $$

Therefore, we find six different critical portraits, see Figs. 1 and 2.

The invariant curves born from the Neimark-Sacker bifurcation cannot exist everywhere. They should disappear through some global bifurcations. To study this, we investigate what happens to the cycles of the approximating vector field $X$ born at the Hopf bifurcation.

**Proposition 3.2.2.** If $a_0, b_0 > 0$ and $\mu_1 < 0$ then the vector field (35) has two saddles, which are always connected by a heteroclinic orbit along the $x_1$-axis. There exists another heteroclinic orbit for

$$ t_J : \mu_2 = \frac{\mu_1}{3+a_0} \left( (a_0+2) \frac{d_0+2b_0}{b} + \frac{c_0-a_0-a_0^2}{a_0} \right) + o(\mu_1). \quad (37) $$

*Proof.* We first shift the $x_1$-coordinate in (35) with

$$ x_1 \to x_1 - \left( \frac{\mu_1}{2} + \frac{\mu_2}{a_0} \right). $$

Then we apply a singular rescaling

$$ x_1 \to \delta x_1, \quad x_2 \to \delta x_2, \quad dt \to \frac{x_2^q}{\delta} dt $$

to obtain

$$ \left\{ \begin{aligned} \dot{x}_1 &= x_2^q \left( \beta_1 + \frac{1}{2} a_0 x_1^2 + \frac{1}{2} b_0 x_2^2 + \delta[d_1 x_1^3 \right. \\ & \qquad \left. + d_2 x_1 x_2^2] \right) \\ \dot{x}_2 &= x_2^q \left( -x_1 x_2 + \delta[\beta_2 x_2 + d_3 x_1^2 x_2 + d_4 x_2^3] \right), \end{aligned} \right. \quad (38) $$

where

$$ \beta_1 = \mu_1 + O(\|\mu\|^2), \quad \beta_2 = \frac{\mu_2}{a_0} + O(\|\mu\|^2). $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Phase portraits of the approximating vector field X at $\mu = 0$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Phase portraits of the normal form N at $\mu = 0$. Compare with the orbits of the vector field to see how the orbits of the map advance.</figcaption></figure>

The system (38) can be rewritten as

$$ \dot{x} = f(x, \beta) + \delta g(x, \beta) \quad (39) $$

with

$$ f(x, \beta) = x_2^q \left( \beta_1 + \frac{1}{2} a_0 x_1^2 + \frac{1}{2} b_0 x_2^2 - x_1 x_2 \right), $$

$$ g(x, \beta) = x_2^q \begin{pmatrix} d_1 x_1^3 + d_2 x_1 x_2^2 \\ \beta_2 x_2 + d_3 x_1^2 x_2 + d_4 x_2^3 \end{pmatrix}. $$

For $\delta = 0$ and $q+1 = a_0$, system (39) is Hamiltonian

with

$$ H(x) = x_2^{q+1} \left( \frac{\beta_1 + \frac{1}{2} a_0 x_1^2}{q+1} + \frac{b_0 x_2^2}{2(q+3)} \right). \quad (40) $$

We have $a_0, b_0 \neq 0$ as nondegeneracy conditions, therefore $q \neq -1$. Level curves of $H$ for several values of $a_0$ are shown in Fig. 3. Now we treat the term $\delta g$ in (39) as a small perturbation of the Hamilton system. We should therefore evaluate the Pontryagin–Melnikov function (see [Guckenheimer & Holmes, 1983, 2002])

$$ \Delta(h, \beta) = \oint_{\Gamma_h} f(x_h(\tau), \beta) \wedge g(x_h(\tau), \beta) d\tau, $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Level curves of H for several $a_0$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Ratio $Q(h) = I_{1,h}/I_{3,h}$.</figcaption></figure>

where $x_h(\tau)$ is a periodic solution of the Hamiltonian system corresponding to a closed regular level set $\Gamma_h = \{x : H(x) = h\}$, while for $a_0 > 0$

$$ \Delta(0, \beta) = \int_{-\infty}^{+\infty} f(x_0(\tau), \beta) \wedge g(x_0(\tau), \beta) d\tau, $$

where $x_0(\tau)$ is the nontrivial heteroclinic solution in the critical level set $H = 0$. Notice that $\lim_{h \to 0^-} \Delta(h, \beta) = \Delta(0, \beta)$, since the integral over the trivial heteroclinic connection equals zero. Then $\Delta(0, \beta) = 0$ defines a linear approximation to a curve on which the heteroclinic connection “survives” in (39) for small $\delta \neq 0$. Our computation is analogous to the one in [Chow et al., 1994].

Using Green's formula, we have

$$ \begin{aligned} \Delta &= \oint_{\Gamma_h} x_2^q (d_1 x_1^3 + d_2 x_1 x_2^2) dx_2 \\ &\quad - x_2^q (\beta_2 x_2 + d_3 x_1^2 x_2 + d_4 x_2^3) dx_1 \\ &= \oint_{\Gamma_h} x_2^q \left( d_1 x_1^3 + d_2 x_1 x_2^2 + (q+1) \left( \beta_2 x_1 + \frac{1}{3} d_3 x_1^3 \right) \right. \\ &\qquad \left. + (q+3) d_4 x_1 x_2^2 \right) dx_2. \end{aligned} $$

Now, using it along $\Gamma_h$ we have $dH = x_2^q(\beta_1 + (1/2)a_0x_1^2 + (1/2)b_0x_2^2)dx_2 + x_2^{q+1}x_1dx_1 = 0$ and continue

$$ \begin{aligned} \Delta &= \oint_{\Gamma_h} x_2^q \left( a_0\beta_2x_1 + \left(d_1 + \frac{1}{3}a_0d_3\right)x_1^3 + (d_2 + (a_0+2)d_4)x_1x_2^2 \right) dx_2 \\ &= \oint_{\Gamma_h} x_2^q \left( x_1 \left(a_0\beta_2 - \frac{2}{b_0}\beta_1(d_2 + (a_0+2)d_4)\right) + x_1^3 \left(d_1 + \frac{1}{3}a_0d_3 - \frac{a_0}{b_0}(d_2 + (a_0+2)d_4)\right) \right) dx_2 \\ &\quad - \oint_{\Gamma_h} \frac{2}{b_0}(d_2 + (a_0+2)d_4)x_1^2 x_2^{q+1} dx_1 \\ &= I_{1,h} \left(a_0\beta_2 - \frac{2}{b_0}\beta_1(d_2 + (a_0+2)d_4)\right) + I_{3,h} \left(d_1 + \frac{a_0}{3}d_3 - \frac{a_0}{3b_0}(d_2 + (a_0+2)d_4)\right). \end{aligned} $$
Here, we define $I_{i,h} = \int_{\Gamma_h} x_2^q x_1^i dx_2$ for $i = 1, 3$.

For $h = 0$ we can evaluate the Pontryagin-Melnikov function as follows. We rewrite our Hamiltonian system as

$$ \begin{cases} \dot{x}_1 = x_2^q \left( \beta_1 + \frac{1}{2} a_0 x_1^2 + \frac{1}{2} b_0 x_2^2 \right) = \frac{b_0}{a_0 + 2} x_2^{q+2} \\ \dot{x}_2 = x_2^q (-x_1 x_2). \end{cases} $$

Here, we used the zero-level set of the Hamiltonian. The initial conditions are

$$ x_1(0) = 0, \quad x_2(0) = \sqrt{\frac{-2(a_0 + 2)\beta_1}{a_0 b_0}}, $$

which are not changed if we reparametrize $d\tau = x_2^{-q} dt$. We get the new system

$$ \begin{cases} \dot{x}_1 = \frac{b_0}{a_0 + 2} x_2^2 \\ \dot{x}_2 = -x_1 x_2, \end{cases} $$

which can be explicitly solved with $x_1(\tau) = \gamma \tanh(\gamma\tau)$, $x_2(\tau) = x_2(0) \cosh^{-1}(\gamma\tau)$, where

$$ \gamma = \sqrt{\frac{-2\beta_1}{a_0}}. $$

With these solutions, we evaluate $I_{i,0}$ defined above:

$$ \begin{aligned} I_{i,0} &= \int_{\Gamma_0} x_2^q x_1^i dx_2 \\ &= - \int_{-\infty}^{\infty} x_2^{a_0} x_1^{i+1} d\tau \\ &= - \int_{-\infty}^{\infty} \frac{\sinh(\gamma\tau)^{i+1}}{\cosh(\gamma\tau)^{i+a_0+1}} d\tau \\ &= -\gamma^i x_2^q(0) (1 + (-1)^{i+1}) \frac{\Gamma\left(\frac{a_0}{2}\right) \Gamma\left(\frac{2+i}{2}\right)}{2\Gamma\left(\frac{1}{2}(2+a_0+i)\right)}. \end{aligned} $$

In this time-parametrization, we have used the explicit forms of $x_1(\tau)$ and $x_2(\tau)$. If we now compute the ratio $Q = I_{3,0}/I_{1,0}$, we find

$$ Q = \frac{-6\beta_1}{a_0(3+a_0)}. $$

Using the expressions for $\beta_k$ and $d_i$, we find that the Melnikov function $\Delta(0, \beta)$ has a zero if

$$ \mu_2 = \frac{\mu_1}{3 + a_0} \left( (a_0 + 2) \frac{d_0 + 2b_0}{b_0} + \frac{c_0 - a_0 - a_0^2}{a_0} \right) + o(\mu_1). \quad (41) $$

This value of $\mu_2$ asymptotically corresponds to the existence of a nontrivial heteroclinic orbit for the perturbed system (38). ■

**Remark 3.2.3.** Taking $a_0 = b_0 = c_0 = d_0 = 1$ we find $\mu_2 = 2\mu_1 + o(\mu_1)$.

**Remark 3.2.4.** We derived the linear approximations to both the Neimark-Sacker bifurcation curve, see Proposition 3.1.1,

$$ \mu_2 = \frac{d_0 + 2b_0}{b_0} \mu_1 \quad (42) $$

and the heteroclinic bifurcation curve, see (41). To analyze their relative position, we compute the difference between their slopes:

$$ \frac{1}{a_0 b_0 (a_0 + 3)} (-a_0 d_0 - 3 a_0 b_0 - b_0 a_0^2 + b_0 c_0). $$

This shows that the curves coincide in the linear approximation if and only if the Lyapunov coefficient $c_{NS}$ vanishes. Thus, changing the relative position of the two curves changes the stability of the closed invariant curve.

**Remark 3.2.5.** Moreover, for the vector field we have the uniqueness of the limit cycle. This can be verified as follows. We should evaluate the Pontryagin function on a level curve of the Hamiltonian with $h \neq 0$. Now $Q(h) = I_{1,h}/I_{3,h}$ cannot be evaluated explicitly, but one can prove the monotonicity of $Q$ as Chow et al. [1994]. This implies the uniqueness of the limit cycle. We include some pictures, we computed $Q(h)$ numerically, which illustrate the monotonicity (see Fig. 4).

### 3.3. Bifurcation diagrams

Although we had six cases for $F_0$, only four bifurcation diagrams will be reported in Figs. 7–14, because two others differ only at the critical parameter values. We start with the bifurcation diagrams of the vector field (35), which are similar to the bifurcation diagrams of the truncated amplitude system for the fold-Hopf bifurcation (see [Guckenheimer & Holmes, 1983, 2002; Kuznetsov, 1998]). In our study, however, we have to take into account that the Neimark-Sacker bifurcation can be either sub- or super-critical, depending on the sign of $c_{NS}$. We included these sketches to indicate the direction of the orbit in the phase portraits for the map. The orbits of the map continuously jump from the lower to the upper half plane and back. This is easily understood from (34), which
<figure><img src="image_6.png" /><figcaption>Fig. 5. Heteroclinic tangencies appearing in the lines $h_{1,2}$ together with a transversal heteroclinic structure between them.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Heteroclinic structures in (30) near the left fixed point for $a = b = d = 1$, $c = 0.5$, $\mu_1 = -0.2$ and (a) $\mu_2 = -0.35345880$; (b) $\mu_2 = -0.35347$; (c) $\mu_2 = -0.35349058$.</figcaption></figure>

implies that (9) can be approximated by the composition of the unit shift along the orbits of X with the reflection R. We go around the origin in the parameter plane and discuss the phase portraits of the truncated normal form (30). The fold and flip bifurcation curves are denoted by $F_{\pm}$ and $P_{\pm}$, respectively. The sub/super-critical Neimark-Sacker bifurcation curve is indicated by $NS_{\pm}$. For the vector field, J is the heteroclinic bifurcation curve.

**Remark 3.3.1.** We took $c = d = 1$ to generate the pictures related to the map, unless stated otherwise. We use a version of DSTOOL by Back *et al.* [1992], incorporating the algorithms from Krauskopf and Osinga [1998a], Krauskopf and Osinga [1998b]. The colors should be interpreted as follow: Orbits are green, while unstable manifolds of saddles are red and stable manifolds are blue. Note that, with the above choice of the coefficients, $c_{NS}$ is negative in case 1 and positive in case 2.

* **Case 1.** In region 1 orbits merely jump to the right. Crossing $F_+$ implies the appearance of two

fixed points on the horizontal axis. In **2** one of these fixed points is totally unstable, while the other is a saddle. While crossing curve $P_+$ from **2** to **3**, the unstable fixed point becomes a saddle and an unstable period-2 cycle appears. If $c_{NS} > 0$, an unstable invariant curve “around” the period-2 cycles appears via the Neimark-Sacker bifurcation on $NS_+$, when we go from **3** to **4+**. The invariant curve disappears through a series of bifurcations associated with the heteroclinic bifurcations near $J_+$, if we come to **5**. The presence of the *J*-curve for the vector field implies for the map the existence of two curves, along which heteroclinic tangencies occur (see sketches in Fig. 5). Between these two curves, a heteroclinic structure is present. Figure 6 shows invariant curve configurations computed with DSTOOL. If $c_{NS} < 0$, a stable closed invariant curve emerges in **4-** through a series of bifurcations associated with the heteroclinic structure. This stable invariant curve exists until we cross $NS_-$, where the stable period-2 cycle becomes attracting in **5**. Next we cross $P_-$ and the period-2 orbit disappears, leaving us with a stable fixed point and
<figure><img src="image_8.png" /><figcaption>Fig. 7. Vector field: Case 1. $a_0 > 0$, $b_0 > 0$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Map: Case 1. $a_0 > 0$, $b_0 > 0$.</figcaption></figure>

a saddle in **6**. These two collide if we return back to **1**.

• **Case 2.** Fix a phase domain near the origin. Now we start with the two fixed points, one stable and

one unstable on the axis in region **1**. Then, crossing the flip curve $P_+$ to **2**, one fixed point exhibits a period doubling and a period-2 cycle appears. The fixed points on the horizontal axis collide at
<figure><img src="image_10.png" /><figcaption>Fig. 9. Vector field: Case 2. $a_0 < 0$, $b_0 > 0$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Map: Case 2. $a_0 < 0$, $b_0 > 0$.</figcaption></figure>

$F_+$ separating region **2** from **3**, where a stable period-2 cycle exists. If $c_{NS} > 0$, then an unstable invariant curve appears when we cross the Neimark-Sacker bifurcation curve $NS_+$. This in-

variant curve grows, until it blows up and disappears from the fixed phase domain at some curve $B_+$. Actually, the invariant curve can lose its smoothness and disappear before touching the
<figure><img src="image_12.png" /><figcaption>Fig. 11. Vector field: Case 3. $a_0 > 0, b_0 < 0$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Map: Case 3. $a_0 > 0, b_0 < 0$.</figcaption></figure>

boundary of the domain. If $c_{NS} < 0$, then we first encounter the “boundary bifurcation” curve $B_-$, where a big stable invariant curve appears in our fixed phase domain. The transition from **4-** to **5** destroys the curve via the Neimark-Sacker bifurcation. Finally, crossing the fold curve $F_-$ produces two fixed points in **6** and through the flip bifurcation on $P_-$ the period-2 cycle disappears again as we are back in **1**.

• **Case 3.** We start with a period-2 saddle cycle in **1**. Entering **2** through the fold curve $F_+$ creates two fixed points on the horizontal axis, a saddle

and a repelling one. Then, while crossing the flip curve $P_+$ to **3**, the period-2 cycle is destroyed and we get two saddles on the $x_1$-axis. Passing $P_-$ one saddle becomes stable and a period-2 orbit in **4** is created. Finally, the fixed points on the horizontal axis collide on $F_-$ and we are in region **1** again.

• **Case 4.** Starting in region **1** we have as in Case 3 a period-two saddle cycle, but also a stable and an unstable fixed point on the $x_1$-axis. The unstable one becomes a saddle when we enter **2** through $P_+$ curve. Then nothing special appears
<figure><img src="image_14.png" /><figcaption>Fig. 13. Vector field: Case 4. $a_0 < 0$, $b_0 < 0$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Map: Case 4. $a_0 < 0$, $b_0 < 0$.</figcaption></figure>

but a “saddle-like flow” in **3** after the saddle and the stable point collided on $F_+$. Going from **3** to **4** we get a saddle and an unstable point through the fold bifurcation on the curve $F_-$. We are back in **1**, when the flip bifurcation creates the period-2 cycle on $P_-$.

The obtained diagrams give a rather detailed description of bifurcations of the truncated normal form (30). However, this description remains incomplete due to the presence of closed invariant curves and heteroclinic tangencies. Indeed, the rotation number on the closed invariant curve can change infinitely many times from rational to irrational and,

moreover, the invariant curve can lose smoothness and disappear. Near a heteroclinic tangency, infinite series of bifurcations happen, including cascades of flips and folds (see [Gavrilov & Shil'nikov, 1972, 1973; Gonchenko et al., 1996]).

<figure><img src="image_16.png" /><figcaption>Fig. 15. A transversal heteroclinic structure near the horizontal axis.</figcaption></figure>
### 3.4. Effect of higher order terms

Adding higher order terms to the truncated normal form (30), i.e. restoring (9), complicates the bifurcation picture further.

Using the Implicit Function Theorem, one can prove that for $||μ||$ sufficiently small, the map (9) has the same local bifurcations as (30). More precisely, Proposition 3.1.1 is valid also for the full normal form (9) with any higher order terms. Therefore, we know what to expect locally. In particular, in Cases 1 and 2 closed invariant curves appear. Moreover, the unit shift along the orbits of the vector field (35) composed of the reflection approximates (9) as good as (30). This implies that (9) also has two bifurcation curves along which heteroclinic tangencies occur. Between these curves, a heteroclinic structure is present. Higher order terms in (9) do affect these curves, but they both remain tangent to the curve (37) from Proposition 3.2.2.

As mentioned in [Gheiner, 1994], there are more differences between the phase portraits of (30) and a generic (9), which are related to other heteroclinic tangencies. Indeed, in the truncated normal form (9) the $x_1$-axis is always invariant. Therefore, in Cases 1 and 3 we have heteroclinic connections between the saddles located on the horizontal axis. However generically, the higher order terms in (9) break down the reflectional symmetry and the heteroclinic connection along the $x_1$-axis is lost. This allows for heteroclinic structures caused by intersections of the invariant manifolds of the saddles near the horizontal axis. These intersections can be either transversal (as in Fig. 15) or tangential. Therefore, in the first three cases, the bifurcation diagrams of (30) and a generic (9) are not locally topologically equivalent.

Gheiner [1994] showed that in Case 3, an additional heteroclinic structure may appear: The unstable manifold of the period-2 cycle could intersect tangentially the stable manifold of a saddle fixed point on the $x_1$-axis.

In Case 4, Gheiner [1994] gave a strong indication that (9) is locally topologically equivalent to (30), where the cubic terms can be omitted.

## 4. Examples

In this section we present two examples of the fold-flip bifurcation. We compute the normal form coefficients and show that the behavior predicted by the normal form is correct.

### 4.1. A note on the generalized Hénon map

Consider the generalized Hénon map

$$ \begin{pmatrix} x \\ y \end{pmatrix} \mapsto \begin{pmatrix} y \\ \alpha - \beta x - y^2 + Rxy + Sy^3 \end{pmatrix}. \quad (43) $$

Setting $R = S = 0$ one obtains the standard Hénon map. We will consider $α$ and $β$ as control parameters and $R$ and $S$ as constants. This map appears as a rescaled first return map in the study of at least two global bifurcations of maps related to codim 2 homoclinic tangencies, when:

(1) a diffeomorphism in $R^2$ has a neutral saddle with a quadratic homoclinic tangency (see [Gonchenko & Gonchenko, 2002; Gonchenko, 2002]);

(2) a diffeomorphism in $R^3$ has a saddle with a generalized homoclinic tangency (i.e. the unstable manifold of the saddle has a quadratic tangency to its stable manifold but is nontransversal to leaves of the strong stable foliation in the stable manifold at homoclinic points) (see [Gonchenko et al., 2001]).

In these studies, codim 2 local bifurcations of fixed points of (43) play an important role, since they allow to predict bifurcations of closed invariant curves. Bifurcations of (43) are well understood for $β > 0$, where the standard Hénon map preserves orientation. In this parameter region, strong resonance points have been found. Much less is known about bifurcations of (43) when $β < 0$. We show below that the map has a fold-flip point in this parameter region, and compute its normal form.

For $α = 0$ the map has the fixed point $(x, y) = (0, 0)$. The Jacobi matrix of (43) at this point is

$$ A = \begin{pmatrix} 0 & 1 \\ -\beta & 0 \end{pmatrix}. $$

Thus, if $β = -1$, we have two multipliers $+1$ and $-1$. So we have found a fold-flip codim 2 point for all values of $R$ and $S$.

First we analyze (43) without the cubic term, i.e. set $S = 0$. It is easy to verify that fixed points of this map can bifurcate at the following curves:

$$ t_{\text{fold}} = \left\{ (\alpha, \beta) : \alpha = \frac{(\beta + 1)^2}{4(R - 1)} \right\}, $$

$$ t_{\text{flip}} = \left\{ (\alpha, \beta) : \alpha = \frac{1}{4}(\beta + 1)^2(3 - R) \right\}. $$
<figure><img src="image_17.png" /><figcaption>Fig. 16. The $(\alpha, \beta)$-plane: (a) Case 1; (b) Case 2; (c) Case 3/4. The thick/dashed line is the fold/flip curve. The dotted one is the Neimark–Sacker curve.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Selected phase portraits of (43).</figcaption></figure>

Note that these two parabolas have a common point: The fold-flip point $(\alpha, \beta) = (0, -1)$. They coincide if $R = 2$, i.e. the flip bifurcation is degenerate for $R = 2$.

If we consider the second iterate of (43) with $S = 0$, a curve can be found, where a period-2 cycle has two multipliers with unit product:

$$t_{\text{NS}} = \left\{ (\alpha, \beta) : \alpha = \frac{(\beta+1)(\beta + R^2 - (1+R))}{R^2} \right\}.$$

To exclude neutral saddles, we should check that the multipliers are indeed complex for some values of $\alpha, \beta$ and $R$. Simple analysis shows that the Neimark–Sacker bifurcation occurs only if $R \neq 0$ and $R > -1$. This means that depending on whether $R < -1$ or $R > -1$ we have either the “difficult” or the “easy” case, respectively. For the difficult case we have a distinction between $R > 1$

and $R < 1$, since the fold and flip curve changes position in the parameter plane when crossing $R = 1$. We sketch in Fig. 16 the bifurcation curves in the $(\alpha, \beta)$-plane for three typical values of $R$. We recognize (a) and (b) as Cases 1 and 2 from Sec. 3.3, respectively, and (c) as Case 3 or 4.

We now compute the critical normal form coefficients using the formulas obtained in Sec. 2.2. The eigenvectors are

$$q = 2q^* = \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \quad p = 2p^* = \begin{pmatrix} 1 \\ -1 \end{pmatrix},$$

while the multilinear functions are given by

$$B(x, y) = \begin{pmatrix} 0 \\ R(x_1 y_2 + x_2 y_1) - 2x_2 y_2 \end{pmatrix},$$

$$C(x, y, z) = \begin{pmatrix} 0 \\ 6Sx_2y_2z_2 \end{pmatrix}.$$
For the quadratic coefficients we find

$$
\begin{align*}
a_1(0) &= -(1-R), & e_1(0) &= -1, \\
b_1(0) &= -(R+1), \\
h_{2000} &= \frac{1}{2}(1-R) \begin{pmatrix} 1 \\ -1 \end{pmatrix}, & h_{1100} &= -\frac{1}{2} \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \\
h_{0200} &= \frac{1}{2}(1+R) \begin{pmatrix} 1 \\ -1 \end{pmatrix},
\end{align*}
$$

then we get

$$
\begin{align*}
c_1(0) &= \frac{3}{2}(1-R), & c_2(0) &= -\frac{1}{2}(1+R), \\
c_3(0) &= -\frac{1}{2}(1-R)^2, & c_4(0) &= \frac{3}{2}(1+R)^2.
\end{align*}
$$

The normalized coefficients of (9) are therefore

$$
\begin{align*}
a(0) &= (1-R), & b(0) &= (1+R), \\
c(0) &= \frac{3}{2}(1-R), & d(0) &= -\frac{1}{2}(5+3R)
\end{align*}
$$

(see Remark 2.2.2). We see that, indeed, depending
on R the following cases occur:

<table>
  <thead>
    <tr>
      <th></th>
      <th>a</th>
      <th>b</th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>R &lt; -1</td>
      <td>+</td>
      <td>-</td>
      <td>case 3</td>
    </tr>
    <tr>
      <td>-1 &lt; R &lt; 1</td>
      <td>+</td>
      <td>+</td>
      <td>case 1</td>
    </tr>
    <tr>
      <td>R &gt; 1</td>
      <td>-</td>
      <td>+</td>
      <td>case 2</td>
    </tr>
  </tbody>
</table>

For Cases 1 and 2 we calculate the leading term of
the Lyapunov coefficient

$c_{\text{NS}} = R^2(1-R),$

so that the closed invariant curves appearing via
the Neimark–Sacker bifurcation will be unstable in
Case 1 and stable in Case 2.

The bifurcation structure does not change much
locally, if $S \neq 0$. For $(\alpha, \beta) = (0, -1)$, we have
the same critical fixed point and the same eigen-
vectors for the center manifold computations. The
critical values for $a(0)$ and $b(0)$ remain unchanged,
therefore we can distinguish the same Cases 1–3,
depending on $R$. The coefficients $c_i(0)$ are different,
namely:

$$
\begin{align*}
c_1(0) &= \frac{3}{2}(1-R) + 3S, & c_2(0) &= -\frac{1}{2}(1+R) + 3S, \\
c_3(0) &= -\frac{1}{2}(1-R)^2 + 3S, & c_4(0) &= \frac{3}{2}(1+R)^2 + 3S.
\end{align*}
$$

From these expressions we obtain the normalized
critical cubic coefficients:

$$
\begin{align*}
c(0) &= \frac{3}{2}(1-R) + 3S, \\
d(0) &= -\frac{1}{2}(5+3R) + S(3+4R).
\end{align*}
$$

The Lyapunov coefficient also gets an extra term:

$$
c_{\text{NS}} = R^2(1-R) + 2RS(1+2R).
$$

Now, while we still have Case 1 for $-1 < R < 1$, the closed invariant curve might be stable or unstable depending on a combination of $R$ and $S$. Similarly, for Case 2 where $R > 1$, the closed invariant curve need not be stable. We excluded $R = 0$ as exceptional. For $(\alpha, \beta, R) = (0, -1, 0)$ we actually have a codim 3 singularity, since the Lyapunov coefficient is then equal to zero.

We present in Fig. 17 some phase portraits to
show the behavior of (43) close to the fold-flip point
for several values of R and S:

<table>
  <thead>
    <tr>
      <th>Subfigure</th>
      <th>R</th>
      <th>S</th>
      <th>&alpha;</th>
      <th>&beta;</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>a</td>
      <td>-2</td>
      <td>0</td>
      <td>0.0004</td>
      <td>-0.999</td>
    </tr>
    <tr>
      <td>b</td>
      <td>-2</td>
      <td>-5</td>
      <td>0.0004</td>
      <td>-0.999</td>
    </tr>
    <tr>
      <td>c</td>
      <td>0.5</td>
      <td>0</td>
      <td>0.0001</td>
      <td>-1.001</td>
    </tr>
    <tr>
      <td>d</td>
      <td>0.5</td>
      <td>-2</td>
      <td>0.001</td>
      <td>-1.001</td>
    </tr>
    <tr>
      <td>e</td>
      <td>3</td>
      <td>0</td>
      <td>0.0004</td>
      <td>-0.9985</td>
    </tr>
    <tr>
      <td>f</td>
      <td>3</td>
      <td>0.5</td>
      <td>0.001</td>
      <td>-0.9985</td>
    </tr>
  </tbody>
</table>

4.2. The extended Lorenz-84 model

As an example of the fold-flip bifurcation in ODEs
we study an extension of a model formulated by
Lorenz [1984]. A bifurcation analysis of this model
was presented by Shil'nikov et al. [1995] and van
Veen [2003]. The latter paper showed, that the
Lorenz-84 model approximates the dynamics of a
low-order Galerkin truncation of an atmospheric
flow model. The atmospheric model describes syn-
optic dynamics, i.e. dynamics on a time scale of
about a week and a length scale of about ten thou-
sand kilometers. The synoptic atmospheric dynam-
ics over the North Atlantic ocean is dominated by
the jet stream, a westerly circulation, and baroclinic
waves, which transport heat and momentum north-
ward. In the Lorenz-84 model the intensity of the
jet stream is given by X and the sine and cosine
coefficients of a baroclinic wave are given by Y and
Z. The dynamical equations for these variables are the first three equations in the system

$$
\begin{equation}
\left\{
\begin{aligned}
\dot{X} &= -Y^2 - Z^2 - \alpha X + \alpha F - \gamma U^2 \\
\dot{Y} &= XY - \beta XZ - Y + G \\
\dot{Z} &= \beta XY + XZ - Z \\
\dot{U} &= -\delta U + \gamma UX + T
\end{aligned}
\right.
\tag{44}
\end{equation}
$$

The damping time scale of the baroclinic wave is about one week and is scaled to unity. As $\alpha = 1/4$ the jet stream is damped more slowly. We extend the Lorenz-84 model by adding the fourth equation in the spirit of Palmer [1995], who studied the influence on the jet stream and the baroclinic waves of external parameters such as the sea surface temperature. Note, that $U$ interacts nonlinearly with the jet stream and that the Lyapunov

function $L = X^2 + Y^2 + Z^2 + U^2$ is conserved in the absence of linear damping and constant forcing. In the following we will fix

$$
\beta = 1, \quad \gamma = 0.987, \quad \delta = 1.04, \quad G = 0.2.
$$

It is known from [van Veen, 2003], that the basic cycle of the Lorenz-84 model, which is created via a Hopf bifurcation of the trivial equilibrium and represents a traveling baroclinic wave, undergoes a period-doubling cascade at certain values of the parameter $F$. By construction, solutions of the Lorenz-84 model are also solutions of the extended model for $U = T = 0$. At a period-doubling bifurcation of the Lorenz-84 model a small perturbation of the extended model can yield a cycle with Floquet multipliers $+1$ and $-1$.

<table><caption>Table 1. The critical fixed point of the Poincaré map associated to (44) with the parameter values and the multipliers.</caption>
  <tr>
    <td>Y = -5.6335141581943.10<sup>-2</sup></td>
    <td>F = 1.7620532879639</td>
    <td>&lambda;<sub>1</sub> = 1 &plusmn; 10<sup>-11</sup></td>
  </tr>
  <tr>
    <td>Z = 4.1293337647981.10<sup>-2</sup></td>
    <td>T = 2.80597685.10<sup>-4</sup></td>
    <td>&lambda;<sub>2</sub> = -1 &plusmn; 10<sup>-9</sup></td>
  </tr>
  <tr>
    <td>U = 0.31352886978279</td>
    <td></td>
    <td>&lambda;<sub>3</sub> = -0.43054026152942</td>
  </tr>
</table>

We study the Poincaré map in the plane $X = 1.05$ at the fold-flip point. In Table 1 the numerical values are listed. Figure 18 shows the bifurcation diagram obtained in a neighborhood of the codim 2 point. On the top left corner there is a generalized Hopf point $GH$ at which a Hopf bifurcation and a fold bifurcation of a cycle meet (see e.g. [Kuznetsov, 1998], Chap. 8.3). Along the fold line the cycle has one multiplier $+1$ and two multipliers within the unit circle. One of them crosses $-1$ at the codim 2 fold-flip point $A$. From the picture and its scale we deduce that $a(0), b(0) > 0$. Secondly we expect that $a(0)$ is small, since the fold and flip curves are very close. To compute the multilin-ear functions, we integrated the variational equa-tions described in Appendix A numerically using a Runge-Kutta-Felbergh scheme of 7-8 order. Then we implemented the formulas for the Poincaré map and its derivatives in MAPLE and applied the for-mulas of Sec. 2.2 for the center manifold reduction. This gives

$$
\begin{align*}
a(0) &= 0.002047, \quad b(0) = 4.4010, \\
c(0) &= -0.02336, \quad d(0) = 232.682.
\end{align*}
$$

These values are in agreement with what we deduced from Fig. 18. For these coefficients the Lyapunov coefficient has the value $c_{NS} = -0.6062$, which indicates that there is a stable invariant curve in the center manifold. A configuration of the stable invariant circles of the Poincaré map, obtained by forward integration, is shown in Fig. 19. A double torus around the period-2 cycle corresponds to the invariant circles of the second iterate of the Poincaré map. Also shown are the intersection points of the unstable cycle that bifurcates on the torus line TR. The period of this cycle is close to that of the motion along the stable torus. The other period of stable solutions on the torus is much longer and in fact goes to infinity at the fold-flip point $A$. Away from TR the torus breaks up and a strange attractor is created in a thin tube around the original invariant circle as shown in Fig. 20.

With $a(0), b(0) > 0$ we expect the torus to become heteroclinic as can be seen from the unfolding in Fig. 7. However, if we calculate the linear approximation of the heteroclinic bifurcation line $J_-$ from Proposition 3.2.2, it turns out to lie extremely close to the period-doubling line PD in Fig. 18. Hence,
<figure><img src="image_19.png" /><figcaption>Fig. 18. Bifurcation diagram near the codim 2 point: Hopf (black), period-doubling (red), fold of periodic orbits (blue), torus bifurcation (darker blue). Obtained using AUTO97 [Doedel et al., 1977].</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. The invariant circles of the second iterate of the Poincaré map to the right of *TR*. The crosses indicate the intersection points of the unstable cycle which bifurcates on *TR*.</figcaption></figure>

we cannot find heteroclinic tangencies numerically. In the Poincaré section of Fig. 19, corresponding to phase portrait 4– in Fig. 8, we do see that the torus is squeezed and looks like the heteroclinic surface that exists on $J_-$.

Summarizing, the stable solutions around the codim 2 point *A* are

* an equilibrium which represents a steady jet stream without wave activity,

* period 1 and 2 cycles which represent traveling baroclinic waves, and are also present in the Lorenz-84 model itself,

* a stable torus, which represents a traveling baroclinic wave with an amplitude that is slowly modulated or

* a strange attractor, which represents a traveling baroclinic wave with an amplitude that is modulated irregularly.
<figure><img src="image_21.png" /><figcaption>Fig. 20. After the destruction of the torus a strange attractor is created. Also shown are the intersection points of the saddle type cycle. The integration was started near the saddle type orbit.</figcaption></figure>

The modulation of the amplitude of the traveling wave is due to interaction with the added mode, $U$.

## 5. Discussion

In this paper we have contributed to the analysis of a codim 2 bifurcation occurring in generic two-parameter families of maps, namely the fold-flip bifurcation. We have applied our results to two examples: a generalized Hénon map and an extended Lorenz-84 ODE. Note that this bifurcation has been found in the period map describing dynamics of a simple robotic arm by Steindl [1990], and recently in a periodically forced Lorenz-84 model by Broer et al. [2002]. In that case the period map is considered as well.

There are open numerical problems appearing while applying the theory of codim 2 bifurcations of fixed points to limit cycles in multi-dimensional ODEs. The approach based on numerical integration of the higher-order variational equations, that we have applied in Sec. 4.2, has obvious limitations. Although it works well for low-dimensional and nonstiff ODEs, a more robust approach to study codim 2 bifurcations of limit cycles has to be developed. Such an approach might combine the center manifold reduction near the cycle with a periodic normalization of the ODE on the center manifold. For first theoretical advances in this direction, see [Iooss, 1988] and [Chow & Wangm, 1986].

Finally, let us note that the bifurcation behavior of the generalized Hénon map (43) is far from being completely understood. Although the complete bifurcation diagram is likely to be never

constructed, its essential global features can be analyzed using existing analytical and numerical tools. In particular, it would be interesting to investigate how the known bifurcation structures existing for $\beta > 0$ match those emanating from the fold-flip point existing in the half-plane $\beta < 0$.

## References

* Arnold, V. I. [1983] *Geometrical Methods in the Theory of Ordinary Differential Equations* (Springer-Verlag, NY).

* Arrowsmith, D. & Place, C. [1990] *An Introduction to Dynamical Systems* (Cambridge University Press).

* Back, A., Guckenheimer, J., Myers, M. R., Wicklin, F. J. & Worfolk, P. A. [1992] “DSTool: Computer assisted exploration of dynamical systems,” *Notices Amer. Math. Soc.* **39**, 303.

* Beyn, W.-J., Champneys, A., Doedel, E., Govaerts, W., Kuznetsov, Yu. A. & Sandstede, B. [2002] *Handbook of Dynamical Systems, Vol. 2*, Numerical Continuation, and Computation of Normal Forms, ed. Fiedler, B. (Elsevier Science, Amsterdam).

* Broer, H. W., Simó, C. & Vitolo, R. [2002] “Bifurcations and strange attractors in the Lorenz-84 climate model with seasonal forcing,” *Nonlinearity* **15**, 1205–1267.

* Chow, S.-N. & Wang, D. [1986] “Normal forms of bifurcating periodic orbits,” *Contemp. Math.* **56**, 9–18.

* Chow, S.-N., Li, C. & Wang, D. [1994] *Normal Forms and Bifurcations of Planar Vector Fields* (Cambridge University Press).

* Doedel, E. J., Champneys, A. R., Fairgrieve, T. F., Kuznetsov, Yu. A., Sandstede, B. & Wang, X. J. [1997] *AUTO97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont)*, Computer Science, Concordia University, Montreal, Canada.
Gavrilov, N. K. & Shil'nikov, L. P. [1972] "On three-dimensional systems close to systems with a structurally unstable homoclinic curve. I," *Math. USSR Sb.* **17**, 467-485.

Gavrilov, N. K. & Shil'nikov, L. P. [1973] "On three-dimensional systems close to systems with a structurally unstable homoclinic curve. II," *Math. USSR Sb.* **19**, 139-156.

Gheiner, J. [1994] "Codimension-two reflection and non-hyperbolic invariant lines," *Nonlinearity* **7**, 109-184.

Gonchenko, S. V., Shil'nikov, L. P. & Turaev, D. V. [1996] "Dynamical phenomena in systems with structurally unstable Poincaré homoclinic orbits," *Int. J. Bifurcation and Chaos* **6**, 1-17.

Gonchenko, S. V., Gonchenko, V. S. & Tatjer, J. C. [2001] "Three-dimensional dissipative diffeomorphisms with codimension two homoclinic tangencies and generalized Hénon maps," in *Proc. Int. Conf. 'Progress in Nonlinear Science' Dedicated to 100th Anniversary of A. A. Andronov*, pp. 63-79.

Gonchenko, S. V. & Gonchenko, V. S. [2002] On Andronov-Hopf bifurcations of two-dimensional diffeomorphisms with homoclinic tangencies," preprint No. 556, WIAAS, Berlin.

Gonchenko, V. S. [2002] "On bifurcations of two-dimensional diffeomorphisms with homoclinic tangency of manifolds of a 'neutral' saddle," *Proc. Stklov. Inst. Math.* **236**, 95-102 (in Russian).

Guckenheimer, J. & Holmes, P. [1983/2002] *Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields* (Springer-Verlag, NY).

Iooss, G. [1988] "Global characterization of the normal form for a vector field near a closed orbit," *J. Diff. Eqs.* **76**, 47-76.

Krauskopf, B. & Osinga, H. [1998a] "Globalizing two-dimensional unstable manifolds of maps," *Int. J. Bifurcation and Chaos* **8**, 483-503.

Krauskopf, B. & Osinga, H. [1998b] "Growing 1D and quasi-2D unstable manifolds of maps," *J. Comput. Phys.* **146**, 404-419.

Kuznetsov, Yu. A. [1998] *Elements of Applied Bifurcation Theory*, 2nd edition (Springer-Verlag, NY).

Kuznetsov, Yu. A. [1999] "Numerical normalization techniques for all codim 2 bifurcations of equilibria in ODE's," *SIAM J. Numer. Anal.* **36**, 1104-1124.

Lorenz, E. N. [1984] "Irregularity: a fundamental property of the atmosphere," *Tellus A36*, 98-110.

Palmer, T. N. [1995] "The influence of north-west Atlantic sea surface temperature: An unplanned experiment," *Weather* **30**, 413-419.

Shil'nikov, A. L., Nicolis, G. & Nicolis, C. [1995] "Bifurcation and predictability analysis of a low-order atmospheric circulation model," *Int. J. Bifurcation and Chaos* **5**, 1701-1711.

Simó, C. [1989] "On the analytical and numerical approximation of invariant manifolds," *Modern Methods in Celestial Mechanics*, Proc. 13th Spring School on Astrophysics in Goutelas, eds. Benest, D. & Froeschlé, C. (Edition Frontières, Gift-sur Yvette, France); available via http://www.maia.ub.es/dsg/2004/index.html.

Steindl, A. [1990] "Bifurcation of codimension 2 for a discrete map," in *Continuation and Bifurcations: Numerical Techniques and Applications* (Leuven, 1989), NATO Adv. Sci. Inst. Ser. C Math. Phys. Sci., Vol. 313, eds. Roosevelt, D., De Dier, B. & Spence, A. (Kluwer Acad. Publ., Dordrecht), pp. 389-396.

Takens, F. [1974] "Forced oscillations and bifurcations," *Comm. Math. Inst. Rijksuniversiteit Utrecht* **3**, 1-59.

Veen, L. van [2003] "Baroclinic flow and the Lorenz-84 model," *Int. J. Bifurcation and Chaos* **13**, 2117-2139.

## Appendix A

### Poincaré Maps and Their Derivatives

In this appendix we present a method to compute numerically the derivatives of a Poincaré map up to and including third order. Let $f(x)$ be a smooth vector field in $\mathbb{R}^n$. We take a local cross-section $\Sigma \subset \mathbb{R}^n$ with dimension $n-1$. The hypersurface $\Sigma$ does not need to be a coordinate plane, but should be chosen transversal to the flow of $f$. This is satisfied if for the normal $n_\Sigma(x)$ for $x \in \Sigma$ we have $\langle f(x), n_\Sigma(x) \rangle \neq 0$. Let $L_0$ be a $\tau_0$-periodic orbit through $\Sigma$ and $U \subset \Sigma$ a subset which contains an intersection point $x_0$ of $L_0 \cap \Sigma$. If $L_0$ intersects $\Sigma$ in multiple points we shrink $U$ until we have one point of intersection, i.e. $x_0$. The Poincaré map $P: U \to \Sigma$ is defined for $x \in U$ by

$$P(x) = \phi(t(x), x),$$

where $t(x)$ is the arrival time after which the orbit intersects $\Sigma$ for the first time again, and $\phi(t, x)$ is the solution to

$$\dot{x} = f(x), \quad x \in \mathbb{R}^n,$$

with the initial condition $\phi(0, x) = x$ (see Fig. 21). Note that $\phi$ is as smooth as $f$. The arrival time $\tau$ depends on $x$ but $t(x_0) = \tau_0$. The Floquet multipliers of $L_0$ can be calculated as the eigenvalues of the monodromy matrix $\partial\phi(\tau_0, x)/\partial x$. It always has a trivial eigenvalue 1. With the Liouville Theorem one can prove that the product of the multipliers is always positive.

We proceed with a method to compute the Poincaré map and its derivatives numerically. We closely follow Simó [1989], but we extend the method to the third-order derivatives. We can write for a displacement $x+h \in \mathbb{R}^n$
<figure><img src="image_22.png" /><figcaption>Fig. 21. The geometry of the Poincaré map.</figcaption></figure>

$$ \phi(t, x+h) = \phi(t, x) + \frac{\partial\phi(t, x)}{\partial x}h + \frac{1}{2}\frac{\partial^2\phi(t, x)}{\partial x^2}h^2 + \dots $$

We assume that $(x + h) \in \Sigma$, or, more precisely, $\partial\phi(t, x)/(\partial x)h = (D_x\phi(t, x))(h)$, $\partial^2\phi(t, x)/(\partial x^2)h^2 = (D_x^2\phi(t, x))(h, h)$, etc. The multilinear functions $\partial^i\phi(t, x)/\partial x^i$ satisfy the variational equations, which we find by successive differentiation of $\dot{\phi}(t, x) = f(\phi(t, x))$. Since $f$ and $\phi$ are smooth we can change the order of the derivatives and find

$$ \begin{aligned} \frac{d}{dt} \frac{\partial \phi(t, x)}{\partial x} &= Df(\phi(t, x)) \frac{\partial \phi(t, x)}{\partial x}, \\ \frac{d}{dt} \frac{\partial^2 \phi(t, x)}{\partial x^2} &= D^2 f(\phi(t, x)) \left( \frac{\partial \phi(t, x)}{\partial x} \right)^2 \\ &\quad + Df(\phi(t, x)) \frac{\partial^2 \phi(t, x)}{\partial x^2}, \\ \frac{d}{dt} \frac{\partial^3 \phi(t, x)}{\partial x^3} &= D^3 f(\phi(t, x)) \left( \frac{\partial \phi(t, x)}{\partial x} \right)^3 \\ &\quad + 3D^2 f(\phi(t, x)) \left( \frac{\partial^2 \phi(t, x)}{\partial x^2}, \frac{\partial \phi(t, x)}{\partial x} \right) + Df(\phi(t, x)) \frac{\partial^3 \phi(t, x)}{\partial x^3}. \end{aligned} $$

The initial conditions are given by

$$ \begin{aligned} \left. \frac{\partial \phi(t, x)}{\partial x} \right|_{t=0} &= I_n, \quad \left. \frac{\partial^2 \phi(t, x)}{\partial x^2} \right|_{t=0} = 0, \\ \text{and} \quad \left. \frac{\partial^3 \phi(t, x)}{\partial x^3} \right|_{t=0} &= 0. \end{aligned} $$

To find an analytic solution to the variational equations is an exception, but numerically we can integrate these equations and we will show how to use

them to compute the derivatives of the Poincaré map. Normally one considers a cross-section $\Sigma$ and the Poincaré map is just the return map to $\Sigma$. This is the case when one looks at periodic orbits, which account for fixed points of the Poincaré map. However, we will set up the structure of the Poincaré map for a flow from a cross-section $\Sigma_1$ to a cross-section $\Sigma_2$ to distinguish the initial and final cross-section. In this way the notation is clearer, but in the end we will set $\Sigma_1 = \Sigma_2$. The sections are defined by equations $g_1(x) = 0$ and $g_2(x) = 0$, respectively and we assume both sections to be transversal to the flow.

Let $x_1 \in \Sigma_1$ be the initial point and we define $P: \Sigma_1 \to \Sigma_2$ by $P(x_1) = x_2 = \phi(t(x_1), x_1)$, such that $g_2(\phi(t(x_1), x_1)) = 0$. Here $t(x_1)$ is the travel time and it depends, of course, on the initial point. We compute the first derivative of $P$ by differentiation with respect to $x_1$

$$ \frac{\partial P}{\partial x_1} = f(x_2) \frac{\partial t}{\partial x_1} + \frac{\partial \phi}{\partial x_1}(t(x_1), x_1). \quad (A.1) $$

**Remark A.0.1.** We add a word of caution. We should make a distinction between the variables of the ODE in $\mathbb{R}^n$ and the variables of the Poincaré map in $\Sigma \sim \mathbb{R}^{n-1}$. The derivatives above are formal w.r.t. $x_1 \in \mathbb{R}^n$, but the LHS of (A.1) has the domain $\Sigma_1$. Therefore the RHS has to be restricted to $\Sigma_1$. Now, if we take $\Sigma_1 = \Sigma_2$ to be given by setting the nth coordinate equal to a constant, then the derivatives of the Poincaré map will contain a block of zeros, namely for the nth component. This result follows easily from (A.2) and (A.5) for the linear and quadratic part and can be extended for higher orders. In this way we can take the variables $\{x_1, \dots, x_n\}$ of $f$ in $\mathbb{R}^n$ as variables for the Poincaré map, excluding the nth.

The first variational equation provides the matrix $\partial\phi/\partial x_1$. The derivative $\partial t/\partial x_1$ is still unknown. We obtain it by differentiating the relation $g_2(\phi(t(x_1), x_1)) = 0$

$$ \left\langle Dg_2(x_2), \left( f(x_2) \frac{\partial t}{\partial x_1} + \frac{\partial \phi}{\partial x_1}(t(x_1), x_1) \right) \right\rangle = 0, $$

or

$$ \frac{\partial t}{\partial x_1} = - \frac{1}{\langle \nabla g_2(x_2), f(x_2) \rangle} Dg_2 \frac{\partial \phi}{\partial x_1}. $$

Note that the transversality condition implies that the denominator is nonzero. Now we suppose that $\Sigma_1$ and $\Sigma_2$ are just coordinate planes, i.e. $g_1$ and $g_2$ are given by taking, for instance the last, a
coordinate equal to a constant. We write $a_{ij} = \partial\phi_i/\partial x_{1,j}(t(x_1), x_1)$. Then, we have $\partial t/\partial x_{1,i} = -a_{ni}/f_n(x_2)$ and in coordinates the $(n-1) \times (n-1)$ matrix $\partial P/\partial x_1(x_1)$ is given by

$$ \left( \frac{\partial P}{\partial x_1} \right)_{ij} = a_{ij} - \frac{f_i(x_2)}{f_n(x_2)} a_{nj}. \quad (A.2) $$

The restriction to the first $n-1$ components is given by taking $1 \le i, j \le n-1$.

We continue with the second derivative. We differentiate another time and find

$$ \frac{\partial^2 P}{\partial x_1^2} = Df(x_2)f(x_2) \left(\frac{\partial t}{\partial x_1}\right)^2 + 2Df(x_2) \frac{\partial \phi}{\partial x_1} \frac{\partial t}{\partial x_1} + f(x_2) \frac{\partial^2 t}{\partial x_1^2} + \frac{\partial^2 \phi}{\partial x_1^2}. $$

Note that the first variational equation is used here to simplify the derivatives. As before $\partial^2\phi/\partial x_1^2$ is obtained from the second variational equation and $\partial^2 t/\partial x_1^2$ by differentiating $g(\phi(t(x_1), x_1))$ once more

$$ D^2 g_2(x_2) \left( f(x_2) \frac{\partial t}{\partial x_1} + \frac{\partial \phi}{\partial x_1} \right)^2 + Dg_2(x_2) \left[ Df(x_2) f(x_2) \left( \frac{\partial t}{\partial x_1} \right)^2 + 2Df(x_2) \frac{\partial \phi}{\partial x_1} \frac{\partial t}{\partial x_1} + f(x_2) \frac{\partial^2 t}{\partial x_1^2} + \frac{\partial^2 \phi}{\partial x_1^2} \right] = 0, \quad (A.3) $$

or more compactly (using (A.1) and (A.3))

$$ D^2 g_2(x_2) \left( \frac{\partial P}{\partial x_1} \right)^2 + Dg_2(x_2) \frac{\partial^2 P}{\partial x_1^2} = 0. \quad (A.4) $$

For notational clarity we drop from now on the subscript 1 in $x_1$. All derivatives are now with respect to the (the components $x_j$ of) $x$ or to $t$. Returning to the specific case of a return map we introduce the notation $b_{ijk} = \partial^2\phi_i/\partial x_j\partial x_k$, $f_i = f_i(x_2)$ and $f_{ij} = D_i f_j(x_2)$. Making these substitutions also for the time derivatives, we obtain

$$ \frac{\partial^2 P_i}{\partial x_j \partial x_k} = b_{ijk} - \frac{f_i}{f_n} b_{njk} + \frac{1}{f_n} \sum_{s=1}^{n} \left( f_{si} - \frac{f_i}{f_n} f_{sn} \right) \times \left( \frac{f_s}{f_n} a_{nj} a_{nk} - a_{sk} a_{nj} - a_{sj} a_{nk} \right). \quad (A.5) $$

As before we have $1 \le i, j, k \le n-1$.

Finally, we do the third derivative as well. Differentiating one more time, we find

$$ \frac{\partial^3 P}{\partial x^3} = D^2 f(x_2) \left( f^2(x_2) \left( \frac{\partial t}{\partial x} \right)^3 + 3f(x_2) \frac{\partial \phi}{\partial x} \left( \frac{\partial t}{\partial x} \right)^2 + 3 \left( \frac{\partial \phi}{\partial x} \right)^2 \frac{\partial t}{\partial x} \right) + Df(x_2)^2 \left( f(x_2) \left( \frac{\partial t}{\partial x} \right)^3 + 3 \frac{\partial \phi}{\partial x} \left( \frac{\partial t}{\partial x} \right)^2 \right) + f(x_2) \frac{\partial^3 t}{\partial x^3} + \frac{\partial^3 \phi}{\partial x^3} + 3Df(x_2) \times \left( f(x_2) \frac{\partial t}{\partial x} \frac{\partial^2 t}{\partial x^2} + \frac{\partial \phi}{\partial x} \frac{\partial^2 t}{\partial x^2} + \frac{\partial \phi}{\partial x^2} \frac{\partial t}{\partial x} \right). $$

More explicitly in coordinates we have

$$ \begin{aligned} \frac{\partial^3 P_i}{\partial x_j \partial x_k \partial x_l} = & \sum_{r,s=1}^{n} D_{rs} f_i(x_2) \left( f_r f_s \frac{\partial t}{\partial x_j} \frac{\partial t}{\partial x_k} \frac{\partial t}{\partial x_l} + f_s \frac{\partial \phi_r}{\partial x_j} \frac{\partial t}{\partial x_k} \frac{\partial t}{\partial x_l} + \frac{\partial \phi_r}{\partial x_j} \frac{\partial \phi_s}{\partial x_k} \frac{\partial t}{\partial x_l} \right) \\ & + \sum_{r,s=1}^{n} D_s f_i(x_2) D_r f_s(x_2) \left( f_r \frac{\partial t}{\partial x_j} \frac{\partial t}{\partial x_k} \frac{\partial t}{\partial x_l} + \frac{\partial \phi_r}{\partial x_j} \frac{\partial t}{\partial x_k} \frac{\partial t}{\partial x_l} \right) \\ & + \sum_{s=1}^{n} D_s f_i(x_2) \left( f_s(x_2) \frac{\partial t}{\partial x_j} \frac{\partial^2 t}{\partial x_k \partial x_l} + \frac{\partial \phi_s}{\partial x_j} \frac{\partial^2 t}{\partial x_k \partial x_l} + \frac{\partial^2 \phi_s}{\partial x_j \partial x_k} \frac{\partial t}{\partial x_l} \right) \\ & + f_i \frac{\partial^3 t}{\partial x_j \partial x_k \partial x_l} + \frac{\partial^3 \phi_i}{\partial x_j \partial x_k \partial x_l}. \end{aligned} \quad (A.6) $$
Here the * means that the term with j, k, l cyclically permuted should be included as well. The above expression is then invariant under changing the order of the differentiation, if f is smooth. The third-order derivative of the time can be found from

$$D^3 g_2(x_2) \left(\frac{\partial P}{\partial x}\right)^3 + 3D^2 g_2^2(x_2) \frac{\partial^2 P}{\partial x^2} \frac{\partial P}{\partial x} \\ + Dg_2(x_2) \frac{\partial^3 P}{\partial x^3} = 0.$$