Gurzadyan's problem 5 and improvement of softenings for cosmological simulations using the PP method.
Among Gurzadyan's "10 key problems in stellar dynamics" , 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 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 Problem 5 appears particularly important.
In this paper one such 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 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 [LAMBDA]CDM ([LAMBDA] 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:
d[s.sup.2] [approximately equal to] [a.sup.2] [(1 + 2[PHI])d[[eta].sup.2] - (1 - 2[PHI])[[delta].sub.[alpha][beta]] d[x.sup.[alpha]] d[x.sup.[beta]]], [alpha], [beta] = 1, 2, 3, (1)
where [alpha]([eta]) is the scale factor depending on the conformal time [eta],
[PHI] = ([eta], r) = [phi](r)/[c.sup.2]a([eta]), [DELTA][phi] = 4[pi]G ([rho] - [bar.[rho])
Here r is the comoving radius-vector, [DELTA] = [[delta].sup.[alpha][beta]][[partial derivative].sup.2]/([partial derivative][x.sup.[alpha]][partial derivative][x.sup.[beta]]) stands for the Laplace operator, G is the gravitational constant, and p represents the rest mass density in the comoving coordinates [x.sup.1] [equivalent] x, [x.sup.2] [equivalent] y, and [x.sup.3] [equivalent] z. This quantity is time-independent within the adopted accuracy (both the nonrelativistic and weak field limits are applied), and [bar.[rho]] denotes its constant average value. Really, p 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 [rho] maybe 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 [PHI](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 a may be dimensionless; then the comoving coordinates [x.sup.[alpha]] have a dimension of length, and vice versa. In order to be specific, let us choose the first option; then the rest mass density p is measured in its standard units (namely, mass/[length.sup.3]).
In complete agreement with [7-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) (the cutoff of the gravitational potential in the special relativity spirit introduced in  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):
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)
where [R.sub.i] = a[r.sub.i] stands for the physical radius-vector of the ith particle and [m.sub.i] 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 m is located at the point R = 0, then its rest mass density in the physical coordinates [[rho].sub.ph] = m[delta](R), and the potential of the produced gravitational field [phi] = -Gm/R 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]:
[[rho].sup.(P).sub.ph] = [3m/4[pi][[epsilon].sup.3]] [(1 + [R.sup.2]/[[epsilon].sup.2]).sup.-5/2], [[phi].sup.(P)] = - Gm/[square root of [R.sup.2] + [[epsilon].sup.2]]. (4)
The same quantities in the Hernquist model [13-16]:
[[rho].sup.(H).sub.ph] = [m[epsilon]/2[pi]R] [1/[(R + [epsilon]).sup.3]], [[rho].sup.(H)] = [Gm/[R + [epsilon]], (5)
where [epsilon] 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 [[phi].sup.(P)] or [[phi].sup.(H)] in modern computer simulations based on the PP method. Thus, the force resolution e should not be attributed to any real physical sense representing just a mathematical trick eliminating singularity. Both Plummer [[phi].sup.(P)] and Hernquist [[phi].sup.(H)] potentials converge to the Newtonian one when [epsilon] [right arrow] 0 (in this limit, as one can also easily verify, both densities in (4) and (5) tend to the same expression m[delta](R) corresponding to a point-like massive particle, as expected). However, for any nonzero value of e 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 -Gm/R when R [right arrow] +[infinity] (R [much greater than] [epsilon]), but in the opposite limit R [right arrow] 0 (R [much less than] [epsilon]) they behave as a constant -Gm/[epsilon], 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 m in the gravitational field of the mass M resting in the origin of coordinates. This trajectory is given by the following functions  (X, Y, and t denote Cartesian coordinates on the plane of motion and time, resp.):
X = A([epsilon] - cosh [xi]), Y = A [square root of ([[epsilon].sup.2] - 1)] sinh [xi], t = [square root of [[A.sup.3]/Gm] ([epsilon] sinh [xi] - [xi]), (6)
where [epsilon] > 1 stands for the eccentricity, [xi] represents the varying parameter, and A is the so-called semiaxis of a hyperbola, being interrelated with the distance to perihelion [R.sub.min]: A([epsilon] - 1) = [R.sub.min]. In what follows, the values [epsilon] = 1.1 and [xi] [member of] [0, 0.15] are used.
The functions (6) satisfy the equations of motion:
[[d.sup.2]X/d[t.sup.2]] = -GM [X/[R.sup.3]], [[d.sup.2]Y/d[t.sup.2]] = -GM [Y/[R.sup.3]], R = [square root of ([X.sup.2] + [Y.sup.2])].
Introducing the normalized quantities
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (8)
one can rewrite (7) in the form being more convenient for the subsequent numerical integration:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (9)
According to (8), if [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. 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 [DELTA][??] = 0.0025.
Orange points correspond to the modified equations of motion:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (10)
that is, to the "Newton-Hernquist" conversion
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)
in the expression for the gravitational potential. The normalized softening length [??] everywhere equals 0.05 (that is, e amounts to 50% of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], 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:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (12)
The idea underlying this modification is simple: at each iteration one can make calculations using both [??] and 2[??] 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 [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] its derivative with respect to [??], being proportional to the gravitation force, behaves as [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] up to the second order of smallness concerning the ratio e/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 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:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (13)
that is, to the "Newton-Plummer" conversion
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)
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 [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] behaves as [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], so the deviation from the pure Newtonian behavior -1/[[??].sup.2] is now four times smaller.
Finally, one gets blue points modifying the Plummer potential :
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (15)
Now for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] the deviation from the pure Newtonian behavior represents a quantity of the fourth order of smallness concerning the ratio [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII. Consequently, precision is really high even despite the fact that the condition [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] holds true as before.
In this paper 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/[square root of [R.sup.2] + [[epsilon].sup.2]] and Hernquist ~1/(R + [epsilon]) potentials. In particular, the proposed ~1/[([R.sup.n] + [[epsilon].sup.n]).sup.1/n] potential with n > 2 gives better approximation since for R > [epsilon] the corresponding gravitation force differs from the standard Newtonian one in a small quantity ~[([epsilon]/R).sup.n]. 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 e 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 > [epsilon] (where owing to this important inequality the expansion into series with respect to the ratio [epsilon]/R < 1 is allowed), the developed scheme still misrepresents the picture for R [less than or equal to] [epsilon] (where the above-mentioned expansion is forbidden). However, if such close approachings seldom happen, 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 Problem 5. In other words, the proposed modifications reducing simulation errors caused by softening can help to bring the phase trajectory of the N-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.
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.
 V. G. Gurzadyan, "10 key problems in stellar dynamics: in retrospect," http://arxiv.org/abs/14070398.
 S. J. Aarseth, Gravitational N-Body Simulations, Cambridge Monographs on Mathematical Physics, Cambridge University Press, Cambridge, UK, 2003.
 J. Makino, T. Fukushige, M. Koga, and K. Namura, "GRAPE-6: massively-parallel special-purpose computer for astrophysical particle simulations," Publications of the Astronomical Society of Japan, vol. 55, no. 6, pp. 1163-1187, 2003.
 M. Eingorn and A. Zhuk, "Hubble flows and gravitational potentials in observable Universe," Journal of Cosmology and Astroparticle Physics, vol. 9, article 026, 2012.
 M. Eingorn, A. Kudinova, and A. Zhuk, "Dynamics of astrophysical objects against the cosmological background," Journal of Cosmology and Astroparticle Physics, vol. 4, article 010, 2013.
 M. Eingorn and A. Zhuk, "Remarks on mechanical approach to observable Universe," Journal of Cosmology and Astroparticle Physics, vol. 5, article 024, 2014.
 P. J. Peebles, The Large-Scale Structure of the Universe, Princeton University Press, Princeton, NJ, USA, 1980.
 V Springel, "The cosmological simulation code GADGET-2," Monthly Notices of the Royal Astronomical Society, vol. 364, no. 4, pp. 1105-1134, 2005.
 L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields, vol. 2 of Course of Theoretical Physics Series, Oxford Pergamon Press, Oxford, UK, 4th edition, 2000.
 M. Eingorn, M. Brilenkov, and B. Vlahovic, "Zero average values of cosmological perturbations as an indispensable condition for the theory and simulations," http://arxiv.org/abs/1407.3244.
 H. C. Plummer, "On the problem of distribution in globular star clusters," Monthly Notices of the Royal Astronomical Society, vol. 71, pp. 460-470, 1911.
 K. Dolag, S. Borgani, S. Schindler, A. Diaferio, and A. M. Bykov, "Simulation techniques for cosmological simulations," Space Science Reviews, vol. 134, no. 1-4, pp. 229-268, 2008.
 F. Iannuzzi and K. Dolag, "Adaptive gravitational softening in GADGET," Monthly Notices of the Royal Astronomical Society, vol. 417, no. 4, pp. 2846-2859, 2011.
 J. E. Barnes, "Gravitational softening as a smoothing operation," Monthly Notices of the Royal Astronomical Society, vol. 425, no. 2, pp. 1104-1120, 2012.
 B. Rottgers, T. Naab, and L. Oser, "Stellar orbits in cosmological galaxy simulations: the connection to formation history and line-of-sight kinematics," Monthly Notices of the Royal Astronomical Society, vol. 445, no. 2, pp. 1065-1083, 2014.
 L. Hernquist, "An analytical model for spherical galaxies and bulges," Astrophysical Journal, vol. 356, no. 2, pp. 359-364, 1990.
 L. D. Landau and E. M. Lifshitz, Mechanics, vol. 1 of Course of Theoretical Physics Series, Oxford Pergamon Press, Oxford, UK, 3rd edition, 2000.
 K. S. Oh, D. N. C. Lin, and S. J. Aarseth, "On the tidal disruption of dwarf spheroidal galaxies around the galaxy," The Astrophysical Journal, vol. 442, no. 1, pp. 142-158, 1995.
CREST and NASA Research Centers, North Carolina Central University, Fayetteville Street 1801, Durham, NC 27707, USA
Correspondence should be addressed to Maxim Eingorn; email@example.com
Received 3 October 2014; Revised 30 November 2014; Accepted 4 December 2014; Published 22 December 2014
Academic Editor: Elias C. Vagenas
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article; particle-particle|
|Publication:||Advances in High Energy Physics|
|Date:||Jan 1, 2014|
|Previous Article:||Radio-wave propagation in salt domes: implications for a UHE cosmic neutrino detector.|
|Next Article:||Present status and future perspectives of the NEXT experiment.|