Singularity avoidance in Bianchi I quantum cosmology

We extend recent discussions of singularity avoidance in quantum gravity from isotropic to anisotropic cosmological models. The investigation is done in the framework of quantum geometrodynamics (Wheeler-DeWitt equation). We formulate criteria of singularity avoidance for general Bianchi class A models and give explicit and detailed results for Bianchi I models with and without matter. We find that the classical singularities can generally be avoided in these models.


I. INTRODUCTION
A major issue in any quantum theory of gravity is the fate of the classical singularities. So far, such a theory is not available in final form, although various approaches exist in which this question can be sensibly addressed [1,2]. It is clear that such an investigation cannot yet be done at the level of mathematical rigor comparable to the singularity theorems in the classical theory (see e.g. [3]). Nevertheless, focusing on concrete approaches and concrete models, one can state criteria of singularity avoidance and discuss their implementation. This is what we shall do here.
We restrict our analysis of singularity avoidance to quantum geometrodynamics, with the Wheeler-DeWitt equation as its central equation [1]. Although this may not be the most fundamental level of quantum gravity, it is sufficient for addressing the issue of singularity avoidance. Quantum geometrodynamics follows directly from general relativity by rewriting the Einstein equations in Hamilton-Jacobi form and formulating quantum equations that yield the Hamilton-Jacobi equations in the semiclassical (WKB) limit. It thus makes as much sense to addressing singularity avoidance here than it does to addressing it in quantum mechanics at the level of the Schrödinger equation. Singularity avoidance has also been discussed in loop quantum gravity [2,4,5], with various results, but we will not consider this here.
Singularity avoidance was already addressed by De-Witt in his pioneering paper on canonical quantum gravity [6]. He suggested to impose the condition Ψ → 0 for the quantum-gravitational wave functional Ψ when approaching the region of a classical singularity. The wave functional is effectively defined on the configuration space of all three-dimensional geometries, also called superspace [1,7]. The "DeWitt criterion" of vanishing wave function then means that Ψ must approach zero when approaching a singular three-geometry (which itself is not part of superspace, but can be envisaged as its boundary). It is important to emphasize that this criterion is a sufficient but not a necessary one: singularities can be avoided for non-vanishing Ψ (recall the solution of the Dirac equation for the hydrogen atom, which is even diverging).
DeWitt had in mind cosmological singularities such as big bang or big crunch. The question arises, of course, also for the singularities that classically arise for gravitational collapse. In simple models of quantum geometrodynamics, their avoidance can be rigorously addressed. One example is the collapse of a null dust shell, which classically develops into a black hole, but quantum gravitationally evolves into a re-expanding shell, with Ψ = 0 in the region of the classical singularity [8]. In general, however, such cases are too difficult to allow a mathematically exact treatment, so most investigations so far were restricted to Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology. The first detailed discussion of singularity avoidance in Wheeler-DeWitt quantum cosmology was performed for the big-rip singularity that occurs in the presence of phantom matter [9]. Other applications followed, which generally focused on singularities occurring in dark-energy models but also on the big bang; see, for example, [10][11][12][13] and the references therein. The question was also investigated for f (R) quantum cosmology [14]. In most cases, the De-Witt criterion was applied.
In the present paper, we make a step forward and discuss the issue of singularity avoidance for anisotropic cosmologies. The simplest case is the Bianchi I model (see e.g. [15]), to which we restrict our investigation here. The more interesting Bianchi IX model is reserved for a future investigation. Anisotropic models are characterized by the fact that the dimension of their configuration space (minisuperspace) is bigger than two already for the pure gravitational case. This will be important for the formulation of the DeWitt criterion.
Our article is organized as follows. In Sec. II, we formulate our criteria for singularity avoidance. In this, a generalization is made that takes into account the conformal structure of minisuperspace. Section III then addresses the vacuum Bianchi I (Kasner) model. Sections IV and V are devoted to Bianchi I models with matter: an effective matter potential is used in Sec. IV, and a dynamical (phantom) scalar field is used in Sec. V. We shall find that singularities can be avoided in all relevant cases. Section VI presents a short conclusion and an outlook.

II. CRITERIA FOR SINGULARITY AVOIDANCE
In this section, we formulate the criteria of singularity avoidance at the level of a general (diagonal) Bianchi class A model. These will be applied in detail to the Bianchi I model in the following sections.
The action for such Bianchi models can be brought into the form [15][16][17] (1) We parametrize the minisuperspace M of these models by the coordinates q = {α, β + , β − , φ}, where α ≡ ln a, β + , and β − are the Misner variables, 1 and φ denotes matter field degrees of freedom; (3) R is the three-dimensional Ricci scalar. Units are chosen such that 3c 6 V 0 /4πG = 1, where V 0 is the volume of three-dimensional space (assumed to be compact here).
Variation with respect to the lapse N yields the Hamiltonian constraint where the p I are the momenta canonically conjugate to the configuration variables q, the G IJ denote the components of the DeWitt metric, G IJ the components of its inverse, and V is the minisuperspace potential which contains contributions from the three-curvature and from the matter part. We remark that the equations of motion can be written as a geodesic equation plus a forcing term ( [17], p. 452). Because of the constraint nature (2) of the Hamiltonian, minisuperspace possesses a natural conformal structure. This can be seen as follows. Let us consider a rescaling of the lapse, N → N = Ω 2 N , with a differentiable function Ω : M → R + . The transformation of the Hamiltonian constraint then follows from the invariance of the total Hamiltonian H according to This rescaling induces a local Weyl transformation of the DeWitt metric, We can thus interpret minisuperspace as a conformal manifold (M, [G IJ dq I ⊗ dq J ]), that is, as a manifold equipped with an equivalence class of metrics, 1 But note that Misner in [17] uses Ω = −α.
Objects of interest on such a manifold are conformally covariant objects, for example tensors that transform like under Weyl transformations. We call k the conformal weight of T and denote it by w(T ) = k. Because of the conformal nature of configuration space, a spacetime which satisfies Einstein's equations can, in fact, be regarded as a sheaf of geodesics on this space [18]. In geometrodynamics, quantization is performed formally by replacing the canonical momenta according to the rule p I → −i ∂ ∂q I and substituting these expressions into the Hamiltonian constraint (2) [1]. This procedure leads to the minisuperspace Wheeler-DeWitt equation where the quotation marks indicate the need for choosing an appropriate factor ordering. The underlying conformal structure of minisuperspace motivates us to choose a factor ordering that makes the Wheeler-DeWitt equation conformally covariant. Following the discussion by Misner ([17], p. 462), 2 this is achieved by where R denotes the Ricci scalar constructed from G IJ and ξ = d−2 4(d−1) , with d = dim(M). If we, in addition, impose the weights w(Ψ) = −(d − 2)/2 and w(V) = −2, the Wheeler-DeWitt equation (7) is indeed conformally covariant. The operator − ξR is called the conformal Laplacian or Yamabe operator. It was shown that, given a compact Riemannian manifold of dimension d ≥ 3, one can find a metric conformal to G IJ with constant scalar curvature [20].
We are now in a position to formulate possible criteria for singularity avoidance. As mentioned in the Introduction, the first one goes back to DeWitt [6], who suggested to take Ψ → 0 in the vicinity of the region where a classical singularity occurs as a sufficient criterion for quantum avoidance. This criterion was successfully applied to cosmological models in a series of recent papers; see, for example, [13] and references therein. These examples deal mostly with two-dimensional minisuperspaces where w(Ψ) = 0 and the usual Laplace-Beltrami operator coincides with the conformal Laplacian. In dimensions d ≥ 3, however, the DeWitt criterion is not conformally invariant, which is why we seek here a generalization to guarantee this invariance.
This leads us to consider conformally invariant objects constructed from Ψ. We first note that we can define a density of conformal weight 0 by where dvol contains the square root of the (absolute value of the) determinant of the DeWitt metric, and ⋆ denotes the Hodge star. Moreover, we address the Klein-Gordon current defined by which is a (d − 1)−form with conformal weight 0. These definitions allow us to propose the following two criteria.
This is the conformally invariant version of the DeWitt criterion [6].
Another criterion, which was introduced in the discussion of the quantum fate of the big-rip singularity in [9], is the following: Criterion 3. A wave packet is said to avoid the singularity if it spreads in the vicinity of the singularity.
The spreading of wave packets indicates the breakdown of the semiclassical approximation. Classical cosmology and in particular the classical singularity theorems then cease to apply. This criterion is indeed fulfilled in the big-rip case of [9]. In a sense, the DeWitt criterion and its conformal generalization also obey criterion 3 because the vanishing of the wave function signals the breakdown of the classical approximation.
We note that the second criterion suffers from the problem that it is not applicable in the case of real wave functions, which often arise as solutions to the (real) Wheeler-DeWitt equation. 3 The Klein-Gordon flux is not positive definite and can thus in general not be interpreted as a probability flux. Exceptions are situations where only one WKB branch is present; this has led to the proposal that the Klein-Gordon current only be applied to such cases [21]. In the following, we will thus mainly concentrate on the first criterion, which is the natural generalization of the DeWitt criterion to higher-dimensional minisuperspaces.
Application of these criteria is based on the idea that the (square of) the wave function is related to probability, as is the case in quantum mechanics. In quantum gravity, this is far from clear [1]. The main reason is the absence of an external time parameter in the Wheeler-DeWitt equation. Only in the semiclassical (Born-Oppenheimer) approximation, where an approximate time parameter emerges, can one impose the usual probability interpretation. Nevertheless, we shall stick heuristically to this idea also in the full theory. Peaks of the wave function have often been interpreted as giving predictions in cosmology; see, for example, [22] and the references therein.
In the semiclassical limit with only one WKB component, an interpretation of the use of the probability interpretation in minisuperspace was suggested in [23]; see also [24], pp. 186-190. It is thus appropriate to end this section with some remarks on the formulation of this proposal in the language of conformal minisuperspace.
Let us consider solutions of the Wheeler-DeWitt equation in the WKB approximation given by Ψ ≈ √ D e i/ S, where S is a solution to the Hamilton-Jacobi equation and D is the van Vleck factor which satisfies the linear transport equation Let now A ⊆ M be a region in minisuperspace and B a thin 'pencil' drawn out by the classical solutions, that is, integral curves of the vector field G IJ (∂ I S)∂ J . It was shown in [23] that where is the conformally invariant and conserved flux through a hypersurface Σ crossing the pencil B. 4 The contribution of B to A∩B ⋆|Ψ| 2 is therefore proportional to the coordinate-time that the classical solutions filling out the pencil B spend in the region A. Note that w(⋆D) = w(⋆|Ψ| 2 ) = 2. This reflects the fact that the integral N dt on the right-hand side of (12) depends on the representation of the lapse. The formula can help us to interpret the behavior of wave packets in regions of minisuperpsace where the WKB approximation is valid.

III. KASNER SOLUTION
The vacuum Bianchi I (Kasner) solution can be written in the form 4 It is assumed that the hypersurface is chosen such as to cross each classical trajectory in the pencil only once [21,24].
The constraints on p x , p y and p z define the so-called Kasner sphere and Kasner plane, respectively. The physical solutions lie on their intersection, which represents a circle in the (p x , p y , p z ) space. The nature of the singularity depends on the value of the coefficients p x , p y , p z . If one of them is equal to 1, the Kasner solution will become the Milne universe, which is diffeomorphic to slices of Minkowski spacetime. The singularity is then only a coordinate singularity. For all other values, the singularity is physical, which is indicated by the divergence of the Kretschmann invariant, Since R µν = 0, the curvature singularity is a pure Weyl singularity. If we use the common parametrization of the Kasner circle, with u ∈ (−∞, ∞), we find that In order to obtain a Hamiltonian for the model, we employ the general symmetry reduced ansatz with the scale factor a = (a x a y a z ) 1/3 =: e α , whose third power describes the volume (which expands as a 3 = t), and with the anisotropy factors β ± , which describe the shape of the universe. Note that the scale factor is chosen here to be dimensionless; the physical length dimension is in the coordinates x, y, and z.
The symmetry reduced Einstein-Hilbert action takes the form The Hamiltonian obtained after the usual Legendre transform reads Choosing for the lapse function the value N = e 3α , it becomes clear that the Hamiltonian is equivalent to the Hamiltonian of a free relativistic particle in 2 + 1 dimensions. We conclude that the solutions represent straight lines in minisuperspace, which can be parametrized as follows: with C ± ∈ R arbitrary constants. The approach to the singularity is called velocity term dominated (VTD); see, for example, [25]. This terminology refers to the dominance of the kinetic over the potential terms, which is trivially fulfilled here (absence of potential). The DeWitt metric on M is given by from which we obtain for the Ricci scalar on M the value R = 9 2 e −3α . The Wheeler-DeWitt equation now reads where the numbers f and ξ parametrize a family of operator orderings. After the transformation Ψ → Ψ = e f α Ψ we obtain The conformal factor ordering is obtained by setting f = 3/4 and ξ = 1/8. We then get In the conformal factor ordering the Wheeler-DeWitt equation is thus identical to the classical wave equation in d = 1 + 2 dimensions. Note that the DeWitt metric is flat in this representation, such that criterion 1 above is equivalent to the DeWitt criterion Ψ → 0 as applied in earlier papers; see, for example, [13].
Let us now turn to the formulation of the criteria for singularity avoidance. There, the minisuperspace dimension d will be crucial. Solutions to the free wave equation in 1+1 dimension can propagate only into two directions. Wave packets are not subject to spreading and their amplitudes do not decay. In higher dimensions, however, the wave can propagate into infinitely many directions. This leads to a spreading and a resulting decay of the amplitude of the wave. The above statement can be made more precise in the form of decay rate estimates.
In d > 2 dimensions, we can apply the following decay rate estimate (see e.g. [31]): Let Φ be a solution to the initial value problem where f and g are smooth functions R d−1 → R with compact support. Then there exist C 1/2 > 0 such that For the situation in question it follows that such wave packets satisfy the above criteria 1, 2, and 3 for singularity avoidance, that is we have, J[Ψ, Ψ] → 0 and ⋆ |Ψ| 6 → 0 as α → ±∞. (24) This is caused by the spreading of the wave packet when approaching the region of the classical singularity. In addition to the avoidance of the initial singularity, there is also an avoidance of the non-singular late stages of the universe for α → +∞, since the wave function goes to zero there, too. This is in conflict with the classical behavior and demonstrates that quantum effects can become relevant at arbitrary scales. Here, they "stop" the classical evolution of the model. Other examples are quantum effects near the turning point of a classically recollapsing universe [26] or the above-mentioned cases of future singularities.
In the next section, we will investigate if and how the situation changes if matter is added to the model.

IV. BIANCHI I MODEL WITH AN EFFECTIVE MATTER POTENTIAL
In this section, we treat matter in a phenomenological way. The representation of matter by a dynamical scalar field φ is relegated to the next section. In anisotropic models, anisotropic pressures can be used, but we address for simplicity the case of a barotropic fluid.
A hypersurface orthogonal (non-tilted) barotropic fluid with an equation of state p = wρ and ρ ∝ a −3(1+w) can be modelled by adding an effective matter potential of the form V(α) = N V 0 e −3(1+w)α ∝ ρ to the Einstein-Hilbert action (18), with V 0 > 0 being constant. The full action then reads We restrict our discussion to w < 1, which excludes the case of a stiff matter fluid. The important cases of a cosmological constant (w = −1), dust (w = 0), and radiation (w = 1/3) are included. The null and weak energy conditions are satisfied for w ≥ −1, while the strong energy condition and the dominant energy condition require w ≥ − 1 3 and −1 ≤ w ≤ 1, respectively. The variables β ± are cyclic and we call their conjugate momenta p ± . Variation of the Lagrangian with respect to N leads tȯ where k := 3(1 − w). We assume that p 2 + + p 2 − = 0 and choose the comoving gauge N = 1. Equation (25) is then solved by where 2 F 1 [a, b; c; z] is the hypergeometric function. The scale factor a(t) is shown for different w in Fig. 1a.
For small a, the hypergeometric function asymptotically equals 1, and we get for a → 0: Thus the universe starts with a big bang at t = 0, independent of the value for the barotropic index w. For large a and w = −1, the hypergeometric function can be simplified, too, and one gets from (26) in the limit a → ∞: For k < 6 (w > −1), the universe expands infinitely, whereas in the phantom case, that is for k > 6 (w < −1), the universe becomes infinitely large already at t = t * and ends with a big rip. We note that (28) is the full solution for the flat FLRW case: for k < 6, there is a big bang, but for k > 6 there is no past singularity. Therefore one can say that the anisotropy introduces the past singularity, leading to a model with big bang and big rip. For the anisotropy factors one has which becomes constant for large a, see also Fig. 1b. Thus in contrast to the vacuum solution this universe isotropizes for late times. For small a, the asymptotic behavior corresponds to (20), which is again independent of the matter content. This property is sometimes called "matter doesn't matter". Since the Kasner behavior is recovered in the limit a → 0, this approach to the singularity is referred to as asymptotically velocity term dominated (AVTD) [25]. We now turn to the quantum version of these models. The Wheeler-DeWitt equation reads where we have set = 1 and skipped the tilde over the wave function. The solutions can be written in the form with the mode functions given by where J ν (z) and Γ(z) denote the Bessel function of the first kind and the gamma function, respectively. Let us now investigate the asymptotic forms of the wave packet. In the limit α → −∞ we can approximate the mode functions by We conclude that the quantum Kasner behavior is recovered in this limit (which follows as a solution of (23)). The discussion of the limit α → ∞ is slightly more complicated, but it turns out that a discussion of the mode functions in the WKB approximation will be sufficient. A solution to the Hamilton-Jacobi equation is given by The corresponding van Vleck factor reads D p+,p− (α) = 1 If we introduce the functions then the approximate wave packet with these coefficients, matches the exact wave packet for large α at the leading order. This follows from the asymptotic expansion of the exact mode functions and an approximation of the WKB modes of the form Then one has We can now draw a clear picture of the behavior of wave packets. In the limit α → −∞, we recover the quantum Kasner behavior. Consequently, we expect a spreading with a resulting decay of amplitudes. The behavior in the limit α → ∞ can be infered from (35): the term in the second line of this equation is just the Fourier transform of B σ and is independent of α. If, for example, we choose B σ to be Gaussian, its Fourier transform will be a Gaussian which is peaked around some particular values of β + and β − . This strongly reflects the classical behavior of isotropization. Most importantly, wave packets do not spread in the region where α is large. The wave packet is modulated by a strongly oscillating factor and an exponentially decaying factor. The exponentially decaying factor comes from the van Vleck factor (33) and can be interpreted as arising from the particular representation of the wave function.
The decay of the mode functions in this representation can be intuitively understood by inspecting the Hawking-Page formula (12): The representation of the wave function Ψ we are working with is related to the gauge N = e 3α by the corresponding representation of the DeWitt metric. In this gauge, classical solutions reach α = ∞ in a finite time t. Hence they spend less and less time t in the region of minisuperspace where α is large. In this sense the decay of the density √ −GD is implied by (12).
For simplicity, we now set B − = 0. Then the large-α limit of the Klein-Gordon current is given by Up to leading order, the current only has an α component given by the Fourier transform of B + (p + , p − ). If we assume that B + is peaked at some particular values p + and p − , we will expect the Fourier transform of B + to be peaked at some particular value of β + and β − . The current thus reflects the classical behavior in the region where α is large (in contrast to the vacuum Kasner case). We have, however, ⋆|Ψ| 6 → 0 as α → ∞. Note that the behavior is qualitatively independent of w, that is, there is no difference between the cases w ≥ −1 and w < −1, although the latter case leads to a big rip. The big rip is thus only avoided by criterion 1.

V. BIANCHI I MODEL WITH A PHANTOM FIELD
In the previous section, we have added matter degrees of freedom through an effective potential V(a). Here, we will instead implement the equation of state p = wρ by a scalar field φ with a potential V (φ), a procedure described in [30]. Matter degrees of freedom are now dynamical. The connection between the kinetic and potential terms and the parameters ρ and p are as follows: Note that l = ±1 depending on whether we consider normal or phantom matter, respectively. Using these relations one finds the same functions for a(t) and β ± (t) as in the previous section. We use here κ ± instead of p ± , because these constants will be used in the construction of V (φ). Combining (36), (25), p = wρ, and ρ(a) we get for the classical solution in configuration space, see Fig. 2a. The scalar field vanishes like a polynomial for small a and diverges logarithmically for large a.
Using the same equations as before we get for the potential (Recall k = 3(1 − w).) After substituting a by φ and using (37) we find compare Fig. 2b. Potentials with sinh-functions also occur frequently in FLRW models [27,28].
As the Wheeler-DeWitt equation will not be analytically solvable for a general potential, we choose here k = 12 (w = −3, l = −1) as a particular example. This is, on the one hand, simply solvable and reflects, on the other hand, the general case. We approximate V (φ) and |φ(a)| for large a and therefore large |φ|; that is, we investigate the limit when approaching the big rip. This Let us now turn to quantum cosmology. Note that the DeWitt metric has here signature (−, −, +, +), since the kinetic term of the (phantom) scalar field has the same sign as the one of the scale factor (cf. Eq. (30) in [9]). The conformally covariant Wheeler-DeWitt equation with the scalar potential in the limit approaching the big rip reads where V 0 := 2ρ 0 (κ 2 + + κ 2 − ). After intermediate steps in which one makes use of the variables u := α + |φ| and v := α − |φ|, we solve the equation using a separation ansatz. Then the full solution is with mode functions where H (1,2) ν (z) are the Hankel functions. Note that the latter assume the WKB form for large arguments, where for the van Vleck determinant D ∝ e −3(α+|φ|) holds. From (39) we see that |φ| ∝ α for large α. Thus the amplitude of the wave function decreases and the wave function vanishes as we approach the big rip singularity. As in the previous section, the DeWitt criterion is fulfilled, and the singularity is avoided if criterion 1 is adopted.
Using the asymptotic WKB form of the Hankel functions, we get a WKB solution for the complete mode wave function. The phase is Using the principle of constructive interference [7,29], we get As expected, we find that |φ| ∝ α + const. and that β ± become constant. Note that we have not recovered here (29), because we have used an approximated form of the potential and an asymptotic expression of the wave function.    In this paper, we have extended previous investigations of singularity avoidance from isotropic to anisotropic models. We have, in particular, adapted the avoidance criterion to the covariant structure of minisuperspace, which becomes relevant for dimensions higher than two. We have found that the DeWitt criterion can, in general, be fulfilled, but not so the vanishing of the Klein-Gordon current. For the reasons mentioned, however, we attribute more relevance to the DeWitt criterion.
While the general criteria were formulated for general Bianchi class A models, detailed investigations were made for the Bianchi I model with and without matter. Bianchi I models admit the prototype of an asymptotically velocity term dominated (AVTD) model. Our results of singularity avoidance should thus be representative for such a kind of singularity with a sufficiently large number of degrees of freedom. Other Bianchi models such as Bianchi VIII and Bianchi IX exhibit an oscillatory behavior when approaching the singularity. The singularity can, however, become AVTD if, for example, a scalar field is added [25].
Because Bianchi IX models are generally considered as reflecting the generic behavior towards a cosmological singularity, future investigations of singularity avoidance should attempt to address these models in as much detail as possible. For this, it would be desirable to have mathematical theorems available such as those discussed here for the Kasner solution. In [32] one finds an existence and uniqueness theorem (Theorem 8.6 there), but for the Bianchi IX potential no decay rate estimates seem to exist. The Wheeler-DeWitt equation for the vacuum Bianchi IX (mixmaster) model was solved numerically in [33] by using the 'hard wall approximation'. The results found in this analysis strongly indicate the decay of wave packet amplitudes.
It is generally believed that the approach to a spacelike singularity at the level of the full Einstein equations can be described by the Belinsky-Khalatnikov-Lifshits (BKL) scenario; see, for example, [25,34] and references therein. This corresponds to a decoupling of spatial points, in which every spatial point exhibits the dynamics of a seperate Bianchi IX model. The eventual goal will be to present a quantum-gravitational analysis of this situation, from which one should be able to draw general conclusions about singularity avoidance. An investigation in the framework of affine quantization was made recently in [35]. Attempts in this direction using the Wheeler-DeWitt equation will be the subject of future investigations.