Study of one-dimensional nonlinear lattice

In this article a brief review of the theory of one-dimensional nonlinear lattice is presented. Special attension is paid for the lattice of particles with exponential interaction between nearest neighbors (the Toda lattice). The historical exposition of findings of the model system, basic equations of motion, special solutions, and the general method of solutions are given as chronologically as possible. Some reference to the Korteweg-de Vries equation is also given. The article consists of three parts. Firstly, the idea of dual system is presented. It is shown that the roles of masses and springs of a harmonic linear chain can be exchanged under certain condition without changing the eigenfrequencies. Secondly, the idea is applied to the anharmonic lattice and an integrable lattice with exponential interaction force between adjacent particles is obtained. Special solutions to the equations of motion and general method of solution are shown. In the last part, some studies on the Yang-Yang’s thermodynamic formalism is given.


Introduction.
It was around 1960 when we organized a group of physicists especially interested in the results of computer experiments just becoming available. 1) It was shown that some perturbation methods rather frequently used in physics failed in leading to correct results obtained by numerical computations. For example, the vibrational states of a lattice studied by computer exhibited sharp localization around light-mass impurities, whereas perturbation method predicted only small spreading of the frequency spectrum. 2)-4) As another example, enhancement of heat flow due to nonlinearlity of the interaction force between crystal atoms was shown by the computer, whereas common theories based on approximate method predicted the opposite effect. 5) Dual lattice. I first examined the eigenfrequencies of the simplest systems, which were onedimensional lattice with linear interaction between particles. I examined lattices which contain mass * ) 5-29-8-108 Yoyogi, Shibuya-ku, Tokyo 151-0053, Japan. * * ) Recipient of the Japan Academy Prize in 2000.
impurities (different masses) with those containing interaction impurities (different force constants), to find that different lattices of the different sort may have the same frequency spectra. I considered the reason, and found the following mechanical theorem. 6) Consider a linear lattice of N particles with masses m 1 , m 2 , ... , m N . If the force constants of the springs are K 1 ,K 2 ,... ,K N , the total Hamiltonian is where x j and p j are canonical coordinate and momentum. We put x 0 = 0, x j = r 1 + r 2 + ··· + r j . Then the relative coordinate r j = x j − x j−1 can be used as the generalized coordinate, together with the conjugate momentum s j = ∂H/∂ṙ j . We see Further we may exchange the roles of the generalized coordinates and momenta writing M. Toda [Vol. 80, r j = P j /a, s j = −aQ j [3] here a is an arbitrary constant. Then the Hamiltonian is transformed to Q N +1 = 0, [4] where 1 m * The lattice thus obtained consists of particles with the masses m * 1 , m * 2 , ... and the spring constants K * 1 ,K * 2 ,.... Lattices expressed by H(p, x) and H(P, Q) have the same spectra, and may be called dual to each other.
Nonlinear lattice. The problem of wave motion in nonlinear media is interesting not only as a purely mechanical problem, but also in connection with many physical phenomena such as shallow water waves, plasma waves, and heat conduction in crystals. Thus nonlinear phenomena have infinite variety compared with linear cases. It seems quite important to find feasible mathematical models to extend our way of maneuvering nonlinearity. Vibration of a system of particles joined by harmonic springs can be described by superposition of normal modes which are mutually independent. If we excite a normal mode, its energy is not transfered to other normal modes. The system of harmonic oscillators never reaches the state of thermal equilibrium, and is non-ergodic. Since the presence of nonlinear terms may make the calculation insurmountably complex, usually it is assumed that the nonlinear terms guarantee the ergodicity and approach to the state thermal equilibrium.
Fermi, Pasta, and Ulam (FPU) intended to verify this expectation by computer experiments. 7) Unfortunately, because of the ill conditions after the world war II, I could not see FPU's paper at that time, because it was printed only as a report of the Los Alamos research center. However I had some information about FPU through the works by J. Ford and some others. 8), 9) It was that, contrary to the expectation of FPU, it was clarified that onedimensional nonlinear lattices marvellously sustained the character of linear lattices.
I happened to believe as follows: There will be some model system, expressible in terms of simple mathematical formula, and admits truely exact analytic solutions. I started to seek out the soluble model, and actually found it before 1967. 10)- 15) The equations of motion of a uniform onedimensional lattice are usually written in the form where V (r) stands for the interaction potential between adjacent particles. I thought it could be more favorable to use the transformation described in the preceding section, which exchanges the role of the interaction terms with that of the momentum. Then the Hamiltonian for the above equations of motion turns to be The canonical equations of motion are then given aṡ We confine ourselves to the case where the last equation is solved for r n in such a way that Then we have the fundamental equation Sometimes it is convenient to use S n = t s n dt [11] to write the above equation in the form Finding of the integrable lattice. We have to find the potential function V (r), or its inverse function χ(ṡ n ), which makes the above equations of motion satisfied for some non-trivial solutions. It is also hoped that the form of V (r) has some similarity with that of the intermolecular potential.
The above equations of motion looks like a recurrence formula for a periodic function by which s n+1 is derived from s n−1 and s n . In the linear case s n can be trigonometrical (sinusoidal) functions. Then, for nonlinear lattice, how about elliptic functions? I tried for some days in vain. Then an idea came about, and I tried the square sn 2 (u) of the Jacobian elliptic function sn(u). I found the addition formula for the Jacobian elliptic function where all of the elliptic functions sn, cn, dn are of the same modulus k, and by definition We define a function ε(u) by [15] and use the formula [16] to have ε (u) = dn 2 u, ε (u) = −2k 2 sn u cn u dn u. [17] Integrating [13] with respect to v, we get . [18] Though ε(u) is not a periodic function, Jacobian zn function and the derivative defined by [19] are periodic functions with the period 2K, where K and E are respectively the complete elliptic integrals of the first and the second kind. Using these functions we rewrite the addition formula as . [20] This is to be compared with [10].
Then we find that we have obtained a periodic wave given as [21] where b is a constant, and u = 2 νt ± n λ K, v = 2K/λ. [22] Here ν (the frequency) and λ (the wavelength) are constants, and since du = 2Kν dt, by comparing [10] with [20] and [21], we get [23] where b and σ are constants. χ(ṡ) = −mr is the inverse function ofṡ = −V (r), which must not contain ν and v. It means that the factor is a constant. Refering to [23] we find that [25] and further by [9], that [26] or, solving forṡ, and referring to [8] we obtain Integrating, we finally obtain V (r) = a b e −b(r−σ) + ar + const. [28] In what follows, we use simpler expression by putting σ = 0 to write . [29] This is called the Toda potential (Fig. 1). For small r we have Thus for sufficiently small oscillation, the lattice looks like a linear lattice with the spring constant κ = ab. For somewhat large motion, the lattice will behave like a system of hard spheres. We shall summarize the results. The equations of motion given by [6] and [29], is called the Toda lattice (exponential lattice) equation, which is also written as In terms of the dual lattice, by [26] we have the equations of motion where by [33] the relation between r n andṡ n is given asṡ n = a e −br n − 1 . [36] The periodic wave solution is given by [21], [36] and [19] as (cnoidal wave, Fig. 2) with the dispersion relation which is given by [24] as 2Kν = ab m where K and E are complete elliptic integrals [39] Thus I found the nonlinear lattice and its periodic solutions at the same time.

Continuous limit.
Sometimes it is seen that the continuum approximation gives good results. In this case it is convenient to use the operator rule e ±d/dn f (n) = f (n ± 1). [40] Then we see the simplified version a = b = m = 1 of the equations of motion [32] can be written as If we straightly expand the left hand side in powers of ∂/∂n and neglect higher powers of r, we can rewrite the above equation as Thus, for the wave advancing to the right, we have and further expanding sinh, If we change the units and signs, we can arrive at an equation of the form This is the famous Korteweg-de Vries (KdV) equation originally derived to describe the shallow water waves. 16) They gave periodic solution which is similar to the periodic wave given in the foregoing solution.
Besides they gave solution of the form (Fig. 3) which is a solitary wave solution or soliton solution. N. Zabusky and M. D. Kruskal 17) studied the KdV equation by computer experiment. They used the cyclic boundary conditions and started from the initial state such as [47] They found the wave split into a group of solitary waves (solitons), each of which proceeded nearly independently and after a while the wave recovered the initial state. If the wave consists of a small number of solitons it will come back to the initial state after a time, which is nearly equal to the least common multiple of the recurrence times of the solitons.

Solitons.
A soliton of the nonlinear lattice can be thought as the limit of infinite wave length. It leads to a soliton solution The velocity of the soliton is which is larger for the soliton of larger height. In a paper P. D. Lax wrote two-soliton state of the KdV equation without any proof. But I found a hint from it to rewrite [48] in the form 12), 18) with β determined by the equations of motion as β = ab m sinh κ. [51] Then I found that two soliton state is given by where A, B, β and γ are determined from the equations of motion. There are two cases. One is the head-on collision of two solutions, given by [53] The other case is for two solitons running in the same direction: [54] Integrability. Ford 19) numerically examined the integrability of the lattice with exponential interaction between particles. He took cyclic lattice of three particles, and applied the Poincaré method of mapping (surface of section) of trajectories in the phase space. He found that the trajectories are smooth and did not become erratic even if the energy is raised extremely high, indicating that the lattice was integrable (Fig. 5).
The Hamiltonian of a three-particle cyclic lattice with exponential interaction can be written in a dimensionless form as He applied the transformation which diagonalizes the corresponding harmonic lattice. Then, by rescaling as he obtained the equations of motion Ford numerically integrated [58] and had the Poincaré mappings of the three-particle exponential lattice. Fig. 5(a) is for the energy E = 1, and Fig. 5(b) is for E =256. He examined mapping up to E =56000, and always he had smooth curves with no indication of stochastic behaviour. Thus the numerical works strongly suggested that the lattice with exponential interaction is integrable. In other words, it admits the so-called third integral besides momentum and energy.
These results were in strong contrast to the behaviour of usual nonlinear system, which are in almost all cases chaotic, non-integrable.
For instance since 1964, we had the so-called Hénon-Heiles system, 20) which is equivalent to a cyclic lattice of three particles with the Hamiltonian [59] The second term on the right hand side is the harmonic interaction and the third is the cubic nonlinear term. After performing the transformations [56] and [57], we obtain the equations of motion [60] As easily seen, 18) if we expand the right hand side of [58] up to the second powers of q 1 and q 2 , we get [60].
But in both cases the trajectories are quite different. In the case of [55] they are entirely smooth showing integrability of the exponential lattice. On the contrary, the trajectories of [60] become chaotic as shown in Fig. 6(a), (b) and (c) with increasing energy E, indicating non-integrability of the lattice with cubic nonlinearlity. This is a very clear example of the transition between integrable and non-integrable systems.
Conserved quantities. The equations of motion of the cyclic three-particle exponential lattice [55] can be written aṡ where [62] From these equations we see that I 1 = P 1 + P 2 + P 3 are conserved quantities. I 1 is the total momentum, I 2 is related to the total energy (E = (1/2)I 2 1 − I 2 ).

M. Toda
[Vol. 80, I 3 is the so-called third integral of motion, which is not interpreted in terms of momentum and energy. The existence of the third integral in this case was proved by the numerical study of the preceding section. M. Hénon 21) and H. Flaschka 22) were stimulated by the results of the numerical study performed by Ford. 19) Hénon proved that uniform cyclic exponential lattice has as many conserved quantities as the number of particles in the lattice. If we write them for the case of three particles, we get I 1 , I 2 and I 3 as mentioned above.
In the same year, Flaschka independently proved the same thing by a method, which we shall explain in the next section. 22) Matrix equation of motion. We continue dealing with cyclic exponential lattice of N particles. We first rewrite the equations of motion 23) [65] If we put we havė a n = a n (b n − b n+1 ),ḃ n = 2(a 2 n−1 − a 2 n ). [67] For a cyclic lattice of N particles forming a ring the equations of motion can be written in the matrix form as This is called the Lax formalism. 24) It was first introduced for the KdV equation, and afterwards used in many cases. It is assumed that B is anti-symmetric. In our case we have Matrix components a n and b n are functions of the displacements and momenta, which are functions of time. Therefore matrices L and B are functions of time. We write them as L(t) and B(t), and the initial values as L(0) and B(0). We also write eigenvalue of L(t) as λ(t) to express the t dependence, though we will find that λ is independent of t. If we write the eigenfunction as ϕ(t) (N × 1 matrix), we have

L(t)ϕ(t) = λ(t)ϕ(t)
[71] where L(t) is a N × N matrix. λ(t) are given as the roots of the equation where I is the N × N unit matrix. Since B is antisymmetric the matrix U defined by is unitary. That is we have Further, by [69], we have so that [76] Thus L(t) and L(0) are unitary equivalent. [78] Comparing this result with det[L(0) − λ(0)I] = 0 we get Therefore the eigenvalues are independent of time.
Thus it is shown that the motion in the lattice conserves its spectrum (isospectral deformation). Expanding [72], we get the eigenvalue equation where c k are functions of a n and b n . From the above simultaneous equations we can write c k as functions of λ j , which are time-independent. Therefore c k are also conserved quantity. It is verified that c k are essentially the same as Hénon's constants of motion. We may also obtain these constants of motion as Thus it is shown that cyclic exponential lattice of N particles has the same number of conserved quantities, which means that it is integrable.
Inverse scattering method. If we write down the eigenvalue equation [71], we have a n−1 ϕ n−1 + b n ϕ n + a n ϕ n+1 = λϕ n . [82] This equation has its counter-part in the case of the KdV equation. Consider a soliton with negative value This is an example of solutions of the KdV equation  25) which solves the initial value problem and called as the inverse scattering method.
For the case of discrete lattice, H. Flaschka 23) applied the inverse scattering method to the Toda lattice. The motion in the infinite Toda lattice was also solved by E. Date and S. Tanaka. 26) Using the theory of hyper elliptic integrals M. Kac and P. van Moerbeke showed similar method to solve the three particle cyclic lattice. 27) Probability distribution. We now turn to another problem. Around 1940, that is when I graduated from the university, there were senior friends in Japan studying thermodynamical problems such as the distribution function of molecules in liquids and possibility of phase transition to gases or solid state, and so on. T. Nagamiya 28) considered one dimensional chain of particles where the total potential is the sum of the potential V (x n − x n−1 ) between adjacent particles. The distribution g n (x) of the n-th particle is governed by for the canonical ensemble (f = pressure). Nagamiya calculated the distribution function g n (x) for some special potentials, such as the harmonic potential and hard sphere potential. The results were to be compared with the distribution function of molecules revealed by the X-ray experiments of these days. 29) invented a simple method for treating the statistical mechanics of one dimensional chain of N molecules. He considered the configurational integral of the system (β = 1/kT )

Partition function. H. Takahashi
Instead of the variables x 1 ,x 2 ,... ,x N , we introduce and observe the Jacobian [91] Then we see that where we have written Q(β) for which is the configurational part of the partition function. By the standard statistical mechanical argument we see that the average length of the system is given as [94] and the energy per particle by We see that the length l is always a single valued function of the pressure, which means that no phase change occurs in one dimensional system. In classical mechanics, the Hamiltonian per particle is [96] If we write the total partition function as ζ(β), it is given by [97] Exponential lattice. In the case of the exponential lattice, we can use the simple expression [98] We have [99] To perform the integration we put e −x = t, and get where Γ(ρ) is the gamma function of the order ρ = βf. [101] Thermodynamics. The Bethe anzatz was first invented for quantum-mechanical one-dimensional spin system, and extended to one-dimensional system of particle interacting via a repulsive delta-function potential. It was further developed by C. N. Yang and C. P. Yang to the thermodynamical system with given temperature and pressure. 30) Following N. Theodorakopoulos, 31) the classical limit of the Yang-Yang's thermodynamical equation written for the chemical potential µ is with the subsidiary condition where f is used for the pressure to discriminate it from momentum. In the above integral equation the factor exp(−β ) is related to some kind of excitation, and θ (p) stands for the derivative of the two-body scattering phase shift θ(p) for the interaction potential, and p denotes the relative momentum. In the classical limit, and for the exponential interaction, we use θ (p) = −2log |p|. [104] then the above equation is rewritten as [105] By setting we have [108] [109] This is nonlinear with respect to φ(x) for ψ(x) implicitly contains φ(x). M. Opper solved this equation. 32) The result can be written as where the constant C is determined as where Comparing the both sides of [110], we get to see [115] Since J = ∞ 0 cos xt e −t 2 /2 t ρ dt [116] by integrating partially we have [117] Similarly we have [118] So that and [109] is satisfied.
Therefore the solution φ(x) and ψ(x) of the Yang-Yang equation [110] for the exponential lattice is given by [114] in terms R(x) and J(x). Then the chemical potential µ for each particle is given by [107].
Chemical potential. However, since in [107] the chemical potential µ is independent of x, we may put x = 0 in the right hand side of this equation to estimate µ, so that with [121] Following the standard theory of classical statistical mechanics, we see that e −βµ is nothing but the partition function ζ(β) calculated as [97] and [100]: On the other hand, we have J(0) = R (0) = 0, and from [113] and [114] [123] Therefore [120] gives [124] Thus we have reproduced the formula ζ(β) = e −βµ . Now, our problem is to clarify the meaning of ψ(0) and φ(0) by identifying them to well-known statistical mechanical quantities such as the phase integral or some part of partition function.
Binary interaction. In [120] the factor e 2ψ(0) looks indicating pairing of two factors of ψ(0), which seems to come from the pair of interaction potentials giving rise to with eξ (η) = cosh(η/2), [126] The transformation from (x,y) to (ξ,η) is a canonical transformation, [127] Integrating we find the formula where (ρ = βf ) . [129] Here I have used a special notation ζ Q , which incidentally appears later again.
In [128] the phase integral at the reciprocal temperature β is related to the square of the integral at β/2, or rather to the integral for one half of the potential, 1 2 (e −x + fx), Γ(ρ/2). [130] From the above consideration, it seems clear that ψ(0) has something to do with Q(β/2) above, though they look quite different from each other.
Final result for ψ(0). Since ψ(0) is related to the phase shift θ(p) of binary collision, we turn to dynamics rather than sticking to the configurational integral Q(β).