Mean Square Geodesic Deviation in the Zeldovich Problem on Light Propagation in a Universe with Inhomogeneities

In 1964, Ya.B. Zeldovich formulated the problem of light propagation in the Universe under the influence of inhomogeneities. It is reduced to describing the divergence of two close geodesics in a Riemannian space and is described by the geodesic deviation equation (Jacobi equation) with the curvature along the geodesic line varying randomly. Assuming the curvature to be constant on segments of small but finite length, the problem is reduced to studying the product of random matrices and makes it possible to apply the appropriate well-developed mathematical theory, which, however, did not allow calculating the mean-square growth rate of the geodesic deviation. In our paper, we propose a way to solve this problem by introducing a bilinear quantity, one component of which coincides with the square of the Jacobi field. The system of first-order differential equations for the bilinear quantity is explicitly written out, and the solution, same as the growth rate, is again expressed through the product of matrices. Such a technique can be used in the study of a wide range of problems and is naturally generalized to higher-order moments.


INTRODUCTION
Back in 1964, Zeldovich drew attention to the fact that light in the Universe propagates in such a way that an observer measuring angular distances between objects by means of standard cosmological tests obtains a curvature of space slightly less than the curvature corresponding to the average density of the Universe [1]. In particular, if the average density of the world is exactly equal to the critical density, the observer should come to the conclusion that the Universe is open. At a critical density, the spatial section is flat; at a lower density, it has a negative curvature. In a space of zero curvature, the distance between close geodesics increases linearly with the lengths of geodesics, while in a space of negative curvature it grows exponentially. This phenomenon arises from the collective effect of small spatial curvature inhomogeneities. The observer should notice this phenomenon, if the observation conditions allow the difference between linear and exponential growth to be detected. The physical effect itself is fully acknowledged by modern cosmology and is mentioned in various contexts in studies on gravitational lensing; however, quantitatively, it is rather small due to the fact that the density fluctuations are minor and the density is close to critical [2].
The effect noted by Zeldovich is geometric in nature and is reduced to the problem of geodesic divergence in a Riemannian space, the curvature of which along the geodesic can be considered a random process (see [3] for more details). The nature of this effect is not related to the four-dimensionality of space-time and the presence of the time coordinate. Since we do not plan to examine direct cosmological consequences of Zeldovich's concepts (they are exhausted in his paper) in this study, we are talking simply about the geodesic divergence on a curved twodimensional Riemannian space. It should be noted that the original work was written by Zeldovich before the formation of modern physics of random media, so he solved the problem using a very specific technique. For its solution by modern regular methods, see [3].
The growth rate of the distance between geodesic lines (geodesic deviation) can be defined in various terms. In random media, we can talk about the Lyapunov exponent (selective growth rate) and the growth rates of the normalized statistical moments of geodesic deviation. In his paper, Zeldovich suggested that these growth rates should not coincide; the growth rate of the statistical moments should exceed the Lyapunov exponent, and the higher moments should increase faster than the lower ones.
Over time, it was realized that Zeldovich's problem is a convenient model example for studying the development of various instabilities in a random media, and the ratio of the growth rates of various characteristics of geodesic deviation noted by Zeldovich is a manifes-tation of the general property of the instability development in random media, which eventually became known as intermittency [4].
In his study, Zeldovich arrived at correct answers using reasoning that did not require explicit calculation of the growth rates. Under certain assumptions about curvature properties as a random process (see below for more details), it is possible to calculate the Lyapunov exponent using the mathematical theory of the product of random matrices; however, these methods do not allow calculating the mean-square growth rate and the growth rate of the highest moments [5]. It was possible only in the model of curvature fluctuations in the form of white noise [6]. This model of specifying a random curvature-defining process is poorly compatible with the concept of the properties of gravitational forces, and it is desirable to be replaced with more physical assumptions. This is the content of the present study.
In his study, Zeldovich expressed a number of fruitful thoughts, which are only partly related to the idea that is of interest to us. They were developed in subsequent years by a number of authors, including Zeldovich himself (see, e.g., [7][8][9][10]). Although a detailed review of the modern scientific perception of Zeldovich's ideas would certainly be of interest, this is not the aim of our study.

BASIC EQUATIONS
One of the important technical discoveries in Zeldovich's work was the realization that the phenomenon under study can be reduced to the behavior of geodesics in a two-dimensional Riemannian space spanned by the projection of two close light rays onto a spatial section of the cosmological model, and the results should be rescaled for a non-stationary multidimensional cosmological model.
Let us consider two geodesics in a two-dimensional Riemannian space, intersecting at one point; the angle between the geodesics at the intersection point is considered a small value. Let us mark the distance (in the metric of a Riemannian manifold) on both geodesics in the same direction from the point of their intersection. The distance (in the same metric) between the obtained points in the first approximation with respect to is , where the value is known in Riemannian geometry as the Jacobi field (in physics is known as the geodesic deviation). It turns out (see, e.g., [11]) that the Jacobi field satisfies the equation (1) which is called the Jacobi equation. Here, the derivatives are taken with respect to the variable , and is the Gaussian curvature of the manifold, the only nonzero component of the four-dimensional Riemann tensor, which is considered a random process. Assum-θ x θ θy y ing that the points move along the geodesics at a constant speed, the distance can also be understood as the time over which the points travel a given distance; thus, we are dealing with the evolution of the value in a random medium . We assume that the random process is arranged as follows. On each segment of the form (these are called update intervals), the curvature of does not depend on and is considered a random variable with zero mean value and finite variance ; on different segments, these random variables are statistically independent and equally distributed (for example, according to the Gaussian law). This model of a stochastic process is called an update model. Other assumptions about the structure of the random process are also possible, and they can be studied by such methods as well (see, e.g., [12,13]). The proposed method does not require a special type of statistical distribution of curvature. In order to obtain a specific result, it is necessary to set some specific distribution, as shown below.
To solve Eq. (1), it is necessary to set the initial conditions. We will assume that , .
Let us rewrite Eq. (1) as a system of first-order equations by introducing a two-dimensional vector , the first coordinate of which is and the second coordinate is (the constant factor equal to the length of the update interval is introduced to match the dimensions of the vector components). Then, in the matrix form Naturally, vector is also a (vector) random process.

EQUATION FOR THE CORRELATION TENSOR
Let us now construct the correlation tensor for the random process . To do this, we introduce the twoindex tensor (2) whose mean value (calculated from the distribution) is the correlation tensor. Differentiating the product and using (2), it is easy to show that satisfies the equation (3) where the tensor components are At each update interval, Eq. (3) is an equation with constant coefficients, so its solution is the value of the tensor at the point multiplied by the exponent from tensor , which we denote by , since it is, of course, a fourth-rank tensor as well. We have to supplement this tensor with the subscript , which indicates the update interval at which the calculation was performed. Naturally, the specific values of the tensor depend on the interval number and are random. Now, we can average Eq. (3). First, let us do this for the first interval . It is sufficient to calculate the mean value of the tensor , which is calculated component by component.
Since the values are replaced by independent values at the ends of the update intervals, we can carry out the same calculation sequentially at the next update intervals; the mean value of the tensor does not depend on since the distributions are the same at different update intervals.
For the convenience of further reasoning, we will renumber the variables. Let us introduce an auxiliary four-dimensional quantity with components , , , and , respectively. The quantity satisfies the equation where The solution to this equation at each update interval is expressed through the matrix . Therefore, the evolution of the averaged value is expressed through the matrix . Moving from one update interval to another, approaches the eigenvector of the matrix corresponding to the highest eigenvalue, which, in turn, determines the growth rate of the components. Since the mean value of the component is equal to the mean value of , the highest eigenvalue of the matrix also determines the growth rate of the second moment of the Jacobi field.

CALCULATION OF THE EIGENVALUE
First, let us find the explicit form of the matrix using the definition of the matrix exponent. Due to the simple structure of the matrix itself, its degrees are written out explicitly. The form of the matrix depends on the sign of . For convenience, we temporarily introduce the dimensionless curvature . Then at , we have (4) and at k < 0 the matrix looks as Note that it would be sufficient to write out only one matrix, for example, for ; the form of the second matrix for negative would then follow from the relations and .
Let us further show how in the case of small , we can obtain approximate estimates for the eigenvalues of the matrix . The exact answer will depend on the distribution of the value. Our goal will be to obtain a series expansion of the eigenvalues in powers of .
Let us start with a simple case in which the matrix is calculated explicitly. To do this, we assume that takes only two values, , with equal probability. Then, is the half-sum of matrices (4) and (5), and an explicit characteristic equation can be written out.
The coefficients consist of combinations of and cosh . We will write out the result of their expansion in powers of . The approximate characteristic equation accurate to the terms of the order of has the form (6) Below, we will explain why it is more correct to view the resulting expansion as a series in powers of and not . For now, we will deal with the solution of Eq. (6). It is easy to see that one of the roots is . The remaining cubic equation has the form (7) Let , then (7) will take the form A standard replacement of variables λ = y + allows us to bring the cubic equation to the canonical form (note that all the transformations are now performed with accuracy to ): Another standard replacement , known as Vieta substitution, with accuracy up to leads to the equation     The roots of the last equation are expressed easily; considering all the replacements, we can write out the solutions of (7): (8) where . The highest eigenvalue is what determines the growth of the solutions of the extended system (3).
Note that it was possible to arrive at Eq. (6) in a slightly different way, which is applicable for a wider class of distributions. The idea is to replace the mathematical expectation of a function of a random variable with the mathematical expectation of the first few terms of its Taylor series. The question of the accuracy of such an approximation is rather complicated, since it is necessary to consider not only the local behavior of the function near the expansion point, but also the character of the growth of the random variable moments. In our situation, the issue is simplified, since we consider the product of a random variable and a small parameter . Then, with accuracy to the terms of the order of for an arbitrary sufficiently smooth function in a neighborhood of zero and symmetric distribution, we can write (9) Using this approximation to calculate the elements of the matrix and considering the fact that the first moment for a symmetric distribution is zero, we arrive exactly at expression (6), in which σ 2 = DK = . Note that the following refining corrections in (9) contain the moments of the random variable of higher orders, which are generally not expressed in terms of its variance . For this reason, (6) is generally not an expansion in powers of .
Thus, only in the first approximation for a symmetric distribution, the highest eigenvalue of the matrix is determined by only two parameters, the variance of the distribution and the length of the update interval .
It remains to explicitly recalculate the highest eigenvalue in terms of the growth rate of the second statistical moment of the Jacobi field (1).

( ) ( ) ( )
, λ = + δ + δ + δ , λ = − δ + δ + δ ± δ − δ , No. 5 2021 5. RELATIONSHIP OF THE EIGENVALUE AND GROWTH RATE Recall (see [4]) that the growth rate of the th statistical moment is the value (10) where is the number of the statistical moment, and the growth rate is normalized to the moment number in the denominator of Eq. (10) to make the results comparable for moments of different orders.
The highest eigenvalue found in (8) characterizes the growth of the mean square of the Jacobi field in construction. Specifically, this means that the value increases by times at a distance of update intervals. On the other hand, according to the definition of the growth rate (10), the same change will be written as . Hence, the relation between the eigenvalue and the growth rate of the second moment: Substituting the expression for from (8), with accuracy to terms of the order of , we obtain (11) Note that proceeding to expansion (11), it is necessary to consider the limited region of convergence of the series for the logarithm, hence, the restrictions on the possible values of the product . The numerical estimates show that the condition should be fulfilled for the logarithm expansion formula to be valid.
Of course, we assume that the highest eigenvalue of our matrix and the mean-square growth rate for more complex distributions are calculated numerically.
To demonstrate the agreement of the estimates with the numerical experiment, we consider three families of distributions: Bernoulli with values , uniform on the interval , , and Gaussian with zero mean and variance . The distribution parameters are selected so that the variance of is exactly . Furthermore, we will consider , 1, and 10, and the values running the interval from to 10. For each and , we will simulate realizations of the random variable , numerically estimate the matrix , and find its highest eigenvalue. The results will be plotted on the graphs of the highest eigenvalue versus the value and compared with the σΔ 2 σΔ < . (8). Figure 1 shows that the results for small are nearly identical and closely agree with the approximation; for large values, a discrepancy appears, since the higher moments of the distribution start to prevail.
Finally, Fig. 2 shows the dependence of the growth rate on the length of the update interval for different variances of the distribution. Note that at small , the growth rate is determined by the first term . The value at which the expansion of the logarithm loses accuracy approximately coincides with the value at which the approximation accuracy for the eigenvalues is lost (see Fig. 1). Apparently, a good approximation for small values of the product is . Note one particular case when and is uniformly distributed on the interval considered earlier in [5]. For it, we obtain the approximate value of the mathematical expectation of the matrix , and the highest eigenvalue (formula (8) also gives a close value ). In terms of the growth rate, we obtain . Note that the value found in [5] for the Lyapunov exponent turns out, as it should, to be smaller than the growth rate of the highest moment . It should also be noted that the direct application of formula (11) gives an overestimated value . This is due to the fact that this case lies outside the limits of applicability regarding the expansion of the logarithm.

RESULTS AND DISCUSSION
We have demonstrated how to calculate the meansquare growth rate of the geodesic deviation in the Zeldovich problem. The difficulty of the calculation consists in the fact that we are interested in a quadratic quantity (the square of the geodesic deviation), for which it was previously not possible to obtain a linear equation. We solve this issue by introducing a bilinear quantity (second-order tensor), one component of which is the square of the geodesic deviation. This technique is widely used in the theory of turbulence. For example, in the study of the magnetic energy evolution in a flow of a conductive fluid, it is unclear how to obtain an equation for this quadratic quantity, but it is possible to construct an equation for the correlation tensor of the magnetic field [14].
The technique we use is in no way limited to the specific form of the geodesic deviation equation, but is .    applicable to the transfer of an arbitrary physical quantity described by a system of ordinary linear differential equations with random coefficients. Many problems that study the development of instabilities in random media by means of a transition to the Lagrangian reference frame are reduced to such equations (see, e.g., [4]).
With the help of this technique, it is also possible to study the growth of statistical moments of a higher order. For example, to study the growth rate of the fourth degree of geodesic deviation, it will be necessary to introduce the fourth-order tensor . Naturally, this will require consideration of matrices of a higher order and will make the calculations more cumbersome, but will not fundamentally change their nature.
With regard to the prospects for the application of the results and the entire body of knowledge about the Zeldovich effect, we note once again that at the stages of the evolution of the Universe close to us in time, this effect is quantitatively very small. In principle, it is possible that the combined action of multiple celestial bodies located near one line of sight, which individually do not lead to gravitational lensing, can collectively cause the occurrence of a gravitational lens. To estimate the probability of the occurrence of such a lens, one must be able to calculate various moments of the Jacobi field. However, this line of research does not appear particularly promising. The Zeldovich effect seems much more important for the study of the earliest Universe, since it clearly says that the influ-= ijkl i j k l Z zz z z ence of inhomogeneities blurs the concept of critical density for the cosmological model and makes one think how to correctly formulate the Einstein equations for fluctuating space-time. It seems that the problem formulated by Zeldovich successfully isolates a fragment of this large problem; this fragment can be researched without consistent mathematical elaboration of the concept of a pseudo-Riemannian manifold with a random metric and, at the same time, indicates a possible direction of the fundamental development of general relativity. The authors hope to make a feasible contribution to this development in the future.   unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.