Ioannis Bakas
email: bakas@ajax.physics.upatras.gr

Francois Bourliot
email: bourliot@cpht.polytechnique.fr

Dieter Lüst
email: dieter.luest@lmu.de

Marios Petropoulos
email: marios@cpht.polytechnique.fr

Mixmaster universe in Hořava-Lifshitz gravity

We consider spatially homogeneous (but generally non-isotropic) cosmologies in the recently proposed Hořava-Lifshitz gravity and compare them to those of general relativity using Hamiltonian methods. In all cases, the problem is described by an effective point particle moving in a potential well with exponentially steep walls. Focusing on the closed-space cosmological model (Bianchi type IX), the mixmaster dynamics is now completely dominated by the quadratic Cotton tensor potential term for very small volume of the universe. Unlike general relativity, where the evolution towards the initial singularity always exhibits chaotic behavior with alternating Kasner epochs, the anisotropic universe in Hořava-Lifshitz gravity (with parameter λ > 1/3) is described by a particle moving in a frozen potential well with fixed (but arbitrary) energy E. Alternating Kasner epochs still provide a good description of the early universe for very large E, but the evolution appears to be non-ergodic. For very small E there are harmonic oscillations around the fully isotropic model. The question of chaos remains open for intermediate energy levels.

Introduction

General relativity predicts the existence of space-time singularities under some general conditions, which in a cosmological context are space-like and correspond to the initial singularity of the universe, [START_REF] Hawking | The singularities of gravitational collapse and cosmology[END_REF], [START_REF] Hawking | The large Scale Structure of Space-Time[END_REF]. Studying the asymptotic behavior of Einstein equations in the vicinity of space-like singularities provides the way that the initial state of the universe is reached classically. A broad framework for this purpose was developed in the seminal work of Belinskii, Khalatnikov and Lifshitz, [START_REF] Belinskii | Oscillatory approach to a singular point in the relativistic cosmology[END_REF], [START_REF] Belinskii | Construction of a general cosmological solution of the Einstein equation with a time singularity[END_REF], [START_REF] Belinskii | A general solution of the Einstein equations with a time singularity[END_REF], who found that the spatial points decouple leading to dimensional reduction of the field equations to one (time) coordinate. Then, the universe is described as a point particle moving in an effective potential well, and, remarkably, the dynamical equations are the same as in spatially homogeneous (but generally non-isotropic) cosmological models, in particular Bianchi IX for having a closed universe with spherical topology. Thus, in this context, it becomes important to investigate the main features of the homogeneous and anisotropic cosmologies in the small volume limit of the universe, where the matter sources are ignored.

The dynamics of Bianchi IX model in vacuum (also known as mixmaster universe) has been thoroughly analyzed in the literature by Belinskii, Khalatnikov and Lifshitz, [START_REF] Belinskii | Oscillatory approach to a singular point in the relativistic cosmology[END_REF], [START_REF] Belinskii | Construction of a general cosmological solution of the Einstein equation with a time singularity[END_REF], [START_REF] Belinskii | A general solution of the Einstein equations with a time singularity[END_REF], and independently by Misner who used Hamiltonian methods, [START_REF] Misner | Mixmaster universe[END_REF], [START_REF] Misner | Quantum cosmology. I[END_REF], [START_REF] Misner | Classical and quantum dynamics of a closed universe[END_REF], [START_REF] Misner | Minisuperspace[END_REF], [START_REF] Misner | The mixmaster cosmological metrics[END_REF] (but see also the textbook [START_REF] Misner | Gravitation[END_REF] and the monograph [START_REF] Ryan | Hamiltonian Cosmology[END_REF]). The results can be summarized in a nut-shell by saying that the evolution consists of alternating Kasner epochs, acting as oscillations that permute the principal axes of the spherical spatial slices, as the universe is approaching the initial singularity. A good picture of the dynamics close to the singularity is then provided by a billiard motion in a finite region of Lobachevsky plane, which turns out to be chaotic. More general studies of the chaotic behavior of the mixmaster universe have been carried out in the literature over the years, [START_REF] Barrow | Chaotic behaviour in general relativity[END_REF], [START_REF] Chernoff | Chaos in the mixmaster universe[END_REF], [START_REF] Rugh | Chaos in the Einstein equations -characterization and importance[END_REF], [START_REF] Cornish | Mixmaster universe: A chaotic Farey tale[END_REF], and the whole subject is now well established and understood for general relativity in four space-time dimensions. Several generalizations have also been considered in detail, in particular in the context of higher dimensional Kaluza-Klein theories of gravity, [START_REF] Elskens | Ergodic theory of the mixmaster universe in higher spacetime dimensions[END_REF], [START_REF] Demaret | Are Kaluza-Klein models of the universe chaotic?[END_REF], where the billiard picture appears to be universal and the criteria for the appearance of chaos can be formulated in Lie algebraic terms that depend on the dimensionality of space-time, [START_REF] Damour | Oscillatory behaviour in homogeneous string cosmology models[END_REF], [START_REF] Damour | Hyperbolic Kac-Moody algebras and chaos in Kaluza-Klein models[END_REF], [START_REF] Damour | Einstein billiards and overextensions of finite-dimensional simple Lie algebras[END_REF] (but see also [START_REF] Damour | Cosmological billiards[END_REF] for a comprehensive review of these matters and references therein). In parallel, there have also been studies of the same problem in higher curvature generalizations of gravity, in particular in four dimensions, by adding R 2 (and possibly other) curvature terms to the gravitational action, where the chaotic behavior was found to be absent, [START_REF] Barrow | Mixmaster cosmological model in theories of gravity with a quadratic Lagrangian[END_REF], [START_REF] Barrow | Chaotic behaviour in higher-order gravity theories[END_REF], [START_REF] Cotsakis | Mixmaster universe in fourthorder gravity theories[END_REF]. Thus, the subject is quite rich and interesting in all generally covariant effective gravitational theories, including those that arise from string theory.

Recently, there has been a rather odd proposal in the literature to replace the rel-ativistic theory of Einstein gravity by a non-relativistic field theory of Lifshitz type that is only applicable to the ultra-violet regime, [START_REF] Hořava | Membranes at quantum criticality[END_REF], [START_REF] Hořava | Quantum gravity at a Lifshitz point[END_REF]. The resulting theory became known as Hořava-Lifshitz gravity and it is by construction a higher derivative modification of ordinary general relativity with anisotropic scaling in the space and time coordinates. As such, its field equations contain second derivatives in time and higher derivatives in space coordinates (actually up to six in four space-time dimensions where the present work will focus). This proposal aims to provide a renormalizable theory of quantum gravity at short distances that flows to ordinary general relativity in the infra-red domain of large distances. It is, however, quite different in nature from the (more conventional) higher derivative generalizations of Einstein gravity that have been considered so far, which remain fully covariant, whereas here the modification by higher curvature terms affects only the spatial dimensions.

It should be said straight from the beginning that there are no general theorems for the existence of singularities in Hořava-Lifshitz gravity, and under which conditions these may be valid, in particular for the existence of an initial space-like singularity in a cosmological context 1 . Furthermore, no analysis has been made so far on how these singularities, if present, will be approached asymptotically at very early times (one may phrase it by simply saying that the analysis "BKL for HL" is still lacking).

Remarkably, it can be seen that all spatially homogeneous (but generally anisotropic) cosmological space-times, including, in particular, the Bianchi IX model, provide consistent mini-superspace truncations of the field equations in Hořava-Lifshitz gravity, as in general relativity. In this paper we begin studying these models in detail, first because they are interesting in their own right, as they can provide a basis for comparing the two different theories of gravity in the classical and (hopefully in the future) in their quantum regime, and, second because they can also play a role in understanding the evolution of the universe close to the initial singularity, as in general relativity. In all cases, if the spatial points decouple close to the singularity, which is a reasonable expectation in general, the closed universe will be effectively described by mixmaster dynamics, viewed as point particle moving in a potential well whose structure depends on the theory. The framework that will be adopted throughout this study is that of Hamiltonian dynamics, since the action of Hořava-Lifshitz gravity is only defined through the 3 + 1 ADM decomposition of the space-time metric.

The rest of this paper is organized as follows: Section 2 provides a brief overview 1 The only case that has already been extensively studied in the literature is the analogue of Friedmann universe, which is an isotropic Bianchi IX model with suitable matter sources, and was found to exhibit a bounce -rather than a singularity -in the past under some technical but rather general conditions on the couplings of the theory, [START_REF] Calcagni | Cosmology of the Lifshitz universe[END_REF], [START_REF] Kiritsis | Hořava-Lifshitz cosmology[END_REF], [START_REF] Brandenberger | Matter bounce in Hořava-Lifshitz cosmology[END_REF]. Such a bounce, however, might only be attributed to the particular model, since the addition of shearing components, due to anisotropies, could circumvent it and render it unstable. We will discuss more extensively this point later, in section 4, offering also some simple anisotropic models that do not exhibit a bounce.

of Hořava-Lifshitz gravity in 3 + 1 dimensions and introduces the necessary notions in self-contained way. Section 3 explains why the homogeneous cosmologies are consistent models in Einstein gravity as well as in Hořava-Lifshitz gravity and describes mixmaster dynamics as an effective particle model moving in a potential well that is applicable to both theories. The effective potentials are derived in each case separately. Section 4 investigates the evolution of the universe close to the initial singularity, where the problem simplifies considerably, but it still exhibits rich structure. Unlike general relativity, where the potential vanishes for generic values of the anisotropy parameters and the evolution towards the initial singularity proceeds in an oscillatory fashion with alternating Kasner epochs, the universe in Hořava-Lifshitz gravity (with parameter λ > 1/3) is described by a particle moving in a frozen potential well with prescribed (but arbitrary) energy. The question of chaos in the corresponding motion is briefly addressed. Finally, section 5 contains our conclusions and discusses some open questions and directions for further research. Two appendices are also included at the end. The first contains useful formulas for the Bianchi IX model of three-geometries, collecting, in particular, the expressions for its Ricci curvature and Cotton tensors. The second contains a derivation of the bounce law from the exponentially steep walls of the potential well that will be used in mixmaster dynamics.

Hořava-Lifshitz gravity

General relativity as well as Hořava-Lifshitz gravity are formulated in a similar fashion using the ADM decomposition of the four-dimensional metric in space-time M 4 , which is assumed to be topologically R × Σ 3 ,

ds 2 = -N 2 dt 2 + γ ij (dx i + N i dt)(dx j + N j dt) .
(2.1)

The three-dimensional slices Σ 3 have metric γ ij and extrinsic curvature tensor

K ij = 1 2N ∂γ ij ∂t -∇ i N j -∇ j N i . (2.2)
The space of all three-dimensional metrics γ ij , which is known as Wheeler-DeWitt superspace, is very important in this study. It is endowed with a metric, often called DeWitt metric, [START_REF] Dewitt | Quantum theory of gravity: 1. The canonical theory[END_REF], which is taken here to depend on a parameter λ in general. The metric in superspace, and its inverse, are defined as usual,

G ijkl = 1 2 γ ik γ jl + γ il γ jk -λγ ij γ kl , (2.3)
and

G ijkl = 1 2 γ ik γ jl + γ il γ jk - λ 3λ -1 γ ij γ kl , (2.4) so that G ijmn G mnkl = 1 2 δ k i δ l j + δ l i δ k j . (2.5)
The DeWitt metric is positive definite for λ < 1/3 and indefinite for λ > 1/3, which includes, in particular, the special value λ = 1 applicable to Einstein gravity.

The gravitational theories that will be considered in the sequel, using the ADM formalism (see, for instance, the textbook [START_REF] Misner | Gravitation[END_REF]), admit a four-dimensional action

S = S K -S V , (2.6)
where S K is the kinetic part of the action with universal form

S K = 2 κ 2 dtd 3 x detγ N K ij G ijkl K kl . (2.7)
The potential part of the action, S V , is given by

S V = dtd 3 x detγ N V , (2.8)
where V is chosen according to the theory.

Ordinary general relativity corresponds to λ = 1 and has potential term

V GR = - 2 κ 2 (R -2Λ) (2.9)
that involves second derivatives in the space coordinates. Here, Λ is the cosmological constant in M 4 , R is the Ricci scalar curvature of the three-dimensional metric γ ij and κ 2 = 32πG is expressed in terms of Newton's constant in four space-time dimensions. On the other hand, Hořava-Lifshitz gravity has a potential that involves higher order terms, thus breaking relativistic invariance of the four-dimensional theory, [START_REF] Hořava | Membranes at quantum criticality[END_REF], [START_REF] Hořava | Quantum gravity at a Lifshitz point[END_REF]. These terms have a specific form, composed of several higher order (quadratic) curvature corrections, which are designed to smooth out the ultra-violet behavior of gravity. Also, the parameter λ is left undetermined in this context and may run with the energy scale in quantum theory.

A particularly simple choice of V in Hořava-Lifshitz gravity, though by no means unique, corresponds to the so called "detailed balance" condition, meaning

V HL = κ 2 2 E ij G ijkl E kl , (2.10)
so that V is derived from a superpotential W in the sense

E ij = 1 2 detγ δW δγ ij .
(2.11)

In 3 + 1 dimensions that will be considered here, the superpotential is taken to be the Euclidean action of three-dimensional topological gravity on Σ 3 with cosmological constant Λ w (other than Λ),

W = 1 w 2 W CS + µW EH , (2.12)
where the first term refers to the gravitational Chern-Simons action, [START_REF] Deser | Three-dimensional massive gauge theories[END_REF],

W CS = Σ 3 d 3 x detγ ijk Γ l im ∂ j Γ m kl + 2 3 Γ m jn Γ n kl , (2.13)
with 123 = 1, and the second term is the corresponding three-dimensional Einstein-Hilbert action

W EH = Σ 3 d 3 x detγ (R -2Λ w ) . (2.14)
Thus, V HL with "detailed balance" follows by computing E ij varying the superpotential functional W with respect to the metric γ ij . The result reads

E ij = 1 w 2 C ij - µ 2 G ij , (2.15)
where C ij is the Cotton tensor of the metric γ ij defined as follows,

C ij = 1 2 detγ δW CS δγ ij = 1 detγ ikl ∇ k R j l - 1 4 δ j l R . (2.16)
This is a symmetric and traceless tensor that vanishes if and only if the three-dimensional metric is conformally flat. The second term is the familiar Einstein tensor on Σ 3 with three-dimensional cosmological constant Λ w ,

G ij = - 1 detγ δW EH δγ ij = R ij - 1 2 Rγ ij + Λ w γ ij . (2.17)
Putting all together, the potential terms of Hořava-Lifshitz gravity in 3 + 1 dimensions satisfying the "detailed balance" condition are

V HL = αC ij C ij + βC ij R ij + γR ij R ij + δR 2 + R + ζ (2.18) with coefficients α = κ 2 2w 4 , β = - µκ 2 2w 2 , γ = µ 2 κ 2 8 , δ = - µ 2 κ 2 (4λ -1) 32(3λ -1) , = µ 2 κ 2 Λ w 8(3λ -1) , ζ = - 3µ 2 κ 2 Λ 2 w 8(3λ -1) . (2.19)
The last two terms are identical to the potential V GR of general relativity, with the appropriate identifications of the coefficients and ζ, whereas the remaining ones are higher order (quadratic) curvature corrections that apparently are suppressed in the infra-red limit of the theory.

More general choices of V HL , other than "detailed balance", are also admissible and correspond to the sum (2.18) with arbitrary coefficients; they are only subject to the restriction that general relativity will emerge in the infra-red regime. The analysis that will be performed in the sequel applies equally well to all such general choices of potential in Hořava-Lifshitz gravity with or without "detailed balance".

Finally, it is important to note that the Hořava-Lifshitz theory of gravity is not invariant under general coordinate transformations in space-time; this should be contrasted with other higher order theories of gravity that remain relativistic. Since M 4 is topologically R × Σ 3 , it is only appropriate here to consider invariance of the action under the restricted class of foliation preserving diffeomorphisms, t = t(t) , xi = xi (t, x) .

(2.20)

Thus, the lapse function N associated to the freedom of time reparametrization is restricted to be a function of t alone, whereas the shift vector N i associated with diffeomorphisms of Σ 3 can depend on all space-time coordinates.

Mixmaster dynamics

The mixmaster universe arises as mini-superspace model of gravity assuming that the the three-dimensional slices Σ 3 are homogeneous geometries with the topology of S 3 and isometry group SU [START_REF] Hawking | The large Scale Structure of Space-Time[END_REF]. Thus, by employing the Bianchi IX ansatz, as explained in Appendix A, the four-dimensional metric is given in diagonal form2 by

ds 2 = -N 2 (t)dt 2 + γ 1 (t)σ 2 1 + γ 2 (t)σ 2 2 + γ 3 (t)σ 2 3 , (3.1)
using the invariant 1-forms σ i of SU [START_REF] Hawking | The large Scale Structure of Space-Time[END_REF]. The metric coefficients are all taken to depend only on the time coordinate t. This class of metrics provides consistent reduction of vacuum Einstein equations to an autonomous system of ordinary non-linear differential equations for γ i (t) that has been studied extensively in the literature in the past fourty years. It provides a simple model of homogeneous, but generally anisotropic, universe, which proves valuable for studying the chaotic behavior of general relativity close to the initial singularity. The same ansatz also works consistently for the Hořava-Lifshitz gravity with or without "detailed balance" and gives rise to another -though more complicated -system of ordinary non-linear differential equations for the coefficients γ i (t).

The purpose of this section is to describe in detail the mini-superspace reduction of the field equations and transform them into an effective point particle model using Hamiltonian methods, as outlined by Misner, [START_REF] Misner | Mixmaster universe[END_REF], [START_REF] Misner | Quantum cosmology. I[END_REF], [START_REF] Misner | Classical and quantum dynamics of a closed universe[END_REF], [START_REF] Misner | Minisuperspace[END_REF], [START_REF] Misner | The mixmaster cosmological metrics[END_REF] (but see also the textbook [START_REF] Misner | Gravitation[END_REF] and the monograph [START_REF] Ryan | Hamiltonian Cosmology[END_REF]). Although our discussion is entirely confined to the Bianchi IX case, it should be noted here that all homogeneous spaces arising in the Bianchi classification of model three-geometries provide consistent reduction of general relativity as well as Hořava-Lifshitz gravity. The details for all other homogeneous cosmologies in Hořava-Lifshitz gravity will not be included here.

Effective particle model

Hamiltonian methods for homogeneous cosmologies are most appropriate to use in the ADM decomposition of space-time and they lead naturally to an effective point particle model with appropriately chosen external potential. The method is applicable to both general relativity and Hořava-Lifshitz gravity because the lapse function N is taken to depend only on t in such cases.

Recall that the canonical momenta conjugate to γ ij are simply given by

π ij = δS δ(∂γ ij /∂t) = 2 κ 2 detγ G ijkl K kl (3.2)
given the dependence of the general gravitational action S = S K -S V upon the extrinsic curvature, whereas the momenta conjugate to N and N i vanish. Then, the Hamiltonian form of the action is

S = dtd 3 x π ij ∂γ ij ∂t -NH -N i H i , (3.3)
where

H = κ 2 2 detγ π ij G ijkl π kl + detγ V (3.4)
and

H i = -2∇ j π ij . (3.5)
In ordinary general relativity H = 0 and H i = 0 are the constraints of the theory associated to general reparametrization invariance in the 3 + 1 decomposition of the metric. In this context, N and N i serve as Lagrange multipliers whose variation yields the constraints. On the other hand, the invariance of Hořava-Lifshitz gravity under the restricted class of foliation preserving diffeomorphisms leaves intact the momentum constraints H i = 0 and replaces the time constraint H = 0 with the much weaker condition

d 3 xH = 0 . (3.6)
For homogeneous cosmologies, H is only a function of t, and so is the lapse function N, and, therefore, there is no difference in the Hamiltonian description of the two theories other than the form of the potential. Also, in these cases, the momentum constraints are satisfied identically, since γ ij and π ij are functions of t only and the covariant derivative in H i with respect to space coordinates reduces to an ordinary derivative. Thus, one may consistently choose N i = 0 and forget altogether the momenta constraints.

Based on these observations, the Hamiltonian form of the gravitational action for homogeneous cosmologies takes the following form

S = 16π 2 dt π ij dg ij dt -NH (3.7)
after performing the integration over space that accounts for the 16π 2 factor. Here,

H = κ 2 2 detγ π ij π ij - λ 3λ -1 (π k k ) 2 + detγ V (3.8)
depends only on t and it is constrained to vanish by varying S with respect to N.

At this point, one may employ the freedom of time reparametrizations t = t(t) to eliminate N from the variational problem. Any choice of N(t) inserted in equation (3.7) leaves the action in canonical Hamiltonian form, but the content of the gauge fixed action will be equivalent to the original one only if it is supplemented by the constraint H = 0, which can no longer be derived from the variational principle. The most convenient choice is

N(t) = 6 κ 2 detγ(t) , (3.9)
which will be adopted from now on. Thus, we arrive at an effective point particle Hamiltonian model for all homogeneous cosmologies; yet another formulation of this point particle model will be mentioned shortly.

Further simplification occurs by introducing the volume and shape moduli of the three-geometry

γ 1 = e 2Ω+β + + √ 3β -, γ 2 = e 2Ω+β + - √ 3β -, γ 3 = e 2Ω-2β + , ( 3.10)
as explained in Appendix A, so that detγ = √ γ 1 γ 2 γ 3 = exp(3Ω). Then, it is appropriate to parametrize the components of the momenta matrix π i j , which is diagonal as is the metric matrix, as follows,

p Ω = 2π k k , p i j = π i j - 1 3 δ i j π k k , (3.11)
and set

p 1 1 = 1 6 (p + + √ 3 p -) , p 2 2 = 1 6 (p + - √ 3 p -) , p 3 3 = - 1 3 p + . (3.12)
It turns out that p Ω is the conjugate momentum to the volume moduli Ω and p ± are the conjugate momenta to the shape moduli β ± . In terms of these variables, the gauge fixed action S becomes

S = 16π 2 dt p Ω dΩ dt + p + dβ + dt + p - dβ - dt -H , (3.13)
where t denotes now the specific time coordinate following the choice (3.9) and

H = 1 2 p 2 + + p 2 -- 1 2(3λ -1) p 2 Ω + V(β + , β -, Ω) = 0 . (3.14)
With the given choice of lapse, the effective potential of the model is

V = 6 κ 2 e 6Ω V . (3.15)
This is the final form of the action that will be used to study the mixmaster universe in general relativity and Hořava-Lifshitz gravity and compare the results for the two theories. Note that in writing H above we left undetermined the parameter λ of the superspace metric in order to apply it to all cases of present interest.

Finally, as promised, we mention for completeness that there is an alternative formulation of the same effective point particle model based on the reduced ADM action. In this approach, which has been mostly used by Misner (and many others) in general relativity, one eliminates Ω by solving the Hamiltonian constraint for p Ω and sets t = Ω, choosing also N appropriately. The procedure is similar to going from the quadratic form of the action of a relativistic particle to its square-root form. Then, the reduced variational problem corresponds to a point particle moving in two dimensions (provided by β + and β -alone) under the influence of a time dependent potential well. The reduced ADM formulation is only adequate for λ > 1/3, but it will not be used here in any case.

Next, we will specialize the discussion to general relativity and Hořava-Lifshitz gravity and derive, in each case separately, the effective potential of the mixmaster universe that enters in the canonical Hamiltonian variational principle (3.13), (3.14). Hamilton's equations follow by varying Ω, β ± , p Ω and p ± and yield

d 2 β ± dt 2 = - ∂V ∂β ± , d 2 Ω dt 2 = 1 2(3λ -1) ∂V ∂Ω . (3.16)
We will not write down explicitly the resulting equations of motion -they cannot be solved in any case -but focus only on the qualitative features of mixmaster dynamics, following from the potential, which make it a valuable tool for exploring the behavior of the universe as it approaches the initial singularity3 .

General relativity

First, we consider the mixmaster dynamics in general relativity and present the form of the effective particle potential for comparison later with Hořava-Lifshitz gravity. It follows by expressing the three-dimensional Ricci scalar R, given in Appendix A, in terms of Ω and β ± and reads, [START_REF] Misner | Mixmaster universe[END_REF], [START_REF] Misner | Quantum cosmology. I[END_REF], [START_REF] Misner | Classical and quantum dynamics of a closed universe[END_REF], [START_REF] Misner | Minisuperspace[END_REF], [START_REF] Misner | The mixmaster cosmological metrics[END_REF],

V GR = 6 κ 4 e 4Ω 2e 2β + cosh(2 √ 3β -) -1 -4e -β + cosh( √ 3β -) + e -4β + + 24Λ κ 4 e 6Ω .
(3.17)

The equations of motion that follow from the variational principle (setting also λ = 1) cannot be solved exactly. However, they have been studied extensively in the literature for many years and found to exhibit some very interesting qualitative features close to the initial singularity, where the volume of the universe vanishes as Ω → -∞. These features will be discussed in some detail later. It will also be helpful in this context to have a good qualitative picture of the potential well.

The effective potential of mixmaster dynamics corresponds to a well shown in Fig. 1 for fixed Ω. It has three canyon lines located at β -= 0 and β -= ± √ 3β + , where any two γ i 's become equal. The potential is bounded from below and exhibits discrete Z 3 symmetry by permuting the principal axes of rotation of S 3 . Thus, it has the shape of an equilateral triangle in the anisotropy space (β + , β -) and exponentially steep walls far away from the origin. Very close to the origin the well is round, as can be seen by expanding V GR up to quadratic order in β ± . As can be seen, they extend symmetrically between the canyon lines β -= 0 and β -= ± √ 3β + , which correspond to a partially anisotropic universe with axial symmetry (known as Taub space-time, [START_REF] Taub | Empty spacetimes admitting a three parameter group of motions[END_REF]). The fully isotropic case corresponds to the origin β + = 0 = β -, where the potential attains its minimum.

The asymptotic form of the potential for very large values of anisotropy is independent of the cosmological constant Λ and looks like

V GR 6 κ 4 e 4Ω e -4β + , as β + → -∞ (3.18) and V GR 72 κ 4 e 4Ω β 2 -e 2β + , as β + → +∞ , |β -| << 1 . (3.19)
Then, the asymptotic structure of the potential is completely characterized by combining these two relations with the triangular symmetry of the model.

The effective point particle can only escape to infinity along the canyon lines where the potential has the shape shown in Fig. 3, keeping Ω fixed. The smallest deviation from axial symmetry will turn the particle against the infinitely steep walls.

Hořava-Lifshitz gravity

Applying the same framework to Hořava-Lifshitz gravity, the effective point particle potential is described by

V HL = 6 κ 2 e 6Ω αC ij C ij + βC ij R ij + γR ij R ij + δR 2 + R + ζ . (3.20)
The coefficients are left arbitrary so that the discussion can be made as general possible without necessarily imposing the "detailed balance" condition. The only restriction we put here is α > 0 for well-definiteness, so that V HL (β + , β -) stays bounded from below.

Using the explicit expressions for C ij , R ij and R found in Appendix A and rewriting them in terms of the variables Ω and β ± , we obtain the following results for the individual terms entering in the effective potential of Bianchi IX cosmology:

C ij C ij = 1 2 e -6Ω 2e 6β + 3cosh(6 √ 3β -) -3cosh(4 √ 3β -) + cosh(2 √ 3β -) -1 - 2e 3β + 3cosh(5 √ 3β -) -cosh(3 √ 3β -) -2cosh( √ 3β -) + 2 cosh(4 √ 3β -) + 2cosh(2 √ 3β -) -3 - 4e -3β + cosh(3 √ 3β -) -cosh( √ 3β -) + e -6β + 2cosh(2 √ 3β -) + 1 - 6e -9β + cosh( √ 3β -) + 3e -12β + , (3.21
)

C ij R ij = -e -5Ω 2e 5β + cosh(5 √ 3β -) -cosh(3 √ 3β -) -2e 2β + cosh(4 √ 3β -) -cosh(2 √ 3β -) + e -4β + -2e -7β + cosh( √ 3β -) + e -10β + , (3.22)
R ij R ij = 1 4 e -4Ω 2e 4β + 3cosh(4 √ 3β -) -4cosh(2 √ 3β -) + 1 - 8e β + cosh(3 √ 3β -) -cosh( √ 3β -) + 4e -2β + cosh(2 √ 3β -) + 1 - 8e -5β + cosh( √ 3β -) + 3e -8β + , (3.23) R 2 = 1 4 e -4Ω 2e 4β + cosh(4 √ 3β -) -4cosh(2 √ 3β -) + 3 - 8e β + cosh(3 √ 3β -) -cosh( √ 3β -) + 4e -2β + 3cosh(2 √ 3β -) + 1 - 8e -5β + cosh( √ 3β -) + e -8β + (3.24) and R = - 1 2 e -2Ω 2e 2β + cosh(2 √ 3β -) -1 -4e -β + cosh( √ 3β -) + e -4β + . (3.25)
The potential V HL (β + , β -) also has the shape of an equilateral triangle with exponentially steep walls when α > 0. Fig. 1 and Fig. 2 still provide a good qualitative picture of it far from the origin in (β + , β -) parameter space. Only the bottom area close to the origin has slightly different shape that depends on the relative coefficients of the individual terms of the potential. The equations of motion that follow from it provide the analogue of mixmaster dynamics in Hořava-Lifshitz gravity. We will not attempt to solve them here but rather confine ourselves to study some qualitative features that make the model useful for early time cosmology, as in general relativity.

For α > 0, which will be assumed from now on, the asymptotic form of the potential V HL for very large values of anisotropy is dominated completely by the quadratic Cotton tensor term, which happens to contain the steepest walls of all terms, and one has

V HL 9 κ 2 αe -12β + , as β + → -∞ . (3.26)
Likewise, we have

V HL 576 κ 2 αβ 2 -e 6β + , as β + → +∞ , |β -| << 1 . (3.27)
Thus, unlike general relativity, we note that the asymptotic form of the potential is independent of Ω.

As before, the effective point particle can only escape to infinity along the canyon lines β -= 0 and β -= ± √ 3β + which arise for general choices of V HL .

Β V Figure 4: The form of the potential V Cotton (β + , β -) along the canyon line β -= 0 Fig. 4 shows only the plot of the quadratic Cotton term of the potential, V Cotton , along one of these lines, say β -= 0; it exhibits a local maximum at β + = (log2)/3, where γ 1 = γ 2 = 2γ 3 , and similar local maxima show up along the other two lines obtained by permuting γ i . As will be seen later V Cotton is the only relevant term of V HL in early time cosmology.

Approach to the initial singularity

In this section we examine the dynamical behavior of the universe close to the initial singularity, where Ω → -∞, using the Bianchi IX model in vacuum. In this context, it is important to have anisotropic models with general parameters β ± , since, otherwise, the universe will not be able close up to S 3 without radiation or matter sources. First, we will make some general -though crude -remarks about the existence of the initial singularity and then study the problem in question for the two different theories of gravitation.

General considerations

The isotropic Bianchi IX case in vacuum corresponds to a closed Robertson-Walker space-time

ds 2 = -N 2 (t)dt 2 + e 2Ω(t) dΩ 2 3 , (4.1)
with N(t) = (6/κ 2 )exp(3Ω). Note, however, that this metric cannot remain isotropic more than instantaneously; consistency of the dynamics also requires the presence of some shearing components provided by space anisotropy in vacuum -beyond the pure dilation -or suitable sources and combinations thereof.

In general relativity, this follows from the form of the potential assuming β + = 0 = β -for all time. Since the potential

V isotropic GR (Ω) = - 6 κ 4 3e 4Ω -4Λe 6Ω (4.2)
is always negative for small volume, irrespective of Λ, it fails to satisfy the Hamiltonian constraint (setting λ = 1),

1 8 p 2 Ω = 2 dΩ dt 2 = V isotropic GR (Ω) . (4.3)
Adding sources, in the form of perfect fluid, remedies the situation and yields the Friedmann universe. Recall in this case that the potential density V GR is modified by adding the (positive) contribution of the energy density ρ of the fluid, so that 4

2

dΩ dt

2 = V GR + µ e ν Ω , (4.4)
where the contribution of non-relativistic matter sources corresponds to ν = 3 (since ρ ∼ V -1 ) and that of radiation to ν = 2 (since ρ ∼ V -4/3 ) with V = exp(3Ω). Clearly, the radiation term provides the dominant contribution in the small volume limit and one obtains an expanding isotropic and homogeneous universe with initial singularity at T = 0. If deviations from isotropy are subdominant in the small volume limit, compared to radiation, the initial singularity will persist in the presence of sources. This argument only applies to the anisotropy terms of the potential V GR that depends on the volume as exp(4Ω) via the curvature. The anisotropy, however, introduces additional (positive definite) terms in the Hamiltonian constraint, namely the kinetic energy of the anisotropy parameters (dβ + /dt) 2 + (dβ -/dt) 2 if one considers the Bianchi IX model, which completely dominate the evolution close to the initial singularityessentially trying to "avoid" it by oscillations -as will be discussed later. In all cases, the universe can come to a singular state at T = 0 satisfying the field equations.

In Hořava-Lifshitz gravity, the isotropic potential does not receive contribution from the Cotton tensor, since C ij vanishes identically. Thus, the potential given in general by equation (3.20), without necessarily assuming "detailed balance" condition, becomes V isotropic HL

(Ω) = 3 2κ 2 3(γ + 3δ)e 2Ω + 6 e 4Ω + 4ζe 6Ω . (4.5)

The first term arises from the combined effect of R ij R ij and R 2 and dominates the dynamics for small volume. The Hamiltonian constraint now reads 1 4(3λ -1)

p 2 Ω = (3λ -1) dΩ dt 2 = V isotropic HL (Ω) (4.6)
and cannot be possibly fulfilled when

(3λ -1)(γ + 3δ) < 0 . (4.7)
This inequality, which is certainly satisfied for in the case of "detailed balance" condition (see the choice of coefficients (2.19)), means that the quadratic curvature terms correspond to "dark radiation" (since they effectively have "ν = 2"), but with negative energy density. It also implies that the universe cannot evolve isotropically in vacuum without turning on some shearing components, as in general relativity.

Adding sources, in the form of perfect fluid, leads to an interesting possibility when inequality (4.7) is fulfilled with λ > 1/3. In analogy with the previous analysis one obtains

(3λ -1) dΩ dt 2 = V HL + µ e ν Ω (4.8)
and the dominant contribution in the small volume limit is provided by the quadratic Ricci curvature terms and the matter sources with suitable ν. Then, isotropic evolution becomes possible, leading to a Friedmann universe in Hořava-Lifshitz gravity. When ν > 2, there can be a bounce in Ω that replaces the initial singularity of the universe, [START_REF] Calcagni | Cosmology of the Lifshitz universe[END_REF], [START_REF] Kiritsis | Hořava-Lifshitz cosmology[END_REF], [START_REF] Brandenberger | Matter bounce in Hořava-Lifshitz cosmology[END_REF], as can be easily seen by neglecting the contribution of the curvature R and the cosmological constant term. This is the only case for which the energy density of "dark radiation" can grow with respect to the regular matter energy5 . Note, however, that possible deviations from isotropy will become dominant in the small volume limit, since the quadratic Cotton tensor term, which is independent of Ω, kicks in V HL and washes away the effect of the previously thought relevant terms. The kinetic energy of the anisotropy parameters also contributes on equal footing. This indicates that the cosmological bounce is unstable against anisotropy, and, generically, the universe can come in a singular state at T = 0 satisfying the field equations. The validity or not of inequality (4.7) becomes irrelevant in the presence of anisotropy. Consistency also requires λ > 1/3, otherwise the Hamiltonian constraint cannot be possibly satisfied in the small volume limit in the presence of anisotropy; by the same token, the universe can only remain still in an isotropic state when λ < 1/3, and, therefore, this possibility will not be considered further.

Although the argument above does not provide a rigorous proof for the existence of an initial singularity in Hořava-Lifshitz cosmology, and under which general conditions this may be possible, it seems sufficient for the purposes of the present work 6 . Thus, in the following, we will use mixmaster cosmology to explore the approach to the initial singularity, as in general relativity.

General relativity

The potential V GR appears to vanish as one approaches the initial singularity. This is true for generic values of the anisotropy parameters β ± implying Kasner behavior of the universe close to the singularity, which is taken to occur at the beginning of cosmic time, [START_REF] Belinskii | Oscillatory approach to a singular point in the relativistic cosmology[END_REF], [START_REF] Belinskii | Construction of a general cosmological solution of the Einstein equation with a time singularity[END_REF], [START_REF] Belinskii | A general solution of the Einstein equations with a time singularity[END_REF], [START_REF] Misner | Mixmaster universe[END_REF], [START_REF] Misner | Quantum cosmology. I[END_REF], [START_REF] Misner | Classical and quantum dynamics of a closed universe[END_REF], [START_REF] Misner | Minisuperspace[END_REF], [START_REF] Misner | The mixmaster cosmological metrics[END_REF]. In fact, since the Ricci scalar curvature of the homogeneous space Σ 3 vanishes in this limit, the space is effectively flat, as in Bianchi type-I cosmology, and it is more appropriate to use Cartesian dx, dy and dz instead of the 1-forms σ i of SU [START_REF] Hawking | The large Scale Structure of Space-Time[END_REF].

More precisely, when the potential vanishes, all momenta are constant satisfying p 2 Ω = 4(p 2 + + p 2 -) by the Hamiltonian constraint (with λ = 1). Then, it is convenient to introduce the following parametrization of the constant momenta,

n 1 = 1 3p Ω p Ω -2p + -2 √ 3p -, n 2 = 1 3p Ω p Ω -2p + + 2 √ 3p -, n 3 = 1 3p Ω (p Ω + 4p + ) , (4.9)
so that

n 1 + n 2 + n 3 = 1 = n 2 1 + n 2 2 + n 2 3 . (4.10)
The remaining equations dβ ± /dt = p ± and dΩ/dt = -p Ω /4 (so that dβ ± /dΩ = -4p ± /p Ω ) can be easily solved to yield the metric coefficients γ i (t) = T 2n i with respect to a time frame T defined by absorbing the lapse function as dT = N(t)dt. Then, the metric takes the familiar Kasner form

ds 2 = -dT 2 + T 2n 1 dx 2 + T 2n 2 dy 2 + T 2n 3 dz 2 , ( 4.11)
which describes an expanding flat universe with linearly varying volume element, detγ = T.

Thus, the mixmaster dynamics close to the initial singularity appears to follow the Kasner evolution with some fixed parameters (n 1 , n 2 , n 3 ). The Kasner universe is anisotropic as it always contains a direction, say z, along which distances contract rather than expand; this follows from the algebraic conditions (4.10), which imply that one of the n i 's, say n 3 , is lying in the interval -1/3 ≤ n 3 ≤ 0. We may order the Kasner exponents as

-1/3 ≤ n 3 ≤ 0 ≤ n 2 ≤ 2/3 ≤ n 1 ≤ 1 (4.12)
without loss of generality. The axially symmetric case corresponds to the choice of parameters n 1 = n 2 = 2/3 and n 3 = -1/3 (and permutations of the axes thereof).

This description is valid at generic points of (β + , β -) parameter space, but it can break down far away from the origin when the effective point particle experiences the exponentially steep walls of the potential V GR . For example, when the particle approaches one of the triangular walls arising at β + → -∞, as it moves within the wedge |β -| < -√ 3β + , the dominant term of the potential is proportional to exp[4(Ωβ + )], as can be seen from the asymptotic behavior (3.18), and will become sufficiently large to influence the motion if dβ + /dΩ > 1 7 . In the simplest case of axially symmetric evolution towards the wall, so that β -stays zero along the trajectory and the Kasner exponents are n 1 = n 2 = 2/3 and n 3 = -1/3, it is clear from the Hamiltonian constraint (dβ + /dΩ) 2 + (dβ -/dΩ) 2 = 4 that dβ + /dΩ = 2 and the inequality for having a bounce is satisfied. More generally, when the particle moves within the wedge |β -| < -√ 3β + against the wall, we have

dβ + dΩ = 1 -3n 3 (4.13)
and the inequality dβ + /dΩ > 1 is always satisfied since n 3 < 0. Similar considerations also apply to all other walls by the triangular symmetry of the model. Thus, the point particle will always bounce against the walls of the potential.

Summarizing, close to the initial singularity, the evolution of the universe is accurately described by Kasner dynamics until the particle hits the walls and enters into a new Kasner phase (with different parameters, in general) after the bounce 8 . The bounce repeats itself again and again, in general, leading to an oscillatory behavior of S 3 that alternates its three principal axes, while the universe is approaching the initial singularity, [START_REF] Belinskii | Oscillatory approach to a singular point in the relativistic cosmology[END_REF], [START_REF] Belinskii | Construction of a general cosmological solution of the Einstein equation with a time singularity[END_REF], [START_REF] Belinskii | A general solution of the Einstein equations with a time singularity[END_REF], [START_REF] Misner | Mixmaster universe[END_REF], [START_REF] Misner | Quantum cosmology. I[END_REF], [START_REF] Misner | Classical and quantum dynamics of a closed universe[END_REF], [START_REF] Misner | Minisuperspace[END_REF], [START_REF] Misner | The mixmaster cosmological metrics[END_REF]. It appears that almost all solutions obtained by successive bounces come arbitrarily close to the corners of the well associated to the special values of Misner parameter s = 0 or s = ±3 (by triangular symmetry). Then, the space-time metric comes close to ds 2 = -dT 2 + dx 2 + dy 2 + T 2 dz 2 (up

n 1 = 2s(s -3) 3(s 2 + 3) , n 2 = 2s(s + 3) 3(s 2 + 3) , n 3 = - (s -3)(s + 3) 3(s 2 + 3
) .

An alternative parametrization has been provided by Belinskii, Khalatnikov and Lifshitz, [START_REF] Belinskii | Oscillatory approach to a singular point in the relativistic cosmology[END_REF], [START_REF] Belinskii | Construction of a general cosmological solution of the Einstein equation with a time singularity[END_REF], [START_REF] Belinskii | A general solution of the Einstein equations with a time singularity[END_REF]. The axially symmetric case n 1 = n 2 = 2/3 and n 3 = -1/3 (corresponding to s = ∞) is special as the particle heads to the corner of the triangular well after the bounce, following the canyon line β -= 0 with Kasner parameters n 1 = n 2 = 0 and n 3 = 1 (corresponding to s = 0); this is also apparent from equations (4.9), since p -= 0 all the time and p + only reverses sign relative to p Ω after the bounce. The fixed points of the bounce law arise for s = ±3 and correspond to Kasner parameters (1, 0, 0) and (0, 1, 0) that describe the evolution of the particle along the canyon lines β -= ± √ 3β + without bounce. The bounce law from the other two walls follows by permutation of the axes, which is equivalent to replacing s by (s ± 3)/(s ∓ 1), and amounts to (s + 3)/3(s -1) → 3(s -1)/(s + 3) or (s -3)/3(s + 1) → 3(s + 1)/(s -3) with fixed points s = 0 and 3 or s = 0 and -3 respectively.

to permutations of the axes), which is equivalent to the flat metric ds 2 = -dη 2 + dξ 2 + dx 2 + dy 2 by the transformation ξ = T sinhz and η = T coshz. A new era of alternating Kasner epochs subsequently starts by permuting the axes, and so on. Actually, this motion can be formulated as billiard in a finite region of hyperbolic two-dimensional space, obtained by appropriate transformation of (β + , β -) parameter space, and, as such, it is chaotic (for an overview and history of the developments in the subject see, for instance, [START_REF] Misner | The mixmaster cosmological metrics[END_REF] and references therein); it also provides the origin of the chaotic behavior exhibited by mixmaster dynamics in general, [START_REF] Barrow | Chaotic behaviour in general relativity[END_REF], [START_REF] Chernoff | Chaos in the mixmaster universe[END_REF], [START_REF] Cornish | Mixmaster universe: A chaotic Farey tale[END_REF].

The intuitive characterization of having chaos in early time cosmology is that when the universe starts with a well-defined state, it will evolve towards the singularity by going through almost all possible anisotropic stages by changing shape, as it does in general relativity.

Hořava-Lifshitz gravity

In this case, as one approaches the initial singularity, where Ω → -∞, the potential does not vanish for generic values of β ± , contrary to what happens in general relativity. Instead, it is well approximated by

V HL = V Cotton ≡ 6α κ 2 e 6Ω C ij C ij , as Ω → -∞ . (4.14)
This term can be alternatively written (in more compact form) as

V Cotton = α κ 2 ∂W ∂β + 2 + ∂W ∂β - 2 , (4.15)
where the corresponding superpotential is

W = e 3β + cosh(3 √ 3β -) -cosh( √ 3β -) -cosh(2 √ 3β -) - e -3β + cosh( √ 3β -) + 1 2 e -6β + , (4.16)
and it is positive definite when α > 0. Clearly, this applies to all models of Hořava-Lifshitz gravity, with or without "detailed balance". Consistency with the Hamiltonian constraint requires that λ > 1/3, which we assume in the following. This is not surprising in retrospect because the quadratic Cotton tensor term is marginal in the gravitational action and it is expected to dominate in the ultra-violet regime of the theory. Furthermore, since V Cotton is independent of Ω, the scale factor of the universe will evolve as a free particle with fixed (but arbitrary) momentum p Ω , so that the volume of space diminishes linearly at early times (in the appropriate time coordinate T). As for β ± that determine the shape of the universe, they will keep changing all the time following the motion of a particle in a frozen (time-independent) well

E = 1 2 dβ + dt 2 + 1 2 dβ - dt 2 + V Cotton (β + , β -) (4.17)
with fixed energy level 9

E = 1 4(3λ -1) p 2 Ω . (4.18)
Thus, the dynamics appears more complicated now, compared to general relativity, and the universe will no be -in general -in a Kasner epoch before bouncing off the walls.

Solving this effective point particle problem is not an easy task, but one can examine some qualitative features of the motion depending on the energy E. When E is very large, the potential can be approximated by zero for generic values of β ± , since the particle looks insensitive to the small bumps at the bottom of the well 10 . Thus, only in this case, which resembles general relativity, the universe will be in a Kasner epoch far away from the walls (recall that the Kasner model is a common solution to the two theories and it is insensitive to the cosmological constant in the small volume limit). Note, however, that the bounce law is modified 11 compared to general relativity and resumes its standard form, as can be seen in Appendix B. But the qualitative picture remains the same: after the bounce the universe will enter into another Kasner epoch and so on, as in general relativity.

The picture of a billiard is also very useful here. For very large E we have a very 9 We must require E > 0 so that the universe is anisotropic. Then, for small E, Ω diminishes slowly towards the initial singularity, whereas for large E it diminishes very fast. 10 Actually, this is an assumption which can be safely made only for short time development of the system, based on intuition. In general, one should also prove that these bumps, no matter how small they are compared to E, do not influence much the long time development of the system after several iterations. Such non-perturbative results, which go back to Poincaré and fall within the classic theory of Kolmogorov, Arnol'd and Moser (KAM theory), [START_REF] Kolmogorov | On conservation of conditionally periodic motions for a small change in Hamilton's function[END_REF] (but see also [START_REF] Arnol | Mathematical Methods of Classical Mechanics[END_REF] and [START_REF] Pöschel | A lecture on the classical KAM theorem[END_REF] for comprehensive exposition), are rather difficult to establish in classical mechanics and they will not be considered in detail in the present work. They should be properly investigated, however, as they might change the physical picture we are about to present after a very long time.

11 A bounce against the steady wall occurring at β + → -∞, follows the rule of ordinary reflections, namely the incidence and reflection angles are equal. It is neatly described as s/ √ 3 → √ 3/s, using the same parametrization of the Kasner parameters in terms of a single variable s as in general relativity. It should be compared to the bounce law s/3 → 3/s that governs the mixmaster universe in general relativity. In the present case, the fixed points arise at s = ± √ 3 with associated Kasner exponents

n 1 = (1 ∓ √ 3)/3, n 2 = (1 ± √
3)/3 and n 3 = 1/3, respectively. As can be seen from the corresponding expressions in Appendix B, p + = 0 in this case and the point particle moves parallel to the wall (the two signs refer to the two directions of motion), until it hits the other walls. The bounce from the other two walls follows by permutation of the axes, as usual, which is equivalent to replacing s by (s ± 3)/(s ∓ 1). Then, the bounce law yields the other fixed points s = 3 ± 2 √

3 and s = -3 ± 2 √ 3 that describe a particle moving parallel to the other two walls in either direction. Clearly, one can have a closed orbit forming equilateral triangle, which does not seem possible in general relativity. large region on the plane (β + , β -) bounded by an equilateral triangle inside which the particle moves freely with very large velocity. By simple rescaling, one can reformulate this problem as a particle moving freely inside an equilateral triangle of finite size with finite energy. The motion takes place on the plane following the standard rule of equal incidence and reflection angles. It is easy to prove that such a billiard is not ergodic. In fact, billiards in any triangular domain on the plane are non-ergodic when the angles of the triangle are rational multiples of 2π, since, then, all possible reflection angles along a given path of the particle assume only a finite set of values (see, for instance, [START_REF] Ya | Introduction to Ergodic Theory[END_REF]).

This should be contrasted with the ergodic behavior of the triangular billiard with moving walls that arises in general relativity, which can be viewed as a fixed triangular billiard with non-standard bounce law. The equivalent formulation of this problem in general relativity, as billiard in a compact domain of Lobachevsky plane with standard reflection rules, is another way of establishing ergodicity in that case, since the geodesic flow on surfaces with negative constant curvature is the prime example of ergodic behavior in Hamiltonian systems, [START_REF] Hopf | Fuchsian groups and ergodic theory[END_REF] (but see also [START_REF] Ya | Introduction to Ergodic Theory[END_REF] for an overview).

In summary, for very large E, the motion of the effective point particle in Hořava-Lifshitz gravity is not chaotic, provided, of course, that the small bumps of the potential do not affect much the long time behavior of the system after many bounces. It is an immediate and important consequence of the KAM theory that small perturbations of non-degenerate Hamiltonian systems, like ours, are not ergodic 12 . Thus, in this case, if the universe starts with a well-defined state, it will evolve towards the singularity by changing shape but without passing through all possible anisotropic stages.

On the other hand, for intermediate E, the landscape of the bottom of the potential becomes visible to the particle and the bumps can no longer be ignored even for short time. Thus, the evolution between the walls becomes more complicated and (unfortunately) it cannot be described in simple terms. Also note that the canyon lines now exhibit a small bump, as shown in Fig. 4, and, therefore, below a threshold,

E = 9α 16κ 2 , ( 4.19)
the particle cannot exit the well and head towards its corners. For E < E the motion remains bounded, provided that the anisotropies are relatively small, and the particle oscillates around the origin (fully isotropic model). Thus, for intermediate E, the dy-namics of the universe close to the singularity is very different and complex and it is not yet clear if it remains non-ergodic. If a second integral of motion exists, the system will be integrable, but we have not been able to find such thing.

Finally, for very small values E and relatively small anistropies, the particle moves around the minimum of the potential as a two-dimensional isotropic oscillator,

E = 1 2 dβ + dt 2 + 1 2
dβ - dt 2 + 81α κ 2 (β 2 + + β 2 -) , (4.20)
which follows by expanding V Cotton up to quadratic order around β + = 0 = β -. In such case the motion is integrable and chaos is obviously absent. The universe still exhibits oscillatory behavior by changing shape in its descent towards the initial singularity, but the evolution is not Kasner-like. The universe passes through the isotropic model periodically with cyclic frequency 9

√ 2α/κ.
The problem is certainly very rich and should be investigated in more detail, including numerical studies, in order to be able to make more conclusive and safe statements about chaos in the motion for general values of E. In fact, a Hamiltonian system like (4.17) can be ergodic at certain energy levels and non-ergodic at other levels. Here, there are also intermediate energies E separating into phases the behavior of mixmaster dynamics close to the initial singularity of the universe.

It is instructive to compare this behavior with the absence of chaos in fully covariant higher curvature generalizations of Einstein gravity, [START_REF] Barrow | Mixmaster cosmological model in theories of gravity with a quadratic Lagrangian[END_REF], [START_REF] Barrow | Chaotic behaviour in higher-order gravity theories[END_REF], [START_REF] Cotsakis | Mixmaster universe in fourthorder gravity theories[END_REF], where the reasoning is different. In the general context of f (R) gravity in four space-time dimensions, there is a well known conformal relationship between the vacuum higher derivative theory and ordinary general relativity coupled to a scalar field

ϕ = log f (R) (4.21)
with potential

V(ϕ) = 1 2 f 2 (R f -f ) , (4.22)
where prime denotes the derivative with respect to the four-dimensional scalar Ricci curvature. In the simplest case, the Einstein-Hilbert Lagrangian is replaced by f (R) = R + αR 2 , but f (R) can also assume more general forms. The mixmaster universe provides a consistent mini-superspace model of f (R) gravity, which is still described by a point particle that bounces off the walls of a triangular potential in the small volume limit. There is also an analogue of the Kasner solution in general relativity coupled to a scalar field, which is appropriate to use in this case. However, the effect of the scalar field ϕ is to slow down the speed of the point particle relative to the moving walls and the particle will bounce back only if it moves not too oblique relative to the walls; it should be contrasted to general relativity without scalar field, where the point particle can hit the moving wall, say the one located at β + → -∞, from anywhere within the wedge |β -| < -√ 3β + . As a result, a few collisions are sufficient to make it so oblique that it will not bounce off another wall. So, the universe will enter quickly in a definite Kasner trajectory and stay there all the time in its approach to the singularity. Thus, the evolution is not chaotic in these theories.

Conclusions

We have shown that the homogeneous cosmologies provide consistent truncations of Hořava-Lifshitz gravity in vacuum and investigated them using Hamiltonian methods with emphasis on the closed space universe of Bianchi IX type. The field equations reduce to an autonomous system of ordinary differential equations describing the motion of a point particle in a potential well with Z 3 symmetry. In general, the potential depends on time, through the volume moduli, but when the universe approaches the initial singularity it freezes, as it becomes independent of time. Then, for λ > 1/3, the universe flows to the singularity by continuous changes of its shape, as in ordinary mixmaster cosmology, rolling like a particle in the well with fixed (but arbitrary) energy E. The main difference from general relativity is that the potential does not vanish for generic values of the anisotropy parameters, and, thus, the evolution of the early universe is not described by the Kasner solution away from the steep walls. The dynamics is more intricate now, but, still, the shape of the potential far away from the origin resembles that of general relativity and the particle can bounce off the walls. In a certain limit (large E), the motion appears to be non-ergodic, and, thus, chaos is absent. The same thing applies to very small values of the energy E, though for a different reason. However, it remains to be seen if the motion is chaotic for more general values of the parameter E, as in general relativity, and compare it further with mixmaster dynamics in fully covariant higher curvature generalizations of Einstein gravity. More work is certainly required in this direction and we hope to return to it elsewhere. This work should be considered as the beginning of a more general investigation in Hořava-Lifshitz gravity. First, the most pressing open question is to revisit the singularity theorems of general relativity and examine under which general conditions they remain valid in theories with anisotropic scaling. In this context, we will also be able to see whether the matter bounce of the Friedmann model has more general value beyond the homogeneous and isotropic case. Assuming, however, that the occurrence of space-like singularities is generic in Hořava-Lifshitz gravity, it will be very important to examine how the singularity is approached by extending the standard analysis of Belinskii, Khalatnikov and Lifshitz. If the spatial points decouple from the dynamics, which is a reasonable expectation even for theories with anisotropic scaling in space and time, then, the homogeneous cosmologies considered here will prove valuable tool for understanding the behavior of the universe near the initial singularity.

Second, it is interesting to consider the canonical quantization of Bianchi IX cosmology (or any other homogeneous model for that matter) as mini-superspace models to Hořava-Lifshitz gravity. This can provide a tractable way to compare it with Einstein gravity in the quantum regime. The Wheeler-DeWitt equation for the "wavefunction" of the universe appears to be more manageable here because the walls of the effective potential are frozen in time in the domain of validity of quantum cosmology. On the contrary, in ordinary quantum gravity, the Bianchi IX model is more difficult to treat and interpret canonically because the corresponding potential is not inert to the evolution, but scales with time. We plan to address these issues in detail elsewhere.

Third, it should be noted that there is an intimate connection between the Euclidean version of Hořava-Lifshitz gravity with "detailed balance" and the theory of geometric flows. Namely, the gradient flow of the metric derived from the functional (superpotential) W yields a continuous deformation of the geometry on Σ 3 that is first order in time and trivially satisfies the higher order equations of motion of the theory. This was first pointed out in the original works, [START_REF] Hořava | Membranes at quantum criticality[END_REF], [START_REF] Hořava | Quantum gravity at a Lifshitz point[END_REF], focussing, in particular, in 2 + 1 dimensional gravity with anisotropic scaling and its connection to theory of Ricci flows on two-dimensional surfaces. In 3 + 1 dimensions, the analogous deformation theory is provided by the so called Cotton flow of three-geometries, since the variation of the Chern-Simons term in the functional W is the Cotton tensor. This, then, provides the leading order term of the flow, which is third order in space, and it should also be augmented with the Ricci curvature terms that come about by varying the three-dimensional Einstein-Hilbert term in W. Remarkably, the Cotton flow admits consistent truncation to an autonomous system of ordinary differential equations for all homogeneous three-geometries, [START_REF] Kisisel | Cotton flow[END_REF], and the same thing applies to the Ricci flow, [START_REF] Isenberg | Ricci flow of locally homogeneous geometries on closed manifolds[END_REF]. In a forthcoming paper we discuss solutions of the combined Cotton-Ricci flow for Bianchi IX geometries, which can be thought as gravitational instanton solutions in the 3 + 1 Hořava-Lifshitz gravity, [START_REF] Bakas | Geometric flows in Hořava-Lifshitz gravity[END_REF]. These configurations might also have important role in quantum gravity, in the spirit of the Hartle-Hawking proposal for the construction of the "wave-function" of the universe, using Euclidean path integrals. Thus, in this context, it is natural to expect that the theory of geometric flows will connect naturally to the problem of quantization of non-relativistic theories of gravitation, in general, and for the Bianchi IX mini-superspace model, in particular.

Finally, among other things, we mention the interesting possibility of higher dimensional generalizations. In such cases, the theory is still defined using the ADM decomposition of the metric, as in 3 + 1 dimensions, but the potential contains even higher curvature terms depending on the space dimension. For example, in 4 + 1 dimensions, we can have spatial derivative terms up to order eight, which follow from a superpotential W that involves in its integrand the square of the Weyl tensor C ijkl C ijkl and the square of Ricci curvature R 2 , assuming the "detailed balance" condition, [START_REF] Hořava | Quantum gravity at a Lifshitz point[END_REF]; one can also add the corresponding Einstein-Hilbert term to W, with cosmological constant, to insure that the theory flows to five-dimensional Einstein gravity at large scales . In such case, the Bach tensor of the four-dimensional spatial geometry, which is obtained by varying the square of the Weyl tensor in W, replaces the Cotton tensor that was featuring earlier in the 3 + 1 dimensional theory; if an R 2 term is also present in W one has to add its contribution, which scales in the same way. Then, the leading term of the potential V HL (at least when the volume of the five-dimensional universe is very small) is provided by the square of the Bach tensor, plus possible additional contributions coming from R 2 in W, and it is scale invariant. In analogy with the previous analysis, this term will dominate the cosmological evolution at early times and lead to a frozen potential well, irrespective "detailed balance", where the effective point particle rolls. Similar considerations apply to all higher dimensional cases. In view of the results obtained in the literature for higher dimensional homogeneous string cosmology models, [START_REF] Elskens | Ergodic theory of the mixmaster universe in higher spacetime dimensions[END_REF], [START_REF] Demaret | Are Kaluza-Klein models of the universe chaotic?[END_REF], [START_REF] Damour | Oscillatory behaviour in homogeneous string cosmology models[END_REF], [START_REF] Damour | Hyperbolic Kac-Moody algebras and chaos in Kaluza-Klein models[END_REF], [START_REF] Damour | Einstein billiards and overextensions of finite-dimensional simple Lie algebras[END_REF], it is also important here to explore the universal features of dynamics close to the singularity and expose their dependence on space-time dimensionality.

terms of a single variable s, following Misner, [START_REF] Misner | Mixmaster universe[END_REF], which leaves pΩ /p -invariant and changes the Kasner exponents accordingly.

n

Hořava-Lifshitz gravity:

The Hamiltonian is now derived using the asymptotic form of the potential V HL , which is independent of Ω when Ω → -∞, and becomes approximately

2H HL = p 2 + + p 2 -- 1 2(3λ -1) p 2 Ω + 18α κ 2 e -12β + . (B.10)
The point particle moves freely before and after the bounce only when p Ω is very large, in which case the Hamiltonian constraint simplifies to p 2 Ω = 2(3λ -1)(p 2 + + p 2 -) for generic values of β ± . This will be implicitly assumed here and also that α > 0 and λ > 1/3.

In the present case, there is no Ω-dependence on the wall located at β + → -∞ and it follows immediately that p Ω p - = constant . (B.11)

Then, this is an ordinary bounce from a steady wall following the standard rule that the incidence and reflection angles are equal asymptotically, namely p + flips sign and p -remains unchanged. To compare with the previous case, it is convenient to describe the free motion of the particle well before and after the bounce using Kasner exponents, which are now defined as Summarizing, the bounce in Hořava-Lifshitz gravity follows the standard reflection rule from a steady wall, whereas in general relativity this rule is modified by the moving walls and it is effectively described by inserting a factor of √ 3 in the corresponding Misner parameter. In general relativity, the standard rule of equal incidence and reflection angles only applies to the transformed Hamiltonian (B.4) from which the bounce law for the original momenta p ± was derived.

n 1 = 1 3p Ω p Ω -

Figure 1 :Figure 2 :

 12 Figure 1: The potential well of mixmaster dynamics and its three canyons

ΒVFigure 3 :

 3 Figure 3: The form of the potential V GR (β + , β -) along the canyon line β -= 0

  2(3λ -1)p + -6(3λ -1)p -, Ω -2(3λ -1)p + + 6(3λ -1)p -, + n 2 3 = 1 by virtue of the Hamiltonian constraint far away from the wall. By employing Misner's parametrization (B.7), as before, we obtain Thus, the bounce law against the wall at β + → -∞ is now described as Ω /p -invariant and changes the Kasner exponents accordingly. It has the same from as in general relativity, s/3 → 3/s, setting

	and satisfy n 2 1 + n 2 2 p Ω p -= 2(3λ -1) 3	2 n 2 -n 1	=	3λ -1 2		s √	3	+	√ s	3	= constant .	(B.13)
							√
			s √	3	→	s	3	,			(B.14)
	which leaves p s GR =	√	3 s HL .			(B.15)
	n 2 = p n 3 = 1 3p Ω 1 3p Ω p Ω + 2 2(3λ -1)p +		(B.12)

More general non-diagonal metrics of the form ds 2 = -N 2 (t)dt 2 + γ ij (t)σ i σ j are also legitimate for investigation, but they will not be considered at all in this paper.

In quantum cosmology one implements the Hamiltonian constraint by postulating the Wheeler-DeWitt equation ĤΨ = 0 for the "wave-function" of the universe with the appropriate factor ordering prescription. Some aspects of the problem have already been studied in the literature for the mixmaster model of ordinary quantum gravity,[START_REF] Hawking | The isotropy of the universe[END_REF],[START_REF] Kodama | Holomorphic wave function of the universe[END_REF],[START_REF] Moncrief | Amplitude-real-phase exact solutions for quantum mixmaster universe[END_REF], and similar considerations can also be applied to the case of Hořava-Lifshitz canonical quantum gravity. We will postpone any further discussion of the quantum aspects of early time cosmology to future publication.

This is one of the Friedmann equation in standard cosmology with p Ω being the Hubble parameter. Also, to compare with the more standard form of these equations, it is appropriate to use another time frame defined as dT = N(t)dt, where the Robertson-Walker metric takes the more familiar form ds 2 = -dT 2 + a 2 (T)dΩ[START_REF] Hawking | The large Scale Structure of Space-Time[END_REF] 3 with a = expΩ. Then, the initial singularity occurs at some finite past proper time, say T = 0, instead of t = -∞.

This condition by itself is quite restrictive since it rules out regular radiation before the bounce.

Another example for having an initial singularity -rather than a bounce -is provided by the anisotropic Bianchi I model (also known as Kasner solution), although the reasoning is slightly different here. This is a common solution to general relativity and Hořava-Lifshitz gravity in vacuum because Σ 3 is a flat three-dimensional space and all components of the Ricci curvature and Cotton tensor vanish. In this case, the only contribution to the Hamiltonian constraint (neglecting the cosmological constant term for small volume) is provided by the kinetic energy of the anisotropy parameters and the universe can evolve towards the initial singularity without reaching a minimum volume. Thus, the bounce in the Friedmann model does not appear to be generic, in particular in the presence of anisotropy.

The rate dβ + /dΩ is positive because β + decreases as Ω decreases while the particle is heading towards the wall in its descent towards the singularity.

The bounce law for the Kasner parameters, n i → n i , has been worked out in the literature (but see also Appendix B). For a bounce against the wall at β + → -∞, it is most easily described as s/3 → 3/s, using Misner's parametrization,[START_REF] Misner | Mixmaster universe[END_REF],

At a given energy level, however, the invariant non-resonant (Kronecker) tori in phase space form a Cantor set, which has no interior points. Therefore, it is impossible to tell with finite precision whether a given initial position falls on an invariant torus or in a gap between such tori. In such cases, according to KAM theory, one can only make probabilistic statements for a chosen orbit to be on an invariant torus, and, hence, be stable.

Acknowledgements

This work was partially supported by the Cluster of Excellence "Origin and the Structure of the Universe" in Munich, Germany, the French Agence Nationale pour la Recherche programme "String Cosmology" under contract 05-BLAN-NT09-573739, the ERC Advanced Grant "Mass Hierarchy and Particle Physics at the TeV Scale" under contract 226371, the ITN programme "Unification in the LHC Era" under contract PITN-GA-2009-237920 and the bilateral programme GRC APIC PICS-Grèce for French-Greek scientific collaboration under contract 3747. The authors I.B., F.B. and M.P. are grateful to the hospitality extended to them at the Arnold Sommerfeld Center for Theoretical Physics in Munich during the early stages of this work. We also thank Anastasios Petkou and Marika Taylor for stimulating discussions during the GGI workshop "Perspectives in String Theory" that triggered parts of the present work as well as Elias Kiritsis, Hermann Nicolai and George Savvidy for fruitful exchanges. Finally, we thank Christos Sourdis for assistance in drawing the figures.

A Bianchi IX model geometry

The Bianchi IX model describes a homogeneous (but generally non-isotropic) threedimensional geometry with the topology of S 3 and isometry group SU [START_REF] Hawking | The large Scale Structure of Space-Time[END_REF]. The line element is constructed using the corresponding left-invariant 1-forms σ i ,

σ 2 = cosψsinθdφsinψdθ , (A.2)

and takes the form

The 1-forms σ i satisfy the defining SU(2) relations

whereas the angles range as 0 ≤ θ ≤ π, 0 ≤ φ ≤ 2π and 0 ≤ ψ ≤ 4π, since ψ is extended to the double covering of the rotation group. The space integration is carried out using

Next, we present the expressions for the Ricci curvature and Cotton tensors of Bianchi IX metrics that will be used in the main text. Proper discussion requires the use of e i = √ γ i σ i and the corresponding connection 1-forms ω i j satisfying the zero torsion relations de i + ω i j ∧ e j = 0. Then, the curvature 2-forms are computed as

and the Ricci 1-forms are (Ric

The Ricci curvature scalar is R = i k (Ric) k . Also, the Cotton 2-form is given by

where Y j are simply given by

With these explanations in mind, the non-vanishing components of the Ricci curvature tensor take the following form

and, therefore, the Ricci scalar curvature is

Also, the non-vanishing components of the Cotton tensor take the form

) 17) It is convenient (and will be used throughout this paper) to parametrize the metric coefficients γ i as follows,

)

The volume of S 3 is parametrized by Ω, whereas β + and β -measure the deviations from the isotropic metric that is associated to β + = 0 = β -. Thus, for γ 1 = γ 2 = γ 3 , the metric on S 3 is homogeneous, but not isotropic, having different circumference on great circles in each of the three mutually orthogonal principal directions. Axially symmetric non-isotropic metrics are obtained by choosing one of the deformation parameters equal to zero, say β -= 0, in which case

Finally, note that the Cotton tensor vanishes in the isotropic case γ 1 = γ 2 = γ 3 , since the metric of the round S 3 is conformally flat.

B The bounce law

In this appendix, we derive the bounce law of the effective point particle as it hits the wall located at β → -∞. Far away from the wall the potential is assumed to vanish and the particle moves freely, where this is applicable.

General relativity:

In this case, the Hamiltonian is derived using the asymptotic form of the potential V GR , as Ω → -∞, and becomes approximately

The Ω-dependence of the potential can be easily transformed away by introducing new variables

and their conjugate momenta

Then, the Hamiltonian takes the simpler form

In terms of these variables, the potential is frozen in time and both pΩ and p -are constants of motion. Dividing by parts, one obtains

relating the original momenta before and after the bounce.

At this point, it is useful to introduce Kasner exponents, which are applicable to the evolution of the universe well before and after the bounce,