Subspace-based non-blind deconvolution

Foreward

The present review has an appropriate place in the Volume to honor Prof. Giorgilli. In fact Antonio Giorgilli has played a major role in the theory of the third integral. I was very happy to have many collaborations with Antonio on many aspects of the theory over many years. Thus the present review is a tribute to his insightful contributions during these years.

1.   Introduction

It is well known that a stationary axisymmetric galaxy has two integrals of motion, the energy and the angular momentum along the symmetry axis. In some exceptional cases it has also a third integral of motion. These are called integrable cases. Such cases were found by Stäckel ^[116,117] and later by several people. A systematic search of integrable cases applicable to galactic dynamics was made by Lynden-Bell ^[100]. Further detailed searches were made by Hietarinta ^[88], Ramani et al. ^[112] and Lakshmanan and Sahadevan ^[93]. For a review see ^[38].

However in general dynamical systems are nonintegrable. A proof of this statement was given by Poincaré in his Méthodes Nouvelles ^[111]. Generalizations of Poincaré's theorem were provided by Fermi ^[73,74]. A most general proof that in general dynamical systems are nonintegrable was provided by Siegel ^[114]. In fact there is an infinity of integrable systems, but they are exceptional in the same way as the rational numbers are exceptional with respect to the irrational numbers.

My own involvement in this problem started during my visit in Stockholm in 1958. At the suggestion of the leading Swedish Astronomer Prof. Bertil Lindblad. I extended his work of first order epicyclic orbits on the plane of a galaxy to all higher orders ^[19]. Then I tried to extend my calculations to the third dimension. Very little work existed at that time on the vertical oscillations. Thus I calculated on the computer BESK of the University of Stockholm two orbits on the meridian plane of a very simple model of the galactic Hamiltonian:

representing the motions of stars near the Sun.

The results were surprising (Figure 1). I expected that the orbits should fill the whole space inside the Curve of Zero Velocity (CZV), but instead they filled two different "curvilinear parallelograms", in a way similar to Lissajous figures. I published these orbits in 1958 ^[20] and I presented my results at a special meeting of the IAU General Assembly in Moscow (1958). At the meeting there was Dr. G. Kuzmin, who commented that my calculations indicated the existence of a third integral in this galactic model. I replied that the existence of a third integral was not probable for two reasons: (a) the theoretical nonexistence theorems of Poincaré and Fermi and (b) because the system I considered was so simple that if there was a third integral someone would have found it already. In fact, later I realized that in fact Whittaker ^[126,127], Cherry ^[15,16,17] and Birkhoff ^[7] had found formal integrals in similar cases.

After the Moscow meeting Dr. Kuzmin wrote me on 7-2-59: "I think that the Poincaré theorem is applicable to the third integral problem and that it proves the non-existence of this integral in general case... Although these considerations are perhaps not fully free from gaps they dispel to a great extent my doubts about the non-existence of the one-valued third integral in general case...".

But meanwhile I could calculate a formal third integral in the form of a series and I presented my results at the Kiel meeting of the Deutsche Astronomische Gesellschaft in 1959. Thus I replied to Dr. Kuzmin on 29-9-59: "Meanwhile I have found a third integral in a special case, and then in a quite general case. I made a short communication about this problem in the Kiel Meeting of the Deutsche Astronomische Gesellschaft. I include herewith a short summary of this communication. I will send you a complete copy of my paper as soon as I have one available".

In fact my paper appeared in the Zeitschrift für Astrophysik in 1960 ^[21] and Dr. Kuzmin wrote me in 15-7-60: "Recently I found in ZAp your very interesting paper about which you spoke in your letter. Your paper inspired me to return to the problem at last".

The third integral of the Hamiltonian (1) is a series of the form:

where

is the energy along the x-axis (in the lowest order). The third order term is:

and higher order terms can be calculated step by step.

At the Kiel meeting there was present Prof. Otto Heckmann who was very much interested in my work and encouraged me to proceed further. In particular he wanted me to prove the convergence of the third integral. However I could not find any such proof and I left the problem open.

Later I realized that in general the third integral is not convergent, but only formal. In fact if we truncate the third integral at successively higher orders we find first a better and better convergence, but beyond a certain order the values of the truncated third integral have larger and larger deviations (Figure 2) ^[51]. As $x$ increases the best truncation occurs at smaller orders. This is consistent with the Nekhoroshev ^[108] theorem for formal series.

Nevertheless formal series may give very useful approximate results. An example was provided by Poincaré ^[111] who considered the series.

This series is divergent, but if we truncate it at high orders smaller than $n = 1000$, we find very constant results. On the other hand the series

is convergent but requires orders higher than $n = 1000$ in order to show its convergent character.

Thus Poincaré states that while the mathematicians consider (correctly) the first series is divergent and the second convergent, the astronomers on the contrary consider the first series convergent and the second divergent, because their numerical calculations never go to very high orders.

One of the first astronomical applications of the third integral was in finding the forms of the velocity ellipsoid in our Galaxy. It was observed that the velocities of the stars near the Sun are distributed around the velocity of the Sun with decreasing densities outwards, the equidensities forming three-axial ellipsoids. If there were only two integrals of motion, the energy $E = \frac{1}{2}(\dot{r}^2 + \Theta^2 + \dot{z}^2) + V(r, z)$ and the z-component of the angular momentum $A = r \Theta$ one would expect a distribution function

that would have the $\dot{r}$ and $\dot{z}$ axes equal. However the observed distribution gives the $\dot{z}$-axis much smaller than the $\dot{r}$-axis. This can be explained if a third integral $\Phi = \frac{1}{2}\dot{z}^2 +...$ exists beyond the two other integrals. Then the $r$ and $z$ axes do not need to be equal. (This was the thesis of B. Barbanis ^[5]).

Much work on the third integral had been done by Whittaker ^[126,127] who called this integral "adelphic" (from the Greek word "adelphos" = brother) and also by Cherry ^[15,16,17] and Birkhoff ^[7]. Whittaker used action-angle variables while Cherry and Birkhoff used complex variables. My own work was using cartesian variables that are very useful for galactic applications. Furthermore I used a computer program to calculate higher order terms of the third integral, that gave very good approximate results (^[24,43]).

Another important development around that time was the celebrated KAM theorem ^{[1,2,3,92,105,106]}. This theorem states that in an autonomous Hamiltonian system, close to an integrable system expressed in action-angle variables, there are N-dimensional invariant tori containing quasi-periodic motions with frequencies $\omega_i$, satisfying a diophantine condition (after the ancient Greek mathematician Diophantos):

where

with $k_1, k_2, \dots, k_N$ integers and $|k| = |k_1|+|k_2|+\dots |k_N|\neq 0$, while $\gamma$ is a small quantity depending on the perturbation. Thus if a system is close to integrable and the frequencies $\omega_i$ are far from all resonances the orbits lie on invariant tori.

On a surface of section of a 2-D system these tori form invariant curves around a stable invariant point which is the intersection of a stable periodic orbit.

It is remarkable that Birkhoff ^[7] did not believe in the existence of such invariant curves and stated that orbits starting close to an invariant point of the stable type would gradually go very far from it.

The KAM theorem had many important applications in various dynamical systems. It was later supplemented by the Nekhoroshev theorem ^[108], that used expansions of the third integral type to approximate the motions in non-integrable systems close to an integrable system. Namely if the deviation from an integrable system is of order $\varepsilon$, the formal integrals are valid over an exponentially long time of the order of

where $t_*, M$ and $\tau$ are positive constants. This time in particular astronomical problems may exceed the age of the universe.

Improvements of the Nekhorosev theory were provided by Giorgilli ^[77,78] and in particular by Morbidelli and Giorgilli ^[102] who found superexponential stability times. In the limiting cases of the KAM tori this time goes to infinity ^[81]. Therefore the KAM theorem is a limit of the Nekhorosev theorem. Further developments of the Nekhorosev theory were provided by Contopoulos, Efthymiopoulos and Giorgilli in ^[51] and by Efthymiopoulos, Giorgilli and Contopoulos in ^[68].

Around 1962 I had some interesting interactions with J. Moser. At first Moser did not believe that asymptotic series like that of the third integral were useful. But he asked me to give a seminar on this subject at the Courant Institute in New York in 1962. When he saw how well the third integral was conserved (I had at that time computer calculations up to the 10th order) he was impressed. Later he wrote in a paper ^[107] "These questions received renewed interest on account of the work of Contopoulos, who searched for and calculated a "third integral" for a system describing a galactic model".

When I was a visiting professor at the Yale University (1962) I gave a series of lectures on the third integral. I summarized my lectures in a paper ^[22] under the title "On the existence of a third integral of motion". I had included in this paper a section devoted to the use of a "surface of section" to visualize the orbits in the phase space of a conservative system of 2 degrees of freedom. My paper was accepted and I had received the proofs, when the referee sent me a belated remark, that the surface of section part of my paper was not so important and it would be good to omit it. Thus I took out this part. That was a mistake because soon afterwards there was the famous paper by Hénon and Heiles ^[87] under the title "The Applicability of the Third Integral of Motion: Some Numerical Experiments" that was based mainly on figures drawn by the computer on a surface of section.

During my stay at Yale I organized a meeting on dynamical systems that was attended by Moser, Hénon, Toomre, Szebehely, Brouwer etc. Moser spoke about the KAM theorem. But when I asked him if he could give an estimate of the applicability of his theorem to a practical problem like the Sun-Earth-Moon system, he refused to give such an estimate.

Then Hénon made some calculations in his hotel and next day he gave an estimate for the convergence of the KAM theorem that was published later ^[85]. This estimate was $10^{-48}$ in some appropriate units, which would imply that the orbit of the Moon would be stable despite the perturbation due to the Sun, if the Earth and the Moon were point masses at a distance less than 1mm. Later, however, it was shown that the KAM theorem was giving realistic results in many real cases. As a consequence people could apply the third integral to prove the stability of the Trojan asteroids at distances from the equilateral equilibria $L_4$ and $L_5$ of the order of 1 astronomical unit ^[64,67].

During my stay at Yale I had a Greek student, Mr. G. Bozis, and I suggested to him to study numerically the restricted three-body problem (RTBP) to find whether a third integral was applicable in this case. Then I left Yale. Bozis asked permission from his supervisor, Prof. V. Szebehely to use the computer to calculate orbits in the restricted problem. But Prof. Szebehely was reluctant. He told him that it was known that Poincaré's theorem is applicable in the RTBP, thus there should be no third integral in this case. Bozis, however insisted and finally Szebehely reluctantly agreed to let Bozis use the computer for these calculations, but he told him "You do that at your own risk". But when Bozis got his first results Szebehely was amazed. Bozis' figures indicated a well established third integral.

After my return to Greece (1964) I organized an IAU Symposium in Thessaloniki on "the theory of orbits", that was attended by many people working in Celestial Mechanics and Stellar Dynamics. Among them there were Hénon and Szebehely. Hénon spoke about the orbits in the restricted three-body problem and showed many figures indicating the existence of a third integral ^[86]. Then Szebehely jumped up at the discussion period saying that he had done similar work with G. Bozis, a student of Prof. Contopoulos. He presented many figures, making essentially a second lecture on the subject ^[118]. Bozis published his results in detail later ^[11].

2.   Resonant cases

The third integral was very useful in calculating the usual types of orbits in a Galaxy ("boxes", i.e., deformed Lissajous figures, like those in Figure 1 or more deformed as in Figure 3) The agreement between the theoretical results and the numerically calculated orbits was perfect. But Torgard and Ollongren ^[119] had found also "tube orbits" i.e., orbits filling tubes around various stable periodic orbits. While I was at Yale (1962) I found also several types of tube orbits (Figure 4) using the IBM 707 computer of the Institute of Space Sciences in New York.

At first people thought that the third integral was not applicable to such tube orbits. However at that time I developed the theory of the resonant forms of the third integral and I concluded that the tube orbits were due to particular resonances.

After my return to Greece (1963) I continued using the computers of the Institute of Space Studies by mail. I was sending the initial conditions of the orbits to a friend of mine at the Institute of Space Studies, in the form of punched cards. After about one month I was receiving the output, mainly in the form of figures, done by the computer. (Today such calculations are done in a few seconds).

Most orbits in the meridian plane of a galaxy were boxes. But I found that some initial conditions should lead to particular tube orbits. I sent these initial conditions to my friend in New York. When the results arrived, in the form of a thick roll of computer figures, the whole Astronomy Department was excited. We were looking with anxiety at the figures, as we developed the roll on the corridor. And behold! There appeared the nice tube orbits that we expected. And we saw also the expected transitions from the tube to the box orbits. These results were published later ^[23,25].

The calculation of the higher order terms of the third integral cannot be done by hand, as there are many thousands or millions of terms to be calculated. Thus I decided to write a computer program to calculate these higher order terms. At that time there were no packages of computer algebra, like Mathematica etc. Thus I wrote a program in Fortran that could calculate the third integral up to order 12. We had an IBM 1620 computer at the University of Thessaloniki, and we used punched cards. Because of limitations in memory I had to split my program in two parts. The first part produced a deck of cards 1 meter long. This was used as an initial condition for the second part, that produced a deck 2 meters long, and this finally was used to calculate our orbits, periodic orbits etc. The results were impressive. We published these results in 1965 ^[43] and in more detail in ^[24]. A little later Gustavson ^[82] prepared a similar paper, also finding higher order terms of the normal form (equivalent to the third integral).

When $\omega_1 = \omega_2$ the form of the third integral changes completely. It is not a part of the energy, as in non resonant cases. The forms of the orbits in this case are very different from distorted Lissajous figures (Figure 5) ^[25,43,44].

I was in Chicago in 1969 when I saw Dr. Chandrasekhar in the corridor carrying a gadget to show to his students. This contained a sphere that could go up and down, or rotate around its vertical axis. Chandrasekhar wanted to show the coupling between the up and down oscillation and the rotational oscillation. If one would set the sphere in the up and down mode, later this up and down motion would gradually stop and all the energy would go in a rotational oscillation clockwise and counterclockwise. Even later the rotation would stop and the energy would return to the up and down motion and so on.

The use of an "experiment" by Chandrasekhar was quite unusual. Then Chandra told me "You see, George, a case where your third integral is not applicable". But, as I explained to him, this was a perfect case of a 1:1 resonance, where the third integral has a very different form.

One important application of the resonant forms of the third integral is in finding periodic orbits. Poincaré has indicated that the periodic orbits are dense in generic dynamical systems, both integrable and nonintegrable.

A periodic orbit in an autonomous Hamiltonian

has two characteristic exponents zero and a pair of exponents $\pm \alpha$ opposite to each other. If there is another integral of motion

then either $\alpha = 0$ or

(where $\Phi_x = \partial{\Phi}/\partial{x}$ etc) which means that the surfaces $H$ and $\Phi$ are tangent along the periodic orbit. We found many cases where the latter relation is satisfied. However for every resonance the form $\Phi$ of the third integral is different.

If we replace $\dot{y}$ from Eq. (2.1) into Eq. (2.2) we have on a surface of section (say $y = 0$), for a given energy $h$, an invariant curve

The periodic orbits are singular points on this surface of section

We found several resonant orbits in this way using our computer program to calculate the third integral ^[24]. In some cases we needed to use high order truncations in order to find periodic orbits with sufficient accuracy ^[27].

However in order to find the third integral up to very high order an effective computer program was necessary. This was done in a way that was quite unexpected. I was invited by two Italian Professors, Scotti and Galgani, to the Cagliari meeting of the Italian Physical Society (1972). I arrived an evening at Cagliari, and Scotti and Galgani were waiting for me at the airport. We went to my hotel and discussed many problems related to the third integral. The discussion continued the next day before and after the meeting of the Physical Society and then they drove me to the airport in the afternoon. I did not see Cagliari or Sardinia at all. But our discussion had one important consequence. When Galgani returned to Milan, he found a student (Giorgilli) who was eager to follow up this problem. And Giorgilli produced the best computer program to calculate the third integral to very high orders ^[76,80]. Giorgilli used a computer program that calculates the "generating function", the normal form and the integrals as series in the original variables. This method was a variant of the Lie series method of Hori ^[89] and Deprit ^[63].

The calculation of high orders was necessary because for low orders the third integral looked exact. E.g., Kaluza and Robnik ^[90] calculated the new integral in the Hénon-Heiles ^[87] problem up to order 14 and they found that the system looked integrable. Only if one goes at least to order 40 one sees definitely that the system is not integrable.

Later Efthymiopoulos calculated also high order integrals and applied them to the Trojan asteroids around the Lagrangian points $L_4$, $L_5$ of the solar system ^[64].

The periodic orbits for relatively small perturbations follow a regular pattern (Figure 6) for the Hamiltonian

with $\epsilon = 4.0$ ($\omega_1^2 = 1.6, \omega_2^2 = 0.9, h = 0.00765$).

Every periodic and nonperiodic orbit has a particular rotation number and on the left the successive rotation numbers seem to vary smoothly in general. Only near $2/3$ on the left of Figure 6 the decrease is abrupt (indicating chaos) and on the right there is an island of stability, where the rotation number is constant. In fact near every stable orbit there is a small plateau and the curve rot as a function of $\overline{x} = \omega_1x$ is a "devil's staircase". In any case nonperiodic orbits form invariant curves with non rational rotation numbers and that set has a large measure. In fact the resonant effects around any rational $n/m$ appear when the difference of the ratio $\omega_1/\omega_2$ from $n/m < 1$ is smaller than $\varepsilon/m^3$ i. e.,

Therefore the resonant effects affect an interval $2\varepsilon/m^3$. The sum of all the resonant intervals for $n < m$ is smaller than

where $c = 1.64$ ^[38]. Thus if $\varepsilon$ is small the total set of near resonant values of $\omega_1/\omega_2$ is small.

But when the perturbation increases (Figure 7a) to $\varepsilon = 4.4$ most invariant curves are destroyed. However the periodic orbits (mostly unstable) remain. Then for even higher perturbations (Figure 7b) $\varepsilon = 4.5$ a new type of periodic orbits appear, which are called "irregular" (vertical lines) ^[28]. These have in general rotation numbers different from those of their surroundings, or do not have any definite rotation number at all. In fact while most periodic orbits are produced as bifurcations of the central periodic orbit (1:1), or of the boundary, at smaller energies, the irregular orbits are produced in couples (one stable and one unstable) at a minimum $\varepsilon$ (this is called a "tangent bifurcation"). The number of irregular orbits increases exponentially with their multiplicity $m$, while the number of regular orbits increases as $O (m^2)$ ^[50]. Thus the irregular orbits dominate the system for large perturbations.

Many irregular orbits are found in the lobes of the asymptotic curves of any unstable periodic orbit $O$. The lobes are formed by the intersection of the stable and unstable manifolds (asymptotic curves) from the periodic orbit (Figure 8). Such orbits cannot be regular, i.e., they are not bifurcations of regular orbits that exist outside the loops. In fact, if there would be a continuous family of orbits inside and outside a loop there should be a case where the orbit would cross the boundary of the loop and this orbit would be asymptotic and not periodic. This was verified later by Simó and Vieiro ^[115].

3.   Resonance overlap, chaos

The third integral is applicable if the perturbation of an integrable system is small. If the perturbation is large most orbits are chaotic and the third integral is destroyed.

Chaos appears in nonintegrable systems near the unstable periodic orbits. E.g., on a surface of section around an equilibrium point we may have a set of 3 islands around 3 points representing a stable periodic orbit and 3 chaotic regions around 3 points representing an unstable periodic orbit (Figure 9b).

In an integrable case there is no chaos and the 3 unstable points are joined by two separatrices (asymptotic curves forming the limits of the islands around the 3 stable points. (Figure 9a)). However if we calculate the asymptotic curves from the 3 unstable points in a nonintegrable case, they are intersected at an infinity of homoclinic points (Figure 9b). These intersections generate chaos, i.e., most nearby orbits form a set of scattered points on a surface of section (Figure 10) close to the unstable periodic orbits and their asymptotic curves.

In a nonintegrable system there is an infinity of stable and unstable periodic orbits at various distances from the equilibrium point. For small perturbations the islands of various orders and the chaotic regions around every unstable orbit are well separated by invariant (KAM) curves closing around the origin (Figure 10a). However when the perturbation increases the various resonances overlap (Figure 10b) and a large degree of chaos is generated. Then the KAM curves around the origin are destroyed and the third integral is no more applicable. \clearpage

A way to find the critical perturbation when large chaos is introduced and the third integral breaks down is by calculating the sizes of the islands using a different resonant form of the third integral for every set of islands. I did such a calculation in 1966 ^[26]. Around the same time a similar calculation was done by Rosenbluth et al. ^[113]. If the theoretical sizes of the islands become larger than the distance between the various types of islands it is obvious that the theoretical results are not valid any more and an overlap of resonances appears.

The theory of resonance overlap was later developed in detail by Chirikov ^[18] and by many authors it is called Chirikov's criterion. However Chirikov refers to my early paper ^[26] and to the paper of Rosenbluth et al. ^[113] when dealing with this problem ^[128]. I had sent my 1966 paper to Chirikov and he wrote me on 23-10-68: "My criterion of stochasticity by resonance overlapping is essentially the same as yours in your paper in Bulletin Astronomique".

A way to find whether an orbit is ordered or chaotic is by calculating its maximal Lyapunov characteristic number (LCN). Namely if $\boldsymbol{\xi}, \boldsymbol{\xi_0}$ are the infinitesimal deviations from this orbit at times $t$ and $t_0 = 0$ we have

where $\chi$ is the "finite time LCN". If $LCN = 0$ the orbit is ordered while if $LCN > 0$ the orbit is chaotic.

In numerical calculations we find that $\log \chi$ has irregular variations for small t but as $\log t$ increases it tends to vary linearly downwards for ordered orbits (orbit (1) in Figure 11, or it tends to a constant value (orbit (2) in Figure 11).

Theoretically $\boldsymbol{\xi}$ was found by Contopoulos, Galgani and Giorgili ^[47] as a function of $\boldsymbol{\xi_0}$ and $t$ by solving the variational equations together with the equations of motion.

As the calculation of the limiting values of $\chi$ requires a long computer time, there have been various more efficient methods for finding ordered and chaotic domains in dynamical systems. In particular in mappings one can use the "stretching number"

i.e., a local LCN, after only one iteration. Then after about 20 successive iterations we can find whether an orbit is ordered or chaotic ^[125]. Similar methods were developed by several other authors. A comparison of the various methods can be found in my book ^[38].

4.   Galactic dynamics

During my stay at the Yerkes Observatory in 1963 Dr. L. Woltjer visited us and we had many discussions about the third integral. Woltjer argued that my integral was applicable only to smooth potentials but could not be applied to more complicated cases, like the spiral arms of a galaxy. Dr. Chandrasekhar was present during our discussions and he sided with Woltjer. But the evening, during the dinner, I could convince Woltjer that I was right. The next morning Dr. Chandrasekhar told me "George, I thought about your problem and I think that Woltjer is right". But I replied "Chandra you are supposed to be an arbiter between me and Woltjer. But now Woltjer agrees with me. Thus you are obliged to agree also". Then Woltjer and I wrote a paper on a rough model of a spiral galaxy, that was published later ^[46]. This was a parallel approach to the density wave theory of spiral galaxies that was developed at that time by Lin and Shu ^[97].

The density wave theory of the spiral galaxies was first developed by Lindblad ^[98]. His work was most remarkable, since it predated the theory of plasma density waves. Lindblad pointed out, correctly, that this theory was applicable both to trailing and to leading spirals ^[99]. However Lindblad was prejudiced in favor of leading spirals ^[99], something that was objected by the majority of astronomers. This fact and the mathematical difficulties of his papers made most people to disregard his work.

When I met C. C. Lin in Boston, during my visit there in 1963, he told me that he wanted to develop a density wave theory of spiral structure. I told him that Lindblad had already done such a theory and we found in the MIT library several related volumes of the Stockholm Observatoriums Annaler. However the next day C. C. Lin told me that he found Lindblad's papers too difficult to follow, thus he would try to develop a theory from scratch (^[97], and later papers).

Later I was involved quite a lot in the density wave theory. While Lin and Shu, and others, developed the linear theory, I developed a nonlinear theory ^[29,31] that is particularly useful near the resonances. In a galaxy we have two basic frequencies, the epicyclic frequency $\omega_1$ and the angular frequency $\omega_2$ in a frame rotating with the spiral arms. We have a resonance when $\omega_1 / \omega_2 = m / n$ (integers). The linear theory is a first order theory away from resonances, but it is not applicable near the main resonances of a galaxy, the inner and outer Lindblad resonances, where $\omega_1 / \omega_2 = \pm 2/1$, and corotation (particle resonance), where $\omega_2 = 0$ ^[30].

I could show that the nonlinear theory is an application of the third integral theory. In particular the third integral gives the main periodic orbits in the galaxy. Away from the resonances the spiral arms can be constructed by the successive periodic orbits, which are like precessing ellipses (Figure 12). The density maxima are along two symmetric spirals, where the successive orbits approach each other.

However near the inner Lindblad resonance there are two populations of near elliptical periodic orbits (Figure 13) perpendicular to each other. Outside this resonance the dominant family of periodic orbits forms the usual maxima along the spiral arms, while inside the resonance the response density is completely out of phase with the imposed density and tends to destroy it. Thus finally the density wave terminates near the inner Lindblad resonance.

On the other hand near corotation we have two couples of families of periodic orbits, close to the Lagrangian points $L_1$, $L_2$ (unstable) and $L_4$, $L_5$ (stable unless there is a very strong bar at the center of the galaxy ($L_3$), in which case $L_4$, $L_5$ become unstable). Around the stable $L_4$, $L_5$ orbits there are vortices formed by banana-like orbits (Figure 14) ^[30].

A secondary but important resonance is the $4/1$ resonance. Near this resonance the periodic orbits have the form of squares with smooth corners (Figure 12) ^[35]. However, while inside the resonance the orbits support the spiral arms (thick curves), outside the resonance the orbits have a quite different orientation and do not support the spirals. As a consequence the relatively strong spirals (of order $10\%$) must terminate near the $4/1$ resonance, unless the amplitude of the spirals is quite small ($< 5\%$), in which case weak extensions of the spirals can extend up to corotation.

On the other hand in barred galaxies the perturbation is of order 50-100%. In such cases the orbits support the bar up to corotation, but tend to destroy it beyond corotation. Thus the bars terminate close but inside corotation ^[32]. (However there are cases where a bar terminates well inside corotation ^[120]). Beyond corotation we found that most orbits in barred galaxies are chaotic ^[91] (Figure 15).

5.   3D problems

The extension of our studies to 3 degrees of freedom opened new important areas of research. First of all we found new stability types of periodic orbits. In fact the eigenvalues (characteristic exponents) $\lambda$ of a 3-d periodic orbit are given by an equation of the form

The Eq. (5.1) can be factorized as follows:

where

with

Thus the eigenvalues are

(hence $\lambda_1\lambda_2 = \lambda_3\lambda_4 = 1$, i.e., the eigenvalues are inverse in pairs). If $|b_1| < 2$ and $|b_2| < 2$ the orbit is stable. Then the eigenvalues are complex conjugate on the unit circle. If $|b_1| < 2, |b_2| > 2$ the orbit is simply unstable. Then the eigenvalues $\lambda_1, \lambda_2$ are on the unit circle and the eigenvalues $\lambda_3, \lambda_4$ are on the real axis. Similarly if $|b_1| > 2, |b_2| < 2$. If $\Delta > 0$ and $|b_1| > 2, |b_2| > 2$ we have double instability and all eigenvalues are on the real axis. Finally if $\Delta < 0$ we have complex instability. Then the eigenvalues are complex of the form

This classification is a simplification of a classification devised by Broucke ^[12].

E.g., in a Hamiltonian

a particular type of periodic orbits (family 1a) in a diagram ($\epsilon$, $\eta$) has domains of stability (S), simple instability (U), double instability (UU) and complex instability ($\Delta$)b (Figure 16)^[42]. The most important type of instability is the complex instability ($\Delta$ in Figure 16). (Details about the complex instability can be found in my book ^[38]). If we change the parameters A, B, C of Eq. (5.7) we can predict theoretically the behaviour of the various periodic orbits ^[34].

The periodic orbits can be found theoretically if we use the formal integrals of the "third integral" type.

If, besides the Hamiltonian (5.7) there is only one (formal) integral

then a periodic orbit has either 4 characteristic exponents zero or

(the subscripts denote partial derivatives) This is similar to Eq. (2.3) in 2 dimensions. Then a family of periodic orbits depends on one parameter (the energy h).

If there are two (formal) integrals $\Phi^{(1)}$ and $\Phi^{(2)}$, then either all 6 characteristic exponents are zero or all the $3\times 3$ determinants of the matrix

are zero. Then a family of periodic orbits depends on 2 parameters.

If we solve Eq. (5.7) for $\dot{z}$ and insert this value in Eq. (5.8) for $z = 0$ we find

This represents a 3-D surface in the 4-D space $(x, y, \dot{x}, \dot{y})$. Then we have the following cases ^[33].

● Case 1. A periodic orbit is represented by a singular point. In this case the Eq. (5.9) are satisfied.

● Case 2. There is a family of periodic orbits for a given value of $h$. Then Eqs. (5.9) are not satisfied and there are 4 characteristic exponents zero. But in this case the $3\times 3$ determinants (5.10) are zero.

● Case 3. This is an exceptional case where all 6 characteristic exponents are zero.

Numerical experiments for the periodic orbits of the Hamiltonian (5.7) and other polynomial Hamiltonians have given good agreement with the theoretical calculations of periodic orbits using integrals $\Phi^{(1)}$ and $\Phi^{(2)}$ truncated at order 6.

I spoke about these new phenomena at a meeting in California (1985). After my talk I was approached by Prof. Lichtenberg, who told me "Usually the speakers at a meeting speak about a subject that they have covered, again and again, in previous meetings. Your talk was the only novel lecture today". After that we had long discussions with Profs. Lichtenberg and Lieberman at their Department. Later Lichtenberg and Lieberman ^[96] wrote a very interesting book on "Regular and Chaotic Dynamics", where they refer extensively to my work.

A numerical study of nonperiodic orbits in the system (5.7) ^[47] has shown a large ordered domain, an important chaotic sea, and smaller chaotic seas. In the ordered domain there are three formal integrals of motion, in the small seas there are two integrals and in the large sea there is only one integral (the energy).

Using the computer program of Giorgilli ^[79] we calculated 3 formal integrals

where $\Phi_{in}$ contains terms of degree $n$. In particular $\Phi_{i0}$ are the harmonic energies along the axes x, y and z. We found that there is an ordered region where all these formal integrals (truncated at order 8) were approximately conserved, a large chaotic region (big sea) where all three formal integrals had large variations and only the total energy $H = h$ was conserved, and some small regions (small seas) where we had one formal integral approximately conserved, besides the total energy.

A check of the validity of the integrals along particular orbits was made by calculating the finite time Lyapunov characteristic number $\chi$ for a large interval of time. In the ordered region the value of $\chi$ decreases linearly in time (Figure 17) thus the limiting value of $\chi$ is zero. In the chaotic regions the value of the "finite time LCN" $\chi$ tends to a constant $LCN\neq 0$. But in the big sea the limiting value of LCN is large (and independent of the initial conditions) while in every small sea the limiting LCN is maller than in the big sea. However the question is whether the small seas communicate with the big sea after a long time. This is probable because of Arnold diffusion (see below). If this happens then the final value of LCN for all chaotic orbits should be the same.

Arnold diffusion refers to a diffusion that appears in systems of three or more degrees of freedom despite the KAM theorem. In fact in a system of two degrees of freedom the KAM tori are 2-dimensional and separate the orbits in the 3-D phase space inside and outside these tori. However in a system of $N\geq 3$ degrees of freedom the KAM tori are N-dimensional in the 2N-1 phase space and do not separate an inner and an outer domain (E.g., if $N = 3$ the KAM tori are 3-dimensional in a 5 dimensional phase space). Thus there can be diffusion from one side of an invariant (KAM) torus to the other. This is called Arnold diffusion ^[4] and appears even if the perturbation from an integrable case is infinitesimal. Of course in an integrable case all orbits lie on invariant surfaces and no diffusion appears.

In order to understand the Arnold diffusion we may reduce all dimensions by 2. Thus the phase space is 3-dimensional and the KAM tori are 1-dimensional, like strings. These strings leave gaps around them and they allow the motions to go all over the 3-D phase space.

In the gaps between invariant tori there are resonant unstable periodic orbits surrounded by layers of chaotic orbits. These layers communicate with each other and generate Arnold diffusion.

The periodic orbits correspond to resonances of the form

where $v_1 = \omega_1/\omega_3, v_2 = \omega_2/\omega_3$ and $m_1, m_2, m_3$ are integers. The lines in the plane $(v_1, v_2)$ form the so called Arnold web (Figure 18). A chaotic orbit starting close to a point A can go close to a distant point B by following approximately successive lines of the web. Close to every intersection of these lines we have a resonances overlap and an orbit can go from one resonance to another.

The estimates for the diffusion rate according to Arnold are of order

where $\varepsilon$ is the perturbation, therefore the diffusion is exponentially slow. However, if $\varepsilon$ increases considerably the approximate lines become very thick and we may have a direct transition from A to B.

The phenomenon of Arnold diffusion was studied in detail by Froeschlé et al. ^[75] and other related works (see, e.g., ^[95]). A more extensive study was given recently by Efthymiopoulos and Harsoula ^[66] that shows clearly the details of Arnold diffusion.

6.   Relativity and cosmology

Dr. Chandrasekhar and I had some rather interesting interactions on Relativity. Our first work was on the post-Newtonian approximations ^[14]. We first found the transformations of coordinates and time that generalize the Lorentz formulae in the first post-Newtonian approximation.

Later Chandrasekhar asked me to find order and chaos in the case of two fixed black holes. In fact Chandrasekhar ^[13] had calculated a few orbits in this case and he wanted me to distinguish between ordered and chaotic orbits. I wrote three papers on this topic ^[36,37,45] and this was the first systematic study of order and chaos in Relativity.

First of all I calculated the main periodic orbits around the two black holes (Figure 19) (orbits around the mass $M_1$ (a), or $M_2$ (b), and orbits surrounding both $M_1$ and $M_2$ (c), or forming a figure eight around $M_1$ and $M_2$ (d)). Then I distinguished orbits of type Ⅰ falling on $M_1$, of type Ⅱ falling on $M_2$ and of type Ⅲ escaping to infinity. There are infinite sets f orbits of type Ⅰ, Ⅱ and Ⅲ. In fact between 2 orbits of two different types (say Ⅰ and Ⅱ) there is an orbit of the third type (in this case Ⅲ) (Figure 20). (Between an orbit of type Ⅱ and an orbit of type Ⅲ there is an orbit of type Ⅰ and so on). Thus the orbits form fractal sets, and are in general chaotic. However there are some intervals of regular orbits. Namely there are some stable periodic orbits, surrounded by quasi-periodic orbits, as seen in Figure 21.

A by-product of this study was that in the classical problem of 2 fixed centers there are no periodic orbits of type (a) and (b) around a single mass $M_1$ or $M_2$. This fact had escaped the notice of previous authors on the two fixed centers, but I gave an analytic proof of the nonexistence of such orbits ^[36].

The case of relativistic order and chaos received much attention in the following years. In particular the question rose whether the mixmaster cosmological model is integrable or not. The mixmaster model ^[6,101] has expansion and contraction along different directions, which change in a chaotic way an infinite number of times as we approach the initial singularity (backwards in time), This problem can be written by the following Hamiltonian

with energy equal to zero. Then the question is whether this Hamiltonian has the Painlevé property ^[112], i.e., whether the solutions can be expressed in Laurent series of time around every singular point.

We found a general solution, depending on six arbitrary constants, which is of Painlevé type ^[48] and this provided an indication of integrability. However later it was realized ^[94], ^[49] that there is another general solution, which is not of Painlevé type. Thus the mixmaster model is not integrable.

But is the mixmaster model chaotic? In a paper by Cushman and Sniatycki ^[55] it is stated that the mixmaster model is not chaotic, while a paper by Cornish and Levin ^[54] states that the mixmaster model is chaotic.

I pointed out that both papers are correct from different points of view. The first paper proves that almost all mixmaster orbits escape to infinity, thus they cannot be ergodic. On the other hand the second paper emphasizes the fractal structure of the orbits. Thus although the mixmaster model is not chaotic in the usual sense (because its phase space is not finite) it is a case of "chaotic scattering".

I spoke about this problem in Moscow in 2000. Then a Russian participant (Dr. Cernin) presented a paper under the title "The Contopoulos' paradox". The paradox was my opinion that the above papers are both correct.

7.   Quantum chaos

An important development in recent years was the study of order and chaos in quantum mechanics. In the Bohmian approach of quantum mechanics ^{[57,58,59,60,61]}, ^[8,9] we can define orbits by using, besides the Schrödinger equation, the pilot wave equation.

where $\Psi$ is a solution of Schrödinger's equation and $\Im$ indicates the imaginary part.

The Bohmian orbits are either ordered or chaotic (Figure 22). We found cases where the orbits are ordered classically and chaotic in the corresponding quantum problem, and vice versa, cases where the orbits are chaotic classically while the corresponding quantum orbits are ordered ^[65].

Then we checked whether the Born rule of quantum mechanics, that gives the probability as the square of the wavefunction

can be derived theoretically. In fact Bohm and Vigier ^[10] claimed that an arbitrary initial particle distribution $P$ tends asymptotically to $|\Psi|^2$.

However we found ^[65] that this is not the case for ordered orbits. On the other hand for chaotic orbits in most cases the value of $P$ approaches $|\Psi|^2$, following an exponential law in time. This defines a quantum analog of the Nekhoroshev time of classical mechanics. However there are also exceptional cases where even some chaotic orbits do not follow the Born rule ^[52].

The solutions of Eq. (7.1) have a singular point when $\Psi = \Psi_R+i\Psi_I = 0$. This is called a nodal point. The orbits around the nodal point are of spiral form (Figure 23). The nodal point is either an attractor or a repellor. But close to the nodal point there is a second singular point (the X-point in Figure 23), which is unstable. If an orbit never approaches the region around the moving nodal point (the nodal point-X-point complex) then it is ordered and not chaotic (red orbit in Figure 22). However if an orbit approaches the region around the nodal point, then it is chaotic (black orbit in Figure 22).

The regular orbits can be given by series expansions involving trigonometric terms ^[69,40]. These series are similar to the classical formulae derived by the third integral theory for periodic time potentials. Chaos is generated when the orbits approach the X-point near the nodal point. This is an important new result in the dynamics of Bohmian quantum mechanics.

We considered also cases with more than one nodal point at every time ^[70] and cases of 3 dimensions ^[122,123]. Furthermore we found the diffraction patterns due to small material targets. In this case we have large numbers of nodal points along a line called "separator" (Figure 24)^[62]. The particles follow nearly parallel orbits up to the separator, but their orbits deviate almost radially from the target beyond the separator. Orbits passing close to the target have small deviations, while orbits further away have larger deviations (Figure 24a). This behaviour is quite different from the classical picture (Figure 24b). where orbits coming close to the target have larger deviations than orbits passing further away.

The geometry of the quantum orbits explains the diffraction patterns observed in scattering experiments, including the appearance of enhanced flow along certain angles called Bragg angles.

An important consequence of this study is that the times of flight of orbits passing close to the target are larger than the times of flight of further away orbits (Figure 24a)^[71]. This is the opposite of what happens in the classical case (Figure 24b) and requires an experimental verification that would support the Bohmian approach of quantum mechanics.

In 3 dimensions we have figures like Figure 23 at every level z. Thus we have nodal lines and X-lines (Figure 25). The orbits are chaotic if they approach these lines. In particular cases we have one or two integrals of motion. In the case of one integral (called partially integrable) the motions take place on a surface ^[53,121,122] (Figure 26) and can be either ordered (blue curve) or chaotic (red curve) on this surface. In the case of two integrals (integrable case) all orbits are ordered. However in general orbits have no integral, but are chaotic filling a volume in 3 dimensions (Figure 27).

8.   Integrals in chaos

An important recent development refers to the use of the third integral in calculating chaotic orbits. As it was shown by Moser ^[103,104] the third integral around an unstable periodic orbit is convergent. This is due to the fact that there are no small divisors in the series of the third integral in this case, while there are infinitely many small divisors in the series around a stable periodic orbit. In fact near a stable periodic orbit the (formal) integral has divisors of the form $m_1\omega_1+m_2\omega_2$ (where $\omega_1$ and $\omega_2$ are real and $m_1, m_2$ are integers) and some of these divisors are very small, generating nonintegrability. However near an unstable periodic orbit one eigenvalue is real and the other is imaginary, therefore the divisors $m_1\omega_1+m_2\omega_2$ never become small. This problem was first discussed by Cherry ^[17] and in more detail by Moser ^[103,104]. The proof of Moser was completed by Giorgilli ^[79]. The method of Moser was applied by da Silva Ritter et al. ^[56], Vieira and Ozorio de Almeida ^[124] and Ozorio de Almeida and Vieira ^[110] in finding the asymptotic curves and nearby orbits in the case of the Hénon hyperbolic map:

Da Silva Ritter, Ozorio de Almeida and Douady ^[56] introduced a new, near to the identity transformation in the variables

such that in the new variables ($\xi, \eta$) the transformation becomes

where

and

($\lambda_1$ being the largest eigenvalue of the mapping for a given $\kappa$ and $w_i$ are constants).

We extended this work by calculating the limits of applicability of the Moser series ^[72]. Namely the mapping (8.1) gives the consequents of an initial point $(\xi_0, \eta_0)$ along a hyperbola $\xi_0 \eta_0 = c$ (Figure 28a). The limiting value of $c$ is found by using the d' Alembert criterion. Namely we calculated the limits of the ratios

along any given direction $\phi$, where $\xi = R \cos \phi, \, \eta = R \sin \phi$.

We found numerically that in this particular case the limits of the ratio (43) depend on $R$ but not on $\phi$.

For $\kappa = 1.43$ the series converge for $|c| \leq 0.49$ ^[83].

Any orbit with initial conditions inside the domain of convergence, has its images and pre-images along a hyperbola passing through this point (Figure 28a).

If we back-transform the hyperbolae of Figure 28a to the original variables $(x, y)$, by inverting Eqs. (8.2), we find the image of the domain of convergence of the Moser series in the plane $(x, y)$ (red in Figure 28b). This domain is limited outwards by the images of the curves $c = \pm 0.49$ of the regions 2, 3, 4 of Figure 28a. However there is also an inner limit, which is the image of the curve $c = 0.49$ of the region 1, in the ($\xi, \eta$) plane. Thus there is an inner region (white in Figure 28b) where the series do not converge. This region is around a stable periodic point $S$.

In Figure 29 we have drawn a number of invariant curves $\xi \eta = c$ as they are mapped in the variables $(x, y)$. We call these curves "Moser invariant curves". In particular the curves $c = 0$ (red) represent the axes $\eta = 0$ and $\xi = 0$. These curves correspond to the stable and unstable asymptotic curves of the point $(0, 0)$ and extend to infinity within the domain of convergence, because the series (39) converge all the way to infinity if $\eta = 0$ or $\xi = 0$. The images of the curves $\xi = 0$ and $\eta = 0$ intersect each other at an unlimited number of homoclinic points that can be found analytically.

Any point in the domain of convergence (red) of Figure 28b has its images inside the same domain. On the other hand orbits with initial conditions outside the domain of convergence have images that approach closer and closer the outer limits of this domain. This is seen in Figure 28b where we have mapped a grid of initial conditions $(-10, 10) \times (-10, 10)$ outside the convergence domain. Their first images are in green and their second images are in blue (covering also the inner green regions). Higher order images are congested even closer to the outer limits of the domain of convergence. Thus the outer limits of the domain of convergence acts as an attractor for the orbits outside this domain ^[41].

The main application of the Moser series regards the orbits starting close to the origin $O$ where chaos is dominant. In this case the successive images of an initial condition close to $O$ look quite random (Figure 30a). However, all these points lie on a particular Moser invariant curve (Figure 30b) that can be given accurately analytically. Thus these chaotic orbits can be given by analytical formulae. The only indication of chaos is that although all the iterates lie on the same Moser curve, the distance between successive iterates increases at every iteration. Furthermore, the Moser invariant curves make many oscillations that extend to large distances, therefore some points on them may be at considerable distances from the center.

A particular application of the Moser series is in finding periodic orbits.

In Figure 31 we give the characteristics of the various families 6, 5, 4, 3, 2 together with the characteristics of the periodic orbit $S$ and of the homoclinic point $H$. In particular the two families of period 4 are generated from the stable family $S$ at $\kappa = 1.317$. The upper family is initially stable while the lower one is unstable. The stable family becomes unstable at a period doubling bifurcation, and all the families that were produced by a cascade of further bifurcations become unstable beyond $\kappa = 1.383$. The same happens with all the families of order higher than 4. All these families are congested close to the homoclinic point $H$ for larger values of $\kappa$ (Figure 31).

Of special interest is the family of period 2 because when this family is bifurcated from the original family $S$ (for $\kappa = 1.76$), the orbit $S$ itself becomes unstable. The bifurcating stable family 2 also becomes unstable at $\kappa = 1.84$ and then follows a cascade of periodic doubling bifurcations so that beyond $\kappa = 1.86$ all the bifurcated families become unstable.

When $S$ becomes unstable it has its own asymptotic curves and new Moser invariant curves close to them. We found a Moser transformation that gives these new asymptotic invariant curves ^[41]. There is a Moser domain of convergence of these transformations (green) which is completely inside the Moser domain of convergence around the origin point $O$ (red in Figure 32a). In this case ($\kappa = 2$) there is no inner limit of the domain of convergence around $S$.

The asymptotic curves from the orbit $S$ intersect the asymptotic curves from the orbit $O$ at heteroclinic points (Figure 32b). All the homoclinic and heteroclinic points are found analytically with a very good accuracy as compared with the numerical results ^[41].

An application of the Moser series is in finding the chaotic orbits that generate the spiral arms outside corotation in barred galaxies ^[84]. In fact the spiral arms in barred galaxies are composed of chaotic orbits that can be given theoretically by using the Moser theory. Thus the use of the Moser series has important applications in galactic dynamics.

Further details about the history of the third integral can be found in the books of Contopoulos ^[38,39].

Conflict of interest

The author declares no conflict of interest.