1 Introduction

The Fermi-Pasta-Ulam (FPU) chain or lattice is an \(n\) degrees-of-freedom (dof) Hamiltonian system that models a chain of oscillators with nearest-neighbour interaction, see [7] and [8]. We will describe the model in Sect. 2. There exists a huge amount of literature on the FPU chain but nearly always regarding the classical case of equal masses, see for recent references [4]. One should note that, although the classical case looks physically natural, its dynamics is non-generic because of symmetries. One of the consequences is that only two dof resonances are effective in such a chain with \(n\) particles; see [13]. In the conclusions we will mention applications and make some observations about the classical FPU chain in the light of the present analysis.

In this paper we will outline a research program to study the case where the masses are different. An inhomogeneous nonlinear lattice with nearest neighbour interaction is studied in [15] with emphasis on energy control. The case of alternating masses was studied in [9] for a FPU chain to obtain insight in the equipartition of energy. A preliminary but important conclusion in [9] is that for the masses considered and on long timescales no equipartition takes place. This has been confirmed by analysis based on symmetries and normal forms in [3].

It is understandable that only a few results were obtained for inhomogeneous lattices as the choice of inhomogeneities, the masses of the lattice, seems to be arbitrary. We will solve this arbitrariness by focusing on the presence of resonances induced by the choice of masses. After referring to some basic material on Hamiltonians and normal forms we formulate in Sect. 2 the periodic FPU \(\alpha \) chain with arbitrary (positive) masses. In such a \(n\) degrees-of-freedom system there exists a momentum integral that enables us to reduce to a \(n-1\) dof system. An inverse problem is considered in Sect. 3: how do we find mass distributions producing prominent resonances in the spectrum induced by \(H_{2}(p, q)\), the quadratic part? This involves the analysis of the inverse map of the vector of mass distribution to the vector of positive eigenvalues of an associated coefficient matrix. This problem is solved in Sect. 3 for the cases of 3 and 4 particles; in the latter case it turns out that of the four 1st order resonances that exist in general (see for the terminology [14]) 3 exist, of the 12 possible 2nd order resonances 10 exist in this FPU chain. In Sect. 4 we focus on the \(1:2:3\) resonance that arises for a one-dimensional variety of mass ratios. It turns out that for one particular combination of mass ratios, the normal form of the nonlinear system is integrable. Moving from this particular case into the variety of mass ratios, one of the periodic solutions shows Hamilton-Hopf bifurcation that corresponds with Shilnikov-Devaney bifurcation in this Hamiltonian system and produces a chaotic normal form.

In the conclusions we draw attention to possible applications of the present paper.

The Appendix contains general statements on the relation between mass ratios and the spectrum induced by \(H_{2}(p, q)\) that can be useful for future research. Table 3 summarises the instructions for the case of 4 particles. It is shown that for a given \(n\)-dimensional eigenvector characterizing the FPU chain, all positive solutions of an \(n\)-dimensional mass distribution are in a compact subset of \({\mathbb{R}}^{n}\). This subset is empty in some cases, for instance the important \(1:1: \ldots :1\) resonance does not arise for the periodic FPU \(\alpha \) chain with four or more particles.

1.1 Hamiltonian Formulation

For an autonomous Hamiltonian system with \(n\) dof, \(n\) independent integrals suffice for integrability, in that case there will be no chaotic motion in such a system. However, in general, Hamiltonian systems with two or more dof are non-integrable. In many cases, this phenomenon was identified with homoclinic chaos as predicted by Poincaré in the nineteenth century, see [11, vol. 3]; for a description see [17, Sects. 5.4 and 9.3].

In the seventies of last century, a number of scientists started with the computation and analysis of normal forms of general Hamiltonian systems near equilibrium. Introductions and surveys of results can be found in [14, Chap. 10] and [18]. We will follow the same terminology.

In the sequel, a periodic solution should be understood as a periodic solution for a fixed value of the energy (iso-energetic solution), so actually it corresponds for the full Hamiltonian system with a family of periodic solutions parameterized by the energy.

Normal form computations for Hamiltonian systems can be carried out in various ways. Apart from efficiency, the main point is to keep the system energy-preserving and preferably canonical. Using for instance polar coordinates like action-angle coordinates or amplitude-phase coordinates one can perform averaging over the angles or explicitly time to obtain a first-order normal form. One may consult [14] or [18] for more details. An introductory text is [16, Chaps. 11 and 12].

In Sect. 4 we will analyze periodic \(\alpha \)-chains (FPU chains where the Hamiltonian is truncated after the cubic terms), containing the \(1:2:3\) resonance with main objective to investigate the stability of the short-periodic solutions on the energy manifold and the integrability of the normal form. This is highly relevant for the characterization of the chaotic dynamics of the system but, as mentioned above, it raises special problems. In the cases of vanishing actions or amplitudes, for instance when studying normal modes, the procedure will be as follows (see also Sect. 4.1). Starting with the equations of motion, we will use co-moving coordinates (see for instance transformation (11.9-10) in Chap. 11 of [16]) to obtain a first order normal form. This normal form is used to localize the short-periodic solutions; the normal form conserves the energy but the transformation is not canonical. We will use averaging-normalization as it yields rigorous approximation results (see [14]), the results are qualitatively and quantitatively precise. The same holds when we use polar coordinates outside the coordinate planes.

In Sect. 4, the short-periodic solutions can be computed explicitly. The next step is then to linearize near the periodic solutions and to determine the Floquet exponents for which we have to study coupled Mathieu equations. This is still a formidable task, but we can obtain a first order approximation of the exponents by normalizing the coupled Mathieu equations. This will give a number of spectral stability results in Sect. 4.

In the Appendix we outline the general procedure to obtain mass ratios that correspond with resonance. We also give indications of the extension to periodic FPU chains with \(n\) particles.

1.2 Outline of a Research Programme

The original Fermi-Pasta-Ulam chain [7] consists of \(n\) oscillators of equal mass with nearest-neighbour interaction. In a neighbourhood of equilibrium, the spectrum of the linear part of the equations of motion plays a crucial part regarding the nature of the ensuing dynamics, see for instance [14] or [18]. Considering inhomogeneous mass distributions in FPU chains, one can produce a great many different spectra induced by \(H_{2}\). Each of these cases may produce different dynamics in the corresponding FPU chain. In Sect. 3 we will consider resonant spectra for the case of three and more extensively four particles with periodic boundary conditions i.e. chains where the first and the last oscillator are connected. For the case of four particles we will focus on the rich dynamics of the \(1:2:3\) resonance. An outline of possible further research follows here:

  1. 1.

    According to Table 1 regarding the case of four particles, we also have to study two first order resonances (\(1:2:1\) and \(1:2:4\)) and ten second order resonances. Also, higher order resonances may be worthwhile to investigate. Special attention should be given to the \(1:1:2\) and \(1:1:3\) cases as only four special mass ratios produce these resonances. In such a case degenerations may arise so that we have to consider detuning phenomena, see [14].

    Table 1 Fibers of resonances
    Full size table
  2. 2.

    Cases of five and more particles will present many more problems.

  3. 3.

    The present study is restricted to so-called periodic \(\alpha \)-chains. Including quartic terms in the Hamiltonian (\(\beta \)-chains) and considering lattices with fixed begin- and end-point will produce new results.

  4. 4.

    The study presented here and possibly future studies will throw light on qualitative and quantitative differences between systems in nearest-neighbour interaction and non-local interaction, a topic that is relevant for plasma physics and stellar dynamics.

2 The Fermi-Pasta-Ulam Chain

For the mono-atomic case of the original periodic FPU-problem (all masses equal) it was shown in [12] for up to six degrees-of-freedom (dof) and much more general in [13], that the corresponding normal forms are governed by \(1:1\) resonances and that these Hamiltonian normal forms are integrable. This explains the recurrence phenomena near equilibrium, it also shows that the classical FPU chain has symmetries to make the problem non-generic.

We will drop the original assumption of identical (mono-atomic) particles to consider the periodic FPU-problem again. For \(n\) particles with mass \(m_{j}>0\), position \(q_{j}\) and momentum \(p_{j}= m_{j} \dot{q}_{j}, j= 1 \ldots n\), \(\varepsilon \geq 0\) a small parameter, the Hamiltonian is of the form:

$$ H(p, q)= \sum_{j=1}^{n} \biggl(\frac{1}{2m_{j}}p_{j}^{2} + V(q_{j+1}-q_{j}) \biggr) \quad \mbox{with}\ V(z)= \frac{1}{2}z^{2}+ \varepsilon \frac{\alpha }{3}z^{3} + \varepsilon^{2} \frac{\beta }{4}z^{4}. $$
(1)

The quadratic part of the Hamiltonian is not in diagonal form; for \(n=3, 4 \ldots \) the linearized equations of motion can be written as:

$$\begin{aligned} \left\{\textstyle\begin{array}{ll} m_{1} \ddot{q}_{1} +2q_{1}-q_{2} -q_{n} & = 0, \\ m_{2} \ddot{q}_{2} +2q_{2}-q_{3} -q_{1} & = 0, \\ m_{3} \ddot{q}_{3} +2q_{3}-q_{4} -q_{2} & = 0, \\ \ldots & = 0, \\ m_{n} \ddot{q}_{n} +2q_{n}-q_{1} -q_{n-1} & = 0. \end{array}\displaystyle \right. \end{aligned}$$
(2)

We can write for the quadratic part of \(H(p, q)\):

$$ H_{2} = \frac{1}{2}p^{T}A_{n} p + \frac{1}{2} q^{T}C_{n} q, $$
(3)

with \(A_{n}\) the \(n\times n\) diagonal matrix with at position \((i,i)\) the value \(m_{i}^{-1} =: a_{i}\), \(C_{n}\) is an \(n \times n\) matrix. For an analysis of the quadratic term \(H_{2}(p, q)\) we need to know the eigenvalues of \(A_{n}C_{n}\). The relation between the eigenvalues of \(A_{n}C_{n}\) and the eigenvalues of the matrix of coefficients of system (2) will be given below. Since the null space of \(C_{n}\) has dimension one, the matrix \(A_{n}C_{n}\) has an eigenvalue 0 corresponding to a (translational) momentum integral. It will turn out that the other eigenvalues of \(A_{n}C_{n}\) are positive, as expected. For a given set of masses, the calculation of the remaining eigenvalues corresponding with the frequencies of the linearized system is easy, but we are faced with another, an inverse problem. To focus ideas, suppose that \(n=4\). The presence of the momentum integral implies that we have to consider a three degrees-of-freedom (dof) Hamiltonian problem. We know, see for instance [14, Chap. 10] or [18], that the first order resonances are \(1:2:1, 1:2:2, 1:2:3\) and \(1:2:4\). The question is then if and how we can choose the masses so that these prominent resonances are present. Of course, this problem will be more formidable if \(n>4\). In the next section we determine for \(n=4\) the ratios of masses that produce the resonance \(1:2:3\). The approach works equally well for other prescribed rations of eigenvalues, as we discuss in the Appendix. Prominent resonances for \(n > 4\) can be found but a systematic study of these cases poses a difficult open algebraic problem.

3 The Spectrum Induced by \(H_{2}\)

After a number of general considerations we will give details for the cases of three and four particles. The first case is rather trivial as far as the spectrum goes, the case of four particles is already quite complicated. Here we mention the main facts that we need in the later sections. In the Appendix we will give more details.

3.1 The Matrix for Inhomogeneous FPU-Lattices and Its Eigenvalues

The linear system (2) can be written as

$$ \left(\textstyle\begin{array}{c} \dot{q} \\ \ddot{q} \end{array}\displaystyle \right) = M \left(\textstyle\begin{array}{c} q \\ \dot{q} \end{array}\displaystyle \right) ,\quad \mbox{where } M= \left(\textstyle\begin{array}{c@{\quad }c} 0 & I_{n} \\ -A_{n}C_{n} & 0 \end{array}\displaystyle \right) , $$
(4)

where the matrix \(A_{n}\) is a diagonal matrix with the inverse masses \(m_{j}^{-1}=:a_{j}\) on the diagonal, and where the matrix \(C_{n}\) has elements 2 on the diagonal, and −1 at positions \((i,i+1)\) and \((i,i-1)\), with the indices taken modulo \(n\). For instance,

$$ C_{5}= \left(\textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} 2&-1&0&0&-1 \\ -1&2&-1&0&0 \\ 0&-1&2&-1&0 \\ 0&0&-1&2&-1 \\ -1&0&0&-1&2 \end{array}\displaystyle \right) . $$

(This matrix turns up elsewhere in mathematics. It is the affine Cartan matrix of the completed root system \(\bar{A}_{n}\). See e.g. [2, Déf. 3 in 1.5 of Chap. 6, and Planche I].)

The \((2n)\times (2n)\) matrix \(M\) has a double eigenvalue 0, corresponding to the momentum integral

$$ \sum_{j=1}^{n} m_{j} \dot{q}_{j} = \mbox{constant}. $$
(5)

In the sequel we will choose the case of vanishing momentum integral which is not a restriction of generality. If \(\lambda \) is a positive eigenvalue of \(A_{n}C_{n}\), then \(i\sqrt{\lambda }\) and \(-i\sqrt{ \lambda }\) are eigenvalues of \(M\), corresponding to frequencies of eigenmodes of the linearized system. So it is useful to collect results concerning the eigenvalues of \(A_{n}C_{n}\).

Proposition 3.1

For \(n= 3, 4 \ldots \) the matrix \(A_{n}C_{n}\) has one eigenvalue 0 and \(n-1\) positive eigenvalues \(\lambda_{1},\ldots , \lambda_{n-1}\), possibly coinciding. If eigenvalues coincide the corresponding eigenspace has maximal dimension.

Proof

Since the \(a_{j}=m_{j}^{-1}\) are positive, the matrix \(A_{n}^{1/2}\) is well-defined. The symmetric matrix \(A_{n}^{1/2}C_{n}A_{n}^{1/2}\) has real eigenvalues, and the algebraic and geometric multiplicities of eigenvalues coincide.

If \(y\) is an eigenvector of \(A_{n}C_{n}\) with eigenvalue \(\lambda \), then \(\lambda A_{n}^{-1}y= C_{n} y\), and the scalar product with \(y\) for both sides of the equation gives the equality

$$ \lambda \sum_{i} a_{i}^{-1} y_{i}^{2} = 2\sum_{i} y_{i}^{2} - 2 \sum_{i} y_{i} y_{i-1} . $$

(Indices taken modulo \(n\).) So

$$ \sum_{i} (2-\lambda /a_{i}) y_{i}^{2} =2 \sum_{i} y_{1} y_{i-1} . $$

With Schwarz’s inequality applied to the vector \((y_{i})\) and the vector \((y_{i-1})\) this implies \(\lambda \geq 0\). Equality occurs only if the vectors \((y_{i})\) and \((y_{i-1})\) are positive multiples of each other, which occurs only for multiples of \((1,1,\ldots ,1)\). □

For the investigation of the linearized problem we need to understand the map \({\mathbb{R}}_{>0}^{n} \rightarrow {\mathbb{R}}_{>0}^{n-1}\), from a vector \((a_{1},\ldots , a_{n})\) of inverse masses to a vector \((\lambda_{1},\ldots ,\lambda_{n-1})\) of positive eigenvalues. The order of the eigenvalues is not determined, so we have, more precisely, a map \(\rho_{n}: {\mathbb{R}}_{>0}^{n} \rightarrow S_{n-1}\backslash {\mathbb{R}} _{>0}^{n-1}\), with the action of the symmetric group \(S_{n-1}\) on the coordinates. For the linearized inhomogeneous FPU-chain described by system (2), the dihedral group \(D_{n}\) with \(2n\) elements permutes the coordinate \(q_{j}\) (generated by a shift and a reflection). This transforms system (2) into an equivalent system. Another symmetry is by scaling: \(\rho ( t(a_{1},\ldots ,a_{n}) ) = t \rho (a_{1},\ldots ,a_{n})\) for \(t>0\).

To investigate the correspondence between eigenvalues and inverse masses we use the equality

$$ \det ( A_{n} C_{n} - \lambda I_{n}) = -\lambda \prod (\lambda _{j}-\lambda ) , $$

for \((\lambda_{1},\ldots ,\lambda_{n-1})=\rho (a_{1},\ldots ,a_{n})\). Comparing the coefficients of this polynomial we obtain the equalities

$$ p_{j} (A_{n}) = e_{n-j}\bigl(\{ \lambda_{1},\ldots ,\lambda_{n-1}\}\bigr) \quad (1\leq j \leq n-1) , $$
(6)

with the elementary symmetric functions \(e_{k}\) and homogeneous polynomials \(p_{j}(A_{n})\) in the \(a_{j}\) of degree \(n-j\). This describes the structure of the set of diagonal matrices \(A_{n}\) for a prescribed spectrum of \(A_{n} C_{n}\). It is the set of points with positive coordinates in an algebraic set in \(\mathbb{C}^{n}\) which is the intersection of \(n-1\) hyperplanes given by equations of degree \(1, 2, \ldots , n-1\).

In Lemma A.2 in the Appendix we’ll show that

$$ p_{n-1}(A_{n}) = 2\sum _{i} a_{i} ,\qquad p_{n-2}(A_{n}) = \sum_{1\leq i < j\leq n} c_{i,j} a_{i}a_{j} , $$
(7)

with \(c_{i,j}=3\) if \(i-j=\pm 1\bmod n\), and \(c_{i,j}=4\) otherwise. All \(p_{j}(A_{n})\) are invariant under the action of the dihedral group \(D_{n}\) on the coordinates \(a_{j}\).

In the Appendix we will also show that all real solutions \((a_{1}, \ldots ,a_{n})\) for a given eigenvalue vector \((\lambda_{1},\ldots , \lambda_{n-1})\) are in a compact subset of \({\mathbb{R}}^{n}\). This subset may be empty. For all \(n\geq 4\) the \(1:1:\cdots :1\) resonance does not occur for any mass distribution.

3.2 The Case of Three Particles

For \(n=3\) the determination of the eigenvalues for given inverse masses amounts to solving the quadratic equation

$$ \lambda^{2} -2(a_{1}+a_{2}+a_{3}) \lambda + 3(a_{1}a_{2}+a_{1}a_{3}+a _{2}a_{3})=0, $$

which has positive solutions.

Conversely, for all choices \((\lambda_{1},\lambda_{2})\) of positive eigenvalues, values of \(a_{1},a_{2},a_{3}\) can be found such that \(A_{3}C_{3}\) has eigenvalues \(\lambda_{1}\), \(\lambda_{2}\) and 0. If \(\lambda_{1}=\lambda_{2}\) there is exactly one solution \(a_{1}=a_{2}=a _{3}=\frac{1}{3} \lambda \) (equal masses). If the eigenvalues have ratio \(\lambda_{1}/\lambda_{2}>1\) then the corresponding points \((a_{1},a _{2},a_{3})\) in \({\mathbb{R}}^{3}\) form an ellipse. This ellipse may or may not be contained in the positive octant. See Fig. 1.

Fig. 1
figure 1

Solutions sets of inverse masses \((a_{1},a_{2},a_{3})\) for the FPU-chain with three particles. For the eigenvalue ratio \(\lambda_{1}/\lambda_{2}=2\) of matrix \(A_{n}C_{n}\) the solution set is compact; it is an ellipse in \({\mathbb{R}}_{>0}^{3}\). For the eigenvalue ratio \(\lambda_{1}/\lambda_{2}=4\) the solutions are on a larger ellipse in \({\mathbb{R}}^{3}\), which intersects \({\mathbb{R}}_{>0}^{3}\) in three open curves

Full size image

3.3 The Case of Four Particles

In the case \(n=4\) we use the scaling to restrict our further investigation to eigenvalues satisfying \(\lambda_{1}+\lambda_{2}+ \lambda_{3}=1\). Working out the polynomials in (6) for the case \(n=4\) we obtain three equations for a given vector \((\lambda_{1}, \lambda_{2},\lambda_{3})\in {\mathbb{R}}_{>0}^{3}\):

$$ \begin{aligned} 4(a_{1}a_{2}a_{3}+a_{2}a_{3}a_{4}+a_{3}a_{4}a_{1}+a_{4}a_{1}a_{2}) &= \lambda_{1}\lambda_{2}\lambda_{3} =: \xi , \\ 3(a_{1}a_{2}+a_{2}a _{3}+a_{3}a_{4}+a_{4}a_{1})+4(a_{1}a_{3}+a_{2}a_{4})&= \lambda_{1} \lambda_{2} +\lambda_{2} \lambda_{3}+ \lambda_{3}\lambda_{1} =: \eta , \\ 2(a_{1}+a_{2}+a_{3}+a_{4})&= 1 . \end{aligned} $$
(8)

We call the set of \((a_{1},\ldots ,a_{4})\in {\mathbb{R}}_{>0}^{4}\) satisfying these relations the fiber of \((\xi ,\eta )\in {\mathbb{R}} _{>0}^{2}\). In Proposition A.5 in the Appendix we will give a precise characterization of the set of \((\xi ,\eta )\) for which the fiber is non-empty.

The resonances deserve special attention. A resonance \(n_{1}:n_{2}:n _{3}\) in the linearized system (2) corresponds to an eigenvalue vector of \(A_{4}C_{4}\) with the ratios \(n_{1}^{2} : n_{2} ^{2} : n_{3}^{2}\). We considered all resonances of order one and two, and obtained the results in Table 1. As noted in Sect. 1.2, the resonances \(1:1:2\) and \(1:1:3\) need special attention.

In Sect. A.2.2 in the Appendix we describe how to determine the fiber explicitly by choosing an additional parameter and solve the system of equations (8). The idea is to determine successively quantities that are invariant under a decreasing sequence of subgroups of the dihedral group \(D_{4}\) acting by permutation on the vector \((a_{1},a_{2},a_{3},a_{4})\).

3.4 The Resonance \(1:2:3\)

Here we consider the resonance that is the subject of study in the next section.

By scaling we arrange \(\lambda_{1}=\frac{9}{14}\), \(\lambda_{2}= \frac{2}{7} \), \(\lambda_{3}=\frac{1}{14}\) to satisfy the last equation in (8). We follow the computational scheme in Table 3 in the Appendix. This leads to

$$ \begin{aligned} a_{1},a_{3} &= \frac{2+u}{56}\mp \frac{\sqrt{2}}{112} \sqrt{ \frac{u(6-u)(16-u)}{5-u}} , \\ a_{2},a_{4}&= \frac{12-u}{56} \mp \frac{1}{56 \sqrt{2}} \sqrt{\frac{(6+u)(4-u)(10-u)}{5-u}} , \end{aligned} $$
(9)

where we take the minus sign for \(a_{1}\) and \(a_{2}\). The parameter \(u\) runs through the interval \([0,u_{1})\) with

$$ u_{1}= \frac{8}{3} - \frac{2}{3}\sqrt[3]{19} \approx 0.887732 . $$
(10)

Figure 2 gives a plot.

Fig. 2
figure 2

One branch of the fiber for the resonance \(1:2:3\) is given by the functions \(a_{1}\leq a_{3}< a_{2}<a_{4}\) on the interval \([0,u_{1})\). (Horizontal axis: parameter \(u\); vertical axis: values of \(a_{j}(u)\).) The three dots on the horizontal axis correspond to the values \(0, 0.534105\) and 0.826713 of the parameter \(u\) for which we will carry out simulations in the next section (the cases 0, 1 and 2)

Full size image

We get all solutions when we let the dihedral group \(D_{4}\) act on the solutions that we constructed. The branch in Fig. 2 and its image under \([1,3]\in D_{4}\) have the point parametrized by \(u=0\) in common.

3.5 Illustration of the Fiber

The equations in (9) describe a curve \(u\mapsto (a_{1}(u),a _{2}(u),a_{3}(u),a_{4}(u) )\) in \({\mathbb{R}}_{>0}^{4}\) corresponding to a one-parameter family of solutions for the inverse masses. To illustrate it we use the second and last equation in (8), which describe an ellipsoid in the hyperplane \(a_{1}+a _{2}+a_{3}+a_{4}=\frac{1}{2}\). In (53) and (54) in the Appendix we’ll describe this ellipsoid in a more explicit way. The first equation in (8) produces an intersection with this ellipsoid in some curves. The points with positive coordinates in this intersection form the fiber.

On this ellipsoid we can use a system of spherical coordinates, mapping the ellipsoid to the rectangle \([-\pi ,\pi ]\times [-\frac{1}{2} \pi ,\frac{1}{2}\pi ]\), with boundary identifications. The image of the fiber under this map is given in Fig. 3.

Fig. 3
figure 3

The fiber for the resonance \(1:2:3\) is contained in an ellipsoid, which we describe in spherical coordinates. (Horizontally the azimuth \(\phi \), and vertically the inclination \(\psi \). See (54).) The thick curve corresponds to the branch of the fiber in Fig. 2. The dotted curves correspond to the translates of this branch under the dihedral group \(D_{4}\). The thin curves indicate the boundary of the region corresponding to coordinates in \({\mathbb{R}} _{>0}^{4}\). The picture illustrates that the fiber for \((1:2:3)\) consists of four open curves, and that (9) describes a fundamental domain for the action of the dihedral group on the fiber

Full size image

3.6 Transformation of the Hamiltonian

We form the diagonal matrix \(A_{4}(u)\) with diagonal elements \(a_{j}(u)\), \(1\leq j \leq 4\). In the proof of Proposition 3.1 we noted that \(A_{4}(u)^{1/2}C_{4} A_{4}(u)^{1/2}\) is a symmetric matrix (as long as \(u\in [0,u_{1})\)), so we can find an orthogonal matrix \(U(u)\) such that \(A_{4}(u)^{1/2}C_{4} A_{4}(u)^{1/2}= U(u) \varLambda U(u)^{T}\), where \(\varLambda \) is the diagonal matrix with diagonal elements \(\frac{9}{14}\), \(\frac{2}{7}\), \(\frac{1}{14}\), and 0. Then the transformation matrices

$$ K(u) = A_{4}(u)^{-1/2} U(u) , \qquad L(u) = A_{4}(u)^{1/2}U(u) $$
(11)

determine a symplectic transformation

$$ p=K(u) y ,\qquad q=L(u)x , $$
(12)

which transforms the quadratic part in (3) of the Hamiltonian into

$$ H_{2} = \frac{1}{2} y^{T} y + \frac{1}{2} x^{T} \varLambda x = \frac{1}{2}\sum _{j=1}^{4} \bigl( y_{j}^{2} + \lambda_{j} x_{j}^{2} \bigr) . $$
(13)

This will produce the so-called quasi-harmonic form of the equations of motion. To see that \(H_{2}\) takes the form (13) we need the existence of an orthogonal matrix \(U(u)\) diagonalizing \(A_{4}(u)^{1/2} C_{4} A_{4}(u)^{1/2}\). We do not need to know \(U(u)\), \(K(u)\) or \(L(u)\) explicitly.

To transform the cubic and higher order terms of the Hamiltonian to coordinates corresponding to the eigenmodes of the linearized system we need to know the transformation matrix \(L(u)\) explicitly. It is no problem to do this numerically with Mathematica or Matlab for any given \(u\in [0,u_{1})\). It is nicer to have \(U(u)\), and hence \(L(u)\) and \(K(u)\), symbolically in terms of the parameter \(u\). Lemma A.6 in the Appendix describes the method that we used to obtain a symbolic description.

For the cubic term we note that (with indices modulo 4)

$$ \frac{1}{3} \sum_{j} (q_{j+1}-q_{j} )^{3} = \sum_{j} (q _{j+1}-q_{j-1} )q_{j}^{2} . $$
(14)

The substitution \((q_{1},\ldots ,q_{4})^{T} = L(u) ( x_{1}, \ldots ,x_{4})^{T}\) gives

$$\begin{aligned} H_{3} =& \varepsilon \bigl( d_{1}(u)x_{1}^{3}+d_{2}(u)x_{1}^{2}x_{2}+d _{3}(u)x_{1}^{2}x_{3} +d_{4}(u)x_{2}^{2}x_{1}+d_{5}(u)x_{3}^{2}x_{1} \\ &{} +d_{6}(u)x_{1}x_{2}x_{3} +d_{7}(u)x_{2}^{3}+d_{8}(u)x_{3} ^{3}+ d_{9}(u)x_{3}^{2}x_{2}+d_{10}(u)x_{2}^{2}x_{3} \bigr) , \end{aligned}$$
(15)

with the functions \(d_{j}\) as indicated in Table 2.

Table 2 Coefficients of the cubic term of the Hamiltonian, transformed to eigenmodes, the so-called quasi-harmonic form. For \(u=0\) only \(d_{3}\), \(d_{6}\) and \(d_{10}\) are non-zero
Full size table

4 The \(1:2:3\)-Resonance for the Periodic \(\alpha \)-Lattice (\(n=4\))

For any possible inhomogeneous FPU \(\alpha \)-chain with four dof we have the system:

$$\begin{aligned} \textstyle\begin{cases} \dot{q_{1}} = v_{1}, \\ \dot{v}_{1} = [-2q_{1}+q_{2} +q_{4} - \varepsilon \alpha ((q_{1}-q _{4})^{2}-(q_{2}-q_{1})^{2}) ] a_{1}, \\ \dot{q_{2}} = v_{2}, \\ \dot{v}_{2} = [-2q_{2}+q_{3} +q_{1} - \varepsilon \alpha ((q_{2}-q _{1})^{2}-(q_{3}-q_{2})^{2} )] a_{2}, \\ \dot{q_{3}} = v_{3}, \\ \dot{v}_{3} = [-2q_{3}+q_{4} +q_{2} - \varepsilon \alpha ((q_{3}-q _{2})^{2} -(q_{4}-q_{3})^{2})] a_{3}, \\ \dot{q_{4}} = v_{4}, \\ \dot{v}_{4} = [-2q_{4}+q_{1} +q_{3}- \varepsilon \alpha ((q_{4}-q _{3})^{2} -(q_{1}-q_{4})^{2})] a_{4}, \end{cases}\displaystyle \end{aligned}$$
(16)

The coefficient \(\alpha \) has been retained for reference to the literature; here we will take \(\alpha =1\). If \(a_{1}= \cdots =a_{4}=1\), we have the classical periodic FPU chain with four particles; it was shown in [12], that in this case the normal form is integrable. The implication is that for \(\varepsilon \) small, chaos is negligible in this classical case. The case of the \(1:2:3\)-resonance is strikingly different as only if we have two masses equal, the normal form is integrable.

Apart from the Hamiltonian we have from (5) as a second (momentum) integral:

$$ m_{1}v_{1}+m_{2}v_{2}+m_{3}v_{3}+m_{4}v_{4} = \mbox{constant}. $$
(17)

The presence of the momentum integral results in two zero eigenvalues of the matrix M in Eq. (4), so by reduction we have to deal essentially with a three dof system. This will be used and explicitly shown in the next three subsections.

According to Table 1 the \(1:2:3\) resonance is present among the possible inhomogeneous FPU lattices. Figure 2 gives one branch of values of inverse masses \(a_{1}, \ldots , a_{4}\) producing this resonance. All vectors \((a_{1}, \ldots , a_{4})\) are obtained by the action of the dihedral group \(D_{4}\) on the coordinates and the scaling \((a_{1}, \ldots , a_{4}) \mapsto (ta_{1}, \ldots , ta_{4})\) with \(t>0\).

Table 1 and Fig. 2 show that the \(1 : 2 : 3\)-resonance appears in one case with relatively well-balanced masses, two of which are equal. We denote this by case 0; it will turn out in Sect. 4.1 that this case is quite special dynamically. The other cases are less balanced regarding the masses. Case 0 corresponds to \(u = 0\); as \(u\) increases (we have \(0 \leq u < u_{1}\) with \(u_{1}=0.887732\)), the masses get less well-balanced, one of them tending to infinity. We study the dynamical behaviour in Sect. 4.2. For numerical simulations we have singled out two more cases indicated in Fig. 2.

The expression for the quadratic part of the Hamiltonian \(H_{2}\) is:

$$ H_{2} = \frac{1}{2} \sum _{i=1}^{4} a_{i}p_{i}^{2} + \frac{1}{2} \bigl[(q _{2}-q_{1})^{2}+(q_{3}-q_{2})^{2}+(q_{4}-q_{3})^{2}+(q_{1}-q_{4})^{2} \bigr]. $$
(18)

\(H_{2}\) is a first integral of the linear system (2), it is also a first integral of the normal form of the full system (16). When using \(H_{2}\) from the solutions of the truncated normal form

$$ \bar{H}(p, q) = H_{2}(p, q) + \varepsilon \bar{H}_{3}(p, q), $$

we obtain an \(O(\varepsilon )\) approximation of the (exact) \(H_{2}(p(t), q(t))\) valid for all time; for a proof see [14, Chap. 10]. Note that in the equations we find it convenient to use the non-canonical velocities instead of the momenta. Using the expression \(H_{2}(p(t), q(t))\) for the solutions of the full system (16) shows the accuracy of the normal form and gives an impression of the nature of the dynamics.

The normal form \(\bar{H}_{3}(p, q)\), written in action-angle coordinates or amplitude-phase coordinates (see below), will contain certain combination angles corresponding with the resonance. If \(\bar{H}_{3}\) contains only one combination angle, we have an additional integral of motion and the normal form \(H_{2}+ \bar{H}_{3}\) is integrable. In the case of two or more independent combination angles, we have to investigate the (non-)integrability of the normal form.

To display the quantitative aspects of the solutions we have the possibility of drawing an energy- or action-simplex or as an alternative to produce a time series for explicit solutions or integrals of the normal forms. Both techniques will be used.

As the short-periodic solutions have constant actions (or constant radii in polar coordinates), the integral \(H_{2}\) of the normal form produces for fixed energy an action-simplex with short-periodic solutions represented by points; the actions \(\tau_{i}\) and the polar coordinates \(r_{i}\) are related. One way of displaying the position of short-periodic solutions and their stability on the 5-dimensional energy manifold is the use of this action-simplex with normal modes at the vertices and solutions in the coordinate planes at the sides. The interior of the faces may contain short-periodic solutions in general position. Their stability is indicated by \(E\) (elliptic i.e. imaginary eigenvalues), \(H\) (hyperbolic i.e.real eigenvalues) and \(C\) (complex eigenvalues with real parts non-zero). See for instance for the action simplices displaying periodic solutions Fig. 6.

4.1 Case 0: The FPU Chain with Well-Balanced Masses

In this case we have the \(1:2:3\) resonance with mass values that are as much as possible similar; we have with \(u = 0\) in (9):

$$ a_{1}= 0.0357143, \quad \quad a_{2}= 0.126804,\quad \quad a_{3}= 0.0357143, \quad \quad a_{4}= 0.301767. $$

Note that \(a_{1}=a_{3}\). We checked numerically that the time series \(H_{2}(p(t), q(t))\) based on the original formulation of system (16) and the time series obtained from the transformed Hamiltonian (19) produce the same result as it should.

To put system (16) in the standard form of quasi-harmonic equations we have to apply the symplectic transformation \(p = K(0)y, q = L(0)x\) in (12). This leads with (15) and Table 2 to the transformed Hamiltonian

$$ H(y, x) = \frac{1}{2} \sum_{i=1}^{4} \bigl(y_{i}^{2}+ \omega_{i}^{2} x_{i} ^{2}\bigr)+ \varepsilon \bigl(d_{3}x_{1}^{2}+d_{10}x_{2}^{2}+d_{6}x_{1}x_{2} \bigr)x _{3}, $$
(19)

with

$$\begin{aligned} &\omega_{1}^{2}= \frac{9}{14},\quad \quad \omega_{2}^{2}= \frac{4}{14},\quad \quad \omega _{3}^{2}= \frac{1}{14},\quad \quad \omega_{4}^{2}=0,\quad \quad d_{3}= -9 \frac{\sqrt{21}}{490},\\ &d_{10}= 2 \frac{\sqrt{21}}{245},\quad \quad d_{6}= -3 \frac{ \sqrt{14}}{490}. \end{aligned}$$

Rescaling time \(t/ \sqrt{14} \rightarrow t\), the equations of motion for the three dof system become:

$$\begin{aligned} \textstyle\begin{cases} \ddot{x}_{1} + 9 x_{1} = - \varepsilon 14(2d_{3}x_{1}x_{3} + d_{6}x _{2}x_{3}), \\ \ddot{x}_{2} + 4 x_{2} = - \varepsilon 14(2d_{10}x_{2}x_{3} + d _{6} x_{1}x_{3}), \\ \ddot{x}_{3} + x_{3} = - \varepsilon 14(d_{3}x_{1}^{2} + d_{10}x _{2}^{2} +d_{6} x_{1}x_{2}). \end{cases}\displaystyle \end{aligned}$$
(20)

According to the Weinstein [20] result there exist at least three families of short-periodic solutions of system (20). Inspection of the equations provides us directly with one family given by:

$$ x_{1}(t)=\dot{x}_{1}(t)=x_{2}(t)= \dot{x}_{2}(t)=0, \quad \quad \ddot{x}_{3} + x _{3}=0. $$
(21)

For fixed energy we refer to this periodic solution as the \(x_{3}\) normal mode; to find such an exact solution explicitly is slightly unusual, the solution is harmonic. Additional periodic solutions are obtained as approximations from normal forms as in [10]. In general, when normalizing a three dof system, one recovers the three actions and one expects to find the angles in combinations according to the actual resonances. For the \(3:2:1\) resonance these are to first order after normalization the so-called combination angles \(\phi_{1}- \phi_{2} - \phi_{3}\) and \(2 \phi_{3} - \phi_{2}\). At second order the combination angle \(\phi_{1} - 3 \phi_{3}\) will arise etc., for details see Sect. 10.2.1 of [14]; for instance the term ‘genuine resonance’ associated with the so-called ‘annihilators’ of \(H_{2}\) can be found in definition 10.2.2 of [14].

Computing the normal form (\(H_{2}+ \varepsilon \bar{H}_{3}\)) of system (20) to \(O(\varepsilon )\) as in [10] or [14] and as we shall explicitly show below, only the \(d_{6}\) term survives in \(\bar{H}_{3}\); this makes the Hamiltonian (19) non-generic. An intermediate normal form of the equations of motion becomes:

$$\begin{aligned} \textstyle\begin{cases} \ddot{x}_{1} + 9 x_{1} = - \varepsilon 14d_{6}x_{2}x_{3}, \\ \ddot{x}_{2} + 4 x_{2} = - \varepsilon 14d_{6}x_{1}x_{3}, \\ \ddot{x}_{3} + x_{3} = - \varepsilon 14d_{6}x_{1}x_{2} . \end{cases}\displaystyle \end{aligned}$$
(22)

System (22) is called intermediate as it has still to be normalized, but the omitted terms play no part in the normalization to first order. There is a lot of freedom in choosing coordinate systems to compute the normal form of the equations of motion. Near the coordinate planes, in particular to study the stability of the normal modes, we will use co-moving coordinates. Away from the coordinate planes (solutions in general position), action-angle variables or polar coordinates are easier to handle than co-moving coordinates. For general position orbits we will use in system (20) transformations \(x_{i}, \dot{x}_{i} \rightarrow r _{i}, \psi_{i}\) of the form:

$$ x_{i} = r_{i} \cos (\omega_{i} t + \psi_{i}), \quad \quad \dot{x}_{i} = - r_{i} \omega_{i} \sin (\omega_{i} t + \psi_{i}). $$
(23)

The actions \(\tau_{i}\) are related to the \(r_{i}^{2}\), the angles \(\phi_{i}\) to the arguments \((\omega_{i} t + \psi_{i})\). Putting \(\chi = \psi_{1}- \psi_{2}- \psi_{3}\) and averaging over time \(t\), the averaging-normal form equations outside the coordinate planes become:

$$\begin{aligned} \textstyle\begin{cases} \dot{r}_{1} = \varepsilon \frac{7}{6}d_{6}r_{2}r_{3} \sin \chi , \\ \dot{r}_{2} = - \varepsilon \frac{7}{4}d_{6}r_{1}r_{3} \sin \chi , \\ \dot{r}_{3} = - \varepsilon \frac{7}{2}d_{6}r_{1}r_{2} \sin \chi , \\ \dot{\chi } = \varepsilon \frac{7}{2}d_{6} \frac{\cos \chi }{r_{1}r _{2}r_{3}} ( \frac{r_{2}^{2}r_{3}^{2}}{3} - \frac{r_{1}^{2}r_{3} ^{2}}{2} - \frac{r_{1}^{2}r_{2}^{2}}{1} ) . \end{cases}\displaystyle \end{aligned}$$
(24)

The integral \(H_{2}\) of the normal form equations becomes:

$$ 9r_{1}^{2}+4r_{2}^{2}+r_{3}^{2} = 2E_{0}, $$
(25)

with \(E_{0}\) a positive (energy) constant. The combination angle \(2 \phi_{3} - \phi_{2}\) is missing; another integral of the normal form (24) is:

$$ 2r_{2}^{2} - r_{3}^{2} = C\quad (\mbox{constant}). $$
(26)

In the original variables this integral is:

$$ 2x_{2}^{2}+ \frac{1}{2}\dot{x}_{2}^{2}-x_{3}^{2}- \dot{x}_{3}^{2} = \mbox{constant}. $$

As we have three independent integrals of the normal form equations (24), the normal form is integrable. Because of the approximative character of the normal form, this means that chaotic motion in the original system (20) is restricted to \(O(\varepsilon )\).

Periodic solutions in general position of system (24) exist if \(\sin \chi =0, t \geq 0\) for certain values of the \(r_{i}\). From the 4th equation of system (24) we find the requirement:

$$ \frac{r_{2}^{2}r_{3}^{2}}{3} - \frac{r_{1}^{2}r_{3}^{2}}{2} - \frac{r _{1}^{2}r_{2}^{2}}{1} =0 . $$

Eliminating \(r_{1}\) by the \(H_{2}\) integral we find after some rearrangements the condition

$$ 2r_{2}^{2}r_{3}^{2}+ \frac{4}{3}r_{2}^{4}+ \frac{1}{6}r_{3}^{4} = \frac{1}{3}E_{0}\bigl(2r_{2}^{2}+r_{3}^{2} \bigr), \quad 0 < r_{2} < \sqrt{\frac{E _{0}}{2}},\ 0 < r_{3} < \sqrt{2E_{0}}. $$
(27)

Both for \(\chi =0\) and for \(\chi = \pi \) we find from condition (27) tori imbedded in the energy manifold. The two tori consist of periodic solutions in general position connecting the \(x_{2}\) and \(x_{3}\) normal modes. Their period is \(O(\varepsilon )\) modulated by their position on the tori. The relation between the presence of a continuous family of periodic solutions on the energy manifold and the existence of another integral (26) is an example of a more general theory on characteristic exponents of periodic solutions developed by Poincaré in [11, vol. 1].

Periodic Solutions in the Coordinate Planes

It is clear from the intermediate normal form (22) that the normalized equations of motion will contain all three normal modes. We will use co-moving coordinates \((y_{1}, y_{2}, z_{1}, z_{2}, u_{1}, u_{2})\) to normalize and to study the stability:

$$\begin{aligned} \textstyle\begin{cases} x_{1} = y_{1} \cos 3t + \frac{1}{3}y_{2} \sin 3t,\quad& \dot{x}_{1} = -3y _{1} \sin 3t + y_{2} \cos 3t, \\ x_{2} = z_{1} \cos 2t + \frac{1}{2}z_{2} \sin 2t, & \dot{x}_{2} = -2z _{1} \sin 2t + z_{2} \cos 2t, \\ x_{3} = u_{1} \cos t + u_{2} \sin t,& \dot{x}_{3} = -u_{1} \sin t + u_{2} \cos t. \end{cases}\displaystyle \end{aligned}$$
(28)

The normalized variables are obtained by averaging over time \(t\) and are satisfying the system:

$$\begin{aligned} \textstyle\begin{cases} \dot{y}_{1} = \varepsilon \frac{7}{6} d_{6} (z_{1}u_{2} + \frac{1}{2}z_{2}u_{1}), \\ \dot{y}_{2} = - \varepsilon \frac{7}{2} d_{6} (z_{1}u_{1} - \frac{1}{2}z_{2}u_{2}), \\ \dot{z}_{1} = \varepsilon \frac{7}{4} d_{6} (-y_{1}u_{2} + \frac{1}{3}y_{2}u_{1}), \\ \dot{z}_{2} = - \varepsilon \frac{7}{2} d_{6} (y_{1}u_{1} + \frac{1}{3}y_{2}u_{2}), \\ \dot{u}_{1} = \varepsilon \frac{7}{2} d_{6} (-\frac{1}{2}y_{1}z _{2} + \frac{1}{3}y_{2}z_{1}), \\ \dot{u}_{2} = - \varepsilon \frac{7}{2} d_{6} (y_{1}z_{1} + \frac{1}{6}y_{2}z_{2}). \end{cases}\displaystyle \end{aligned}$$
(29)

The generic picture for the existence of short-periodic solutions in the Hamiltonian \(1:2:3\) resonance is given in [10]. As stated above we recover three normal modes instead of generically two; this is caused by the already mentioned degenerate form of Hamiltonian (19).

The three normal modes of the normalized system are harmonic functions:

$$ A \cos mt + B \sin mt,\quad m=3,2,1,\ A^{2}+B^{2} >0. $$

To study their stability we linearize around the three normal modes of system (22) to obtain coupled Mathieu-equations; we approximate the characteristic exponents by normalizing the coupled systems. We find in these three cases after some calculations:

1. Normal mode \(x_{1}\): put \(x_{1}= A\cos 3t + B \sin 3t + w_{1}, x _{2}=w_{2},x_{3}=w_{3}\).

Transforming in the linearized system by (28) and normalization we find:

$$\begin{aligned} \dot{z}_{1} =& - \varepsilon \frac{7}{4} d_{6}(Bu_{1} - Au_{2}), \\ \dot{z}_{2} =& \varepsilon \frac{7}{2} d_{6}(Au_{1} + Bu_{2}), \\ \dot{u}_{1} =& - \varepsilon \frac{7}{2} d_{6} \biggl(Bz_{1} - \frac{1}{2}Az _{2}\biggr), \\ \dot{u}_{2} =& \varepsilon \frac{7}{2} d_{6} \biggl(Az_{1} + \frac{1}{2}Bz _{2}\biggr). \end{aligned}$$

The eigenvalues of the matrix describing this linear system have multiplicity 2 and are multiples of:

$$ \pm \sqrt{A^{2}+B^{2}}. $$

In the nomenclature of [14, Sect. 10.7.3] this is the unstable case HH.

It is interesting to consider the action-simplex with a number of initial conditions near the \(x_{1}\) normal mode, see Fig. 4. The unstable manifold of the normal mode is two-dimensional but the solutions, displayed by dots in the simplex, remain in a narrow strip extending to the edge where \(x_{1}=0\). This is caused by the third integral (26) of the normal form which tells us that the action corresponding with \(x_{1}\) is proportional to the action of \(x_{2}\).

Fig. 4
figure 4

The \(\omega =3\) normal mode (\(x_{1}\)) exists in the case 0 and is unstable (see also Fig. 6). We consider the time evolution of 98 initial positions near this normal mode by displaying the actions in the action-simplex at \(t=0, 225, 450\). The evolution is based on Hamiltonian (19); \(\varepsilon = 0.2\). The unstable manifold is two-dimensional after which the action points remain near a line in the action simplex. The inclination is explained by the expression of the third integral (26) of the normal form

Full size image

2. Normal mode \(x_{2}\): put \(x_{1}= w_{1}, x_{2}= A\cos 2t + B \sin 2t +w_{2},x_{3}=w_{3}\).

Transforming in the linearized system by (28) and normalization by averaging we find:

$$\begin{aligned} \dot{y}_{1} = - \varepsilon \frac{7}{6} d_{6}(Bu_{1} + Au_{2}), \\ \dot{y}_{2} = \varepsilon \frac{7}{2} d_{6}(Au_{1} - Bu_{2}), \\ \dot{u}_{1} = \varepsilon \frac{7}{2} d_{6} \biggl(By_{1} - \frac{1}{3}Ay _{2}\biggr), \\ \dot{u}_{2} = \varepsilon \frac{7}{2} d_{6} \biggl(Ay_{1} + \frac{1}{3}By _{2}\biggr). \end{aligned}$$

The eigenvalues have multiplicity 2 and are multiples of:

$$ \pm i \sqrt{A^{2}+B^{2}}. $$

In the nomenclature of [14] this is the spectrally stable case EE, but with both positive and negative imaginary eigenvalues coincident. A numerical calculation confirms the stability in the sense that the solutions remain near the normal mode during a finite time.

When varying \(u\), this will produce a Hamiltonian-Hopf bifurcation, see the next subsection.

As the normal mode is spectrally stable, it is of interest to display the behaviour of the actions of solutions starting near this normal mode. In Fig. 5 we show that for a limited time interval, the actions stay nearby.

Fig. 5
figure 5

The \(\omega =2\) normal mode (\(x_{2}\)) is stable in the case 0. We consider the time evolution based on Hamiltonian (19) of 98 initial positions near this normal mode by displaying the actions in the action-simplex at \(t=0, 225, 450\); \(\varepsilon = 0.2\)

Full size image

3. Normal mode \(x_{3}\): put \(x_{1}= w_{1}, x_{2}= w_{2},x_{3}=A\cos t + B \sin t +w_{3}\).

Transforming in the linearized system by (28) and normalization we find:

$$\begin{aligned} \dot{y}_{1} =& - \varepsilon \frac{7}{6} d_{6} \biggl(Bz_{1} + \frac{1}{2}Az _{2}\biggr), \\ \dot{y}_{2} =& \varepsilon \frac{7}{2} d_{6} \biggl(Az_{1} - \frac{1}{2}Bz _{2}\biggr), \\ \dot{z}_{1} =& \varepsilon \frac{7}{4} d_{6} \biggl(By_{1} - \frac{1}{3}Ay _{2}\biggr), \\ \dot{z}_{2} =& \varepsilon \frac{7}{2} d_{6} \biggl(Ay_{1} + \frac{1}{3}By _{2}\biggr). \end{aligned}$$

The eigenvalues have multiplicity 2 and are multiples of:

$$ \pm i \sqrt{A^{2}+B^{2}}. $$

In the nomenclature of [14, Sect. 10.7.3] this is the spectrally stable case EE, but again with both positive and negative imaginary eigenvalues coincident. The numerical behaviour (not shown) looks similar to Fig. 5.

Our choice of well-balanced masses involves the symmetry \(a_{1}=a_{3}\). In the sequel we will see that other choices of masses producing \(1:2:3\) resonance give qualitatively different results. It is interesting to compare the dynamics of case 0 (\(u=d_{9}=0\)) with the dynamics for \(u>0\). Such a comparison will be given in the next subsections.

4.2 The Hamiltonian-Hopf Bifurcation

In the preceding subsection we considered a rather symmetric case, \(a_{1}=a_{3}\), corresponding with \(u=0\), producing an integrable normal form; see Sect. 3.6 and Table 2. We will now consider the cases \(0< u< u_{1} (= 0.887732\ldots )\); as \(u\) increases through the interval \((0, u_{1})\) the masses will differ more and more, producing generic Hamiltonians. To put system (16) in the standard form of perturbed harmonic equations we have to apply again a symplectic transformation, i.e. (12) from Sect. 3.6. This leads to a transformed Hamiltonian (with rescaled frequencies) of the form \(H_{2}+ \varepsilon H_{3}\) with

$$ H_{2} = \frac{1}{2} \biggl(\dot{x}_{1}^{2}+ \frac{9}{14}x_{1}^{2} + \dot{x} _{2}^{2}+ \frac{4}{14}x_{2}^{2} +\dot{x}_{3}^{2}+ \frac{1}{14}x_{3} ^{2}\biggr) $$

and

$$\begin{aligned} \textstyle\begin{cases} H_{3} = d_{1}x_{1}^{3}+ d_{2}x_{2}x_{1}^{2} +d_{3}x_{3}x_{1}^{2}+d _{4}x_{2}^{2}x_{1}+d_{5}x_{3}^{2}x_{1}+d_{6}x_{1}x_{2}x_{3}+ d_{7}x _{2}^{3} + d_{8}x_{3}^{3}\\ \phantom{H_{3} =\,\,} {}+ d_{9}x_{2}x_{3}^{2} + d_{10}x_{2}^{2}x_{3}, \end{cases}\displaystyle \end{aligned}$$
(30)

with all coefficients non-zero, see Table 2. After rescaling time \(t \rightarrow t/ \sqrt{14}\), the equations of motion for the three dof system can be written as:

$$\begin{aligned} \textstyle\begin{cases} \ddot{x}_{1} + 9x_{1} = - \varepsilon 14(3d_{1}x_{1}^{2}+2d_{2}x _{1}x_{2} +2d_{3}x_{1}x_{3}+d_{4}x_{2}^{2}+ d_{5}x_{3}^{2}+ d_{6}x _{2}x_{3}), \\ \ddot{x}_{2} + 4 x_{2} = - \varepsilon 14(d_{2}x_{1}^{2}+2d_{4}x _{2}x_{1}+d_{6}x_{1}x_{3}+3d_{7}x_{2}^{2}+ d_{9}x_{3}^{2}+ 2d_{10}x _{2}x_{3}), \\ \ddot{x}_{3} + x_{3} = - \varepsilon 14(d_{3}x_{1}^{2}+2d_{5}x_{3}x _{1}+d_{6}x_{1}x_{2} +3d_{8}x_{3}^{2}+ 2d_{9}x_{2}x_{3} + d_{10}x_{2} ^{2}). \end{cases}\displaystyle \end{aligned}$$
(31)

The size of the coefficients of \(H_{3}\) are comparable with the size of \(d_{6}\) or smaller, we will give them explicitly as examples for the cases 1 and 2 in Sect. 4.3 with less balanced masses.

In the cubic part of the normalized Hamiltonian we retain of the cubic part only the terms with \(d_{6}\) and \(d_{9}\); the other terms are, after normalization, active only at higher order. So, anticipating this, an intermediate normal form of the equations of motion becomes:

$$\begin{aligned} \textstyle\begin{cases} \ddot{x}_{1} + 9 x_{1} = - \varepsilon 14d_{6}x_{2}x_{3}, \\ \ddot{x}_{2} + 4 x_{2} = - \varepsilon 14(d_{6}x_{1}x_{3}+d_{9}x _{3}^{2}), \\ \ddot{x}_{3} + x_{3} = - \varepsilon 14(d_{6}x_{1}x_{2} + 2d_{9}x _{2}x_{3}), \end{cases}\displaystyle \end{aligned}$$
(32)

The Normal Form and Periodic Solutions Outside the Coordinate Planes

Using transformation (23) and putting \(\phi_{1}- \phi_{2} - \phi _{3} = \chi_{1}\), \(2 \phi_{3} - \phi_{2} = \chi_{2}\), we find after averaging-normalization:

$$\begin{aligned} \textstyle\begin{cases} \dot{r}_{1} = \varepsilon \frac{7}{6}d_{6}r_{2}r_{3} \sin \chi_{1}, \\ \dot{r}_{2} = - \varepsilon \frac{7}{4} ( d_{6}r_{1}r_{3} \sin \chi_{1} +d_{9} r_{3}^{2} \sin \chi_{2} ) , \\ \dot{r}_{3} = - \varepsilon \frac{7}{2} ( d_{6}r_{1}r_{2} \sin \chi_{1} - 2d_{9}r_{2}r_{3} \sin \chi_{2} ) , \\ \dot{\chi _{1}} = \varepsilon \frac{7}{2} [ d_{6} \frac{\cos \chi_{1}}{r_{1}r_{2}r_{3}} ( \frac{r_{2}^{2}r_{3}^{2}}{3} - \frac{r _{1}^{2}r_{3}^{2}}{2} - \frac{r_{1}^{2}r_{2}^{2}}{1} ) - d_{9} \frac{ \cos \chi_{2}}{r_{2}} ( \frac{1}{2}r_{3}^{2}+2r_{2}^{2} ) ] , \\ \dot{\chi _{2}} = \varepsilon \frac{7}{4} ( d_{6} \frac{r_{1} \cos \chi_{1}}{r_{2}r_{3}} (4r_{2}^{2}-r_{3}^{2}) + d_{9} \frac{ \cos \chi_{2}}{r_{2}} (8r_{2}^{2} - r_{3}^{2}) ) . \end{cases}\displaystyle \end{aligned}$$
(33)

The integral \(H_{2}\) of the normal form equations becomes again:

$$ 9r_{1}^{2}+4r_{2}^{2}+r_{3}^{2} = 2E_{0}, $$
(34)

Periodic solutions in general position with constant amplitude have to satisfy \(\sin \chi_{1} = \sin \chi_{2} =0\) or \(\chi_{1}= 0, \pi \) and \(\chi_{2}= 0, \pi \). We have

$$ \cos \chi_{1} \cos \chi_{2} = \pm 1, \quad \quad q = \frac{d_{9}}{d_{6}} >0. $$

From the last two equations of system (33) we have the conditions:

$$\begin{aligned} & \frac{r_{2}^{2}r_{3}^{2}}{3} - \frac{r_{1}^{2}r_{3}^{2}}{2} - \frac{r _{1}^{2}r_{2}^{2}}{1}= \pm qr_{1}r_{3} \biggl( \frac{1}{2}r_{3}^{2}+2r _{2}^{2} \biggr) , \end{aligned}$$
(35)
$$\begin{aligned} & 4r_{2}^{2}-r_{3}^{2} = \pm q \frac{r_{3}}{r_{1}} \bigl(r_{3}^{2}-8r_{2} ^{2}\bigr). \end{aligned}$$
(36)

Eliminating \(r_{1}\) from (35) using (36) we obtain two equations that are quadratic in \(r_{2}^{2}\) and \(r_{3}^{2}\). Eliminating \(r_{1}\) from the \(H_{2}\) integral we find one equation that is quadratic in \(r_{2}^{2}\) and \(r_{3}^{2}\). These expressions have to be handled for the range of \(q\) determined by \(u \in (0, u_{1})\). Using Mathematica and corresponding plots we find four positive solutions corresponding with four periodic solutions characterized by two different phases.

We omit the stability analysis, but note that the generic case of the \(1:2:3\) resonance was studied in [10] that produces four general position periodic solutions with the stability types \(EE\) and \(EH\).

Periodic Solutions in the Coordinate Planes

Inspection of the intermediate normal form system (32) shows that the \(x_{1}\) and \(x_{2}\) normal modes exist as solutions of this system, the \(x_{3}\) normal mode does not. It is shown in [10] that the normal mode \(x_{2}\) is unstable. If the instability is of class \(C\) (complex eigenvalues), a Shilnikov-Devaney bifurcation [6] may take place resulting in chaotic dynamics originating from a neighbourhood of the complex unstable normal mode. To avoid singularities near the normal modes we use again the co-moving variables from transformation (28). The normalized variables satisfy the system:

$$\begin{aligned} \textstyle\begin{cases} \dot{y}_{1} = \varepsilon \frac{7}{6} d_{6} (z_{1}u_{2} + \frac{1}{2}z_{2}u_{1}), \\ \dot{y}_{2} = - \varepsilon \frac{7}{2} d_{6} (z_{1}u_{1} - \frac{1}{2}z_{2}u_{2}), \\ \dot{z}_{1} = \varepsilon \frac{7}{2} [ \frac{1}{2}d_{6} (-y_{1}u _{2} + \frac{1}{3}y_{2}u_{1})+d_{9}u_{1}u_{2}], \\ \dot{z}_{2} = - \varepsilon \frac{7}{2} [d_{6} (y_{1}u_{1} + \frac{1}{3}y_{2}u_{2})+d_{9}(u_{1}^{2}-u_{2}^{2})], \\ \dot{u}_{1} = \varepsilon \frac{7}{2}[ d_{6} (-\frac{1}{2}y_{1}z _{2} + \frac{1}{3}y_{2}z_{1})+ d_{9}(-2z_{1}u_{2}+z_{2}u_{1})], \\ \dot{u}_{2} = - \varepsilon \frac{7}{2}[ d_{6} (y_{1}z_{1} + \frac{1}{6}y_{2}z_{2})+d_{9}(2z_{1}u_{1}+z_{2}u_{2})]. \end{cases}\displaystyle \end{aligned}$$
(37)

We find three families of short-periodic solutions; the constants \(A, B\) are real, \(A^{2}+B^{2}>0\).

  1. 1.

    \(x_{1}(t)=A \cos 3t + B \sin 3t, x_{2}=x_{3}=0\).

  2. 2.

    \(x_{2}(t)= A \cos 2t + B \sin 2t, x_{1}=x_{3}=0\).

  3. 3.

    If \(x_{2}(t)=0, d_{6} \neq 0\):

    $$\begin{aligned} \textstyle\begin{cases} x_{1}(t) = \frac{d_{9}}{d_{6}} ( \frac{A}{A^{2}+B^{2}}(3B^{2}-A ^{2}) \cos 3t - \frac{B}{A^{2}+B^{2}}(3A^{2}-B^{2}) \sin 3t ) , \\ x_{3}(t) = A \cos t + B \sin t. \end{cases}\displaystyle \end{aligned}$$
    (38)

    If \(d_{9}\) differs from zero, this family of periodic solutions moves along the \(x_{2}=0\) edge of the simplex in Fig. 6 starting from the \(x_{3}\) normal mode that exists if \(d_{9}=0\).

    Fig. 6
    figure 6

    The action simplices of the cases \(u=0\) and \(0< u< u_{1}\); the cases 1 and 2 are typical for the family of Hamiltonians where \(u \neq 0\). The actions \(\tau_{i}\) (related to \(r_{i}^{2}\)) form a triangle for fixed values of \(H_{2}\) which is an integral of the normal forms. The frequencies have been normalized to \(1, 2, 3\) to indicate the \(x_{3}, x_{2}, x_{1}\) normal mode positions at the vertices. The black dots indicate periodic solutions, the indicated stability types are HH (hyperbolic-hyperbolic), EE (elliptic-elliptic) and C (complex). The two (roughly sketched) curves connecting the \(x_{2}\) and \(x_{3}\) normal modes in the left simplex correspond with two tori consisting of periodic solutions, respectively with combination angles \(\chi =0\) and \(\pi \). The tori break up into 4 general position periodic solutions if \(u>0\) (cases 1 and 2)

    Full size image

To evaluate the stability of the periodic solutions we will linearize system (32) near these solutions; this produces coupled Mathieu equations which we will analyze by normalization.

The \(x_{2}\) Normal Mode

Put:

$$ x_{1}=w_{1}, \quad \quad x_{2} = A \cos 2t + B \sin 2t +w_{2}, \quad \quad x_{3}=w_{3}, $$

with real constants \(A, B, A^{2}+B^{2} >0\) and corresponding expressions for the derivatives. We find after linearization

$$\begin{aligned} \textstyle\begin{cases} \ddot{w}_{1} + 9 w_{1} = - \varepsilon 14 d_{6} (A \cos 2t + B \sin 2t)w_{3}, \\ \ddot{w}_{2} + 4 w_{2} = 0, \\ \ddot{w}_{3} + w_{3} = - \varepsilon 14 [d_{6}w_{1}(A \cos 2t + B \sin 2t) + 2d_{9} (A \cos 2t + B \sin 2t)w_{3}], \end{cases}\displaystyle \end{aligned}$$
(39)

We study the stability of this system by normalization to find the eigenvalues of the matrix (omitting the factor \(7 \varepsilon /2\))

$$ \left ( \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c} 0 & 0 & \frac{d_{6}}{3}B & \frac{d_{6}}{3}A \\ 0 & 0 & -d_{6}A & d_{6}B \\ -d_{6}B & \frac{d_{6}}{3}A & 2d_{9}B & - 2d_{9}A \\ -d_{6}A & -\frac{d_{6}}{3}B & - 2d_{9}A & -2 d_{9}B \end{array}\displaystyle \right ) $$

produce first order approximations of the characteristic exponents of system (39). For the eigenvalues we find apart from the factor \(7 \varepsilon /2\):

$$ \lambda^{2} = -\bigl(A^{2}+B^{2}\bigr) \biggl[ \biggl(\frac{1}{3}d_{6}^{2}-2d_{9}^{2} \biggr) \pm 2d_{9} \sqrt{d_{9}^{2}- \frac{1}{3}d_{6}^{2}} \biggr] . $$

A sufficient condition for the complex case \(C\) to arise is

$$ d_{6}^{2} > 6d_{9}^{2}. $$
(40)

This condition corresponds with the condition in Table 1 of [10]. Condition (40) is satisfied for \(0< u < u_{1}\) so that the complex case \(C\) arises for \(u>0\).

Another view of the eigenvalues is obtained by realizing that in Sect. 4.1 we had \(u=0\) resulting in \(d_{6} \neq 0, d_{9}=0\); \(u=0\) gives for the \(x_{2}\) normal mode purely imaginary eigenvalues with multiplicity two. As \(u\) increases (\(d_{9} \neq 0\)), the eigenvalues move from the imaginary axis into the complex domain. This is part of the Hamiltonian-Hopf bifurcation, see Fig. 7.

Fig. 7
figure 7

The Hamiltonian-Hopf bifurcation of a periodic solution in a three dof system as takes place for the \(x_{2}\) normal mode in \(1:2:3\) resonance of [10]. In our problem we have only the transition from the case of double imaginary eigenvalues to complex eigenvalues as the double eigenvalues are generated by the symmetry at the start of the interval \(0 < u < u_{1}\)

Full size image

For case 2 (see Sect. 4.3) we show in the action-simplex of Fig. 8 the behaviour of solutions starting near this complex unstable normal mode.

Fig. 8
figure 8

The \(\omega =2\) normal mode (\(x_{2}\)) exists in the case 2 and is complex unstable (see also Fig. 6). We consider the time evolution of 98 initial positions near this normal mode by displaying the actions in the action-simplex at \(t=0, 225, 450\); \(\varepsilon = 0.2\)

Full size image

The \(x_{1}\) Normal Mode

For \(d_{9}=0\) we have found in the preceding subsection the case HH. This is a generic case of eigenvalues, so for \(d_{9}\) small enough and \(u \in (0, u_{1})\) the nature of the instability will not change but the dynamics is very different as the normal form is not integrable.

For case 2 (see Sect. 4.3) we show in the action-simplex of Fig. 9 the behaviour of solutions starting near this unstable normal mode.

Fig. 9
figure 9

The \(\omega =3\) normal mode (\(x_{1}\)) exists in the case 2 and is unstable (see also Fig. 6). We consider the time evolution of 98 initial positions near this normal mode by displaying the actions in the action-simplex at \(t=0, 225, 450\), \(\varepsilon = 0.2\). The behaviour is different from the case 0, see Fig. 4, as in this case the normal form is not integrable

Full size image

The Periodic Solution for \(x_{2}(t)=0\)

For the periodic solution (38) we put:

$$ x_{1}= C \cos 3t + D \sin 3t, \quad \quad x_{3} = A \cos t + B \sin t. $$

Transforming

$$ x_{1}= C \cos 3t + D \sin 3t + w_{1},\quad \quad x_{2}=w_{2},\quad \quad x_{3}= A \cos t + B \sin t + w_{3}, $$

and substitution into system (32), we find after linearization:

$$\begin{aligned} \textstyle\begin{cases} \ddot{w}_{1} + 9 w_{1} = - \varepsilon 14d_{6} (A \cos t + B \sin t)w _{2}, \\ \ddot{w}_{2} + 4 w_{2} = - \varepsilon 14[d_{6}(C \cos 3t +D \sin 3t)w_{1} + d_{6} ( A \cos t + B \sin t)w_{1}\\ \phantom{\ddot{w}_{2} + 4 w_{2} =\,\,}{} +2d_{9}( A \cos t + B \sin t)w _{3}], \\ \ddot{w}_{3} + w_{3} = - \varepsilon 14[d_{6}(C \cos 3t +D \sin 3t)w _{2} +2d_{9} (A \cos t + B \sin t)w_{2}]. \end{cases}\displaystyle \end{aligned}$$
(41)

To investigate stability we normalize near the periodic solution; apart from a factor \(7 \varepsilon /2\), this produces the matrix:

$$ \left ( \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} 0 & 0 & \frac{d_{6}}{3}B & \frac{d_{6}}{6}A & 0 & 0 \\ 0 & 0 & -d_{6}A & \frac{d_{6}}{2}B & 0 & 0 \\ -\frac{d_{6}}{2}B & \frac{d_{6}}{6}A & 0 & 0 & d_{6} \frac{D}{2} +d _{9}B & - d_{6} \frac{C}{2}+d_{9}A \\ -d_{6}A & - \frac{d_{6}}{3}B & 0 & 0 & -d_{6}C- 2d_{9}A & -d_{6}D+2 d _{9}B \\ 0 & 0 & d_{6}D-2d_{9}B & - d_{6} \frac{C}{2}+d_{9}A & 0 & 0 \\ 0 & 0 & -d_{6}C - 2d_{9}A & -d_{6} \frac{D}{2}-d_{9}B & 0 & 0 \end{array}\displaystyle \right ) . $$

Using the values of \(C\) and \(D\) given in (38), we find purely imaginary eigenvalues with multiplicity two. The results have been summarized in Fig. 6.

4.3 Experiments for Two Cases with \(u>0\)

The first order normal form of system (31) produces qualitatively the same dynamics for \(u \in (0, u_{1})\). We consider a few experiments for two cases that are typical for this dynamics.

Case 1 with Less-Balanced Masses

We choose for \(u = 0.534105\) from Eqs. (9)):

$$ a_{1}= 0.00510292, \quad \quad a_{2}= 0.117265,\quad \quad a_{3}= 0.0854008,\quad \quad a_{4}= 0.292231 $$

In this case we have \(m_{1} > m_{3} > m_{2}> m_{4}\). With these mass (\(a _{i}\)) values the symplectic transformation of Sect. 3.6 to system (31) produces the expression:

$$\begin{aligned} H_{3} =& 0.0281999 x_{1}^{3} -0.0258437 x_{1}^{2}x_{2} -0.0777574x _{1}^{2}x_{3}- 0.0275058x_{1}x_{2}^{2}\\ &{}- 0.00252349x_{1}x_{3}^{2} -0.0306229x_{1}x_{2}x_{3}+ 0.0157538x_{2}^{3} + 0.000502655x_{3} ^{3}\\ &{}- 0.0089438x_{2}x_{3}^{2} + 0.028527x_{2}^{2}x_{3}. \end{aligned}$$

We have the case:

$$ d_{6}= -0.0306229, \quad \quad d_{9}= -0.0089438. $$

so that the \(x_{2}\) normal mode is complex unstable; see Fig. 6. The \(H_{2}(t)\) time series are shown in Fig. 10.

Fig. 10
figure 10

The \(H_{2}(t)\) time series (sum of the actions) based on system(32) of case 1 for two sets of initial values. (Left) Starting near the \(x_{1}\) normal mode plane with \(x_{1}=1, x_{2}=0.1, x_{3}= 0.1, \dot{x}_{i}=0, i=1, 2, 3\); \(\varepsilon = 0.5, H_{2}(0) \approx 4.52\). The variations of \(H_{2}\) start with 0.4 and are near \(t =1000\) 0.15. Energy is pumped into \(H_{3}\). (Right) The \(H_{2}(t)\) time series starting near the complex unstable normal mode \(x_{2}\) with \(x_{1}=0.1, x_{2}=1.5, x_{3}= 0.1, \dot{x}_{i}=0, i=1, 2, 3\); \(\varepsilon = 0.5, H_{2}(0) \approx 4.52\). The variations of \(H_{2}\) are 0.4. Horizontal scales: time in \([0,1000]\), vertical scales: \(H_{2}\) in \([4.25,4.75]\)

Full size image

Note that \(d_{9}\) is still fairly small with the implication that the expansion of the flow near the \(x_{2}\) normal mode will not be very explosive. This may reduce the amount of chaos present in the system. We compare with case 0 and give a few more details for different initial conditions based on integration of system (22) and system (32). We established that in all cases the \(x_{1}\) normal mode is unstable (\(HH\)), see also Fig. 6. Starting near the \(x_{1}\) normal mode in case 0, the solutions move away, guided by the two-dimensional unstable manifold of the normal mode; the integrability of the normal form produces a fairly regular \(H_{2}(t)\) with variations 0.05. In Fig. 10 we display \(H_{2}(t)\) for case 1 with the same initial conditions showing strong variations of \(H_{2}(t)\); its behaviour is influenced by the chaotic character of the normal form. On this interval of time \([0, 1000]\), energy is clearly pumped into \(H_{3}\) but the recurrence of the Hamiltonian system will return this on a much longer timescale.

The chaos in case 1 (and 2) is strongly influenced by the complex instability of the \(x_{2}\) normal mode. In case 0 this mode is stable so that \(H_{2}(t)\) varies again with 0.05 only. Using the same initial conditions for case 1 we find strong variations (0.4) of \(H_{2}(t)\), but always within the limits of the error estimates; see Fig. 10, right.

Case 2 with Less-Balanced Masses

We choose for \(u = 0.826713\) from Eqs. (9) a case with even less balanced masses; in this case \(m_{1}\) is quite massive. We have:

$$ a_{1}= 0.000685158,\quad \quad a_{2}= 0.11239,\quad \quad a_{3}= 0.100269,\quad \quad a_{4}= 0.286656. $$

With these mass (\(a_{i}\)) values the symplectic transformation of Sect. 3.6 to system (31) produces the expression:

$$\begin{aligned} H_{3} =& 0.0352657 x_{1}^{3} -0.0272316 x_{1}^{2}x_{2} -0.0743155x _{1}^{2}x_{3}- 0.0366184x_{1}x_{2}^{2}\\ &{}- 0.00260064x_{1}x_{3}^{2} -0.0337877x_{1}x_{2}x_{3}+ 0.0181144x_{2}^{3} + 0.000760425x_{3} ^{3}\\ &{}- 0.0105601x_{2}x_{3}^{2} + 0.023904x_{2}^{2}x_{3}. \end{aligned}$$

We have the case:

$$ d_{6}= -0.0337877,\quad \quad d_{9}= -0.0105601 $$

If \(d_{9} \neq 0\) (the cases 1 and 2), the \(x_{3}\) normal mode does not exist. In Fig. 11 we show the action-simplex for solutions starting near the \(x_{1}=x_{2}=0\) position, so near the \(\omega =1\) vertex.

Fig. 11
figure 11

We consider for case 2 the time evolution of 98 starting points near the \(\omega =1\) vertex by displaying the action-simplex at various times

Full size image

We computed the \(H_{2}(t)\) time series of case 2 based on system (32); we omit the picture as it is similar to Fig. 10, right.

4.4 Comparison with Another Hamiltonian System in \(1:2:3\) Resonance

It is instructive to discuss our results for the inhomogeneous FPU chain with another Hamiltonian system in \(1:2:3\) resonance, and compare the instability types of the \(x_{2}\) normal mode. In this system we have a full Hamiltonian-Hopf bifurcation as sketched in Fig. 7.

For the inhomogeneous FPU lattice in \(1:2:3\) resonance we found complex instability (\(C\)) of the \(x_{2}\) normal mode and no cases of \(HH\) instability. Both cases, \(HH\) and \(C\) lead to a non-integrable normal form but the dynamics is different. See [5].

To illustrate the different dynamics consider the Hamiltonian presented as an example in [18]:

$$ H(p, q)= \frac{1}{2}\bigl(p_{1}^{2}+q_{1}^{2} \bigr)+\bigl(p_{2}^{2}+q_{2}^{2}\bigr)+ \frac{3}{2}\bigl(p_{3}^{2}+q_{3}^{2} \bigr)- \varepsilon q_{1}^{2}(a_{2}q_{2}+a _{3}q_{3})- \varepsilon b q_{1}q_{2}q_{3}. $$
(42)

This system is in \(1:2:3\) resonance but it is not derived from a FPU chain. We present \(H_{2}(t)\) for both instability cases in Fig. 12. The dynamics is chaotic but in the case left, the \(q_{2}\) normal mode is unstable with real eigenvalues (HH); transverse homoclinic intersections produce chaotic motion. On the right the \(q_{2}\) normal mode is complex unstable (C) which produces the Hamiltonian Devaney-Shilnikov phenomenon. This involves a homoclinic orbit surrounded by an infinite number of unstable periodic solutions producing more violent chaotic motion as predicted in [6].

Fig. 12
figure 12

Two \(H_{2}(t)\) time series based on Hamiltonian (42). On the left \(x_{1}(0)=0.1\), \(x_{2}(0)=1\), \(x_{3}(0)=0.5\), and on the right \(x_{1}(0)=2\), \(x_{2}(0)=1 x_{3}(0)=-0.05\). For both time series we use \(\varepsilon = 0.5\), \(a_{2}=3\), \(a _{3}=1\), \(b=1\), and \(\dot{x}_{1}(0)=\dot{x}_{2}(0)=\dot{x}_{3}(0)=0\). For the \(x_{2}\) normal mode we have instability \(HH\) on the left and instability \(C\) on the right. In both cases the Hamiltonian flow is chaotic but in the right picture the system has undergone Devaney-Shilnikov bifurcation. Horizontal scales: time in \([0,500]\), vertical scales: energy in \([1.15,1.65]\) (left) and in \([3,4.5]\) (right)

Full size image

5 Conclusions

General

  • For an inhomogeneous periodic FPU-chain with four particles, most frequency ratios occur for a one-dimensional variety of mass ratios. The frequency ratios \(1:2:1\) and \(1:1:3\) arise for a finite number of mass ratios, the ratios \(1:2:2\), \(1:1:1\) and \(1:3:3\) do not occur at all in this FPU-chain. See Table 1.

  • For any number of particles \(n\geq 3\) the set of mass distributions for a given frequency distribution has a relatively simple algebraic structure. For \(n=4\) we describe algorithmically how to determine this set for a given frequency distribution. For \(n\geq 4\) there are frequency distributions that do not correspond to any mass distribution.

  • Applications.

    There are many applications of chains of particles. A large number focuses on molecular dynamics; in fact normalization was introduced and used in one of the classics on atomic dynamics [1]. A recent result is to construct a chain of FPU cells where each cell is a 4-particles FPU system, weakly connected with identical cells; see [19]. Another application based on normal forms and related symmetry considerations is to consider a chain of \(2n\) particles with alternating masses; see [9] and [3].

The Case of Four Particles in \(1:2:3\) Resonance

  • The \(\alpha \)-chain with four particles in \(1:2:3\) resonance does not contain the case of the second normal mode with \(HH\) instability. For nearly all mass ratios it contains the second normal mode with complex \(C\) instability showing the generic features of the general \(1:2:3\) resonance described in [10].

    In this sense it is very different form the classical FPU system with equal masses where many resonances arise in the spectrum that are not effective, see [13].

  • A special case of the resonance \(1:2:3\) has the symmetry of two opposite equal masses and two quite different masses. Along the variety of mass ratios as a limit case one of the masses tends to infinity.

  • The symmetric case of two equal masses differs dynamically from the other cases. The transition to four different masses corresponds to half of a Hamiltonian-Hopf bifurcation with a Shilnikov-Devaney bifurcation producing chaotic dynamics. In a more general context such behaviour of the \(1:2:3\) resonance was described in [10].

  • The normalized system for the symmetric case of two equal masses is integrable and has periodic solutions for each of the three eigenmodes (the normal modes). Moreover, there are on the energy manifold two families of periodic solutions connecting the second and the third eigenmode. This is a degeneration in the sense described by Poincaré [11, vol. 1]. The symmetry in this special case makes the \(1:2:3\) resonance more similar to the classical FPU problem.

  • Under the transition away from the symmetric case, the eigenmodes \(x_{1}\) (associated with frequency 3) and \(x_{2}\) (associated with frequency 2) produce a periodic solution (normal mode) in the nonlinear system. The periodic solution that was associated to the third eigenmode in the symmetric case moves away along an edge of the action simplex. The two continuous families of periodic solutions of the symmetric case break up into four periodic solutions.

  • To summarize: the inhomogeneous periodic FPU \(\alpha \)-chain with four particles is characterized by a non-integrable normal form, except in the symmetric case of two equal masses. The implication is that near stable equilibrium its chaotic behaviour is not restricted to exponentially small sets as in the case of two dof systems and as in the case of the classical FPU \(\alpha \)-chain.

Considering Again the Classical FPU Chain

  • We have shown that in the case of four particles the presence of two equal masses produces a symmetry in the dynamical system that makes the system structurally unstable. In this perspective the model of the classical FPU chain with all masses equal is also structurally unstable and misleading as a model.

  • The averaging-normal form technique we have used is valid for an arbitrary number of particles as long as the total energy of the chain is finite and small. This enables us to extend the analysis to chains with many particles as was shown in [13].

  • The dynamics on the energy manifold is structured by approximate invariant manifolds, some of them valid for all time, some with finite but long validity (\(1/ \varepsilon^{m}\) intervals for some positive \(m\)). At the same time the Poincaré recurrence theorem produces relatively short recurrence times, see [19]. Altogether this suggests that the classical FPU chain for low energy values does not lead to equipartition of energy and is not a good model for statistical mechanics.