Curvature invariants for the Bianchi IX spacetime filled with tilted dust

We present an analysis of the Kretschmann and Weyl squared scalars for the general Bianchi IX model filled with tilted dust. Particular attention is given to the asymptotic regime close to the singularity for which we provide heuristic considerations supported by numerical simulations. The present paper is an extension of our earlier publication [1].


I. INTRODUCTION
Einstein's theory of general relativity, GR, suffers from singularities (see, e.g. [2][3][4][5]). These are pathologies of spacetime, which determine the limits of the validity of GR. They also signal the occurrence of interesting phenomena like black holes and the cosmological big bang, which seem to be supported by observational data.
The present paper is an extension of our recently published paper [1] on the evolution towards the singularity of the general Bianchi IX spacetime. The dynamics of this model underlies the dynamics of the Belinski-Khalatnikov-Lifshitz (BKL) scenario (for an overview see e.g. [6]), which is conjectured to describe the approach towards a generic spacelike singularity of GR. As far as we know, the time evolution of particular curvature invariants of the general BIX has not been considered yet.
The commonly known singularity theorems (see e.g. [2][3][4]) predict the existence of incomplete geodesics, but say little about other features of singularities. Studying particular curvature invariants gives more characteristics of these pathologies. The blowing up of the Kretschmann scalar was proved rigorously by Ringstrom [7] in the vacuum case of the BIX model (mixmaster universe). Barrow and Hervik [8] have studied the Weyl tensor and provided asymptotic expression for homogeneous cosmological models filled with non-tilted perfect fluids. The behaviour of the Weyl squared scalar is of particular interest since it is conjectured that this invariant acts as a measure of gravitational entropy [8,9].
In our paper we examine the time evolution of the Kretschmann scalar for the non-diagonal (general) BIX spacetime which requires the coupling, for instance, to a tilted matter field. The latter is chosen to be dust (for simplicity). The difference in the dynamics between the diagonal and non-diagonal cases has been considered in [10].
Our paper is organized as follows: We first review in section II the Hamiltonian formulation of the dust filled Bianchi IX model which was employed for the numerical simulation of the dynamics in [1]. In section III we compute an expression for the Kretschmann scalar which allows for a numerical evaluation via the framework based on the Hamiltonian formulation.
Section IV is devoted to the asymptotic regime close to the singularity. We provide heuristic considerations predicting the behaviour of the Kretschmann scalar and support our result with a numerical analysis. We conclude in Sec. V.

II. HAMILTONIAN FORMULATION OF THE DUST FILLED BIANCHI IX MODEL
An examination of spatially homogeneous models with tilted fluids is given in [11]. The Bianchi IX model with moving matter had been studied in [12,13]. The Hamiltonian formulation of the Bianchi IX model filled with tilted dust was first derived in a series of papers by Ryan [14,15]. For a similar formulation concerning the other Bianchi models see, for example, [16][17][18]. The metric in the ADM form is given by The basis one-forms 1 satisfy the relation with C i jk = ε ijk being the structure constants of the Lie algebra so(3, R). We parametrize the metric coefficients in this frame as follows: whereh The variables α, β + , and β − are known as Misner variables. 2 The scale factor exp(α) is related to the volume, while the anisotropy factors β + and β − describe the shape of the universe. The variables Γ 1 , Γ 2 and Γ 3 were used by BKL in their original analysis [19]. In addition, we have introduced the SO(3, 1 Throughout this work we choose the units such that 3 4πG σ 1 ∧ σ 2 ∧ σ 3 = 1, which correspond the setting κ = 8πG/c 4 = 6. 2 Note that Misner originally used the variable Ω = −α.
The Euler angles θ, φ, and ψ are dynamical quantities and describe nutation, precession, and pure rotation of the principal axes, respectively. The Hamiltonian and momentum constraints are given by The l i denote (non-canonical) angular momentum-like variables: where we defined the 'angular velocities' ω = {ω i j } ≡ O TȮ and the 'moments of inertia' The l i obey the Poisson bracket algebra {l i , l j } = −C k ij l k . The variable p T , which is in fact a constant of motion, denotes the momentum canonically conjugate to the 'dust time' T , that is, the proper time measured by observes co-moving with the dust particles. The three-curvature scalar on spatial hypersurfaces of constant coordinate time t is given by The formalism described so far is not entirely canonical and must be complemented by the geodesic equation for the dust particles where is a constant of motion and "×" denotes the usual cross product in 3 dimensional Euclidean space. The angular velocities ω can be eliminated from the equations of motion by using the momentum constraints. The equations of motion have a lengthy form and were given in [1]. Moreover, the Hamiltonian formulation can be employed to obtain a qualitative picture of the dynamics (see e.g. [1,[14][15][16]19]).
Our numerical simulations were performed by using the gauge N = e 3α and N i = 0. Setting N = e 3α "moves" the singularities to t = ±∞. This is required to resolve the oscillations when evolving the system towards the singularity. Using the shooting method described in [1] we numerically obtain the solution in terms of the variables as functions of t over some finite time interval [t 1 , t 2 ]. We restrict our attention to the so called tumbling case, that is, all v i are non-zero initially. This case might be considered as the most generic one in the context of the model under consideration. For convenience, we restrict ourselves in this work to the numerical solution which was already considered in [1]. The part of the solution plotted in Fig. 2b extends over one Kasner era. We regard this solution to occur at the transition into the asymptotic regime. The numerical accuracy of our solution has been confirmed by checking the (approximate) preservation of the Hamiltonian constraint H = 0, and the constant of motion C 2 − 1 = 0. Both H and C 2 − 1 stay at order 10 −15 . We believe that the effect of chaos is negligible in a single Kasner era. For further details see [1].

III. CALCULATION OF THE KRETSCHMANN SCALAR
We are interested in the temporal evolution of curvature invariants when approaching the singularity. Particular interest lies in the evolution of the Kretschmann scalar, which can be decomposed as where C µ νλσ is the Weyl tensor, R µν is the Ricci tensor and R the Ricci scalar. For our purposes it is convenient to make use of the constraints and the equations of motion to simplify the expressions such that they are suited for a numerical evaluation. We will do so throughout the calculation in this section and bring our expression into a form that is ready for a numerical evaluation. This means that all expressions should only involve the variables (11) as well as the constants p T ≡ 12p T and C. Furthermore we shall use the quasi-Gaussian gauge N i = 0 while keeping the lapse N unspecified. We now proceed by calculating the terms on the right-hand side of equation (12).
From the Einstein field equations R µν − 1 2 Rg µν = κT µν = κρu µ u ν it follows that we can write Recall that in the model under consideration .
We conclude that the Ricci part 2R µν R µν − 1 3 R 2 of the Kretschmann scalar blows up as when approaching the singularity. The calculation of the Weyl part of the Kretschmann scalar, however, is less trivial. The 3+1 split allows for a decomposition of the Weyl tensor into electric and magnetic part (see e.g. [20]) according to with ikl being the Levi-Civita tensor. The other objects involved in the decomposition are explained below. Now let P µ ν = δ µ ν + n µ n ν denote the projector onto spatial hypersurfaces orthogonal to the normal vector where , S ij and j i are the energy density, the shear density and the momentum density as measured by Eulerian observers (observers with four velocity n µ ).
A direct calculation yields Let us now turn to the computation of E ij E ij , which can be written out as We now evaluate the single terms. The term in the first line of (19) right after the equal sign can be written as We denote the three-dimensional Ricci tensor of the diagonal model by (3)R ij . The three-dimensional Ricci tensor of the non-diagonal model can then obtained via rotation according to (3) kl . The only non-vanishing components of (3)R ij are given by 8 The first term in the second line of (19) reads The three-dimensional Ricci squared scalar can be written as The term in the third line of (19) becomes We can use the Hamiltonian constraint to simplify We therefore obtain a simple expression for the term in the fourth line of (19). Since we have direct numerical access to the quantities in the fourth and fifth line of (19), we will not manipulate them further. It is well known that the Weyl squared scalar vanishes for the Friedmann models. The dust filled closed Friedmann universe is included in the model under consideration as the particular case for which Γ 1 = Γ 2 = Γ 3 and C = 0. As a consistency check of our calculation we convinced ourselves that the Weyl squared scalar vanishes for these restrictions. We find that B ij B ij and E ij E ij vanish separately and hence C µνλσ C µνλσ = 0 as expected.
A numerical evaluation "close" to the point of recollapse 3 is shown in Fig. 1. We note that the bare Kretschmann scalar appears to roughly blow up exponentially in t and it rapidly exceeds the range of numbers that are accessible in Matlab. This is why from now on we turn to numerically evaluating the so-called Hubble normalized Kretschmann scalar R µνλσ R µνλσ /|K i i |. This quantity has the virtue of being dimensionless and numerically well behaved. The expansion scalar is given by   1: Plot (a) shows the result of the application of the shooting method described in [1] "close" to the recollapse (which is roughly at t ≈ −1000). The plot (b) was obtained from a numerical calculation based on the result of section III.

IV. THE ASYMPTOTIC REGIME CLOSE TO THE SINGULARITY
We remark again that we are considering the tumbling case. We expect, however, that a similar discussion also holds for the non-tumbling and non-rotating cases.
In order to simplify the dynamics of Bianchi IX (see, e.g. the review in [19]) we make two assumptions. The first assumption states that the anisotropy of space grows without bound. This means that the solution enters the regime where the ordering of indices depends on the initial conditions and is mostly irrelevant for our discussion. The second assumption 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. For our variables this means that the dust velocities v i assume constant values. The main purpose of the article [1] was to provide a numerical verification for the two assumptions. According to the phrase "matter does not matter" we expect the matter terms in the Kretschmann scalar to be negligible in the asymptotic regime, that is, the Weyl part should dominate over the Ricci part.
During Kasner epochs (i.e. between two successive bounces in the asymptotic regime) we expect the most relevant term to be the term in the first two lines of (20) right after the equality-sign. We therefore assume now that the Kretschmann scalar can be approximated by This claim is confirmed by our numerical simulations. We remark that this term corresponds exactly to the Weyl squared scalar of the Bianchi I model with the metric The Weyl tensor of the Bianchi I model has only an electric part, and the magnetic part vanishes (in the quasi-Gaussian gauge). During Kasner epochs the time evolution can be parameterized using the Lifshitz -Khalatnikov parameter u following the considerations in [19]. Doing so and using the assumption (29) we obtain that the Hubble-normalized Kretschmann scalar can be approximated by Consequently the Kretschmann scalar blows up like the expansion K i i to the power 4 during Kasner eras. In order to understand the temporal evolution of the Kretschmann scalar over the course of one epoch we have plotted the expression on the right hand side of (31) as a function of u in Fig. 2a. It is important that the function has a maximum in u = 1. 4 The BKL refers to bounces from the curvature potential as transformations of the first kind while they call bounces from centrifugal walls transformations of the second kind. Transformations of the first kind change the Lifshitz-Khalatnikov parameter according to u 1 → u − 1. Transformations of the second kind interchange the values of the velocities according to (log Γ 1 ) · 2 → (log Γ 2 ) · , (log Γ 2 ) · 2 → (log Γ 1 ) · and leave the value of u unchanged, i.e. u 2 → u. It follows that 1 → changes the value of the Hubble normalized Kretschmann scalar (31) while 2 → does not. According to the analysis in [19] a typical Kasner era can be expressed as a sequence of n Kasner epochs which starts with an epoch that has a maximum u-value larger than 1 when evolving towards the singularity. The value of u decreases with each transformation of the first kind and ends with the epoch for which u becomes smaller than 1 for the first time, e.g.
It should be remarked at this point that the u-map was found to be asymptotically exact for particular cases (for a collection of rigorous results concerning the u-map see [5]). A solution of the discrete mixmaster map and a detailed study of its chaotic nature for the vacuum Bianchi IX case can be found in [22]. We are now in the position to provide a picture of the behaviour of the Kretschmann scalar over the course of one Kasner era: According to the formula (31) plotted in Fig. 2a and (32) we expect the Hubble normalized Kretschmann scalar to increase its value with each transformation of the first kind before it hits the value u min < 1 for which it decreases again. This is apart from the behaviour in the vicinity of the bounces precisely what we observe in the numerically evaluated Hubble normalized scalar plotted in Fig. 2. We remark that we regard the part of the solution plotted in Fig. 2b to be not quite in the asymptotic regime but rather at the tansition into the asymptotic regime.
It was helpful, in this paper, to support the computations by using the tensor algebra package xAct [23]. Numerical calculations were performed using MATLAB R2016b.

V. SUMMARY
The main purpose of this paper is to provide a description of the temporal behaviour of the Kretschmann scalar in the asymptotic regime. In this regime the volume density, being proportional to the product of the three directional scale factors, evolves towards zero [19], but this is not a satisfactory indication of the singularity as it depends on the choice of coordinates. The blowing up of curvature invariants, on the other hand, shows that we are dealing with a genuine curvature singularity.
During Kasner epochs, R µνλσ R µνλσ increases like the expansion to the power four. Over the course of a single Kasner era the value of the Hubble normalized Kretschmann scalar increases until it drops down to a finite value when it ends. This process will repeat itself with the beginning of the next Kasner era until the system approaches the singularity.
The present paper is supposed to be an extension of our previous paper [1], which considers, for simplicity, only the tilted dust field as a source. The discussion of other tilted fluids goes beyond the scope of our present programme. The effect of tilted radiation, which has been studied analytically and numerically in the recent paper [24], is particularly interesting.
The asymptotic regimes of the Bianchi IX and BVIII models are quite similar [25]. Both models have been used to derive the BKL scenario [26].
The asymptotic regime approximates well the dynamics near the singularity. This is why it was recently used in the struggle for removing the singularity of the BKL scenario by quantization [27].   FIG. 2: The plot (b) was obtained by the shooting method described in [1]. One can see a typical Kasner era composed of epochs that are approximately connected by transformations of the first and second kind. The output of the shooting method was used to obtain the plots (c) and (d) based on the calculation in section III of this paper.