Advances in High Energy Physics

Volume 2014 (2014), Article ID 903642, 4 pages

http://dx.doi.org/10.1155/2014/903642

## Gurzadyan’s Problem 5 and Improvement of Softenings for Cosmological Simulations Using the PP Method

CREST and NASA Research Centers, North Carolina Central University, Fayetteville Street 1801, Durham, NC 27707, USA

Received 3 October 2014; Revised 30 November 2014; Accepted 4 December 2014; Published 22 December 2014

Academic Editor: Elias C. Vagenas

Copyright © 2014 Maxim Eingorn. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The publication of this article was funded by SCOAP^{3}.

#### 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, masslength^{3}).

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 .