Received June 20, 1966
Computer experiments with a one-dimensional model have been performed to trace the evolution of stellar systems. The non-stationary equilibrium state for a one-dimensional system of stars computed for an interesting class of initial conditions is found to correspond to a minimum energy configuration. The results of the numerical experiments are compared with theory. For initial energies far from the minimum equilibrium energy the system becomes unstable and breaks up into smaller clusters. The one-dimensional model is of interest as an approximation to the distribution of velocity and mass normal to the galactic plane of a highly flattened galactic system.
This paper is concerned with the time behavior of a system of stars. Since the model to be investigated is one-dimensional, the stars can be represented by a large number of mass sheets. The motion of such sheets is sufficiently simple and the trajectories of several thousand stars can be followed on an electronic computer. We can thus perform numerical experiments to investigate the evolution of a stellar system and compare the results with the equilibrium conditions obtained from the self-consistent collision-less Boltzmann or Liouville equation.
As indicated by Michie (1964) there is a need to investigate the importance of orbital (or phase) mixing in the initial evolution of a system of stars. As the system of stars evolves, the gravitational field will change with time and the stars will follow complicated trajectories along which the individual stellar energies are not conserved. Henon (1966) has recently performed numerical experiments with a system of concentric spherical shells to study the relaxation of the mean gravitational field and the resulting approach to equilibrium for a spherically symmetric star cluster. The analogies between models where the material is stratified in parallel planes and where the material is distributed in concentric spherical shells have been examined by Woolley (1957).
According to Oort (1932) the velocities of stars normal to the galactic plane are decoupled from the other velocity components. Since the force on a star is approximately normal to the galactic plane the stars will oscillate perpendicular to the galactic plane with a period independent of the revolution of the galactic system. A one-dimensional model representing a system stratified in infinite parallel planes can therefore be used as an approximation to the distribution of velocity and density of stars normal to the galactic plane of a greatly flattened galactic system. Camm (1950) has considered steady solutions of the collisionless Boltzmann and the Poisson equations for such a model, and his results compared favorably with the observed densities.
The present paper is concerned primarily with investigating the conditions under which a system will approach a non-stationary equilibrium state and the conditions for which a system will become unstable and break up into smaller systems. The relaxation or characteristic time for a system of stars is given by
where G is the gravitational constant and ρ is the mass density. The characteristic length of interest is the Debye length D which is defined by
where VT τ is the velocity dispersion of the system of stars. The dimension of a system of stars near equilibrium is of the order of a Debye length.
The phrase non-stationary equilibrium refers to the equilibrium state which is established by the interaction of the stars with the smoothed potential of the system in a time of the order of 2π τc. Even though the rate of change of the system due to close encounters of stars is small, since changes do occur the system can not be regarded as steady. The binary encounters eventually cause the system to approach a state of statistical equilibrium. An actual stellar system can never completely attain a state of statistical equilibrium since this would involve a Gaussian velocity distribution. The absence of a potential barrier outside an actual stellar system permits the escape from the system of all stars with positive energies. Thus the leakage of stars will prevent the establishment of a Gaussian distribution necessary for a state of statistical equilibrium. Nevertheless, the time constant with which the three-dimensional system will continually tend toward a state of statistical equilibrium (even though it can never be completely attained) is of the order nD3τc, where n is the density of stars. This time is of the same order of magnitude as the time constant involved in Chandrasekhar's (1942) calculation of the dynamical friction. For the one-dimensional model the behavior is different because of the potential barrier outside the system. Fluctuations may now cause the establishment of a state of statistical equilibrium in a time which is at least of the order of (nD)2τc.
Lecar (1966a) has performed numerical integrations of a one-dimensional system of stars and has shown that no thermalization exists to order nDτc. For the case of a plasma this has been shown analytically by Eldridge and Feix (1963). Presently only times of order τc, during which collisional effects will be negligible, are of interest.
The model consists of a system of N stars which are represented by mass sheets. These sheets are of infinite extent along the y- and z-axes and the sheets are constrained to move along the x-axis. The equations of motion of all the N stars are solved simultaneously by computing the gravitational field at the position of each star and then integrating the equations of motion for each star over a small time interval Δt and repeating the process. When two sheets cross, they are allowed to pass freely through each other.
The evolution of the system of stars is studied by using various non-equilibrium initial distributions as input to the computer and by then following the time development of the system on the computer.
The equation of motion for the ith sheet with a mass mi per unit area is given by the expression
where xi is the position of the ith sheet and E(x,t) is obtained from
The position for the field is then
where
The IBM 7094 data-processing system was used to calculate the self-consistent motion of systems containing several thousand stars. The position and velocity of each sheet are computed at successive times t1, t2, t3, ... tr. For each time step Δt = tn+1−tn the new positions and velocities of the sheets are computed from the equations
and
where
and
Introduction of the term d3xi/dt2 allows us to take a larger time interval Δt in the computations while keeping the error very small. The accuracy of the computer program was checked by comparing results for various Δt's and also by checking the reversibility of the motion. For example, the fluctuations in the computed total energy were less than 0.005 per cent for the systems investigated. If n is the average density of stars and D the Debye length or dimension of the system, then the numerical experiments using the model described above are very nearly exact for nD ≫ 1 and actually give the time-dependent solution of the non-linear Liouville equation. The motion of the mass sheets can also be computed by a more exact method. That is, the times at which neighboring pairs of sheets cross are computed and by use of the shortest of these crossing times the new positions and velocities for the sheets affected are recomputed. The next pair of sheets that crosses is then found, and the process is repeated. The accelerations are constant until a pair of sheets crosses. This very accurate program has been used by Lecar (1966) to study certain invariants of the system. However, the exact program takes a large amount of machine time and is not suitable for investigating the motion of large numbers of sheets.
The potential energy per unit area for a system of N stars, each of mass m per unit area, is given by the equation
where Ej(x) is defined by equation (5) such that for each j, xj-1 <x<xj. The above definition of the potential energy is such that, if all stars are at the same point, then the potential energy of the system is zero. The kinetic energy of the system is given by the usual expression
In problems concerning the structure and evolution of stellar systems we can as a first approximation neglect the difference in mass and structure of individual stars. We therefore set all the star masses equal to m so that a distribution function f(x,v,t) will completely define the state of the system. When the effects of binary collisions between stars are negligible the distribution function satisfies the usual equation of stellar dynamics, that is, the Liouville equation (mass conservation)
where E=d2x/dt2 = −∂φ/∂x is given by the Poisson equation
The limitations of the Liouville equation in describing stellar systems are discussed by Kurth (1957).
To study the time development of a system of stars one could use the distribution function at t=0, that is, f(x,v,t = 0), as the initial condition and then solve the timedependent non-linear Liouville equation. Since in general this problem cannot be solved analytically a high-speed computer is used to determine the evolution of the system, and equations (12) and (13) are used only to obtain the steady-state solution for the system.
The virial theorem can be applied to the one-dimensional system to obtain a relation between the potential and kinetic energy of the system in equilibrium. Consider the identity
if this expression is multiplied by the distribution function and integrated over phase space, the following result is obtained:
Since the system is in equilibrium, the left-hand side of equation (15) is equal to zero. Also
where T is the total kinetic energy of the system. Using ρ(x) = ∫mfdv and d2x/dt2=-dφ/dx in equation (15), one obtains the expression
The Poisson equation can be written as
Therefore, using an integration by parts, one obtains the result
The last equation is simply the potential energy P of the system. Equation (15) then becomes
irrespective of the distribution function. In most of our numerical studies we found that the system very quickly approached conditions such that equation (20) was satisfied.
In the steady state ∂f/∂t=0, and the Liouville equation tells us that the system is described by a distribution function which is a function only of U=½mv2+mφ, the total energy of a star. That is, if a time-independent equilibrium state exists, then
where F(U) may be any function of U. Of course we are interested only in functions F(U) that are stable. In general, the form of F(U) depends on the initial distribution and must be obtained by following the time development of the non-linear Liouville equation. However, there is one type of initial distribution where F(U) is known without actually solving the Liouville equation. For this distribution the initial f is taken to be constant over a certain region of phase space and is zero outside this region. According to the Liouville theorem f remains constant along the different trajectories so that the region can only change its shape with time while keeping its area constant. For this reason the distribution function just described has been called the water-bag model by DePackh (1962). The water-bag model is of interest to us primarily because it allows us to calculate exactly the equilibrium configuration of the one-dimensional star gas for comparison with the computer results.
Two initial distributions considered in this paper are the water-bag model and a distribution which has a constant density over a region of the x-axis and has a Maxwellian velocity distribution. For the water-bag model we know F(U); it is constant for 0 ≤U≤mε and is zero for U outside of this region; mε is some maximum energy to be determined later. For the second distribution we assume that F(U) is Maxwellian.
The equilibrium solution for the water-bag model will now be obtained. The initial shape of the water bag is taken to be a rectangle defined by the area in phase space between ±x0 and ±v0. Using f(x,v,t→∞) = F(U) in the Poisson equation (13) and changing to an integration over dU, one obtains
where mε is the energy such that F(U)=0 for U>mε and F(U)=A for 0≤U≤mε. Since the area of the system in phase space remains constant, the value of A is obtained from the relation
or
A first integration of equation (22) gives the result
where we have chosen φ=dφ/dx=0 at x=0. A second integration gives the final result
Figure 1 shows the resulting Φ and E where the normalizations ψ = (2πGmN/ε)x and Φ=φ/ε have been used.
Let xs be the coordinate defining the boundary of the one-dimensional system, then φ(xs)=ε and (dφ/dx)x=xs)= 2πGNm. 5">equation (25), one obtains the result
The value of x is obtained from equation (26) by taking ε as the upper limit of the integral, thus
The value of the maximum velocity is given by the expression
Equation (22) shows that the number density is given by
where we have made use of equation (24).
The steady-state solution for the water-bag distribution was obtained by using only conservation of area in the two-dimensional phase space. However, if the equilibrium state is to be reached, then the initial and final energy of the system must also be equal. For the initial rectangular water bag the total energy W is given by
As before ±v0 and ±x0 define the rectangular area of the initial distribution. The equilibrium state as defined by equation (26) has a total energy given by the expression
where
It is easily shown that for any given x0 and v0 the value of equation (32) is always less than that of equation (31). It is thus likely that the equilibrium distribution for the water-bag model has a minimum energy property. To check this property the analysis of DePackh (1962) is first applied to the water-bag model. Let the function V±(x,t) describe the contour of the water bag in phase space. The subscripts ± refer to the upper and lower branches of the contour. The Liouville equation can then be written as
and δ(ζ)is the Dirac delta function. Equation (33) can be written in the following form
After integrating equation (34) over v one obtains the result
or
The equilibrium contour is obtained from the equation
In order to apply the variational principle a generation function g [ x,v+(x), θ(x) ] with
is introduced such that the total energy of the system is given by
For the sake of simplicity the contour of the water bag has been taken to be symmetric such that v+=−v−. Following Courant and Hilbert (1953) the extremum of equation (38) is found by the usual methods of variational calculus. The end points ±xs of the contour are held fixed, and the function v+(x) receives a variation aη(x), where a is an arbitrary constant and η(x) is an arbitrary function which vanishes at the end points. One then finds that
Since n is an arbitrary function, the condition for obtaining an extremum is
Using the first integral in equation (19) for the potential energy one can easily show that for the present problem g has the form
where 4πGm(N/2 − 2A θ) = E and A is the magnitude of f. If the expression for g given by equation (41) is used in equation (40), one obtains the result
which is identical to equation (37). Thus the equilibrium solution given by equation (26) is a minimum energy configuration. The consequence of the minimum energy property is that, starting from any non-equilibrium state, energy conservation will prevent the system from reaching the steady state described by equations (26) and (30). This result was to be expected from the analysis of DePackh (1962), which shows that for a plasma small perturbations of the equilibrium state are not damped. Nevertheless, the interesting point is that our numerical results show that the system does its best within the limitation of energy conservation to approach the equilibrium state. In general, we can hope to approach the equilibrium state very closely whenever the initial energy is not too different from the energy of the equilibrium state given by equation (32).
Two types of initial water-bag distributions have to be separately considered. The first type is a fat water bag where the total energy is near the energy of the equilibrium distribution. For this case the system quickly approaches the equilibrium state with the exception of the arms in phase space which develop to accommodate the excess energy. The second type of initial distribution is a thin water bag which has a total energy much larger than the energy of the equilibrium distribution. For this second case the system cannot even come near its equilibrium state, and an instability develops which causes the system to break up into smaller clusters.
For a system with an initially constant density over a region and a Gaussian velocity distribution it is assumed the final state is described by an F(U) of the form A exp (-kU). Following the same method used in obtaining the solution for the water-bag model one obtains the results
The one-dimensional sheet model includes individual or star-star interactions as well as the interaction of the stars with the smoothed potential of the system since the computer simply solves the equations of motion of the stars in the system. The individual behavior is due to the graininess and will be affected by going to the fluid limit as implied by the Liouville equation.
The effects of graininess on the time development of a system can be determined by varying the number of sheets while keeping the total mass of the system constant. The effects of graininess must be checked to determine whether the Liouville equation adequately describes the system. Figure 2 shows the time development of two equivalent systems with different graininess. The curves indicate the variation in time of the kinetic energy for a water-bag model. The initial conditions are obtained by using a random-number generator to give a nearly uniform distribution over a rectangular region in phase space. For all the numerical calculations the gravitational constant was normalized such that 4πG = 1. Also all time scales are normalized to τco= τc (t = 0). The initial ratio of kinetic to potential energy for both curves in Figure 2 is 2 whereas the equilibrium value is 0.5. The upper curve is for a system of 1000 stars with nD ∼ 1000 and m=2. The lower curve is for a system of 2000 stars with nD > 2000 and m = 1. The oscillations of the kinetic energy are identical for the two systems, indicating that the graininess effects, that the effects of binary collisions between stars, are negligible for any time scales of present interest. We also note that after only a few periods non-linear damping or phase mixing has almost eliminated the oscillations in the kinetic energy and the kinetic energy remains near the equilibrium value indicated by the dashed line. We can thus be assured that the Liouville equation correctly describes the system.
To determine when the effects of graininess become important, the evolution of several equivalent systems with varying numbers of stars is compared. Figure 3 shows the variation of the kinetic energy for equivalent systems which contain from 20 to 500 stars. For the system with 500 stars the behavior is still very similar to that shown in Figure 2. The time behavior of the energy distribution function and the density for equivalent systems which contain more than 500 stars is also very similar. As can be seen from Figure 3, when the number of stars of the system becomes 200 or less the time behavior of the kinetic energy differs markedly from that of Figure 2. Also, the time development of the energy distribution function and of the density is different for the systems shown in Figure 3. In general, it was found that oscillations in the various parameters of the system persisted much longer when the number of stars of the system was less than 500.
Figure 4 shows the time development of the density of the water-bag model for a system of 2000 stars with a mass m per unit area equal to 1. The ratio of initial to equilibrium energy (as given by eq. [32]) for this system is 1.33, and the initial ratio of kinetic to potential energy is 2. The corresponding variation of the kinetic energy is shown by the lower curve of Figure 2. The dashed curve in Figure 4 is the theoretical density as obtained from equation (30). After only a few periods 2πτc, the experimental density is close to the calculated value. The time development of the energy distribution function for the same system is shown in Figure 5. The dashed curve represents the theoretical distribution. Again, after only about six periods, 2πτc, the experimental distribution closely approximates the theoretical distribution. Thus, even for values of the ratio of initial to equilibrium energy of 1.33, the system approaches its non-stationary equilibrium state closely, in a time of the order of 2πτc.
The best way to see the approach to equilibrium is to take a sequence of pictures showing the time development of the system in phase space. Figure 6 shows such a sequence for the system of 2000 stars. Each star is represented by a small circle. While approaching the non-stationary equilibrium state indicated by the area inside the oval in Figure 6, d, the system rotates in phase space with a frequency near 1/τc. The relatively few particles in the long arms are required to accommodate the excess energy, and they correspond to the high-energy peak in Figure 5. The arms develop because the period of oscillation of the sheets increases with energy. Figure 7 shows the period as a function of energy for the equilibrium state.
The results of Henon (1966) show that the oscillations for a system of spherically concentric mass shells are damped more rapidly than for the one-dimensional system shown in Figure 2. The reason for this is that the period of oscillation of a spherical shell increases much more rapidly with increasing energy than is the case for a sheet in the one-dimensional system. Therefore, the effectiveness of phase mixing, which depends on the change in
In the equilibrium state the gravitational potential is time-independent and the trajectories of individual stars form fixed closed contours in phase space as shown in Figure 8. Figure 9 shows the actual trajectory of a mass sheet for the system shown in Figure 6. The phase-space trajectory shown in Figure 9 is typical for systems which are initially near the equilibrium state. Initially, while the system is approaching the equilibrium configuration the orbits are very much perturbed. This causes the redistribution of the orbits required to approach equilibrium. After only about two orbits the trajectories are nearly the same as the theoretical orbits.
The time behavior described so far is characteristic of the water-bag model whenever the initial energy given by equation (31) is near the equilibrium energy given by equation (32). Since the system approaches the equilibrium configuration within the limits of energy conservation, it appears that the steady state given by equation (26) does have a physical meaning. The stability of the steady state has been checked using an initial distribution which satisfies equation (26). Even after many periods 22πτc the system is found to remain in the equilibrium state. The effect of a perturbation of the equilibrium state has also been checked. Figure 10 shows the time development of a sinusoidal perturbation of the equilibrium state. The perturbation results in non-linear effects, and because the outermost stars are now outside the main system they have a longer period and cause the development of arms. After this small perturbation the non-stationary equilibrium state can never be completely reached.
>
Figure 11 shows the time development for a water-bag distribution of a system of 1000 stars with a ratio of initial to equilibrium energy equal to 6. The initial ratio of kinetic to potential energy for this system is 32. The system becomes unstable and quickly breaks up into four smaller clusters. Each of these clusters individually tends to approach an equilibrium configuration. In Figure 12 the position of every 20th star for this system is plotted as a function of time. After only one complete oscillation the bunching of the stars into four distinct regions is clearly visible. Calculations for systems where the initial distribution of stars in phase space is obtained by using different random-number sequences were made which indicated that the number of clusters into which a system breaks up depends on the inhomogeneities in the initial distribution of stars in phase space.
A system with an initially constant density between ±x0 and with a Maxwellian velocity distribution has a time development similar to that for the water-bag distribution. For example, if the initial ratio of kinetic to potential energy is of the order of the equilibrium value, then the system shows a time behavior very similar to that shown in Figure 6. However, the arms in phase space develop more rapidly because some of the high-velocity mass sheets are now able to stay outside the main system longer and therefore have a longer period. Also the density and energy distribution function tend to approach conditions satisfied by equations (43)-(45). If the initial ratio of kinetic to potential energy is far from the equilibrium value, the systems with a Maxwellian velocity distribution will break up into smaller clusters. However, the breakup now does not occur as readily as for the water-bag distribution.
In comparing his theoretical results with the force E(z) normal to the galactic plane, Camm (1950) used the results of Oort (1932), who has since revised them. The revised values of E(z) have been tabulated by Lindblad (1959) and are compared in Figure 13 with the two distributions studied in this paper. The variation of E(x) for the Gaussian distribution is in good agreement with the experimentally obtained variation. The difference for small values of x is probably due to the central core of the Galaxy. As expected, for the water-bag distribution the force E(x) differs markedly from the experimentally obtained variation for large x.
For an initial distribution of the water-bag type it was shown that the equilibrium state for a given area in phase space is a minimum energy configuration. Starting from any non-equilibrium state, the stationary state described by the time-independent Liouville equation can never be completely reached irrespective of any graininess effect in the system. This is because the arms which develop to accommodate the excess energy wind around the main system and become longer and thinner as time progresses. The beginning of this process is illustrated in Figure 6. However, after a sufficiently long time it will become progressively more difficult to detect changes in f(x,v,t) since successive turns of the arms approach each other. It may then be possible to obtain an averaged distribution function by averaging over φ. Such an analysis has been performed by Lynden-Bell (1966). Lynden-Bell suggests that the encounterless relaxation of an unsteady system will lead to an equilibrium related to Fermi-Dirac statistics. However, in our investigation of encounterless relaxation with a one-dimensional model we have not been able to confirm Lynden-Bell's theory.
Our computer calculations show that for initial distributions other than the water-bag type the results are nevertheless very similar to those for the water-bag distribution. It should therefore be investigated whether the minimum energy principle is more general than investigated in this paper.
It is a pleasure to thank Barbara Weigel for her help in carrying out the computer calculation and Leo Staton for suggesting the use of the variational method in determining the minimum energy property.
Camm, G. L. 1950, M.N., 110, 305.
Chandrasekhar, S. 1942, Principles of Stellar Dynamics (Chicago: University of Chicago Press).
Courant, R., and Hilbert, D. 1953, Methods of Mathematical Physics (New York: Interscience Publishers).
DePackh, D. C. 1962, Journal of Electronics and Control, 13, 417.
Eldridge, 0. C., and Feix, M. R. 1963, Phys. Fluids, 6, 398.
Henon, M. 1966, paper presented at the 14th International Astrophysical Symposium, University of Liege, Belgium, June, 1966.
Kurth, R. 1957, Introduction to the Mechanics of Stellar Systems (New York: Pergamon Press).
Lecar, M. 1966a, paper presented at the 14th International Astrophysical Symposium, University of Liege, Belgium, June, 1966.
Lecar, M. 1966b, private communication.
Lindblad, B. 1959, Encyclopedia of Physics, Vol. 53 (Berlin: Springer-Verlag).
Lynden-Bell, D. 1966, paper presented at the 14th International Astrophysical Symposium, University of Liege, Belgium, June, 1966.
Michie, R. W. 1964, "The Dynamics of Star Clusters," in Annual Review of Astronomy and Astrophysics, vol. 2, p.49, 1964
Oort, J. H. 1932, Bulletin of the Astronomical Institutes of the Netherlands, Vol. 6, p.249,
Woolley, R. v. d. R. 1957, Mon.Not.R., astr Soc 117, 198,