Existence and implications of a pitchfork-Hopf bifurcation in a continuous-time two-sector growth model

Abstract

This paper shows that the dynamics of the Lucas (J Monet Econ, 22:3–42, 1988) endogenous growth model with flow externalities may give rise to a 2-torus, a compact three-dimensional manifold enclosed by a two-dimensional surface. The implications of this result are relevant for many fields of economic theory. It is first of all clear that if we choose to initialize the dynamics in the basin of attraction of this trapping region, a continuum of perfect foresight solutions may be observed. A simple econometric exercise, linking the physical-to-human capital ratio (state variable) to the 5-years forward variance of the growth rate of an unbalanced sample of 183 countries, seems to provide empirical backing for the phenomenon. Other important consequences, relevant from the point of view of endogenous cycles theory, are also scrutinized in the paper.

1 Introduction

A large body of literature has clearly established that the equilibrium conditions of two-sector continuous-time endogenous growth models with market imperfections do not always determine a unique perfect foresight path. If the number of stable roots of the linearization matrix (evaluated at a hyperbolic stationary point, which may correspond to a balanced growth path) exceeds the number of state (predetermined) variables, then a continuum of perfect foresight paths in a close neighborhood of the stationary point exists. If this happens, the equilibrium is said to be locally indeterminate. More recently, researchers have shown that the indeterminacy property is more pervasive than originally thought in the sense that a continuum of perfect foresight solutions of a dynamic problem may also be observed when the system is initialized far away from a steady state. This is the idea of global indeterminacy, which dates back, at least, to the seminal works of Krugman (1991) and Matsuyama (1991).

To discuss this literature, let us recall the dynamic problem associated with the following variant of the Lucas (1988) model:

$$ \begin{gathered} ~\mathop {\max }\limits_{{C\left( t \right),~u\left( t \right)}} \mathop \smallint \limits_{0}^{\infty } \frac{{C^{{1 - \sigma }} - 1}}{{1 - \sigma }}e^{{ - \rho t}} dt\,\,\,\left( {\mathcal{P}} \right) \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,{\text{subject}}\,\,\,{\text{to}} \hfill \\ \end{gathered} $$

$$ \dot{K} = AK^{\beta } \left( {uH} \right)^{{1 - \beta }} H_{a}^{\gamma } - C $$

$$ \dot{H} = \delta H\left( {1 - u} \right)\left( {1 - u} \right)_{a}^{\eta } $$

$$ \begin{gathered} K\left( 0 \right) = K_{0} \hfill \\ H\left( 0 \right) = H_{0} \hfill \\ \end{gathered} $$

which is extended to account for the presence of flow externalities, as in Chamley (1993), that depend on human capital interactions across individuals (cf. Mattana et al. 2009; Bella and Mattana 2014). The notations are as follows. $C$ is consumption, $K$ is physical capital, $H$ is human capital, and $u$ is the fraction of labor allocated to the production of physical capital. The feasibility requirement is $0 \le u \le 1$. $H_{a}$ and $(1 - u)_{a}$ are externality factors. $\left( {A,\beta ,\delta ,\gamma ,\eta ,\rho ,\sigma } \right)$ is the set of parameters. $\left( {A,\delta } \right) \in \mathbb{R}_{{ + + }}^{2}$ measure the technological level in the physical and human capital sectors, respectively. $\left( {\gamma ,\eta } \right) \in \mathbb{R}_{ + }^{2}$ are numerical factors that tune the effects of the externalities. Finally, $\beta \in \left( {0,1} \right)$ is the share of physical capital in the goods sector, $\rho \in \mathbb{R}_{ + }$ is the time preference rate, and $\sigma \in \mathbb{R}_{ + } - \left\{ 1 \right\}$ is the inverse of the intertemporal elasticity of substitution.

The market solution of problem ${\mathcal{P}}$ implies the following non-linear three-dimensional system of first-order differential equations:

$$ \begin{aligned} \dot{X} & = AX^{\beta } u^{{1 - \beta }} + \frac{{\delta \left( {1 - \beta + \gamma } \right)}}{{\beta - 1}}X(1 - u)^{{1 + \eta }} - QX \\ \dot{u} & = \frac{{u\left( {1 - u} \right)}}{{\beta \left( {1 - u} \right) - \eta u}}\left[ {\delta (1 - u)^{\eta } - \delta \left( {\beta - \gamma } \right)(1 - u)^{{1 + \eta }} - \beta Q} \right]\,\,\,\,\,\,\,\,\left( {\mathcal{S}} \right) \\ \dot{Q} & = - \frac{\rho }{\sigma }Q + \frac{{\beta - \sigma }}{\sigma }AX^{{\beta - 1}} u^{{1 - \beta }} Q + Q^{2} \\ \end{aligned} $$

in the stationary variables $X \equiv K^{{\frac{{1 - \beta + \gamma }}{{1 - \beta }}}} H^{{ - 1}}$, $u$ and $Q \equiv \frac{C}{K}$ (cf. Mattana et al. 2009, for the derivation). To prove global indeterminacy, scholars have typically searched for conditions in the parameter space that imply critical eigenvalues. There are in fact well-established theoretical results about the global properties of systems of ODEs that are close to linear degeneracies. A clear picture can be obtained by discussing the unfolding of the hypernormal form associated with system ${\mathcal{S}}$.¹ The underlying idea behind hypernormal form theory is to use near-identity transformations to remove the terms that are not essential in the analytical expression of the vector field for a complete understanding of the local dynamic behavior. Of course, by topological equivalence, all results referring to the hypernormal form can be carried on to the original system.

Let

$$ \left( {\begin{array}{*{20}c} {\dot{w}_{1} } \\ {\dot{w}_{2} } \\ {\dot{w}_{3} } \\ \end{array} } \right) = \left[ {\begin{array}{*{20}c} 0 & 1 & 0 \\ 0 & 0 & 1 \\ {\varepsilon _{1} } & {\varepsilon _{2} } & {\varepsilon _{3} } \\ \end{array} } \right]\left( {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {w_{3} } \\ \end{array} } \right) + \left( {\begin{array}{*{20}c} 0 \\ 0 \\ {F_{3} \left( {w_{1} ,w_{2} ,w_{3} } \right)} \\ \end{array} } \right)\,\,\,\,\,\,\left( {\mathcal{N}} \right) $$

be the hypernormal transformation of a generic 3-dimensional system of ODEs. In which ${\mathcal{N}}$ and $w_{i}$,$i = 1,2,3$ represent the new coordinates and

$$ F_{3} \left( {w_{1} ,w_{2} ,w_{3} } \right) = \mathop \sum \limits_{{k \ge 2}} \left\{ {\mathop {\mathop \sum \limits^{m} }\limits_{{j = 0}} \left( {a_{j}^{k} w_{1}^{{k - j - 1}} w_{3}^{{j + 1}} + b_{j}^{k} w_{1}^{{k - 2j - 1}} w_{2}^{{2j + 1}} } \right) + \mathop {\mathop \sum \limits^{m} }\limits_{{j = 0}} c_{j}^{k} w_{1}^{{k - 2j}} w_{2}^{{2j}} } \right\} $$

is the reduced non-linear part up to order $m$.² In $F_{3} \left( {w_{1} ,w_{2} ,w_{3} } \right)$, $a_{j}$, $b_{j}$ and $c_{j}$ are coefficients that depend on the original parameters of the model, and $k$ denotes the polynomial degree. $\varepsilon _{1}$, $\varepsilon _{2}$ and $\varepsilon _{3}$ are the unfolding parameters. The origin is stable when $\varepsilon _{2} < 0$, $\varepsilon _{3} < 0$ and $0 < \varepsilon _{1} < \varepsilon _{2} \varepsilon _{3}$. Additionally, the model undergoes a pitchfork bifurcation at $\varepsilon _{1} = 0$.³ The following properties are more relevant from the perspective of global indeterminacy. Furthermore, ${\mathcal{N}}$ undergoes the following non-linear degeneracies:

• codimension 1 bifurcation of the Hopf type at $\varepsilon _{1} = \varepsilon _{2} \varepsilon _{3}$; that is, ${\mathbf{J}}$ at the origin has one real eigenvalue equal to $\varepsilon _{3} \ne 0$ and two purely imaginary eigenvalues equal to ($\pm i\sqrt { - \varepsilon _{2} }$);

• codimension 2 bifurcation for $\varepsilon _{1} = \varepsilon _{2} = 0$; that is, ${\mathbf{J}}$ at the origin has a double-zero eigenvalue and a third real eigenvalue equal to $\varepsilon _{3} \ne 0$;

• codimension 2 pitchfork–Hopf interaction for $\varepsilon _{1} = \varepsilon _{3} = 0$ and $\varepsilon _{2} < 0$; that is, ${\mathbf{J}}$ at the origin has one zero eigenvalue and two complex eigenvalues given by $\pm i\sqrt { - \varepsilon _{2} }$, where $\varepsilon _{2} < 0$.

Mattana and Venturi (1999), as well as Nishimura and Shigoka (2006), exploited a singularity of type (i) to show the possibility of Hopf cycles for system ${\mathcal{S}}$ for the case of $\eta = 0$. Then, given the initial condition $X\left( 0 \right)$, there exists a continuum of initial values of the jump variables located inside the Hopf orbit that gives rise to perfect foresight equilibria. Since a region bounded by a Hopf orbit may exceed the small neighborhood of the local analysis, the equilibrium can be considered globally indeterminate. The properties of a singularity of type (ii) have been used by Mattana et al. (2009) and by Bella and Mattana (2014) to achieve global indeterminacy results for a variant of the model with $\gamma = 0$. In the former contribution, the Kopell and Howard (1975) theorem is used to show that there are regions of the parameter space in which there exists either a closed loop or a homoclinic orbit in a well-located reduced manifold. Again, given $X\left( 0 \right)$, the possibility of a continuum of valid initial values of the jump variables located inside the closed loop of the homoclinic orbit implies global indeterminacy. Bella and Mattana (2014) derived similar results by means of the Bogdanov-Takens bifurcation theorem. Implications for global indeterminacy of case (iii) remain, to our knowledge, totally unexplored. However, as shown, for instance, in Kuznetsov (2004), Wiggins (1991), and Guckenheimer and Holmes (1983), the simultaneous occurrence of both linear degeneracies establishes the possibility of interesting three-dimensional dynamic behavior. Of particular interest for this paper is the possibility that system ${\mathcal{S}}$ gives rise to a 2-torus, a compact three-dimensional manifold enclosed by a two-dimensional surface. Then, if the initial state of the economy belongs to this compact set, there exists, outside the small neighborhood relevant for the local analysis, a continuum of initial values of the jump variables, which give rise to recurrent orbits.⁴

The investigation of the conditions under which the uniqueness of an equilibrium is not warranted for the solution of problem ${\mathcal{P}}$ is of crucial importance for growth theory, since the Lucas model (1988) and its variants remain the standard benchmarks for ongoing research (cf., inter al., Tsuboi 2018; Bella et al. 2019; Petrakis 2020; Sasakura 2020; Gomez 2021; Neves-Sequeira et al. 2021) and for devising policies (cf., inter al., Bretschger et al. 2017; Broitman and Czamansky 2020).

The onset of a 2-torus around a balanced growth path (BGP) also has significance for the theory of business cycles. Since a standard 2-torus possesses two distinct 1-dimensional cycles (or closed loops), the mere existence of such phenomenon demonstrates that business oscillations may not be efficient responses of rational agents to external technological shocks; instead, they may be the result of inefficient fluctuations in the variables of a system, caused by shocks to self-fulfilling beliefs of households and firms. There are many innovative aspects involved. First, orbits that wrap around a standard 2-torus will never come back exactly onto themselves because of the typical a-periodicity of their trajectories. This is a major improvement in terms of the plausibility of cycles over other types of endogenous oscillating solutions (Hopf cycles in primis), where unperturbed orbits are bound to regularly visit the same coordinates in the phase space. In this regard, to generate a-periodicity of a cycle, recent literature has added some form of stochasticity to the model (cf., inter al. Beaudry et al. 2020). Finally, it is clear that a 2-torus is generically a fast/slow system; namely, it is supposed to generate two-frequency oscillations (conventionally termed drift and gyro oscillations). Moreover, it is also clear that a drift motion on a standard torus may be much slower than a gyro motion. Our paper contributes to this new line of research and presents a mathematical structure that is able to identify the presence of slow and fast oscillating variables and, possibly, the emergence of lag-lead structures that can be inferred from real-world economic data.

The paper develops as follows. In the first section, we present the set of conditions necessary to represent system ${\mathcal{S}}$ in its hypernormal form. In Sect. 2, by exploiting the Gavrilov-Guckenheimer bifurcation theorem, we construct a correspondence with a simple topological-equivalent truncated planar system in cylindrical coordinates. The relevance of this result for business cycle theory is explored in Section three. The fourth section is devoted to the discussion regarding the implications for the indeterminacy of the equilibrium. A brief conclusive section reassesses the main findings of the paper. The Appendix provides the proofs of the main propositions.

2 The onset of a 2-torus

Recall system ${\mathcal{S}}$, and set $\gamma = 0$. Let $\left( {X^{*} ,u^{*} ,Q^{*} } \right)$ be the values of $\left( {X,u,Q} \right)$ such that $\dot{X} = \dot{u} = \dot{Q} = 0$. Then, along the balanced growth path (BGP), the dynamics satisfy the following equations:

$$ X^{*} = u^{*} \left( {\beta A} \right)^{{\frac{1}{{1 - \beta }}}} \left[ {\delta \left( {1 - u^{*} } \right)^{\eta } } \right]^{{\frac{1}{{\beta - 1}}}} $$

(1.1)

$$ 0 = \left( {1 - u^{*} } \right)^{\eta } - \sigma \left( {1 - u^{*} } \right)^{{1 + \eta }} - \frac{\rho }{\delta }~ $$

(1.2)

$$ Q^{*} = \frac{{\delta \left( {1 - u^{*} } \right)^{\eta } }}{\beta }\left[ {1 - \beta \left( {1 - u^{*} } \right)} \right]~ $$

(1.3)

Mattana et al. (2009) used the discount rate $\rho$ as a bifurcation parameter to show that the linearization matrix of the right-hand side of system ${\mathcal{S}}$, evaluated at a steady state, has a zero eigenvalue. More specifically, if

$$ \rho = \tilde{\rho }\left( {\beta ,\delta ,\eta ,\sigma } \right) \equiv \delta \left( {\frac{\eta }{\sigma }} \right)^{\eta } \left( {\frac{1}{{1 + \eta }}} \right)^{{1 + \eta }} $$

(2)

${\mathbf{Det}}\left( {\mathbf{J}} \right)$ in (A1.3) vanishes and ${\mathbf{J}}$ has at least one eigenvalue equal to zero.

Now, consider the formula for ${\mathbf{Tr}}\left( {\mathbf{J}} \right)$ in (A1.2). Then, if

$$ \eta = \tilde{\eta }\left( {\beta ,\sigma } \right) \equiv \frac{{2\left( {\sigma - \beta } \right)}}{{1 + \beta - 2\sigma }} $$

(3)

${\mathbf{Tr}}\left( {\mathbf{J}} \right) = 0$. Then, if $\eta = \mathop \eta$ and $\rho = \mathop \rho$, ${\mathbf{J}}$ has two purely imaginary eigenvalues.

Let us now substitute $\eta = \mathop \eta \left( {\beta ,\sigma } \right)$ and $\rho = \mathop \rho \left( {\beta ,\delta ,\sigma } \right)$ into the formula for ${\mathbf{B}}\left( {\mathbf{J}} \right)$ in (A1.4). We find that Sign $\left( {{\mathbf{B}}\left( {\mathbf{J}} \right)} \right)$ only depends on two left parameters according to the following formula:

$$ {\text{Sign}}\left( {{\mathbf{B}}\left( {\mathbf{J}} \right)} \right) = {\text{Sign}}\left( {\sigma - \frac{{1 - \left( {3\beta - 1} \right) - \beta ^{3} - 5\beta ^{2} + \sqrt {\left( {1 + \beta } \right)\left( {3\beta - 1} \right)\left( {\beta - 1} \right)^{4} } }}{{\beta ^{3} - \beta ^{2} + 11\beta - 3}}} \right) $$

In Fig. 1, we spotlight the region of the parameter space spanned by the pair $\left( {\beta ,\sigma } \right)$, where ${\mathbf{B}}\left( {\mathbf{J}} \right) < 0$ and ${\mathbf{Tr}}\left( {\mathbf{J}} \right) = {\mathbf{Det}}\left( {\mathbf{J}} \right) = 0$ simultaneously. For notational convenience, we shall denote this region as ${\mathbf{\Omega }} \in {\mathbb{R}}^{2}$. As shown in Fig. 1, ${\mathbf{\Omega }} \notin \emptyset$.

Fig. 1

The Ω set in the (β,σ) space

The following formal result follows.

Lemma 1 (Gavrilov-Guckenheimer bifurcation)

Let $\left( {\beta ,\sigma } \right) \in {\mathbf{\Omega }}$ .; furthermore, let $\eta = \mathop \eta \left( {\beta ,\sigma } \right)$ and $\rho = \mathop \rho \left( {\beta ,\delta ,\sigma } \right) = \delta \mathop \rho \left( {\beta ,\sigma } \right)$ . Then, ${\mathbf{J}}$ has one zero eigenvalue and two purely imaginary eigenvalues given by $\pm i\sqrt \omega$ , where $\omega = - {\mathbf{B}}\left( {\mathbf{J}} \right)$ .

The type of linear degeneracy implied by Lemma 1 has important global implications. There exist several formal classifications of these effects. One has been developed in standard textbooks, such as Wiggins (1991), Kuznetsov (2004) and Seydel (1994), where the original system is put in cylindrical coordinates. Then, one of the coordinates is decoupled from the remaining two so that the classification can be derived with phase plane techniques. In some sense, the dynamics in this phase plane can be viewed as an approximation to a Poincaré map of the full three-dimensional normal form.

As shown in Appendix A.3, the decoupled topologically equivalent system has the form

$$ \begin{array}{*{20}r} \hfill {\dot{r} = \mu _{1} r + arz~~~~~~~~~~~~~~~~~~~~~~} \\ \hfill {\dot{z} = \mu _{2} + kr^{2} - bz^{2} ~~~~k = \pm 1} \\ \end{array} ~~~~~\left( {\mathcal{C}} \right) $$

where $\mu _{1} = \frac{{\hat{\varepsilon }_{3} }}{\kappa }$ and $\mu _{2} = \frac{{\hat{\varepsilon }_{1} }}{\kappa }\left| {\frac{b}{a}} \right|$ are the standard unfolding parameters, given $\kappa = \frac{{a^{2} }}{{\left| b \right|}}$ (cf. Harlim and Langford, 2007; Freire et al. 2002).

Wiggins (1991) shows that the onset of a 2-torus embedded in the full original space described in (A3.5) implies a Hopf cycle in the dimensionally reduced planar system ${\mathcal{C}}$. In turn, this topological case occurs in two different parametric circumstances:

• 1. $b = + 1;\,\,\,\,a > 0$ (case II a, b);

• 2. $b = - 1;\,\,\,a > 0$ (case III).

We computed $b$ and $a$ as functions of $\left( {\beta ,\sigma } \right)$. We found that $> 0~\forall$ $\left( {\beta ,\sigma } \right) \in {\mathbf{\Omega }}$. In contrast, $b$ changes signs in subregion ${\mathbf{\Omega }}$. Therefore, consider Fig. 2, where we superimpose the critical $b = 0$ curve onto the ${\mathbf{\Omega }}$ subregion in Fig. 1. In the dark gray area, the parameters imply $b = + 1$ and $a > 0$. Conversely, $b = - 1$ and $a > 0$ in the light gray area.

Fig. 2

The (β,σ) parameters such that b=+1 and b= − 1

Now, let $P \subset {\mathbb{R}}^{2}$ be a small open neighborhood of $0$. Furthermore, let $\nu = \rho - \mathop \rho$ and $\tau = \eta - \mathop \eta$, which are small deviations of the two control parameters from their respective bifurcation values. Then, recalling that $a$ is always positive in the feasible space of parameters, the following statement can be proven.

Proposition 1 (Existence of a 2-torus)

Recall Lemma $1$ . Furthermore, let $\left( {\nu ,\tau } \right) \subset P$ . First, it is assumed that $b = + 1$ . Then, if $\mu _{1} \left( {\nu ,\tau } \right) = 0$ and $\mu _{2} \left( {\nu ,\tau } \right) < 0$ , reduced system ${\mathcal{C}}$ has a family of periodic solutions. As a consequence, a $2$ -torus attractor emerges in (A3.5) and, by topological equivalence, in system ${\mathcal{S}}$ . The same occurs when $b = - 1$ , provided that $\mu _{1} \left( {\nu ,\tau } \right) = 0$ and $\mu _{2} \left( {\nu ,\tau } \right) > 0$ . The solution trajectories that evolve on the torus satisfy the TVC if $\left( {\nu ,\tau } \right)$ are chosen to be sufficiently small.

Proof

We have provided enough discussion to prove the main statements in the proposition. To show that the TVC is satisfied, consider the following arguments. Mattana et al. (2009) have shown that, at a steady state, the TVC is satisfied for $0 < u^{*} < 1$. Since $u^{*}$ is smooth with regard to all its arguments, choosing small deviations of the bifurcation parameters from their critical values prevents $u^{*}$ from becoming negative or exceeding $1$.

A discussion regarding the relevance for the economic literature of the range of parameters in Fig. 2 is presented here. Notice first that the Gavrilov-Guckenheimer bifurcation may only occur at the bottom-left region of the parameter space spanned by $\beta$ and $\sigma$. More specifically, bifurcation can only occur for $\beta { \lesssim }0.37$ and $\sigma { \lesssim }0.65$.⁵ Therefore, we can state the following.

Remark 1

The Gavrilov-Guckenheimer bifurcation may only occur when economic agents have a preference against smoothing their consumption over time ( $\sigma$ is low with regard to its domain of existence $)$ and the elasticity of skilled labor to good production is high ( $\beta$ is also low with regard to its domain of existence).

Remark 1 introduces an important departure with regard to the parameter space supporting other types of bifurcations in the Lucas (1988) model and its variants. For example, low values of $\sigma$ and $\beta$ do not seem to be critical ingredients for the onset of a homoclinic bifurcation (cf. Mattana et al. 2009; Bella and Mattana 2014). The same applies for the onset of a Hopf bifurcation (cf. Mattana and Venturi 1999; Nishimura and Shigoka 2009). Low values of $\sigma$ and $\beta$ conversely come back into sight in Bella, Mattana and Venturi (2019), where the focus is on a limit cycle appearing in the full three-dimensional ambient space after the rupture of a homoclinic orbit connecting the unique steady state to itself. Therefore, it seems that for bifurcations giving rise to cycles in well-located planar manifolds,⁶ the parameter space does not need to be restricted to regions where $\sigma$ and $\beta$ are low with regard to their domain of existence. The contrary happens when cycles have effective three-dimensional properties. These elements also provide useful hints regarding the economics at work behind the Gavrilov-Guckenheimer bifurcation in a two-sector continuous time growth model with market imperfections.

3 Implications of toroidal motion for business cycles

An obvious domain where the implications or our results are worth discussing is that of business cycle theory (see Farmer 2014, for a recent appraisal).⁷ As discussed in the Introduction, although the observation of endogenous cycles in a two-sector growth model does not come as a surprise, the specific properties of toroidal motion add a potentially large set of innovative elements to the theory.⁸

Before starting the discussion, consider the following example.⁹

Example 1

Set $\left( {\beta ,\sigma } \right) \in {\mathbf{\Omega }} = \left( {0.3,0.4} \right)$. Additionally, set $A = 1$ and $\delta = 0.05$. It follows that $\mathop \eta = 0.4$ and $\mathop \rho \cong 0.03122$. Then, if $\eta = \mathop \eta$ and $\rho = \mathop \rho$, ${\mathbf{Tr}}\left( {\mathbf{J}} \right) = {\mathbf{Det}}\left( {\mathbf{J}} \right) = 0$ and ${\mathbf{B}}\left( {\mathbf{J}} \right) = - 0.01659$. This economy has $\left( {X^{*} ,u^{*} ,Q^{*} } \right) \cong \left( {4.47794,0.28571,0.11446} \right)$. Since the parameters fall in the green-shaded area of Fig. 1, by Lemma $1$, we know that system ${\mathcal{S}}$ undergoes the Gavrilov-Guckenheimer bifurcation. Now, consider $\left( {\rho ,\eta } \right) = \left( {0.035,0.06} \right)$, which deviates from the critical pair ($\mathop \eta ,\mathop \rho$). Then, ${\mathbf{Tr}}\left( {\mathbf{J}} \right) = 0$, ${\mathbf{Det}}\left( {\mathbf{J}} \right) \cong 0.0012$ and ${\mathbf{B}}\left( {\mathbf{J}} \right) \cong - 0.01496$. By Proposition $1$, this implies that a family of Hopf cycles emerges in the truncated ${\mathcal{C}}$ system. This corresponds to a $2$-torus in the full-dimensional system in (A3.5) and, by topological equivalence, in system ${\mathcal{S}}$.

Figure 3 depicts the solution trajectories, which evolve on a $2$-torus, implied by the coordinates of system ${\mathcal{S}}$.

Fig. 3

The 2-torus represented in the original coordinates $\left( {\user2{X},\user2{u},\user2{Q}} \right)$

In Figs. 4 and 5, we also present the time profiles of the pairs $\left( {K,C} \right)$ and $\left( {H,u} \right)$ over $350$ iterations.¹⁰

Fig. 4

Simulated fluctuations of $\user2{K}$ and $\user2{C}$ around the BGP (350 iterations)

Fig. 5

Simulated fluctuations of $\user2{u}$ and $\user2{~H}$ around the BGP (350 iterations)

It is quite clear that physical capital and consumption are associated with drift dynamics, which give rise to long, slow swings. The opposite occurs for human capital and working time, whose fluctuations are more rapid $($ gyro dynamics$)$. It is also notable that human capital is, by no means, the least volatile variable. Finally, it appears that both $H$ and $u$ alternate between periods of expanding and slowing volatility around the BGP. Interestingly, the period in which the volatility of $u$ increases corresponds to the indents in the physical capital cycle; this phenomenon does not conversely have an evident influence on the smoothness of the cycle for the profile of intertemporal consumption.

How many of these characteristics are matched by real data? Consider the cyclical behavior of $K$ (solid line) and $C$ (dashed line) in Fig. 6 and that of $H$ (dashed line) and $u$ (solid line) in Fig. 7, which were computed with US data over the period 1950–2019.¹¹ To obtain the time series presented hereafter, we removed the quadratic trend present in physical capital, consumption and human capital and then applied an HP filter to all variables to eliminate possible noise effects.

Fig. 6

US real-world fluctuations of $\user2{K}$ and $\user2{C}$ (1950–2019)

Fig. 7

US real-world fluctuations of $\user2{u}$ and $\user2{H}$ (1950–2019)

It is interesting to observe the emerging similarities between the time profiles of the variables of our model (shown in Figs. 4 and 5) and the real-world observations (presented in Figs. 6 and 7). However, some crucial differences are also self-evident. First, in both model series and real-world cases, consumption and physical capital give rise to aperiodic slow swings. There is also a strict correlation between the oscillations of consumption and physical capital and a clear disconnection with the oscillations of $H$ and $u$. A further common characteristic is that human capital seems to be the least volatile variable of the group. Considering these differences, it is clear that real-world human capital and working time do not seem to display the kind of fast oscillations depicted in Fig. 5. It is noteworthy to remember that Penn World Tables (PWT10) only provide a very narrow index of human capital, which is based only on years of schooling and returns with respect to education.

4 A global indeterminacy result in ${\mathbb{R}}^{3}$

Proposition 1 implies that when parameters have the right magnitude, system ${\mathcal{S}}$ has a three-dimensional bounded region enclosed by a two-dimensional surface. Let ${\mathcal{T}}_{{\left( {\nu ,\tau } \right)}} \subset {\mathbb{R}}^{3}$ denote this three-dimensional bounded region that arises for a given $\left( {\nu ,\tau } \right) \subset P$, and let ${\mathcal{T}}_{{\left( {\nu ,\tau } \right){{\%}}}}^{{int}}$ and ${\mathcal{T}}_{{\left( {\nu ,\tau } \right)}}^{{bd}}$ denote the set of all interior points on this three-dimensional bounded region and its boundary (i.e., the 2-torus itself), respectively. It is useful to introduce the following topological definition.

Definition 1

The set of points ${\mathcal{T}}_{{\left( {\nu ,\tau } \right)~}}^{{int}} \cup {\mathcal{T}}_{{\left( {\nu ,\tau } \right)~}}^{{bd}}$ forms the so-called “solid torus”, a connected, compact, and orientable $3$ -dimensional manifold with a boundary.

Now, let

$$ {\mathcal{E}}_{{\left( {\nu ,\tau } \right)}} = \left\{ {\left( {X,u,Q} \right) \in {\mathbb{R}}^{3} :\left( {X,u,Q} \right) \in {\mathcal{T}}_{{\left( {\nu ,\tau } \right)}}^{{int}} } \right\} $$

be a (family of) three-dimensional manifold containing the set of all possible paths starting on ${\mathcal{T}}_{{\left( {\nu ,\tau } \right)}}^{{int}}$. Then, by Definition 1, all paths starting on the (compact) set ${\mathcal{E}}_{{\left( {\nu ,\tau } \right)}}$ are bound to stay in ${\mathcal{E}}_{{\left( {\nu ,\tau } \right)}}$.

Finally, let

$$ {{\Psi }} = \left\{ {\left( {A,\beta ,\delta ,\eta ,\rho ,\sigma } \right):\left({\dot{X},\dot{u},\dot{Q}} \right) = {\mathcal{S}}{\text{~implies~~a~}}2{\text{ - Torus}}} \right\} $$

be the set of all parametric values that give rise to a 2-torus.

We are now ready to prove the main proposition of the paper.

Proposition 2 (Global indeterminacy of the equilibrium)

Recall Proposition $1$ and let $\left( {A,\beta ,\delta ,\eta ,\rho ,\sigma } \right) \in \Psi$ . Assume $\left( {\nu ,\tau } \right) \subset P$ . Consider the case in which the initial condition of the aggregate capital stock is chosen such that $X\left( 0 \right) \in {\mathcal{E}}_{{\left( {\nu ,\tau } \right)}}$ . Then, system ${\mathcal{S}}$ exhibits global indeterminacy of the equilibrium.

Proof

In Proposition 1, we have shown that there exist regions in the parameter space such that a 2-torus emerges in the full $\left( {X,u,Q} \right)$ space. We have also discussed that this implies the existence of a compact, $3$-dimensional manifold with a boundary. Consider the case in which we choose to initialize the dynamics at $X\left( 0 \right) \in {\mathcal{E}}_{{\left( {\nu ,\tau } \right)}} = \left\{ {\left( {X,u,Q} \right) \in {\mathbb{R}}^{3} :\left( {X,u,Q} \right) \in {\mathcal{T}}_{{\left( {\nu ,\tau } \right)}}^{{int}} } \right\}$. Then, there exists a continuum of initial values $\left( {u\left( 0 \right),Q\left( 0 \right)} \right)$ of the jump variables, which give rise to recurrent orbits, namely, perfect-foresight equilibrium trajectories that are bound to stay in the three-dimensional manifold ${\mathcal{E}}_{{\left( {\nu ,\tau } \right)}}$ at all times.

Remark 2

The ability of a three-dimensional object composed of ${\mathcal{T}}_{{\left( {\nu ,\tau } \right)~}}^{{int}} \cup {\mathcal{T}}_{{\left( {\nu ,\tau } \right)~}}^{{bd}}$ to attract orbits wandering in the outside domain depends on the stability properties of the Hopf cycle that arises in system ${\mathcal{C}}$ . For the sake of a simple discussion, we do not formally investigate the stability properties of a $2$ -torus. ¹²

As is clear, Proposition 2 establishes a global indeterminacy property for a two-sector continuous time growth model for reasons different than the mere appearance of closed orbits in well-behaved planar submanifolds. At this point, an interesting but subtle detail involves the empirical relevance of this type of global indeterminacy. To highlight this issue, we construct in Fig. 8 the basin of attraction of the torus implied by the parameter values in Example 1 projected onto the $\left( {X,u} \right)$ plane.¹³ A note of caution is due here. Since there is no unambiguous definition and measurement of human capital investment, this type of perspective can only be considered in relative terms.

Fig. 8

The basin of attraction of a 2-torus

The continuous line represents the boundary of the gray-shaded area containing the set of initial conditions that converge to the torus.¹⁴ The darker-colored object inside the shaded area, containing the coordinates of the steady state $\left( {X^{*} ,u^{*} } \right) \cong \left( {4.47794,0.28571} \right)$, is the projection of the torus. Note that since the gray-shaded area extends outside the projection of the torus, we observe a stable object, meaning that the 2-torus is attractive. The position of the gray-shaded area in the phase space provides interesting hints on the characteristics of real economies in an indeterminacy trap. At least for the parameter configuration in Example 1, it seems that given the parameter configuration, only countries with a fraction of the total time devoted to work $\left( {u{ \lesssim }0.45} \right)$ are prone to multiplicity; the phenomenon occurs for a wide range of $X = K/H$ ratio values that are above a threshold level of approximately $X \simeq 4.15$; however, for values of $X{ \gtrsim }4.7$ the area shrinks considerably.

To put these findings in perspective, we also run a similar procedure for the case in which system ${\mathcal{S}}$ presents global indeterminacy due to the emergence of a homoclinic bifurcation. As discussed in the Introduction, the onset of the homoclinic bifurcation is studied in Mattana et al. (2009) and in Bella and Mattana (2014) by the application of the Kopell-Howard and Bogdanov-Takens theorems, respectively. To keep the analysis qualitatively similar, we numerically study the basin of attraction of the homoclinic loop by using the parameter same values as those in Example 3 in Mattana et al. (2009) in the transformed $\left( {X,u} \right)$ planar space given by the center manifold. The two steady states of system ${\mathcal{S}}$ coalesce at a unique fixed point at the critical values of the bifurcation parameters $\left( {\tilde{\rho },\tilde{\sigma }} \right) \cong \left( {0.03852,0.24696} \right)$. Note that in this parametric configuration, $\sigma$ and $A$ are small with regard to their domain of existence, whereas $\beta$ is quite large. Now, setting $\left( {\rho ,\sigma } \right) = \left( {0.0385,0.244} \right)$, which is slightly away from the critical values, system ${\mathcal{S}}$ presents two steady states; one of which is a saddle, while the other is a nonsaddle. This scenario is depicted in Fig. 9, where the limit homoclinic orbit and the Hopf cycle are reported.

Fig. 9

The basin of attraction of a Hopf cycle

The light gray-shaded area is the basin of attraction of the closed Hopf loop, whereas the dark gray-shaded area is the basin of attraction of the homoclinic limit cycle in the interior of the bundle of unstable orbits. Moreover, the dots reported represent the two steady states. Finally, note that we have chosen a basis of eigenvectors so that the oscillations of the $u$ variable roughly match those observed in real-world data of the US economy (cf. Figure 5). The union of the two gray-shaded basins of attraction is therefore the locus of initial conditions, which give rise to recurrent orbits for this specific choice of parameters. Again, countries prone to indeterminacy problems are those with a relatively low fraction of total time devoted to work $\left( {u{ \lesssim }0.32} \right)$; furthermore, the phenomenon occurs for a specific range of $X = K/H$ ratio values, which is given by $0.15{ \lesssim }X{ \lesssim }0.4$.

Although it is not possible to directly compare the regions in the phase space at which economies find themselves in the indeterminacy trap, some interesting similarities are self-evident. First, the range of $u$ for which there is indeterminacy roughly overlaps in the two examples. Furthermore, there is a region to the left of both basins of attraction in Figs. 8 and 9 for which no indeterminacy problems arise. It seems that countries with physical-to-human capital ratios at the left of a given threshold are more likely to face unique growth paths. The comparison when we consider the indeterminacy implications when the physical-to-human capital ratio tends to become higher is not as straightforward. In this case, when we look at indeterminacy problems arising because of the onset of a Hopf cycle, countries at the right of the basin of attraction are again free of indeterminacy considerations; in the case when indeterminacy occurs due to the existence of a 2-torus, the indeterminacy area extends further right.

Gathering these traits of economies that are prone to indeterminacy problems helps to determine the mechanics behind the onset of this phenomenon. Assume an economy needs to increase an otherwise low $X = K/H$ ratio. As discussed in Mulligan and Sala-i-Martin (1993), a low value of $X$ can be raised either through an increase in savings (a decrease in $Q$, which is the consumption to capital ratio) or through an increase in time allocation, $u$, to the production of physical output. Since $\sigma$ is low, agents have a preference against smoothing their consumption; hence, they prefer the second option (wage rate effect). Now, as long as the system is sufficiently far from the steady state (strong unbalance) or the externality factor (which, as discussed in Krugman (1991) and Matsuyama (1991), also measures the degree of complementarity across agents’ decisions), is sufficiently low, there is a unique perfect-foresight choice of $u$ that is able to lead the economy to a long-run equilibrium. There is, however, a threshold level of $X$ beyond which the complementarity effects dominate, and there appear multiple (inefficient) choices of $u$, each of which converges to the steady state. This phenomenon can typically be categorized as a market coordination failure when individual decisions are called strategic complements. A complementary discussion applies when $X$ is high.

Does this theoretically derived information have a counterpart in real world data? In other words, are there signs that an economy might evolve, given some regions in the domain of the state (fundamental) variable, under the influence of indeterminacy? A perfect-foresight equilibrium is determinate if it is locally unique; it is indeterminate if many other equilibria are arbitrarily close to the first equilibrium, depending on specific choices of the control variables. As widely recognized in the specialized literature (cf., inter al. Lubik and Schorfheide 2004; Beyer and Farmer 2007), shocks to these nonfundamental variables may then contribute to the total variance of an economy.

Therefore, to find evidence of indeterminacy of the equilibrium in real economies, we again draw information from Penn World Tables (PW10) and propose the following econometric exercises. The PW10 dataset provides an index of human capital accumulation based on years of schooling and returns with respect to education (not monetary values) and the physical capital stock at constant 2017 national prices (in millions of 2017 USD) for 183 countries. The real GDP at constant 2017 national prices (in millions of 2017 USD) is also available in the dataset. The period covered includes 1950–2019; however, since there appears to be clear overall erratic variable behavior for the period 1950–1955, we drop the pre-1955 observations in the construction of our sample. Then, for each country in the panel, we compute the following: (i) the ratio between physical capital (in monetary terms) and the index of human capital (variable $I_{{K/H}}$ below),¹⁵ and (ii) a proxy for working time ($P_{u} )$ (cf. footnote 14, for the computation). Following the cited literature, we also need a proxy of the volatility of (future) growth rates. Thus, we compute the five-year forward variance $\left( {\sigma _{{g\_y}}^{5} } \right)$ from time $t_{i}$, where $i = 5$. If indeterminacy occurs only for specific ranges of the domain of the $X$ variable, we expect to observe non-linearity in the empirical relationship between $X$ and the variance of (future) growth rates.

There are several possibilities for empirically modeling indeterminacy-induced non-linearity in the forward variance of the growth rate. We first consider the following quadratic specification

$$ \sigma _{{g_{{\_y}} ,t_{{i - 5}} }}^{5} = \alpha _{0} + \alpha _{1} i_{{K/H}} + \alpha _{2} i_{{K/H}}^{2} + \sum {\mathbf{\varphi Z}} + e~ $$

(4)

where we have dropped the time and unit foots for a simpler appraisal. In Eq. (4), $i_{{K/H}}$ is the logarithm of $I_{{K/H}}$ and ${\mathbf{Z}}$ is a vector of control variables. $\alpha _{0}$ is the (unit-specific) constant. $\alpha _{1}$ and $\alpha _{2}$ are the parameters of interest, whereas ${\mathbf{\varphi }}$ is a vector of other parameters. $e$ is a stochastic error term. The (joint) null hypothesis is $H_{0} :\alpha _{1} = \alpha _{2} = 0$, which implies that time $t$ observations of the $i_{{K/H}}$ (state) variable carry no explanatory power over the growth rate forward volatility. Conversely, the rejection of the null hypothesis leads to a world where uncertainty on growth prospects (proxied by $\sigma _{{g\_y}}^{5} )$ depends on the initial position of an economy in terms of the state variable (indeterminacy)). Note that we are specifically interested in the case where $\alpha _{2} < 0$, which implies the existence of a localized $I_{{K/H}}$ ratio at which the internal mechanisms of the economy maximize uncertainty with regard to future growth. The panel is strongly unbalanced, and the resulting observations are 5.561. As a control vector, we use time $t$ growth rates of the GDP, population and hours worked.

A robust least-squares estimation of Eq. (4) leads to strongly significant values of $\alpha _{1} > 0$ and $\alpha _{2} < 0$ (p values lower than $0.01$); the additional variance explained by the variables $i_{{K/H}}$ and $i_{{K/H}}^{2}$ is, however, quite small (only slightly above 1%).¹⁶ Significant coefficients are also obtained when we re-estimate Eq. (4) with $P_{u}$ and $P_{u}^{2}$ instead of $i_{{K/H}}$ and $i_{{K/H}}^{2}$. With the estimated parameters, it is easy to verify that the values of $I_{{K/H}}$ at which the growth uncertainty is maximized fall within the 7^th decile of the frequency distribution of $I_{{K/H}}$. For the $P_{u}$ variable, the same happens at the 5^th decile. These findings provide remarkable support for the pervasiveness of global indeterminacy in the growth process.

Confirmation of these results is also obtained by estimating the following (discrete) model with interaction terms:

$$ \sigma _{{g_{{\_y}} ,t_{{i - 5}} }}^{5} = \beta _{0} + \beta _{1} i_{{K/H}} + \mathop \sum \limits_{1}^{9} \beta _{{2i}} d_{i} i_{{K/H}} + \sum {\mathbf{\varphi Z}} + e~ $$

(5)

where the $d_{i}$ dummies are coded as 0 or 1 according to whether the observations of $i_{{K/H}}$ are in a particular decile of their frequency distribution. The robust least squares estimation method provides evidence of significant $\beta _{{2i}}$ coefficients for observations lying in the 7th, 8th and 9th deciles of the frequency distribution of $i_{{K/H}}$ (p values lower than $0.05$). By repeating the estimation of Eq. (5) with the $P_{u}$ variable instead, the results are less straightforward and denote weak significance (p values lower than $0.1$) for the coefficients associated with the 2nd, 3rd and 4th deciles of the frequency distribution of $P_{u}$.

5 Conclusions

This paper shows that the Lucas (1988) growth model with flow externalities à la Chamley (1993) may give rise to global indeterminacy of the equilibrium for reasons different than the existence of closed loops in well-located two-dimensional manifolds. The mathematical underpinning of our results lies in the properties of a pitchfork–Hopf interaction, which is a simple linear degeneracy that can be associated with the onset of a 2-torus trapping region in ${\mathbb{R}}^{3}$ enclosed by a two-dimensional surface. By studying the basin of attraction of the torus, given a specific parametrization, we show that economies enmeshed in indeterminacy problems typically have a low fraction of the total time devoted to work. Interestingly, the same happens in the case of indeterminacy caused by the birth of closed loops in well-located planar submanifolds. Using data from PWT10, we find an empirical counterpart to these results. Indeed, only in specific ranges in their frequency of distribution do our proxies of time devoted to work and physical-to-human capital ratios have predictive power with regard to future growth uncertainty. Our results also have significance for the field of endogenous business cycles. In particular, without introducing any stochasticity, a solution orbit wrapping around a torus is naturally aperiodic; namely, it never comes back exactly onto itself. This is an important difference with regard to other cases of cyclical solutions (such as Hopf cycles) where unperturbed orbits are bound to intersect the same coordinates in the phase space forever. Furthermore, we have shown that physical capital and consumption are highly correlated and give rise to long, slow swings, very similar to those observed in real data of the US economy (data from the PWT10). The opposite occurs for human capital and working time, whose fluctuations appear more rapid. In our examples, human capital appears to be the least volatile variable of the model.

Appendix

1.1 Derivation of system ${\mathbf{\mathcal{S}}}$, its Jacobian and characteristic polynomial

Recall problem ${\mathcal{P}}$ in Sect. 1.

$$ \begin{gathered} ~\mathop {\max }\limits_{{C\left( t \right),~u\left( t \right)}} \mathop \smallint \limits_{0}^{\infty } \frac{{C^{{1 - \sigma }} - 1}}{{1 - \sigma }}e^{{ - \rho t}} dt \hfill \\ \,\,\,\,\,\,\,\,\,{\text{subject to}} \hfill \\ \end{gathered} $$

$$ \dot{K} = AK^{\beta } \left( {uH} \right)^{{1 - \beta }} H_{a}^{\gamma } - C $$

$$ \dot{H} = \delta H\left( {1 - u} \right)\left( {1 - u} \right)_{a}^{\eta } $$

$$ \begin{gathered} K\left( 0 \right) = K_{0} \hfill \\ H\left( 0 \right) = H_{0} \hfill \\ \end{gathered} $$

The associated current-value Hamiltonian is

$$ {\mathcal{H}} = \frac{{C^{{1 - \sigma }} - 1}}{{1 - \sigma }} + \lambda \left[ {AK^{\beta } (uH)^{{1 - \beta }} H_{a}^{\gamma } - C} \right] + \mu \left[ {\delta H\left( {1 - u} \right)(1 - u)_{a}^{\eta } } \right] $$

First-order necessary conditions imply

$$ \begin{aligned} \partial {\mathcal{H}}/\partial C & = C^{{ - \sigma }} - \lambda = 0 \\ \partial {\mathcal{H}}/\partial u & = \lambda \left( {1 - \beta } \right)AK^{\beta } u^{{ - \beta }} H^{{1 - \beta + \gamma }} - \mu \delta H(1 - u)^{\eta } = 0 \\ \dot{\lambda } & = \lambda \rho - \frac{{\partial {\mathcal{H}}}}{{\partial K}} = \lambda \rho - \lambda \beta AK^{{\beta - 1}} (uH)^{{1 - \beta + \gamma }} \\ \dot{\mu } & = \mu \rho - \frac{{\partial {\mathcal{H}}}}{{\partial H}} = \mu \rho - \lambda \left( {1 - \beta } \right)AK^{\beta } u^{{1 - \beta }} H^{{ - \beta + \gamma }} - \mu \delta (1 - u)^{{1 + \eta }} \\ \end{aligned} $$

Taking logs, differentiating with respect to time, we obtain

$$ \begin{aligned} \dot{C}/C & = - \frac{1}{\sigma }\left[ {\rho - \beta AK^{{\beta - 1}} (uH)^{{1 - \beta + \gamma }} } \right] \\ \dot{K}/K & = AK^{{\beta - 1}} u^{{1 - \beta }} H^{{1 - \beta + \gamma }} - \frac{C}{K} \\ \dot{H}/H & = \delta (1 - u)^{{1 + \eta }} \\ \dot{u} & = \frac{{u\left( {1 - u} \right)}}{{\beta \left( {1 - u} \right) - \eta u}}\,\,\left[ {\delta (1 - u)^{\eta } - \delta \left( {\beta - \gamma } \right)(1 - u)^{{1 + \eta }} - \beta \frac{C}{K}} \right] \\ \end{aligned} $$

Consider now the stationarizing variables, $X \equiv KH^{{\frac{{1 - \beta + \gamma }}{{\beta - 1}}}}$ and $Q \equiv \frac{C}{K}$, we end up with system ${\mathcal{S}}$ in the text

$$ \begin{aligned} \dot{X} & = AX^{\beta } u^{{1 - \beta }} + \frac{{\delta \left( {1 - \beta + \gamma } \right)}}{{\beta - 1}}X(1 - u)^{{1 + \eta }} - QX \\ \dot{u} & = \frac{{u\left( {1 - u} \right)}}{{\beta \left( {1 - u} \right) - \eta u}}\left[ {\delta (1 - u)^{\eta } - \delta \left( {\beta - \gamma } \right)(1 - u)^{{1 + \eta }} - \beta Q} \right] \\ \dot{Q} & = - \frac{\rho }{\sigma }Q + \frac{{\beta - \sigma }}{\sigma }AX^{{\beta - 1}} u^{{1 - \beta }} Q + Q^{2} \\ \end{aligned} $$

Notice that, since ${\mathcal{H}}$ is concave in the state variables, Focs are also sufficient to solve problem ${\mathcal{P}}$.

Let $\left( {X^{*} ,u^{*} ,Q^{*} } \right)$ be values of $\left( {X,u,Q} \right)$ such that $\dot{X} = \dot{u} = \dot{Q} = 0$. Solving system ${\mathcal{S}}$ in the text we obtain that, along a BGP

$$ \begin{aligned} 0 & = \left( {1 - u^{*} } \right)^{\eta } - \sigma \left( {1 - u^{*} } \right)^{{1 + \eta }} - \frac{\rho }{\delta } \\ X^{*} & = u^{*} \left( {\beta A} \right)^{{\frac{1}{{1 - \beta }}}} \left[ {\delta \left( {1 - u^{*} } \right)^{\eta } } \right]^{{\frac{1}{{\beta - 1}}}} \\ Q^{*} & = \frac{{\delta \left( {1 - u^{*} } \right)^{\eta } }}{\beta }\left[ {1 - \beta \left( {1 - u^{*} } \right)} \right] \\ \end{aligned} $$

Let ${\mathbf{J}}$ be the Jacobian of the right-hand side of system ${\mathcal{S}}$ evaluated along the BGP. The single elements of ${\mathbf{J}}$ are as follows

$$ j_{{11}}^{*} = \frac{{\partial \dot{X}}}{{\partial X}} = \beta AX^{{*\beta - 1}} u^{{*1 - \beta }} - Q^{*} - \delta \left( {1 - u^{*} } \right)^{{1 + \eta }} = \frac{{\dot{X}}}{X} - \left( {1 - \beta } \right)AX^{{*\beta - 1}} u^{{*1 - \beta }} $$

$$ j_{{12}}^{*} = \frac{{\partial \dot{X}}}{{\partial u}}{\text{ = ~}}\left( {1 - \beta } \right)AX^{{*\beta }} u^{{* - \beta }} + \delta \left( {1 + \eta } \right)X^{*} \left( {1 - u^{*} } \right)^{\eta } $$

$$ j_{{13}}^{*} = \frac{{\partial \dot{X}}}{{\partial Q}} = - X^{*} $$

$$ j_{{21}}^{*} = \frac{{\partial \dot{u}}}{{\partial X}} = 0 $$

$$ \begin{aligned} j_{{22}}^{*} & = \frac{{\partial \dot{u}}}{{\partial u}} = \frac{{\left( {1 - 2u^{*} } \right)\left( {\beta \left( {1 - u^{*} } \right) - \eta u^{*} } \right) - u^{*} \left( {1 - u^{*} } \right)\left( { - \beta - \eta } \right)}}{{\left[ {\beta \left( {1 - u^{*} } \right) - \eta u^{*} } \right]^{2} }}\left[ {\delta \left( {1 - u^{*} } \right)^{\eta } - \delta \beta \left( {1 - u^{*} } \right)^{{1 + \eta }} - \beta Q^{*} } \right] \\ & \,\,\,\,\, + \frac{{\delta u^{*} \left( {1 - u^{*} } \right)^{\eta } }}{{\beta \left( {1 - u^{*} } \right) - \eta u^{*} }}\left[ {\beta \left( {1 + \eta } \right)\left( {1 - u^{*} } \right) - \eta } \right] \\ \end{aligned} $$

$$ j_{{23}}^{*} = \frac{{\partial \dot{u}}}{{\partial Q}} = - \frac{{\beta u^{*} \left( {1 - u^{*} } \right)}}{{\beta \left( {1 - u^{*} } \right) - \eta u^{*} }} $$

$$ j_{{31}}^{*} = \frac{{\partial \dot{Q}}}{{\partial X}} = - \left( {1 - \beta } \right)\frac{{\beta - \sigma }}{\sigma }AX^{{*\beta - 2}} u^{{*1 - \beta }} Q^{*} $$

$$ j_{{32}}^{*} = \frac{{\partial \dot{Q}}}{{\partial u}} = \frac{{\beta - \sigma }}{\sigma }\left( {1 - \beta } \right)AX^{{*\beta - 1}} u^{{* - \beta }} Q^{*} $$

$$ j_{{33}}^{*} = \frac{{\partial \dot{Q}}}{{\partial Q}} = \frac{{\beta - \sigma }}{\sigma }AX^{{*\beta - 1}} u^{{*1 - \beta }} - \frac{\rho }{\delta } + 2Q^{*} $$

Therefore, along the BGP, using $Z^{*} = AX^{{*\beta - 1}} u^{{*1 - \beta }}$, we have

$$ {\mathbf{J}} = \left[ {\begin{array}{*{20}c} { - \left( {1 - \beta } \right)Z^{*} } & {{\text{~}}j_{{12}}^{*} } & { - X^{*} } \\ 0 & {j_{{22}}^{*} } & { - \frac{{\beta u^{*} \left( {1 - u^{*} } \right)}}{{\beta \left( {1 - u^{*} } \right) - \eta u^{*} }}} \\ { - \frac{{\beta - \sigma }}{\sigma }\frac{{Z^{*} }}{{X^{*} }}Q^{*} } & {\frac{{\beta - \sigma }}{\sigma }\frac{{Z^{*} }}{{u^{*} }}Q^{*} } & {Q^{*} } \\ \end{array} } \right] $$

(A1.1)

where $j_{{12}}^{*} = \frac{{\left( {1 - \beta } \right)Z^{*} X^{*} }}{{u^{*} }} + \delta \left( {1 + \eta } \right)X^{*} (1 - u^{*} )^{\eta }$ and $j_{{22}}^{*} = \frac{{\delta u^{*} (1 - u^{*} )^{\eta } }}{{\beta \left( {1 - u^{*} } \right) - \eta u^{*} }}\left[ {\beta \left( {1 + \eta } \right)\left( {1 - u^{*} } \right) - \eta } \right]$.

The eigenvalues of ${\mathbf{J}}$ are the solutions of the characteristic equation

$$ {\text{det}}\left( {\lambda {\mathbf{I}} - {\mathbf{J}}} \right) = \lambda ^{3} - {\text{Tr}}\left( {\mathbf{J}} \right)\lambda ^{2} + {\text{B}}\left( {\mathbf{J}} \right)\lambda - {\text{Det}}\left( {\mathbf{J}} \right) $$

where I is the identity matrix and Tr $\left( {\mathbf{J}} \right)$, Det $\left( {\mathbf{J}} \right)$ and B $\left( {\mathbf{J}} \right)$ are Trace, Determinant and Sum of Principal Minors of ${\mathbf{J}}$, respectively. We obtain

$$ {\mathbf{Tr}}\left( {\mathbf{J}} \right) = \frac{{\delta u^{*} (1 - u^{*} )^{\eta } }}{{\beta \left( {1 - u^{*} } \right) - \eta u^{*} }}\left[ {\beta \left( {2 + \eta } \right)\left( {1 - u^{*} } \right) - \eta \left( {1 + u^{*} } \right)} \right] $$

(A1.2)

$$ {\mathbf{Det}}\left( {\mathbf{J}} \right) = - \frac{{\beta Z^{*} Q^{*} }}{{\beta - \left( {\beta + \eta } \right)u^{*} }}\frac{{\delta u^{*} (1 - u^{*} )^{\eta } }}{\sigma }\left[ {\sigma \left( {1 + \eta } \right)\left( {1 - u^{*} } \right) - \eta } \right] $$

(A1.3)

$$ {\mathbf{B}}\left( {\mathbf{J}} \right) = j_{{22}}^{*} \left[ {Q^{*} - \left( {1 - \beta } \right)Z^{*} } \right] + Z^{*} Q^{*} \frac{{\beta \left( {1 - \sigma } \right)\eta u^{*} - \beta \sigma \left( {1 - \beta } \right)\left( {1 - u^{*} } \right)}}{{\sigma \left[ {\beta \left( {1 - u^{*} } \right) - \eta u^{*} } \right]}} $$

(A1.4)

1.2 Existence of the Hypernormal form

1.2.1 Translation to the origin

Substitute $\mathop X \equiv X - \hat{X}^{*}$, $\mathop u \equiv u - \hat{u}^{*}$, $\mathop Q = Q - \hat{Q}^{*}$, $v = \rho - \hat{\rho },\mu = \sigma - \hat{\sigma },\tau = \eta - \hat{\eta }$ in system ${\mathcal{S}}$. We obtain

$$ \begin{gathered} \mathop {\tilde{X}}\limits^{ \cdot } = A\left( {\hat{X}^{*} + \mathop X } \right)^{\beta } \left( {\hat{u}^{*} + \mathop u } \right)^{{1 - \beta }} - \delta \left( {\hat{X}^{*} + \mathop X } \right)(1 - \hat{u}^{*} - \mathop u )^{{1 + \hat{\eta }{\text{~}} + \tau }} - \left( {\hat{Q}^{*} + \mathop Q } \right)\left( {\hat{X}^{*} + \mathop X } \right) \hfill \\ \mathop {\tilde{u}} = \frac{{\left( {\hat{u}^{*} + \mathop u } \right)\left( {1 - \hat{u}^{*} - \mathop u } \right)\left[ {\delta (1 - \hat{u}^{*} - \mathop u )^{{\hat{\eta }{\text{~}} + \tau }} - \beta \delta (1 - \hat{u}^{*} - \mathop u )^{{1 + \hat{\eta }{\text{~}} + \tau }} - \beta \left( {\hat{Q}^{*} + \mathop Q } \right)} \right]}}{{\beta \left( {1 - \hat{u}^{*} - \mathop u } \right) - \left( {\hat{\eta }{\text{~}} + \tau } \right)\left( {\hat{u}^{*} + \mathop u } \right)}} \hfill \\ \mathop {\tilde{Q}}\limits^{ \cdot } = - \frac{{\hat{\rho } + v}}{{\hat{\sigma } + \mu }}\left( {\hat{Q}^{*} + \mathop Q } \right) + \left( {\frac{\beta }{{\hat{\sigma } + \mu }} - 1} \right)A\left( {\hat{X}^{*} + \mathop X } \right)^{{\beta - 1}} \left( {\hat{u}^{*} + \mathop u } \right)^{{1 - \beta }} \left( {\hat{Q}^{*} + \mathop Q } \right) + \left( {\hat{Q}^{*} + \mathop Q } \right)^{2} \hfill \\ \end{gathered} $$

(A2.1)

1.2.2 Second-order Taylor expansion

Consider the second-order expansion of system in (A2.1) with respect to the variables $\left( {\mathop X ,\mathop u ,\mathop Q ,v,\mu ,\tau } \right)$

$$ \left( {\begin{array}{*{20}c} {\mathop {\tilde{X}} } \\ {\mathop {\tilde{u}} } \\ {\mathop {\tilde{Q}} } \\ \end{array} } \right) = {\mathbf{J}}\left( {\begin{array}{*{20}c} {\mathop X } \\ {\mathop u } \\ {\mathop Q } \\ \end{array} } \right) + {\mathbf{M}}\left( {\mu ,\nu ,\tau } \right)\left( {\begin{array}{*{20}c} {\mathop X } \\ {\mathop u } \\ {\mathop Q } \\ \end{array} } \right) + \left( {\begin{array}{*{20}c} {\mathop f _{1} \left( {\mathop X ,\mathop u ,\mathop Q } \right)} \\ {\mathop f _{2} \left( {\mathop X ,\mathop u ,\mathop Q } \right)} \\ {\mathop f _{3} \left( {\mathop X ,\mathop u ,\mathop Q } \right)} \\ \end{array} } \right) $$

where ${\mathbf{J}}$ is as in (A1.1) and ${\mathbf{M}}\left( {\mu ,\nu ,\tau } \right)$ is as follows

$$ {\mathbf{M}}\left( {\mu ,\nu ,\tau } \right) = \left[ {\begin{array}{*{20}c} { - \delta (1 - \hat{u}^{*} )^{{1 + \eta }} N\tau } & {\delta \hat{X}^{*} (1 - \hat{u}^{*} )^{\eta } \left( {1 + N^{{1 + \eta }} } \right)\tau } & 0 \\ 0 & { - \delta \hat{u}^{*} (1 - \hat{u}^{*} )^{\eta } N\tau } & { - \frac{{\beta \hat{u}^{{*2}} \left( {1 - \hat{u}^{*} } \right)\tau }}{{\left( {\beta \left( {1 - \hat{u}^{*} } \right) - \eta \hat{u}^{*} } \right)^{2} }}} \\ {\frac{{\beta \left( {1 - \beta } \right)Z^{*} }}{{\sigma ^{2} \hat{X}^{{*2}} }}\mathop Q \mu } & { - \frac{{\beta \left( {1 - \beta } \right)Z^{*} }}{{\sigma ^{2} \hat{u}^{*} }}\mathop Q \mu } & {\frac{{\rho - \beta Z^{*} }}{{\sigma ^{2} }}\mu - \frac{\rho }{{\sigma + \mu }}v} \\ \end{array} } \right] $$

where $N = {\text{ln}}\left( {1 - \hat{u}^{*} } \right)$ and

$$ \mathop f _{1} \left( {\mathop X ,\mathop u ,\mathop Q } \right) = A_{1} \mathop X ^{2} + A_{2} \mathop u ^{2} + A_{3} \mathop X \mathop u - \mathop X \mathop Q $$

$$ \mathop f _{2} \left( {\mathop X ,\mathop u ,\mathop Q } \right) = B_{1} \mathop u ^{2} + B_{2} \mathop u \mathop Q $$

$$ \mathop f _{3} \left( {\mathop X ,\mathop u ,\mathop Q } \right) = C_{1} \mathop X ^{2} + C_{2} \mathop u ^{2} + \mathop Q ^{2} + C_{4} \mathop X \mathop Q + C_{5} \mathop u \mathop Q + C_{6} \mathop X \mathop u $$

All relevant coefficients have been calculated but not reported for the sake of a simple representation.

1.2.3 (The organizing center) Calculation of the Eigenvectors in correspondence of $\lambda _{1} = \lambda _{2} = \lambda _{3} = 0$.

We need to solve

$$ \begin{array}{*{20}c} {{\mathbf{Ju}} = 0} \\ {{\mathbf{Jv}} = {\mathbf{u}}} \\ {{\mathbf{Jz}} = {\mathbf{v}}} \\ \end{array} $$

Possible candidates are the following

$$ {\mathbf{u}} = \left[ {\begin{array}{*{20}c} {\frac{{j_{{12}}^{*} j_{{23}}^{*} - j_{{22}}^{*} j_{{13}}^{*} }}{{j_{{11}}^{*} j_{{22}}^{*} }}} \\ { - \frac{{j_{{23}}^{*} }}{{j_{{22}}^{*} }}} \\ 1 \\ \end{array} } \right]~;~{\mathbf{v}} = \left[ {\begin{array}{*{20}c} {\frac{{(j_{{22}}^{*} )^{2} + j_{{23}}^{*} j_{{32}}^{*} }}{{j_{{31}}^{*} (j_{{22}}^{*} )^{2} }}} \\ { - \frac{{j_{{23}}^{*} }}{{(j_{{22}}^{*} )^{2} }}} \\ 0 \\ \end{array} } \right]~;~{\mathbf{z}} = \left[ {\begin{array}{*{20}c} {\frac{{(j_{{22}}^{*} )^{2} + j_{{23}}^{*} j_{{32}}^{*} + j_{{13}}^{*} j_{{31}}^{*} }}{{j_{{11}}^{*} j_{{31}}^{*} (j_{{22}}^{*} )^{2} }}} \\ 0 \\ { - \frac{1}{{(j_{{22}}^{*} )^{2} }}} \\ \end{array} } \right]~ $$

We can therefore form the following transformation matrix

$$ {\mathbf{T}} = \left[ {\begin{array}{*{20}c} {u_{1} } & {v_{1} } & {z_{1} } \\ {u_{2} } & {v_{2} } & 0 \\ 1 & 0 & {z_{3} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\frac{{j_{{12}}^{*} j_{{23}}^{*} - j_{{22}}^{*} j_{{13}}^{*} }}{{j_{{11}}^{*} j_{{22}}^{*} }}} & {\frac{{(j_{{22}}^{*} )^{2} + j_{{23}}^{*} j_{{32}}^{*} }}{{j_{{31}}^{*} (j_{{22}}^{*} )^{2} }}} & {\frac{{(j_{{22}}^{*} )^{2} + j_{{23}}^{*} j_{{32}}^{*} + j_{{13}}^{*} j_{{31}}^{*} }}{{j_{{11}}^{*} j_{{31}}^{*} (j_{{22}}^{*} )^{2} }}} \\ { - \frac{{j_{{23}}^{*} }}{{j_{{22}}^{*} }}} & { - \frac{{j_{{23}}^{*} }}{{(j_{{22}}^{*} )^{2} }}} & 0 \\ 1 & 0 & { - \frac{1}{{(j_{{22}}^{*} )^{2} }}} \\ \end{array} } \right] $$

1.2.4 Hypernormal form calculation (Gamero et al. 1991)

Using the change of coordinates

$$ \left( {\begin{array}{*{20}c} {\mathop X } \\ {\mathop u } \\ {\mathop Q } \\ \end{array} } \right) = {\mathbf{T}}\left( {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {w_{3} } \\ \end{array} } \right) $$

we can proceed in the computation of the normal form. Considering that

$$ {\mathbf{T}}^{{ - 1}} = 1/D\left[ {\begin{array}{*{20}c} { - v_{2} z_{3} } & {v_{1} z_{3} } & {v_{2} z_{1} } \\ {u_{2} z_{3} } & {z_{1} - u_{1} z_{3} } & { - u_{2} z_{1} } \\ {v_{2} } & { - v_{1} } & { - \mathop u _{1} v_{2} + \mathop u _{2} v_{1} } \\ \end{array} } \right] $$

where $D = v_{2} z_{1} - u_{1} v_{2} z_{3} + u_{2} v_{1} z_{3}$, we are able to put system (A2.1) in the following third-order normal form

$$ \left( {\begin{array}{*{20}c} {\mathop {w_{1} } } \\ {\mathop {w_{2} } } \\ {\mathop {w_{3} } } \\ \end{array} } \right) = {\mathbf{L}}\left( {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {w_{3} } \\ \end{array} } \right) + {\mathbf{M}}\left( {\mu ,\nu ,\tau } \right)\left( {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {w_{3} } \\ \end{array} } \right) + \left( {\begin{array}{*{20}c} {\bar{F}_{1} \left( {w_{1} ,w_{2} ,w_{3} } \right)} \\ {\bar{F}_{2} \left( {w_{1} ,w_{2} ,w_{3} } \right)} \\ {\bar{F}_{3} \left( {w_{1} ,w_{2} ,w_{3} } \right)} \\ \end{array} } \right) + {\mathcal{O}}\left( {w_{1} ,w_{2} ,w_{3} ^{4} } \right)~ $$

(A2.2)

where

$$ {\mathbf{L}} = \left[ {\begin{array}{*{20}c} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ \end{array} } \right]~ $$

(A2.3)

and $\bar{F}_{i} ,$ $i = 1,2,3$ contain the second-order non-linear terms of the analytic expression of the vector field. All relevant coefficients involved in matrix M and in the non-linear parts contained in $\bar{F}_{i} ,$ $i = 1,2,3$ have been calculated but not reported for the sake of a simple representation.

1.2.5 Hypernormal form simplification

Gamero et al. (1991) detail the procedure to be used to further simplify the non-linear part of system (A2.2) by removing the terms which are not essential for a complete understanding of the local dynamic behavior. The procedure leads to the following second-order hypernormal form

$$ \left( {\begin{array}{*{20}c} {\dot{w}_{1} } \\ {\dot{w}_{2} } \\ {\dot{w}_{3} } \\ \end{array} } \right) = {\mathbf{L}}\left( {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {w_{3} } \\ \end{array} } \right) + {\mathbf{M}}\left( {\mu ,\nu ,\tau } \right)\left( {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {w_{3} } \\ \end{array} } \right) + \left( {\begin{array}{*{20}c} 0 \\ 0 \\ {\bar{a}w_{1} w_{3} + \bar{b}w_{1} w_{2} + \bar{c}w_{1}^{2} } \\ \end{array} } \right) $$

where $\bar{a}$, $\bar{b}$ and $\bar{c}$ are intricate combinations of the original parameters. All necessary coefficients have been calculated but are not reported.

1.2.6 Versal deformation

Consider the linear part of system (A2.2)

$$ {\mathbf{V}}\left( {\mu ,\nu ,\tau } \right) = {\mathbf{L}} + {\mathbf{M}}\left( {\mu ,\nu ,\tau } \right) $$

A candidate for the versal deformation occurs at

$$ \left( {\begin{array}{*{20}c} 0 & 1 & 0 \\ 0 & 0 & 1 \\ {\hat{\varepsilon }_{1} \left( {\mu ,\nu ,\tau } \right)} & {\hat{\varepsilon }_{2} \left( {\mu ,\nu ,\tau } \right)} & {\hat{\varepsilon }_{3} \left( {\mu ,\nu ,\tau } \right)} \\ \end{array} } \right)~ $$

(A2.4)

where $\hat{\varepsilon }_{1} \left( {\nu ,\mu ,\tau } \right) \equiv {\mathbf{Det}}\left( {\mathbf{V}} \right)$, $\hat{\varepsilon }_{2} \left( {\nu ,\mu ,\tau } \right) \equiv - {\mathbf{B}}\left( {\mathbf{V}} \right)$ and $\hat{\varepsilon }_{3} \left( {\nu ,\mu ,\tau } \right) \equiv {\mathbf{Tr}}\left( {\mathbf{V}} \right)$ are the unfolding parameters, and ${\mathbf{V}}$ is topologically equivalent to ${\mathbf{J}}$ (cf. Algaba et al. 2003). It is easy to verify that the transversality condition

$$ \left. {\frac{{\partial \left( {\hat{\varepsilon }_{1} ,\hat{\varepsilon }_{2} ,\hat{\varepsilon }_{3} } \right)}}{{\partial \left( {\tau ,\mu ,\nu } \right)}}} \right|_{{\tau = \mu = \nu = 0}} = \left| {\begin{array}{*{20}c} {\frac{{\partial \hat{\varepsilon }_{1} }}{{\partial \tau }}} & {\frac{{\partial \hat{\varepsilon }_{1} }}{{\partial \mu }}} & {\frac{{\partial \hat{\varepsilon }_{1} }}{{\partial v}}} \\ {\frac{{\partial \hat{\varepsilon }_{2} }}{{\partial \tau }}} & {\frac{{\partial \hat{\varepsilon }_{2} }}{{\partial \mu }}} & {\frac{{\partial \hat{\varepsilon }_{2} }}{{\partial v}}} \\ {\frac{{\partial \hat{\varepsilon }_{3} }}{{\partial \tau }}} & {\frac{{\partial \hat{\varepsilon }_{3} }}{{\partial \mu }}} & {\frac{{\partial \hat{\varepsilon }_{3} }}{{\partial v}}} \\ \end{array} } \right| \ne 0.~ $$

(A2.5)

holds.

1.3 Transformation into phase/amplitude equations

Algaba et al. (2003) shows how a codimension 2 pitchfork–Hopf interaction can be exploited in order to put system ${\mathcal{S}}$ in a form that allows to decuple the hypernormal form into phase/amplitude equations (cf. Wiggins 1991). Consider the following variable transformation

$$ \left( {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ {w_{3} } \\ \end{array} } \right) = \left[ {\begin{array}{*{20}c} 0 & { - 1} & s \\ { - \omega } & 0 & 0 \\ 0 & {\omega ^{2} } & 0 \\ \end{array} } \right]\left( {\begin{array}{*{20}c} x \\ y \\ z \\ \end{array} } \right)~ $$

(A3.1)

where $s$ is a scale parameter, and $\omega = - \hat{\varepsilon }_{2}$. That is

$$ \left( {\begin{array}{*{20}c} {\dot{x}} \\ {\dot{y}} \\ {\dot{z}} \\ \end{array} } \right) = \left[ {\begin{array}{*{20}c} 0 & { - \frac{1}{\omega }} & 0 \\ 0 & 0 & {\frac{1}{{\omega ^{2} }}} \\ {\frac{1}{s}} & 0 & {\frac{1}{{s\omega ^{2} }}} \\ \end{array} } \right]\left( {\begin{array}{*{20}c} 0 & 1 & 0 \\ 0 & 0 & 1 \\ {\hat{\varepsilon }_{1} \left( {\mu ,\nu ,\tau } \right)} & {\hat{\varepsilon }_{2} \left( {\mu ,\nu ,\tau } \right)} & {\hat{\varepsilon }_{3} \left( {\mu ,\nu ,\tau } \right)} \\ \end{array} } \right)\left[ {\begin{array}{*{20}c} 0 & { - 1} & s \\ { - \omega } & 0 & 0 \\ 0 & {\omega ^{2} } & 0 \\ \end{array} } \right]\left( {\begin{array}{*{20}c} x \\ y \\ z \\ \end{array} } \right) + n.l.t.~ $$

(A3.2)

At the Gavrilov-Guckenheimer bifurcation ($\hat{\varepsilon }_{1} = \hat{\varepsilon }_{3} = 0$; $\hat{\varepsilon }_{2} < 0$) system (A3.2) becomes

$$ \left( {\begin{array}{*{20}c} {\dot{x}} \\ {\dot{y}} \\ {\dot{z}} \\ \end{array} } \right) = \left( {\begin{array}{*{20}c} 0 & { - \omega } & 0 \\ \omega & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} } \right)\left( {\begin{array}{*{20}c} x \\ y \\ z \\ \end{array} } \right) + \left( {\begin{array}{*{20}c} {a_{1} x^{2} + a_{2} y^{2} + a_{3} z^{2} + a_{4} xy + a_{5} xz + a_{6} yz} \\ {b_{1} x^{2} + b_{2} y^{2} + b_{3} z^{2} + b_{4} xy + b_{5} xz + b_{6} yz} \\ {c_{1} x^{2} + c_{2} y^{2} + c_{3} z^{2} + c_{4} xy + c_{5} xz + c_{6} yz} \\ \end{array} } \right) $$

(A3.3)

All parameters have been computed but are not displayed for the sake of a contained presentation.

We can finally translate system (A3.3) into cylindrical coordinates $\left( {r,z,\theta } \right)$. Standard arguments (cf., inter al., Wiggins 1991) imply

$$ \begin{aligned} \dot{r} & & & = \hat{a}_{1} rz + \hat{a}_{2} r^{3} + \hat{a}_{3} rz^{2} \\ \dot{z} & & = \hat{b}_{1} r^{2} + \hat{b}_{2} z^{2} + \hat{b}_{3} z^{3} + \hat{b}_{4} r^{2} z \\ \dot{\theta } & = \omega + \hat{c}_{1} z + \hat{c}_{2} r^{2} + \hat{c}_{3} z^{2} \\ \end{aligned} $$

(A3.4)

where $\hat{a}_{i} ,\,i = 1,2,3$, $\hat{b}_{j}$, $j = 1,2,3,4$ and $\hat{c}_{n}$, $n = 1,2,3$ are the so-called resonant non-linear coefficients (e.g., Guckenheimer and Holmes 1983).

Finally, given the $\theta$-independence in the ($\dot{r},\dot{z}$) equations, we can use a standard rescaling technique and derive the decoupled topologically equivalent system

$$ \begin{array}{*{20}r} \hfill {\dot{r} = \mu _{1} r + \tilde{a}rz~~~~~~~~~~~~~~~~~~~~~~} \\ \hfill {\dot{z} = \mu _{2} + kr^{2} - \tilde{b}z^{2} ~~~~k = \pm 1} \\ \end{array} ~ $$

(A3.5)

where $\mu _{1} = \frac{{\hat{\varepsilon }_{3} }}{\kappa }$, $\mu _{2} = \frac{{\hat{\varepsilon }_{1} }}{\kappa }\left| {\frac{b}{a}} \right|$, are the standard unfolding parameters, given $\kappa = \frac{{a^{2} }}{{\left| b \right|}}$ (cf. Harlim and Langford 2007; and Freire et al. 2002). All necessary mappings between the parameters of the planar system (A3.6) and the parameters of the original system ${\mathcal{S}}$ have been computed. Results are not displayed for the sake of a contained presentation but are available upon request.