Supported by the Office of Naval Research, contract NONR-225(83) and AFOSR AF 49(638)1321.
Author's current address: Present address: NASA Langley Research Center, Hampton, Va.
A collisionless computer model is used to perform gravitational experiments with a two-dimensional galaxy. The motion of 2000 rods of mass is computed step-wise in time as they move according to Newton's laws of motion in their mutual gravitational field. The results are displayed as a movie and show the types of galactic structure that can arise from purely gravitational effects. A galaxy with insufficient rotation to balance gravitational attraction is seen to collapse and form filamentary condensations, then to change through structures similar to the spirals and barred spirals, and finally to settle down to a cylindrical form reminiscent of an elliptical galaxy.
Computer models have been extensively used in the study of plasmas and charged particle flow (Buneman 1959; Dawson 1962; Eldridge and Feix 1962, 1963; Burger 1964, 1965; Hockney 1966), and in this paper we apply some of the techniques developed for plasma applications to gravitational problems.
Computer simulation of stellar motions can be divided into two categories. In the first the exact motion of a small number of stars, of the order of 50, is followed and great care is taken to integrate accurately through close binary encounters (von Hoerner 1960, 1963; Aarseth 1963, Hayli 1967, Szebehely, Standish, and Peters 1967). Such a simulation is appropriate in the study of open star clusters and in understanding the formation of binary stars. On the other hand there are the models for galaxies of stars which were pioneered by P. O. Lindblad (1960a, b, 1962) and B. Lindblad (1961). In these models it is important to eliminate the effect of binary encounters since these are negligible in a real galaxy. It is in this latter application that the techniques which have been developed to simulate plasmas may be used, since in this case also collisional effects due to binary encounters are unimportant. A one-dimensional model of interacting sheets of mass has been used by Hohl and Feix (1966), as a direct development from the plasma model of Eldridge and Feix (1962, 1963). In this model 1000-2000 sheet stars are moved. Hohl (1967) has also used a two-dimensional model in which 500 rod stars are moved. In this paper we make use of the techniques developed by Hockney (1966) in a two-dimensional rod model of a plasma to make a simulation of a galaxy with 2000 rod stars. The model is used to perform some simple experiments concerning the stability and condensation of a rotating cylinder of gravitating masses. Some of the structures observed show a remarkable similarity to structures seen in real galaxies. There is evidence for an evolution in form from the irregular through the spiral and barred spiral to the elliptical, due to purely gravitational effects.
These calculations are necessarily exploratory, but it is hoped that they will aid our understanding of galactic structure and evolution. In particular it must be remembered that the results are two-dimensional and correspond to the motion in a plane of infinitely long rods of mass. The more realistic case of mass points moving in a plane seems, unfortunately, to be a much more difficult problem, at least in the case when the gravitational force is obtained from a solution of Poisson's equation as it is in this model.
In the computer model of the galaxy we store the positions and velocities of a large number of rod stars and advance them step-wise in time according to Newton's law of motion. The principal difficulty in such a simulation is the calculation of the force on each star. For example, if there are N stars, then we may find the force on any one of them by adding up the contributions from each of the other (N-1) stars. Since this must be repeated for each of the N stars, the total computer time for the force calculation by this direct-action technique is, for large N, proportional to N2. In practice we find that the direct-action technique is limited to a few hundred stars on computers such as the IBM 7090 and 7094 (see, e.g., Hohl 1967).
An alternative technique - which perhaps we might call the Poisson technique - is to solve for the gravitational potential using Poisson's equation and obtain the force on each star by differentiating the local gravitational field. If use is made of the special techniques that are available for the solution of Poisson's equation, it is possible to move several thousand stars by this technique because the force calculation is now proportional to N log N instead of N2.
In this method we limit the calculation of the fields to a square region in space and divide this into a square array of M cells (see Fig. 1). At the center of each cell is a mesh point at which the gravitational potential is to be calculated. A time step of the calculation procedes in stages as follows:
1. Mass distribution - The coordinates and velocities of each star are stored in the computer to the full precision of the machine (8 decimals on the IBM 7090). The coordinates of each star are examined and the mass of the star is spread uniformly over the cell in which the star resides. This amounts to associating the mass of the star with the potential mesh point at the center of the cell in which it resides. When the mass of all the stars has been apportioned in this way, one is left with a mass distribution over the M mesh points.
2. Poisson's equation - Poisson's equation.The five-point finite-difference form of Poisson's equation is now solved, giving as its solution the potential at each of the M mesh points, i.e., at the center of each cell. Poisson's equation is solved by the fast direct method developed by Hockney (1965), which obtains the solution on a 48 × 48 mesh (M = 2304) in 0.88 sec on an IBM 7090, with an error in the solution of the difference equations of less than 1 part in a million. The method is described in more detail in the next section when the problem of boundary conditions is discussed, but for the moment we note that the number of operations is proportional to M log M.
3. Field calculation - The force on each star is obtained by differencing the potential as follows. If a star lies anywhere within a cell it experiences the same field. This is taken, in the case of the horizontal field, to be the difference between the potential at the mesh point immediately to the right and the mesh point immediately to the left, divided by twice the mesh spacing, H. The use of a finite mesh and the above simple approximation to the field leads to a force law between masses in the model as shown in Figure 2. We note that for two masses which are separated by greater than 6H the force computed by the model is within 2 per cent of the r-1 force between line masses and rapidly gets closer as the distance increases. At closer separations the two force laws have entirely different characteristics. In the model the force reaches a maximum at a separation of one mesh distance and is zero if two mass rods are in the same cell. This behavior contrasts with the interaction of line masses which goes to infinity monotonically as the separation decreases. The law of interaction in the model is, however, similar to the law of interaction of two clouds of mass with a diameter equal to a mesh interval.
4. Newton's law - Newton's law of motion is integrated for a short time interval DT, during which time the force on each star is held constant. A new position and velocity is computed for each star and the cycle repeats at step 1.
The computer time taken for a complete cycle with 2000 rods is 3.0 sec on an IBM 7090. The cycle time increases to 3.3 sec if a movie film of the star positions is made.
In a model such as that just described one must make a balance between the resolution in space, which is determined by the number of space cells M, and N the number of stars that are moved. In models of star clusters such as those of von Hoerner (1960) and Aarseth (1963) in which it is important to follow the motion of individual stars as they experience close binary encounters with their neighbors, obviously the space resolution M must be much greater than N. However, in our case we wish to simulate the motion of stars in a galaxy. In this case it is known that the effect of close binary encounters between individual stars is quite negligible. A computer model which must perforce condense 108 or so stars into one representative rod would exaggerate the effect of close binary encounters by about the same factor if the exact r-1 force between line rods was used. However, if we eliminate the force of interaction between rods in the model when the rods are closer than the average interrod separation, close binary encounters can be eliminated. For the same reason Lindblad (1962) in his direct action galactic model also found it necessary to introduce a safety sphere around each star, inside of which there was no interaction with other stars. In the Poisson method such a cutoff of the force is achieved naturally by choosing the number of space cells to be approximately equal to the number of particles. If this is done the number of computer operations required to move N stars is proportional to N log N for large N. Since the long-range force of interaction is correctly represented in the model, there will be relaxation effects due to the cumulative effects of many small angle collisions. The effect of these collisions is discussed in IIb.
The key to the successful use of the Poisson technique is the existence of a fast direct method for solving Poisson's equation in simple geometries. For complete details reference must be made to the original paper (Hockney 1965), but we give here the salient features of the method.
The technique works best for the interior of a square or rectangle in (x,y) cartesian coordinates, in which case the mass distribution is first Fourier-analyzed in one of these directions, say, in x. If the boundary conditions on the potential are simple, such as zero value or zero slope specified on the boundaries parallel to the y-axis, then Poisson's equation becomes converted to a set of uncoupled ordinary differential equations for the harmonic amplitudes of the potential with right-hand sides containing the known values of the harmonic amplitudes of the mass distribution. These equations, or the tridiagonal matrix to which they are analogous, are readily solved. The actual potential is obtained by performing a Fourier synthesis of the harmonic amplitudes of the potential. Obviously one could also solve the problem by Fourier analyzing in both the x- and y-directions, but a count of the computer operations shows that a double analysis would take longer on the computer than the method advocated here. Furthermore a double Fourier analysis would not be applicable to (r,z) or (r,θ) coordinate systems whereas a single analysis may be done in the z- and θ-directions. The resulting ordinary differential equations in the r-coordinate, while more complex than those in the y-coordinate, are still quite tractable and readily solved.
This technique is only fast in practice if the Fourier analysis is performed in an efficient manner. If the number of mesh points is restricted to be a power of 2 or to be 3 times a power of 2, then the symmetry of the sine and cosine functions can be exploited to reduce drastically the number of arithmetic operations required. Such methods were first introduced by Runge (1903) and have been described by Whittaker and Robinson (1944). Recently a new and more general formulation has been given by Cooley and Tukey (1965). In any event a finite Fourier analysis of n data points need only take about 2.5 n log2 n operations. The total number of operations for the complete solution of Poisson's equation on a square mesh with a total of M mesh points then becomes 1.25 M log2M for large M. Iterative methods of solution such as successive overrelaxation require at least 5 M operations per iteration. Hence, for the case of 1024 mesh points, the iterative method must converge in less than three iterations if it is to be faster than the direct method. Even with a good guess available from the last time step it is not really credible that such a rate of convergence is possible. Since the number of iterations for a given convergence goes up as √M for relaxation methods or at best as log M for alternating direction methods (see, e.g., Varga 1962), this comparison will be no more favorable to the iterative scheme as M increases. We expect therefore to find the direct method using Fourier analysis to be quicker than alternative iterative methods wherever it can be applied. In one case where the same problem was solved both ways the direct method was found to be five times faster than successive overrelaxation (J. Gary, private communication).
The method just described will solve Poisson's equation in the interior of a square if the value of the potential is given round the boundary. To determine the values on the boundary we make use of the expansion:
where Z = X + iY are coordinates on the boundary and R = (X2+ Y2)½. L is an arbitrary constant conveniently taken to be the width of the square.
This expansion is valid provided none of the stars reach the boundary, and in practice we limit the summation to the first twelve terms. The coefficients in the expansion are obtained by taking the various moments of the mass distribution according to the definition
where zj = xj + iyj is the coordinate of the jth mass and there are a total of N rods of mass.
The success of this adjustment can be gauged by the extent to which angular momentum is conserved. In all the experiments reported here the maximum relative change in angular momentum was 2 per cent. Lack of energy conservation arises from the finite size of the time step as well as the boundary adjustment. The maximum relative change in total energy was 3 per cent. These changes occurred over runs of 500 time steps or approximately five galactic rotations.
The motion of stars in a galaxy is collisionless in the sense that the distribution of stars can be computed from the collisionless Boltzmann (or Vlasov) equation. In this case the motion of stars is determined by the smoothed-out potential distribution obtained from Poisson's equation. For example, Chandrasekhar (1942, chap. iii) estimates that the time between collisions (or relaxation time) is typically several hundred galactic rotations. It is important, therefore, to be sure that collisions are unimportant in the computer model even though there are only 2000 representative rods instead of the 1011 stars that might make up a galaxy.
We have already explained that the computer model prevents close binary encounters, but it remains to be shown that the long-range encounters are unimportant. We shall adopt a stochastic view of the collision process as was done by Hockney (1966).
Consider the gravitational field experienced by a test star at the center of a cylindrical distribution of model rods that are distributed with uniform average density. Since the macroscopic field on the test charge is zero by symmetry, any field experienced by the test star will be due to fluctuations caused by the discrete nature of the masses.
The field on the star is
where Ei is the gravitational field at the test star due to the ith star. The summation extends over all the N stars of the system.
The mean square fluctuating field is then
Expressing this as an average over space we have
where n is the density of rods per cm2; m is the mass per unit length of rod; L is the diameter of the system; H is the mesh distance. The lower limit of the integration recognizes the fact that the force between stars is zero if the stars are in the same space cell.
To estimate the collision rate we also need to know the correlation time, τc, of the field fluctuations. This is certainly greater than or equal to the time step DT since the field is held constant (i.e., perfectly correlated) during a time step. We also expect the fields to be uncorrelated after a time lag sufficient to allow all rods to traverse the system. This is a time characterized by the dynamic period of the system
We shall take the correlation time to be the dynamic period of the system.
The collision rate <v> due to fluctuations is
where v is a characteristic velocity of motion, which we take as L/(2τD).
The importance of collision may be assessed by comparing the collision time (τcoll = <v>-1) with the dynamic time τD of the system which is the same as the rotation period τrot for a galaxy of uniform density in which rotation balances gravitational attraction.
After some manipulation we get
This relation not only gives the familiar 1/N dependence of other estimates but also shows the smoothing effect caused by taking a finite cell size, H, in the solution of Poisson's equation.
In our case we have N=2000, L=24H giving
We shall be studying here phenomena over a time span between 1 and 10 galactic rotations and for these intervals the model may be considered collisionless.
The experiments to be described concern the stability and evolution of a cylindrical collection of mass rods (henceforth referred to as a galaxy of stars) which initially rotates clockwise as a solid body with angular velocity independent of radius. We consider the cases where the rotation is 0, 0.5, 1.0, 1.5 times that required to balance centrifugal force against gravitational attraction. There are no random velocities inserted in the initial condition. However, a slight initial randomness is introduced through the initial positions of the stars, which are uniformly dense on average but not in detail. The initial velocities are chosen such that there is no initial linear momentum, and this remains true since linear momentum is accurately conserved by the model. Such initial conditions are chosen to represent a vortex in the protogalactic medium, which has condensed sufficiently to be considered as an independent gravitating system.
The most dramatic display of the computed results is in the form of a 16-mm computer-made movie which shows the space and time development of the star system. In Figures 3-7 (Plates 2-6) we have attempted to show the essential features of the movie by printing a fairly large number of frames, most of which are taken at equal intervals of time. However certain features such as individual star velocities and orbits do not show up in these figures and these must be left to the description.
Figure 3 (Plate 2) shows the time development of a rotating cylindrical galaxy of stars in which the initial rotation is such as to balance the gravitational attraction. For the case of constant density, this condition is satisfied if the cylinder is given a solid-body rotation with angular frequency:
where ρ is the constant mass density of the cylinder. A characteristic time of the problem is then the rotation period
which is used as the unit for time measurements. The time step is taken to be approximately one hundredth of this period. In the initial condition (t = 0) mass rods are distributed at random within a circle and given an angular rotation of ω*rot. The small deviations from a constant density, due to the random distribution of discrete rods, provide small initial disturbances to a perfectly balanced state from which instabilities may grow. Between t=0 and t=0.8 surface waves grow, leading to a distinctly fluted shape at t=0.8. If we describe the surface perturbation in radial coordinates as ∞ cos kθ, where θ is the azimuthal angle, then these flutes correspond to a wavenumber k, of 4 or 5. After t=1.0 the wavenumber k=2 begins to dominate the disturbance. The cylinder now resembles an egg rotating at ω*rot. The amplitude of the k=2 mode varies periodically, resulting in the galaxy having maximum eccentricity at t≈3.0 and t≈7.0 and becoming almost circular again at t≈5.0 and t=9.0.
The k=1 harmonic remains at constant amplitude since this follows from the absence of linear momentum in the initial condition. Harmonic k=2 appears as a standing wave with respect to a frame rotating with angular velocity ω*rot. Similarly harmonics k=3 and k=4 behave as standing waves with respect to frames rotating at 2ω*rot and 3ω*rot, respectively. The measured frequencies of the standing waves have been fitted empirically to the approximate dispersion relationship.
The balanced rotating cylinder is thus unstable, but non-linear effects limit the amplitudes of the disturbances to about a fifth of the initial radius. There is no evidence (at least over the time scale of ten rotations) of any tendency for the cylinder to fission and divide into two lumps. This confirms, in two dimensions, the findings of Lyttleton (1953), who concludes that such a fission process will not occur in a rotating gravitating fluid and therefore cannot be the origin of double stars, as was proposed by Jeans (1919, p102) and others.
Attempts to predict this surface instability and dispersion relation from cold fluid perturbation theory have proved unsuccessful and we present this result here as an experimental observation. It is reassuring that similar surface flutes are also observed by Hohl (1967) in a model which uses an entirely different method for the calculation of the fields. It is most likely therefore that these flutes represent a physical phenomenon and are not due to some mathematical instability.
Figure 4 (Plate 3) shows the time development of the collapse of a cylindrical galaxy of stars with no initial angular momentum. Between t=0 and t=0.2 there is a radially directed implosion. The collisionless nature of the model, however, allows the stars to pass through each other and expand again to about the initial radius, during which time condensation occurs. At t=0.4 five distinct condensations are visible which after a further implosion and expansion have, at t=0.8, become concentrated into three globular lumps. By t=2.0 these three lumps have become fused into a single condensation which represents a steady state. This state is one in which the gravitational attraction is balanced by the temperature (or random motion of the stars about the center).
In the movie one can distinguish two classes of stars. The vast majority of the stars are describing oscillatory motion with a short period in the main mass of the condensation, or nucleus. Some fifty, which might be described as the halo stars, have a higher energy and describe oscillatory motion out to about four or five times the radius of the nucleus. The period of oscillation of the halo stars is some ten times that of the stars in the nucleus. The separation of the stars into two classes is clearly seen even as early as t=0.6, when the majority of stars have imploded for the second time into a small radius while the halo stars are left scattered out as far as twice the radius of the initial cylinder.
Figure 5 (Plate 4) shows the collapse of a cylindrical galaxy which has initially half the angular velocity necessary to balance gravitational attraction. In the time development we see a mixture of the effects noted in IIIa, IIIb. There is now a tendency for the cylinder to pulsate between a maximum and minimum radius as seen at t=0.4 and t=0.8, respectively. The tendency to condense into smaller pieces is most pronounced when the maximum radius is reached and the velocities are small. These subsidiary condensations as seen at t=0.6 and t=1.0 are filamentary in shape rather than the globular condensations seen in the absence of rotation. Between t=1.4 and t=1.8 barred spiral structures are evident, but these seem to be shortlived; from t=2.0 to t=4.4 there is no evidence of spiral structure, and a steady state seems to have been set up which is roughly circular and is the two-dimensional analogue of an elliptical galaxy.
In order to show the structures in more detail we have performed the case of half the necessary rotation with twice the resolution in space (there are now 32 space cells across the cylinder as compared with 16 previously) and with half the time step. This is shown in Figure 6 (Plate 5). The motion is followed from t=0.0 to t=2.4 and follows the same general development as that of Figure 5, but the greater space resolution enables more subsidiary condensations to occur (about ten at t=0.6), and the spiral structures are more pronounced. This is because the model divides the square region into a 48 × 48 mesh or 2304 space cells and is such that the force between mass rods vanishes if both rods are in the same space cell. Thus no condensation can be held together that is smaller than three or four cells in diameter, which is about the size of the smaller condensations at t=0.6. The filamentary nature of these condensations appears to indicate the presence of significant rotation and is reminiscent of the condensations seen in red-light photographs of the Crab Nebula, NGC 1952. Filamentary condensations are also seen in some peculiar galaxies; for example, compare Figure 6 at t=1.1 or Figure 5 at t=1.0 with the galaxy NGC 1741 (Arp 1966, No. 259, pl. 44).
Between t=2.0 and t=2.2 one can see the evolution of a spiral structure into a barred structure by purely gravitational effects. The structure at t=2.0 resembles, for example, the Sc galaxy NGC 5457 (M101, Sandage 1961, pl. 27) and that at t=2.2 the SBb galaxy NGC 5850 (frontispiece of Chandrasekhar 1942, pl. II; Hubble 1936). Examination of the movie shows that the motion of the stars in the spiral arms at t=2.0 is predominatly longitudinally along the arm, from the outside inward. Thus the arm is roughly the orbit of the stars which compose the arm. One can also distinguish a halo of stars which execute eccentric orbits with a long period. These orbits pass close to the center of the nucleus and reach out to the edge of the square region. The majority of the stars move in much more circular orbits in the nucleus of the condensation.
Figure 7 (Plate 6) shows the case of a small cylinder which is given initially 1.5 times the angular velocity necessary to balance gravity. This is similar to the case of Figure 5 if we take t=0.2 as the origin of time, and shows the same general evolution from filamentary condensations, through spirals and barred spirals to an elliptical steady state. One can see a resemblance between the S-shaped condensation at t=4.8 and the SBc galaxy NGC 7479 (Chandrasekhar frontispiece, Hubble pl. II) or NGC 3359 (Sandage 1961, pl. 49). This form is also seen in Figure 5 at t=1.4.
The normalized time measurements change rapidly in this case due to the higher density in the initial state and the resulting small unit of time.
A number of simple gravitational experiments have been performed with a computer model in which the motion of 2000 rods of mass is calculated in their mutual gravitational fields. The model is shown to be collisionless and may therefore be taken to represent a two-dimensional galaxy of stars. An initial cylindrical distribution, with a solid-body rotation just sufficient to balance gravitational attraction, is observed to be unstable; however, the instability does not reach sufficient amplitude to cause a fission of the cylinder. The case of no initial rotation shows several globular condensations which finally merge into a single condensation in which temperature balances gravitational attraction. The case in which either too much or too little rotation is given to the cylinder shows an evolution of forms from filamentary condensations, through spiral and barred spiral to the two-dimensional equivalent of an elliptical galaxy. Thus if gravitational effects dominate the evolution of galaxies, these experiments would indicate that the barred galaxies are young and the elliptic galaxies old, as is generally believed from considerations of stellar content (Sandage 1961). There is evidence for the stars to divide into two classes; a small number of halo stars with eccentric orbits of large amplitude and long period, and the majority of stars forming a compact nucleus with relatively circular orbits with a short period.
The author would like to thank Professors P.A. Sturrock and O. Buneman for encouraging the application of the computer model to gravitational problems and to the Office of Naval Research for providing the funds for the initial work under contract Nonr-225(83). Thanks are also due to the staff at the Stanford University Computation Center and the Lockheed Stromberg-Carlson facility, where the calculations were performed.
Aarseth, S J, 1963, Dynamical evolution of clusters of galaxies, Monthly Notices of the Royal Astronomical Society, Vol 126, p223
Arp, H, 1966, Atlas of Peculiar Galaxies, California Institute of Technology.
Buneman, 0, 1959, Dissipation of Currents in Ionized Media, Phys Rev 115, 503.
Burger, P, 1964, Non-existence of dc States in Low-Pressure Thermionic Converters, J Appl Phys, 35, 3048.
Burger, P, 1965, Theory of Large-Amplitude Oscillations in the One-Dimensional Low-Pressure Cesium Thermionic Converter, J Appl Phys, 36, 1938.
Chandrasekhar, S, 1942, Principles of Stellar Dynamics, University of Chicago Press.
Cooley, J W, Tukey, J W, 1965, An algorithm for the machine calculation of complex Fourier series, Math Comp, 19, 297-301.
Dawson, J M. 1962, One-Dimensional Plasma Model, Phys Fluids, 5, 445.
Eldridge, O C, Feix, M R, 1962, One-Dimensional Plasma Model at Equilibrium, Phys Fluids, vol 5, 1076.
Eldridge, O C, Feix, M R, 1963, Numerical Experiments with a Plasma Model, Phys Fluids, vol 6, 398.
Hayli, A, 1967, Proc Symp Computer Simulation of Plasma and Many-Body Problems, Williamsburg, Va, April.
Hockney, R W, 1965, A Fast Direct Solution of Poisson's Equation Using Fourier Analysis, JACM, vol 12, 95.
Hockney, R W, . 1966, Computer Experiment of Anomalous Diffusion, Phys Fluids, 9, 1826.
von Hoerner, S, 1960, Die numerische Integration des n-Korper-Problemes fur Sternhaufen I, Zs f Ap, vol 50, 184-214.
von Hoerner, S, 1963, Die numerische Integration des n-Korper-Problemes fur Sternhaufen II, Zs f Ap, vol 57, 47-82.
Hohl, F, 1967, Collective Effects in Stellar Dynamics and Plasma Physics, Ph.D. thesis, College of William and Mary.
Hohl, F, Feix, M R. 1966, A one-dimensional model to study collective effects in stellar dynamics, Phys. Letters, 22, 432-434.
Hubble, E, 1936, The Realm of the Nebulae, Yale University Press.
Jeans, J H, 1919, Problems of Cosmogony and Stellar Dynamics, Cambridge University Press.
Lindblad, B, 1961, On the formation of dispersion rings in the central layer of a galaxy, Stockholm Obs Ann, 21, 8.
Lindblad, P O, 1960, The development of spiral structure in a galaxy approached by numerical computations, Stockholm Obs Ann, 21, 4.
Lindblad, P O, 1960b, Stockholm Obs Ann, 21, 3.
Lindblad, P O, 1962, in The Distribution and Motion of Interstellar Matter in Galaxies, ed. L. Woltjer (New York: W. A. Benjamin, Inc.), pp. 222-233.
Lyttleton, R A. 1953, The Stability of Rotating Liquid Masses, Cambridge University Pres.
Runge, C, 1903, Uber die Zerlegung empirisch gegebener Periodischer Funktionen in Sinuswellen, Zs f Math u Phys, 48, 443.
Sandage, A, 1961, The Hubble Atlas of Galaxies, Carnegie Institute of Washington Pub, 618. Washington, D.C.
Szebehely, V G, Standish, E M, Peters, C F, 1967, Proc Symp Computer Simulation of Plasma and Many-Body Problems, Williamsburg, Va., April.
Varga, R S, 1962, Matrix Iterative Analysis, Englewood Cliffs, N.J.: Prentice-Hall).
Whittaker, E T, Robinson, G, 1944, Calculus of Observations, Glasgow, Scotland: Blackie & Sons).