Effect of the Cubic Torus Topology on Cosmological Perturbations

Abstract

We study the effect of the cubic torus topology of the Universe on scalar cosmological perturbations which define the gravitational potential. We obtain three alternative forms of the solution for both the gravitational potential produced by point-like masses, and the corresponding force. The first solution includes the expansion of delta-functions into Fourier series, exploiting periodic boundary conditions. The second one is composed of summed solutions of the Helmholtz equation for the original mass and its images. Each of these summed solutions is the Yukawa potential. In the third formula, we express the Yukawa potentials via Ewald sums. We show that for the present Universe, both the bare summation of Yukawa potentials and the Yukawa-Ewald sums require smaller numbers of terms to yield the numerical values of the potential and the force up to desired accuracy. Nevertheless, the Yukawa formula is yet preferable owing to its much simpler structure.

1. Introduction

Spatial topology of the Universe is the fundamental problem of modern cosmology. Is the Universe spatially flat, open, or closed? Moreover, is it simply or multiply connected? This issue was and is the subject of debate in many scientific articles. Theory, e.g., General Relativity, does not provide direct answers to these questions, hence it is observation instead that plays the decisive role here. For instance, detection of multiple images of the same object would directly indicate that the space is multiply connected. Furthermore, such an extended object as the last scattering surface can self-intersect along pairs of circles, the so-called circles-in-the-sky [1,2,3,4]. These pairs of matched circles have the same temperature fluctuation distribution. While nearly antipodal circles-in-the-sky have not yet been revealed in the CMB radiation maps, the analysis of CMB anisotropies for repeated patterns is very promising, and even detection of a single pair of matched circles could confirm the flatness of the Universe together with its multi-connectedness [5,6].

At large angular scales, there are observable CMB anomalies in the form of quadrupole moment suppression and the quadrupole and octopole alignment. From the topological point of view, it is natural to explain the suppression by the absence of long wavelengths in sufficiently compact spaces. Such spaces should also have all dimensions of the same order of magnitude, thus being well-proportioned [7,8]. A cubic torus $T^3$ represents a glaring example of well-proportioned compact space, and a more general rectangular type $T\times T\times T$ (with unequal periods of tori) can bring forth a symmetry plane or a symmetry axis in the CMB pattern [9].

A multiply connected flat space with toroidal topology $T^3$ is compact in all directions and has a finite volume. There are observational limits on the size of such Universe, including the ones coming from the analysis of 7-year and 9-year WMAP temperature maps [8,10]. For example, according to the 7-year WMAP data, the lower bound on the size of the fundamental topological domain is 27.9 Gpc [11]. In the case of $T^3$ topology, more recent Planck mission data give $R_i>0.92{\chi }_{\mathrm{rec}}$ and $R_i>0.97{\chi }_{\mathrm{rec}}$ for Planck 2013 and 2015 results, respectively, where $R_i$ is the radius of the largest sphere which can be inscribed in the topological domain, and ${\chi }_{\mathrm{rec}}\sim 14$ Gpc is the distance to the recombination surface [12,13]. Based on these results, the Planck Collaboration has reported that currently, there is no detection of compact topology with a characteristic scale being less than the last scattering surface diameter. Meanwhile, as pointed out in [14], it is quite possible that the Universe does have compact topology, detectable through the values of observable parameters which lie outside the ranges covered by the WMAP and Planck missions (at least with respect to the circles-in-the-sky search).

Symmetries associated with the cubic torus $T^3$ are inherent in cosmological simulations of the large scale structure formation. Indeed, the cosmological N-body problem is almost always numerically solved in a cubic box with periodic boundary conditions [15,16,17,18,19,20,21,22,23,24,25]. In view of computational limitations, the edge of the simulation box is smaller than the lower experimental bound on the torus period and ranges from hundreds to thousands of Mpc. To perform such simulations, we need to know the form of metric perturbations, in particular, the expression for the gravitational potential generated by discrete masses. Such an investigation was performed, e.g., in [26,27] where the Authors considered a toroidal lattice with period L and equal masses M placed at the center of each cell. They found a solution of the Einstein equations, expanded into series in powers of the small parameter ${\left(M/L\right)}^{1/2}$. It turns out that in this case, the discrete mass distribution is characterized by non-convergent series. The inherent ultraviolet (UV) divergence is related to the point-like nature of the investigated matter sources. In order to avoid this problem, the Authors provided the masses with a small finite extension, thus introducing a UV cutoff scale. If the Schwarzschild radius of the masses is chosen as this cutoff scale, then one needs the first ${10}^9$ summands in the series. Meanwhile, if a typical galaxy dimension is chosen instead, then only first 200 summands are needed for an accurate description of the exterior solution. This cosmological model is characterized by a number of apparent limitations. First of all, clustering in the real Universe is much more complicated and irregular. In the second place, it is not difficult to show that the obtained solution with point-like gravitating masses has no definite values on the straight lines joining identical masses in neighboring cells, i.e., at points where masses are absent. The only way to avoid this problem and get a regular solution at any point of the cell is to perform the smearing of these masses over some region, i.e., to employ again a UV cutoff. Exactly the same situation takes place for the gravitational potential as a solution of the Poisson equation with periodic boundary conditions [28].

The situation is drastically changed if the gravitational potential satisfies the Helmholtz-type equation, as it takes place within the cosmic screening approach [29,30,31,32,33]. Careful analysis of the perturbed Einstein equations reveals that first-order cosmological perturbations (e.g., the gravitational potential) satisfy the Helmholtz equation. This relativistic effect arises due to the interaction between the gravitational field and the nonzero cosmological background. In the present paper, we analyze this equation in the case of periodic boundary conditions usually assumed for cosmological N-body simulations. In other words, we investigate the impact of the cubic torus topology on the shape of the gravitational potential. We present three alternative expressions for the potential. The main purpose of this paper is to determine among these solutions the one that is most advantageous with respect to numerical applications. Namely, to find which of the solutions requires less terms in series to attain the necessary precision. Our investigation shows that the solution based on Yukawa-type potentials is preferable, provided that the screening length is smaller than the period of the torus. This condition is in agreement with observational bounds.

The paper is structured as follows. In Section 2, we introduce the general setup of the model and present three alternative solutions for the gravitational potential for cubic torus topology. Section 3 and Section 4 are devoted to the detailed study of these potentials and the corresponding forces, respectively, in view of their usefulness for numerical computations. In Section 5, we summarize the obtained results.

2. The Model and Alternative Solutions

We consider the $\mathsf{\Lambda }$CDM model where matter (cold dark and baryonic) is taken in the form of point-like gravitating masses $m_n$. These inhomogeneities perturb the background Friedmann-Lemaître-Robertson-Walker metric. In the conformal Newtonian gauge, the perturbed metric reads [34,35]

$$d{s}^2=a^2\left[(1+2\mathsf{\Phi })d{\eta }^2-(1-2\mathsf{\Phi }){\delta }_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }\right].$$

(1)

The first-order scalar perturbation $\mathsf{\Phi }$, $\left|\mathsf{\Phi }\right|\ll 1$, defines the total gravitational potential of the system and satisfies the equation [33]

$$\mathsf{\Delta }\mathsf{\Phi }-\frac{a^2}{{\lambda }_{\mathrm{eff}}^2}\mathsf{\Phi }=\frac{\kappa c^2}{2a}\left(\rho -\overline{\rho }\right),$$

(2)

where $\kappa \equiv 8\pi G_N/c^4$ ($G_N$ is the Newtonian gravitational constant and c represents the speed of light), $a\left(\eta \right)$ is the scale factor, and $\mathsf{\Delta }$ denotes the Laplace operator in comoving coordinates. In addition,

$$\rho =\sum _n{m}_n\delta (\mathbf{r}-{\mathbf{r}}_n)$$

(3)

is the comoving mass density and $\overline{\rho }=\mathrm{const}$ is its average value. The Helmholtz Equation (2) was derived in [33] within the cosmic screening approach [29,30,31,32] and effectively takes into account peculiar velocities of inhomogeneities. The effective screening length

$${\lambda }_{\mathrm{eff}}=\sqrt{\frac{c^2{a}^{2}H}{3}\int \frac{da}{a^3{H}^3}},$$

(4)

where $H=\left(c/a^2\right)da/d\eta$ is the Hubble parameter. It can be easily seen that ${\lambda }_{\mathrm{eff}}$ admits the time dependence. If we substitute the cosmological parameters according to the Planck 2018 data [36], i.e., $H_0=67.4\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, ${\mathsf{\Omega }}_{\mathrm{M}}=0.315$, ${\mathsf{\Omega }}_{\mathsf{\Lambda }}=0.685$, at the present time we get ${\left({\lambda }_{\mathrm{eff}}\right)}_0=2.57$ Gpc [33].

It is convenient to introduce a shifted gravitational potential

$$\widehat{\mathsf{\Phi }}\equiv \mathsf{\Phi }-{\lambda }_{\mathrm{eff}}^2\frac{\kappa c^2}{2a^3}\overline{\rho },$$

(5)

which satisfies

$$\mathsf{\Delta }\widehat{\mathsf{\Phi }}-\frac{a^2}{{\lambda }_{\mathrm{eff}}^2}\widehat{\mathsf{\Phi }}=\frac{\kappa c^2}{2a}\rho .$$

(6)

The superposition principle allows us to solve this equation for a selected particle and then, to simply re-express the solution for a system of randomly distributed particles.

Herein we intend to solve this equation in the case of the three-torus topology $T\times T\times T$ with periods $l_1,l_2$ and $l_3$. Obviously, $\overline{\rho }={\sum }_n{m}_n/\left(l_1{l}_2{l}_3\right)$. It is worth noting that the form of the equation is determined by General Relativity with the appropriate choice of the metric and energy-momentum tensor of matter. Therefore, we have the same equation for both flat simply-connected and multiply-connected topologies. Evidently, the solution of this equation depends on the boundary conditions. First, we find the solution for the selected particle m, chosen to be, without loss of generality, at the center of Cartesian coordinates. For the given source, Equation (6) becomes

$$\mathsf{\Delta }\widehat{\mathsf{\Phi }}-\frac{a^2}{{\lambda }_{\mathrm{eff}}^2}\widehat{\mathsf{\Phi }}=\frac{\kappa c^2}{2a}m\delta \left(x\right)\delta \left(y\right)\delta \left(z\right).$$

(7)

Toroidal topology also implies periodic boundary conditions. Therefore, the expansion of the delta-function into Fourier series reads

$$\delta \left(x\right)=\frac{1}{l_1}\sum _{k_1=-\infty }^{+\infty }cos\left(\frac{2\pi k_1}{l_1}x\right)$$

(8)

and similar expressions apply for $\delta \left(y\right)$ and $\delta \left(z\right)$. Employing this delta-function presentation, it can be easily verified that the solution of Equation (7) is

$$\begin{array}{ccc}\hfill \widehat{\mathsf{\Phi }}& =& -\frac{\kappa c^2}{2a}\frac{m}{l_1{l}_2{l}_3}\sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }{\left[4{\pi }^2\left(\frac{k_1^2}{l_1^2}+\frac{k_2^2}{l_2^2}+\frac{k_3^2}{l_3^2}\right)+\frac{a^2}{{\lambda }_{\mathrm{eff}}^2}\right]}^{-1}\hfill \\ & \times & cos\left(\frac{2\pi k_1}{l_1}x\right)cos\left(\frac{2\pi k_2}{l_2}y\right)cos\left(\frac{2\pi k_3}{l_3}z\right).\hfill \end{array}$$

(9)

Thus, for a system of arbitrarily located massive particles in a cell, the total gravitational potential is

$$\begin{array}{ccc}\hfill \widehat{\mathsf{\Phi }}& =& -\frac{\kappa c^2}{2a}\frac{1}{l_1{l}_2{l}_3}\sum _n{m}_n\left\{\sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }{\left[4{\pi }^2\left(\frac{k_1^2}{l_1^2}+\frac{k_2^2}{l_2^2}+\frac{k_3^2}{l_3^2}\right)+\frac{a^2}{{\lambda }_{\mathrm{eff}}^2}\right]}^{-1}\right.\hfill \\ & \times & \left.cos\left[\frac{2\pi k_1}{l_1}(x-x_n)\right]cos\left[\frac{2\pi k_2}{l_2}(y-y_n)\right]cos\left[\frac{2\pi k_3}{l_3}(z-z_n)\right]\right\}.\hfill \end{array}$$

(10)

The obtained solutions satisfy two important natural conditions. First, Equation (9) yields the correct Newtonian limit in the close vicinity of the source particle. Second, using the relation (5), it can be demonstrated that the average value of $\mathsf{\Phi }$ is equal to zero, as is required of fluctuations at the first-order level. It is worth noting that the sum of Newtonian potentials does not satisfy this condition (see also the reasoning in [37]).

The solution of Equation (7) can also be found in the alternative way. Owing to periodic boundary conditions, each mass in the fundamental cell has its counterparts shifted by multiples of tori periods $l_1,l_2$ and $l_3$. Therefore, we may solve Equation (7) by merely counting the distinct contributions of these images. Since this is a Helmholtz-type equation, the solution is the sum of the corresponding Yukawa potentials:

$$\begin{array}{ccc}\hfill \widehat{\mathsf{\Phi }}& =& -\frac{\kappa c^{2}m}{8\pi a}\sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }\frac{1}{\sqrt{{(x-k_1{l}_1)}^2+{(y-k_2{l}_2)}^2+{(z-k_3{l}_3)}^2}}\hfill \\ & \times & exp\left(-\frac{a\sqrt{{(x-k_1{l}_1)}^2+{(y-k_2{l}_2)}^2+{(z-k_3{l}_3)}^2}}{{\lambda }_{\mathrm{eff}}}\right).\hfill \end{array}$$

(11)

We rewrite the alternative solutions (9) and (11) as

$$\begin{array}{ccc}\hfill {\tilde{\mathsf{\Phi }}}_{cos}\equiv {\left(-\frac{G_{N}m}{c^{2}al}\right)}^{-1}{\widehat{\mathsf{\Phi }}}_{cos}& =& \sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }{\left[\pi \left(k_1^2+k_2^2+k_3^2\right)+\frac{1}{4\pi {\tilde{\lambda }}_{\mathrm{eff}}^2}\right]}^{-1}\hfill \\ & \times & cos\left(2\pi k_1\tilde{x}\right)cos\left(2\pi k_2\tilde{y}\right)cos\left(2\pi k_3\tilde{z}\right)\hfill \end{array}$$

(12)

and

$$\begin{array}{ccc}\hfill {\tilde{\mathsf{\Phi }}}_{exp}\equiv {\left(-\frac{G_{N}m}{c^{2}al}\right)}^{-1}{\widehat{\mathsf{\Phi }}}_{exp}& =& \sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }\frac{1}{\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}}\hfill \\ & \times & exp\left(-\frac{\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}}{{\tilde{\lambda }}_{\mathrm{eff}}}\right),\hfill \end{array}$$

(13)

where, for simplicity, we consider an equal-sided cubic torus with $l_1=l_2=l_3\equiv l$ and introduce the notation

$$x=\tilde{x}l,y=\tilde{y}l,z=\tilde{z}l,{\lambda }_{\mathrm{eff}}={\tilde{\lambda }}_{\mathrm{eff}}al.$$

(14)

Yukawa potentials with periodic boundaries can also be expressed in the form of Ewald sums, i.e., as two distinct rapidly converging series, each of which exists in one of the real and Fourier spaces. This type of presentation is usually used in depicting electrostatic interactions in plasma, colloids etc., and for this purpose, the Yukawa potential for systems with three-dimensional periodicity was obtained earlier in [38]. In the cosmological framework, the Yukawa-Ewald potential for gravitational interactions takes the form

$$\begin{array}{ccc}\hfill {\tilde{\mathsf{\Phi }}}_{\mathrm{mix}}& \equiv & {\left(-\frac{G_{N}m}{c^{2}al}\right)}^{-1}{\widehat{\mathsf{\Phi }}}_{\mathrm{mix}}\hfill \\ & =& \sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }\left\{\frac{D\left(\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2};\alpha ;{\tilde{\lambda }}_{\mathrm{eff}}\right)}{2\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}}\right.\hfill \\ & +& \left.4\pi cos\left[2\pi \left(k_1\tilde{x}+k_2\tilde{y}+k_3\tilde{z}\right)\right]\frac{exp\left[-\left(4{\pi }^2{k}^2+{\tilde{\lambda }}_{\mathrm{eff}}^{-2}\right)/\left(4{\alpha }^2\right)\right]}{4{\pi }^2{k}^2+{\tilde{\lambda }}_{\mathrm{eff}}^{-2}}\right\},\hfill \end{array}$$

(15)

where $k^2\equiv k_1^2+k_2^2+k_3^2$,

$$\begin{array}{ccc}\hfill & & D\left(\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2};\alpha ;{\tilde{\lambda }}_{\mathrm{eff}}\right)\hfill \\ & \equiv & exp\left(\frac{\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}}{{\tilde{\lambda }}_{\mathrm{eff}}}\right)\hfill \\ & \times & \mathrm{erfc}\left(\alpha \sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}+\frac{1}{2\alpha {\tilde{\lambda }}_{\mathrm{eff}}}\right)\hfill \\ & +& exp\left(-\frac{\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}}{{\tilde{\lambda }}_{\mathrm{eff}}}\right)\hfill \\ & \times & \mathrm{erfc}\left(\alpha \sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}-\frac{1}{2\alpha {\tilde{\lambda }}_{\mathrm{eff}}}\right).\hfill \end{array}$$

(16)

In these expressions, erfc represents the complementary error function and $\alpha$, the free parameter, is to be assigned the optimal value to save computational effort while operating with adequate accuracy. Below we will test a number of values of $\alpha$. Our research demonstrates that for the chosen range of ${\tilde{\lambda }}_{\mathrm{eff}}$, the optimal one is around 2.

We note that alternative expressions for the gravitational potential were also found in the cases of slab and chimney topologies [39,40]. Direct comparison of the obtained formulas for different topologies shows that only the Formula (13) above (Yukawa potentials) may be interpreted as a simple extension of the “slab” and “chimney” counterparts (2.28) in [39] and (18) in [40], respectively. However, obviously, both Formulas (12) and (15) are quite different from their counterparts, and this is a nontrivial task to derive them from the previous papers. Moreover, the Yukawa-Ewald potential is not presented in the case of slab topology. Thus, formulas found in the present paper are new and otherwise absent in the literature. As regards the gravitational forces derived below, again, the similarity is present only in the case of the Yukawa formula, but both alternatives are new and cannot be easily derived from the previous results on different topologies.

The obtained solutions (12), (13) and (15) depend on time. It is important to note that they all satisfy the same Helmholtz equation and are three representations of the solution. In our papers [29,30,31,33] it was investigated in detail that such a solution satisfies the complete system of perturbed Einstein equations.

3. Gravitational Potentials

In the previous section, we have obtained three alternative formulas for the gravitational potential created by a point-like particle placed at the center of Cartesian coordinates $(x,y,z)=(0,0,0)$ and by its infinitely many images placed at points $(x,y,z)=(k_{1}l,k_{2}l,k_{3}l)$ where $k_{1,2,3}=0,\pm 1,\pm 2,\dots$ Due to periodic boundary conditions, these formulas include infinite series. We aim to find out which of these expressions requires fewer terms in the series sum to yield the value of the potential to given accuracy. The less the number n of these terms, the more advantageous the corresponding formula for numerical calculations. This number n is defined via the condition that the ratio $|\mathrm{exact}\tilde{\mathsf{\Phi }}-\mathrm{approximate}\tilde{\mathsf{\Phi }}|/|\mathrm{exact}\tilde{\mathsf{\Phi }}|$ is either equal to or less than 0.001. This is our demanded level of precision in determining the approximate value of $\tilde{\mathsf{\Phi }}$. Each of the alternative expressions (12), (13) and (15) has its own number n designated as $n_{\cos },n_{\exp }$ and $n_{\mathrm{mix}}$, correspondingly. Evidently, the formula with the smallest number is, all other things being equal, the most convenient for numerical computations. All three formulas to be compared contain triple series. Consequently, the sought values n correspond to the minimum number of triplets $(k_1,k_2,k_3)$ included in series for which the required precision is achieved. To find this number, we generate a sequence in increasing order of $\sqrt{k_1^2+k_2^2+k_3^2}$ in Mathematica [41] and count the number n of terms involved in it.

We calculate the potentials (12), (13) and (15) in Mathematica [41] up to the adopted accuracy at a selection of points in the cell and display the results in

Table 1. Rescaled potential $\tilde{\mathsf{\Phi }}$ and the corresponding numbers $n_{exp}$ of terms in the series sum at nine points in the cell for ${\tilde{\lambda }}_{\mathrm{eff}}=0.01$ (left chart) and ${\tilde{\lambda }}_{\mathrm{eff}}=0.1$ (right chart).

${\tilde{\mathit{\lambda }}}_{\mathbf{eff}}\mathbf{=}\mathbf{0.01}$						${\tilde{\mathit{\lambda }}}_{\mathbf{eff}}\mathbf{=}\mathbf{0.1}$
	$\tilde{\mathit{x}}$	$\tilde{\mathit{y}}$	$\tilde{\mathit{z}}$	$\tilde{\mathsf{\Phi }}$	${\mathit{n}}_{exp}$		$\tilde{\mathit{x}}$	$\tilde{\mathit{y}}$	$\tilde{\mathit{z}}$	$\tilde{\mathsf{\Phi }}$	${\mathit{n}}_{exp}$
$A_1$	0.5	0	0.5	$1.105\times {10}^{-30}$	9	$A_1$	0.5	0	0.5	$4.837\times {10}^{-3}$	20
$A_2$	0.5	0	0.1	$2.810\times {10}^{-22}$	2	$A_2$	0.5	0	0.1	$2.406\times {10}^{-2}$	9
$A_3$	0.5	0	0	$7.715\times {10}^{-22}$	2	$A_3$	0.5	0	0	$2.705\times {10}^{-2}$	9
$B_1$	0.1	0	0.1	$5.101\times {10}^{-6}$	1	$B_1$	0.1	0	0.1	$1.719$	1
$B_2$	0.1	0	0	$4.540\times {10}^{-4}$	1	$B_2$	0.1	0	0	$3.679$	1
$C_1$	0.5	0.5	0.5	$2.262\times {10}^{-37}$	20	$C_1$	0.5	0.5	0.5	$1.602\times {10}^{-3}$	20
$C_2$	0.5	0.5	0.1	$5.413\times {10}^{-31}$	8	$C_2$	0.5	0.5	0.1	$4.478\times {10}^{-3}$	20
$C_3$	0.1	0.5	0.1	$1.044\times {10}^{-22}$	3	$C_3$	0.1	0.5	0.1	$2.146\times {10}^{-2}$	12
$C_4$	0.1	0.1	0.1	$1.735\times {10}^{-7}$	1	$C_4$	0.1	0.1	0.1	$1.022$	1

and

Table 2. Rescaled potential $\tilde{\mathsf{\Phi }}$ and the corresponding numbers $n_{exp}$ and $n_{\mathrm{mix}}$ of terms in the series sum at nine points for ${\tilde{\lambda }}_{\mathrm{eff}}=1$ (top chart) and ${\tilde{\lambda }}_{\mathrm{eff}}=5$ (bottom chart).

	$\tilde{\mathit{x}}$	$\tilde{\mathit{y}}$	$\tilde{\mathit{z}}$	$\tilde{\mathsf{\Phi }}$	${\mathit{n}}_{exp}$	${\mathit{n}}_{\mathbf{mix}}^{\mathit{\alpha }\mathbf{=}\mathbf{1}}$	${\mathit{n}}_{\mathbf{mix}}^{\mathit{\alpha }\mathbf{=}\mathbf{2}}$	${\mathit{n}}_{\mathbf{mix}}^{\mathit{\alpha }\mathbf{=}\mathbf{3}}$
$A_1$	0.5	0	0.5	$12.00$	3449	64	9	25
$A_2$	0.5	0	0.1	$12.42$	3352	48	7	22
$A_3$	0.5	0	0	$12.47$	3345	47	7	22
$B_1$	0.1	0	0.1	$16.77$	2965	37	7	23
$B_2$	0.1	0	0	$19.66$	2794	33	7	24
$C_1$	0.5	0.5	0.5	$11.79$	3510	64	20	24
$C_2$	0.5	0.5	0.1	$11.98$	3451	62	8	24
$C_3$	0.1	0.5	0.1	$12.37$	3370	54	7	21
$C_4$	0.1	0.1	0.1	$15.50$	3058	39	7	22
	$\tilde{\mathit{x}}$	$\tilde{\mathit{y}}$	$\tilde{\mathit{z}}$	$\tilde{\mathsf{\Phi }}$	${\mathit{n}}_{\mathrm{mix}}^{\mathbf{\alpha }=\mathbf{1}}$	${\mathit{n}}_{\mathrm{mix}}^{\mathbf{\alpha }=\mathbf{2}}$	${\mathit{n}}_{\mathrm{mix}}^{\mathbf{\alpha }=\mathbf{2}.\mathbf{5}}$	${\mathit{n}}_{\mathrm{mix}}^{\mathbf{\alpha }=\mathbf{5}}$
$A_1$	0.5	0	0.5	$313.6$	14	1	1	17
$A_2$	0.5	0	0.1	$314.0$	11	2	1	16
$A_3$	0.5	0	0	$314.1$	11	2	1	17
$B_1$	0.1	0	0.1	$318.4$	12	1	3	23
$B_2$	0.1	0	0	$321.3$	12	1	4	29
$C_1$	0.5	0.5	0.5	$313.4$	20	1	1	12
$C_2$	0.5	0.5	0.1	$313.6$	13	1	1	6
$C_3$	0.1	0.5	0.1	$314.0$	13	2	1	11
$C_4$	0.1	0.1	0.1	$317.2$	12	1	3	17

. The number $n_{\exp }$ is defined employing Equation (13): for any $n>n_{\exp }$, the approximate ${\tilde{\mathsf{\Phi }}}_{exp}$ will be determined with better accuracy than a tenth of a percent. The values $n_{\cos }$ and $n_{\mathrm{mix}}$ follow from the Formulas (12) and (15) under the condition that the gravitational potential is calculated with the same accuracy at the point of interest. We have found that the Yukawa-Ewald formula (15) works well (i.e., requires the smallest number of terms) both for small and large selected values of the screening length ${\tilde{\lambda }}_{\mathrm{eff}}$. Therefore, we evaluate the exact $\tilde{\mathsf{\Phi }}$ by the Formula (15) for $n\gg n_{\mathrm{mix}}$. Additionally, we have observed that the use of Equation (12) to get $n_{\cos }$ yields faulty outputs due to problematic aspects of the computational process. The “trigonometric” Formula (12) contains an alternating series. The summation of such a series is accompanied by significant round-off errors and to reach the required accuracy, when possible, it is necessary to take into account a very large number of terms (more than ${10}^5$ in our case). Therefore, the use of this formula looks absolutely unreasonable in comparison with the rapidly converging expressions (13) and (15). Hence, this trigonometric formula is not suitable for numerical calculations, and the related values are excluded from the tables.

As follows from the Formulas (12), (13) and (15), the resulting values of the potential are sensitive to the choice of ${\tilde{\lambda }}_{\mathrm{eff}}$. For our calculations, we choose four different values, that are ${\tilde{\lambda }}_{\mathrm{eff}}=0.01,0.1,1$ and 5. The rescaled screening length ${\tilde{\lambda }}_{\mathrm{eff}}$ is the ratio of the physical effective screening length ${\lambda }_{\mathrm{eff}}$ to the physical size of the period $al$. As we have mentioned previously, today ${\lambda }_{\mathrm{eff}}\sim 2.6$ Gpc for the $\mathsf{\Lambda }$CDM model [33], and the size of the fundamental domain for the cubic torus topology is restricted by Planck 2015 results to no less than $al\sim$ 27 Gpc [13], i.e., the observational data require that ${\tilde{\lambda }}_{\mathrm{eff}}\ll 1$. Nevertheless, taking into account that many N-body simulations are indeed performed in boxes with sizes less than 1 Gpc, we also consider ${\tilde{\lambda }}_{\mathrm{eff}}\ge 1$.

The Yukawa-Ewald potential (15) is sensitive to the free parameter $\alpha$ as well. Therefore, choosing different values of this quantity, we also seek those at which $n_{\mathrm{mix}}$ will be minimum. Our calculations demonstrate that for the chosen range of $0.01\le {\tilde{\lambda }}_{\mathrm{eff}}\le 5$, the optimal value of $\alpha$ is around 2.

In Table 1, we give the minimum numbers $n_{exp}$ of terms in the series (13) that return the value of the gravitational potential up to the adopted accuracy. The left and right charts correspond to ${\tilde{\lambda }}_{\mathrm{eff}}=0.01$ and ${\tilde{\lambda }}_{\mathrm{eff}}=0.1$, respectively. As for the values $n_{\mathrm{mix}}$, these numbers are the same as $n_{exp}$ for both left and right charts as long as $\alpha$ lies between ${10}^{-3}$ and 2. However, outside this interval, the Yukawa-Ewald formula may require more summands. For example, if ${\tilde{\lambda }}_{\mathrm{eff}}=0.1$ and $\alpha =2.5$, then $n_{\mathrm{mix}}=29,23,23,4,1,48,28,23,6$ for the points $A_1,A_2,\dots ,C_4$, respectively. All in all, when ${\tilde{\lambda }}_{\mathrm{eff}}\ll 1$ (in accordance with observational bounds), both the Yukawa (13) and Yukawa-Ewald (15) formulas demonstrate good results since they need less terms in the series sum. The potential expression in (13) is much simpler, though. Thus, from that aspect, the Yukawa formula is a more practical tool for computational purposes in the case of small ${\tilde{\lambda }}_{\mathrm{eff}}$.

In Table 2, we present the results of similar calculations for ${\tilde{\lambda }}_{\mathrm{eff}}=1$ and ${\tilde{\lambda }}_{\mathrm{eff}}=5$. In the case ${\tilde{\lambda }}_{\mathrm{eff}}=1$ (top chart), $n_{exp}\gg n_{\mathrm{mix}}$ for all selected points, and the optimal choice for the parameter $\alpha$ in the Yukawa-Ewald formula (15) is 2. For ${\tilde{\lambda }}_{\mathrm{eff}}=5$ (bottom chart), as regards the Yukawa Formula (13), the $n_{exp}$ values are ill-suited (they are extremely large). Therefore, we do not show them here. The value $\alpha =2$ is again the optimal choice for the Yukawa-Ewald potential. Hence, when ${\tilde{\lambda }}_{\mathrm{eff}}\ge 1$, the Yukawa-Ewald formula (15) delivers the best performance in numerical calculations.

We demonstrate in Figure 1 and Figure 2 (plotted in Mathematica [41]) the $\tilde{z}=0$ sections of the rescaled potential $\tilde{\mathsf{\Phi }}$ for different values of the rescaled effective screening length ${\tilde{\lambda }}_{\mathrm{eff}}$ considered in Table 1 and Table 2. To plot both figures, we use the Yukawa-Ewald formula (15) for $\alpha =2$ and $n\gg n_{\mathrm{mix}}$.

4. Gravitational Forces

In this section we provide the gravitational force formulas associated with the alternative expressions (12), (13) and (15) for the gravitational potential. Clearly, it is sufficient to consider the force projection onto one of the axes in Cartesian coordinates. Let it be the x-axis. Then, the projections read

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial \tilde{x}}\left({\tilde{\mathsf{\Phi }}}_{cos}\right)& =& -2\pi \sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }{\left[\pi \left(k_1^2+k_2^2+k_3^2\right)+\frac{1}{4\pi {\tilde{\lambda }}_{\mathrm{eff}}^2}\right]}^{-1}\hfill \\ & \times & k_{1}sin\left(2\pi k_1\tilde{x}\right)cos\left(2\pi k_2\tilde{y}\right)cos\left(2\pi k_3\tilde{z}\right),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill & & \frac{\partial }{\partial \tilde{x}}\left({\tilde{\mathsf{\Phi }}}_{exp}\right)=-\sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }exp\left(-\frac{\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}}{{\tilde{\lambda }}_{\mathrm{eff}}}\right)\hfill \\ & \times & \left\{\frac{\tilde{x}-k_1}{{\left[{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2\right]}^{3/2}}+\phantom{\frac{\tilde{x}-k_1}{{\left[{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2\right]}^{3/2}}}\frac{\tilde{x}-k_1}{{\tilde{\lambda }}_{\mathrm{eff}}\left[{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2\right]}\right\},\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill & & \frac{\partial }{\partial \tilde{x}}\left({\tilde{\mathsf{\Phi }}}_{\mathrm{mix}}\right)=\hfill \\ & & -\frac{1}{2}\sum _{k_1=-\infty }^{+\infty }\sum _{k_2=-\infty }^{+\infty }\sum _{k_3=-\infty }^{+\infty }\left\{\frac{\left(\tilde{x}-k_1\right)D\left(\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2};\alpha ;{\tilde{\lambda }}_{\mathrm{eff}}\right)}{{[{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2]}^{3/2}}\right.\hfill \\ & +& C_{-}\frac{\left(\tilde{x}-k_1\right)exp\left(-\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}/{\tilde{\lambda }}_{\mathrm{eff}}\right)}{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}\hfill \\ & +& C_{+}\frac{\left(\tilde{x}-k_1\right)exp\left(\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}/{\tilde{\lambda }}_{\mathrm{eff}}\right)}{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}\hfill \\ & +& \left.16{\pi }^2{k}_{1}sin\left[2\pi \left(k_1\tilde{x}+k_2\tilde{y}+k_3\tilde{z}\right)\right]\frac{exp\left[-\left(4{\pi }^2{k}^2+{\tilde{\lambda }}_{\mathrm{eff}}^{-2}\right)/\left(4{\alpha }^2\right)\right]}{4{\pi }^2{k}^2+{\tilde{\lambda }}_{\mathrm{eff}}^{-2}}\right\},\hfill \end{array}$$

(19)

where

$$\begin{array}{ccc}\hfill C_{\mp }& =& C_{\mp }\left(\sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2};\alpha ;{\tilde{\lambda }}_{\mathrm{eff}}\right)\hfill \\ & \equiv & \frac{2\alpha }{\sqrt{\pi }}exp\left[-{\left(\alpha \sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}\mp \frac{1}{2\alpha {\tilde{\lambda }}_{\mathrm{eff}}}\right)}^2\right]\hfill \\ & \pm & \frac{1}{{\tilde{\lambda }}_{\mathrm{eff}}}\mathrm{erfc}\left(\alpha \sqrt{{(\tilde{x}-k_1)}^2+{(\tilde{y}-k_2)}^2+{(\tilde{z}-k_3)}^2}\mp \frac{1}{2\alpha {\tilde{\lambda }}_{\mathrm{eff}}}\right).\hfill \end{array}$$

(20)

The x-components of the force are zero at the points $A_1,A_2,A_3,C_1$ and $C_2$, so we calculate these components up to the adopted accuracy only at the points $B_1,B_2,C_3$ and $C_4$.

The results we arrive at are as follows: first, as also is the case for potentials, the trigonometric Formula (17) does not provide acceptable values because of the complications that arise during the computational stage. Next, if ${\tilde{\lambda }}_{\mathrm{eff}}=0.01$, then $n_{exp}=1,1,3,1$ for the points $B_1,B_2,C_3,C_4$, respectively, and the numbers $n_{\mathrm{mix}}$ are exactly the same (assuming here and in what follows that $\alpha =2$). In the case ${\tilde{\lambda }}_{\mathrm{eff}}=0.1$, we get $n_{exp}=1,1,12,1$ for these points, and the corresponding numbers $n_{\mathrm{mix}}$ are, again, identical. However, for ${\tilde{\lambda }}_{\mathrm{eff}}=1$, we have $n_{exp}=258,82,987,486$ while $n_{\mathrm{mix}}=7,7,21,9$. Finally, if ${\tilde{\lambda }}_{\mathrm{eff}}=5$, then still $n_{\mathrm{mix}}=7,7,21,9$, but $n_{exp}$ acquires unreasonably large values.

Our calculations demonstrate that the numbers $n_{exp}$ start to grow once ${\tilde{\lambda }}_{\mathrm{eff}}$ exceeds 0.1 and they acquire large values as ${\tilde{\lambda }}_{\mathrm{eff}}$ approaches 1, while $n_{\mathrm{mix}}$ remain small throughout (provided that the value of the parameter $\alpha$ is optimal, e.g., for $\alpha =2$ in the above cases). The Yukawa formula is, after all, an attractive option in view of its simpler structure. The Yukawa-Ewald formula is much more complicated and consequently, it takes longer to numerically calculate the potentials and forces for comparable values of $n_{\exp }$ and $n_{\mathrm{mix}}$. There exists a moment, though, when the execution time for the Yukawa formula, with increased number of terms in the sum, is approximately equal to the one for the more complex Yukawa-Ewald formula. We have seen that with respect to gravitational forces, this takes place when $n_{\exp }$ is about 6 times larger than $n_{\mathrm{mix}}$.

In Figure 3 and Figure 4 (plotted in Mathematica [41]), we depict the $\tilde{z}=0$ sections of the x-components of gravitational forces for ${\tilde{\lambda }}_{\mathrm{eff}}=0.01$, $0.1$ and 1. We employ the Yukawa-Ewald formula (19) for $n\gg n_{\mathrm{mix}}$ and $\alpha =2$. For ${\tilde{\lambda }}_{\mathrm{eff}}=5$, the picture is similar to one in the case ${\tilde{\lambda }}_{\mathrm{eff}}=1$.

5. Conclusions

In this paper we have analyzed the influence of topology on the gravitational interaction in the Universe. We have considered a model in which the fundamental domain is a three-torus $T\times T\times T$. In such a space, gravitating masses are subject to periodic boundary conditions along three coordinate axes, that is, every mass in the fundamental domain has its counterparts in infinitely many cells shifted along each axis by multiples of tori periods. For this lattice Universe, we have obtained three alternative forms of the expression for the gravitational potential produced by a point-like mass. The first one (see Equation (12)) exploits the periodic structure of space: it involves the expansion of delta-functions into Fourier series. This solution is a trigonometric series, thus we named it trigonometric. The second one (see Equation (13)), the Yukawa solution, was obtained by directly summing the fields produced by the original mass and its images. Since the summed potentials satisfy the Helmholtz equation, they are the Yukawa potentials. In the third formula (see Equation (15)), we have expressed them via Ewald sums (the Yukawa-Ewald formula) and shown that in some cases, such a trick facilitates (despite the complex form of the resulting expression) numerical calculations. We have also presented the corresponding formulas for the x-component of the gravitational force (see Equations (17)–(19)). All these formulas, both for the potentials and forces, can be easily generalized for a system of arbitrarily located massive particles (see, e.g., (10)). It is well known that the gravitational potential of a system of masses distributed in the Universe defines the scalar perturbations of the metric.

A reasonable question to ask is, then, which of these formulas is indeed preferable for numerical applications for the given accuracy. Since all three expressions are sensitive to the rescaled effective screening length ${\tilde{\lambda }}_{\mathrm{eff}}={\lambda }_{\mathrm{eff}}/\left(al\right)$, we analyzed both small values and values equal to or greater than 1: $0.01\le {\tilde{\lambda }}_{\mathrm{eff}}\le 5$. The Yukawa-Ewald formula additionally admits a free parameter $\alpha$, and we have revealed that for the given range of ${\tilde{\lambda }}_{\mathrm{eff}}$, the optimal value of $\alpha$ is 2. Our calculations show that the trigonometric formula does not provide reasonable results for the potentials or forces as complications arise in the computational process. In the case of small ${\tilde{\lambda }}_{\mathrm{eff}}$ (as also demanded by the observational bounds) both the Yukawa and Yukawa-Ewald formulas deliver good results since they require rather small numbers of terms $n_{\exp }$ and $n_{\mathrm{mix}}$ to yield the values for the potentials and forces up to the given accuracy. Nevertheless, employing the Yukawa formula is a more convenient choice owing to its notably simpler structure. The situation is altered when ${\tilde{\lambda }}_{\mathrm{eff}}>0.1$: $n_{\exp }$ begins to increase quickly while $n_{\mathrm{mix}}$ still takes on rather small values. Therefore, for such ${\tilde{\lambda }}_{\mathrm{eff}}$, the Yukawa-Ewald formula is preferable instead.

Finally, we emphasize two important points. First, our results directly confirm that the undesirable impact of periodicity on simulation outputs can be weakened if the edge of the box (cubic torus period $al$) is set to be larger than the predicted Yukawa interaction range ${\lambda }_{\mathrm{eff}}$ (see Table 1). The Yukawa formula reflects the contribution of images to the value of the potential via the number $n_{\exp }$ (needed to reach the required accuracy). We can easily see that $n_{\exp }$ gets smaller (i.e., goes to 1) with decreasing ${\tilde{\lambda }}_{\mathrm{eff}}$. Second, operating with summed Yukawa potentials, we provide a reliable description of the inhomogeneous gravitational field generated by a toroidal lattice of point-like masses, avoiding non-convergent series. The obtained series converge at all points only except those where discrete masses themselves are located.