A two-dimensional model is used to perform computer experiments for a collisonless self-gravitating system. The gravitational field is obtained by solving the Poisson equation and the system is advanced stepwise in time. Computer simulations have been performed for systems with an initially uniform distribution over a circular or elliptical region in x-y space, zero thermal velocity, and with given values of initial solid-body rotation. Up to 10,000 stars are used in the calculations. As such systems evolve in time they display the various shapes observed in actual galaxies. It is found that some of the calculations recently reported cannot be considered collisionless.
Actual stellar systems, such as galaxies, contain about 1011 stars of sufficiently large average separation so that binary encounters between stars can be neglected. Stellar systems can therefore be described by the collisionless Boltzmann equation. In order to simulate the evolution of such stellar systems on a computer one would like to follow the motion of at least several thousand masses or stars. Hohl and Feix (ref1 and ref2) and Lecar (ref3) used one-dimensional sheet models to study the evolution of self-gravitating systems. A similar model where the motion of a large number of concentric spherical mass shells is followed has been used by Henon (ref4 and ref5). Recently Hockney (ref6) and Hohl (ref7 and ref8) introduced a two-dimensional model where the stars are represented by infinitely long mass rods. In the present paper the results of some simple experiments with the two-dimensional model are given. The calculations show that spiral and other structures similar to those of some actual stellar systems can be obtained by purely gravitational effects.
The model consists of a large number of rods of mass m per unit length and the rods are of infinite extent in the z-direction. These rods move in the x-y plane under the action of their mutual gravitational attraction. The system of mass rods is advanced in time in the following manner. First, the distribution of mass (x,y) is used to obtain the gravitational potential p(x,y) by numerically solving the Poisson equation. Second, the gravitational field at the position of the particles is computed from the potential p(x,y). Third, using Newton's laws the motion of all the mass rods is advanced for a constant time step δt. These three steps represent one cycle and they are repeated until the desired evolution of the system is achieved.
The crucial point in the computations is the solution of the Poisson equation. It is desirable that the time required for this process be only a small fraction of the cycle time. If the system is advanced for a small time step δt the mass distribution (x,y) will not change very much. The change in the gravitational potential will then also be very small. Thus, the solution of the finite difference form of the Poisson equation by an iteration method which uses the potential from the previous cycle as an initial guess will converge very rapidly. The accuracy of the iterative solution of the Poisson equation is continually checked during the calculations. This is done by obtaining the values of the potential and the field at several selected points by summing directly the contribution from each star. The values so obtained agree with those obtained from the solution of the Poisson equation to at least the first three digits. For a 10,000 star system with a 101 × 101 grid to solve the Poisson equation the cycle time is 7 seconds on the Langley Research Center's CDC 6600 data processing system. This time includes such operations as the checking of the gravitational potential, the calculation of the energy and angular momentum, and writing the positions and velocities of all stars on tape.
To solve the Poisson equation one first requires the boundary conditions around the rectangular mesh. The potential at an arbitrary boundary point z = x + iy can be written as
where zn,m=xn,m+iyn,m is the coordinate of the cell n,m of the rectangular mesh, ρn,m is the mass density in that cell and G is the gravitational constant. Equation (1) can be written as
where N is the total number of stars in the system. Since ρn,m is nonzero only for | zn,m / z | < 1 equation (2) can be written as
where
and the series expansion for ln(1 - zn,m/z) has been truncated after 15 terms.
The Poisson equation
is solved by using the standard five-point difference equation
where Δx = Δy = 1 so that ρn,m represents the number of mass rods in the cell n,m. This set of simultaneously equations is solved by an iteration mechanism on the CDC 6600 data processing system in the form
The superscript r refers to the rth iteration and the parameter γ is adjusted to give the maximum rate of convergence.
The differential equation of motion for the stars are
The evolution in time of the particle trajectories is given by the finite difference approximation
Equations which were found to be as accurate as equations (10) and (11) but required fewer computer operations are
Similarly, the use of the new velocity in equation (14) removes an unconditional instability.
The kinetic energy of the system is
and the potential energy is
where the summation goes to N, the number of stars in the system. The total energy of the system U remains constant and is
U = T + P (17)
The relaxation time (ref 9) for the two-dimensional model is of interest to determine for how long the system can be treated as collisionless.
Consider an encounter between a rod of mass m and velocity V with a stationary rod of mass M. If D is the distance of closest approach (impact parameter) then the transverse force acting on the moving star during a time 2D/V can be approximated as 2GmM/D. The moving star has therefore acquired a transverse velocity
δV ≈ 4GM/V (18)
For a star interacting with a system containing many stars the effect of many individual encounters must be summed. The number of encounters in a time t with impact parameter between D and D + dD is tVn dD where n is the density of stars in the system. The expectation value of Vi2 is then given by
since in general Dmax ≫ Dmin where Dmax is the dimension of the system. The relaxation time τc is defined to be the time required for <VT2> to be of the same order of magnitude as v2. The velocity V is taken to be the characteristic velocity Vg = Dmax/τg where τg = 2π / √(2πGmn) is the rotation period of a uniform circular distribution of mass rods with the centrifugal force equal to the gravitational attraction. Equation (19) then becomes
where M = m. The ratio of the relation time to the rotation period is then given by
Equation (21) is only a rough estimate but it nevertheless indicates that for N equal to a few thousand, collisional effects can be expected to become important after only a few rotations.
Computer simulations have been performed for a number of systems with an initially uniform distribution over a circular or elliptical region in x-y coordinate space, zero thermal velocity, and various values of initial solidbody rotation. In particular, systems with a ratio of minor axis a to major axis b equal to 3, and with an initial solid-body rotation given by
have been investigated in detail. It was found that for a system containing 2000 stars with a 51 × 51 grid to solve the Poisson equation, collisional effects caused sufficient randomization to destroy the spiral structure after only 2 to 3 rotations. The results are shown in figure 1, figure 2 and figure 3. The normalization G = m = 1 were used in the calculations and the time is shown in rotation periods τg = 2π/ωg. Figure 1 shows the evolution of the 2000-star system in x-y coordinate space. The evolution of the velocity distribution in Vx - Vy space is shown in figure 2. Figure 3 shows the evolution of the velocity distribution after the initial solid body rotation has subtracted.
The radial velocity Vr is plotted versus (Vθ - rωg) where Vθ is the azimuthal velocity and r is the radius from the center of the system to a star. Initially the thermal (or random) velocity of all mass rods is zero. As the system evolves in time the thermal velocity builds up and the system expands in Vr - (Vθ- rω) space. After a time t = 1.75τg there is little further change in the velocity distribution. Thus, the calculations reported by Hockney (ref 6) with 2000 stars and a 48 × 48 grid cannot be considered collisionless for more than 2 rotations. The same holds for the results given by Hohl (ref 8) for 2000 stars (with a 51 × 51 grid). To reduce collisional effects the calculations were repeated with an identical initial elliptical distribution but with 4000 stars and a 101 × 101 grid to solve the Poisson equation. The results are shown in figure 4. and figure 5. Figure 4 shows that the spiral structure now has a somewhat greater stability but it still tends to randomize after about 4 rotations. The corresponding evolution in Vx - Vy velocity space is shown in figure 5. The number of stars was increased to 10,000 and the resulting evolution in x-y coordinate space is shown in figure 6. The evolution of yet another system with 20,000 stars was computed up to four rotations and it resulted in a structure nearly identical to that shown in figure 6. Thus, it appears that in order to simulate collisionless systems for the first n rotations, nearly 2n×103 stars are required. Figure 1 shows that at t = 6τg the spiral structure begins to randomize. To keep the system collisionless for a longer time calculations are now being performed with systems containing up to 200,000 stars.
An interesting point is that all systems investigated which have an initially elongated distribution and a given solid-body rotation go through an initial evolution similar to that shown in figure 6 up to t = 1.0τg. At t = 0.75τg the system shows a structure nearly identical to that of the peculiar galaxy (GB1) discussed by Burbidge, Burbidge, and Shelton (ref 10). Another such integral-sign shaped galaxy is NGC 4656/7 reproduced in the Hubble Atlas (ref 11). Figure 6 shows that at t = 3.0τg the two spiral arms begin to overlap and give the appearance of a circular ring separate from the nucleus, presenting a shape similar to NGC 3145 (ref 11). At t = 4.0τg the system shows a structure similar to NGC 5194 (ref 11). An interesting point is that at t = 3.0τg the spiral structure has almost disappeared and at t = 3.75τg, after less than one rotation, the spiral structure is again very pronounced. This phenomena appears to reoccur from t = 4.5τg to t = 6.00τg. However, at t = 6.0τg collisional effects (or computer noise) have already appreciably affected the evolution of the system. The author expects that for the calculations presently in progress with 100,000 stars, collisional effects will be small enough so that several such transitions from a ring shape to a spiral structure can be observed. Such repeating density waves resulting in spiral structure will prevent the arms from winding up. This overall picture is in agreement with the theory of Lin and Shu (ref 12) where the stars pile up for a certain time in the spiral-shaped potential wells. The evolution of the velocity distribution in Vx - Vy space is shown in figure 7. At t = 0.25τg, the system shows an interesting structure. Condensation in velocity space has occurred and striation appears in the velocity distribution. The corresponding evolution in Vr - (Vθ - τω) velocity space is shown in figure 8. This distribution shows little change after the first rotations.
One might expect that the appearance of the two-arm spiral structure shown in figure 6 depends on the initially elongated distribution. This is not the case as is illustrated in figure 9. Figure 9 shows the evolution of a 10,000-star system with an initially uniform distribution over a circular region and with an initial solid body rotation equal to that required to balance gravitation, that is, ω = ωg = 4ωGNm For the first two rotations the evolution of the balanced cylinder is very similar to that of the 2000-star system presented by Hockney (ref 6) and Hohl (ref 8). Figure 10 shows the evolution of such a 2000-star system in x-y coordinate space. For the 2000-star system it is found that the relatively small fourth mode surface wave around the periphery of the cylinder quickly change to a second mode surface wave (egg shape), then all structure is lost and the system finally assumes a circular shape. Figure 9 shows that the evolution is quite different when collisional effects are reduced. The fourth mode surface wave now continues to grow and after five rotations a four-arm spiral begins to form. After about 8 rotations two of the arms have disappeared and the system has assumed a bar shape. One would now expect that a two-arm spiral structure would form. However, figure 11 and figure 12 shows that after about nine rotations collisional effects have increased the random or thermal velocity to such an extent that they are of the same magnitude as the rotational velocities and the model no longer simulates a collisionless system. Nevertheless, figure 9 shows that the second mode surface wave will finally dominate. figure 11 shows that there is little change of the distribution in Vx - Vy, velocity space for the first 5 rotations. After that time condensation in velocity space begin to occur. The evolution of the random velocity shown in figure 12 shows that there is an almost linear increase in the random or thermal motion of the mass rods.
It would be more appropriate to simulate spiral galaxies by means of a model using finite length mass rods or even point masses moving in a plane. The two-dimensional model does simulate spiral structure and shows that two-arm spirals are very likely to occur. However, whether the two-dimensional model correctly simulates some of the more refined properties of spiral galaxies can only be determined by repeating the calculations with a model of finite length rods or point masses. The reason for using the two-dimensional model is that the solution of the Poisson equation is relatively simple. Since it was found that the number of stars required to simulate a collisonless system for a sufficiently long time is of the order of 100,000, the time used for solving the Poisson equation represents only a very small portion of the cycle time. The additional time required for obtaining the potential distribution for point masses will, therefore, no longer be an important factor. For this reason future calculations will be performed for systems of finite length rods and point masses moving in a plane.
[1] Hohl, F, Feix, M R, A One-Dimensional Plasma Model for a Gas of Stars, Bull. of the American Physical Soc., Vol 11, No 4, p. 560, 1966.
[2] Hohl, F, Feix, M R, Numerical Experiments with a One-Dimensional Model for Self-Gravitating Star System, Astrophys J, Vol 147, No 3, p 116l4, 1967.
[3] Lecar, M, Relaxation of a One-Dimensional Self-Gravitating Gas, Symp. on Computer Simulation of Plasma and Many-Body Problems, Williamsburg, Va, April 19-21, 1967. NASA SP-153, 1967.
[4] Henon, M, Initial Collapse and Dynamical Mixing of a Spherical Cluster. Gravitational Instability and the Formation of Stars, Galaxies and Their Characteristic Structures, 14th International Astrophysical Symp, Univ. of Liege, Belgium, June 1966.
[5] Henon, M, Collective Motions in a Spherical Star Cluster. Symp. on Computer Simulation of Plasma and Many-Body Problems, Williamsburg, Va, April 19-21, NASA SP-153. 1967.
[6] Hockney, R W, Gravitational Experiments with a Cylindrical Galaxy, Astrophys J, Vo1 150, No 3, p 797, 1967.
[7] Hohl, F, One- and Two-Dimensional Models to Study the Evolution of Stellar Systems, Symp. on Computer Simulation of Plasma and Many-Body Problems, Williamsburg, Va, April 19-21, NASA SP-153. 1967.
[8] Hohl, F, Computer Solutions of the Gravitational N-Body Problem, Proc of the International Astronomical Union Colloquium on the Gravitational Problem of N-Bodies, Paris, France, Aug. 1967.
[9] Chandrasekhar, S, Principles of Stellar Dynamics, Dover Publications, Inc 1942.
[10] Burbidge, E, M, Burbidge, G, R, Shelton, J W, A Peculiar Galaxy. That May be at a Young Evolutionary Stage, Astrophys J Vol 150, No 3, p 783, 1967.
[11] Sandage, A, The Hubble Atlas of Galaxies, Carnegie Institute of Washington Pub 618, Washington, DC, 1961.
[12] Lin, C C, Shu, F H, On the Spiral Structure of Disk Galaxies, Astrophys. J, vol 1k40, p 646, 1964.
