text:
# Stable regions and singular trajectories in chaotic soft wall billiards

## 1 Introduction

Kolmogorov-Arnold-Moser (KAM) theory provides a framework for the understanding of an integrable Hamiltonian system which is subject to a perturbation . As the strength of the perturbation is increased areas of phase-space became chaotic and eventually only islands of stability, surrounded by a chaotic ”sea”, survive. The opposite question is of interest as well: What will happen to a completely chaotic system when a perturbation is applied? will it remain chaotic or will islands of stability appear? The appearance of KAM islands is of great importance: First, for initial conditions falling inside the island, there will be finite range oscillations in the the tail of the correlation function. In addition, these islands are usually accompanied by areas of ”stickiness” where particles belonging to the chaotic region of phase space spend an anomalous amount of time, and hence they affect the temporal correlation function even for initial conditions in the chaotic ”sea”. Indeed, it was argued that this stickiness can cause a power-law decay of correlations, instead of an exponential one, and thus affect the transport properties of the whole system even for a relatively small island size.

For the billiard system (i.e. a particle moving freely in a bounded region and reflecting elastically from the boundary), one such perturbation is making the walls ”soft”, as opposed to the infinitely steep potential of the ideal case. Since the billiard is a widely used paradigm for understanding the microscopic foundations of statistical mechanics, and noting that real classical molecules move in smooth potentials, we conclude that the study of soft-wall billiards is an important step towards sustaining these foundations and, in particular, Boltzmann’s ergodic hypothesis.

For a certain kind of dispersing billiards it was theoretically proven that when the wall becomes soft, an island appears near a singular, tangent trajectory (see details in section 3). More recently, and partly inspired by the numerical results presented herein, the appearance of an additional island, this time near a corner, was also shown. Experimentally, some conjectures were made about the role played by the softness of walls in semi-conductor quantum dots. It was considered to be responsible for the algebraic tail in the distribution of the length of trajectories, which gives rise to the fractal nature of conductance fluctuations. Observations of the formation of stable islands in a stadium shaped soft atom-optics billiard were reported in Ref. .

In this article, we present additional results of numerical and experimental studies of soft wall atom-optics billiards. Three unifying concepts are found in our results: First, islands of stability arise around singular trajectories (tangencies, corners) as predicted by Refs. . Next, the size of these islands depend on the wall’s softness in a very sensitive way. Finally, ”stickiness” regions surround the stable islands and affect the decay properties of the billiards. In section 2, we present our numerical and experimental studies of the stability induced in an atom-optics billiard of the Bunimovich type by making the walls soft. In section 3, we present results of numerical simulations for the dynamics of atoms inside totally dispersing billiards, which exhibit softness induced stability around singular trajectories.

## 2 Islands of stability in converging billiards

The atom-optics billiard is a recently developed experimental system in which billiard dynamics can be studied. As described in detail in Ref. , the billiard is realized by the use of a laser beam which is rapidly ($`100`$ kHz) scanned using two perpendicular acousto-optic scanners (AOSs). The laser frequency is tuned above the atomic resonance (the $`D_2`$ line of 85Rb in our experiment), hence it applies a repulsive force on the atoms. By controlling the deflection angle of both AOSs, we create the required billiard shapes which confine the atoms in the transverse direction. The instantaneous potential is given by the dipole potential of the laser gaussian beam:

$$U(x,y,t)=U_0\mathrm{exp}\left[2\left(\left(xx_0\left(t\right)\right)^2+\left(yy_0\left(t\right)\right)^2\right)/w_0^2\right],$$
(1)
where the curve $`(x_0\left(t\right),y_0\left(t\right))`$ is the shape of the ideal billiard, along which the center of the Gaussian laser beam scans. $`w_0`$ is the laser beam waist and $`U_0`$ is the instantaneous potential height. The time averaged potential is kept constant for different values of $`w_0`$ by changing the laser power. In such way, $`w_0`$ serves as the softness control parameter, and is experimentally controlled by the use of a a telescope with a variable magnification, located prior to the AOS’s such that $`w_0`$ can be changed without affecting the billiard’s size and shape. Fast enough scanning of the beam results in an effective time-averaged potential wall . Photon scattering from the laser beam is greatly suppressed due to the large detuning, and collisions between the atoms are very rare for the range of densities used, hence motion of the atoms between reflections from the light walls can be regarded as strictly ballistic.

The loading scheme of cold atoms into the billiard was described in . Briefly, laser cooled 85Rb atoms are loaded from a magneto optical trap into a red-detuned one-dimensional optical lattice which is laterally shifted (250 $`\mu `$m) from the billiard’s location. The atoms are then transferred from the lattice into the billiard by pushing them with a pulse of a strong on-resonance beam which is perpendicular to both the lattice and the billiard beams. Simultaneously with the pushing, the lattice beams are adiabatically switched-off, resulting in additional cooling. Further reduction of the radial velocity spread and especially a decrease in the number of very slow atoms is achieved by capturing only a central velocity group of atoms into the billiard, which is turned on at the proper time after the push. Typically, $`310^5`$ cold 85Rb atoms are loaded into the billiard, with a typical velocity of $`100`$ mm/s and a velocity spread of $`20`$ mm/s. The cooling in the longitudinal direction ensures that the system can be approximated as a two dimensional one.

As a way to probe the nature of the dynamics of the atoms confined by such an ”atom-optics” billiard we measure the decay of the number of confined atoms through a small hole on the boundary. A purely exponential decay is expected for chaotic motion, with a decay time given by $`\tau _c=\pi A/vL`$, where $`A`$ is the billiard’s area, $`v`$ the atomic velocity and $`L`$ the length of the hole . A slower algebraic decay usually characterizes a stable or partly stable phase space. In our system, the hole is produced by switching off one of the AOSs for $`1`$ $`\mu `$s every scan cycle, synchronously with the scan. The number of atoms remaining in the trap is measured using fluorescence detection. The ratio of the number of trapped atoms with and without the hole, as a function of time, is the main data of our experiments.

First, we study a tilted Bunimovich stadium (see Fig. 3), which is chaotic and ergodic for the ideal hard-wall case. We used a ”tilted” stadium, to avoid the long segments of regular motion associated with ”bouncing ball” trajectories in the original stadium. It is composed of two semicircles of different radii ($`64`$ $`\mu `$m and $`31`$ $`\mu `$m), connected by two non-parallel straight lines ($`192`$ $`\mu `$m long). The results of numerical simulations for the decay in the number of confined atoms, through a small hole in the boundary, are presented in Fig. 1. For the stadium with $`w_0=12`$ $`\mu `$m a nearly exponential decay is seen, as expected for the ideal hard-wall billiard. As the walls are made softer, a monotonic slowing down in the decay is observed. We compare this behavior with that of a circle with the same area and hole size, which exhibits integrable motion. In the ideal hard-wall case the existence of nearly periodic trajectories results in many time scales for the decay through a small hole on the boundary, and then yield an algebraic decay. Simulations for the decay from a circular billiard (dotted lines in Fig. 1) show that the decay is almost unaffected by the change in softness parameter, as expected.

Figure 2 shows the Poincaré surface of section obtained from the results of numerical simulations for classical trajectories of Rb atoms inside the tilted-stadium billiard, with $`w_0=24\mu `$m. For clarity, we assume a monoenergetic ensemble (with $`v=120`$ mm/s), a two-dimensional system, and no gravity. However, including the experimental velocity spread, three dimensional motion and gravity have only a small effect on the results. In analogy to our experimental setup, we use for our simulation a scanning gaussian potential, with a scan frequency of 100 kHz. The dimensions of the billiard are equal to the experimental ones. The Poincaré surface of section shows $`v_x`$ versus $`x`$ at every trajectory intersection with the billiard’s symmetry axis ($`y=0`$), provided that $`v_y>0`$. The smoothness of the wall results in a stability region that appears around the singular trajectory which connects the points where the big semicircle joins the straight lines. This fact is also confirmed by analyzing some typical trajectories (see Fig. 3). Additional islands (smaller than the central one) appear around it, and correspond to trajectories with higher periodicity, like the one shown in Fig. 3C. Around these islands there is a large area of ”stickiness”, where the trajectories spend a long time. Such a ”sticky” trajectory is presented in Fig. 3D. The exact structure of the island and its vicinity depends on the softness parameter $`w_0`$ in a sensitive way. In general, the size of the island and the stickiness region increase with the softness parameter. The slowing down in the decay is attributed mostly to the sticky trajectories, which in general cover a phase space volume greater than that of the island itself (see Fig. 2 of ).

In Fig. 4, experimental results for the decay from a tilted stadium shaped atom-optics billiard are presented, for two different values of the softness parameter ($`w_0=14.5`$ $`\mu `$m and $`w_0=24`$ $`\mu `$m). The experimental results confirm that the soft wall causes an increased stability, and a slowing down in the decay curve, in good agreement with the numerical results. We also observed that when the hole includes the singular point where the semicircle meets the straight line, no effect for the change in $`w_0`$ is seen. In this case the stable island that is formed around the singular trajectory and the region of ”stickiness” around it (see Fig. 3) are destroyed by the hole. By moving the hole we can thus experimentally map the location of the stable island.

Similar decay measurements for a circular atom-optics billiard (see Fig. 5) showed no dependence on $`w_0`$ in the range $`14.524`$ $`\mu `$m. We were not able to see any dependence on the hole position either.

## 3 Islands of stability in dispersive billiards

We studied the effect of a soft potential wall also in a dispersing billiard. The billiard that we studied was formed from the intersection of four arcs with different radii and it had a tangency, namely, a specific periodic orbit is tangent to the billiard’s boundary at one point (see Fig. 7). In this billiard, the periodic orbit had a period of two. More generally, we consider a family of such billiards, with a variable parameter $`\zeta `$, which is the vertical distance between the lower arc and the two-periodic orbit. For $`\zeta =0`$ the periodic orbit is tangent to the lower arc. For $`\zeta >0`$ the lower arc is below the periodic orbit, while for $`\zeta <0`$ it is above and actually intercepts the orbit. A billiard with a two-periodic tangent orbit was chosen, to increase the size of the island. Every hard-wall billiard of this family is fully chaotic. This specific billiard shape was motivated by a theoretical study which predicts the formation of an elliptic island around such a tangency .

In Fig. 6, the Poincaré surface of section in the vicinity of the two-periodic orbit is presented, for a billiard with $`\zeta =40`$ $`\mu `$m, and $`w_0=24`$ $`\mu `$m. An elliptic island is clearly observed. Typical trajectories for this case are shown in Fig. 7. In Fig. 7A, a typical chaotic trajectory is shown, corresponding to the areas in phase space marked as ”A” in Fig. 6. A periodic trajectory of period 2 (corresponding to the central island in Fig 6, marked as ”B”) is shown in Fig. 7B. Also shown in Fig. 7 are a high-period periodic trajectory and a ”sticky” one. We repeated these calculations for billiards with different values of $`\zeta `$ and $`w_0`$, and plotted the size of the island as a function of these two parameters in Fig. 8. It can be seen that the island becomes observable above some softness of the potential, and that its size initially increases with an increased softness. However, further softening reduces the size of the island. In general, higher values of $`\zeta `$ require higher values of $`w_0`$ to induce the stability. These results, although obtained for a strongly perturbed case, are in agreement with Refs. . Some intuition can be gained by noting that the stable trajectory is slightly reflected by the lower arc. Since the stable trajectory is almost tangent to the lower arc, its perpendicular velocity is very low, and a small potential is already sufficient to reflect it. This small potential is provided at a certain distance from the lower arc, depending on $`w_0`$. If $`\zeta `$ is increased, the trajectory will ”feel” the lower arc (and hence will be stabilized) only for a larger value of $`w_0`$. However, if $`w_0`$ is further increased, the potential at the trajectory would become too high, thus ”blocking” it and destroying the stable island. In additional simulations we observed that the center of the island moves towards larger $`y`$ values when $`w_0`$ is increased, in good agreement with the above intuitive model.

We observed also another kind of stable islands in these billiards, around a different type of singular trajectories, i.e. trajectories that are supported by the corners of the billiard. An example for such a trajectory, and a sticky trajectory in the area surrounding this new island, are shown in Fig. 9. Our observations are consistent with the theoretical proof for the existence of this large islands, provided in Ref. . There it was shown that, for a small perturbation, its existence depends both on the geometry of the corner and on the smooth potential behavior at the corners. Our results show the existence of this island also for a strong perturbation. In all cases where an island was found, we also found a sticky region near it with atoms spending arbitrary long times near the island.

We performed decay simulations for the scattering billiard, by opening a hole on the boundary, and counting the number of remaining atoms as a function of time. Here, as in the case of the Stadium billiard, an exponential decay over 3 orders of magnitude was observed for $`w_0=8`$ $`\mu `$m. In cases where an island develops, a slowing in the decay was observed when the hole did not include a reflection point of the periodic orbit. This slowing down in the decay was similar to that for the stadium billiard of Fig. 1.

## 4 Summary

In this paper we showed, numerically and experimentally, that stable islands are created in the phase space of chaotic billiards, either dispersive or focusing, when the walls are softened. Using cold atoms in atom-optics billiards, an island of stability was demonstrated in a tilted Bunimovich stadium billiard, around the periodic trajectory that connects the points where the big semicircle joins the straight lines. For a dispersive billiard, we confirmed numerically the existence of a stable island around a periodic orbit which is tangent to the billiard’s wall. In addition, we found numerical evidence, for the first time, for the formation of a stable island near a ”corner” trajectory. In all three cases, the elliptic islands appear around trajectories which are singular in the ”ideal” infinitely hard wall billiard. The orbits are singular in the sense that they do not have a smooth limit when approaching the singularity from both sides.

Our numerical results show that the appearance of a stable island in a soft-wall billiard is in general accompanied by areas of ”stickiness” surrounding it. Although the size of the stable island is a sensitive function of softness, not always showing a monotonic behavior, we confirmed numerically and experimentally that the ”sticky” trajectories modify this behavior and, in general, introduce a monotonic slowing down in the decay of the number of trapped particles, as the walls are made softer.

An interesting question in this context is wether the existence of a singularity in the hard-wall billiard is a necessary and sufficient condition for stability in the soft-wall one. More generally, one can ask to what extent is the relation between the existence of a singularity in a chaotic system and the emergence of stability when this system is perturbed, a universal feature. From the practical point of view, the softness of the wall can be used as a control parameter for the structure of phase space, in particular in the context of controlling chaos. Finally, we note that since stable and ”sticky” regions follow the softness of the walls, and since any physically realizable potential is inherently soft, the physical realization of a truly chaotic and ergodic billiard is questioned.

## 5 Acknoledgements

We acknowledge helpful discussions with V. Rom-Kedar and support from the Israel Science Foundation, the Minerva Foundation, the United States-Israel Binational Science Foundation and Foundation Antorchas. M. F. Andersen acknowledges support from the Nachemsohn Dansk-Israelsk Studienfond.


format: text

text:
# Anomalous Diffusion in Infinite Horizon Billiards

## I Introduction

Diffusion of particles in an infinite domain billiard is a well studied problem \[1-5,9,14\]. By a billiard we refer to the motion of a point particle in a two dimensional domain in which the particle moves with constant velocity in straight line orbits executing specular reflection (i.e., angle of incidence equals angle of reflection) from fixed boundaries. By an infinite domain we refer to an unbounded two dimensional region. An early consideration of diffusion in a billiard as a model in physics was made by Lorentz  to model electrons in a metal. In this model (called the Lorentz gas) particles move freely and reflect specularly from fixed, randomly placed, hard-wall scatterers. A modification of the two dimensional Lorentz gas in which there are circular scatterers on a square lattice is an example of an infinite horizon billiard, called the Sinai billiard , and is illustrated in Fig. 1(a).

Infinite horizon billiards are the subset of infinite domain billiards that contain channels through which a particle may pass without ever reflecting off a billiard wall.

In this paper we consider diffusion in infinite horizon billiards. The examples that we study numerically are shown in Figs. 1(a)-(d). The billiards in Fig. 1 include: (a) the Sinai billiard, composed of circular, hard wall scatterers arranged on a square lattice such that the scatterers do not touch each other; (b) a modification of model (a) in which the circular scatterers are randomly displaced (random in direction and magnitude) by at most $`\mathrm{\Delta }<L/2R`$, so that there are channels of width $`L2(R+\mathrm{\Delta })`$ accommodating free motion; (c) randomly oriented square scatterers on a square lattice; and (d) the scalloped channel, in which the domain is infinite in the $`y`$ direction and bounded in the $`x`$ direction by circular arc segments, each subtending an angle less than or equal to $`180^o`$. Figure 1(d1) shows the case of the scalloped channel where the circular arc segments are semi-circles, while Fig. 1(d2) shows the case where the arcs subtend an angle less than $`180^o`$. Particle motion for the situation in Fig. 1(d1) is equivalent to particle motion for a stadium-type billiard \[see Fig. 2(a)\]; the particle motion within the scalloped channel can be folded into the stadium billiard via reflection of the particle at a straight wall as it passes to the next cell.

By cells we mean each portion of the scalloped channel domain between the dotted lines of Fig. 1(d). In a similar manner particle motion in the bounded billiard of Fig. 2(b) can be thought of as equivalent to motion in the infinite billiard of Fig. 1(a).

One important means of characterizing transport in an infinite domain two dimensional billiard is through the phase space probability distribution function (pdf), $`P(x,y,\theta ,t)`$, where $`\theta `$ is the angle of the particle velocity and we take all particle velocities to have magnitude one. From the pdf one can calculate the displacement moments of the distribution. The $`q`$th moment at time $`t`$ is

$$r^q=(x^2+y^2)^{q/2}=_0^{2\pi }_{}(x^2+y^2)^{q/2}P(x,y,\theta ,t)dxdy\frac{\mathrm{d}\theta }{2\pi },$$
(1)
where $`r=\sqrt{x^2+y^2}`$ and $``$ is the (infinite) spatial domain of the billiard and where throughout this paper we take the initial pdf, $`P(x,y,\theta ,0)`$, to be bounded, $`|P(x,y,\theta ,0)|<K`$, and to be zero outside some finite region. For the infinite horizon billiards shown in Fig. 1 we find that the moments of the displacement have a time dependence

$$r^qt^{\gamma _q},$$
(2)
which we use as shorthand for

$$\gamma _q=\underset{t\mathrm{}}{lim}\frac{\mathrm{log}r^q}{\mathrm{log}t}.$$
(3)
For all cases in Fig. 1 we find the exponent $`\gamma _q`$ to be

$$\gamma _q=\{\begin{array}{cc}q/2\hfill & q<2\hfill \\ q1\hfill & q>2\hfill \end{array}.$$
(4)
Results of the form (2) with $`\gamma _q`$ composed of piecewise linear functions (different from (3)) have also been obtained in other situations of Hamiltonian transport . The occurrence of an exponent $`\gamma _qq/2`$ is commonly referred to as anomalous diffusion, and $`\gamma _q>q/2`$ $`(\gamma _q<q/2)`$ is called superdiffusion (subdiffusion).

Figure 3 shows the results of numerical experiments testing (4). In these numerical experiments we start with a cloud of many initial conditions distributed uniformly in the accessible space occupied by one cell \[the region outside of the scatterer and within $`0yL,0xL`$ for Figs. 1(a)-(c), and $`0yL`$ for Figs. 1(d)\] and uniformly in $`0\theta 2\pi `$. We then evolve the orbit of each particle in the cloud forward in time  and obtain $`r^q`$. In all cases the results we obtain are consistent with (4). The motivation for choosing an initial distribution uniform in space and angle is that motion of an orbit from a typical initial condition in the billiards of Fig. 2 is known to be ergodic, generating an invariant density uniform in space and in angle $`0\theta 2\pi `$. In another set of numerical experiments we used initial conditions distributed uniformly in intervals of $`\theta `$ such that the lengths of the initial flights are bounded \[e.g., for Figs. 1(a)-(c), $`\theta _0\theta \pi /2\theta _0`$, $`\pi /2+\theta _0\theta \pi \theta _0`$, $`\pi +\theta _0\theta 3\pi /2\theta _0`$, $`3\pi /2+\theta _0\theta 2\pi \theta _0`$ where $`\theta _0<\pi /4`$\]. Again, agreement with (3) and (4) was found. This indicates that the angular particle distribution relaxes to the uniform distribution sufficiently fast that the results (3) and (4) are not modified. Explaining the reason for this insensitivity to the initial distribution is one of the main contributions of this paper (Sec. 2C). Our other main contribution is the result that, for long time, initial distributions with no particles in a $`\theta `$ interval about a direction of infinitely long flight (say $`\theta =0`$) lead to long-time distributions with an invariant scaling form. Specifically, $`\widehat{P}(\theta ,t)=_{}P(x,y,\theta ,t)dxdy`$ approaches $`\stackrel{~}{P}(\varphi )`$, where $`\varphi =\theta t`$ and $`\stackrel{~}{P}(\varphi )`$ is independent of $`t`$. Furthermore, the same $`\stackrel{~}{P}(\varphi )`$ universally applies for the billiards of Figs. 1(a), 1(b), 1(c), 1(d2), but $`\stackrel{~}{P}(\varphi )`$ is different for the billiard of Fig. 1(d1). We use this scaling result to show the insensitivity of (4) to the initial particle distribution.

Bleher  showed for infinite horizon billiards that the limit in distribution as $`t\mathrm{}`$ of the pdf of the particle displacements in the billiard is Gaussian with a width that increases with time as $`\sqrt{t\mathrm{log}t}`$. (For ordinary diffusion, the result is the same except that the width increases as $`\sqrt{t}`$.) However, the asymptotic $`t`$ dependence of the moments cannot be calculated from Bleher’s result. This is discussed in Sec. 3. Note that the definition of the symbol $``$ given in Eqs. (2) and (3) is such that logarithmic corrections to the scaling of $`r^q`$ with time are not included \[e.g., if $`r^2Ct\mathrm{log}t`$, where $`C`$ is a constant, as suggested by , then $`\gamma _2`$ from (3) is one, consistent with (4)\].

## II Theory

The difference between particle transport in an infinite horizon billiard and particle transport in the case of normal diffusion is due to the arbitrarily long trajectories found in the infinite horizon billiard. These long flights occur in the channels between the scattering boundaries of the billiard. When a particle is traveling nearly parallel to the axis of the channel \[for Figs. 1(a)-(c) there are many channels parallel to both the $`x`$ and $`y`$ axes whereas for Fig. 1(d) there is only one channel, which is parallel to the $`y`$ axis\], it will travel long distances between reflections off the billiard wall. Let $`\theta _n`$ be the angle between the trajectory of a particle and the axis of the channel after the $`n`$th reflection of that particle with a billiard wall. The length $`r_n`$ of the flight is of the order of $`1/|\theta _n|`$ for $`|\theta _n|1`$.

For some of the billiards we consider there are strong correlations between the $`\theta `$ values from one reflection to the next. This is especially true of the scalloped channel with semicircular arcs (stadium billiard). Here, upon reflection, the angle $`\theta _n`$ can change by, at most, a factor of three ; $`|\theta _n|/3|\theta _{n+1}|3|\theta _n|`$. Thus if $`|\theta _n|`$ is small, $`|\theta _{n+1}|`$ is also small, and both represent long flights. This suggests an extreme model for the particle transport in which the length of a flight for each particle does not change from reflection to reflection; the length of each flight is completely correlated with the previous flight, but different particles have different flight lengths.

### II.1 The Completely Correlated Model

Consider a one dimensional system in which an ensemble of particles executes a random walk. The particles in the ensemble differ from each other in the length of the step, $`\mathrm{\Delta }r_\xi `$, each particle takes.

$$\mathrm{\Delta }r_\xi =1/\xi ,\xi [0,1].$$
(5)
The random variable $`\xi `$ is distributed uniformly in $`[0,1]`$ ($`\xi `$ is inspired by $`\theta _n`$ defined above, but is chosen to occupy the interval $`[0,1]`$ for simplicity), and $`\xi `$ for a particle does not change from step to step (complete correlation). Since every billiard particle has the same speed we normalize the magnitude of the velocity to 1 and thus take the time between steps in our one dimensional random walk model to be

$$\mathrm{\Delta }t_\xi =1/\xi .$$
(6)
Since the random walk behaves like normal diffusion when the number of steps is large, the $`q`$th moment for each $`\xi `$, $`|r|^q_\xi `$, will be

$$|\frac{r}{\mathrm{\Delta }r_\xi }|^q_\xi (\frac{t}{\mathrm{\Delta }t_\xi })^{q/2}\text{, for }t\mathrm{\Delta }t_\xi .$$
(7)
Substitution for $`\mathrm{\Delta }r_\xi `$ and $`\mathrm{\Delta }t_\xi `$ yields, for long time $`t`$, $`(t\mathrm{\Delta }t_\xi )`$

$$|r|^q_\xi (t/\xi )^{q/2}$$
(8)
for each $`\xi `$.

To find the $`q`$th moment for the entire ensemble one needs to take an average over all particles:

$$|r|^q=_0^1|r|^q_\xi d\xi .$$
(9)
This average has two main contributions:

$$|r|^q_{N/t}^1(t/\xi )^{q/2}d\xi +_0^{1/t}t^qd\xi $$
(10)
where $`N`$ is chosen to be large, so that the random walks have made many steps. The first term in Eq. (10) represents the fraction of particles that have executed at least $`N`$ steps and so can be described by Eq. (8). The second term comes from those particles that are still in their first flight (with velocity = 1). The contribution from the interval $`N/t\xi 1/t`$ \[omitted in (10)\] has neither the majority of walkers (in the limit of long time) nor the most extreme displacements and so does not dominate the other terms. For long time Eq. (10) yields

$$|r|^q\{\begin{array}{cc}t^{q/2}\hfill & q2\hfill \\ t^{q1}\hfill & q2\hfill \end{array}.$$
(11)
This simple toy model shows that it is the very long flights allowed by the open channels that give rise to the non-Gaussian behavior for moments greater than 2. One might then think that for the real billiard systems of Fig. 1 the behavior for $`q>2`$ is a trivial result of the initial distribution of angles. However, this is not the case and the behavior (11) is robust to changes in the initial distribution of particles. As discussed in Sec. 2C, non-uniform initial angular distributions scatter rapidly enough that (11) still holds even if the initial distribution has no particles traveling nearly parallel to the channels.

### II.2 Moment Equation for Infinite Horizon Billiards with a Uniform Initial Angular Distribution

We now consider a uniform initial spatial (within a cell) and angular distribution and for the case of the ergodic billiards of Figs. 1(a) (Sinai billiard) and 1(d) (scalloped billiard). In these cases such a distribution is stationary when $`(x,y)`$ is taken modulo the appropriate cell period. We show that $`\gamma _q`$ is given by Eq. (4) for these real billiard systems. The particle transport in an infinite horizon billiard must proceed at least as quickly as a random walk process (i.e., as fast as normal diffusion). While there exist mechanisms for faster than diffusive transport (to be discussed below), there is no stable mechanism to stop or trap a billiard particle. The periodic orbits that might trap the particle in these systems are all exponentially unstable. That a random walk is a lower bound on the transport in the billiard system implies that $`\gamma _qq/2`$.

In addition to diffusive like behavior, particle trajectories that consist of a single long flight (“ballistic” flight) also participate in particle transport. If every particle were to move ballistically, we would find $`\gamma _q=q`$. This provides an upper bound on $`\gamma _q`$; not every particle will execute a single uninterrupted flight. We can, however, place a lower bound on the fraction of particle that do. Since the distribution of particle velocities is uniform in orientation, the fraction of particles executing flights of distance $`vt`$ or more (where $`v`$ is the velocity of the particle) involve a fraction of order of $`W/t`$ of the particles occupying a channel at any one time, where $`W`$ is the channel width defined in Fig. 1. (Although we have defined $`v=1`$, we retain $`v`$ in this section to clarify when we are speaking of distances and when we are speaking of times.) The contribution to $`|r|^q`$ from these particles is

$$|r|^q(vt)^q(C/t)=Ct^{q1},$$
(12)
where $`C`$ is a constant depending on the size of the channel.

There are two points that can be fixed on the graph of $`\gamma _q`$ vs. $`q`$. The first is the zeroth moment, which by the conservation of probability is identically equal to one. Thus $`\gamma _0=0`$. The second point is fixed by the diffusion coefficient, which relates the second moment to time. It has been found  that $`r^2`$ is of the order $`t\mathrm{log}t`$ and so, consistent with the definition (3), ignoring factors of $`\mathrm{log}t`$ we have $`\gamma _2=1`$.

Next we argue that $`\gamma _q`$ is a concave up function of $`q`$. Invoking the Cauchy-Schwartz inequality $`((xy)x\times y)`$, for $`0ϵq`$

$$|r|^q=|r|^{(q+ϵ)/2}|r|^{(qϵ)/2}|r|^{q+ϵ}^{1/2}|r|^{qϵ}^{1/2}.$$
(13)
Thus

$$\mathrm{log}|r|^q\frac{\mathrm{log}|r|^{q+ϵ}+\mathrm{log}|r|^{qϵ}}{2}.$$
(14)
From (2), (3), and (14) we have that the graph of $`\gamma _q`$ versus $`q`$ is indeed concave up

$$\gamma _q\frac{\gamma _{q+ϵ}+\gamma _{qϵ}}{2}.$$
(15)
Finally, we show that Eq. (4) holds, i.e., that $`\gamma _q=q/2`$ for $`q2`$ and $`\gamma _q=q1`$ for $`q2`$. We have given lower bounds on $`\gamma _q`$, i.e. $`\gamma _qq/2`$ and $`\gamma _qq1`$. We have also pinned the value of $`\gamma _q`$ at two values: $`\gamma _0=0`$ and $`\gamma _2=1`$. The first lower bound, the two known values of $`\gamma _q`$, and the concavity condition (15) force $`\gamma _q=q/2`$ for $`q2`$. There is a trivial upper bound on $`\gamma _q`$, $`\gamma _qq`$. This upper bound, along with Eq. (14) implies that the slope of $`\gamma _q`$ never exceeds one. Since $`\gamma _2=1`$, $`\gamma _qq1`$ for $`q2`$. Therefore, for $`q2`$ both the upper bound and the lower bound coincide resulting in $`\gamma _q=q1`$ for $`q2`$.

### II.3 Non-uniform Initial Angular Distributions

The discussions of the previous two subsections relied on the existence of a uniform distribution for $`\theta _n`$ in the particle ensemble. This applies if one starts with a distribution uniform in angle and uniform in the space within a cell, in which case the particle ensemble will retain the uniform angular distribution. On the other hand, even if we have an initial distribution that is non-uniform, it will (under very general conditions) relax to the uniform distribution as time increases. The question we now address is whether this relaxation is fast enough to yield the same result for the temporal scaling of $`|r^q|`$, Eq. (11), as the initially uniform distribution.

To consider this question, we first study the relaxation of an initial particle distribution with no particles in a finite gap around the channel direction. For example, for the scalloped channel we consider the case where the initial particle distribution is uniform in the space within a cell $`𝒞`$ and

$$P(x,y,\theta ,t=0)=\{\begin{array}{cc}0\hfill & \text{for }|\theta |<\theta _o\text{ or }|\theta \pi |<\theta _o,\hfill \\ K\hfill & \text{ otherwise}\hfill \end{array},$$
(16)
where $`0\theta 2\pi `$ is the angle of the particle velocity vector with the $`y`$-axis, $`K`$ and $`\theta _o`$ are constants, and $`(x,y)`$ is in $`𝒞`$. Similarly, for the cases of Figs. 1(a)-1(c) we consider an initial distribution,

$$P(x,y,\theta ,t=0)=\{\begin{array}{cc}0\hfill & \text{for }|\theta |,|\theta \pi /2|,|\theta \pi |,|\theta 3\pi /2|<\theta _o,\hfill \\ K\hfill & \text{otherwise}\hfill \end{array}.$$
(17)
In particular, we focus on the behavior of $`P(x,y,\theta ,t)`$ for $`|\theta |`$ small, $`t`$ large and $`(x,y)`$ in a channel. (Similar results apply for (17) with $`|\theta \pi /2|`$ or $`|\theta \pi |`$ or$`|\theta 3\pi /2|`$ small.) For all cases shown in Fig. 1, we find a remarkable scaling behavior. Let $`\widehat{P}(\theta ,t)=Pdxdy`$ where the spatial integral is over a channel. Then, if we introduce the scaled variable $`\varphi =\theta t`$, we find that, in all the cases we have tested, the angular distribution function $`\widehat{P}(\theta ,t)`$ approaches a stationary form. That is,

$$\widehat{P}(\theta ,t)\stackrel{~}{P}(\varphi )\text{ as }t\mathrm{}.$$
(18)
We illustrate this with a numerical calculation on the scalloped billiard with $`180^o`$ arcs (Fig. 1(d1)) in Fig. 4. In generating this figure we follow the evolution of a large number of orbits initialized according to Eq. (16), and form the distribution $`\widehat{P}(\theta ,t)`$ using a histogram approximation. As shown, at successively larger $`t`$ the distribution approaches a time independent form $`\stackrel{~}{P}(\varphi )`$ where $`\stackrel{~}{P}(\varphi )`$ is zero for $`\varphi <\varphi _{\text{min}}`$, and, as $`\varphi `$ increases past $`\varphi _{\text{min}}`$, $`\stackrel{~}{P}(\varphi )`$ increases, asymptoting to a constant for large $`\varphi `$.

A similar numerical experiment for the Sinai billiard (Fig. 1(a)) yields the result in Fig. 5. In this case the large time distribution assumes the form,

$$\stackrel{~}{P}(\varphi )=\{\begin{array}{cc}C(\varphi /\varphi _o)\hfill & \text{ for }0<\varphi <\varphi _o,\hfill \\ C\hfill & \text{ for }\varphi >\varphi _o,\hfill \end{array},$$
(19)
where $`C`$ and $`\varphi _o`$ are constants. Furthermore, we numerically obtain this same form for all the other cases of Fig. 1, \[except for the case of the scalloped billiard with $`180^o`$ arcs (Fig. 1(d1)) which gives the result in Fig. 4\]. We explain the reason for the result (19) and why it does not apply for the billiard of Fig. 1(d1) subsequently, but before doing that we first show that these results for $`\stackrel{~}{P}(\varphi )`$ imply the applicability of Eq. (4) for the large time behavior of $`|r|^q`$.

In order to show that (4) applies consider the fraction $`\mathrm{\Xi }(t_o)`$ of particles that, at some large time $`t_o`$, have $`\varphi <\varphi _{}`$, where we take $`\varphi _{}=2\varphi _{\text{min}}`$ for the case of Fig. 1(d1) and $`\varphi _{}=\varphi _o`$ for the other cases. Noting that $`\varphi <\varphi _{}`$ implies $`\theta <\varphi _{}/t`$, we have that for large $`t_o`$

$$\mathrm{\Xi }(t_o)_0^{\varphi _{}/t_o}\stackrel{~}{P}(\theta t)d\theta =(t_o)^1_0^\varphi _{}\stackrel{~}{P}(\varphi )d\varphi =Kt_o^1,$$
(20)
where $`K`$ is a constant. Between times $`t_o`$ and $`2t_o`$ these particles experience flights of length $`vt_o`$. Hence for $`t=2t_o`$ these flights give a contribution to $`r^q`$ that is approximately $`(2K/t)(vt)^qt^{q1}`$. Thus the lower bound (12) still applies, and, by the reasoning in Sec. 2B, we again obtain (4).

We note that the asymptotic time dependence $`\stackrel{~}{P}(\varphi )`$ with $`\varphi \theta t`$ is marginal in the sense that, if the repopulation of an initially empty channel were slower (in the sense below), then (4) would not be recovered. In order to see this, consider the hypothetical case where an asymptotic distribution $`P(\varphi )`$ of the form in Fig. 4 or Fig. 5 is still approached, but with a self-similar $`\varphi `$ scaling given by $`\varphi _\alpha \theta t^\alpha `$. We have already considered the case $`\alpha =1`$, while $`\alpha <1`$ ($`\alpha >1`$) corresponds to slower (faster) filling in of the channel. If, for $`0<\alpha <1`$ (i.e. slow repopulation) we pursue the same reasoning as above for the $`\alpha =1`$ case, then we obtain the bound $`\gamma _q\alpha (q1)`$, and we can no longer conclude that (4) holds. For $`\alpha >1`$, we consider at time $`t_0`$ a range of angles $`\theta _+1/t_o>\theta >\theta _{}\varphi _{}/t_o^\alpha `$. Since $`\alpha >1`$, we have that $`\theta _+\theta _{}`$, and replacing $`\theta _{}`$ by zero does not alter the estimate for the contribution to $`r^q`$ from $`\theta _+>\theta >\theta _{}`$. Thus the lower bound estimate of Sec. 2B still applies for $`0<\alpha <1`$.

We now discuss the asymptotic forms $`\stackrel{~}{P}(\varphi )`$ illustrated in Fig. 4 for the scalloped billiard with $`180^o`$ arcs, and in Fig. 5 and Eq. (19) for the other cases.

First we discuss the scalloped billiard with $`180^o`$ arcs. A full treatment of the theory yielding $`\stackrel{~}{P}(\varphi )`$ in this case will be given elsewhere ; in the present paper, we will limit our discussion to the basic reason for the difference between this case and the other cases. In particular, we discuss why, for this case, $`\stackrel{~}{P}(\varphi )=0`$ in a finite interval about zero, $`\varphi <\varphi _{\text{min}}`$. From Fig. 1(d1) we see that after a long flight in the channel a particle will collide with the channel wall close to one of the cusp points where two arcs touch. For the case of $`180^o`$ arcs, the tangent to such a section of the channel wall is nearly horizontal. Thus, upon reflection $`\widehat{\theta }\text{min}(|\theta |,|\theta \pi |)`$ will still be small. In fact as shown in  by consideration of the geometry, $`\widehat{\theta }_{n+1}`$ on the $`(n+1)`$st reflection cannot change by more than a factor of 3 from $`\widehat{\theta }_n`$,

$$\widehat{\theta }_n/3\widehat{\theta }_{n+1}3\widehat{\theta }_n,\text{ for }\widehat{\theta }_n1.$$
(21)
Thus, if $`\widehat{\theta }_n`$ is small, $`\widehat{\theta }_{n+1}`$ is still relatively small. We can obtain the lower bound on $`\varphi `$ by considering the most extreme case where $`\widehat{\theta }`$ always decreases by 3 on every bounce,

$$\widehat{\theta }_{n+1}=\frac{1}{3}\widehat{\theta }_n,$$
(22)
$$t_{n+1}=t_n+W/(v\widehat{\theta }_n),$$
(23)
where $`W`$ is the channel width defined in Fig. 1(d1), $`v`$ is the particle velocity, and $`t_n`$ is the time of the $`n`$th reflection. Multiplying (22) by (23) we have

$$\widehat{\varphi }_{n+1}=\frac{1}{3}\widehat{\varphi }_n+W/v,$$
(24)
which, for large $`n`$, asymptotes to the solution $`\widehat{\varphi }=3W/(2v)`$. Thus $`\varphi _{\text{min}}=3W/(2v)`$, which agrees with our numerical solution Fig. 4 (see also ).

In contrast to the case of the scalloped channel with $`180^o`$ arcs, in the other cases shown in Fig. 1, the scattering of a long flight upon reflection from a channel wall leads to a much more drastic change in the angle of a particle’s velocity vector with respect to the channel axis. For example, for the case of the Sinai billiard, the angular deflection is typically of order $`\widehat{\theta }^{1/2}`$ which , for small $`\widehat{\theta }`$, is much larger than $`\widehat{\theta }`$. For the case of the scalloped channel with arcs of less than $`180^o`$, a particle moving nearly parallel to a channel axis is scattered by an angle of order one. Furthermore, after a large deflection, the orientation of the particle’s velocity vector is rapidly randomized by a succession of many reflections which, since the particle is no longer in a long flight, can occur in a relatively short time. These considerations lead us to a model for the cases in Figs. 1(a)-1(c) and 1(d2) in which we adopt the model hypothesis that, when a particle in a long flight suffers a collision with a billiard wall, the orientation of its velocity vector is randomly scattered with uniform probability density in $`[0,2\pi ]`$. We wish to determine the evolution from the initial condition (a) in Fig. 6 for the case $`\theta _{\text{max}}1`$.

This initial distribution is equal to the initial distribution (b) in Fig. 6 minus the initial distribution (c) in Fig. 6. The initial condition (b), which is uniformly distributed in angle, remains unchanged when it is evolved forward in time $`[\widehat{P_b}(\theta ,t)=\widehat{P_b}(\theta ,0)]`$, since it is an invariant distribution. Thus, to find the evolution from initial condition (a), we can determine the evolution from (c), and then subtract it from (b). The long time evolution from (c) can be found by considering the time at which particles are scattered. Consider, for example, the scalloped channel, Fig. 1(d2), and a particle with a small initial $`\theta _o`$. Suppose the particle is located in the channel at a distance $`\mathrm{\Delta }x`$ from the boundary of the channel with which it will collide \[left or right vertical dashed line in Fig. 1(d2)\]. If $`\mathrm{\Delta }x<vt\mathrm{sin}\theta _ovt\theta _o`$, the particle scatters; if $`\mathrm{\Delta }x>vt\theta _o`$, it does not scatter. Since the particles we are considering are in the channel, $`\mathrm{\Delta }x<W`$, every particle with $`\theta _o>W/(vt)`$ must have scattered at least once. We assume that $`t>W/(v\theta _{max})t_o`$. Since $`\theta _{\text{max}}`$ is small, the scattered particles contribute a small positive value of order $`\theta _{\text{max}}`$ to $`\widehat{P}_c(\theta ,t)`$ in $`0\theta 2\pi `$. Thus $`\widehat{P}_c(\theta ,t)`$ is small (i.e., of order $`\theta _{\text{max}}`$) for $`\theta >W/vt`$. For $`\theta _o<\mathrm{\Delta }x/(vt)`$, $`t>t_o`$, the particle has not yet scattered. Assuming that the initial spatial distribution of particles in the channel is uniform, the fraction of particles with initial angle $`\theta _o`$ that have scattered is $`\theta _ovt/W`$. Thus

$$\widehat{P}_c(\theta ,t)\{\begin{array}{cc}\widehat{P}_c(\theta ,0)(1\theta vt/W)\hfill & \text{ for }\theta <W/(vt),\hfill \\ 0\hfill & \text{ for }\theta >W/(vt)\hfill \end{array},$$
(25)
where we have neglected the small, order $`\theta _{max}`$, contribution to $`\widehat{P}_c(\theta ,t)`$ from scattered particles. Subtracting $`\widehat{P}_c`$ from $`\widehat{P}_b`$ as illustrated in Fig. 7, we obtain the time asymptotic form in Fig. 5 and Eq. (19).

## III Discussion

As mentioned in the introduction, Bleher  proved, up to some natural conjectures that “For any periodic configuration of scatterers with an infinite horizon the limit in distribution

$$\underset{t\mathrm{}}{lim}\frac{r(t)r(0)}{(t\mathrm{log}t)^{1/2}}=\eta $$
(26)
exists and $`\eta `$ is a Gaussian random variable.” (We have substituted our $`r(t)`$ for Bleher’s $`𝐱(t)`$ for the sake of notational consistency).

We take no issue with Bleher’s proof referred to above, but we emphasize its use of limit in distribution convergence. This kind of convergence is the least restrictive kind of convergence considered within the context of probability and statistics. A proof of convergence in distribution implies only the ability to calculate the expectation value of functions that remain bounded. The convergence required by a limit in distribution is only strong enough to allow calculation of the expectations of functions that remain bounded . Therefore the convergence of Bleher’s pdf is not strong enough to allow the moments of the distribution to be calculated. As a result of the strong weight $`|r|^q`$ puts on the tails of the distribution $`(|r|\mathrm{})`$, two different distributions with the same limit in distribution can have very different moments.

The fact that Bleher is able to accurately (according to our simulations) calculate the second and lower moments of the displacement distribution suggests that his result can be strengthened to “convergence in $`q^{th}`$ mean” which is satisfied for a sequence $`X_n`$ if the expectation value of $`|X_nX|^q0`$ as $`n\mathrm{}`$. Convergence in $`q^{th}`$ mean also implies that the expectation value of $`|X_n|^p`$ limits to the expectation value of $`|X|^p`$ for $`1pq`$ . Thus our results are consistent with convergence in $`q`$th mean to Bleher’s distribution for $`q=2`$, but rule out convergence for any higher value of $`q`$.

The inapplicability of Bleher’s result explains the discrepancy between Eqs. (2)-(4) and the result $`|r|^q(t\mathrm{log}t)^{q/2}`$ one would find by mistakenly calculating moments using Bleher’s pdf. Equations (2)-(4) also differ from normal diffusion, $`|r|^qt^{q/2}`$, as well as from the result suggested in .

In conclusion, our two main results are as follows:

(a) The moments $`r^q`$ scale as $`t^{\gamma _q}`$ with $`\gamma _q`$ given by Eq. (4) for any initial bounded distribution,$`|P(x,y,\theta ,0)|<K`$, that is zero outside some finite region (in particular, Eq. (4) still applies if the initial distribution has no particles with infinite flights).

(b) If the initial distribution has no particles in a $`\theta `$interval about a direction of infinitely long flight (say $`\theta =0`$), then $`\widehat{P}(\theta ,t)=_{}P(x,y,\theta ,t)`$ approaches a time-invariant scaling form $`\stackrel{~}{P}(\varphi )`$, where $`\varphi =\theta t`$ and $`\stackrel{~}{P}(\varphi )`$ is universally the same (Fig. 5) for the billiards of Figs. 1(a),1(b), 1(c), 1(d2), but is different (Fig. 4) for the billiard of Fig. 1(d1).

This work was supported by the Office of Naval Research (Physics) and by the National Science Foundation (award DMS 0104087). We thank L. Bunimovich and J. R. Dorfman for discussion.

## IV References

a. Department of Physics and Institute for Research in Electronics and Applied Physics.

b. Department of Mathematics and Institute for Physical Science and Technology.

c. Department of Electrical and Computer Engineering.

1. H. A. Lorentz, Proc. Amst. Acad. 7:438, 585, 604 (1905)

2. E. H. Hauge, Lecture Notes in Physics, Vol. 31, p. 337. (1974).

3. G. Gallavotti, Lecture Notes in Physics, Vol. 38, p 236 (1975).

4. Ya. G. Sinai, Funkts. Anal. Ego Prilozh. 13:46 (1979).

5. L. A. Bunimovich and Ya. G. Sinai, Commun. Math. Phys. 78:479 (1981).

6. P. Castiglione, A. Mazzino, P. Muratore-Ginanneschi, A. Vulpiani, Physica D 134:75 (1999).

7. R. Ferrari, A. J. Manfroi, W. R. Young, Physica D 154:111 (2001).

8. The numerical simulations of the four examples were carried out in a single cell, such as the cells shown in Fig. 2. The path of a particle in the infinite domain billiard is extracted by reflection of the particle trajectory at each straight wall of the cell. For infinite domain billiards with random displacements, the displacement from the center of the cell was altered every time the particle reflected from a straight wall (i.e. entered a new cell). For simplicity we altered the displacement even when the particle entered a previously visited cell; thus the simulation differs slightly from the fixed scatterer scenario, but the discrepancy affects only lower order terms. Likewise the orientation of the square was altered every time the particle reflected from a straight wall.

9. P. M. Bleher J. Stat. Phys. 66:315 (1992).

10. A. Papoulis, Probability, Random Variables, and Stochastic Process, 4th ed. (McGraw-Hill, 1994).

11. G.M. Zaslavsky and M. Edelman, Phys. Rev. E 56:5310 (1997).

12. K.C. Lee, Phys. Rev. Lett. 60:1991 (1988).

13. D. N. Armstead B. R. Hunt and E. Ott, to be published (2002).

14. P. Dahlqvist, J. Stat. Phys. 84:773 (1996).

15. P. Pollett MS308 Probability Theory web page,
http://www.maths.uq.edu.au/ pkp/teaching/ms308/summaries/conv/conv.html.


text:
# A generalized Numerov method for linear second-order differential equations involving a first derivative term

## 1 Introduction

The linear second-order differential equations of the type

$$y^{\prime \prime }(x)+g(x)y^{}(x)+f(x)y(x)=0$$
(1)
occur in various fields of physics. Here we shall mean the problem of the numerical self-consistent solution of the Hartree-Fock (HF) equations of motion deduced from the Skyrme energy functional  describing ground-state properties of atomic nuclei. For spherical nuclei these equations can be reduced to the form (1), the first derivative term arising due to the radial dependence of the nucleon effective mass. Usually the task is solved by the Runge-Kutta method. But this method is not the best one for the HF self-consistency procedure because it requires the interpolation of functions $`g(x)`$ and $`f(x)`$ in Eq. (1) between the grid points where function $`y(x)`$ is not calculated. Another well-known method (see, for example, Refs. ), which was proposed by B. V. Numerov in 1923, is free of the pointed difficulty but the first derivative term in Eq. (1) precludes from its immediate application. Several modifications of the Numerov method (NM) were developed  in order to include equations of the type (1) and more general nonlinear equations. Here another generalization of the NM is presented which is most suitable for the HF calculations  and yields the same degree of accuracy as the initial Numerov method.

## 2 Generalized linear Numerov method

Let us introduce notations: $`y_0=y(x_0)`$, $`y_\pm =y(x_0\pm h)`$ and analogously for $`g(x)`$ and $`f(x)`$, where $`x_0`$ is the fixed grid point, $`h`$ is the step length. Developing quantities $`y_\pm `$, $`y_\pm ^{}`$, $`y_\pm ^{\prime \prime }`$ in powers of $`h`$, we obtain

$`y_++y_{}`$ $`=`$ $`2y_0+h^2y_0^{\prime \prime }+{\displaystyle \frac{h^4}{12}}y_0^{(4)}+{\displaystyle \frac{h^6}{360}}y_0^{(6)}+O(h^8),`$ (2)
$`y_+y_{}`$ $`=`$ $`2hy_0^{}+{\displaystyle \frac{h^3}{3}}y_0^{\prime \prime \prime }+{\displaystyle \frac{h^5}{60}}y_0^{(5)}+O(h^7),`$ (3)
$`y_+^{}+y_{}^{}`$ $`=`$ $`2y_0^{}+h^2y_0^{\prime \prime \prime }+{\displaystyle \frac{h^4}{12}}y_0^{(5)}+O(h^6),`$ (4)
$`y_+^{}y_{}^{}`$ $`=`$ $`2hy_0^{\prime \prime }+{\displaystyle \frac{h^3}{3}}y_0^{(4)}+{\displaystyle \frac{h^5}{60}}y_0^{(6)}+O(h^7),`$ (5)
$`y_+^{\prime \prime }+y_{}^{\prime \prime }`$ $`=`$ $`2y_0^{\prime \prime }+h^2y_0^{(4)}+{\displaystyle \frac{h^4}{12}}y_0^{(6)}+O(h^6),`$ (6)
$`y_+^{\prime \prime }y_{}^{\prime \prime }`$ $`=`$ $`2hy_0^{\prime \prime \prime }+{\displaystyle \frac{h^3}{3}}y_0^{(5)}+O(h^5).`$ (7)
In addition, Eq. (1) yields:

$`y_0^{\prime \prime }`$ $`=`$ $`g_0y_0^{}f_0y_0,`$ (8)
$`y_+^{\prime \prime }`$ $`=`$ $`g_+y_+^{}f_+y_+,`$ (9)
$`y_{}^{\prime \prime }`$ $`=`$ $`g_{}y_{}^{}f_{}y_{}.`$ (10)
Let us consider Eqs. (3)–(10) as a system of eight linear equations for eight unknown quantities: $`y_0^{}`$, $`y_0^{\prime \prime }`$, $`y_0^{\prime \prime \prime }`$, $`y_0^{(4)}`$, $`y_\pm ^{}`$, $`y_\pm ^{\prime \prime }`$. Solving these equations and substituting the found quantities $`y_0^{\prime \prime }`$ and $`y_0^{(4)}`$ in Eq. (2) we get

$$T_0y_0=T_+y_++T_{}y_{}+R^{(6)},$$
(11)
where

$`T_0`$ $`=`$ $`2a{\displaystyle \frac{5h^2}{6}}b_0f_0,`$ (12)
$`T_\pm `$ $`=`$ $`a\pm {\displaystyle \frac{h}{24}}(10cg_0+g_++g_{})+{\displaystyle \frac{h^2}{12}}b_\pm f_\pm ,`$ (13)
$`a`$ $`=`$ $`\left(1+{\displaystyle \frac{h}{3}}g_+\right)\left(1{\displaystyle \frac{h}{3}}g_{}\right)+{\displaystyle \frac{h^2}{18}}g_0(g_++g_{}),`$ (14)
$`b_0`$ $`=`$ $`\left(1+{\displaystyle \frac{4h}{15}}g_+\right)\left(1{\displaystyle \frac{4h}{15}}g_{}\right)+\left({\displaystyle \frac{h}{15}}\right)^2g_+g_{},`$ (15)
$`b_\pm `$ $`=`$ $`\left(1\pm {\displaystyle \frac{5h}{6}}g_0\right)\left(1{\displaystyle \frac{h}{3}}g_{}\right)+\left({\displaystyle \frac{h}{3}}\right)^2g_0g_{},`$ (16)
$`c`$ $`=`$ $`\left(1+{\displaystyle \frac{7h}{20}}g_+\right)\left(1{\displaystyle \frac{7h}{20}}g_{}\right)+\left({\displaystyle \frac{3h}{20}}\right)^2g_+g_{},`$ (17)
$`R^{(6)}`$ $`=`$ $`{\displaystyle \frac{h^6}{240}}(y_0^{(6)}+3g_0y_0^{(5)})+O(h^8).`$ (18)
In more detail this result can be obtained by the following way. Making use of Eqs. (8)–(10), we get from Eqs. (7) and (6)

$`y_0^{\prime \prime \prime }`$ $`=`$ $`{\displaystyle \frac{g_{}y_{}^{}+f_{}y_{}g_+y_+^{}f_+y_+}{2h}}{\displaystyle \frac{h^2}{6}}y_0^{(5)}+O(h^4),`$ (19)
$`y_0^{(4)}`$ $`=`$ $`{\displaystyle \frac{2(g_0y_0^{}+f_0y_0)g_+y_+^{}f_+y_+g_{}y_{}^{}f_{}y_{}}{h^2}}`$ (20)
$`{\displaystyle \frac{h^2}{12}}y_0^{(6)}+O(h^4).`$
Substituting these equalities into Eqs. (3)–(5) we obtain

$`y_+y_{}`$ $`=`$ $`2hy_0^{}+{\displaystyle \frac{h^2}{6}}(g_{}y_{}^{}+f_{}y_{}g_+y_+^{}f_+y_+)`$ (21)
$`{\displaystyle \frac{7h^5}{180}}y_0^{(5)}+O(h^7),`$
$`y_+^{}+y_{}^{}`$ $`=`$ $`2y_0^{}+{\displaystyle \frac{h}{2}}(g_{}y_{}^{}+f_{}y_{}g_+y_+^{}f_+y_+)`$ (22)
$`{\displaystyle \frac{h^4}{12}}y_0^{(5)}+O(h^6),`$
$`y_+^{}y_{}^{}`$ $`=`$ $`{\displaystyle \frac{h}{3}}(4g_0y_0^{}+4f_0y_0+g_+y_+^{}+f_+y_++g_{}y_{}^{}+f_{}y_{})`$ (23)
$`{\displaystyle \frac{h^5}{90}}y_0^{(6)}+O(h^7).`$
It is useful to rewrite Eq. (21) in the form:

$`y_0^{}`$ $`=`$ $`{\displaystyle \frac{y_+y_{}}{2h}}+{\displaystyle \frac{h}{12}}(g_+y_+^{}+f_+y_+g_{}y_{}^{}f_{}y_{})`$ (24)
$`+{\displaystyle \frac{7h^4}{360}}y_0^{(5)}+O(h^6).`$
Substitution for this formula into Eqs. (22) and (23) leads to the following system of two equations for quantities $`y_\pm ^{}`$

$`\alpha _+y_+^{}+\alpha _{}y_{}^{}`$ $`=`$ $`u,`$ (25)
$`\beta _+y_+^{}\beta _{}y_{}^{}`$ $`=`$ $`v,`$ (26)
where

$`\alpha _\pm `$ $`=`$ $`1\pm {\displaystyle \frac{h}{3}}g_\pm ,\beta _\pm =\alpha _\pm +{\displaystyle \frac{h^2}{9}}g_0g_\pm ,`$ (27)
$`u`$ $`=`$ $`{\displaystyle \frac{1}{h}}(y_+y_{})+{\displaystyle \frac{h}{3}}(f_{}y_{}f_+y_+){\displaystyle \frac{2h^4}{45}}y_0^{(5)}+O(h^6),`$ (28)
$`v`$ $`=`$ $`{\displaystyle \frac{2}{3}}g_0(y_{}y_+){\displaystyle \frac{h}{3}}(4f_0y_0+f_+y_++f_{}y_{})`$ (29)
$`+{\displaystyle \frac{h^2}{9}}g_0(f_{}y_{}f_+y_+){\displaystyle \frac{7h^5}{270}}g_0y_0^{(5)}{\displaystyle \frac{h^5}{90}}y_0^{(6)}+O(h^7).`$
The solution of the system (25), (26) is as follows:

$$y_\pm ^{}=\frac{\beta _{}u\pm \alpha _{}v}{2a},$$
(30)
where $`a=\frac{1}{2}(\alpha _+\beta _{}+\alpha _{}\beta _+)`$, that coincides with the definition (14).

Substituting the formulas (30) into right-hand side of Eq. (24) we obtain after some algebra with taking into account Eqs. (27)–(29)

$$y_0^{}=\frac{S_+y_+S_{}y_{}S_0y_0}{2ah}+R^{(4)},$$
(31)
where

$`S_0`$ $`=`$ $`{\displaystyle \frac{h^3}{9}}(g_++g_{})f_0,`$ (32)
$`S_\pm `$ $`=`$ $`\left(1+{\displaystyle \frac{5h}{12}}g_+\right)\left(1{\displaystyle \frac{5h}{12}}g_{}\right)+\left({\displaystyle \frac{h}{12}}\right)^2g_+g_{}+{\displaystyle \frac{h^2}{6}}\left(1{\displaystyle \frac{h}{3}}g_{}\right)f_\pm ,`$ (33)
$`R^{(4)}`$ $`=`$ $`{\displaystyle \frac{7h^4}{360}}y_0^{(5)}+O(h^6).`$ (34)
Finally, substitution for the found solutions $`y_\pm ^{}`$, $`y_0^{}`$ (Eqs. (30), (31)) into Eqs. (8) and (20) yields the explicit formulas for the quantities $`y_0^{\prime \prime }`$ and $`y_0^{(4)}`$ in terms of $`y_0`$, $`y_\pm `$. After substitution for these formulas into Eq. (2) and a series of lengthy but straightforward algebraic transformations we arrive at the result (11)–(18). Omitting the term $`R^{(6)}`$ in Eq. (11) we obtain the recurrence three-point formula of the generalized Numerov method for linear second-order differential equations or, for brevity, of the generalized linear NM (GLNM). Clearly this method reduces to the ordinary NM if $`g(x)=0`$ in Eq. (1). The local truncation error of the GLNM, which is contained in the term $`R^{(6)}`$, is one of the same order $`h^6`$ as the error of the ordinary NM.

The formula (31) enables one to calculate the first derivative if the function $`y(x)`$ is known at the grid points. The local truncation error of order $`h^4`$ is determined by the term $`R^{(4)}`$. The more precise formula follows immediately from Eq. (23)

$$y_\pm ^{}=\frac{\left(3hg_{}\right)y_{}^{}h\left(4g_0y_0^{}+4f_0y_0+f_+y_++f_{}y_{}\right)3R^{(5)}}{3\pm hg_\pm },$$
(35)
where

$$R^{(5)}=\frac{h^5}{90}y_0^{(6)}+O(h^7).$$
(36)
The presence of the derivatives in the right-hand side is a shortcoming of this formula, nevertheless Eq. (35) is practical for the evaluation of $`y^{}(x)`$ at the endpoints of the grid.

## 3 Application to the Hartree-Fock calculations

Consider a general scheme within which the method proposed can be applied to the numerical solution of the HF equations. The HF approximation is a basis of numerous microscopic physical theories describing the quantum many-body systems. In general formulation, the HF method leads to a system of nonlinear integrodifferential equations. We shall consider a special case of the HF equations of motion deduced from the energy functional constructed on the base of the zero-range Skyrme forces  which are widely used for the description of atomic nuclei properties (see, for example, Ref. ). The variational principle applied to the Skyrme energy functional leads in the case of spherical nuclei to the following system of equations (in proper units):

$$z_\lambda ^{\prime \prime }(x)\frac{m_q^{}(x)}{m_q(x)}z_\lambda ^{}(x)+2m_q(x)[e_\lambda V_\lambda (x)]z_\lambda (x)=0,$$
(37)
where $`q`$ denotes a sort of nucleon (proton or neutron), the index $`\lambda `$ stands for the set of orbital quantum numbers (including $`q`$), $`x`$ denotes the radial coordinate, $`z_\lambda (x)`$ is the radial wave function, $`e_\lambda `$ is the eigenvalue playing the role of single-particle energy, $`V_\lambda (x)`$ is the state-dependent mean-field potential, $`m_q(x)`$ is the nucleon effective mass.

Comparing Eqs. (1) and (37), we see that they have the same form and would be identical if we put

$$y(x)=z_\lambda (x),g(x)=\frac{m_q^{}(x)}{m_q(x)},f(x)=2m_q(x)[e_\lambda V_\lambda (x)].$$
(38)
The essential difference consists in the following: actually formula (37) stands for the system of $`N`$ coupled nonlinear integrodifferential equations because the quantities $`m_q(x)`$ and $`V_\lambda (x)`$ are functionals of the densities which depend in turn on the set $`\{z_\lambda ,z_\lambda ^{}\}`$ of all the wave functions and their derivatives with $`\lambda \{\lambda _1,\mathrm{},\lambda _N\}`$ (see ). In practice, the system of equations (37) is solved by making use of some iteration procedure. The convergence is achieved by averaging of the densities calculated on two successive iterations. The description and the analysis of the procedure in more detail are outside the scope of the present paper (see, e. g., Ref. , and references therein where some relevant methods are discussed). Here it is important only that on each fixed HF iteration one has to solve the set of $`N`$ uncoupled linear differential equations which have the same form (37) but with already known functions $`m_q(x)`$ and $`V_\lambda (x)`$ determined by the results of previous iterations.

The numerical integration of Eq. (37) in this case can be performed by means of the GLNM described above. For the sake of simplicity it is convenient to come back to the notations of the preceding section taking into account Eqs. (38). Setting

$`Y_0^{out}`$ $`=`$ $`y_0/y_+,Y_{}^{out}=y_{}/y_0,`$ (39)
$`Y_0^{in}`$ $`=`$ $`y_0/y_{},Y_+^{in}=y_+/y_0,`$ (40)
we obtain the following recurrence relations from Eq. (11) (omitting the term $`R^{(6)}`$):

$`Y_0^{out}`$ $`=`$ $`T_+/(T_0T_{}Y_{}^{out}),`$ (41)
$`Y_0^{in}`$ $`=`$ $`T_{}/(T_0T_+Y_+^{in}).`$ (42)
Note that these transformations of Eq. (11) are similar but not identical to ones of the renormalized Numerov method developed in Ref. .

In the proposed scheme Eqs. (41) and (42) are used for the outward and the inward integrations, respectively. The outward integration starts from a point near $`x=0`$. The initial value of $`Y_{}^{out}`$ in Eq. (41) is calculated using analytic expansions of the regular solutions $`z_\lambda (x)`$ about $`x=0`$. To improve the accuracy of calculations it is practical to decrease the step length of the grid near this point. The inward integration starts from an outside endpoint of the grid where the irregular solutions of Eq. (37) are known analytically. The initial value of $`Y_+^{in}`$ in Eq. (42) is calculated as the ratio of the Whittaker functions for protons and of the spherical Hankel ones for neutrons. The eigenvalue $`e_\lambda `$ is found from the condition $`Y_0^{out}=1/Y_+^{in}`$ at some matching point $`x_0=x_m`$. The reasonable choice for this point is (see, e. g., Ref. ) the approximate position of the last extremum of the wave function $`z_\lambda (x)`$. In the end of this procedure, which is performed for each $`z_\lambda (x)`$ separately, the function $`z_\lambda (x)`$ is calculated at the grid points using the ratios $`Y_0^{out}`$, $`Y_0^{in}`$, and the normalization condition:

$$_0^{\mathrm{}}𝑑xz_\lambda ^2(x)=1.$$
(43)
Finally, the derivatives $`z_\lambda ^{}(x)`$ are calculated employing Eq. (31) at the inside points of the grid and Eq. (35) at the endpoints. After this, new approximations to the functions $`m_q(x)`$, $`V_\lambda (x)`$, which are used in the next HF iteration, are calculated.

The algorithm, that was briefly outlined above, has been realized in the computer code  intended for the Skyrme-Hartree-Fock calculations. The convergence and stability of the procedure described were tested in the calculations of ground-state properties of all doubly magic atomic nuclei using most of the present Skyrme-force parametrizations. It was obtained that the algorithm based on the GLNM reproduce the known reference results within their accuracy.


text:
# 1 In Euler problem as well as in the generalizations we present, a particle at the position (𝑥,𝑦) evolves under the action of two centers A and B, with respective coordinates (1,0) and (-1,0).

QUASI-INTEGRABILITY IN A CLASS OF SYSTEMS

GENERALIZING THE PROBLEM OF TWO FIXED CENTERS

A. ALBOUY

IMCCE

77, av. Denfert-Rochereau, 75014 Paris, France

albouy@bdl.fr

T.J. STUCHI

Instituto de Física, Universidade Federal do Rio de Janeiro

Caixa Postal 68528

21941-972, Rio de Janeiro, Brazil

tstuchi@if.ufrj.br

Abstract. The problem of two fixed centers is a classical integrable problem, stated and integrated by Euler in 1760. The integrability is due to the unexpected first integral $`G`$. Some straightforward generalizations of the problem still have the generalization of $`G`$ as a first integral, but do not possess the energy integral. We present some numerical integrations suggesting that in the domain of bounded orbits the behavior of these a priori non hamiltonian systems is very similar to the behavior of usual quasi-integrable systems.

The equations. Euler’s problem in the plane (see Figure 1) is defined by the system of differential equations

$$\ddot{x}=a(x_A,y)x_Ab(x_B,y)x_B,\ddot{y}=a(x_A,y)yb(x_B,y)y.$$
$`(1)`$
The two fixed centers are the points $`(1,0)`$ and $`(1,0)`$, and the moving particle is the point $`(x,y)`$. We have set $`x_A=x1`$, $`x_B=x+1`$,

$$a(\xi ,\eta )=m_A(\xi ^2+\eta ^2)^{3/2},b(\xi ,\eta )=m_B(\xi ^2+\eta ^2)^{3/2}.$$
$`(2)`$
The problem can be defined in the 3-dimensional space in the same way, and is also integrable, as was noticed by Euler. However, we will restrict ourselves to the planar case.

The first step in Euler’s integration was to exhibit two independent first integrals of the motion. One is the energy

$$H=(\dot{x}^2+\dot{y}^2)/2\stackrel{~}{a}(x_A,y)\stackrel{~}{b}(x_B,y),$$

with $`\stackrel{~}{a}(\xi ,\eta )=m_A(\xi ^2+\eta ^2)^{1/2}`$, $`\stackrel{~}{b}(\xi ,\eta )=m_B(\xi ^2+\eta ^2)^{1/2}`$. We will call the second one Euler’s integral:

$$G=C_AC_B2\stackrel{~}{a}(x_A,y)x_A+2\stackrel{~}{b}(x_B,y)x_B,$$

with $`C_A=x_A\dot{y}y\dot{x}`$, $`C_B=x_B\dot{y}y\dot{x}`$. Euler continued the integration, eliminating the second derivatives in (1) using the first integrals, and separating the variables.

Our generalization is simply to consider system (1) in the case where $`a`$ and $`b`$ are any homogeneous functions of degree $`3`$. Indeed, we want to put a little restriction on these homogeneous functions. We will suppose that both differential forms $`\xi a(\xi ,\eta )(\xi d\eta \eta d\xi )`$ and $`\eta a(\xi ,\eta )(\xi d\eta \eta d\xi )`$ are exact forms on the plane minus the origin, and that the same is true when we change $`a`$ in $`b`$. This hypothesis comes from the study of the problem of one fixed center (see  and ). It is not a strong restriction: the forms are already closed, so the condition on each function is just the cancellation of two scalar quantities, namely the integrals of both forms on a closed path around the origin.

It can be shown that any function $`a(\xi ,\eta )`$ satisfying the above conditions comes from a function $`A(\xi ,\eta )`$ homogeneous of degree $`1`$ as follows. Let us denote by $`A_\xi `$, $`A_\eta `$ the first derivatives of the function $`A`$ and by $`A_{\xi \xi }`$, $`A_{\xi \eta }`$, $`A_{\eta \eta }`$ the second derivatives. Then

$$a(\xi ,\eta )=\eta ^2A_{\xi \xi }=\xi ^1\eta ^1A_{\xi \eta }=\xi ^2A_{\eta \eta }.$$

The function $`a(\xi ,\eta )`$ in Euler’s case $`(2)`$ is obtained in this way from the function $`A=m_A(\xi ^2+\eta ^2)^{1/2}`$. We also have $`A_\xi =\xi \stackrel{~}{a}`$. It was discovered by the first author that $`G`$ persists in the form

$$G=C_AC_B2A_\xi (x_A,y)+2B_\xi (x_B,y),$$

where $`B(\xi ,\eta )`$ and $`B_\xi (\xi ,\eta )`$ are associated to the function $`b`$ in the same way as $`A`$ and $`A_\xi `$ are associated to $`a`$. In general, no integral takes the place of the energy integral.

Quasi-integrability. We report our numerical exploration of these generalized Euler’s problems, showing three examples that seem to us significant. In all cases we met, the result is either escape or quasi-integrable behavior. The third experiment displays some islands suggesting non integrability. Magnifying the neighborhood of a saddle point, a domain of irregular dynamics can be observed.

The obvious choice for a Poincaré section is to fix the integral $`G`$ and take, for example, $`y=0`$ $`(\dot{y}>0)`$. In each case, we show the iterates of some points of this Poincaré mapping, the central orbit in the section and some typical quasiperiodic orbit. All the orbits in a given field of forces have the same value of Euler’s integral. Since the examples have very large orbits, we have taken throughout the numerical experiments a somewhat arbitrary cut-off criterion given by the value of any coordinate or velocity greater than a thousand. In all the examples, we have taken $`A(\xi ,\eta )=(5\xi ^2\xi \eta +5\eta ^2)/10r`$, where $`r=\sqrt{\xi ^2+\eta ^2}`$, which corresponds to $`a(\xi ,\eta )=(5\xi ^2+3\xi \eta +5\eta ^2)/10r`$.

First example. Figures 2 to 4 correspond to $`B(\xi ,\eta )=(\xi ^33\xi \eta ^2)/4r^2+r`$, and thus $`b(\xi ,\eta )=2(3\xi \eta ^2\xi ^3)/r^6+r^3`$. The Poincaré map is close to a linear map, on a whole domain delimited by the escape criterion. This is quite strange. An explanation for this phenomena comes from geometrical considerations. We have chosen the plane as the domain for the motion, but there is a natural bigger domain for this kind of systems. It is the manifold of half lines drawn from the origin in a $`3`$dimensional vector space. Our plane is from this point of view just one half of the natural domain, a hemisphere chosen arbitrarily. Escaping orbits appear as orbits cut by the boundary of the hemisphere (in the classical Kepler problem, hyperbolas appear in the same way as cut ellipses). The theoretical grounds for this remark may be found in .

Second example. In Figures 5 to 7, we have chosen $`B=4(\xi ^4+\eta ^4)^{1/4}`$, $`b=12\xi ^2\eta ^2(\xi ^4+\eta ^4)^{7/4}`$. Here the section displays a wide domain with strong torsion but we are still very close from an integrable system. This rises the question: what are the integrable systems nearby? We know very few cases where our generalized Euler problem is integrable, namely the classical case and its projective transformations defined in  (which correspond for example to replace $`\xi ^2+\eta ^2`$ in Eq. (2) for $`b`$ by any homogeneous quadratic expression in $`\xi `$, $`\eta `$, and leave $`a`$ as it is.) Because we needed to get sufficiently many bounded orbits we were probably forced to stay close from integrable cases.

Third example and final comments. In Figures 8 to 11, we have chosen $`B=(3\xi ^2\xi \eta +3\eta ^2)/3r`$, $`b=(\xi ^2+\xi \eta +\eta ^2)/r^5`$. Here the system behaves as a typical conservative system close to an integrable system. We are more accustomed to observe this in the class of hamiltonian systems, and one can argue that maybe the system is hamiltonian for some symplectic form. We do not believe so, and rather relate this quasi-integrability to KAM theory applied to reversible systems (see , Theorem 2.9). Our systems are clearly reversible.

Acknowledgments. We thank Prof. H. Cabral and the Departamento de Matemática da Universidade Federal de Pernambuco for their hospitality, Prof. A. Neishtadt for stimulating comments, and Prof. Carles Simó for his rk78 code and his very useful comments.

References

 L. Euler, De motu corporis ad duo centra virium fixa attracti, Problème. Un corps étant attiré en raison réciproque quarrée des distances vers deux points fixes donnés trouver les cas où la courbe décrite par ce corps sera algébrique. Opera Omnia, S. 2, vol 6. (1764, 1765, 1760) pp. 209–246, 247–273, 274–293

 A. Albouy, Lectures on the two-body problem, Classical and Celestial Mechanics: The Recife Lectures, H. Cabral F. Diacu editors, volume “in press”, Princeton University Press

 G. Darboux, Sur une loi particulière de la force signalée par Jacobi, Note 11 au cours de mécanique de Despeyrous, tome premier, A. Hermann, Paris (1884)

 P. Appell, Sur les lois de forces centrales faisant décrire à leur point d’application une conique quelles que soient les conditions initiales, American Journal of Mathematics 13 (1891) pp. 153–158

 J. Moser, Stable and Random Motions in Dynamical Systems, With Special Emphasis on Celestial Mechanics, Princeton University Press (1973)