Gurzadyan's Problem 5 and improvement of softenings for cosmological simulations using the PP method

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.


Introduction
Among Gurzadyan's "10 key problems in stellar dynamics" [1], any successful step towards solving the Problem 5 will exert powerful influence on the state of affairs in stellar and galactic dynamics. The statement consists in creation of a computer code describing the N -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 the Problem 5 appears particularly important.
In this Letter one such a step is proposed. Obviously, the higher precision can be achieved in N -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 N -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 Letter is organized as follows. In the next section the equations of motion are written down and two standard softenings, namely, the Plummer and Hernquist ones, are introduced. In the subsequent section various modifications are proposed and their efficiency is compared with respect to the same illustrative example. Main results are discussed briefly in Summary.

Plummer and Hernquist softenings
According to the mechanical approach/nonrelativistic discrete cosmology, developed recently in a series of papers [4,5,6] in the framework of the conventional ΛCDM model, the scalar cosmological perturbations in the late Universe with flat spatial topology can be described by the following perturbed FLRW metric: where a(η) is the scale factor, Here △ = δ αβ ∂ 2 /(∂x α ∂x β ) stands for the Laplace operator, G is the gravitational constant, and ρ represents the rest mass density in the comoving coordinates. This quantity is time-independent within the adopted accuracy (both the nonrelativistic and weak field limits are applied), and ρ denotes its constant average value.
In complete agreement with [7,8,9], the following equations of motion, which describe dynamics of the N -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) 1 : where R i = ar i stands for the physical radius-vector of the i-th particle, and m i represents its mass. These equations are also commonly used for simulations at astrophysical (i.e. non-cosmological) 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 m is located at the point R = 0, then its rest mass density in the physical coordinates ρ ph = mδ(R), and the potential of the produced gravitational field φ = −Gm/R is singular at the location point. In order to suppress this singularity for reasons enumerated briefly in Introduction, 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,12,13,14,15] while the same quantities in the Hernquist model [13,14,15,16] ρ where ε is the softening length/parameter (called the force resolution as well) typically amounting to a few percent of the mean interparticle distance. While still dealing with point-like particles, one usually takes into account the gravitational interaction by means of φ (P) or φ (H) in modern computer simulations based on the PP method. Thus, the force resolution ε should not be attributed any real physical sense representing just a mathematical trick eliminating singularity. Both Plummer φ (P) and Hernquist φ (H) potentials converge to the Newtonian one when ε → 0. 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. The next section is entirely devoted to controllable elimination of this defect of the above-mentioned softenings.

Illustration of accuracy improvement
For illustration purposes let us consider a hyperbolic trajectory of a test particle with the mass m in the gravitational field of the mass M resting in the origin of coordinates. This trajectory is given by the following functions [17]: where ǫ > 1 stands for the eccentricity, ξ represents the varying parameter, and A is interrelated with the distance to perihelion R min : A(ǫ − 1) = R min . In what follows, the values ǫ = 1.1 and ξ ∈ [0, 0.15] 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 ξ = 0, thent = 0,X = ǫ − 1,Ỹ = 0, besides, dX/dt = 0, The enumerated values will serve as initial conditions for all simulations being compared with each other hereinafter.
The exact solution (8) is depicted on Fig. 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 ∆t = 0.0025. Orange points correspond to the modified equations of motion i.e. to the "Newton -Hernquist" conversion in terms of gravitational potentials. The normalized softening lengthε everywhere equals 0.05 (that isε amounts to 50 % ofR min :ε =R min /2, 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 2ε softenings and then interpolating 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 forR ≫ε its derivative with respect toR, being proportional to the gravitation force, behaves as −1/R 2 +6/R 2 ·(ε/R) 2 up to the second order of smallness concerning the ratiõ ε/R. 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, then 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 i.e. to the "Newton -Plummer" conversion in terms of gravitational potentials. Precision is higher than in the previous case because the derivative of the expression in the right hand side of (14) with respect toR forR ≫ε behaves as −1/R 2 + 1.5/R 2 · (ε/R) 2 , so the deviation from the pure Newtonian behavior −1/R 2 is now four times smaller. Finally, one gets blue points modifying the Plummer potential: Now forR ≫ε the deviation from the pure Newtonian behavior represents a quantity of the forth order of smallness concerning the ratioε/R. Consequently, precision is really high even despite the fact that the conditionε =R min /2 holds true as before.

Summary
In this Letter a promising opportunity of increasing the accuracy of computer N -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/ √ R 2 + ε 2 and Hernquist ∼ 1/(R + ε) potentials. In particular, the proposed ∼ 1/ (R n + ε n ) 1/n potential with n > 2 gives better approximation. This is demonstrated explicitly for n = 4 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 R > ε, the developed scheme still misrepresents the picture for R ε. However, if such close approachings happen seldom, this misrepresentation is not significant for the whole N -body system behavior description. Thus, this scheme can really play an important role in astrophysical/cosmological modeling and solving the above-mentioned Problem 5.