Effect of the cubic torus topology on cosmological perturbations

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.


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 socalled 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 × T × 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χ rec and R i > 0.97χ 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 χ rec ∼ 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 (M/L) 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 Helmholtztype 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. Sections 3 and 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.

The Model and Alternative Solutions
We consider the Λ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] (1) The first-order scalar perturbation Φ, |Φ| 1, defines the total gravitational potential of the system and satisfies the equation [33] ∆Φ − a 2 where κ ≡ 8πG N /c 4 (G N is the Newtonian gravitational constant and c represents the speed of light), a(η) is the scale factor, and ∆ denotes the Laplace operator in comoving coordinates. In addition, is the comoving mass density andρ = 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 where H = c/a 2 da/dη is the Hubble parameter. It can be easily seen that λ eff admits the time dependence. If we substitute the cosmological parameters according to the Planck 2018 data [36], i.e., H 0 = 67.4 km s −1 Mpc −1 , Ω M = 0.315, Ω Λ = 0.685, at the present time we get (λ eff ) 0 = 2.57 Gpc [33]. It is convenient to introduce a shifted gravitational potential which satisfies 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 × T × T with periods l 1 , l 2 and l 3 . Obviously,ρ = ∑ n m n /(l 1 l 2 l 3 ). 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 Toroidal topology also implies periodic boundary conditions. Therefore, the expansion of the delta-function into Fourier series reads and similar expressions apply for δ(y) and δ(z). Employing this delta-function presentation, it can be easily verified that the solution of Equation (7) is Thus, for a system of arbitrarily located massive particles in a cell, the total gravitational potential is 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 Φ 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: We rewrite the alternative solutions (9) and (11) as and where, for simplicity, we consider an equal-sided cubic torus with l 1 = l 2 = l 3 ≡ l and introduce the notation 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 threedimensional periodicity was obtained earlier in [38]. In the cosmological framework, the Yukawa-Ewald potential for gravitational interactions takes the form where In these expressions, erfc represents the complementary error function and α, 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 α. Our research demonstrates that for the chosen range ofλ 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.

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, ±1, ±2, . . . 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 |exactΦ − approximateΦ|/|exactΦ| is either equal to or less than 0.001. This is our demanded level of precision in determining the approximate value ofΦ. Each of the alternative expressions (12), (13) and (15) has its own number n designated as n cos , n exp and n 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 k 2 1 + k 2 2 + k 2 3 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 Tables 1 and 2. The number n exp is defined employing Equation (13): for any n > n exp , the approximateΦ exp will be determined with better accuracy than a tenth of a percent. The values n cos and n 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λ eff . Therefore, we evaluate the exactΦ by the Formula (15) for n n 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λ eff . For our calculations, we choose four different values, that areλ eff = 0.01, 0.1, 1 and 5. The rescaled screening lengthλ eff is the ratio of the physical effective screening length λ eff to the physical size of the period al. As we have mentioned previously, today λ eff ∼ 2.6 Gpc for the Λ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 ∼ 27 Gpc [13], i.e., the observational data require thatλ eff 1. Nevertheless, taking into account that many N-body simulations are indeed performed in boxes with sizes less than 1 Gpc, we also considerλ eff ≥ 1.
The Yukawa-Ewald potential (15) is sensitive to the free parameter α as well. Therefore, choosing different values of this quantity, we also seek those at which n mix will be minimum. Our calculations demonstrate that for the chosen range of 0.01 ≤λ eff ≤ 5, the optimal value of α 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λ eff = 0.01 andλ eff = 0.1, respectively. As for the values n mix , these numbers are the same as n exp for both left and right charts as long as α lies between 10 −3 and 2. However, outside this interval, the Yukawa-Ewald formula may require more summands. For example, ifλ eff = 0.1 and α = 2.5, then n mix = 29, 23, 23, 4, 1, 48, 28, 23, 6 for the points A 1 , A 2 , . . . , C 4 , respectively. All in all, whenλ eff 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λ eff .
In Table 2, we present the results of similar calculations forλ eff = 1 andλ eff = 5. In the caseλ eff = 1 (top chart), n exp n mix for all selected points, and the optimal choice for the parameter α in the Yukawa-Ewald formula (15) is 2. Forλ 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 α = 2 is again the optimal choice for the Yukawa-Ewald potential. Hence, whenλ eff ≥ 1, the Yukawa-Ewald formula (15) delivers the best performance in numerical calculations.
We demonstrate in Figures 1 and 2

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 where 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λ 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 mix are exactly the same (assuming here and in what follows that α = 2). In the caseλ eff = 0.1, we get n exp = 1, 1, 12, 1 for these points, and the corresponding numbers n mix are, again, identical. However, forλ eff = 1, we have n exp = 258, 82, 987, 486 while n mix = 7, 7, 21, 9. Finally, ifλ eff = 5, then still n mix = 7, 7, 21, 9, but n exp acquires unreasonably large values.
Our calculations demonstrate that the numbers n exp start to grow onceλ eff exceeds 0.1 and they acquire large values asλ eff approaches 1, while n mix remain small throughout (provided that the value of the parameter α is optimal, e.g., for α = 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 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 mix .

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 × T × 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λ eff = λ eff /(al), we analyzed both small values and values equal to or greater than 1: 0.01 ≤λ eff ≤ 5. The Yukawa-Ewald formula additionally admits a free parameter α, and we have revealed that for the given range ofλ eff , the optimal value of α 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λ 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 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λ eff > 0.1: n exp begins to increase quickly while n mix still takes on rather small values. Therefore, for suchλ 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 λ 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λ 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.