Abstract
This paper is devoted to different modifications of two standard softenings of the gravitational attraction (namely, the Plummer and Hernquist softenings), which are commonly used in cosmological simulations based on the particle-particle (PP) method, and their comparison. It is demonstrated that some of the proposed alternatives lead to almost the same accuracy as in the case of the pure Newtonian interaction, even despite the fact that the force resolution is allowed to equal half the minimum interparticle distance. The revealed way of precision improvement gives an opportunity to succeed in solving Gurzadyan’s Problem 5 and bring modern computer codes up to a higher standard.
1. Introduction
Among Gurzadyan’s “10 key problems in stellar dynamics” [1], any successful step towards solving Problem 5 will exert powerful influence on the state of affairs in stellar and galactic dynamics. The statement of this problem consists in creation of a computer code describing the -body system with the phase trajectory being close to that of the real physical system for long enough time scales. In view of primary importance of cosmological simulations in analyzing the large scale structure formation and comparing results with predictions of different theories of the Universe evolution, the goal of Problem 5 appears particularly important.
In this paper one such step is proposed. Obviously, the higher precision can be achieved in -body computer simulations, the more rigorous constraints can be imposed on parameters of a concrete cosmological model. The well known PP method computing forces according to the Newton law of gravitation represents an accurate -body technique (see, e.g., [2, 3]). At the same time, it suffers from an evident shortcoming. Since the Newtonian gravitational potential is singular at the particles’ positions, softening is required in order to avoid divergences of forces when the interparticle separation distances are very small. Introduction of softening ensures collisionless behavior of the system and simplifies numerical integration of its equations of motion essentially. However, a high price to pay is noticeable reduction of precision. Below an attempt is made to modify two generally accepted softenings in such a way that the accuracy of computer simulations becomes improved dramatically without any unjustified complication of the equations of motion or the integration technique.
The paper is organized as follows. In Section 2 the equations of motion are written down and two standard softenings, namely, the Plummer and Hernquist ones, are introduced. In Section 3 various modifications are proposed and their efficiency is compared with respect to the same illustrative example. Main results are discussed briefly in Section 4.
2. Plummer and Hernquist Softenings
According to the mechanical approach/nonrelativistic discrete cosmology, developed recently in a series of papers [4–6] in the framework of the conventional CDM ( Cold Dark Matter) model, the scalar cosmological perturbations in the Universe with flat spatial topology can be described by the following perturbed FLRW metric in both nonrelativistic matter dominated and dark energy dominated eras: where is the scale factor depending on the conformal time ,
Here is the comoving radius-vector, stands for the Laplace operator, is the gravitational constant, and represents the rest mass density in the comoving coordinates , and . This quantity is time-independent within the adopted accuracy (both the nonrelativistic and weak field limits are applied), and denotes its constant average value. Really, changes with the lapse of time in view of peculiar motion of cosmic bodies; however, the corresponding velocities are small enough at the considered nonrelativistic matter dominated and dark energy dominated stages of the Universe evolution. Consequently, this temporal change of may be disregarded when determining the gravitational potential from (2). In other words, it is determined by the positions of cosmic bodies but not by their peculiar velocities, as it certainly should be in the nonrelativistic and weak field limits. Naturally, the introduced function (2) satisfies the linearized Einstein equations for the metric (1) within the adopted accuracy (see [4, 6]).
It is worth noting that in the considered flat spatial topology case the scale factor may be dimensionless; then the comoving coordinates have a dimension of length, and vice versa. In order to be specific, let us choose the first option; then the rest mass density is measured in its standard units (namely, masslength3).
In complete agreement with [7–9], the following equations of motion, which describe dynamics of the -body system experiencing both the gravitational attraction between its constituents and the global cosmological expansion of the Universe, can be immediately derived from (1) and (2) (the cutoff of the gravitational potential in the special relativity spirit introduced in [10] in order to resolve the problem of nonzero average values of first-order scalar perturbations is not considered here being irrelevant for the given investigation. Really, this cutoff is supposed to take place at great distances of the order of the particle horizon and certainly does not affect dynamics on smaller scales described by the written down standard equations of motion): where stands for the physical radius-vector of the th particle and represents its mass. These equations are also commonly used for simulations at astrophysical (i.e., noncosmological) scales when the second term in the left hand side of (3) is irrelevant and may be neglected.
Apparently, the right hand side of (3) is a superposition of forces which originate from Newtonian gravitational potentials of single point-like particles. If such a particle possessing the mass is located at the point , then its rest mass density in the physical coordinates , and the potential of the produced gravitational field is singular at the location point. In order to suppress this singularity for reasons enumerated briefly in Section 1 (namely, avoiding divergences of interparticle forces, ensuring collisionless dynamics, and simplifying numerical integration of equations of motion), let us consider two different models dealing with non-point-like gravitating masses. The density and the potential of a single body in the Plummer model [11–15]: The same quantities in the Hernquist model [13–16]: where is the softening length/parameter (called the force resolution as well) typically amounting to a small percent of the mean interparticle distance. While still dealing with point-like particles, one usually takes into account the gravitational interaction by means of or in modern computer simulations based on the PP method. Thus, the force resolution should not be attributed to any real physical sense representing just a mathematical trick eliminating singularity. Both Plummer and Hernquist potentials converge to the Newtonian one when (in this limit, as one can also easily verify, both densities in (4) and (5) tend to the same expression corresponding to a point-like massive particle, as expected). However, for any nonzero value of the attraction between every two bodies is changed with respect to the Newton law of gravitation at all separations. In particular, both analyzed potentials (4) and (5) tend to zero as when (), but in the opposite limit () they behave as a constant , so Newtonian singularity is absent. The next section is entirely devoted to controllable elimination of this defect of the above-mentioned softenings.
3. Illustration of Accuracy Improvement
For illustration purposes let us consider a hyperbolic trajectory of a test particle with the mass in the gravitational field of the mass resting in the origin of coordinates. This trajectory is given by the following functions [17] (, , and denote Cartesian coordinates on the plane of motion and time, resp.): where stands for the eccentricity, represents the varying parameter, and is the so-called semiaxis of a hyperbola, being interrelated with the distance to perihelion : . In what follows, the values and are used.
The functions (6) satisfy the equations of motion:
Introducing the normalized quantities one can rewrite (7) in the form being more convenient for the subsequent numerical integration:
According to (8), if , then , , and ; besides, . The enumerated values will serve as initial conditions hereinafter.
The exact solution (8) is depicted on Figure 1 (the black curve) together with the numerical solution of (9) (red points). The leapfrog “drift-kick-drift” numerical integration scheme is applied here and below with the fixed time step .
Orange points correspond to the modified equations of motion: that is, to the “Newton-Hernquist” conversion in the expression for the gravitational potential. The normalized softening length everywhere equals (that is, amounts to of : , meaning quite close approaching). Obviously, the orange points lie rather far from the red ones. Consequently, the error is significant.
Further, one obtains purple points modifying the Hernquist potential:
The idea underlying this modification is simple: at each iteration one can make calculations using both and softenings and then interpolate results to the zero softening parameter. In other words, the expression in the right hand side of (12) is constructed purposely in such a way that for its derivative with respect to , being proportional to the gravitation force, behaves as up to the second order of smallness concerning the ratio . The term of the first order is missing; therefore, the actual superposition of two Hernquist potentials with different softenings reduces the simulation error in comparison with the previous case. Really, the purple points are noticeably closer to the red ones than the orange points. However, precision is still low. Of course, one can increase a number of terms in the superposition and apply a higher order interpolation, but introduction of each additional term requires more computational time and, consequently, is not reasonable.
Green points correspond to the modified equations of motion: that is, to the “Newton-Plummer” conversion in the expression for the gravitational potential. Precision is higher than in the previous case because the derivative of the expression in the right hand side of (14) with respect to for behaves as , so the deviation from the pure Newtonian behavior is now four times smaller.
Finally, one gets blue points modifying the Plummer potential [18]:
Now for the deviation from the pure Newtonian behavior represents a quantity of the fourth order of smallness concerning the ratio . Consequently, precision is really high even despite the fact that the condition holds true as before.
4. Conclusion
In this paper a promising opportunity of increasing the accuracy of computer -body simulations based on the PP method is addressed. Namely, the inevitable error arising from gravitational softening is reduced considerably by modifying the commonly used Plummer ~1 and Hernquist ~1 potentials. In particular, the proposed ~1 potential with gives better approximation since for the corresponding gravitation force differs from the standard Newtonian one in a small quantity ~(. This is demonstrated explicitly for with the help of the concrete illustrative example of one particle moving along the hyperbolic trajectory in the softened gravitational field of another one. The force resolution is taken amounting to half the minimum separation distance, but despite this fact the suggested alternative softening is characterized by much higher precision being much closer to the pure Newtonian picture than the standard ones.
Apparently, while improving numerical integration in the region (where owing to this important inequality the expansion into series with respect to the ratio is allowed), the developed scheme still misrepresents the picture for (where the above-mentioned expansion is forbidden). However, if such close approachings seldom happen, this misrepresentation is not significant for the whole -body system behavior description. Thus, this scheme can really play an important role in astrophysical/cosmological modeling and solving Problem 5. In other words, the proposed modifications reducing simulation errors caused by softening can help to bring the phase trajectory of the -body system in a corresponding computer code much closer to that of the real physical one.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by NSF CREST award HRD-1345219 and NASA Grant NNX09AV07A. The author would like to thank S. J. Aarseth for valuable comments and the referee for critical remarks which have considerably improved the presentation of the obtained results.