# BKL singularity

Figure 1. A spherical body undergoing a chaotic BKL (Mixmaster) dynamics close to singularity according to rules eq. 35 with initial ${\displaystyle u={\sqrt {13}}}$

A Belinskii-Khalatnikov-Lifshitz (BKL) singularity is a model of the dynamic evolution of the Universe near the initial singularity, described by an anisotropic, chaotic solutions of the Einstein field equations of gravitation.[1] According to this model, the Universe is chaotically oscillating around a gravitational singularity in which time and space become equal to zero. This singularity is physically real in the sense that it is a necessary property of the solution, and will appear also in the exact solution of those equations. The singularity is not artificially created by the assumptions and simplifications made by the other special solutions such as the Friedmann–Lemaître–Robertson–Walker, quasi-isotropic, and Kasner solutions.

The picture developed by BKL has several important elements. These are:

• Near the singularity the evolution of the geometry at different spatial points decouples so that the solutions of the partial differential equations can be approximated by solutions of ordinary differential equations with respect to time for appropriately defined spatial scale factors. This is called the BKL conjecture.
• For most types of matter the effect of the matter fields on the dynamics of the geometry becomes negligible near the singularity. Or, in the words of John Wheeler, "matter doesn't matter" near a singularity. The original BKL work posed a negligible effect for all matter but later they theorized that "stiff matter" (equation of state p = ε) equivalent to a massless scalar field can have a modifying effect on the dynamics near the singularity.
• The ordinary differential equations which describe the asymptotics are those which come from a class of spatially homogeneous solutions which constitute the Mixmaster dynamics: a complicated oscillatory and chaotic model that exhibits properties similar to those discussed by BKL.

The study of the dynamics of the universe in the vicinity of the cosmological singularity has become a rapidly developing field of modern theoretical and mathematical physics. The generalization of the BKL model to the cosmological singularity in multidimensional (Kaluza–Klein type) cosmological models has a chaotic character in the spacetimes whose dimensionality is not higher than ten, while in the spacetimes of higher dimensionalities a universe after undergoing a finite number of oscillations enters into monotonic Kasner-type contracting regime.[2][3][4]

The development of cosmological studies based on superstring models has revealed some new aspects of the dynamics in the vicinity of the singularity.[5][6][7] In these models, mechanisms of changing of Kasner epochs are provoked not by the gravitational interactions but by the influence of other fields present. It was proved that the cosmological models based on six main superstring models plus D = 11 supergravity model exhibit the chaotic BKL dynamics towards the singularity. A connection was discovered between oscillatory BKL-like cosmological models and a special subclass of infinite-dimensional Lie algebras – the so called hyperbolic Kac–Moody algebras.[8]

## Introduction

The basis of modern cosmology are the special solutions of the Einstein field equations found by Alexander Friedmann in 1922–1924. The Universe is assumed homogeneous (space has the same metric properties (measures) in all points) and isotropic (space has the same measures in all directions). Friedmann's solutions allow two possible geometries for space: closed model with a ball-like, outwards-bowed space (positive curvature) and open model with a saddle-like, inwards-bowed space (negative curvature). In both models, the Universe is not standing still, it is constantly either expanding (becoming larger) or contracting (shrinking, becoming smaller). This was confirmed by Edwin Hubble who established the Hubble redshift of receding galaxies. The present consensus is that the isotropic model, in general, gives an adequate description of the present state of the Universe; however, isotropy of the present Universe by itself is not a reason to expect that it is adequate for describing the early stages of Universe evolution. At the same time, it is obvious that in the real world homogeneity is, at best, only an approximation. Even if one can speak about a homogeneous distribution of matter density at distances that are large compared to the intergalactic space, this homogeneity vanishes at smaller scales. On the other hand, the homogeneity assumption goes very far in a mathematical aspect: it makes the solution highly symmetric which can impart specific properties that disappear when considering a more general case.

Another important property of the isotropic model is the inevitable existence of a time singularity: time flow is not continuous, but stops or reverses after time reaches some (very large or very small) value. Between singularities, time flows in one direction: away from the singularity (arrow of time). In the open model, there is one time singularity so time is limited at one end but unlimited at the other, while in the closed model there are two singularities that limit time at both ends (the Big Bang and Big Crunch).

The only physically interesting properties of spacetimes (such as singularities) are those which are stable, i.e., those properties which still occur when the initial data is perturbed slightly. It is possible for a singularity to be stable and yet be of no physical interest: stability is a necessary but not a sufficient condition for physical relevance. For example, a singularity could be stable only in a neighbourhood of initial data sets corresponding to highly anisotropic universes. Since the actual universe is now apparently almost isotropic such a singularity could not occur in our universe. A sufficient condition for a stable singularity to be of physical interest is the requirement that the singularity be generic (or general). Roughly speaking, a stable singularity is generic if it occurs near every set of initial conditions and the non-gravitational fields are restricted in some specified way to "physically realistic" fields so that the Einstein equations, various equations of state, etc., are assumed to hold on the evolved spacetimes. It might happen that a singularity is stable under small variations of the true gravitational degrees of freedom, and yet it is not generic because the singularity depends in some way on the coordinate system, or rather on the choice of the initial hypersurface from which the spacetime is evolved.

For a system of non-linear differential equations, such as the Einstein equations, a general solution is not unambiguously defined. In principle, there may be multiple general integrals, and each of those may contain only a finite subset of all possible initial conditions. Each of those integrals may contain all required independent functions which, however, may be subject to some conditions (e.g., some inequalities). Existence of a general solution with a singularity, therefore, does not preclude the existence of other additional general solutions that do not contain a singularity. For example, there is no reason to doubt the existence of a general solution without a singularity that describes an isolated body with a relatively small mass.

It is impossible to find a general integral for all space and for all time. However, this is not necessary for resolving the problem: it is sufficient to study the solution near the singularity. This would also resolve another aspect of the problem: the characteristics of spacetime metric evolution in the general solution when it reaches the physical singularity, understood as a point where matter density and invariants of the Riemann curvature tensor become infinite.

## Existence of physical time singularity

One of the principal problems studied by the Landau group (to which BKL belong) was whether relativistic cosmological models necessarily contain a time singularity or whether the time singularity is an artifact of the assumptions used to simplify these models. The independence of the singularity on symmetry assumptions would mean that time singularities exist not only in the special, but also in the general solutions of the Einstein equations. It is reasonable to suggest that if a singularity is present in the general solution, there must be some indications that are based only on the most general properties of the Einstein equations, although those indications by themselves might be insufficient for characterizing the singularity.

A criterion for generality of solutions is the number of independent space coordinate functions that they contain. These include only the "physically independent" functions whose number cannot be reduced by any choice of reference frame. In the general solution, the number of such functions must be enough to fully define the initial conditions (distribution and movement of matter, distribution of gravitational field) at some moment of time chosen as initial. This number is four for an empty (vacuum) space, and eight for a matter and/or radiation-filled space.[9][10]

Previous work by the Landau group[11][12][13] (reviewed in[9] ) led to the conclusion that the general solution does not contain a physical singularity. This search for a broader class of solutions with a singularity has been done, essentially, by a trial-and-error method, since a systematic approach to the study of the Einstein equations was lacking. A negative result, obtained in this way, is not convincing by itself; a solution with the necessary degree of generality would invalidate it, and at the same time would confirm any positive results related to the specific solution.

At that time, the only known indication was related to the form of the Einstein equations written in a synchronous frame, that is, in a frame in which the proper time x0 = t is synchronized throughout the whole space; in this frame the space distance element dl is separate from the time interval dt.[14] The Einstein equation

${\displaystyle R_{0}^{0}=T_{0}^{0}-{\tfrac {1}{2}}T}$

(eq. 1)

written in synchronous frame gives a result in which the metric determinant g inevitably becomes zero in a finite time irrespective of any assumptions about matter distribution.[9][10]

This indication, however, was dropped after it became clear that it is linked with a specific geometric property of the synchronous frame: the crossing of time line coordinates. This crossing takes place on some encircling hypersurfaces which are four-dimensional analogs of the caustic surfaces in geometrical optics; g becomes zero exactly at this crossing.[13] Therefore, although this singularity is general, it is fictitious, and not a physical one; it disappears when the reference frame is changed. This, apparently, removed the incentive among the researchers for further investigations along these lines.

However, the interest in this problem waxed again in the 1960s after Penrose published his theorems[15] that linked the existence of a singularity of unknown character with some very general assumptions that did not have anything in common with a choice of reference frame. Other similar theorems were found later on by Hawking[16][17] and Geroch[18] (see Penrose–Hawking singularity theorems). This revived interest in the search for singular solutions.

## Generalized homogeneous solution

In a space that is both homogeneous and isotropic the metric is determined completely leaving free only the sign of the curvature. Assuming only space homogeneity with no additional symmetry such as isotropy leaves considerably more freedom in choosing the metric. The following pertains to the space part of the metric at a given instant of time t assuming a synchronous spacetime reference system so that t is the same synchronized time for the whole space.

Homogeneity implies identical metric properties at all points of the space. An exact definition of this concept involves considering sets of coordinate transformations that transform the space into itself, i.e. leave its metric unchanged: if the line element before transformation is

${\displaystyle dl^{2}=\gamma _{\alpha \beta }\left(x^{1},x^{2},x^{3}\right)dx^{\alpha }dx^{\beta },}$

then after transformation the same line element is

${\displaystyle dl^{2}=\gamma _{\alpha \beta }\left(x^{\prime 1},x^{\prime 2},x^{\prime 3}\right)dx^{\prime \alpha }dx^{\prime \beta },}$

with the same functional dependence of γαβ on the new coordinates. (For a more theoretical and coordinate-independent definition of homogeneous space see homogeneous space). A space is homogeneous if it admits a set of transformations (a group of motions) that brings any given point to the position of any other point. Since space is three-dimensional the different transformations of the group are labelled by three independent parameters.

Figure 2. The triad e(a) (e(1), e(2), e(3)) is an affine coordinate system (including as a special case Cartesian coordinate system) whose coordinates are functions of the curvilinear coordinates xα (x1, x2, x3).

In Euclidean space the homogeneity of space is expressed by the invariance of the metric under parallel displacements (translations) of the Cartesian coordinate system. Each translation is determined by three parameters — the components of the displacement vector of the coordinate origin. All these transformations leave invariant the three independent differentials (dx, dy, dz) from which the line element is constructed. In the general case of a non-Euclidean homogeneous space, the transformations of its group of motions again leave invariant three independent linear differential forms, which do not, however, reduce to total differentials of any coordinate functions. These forms are written as ${\displaystyle e_{\alpha }^{(a)}dx^{\alpha }}$  where the Latin index (a) labels three independent vectors (coordinate functions); these vectors are called a frame field or triad. The Greek letters label the three space-like curvilinear coordinates. A spatial metric invariant is constructed under the given group of motions with the use of the above forms:

${\displaystyle dl^{2}=\eta _{ab}\left(e_{\alpha }^{(a)}dx^{\alpha }\right)\left(e_{\beta }^{(b)}dx^{\beta }\right)}$

(eq. 6a)

i.e. the metric tensor is

${\displaystyle \gamma _{\alpha \beta }=\eta _{ab}e_{\alpha }^{(a)}e_{\beta }^{(b)}}$

(eq. 6b)

where the coefficients ηab, which are symmetric in the indices a and b, are functions of time. The choice of basis vectors is dictated by the symmetry properties of the space and, in general, these basis vectors are not orthogonal (so that the matrix ηab is not diagonal).

The reciprocal triple of vectors ${\displaystyle e_{(a)}^{\alpha }}$  is introduced with the help of Kronecker delta

${\displaystyle e_{(a)}^{\alpha }e_{\alpha }^{(b)}=\delta _{a}^{b};\quad e_{(a)}^{\alpha }e_{\beta }^{(a)}=\delta _{\alpha }^{\beta }}$

(eq. 6c)

In the three-dimensional case, the relation between the two vector triples can be written explicitly

${\displaystyle \mathbf {e} _{(1)}={\frac {1}{v}}\mathbf {e} ^{(2)}\times \mathbf {e} ^{(3)},\quad \mathbf {e} _{(2)}={\frac {1}{v}}\mathbf {e} ^{(3)}\times \mathbf {e} ^{(1)},\quad \mathbf {e} _{(3)}={\frac {1}{v}}\mathbf {e} ^{(1)}\times \mathbf {e} ^{(2)},}$

(eq. 6d)

where the volume v is

${\displaystyle v=\left\vert e_{\alpha }^{(a)}\right\vert =\mathbf {e} ^{(1)}\cdot \mathbf {e} ^{(2)}\times \mathbf {e} ^{(3)},}$

with e(a) and e(a) regarded as Cartesian vectors with components ${\displaystyle e_{(a)}^{\alpha }}$  and ${\displaystyle e_{\alpha }^{(a)}}$ , respectively. The determinant of the metric tensor eq. 6b is γ = ηv2 where η is the determinant of the matrix ηab.

The required conditions for the homogeneity of the space are

${\displaystyle \left({\frac {\partial e_{\alpha }^{(c)}}{\partial x^{\beta }}}-{\frac {\partial e_{\beta }^{(c)}}{\partial x^{\alpha }}}\right)e_{(a)}^{\alpha }e_{(b)}^{\beta }=C_{ab}^{c}.}$

(eq. 6e)

The constants ${\displaystyle C_{ab}^{c}}$  are called the structure constants of the group.

Multiplying by ${\displaystyle e_{(c)}^{\gamma }}$ , eq. 6e can be rewritten in the form

${\displaystyle e_{(a)}^{\alpha }{\frac {\partial e_{(b)}^{\gamma }}{\partial x^{\alpha }}}-e_{(b)}^{\beta }{\frac {\partial e_{(a)}^{\gamma }}{\partial x^{\beta }}}=C_{ab}^{c}e_{(c)}^{\gamma }.}$

(eq. 6f)

Equation 6e can be written in a vector form as

${\displaystyle \mathbf {e} _{(a)}\times \mathbf {e} _{(b)}{\text{curl }}\mathbf {e} ^{(c)}=-C_{ab}^{c},}$

where again the vector operations are done as if the coordinates xα were Cartesian. Using eq. 6d, one obtains

${\displaystyle {\frac {1}{v}}\left(\mathbf {e} ^{(1)}{\text{curl }}\mathbf {e} ^{(1)}\right)=C_{32}^{1},\quad {\frac {1}{v}}\left(\mathbf {e} ^{(2)}{\text{curl }}\mathbf {e} ^{(1)}\right)=C_{13}^{1},\quad {\frac {1}{v}}\left(\mathbf {e} ^{(3)}{\text{curl }}\mathbf {e} ^{(1)}\right)=C_{21}^{1}}$

(eq. 6g)

and six more equations obtained by a cyclic permutation of indices 1, 2, 3.

The structure constants are antisymmetric in their lower indices as seen from their definition eq. 6e: ${\displaystyle C_{ab}^{c}=-C_{ba}^{c}}$ . Another condition on the structure constants can be obtained by noting that eq. 6f can be written in the form of commutation relations

${\displaystyle \left[X_{a},X_{b}\right]\equiv X_{a}X_{b}-X_{b}X_{a}=C_{ab}^{c}X_{c}}$

(eq. 6h)

for the linear differential operators

${\displaystyle X_{a}=e_{(a)}^{\alpha }{\frac {\partial }{\partial x^{\alpha }}}}$

(eq. 6i)

In the mathematical theory of continuous groups (Lie groups) the operators Xa satisfying conditions of the form eq. 6h are called the generators of the group. However, to avoid confusion when comparing with other presentations, it should be mentioned that the systematic theory usually starts from operators defined using the Killing vectors ${\displaystyle \xi _{(a)}^{\alpha }}$  (since in the synchronous metric none of the γαβ components depends on time, the Killing vectors are time-like):

${\displaystyle X_{a}=\xi _{(a)}^{\alpha }{\frac {\partial }{\partial x^{\alpha }}}.}$

The condition mentioned above follows from the Jacobi identity

${\displaystyle [[X_{a},X_{b}],X_{c}]+[[X_{b},X_{c}],X_{a}]+[[X_{c},X_{a}],X_{b}]=0}$

and has the form

${\displaystyle C_{ab}^{e}C_{ec}^{d}+C_{bc}^{e}C_{ea}^{d}+C_{ca}^{e}C_{eb}^{d}=0}$

(eq. 6j)

It is a definite advantage to use, in place of the three-index constants ${\displaystyle C_{ab}^{c}}$ , a set of two-index quantities, obtained by the dual transformation

${\displaystyle C_{ab}^{c}=e_{abd}C^{dc}}$

(eq. 6k)

where eabc = eabc is the unit antisymmetric symbol (with e123 = +1). With these constants the commutation relations eq. 6h are written as

${\displaystyle e^{abc}X_{b}X_{c}=C^{ad}X_{d}}$

(eq. 6l)

The antisymmetry property is already taken into account in the definition eq. 6k, while property eq. 6j takes the form

${\displaystyle e_{bcd}C^{cd}C^{ba}=0}$

(eq. 6m)

Figure 3. The parameter space as a 3-plane (class A) and an orthogonal half 3-plane (class B) in R4 with coordinates (n(1), n(2), n(3), a), showing the canonical representatives of each Bianchi type.

The choice of the three frame vectors in the differential forms ${\displaystyle e_{\alpha }^{(a)}dx^{\alpha }}$  (and with them the operators Xa) is not unique. They can be subjected to any linear transformation with constant coefficients:

${\displaystyle \mathbf {e} _{(a)}=A_{a}^{b}\mathbf {e} _{(b)}.}$

(eq. 6n)

The quantities ηab and Cab behave like tensors (are invariant) with respect to such transformations.

The conditions eq. 6m are the only ones that the structure constants must satisfy. But among the constants admissible by these conditions, there are equivalent sets, in the sense that their difference is related to a transformation of the type eq. 6n. The question of the classification of homogeneous spaces reduces to determining all nonequivalent sets of structure constants. This can be done, using the "tensor" properties of the quantities Cab, by the following simple method (C. G. Behr, 1962).

The unsymmetric "tensor" Cab can be resolved into a symmetric and an antisymmetric part. The first is denoted by nab, and the second is expressed in terms of its "dual vector" ac:

${\displaystyle C^{ab}=n^{ab}+e^{abc}a_{c}.}$

(eq. 6o)

Substitution of this expression in eq. 6m leads to the condition

${\displaystyle n^{ab}a_{b}=0.}$

(eq. 6p)

By means of the transformations eq. 6m the symmetric "tensor" nab can be brought to diagonal form with eigenvalues n1, n2, n3. Equation 6p shows that the "vector" ab (if it exists) lies along one of the principal directions of the "tensor" nab, the one corresponding to the eigenvalue zero. Without loss of generality one can therefore set ab = (a, 0, 0). Then eq. 6p reduces to an1 = 0, i.e. one of the quantities a or n1 must be zero. The Jacobi identities take the form:

${\displaystyle [X_{1},X_{2}]=-aX_{2}+n_{3}X_{3},\quad [X_{2},X_{3}]=n_{1}X_{1},\quad [X_{3},X_{1}]=n_{2}X_{2}+aX_{3}.}$

(eq. 6q)

The only remaining freedom is a change of sign of the operators Xa and arbitrary scale transformations of them (multiplication by constants). This permits us simultaneously to change the sign of all the na and also to make the quantity a positive (if it is different from zero). Also all structure constants can be made equal to ±1, if at least one of the quantities a, n2, n3 vanishes. But if all three of these quantities differ from zero, the scale transformations leave invariant the ratio h = a2(n2n3)−1.

Thus one arrives at the Bianchi classification listing the possible types of homogeneous spaces classified by the values of a, n1, n2, n3 which is graphically presented in Fig. 3. In the class A case (a = 0), type IX (n(1)=1, n(2)=1, n(3)=1) is represented by octant 2, type VIII (n(1)=1, n(2)=1, n(3)=–1) is represented by octant 6, while type VII0 (n(1)=1, n(2)=1, n(3)=0) is represented by the first quadrant of the horizontal plane and type VI0 (n(1)=1, n(2)=–1, n(3)=0) is represented by the fourth quadrant of this plane; type II ((n(1)=1, n(2)=0, n(3)=0) is represented by the interval [0,1] along n(1) and type I (n(1)=0, n(2)=0, n(3)=0) is at the origin. Similarly in the class B case (with n(3) = 0), Bianchi type VIh (a=h, n(1)=1, n(2)=–1) projects to the fourth quadrant of the horizontal plane and type VIIh (a=h, n(1)=1, n(2)=1) projects to the first quadrant of the horizontal plane; these last two types are a single isomorphism class corresponding to a constant value surface of the function h = a2(n(1)n(2))−1. A typical such surface is illustrated in one octant, the angle θ given by tanθ = |h/2|1/2; those in the remaining octants are obtained by rotation through multiples of π/2, h alternating in sign for a given magnitude |h|. Type III is a subtype of VIh with a=1. Type V (a=1, n(1)=0, n(2)=0) is the interval (0,1] along the axis a and type IV (a=1, n(1)=1, n(2)=0) is the vertical open face between the first and fourth quadrants of the a = 0 plane with the latter giving the class A limit of each type.

### The BKL conjecture

In their 1970 work,[1] BKL stated that as one approaches a singularity, terms containing time derivatives in Einstein’s equations dominate over those containing spatial derivatives. This has since been known as the BKL conjecture and implies that Einstein’s partial differential equations (PDE) are well approximated by ordinary differential equations (ODEs), whence the dynamics of general relativity effectively become local and oscillatory. The time evolution of fields at each spatial point is well approximated by the homogeneous cosmologies in the Bianchi classification.

By separating the time and space derivatives in the Einstein equations, for example, in the way used above for the classification of homogeneous spaces, and then setting the terms containing space derivatives equal to zero, one can define the so-called truncated theory of the system (truncated equations).[19] Then, the BKL conjecture can be made more specific:

Weak conjecture: As the singularity is approached the terms containing space derivatives in the Einstein equations are negligible in comparison to the terms containing time derivatives. Thus, as the singularity is approached the Einstein equations approach those found by setting derivative terms to zero. Thus, the weak conjecture says that the Einstein equations can be well approximated by the truncated equations in the vicinity of the singularity. Note that this does not imply that the solutions of the full equations of motion will approach the solutions to the truncated equations as the singularity is approached. This additional condition is captured in the strong version as follows.

Strong conjecture: As the singularity is approached the Einstein equations approach those of the truncated theory and in addition the solutions to the full equations are well approximated by solutions to the truncated equations.

In the beginning, the BKL conjecture seemed to be coordinate-dependent and rather implausible. Barrow and Tipler,[20][21] for example, among the ten criticisms of BKL studies, include the inappropriate (according to them) choice of synchronous frame as a means to separate time and space derivatives. The BKL conjecture was sometimes rephrased in the literature as a statement that near the singularity only the time derivatives are important. Such a statement, taken at face value, is wrong or at best misleading since, as shown in the BKL analysis itself, spacelike gradients of the metric tensor cannot be neglected for generic solutions of pure Einstein gravity in four spacetime dimensions, and in fact play a crucial role in the appearance of the oscillatory regime. However, there exist reformulations of Einstein theory in terms of new variables involving the relevant gradients, for example in Ashtekar-like variables, for which the statement about the dominant role of the time derivatives is correct.[19] It is true that one gets at each spatial point an effective description of the singularity in terms of a finite dimensional dynamical system described by ordinary differential equations with respect to time, but the spatial gradients do enter these equations non-trivially.

Subsequent analysis by a large number of authors has shown that the BKL conjecture can be made precise and by now there is an impressive body of numerical and analytical evidence in its support.[22] It is fair to say that we are still quite far from a proof of the strong conjecture. But there has been outstanding progress in simpler models. In particular, Berger, Garfinkle, Moncrief, Isenberg, Weaver, and others showed that, in a class of models, as the singularity is approached the solutions to the full Einstein field equations approach the "velocity term dominated" (truncated) ones obtained by neglecting spatial derivatives.[22][23][24][25][26] Andersson and Rendall [27] showed that for gravity coupled to a massless scalar field or a stiff fluid, for every solution to the truncated equations there exists a solution to the full field equations that converges to the truncated solution as the singularity is approached, even in the absence of symmetries. These results were generalized to also include p-form gauge fields in.[28] In these truncated models the dynamics are simpler, allowing a precise statement of the conjecture that could be proven. In the general case, the strongest evidence to date comes from numerical evolutions. Berger and Moncrief began a program to analyze generic cosmological singularities.[29] While the initial work focused on symmetry reduced cases,[30] more recently Garfinkle[31] has performed numerical evolution of space-times with no symmetries in which, again, the mixmaster behavior is apparent. Finally, additional support for the conjecture has come from a numerical study of the behavior of test fields near the singularity of a Schwarzschild black hole.[32]

The Einstein equations for a universe with a homogeneous space can reduce to a system of ordinary differential equations containing only functions of time with the help of a frame field. To do this one must resolve the spatial components of four-vectors and four-tensors along the triad of basis vectors of the space:

${\displaystyle R_{(a)(b)}=R_{\alpha \beta }e_{(a)}^{\alpha }e_{(b)}^{\beta },\quad R_{0(a)}=R_{0\alpha }e_{(a)}^{\alpha },\quad u^{(a)}=u^{\alpha }e_{\alpha }^{(a)},}$

where all these quantities are now functions of t alone; the scalar quantities, the energy density ε and the pressure of the matter p, are also functions of the time.

The Einstein equations in vacuum in synchronous reference frame are[9][10]

${\displaystyle R_{0}^{0}=-{\frac {1}{2}}{\frac {\partial \varkappa _{\alpha }^{\alpha }}{\partial t}}-{\frac {1}{4}}\varkappa _{\alpha }^{\beta }\varkappa _{\beta }^{\alpha }=0,}$

(eq. 11)

${\displaystyle R_{\alpha }^{\beta }=-{\frac {1}{2{\sqrt {-g}}}}{\frac {\partial }{\partial t}}\left({\sqrt {-g}}\varkappa _{\alpha }^{\beta }\right)-P_{\alpha }^{\beta }=0,}$

(eq. 12)

${\displaystyle R_{\alpha }^{0}={\frac {1}{2}}\left(\varkappa _{\alpha ;\beta }^{\beta }-\varkappa _{\beta ;\alpha }^{\beta }\right)=0,}$

(eq. 13)

where ${\displaystyle \varkappa _{\alpha }^{\beta }}$  is the 3-dimensional tensor ${\displaystyle \varkappa _{\alpha }^{\beta }={\frac {\partial \gamma _{\alpha }^{\beta }}{\partial t}}}$ , and Pαβ is the 3-dimensional Ricci tensor, which is expressed by the 3-dimensional metric tensor γαβ in the same way as Rik is expressed by gik; Pαβ contains only the space (but not the time) derivatives of γαβ. Using triads, for eq. 11 one has simply

${\displaystyle \varkappa _{(a)(b)}={\frac {\partial \eta _{ab}}{\partial t}},\quad \varkappa _{(a)}^{(b)}={\frac {\partial \eta _{ac}}{\partial t}}\eta ^{cb}.}$

The components of P(a)(b) can be expressed in terms of the quantities ηab and the structure constants of the group by using the tetrad representation of the Ricci tensor in terms of quantities ${\displaystyle \lambda _{abc}=\left(e_{(a)i,k}-e_{(a)k,i}\right)e_{(b)}^{i}e_{(c)}^{k}}$ [33]

${\displaystyle R_{(a)(b)}=-{\frac {1}{2}}\left(\lambda _{ab,c}^{c}+\lambda _{ba,c}^{c}+\lambda _{ca,b}^{c}+\lambda _{cb,a}^{c}+\lambda _{b}^{cd}\lambda _{cda}+\lambda _{b}^{cd}\lambda _{dca}-{\frac {1}{2}}\lambda _{b}^{cd}\lambda _{acd}+\lambda _{cd}^{c}\lambda _{ab}^{d}+\lambda _{cd}^{c}\lambda _{ba}^{d}\right).}$

After replacing the three-index symbols ${\displaystyle \lambda _{bc}^{a}=C_{bc}^{a}}$  by two-index symbols Cab and the transformations:

${\displaystyle \eta _{ad}\eta _{bc}\eta _{cf}e^{def}=\eta e_{abc},\quad e_{abf}e^{cdf}=\delta _{a}^{c}\delta _{b}^{d}-\delta _{a}^{d}\delta _{b}^{c}}$

one gets the "homogeneous" Ricci tensor expressed in structure constants:

${\displaystyle P_{(a)}^{(b)}={\frac {1}{2\eta }}\left\{2C^{bd}C_{ad}+C^{db}C_{ad}+C^{bd}C_{da}-C_{d}^{d}\left(C_{a}^{b}+C_{a}^{b}\right)+\delta _{a}^{b}\left[\left(C_{d}^{d}\right)^{2}-2C^{df}C_{df}\right]\right\}.}$

Here, all indices are raised and lowered with the local metric tensor ηab

${\displaystyle C_{a}^{b}=\eta _{ac}C^{cb},\quad C_{ab}=\eta _{ac}\eta _{bd}C^{cd}.}$

The Bianchi identities for the three-dimensional tensor Pαβ in the homogeneous space take the form

${\displaystyle P_{b}^{c}C_{ca}^{b}+P_{a}^{c}C_{cb}^{b}=0.}$

Taking into account the transformations of covariant derivatives for arbitrary four-vectors Ai and four-tensors Aik

${\displaystyle A_{i;k}e_{(a)}^{i}e_{(b)}^{k}=A_{(a)(b)}-A^{(d)}\gamma _{dab},}$
${\displaystyle A_{ik;l}e_{(a)}^{i}e_{(b)}^{k}e_{(c)}^{l}=A_{(a)(b)(c)}-A_{(b)}^{(d)}\gamma _{dac}+A_{(a)}^{(d)}\gamma _{dbc},}$

the final expressions for the triad components of the Ricci four-tensor are:

${\displaystyle R_{0}^{0}=-{\frac {1}{2}}{\frac {\partial \varkappa _{(a)}^{(a)}}{\partial t}}-{\frac {1}{4}}\varkappa _{(a)}^{(b)}\varkappa _{(b)}^{(a)},}$

(eq. 11a)

${\displaystyle R_{(a)}^{(b)}=-{\frac {1}{2{\sqrt {\eta }}}}\left({\sqrt {\eta }}\varkappa _{(a)}^{(b)}-P_{(a)}^{(b)}\right),}$

(eq. 12a)

${\displaystyle R_{(a)}^{0}=-{\frac {1}{2}}\varkappa _{(b)}^{(c)}\left(C_{ca}^{b}-\delta _{a}^{b}C_{dc}^{d}\right).}$

(eq. 13a)

It should be emphasized that in setting up the Einstein equations there is thus no need to use explicit expressions for the basis vectors as functions of the coordinates.

### Kasner solution

Figure 3. Dynamics of Kasner metrics eq. 2 in spherical coordinates towards singularity. The Lifshitz-Khalatnikov parameter is u=2 (1/u=0.5) and the r coordinate is 2pα(1/u)τ where τ is logarithmic time: τ = ln t.[34] Shrinking along the axes is linear and anisotropic (no chaoticity).

Much more general solutions are obtained by a generalization of an exact particular solution derived by Edward Kasner[35] for a field in vacuum, in which the space is homogeneous and has a Euclidean metric that depends on time according to the Kasner metric

${\displaystyle dl^{2}=t^{2p_{1}}dx^{2}+t^{2p_{2}}dy^{2}+t^{2p_{3}}dz^{2}}$

(eq. 2)

(dl is the line element; dx, dy, dz are infinitesimal displacements in the 3 spatial dimensions, and t is time period passed since some initial moment t0 = 0). Here, p1, p2, p3 are any 3 numbers that satisfy the following Kasner conditions

${\displaystyle p_{1}+p_{2}+p_{3}=p_{1}^{2}+p_{2}^{2}+p_{3}^{2}=1.}$

(eq. 3)

Because of these relations, only 1 of the 3 numbers is independent (2 equations with 3 unknowns). All 3 numbers are never the same; 2 numbers are the same only in the sets of values ${\textstyle (-{\frac {1}{3}},{\frac {2}{3}},{\frac {2}{3}})}$  and (0, 0, 1).[36] In all other cases the numbers are different, one number is negative and the other two are positive. This is partially proved by squaring both sides of the first condition eq. 3 and developing the square:

${\displaystyle \left(p_{1}+p_{2}+p_{3}\right)^{2}=\left(p_{1}^{2}+p_{2}^{2}+p_{3}^{2}\right)+\left(2p_{1}p_{2}+2p_{2}p_{3}+2p_{1}p_{3}\right)=1}$

The term ${\displaystyle \left(p_{1}^{2}+p_{2}^{2}+p_{3}^{2}\right)}$  is equal to 1 by dint of the second condition eq. 3 and therefore the term with the mixed products should be zero. This is possible if at least one of the p1, p2, p3 is negative.

If the numbers are arranged in increasing order, p1 < p2 < p3, they change in the ranges (Fig. 4)

${\displaystyle {\begin{matrix}-{\tfrac {1}{3}}\leq p_{1}\leq 0,\\\ 0\leq p_{2}\leq {\tfrac {2}{3}},\\{\frac {2}{3}}\leq p_{3}\leq 1.\end{matrix}}}$

(eq. 4)

Figure 4. Plot of p1, p2, p3 with an argument 1/u. The numbers p1(u) and p3(u) are monotonously increasing while p2(u) is a monotonously decreasing function of u.

The Kasner metric eq. 2 corresponds to a flat homogenous but anisotropic space in which all volumes increase with time in such a way that the linear distances along two axes y and z increase while the distance along the axis x decreases. The moment t = 0 causes a singularity in the solution; the singularity in the metric at t = 0 cannot be avoided by any reference frame transformation. At the singularity, the invariants of the four-dimensional curvature tensor go to infinity. An exception is the case p1 = р2 = 0, р3 = 1; these values correspond to a flat spacetime: the transformation t sh z = ζ, t ch z = τ turns the metric eq. 2 into Galilean.

BKL parametrize the numbers p1, p2, p3 in terms of a single independent (real) parameter u (Lifshitz-Khalatnikov parameter[37]) as follows

${\displaystyle p_{1}(u)={\frac {-u}{1+u+u^{2}}},\ p_{2}(u)={\frac {1+u}{1+u+u^{2}}},\ p_{3}(u)={\frac {u(1+u)}{1+u+u^{2}}}}$

(eq. 5)

The Kasner index parametrization appears mysterious until one thinks about the two constraints on the indices eq. 3. Both constraints fix the overall scale of the indices so that only their ratios can vary. It is natural to pick one of those ratios as a new parameter, which can be done in six different ways. Picking u = u32 = p3 / p2, for example, it is trivial to express all six possible ratios in terms of it. Eliminating p3 = up2 first, and then using the linear constraint to eliminate p1 = 1 − p2up2 = 1 − (1 + u)p2, the quadratic constraint reduces to a quadratic equation in p2

${\displaystyle \left[1-\left(1+u\right)p_{2}\right]^{2}+p_{2}^{2}+\left(up_{2}\right)^{2}=1}$

(eq. 5a)

with roots p2 = 0 (obvious) and p2 = (1 + u) / (1 + u + u2), from which p1 and p3 are then obtained by back substitution. One can define six such parameters uab = pa / pb, for which pcpbpa when (c, b, a) is a cyclic permutation of (1, 2, 3).[38]

All different values of p1, p2, p3 ordered as above are obtained with u running in the range u ≥ 1. The values u < 1 are brought into this range according to

${\displaystyle p_{1}\left({\frac {1}{u}}\right)=p_{1}(u),\ p_{2}\left({\frac {1}{u}}\right)=p_{3}(u),\ p_{3}\left({\frac {1}{u}}\right)=p_{2}(u)}$

(eq. 6)

In the generalized solution, the form corresponding to eq. 2 applies only to the asymptotic metric (the metric close to the singularity t = 0), respectively, to the major terms of its series expansion by powers of t. In the synchronous reference frame it is written in the form of eq. 1 with a space distance element

${\displaystyle dl^{2}=\left(a^{2}l_{\alpha }l_{\beta }+b^{2}m_{\alpha }m_{\beta }+c^{2}n_{\alpha }n_{\beta }\right)dx^{\alpha }dx^{\beta }}$

(eq. 7)

where

${\displaystyle a=t^{p_{l}},\ b=t^{p_{m}},\ c=t^{p_{n}}}$

(eq. 8)

The three-dimensional vectors l, m, n define the directions at which space distance changes with time by the power laws eq. 8. These vectors, as well as the numbers pl, pm, pn which, as before, are related by eq. 3, are functions of the space coordinates. The powers pl, pm, pn are not arranged in increasing order, reserving the symbols p1, p2, p3 for the numbers in eq. 5 that remain arranged in increasing order. The determinant of the metric of eq. 7 is

${\displaystyle -g=a^{2}b^{2}c^{2}v^{2}=t^{2}v^{2}}$

(eq. 9)

where v = l[mn]. It is convenient to introduce the following quantities [39]

${\displaystyle \lambda ={\frac {\mathbf {l} \ \mathrm {curl} \ \mathbf {l} }{v}},\ \mu ={\frac {\mathbf {m} \ \mathrm {curl} \ \mathbf {m} }{v}},\ \nu ={\frac {\mathbf {n} \ \mathrm {curl} \ \mathbf {n} }{v}}.}$

(eq. 10)

The space metric in eq. 7 is anisotropic because the powers of t in eq. 8 cannot have the same values. On approaching the singularity at t = 0, the linear distances in each space element decrease in two directions and increase in the third direction. The volume of the element decreases in proportion to t.

The Kasner metric is introduced in the Einstein equations by substituting the respective metric tensor γαβ from eq. 7 without defining a priori the dependence of a, b, c from t:

${\displaystyle \varkappa _{\alpha }^{\beta }={\frac {2{\dot {a}}}{a}}l_{\alpha }l^{\beta }+{\frac {2{\dot {b}}}{b}}m_{\alpha }m^{\beta }+{\frac {2{\dot {c}}}{c}}n_{\alpha }n^{\beta }}$

where the dot above a symbol designates differentiation with respect to time. The Einstein equation eq. 11 takes the form

${\displaystyle -R_{0}^{0}={\frac {\ddot {a}}{a}}+{\frac {\ddot {b}}{b}}+{\frac {\ddot {c}}{c}}=0.}$

(eq. 14)

All its terms are to a second order for the large (at t → 0) quantity 1/t. In the Einstein equations eq. 12, terms of such order appear only from terms that are time-differentiated. If the components of Pαβ do not include terms of order higher than 2, then

${\displaystyle -R_{l}^{l}={\frac {({\dot {a}}bc){\dot {}}}{abc}}=0,\ -R_{m}^{m}={\frac {(a{\dot {b}}c){\dot {}}}{abc}}=0,\ -R_{n}^{n}={\frac {(ab{\dot {c}}){\dot {}}}{abc}}=0}$

(eq. 15)

where indices l, m, n designate tensor components in the directions l, m, n.[9] These equations together with eq. 14 give the expressions eq. 8 with powers that satisfy eq. 3.

However, the presence of 1 negative power among the 3 powers pl, pm, pn results in appearance of terms from Pαβ with an order greater than t−2. If the negative power is pl (pl = p1 < 0), then Pαβ contains the coordinate function λ and eq. 12 become

{\displaystyle {\begin{aligned}-R_{l}^{l}&={\frac {({\dot {a}}bc){\dot {}}}{abc}}+{\frac {\lambda ^{2}a^{2}}{2b^{2}c^{2}}}=0,\\-R_{m}^{m}&={\frac {(a{\dot {b}}c){\dot {}}}{abc}}-{\frac {\lambda ^{2}a^{2}}{2b^{2}c^{2}}}=0,\\-R_{n}^{n}&={\frac {(ab{\dot {c}}){\dot {}}}{abc}}-{\frac {\lambda ^{2}a^{2}}{2b^{2}c^{2}}}=0.\\\end{aligned}}}

(eq. 16)

Here, the second terms are of order t−2(pm + pnpl) whereby pm + pnpl = 1 + 2 |pl| > 1.[40] To remove these terms and restore the metric eq. 7, it is necessary to impose on the coordinate functions the condition λ = 0.

The remaining 3 Einstein equations eq. 13 contain only first order time derivatives of the metric tensor. They give 3 time-independent relations that must be imposed as necessary conditions on the coordinate functions in eq. 7. This, together with the condition λ = 0, makes 4 conditions. These conditions bind 10 different coordinate functions: 3 components of each of the vectors l, m, n, and one function in the powers of t (any one of the functions pl, pm, pn, which are bound by the conditions eq. 3). When calculating the number of physically arbitrary functions, it must be taken into account that the synchronous system used here allows time-independent arbitrary transformations of the 3 space coordinates. Therefore, the final solution contains overall 10 − 4 − 3 = 3 physically arbitrary functions which is 1 less than what is needed for the general solution in vacuum.

The degree of generality reached at this point is not lessened by introducing matter; matter is written into the metric eq. 7 and contributes 4 new coordinate functions necessary to describe the initial distribution of its density and the 3 components of its velocity. This makes possible to determine matter evolution merely from the laws of its movement in an a priori given gravitational field which are the hydrodynamic equations

${\displaystyle {\frac {1}{\sqrt {-g}}}{\frac {\partial }{\partial x^{i}}}\left({\sqrt {-g}}\sigma u^{i}\right)=0,}$

(eq. 17)

${\displaystyle (p+\varepsilon )u^{k}\left\{{\frac {\partial u_{i}}{\partial x^{k}}}-{\frac {1}{2}}u^{l}{\frac {\partial g_{kl}}{\partial x^{i}}}\right\rbrace =-{\frac {\partial p}{\partial x^{i}}}-u_{i}u^{k}{\frac {\partial p}{\partial x^{k}}},}$

(eq. 18)

where ui is the 4-dimensional velocity, ε and σ are the densities of energy and entropy of matter.[41] For the ultrarelativistic equation of state p = ε/3 the entropy σ ~ ε1/4. The major terms in eq. 17 and eq. 18 are those that contain time derivatives. From eq. 17 and the space components of eq. 18 one has

${\displaystyle {\frac {\partial }{\partial t}}\left({\sqrt {-g}}u_{0}\varepsilon ^{\frac {3}{4}}\right)=0,\ 4\varepsilon \cdot {\frac {\partial u_{\alpha }}{\partial t}}+u_{\alpha }\cdot {\frac {\partial \varepsilon }{\partial t}}=0,}$

resulting in

${\displaystyle abcu_{0}\varepsilon ^{\frac {3}{4}}=\mathrm {const} ,\ u_{\alpha }\varepsilon ^{\frac {1}{4}}=\mathrm {const} ,}$

(eq. 19)

where 'const' are time-independent quantities. Additionally, from the identity uiui = 1 one has (because all covariant components of uα are to the same order)

${\displaystyle u_{0}^{2}\approx u_{n}u^{n}={\frac {u_{n}^{2}}{c^{2}}},}$

where un is the velocity component along the direction of n that is connected with the highest (positive) power of t (supposing that pn = p3). From the above relations, it follows that

${\displaystyle \varepsilon \sim {\frac {1}{a^{2}b^{2}}},\ u_{\alpha }\sim {\sqrt {ab}}}$

(eq. 20)

or

${\displaystyle \varepsilon \sim t^{-2(p_{1}+p_{2})}=t^{-2(1-p_{3})},\ u_{\alpha }\sim t^{\frac {(1-p_{3})}{2}}.}$

(eq. 21)

The above equations can be used to confirm that the components of the matter stress-energy-momentum tensor standing in the right hand side of the equations

${\displaystyle R_{0}^{0}=T_{0}^{0}-{\frac {1}{2}}T,\ R_{\alpha }^{\beta }=T_{\alpha }^{\beta }-{\frac {1}{2}}\delta _{\alpha }^{\beta }T,}$

are, indeed, to a lower order by 1/t than the major terms in their left hand sides. In the equations ${\displaystyle R_{\alpha }^{0}=T_{\alpha }^{0}}$  the presence of matter results only in the change of relations imposed on their constituent coordinate functions.[9]

The fact that ε becomes infinite by the law eq. 21 confirms that in the solution to eq. 7 one deals with a physical singularity at any values of the powers p1, p2, p3 excepting only (0, 0, 1). For these last values, the singularity is non-physical and can be removed by a change of reference frame.

The fictional singularity corresponding to the powers (0, 0, 1) arises as a result of time line coordinates crossing over some 2-dimensional "focal surface". As pointed out in,[9] a synchronous reference frame can always be chosen in such a way that this inevitable time line crossing occurs exactly on such surface (instead of a 3-dimensional caustic surface). Therefore, a solution with such simultaneous for the whole space fictional singularity must exist with a full set of arbitrary functions needed for the general solution. Close to the point t = 0 it allows a regular expansion by whole powers of t.[42]

## Oscillating mode towards the singularity

The four conditions that had to be imposed on the coordinate functions in the solution eq. 7 are of different types: three conditions that arise from the equations ${\displaystyle R_{\alpha }^{0}}$ = 0 are "natural"; they are a consequence of the structure of Einstein equations. However, the additional condition λ = 0 that causes the loss of one derivative function, is of entirely different type.

The general solution by definition is completely stable; otherwise the Universe would not exist. Any perturbation is equivalent to a change in the initial conditions in some moment of time; since the general solution allows arbitrary initial conditions, the perturbation is not able to change its character. In other words, the existence of the limiting condition λ = 0 for the solution of eq. 7 means instability caused by perturbations that break this condition. The action of such perturbation must bring the model to another mode which thereby will be most general. Such perturbation cannot be considered as small: a transition to a new mode exceeds the range of very small perturbations.

The analysis of the behavior of the model under perturbative action, performed by BKL, delineates a complex oscillatory mode on approaching the singularity.[1][43][44][45] They could not give all details of this mode in the broad frame of the general case. However, BKL explained the most important properties and character of the solution on specific models that allow far-reaching analytical study.

These models are based on a homogeneous space metric of a particular type. Supposing a homogeneity of space without any additional symmetry leaves a great freedom in choosing the metric. All possible homogeneous (but anisotropic) spaces are classified, according to Bianchi, in 9 classes.[46] BKL investigate only spaces of Bianchi Types VIII and IX.

If the metric has the form of eq. 7, for each type of homogeneous spaces exists some functional relation between the reference vectors l, m, n and the space coordinates. The specific form of this relation is not important. The important fact is that for Type VIII and IX spaces, the quantities λ, μ, ν eq. 10 are constants while all "mixed" products l rot m, l rot n, m rot l, etc. are zeros. For Type IX spaces, the quantities λ, μ, ν have the same sign and one can write λ = μ = ν = 1 (the simultaneous sign change of the 3 constants does not change anything). For Type VIII spaces, 2 constants have a sign that is opposite to the sign of the third constant; one can write, for example, λ = − 1, μ = ν = 1.[47]

The study of the effect of the perturbation on the "Kasner mode" is thus confined to a study on the effect of the λ-containing terms in the Einstein equations. Type VIII and IX spaces are the most suitable models exactly in this connection. Since all 3 quantities λ, μ, ν differ from zero, the condition λ = 0 does not hold irrespective of which direction l, m, n has negative power law time dependence.

The Einstein equations for the Type VIII and Type IX space models are[48]

{\displaystyle {\begin{aligned}-R_{l}^{l}&={\frac {\left({\dot {a}}bc\right){\dot {}}}{abc}}+{\frac {1}{2}}\left(a^{2}b^{2}c^{2}\right)\left[\lambda ^{2}a^{4}-\left(\mu b^{2}-\nu c^{2}\right)^{2}\right]=0,\\-R_{m}^{m}&={\frac {(a{\dot {b}}c){\dot {}}}{abc}}+{\frac {1}{2}}\left(a^{2}b^{2}c^{2}\right)\left[\mu ^{2}b^{4}-\left(\lambda a^{2}-\nu c^{2}\right)^{2}\right]=0,\\-R_{n}^{n}&={\frac {\left(ab{\dot {c}}\right){\dot {}}}{abc}}+{\frac {1}{2}}\left(a^{2}b^{2}c^{2}\right)\left[\nu ^{2}c^{4}-\left(\lambda a^{2}-\mu b^{2}\right)^{2}\right]=0,\\\end{aligned}}}

(eq. 22)

${\displaystyle -R_{0}^{0}={\frac {\ddot {a}}{a}}+{\frac {\ddot {b}}{b}}+{\frac {\ddot {c}}{c}}=0}$

(eq. 23)

(the remaining components ${\displaystyle R_{l}^{0}}$ , ${\displaystyle R_{m}^{0}}$ , ${\displaystyle R_{n}^{0}}$ , ${\displaystyle R_{l}^{m}}$ , ${\displaystyle R_{l}^{n}}$ , ${\displaystyle R_{m}^{n}}$  are identically zeros). These equations contain only functions of time; this is a condition that has to be fulfilled in all homogeneous spaces. Here, the eq. 22 and eq. 23 are exact and their validity does not depend on how near one is to the singularity at t = 0.[49]

The time derivatives in eq. 22 and eq. 23 take a simpler form if а, b, с are substituted by their logarithms α, β, γ:

${\displaystyle a=e^{\alpha },\ b=e^{\beta },\ c=e^{\gamma },}$

(eq. 24)

substituting the variable t for τ according to:

${\displaystyle dt=abc\ d\tau .}$

(eq. 25)

Then (subscripts denote differentiation by τ):

{\displaystyle {\begin{aligned}2\alpha _{\tau \tau }&=\left(\mu b^{2}-\nu c^{2}\right)^{2}-\lambda ^{2}a^{4}=0,\\2\beta _{\tau \tau }&=\left(\lambda a^{2}-\nu c^{2}\right)^{2}-\mu ^{2}b^{4}=0,\\2\gamma _{\tau \tau }&=\left(\lambda a^{2}-\mu b^{2}\right)^{2}-\nu ^{2}c^{4}=0,\\\end{aligned}}}

(eq. 26)

${\displaystyle {\frac {1}{2}}\left(\alpha +\beta +\gamma \right)_{\tau \tau }=\alpha _{\tau }\beta _{\tau }+\alpha _{\tau }\gamma _{\tau }+\beta _{\tau }\gamma _{\tau }.}$

(eq. 27)

Adding together equations eq. 26 and substituting in the left hand side the sum (α + β + γ)τ τ according to eq. 27, one obtains an equation containing only first derivatives which is the first integral of the system eq. 26:

${\displaystyle \alpha _{\tau }\beta _{\tau }+\alpha _{\tau }\gamma _{\tau }+\beta _{\tau }\gamma _{\tau }={\frac {1}{4}}\left(\lambda ^{2}a^{4}+\mu ^{2}b^{4}+\nu ^{2}c^{4}-2\lambda \mu a^{2}b^{2}-2\lambda \nu a^{2}c^{2}-2\mu \nu b^{2}c^{2}\right).}$

(eq. 28)

This equation plays the role of a binding condition imposed on the initial state of eq. 26. The Kasner mode eq. 8 is a solution of eq. 26 when ignoring all terms in the right hand sides. But such situation cannot go on (at t → 0) indefinitely because among those terms there are always some that grow. Thus, if the negative power is in the function a(t) (pl = p1) then the perturbation of the Kasner mode will arise by the terms λ2a4; the rest of the terms will decrease with decreasing t. If only the growing terms are left in the right hand sides of eq. 26, one obtains the system:

${\displaystyle \alpha _{\tau \tau }=-{\frac {1}{2}}\lambda ^{2}e^{4\alpha },\ \beta _{\tau \tau }=\gamma _{\tau \tau }={\frac {1}{2}}\lambda ^{2}e^{4\alpha }}$

(eq. 29)

(compare eq. 16; below it is substituted λ2 = 1). The solution of these equations must describe the metric evolution from the initial state, in which it is described by eq. 8 with a given set of powers (with pl < 0); let pl = р1, pm = р2, pn = р3 so that

${\displaystyle a\sim t^{p_{1}},\ b\sim t^{p_{2}},\ c\sim t^{p_{3}}.}$

(eq. 30)

Then

${\displaystyle abc=\Lambda t,\ \tau =\Lambda ^{-1}\ln t+\mathrm {const} }$

(eq. 31)

where Λ is constant. Initial conditions for eq. 29 are redefined as[50]

${\displaystyle \alpha _{\tau }=\Lambda p_{1},\ \beta _{\tau }=\Lambda p_{2},\ \gamma _{\tau }=\Lambda p_{3}\ \mathrm {at} \ \tau \to \infty }$

(eq. 32)

Equations eq. 29 are easily integrated; the solution that satisfies the condition eq. 32 is

${\displaystyle {\begin{cases}a^{2}={\frac {2|p_{1}|\Lambda }{\operatorname {ch} (2|p_{1}|\Lambda \tau )}},\\b^{2}=b_{0}^{2}e^{2\Lambda (p_{2}-|p_{1}|)\tau }\operatorname {ch} (2|p_{1}|\Lambda \tau ),\\c^{2}=c_{0}^{2}e^{2\Lambda (p_{2}-|p_{1}|)\tau }\operatorname {ch} (2|p_{1}|\Lambda \tau ),\end{cases}}}$

(eq. 33)

where b0 and c0 are two more constants.

It can easily be seen that the asymptotic of functions eq. 33 at t → 0 is eq. 30. The asymptotic expressions of these functions and the function t(τ) at τ → −∞ is[51]

${\displaystyle a\sim e^{-\Lambda p_{1}\tau },\ b\sim e^{\Lambda (p_{2}+2p_{1})\tau },\ c\sim e^{\Lambda (p_{3}+2p_{1})\tau },\ t\sim e^{\Lambda (1+2p_{1})\tau }.}$

Expressing a, b, c as functions of t, one has

${\displaystyle a\sim t^{p'_{l}},b\sim t^{p'_{m}},c\sim t^{p'_{n}}}$

(eq. 34)

where

${\displaystyle p'_{l}={\frac {|p_{1}|}{1-2|p_{1}|}},p'_{m}=-{\frac {2|p_{1}|-p_{2}}{1-2|p_{1}|}},p'_{n}={\frac {p_{3}-2|p_{1}|}{1-2|p_{1}|}}.}$

(eq. 35)

Then

${\displaystyle abc=\Lambda 't,\ \Lambda '=(1-2|p_{1}|)\Lambda .}$

(eq. 36)

The above shows that perturbation acts in such a way that it changes one Kasner mode with another Kasner mode, and in this process the negative power of t flips from direction l to direction m: if before it was pl < 0, now it is p'm < 0. During this change the function a(t) passes through a maximum and b(t) passes through a minimum; b, which before was decreasing, now increases: a from increasing becomes decreasing; and the decreasing c(t) decreases further. The perturbation itself (λ2a in eq. 29), which before was increasing, now begins to decrease and die away. Further evolution similarly causes an increase in the perturbation from the terms with μ2 (instead of λ2) in eq. 26, next change of the Kasner mode, and so on.

It is convenient to write the power substitution rule eq. 35 with the help of the parametrization eq. 5:

${\displaystyle {\begin{matrix}\mathrm {if} &p_{l}=p_{1}(u)&p_{m}=p_{2}(u)&p_{n}=p_{3}(u)\\\mathrm {then} &p'_{l}=p_{2}(u-1)&p'_{m}=p_{1}(u-1)&p'_{n}=p_{3}(u-1)\end{matrix}}}$

(eq. 37)

The greater of the two positive powers remains positive.

BKL call this flip of negative power between directions a Kasner epoch. The key to understanding the character of metric evolution on approaching singularity is exactly this process of Kasner epoch alternation with flipping of powers pl, pm, pn by the rule eq. 37.

The successive alternations eq. 37 with flipping of the negative power p1 between directions l and m (Kasner epochs) continues by depletion of the whole part of the initial u until the moment at which u < 1. The value u < 1 transforms into u > 1 according to eq. 6; in this moment the negative power is pl or pm while pn becomes the lesser of two positive numbers (pn = p2). The next series of Kasner epochs then flips the negative power between directions n and l or between n and m. At an arbitrary (irrational) initial value of u this process of alternation continues unlimited.[52]

In the exact solution of the Einstein equations, the powers pl, pm, pn lose their original, precise, sense. This circumstance introduces some "fuzziness" in the determination of these numbers (and together with them, to the parameter u) which, although small, makes meaningless the analysis of any definite (for example, rational) values of u. Therefore, only these laws that concern arbitrary irrational values of u have any particular meaning.

The larger periods in which the scales of space distances along two axes oscillate while distances along the third axis decrease monotonously, are called eras; volumes decrease by a law close to ~ t. On transition from one era to the next, the direction in which distances decrease monotonously, flips from one axis to another. The order of these transitions acquires the asymptotic character of a random process. The same random order is also characteristic for the alternation of the lengths of successive eras (by era length, BKL understand the number of Kasner epoch that an era contains, and not a time interval).

To each era (s-th era) correspond a series of values of the parameter u starting from the greatest, ${\displaystyle u_{\max }^{(s)}}$ , and through the values ${\displaystyle u_{\max }^{(s)}}$  − 1, ${\displaystyle u_{\max }^{(s)}}$  − 2, ..., reaching to the smallest, ${\displaystyle u_{\min }^{(s)}}$  < 1. Then

${\displaystyle u_{\min }^{(s)}=x^{(s)},\ u_{\max }^{(s)}=k^{(s)}+x^{(s)},}$

(eq. 41)

that is, k(s) = [${\displaystyle u_{\max }^{(s)}}$ ] where the brackets mean the whole part of the value. The number k(s) is the era length, measured by the number of Kasner epochs that the era contains. For the next era

${\displaystyle u_{\max }^{(s+1)}={\frac {1}{x^{(s)}}},\ k^{(s+1)}=\left[{\frac {1}{x^{(s)}}}\right].}$

(eq. 42)

In the limiteless series of numbers u, composed by these rules, there are infinitesimally small (but never zero) values x(s) and correspondingly infinitely large lengths k(s).

The era series become denser on approaching t = 0. However, the natural variable for describing the time course of this evolution is not the world time t, but its logarithm, ln t, by which the whole process of reaching the singularity is extended to −∞.

According to eq. 33, one of the functions a, b, c, that passes through a maximum during a transition between Kasner epochs, at the peak of its maximum is

${\displaystyle a_{\max }={\sqrt {2\Lambda |p_{1}(u)|}}}$

(eq. 38)

where it is supposed that amax is large compared to b0 and c0; in eq. 38 u is the value of the parameter in the Kasner epoch before transition. It can be seen from here that the peaks of consecutive maxima during each era are gradually lowered. Indeed, in the next Kasner epoch this parameter has the value u' = u − 1, and Λ is substituted according to eq. 36 with Λ' = Λ(1 − 2|p1(u)|). Therefore, the ratio of 2 consecutive maxima is

${\displaystyle {\frac {a'_{\max }}{a_{\max }}}=\left[{\frac {p_{1}(u-1)}{p_{1}(u)}}\left(1-2|p_{1}(u)|\right)\right]^{\frac {1}{2}};}$

and finally

${\displaystyle {\frac {a'_{\max }}{a_{\max }}}={\sqrt {\frac {u-1}{u}}}\equiv {\sqrt {\frac {u'}{u}}}.}$

(eq. 39)

The above are solutions to Einstein equations in vacuum. As for the pure Kasner mode, matter does not change the qualitative properties of this solution and can be written into it disregarding its reaction on the field. However, if one does this for the model under discussion, understood as an exact solution of the Einstein equations, the resulting picture of matter evolution would not have a general character and would be specific for the high symmetry imminent to the present model. Mathematically, this specificity is related to the fact that for the homogeneous space geometry discussed here, the Ricci tensor components ${\displaystyle R_{\alpha }^{0}}$  are identically zeros and therefore the Einstein equations would not allow movement of matter (which gives non-zero stress energy-momentum tensor components ${\displaystyle T_{\alpha }^{0}}$ ). In other words, the synchronous frame must also be co-moving with respect to matter. If one substitutes in eq. 19 uα = 0, u0 = 1, it becomes ε ~ (abc)−4/3 ~ t−4/3.

This difficulty is avoided if one includes in the model only the major terms of the limiting (at t → 0) metric and writes into it a matter with arbitrary initial distribution of densities and velocities. Then the course of evolution of matter is determined by its general laws of movement eq. 17 and eq. 18 that result in eq. 21. During each Kasner epoch, density increases by the law

${\displaystyle \varepsilon =t^{-2(1-p_{3})},}$

(eq. 40)

where p3 is, as above, the greatest of the numbers p1, p2, p3. Matter density increases monotonously during all evolution towards the singularity.

## Metric evolution

Very large u values correspond to Kasner powers

${\displaystyle p_{1}\approx -{\frac {1}{u}},\ p_{2}\approx {\frac {1}{u}},\ p_{3}\approx 1-{\frac {1}{u^{2}}},}$

(eq. 43)

which are close to the values (0, 0, 1). Two values that are close to zero, are also close to each other, and therefore the changes in two out of the three types of "perturbations" (the terms with λ, μ and ν in the right hand sides of eq. 26) are also very similar. If in the beginning of such long era these terms are very close in absolute values in the moment of transition between two Kasner epochs (or made artificially such by assigning initial conditions) then they will remain close during the greatest part of the length of the whole era. In this case (BKL call this the case of small oscillations), analysis based on the action of one type of perturbations becomes incorrect; one must take into account the simultaneous effect of two perturbation types.

### Two perturbations

Consider a long era, during which 2 out of the 3 functions a, b, c (let them be a and b) undergo small oscillations while the third function (c) decreases monotonously. The latter function quickly becomes small; consider the solution just in the region where one can ignore c in comparison to a and b. The calculations are first done for the Type IX space model by substituting accordingly λ = μ = ν = 1.[44]

After ignoring function c, the first 2 equations eq. 26 give

${\displaystyle \alpha _{\tau \tau }+\beta _{\tau \tau }=0,\,}$

(eq. 44)

${\displaystyle \alpha _{\tau \tau }-\beta _{\tau \tau }=e^{4\beta }-e^{4\alpha },\,}$

(eq. 45)

and eq. 28 can be used as a third equation, which takes the form

${\displaystyle \gamma _{\tau \tau }\left(\alpha _{\tau \tau }+\beta _{\tau \tau }\right)=-\alpha _{\tau }\beta _{\tau }+{\frac {1}{4}}\left(e^{2\alpha }-e^{2\beta }\right)^{2}.}$

(eq. 46)

The solution of eq. 44 is written in the form

${\displaystyle \alpha +\beta ={\frac {2a_{0}^{2}}{\xi _{0}}}\left(\tau -\tau _{0}\right)+2\ln a_{0},}$

where α0, ξ0 are positive constants, and τ0 is the upper limit of the era for the variable τ. It is convenient to introduce further a new variable (instead of τ)

${\displaystyle \xi =\xi _{0}\exp \left[{\frac {2a_{0}^{2}}{\xi _{0}}}\left(\tau -\tau _{0}\right)\right].}$

(eq. 47)

Then

${\displaystyle \alpha +\beta =\ln {\frac {\xi }{\xi _{0}}}+2\ln a_{0}.}$

(eq. 48)

Equations eq. 45 and eq. 46 are transformed by introducing the variable χ = α − β:

${\displaystyle \chi _{\xi \xi }={\frac {\chi _{\xi }}{\xi }}+{\frac {1}{2}}\operatorname {sh} 2\chi =0,}$

(eq. 49)

${\displaystyle \gamma _{\xi }=-{\frac {1}{4}}\xi +{\frac {1}{8}}\xi \left(2\chi _{\xi }^{2}+\operatorname {ch} 2\chi -1\right).}$

(eq. 50)

Decrease of τ from τ0 to −∞ corresponds to a decrease of ξ from ξ0 to 0. The long era with close a and b (that is, with small χ), considered here, is obtained if ξ0 is a very large quantity. Indeed, at large ξ the solution of eq. 49 in the first approximation by 1/ξ is

${\displaystyle \chi =\alpha -\beta ={\frac {2A}{\sqrt {\xi }}}\sin \left(\xi -\xi _{0}\right),}$

(eq. 51)

where A is constant; the multiplier ${\displaystyle {\tfrac {1}{\sqrt {\xi }}}}$  makes χ a small quantity so it can be substituted in eq. 49 by sh 2χ ≈ 2χ.[53]

From eq. 50 one obtains

${\displaystyle \gamma _{\xi }={\frac {1}{4}}\xi \left(2\chi _{\xi }^{2}+\chi ^{2}\right)=A^{2},\ \gamma =A^{2}\left(\xi -\xi _{0}\right)+\mathrm {const} .}$

After determining α and β from eq. 48 and eq. 51 and expanding eα and eβ in series according to the above approximation, one obtains finally:[54]

${\displaystyle {\begin{cases}a\\b\end{cases}}=a_{0}{\sqrt {\frac {\xi }{\xi _{0}}}}\left[1\pm {\frac {A}{\sqrt {\xi }}}\sin \left(\xi -\xi _{0}\right)\right],}$

(eq. 52)

${\displaystyle c=c_{0}e^{-A^{2}\left(\xi _{0}-\xi \right)}.}$

(eq. 53)

The relation between the variable ξ and time t is obtained by integration of the definition dt = abc dτ which gives

${\displaystyle {\frac {t}{t_{0}}}=e^{-A^{2}\left(\xi _{0}-\xi \right)}.}$

(eq. 54)

The constant c0 (the value of с at ξ = ξ0) should be now c0 ${\displaystyle \ll }$  α0·

Figure 5. Bianchi type VIII (open) space undergoing a chaotic BKL (Mixmaster) dynamics close to singularity according to rules eq. 35 with initial ${\displaystyle u={\sqrt {13}}}$ . The singularity is at the central pinch of the hyperboloid surface.

Let us now consider the domain ξ ${\displaystyle \ll }$  1. Here the major terms in the solution of eq. 49 are:

${\displaystyle \chi =\alpha -\beta =k\ln \xi +\mathrm {const} ,\,}$

where k is a constant in the range − 1 < k < 1; this condition ensures that the last term in eq. 49 is small (sh 2χ contains ξ2k and ξ−2k). Then, after determining α, β, and t, one obtains

${\displaystyle a\sim \xi ^{\frac {1+k}{2}},\ b\sim \xi ^{\frac {1-k}{2}},\ c\sim \xi ^{-{\frac {1-k^{2}}{4}}},\ t\sim \xi ^{\frac {3+k^{2}}{4}}.}$

(eq. 55)

This is again a Kasner mode with the negative t power coming into the function c(t).[55]

These results picture an evolution that is qualitatively similar to that, described above. During a long period of time that corresponds to a large decreasing ξ value, the two functions a and b oscillate, remaining close in magnitude ${\displaystyle {\tfrac {a-b}{a}}\sim {\tfrac {1}{\sqrt {\xi }}}}$ ; in the same time, both functions a and b slowly (${\displaystyle \sim {\sqrt {\xi }}}$ ) decrease. The period of oscillations is constant by the variable ξ : Δξ = 2π (or, which is the same, with a constant period by logarithmic time: Δ ln t = 2πΑ2). The third function, c, decreases monotonously by a law close to c = c0t/t0.

This evolution continues until ξ ≈1 and formulas eq. 52 and eq. 53 are no longer applicable. Its time duration corresponds to change of t from t0 to the value t1, related to ξ0 according to

${\displaystyle A^{2}\xi _{0}=\ln {\frac {t_{0}}{t_{1}}}.}$

(eq. 56)

The relationship between ξ and t during this time can be presented in the form

${\displaystyle {\frac {\xi }{\xi _{0}}}={\frac {\ln {\tfrac {t}{t_{1}}}}{\ln {\tfrac {t_{0}}{t_{1}}}}}.}$

(eq. 57)

After that, as seen from eq. 55, the decreasing function c starts to increase while functions a and b start to decrease. This Kasner epoch continues until terms c2/a2b2 in eq. 22 become ~ t2 and a next series of oscillations begins.

The law for density change during the long era under discussion is obtained by substitution of eq. 52 in eq. 20:

${\displaystyle \varepsilon \sim \left({\frac {\xi _{0}}{\xi }}\right)^{2}.}$

(eq. 58)

When ξ changes from ξ0 to ξ ≈1, the density increases ${\displaystyle \xi _{0}^{2}}$  times.

It must be stressed that although the function c(t) changes by a law, close to c ~ t, the metric eq. 52 does not correspond to a Kasner metric with powers (0, 0, 1). The latter corresponds to an exact solution (found by Taub[56]) which is allowed by eqs. 26'–'27 and in which

${\displaystyle a^{2}=b^{2}={\frac {p}{2}}{\frac {\mathrm {ch} (2p\tau +\delta _{1})}{\mathrm {ch} ^{2}(p\tau +\delta _{2})}},\;c^{2}={\frac {2p}{\mathrm {ch} (2p\tau +\delta _{1})}},}$

(eq. 59)

where p, δ1, δ2 are constant. In the asymptotic region τ → −∞, one can obtain from here a = b = const, c = const.t after the substitution ерτ = t. In this metric, the singularity at t = 0 is non-physical.

Let us now describe the analogous study of the Type VIII model, substituting in eqs. eqs. 26'–'28 λ = −1, μ = ν = 1.[45]

If during the long era, the monotonically decreasing function is a, nothing changes in the foregoing analysis: ignoring a2 on the right side of equations 26 and 28, goes back to the same equations 49 and 50 (with altered notation). Some changes occur, however, if the monotonically decreasing function is b or c; let it be c.

As before, one has equation 49 with the same symbols, and, therefore, the former expressions eq. 52 for the functions a(ξ) and b(ξ), but equation 50 is replaced by

${\displaystyle \gamma _{\xi }=-{\frac {1}{4}}\xi +{\frac {1}{8}}\xi \left(2\chi _{\xi }^{2}+\mathrm {ch} 2\chi +1\right).}$

(eq. 60)

The major term at large ξ now becomes

${\displaystyle \gamma _{\xi }\approx {\frac {1}{8}}\xi \cdot 2,\quad \gamma \approx {\frac {1}{8}}\left(\xi ^{2}-\xi _{0}^{2}\right),}$

so that

${\displaystyle {\frac {c}{c_{0}}}={\frac {t}{t_{0}}}=e^{-{\frac {1}{8}}\left(\xi _{0}^{2}-\xi ^{2}\right)}.}$

(eq. 61)

The value of c as a function of time t is again c = c0t/t0 but the time dependence of ξ changes. The length of a long era depends on ξ0 according to

${\displaystyle \xi _{0}={\sqrt {8\ln {\frac {t}{t_{0}}}}}.}$

(eq. 62)

On the other hand, the value ξ0 determines the number of oscillations of the functions a and b during an era (equal to ξ0/2π). Given the length of an era in logarithmic time (i.e., with given ratio t0/t1) the number of oscillations for Type VIII will be, generally speaking, less than for Type IX. For the period of oscillations one gets now Δ ln t = πξ/2; contrary to Type IX, the period is not constant throughout the long era, and slowly decreases along with ξ.

### The small-time domain

As shown above, long eras violate the "regular" course of evolution; this fact makes it difficult to study the evolution of time intervals, encompassing several eras. It can be shown, however, that such "abnormal" cases appear in the spontaneous evolution of the model to a singular point in the asymptotically small times t at sufficiently large distances from a start point with arbitrary initial conditions. Even in long eras both oscillatory functions during transitions between Kasner epochs remain so different that the transition occurs under the influence of only one perturbation. All results in this section relate equally to models of the types VIII and IX.[57]

During each Kasner epoch abc = Λt, i. e. α + β + γ = ln Λ + ln t. On changing over from one epoch (with a given value of the parameter u) to the next epoch the constant Λ is multiplied by 1 + 2p1 = (1 - u + u2)/(1 + u + u2) < 1. Thus a systematic decrease in Λ takes place. But it is essential that the mean (with respect to the lengths k of eras) value of the entire variation of ln Λ during an era is finite. Actually the divergence of the mean value could be due only to a too rapid increase of this variation with increasing k. For large value of the parameter u, ln(1 + 2p1) ≈ −2/u. For a large k the maximal value u(max) = k + x ≈ k. Hence the entire variation of ln Λ during an era is given by a sum of the form

${\displaystyle \sum \ln \left(1+2p_{1}\right)=\dots +{\frac {1}{k-2}}+{\frac {1}{k-1}}+{\frac {1}{k}}}$

with only the terms that correspond to large values of u written down. When k increases this sum increases as ln k. But the probability for an appearance of an era of a large length k decreases as 1/k2 according to eq. 76; hence the mean value of the sum above is finite. Consequently, the systematic variation of the quantity ln Λ over a large number of eras will be proportional to this number. But it is seen in eq. 85 that with t → 0 the number s increases merely as ln |ln t|. Thus in the asymptotic limit of arbitrarily small t the term ln Λ can indeed be neglected as compared to ln t. In this approximation [58]

${\displaystyle \alpha +\beta +\gamma =-\Omega ,\,}$

(eq. 63)

where Ω denotes the "logarithmic time"

${\displaystyle \Omega =-\ln t.\,}$

(eq. 64)

and the process of epoch transitions can be regarded as a series of brief time flashes. The magnitudes of maxima of the oscillating scale functions are also subject to a systematic variation. From eq. 39 for u ≫ 1 we find that ${\displaystyle a_{\max }^{\prime }-a_{\max }\approx -1/2u}$ . In the same way as it was done above for the quantity ln Λ, one can hence deduce that the mean decrease in the height of the maxima during an era is finite and the total decrease over a large number of eras increases with t → 0 merely as ln Ω. At the same time the lowering of the minima, and by the same token the increase of the amplitude of the oscillations, proceed (eq. 77) proportional to Ω. In correspondence with the adopted approximation the lowering of the maxima is neglected in comparison with the increase of the amplitudes so αmax = 0, βmax = 0, γmax = 0 for the maximal values of all oscillating functions and the quantities α, β, γ run only through negative values that are connected with one another at each instant of time by the relation eq. 63.

Figure 4. Variation of α, β, and γ as functions of the logarithmic time Ω during one era. The vertical dash lines denote alterations of the Kasner epochs, corresponding to the linear segments of the curves. On the top are indicated the values of the parameter u that determine the Kasner exponents. The last epoch has a longer duration if x is small. In the first epoch of the next era, γ begins to increase and α becomes a monotonically decreasing function.

Considering such instant change of epochs, the transition periods are ignored as small in comparison to the epoch length; this condition is actually fulfilled.[59] Replacement of α, β, and γ maxima with zeros requires that quantities ln (|p1|Λ) be small in comparison with the amplitudes of oscillations of the respective functions. As mentioned above, during transitions between eras |p1| values can become very small while their magnitude and probability for occurrence are not related to the oscillation amplitudes in the respective moment. Therefore, in principle, it is possible to reach so small |p1| values that the above condition (zero maxima) is violated. Such drastic drop of αmax can lead to various special situations in which the transition between Kasner epochs by the rule eq. 37 becomes incorrect (including the situations described above). These "dangerous" situations could break the laws used for the statistical analysis below. As mentioned, however, the probability for such deviations converges asymptotically to zero; this issue will be discussed below.

Consider an era that contains k Kasner epochs with a parameter u running through the values

${\displaystyle u_{n}=k+x-1-n,\quad n=0,1,\cdots ,k-1,}$

(eq. 65)

and let α and β are the oscillating functions during this era (Fig. 4).[60]

Initial moments of Kasner epochs with parameters un are Ωn. In each initial moment, one of the values α or β is zero, while the other has a minimum. Values α or β in consecutive minima, that is, in moments Ωn are

${\displaystyle \alpha _{n}=-\delta _{n}\Omega _{n}\,}$

(eq. 66)

(not distinguishing minima α and β). Values δn that measure those minima in respective Ωn units can run between 0 and 1. Function γ monotonously decreases during this era; according to eq. 63 its value in moment Ωn is

${\displaystyle \gamma _{n}=-\Omega _{n}(1-\delta _{n}).\,}$

(eq. 67)

During the epoch starting at moment Ωn and ending at moment Ωn+1 one of the functions α or β increases from −δnΩn to zero while the other decreases from 0 to −δn+1Ωn+1 by linear laws, respectively:

${\displaystyle \mathrm {const} +|p_{1}(u_{n})|\Omega \,}$  and ${\displaystyle \mathrm {const} -p_{2}(u_{n})\Omega \,}$

resulting in the recurrence relation

${\displaystyle \delta _{n+1}\Omega _{n+1}={\frac {1+u_{n}}{u_{n}}}\delta _{n}\Omega _{n}={\frac {1+u_{0}}{u_{n}}}\delta _{0}\Omega _{0}}$

(eq. 68)

and for the logarithmic epoch length

${\displaystyle \Delta _{n+1}\equiv \Omega _{n+1}-\Omega _{n}={\frac {f(u_{n})}{u_{n}}}\delta _{n}\Omega _{n}={\frac {f(u_{n})(1+u_{n-1})}{f(u_{n-1})u_{n}}}\Delta _{n},}$

(eq. 69)

where, for short, f(u) = 1 + u + u2. The sum of n epoch lengths is obtained by the formula

${\displaystyle \Omega _{n}-\Omega _{0}=\left[n(n-1)+{\frac {nf(u_{n-1})}{u_{n-1}}}\right]\delta _{0}\Omega _{0}.}$

(eq. 70)

It can be seen from eq. 68 that |αn+1| > |αn|, i.e., the oscillation amplitudes of functions α and β increase during the whole era although the factors δn may be small. If the minimum at the beginning of an era is deep, the next minima will not become shallower; in other words, the residue |α — β| at the moment of transition between Kasner epochs remains large. This assertion does not depend upon era length k because transitions between epochs are determined by the common rule eq. 37 also for long eras.

The last oscillation amplitude of functions α or β in a given era is related to the amplitude of the first oscillation by the relationship |αk−1| = |α0| (k + x) / (1 + x). Even at k 's as small as several units x can be ignored in comparison to k so that the increase of α and β oscillation amplitudes becomes proportional to the era length. For functions a = eα and b = eβ this means that if the amplitude of their oscillations in the beginning of an era was A0, at the end of this era the amplitude will become ${\displaystyle A_{0}^{k/(1+x)}}$ .

The length of Kasner epochs (in logarithmic time) also increases inside a given era; it is easy to calculate from eq. 69 that Δn+1 > Δn.[61] The total era length is

${\displaystyle \Omega _{0}^{\prime }-\Omega _{0}\equiv \Omega _{k}-\Omega _{0}=k\left(k+x+{\frac {1}{x}}\right)\delta _{0}\Omega _{0}}$

(eq. 71)

(the term with 1/x arises from the last, k-th, epoch whose length is great at small x; cf. Fig. 2). Moment Ωn when the k-th epoch of a given era ends is at the same time moment Ω'0 of the beginning of the next era.

In the first Kasner epoch of the new era function γ is the first to rise from the minimal value γk = − Ωk (1 − δk) that it reached in the previous era; this value plays the role of a starting amplitude δ'0Ω'0 for the new series of oscillations. It is easily obtained that:

${\displaystyle \delta _{0}^{\prime }\Omega _{0}^{\prime }=\left(\delta _{0}^{-1}+k^{2}+kx-1\right)\delta _{0}\Omega _{0}.}$

(eq. 72)

It is obvious that δ'0Ω'0 > δ0Ω0. Even at not very great k the amplitude increase is very significant: function c = eγ begins to oscillate from amplitude ${\displaystyle A_{0}'\sim A_{0}^{k^{2}}}$ . The issue about the abovementioned "dangerous" cases of drastic lowering of the upper oscillation limit is left aside for now.

According to eq. 40 the increase in matter density during the first (k − 1) epochs is given by the formula

${\displaystyle \ln \left({\frac {\varepsilon _{n+1}}{\varepsilon _{n}}}\right)=2\left[1-p_{3}(u_{n})\right]\Delta _{n+1}.}$

For the last k epoch of a given era, it should be taken into account that at u = x < 1 the greatest power is p2(x) (not p3(x) ). Therefore, for the density increase over the whole era one obtains

${\displaystyle \ln \left({\frac {\varepsilon _{k}}{\varepsilon _{0}}}\right)\equiv \ln \left({\frac {\varepsilon _{0}'}{\varepsilon _{0}}}\right)=2(k-1+x)\delta _{0}\Omega _{0}.}$

(eq. 73)

Therefore, even at not very great k values, ${\displaystyle \varepsilon _{0}'/\varepsilon _{0}\sim A_{0}^{2k}}$ . During the next era (with a length k ' ) density will increase faster because of the increased starting amplitude A0': ${\displaystyle \varepsilon _{0}''/\varepsilon _{0}'\sim A_{0}'^{2k''}\sim A_{0}^{2k^{2}k'}}$ , etc. These formulae illustrate the steep increase in matter density.

### Statistical analysis near the singularity

The sequence of era lengths k(s), measured by the number of Kasner epochs contained in them, acquires asymptotically the character of a random process. The same pertains also to the sequence of the interchanges of the pairs of oscillating functions on going over from one era to the next (it depends on whether the numbers k(s) are even or odd). A source of this stochasticity is the rule eqs. 4142 according to which the transition from one era to the next is determined in an infinite numerical sequence of u values. This rule states, in other words, that if the entire infinite sequence begins with a certain initial value ${\displaystyle u_{\max }^{(0)}=k^{(0)}+x^{(0)}}$ , then the lengths of the eras k(0), k(1), ..., are the numbers in the continued fraction expansion

${\displaystyle k^{(0)}+x^{(0)}=k^{(0)}+{\frac {1}{k^{(1)}+{\frac {1}{k^{(2)}+\dots }}}}.}$

(eq. 73a)

This expansion corresponds to the mapping transformation of the interval [0, 1] onto itself by the formula Tx = {1/x}, i.e., xs+1 = {1/xs}. This transformation belongs to the so-called expanding transformations of the interval [0, 1], i.e., transformations xf(x) with |f′(x)| > 1. Such transformations possess the property of exponential instability: if we take initially two close points their mutual distance increases exponentially under the iterations of the transformations. It is well known that the exponential instability leads to the appearance of strong stochastic properties.

It is possible to change over to a probabilistic description of such a sequence by considering not a definite initial value x(0) but the values x(0) = x distributed in the interval from 0 to 1 in accordance with a certain probabilistic distributional law w0(x). Then the values of x(s) terminating each era will also have distributions that follow certain laws ws(x). Let ws(x)dx be the probability that the s-th era terminates with the value ${\displaystyle u_{\max }^{(s)}=x}$  lying in a specified interval dx.

The value x(s) = x, which terminates the s-th era, can result from initial (for this era) values ${\displaystyle u_{\max }^{(s)}=x+k}$ , where k = 1, 2, ...; these values of ${\displaystyle u_{\max }^{(s)}}$  correspond to the values x(s–1) = 1/(k + x) for the preceding era. Noting this, one can write the following recurrence relation, which expresses the distribution of the probabilities ws(x) in terms of the distribution ws–1(x):

${\displaystyle w_{s}(x)dx=\sum _{k=1}^{\infty }w_{s-1}\left({\frac {1}{k+x}}\right)\left\vert d{\frac {1}{k+x}}\right\vert }$

or

${\displaystyle w_{s}(x)=\sum _{k=1}^{\infty }{\frac {1}{\left(k+x\right)^{2}}}w_{s-1}\left({\frac {1}{k+x}}\right).}$

(eq. 73c)

If the distribution ws(x) tends with increasing s to a stationary (independent of s) limiting distribution w(x), then the latter should satisfy an equation obtained from eq. 73c by dropping the indices of the functions ws−1(x) and ws(x). This equation has a solution

${\displaystyle w(x)={\frac {1}{\left(1+x\right)\ln 2}}}$

(eq. 74)

(normalized to unity and taken to the first order of x).[62]

In order for the s-th era to have a length k, the preceding era must terminate with a number x in the interval between 1/(k + 1) and 1/k. Therefore, the probability that the era will have a length k is equal to (in the stationary limit)

${\displaystyle W(k)=\int \limits _{1/(k+1)}^{1/k}w(x)\,dx={\frac {1}{\ln 2}}\ln {\frac {(k+1)^{2}}{k(k+2)}}.}$

(eq. 75)

At large values of k

${\displaystyle W(k)\approx {\frac {1}{k^{2}\ln 2}}.}$

(eq. 76)

In relating the statistical properties of the cosmological model with the ergodic properties of the transformation xs+1 = {1/xs} an important point must be mentioned. In an infinite sequence of numbers x constructed in accordance with this rule, there will be observed arbitrary small (but never vanishing) values of x and accordingly arbitrarily large lengths k. Such cases can (by no means necessarily!) give rise to certain specific situations when the notion of eras, as of sequences of Kasner epochs interchanging each other according to the rule eq. 37, loses its meaning (although the oscillatory mode of evolution of the model still persists). Such an "anomalous" situation can be manifested, for instance, in the necessity to retain in the right-hand side of eq. 26 terms not only with one of the functions a, b, c (say, a4), as is the case in the "regular" interchange of the Kasner epochs, but simultaneously with two of them (say, a4, b4, a2b2).

On emerging from an "anomalous" series of oscillations a succession of regular eras is restored. Statistical analysis of the behavior of the model which is entirely based on regular iterations of the transformations eq. 42 is corroborated by an important theorem: the probability of the appearance of anomalous cases tends asymptotically to zero as the number of iterations s → ∞ (i.e., the time t → 0) which is proved at the end of this section. The validity of this assertion is largely due to a very rapid rate of increase of the oscillation amplitudes during every era and especially in transition from one era to the next one.

The process of the relaxation of the cosmological model to the "stationary" statistical regime (with t → 0 starting from a given "initial instant") is less interesting, however, than the properties of this regime itself with due account taken for the concrete laws of the variation of the physical characteristics of the model during the successive eras.

An idea of the rate at which the stationary distribution sets in is obtained from the following example. Let the initial values x(0) be distributed in a narrow interval of width δx(0) about some definite number. From the recurrence relation eq. 73c (or directly from the expansion eq. 73a) it is easy to conclude that the widths of the distributions ws(x) (about other definite numbers) will then be equal to

${\displaystyle \delta x^{(s)}\approx \delta x^{(0)}\cdot k^{(1)2}k^{(2)2}\dots k^{(s)2}}$

(eq. 76a)

(this expression is valid only so long as it defines quantities δx(s) ≪ 1).

The mean value ${\displaystyle {\bar {k}}}$ , calculated from this distribution, diverges logarithmically. For a sequence, cut off at a very large, but still finite number N, one has ${\displaystyle {\bar {k}}\sim \ln N}$ . The usefulness of the mean in this case is very limited because of its instability: because of the slow decrease of W(k), fluctuations in k diverge faster than its mean. A more adequate characteristic of this sequence is the probability that a randomly chosen number from it belongs to an era of length K where K is large. This probability is lnK / lnN. It is small if ${\displaystyle 1\ll K\ll N}$ . In this respect one can say that a randomly chosen number from the given sequence belongs to the long era with a high probability.

It convenient to average expressions that depend simultaneously on k(s) and x(s). Since both these quantities are derived from the same quantity x(s–1) (which terminates the preceding era), in accordance with the formula k(s) + x(s) = 1/x(s–1), their statistical distributions cannot be regarded as independent. The joint distribution Ws(k,x)dx of both quantities can be obtained from the distribution ws–1(x)dx by making in the latter the substitution x → 1/(x + k). In other words, the function Ws(k,x) is given by the very expression under the summation sign in the right side of eq. 73c. In the stationary limit, taking w from eq. 74, one obtains

${\displaystyle W(k,x)={\frac {k+x+1}{\left(k+x\right)\ln 2}}.}$

(eq. 76b)

Summation of this distribution over k brings us back to eq. 74, and integration with respect to dx to eq. 75.

The recurrent formulas defining transitions between eras are re-written with index s numbering the successive eras (not the Kasner epochs in a given era!), beginning from some era (s = 0) defined as initial. Ω(s) and ε(s) are, respectively, the initial moment and initial matter density in the s-th era; δ(s)Ω(s) is the initial oscillation amplitude of that pair of functions α, β, γ, which oscillates in the given era: k(s) is the length of s-th era, and x(s) determines the length (number of Kasner epochs) of the next era according to k(s+1) = [1/x(s)]. According to eqs. 7173

${\displaystyle \Omega ^{(s+1)}/\Omega ^{(s)}=1+\delta ^{(s)}k^{(s)}\left(k^{(s)}+x^{(s)}+{\frac {1}{x^{(s)}}}\right)\equiv \exp \xi ^{(s)},}$

(eq. 77)

${\displaystyle \delta ^{(s+1)}=1-{\frac {\left(k^{(s)}/x^{(s)}+1\right)\delta ^{(s)}}{1+\delta ^{(s)}k^{(s)}\left(1+x^{(s)}+1/x^{(s)}\right)}},}$

(eq. 78)

${\displaystyle \ln \left({\frac {\varepsilon ^{(s+1)}}{\varepsilon ^{(s)}}}\right)=2\left(k^{(s)}+x^{(s)}-1\right)\delta ^{(s)}\Omega ^{(s)}}$

(eq. 79)

(s) is introduced in eq. 77 to be used further on).

The quantities δ(s) have a stable stationary statistical distribution P(δ) and a stable (small relative fluctuations) mean value. For their determination BKL used [57] (with due reservations) an approximate method based on the assumption of statistical independence of the random quantity δ(s) and of the random quantities k(s), x(s). For the function P(δ) an integral equation was set up which expressed the fact that the quantities δ(s+1) and δ(s) interconnected by the relation eq. 78 have the same distribution; this equation was solved numerically. In a later work,[63] KL et al. showed that the distribution P(δ) can actually be found exactly by an analytical method.

For statistical properties in the stationary limit, it is reasonable to introduce the so-called natural extension of the transformation Tx = {1/x} by continuing it without limit to negative indices. Otherwise stated, this is a transition from a one-sided infinite sequence of the numbers (x0, x1, x2, ...), connected by the equalities Tx = {1/x}, to a "doubly infinite" sequence X = (..., x−1, x0, x1, x2, ...) of the numbers which are connected by the same equalities for all –∞ < s < ∞. Of course, such expansion is not unique in the literal meaning of the word (since xs–1 is not determined uniquely by xs), but all statistical properties of the extended sequence are uniform over its entire length, i.e., are invariant with respect to arbitrary shift (and x0 loses its meaning of an "initial" condition). The sequence X is equivalent to a sequence of integers K = (..., k−1, k0, k1, k2, ...), constructed by the rule ks = [1/xs–1]. Inversely, every number of X is determined by the integers of K as an infinite continuous fraction

${\displaystyle x_{s}={\frac {1}{k_{s+1}+{\frac {1}{k_{s+2}+\dots }}}}\equiv x_{s+1}^{+}}$

(eq. 79a)

(the convenience of introducing the notation ${\displaystyle x_{s+1}^{+}}$  with an index shifted by 1 will become clear in the following). For concise notation the continuous fraction is denoted simply by enumeration (in square brackets) of its denominators; then the definition of ${\displaystyle x_{s}^{+}}$  can be written as

${\displaystyle x_{s}^{+}=\left[k_{s},k_{s+1},\dots \right]}$

(eq. 79b)

Reverse quantities are defined by a continuous fraction with a retrograde (in the direction of diminishing indices) sequence of denominators

${\displaystyle x_{s}^{-}=\left[k_{s-1},k_{s-2},\dots \right]}$

(eq. 79c)

The recurrence relation eq. 78 is transformed by introducing temporarily the notation ηs = (1 − δs)/δs. Then eq. 78 can be rewritten as

${\displaystyle \eta _{s+1}={\frac {1}{\eta _{s}x_{s-1}+k_{s}}}}$

By iteration an infinite continuous fraction is obtained

${\displaystyle \eta _{s+1}x_{s}=\left[k_{s},k_{s-1},\dots \right]=x_{s+1}^{-}}$

Hence ${\displaystyle \eta _{s}=x_{s}^{-}/x_{s}^{+}}$  and finally

${\displaystyle \delta _{s}={\frac {x_{s}^{+}}{x_{s}^{+}+x_{s}^{-}}}}$

(eq. 79d)

This expression for δs contains only two (instead of the three in [57]) random quantities ${\displaystyle x_{s}^{+}}$  and ${\displaystyle x_{s}^{-}}$ , each of which assumes values in the interval [0, 1].

It follows from the definition eq. 79c that ${\displaystyle 1/x_{s}^{-}=x_{s}^{-}+k_{s}=x_{s}^{-}+\left[1/x_{s}^{+}\right]}$ . Hence the shift of the entire sequence X by one step to the right means a joint transformation of the quantities ${\displaystyle x_{s}^{+}}$  and ${\displaystyle x_{s}^{-}}$  according to

${\displaystyle x_{s+1}^{+}={\frac {1}{x_{s}^{+}}},\quad x_{s+1}^{-}={\frac {1}{{\frac {1}{x_{s}^{+}}}+x_{s}^{-}}}}$

(eq. 79e)

This is a one-to-one mapping in the unit square. Thus we have now a one-to-one transformation of two quantities instead of a not one-to-one transformation Tx = {1/x} of one quantity.

The quantities ${\displaystyle x_{s}^{+}}$  and ${\displaystyle x_{s}^{-}}$  have a joint stationary distribution P(x+, x). Since eq. 79e is a one-to-one transformation, the condition for the distribution to be stationary is expressed simply by a function equation

${\displaystyle P\left(x_{s}^{+},x_{s}^{-}\right)=P\left(x_{s+1}^{+},x_{s+1}^{-}\right)J}$

(eq. 79f)

where J is the Jacobian of the transformation.

A shift of the sequence X by one step gives rise to the following transformation T of the unit square:

${\displaystyle x^{\prime }={\frac {1}{x}},\quad y^{\prime }={\frac {1}{{\frac {1}{x}}+y}}}$

(with ${\displaystyle x\equiv x_{0}^{+}}$ , ${\displaystyle y\equiv x_{0}^{-}}$ , cf. eq. 79e). The density P(x, y) defines the invariant measure for this transformation. It is natural to suppose that P(x, y) is a symmetric function of x and y. This means that the measure is invariant with respect to the transformation S(x, y) = (y, x) and hence with respect to the product ST with ST(x, y) = (x″, y″) and

${\displaystyle x''={\frac {1}{{\frac {1}{x}}+y}},\quad y''={\frac {1}{x}}}$

Evidently ST has a first integral H = 1/x + y. On the line H = const ≡ c the transformation has the form

${\displaystyle {\frac {1}{x''}}=\left[{\frac {1}{x}}\right]+y=\left[{\frac {1}{x}}\right]+c-{\frac {1}{x}}=c-\left\{{\frac {1}{x}}\right\}}$

Hence the invariant measure density of ST must be of the form

${\displaystyle f(c)\ dc\ d{\frac {1}{x}}=f\left({\frac {1}{x}}+y\right){\frac {1}{x^{2}}}dx\ dy}$

With the account taken of the symmetry P(x, y)= P(y, x) this becomes f(c)= c−2 and hence (after normalization)

${\displaystyle P\left(x_{s}^{+},x_{s}^{-}\right)={\frac {1}{\left(1+x^{+}x^{-}\right)^{2}\ln 2}}}$

(eq. 79g)

(its integration over x+ or x yields the function w(x) eq. 74). The reduction of the transformation to the one-to-one mapping was used already by [64] and they obtained a formula of the form of eq. 79g but for other variables; their paper does not contain applications to the problems which are considered in.[63]

The correctness of eq. 79g be verified also by a direct calculation; the Jacobian of the transformation eq. 79e is

${\displaystyle J={\frac {\partial \left(x_{s+1}^{+},x_{s+1}^{-}\right)}{\partial \left(x_{s}^{+},x_{s}^{-}\right)}}={\frac {\partial x_{s+1}^{+}}{\partial x_{s}^{+}}}{\frac {\partial x_{s+1}^{-}}{\partial x_{s}^{-}}}=\left({\frac {x_{s+1}^{+}}{x_{s}^{+}}}\right)^{2}}$

(in its calculation one must note that ${\displaystyle \left[1/x_{s}^{+}\right]+\left\{1/x_{s}^{+}\right\}=1/x_{s}^{+}}$ ).

Figure 5. The probability distribution function P(δ). Red line: exact function eq. 79h. Blue line: approximate solution of the integral equation in.[57] Both curves appear to be strikingly similar and the means of both distributions are 0.50.[65]

Since by eq. 79d δs is expressed in terms of the random quantities x+ and x, the knowledge of their joint distribution makes it possible to calculate the statistical distribution P(δ) by integrating P(x+, x) over one of the variables at a constant value of δ. Due to symmetry of the function eq. 79g with respect to the variables x+ and x, P(δ) = P(1 − δ), i.e., the function P(δ) is symmetrical with respect to the point δ = 1/2. Then

${\displaystyle P(\delta )\ d\delta =d\delta \int _{0}^{1}P\left(x^{+},{\frac {x^{+}\delta }{1-\delta }}\right)\left({\frac {\partial x^{-}}{\partial \delta }}\right)_{x^{+}}dx^{+}}$

On evaluating this integral (for 0 ≤ δ ≤ 1/2 and then making use of the aforementioned symmetry), finally

${\displaystyle P(\delta )={\frac {1}{\left(|1-2\delta |+1\right)\ln 2}}}$

(eq. 79h)

The mean value ${\displaystyle {\bar {\delta }}}$ = 1/2 already as a result of the symmetry of the function P(δ). Thus the mean value of the initial (in every era) amplitude of oscillations of the functions α, β, γ increases as Ω/2.

The statistical relation between large time intervals Ω and the number of eras s contained in them is found by repeated application of eq. 77:

${\displaystyle {\frac {\Omega ^{(s)}}{\Omega ^{(0)}}}=\exp \left(\sum _{p=0}^{s-1}\xi ^{(p)}\right).}$

(eq. 80)

Direct averaging of this equation, however, does not make sense: because of the slow decrease of function W(k) eq. 76, the average values of the quantity exp ξ(s) are unstable in the above sense – the fluctuations increase even more rapidly than the mean value itself with increasing region of averaging. This instability is eliminated by taking the logarithm: the "doubly-logarithmic" time interval

${\displaystyle \tau _{s}\equiv \ln \left({\frac {\Omega ^{(s)}}{\Omega ^{(0)}}}\right)=\ln |\ln t_{s}|-\ln |\ln t_{0}|=\sum _{p=0}^{s-1}\xi ^{(p)}}$

(eq. 81)

is expressed by the sum of quantities ξ(p) which have a stable statistical distribution. The mean value of τ is ${\displaystyle {\bar {\tau }}=s{\bar {\xi }}}$ . To calculate ${\displaystyle {\bar {\xi }}}$  note that eq. 77 can be rewritten as

${\displaystyle \xi _{s}=\ln {\frac {\left(k_{s}+x_{s}\right)\delta _{s}}{x_{s}\left(1-\delta _{s+1}\right)}}=\ln {\frac {\delta _{s}}{x_{s}x_{s-1}\left(1-\delta _{s+1}\right)}}}$

(eq. 81a)

For the stationary distribution ${\displaystyle {\overline {\ln x_{s}}}={\overline {\ln x_{s-1}}}}$ , and in virtue of the symmetry of the function P(δ) also ${\displaystyle {\overline {\ln \delta _{s}}}={\overline {\ln \left(\delta _{s+1}\right)}}}$ . Hence

${\displaystyle {\bar {\xi }}=-2{\overline {\ln x}}=-2\int _{0}^{1}w(x)\ln x\ dx={\frac {\pi ^{2}}{6\ln 2}}=2.37}$

(w(x) from eq. 74). Thus

${\displaystyle {\overline {\tau _{s}}}=2.37s,}$

(eq. 82)

which determines the mean doubly-logarithmic time interval containing s successive eras.

For large s the number of terms in the sum eq. 81 is large and according to general theorems of the ergodic theory the values of τs are distributed around ${\displaystyle {\overline {\tau _{s}}}}$  according to Gauss' law with the density

${\displaystyle \rho (\tau _{s})=\left(2\pi D_{\tau }\right)^{-1/2}\exp \left[-\left(\tau _{s}-{\overline {\tau _{s}}}\right)^{2}/2D_{\tau }\right]}$

(eq. 82a)

Calculation of the variance Dτ is more complicated since not only the knowledge of ${\displaystyle {\bar {\xi }}}$  and ${\displaystyle {\overline {\xi ^{2}}}}$  are needed but also of the correlations ${\displaystyle {\overline {\xi _{p}\xi _{p\prime }}}}$ . The calculation can be simplified by rearranging the terms in the sum eq. 81. By using eq. 81a the sum can be rewritten as

${\displaystyle \sum _{p=1}^{s}\xi _{p}=\ln \prod _{p=1}^{s}{\frac {\delta _{p}}{\left(1-\delta _{p+1}\right)x_{p}x_{p-1}}}=\ln \prod _{p=1}^{s}{\frac {\delta _{p}}{\left(1-\delta _{p}\right)x_{p-1}^{2}}}+\ln {\frac {x_{0}}{x_{s}}}+\ln {\frac {1-\delta _{1}}{1-\delta _{s+1}}}}$

The last two terms do not increase with increasing s; these terms can be omitted as the limiting laws for large s are dominating. Then

${\displaystyle \sum _{p=1}^{s}\xi _{p}=\sum _{p=1}^{s}\ln {\frac {x_{p}^{-}}{x_{p}^{+}}}}$

(eq. 82b)

(the expression eq. 79d for δp is taken into account). To the same accuracy (i.e., up to the terms which do not increase with s) the equality

${\displaystyle \sum _{p=1}^{s}\xi _{p}^{+}=\sum _{p=1}^{s}\xi _{p}^{-}}$

(eq. 82c)

is valid. Indeed, in virtue of eq. 79e

${\displaystyle x_{p+1}^{+}+{\frac {1}{x_{p+1}^{-}}}={\frac {1}{x_{p}^{+}}}+x_{p}^{-}}$

and hence

${\displaystyle \ln \left(1+x_{p+1}^{+}x_{p+1}^{-}\right)-\ln x_{p+1}^{-}=\ln \left(1+x_{p}^{+}x_{p}^{-}\right)-\ln x_{p}^{+}}$

By summing this identity over p eq. 82c is obtained. Finally again with the same accuracy ${\displaystyle x_{p}^{+}}$  is changed for xp under the summation sign and thus represent τs as

${\displaystyle \tau _{s}=\sum _{p=1}^{\infty }\eta _{p},\quad \eta _{p}=-2\ln x_{p}}$

(eq. 83)

The variance of this sum in the limit of large s is

${\displaystyle D_{{\tau }_{s}}={\overline {\left(\tau _{s}-{\overline {\tau _{s}}}\right)^{2}}}\approx s\left\{{\overline {\eta ^{2}}}-{\bar {\eta }}^{2}+2\sum _{p=1}^{\infty }\left({\overline {\eta _{0}\eta _{p}}}-{\bar {\eta }}^{2}\right)\right\}}$

(eq. 84)

It is taken into account that in virtue of the statistical homogeneity of the sequence X the correlations ${\displaystyle {\overline {\eta _{p}\eta _{p\prime }}}}$  depend only on the differences |pp′|. The mean value ${\displaystyle {\bar {\eta }}={\bar {\xi }}}$ ; the mean square

${\displaystyle {\overline {\eta ^{2}}}=4\int _{0}^{1}w(x)\ln ^{2}x\ dx={\frac {6\xi (3)}{\ln 2}}=10.40}$

By taking into account also the values of correlations ${\displaystyle {\overline {\eta _{0}\eta _{p}}}}$  with p = 1, 2, 3 (calculated numerically) the final result Dτs = (3.5 ± 0.1)s is obtained.

With increasing s the relative fluctuation ${\displaystyle D_{{\tau }_{s}}/{\overline {\tau _{s}}}}$  tends to zero as s−1/2. In other words, the statistical relation eq. 82 becomes almost certain at large s. This makes it possible to invert the relation, i.e., to represent it as the dependence of the average number of the eras sτ that are interchanged in a given interval τ of the double logarithmic time:

${\displaystyle {\overline {s_{\tau }}}=0.47\tau .}$

(eq. 85)

The statistical distribution of the exact values of sτ around its average is also Gaussian with the variance

${\displaystyle D_{s_{\tau }}=3.5{\frac {{\overline {s_{\tau }}}^{3}}{\tau ^{2}}}=0.26\tau }$

The respective statistical distribution is given by the same Gaussian distribution in which the random variable is now sτ at a given τ:

${\displaystyle \rho (s_{\tau })\propto \exp \left\{-\left(s_{\tau }-0.47\tau \right)^{2}/0.43\tau \right\}.}$

(eq. 86)

From this point of view, the source of the statistical behavior is the arbitrariness in the choice of the starting point of the interval τ superimposed on the infinite sequence of the interchanging eras.

Respective to matter density, eq. 79 can be re-written with account of eq. 80 in the form

${\displaystyle \ln \ln {\frac {\varepsilon ^{(s+1)}}{\varepsilon ^{(s)}}}=\eta _{s}+\sum _{p=0}^{s-1}\xi _{p},\quad \eta _{s}=\ln \left[2\delta ^{(s)}\left(k^{(s)}+x^{(s)}-1\right)\Omega ^{(0)}\right]}$

and then, for the total energy change during s eras,

${\displaystyle \ln \ln {\frac {\varepsilon ^{(s)}}{\varepsilon ^{(0)}}}=\ln \sum _{p=0}^{s-1}\exp \left\{\sum _{q=0}^{p}\xi _{q}+\eta _{p}\right\}.}$

(eq. 87)

The term with the sum by p gives the main contribution to this expression because it contains an exponent with a large power. Leaving only this term and averaging eq. 87, one gets in its right hand side the expression ${\displaystyle s{\bar {\xi }}}$  which coincides with eq. 82; all other terms in the sum (also terms with ηs in their powers) lead only to corrections of a relative order 1/s. Therefore,

${\displaystyle {\overline {\ln \ln \left({\frac {\varepsilon ^{(s)}}{\varepsilon ^{(0)}}}\right)}}={\overline {\ln \left({\frac {\Omega ^{(s)}}{\Omega ^{(0)}}}\right)}}.}$

(eq. 88)

By virtue of the almost certain character of the relation between τs and s eq. 88 can be written as

${\displaystyle {\overline {\ln \ln \left(\varepsilon _{\tau }/\varepsilon ^{(0)}\right)}}=\tau \quad {\text{or}}\quad {\overline {\ln \ln \left(\varepsilon ^{(s)}/\varepsilon ^{(0)}\right)}}=2.1s,}$

which determines the value of the double logarithm of density increase averaged by given double-logarithmic time intervals τ or by a given number of eras s.

These stable statistical relationships exist specifically for doubly-logarithmic time intervals and for the density increase. For other characteristics, e.g., ln (ε(s)(0)) or Ω(s) / Ω(0) = exp τs the relative fluctuation increase exponentially with the increase of the averaging range thereby voiding the term mean value of a stable meaning.

The origin of the statistical relationship eq. 88 can be traced already from the initial law governing the variation of the density during the individual Kasner epochs. According to eq. 21, during the entire evolution we have

${\displaystyle \ln \ln \varepsilon (t)={\text{const}}+\ln \Omega +\ln 2(1-p_{3}(t)),}$

with 1 − p3(t) changing from epoch to epoch, running through values in the interval from 0 to 1. The term ln Ω = ln ln (1/t) increases monotonically; on the other hand, the term ln2(1 − p3) can assume large values (comparable with ln Ω) only when values of p3 very close to unity appear (i.e., very small |p1|). These are precisely the "dangerous" cases that disturb the regular course of evolution expressed by the recurrent relationships eqs. 7779.

It remains to show that such cases actually do not arise that in the asymptotic limiting regime. The spontaneous evolution of the model starts with a certain instant at which definite initial conditions are specified in an arbitrary manner. Accordingly, by "asymptotic" is meant a regime sufficiently far away from the chosen initial instant.

Dangerous cases are those in which excessively small values of the parameter u = x (and hence also |p1| ≈ x) appear at the end of an era. A criterion for selection of such cases is the inequality

${\displaystyle x^{(s)}\exp \left|\alpha ^{(s)}\right|\lesssim 1,}$

(eq. 89)

where | α(s) | is the initial minima depth of the functions that oscillate in era s (it would be more appropriate to choose the final amplitude, but that would only strengthen the selection criterion).

The value of x(0) in the first era is determined by the initial conditions. Dangerous are values in the interval δx(0) ~ exp ( − |α(0)| ), and also in intervals that could result in dangerous cases in the next eras. In order for x(s) to fall in the dangerous interval δx(s) ~ exp ( − | α(s) | ), the initial value x(0) should lie into an interval of a width δx(0) ~ δx(s) / k(1)^2 ... k(s)^2.[66] Therefore, from a unit interval of all possible values of x(0), dangerous cases will appear in parts λ of this interval:

${\displaystyle \lambda =\exp \left(\left|-\alpha ^{(s)}\right|\right)+\sum _{s=1}^{\infty }\sum _{k}{\frac {\exp \left(\left|-\alpha ^{(s)}\right|\right)}{k^{(1)^{2}}k^{(2)^{2}}\dots k^{(s)^{2}}}}}$

(eq. 90)

(the inner sum is taken over all the values k(1), k(2), ... , k(s) from 1 to ∞). It is easy to show that this era converges to the value λ ${\displaystyle \ll }$  1 whose order of magnitude is determined by the first term in eq. 90. This can be shown by a strong majoration of the era for which one substitutes | α(s) | = (s + 1) | α(0) |, regardless of the lengths of eras k(1), k(2), ... (In fact | α(s) | increase much faster; even in the most unfavorable case k(1) = k(2) = ... = 1 values of | α(s) | increase as qs | α(0) | with q > 1.) Noting that

${\displaystyle \sum _{k}{\frac {1}{k^{(1)^{2}}k^{(2)^{2}}\dots k^{(s)^{2}}}}=\left(\pi ^{2}/6\right)^{s}}$

one obtains

${\displaystyle \lambda =\exp \left(\left|-\alpha ^{(0)}\right|\right)\sum _{s=0}^{\infty }\left[\left(\pi ^{2}/6\right)\exp \left(\left|-\alpha ^{(0)}\right|\right)\right]^{s}\approx \exp \left(\left|-\alpha ^{(0)}\right|\right).}$

If the initial value of x(0) lies outside the dangerous region λ there will be no dangerous cases. If it lies inside this region dangerous cases occur, but upon their completion the model resumes a "regular" evolution with a new initial value which only occasionally (with a probability λ) may come into the dangerous interval. Repeated dangerous cases occur with probabilities λ2, λ3, ... , asymptotically converging to zero.

## General solution with small oscillations

In the above models, metric evolution near the singularity is studied on the example of homogeneous space metrics. It is clear from the characteristic of this evolution that the analytic construction of the general solution for a singularity of such type should be made separately for each of the basic evolution components: for the Kasner epochs, for the process of transitions between epochs caused by "perturbations", for long eras with two perturbations acting simultaneously. During a Kasner epoch (i.e. at small perturbations), the metric is given by eq. 7 without the condition λ = 0.

BKL further developed a matter distribution-independent model (homogeneous or non-homogeneous) for long era with small oscillations. The time dependence of this solution turns out to be very similar to that in the particular case of homogeneous models; the latter can be obtained from the distribution-independent model by a special choice of the arbitrary functions contained in it.[67]

It is convenient, however, to construct the general solution in a system of coordinates somewhat different from synchronous reference frame: g = 0 as in the synchronous frame, but instead of g00 = 1 it is now g00 = −g33. Defining again the space metric tensor γαβ = −gαβ one has, therefore

${\displaystyle g_{00}=\gamma _{33},\quad g_{0\alpha }=0.}$

(eq. 91)

The special space coordinate is written as x3 = z and the time coordinate is written as x0 = ξ (as different from proper time t); it will be shown that ξ corresponds to the same variable defined in homogeneous models. Differentiation by ξ and z is designated, respectively, by dot and prime. Latin indices a, b, c take values 1, 2, corresponding to space coordinates x1, x2 which will be also written as x, y. Therefore, the metric is

${\displaystyle ds^{2}=\gamma _{33}\left(d\xi ^{2}-dz^{2}\right)-\gamma _{ab}dx^{a}dx^{b}-2\gamma _{a3}dx^{a}dz.}$

(eq. 92)

The required solution should satisfy the inequalities

${\displaystyle \gamma _{33}\ll \gamma _{ab},}$

(eq. 93)

${\displaystyle \gamma _{a3}^{2}\ll \gamma _{aa}\gamma _{33}}$

(eq. 94)

(these conditions specify that one of the functions a2, b2, c2 is small compared to the other two which was also the case with homogeneous models).

Inequality eq. 94 means that components γa3 are small in the sense that at any ratio of the shifts dxa and dz, terms with products dxadz can be omitted in the square of the spatial length element dl2. Therefore, the first approximation to a solution is a metric eq. 92 with γa3 = 0:[68]

${\displaystyle ds^{2}=\gamma _{33}\left(d\xi ^{2}-dz^{2}\right)-\gamma _{ab}dx^{a}dx^{b}.}$

(eq. 95)

One can be easily convinced by calculating the Ricci tensor components ${\displaystyle R_{0}^{0}}$ , ${\displaystyle R_{3}^{0}}$ , ${\displaystyle R_{3}^{3}}$ , ${\displaystyle R_{a}^{b}}$  using metric eq. 95 and the condition eq. 93 that all terms containing derivatives by coordinates xa are small compared to terms with derivatives by ξ and z (their ratio is ~ γ33 / γab). In other words, to obtain the equations of the main approximation, γ33 and γab in eq. 95 should be differentiated as if they do not depend on xa. Designating

${\displaystyle \gamma _{33}=e^{\psi },\quad {\dot {\gamma }}_{ab}=\varkappa _{ab},\quad \gamma _{ab}^{\prime }=\lambda _{ab},\quad |\gamma _{ab}|=G^{2},}$

(eq. 96)

one obtains the following equations:[69]

${\displaystyle 2e^{\psi }R_{a}^{b}=G^{-1}\left(G\lambda _{a}^{b}\right)^{\prime }-G^{-1}\left(G\varkappa _{a}^{b}\right){\dot {}}=0,}$

(eq. 97)