The Inhomogeneous Fermi-Pasta-Ulam Chain, a Case Study of the 1 : 2 : 3 Resonance

A 4-particles chain with different masses represents a natural generalization of the classical Fermi-Pasta-Ulam chain. It is studied by identifying the mass ratios that produce prominent resonances. This is a technically complicated problem as we have to solve an inverse problem for the spectrum of the corresponding linearized equations of motion. In the case of such an inhomogeneous periodic chain with four particles each mass ratio determines a frequency ratio for the quadratic part of the Hamiltonian. Most prominent frequency ratios occur but not all. In general we ﬁnd a one-dimensional variety of mass ratios for a given frequency ratio. A detailed study is presented of the resonance 1 : 2 : 3. A small cubic term added to the Hamiltonian leads to a dynamical behaviour that shows a difference between the case that two opposite masses are equal and a striking difference with the classical case of four equal masses. For two equal masses and two different ones the normalized system is integrable and chaotic behaviour is small-scale. In the transition to four different masses we ﬁnd a Hamiltonian-Hopf bifurcation of one of the normal modes leading to complex instability and Shilnikov-Devaney bifurcation. The other families of short-periodic solutions can be localized from the normal forms together with their stability characteristics. For illustration we use action simplices and examples of behaviour with time.

. 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 α 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 R n . This subset is empty in some cases, for instance the important 1 : 1 : . . . : 1 resonance does not arise for the periodic FPU α chain with four or more particles.

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 amplitudephase 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 α-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 comoving 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.

Outline of a Research Programme
The original Fermi-Pasta-Ulam chain [7] consists of n oscillators of equal mass with nearestneighbour 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. 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]. 2. Cases of five and more particles will present many more problems. 3. The present study is restricted to so-called periodic α-chains. Including quartic terms in the Hamiltonian (β-chains) and considering lattices with fixed begin-and end-point will produce new results. 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.

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 jqj , j = 1 . . . n, ε ≥ 0 a small parameter, the Hamiltonian is of the form: The quadratic part of the Hamiltonian is not in diagonal form; for n = 3, 4 . . . the linearized equations of motion can be written as: We can write for the quadratic part of H (p, q): with A n the n × n diagonal matrix with at position (i, i) the value m −1 i =: a i , C n is an n × 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.

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.

The Matrix for Inhomogeneous FPU-Lattices and Its Eigenvalues
The linear system (2) can be written as where the matrix A n is a diagonal matrix with the inverse masses m −1 j =: 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, (This matrix turns up elsewhere in mathematics. It is the affine Cartan matrix of the completed root systemĀ n . See e.g.
In the sequel we will choose the case of vanishing momentum integral which is not a restriction of generality. If λ is a positive eigenvalue of A n C n , then i √ λ and −i √ λ 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 . If y is an eigenvector of A n C n with eigenvalue λ, then λA −1 n y = C n y, and the scalar product with y for both sides of the equation gives the equality (Indices taken modulo n.) So With Schwarz's inequality applied to the vector (y i ) and the vector (y i−1 ) this implies λ ≥ 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, . . . , 1).
For the investigation of the linearized problem we need to understand the map R n >0 → R n−1 >0 , from a vector (a 1 , . . . , a n ) of inverse masses to a vector (λ 1 , . . . , λ n−1 ) of positive eigenvalues. The order of the eigenvalues is not determined, so we have, more precisely, a map ρ n : R n >0 → S n−1 \R n−1 >0 , 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: ρ(t (a 1 , . . . , a n )) = tρ(a 1 , . . . , a n ) for t > 0.
To investigate the correspondence between eigenvalues and inverse masses we use the equality det(A n C n − λI n ) = −λ (λ j − λ), for (λ 1 , . . . , λ n−1 ) = ρ(a 1 , . . . , a n ). Comparing the coefficients of this polynomial we obtain the equalities 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 C n which is the intersection of n − 1 hyperplanes given by equations of degree 1, 2, . . . , n − 1.
In Lemma A.2 in the Appendix we'll show that with c i,j = 3 if i − j = ±1 mod 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 , . . . , a n ) for a given eigenvalue vector (λ 1 , . . . , λ n−1 ) are in a compact subset of R n . This subset may be empty. For all n ≥ 4 the 1 : 1 : · · · : 1 resonance does not occur for any mass distribution. For the eigenvalue ratio λ 1 /λ 2 = 2 of matrix A n C n the solution set is compact; it is an ellipse in R 3 >0 . For the eigenvalue ratio λ 1 /λ 2 = 4 the solutions are on a larger ellipse in R 3 , which intersects R 3 >0 in three open curves

The Case of Three Particles
For n = 3 the determination of the eigenvalues for given inverse masses amounts to solving the quadratic equation λ 2 − 2(a 1 + a 2 + a 3 )λ + 3(a 1 a 2 + a 1 a 3 + a 2 a 3 ) = 0, which has positive solutions.
Conversely, for all choices (λ 1 , λ 2 ) of positive eigenvalues, values of a 1 , a 2 , a 3 can be found such that A 3 C 3 has eigenvalues λ 1 , λ 2 and 0. If λ 1 = λ 2 there is exactly one solution a 1 = a 2 = a 3 = 1 3 λ (equal masses). If the eigenvalues have ratio λ 1 /λ 2 > 1 then the corresponding points (a 1 , a 2 , a 3 ) in R 3 form an ellipse. This ellipse may or may not be contained in the positive octant. See Fig. 1.
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 2 1 : n 2 2 : n 2 3 . 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 ).

The Resonance 1 : 2 : 3
Here we consider the resonance that is the subject of study in the next section.
By scaling we arrange λ 1 = 9 14 , λ 2 = 2 7 , λ 3 = 1 14 to satisfy the last equation in (8). We follow the computational scheme in Table 3 in the Appendix. This leads to where we take the minus sign for a 1 and a 2 . The parameter u runs through the interval [0, u 1 ) with  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)

Illustration of the Fiber
The equations in (9) describe a curve u → (a 1 (u), a 2 (u), a 3 (u), a 4 (u)) in R 4 >0 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 = 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 [−π, π] × [− 1 2 π, 1 2 π], with boundary identifications. The image of the fiber under this map is given in Fig. 3.

Transformation of the Hamiltonian
We form the diagonal matrix A 4 (u) with diagonal elements a j (u), 1 ≤ j ≤ 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 ∈ [0, u 1 )), so we can find an orthogonal matrix U(u) such that A 4 (u) 1 where Λ is the diagonal matrix with diagonal elements 9 14 , 2 7 , 1 14 , and 0.
Then the transformation matrices determine a symplectic transformation which transforms the quadratic part in (3) of the Hamiltonian into 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 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 ∈ [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) The substitution (q 1 , . . . , with the functions d j as indicated in Table 2.

The 1 : 2 : 3-Resonance for the Periodic α-Lattice (n = 4)
For any possible inhomogeneous FPU α-chain with four dof we have the system: The coefficient α has been retained for reference to the literature; here we will take α = 1. If a 1 = · · · = 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 ε 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: 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 , . . . , a 4 producing this resonance. All vectors (a 1 , . . . , a 4 ) are obtained by the action of the dihedral group D 4 on the coordinates and the scaling (a 1 , . . . , a 4 ) → (ta 1 , . . . , 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 ≤ 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 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 we obtain an O(ε) 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 formH 3 (p, q), written in action-angle coordinates or amplitude-phase coordinates (see below), will contain certain combination angles corresponding with the resonance. IfH 3 contains only one combination angle, we have an additional integral of motion and the normal form H 2 +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 τ 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.

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): 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 with .
Rescaling time t/ √ 14 → t , the equations of motion for the three dof system become: 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: 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 φ 1 − φ 2 − φ 3 and 2φ 3 − φ 2 . At second order the combination angle φ 1 − 3φ 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 + εH 3 ) of system (20) to O(ε) as in [10] or [14] and as we shall explicitly show below, only the d 6 term survives inH 3 ; this makes the Hamiltonian (19) non-generic. An intermediate normal form of the equations of motion becomes: 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 ,ẋ i → r i , ψ i of the form: The actions τ i are related to the r 2 i , the angles φ i to the arguments (ω i t + ψ i ). Putting χ = ψ 1 − ψ 2 − ψ 3 and averaging over time t , the averaging-normal form equations outside the coordinate planes become: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ṙ 1 = ε 7 6 d 6 r 2 r 3 sin χ, r 2 = −ε 7 4 d 6 r 1 r 3 sin χ, r 3 = −ε 7 2 d 6 r 1 r 2 sin χ, The integral H 2 of the normal form equations becomes: with E 0 a positive (energy) constant. The combination angle 2φ 3 − φ 2 is missing; another integral of the normal form (24) is: In the original variables this integral is: 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)   Eliminating r 1 by the H 2 integral we find after some rearrangements the condition Both for χ = 0 and for χ = π we find from condition (27)
The normalized variables are obtained by averaging over time t and are satisfying the system: 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: 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: Transforming in the linearized system by (28) and normalization we find: The eigenvalues of the matrix describing this linear system have multiplicity 2 and are multiples of: 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 .
2. Normal mode x 2 : put x 1 = w 1 , x 2 = A cos 2t + B sin 2t + w 2 , x 3 = w 3 . Fig. 4 The ω = 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); ε = 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 Transforming in the linearized system by (28) and normalization by averaging we find: The eigenvalues have multiplicity 2 and are multiples of: 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.
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: The eigenvalues have multiplicity 2 and are multiples of: 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.

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 . . .); 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)  and with all coefficients non-zero, see Table 2. After rescaling time t → t/ √ 14, the equations of motion for the three dof system can be written as: 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: The Normal Form and Periodic Solutions Outside the Coordinate Planes Using transformation (23) and putting φ 1 − φ 2 − φ 3 = χ 1 , 2φ 3 − φ 2 = χ 2 , we find after averagingnormalization: 4 (d 6 r 1 r 3 sin χ 1 + d 9 r 2 3 sin χ 2 ), r 3 = −ε 7 2 (d 6 r 1 r 2 sin χ 1 − 2d 9 r 2 r 3 sin χ 2 ), cos χ 1 r 1 r 2 r 3 ( The integral H 2 of the normal form equations becomes again: Periodic solutions in general position with constant amplitude have to satisfy sin χ 1 = sin χ 2 = 0 or χ 1 = 0, π and χ 2 = 0, π. We have cos χ 1 cos From the last two equations of system (33) we have the conditions: Eliminating r 1 from (35) using (36) we obtain two equations that are quadratic in r 2 2 and r 2 3 . Eliminating r 1 from the H 2 integral we find one equation that is quadratic in r 2 2 and r 2 3 .
These expressions have to be handled for the range of q determined by u ∈ (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 comoving variables from transformation (28). The normalized variables satisfy the system: We find three families of short-periodic solutions; the constants A, B are real, A 2 + B 2 > 0. 1.
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.
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: with real constants A, B, A 2 + B 2 > 0 and corresponding expressions for the derivatives. We find after linearization ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ẅ 1 + 9w 1 = −ε14d 6 (A cos 2t + B sin 2t)w 3 , We study the stability of this system by normalization to find the eigenvalues of the matrix (omitting the factor 7ε/2) produce first order approximations of the characteristic exponents of system (39). For the eigenvalues we find apart from the factor 7ε/2: A sufficient condition for the complex case C to arise is 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 = 0, d 9 = 0; u = 0 gives for the x 2 normal mode purely imaginary eigenvalues with multiplicity two. As u increases (d 9 = 0), the eigenvalues move from the imaginary axis into the complex domain. This is part of the Hamiltonian-Hopf bifurcation, see Fig. 7.
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.
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 ∈ (0, u 1 ) the nature of Fig. 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 Fig. 8 The ω = 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; ε = 0.2 Fig. 9 The ω = 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, ε = 0.2. The behaviour is different from the case 0, see Fig. 4, as in this case the normal form is not integrable 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.
The Periodic Solution for x 2 (t) = 0 For the periodic solution (38) we put: Transforming and substitution into system (32), we find after linearization: To investigate stability we normalize near the periodic solution; apart from a factor 7ε/2, this produces the matrix: ⎛ 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.

Experiments for Two Cases with u > 0
The first order normal form of system (31) produces qualitatively the same dynamics for u ∈ (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)): 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: We have the case: so that the x 2 normal mode is complex unstable; see Fig. 6. The H 2 (t) time series are shown in Fig. 10. 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 (H H ), 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: With these mass (a i ) values the symplectic transformation of Sect. 3.6 to system (31) produces the expression: If d 9 = 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 ω = 1 vertex.
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.

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 H H instability. Both cases, H H 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]: 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].

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 ≥ 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 ≥ 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 α-chain with four particles in 1 : 2 : 3 resonance does not contain the case of the second normal mode with H H 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 α-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 α-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/ε 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.
In Sect. A.1 we discuss a few facts that are valid for any number of particles n ≥ 3. This enables us to start studying resonances in arbitrary long chains. In Sect. A.2 we give an overview of the relations between frequency ratios and mass ratios valid for all systems with four particles. This might be useful for studying network models as considered in [19]. In Sect. A.3 we discuss some details relevant for the analysis of the 1 : 2 : 3 ratio of frequencies.

A.1 Results for n Particles
Our main general result is the following.
Proposition A.1 If the number of particles is equal to 3 then each choice of eigenvalues λ 1 ≥ λ 2 > λ 3 = 0 occurs for some positive diagonal matrix A 3 .
If n ≥ 4, there are choices of eigenvalues λ 1 ≥ λ 2 ≥ · · · ≥ λ n−1 > λ n = 0 for which there are no positive diagonal matrices A n such that A n C n has these eigenvalues.
To show this we use Lemma A.2, which gives the description of the top coefficients in (7), and Proposition A.3, which restricts the fibers to compact sets in R n .
Lemma A. 2 The polynomials p n−1 and p n−2 have the form indicated in (7).
Proof If we replace the entries −1 at positions (1, n) and (n, 1) in C n by 0 we obtain the Cartan matrix C n for the root system of type A n . (See, e.g., [2, Déf. 3 in 1.5 of Chap. 6, and Planche I].) The determinant of C n is known to be n + 1.
If all a j are non-zero, the characteristic equation is equivalent to det(C n − λA −1 n ) = 0. We determine first the factor of (−λ) n−1 . In the expansion of the determinant the term with λ at all diagonal positions except at (j, j ) is equal to 2 i =j −λa −1 i = 2(−λ) n−1 a j /(a 1 a 2 · · · a n ).
For the factor of (−λ) n−2 we have contributions of two types: Two diagonal positions j and j + 1 (modulo n) lead to a contribution of the form det(C 2 ) i =j,j +1 (−λ)/a i . Two non-adjoining diagonal positions j 1 , j 2 contribute 2 · 2 i =j 1 ,j 2 (−λ)/a i . This leads to the description of p n−2 (a).
Then (a) If η < 1 2 − 3 4n , then Q η is a compact quadric in the hyperplane a 1 + · · · + a n = 1 2 in R n with a non-empty intersection with R n >0 .
Proof Let P = C n − 6I n + 4E, with E the n × n-matrix with all elements equal to 1. Then, considering a = (a 1 , . . . , a n ) as a row vector, we have There are orthogonal matrices U such that U T C n U = Λ, where Λ is the diagonal matrix with the eigenvalues λ j of C n on the diagonal. We put the eigenvalue 0, with eigenvector e as the last one. Then e T = U(0, . . . , 0, √ n) T , and the last column of U is 1 √ n e T . This gives Because of (44) the points in the hyperplane p n−1 (a) = 1 can be described as We write x = (x 1 , . . . , x n−1 ). We find the equation So the points run through a quadratic set in the hyperplane p n−1 (a) = 1. The eigenvectors of C n can be chosen as (ζ k , ζ 2k , . . . , ζ nk ) with ζ = e 2πi/n , which leads to eigenvalues 2 − 2 cos 2πk/n ∈ [0, 4]. So the 6 − λ j are strictly positive. The equation becomes In case (b) in the proposition the single point x = 0 corresponds to 1 2n e ∈ R n >0 . As η decreases the quadric expands in all directions, some of these stay inside R n >0 .

A.2 Results for systems consisting of four particles
For n = 4 Eqs. (8) determine whether points of the fibers exist. In particular, a (scaled) choice of eigenvalues leads to parameters ξ, η > 0 which determine the equations for the fiber. We first consider the values of (ξ, η) that can occur: where Illustration in Fig. 13.
The image X is the region enclosed by these boundary curves.
The points (ξ, η) for which the fiber is non-empty form a subset of the region in Proposition A.4. The fiber is empty for ( 1 27 , 1 3 ). We give a description of the set of (ξ, η) corresponding to non-empty fibers. We omit the proof. (It can be given along the same lines as  (ξ, η) correspond to a non-empty fiber. On the right is again the region in (51), with the subregion in (52) indicated by the dashed line. The points strictly to the right of the dashed line correspond to compact fibers that of Proposition A.4, but takes much more work.) Below we give a computational scheme for the determination of the fiber for given eigenvalues λ 1 , λ 2 , λ 3 . Following this scheme it becomes clear whether the fiber is empty in any individual case.

A.2.1 Spherical coordinates on ellipsoid
In the case n = 4 we may take the orthogonal matrix in the proof of Proposition A.3 in the form  Instructions to compute fibers for the case n = 4. In these instructions we assume that η = 4ξ + 3 16 . Otherwise we also have to consider p 13 ∈ (0, ξ) and investigate whether this leads to further solutions This approach leads to the resonances in Table 3. We discuss the application to the resonance 1 : 2 : 3, and as an example also to the resonances 1 : 1 : 2, 1 : 3 : 6, and 2 : 3 : 4, which are typical for the different cases in Table 1.
Proof We have We try to solve (C − λA −1 )v = 0 with v = (p, x, q, y). The first and third lines give x + y = μ −1 1 p = μ −1 3 q. Similarly, we get p + q = μ −1 2 x = μ −1 4 y. Since λ is an eigenvalue of A 4 C 4 there are non-zero solutions, for which x and y both have to be non-zero. So there is a solution with x = μ 2 . Then we obtain the vector in the lemma. Now we take for a j the expressions in (9). It is clear that 2a j (u) − λ i is not identically zero in u for any of the eigenvalues λ 1 , λ 2 , λ 3 , and λ 4 = 0, and any a j . So we obtain vectors v i (u), 1 ≤ i ≤ 4, that are eigenvectors of A 4 (u)C 4 for the eigenvalue λ i for generic values of u.
The vectors w i = A 4 (u) −1/2 v i are eigenvectors of A 2 (u) 1/2 C 4 A 4 (u) 1/2 . Since the four eigenvalues are different, the w i are orthogonal. We takew i = n −1 i w i with n i = √ w i · w i to get an orthonormal basis. There is the freedom to choose the sign. We multiplyw 1 with −1, to get consistency with our earlier computations. Thew i can be chosen as the columns of the orthogonal matrix U(u). Then the vectors v i = A 4 (u) 1/2w i = n −1 i v i (u) are the columns of the transformation matrix L(u) = A 4 (u) 1/2 U(u). The construction of the v i allows the components to have singularities. The orthonormalization removes all singularities, so the matrix elements of L(u) are continuous functions on