J. Phys. Chem. 1988, 92. 3154-3163
3154
V. Conclusion We have demonstrated that generalized A Q is a useful technique for finding accurate vibrational eigenvalues. For all the systems considered (except the anharmonic oscillator), including just the first correction in h leads to a significant improvement in energy eigenvalues. Including corrections of high enough order to generate fully quantum results takes little more time than calculating the first correction. As discussed above, this is because high-order corrections in h are polynomials of relatively low degree. We are currently investigating ways of incorporating corrections in h with less effort. For systems with apparently convergent perturbation expansions, it would be consistent to calculate low-order corrections in 6 quantum mechanically, while finding high-order corrections in 6 classically. Such a scheme would probably take only slightly longer than an entirely classical calculation, while generating much better results. Furthermore, generalized AQ could profitably be used to calculate the per-
turbation expansions of Hamiltonians with complicated functional forms. Since the Moyal bracket is expressed entirely in terms of derivatives, it is much easier to calculate the Moyal bracket of two complicated functions to a given order in h than it is to calculate the commutator exactly. Our motivation for exploring the technique, however, transcends the purely computational. Generalized AQ offers a framework within which both classical and quantum perturbation theory can be studied on the same footing. This helps to illuminate the close relationship between the theories, and enables h-dependent effects to be introduced in a controlled manner. Mock phase space is conceivably the appropriate setting for deriving a “quantum KAM theorem”, which would be, to lowest order in h , the classical KAM theorem. Finally, we expect generalized AQ to be of use in examining the summability of classical vs quantum perturbation theory. An examination of these problems, as well as the construction of wavefunctions with generalized AQ, is under way.
(93) BaEiE, Z.; Light, J. C. J . Chem. Phys. 1987, 86, 3065. (94) Guckenheimer, J.; Holmes, P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields; Springer-Verlag: New York, 1983. (95) Voros, A. Ann. Inst. Henri Poincare 1977, 26, 343. (96) Voros, A . In Lecture Notes in Physics; Casati, G., Ford, J., Eds.; Springer-Verlag: New York, 1977. No. 93.
Acknowledgment. We are very grateful to John Frederick for communicating the results of his variational calculations. This work was supported by NSF Grant CHE-8410685. Computations reported here were performed in part on the Cornel1 National Supercomputer Facility, which is supported in part by the NSF and IBM Corp. G.S.E. acknowledges partial support of the Alfred P. Sloan Foundation.
From Classical to Quantum Mechanics with Hard Chaos Martin C. Gutzwiller IBM T. J . Watson Research Center, Yorktown Heights, New York 10598 (Received: June 15, 1987; In Final Form: October 14, 1987)
Most of the interesting dynamical systems in physics and chemistry are characterized by an intimate mixture of two extremes: integrability on one hand, and “hard chaos” on the other. The first of these is well understood, particularly as one tries to connect classical and quantum mechanics. The second, however, has been the exclusive preserve of the mathematicians and suffers, moreoever, from a dearth of good physical examples, although it is probably more frequent in nature than its integrable counterpart. The energy surface in the phase space of a harshly chaotic system foliates into two sets of smooth manifolds, the stable and the unstable ones, each of dimension equal to the number of degrees of freedom. A classical trajectory is the intersection of a stable and an unstable manifold. Its physical properties can be described by a simple coding scheme that gives a unique account of its whole history. As an example, the anisotropic Kepler problem (AKP), isotropic Coulomb potential with anisotropic mass tensor, is discussed in some detail, because it allows the two sets of manifolds to be constructed by varying only a single parameter and because its code is the simplest imaginable, consisting of all binary sequences. The transition to quantum mechanics, using Green’s function and its trace, leads to expressions that resemble very closely the grand-canonical partition function of a one-dimensional king chain. The energy levels are the zeros in the complex ( l / k T ) plane and should be located on the imaginary axis. This analogy is directly related to the underlying hard chaos and is followed up by studying the locations of these zeros for various coupling parameters in the Ising chain. It appears that some further restrictions have to be invoked, in order to guarantee a reasonable transition from classical to quantum mechanics in the case of hard chaos.
1. Introduction The term “hard chaos” is meant to describe one of the two extreme conditions that are encountered in the study of phasespace dynamics. The first of these is the regular behavior which all of us have been taught in school as if there was nothing else to be expected in a mechanical system. It is characterized by the occurrence of invariant tori (cf. ref 1) to which the trajectories in phase space are bound, a situation typical of all integrable systems such as the electron attracted by two fixed nuclei. While this property of a dynamical system allows it to be effectively solved, both classically and quantum mechanically, it turns out to be highly exceptional. ( 1 ) Arnold, V. I.; Avez, A. Ergodic Problems of Classical Mechanics; Benjamin: New York, 1968. Arnold, V. I. Mathematical Methods of Classical Mechanics; Springer: New York, 1978. Lichtenberg, A. J., Liebermann, M. A. Regular and Stochastic Morion; Springer: New York, 1983.
0022-3654/88/2092-3154$01.50/0
A similarly unusual, but basically simple behavior forms the opposite extreme of utter randomness, no better than the throw of a dice although completely deterministic. This condition, to be called “hard chaos” with the adjective “harshly chaotic”, will be discussed in the present paper because it is not as well-known as the first. Most mechanical systems show a mixture between integrability and hard chaos which I like to call “soft chaos”. It is like a very intimate mix of the two extremes which seem to coexist in very close, complementary regions of phase space. The result is most often understood as a perturbation of the regular behavior, but it could equally well be viewed as a corruption of hard chaos. Before going into any detailed discussion, however, some natural boundaries have to be mentioned within which these ideas can be applied. They may have a much wider domain, but there is no point in claiming more territory than can be staked out with some well-understood examples. There will be no dissipation; the 0 1988 American Chemical Society
From Classical to Quantum Mechanics with Hard Chaos dynamical systems are Hamiltonian. They will also be conservative; their Hamiltonian is time-independent. Finally, the number of degrees of freedom is minimal, namely, only two, although many of the concepts can easily be extended to any finite number of degrees of freedom. As a last voluntary restriction to be observed in this paper, I will talk about motions in a compact space or, in other words, about bound-state problems. A number of things change when a particle moves in a noncompact space with hard chaos, Le., in a scattering problem. The quasi-random nature of the events remains the same, however, and is described in a similar manner. Sections 3, 5, and 7 give a brief review of the anisotropic Kepler problem, the trace formula to calculate the spectrum, and the discussion of symmetries and stability exponents. The remaining sections present new material, the idea of hard chaos in section 2, the construction of the foliation in phase space in section 4, the reduction to a grand canonical partition function in section 6, the use of transfer matrices in section 8 and the discussion of the possible spectra in section 9.
2. Hard Chaos A conservative Hamiltonian system with two degrees of freedom has a four-dimensional phase space, but each trajectory belongs to a three-dimensional manifold of fixed total energy E . Classical mechanics takes place on this “surface” of constant energy where the equations of motion define a vector field, and the trajectories are obtained by tying these vectors into smooth curves. Even a reasonable-looking vector field can result in an amazing tangle of trajectories so that we face the central issue: How can we represent what is going on in some understandable fashion? The secret is, of course, Poincare’s surface of section. Rather than surveying a whole three-dimensional manifold, we take a two-dimensional slice which is everywhere transverse to the flow of the trajectories. As we compare two neighboring trajectories in their consecutive intersections with the surface of section, we get an idea of the flow although we lose the information concerning the length of time between these consecutive intersections. The mapping from one to the next is well-known to be area-preserving provided we choose canonically conjugate coordinates in the surface of section. The main question for all the subsequent development is the following. Can the energy surface be subdivided into smooth submanifolds with simple shapes so that each trajectory stays within one of them for its whole duration? If so, these submanifolds will intersect the surface of section in smooth lines, and the consecutive intersections of any given trajectory will be on one and the same of these lines. The phase space is then said to be foliated into twedimensional leaves, and the problem of finding the trajectories has been reduced almost to the point where the complete solution is in sight. Anybody who has numerically integrated a Hamiltonian system knows how clearly one can immediately detect the presence or absence of such a foliation. The hydrogen atom in a magnetic field is a case where the foliation gets progressively destroyed as the field becomes stronger, but ultimately the foliation returns as the electron settles into Landau levels (cf. Gay2). If the system is integrable, one finds the famous invariant tori which show up as concentric deformed rings, commonly called islands. The Kolmogoroff-Arnold-Moser (KAM) theorem tells us how much of this foliation survives when one allows for a sufficiently weak as well as well-behaved perturbation. There are, however, weak but not sufficiently well-behaved perturbations of a regular system which prevent the conclusions of the KAM theorem from being valid. The anisotropic Kepler problem (AKP) to be discussed at some length has this property which leads us straight into hard chaos. Hard chaos is again characterized by a foliation of the energy surface, but it looks quite different, and its function in organizing the flow on the energy surface is entirely novel. There are two sets of leaves which form something like a deformed square pattern (2) Gay, J. C. In Atoms in Unusual Situations; Briand, J. P., Ed.; Plenum: New York, 1986.
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
3155
in the surface of section. For reasons to be explained shortly, one set is called the stable manifolds and the other the unstable manifolds in the dynamical system. Each trajectory r is the intersection of a stable with an unstable manifold. A neighboring trajectory in the unstable manifold drifts away from r at an exponential rate with increasing time, but it approaches r in the same manner as time goes to -a. A neighboring trajectory in the stable manifold has the opposite behavior with respect to time. Such a situation seems rather unlikely, but it is actually quite well understood mathematically. Moreover, even a casual look at the standard examples seems to suggest a very great variety of possible mechanical systems in this category. For comparison’s sake recall that the prime examples of integrable systems are the hydrogen molecular ion, the geodesics on a triaxial ellipsoid, the Toda lattice, and the heavy symmetric gyroscope. While I do not have anything harshly chaotic to compare with these examples except the AKP, it is instructive to look at the super-regular case of geodesics on a sphere (mechanically, the geodesics are the trajectories of a free particle constrained to move on a given surface). Apart from the scaling by the radius of the sphere, this is a unique case in geometry which is characterized by its constant positive Gaussian curvature. On the other hand, there is a tremendous variety of compact two-dimensional Riemannian surfaces of constant negative Gaussian curvature, not mutually isometric, and requiring several continuous parameters for their definition, apart from their radius of curvature (cf. ref 3). The physical interpretation of these purely geometric creatures and their hard chaos will be discussed elsewhere. The purpose of making up this new category of harshly chaotic systems is to examine whether their quantum-mechanical analogues have characteristic features as well defined as the integrable ones. Their classical behavior is based on a very special foliation of phase space which must have some consequences for the wave functions. There is still a lot of progress to be made. The present report deals only with the transition from classical to quantum mechanics when hard chaos prevails in the classical limit.
3. The Anisotropic Kepler Problem This model Hamiltonian deals with the electron in the neighborhood of a donor impurity in a silicon or germanium crystal. The author was the first to study its classical mechanics and to discover its ideally chaotic nature ( G u t ~ w i l l e r ~ ~Most ~ ~ ) of . the material in this and the next section is a review of these early papers. Devaneys has given a more strictly mathematical treatment of the AKP. In order to make contact with the experiments one needs the physical quantities charge of the electron e, dielectric constant of the crystal K , longitudinal and transverse effective mass at the conduction band minimum ml and m2,the geometric mean mo = (mlm2)’12, and Planck’s quantum h . The natural energy scale in this system is the “Rydberg” Eo = m o e 4 / 2 ~ 2 hWe 2 . will use the appropriate combinations of these quantities as units in all the subsequent discussion. The Hamiltonian H is written in terms of the momentum (u,u,w) and the position (x,y,z),where the x axis is pointed in the direction of the longitudinal (large) mass m l . The mass ratios 1 = (ml/m2)’I2 and v = (m2/m1)’12 will be used because they have the nice property that pv = 1. The crucial parameter in the whole model is the mass ratio m l / m 2= p / v , which turns out to be 115 for Si and 1120 for Ge:
(3) Selberg, A. J. Indian Math. Soc. 1956,20,47-87. For a more readable account cf. McKean, H. P. Comm. Pure Math. 1972,25,225, and for a very complete discussion cf. Hejhal, D. Lect. Notes Math. 1976, 548. The same ground is covered from a more physical viewpoint very thoroughly by Balas, N. L.; Voros. A. Phys. Rep. 1986, 143(3), 109-240, and with a slightly different emphasis by Gutmiller, M. C. Phys. Scr. 1985, 79,184-192, as well as Contemp. Math. 1986, 53, 215-251. (4) (a) Gutzwiller, M. C. J . Math. Phys. 1971,12, 343-358; (b) J. Math. Phys. 1973, 14, 139-152; (c) J. Math. Phys. 1977, 18, 806-823. (5) Devaney, R.L. J. Differ. Equations 1978,29, 253-268; Invent. Math. 1978, 45, 221-251; Lect. Notes Math. 1978, 668.
Gutzwiller
3156 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
The energy of the electron is E < 0, but since the kinetic energy is homogeneous of degree 2 in (u,v,w),and the potential energy is homogeneous of degree -1 in (x,y,z), the trajectories can be scaled with the appropriate powers of /El. In other words, the structure of phase space is the same for all E < 0, and we have chosen the value as the most convenient to work with. Most of the physical information in the classical behavior and in the transition toward quantum mechanics is contained in the action integral S(q”,q’,E) from a position q’to a position q ” a t the energy E:
S(q”,q’,E) = @=
(*)
x,”‘
( u dx
-2~’.??
15
10
05
>
@
00
+ u dy + w dz)
-0 5
which defines a somewhat awkward looking domain in the (u,x) plane. A simple area-preserving map into a (XU)plane with the formu 1as
transforms the awkward domain in the ( x , u ) plane into the rectangle
Ivl 5 k1’2(7/2)
-1 0
(2)
(3)
2,
i n i
I/2
where the integral is taken along a classical trajectory from q’ to q” at the energy -’/’ as specified in (1). The normalized integral @ takes the value 27r in the ordinary Kepler problem (k = 1) when taken over any closed Kepler ellipse. In the anisotropic Kepler problem (AKP), however, the possible values of @ are the heart of the matter, and they can only be understood provided the structure of the trajectories in phase space is known. Since the two transverse masses in (1) are equal, there is a conserved angular momentum around the x axis, M = yw - zw. When M # 0, the centrifugal potential stabilizes the trajectories to the point where the chaos becomes soft. We shall, therefore, restrict ourselves to M = 0. The trajectories then lie in planes through the x axis, and we can restrict our attention to the plane z = 0; the conjugate pair (w,z) vanishes from the scene, and the Hamiltonian (1) will always be understood with w = 0 and z = 0. Nevertheless, the basic cylindrical symmetry of the problem, even when M = 0, will stay with us. Indeed, the spectrum of the three-dimensional AKP with M = 0 is different from the spectrum of the two-dimensional AKP. The surface of section to be used throughout is the plane y = 0. As a consequence of (1) the coordinate pair (u,x) satisfies the inequality
I4 5
2.0
(5)
This latter domain will be used henceforth. In contrast to my earlier discussions of the AKP, the foliation of the phase space will now be constructed explicitly, without trying to motivate the necessary steps. The results will look like magic, but I believe that this is exactly what one can expect from harsh chaos. The equations of motion are deceptively simple:
where r is the distance from the center; r = (x2 + y 2 ) ’ i 2 . Notice that the force is still radial as in the ordinary Kepler problem, but the acceleration is always turned away from they axis. We will take advantage of this feature. There is a one-parameter family of trajectories which end up in a collision with the origin as in a black hole or, by time-reversal symmetry, come out of the origin as in a white hole of general relativity (Gutzwiller“). As they approach the origin, they start
I
-1.5
‘
I 1
-2.0 -2.0
-1.5
-1.0
-0.5
Wi
1 0.0
0.5
I ~
I
I
I
10
1.5
2.0
X
Figure 1. Typical AKP trajectory coming out of a collision and making 10 intersections with the (heavy) x axis; mass ratio = 5 as for silicon.
hugging they axis, as can be seen from the approximate solutions of (6): k3i2
u
-A,P-3/z, 2
N
-B(2a 2
- 1){2*-5/2
y
3
N
x=AP, 2{
+ ,pa-’
(7)
in terms of the parameter { > 0 along the trajectory. Moreover, one has the relations
(i
a = 4
1
+
( -$), 1-
3A2 B = 4(2a - 3)(4. - 1)
(8)
which show that (Y is not far below 3 / 2 for Si or Ge with l2= 5 or 20. The different collision trajectories are distinguished by the real parameter A which varies from -a to +m and fixes the value of B by the last formula. This explicit expression (7) for a particular set of trajectories in the neighborhood of the origin constitutes only the lowest term in an expansion in powers of {. Further details can be found in the above reference and do not really matter in the present discussion except to note that the larger A the smaller the upper limit of validity for {. When A = 0, the trajectory simply goes up along t h e y axis with y = 2p and x = 0. When A # 0, the trajectory gets deflected from they axis to the same side as the sign of A . In practice the formulas 7 and 8 serve the purpose of providing initial conditions for the numerical integration of the equations of motion ( 6 ) by the standard routines, when the trajectory comes out of or goes into a collision. We will use these special trjactories to construct the foliation of the phase space. 4. The Foliation of Phase Space in the AKP The purpose of this section is to demonstrate how the phase space in a harshly chaotic system like the anisotropic Kepler problem (AKP) carries a rather simple foliation, just as the phase space of an integrable system is foliated by the invariant tori. What seems even more remarkable is the method by which this foliation can be constructed in spite of the harsh chaos. The whole surface of constant energy is covered by an ever finer net of leaves by varying a single parameter, the collision parameter A of the preceding section. In order to explain the process, consider Figure 1, which shows a trajectory for the mass ratio 5 . After getting to the point at y = 0.02 with the help of the formula 7, an ordinary fourth-order Runge-Kutta routine was used to continue the solution of the equations of motion (6) until the x axis had been intersected 11 times. With A > 0 the first intersection occurs with XI > 0, but
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3157
From Classical to Quantum Mechanics with Hard Chaos the subsequent intersections X 2 ,X,,... have arbitrary signs. The first question to be answered is whether these signs come up in any particular order as we let A increase from 0. Without trying to prove this central proposition, the following situation is found: Let us label each trajectory by the sequence Let us then of 10 binaries a,, a2, ..., alo, where ai = sign (Xi). number all these binary sequences by the obvious device of forming the positive integer
+
+ ... + 2 ~ +9
L = 512~1 256~2
(9 1
Since a, = 1, the values of L go from 1 to 1023 for A > 0. It turns out that L increases monotonically from 1 to 1023 as A goes from O+ to 150; each odd value of L is realized in some open interval of A. For practical purposes it is convenient to let A = 0.4 tan ( ~ ~ 1 and 2 ) have a go from 0 to 1, because the intervals for the various values of L are then fairly evenly distributed. The intervals in A (or in a) for two consecutive values of L are separated by a parameter value A which ends in a collision at the tenth intersection or earlier. One can think of such an early mishap as ending the binary sequence a ] , a2, ..., a, prematurely with n C 10. In this manner we find a discrete set of 5 11 trajectories which are labeled by the even numbers from 1 to 1023, and which are generated by well-defined values of A in increasing order. By the same token, we can find the discrete set of trajectories which have a collision at the 1 lth intersection with the x axis and label them with the odd values of L from 1 to 1023. They are somewhere in the middle of the open interval with the same odd number L.
This scheme can be completed, if we add the trajectory for A = 0 which goes up they axis and falls right back down again into a collision. It has no intersection with the x axis outside the origin and gets the label L = 0. Since we have symmetry with respect to t h e y axis, all trajectories with a negative starting parameter A have a unique companion with the parameter -A and a positive label L. In this manner we extend L to the negative values from -1023 to -1, and the monotonic ordering with respect to A is preserved. What happens inside the open interval with a particular label L odd? The first to tenth intersections have always the same sign for X as A increases. The 1l t h intersection starts out with XI, C 0 and ends up with X,, > 0, with a collision at the transition from one to the other regime. This behavior is quite in keeping with the idea that increasing the parameter A pushes the trajectory from left to right. At the same time the value of U,, goes from its minimum - p 1 J 2 ~ / 2monotonically to its maximum +p112r/2 when the collision at the 11th intersection takes place; and after , that collision U,, starts out again at the minimum - p 1 / 2 ~ / 2and increases monotonically to reach the maximum +pll27r/2 at the end of the interval with the particular label L. As the collision parameter A runs through the open interval with the odd label L, the point with coordinates (Xll,UII)in the surface of section traces out two smooth curves, one in each of the two halves, XI C 0 and XII > 0. The coordinate U increases monotonically from -p1I2r/2 to + p L ’ / 2 ~ for / 2 each curve, while the coordinate X stays strictly away from the limiting values of -2,0, and +2. In this manner one gets a family of curves in each half of the surface of section, which are labeled by the odd numbers L from 1 to 1023. These curves do not intersect, although some of them converge at the upper or the lower limit of U. A trajectory which starts with an initial condition on one of these curves has a well-defined past which is described by the binary sequence a,, a2, ..., alowith a collision at the beginning. Since these curves in the surface of section do not intersect, and all of them go from the bottom at U = -p%/2 to the top at U = + ~ ‘ / ~ ? r /they 2 , can be labeled from left to right. But what is the relation of such a label to the number L which goes with the binary sequence a,, az,..., a,,? The curves definitely do not arrange themselves in the order of their labels L but rather in the order of their “backward number” K, which arises as follows. Let us add the binary a,, to the 10 others; a l l = -1 when XI, < 0, and a l l = +1 when XI, > 0. Obviously, a l l determines the half (left or right) of the surface of section in which the curve is located.
,
In each half, the curves with a I o= -1 are to the left of the curves with a l o = + l ; the binary aIodivides the surface of section into quarters. The further construction is now quite clear: the ”backward number” K is defined by the sum K = 1O24UlI + 512U10+
... + 2 ~ + 2
(10)
which takes on all the odd integers from -2047 to 2047 provided we let a, be both -1 and + l . The curves in the surface of section are ordered by their number K. The curves in the (X,U) plane coming from the 1 lth intersection can be complemented by the curves coming from the loth, the 9th, etc. intersections. But in trying to determine their relative locations, we have again to look at them “backward”, so that the last intersection determines in which half of the surface of section they are, the next to last determines the quarter in each half, and so on. The collision sometimes in the past becomes somewhat irrelevant, and we will therefore define a label which is independent of this particular occurrence. Suppose that the equations of motion (6) for a particular trajectory have been integrated forward in time for a number of intersections until the intersection at (Xo,Uo)was reached. The ... going earlier intersections are called (X-l,U-l),(X-2,U-2), backward in time. A binary sequence is defined for the past by setting a_, = sign (Xi) with the label 4 for the past given by m
c ai(1 /2’+’)
E = i=O
(1 1)
If the trajectory came out of a collision at the intersection with number -J, then we set a_, = 0 for i > j . The label is then a rational number with the denominator 2J. Each of the curves in the above discussion has a label, and the curves in the (X,U) plane are ordered according as this label. Let us call these curves the 5 curves because each is characterized by a fixed number 5 which describes its past. We can construct an exactly similar set of curves which are characterized by their future behavior. If we integrate the equations of motion (6) forward in time with the initial conciirion (Xo,Uo),we find the consecutive intersections (X,,U,), (X2,U2),.... from the correwe get the label sponding binaries a, = sign (A’,), m
1 = Xa,(1/2‘)
(12)
I=,
where the sum is simply truncated in the event of a collision. The initial conditions (Xo,Uo)which belong to the same label 7 again lie on a smooth curve which goes monotonically from the bottom value U = -p1I2r/2 to the top value U = +p1I2r/2 while staying away from the vertical lines X = -2,0, +2. Because of the time invariance of the equations (6), these 7 curves are simply the mirror image with respect to the X axis of the curves. The discussion at the beginning of this section shows how one can compute explicitly the curves of constant [ or constant 1, provided the value of [ or of 7 is a rational number whose denominator is a power of 2. These special values are obtained by varying a single parameter A , and the corresponding smooth curves establish an arbitrarily fine network in the surface of section. Figure 2 shows the network down to the multiples 1/ 16 in E. The curves E = constant define the so-called unstable manifolds, while the curves 7 = constant define the stable manifolds of phase space. What does that mean? Consider some point (Xo,U0)which is the intersection of one of the [ curves, say Eo, with one of the 7 curves, qo, in Figure 2. Its predecessor was ([-,,v-J, and its successor will be ( [ , , q , ) . Let us now change toto the nearest p curve, changing 9 by 2-“ while keeping Eo fixed, and see where the predecessor and the successor of this second point lie with respect to the predecessor and the successor of the first one. The change 2-” in vo implies a change in the binary a, according to (12), but not earlier. When going backward to the predecessor, this change of the nth binary will become a change in the (n+l)th binary which causes q.., to change by only 2-(”+,); the predecessor points are closer. On the other hand, in the calculation of q,, the change in the binary a, expresses
3158
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
100
-0 50
I00
C ~ l l i i i o nTrrliclorir\
Figure 2. Stable (9 = constant) and unstable (( = constant) manifolds in the surface of section y = 0 with the coordinates (X,v).
itself in a change of q1 by 2-(“’); the successor points move away from one another. The convergence of the predecessors and the divergence of the successors is exponential with the number of intersections. If the roles of ,$ and q are switched in the above arguments, the predecessors will be found to diverge, while the successors converge exponentially. Trajectories with initial conditions on the same q curve march toward the same future; trajectories with initial conditions on the same ,$ curve arise out of the same past. As one keeps filling up the surface on section with collision trajectories, each point becomes the intersection of such a stable manifold (q = constant) with an unstable manifold (( = constant). These manifolds as shown in the surface of section are remarkably smooth, although closer inspection reveals that their spacing is coordinates quite irregular. The transformation from the (X,v) to the ((J)coordinates is continuous, but not differentiable as a more detailed argument would show. A small change in (X,rr) can lead to a large change in ((,q) in an almost unpredictable manner. In a somewhat simplified language, if the coordinates ( X , w in the surface of section are changed by 6, then the corresponding values of ((,v) change by e where 6 = # I 2 . Trying to interpolate the initial conditions (X,U) for a trajectory with given values (67) becomes very difficult under these circumstances. The main theme of this report is the following idea: the physical quantities of importance in the connection between classical and quantum mechanics are directly related to the coding of the trajectories rather than to the somewhat arbitrary coordinates in the surface of section. This principle will become clearer once the various uses of Feynman’s path integral have been described in the next section. 5. The Sum over All Orbits Investigating the idiosyncracies of a chaotic dynamical system is a praiseworthy endeavor, but I have tried for some 20 years to go beyond the problems in classical mechanics, in order to establish the connection with quantum mechanics. The primary tool in this enterprise is the path integral (cf. Schulmad) and its approximation by the classical trajectories, including the quadratic fluctuations of the admissible paths around these trajectories. Nobody has been able to come up with a significant improvement ( 6 ) Schulman, L. S . Techniques and Applications of Path Integration; Wiley: New York, 1981.
Gutzwiller of this simple recipe, except in one degree of freedom where there is no difference between regular and chaotic behavior. We have no error estimates for this particular approximation in quantum mechanics. On the other hand, even this simple general recipe leads quickly to some rather difficult mathematics, and it is important to keep the grand perspective from classical to quantum mechanics in mind and not to be put off by some unfamiliar technical problems. The path-integral formulation of quantum mechanics covers all the issues that come up in physics and chemistry, but I want to concentrate immediately on the spectrum of bound states. The search for the relevant wave functions is avoided because it involves technical difficulties in the classical approximation to the path integral which have not received a thorough examination thus far. I will mention, however, at which point in the development these additional problems arise, and why they cause further headaches. The starting point in any path-integral formulation is not a solution of the homogeneous Schroedinger equation, time dependent or not, which is discussed in most textbooks, but rather the inhomogeneous Schroedinger equation. Its solution describes the probability amplitude, F(q’:q:t) or G(q”,q’,E),for the system to go from position q’to position q”in the allotted time t or with the given energy E. Feynman’s path integral applies to F , not to G; the reason is that Huygen’s principle applies to the time domain but most emphatically not to the frequency (energy/ h ) domain. If we want to calculate G rather than F because it will yield the spectrum, then we have to stretch Feynman’s original idea from the paths in position space to the paths in phase space, where they are not nearly a: convincing. Nevertheless, when the classical approximation G to G has been made including the quadratic fluctuations around_the classical trajectories, G does not look much different from F, the corresponding approximation to F. -We seal1 not discuss F and G any further but will proceed to F and G. The tilde above F and G will be omitted henceforth. At this point we need the classical action R(q’:q:t), which corresponds to the propagator F(q’:q’,t), and the classical action S(q”,q’,E) which corresponds to Green’s function G(q”,q’,E)
where we have used the Lagrangian L(q,q) of the system, and the integral is taken along a classical trajectory from q” to q’in the time t . In close analogy one has
S(q”,q’,E) =
JJ,pil
(14)
where ptj is to be interpreted as a scalar product between the vestors p and q, and the classical trajectory from q”to q’is defined for the given energy E. There may be many trajectories from q’to q”in the allotted time t or at the given energy E. If the same trajectory is applicable to both (13) and (14), then one has the relation S(q”,q’,E) = R(q’:q’,t) + Et. The classical approximations, F for the propagator and G for Green’s function, are given for two degrees of freedom by the simple expressions F(q”,q’,t) =
1
-C C(q‘:q:t)]’’* exp 27nh CI.tI.
(15)
where C(q”,q’,t) is the determinant (due to van Vleck’) of the mixed second derivatives with respect to q”and q’ of R(q”,q’,t)
where D = D(q”,q’,E) is a determinant of mixed second derivatives (7) van Vleck, J. H. Proc. Natl. Acad. Sci. U.S.A. 1928.
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3159
From Classical to Quantum Mechanics with Hard Chaos with respect to q”,q’, and E of S(q”,q’,E). The last term in the exponent of both (15) and (16) introduces a phase loss of ?r/2 for each conjugate point along the trajectory from q’ to q”. Neighboring trajectories starting in q’ and going toward q” intersect one another at a conjugate point (in the caustic); the determinants C in (1 5 ) and D in (16) become infinite there, like semiclassical wave functions at the boundary of the classically allowed region. The expressions (1 5 ) and (16) have been known for some time although their generality and applicability has not been sufficiently appreciated. Our purpose for quoting them here is to examine whether they have any chance of being numerically evaluated for a harshly chaotic system. Our contention is that this formidable task can be carried out if the foliation of phase space into stable and unstable manifolds is understood. As an illustration for this claim we will discuss the anisotropic Kepler problem (AKP), which seems to be the simplest by far of the few such systems known. If we renounce the knowledge of the approximate stationary wave functions and are satisfied with the values of the corresponding energy levels in a bound-state problem, then we can take the trace on both sides of (16), Le., put q’ = q ” = q and integrate over all of position space. If we had still the original Green’s function on the left, its trace g ( E ) would be the sum over 1/(E - E,) for all eigenstates j . The right-hand side of (16) was first derived by G ~ t z w i l l e rand ~ ~ yields the “trace formula“ g(E) =
1
i
E-E - Ej - --E h,.,.2 J
TO sinh ( a / 2 ) (17)
A periodic orbit (p.0.) is a solution of the equations of motion (6) with the property that not only initial and final position agree, q‘ = q”, but also initial and final momentum, p’and p”. The system is allowed to run more than once over the closed track, but the period To measures the time for only one round. The stability (Le., lack thereof) enters through the exponent a,which tells us that neighboring trajectories increase their distance by a factor exp(a) after one periodic orbit. (Formula 17 applies to the examples to be discussed here, but one could have cosh ( a / 2 ) instead of sinh ( a / 2 )in harsh chaos; cf. Tabors.)
6. The Reduction to a Grand Canonical Partition Function Over the periodic Orbits in (17) In order to evaluate the for the AKP, we have to enumerate Of them and then find the values of the normalized action integral 9 of (21, of the stability exponent a,and the number Of conjugate points p for each Of them. Since there is a one-to-one correspondence between trajectories and binary sequences, it is enough to enumerate all the periodic binary sequences, i.e., finite chains o f f 1 which repeat themselves over and over again. Of as defined The main job is to find the corresponding in (2). Getting the Of however, us Only a few Pangs Of the number P the number of bits in the chain. Nevertheless, these tasks have to be done with some care and will be relegated to the next section. If one is interested in a particular chain of binaries (a,,a2,...,an), it is easy to work out the corresponding values of [ and q using condition ‘i+n = If then One ( l l ) and (12) with the looks at the grid in Figure 2, one can locate approximately the initial conditions (X,U) for the periodic orbit. It would seem easy to refine these initial conditions in an effort to close the trajectory smoothly. Due to the instability of all the trajectories, however, this task turns out to be extremely delicate (G~tzwiller~). It was Of carried Out by the author a”ost a decade ago with the Gregory Baran?but Only to = lo, using in many instances to get Some convergence’ A completely different approach will be discussed now, which is designed to get away from actually finding the initial conditions, but rather to aim action for the which is the quantity of interest in the trace formula (17). (8) Tabor, M . Physica D (Amsrerdam) 1983, 6D, 195-210.
The beginning of section 4 explains how one can find systematically, by varying only one parameter, all the trajectories that come out of a collision and end in a collision no later than at the 1lth intersection with the x axis. While integrating the equations of motion, one can collect the values of the action integral (2) between any two consecutive intersections. Let us call any trajectory between two consecutive intersections an arc and the corresponding value of 9 its length. Let us further limit ourselves to those arcs whose trajectory ends in a collision no further away than at the fifth intersection either before or after the arc. These arcs form the bulk of the collision trajectories which were described earlier. According to (1 1) and (12), the values of [ and q for these arcs form a grid of mesh size 1/32 inside the square -1 C (, q C + l . The arc length is a relatively smooth function which is defined on those grid points. This function is extended to the boundary of the square by setting it equal to zero, its natural limiting value. Let us represent this function in a slightly different manner by considering only the centers of each mesh and attaching to it the mean of the four values at the corners. These centers are defined by coordinates (&q) whose binary expansion has exactly six bits, f l , no more and no less. In this fashion we end up with an arc length for the center arc of a chain which consists of exactly 12 bits. We could designate the result of this little manipulation as the P(a,,a2,...,a12)and forget completely that we had to integrate various equations of motion for different values of the collision parameter A . The function \k depends only on the 12 binary variables, as well as the assumed mass ratio; its values can be obtained by representing \k as a polynomial in these 12 variables. There is no power higher than 1 in any one of them because one has always a’ = 1. The coefficients in such an expansion are not hard to find as is explained elsewhere (cf. Gutzwiller9), and there are exactly as many such coefficients as known values of P.The interest in such an expansion lies in finding out how complicated a function P really is. The symmetry of P,due to the time invariance of the equations of motion, implies that the coefficients vanish for the products of an odd number of a’s. P is always positive; the constant term in the polynomial expansion simply gives the average value of 9.The remaining terms are like the energy ~. of a chain consisting of 12 coupled spins. The calculation for the mass ratio 5 shows that this is not a simple chain; there are substantial terms which involve the product of four spins, in addition to the expected energies for the pairwise interactions. But now a remarkable thing happens when we calculate the lengths 9 for the periodic orbits with the help of the function \k. Their lengths are the of the arc lengths of 12 individual arcs, each belonging to a spin chain which is a cyclic permutation of the others, Therefore, the function \k has to be summed over cyclic permutations when it is used in the trace formula (17). This new function, call it 9 from now on, has the remarkable property which was noted in the earlier papers, that four-spin interactions and higher terms are quite small compared to the pairwise interactions. M ~the latter~ have an~ almost ~ purely exponential decay with the distance between the individuals in the pair. The numbers are instructive: The average length of a periodic orbit is 26,5054 in this calculation, while the maximum is 34,4563 as obtained in the earlier work for the repetition of the orbit (+-), and the minimum is for the of the orbit (+++...), are -1,061 1, -o.5782, -o.2989, The pairwise coupling to fifth neighbors, with -o.1478, and -o,0687 for the consecutive ratios 0.5449, 0.5169,0.4975, and 0.4648; the sixth neighbor interaction is much smaller at -0.0149; but it can be said to play an exceptional role in a chain of 12. The expansion of the function @ as a power series in the 12 spin variables al, ..., a12has all the properties of an expansion in a complete set of orthogonal polynomials. In particular, the average of the quare of @ equals the sum of the squares of the expansion coefficients, which turns out to be 721.5317. Most of thisnumber comes from the square of the constant term; the 6 two-body terms (with
~
3160
The Journal of Physical Chemistry, Vol. 92, No. I 1 1988 ~
multiplicities 12 or 6) account for 18.9161, while the remaining 1981 four-body and higher terms make up the remaining 0.0790. If we think of the latter as generating some noise on top of the pairwise interactions, its amplitude will amount to (0.0790)’/* = 0.271 1, while the constant term together with the pairwise interactions vary from a minimum of 0.7380 to a maximum of 35.0269. These last two values are too high for two reasons: the higher order terms are missing, and our calculation on the basis of collision orbits yields higher numbers than the earlier calculation which started directly from the periodic orbits. From all these considerations, one arrives finally at a simple expression for the action integral of a periodic orbit with the binary sequence a , , ..., a,: 9 = n(A
+ B + C + ...) -
n
+ B u , , ~+
+ ...)
1=1
(18) where one has the closure condition a,+,, = a,. The values of A , B, C, etc. are positive, and the constant term in (18) was chosen so as to make Q, = 0 for the ferromagnetic chain, a, = 1. Both the previous and the present calculations represent the action integrals for the periodic orbits in terms of spin variables as if one was dealing with a one-dimensional Ising model. But the previous work started from a knowledge of the action integrals of all the periodic orbits of length 10, whereas in this work the starting point is the collection of collision trajectories with no more than 10 intersections between the two collisions. The latter method is easy to automatize, but its results are not as good if one compares with the results of the former method using the same number (10) of intersections with the surface of section. In particular, the approximation with only pairwise interactions between the spins is even better when one uses directly the periodic orbits. In addition to a feasible algorithm, however, the new method provides something close to an argument why the Ising-like expression with exponentially decaying pair interaction gives such a good representation for the action integral of the periodic orbits. That remarkable feature seems to go back to the relatively smooth relation between the coordinates in the surface of section and the natural coding of the trajectories.
7. Stability Exponents and Symmetries in the Spectrum The stability exponents a for the periodic orbits can be obtained from the collision trajectories by a similar procedure of adding up the contributions from individual arcs. This task has not been carried out as yet; instead, the results from the earlier work will be adopted. The values of a do not vary at all as the action integrals of the corresponding periodic orbits; some of the short ones are very unstable, and the longest is also the most stable. The expansion of a into polynomials of the binary variables has none of the simple properties that were explained in the preceding paragraphs. There seems to be a great difference between the AKP and the motion on a hyperbolic manifold, where a equals the geometric length divided by the radius of curvature. The stability exponent in a harshly chaotic system is subject to a very general limitation which will be used to fix its value. In the normalization (1) of the Hamiltonian, a / @is the Lyapounoff exponent X, whose average is also called the metric entropy because it measures the rate at which phase space gets distorted in the neighborhood of a trajectory. A rather different quantity is the topological entropy T, which measures the proliferation of periodic orbits as a function of their length 9; their number N ( 9 < L ) = exp(rL)/L for large values of L. It can be shown that the average value of X equals T under rather general assumptions for the dynamical system. The statistics of the lengths 9 for periodic orbits, even according to the simple formula (1 8) in terms of one-dimensional spin chains, is not trivial, because one has to take into account spin chains of different lengths. But if one ignores this complication, and takes only the chains consisting of n binaries, their total number is 2n, and the topological entropy 7 = n log ( 2 ) / 9 . If this quantity is set equal to the metric entropy X = a/@,one finds that a = n log (2). The calculation of periodic orbits for n = 10 gave an
Gutzwiller average of a equal to 7.36, whereas 10 log (2) = 6.93. The agreement is tantalizing, and we will simply give all periodic orbits the same value a = n log (2). The number of conjugate points fl is basically simple, but it depends on the discrete symmetries of the eigenfunctions. The details are explained elsewhere (GutzwillerlOb);only the results will be quoted. When the Green’s function is made symmetric (or antisymmetric) with respect to reflexions on t h e y axis, and the trace is then taken, a new kind of periodic orbit appears in addition to the usual ones. Instead of ending in the same point of the surface of section, (Xn,Un)= (X,,,Uo), one now has (Xn,Un) = (-Xo,-Uo). The binary sequence (a,, ..., a,) is, therefore, followed (and preceded) by its negative (-al, ..., -a,), alternating back and forth instead of repeating itself. In applying (1 8j one has the closure condition ai+, = --al. The trace g ( E ) for the symmetric states with respect to reflexion on the y axis, called X-even, equals half the sum of (17) for the two kinds of periodic orbits, whereas the trace g ( E ) for the antisymmetric states, X-odd, is half the difference. Notice that n has to be an even number so far, because the trajectory closes only when the momentum in the y direction is the same at the beginning as at the end of the periodic orbit, even if (X,U) has changed sign. The AKP is also symmetric with respect to the x axis, and its Green’s function can be symmetrized (or antisymmetrized) accordingly. When taking the trace, another kind of closed orbit appears; it has (Xn,U,) = (Xo,Uo)in the surface of section, but the sign of the momentum u in the y direction is reversed. The number n is odd, but the sequence of binaries repeats itself without change of signs. Again, the trace for the symmetric states with respect to reflection on the x axis, called Y-even, is half the sum for these two kinds of periodic orbits, whereas the trace for the antisymmetric states, Y-odd, is half the difference. There are altogether four different combinations of sums (17) over the different periodic orbits. But now the original number of degrees of freedom enters into the picture. If the donor impurity were investigated in a two-dimensional crystal, there would simply be as many conjugate points as there are intersections with the surface of section. Also, all four combinations of symmetries with respect to the x axis and the y axis would be realized, yielding four different spectra. But in the three-dimensional crystal, our two-dimensional model applies only to the states which have cylindrical symmetry around the (longitudinal) x axis; the Y-odd states, whether X-even or X-odd, are excluded. Moreover, to each periodic orbit, whatever its kind, there corresponds a one-parameter family which is obtained by rotation around the x axis. Each crossing of the x axis yields an additional conjugate point so that altogether one has two conjugate points for each intersection with the surface of section. Henceforth, the trace will always be defined for a particular combination of symmetries with respect to X and Y and for a fixed dimension of the original crystal. The length n of the periodic orbits will, therefore, vary over all integers, odd as well as even. Also, both closures of the binary sequence, with same or with opposite sign, will occur. The reader will have to check whether the rules have been applied correctly. Only rotationally invariant states in three dimensions will be discussed in detail; since they are Y-even, the periodic orbits with n odd will come in with the same sign as the periodic orbits with n even. But since a periodic orbit in three dimension has 2n conjugate points, each causing , terms from the periodic orbits with odd a phase loss of ~ / 2 the n have an overall negative sign, while the terms from periodic orbits with even n are positive. That still leaves the two possible kinds of closure for the binary sequences. 8. Energy Levels as Zeros of a Grand Partition Function The left-hand side of (17) can be integrated with respect to E , starting at -m, and then exponentiated. Leaving out a constant (9) Gutzwiller, M. C . Classical Dynamics and Dynamical Systems; Devaney, R. L., Nitecki, 2. H., Eds.; Lecture Notes in Pure and Applied Mathematics 70; Marcel Dekker: New York and Basel, 1981; pp 69-90. (10) (a) Gutzwiller, M. C. Phys. Rev. Lett. 1980,45, 150-153; (b) Physica D (Amsterdam) 1982, 5D,182-207.
From Classical to Quantum Mechanics with Hard Chaos factor, one obtains formally the infinite product of ( E - Ej) over all energy levels Ey On the right-hand side, one notices that the period To= aSo/aE and that the stability index cy is independent of E . The integration with respect to E can, therefore, be carried out, leaving a factor I / m in each term because S ( E ) = mSo(E). So far the sum on the right-hand side of (1 7) goes over all periodic orbits, including retracks of primitive periodic orbits. The independent variable in the sum is now switched from the periodic orbits to the binary variables of finite closed spin chains. In doing so, each primitive periodic orbit is overcounted n times, where n is the number of intersections with the surface of section. A further factor l / n is needed when the new variables are used in the summation, and this factor combines with the earlier 1/ m to a common l / n without any distinction between primitive and nonprimitive periodic orbits. Moreover, it is helpful to replace 1/2 Sinh (a/2) by exp(-cy/2) = 1/21i2. To be more conscientious, one should expand Sinh in the denominator; but the additional terms would only generate additional energy levels which are broadened by large imaginary parts. Up to an additive constant, the trace formula (17) becomes
The Journal of Physical Chemistry, Vol. 92, No. 1 1 . 1988 3161
1
This matrix, as well as all the ones with 1 > 2, can be broken up into two diagonal blocks of equal size by a similarlity transform with a matrix SI which has only entries 0, + 1, and -1. As a simple example one can write 1
0 0
0 1 0 -ab
0
1
1 0 - 1 0
-a
(23)
This reduction of the transfer matrices is possible because the matrix element (20) remains the same when all binaries change their sign simultaneously. The two reduced blocks for thirdneighbor interaction show the characteristic structure which is found in all the higher order reduced transfer matrices 1
c
abc
ab
o
o
c
-abc
o
-ab
where @ is given by (18). The X-even spectrum results from adding the sums on the right for the even closure a,+, = a, and the odd closure a,+, = -a,, while the X-odd spectrum results from the difference of these two sums. Since E < 0, the variable z is purely imaginary; but if one adds a small positive imaginary part to E , then z acquires a small positive real part. The right-hand side of (19) is, therefore, convergent sufficiently far up above the real axis of the complex E plane. The sum on the right-hand side of (19) is certainly not absolutely convergent on the real E axis. On one hand, a crude estimate says that the number of terms for a given n is 2", and the absolute value of each term is 2"/'.On the other hand, one does not want any singularities above the real E axis; the only salvation seems to be the destructive interference between the phases -z@ at a fixed value of n. If they are distributed randomly , sum over the spin chains of a given length between 0 and 2 ~the n has a most probable absolute value of 2"/', exactly what it takes to offset the denominator which comes from the stability exponent. This primitive argument gives us a fighting chance to get a reasonable spectrum out of the great proliferation of periodic orbits in hard chaos. Let us compute the right-hand side of (19) of @ as in (18), with pairwise interactions out to the lth neighbor. For this purpose, we write first the so-called transfer matrix, starting with the expression
Notice the occurrence of all possible combinations of a , b, and c, as well as the opposite signs in the lower halves of the two matrices. The idea of the transfer matrix is quite simple, but it takes a little while to get familiar with it. Each multiplication of two transfer matrices takes care of summing over one binary variable. The size of the transfer matrix depends only on the range of the interaction I ; a four-spin interaction, -Ma,+,a,+2a,+lal,can be accommodated in the matrices (24) by multiplying half of the nonvanishing matrix elements with exp(-2zM). When (n - I) transfer matrices have been multiplied, there are two ways of closing the spin chain. The even closure condition, a,+, = a,, does not require anything further to be inserted before taking the trace of the nth power of TI.The odd closure condition, a,+, = -al, is enforced by inverting all the spins in the range of interaction at some point before taking the trace of the nth power of TI.This inversion matrix has 1's in the upper right to lower left diagonal, and 0's elsewhere. It becomes regularly diagonal in the similarity transformation by SI, as in (23), with +1 in the upper half and -1 in the lower half. When calculating the trace for the X-even spectrum by adding the contributions from the two kinds of closure, one ends up with only the trace of the upper half of the reduced transfer matrix. The X-odd spectrum, on the other hand, requires only the trace of the nth power of the lower half. These results can now be inserted into (19). The sum over the exp[-zA(1 - a,+,a,) - zB(1 - a l + l u l ~...l )- zL(1 - a l + l a l ~ ~ + l ) ] set of binaries ( a l , ..., a,) becomes the trace over the nth power (20) of the reduced T,; and the trace becomes the sum of the nth powers of the eigenvalues A, of the reduced TI.At this point the order This exponential is the matrix element in a 2' by 2' matrix, where of summation in (19) is inverted, with the sum over n done first the rows are labeled by the set of 1 binaries (aL+], a,, ..., a,-,+z), and the sum over m done last. The sum over n is the expansion and the columns by the set of binaries (al, a,-,, ..., a,-/+l). The of log (1 x), so that one has now a logarithm on both sides of author learned this way of defining the transfer matrix from T. (19). Up to a factor which is independent of the coupling constants S. Schultz (private communication). Each set of binaries is A , B, C, ..., one finds with k = 2I-l that enumerated in lexicographic order; Le., each set is interpreted as k a word whose alphabet consists of two letters, first a = +1 and ( E - E,) = (1 + 1i,/(2)'/~) second a = -1. Each matrix element is the product of the factors (25) ,=1 m=l a = exp(-2zA), b = exp(-2zB), c = exp(-2zC), .... For only nearest neighbors, Le., 1 = 1 and B = C = ... = 0, the The approximate energy levels E, are the zeros of the right-hand transfer matrix becomes simply side. Such a zero obviously requires that A, = -(2)]i2. This condition can be written more directly in terms of the reduced transfer matrix T Iin the form
+
m
n
For next-nearest neighbors, Le., I = 2 and C = ... = 0, the transfer matrix is
n
det IT, + 2'/*l = 0
(26)
Since the elements of the transfer matrix are exponentials of the
3162
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
variable z in (19), the last equation is transcendental in spite of its algebraic appearance. In order to be physically reasonable, the solution with respect to z should lie on the imaginary axis. If such is indeed the case, it would also vindicate the primitive reasoning concerning the destructive interference of the phases in (19).
9. Pseudospectra The bound states of an electron at a donor impurity in silicon with vanishing angular momentum around the longitudinal axis were calculated with the help of the trace formula some years ago (Gutzwiller’oa,b). Expression 18 was used where the pairwise interactions were decreasing exponentially with distance, and the range was infinite. The transfer matrix was infinite and had to be truncated judiciously. The results were excellent for the X-even spectrum, but some of the X-odd energy levels had imaginary parts which were a good fraction of the level spacings (cf. Wintgen et aI.’l). This question of the level broadening in the trace formula of a harshly chaotic system seems worth investigating in a more abstract setting. The transfer matrices of the preceding section allow us to construct spectra of the kind that appear in the donor impurity, but with a simpler dependence of the action integral on the binary code of the periodic orbits. The simplest case is the nearestneighbor coupling of (21). It should be emphasized that this approximation to the AKP is poor, because without any interaction beyond the nearest neighbor the mean value of the action integral for periodic orbits of a given length n is exactly half its maximum value, in contrast to the numerical data of section 7 . The solution of (26) with the reduced transfer matrix (21) is trivial and yields for integer k X-even: X-odd:
+ 1)ia - log (2Il2 + 1 ) 2zA = -2kia - log (2’1’ + 1)
2zA = -(2k
(27)
(11) Wintgen, D.; Marxer, M.; Briggs, J. S . J . Phys. A, in press.
55
I
I
I
I
I
I
-1.5
-1 0
-0 5
0.0
I
I
30 -
25
I
I
I
-2.5
-2 0
20 -3.0
05
Figure 3. The zero in the complex plane which corresponds to the ninth X-even energy level; 2A = 1, 2B increasing from 0 to 0.623, and 2C increasing from 0 to 0.389, and finally 20 increasing from 0 to 0.244.
(27)
A Balmer-like formula for E is obtained immediately when expression 19 for z is used to get E . In order to make a closer comparison with the hydrogen spectrum, we should choose the value of A to yield na for the “alternating” orbit with n intersections. This condition implies A = a / 2 and gives the right expression, probably for the wrong reason. More serious is the real part in 2zA whose value 0.880 is a sizable fraction of the separations 2a in the imaginary parts. The solution of (26) with the reduced transfer matrices of (23) already presents analytical difficulties; but it can be solved numerically if one starts with the solutions (27) for B = 0 and increases B to the required final value. The action integral for the alternating orbit with n intersections is still 2nA, while the mean action integral for all orbits with n intersections is n(A + B). Therefore, one can at least approximate the real distribution in its most obvious aspects. Figure 3 plots the values of z = x + iy for 2A = 1 and 2B going from 0 to cos ( a / 7 ) = 0.623 489 8, which corresponds approximately to the ratio of these parameters for the AKP with mass ratio 5 as in Si and is also not close to any rational number. The X-even level k = 9 in (27) was chosen to start the solution of (26) at z = 17a -log (2Il2 1). As 8 increases, z goes through an amazing set of gyrations in the complex plane, first staying roughly at the level of y = 5 5 7 , but then going down rather precipitously while swinging back and forth in its real part between 0 and -2.5. The upper limit 2 8 = 0.623 49 is reached at about z = 31i - 1 .O. Until then, however, z has never crossed the imaginary axis, although its real value has gotten occasionally all the way up to 0. This last statement can actually be demonstrated analytically by using ( 2 6 ) with the reduced transfer matrices (23). One finds that the real axis is reached for certain simple rational values of the ratio B / A . The solution of (26) in Figure 3 was continued by allowing for the third-neighbor interaction 2C in (24) to increase from 0 to 0.388 739 5. a factor 0.623 489 8 down again from 2B. The motion
+
Gutzwiller
of the solution z through its complex plane loses some of its systematic behavior, but it keeps on swinging back and forth in the x direction until it reaches z = 24i - 1.0. It crosses the imaginary axis, but only for a small excursion into positive values of its real part. There seems to be no change in this strange behavior, as a fourth-neighbor interaction 2 0 is allowed to increase from 0 to 0.242 375 1 on top of the already existing first, second and third neighbors, again a factor 0.623 489 8 below 2C. The same kind of behavior for the solution of (26) is found with other starting values of z , as well as for the X-odd energy levels, always for the three-dimensional AKP. There is no qualitative change in two dimensions, either. Since the number of conjugate points then equals the number of intersections with the surface of section, the phase loss is na/2 for a spin chain of length n. The condition (26) for the reduced transfer matrix TI becomes det (TIf i2Ii2l = 0
(28)
where the double sign refers to the two symmetries with respect to reflections on the x axis, Y-even and Y-odd. Again the solution is elementary for 1 = 1 and serves as starting point for the solutions of the more complicated equations for 1 > 1. The imaginary parts of z are predominantly negative as before. This last statement is in contradiction with the earlier crude agrument which showed that the zeros of the right-hand side in (25) should be purely imaginary, if the number of periodic orbits (topological entropy) is directly related to their stability (metric entropy) and the values of the action integral are randomly distributed. One would have thought that the model in terms of spin chains with rather long range pairwise interactions easily satisfies the last requirement. Now it seems, however, that the trace formula (17) is a much more subtle test of the statistical properties for the collection of periodic orbits in a harshly chaotic system. 10. Some General Conclusions Classical dynamical systems with hard chaos form a well-defined special class, somewhat analogous in their simplicity to the integrable systems, but much more numerous. The history of a trajectory in a harshly chaotic system can be described uniquely
J . Phys. Chem. 1988,92, 3163-3173 in terms of a simple coding scheme. For the anisotropic Kepler problem (AKP) in two dimensions, infinite sequences of binary numbers are the natural scheme. The characteristic features and the physical properties of such a dynamical system seem all to be contained in the code, not in the usual surface of section. This principle applies in particular when one tries to understand how a harshly chaotic system is connected to its quantum-mechanical ancestor. The classical limit of Feynman's path integral and, in particular, the trace formula appear to be the only general guides available in such an endeavor. This approach works beautifully for the free motion on a surface of constant negative curvature, where Selberg's trace formula (which states that the two sides of (1 7) are strictly equal) establishes a direct relation between the spectrum of the Laplacian and the collection of all periodic orbits (some useful introductions are quoted in ref 3). The mathematical proofs in this case of hard chaos, however, take advantage of the high degree of symmetry in the underlying space
3163
and cannot be generalized. Nor can anything as strong as Selberg's trace formula be valid for the many harshly chaotic systems. The two sides of (17) are nevertheless expected to be very close quantitatively, with the singularities with respect to E on the right-hand side very near the poles on the left-hand side. Any possible deviations have to be small compared to the distance between the poles. In particular, the imaginary parts of the singularities on the right-hand side have to be small compared to the distance between them in the complex plane. The example of a mathematical model for the AKP in terms of spin chains with pairwise interactions shows that this last requirement imposes restrictions on the dynamical system which are not understood. They have to do with some delicate statistical properties of the periodic orbits which do not seem to hold for some values of the parameters in the spin-chain model. This mysterious feature in a harshly chaotic system can be called its third entropy, which measures something beyond its toplogical and its metric entropy.
Dynamics for Nonconservative Systems: Ergodicity beyond the Microcanonical Ensemble Julius Jellinek Chemistry Division, Argonne National Laboratory, Argonne. Illinois 60439 (Received: August 1 7 , 1987; In Final Form: December 3, 1987)
A generalization of the procedure suggested recently by Nos6 for dynamical simulation of a canonical ensemble is presented. It is shown that there are infinitely many different dynamics capable of mimicking the canonical ensemble for a physical system provided the corresponding dynamics for an extended system consisting of the physical system and a "bath" is ergodic in a particular sense. An invariance property of ergodicity under scaling of the time is established. It is shown how the time scaling can be used to "tune the dynamics into ergodicity" with respect to any desired distribution function in a given phase space. The dynamics of the extended system implies an apparent trajectory for the physical system which is quite different from the one generated by the Hamiltonian dynamics for a closed physical system. A general unified scheme for dynamical simulation of any statistical mechanical ensemble is formulated.
I. Introduction: Background and Overview Use of computers as a tool of scientific inquiry resulted in a remarkable extension of our practical abilities to study detailed dynamical as well as statistical properties of physical systems. The systems treated vary in size from an aggregation of a few particles, such as clusters with free or constrained boundary conditions,' to systems with enough degrees of freedom to mimic bulk matter.2 The two most important computational techniques for this purpose are those of molecular dynamics3 (MD) and Monte Carlo4 (MC). The M D procedure follows the time evolution of a system while the M C method simulates its ensemble behavior. In an attempt to keep the results of these two approaches in harmony, the fundamental issue of ergodicity, the relation between the time-averaged and ensemble-averaged properties of a system, emerges as a practical problem. The fundamental question in the general theory of dynamical systems can be formulated as follows5 (we consider only classical systems): does the dynamical system of interest generate an invariant measure, Le., a distribution function in the phase space, and if it does, what is this measure? (1) See, e&; Jellinek, J.; Beck, T. L.; Berry, R. S. J . Chem. Phys. 1986, 84,2783. Davis, H. L.; Jellinek, J.; Berry, R. S. Ibid. 1987,86,6456. Beck, T. L.; Jellinek, J.; Berry, R. S. Ibid. 1987, 87, 545 and references therein. (2) See, e&; Nost, S.;Yonezawa, F. J . Chem. Phys. 1986,84, 1803 and references therein. (3) Kusick, J.; Berne, B. In Statistical Mechanics, Part E Berne, B., Ed.; Plenum: New York, 1977; Vol. 6, Chapter 2. (4) Valleau, J. P.;Torrie, G. M . In Staristical Mechanics, Part A ; Berne, B., Ed.; Plenum: New York, 1977, Vol. 5, Chapters 4 and 5. (5) Eckmann, J.-P.; Ruelle, D. Rev. Mod. Phys. 1985, 57, 617.
0022-3654/88/2092-3 163$01.50/0
It is of at least as much interest to flip the question and to ask what types of dynamics (more explicitly, what dynamical equations of motion, if any) are implied by a given statistical ensemble or distribution. The reason is that one would like, ultimately, to be able to select from all those dynamics yielding the desired equilibrium properties a particular dynamics that will represent in a simple way the nonequilibrium properties of the system in interaction with a specific, complicated environment. For example, if one replaces the real heat bath coupled to a system by a very simplified model, one would like to know how to select the most accurate representation of that bath, in order to compute not only the macroscopic equilibrium properties but also time-dependent properties such as time-correlation functions and response functions. A conceivable way to establish a correlation between a dynamics and a statistical ensemble is to capitalize on the ergodicity property, either proven or assumed. If a description is restricted to the conventional phase space I'(ij,g),every point of which represents a mechanical state of the system, only two of the traditional statistical mechanical ensembles, the microcanonical and the canonical, can be defined through their phase-space density distributions. Other ensembles can also be represented by distribution functions but in an augmented phase space obtained from the conventional "mechanical" phase space by adding additional quantities, e.g., ~ o l u m eas , ~dynamical variables. (6) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976. (7) Andersen, H. C. J . Chem. Phys. 1980, 72, 2384
0 1988 American Chemical Society
