On the dynamics of the general Bianchi IX spacetime near the singularity

We show that the complex dynamics of the general Bianchi IX universe in the vicinity of the spacelike singularity can be approximated by a simplified system of equations. Our analysis is mainly based on numerical simulations. The properties of the solution space can be studied by using this simplified dynamics. Our results will be useful for the quantization of the general Bianchi IX model.


I. INTRODUCTION
The problem of spacetime singularities is a central one in classical and quantum theories of gravity. Given some general conditions, it was proven that general relativity leads to singularities, among which special significance is attributed to big bang and black hole singularities [1].
The occurrence of a singularity in a physical theory usually signals the breakdown of that theory. In the case of general relativity, the expectation is that its singularities will disappear after quantization. Although a theory of quantum gravity is not yet available in finite form, various approaches exist within which the question of singularity avoidance can be addressed [2]. Quantum cosmological examples for such an avoidance can be found, for example, in [3][4][5][6] and the references therein.
Independent of the quantum fate of singularities, the question of their exact nature in the classical theory, and in particular for cosmology, is of considerable interest and has a long history; see, for example, [7] and [8] for recent reviews. This is also the topic of the present paper.
Already in the 1940s, Evgeny Lifshitz investigated the gravitational stability of non-stationary isotropic models of universes. He found that the isotropy of space cannot be retained in the evolution towards singularities [10] (see [11] for an extended physical interpretation). This motivated the activity of the Landau Institute in Moscow to examining the dynamics of homogeneous spacetimes [12]. A group of relativists inspired by Lev Landau, including Belinski, Khalatnikov and Lifshitz (BKL), started to investigating the dynamics of the Bianchi VIII and IX models near the initial spacelike cosmological singularity [13]. After several years, they found that the dynamical behaviour can be generalized to a generic solution of general relativity [14]. They did not present a mathematically rigorous proof, but rather a conjecture based on deep analytical insight. It is called the BKL conjecture (or the BKL scenario if specialized to the Bianchi-type IX model). The BKL conjecture is a locality conjecture stating that terms with temporal derivatives dominate over terms with spatial derivatives when approaching the singularity (with the exception of possible 'spikes' [8,9]). Consequently, points in space decouple and the dynamics then turn out to be effectively the same as those of the (non-diagonal) Bianchi IX universe. (In canonical gravity, this is referred to as the strong coupling limit, see e.g. [2], p. 127.) The dynamics of the Bianchi IX towards the singularity are characterized by an infinite number of oscillations, which give rise to a chaotic character of the solutions (see e.g [15]). Progress towards improving the mathematical rigour of the BKL conjecture has been made by several authors (see e.g. [16]), while numerical studies giving support to the conjecture have been performed (see e.g. [17]).
The dynamics of the diagonal Bianchi IX model, in the Hamiltonian formulation, were studied independently from BKL by Misner [18,19]. Misner's intention was to search for a possible solution to the horizon problem by a process that he called "mixing". 1 Ryan generalized Misner's formalism to the non-diagonal case in [20,21]. A qualitative treatment of the dynamics for all the Bianchi models may be found in the review article by Jantzen [22], and we make reference to it whenever we get similar results.
Part of the BKL conjecture is that non-gravitational ('matter') terms can be neglected when approaching the singularity. An important exception is the case of a massless scalar field, which has analogies with a stiff fluid (equation of state p = ρ) and, in Friedmann models, has the same dependence of the density on the scale factor as anisotropies (ρ ∝ a −6 ). As was rigorously shown in [23], such a scalar field will suppress the BKL oscillations and thus is relevant during the evolution towards the singularity. Arguments for the importance of stiff matter in the early universe were already given by Barrow [24].
In our present work, we shall mainly address the general (non-diagonal) Bianchi IX model near its singularity. Our main motivation is to provide support for a rather simple asymptotic form of the dynamics that can suitably model its exact complex dynamics. We expect this to be of relevance in the quantization of the general Bianchi IX model, which we plan to investigate in later papers; see, for example, [25]. Apart from a few particular solutions that form a set of measure zero in the solution space, no general analytic solutions to the classical equations of motion are known. Therefore, we will restrict ourselves to qualitative considerations which will be supported by numerical simulations. The examination of the non-diagonal dynamics presented in [26], though it is mathematically satisfactory, is based on the qualitative theory of differential equations, which is of little use for our purpose.
Our paper is organized as follows. Section II contains the formalism and presents our main results for a general Bianchi IX model. We first specify the kinematics and dynamics. We then consider a matter field in the form of (tilted) dust. This is followed by investigating the asymptotic regime of the dynamics near the singularity. Our conclusions are presented in Sec. III. The numerical methods used in our numerical simulations are described in the Appendix.

A. Kinematics
The general non-diagonal case describes a universe with rotating principal axes. The metric in a synchronous frame can be given as follows (see for the following e.g. [27] and [28]): where N is the lapse function. Spatial hypersurfaces in the spacetime are regarded topologically as S 3 (describing closed universes), which can be parametrized by using three angles θ ,φ,ψ ∈ [0, π] × [0, 2π] × [0, π]. The basis one-forms read σ 1 = − sin(ψ)dθ + cos(ψ) sin(θ)dφ , σ 2 = cos(ψ)dθ + sin(ψ) sin(θ)dφ , The σ i , i = 1, 2, 3, are dual to the vector fields which together with X 0 = ∂ t form an invariant basis of the Bianchi IX spacetime. The X i are constructed from the Killing vectors that generate the isometry group SO(3, R) (see [27] for more details). The basis one-forms satisfy the relation with C i jk = ε ijk being the structure constants of the Lie algebra so(3, R). The X i obey the algebra [X i , X j ] = −C k ij X k . We parametrize the metric coefficients in this frame as follows whereh The variables α, β + , and β − are known as the Misner variables. The scale factor exp(α) is related to the volume, while the anisotropy factors β + and β − describe the shape of this model universe. The variables Γ 1 , Γ 2 and Γ 3 were used by BKL in their original analysis [31].
The Euler angles θ, φ, and ψ are now dynamical quantities and describe nutation, precession, and pure rotation of the principal axes, respectively. In the case of Bianchi IX spacetime, the group SO(3, R) is the canonical choice for the diagonalization of the metric coefficients. For a treatment of other Bianchi models, see [22].

B. Dynamics
In the following, we shall discuss the Hamiltonian formulation of this model. In order to keep track of the diffeomorphism (momentum) constraints, we replace the metric (1) by the ansatz where N i are the shift functions. The Hamiltonian formulation was first derived in a series of papers by Ryan: the symmetric (non-tumbling) case obtained by constraining ψ, φ to be constant and keeping θ dynamical is discussed in [20], and the general case can be found in [21]. We write the Einstein-Hilbert action in the well known ADM form, where is the extrinsic curvature, and D i is the spatial covariant derivative in the noncoordinate basis {X i }. We will set 3 4πG σ 1 ∧ σ 2 ∧ σ 3 = 1 for simplicity. The threedimensional curvature (3) R on spatial hypersurfaces of constant coordinate time is given by We now turn to the calculation of the kinetic term and the diffeomorphism constraints. For this purpose, we define an antisymmetric angular velocity tensor ω i j by the matrix equation An explicit calculation of the right-hand side gives, using (7), The Lagrangian in the gauge N i = 0 then takes the form where the 'moments of inertia' are given by correspond to the rotational energy of a rigid body if the moments of inertia were constant. The canonical momenta conjugate to the Euler angles are given by It is convenient to introduce the following (non-canonical) angular momentum like variables: The relation to the canonical momenta can now explicitely be given by It is readily shown that the variables l i obey the Poisson bracket algebra {l i , l j } = −C k ij l k . After the usual Legendre transform, we obtain the Hamiltonian constraint, From (9), we find that the diffeomorphism constraints (∂L/∂N i = 0) can be written as where is the ADM momentum. From this expression we can finally compute the diffeomorphism constraints in terms of the angular momentum-like variables and obtain that is, we can identify the diffeomorphism constraints with a basis of the generators of SO(3, R). The full gravitational Hamiltonian then reads From the diffeomorphism constraints (22) we conclude that in the vacuum case l i = 0 and that therefore no rotation is possible, that is, we recover the diagonal case. If we want to obtain a Bianchi IX universe with rotating principal axes, we are thus forced to add matter to the system. A formalism for obtaining equations of motion for general Bianchi class A models filled with fluid matter was developed by Ryan [28]. For simplicity, we will only consider the case of dust as discussed by Kuchař and Brown in [29]. If we were, for example, interested in the study of the quantum version of this model, it would be desirable to introduce a fundamental matter field instead of an ideal fluid. Standard scalar fields alone cannot lead to a rotation for Bianchi IX models. The easiest way to achieve this is, to our knowledge, the introduction of a Dirac field [30].
C. Adding dust to the system The energy momentum tensor for dust reads T µν = ρu µ u ν . The local energy conservation ∇ µ T µν = 0 leads to a geodesic equation for the positions of the dust particles. Let us start therefore by considering the geodesic equation for a single dust particle, whose four-velocity we can express in the non-coordinate frame σ i used above by the Pfaffian form We partially fix the gauge by setting N i = 0. The normalization condition implies We have chosen here the minus sign because this guarantees that the proper time in the frame of the dust particle has the same orientation as the coordinate time t. The geodesic equation for the spatial components of the four velocity can then be written as The geodesic equation implies the existence of a constant of motion. To see this explicitly, we compute the expression i=1,2,3u i u i and convince ourselves that it vanishes identically. Thus the Euclidean sum is a constant of motion. Defining u ≡ (u 1 , u 2 , u 3 ) T , the geodesic equation (26) can be rewritten in vector notation, where "×" denotes the usual cross product in the three dimensional Euclidean space. Defining for convenience v ≡ O T u/C, we have (v 1 ) 2 + (v 2 ) 2 + (v 3 ) 2 = 1, and the geodesic equation simplifies to Note that we can also write ω It will thus be possible to eliminate ω from the geodesic equation by using the diffeomorphism constraints.
We now add homogeneous dust to the system. The formalism developed in [29] leads to the following form of the Hamiltonian and diffeomorphism constraints for dust in a Bianchi IX universe, where p T denotes the momentum canonically conjugate to T , where T is the global 'dust time'. Since the Hamiltonian does not explicitly depend on T , the momentum p T is a constant of motion. The fact that l 2 1 + l 2 2 + l 2 3 commutes with H implies that l 2 1 + l 2 2 + l 2 3 = (Cp T ) 2 is a conserved quantity. This is consistent with (27). We note that a similar form of the constraints (30) was already presented in [20,21]. The formalism is not entirely canonical and must be complemented by the geodesic equation (26).
For our numerical purposes, it will be convenient to rewrite the equations of motion in the variables Γ i introduced in (6). We find that the choice of the variables log Γ i allow for a better control over the error in the Hamiltonian constraint. Moreover, we pick the quasi-Gaussian gauge N = e 3α = √ Γ 1 Γ 2 Γ 3 , N i = 0. Recall that the singularity is reached in a finite amount of comoving time (corresponding to the gauge N = 1). The choice N = e 3α allows to resolve the oscillations in the approach towards the singularity. With these choices the Hamiltonian constraint becomes where the moments of inertia are The diffeomorphism constraints read l i = p T Cv i and can be used to eliminate the angular momentum variables from the equations of motion. These equations can then be written as where we have set p T ≡ 12p T for convenience. Note that these equations are exact.
(In [13], the matter terms were neglected.) We use the diffeomorphism constraints to eliminate ω from the geodesic equation (29). If expressed in the gauge N = e 3α and using the Γ i , the geodesic equation (29) can be written aṡ Together with the constraint v 2 1 + v 2 2 + v 2 3 = 1 this is all we need for numerical integration. Note that all dependence on the Euler angles and their momenta has dropped out from the equations of motion (33,34). The numerical method we use is described in Appendix A.

D. The tilted dust case
A qualitative picture for the dynamics of the universe can be obtained by considering the Hamiltonian (30) in the quasi-Gaussian gauge N = e 3α , N i = 0. The Hamiltonian can then be interpreted as the "relativistic energy" of a point particle (called the universe point) with "spacetime coordinates" (α, β + , β − ). The universe point is subject to the forces generated by the dynamical potential which is depicted in Fig. 1. The contour lines of the curvature potential − e 6α

12
(3) R are represented by solid black lines. The curvature potential is exponentially steep and takes its minimum at the origin β ± = 0. When the universe evolves towards the singularity (α → −∞), the curvature potential walls move away from the origin while becoming effectively hard walls in the vicinity of the singularity. The term e 3α p T 1 + h ij u i u j can be interpreted as three rotational potential walls. These potentials are rather unimportant in the examination of the asymptotic dynamics, since they move away from the origin with unit speed. The term I 3 can be interpreted as three singular centrifugal potential walls. They are represented by the dashed red lines. Asymptotically close to the singularity, these walls are expected to become static. In general, however, the centrifugal potential walls are dynamical and change in a complicated manner dictated by the geodesic equation (34). The centrifugal walls will prevent the universe point from penetrating certain regions of the configuration space. Misner [19] and Ryan [20,21] employed these facts to obtain approximate solutions in a diagrammatic form. The other Bianchi models can be treated in a similar way [22].

Special classes of solutions
Before doing the numerics, we comment on particular classes of solutions: one class of solutions is obtained if we choose, for example, the initial conditions v 1 = 0 = v 2 and v 3 = 1. The geodesic equation (34) implies now that the velocities stay constant in time. This implies that at all times l 1 = 0 = l 2 and l 3 = p ψ = p T C. This class of solutions is known as the non-tumbling case. Furthermore, there are classes of solutions which are rotating versions of the Taub solution. These solutions should be divided into two subclasses: one class that oscillates between the centrifugal walls and the curvature potential and one class that runs through the valley straight into the singularity. We set For the Γ i variables it means that Γ 1 = e 2α e 2β + = Γ 2 and Γ 3 = e 2α e −4β + . With this choice we obtain I 3 = 0 and 3I 1 = 3I 2 = sinh 2 (3β + ). Most importantly, the geodesic equation (34) is trivially satisfied, that is, v 1 = v 2 = 1/2 and v 3 = 0 for all times. When setting C = 0 we obtain the diagonal case which contains the isotropic case of a closed Friedmann universe. The simulation plotted in Fig. 2 was performed for the tumbling case, that is, the v i are chosen to be non-zero.

The asymptotic regime close to the singularity
In order to simplify the dynamics of the general case, BKL made two assumptions based on qualitative considerations of the equations of motion. The first assumption states that anisotropy of space grows without bound. This means that the solution enters the regime The ordering of indices is irrelevant. In fact, there are six possible orderings of indices which each correspond to the universe point being constrained to one of the six regions bounded by the rotation and centrifugal walls sketched in Fig. 1. The region Γ 1 > Γ 2 > Γ 3 corresponds to the right region above the line β − = 0 in Fig. 1. More precisely, the inequality (37) means that Our numerical simulations support the validity of this assumption (see the plot of the ratios Γ 2 /Γ 1 and Γ 3 /Γ 2 in Fig. 3). We provide plots of the two ratios Γ 2 /Γ 1 , Γ 3 /Γ 2 and the velocities v in order to provide a sanity check of the approximation we will perform later on. The second assumption made by BKL states that the Euler angles assume constant values: that is, the rotation of the principal axes stops for all practical purposes and the metric becomes effectively diagonal. The analysis of BKL [31] supports the consistency of making both assumptions at the same time. Similar heuristic considerations can possibly be applied to other Bianchi models as well [22]. In the dust model under consideration, this assumption is equivalent to the statement that the dust velocities v assume constant values v → v (0) . Our numerical results indicate that this is in fact true (see Fig. 3). BKL then arrive at the simplified effective set of equations.
Let us now carry out the approximation and apply it to our equations of motion.  Fig. 2. The peaks in the ratios appear during bounces of the universe point with the centrifugal walls. As we can see, the height of these peaks decreases in the evolution towards the singularity.
The kinetic term stays untouched during the approximation. The first step in the approximation is to ignore the rotational potential. In view of the strong inequality (37), we approximate the curvature potential via Furthermore, we approximate the centrifugal potential by Note that one centrifugal wall was ignored completely. Having Fig. 1 in mind, this approximation is well motivated since only two of the centrifugal walls are expected to have a significant influence on the dynamics of the universe point. After defining the new variables we arrive at a simplified Hamiltonian constraint and equations of motion, (log a) · (log b) · + (log a) · (log c) · + (log b) · (log c) · = a 2 + b/a + c/b, which coincides with the asymptotic form of equations obtained in [31]. Equations (43) can now be treated by the numerical methods which we have used in the previous sections. One must ensure that initial conditions are chosen such that the simulation starts close to the asymptotic regime (37).

III. CONCLUSIONS
The numerical simulations indicate that the non-diagonal Bianchi IX solutions, with tilted dust, evolve into the regime where Γ 1 Γ 2 Γ 3 and v i ≈const. The results motivate us to formulate the conjecture: Given a tumbling solution to the general Bianchi IX model filled with pressureless tilted matter, there exists t 0 ∈ R such that the solution is well approximated by a solution to the asymptotic equations of motion for all times t > t 0 describing the vicinity of the singularity.
To make the notion of "approximation" mathematically more precise, a suitable measure of the "distance" on the set of solutions is needed. For this purpose, we propose to use the following simple measure: where {Γ 1 , Γ 2 , Γ 3 } denotes the numerical solution to the exact equations of motion (33)- (34), and denote the numerical solution to the asymptotic equations of motion (43). We have evolved the exact system of equations from t = 0 forward in time until t = 3 × 10 6 . There we used the same initial conditions as the ones we used to obtain the solution shown in Fig. 2. We then took the final state at t = 3 × 10 6 as an initial condition for the asymptotic system of equations and evolved it backwards in time towards the re-bounce until t = −980. Fig. 4 presents the measure (44) as a function of time. We can see fast decrease of ∆ with increasing time (evolution towards the singularity) and fast increase of ∆ with decreasing time (evolution away from the singularity). Our numerical simulations give strong support to the conjecture concerning the asymptotic dynamics of the general Bianchi IX spacetime put forward long ago by Belinski, Khalatnikov, and Ryan [31]. We remark that approximating the diagonal or non-tumbling case by the asymptotic dynamics (43) is invalid (see [32] for more details).
It is sometimes stated that "matter does not matter" in the asymptotic regime, but this does not mean that one is allowed to use the dynamics of the purely diagonal case. One only encounters an effectively diagonal case, which is expressed in terms of the directional scale factors {a, b, c}. So there exist serious differences between the purely diagonal and effectively diagonal cases (see [32] for more details).
Employing the asymptotic form of the equations of motion may enable one to study the chaotic behaviour and other properties of the solution space for the general model. This is also important for quantizing the general Bianchi IX model, where the quantization of the exact dynamics seems to be quite difficult, whereas the quantization of the asymptotic case seems to be feasible [25].

ACKNOWLEDGMENTS
We are grateful to Vladimir Belinski and Claes Uggla for helpful discussions. This work was supported by the German-Polish bilateral project DAAD and MNiSW, No. 57391638, "Model of stellar collapse towards a singularity and its quantization".

Appendix A: Numerical analysis
Numerical simulations of the diagonal Bianchi IX model were already carried out in the late 1980s and early 1990s (see e.g. [33,34]; for a modern account, see [7]). Our main interest here is in the nondiagonal case. With given initial conditions (respecting the constraints H = 0), the system (33) and (34) can be integrated by using a suitable numerical method. In this work we employ the MATLAB R2016b solver ode113 [35]. This code is an implementation of an Adams-Bashforth-Moulton method. It turns out to lead to the best results when compared to other MATLAB solvers. The relative error tolerance of the solver was chosen to be of the order 10 −14 . We can integrate the equations of motion together with the geodesic equation to obtain a numerical solution to the system. We set up initial conditions at t = 0 and evolve the system forward in time towards the final singularity and away from the rebounce.
A major problem in numerical relativity is that the Hamiltonian constraint is not preserved exactly by the numerical procedure. Similar to [33,34], we find that the error in the Hamiltonian constraint varies strongest after the start of the simulation. Furthermore, it varies strongly when the evolution of the universe approaches the point of maximal expansion. While approaching the singularity, the error approaches an approximately constant value. This can be seen in Fig. 5. Therefore we can minimize the error when we choose the initial conditions far away from the point of maximal expansion. Moreover, it turned out that the error can be further reduced when constraining the solver's maximally allowed time step size. This time step size should, however, not be chosen too small since small time step sizes can drive the propagation of round off errors. Small step sizes are, of course, also numerically more expensive. By manually fine tuning the initial conditions and the maximally allowed time step size it was possible to keep the order of the error lower than 10 −15 .
Recall that the dynamics of Bianchi IX are chaotic, that is, slightly changing initial conditions have a large effect on the long time behaviour of solutions. Since the propagation of random numerical errors cannot be avoided, we will be dealing with a "butterfly effect" and it should in general not be expected that our numerical solution is an actual approximation of some exact solution of the equations of motion shows the speed of the universe point in the beta plane as measured in "α-time". This plot can be viewed as another check of the numerics. Between two successive bounces from the potential walls this quantity should be close to one. If this ceased to be true it would indicate that the error in the Hamiltonian constraint becomes relevant and cannot be neglected in the approach towards the singularity (see [7] for a more detailed discussion). Both plots correspond to the numerical solution shown in Fig. 2.
when considering large time intervals.