Numerical Techniques in Loop Quantum Cosmology

In this article, we review the use of numerical techniques to obtain solutions for the quantum Hamiltonian constraint in loop quantum cosmology (LQC). First, we summarize the basic features of LQC, and describe features of the constraint equations to solve - generically, these are difference (rather than differential) equations. Important issues such as differing quantization methods, stability of the solutions, the semi-classical limit, and the relevance of lattice refinement in the difference equations are discussed. Finally, the cosmological models already considered in the literature are listed, along with typical features in these models and open issues.


Introduction
One of the main open questions for any feasible theory of quantum gravity is that of resolution of classical singularities. Due to the lack of experimental verification in this context currently, the resolution of classical singularities as well as the recovery of the appropriate semi-classical behavior are the key check that any acceptable theory of quantum gravity must overcome. Therefore, testing these two issues has been the main goal of all the numerical work performed in the symmetry reduced context of loop quantum cosmology (LQC).
The geometrodynamical quantization of general relativity [48] is assumed to be the correct semi-classical limit of any quantum theory of gravity. In this context, the metric (or in the symmetry reduced cosmological context, the scale factor) is understood as the fundamental variable. Nevertheless, this quantization can not be fundamental because the quantum states completely follow the classical trajectories of the universe in the phase space. More precisely, if a state representing the current classical universe is evolved backwards making use of the Wheeler-DeWitt (WdW) equation, it ends up in the big bang singularity. As we will explain below, the situation in LQC is completely different and, in all analyzed models, the resolution of the corresponding singularity has been reported.
This article is organized as follows: In Section 2, we provide a brief introduction to LQC with emphasis on the different quantization schemes used. Section 3 discusses (numerical analysis This paper is a contribution to the Special Issue "Loop Quantum Gravity and Cosmology". The full collection is available at http://www.emis.de/journals/SIGMA/LQGC.html inspired) stability properties of the equations that arise in LQC models. In Section 4, we survey a number of isotropic and non-isotropic models. We end this article with a brief summary in Section 5.

Essentials of loop quantum cosmology 2.1 Phase space variables and quantum operators
In this section, we will develop the general setting of LQC models -in particular, the origin of the Hamiltonian constraint difference equation. This will allow us to make some generic comments about solving these constraints using numerical techniques. In LQG, an oriented spatial manifold M is chosen, with a fiducial set of spatial triad vectors e a i mapping from the internal vector space (with a positive definite metric q ij ) to the tangent space (with indices a, b, . . . ) of the manifold M 1 . The dual of the triads are the 1-forms ω i a , which give a positive definite metric q ab = q ij ω i a ω j b on M . For the Bianchi models, these 1-forms satisfy the Maurer-Cartan relation dω i = (1/2)C i jk ω j ∧ ω k , where C i jk are the structure constants of the Lie algebra for the particular Bianchi system.
For the following, we consider specifically the case of a type I Bianchi model [9], but similar results hold for other systems as well. In particular, we use a diagonal form of Bianchi I, so that the physical metric is given in terms of three scale factors a i (τ ) and the lapse N by ds 2 = −N 2 dτ 2 + a 2 1 (τ )dx 2 1 + a 2 2 (τ )dx 2 2 + a 2 3 (τ )dx 2 3 .
As with all Bianchi cosmologies, the spatial slices are non-compact and all physical quantities are homogeneous. This raises the need to restrict all integrations to a fiducial cell V, where all edges of the cell align with the coordinate axes x i and with a chosen metric q ab . For Bianchi I, this metric can be chosen to be flat, and we designate the sides of the cell along the x i coordinate to be L i ; thus, the volume of the cell V is given by V 0 = L 1 L 2 L 3 . Then, we can write the triads and connection components as where q is the determinant of the fiducial metric q ab . Throughout this paper, there is assumed to be no summation if all repeated indices are covariant or contravariant; the usual Einstein summation convention will hold when there are mixed indices. Thus, there is no summation in either equation listed in (2.1). The phase space of the reduced symmetry system is coordinatized by six variablesp i ,c i , satisfying the Poisson bracket relation Here γ is the Barbero-Immirzi parameter, a quantization ambiguity which has been calculated to be γ 0.2735 by matching with the Bekenstein-Hawking black hole entropy formula [55]. However, there is a freedom to re-scale 2 the coordinates independently in Bianchi I as x i → α i x i , under which ω i a → α i ω i a and e a i → (α i ) −1 e a i (no summation on i. Under this rescaling, the variablesp i ,c j will similarly change -an undesirable feature, since the quantization 1 In most works on LQC, these fiducial inputs are designed as, e.g. o e a i , to differentiate them from the physical analogues e a i , which includes measurable functions such as the scale factors ai(τ ). However, we do not need this distinction, so we simplify our notation accordingly. For more discussion, see [5,9]. 2 Although we do not address it explicitly here, there is also the issue of the behavior of the system under parity transformations [6,10]. The appropriate choice (when there are no fermions) allows an effective reduction in the parameter space m over which the wave function is defined, e.g. only mi > 0 for all three parameters mi. of the model should not depend on the fiducial structure chosen. Thus, we normalize the variablesp i ,c j such that p 1 = L 2 L 3p1 , p 2 = L 1 L 3p2 , p 3 = L 1 L 2p3 , and c 1 = L 1c 1 , c 2 = L 2c 2 , c 3 = L 3c 3 .
Another way to say this is that the presence of the fiducial volume V 0 gives a one-parameter ambiguity in the symplectic structure (2.2), which is now fixed by the choice of c i , p j . Therefore, the Poisson brackets become If the structure constants C i jk are non-zero -as in the other Bianchi cosmologies -there is the need to normalize them accordingly [10].
In the full theory, the variables are the fluxes of the triads E a i and the holonomies h e , defined by a connection A i a along an edge e, where the latter is given by h e (A) = P exp( e A i a τ iė a dt). Here, τ i are the generators of the SU(2) Lie algebra, such that τ i τ j = (1/2) k ij τ k −(1/4)δ ij I (these are related to the Pauli matrices σ k as τ k = −iσ k /2), and I is the unit 2 × 2 matrix. The fact there are no quantum operators corresponding to the connection components c i in LQG is an indicator of the diffeomorphism invariance of this theory. To match LQG as much as possible, we exponentiate the connection components appropriately to get the holonomy h (δ) k along the coordinate axis x k with an edge of length δL k , given by Here δ is a constant value and the holonomy is given in terms of the almost periodic functions exp(iδc k ) -"almost" because δ can be any real number, not just an integer. We consider a Hilbert space basis of eigenstates |m 1 , m 2 , m 3 ≡ | m of the momentum operator p i . Here, the values m i ∈ R may either be the momentum eigenvalues themselvesi.e. p i ( m) = m i , and therefore have the appropriate physical dimensions -or dimensionless quantities, so that the dimensions are carried within the function p i ( m), depending on the choice made in a particular physical model. These eigenstates satisfy Because there is a Kronecker delta on the right-hand side, rather than a Dirac delta distribution, the wave function is made of a countable number of these orthonormal basis states, in the following form where the usual finiteness condition is assumed. This quantization is in-equivalent to the standard Schrödinger form, and thus to the usual quantum cosmology using the Wheeler-DeWitt equation as a starting point. When the triad and holonomy are promoted to quantum operators, they act on the states aŝ for a constant k and similarly for the other holonomy operators. The function p i ( m) and k must combine to satisfy the Poisson bracket relation (2.3); for Bianchi I, the triad eigenvalue p i ( m) = m i and k = −8πGγ . Thus, we see that when the Hamiltonian constraintĈ =Ĉ grav +Ĉ matter acts on the wave function, since the operators either give functions of m or shift the basis state | m → | m + k δ , the equationĈψ = 0 gives a difference equation for the coefficients ψ m 1 ,m 2 ,m 3 . Up to now, we have considered δ to be a constant; however, when we consider quantization methods in the following section, we will see it is physically more viable to use a functional form δ = δ( m).

Quantization methods
We consider here only the gravitational piece C grav of the Hamiltonian constraint, coming from integrating the full constraint over the fiducial cell V, The Hamiltonian H grav depends on the curvature tensor F i ab , so to quantize this tensor, we write it in terms of our classical phase space variables, and then promote them to quantum operators. In the classical theory, this is done by considering a plaquette 2 ij -a closed loop whose sides are parallel to the coordinate axes x i and x j -and defining Here, Ar 2 is the area of the plaquette and the holonomy h 2 ij is taken as where δ i L i is the length of the ith edge of the plaquette, measured in the fiducial metric q ab . In this scheme, the area of the plaquette is shrunk to zero, so the particular choice of loop is irrelevant. However, this procedure cannot be done in LQG, due to the discreteness of the spectrum for the area operator Ar; instead, the plaquette can only be shrunk until its size reaches the minimum area eigenvalue ∆ 2 P of the area operator, where P is the Planck length 3 . There appears in the literature three separate methods for doing this, with various physical rationales; we discuss these here, and in particular their impact on the eventual difference equation.
• µ 0 scheme: In this case, the (matrix elements of the) holonomy h (µ 0 ) j is thought of as an eigenstate of the area operator Ar j = |p j |; the accompanying eigenvalue is set equal to the area gap ∆ 2 P . In other words, with for some constant C, we let This gives a constant value for µ 0 ; this scheme was used by Ashtekar, Pawlowski and Singh for the k = 0 isotropic case [5], and Chiou for Bianchi I [34]. The advantage of 3 From the full theory, we have that the dimensionless number ∆ = 4πγ √ 3 although the value ∆ = 2πγ √ 3 is sometimes used in earlier works; see the discussion in Appendix B of [1] for further detail. this method is that the difference equation has constant steps, i.e. defined on a lattice with even spacing, and therefore is relatively easy to find solutions. However, this scheme was realized to be problematic in terms of physical predictions [1]. For instance, as we will explain in Section 4, although the singularity is resolved -i.e. the cosmological model exhibits a "bounce" if a scalar field is chosen to be the time parameter -the critical energy density ρ crit at which this occurs is dependent on the momentum p φ of the scalar field, with ρ crit ∼ p −1 φ . Thus, a semi-classical state for the field, with large p φ , would have a correspondingly low value for ρ crit ; indeed, there is no lower bound for this density, so even densities considered far in the classical regime could be critical densities for the appropriate system.
To mitigate these issues, the constants δ i associated with each edge of the plaquette can be generalized into functions δ i ( m) of the parameters. Two possible ways to do this are considered below; we use here the notation of Chiou [35]. Both of these are dependent on the LQC spin network state | p -parameterized by the triad eigenvalues p i ( m) -and its relation to the plaquette chosen.
•μ scheme: This particular scheme was first suggested by Chiou [34], and it is inspired by the fact that physical areas in LQG are associated with edges between vertices in the spin network. Thus, we imagine the side of the fiducial cell V normal to the x 3 axis to be pierced by N 3 edges, each of which is associated with an area and giving a plaquette 2 12 .
Since for a state |p 1 , p 2 , p 3 on the kinematical Hilbert space, the area of the face of the fiducial cell, orthogonal to the direction x i , is given by p i , we have [9] Since we are assigning each of these edges in the x 3 direction to an area parallel to the 1-2 plane, this gives a plaquette with sides of coordinate size δ 3 ( p)L 1 and δ 3 ( p)L 2 , as in Fig. 1 [35]. The total area of the 1-2 side must be L 1 L 2 in the fiducial metric, so that N 3 δ 2 3 ( p)L 1 L 2 = L 1 L 2 . Solving for δ 3 ( p), and making the same argument for the 1-3 and 2-3 sides, we have that the functions δ i ( p) in the Bianchi I model are In Bojowald, Cartin and Khanna [21], this is the case where the number of vertices is proportional to the transverse area. This method is similar to that of the µ 0 scheme, in that the difference equation for the Bianchi I model is defined on a constant spacing lattice, after the appropriate redefinition of variables [34]. However, like the µ 0 scheme, working through the effective dynamics of the system shows that its physical predictions are dependent on the choice of fiducial cell V -in particular, the densities at which a bounce occurs in one of the three spatial directions [35].
•μ scheme: This method was briefly considered by Chiou [34], and developed in detail by Ashtekar and Wilson-Ewing [9]. The first step in this scheme is similar to that of theμ method, namely considering the edges passing through a particular face of the fiducial cell V. Thus, we again have the relation (2.4) for the relation between the number of edges passing through the 1-2 face and the corresponding value of the triad parameter p 3 . However, we now divide this face into plaquettes with lengths δ 1 ( p)L 1 and δ 2 ( p)L 2 in the fiducial metric as shown in Fig. 1, with the total area equaling L 1 L 2 , so that or If we iterate this for the other two planes, we find for Bianchi I that This is where the number of vertices in a particular direction is proportional to the extension in that direction [21]. It can be shown that the physical predictions of this version do not depend on the original fiducial cell 4 , so the critical density is invariant under changes in any scalar field used as a time parameter. However, it also results in a difference equation for the Hamiltonian constraint where the step sizes are parameter-dependent, so the lattice is not made of constant sized steps. From the analogous situation in computational physics, we denote this as a "lattice refinement" model.
Note that theμ andμ schemes are equivalent when there is only a single triad parameter, as in the isotropic cases. In the literature, the µ 0 method is the "original" quantization scheme, whileμ is called the "improved" quantization. Nevertheless, other schemes are possible in principle; see for instance [62], where isotropic embeddings of consistent quantizations of anisotropic models are considered.
computation. The situation is the opposite for LQC, where the difference equation is fundamental, and cannot be altered. However, there is still some ambiguity in the model, dependent on the quantization method used to arrive at the final equation; the suitability of these equations depends not only on mathematical arguments [40] (such as existence of an inner product in the Hilbert space), but also numerical factors, much as the case with differential equations. In this section, we consider several issues important in the process of numerically solving the Hamiltonian constraint for LQC. First, we must ensure that the difference equation is stable, meaning that solutions "close" in the appropriate sense do not exponentially separate as the system is iterated by increasing one of the model parameters. This leads to the notion of von Neumann stability analysis. Next, the quantum constraint equation must have an appropriate semi-classical limit. One must be careful that this limit is unique, which has not been the case for all models -the Schwarzschild interior will serve as our example here. Finally, it is helpful to consider changes in variables which decrease the difficulty in solving the difference equation. A particular problem, coming from the improved quantizations detailed earlier, is that the lattice spacings are no longer constant, but depend on the model parameters themselves. In some cases, this can be alleviated by choosing a new parameter as a function of the old, where the lattice spacing becomes constant up to higher order terms.

von Neumann stability analysis
In the latter half of the 20th century, techniques were developed to numerically compute solutions for important differential equations by discretizing the equation and evolving in time on a lattice where each point represented a specific space-time location. It was quickly realized that such solutions could grow without bound if a bad choice of discretization or time step was made. Thus was born the idea of stability analysis, to check the algorithm used and ensure solutions remain of finite size. Since the difference equation in LQC models is fixed for a particular quantization scheme, the role of stability analysis here is to provide input on the viability of the quantization chosen; the assumption is that an unphysical quantum Hamiltonian constraint would have several problems indicating its lack of consistency -not only on the basis of its predictions, but also from the characteristics of its numerical solutions. This notion has been born out for several models displaying many undesirable features, such as the isotropic subspace of Bianchi IX [22], the k = 1 Friedmann-Robertson-Walker (FRW) model [44], and both locally rotationally symmetric (LRS) Bianchi I and the Schwarzschild interior [64].
We summarize here the discussion of Bojowald and Date [22]. In their work, they ensure the notion of a stable difference equation by checking whether two solutions with "nearby" initial data will remain close to each other, rather than diverging. For simplicity, we start with a difference equation for a sequence ψ m of one parameter m, or The analogous analysis can be done for multiple parameters as well. We look at small regions in the parameter space, so that the coefficients A i (m) of the equation can be approximated by constant values A i . In general, solutions of difference equations with constant coefficients are found by writing down the characteristic equation and solving for the roots λ a [42]. If there are no repeated roots, the solutions are linear combinations of the form ψ m = β 1 λ m 1 + · · · + β 2k+1 λ m 2k+1 , with complex constants β i . For a single repeated root λ a of multiplicity d a , the basis is modified to {λ m 1 , λ m 2 , . . . , λ m a , mλ m a , . . . , m da−1 λ m a , . . . , λ m 2k+1 }, and similarly for more than one repeated root. As shown by Bojowald and Date, the physical requirement that the difference equation has the appropriate classical limit means that λ = 1 is a double root of the characteristic equation.
The need for a stable difference equation means that two sequences ψ m and ψ m , starting at a parameter value M with initial data {ψ M , ψ M +1 , . . . , ψ M +2k−1 } and {ψ M , ψ M +1 , . . . , ψ M +2k−1 }, will remain close together if the initial data is close as well. These sets of values and the difference equation (3.1) allow one to solve for ψ M +2k , and all subsequent values for m > M + 2k. If we think of each set of initial data as vectors S(M ) and S (M ), respectively, in C 2k and use the Euclidean norm || · ||, then the initial data are close to each other if || S − S || ≤ δ M . Then we require that the sequence δψ m ≡ ψ m −ψ m remains small in size as well; this places conditions on the roots of the characteristic equation (3.2). We see this as follows. For small ranges of the parameter δ M , δψ m satisfies the same constant coefficient difference equation (3.1) that the individual sequences ψ m and ψ m do. Thus, we can write δψ m as a linear combination of the basis elements as where a ranges over the separate characteristic roots, and C a,ra are constant coefficients. For a given sequence ψ m and set of initial data S(M ), we chose a nearby initial data set S (M ) such that only one of the values C a,ra is non-zero. Then as the sequence δψ m goes from parameter values M to M + ∆M , we have the norm of the difference between the two sequences goes as In order to keep δψ M +∆M within the same bound on the norm as δψ M , we must have all the roots |λ a | ≤ 1.
In Appendix A, we examine the stability of two versions of the k = 1 FRW model, that of Green and Unruh [44], the other by Ashtekar, Pawlowski, Singh and Vandersloot [7]. Although the Green and Unruh model is unsatisfactory due to its unphysical predictions, the stability analysis shows that the problem lies with the quantization method and the resulting Hamiltonian constraint, an issue repaired in the work of Ashtekar et al.

Semi-classical limits of the dif ference equation
When using these discrete equations to obtain predictions for semi-classical effects, one must be careful about how the limit is taken. We assume that the model difference equation is essentially continuous in the semi-classical limit, in order to match with the existing theory of general relativity. In other words, we should obtain the WdW equation as the limiting case of the difference equation. For LQC, the (dimensionless) quantum parameter m is related to the continuous variable p(m) by where γ is again the Immirzi parameter, and is proportional to the square of the Planck length. Since p(m) is the eigenvalue of the triad operator, it is related to the scale factor functions appearing in the classical metric. One possible way to take the semi-classical limit is letting γ → 0; however, since γ is a small but fixed constant, we must assume that γ → 0, m → ∞ such that p(m) remains constant. This is what is called the "pre-classical limit" in Bojowald and Date [22]. This respects the physical input that γ represents. Then, we can approximate the discrete evolution equation (3.1) with a differential equation for a smooth wave function ψ(p(m)). This can be obtained by writing the wave function ψ m and the coefficients A m i (γ) in a power series in γ. Putting these expansions back into the original discrete relation, we find a continuous differential equation as the leading order contribution. For LQC, the WdW equation should emerge as the O(γ 0 ) term, with the higher order γ terms corresponding to quantum corrections, and all terms with inverse powers of γ should cancel out. The conditions to obtain such a physically reasonable situation are derived by Bojowald and Date [22].
In some sense, taking the limit as described above is an example of finding the semi-classical equation first, then solving the resulting WdW equation for the wave function. However, we could solve for the wave function first without appealing to the physical nature of the parameters -using generating function techniques [30,31,32,33] or other methods -and then pick out those solutions that have the appropriate behavior. This corresponds to taking the asymptotic limit m → ∞ without reference to the parameter γ. This is the "pre-classical limit" as it has been used by Cartin and Khanna [31,32,33]. The leading order term in a m power series will remain a difference equation. Generally, solutions of difference equations can have many undesirable properties, such as oscillation in sign with each succeeding step (i.e. ψ m (−λ) m for a positive constant λ in the large m limit) or unstable solutions that increase without bound. Some control over these solutions can be obtained by choosing the proper initial data. However, it is possible that there is not enough free parameters to avoid unphysical behavior.
As an example, we look at the original quantization of the Schwarzschild interior [2,58]. If we write it in terms of the two triad components p b , p c (because of spherical symmetry), the Hamiltonian constraint is of the form For most of the terms and coefficients of the difference equation, taking the p b , p c → ∞ limit is equivalent to that of γ 2 P δ → 0. However, there is one problematic term in the equation, namely, Specifically, if we look at the quantum gravity limit γδ → 0, but allowing p b , p c to remain the same size by requiring m 1 , m 2 → ∞ accordingly, then while in the large momentum limit, then we have The difference in the limits means this model has the undesirable feature that the nature of the semi-classical limit changes, depending on how this limit is taken. Thus, the best case is to find a Hamiltonian constraint and difference equation where the semi-classical regime is the same, regardless of which limit is used to reach it.

Lattice ref inement and model reparametrization
The advent of these improved quantization schemes unfortunately leads to difficulties in solving the difference equation, which now involves lattice refinement. As we saw above, theμ andμ schemes both involve lattice steps that are functions of the model parameters m. If we have a difference equation of two parameters, say m 1 and m 2 , then generic terms in the Hamiltonian constraint for the wave function ψ(m 1 , m 2 ) will be of the form where the step functions δ i (m 1 , m 2 ) are non-linear functions of the parameters. These difference equations are more difficult to solve, since unlike relations with fixed steps, the equation for a given choice { m} depends on values of ψ( m) not calculated from any other equations. An example of this is the difference equation arising from a lattice-refined Schwarzschild model [21], where Note that δ 0 is a fixed parameter related to the minimum area eigenvalue (i.e. δ 0 ∼ √ ∆). In this case, the Hamiltonian constraint is of the form To solve these types of lattice refined equations, we consider a change in the parameters used in the model, where the equation in the new variables now has constant size steps [21,59]. This is possible always for one-parameter models, i.e. isotropic cosmologies, but may not be feasible for equations with multiple parameters. This can be seen as follows. If the wave function depends on only a single variable m, then with non-constant steps δ(m), we have terms in the difference equation of the form ψ[m ± δ(m)]. We assume the step function δ(m) is linear in a parameter δ 0 which we use as our lattice spacing. To get constant step sizes, we must write the wave function ψ(m) in terms of a new function N (m), such that This would imply that so the values N (m) form a lattice with a constant spacing δ 0 . Taking the Taylor series expansion of the left-hand side of (3.4) gives that this condition requires If the function N (m) is such that higher order terms in (3.4) correspond to higher order in as well, then they can be dispensed with using the appropriate factor ordering of the Hamiltonian constraint operator [21]. Nelson and Sakellariadou showed the issue of factor ordering is intimately related to that of lattice refinement by showing that, for the flat isotropic model, there is a unique choice of quantization scheme (the "improved" version) where the factor ordering ambiguities disappear at the semi-classical level [60].
An example of this reparametrization procedure is seen for the flat isotropic model [6]. The authors of that work use a geometric argument -considering integral curves of operators related to the holonomy operator -to find the new function N (m) such that ψ(N (m)) has a constant size lattice. Here we show the same result is obtained by the method described above. In particular, in this isotropic model, the need to obtain valid physical predictions (much like for the Bianchi I model) makes it natural to have a step function δ(m) = 2/(3K √ m), where K is a purely numerical constant. Solving the equation (3.5) for N (m) with δ 0 = 1 gives N (m) is proportional to the volume eigenvalue associated with m, that is, v(m) = K sgn(m)|m| 3/2 . Therefore, if we redefine the wave function asψ(v) = ψ(m), then the action of the operator with the improved quantization functionμ(m) gives equidistance steps.
This method becomes problematic for multiple parameter models. As we have seen with the functions δ i ( m) -for Bianchi I, given for theμ scheme in (2.5), while in the Schwarzschild interior, shown in (3.3) -the step functions δ i can depend on all parameters m i . Suppose for specificity we have two parameters. Then, we seek to find two independent functions N (m 1 , m 2 ) such that for some constants k i , and k a function of k i . This requires that For this is to be true for all ±k 1 , ±k 2 , we need specifically for the Schwarzschild case that 2 for constants C i ; these equations only have one independent solution, namely N (m 1 , m 2 ) ∼ m 1 |m 2 |, i.e. the eigenvalue of the volume operator for the model. Therefore we see that it is possible to have constant step sizes in one of the two parameters, but not in both [21]. This simplification, however, is useful in putting the difference equation into a form more tractable for numerical computation [46] 5 . Beyond this analytical technique to place the original difference equation into a more tractable form, there are also various methods of numerically solving the equation along the lines of the Taylor expansions outlined above. One possibility is to interpolate the solution between lattice points already calculated, using a least-squares fit [65] (more generally, the wave function is expanded in a Taylor series of the parameters around these previously computed values [59]) and then use this interpolant to compute the values needed to step the evolution forward.

Examples of LQC models
In this section we will mainly restrict our attention to those studies that have faced the purely quantum difference equation via numerical techniques. This has been performed for a number of homogeneous models: FRW with spatially flat k = 0, closed k = 1 and open k = −1 topologies, with and without the presence of a cosmological constant and, less systematically, for anisotropic Bianchi I models and Schwarzschild interior. The quantization of Bianchi II [10] and Bianchi IX [75] models have also satisfactorily been carried out in LQC, but the numerical analysis is still missing. In a last subsection we will mention several studies where the effective analysis has been done, since we consider this is a very important research area of LQC. On the one hand, many of the main features of the LQC evolution has been first derived via this effective treatment. On the other hand, this effective analysis have been applied to a wide variety of scenarios, such as inflation as well as to some anisotropic and inhomogeneous model, which permitted to extract physical information through systematic numerical evolutions that would be difficult to perform in the full quantum setting. Of course, since to a certain extent, heuristics is used within these approaches, care is needed when interpreting the validity of the physical results.

Flat Friedmann-Robertson-Walker
Because of its simplicity and the fact that it describes to a very good approximation the evolution of our universe, the FRW solution has been the most studied model in LQC literature. In fact, the first positive results regarding the singularity resolution in the context of LQC were due to Bojowald [14,16] for this model. The singularity is indeed resolved not simply because the evolution equation is discrete and, hence, the state jumps over the zero point. In fact, the evolution through the classical singular point is also well defined. Therefore, any initial data can be in principle deterministically evolved through the singularity. But, in order to know precisely what is "on the other side of the singularity" and what is the precise behavior of the physical system, one needs to construct Dirac observables and solve the difference evolution equation. Here is where one has to resort to numerical methods. The seminal work in this aspect was performed by Ashtekar, Pawlowski, and Singh (APS) [5]. Let us comment more precisely about their approach. In their work, the spatially flat FRW model was analyzed with a matter source of massless scalar field φ. This model is particularly interesting for the study of singularity resolution since classically a (big-bang) singularity is unavoidable. The scalar field is a monotonic function of the cosmological time, so it can itself be interpreted as an internal clock instead of the scale factor a. On the one hand, this is a widely used approach in cosmology to describe the evolution of the system without making reference to any fiducial structure or coordinates. On the other hand, in the closed (k = 1) model, although the universe undergoes a recollapse, the scalar field φ, contrary to the scale factor a, is still a monotonic function and can be interpreted as a time variable.
In this model, the phase space is coordinatized by the gravitational degree of freedom (c, p), with Poisson bracket {c, p} = 8πG/3, by the scalar field φ and by its conjugate momentum p φ with bracket {φ, p φ } = 1. The only constraint left is the Hamiltonian that, in terms of these variables, takes the the form, C grav := − 6 γ 2 c 2 |p|, and C matter := 8πG p 2 φ |p| 3/2 .
For this system, a complete set of classical Dirac observables is provided by the constant of motion p φ and by p| φ 0 . The latter determines the value of p at the instant of time at which φ = φ 0 . This set is complete because they univocally define a trajectory on the phase space. Therefore, in order to extract physical information of the quantum system, the evolution of these observables must be analyzed. The quantization of the geometric degrees of freedom is performed by following the methods of full LQG and therefore, as has been explained in Section 2.1, making use of holonomies N µ and triads p as fundamental variables. In this way, the gravitational kinematical Hilbert space is H grav kin = L 2 (R B , dµ B ), being R B the Bohr compactification of the real line with the Bohr measure dµ B . This Hilbert space is properly spanned by the eigenfunctions |µ of the triad operator:p|µ = 8πγ 2 Pl µ/6|µ . On the other hand, the matter degree of freedom φ is quantized via a standard Schrödinger representation on a Hilbert space H matt kin := L 2 (R, dφ). Here, the operator associated to the conjugate momentum of the matter fieldp φ acts as a derivative In this first study by Ashtekar, Pawlowski and Singh, the lapse was taken equal to unity and the µ 0 scheme was used. The physical states Ψ(µ, φ) were then defined as those that are annihilated by the Hamiltonian constraint. On the full kinematical Hilbert space, given by the product between the geometric and matter spaces L 2 (R B , B(µ)dµ B ) ⊗ L 2 (R, dφ), the constraint can be written as an evolution equation in terms of the internal emergent time φ, where the difference operator is given by, for certain functions f (µ) and B(µ). In particular, B(µ) is proportional to the eigenvalues of the operator 1/|p| 3/2 , which is the only geometric factor that appears on the matter part of the Hamiltonian. Due to the difference operator, superselection sectors appear. Let us denote by L ε the lattice of points {ε + 4µ 0 } on the µ axis. Therefore, the subspaces H ε ∈ H kin , defined as those states with support only on L ε , are left invariant under the action of the operator Θ 0 and the Dirac observables. The positive-frequency solutions to this equation form the physical Hilbert space and they can be written in an integral form as, In this expression,Ψ(k) is an arbitrary function, ω(k) 2 = πG(16k 2 + 1)/3, whereas e (s) k is a symmetric linear combination of eigenfunctions of the operator Θ 0 with eigenvalue k. These combinations have support on the lattice L ε and are, by construction, invariant under the triad orientation reversal operation. The physical inner product for a given value of φ = φ 0 is, and the action of the Dirac observables is given as, In this way, there are two possible methods in order to solve the constraint equation. On the one hand, the constraint can be understood as an evolution equation in φ. Hence, choosing initial data Ψ(µ, φ 0 ) at a given value of the scalar field φ = φ 0 , the physical state is obtained by evolving them via equation (4.1). On the other hand, one could obtain the eigenfunctions of the operator Θ 0 , construct the appropriate symmetric linear combinations e (s) k , pick up a functionΨ(k), and perform the integral on the right-hand side of equation (4.3). The advantage of the first method is that it is not necessary to deal with the eigenfunctions of Θ 0 but, from a numerical point of view, it is more difficult since a large number of differential equations must be solved. The choice of the free data in each procedure (either as initial data or as the free functionΨ) will introduce the physical input on the system. Obviously, the most interesting set of initial data is given by those that describe a large classical universe at late times, since this is the current state we observe in our universe. In [5] these both procedures were followed to numerically solve the behavior of the system. Let us comment on those with some more detail.
For the evolution procedure, three different methods were used to pick up initial conditions. In the first one a Gaussian (with respect to the Wheeler-DeWitt inner product) peaked on a classical universe was chosen. In the other two, an analytical solution of the Wheeler-DeWitt equation, that is semi-classical at late times, was used. Since equation (4.1) represents an infinite number of coupled differential equation, one needs to restrict to a finite domain in order to perform the numerical evolution. The boundary was chosen far away from the peak of the initial wave packet, and purely outgoing boundary data were chosen there. In addition, for the discretization of the time derivatives of (4.1), the fourth-order adaptive Runge-Kutta method was applied. Regarding the direct evaluation of the integral (4.3), the functionΨ(k) was chosen as a Gaussian (with the usual measure) peaked on certain value k ( − 1) so that it describes a semi-classical state at a time φ 0 . The eigenfunctions of Θ 0 were numerically obtained and the integral (4.3) was evaluated using fast Fourier transform.
Even if, obviously, the considered different initial data sets led to different quantitative results, the generic behavior and the main qualitative results were robust independently of the method of choosing the initial data or the value of the free parameters within those methods. Evolving backwards, the expectation value of the observables follow the trajectory of their classical counterparts until a certain minimum value. At this point they depart completely from the classical behavior and they bounce, keeping always a finite value, and connecting an expanding with a collapsing branch of the universe. In addition, the state is sharply peaked during the whole evolution. In this way, the classical singularity is resolved and the LQC evolution equation connects two classical branches of the universe through a quantum phase. The bounce happens because the system reaches a point where quantum gravity effects are appreciable and make the gravity repulsive. In principle, one would assume that this must happen for high values of the matter density, since we do not observe any quantum gravity in our daily experience. Indeed, this fact is one of the drawbacks of the analysis presented in [5] since the matter density at the bounce happens to be inversely proportional to the expectation value p φ , as discussed in Section 2.2. Therefore, one could pick up a value of p φ so that the bounce would occur at very low values of the matter density (and, hence, curvatures). This drawback was fixed in [6] by introducing the so-called improved dynamics.
In this analysis, theμ quantization scheme was used. Note that in this homogeneous scenario theμ andμ schemes presented in Section 2.2 coincide. The mathematical structure of the formalism is essentially the same as for the µ 0 case explained above, but physically the main flaw of that approach was cured. Within the improved-dynamics scheme, for all the considered cases, it was obtained that in a backward evolution the expectation values of the Dirac observables remain peaked on their classical trajectories (corresponding to an expanding branch of the universe) until the matter density reaches a critical value of of ρ crit = 0.82ρ P l , where the system undergoes a bounce. After that, the energy density starts decreasing and, once its ratio with ρ crit is small, the system evolves again peaked on a (contracting) classical trajectory. As already mentioned, the value ρ crit is the same for all analyzed simulations. In fact, from the effective analysis of this system, it is possible to write down a modified Friedmann equation which provides an analytical expression for this object: ρ crit = 3/(16π 2 γ 3 G 2 ). This value is in very good agreement with the one found in the full quantum simulations.
Regarding the mathematical issues of these second analysis [6], in order to label the eigenfunctions of the triad operator, the real parameter µ is replaced by the affine one v := Jsgn(µ)|µ| 3/2 , with J = 2 √ 2/(3 3 √ 3). In this representation, the action of the components of the holonomies, e iμ/2 , simply shift in a constant step the states, e iμ/2 |v = |v + 1 .
This basis is best suited to the volume operator (associated to the elementary cell): The quantized Hamiltonian constraint, written as an evolution equation, takes the same form as in equation (4.1), with the same structure for the operator, Here a brief comment about the work by Laguna [49] is in order, which has been a source of some confusion (for further discussion, see Section VIB of [7]). The paper by Laguna first deals with the evolution of wave packets under the LQC Hamiltonian constraint for this isotropic model, and then compares it to the WdW equation for that system. Finding a dispersion equation for plane wave solutions of these models, it is shown that the group velocity goes to zero in the LQC model -and thus the quantum bounce occurs as the contraction of the spatial volume becomes an expansion -precisely at the critical maximum energy density ρ crit shown previously. However, Laguna then goes on to consider modifications of the original LQC equation, which has the form of the shallow water equation, that would eliminate the bounce. In this context, of interest to those working in the field of computational physics, he essentially adds extra terms to the equation, none of which are motivated from the original LQC quantization scheme. In this sense, the elimination of the bounce as a spurious effect is not applicable to LQC directly, since it involves going outside the quantization method used to arrive at the Hamiltonian constraint operator.
Finally, even if this is the simplest model of LQC, there is still no complete consensus on the precise details of the quantization. Even inside the framework of the improved dynamics (choosing theμ scheme) there is freedom in the construction of the quantum Hamiltonian since, as is usual in any quantization process, in order to represent the classical constraint it is possible to choose different factor orderings for the (non-commuting) operators. In addition, the fact that one could choose arbitrarily the lapse function also introduces another ambiguity and gives rise to different quantization schemes and, hence, to different discrete operators Θ. In a recent paper [56] four of these prescriptions have been analyzed and compared. The first prescription is the one due to APS and has already been presented above. In this case, the lapse is taken equal to unity and the quantum constraint is symmetrized with respect to holonomy operators. In [3], a simpler version of this scheme was presented in the sense that it allows to carry out the complete analytical resolution of the system. This is usually known as solvable Loop Quantum Cosmology (sLQC). Its key feature is that the lapse function is chosen to be V /(8πG), what produces a quantum Hamiltonian constraint free of inverse of triads (volume) operator. The origin of the third analyzed prescription lies on the analysis of Bianchi I models [9]. As in the APS case, the lapse is also taken equal to one but, contrary to the previous prescriptions, the operator sgn(v) appears in the quantum constraint. This prescription is usually denoted as MMO, due to the names of the authors of the article [51], where it was first introduced. Finally, in the same paper [56] a new prescription (sMMO) is introduced that combines different aspects of the sLQC and MMO ones.
The system has been analyzed for all the mentioned prescriptions following the route of the construction of eigenfunctions of the corresponding discrete operator Θ. From a numerical point of view, it was found that the MMO and sMMO prescriptions are much more efficient. This is due to the fact that, for these cases, the spectrum of Θ is non-degenerate: the state is required to be symmetric with respect to a change µ → −µ and, hence, the mentioned eigenfunctions are determined once their value at a given lattice point is fixed. In addition, nontrivial differences between prescriptions have been found, although they might be non observable in practice. For instance, three observables (the volume, the Hubble parameter, and the energy density of the scalar field) have been analyzed and, in all the cases, the difference between their expectation values for different schemes turn out to be smaller than their corresponding dispersion.

Closed Friedmann-Robertson-Walker
Quantum gravity effects are supposed to modify the behavior of the system only at Planck scales in such a way that the classical singularity is resolved. In addition, they are expected to be suppressed at low curvatures in order to recover the classical trajectory. Since classically the universe undergoes a recollapse, the closed FRW model represents a more restrictive testbed model for LQC methods than the flat case: this recollapse should not be affected by quantum geometry modifications. One of the first proposals for the quantization of this model was presented in [16] in combination with the flat (k = 0) model. The problem, as has been explained above, is that this quantization was proven to be unstable [22] for the closed case. Nevertheless, since this isotropic model can be understood as a particular case of Bianchi IX [28], Green and Unruh [44] made use of techniques developed for the BIanchi IX model [18] in order to construct a quantum hamiltonian constraint within the µ 0 scheme. In this analysis, the numerical resolution of this constraint was also obtained. Even if an inner product and observables, that would provide a neat physical interpretation of the system, were missing, the divergent behavior of the solutions at large scales already indicated that the classical recollapse could not recovered.
In [7], theμ scheme was used to quantize the closed model. In this article some important developments were made that shed light on the concerns explained above. Physically, the most important result was that, due to the quantum effects, the classical (big-bang and big-crunch) singularities are both resolved and replaced by a quantum bounce. This led to a cyclic universe picture, where the quantum wave function evolves through infinitely many classical cycles. The value of the matter density at the bounce is in agreement with the critical value ρ crit , introduced above for the flat model, (provided that the volume of the universe reaches a macroscopic size). More importantly, the state remains sharply peaked on the classical trajectory (within the corresponding cycle) and, thus, undergoes the recollapse at a maximum volume that is in complete concordance with the predictions of the classical theory. Therefore, the infrared limit of the theory turns out to be right. From a technical point of view, the analysis of [7] is formulated in terms of connections. Contrary to previous quantizations which, for technical reasons, constructed holonomies in terms of the extrinsic curvature, in this approach the connection was used. This procedure mimics better the approach followed in full LQG.
The overall structure and analysis is completely analogous to the flat case [6]. The kinematical Hilbert spaces are again given by H kin . The Hamiltonian constraint can be written as, with a discrete operatorΘ 1 Ψ(v, φ) :=Θ 0 Ψ(v, φ) + Ω(v,μ, r 0 )Ψ(v, φ), so that the scalar field φ serves as emergent time. The function Ω depends on the radius of the 3-sphere r 0 with respect to the fiducial metric. The flat k = 0 case is recovered by letting r 0 tend to infinity, which makes the function Ω vanishing. This discrete operator also produces superselection sectors.
The key difference with theΘ 0 operator is that the eigenvalues ofΘ 1 are discrete (see [70] for an analytical proof). More precisely, eigenfunctions exist for all values of ω but, in order to be normalizable, they must decay exponentially in both limits v → ±∞. And this only happens for certain discrete values of ω. Numerically these eigenvalues were found following a bisection method (taking advantage of the fact that the normalizable eigenfunction is always a critical solution between functions exponentially diverging to ±∞ as v → ∞).

Open Friedmann-Robertson-Walker
The study of the open k = −1 isotropic case has been carried out in [74] within theμ scheme. The discrete evolution operator Θ −1 is constructed and the scalar field is understood as emergent time. Following the procedure introduced in [5,6], eigenfunctions of the Θ −1 operator were obtained and then Fourier transformed, in the analogous equations to (4.3), in order to get the physical states. The qualitative behavior of the system is completely similar to the plane k = 0 and closed k = 1 models described above. Regarding the energy density of the scalar field at the bounce, an analytical expression is obtained by effective means, which coincides with the results of the full quantum system. This expression depends on the momentum of the scalar field p φ , but it tends to ρ crit for large values of p φ . In [69] some issues, like the choice of the loop to define the holonomies, were clarified, but this analysis is not yet totally satisfactory. First, the defined discrete operator is not essentially self-adjoint. Also, as in the first analysis of the closed model, the extrinsic curvature, instead of the connection, is used to construct the holonomies. This route is taken in order to avoid the technical difficulty of constructing an appropriate loop. In principle this last drawback could be cured following the proposal presented in [10] for the Bianchi II model, but a complete analysis is still missing.

Cosmological constant
The inclusion of the cosmological constant has only been considered for FRW spatially flat models with a matter content of a massless scalar field that, as in previous models, also serves as internal time. In particular, in Appendix A of [6], the main techniques for theμ scheme were already put forward and some numerical evolutions of this system were presented. In all cases the the classical singularity happened to be replaced by a quantum bounce and, contrary to previous treatments [11,63], the correct semi-classical behavior of the system was recovered. Moreover, the total energy density (defined as the sum between the energy density corresponding to the matter and to the cosmological constant: ρ := ρ φ + Λ/(8πG)) at the bounce was found to be the same as in the Λ = 0 cases and, hence, independent of the value of the cosmological constant.
Already at a classical level, the behavior of the system is completely different depending on the sign of the cosmological constant. Classically, the universe with Λ < 0 begins at a big bang singularity at a "time" φ = −∞, expands until the total energy density reaches a minimum (in fact it vanishes) where, even if we are in the spatially flat case, undergoes a recollapse and finishes at a big-crunch at φ = ∞. In [12] a concrete and detail analysis of the quantization of this model was carried out and the results were similar to the closed k = 1 FRW with a vanishing Λ. The spectrum of the corresponding discrete operator turns out to be discrete and, therefore, its eigenvalues are isolated. Hence, a similar bisection method to the k = 1 case was applied. The quantum system also resolves both initial and final singularities, replacing them by quantum bounces happening at a value ρ = 0.82ρ Pl , which results in a cyclic universe model. In addition, the classical recollapse was also recovered in complete agreement with the classical results.
In the classical Λ > 0 model, after the big-bang, the scale factor of the universe reaches infinity at a finite value of φ = φ 0 . Therefore, contrary to the Λ < 0 model, in this case the values of the scalar field are constrained to a finite interval. In [47] Kaminski and Pawlowski constructed the discrete evolution operator making use of the three quantization prescription commented above (APS, sLQC, and MMO) and studied its possible self-adjoint extensions. The properties of this operator happened to be completely different depending on the value of the cosmological constant. For 0 < Λ < Λ crit , where Λ crit = 8πGρ crit is the value for which the energy density of the cosmological constant equals the maximum total density ρ crit , the operator admits many extensions; whereas for Λ ≥ Λ crit , the operator turns out to be essentially selfadjoint. Observations favor a positive but small value of the cosmological constant Λ Λ crit and the different self-adjoint extensions could give rise to non-equivalent unitary evolutions. Even so, a recent numerical analysis [4] shows that for the relevant states (which are initially sharply peaked at a classical point of the phase space) the evolution is quite insensitive to the chosen extension. Furthermore, the classical initial singularity is resolved. The self-adjoint extensions permit to follow unitarily the evolution in φ after the classical divergence at φ 0 and give rise to a recollapse. Thus, a cyclic universe picture is again obtained.
In order to finalize this section, we would like to briefly comment on another approach that has been proposed in order to include the cosmological constant in LQC through the so-called unimodular gravity [39]. At classical level, unimodular gravity is equivalent to general relativity with the only difference that allows for a dynamical cosmological constant. In fact, the momentum conjugate to the cosmological constant can serve as emergent time and it seems that the quantum theory one obtains is simpler than in the cases with the scalar field clock. Remarkably, the main results of singularity resolution and the critical value of the energy density coincide for both cases.

Bianchi I
Bianchi I model is the simplest non-isotropic model one could consider. Even though, most of the analysis for this model, as a test case for anisotropic cosmologies, has been performed at the level of the effective dynamics rather than solutions of the Hamiltonian constraint difference equation. In part, this is due to the greater complexity of the constraint, with a lattice having non-constant steps in theμ quantization scheme. Complete quantizations have been performed both for a massless scalar field matter content [34,9] and for vacuum [52,54]. In this section, we will just mention these two quantizations and the only numerical analysis [53], even if not completely systematic, that has been carried out to analyze the full quantum evolution of this model in theμ scheme.
As in the previous cases, the scalar field serves as an internal clock and can be treated as an emergent time. Making use of this fact, the quantization of this system was first carried out in [34] for both µ 0 andμ schemes. Nonetheless, this approach suffered from several drawbacks, like the incorrect infrared limit and dependence of the physical results on the choice of cell [71], that were fixed by considering theμ quantization [9]. The particular form of the Hamiltonian constraint given in [9] is non-symmetric in its three parameters due to a variable transformation chosen by the authors; instead of using all three triad parameters p i to serve as variables in the difference equation, one of the three is transformed (as discussed earlier in Section 3.3) to a volume v. Thus, the model is given in terms of two directional triad parameters and the total volume. In addition, it was shown that the difference equation defined by that Hamiltonian is unstable to growing modes [61], although it is unclear whether these modes are eliminated by the inner product. An open area for further work is whether this construction is the best choice, or if other schemes (including the symmetric one) are more viable in a numerical sense.
On the other hand, in [52] a complete quantization of the vacuum Bianchi I model was developed within theμ scheme. Subsequently, the analysis was also extended to theμ quantization prescription [54]. In this approach, in order to have a concept of evolution, the role of the clock must be fulfilled by some geometric degree of freedom. In the context of the quantization given in [52], two different variables (a triad coefficients and its conjugate momentum) have been tested numerically [53]. This numerical analysis has been restricted to Gaussian states and, due to complicate numerical problems, long term evolutions could not be checked. Therefore, further numerical analysis are needed to conclude, for instance, that all semi-classical states remain so after the bounce. Even though, the bounce picture was proven to be robust under the new interpretation of time.

Ef fective models
The approach given by the effective equations is based on the so-called geometric quantum mechanics. In this context, the quantum Hilbert space is considered to be an infinite dimensional phase space with the symplectic form given by the imaginary part of the inner product. Taking expectation values one can relate both spaces: each operator on the Hilbert space is mapped to a function on the phase space. Then, it can be seen that the commutator operation of the Hilbert space is directly related to the Poisson brackets of the phase space. Since the equations of motion need to be solved only for the classical degrees of freedom (that are interpreted as the expectation values of the observables in the quantum theory) with given quantum geometric corrections, one deals with differential equations and not with discrete equations. In particular, for the case of homogeneous systems, we only need to deal with ordinary differential equations since the unique dependence of all objects will be the time variable. This fact makes easier, specially from a computational point of view, to analyze the different models effectively rather than facing the full quantum evolution equations.
As far as we know, there is only one analysis [72] in the literature, corresponding to the quantization of FRW model, that has systematically derived the effective Hamiltonian and equations of motion by computing the expectation value of the quantum Hamiltonian operator for Gaussian coherent states and then performing an asymptotic expansion. For the rest of the cases, the effective equation program has essentially consisted on replacing the operators in the quantum Hamiltonian constraint by their expectation values. This provides a first order approach inasmuch as one disregards all state dependent information. But it gives rise to a modified classical Hamiltonian that contains corrections by LQG effects, from which the effective equations of motion for different variables can be obtained. The effective evolution of these variables is expected to mimic the behavior of the quantum system. In fact, in all the cases presented in this section, that the quantum discrete equation has been numerically solved, the pure quantum results have been compared with the effective ones, and a very precise agreement has been obtained.
Effective descriptions of quantum FRW are used to intuitively understand different aspects of the system and generalize the obtained results through a purely quantum analysis. For instance, the analytical form of ρ crit has been obtained [5]. In addition, the presence of a short phase of superinflation immediately after the bounce was first shown for a scalar matter field in [15]. This result was afterwards generalized to other matter models [67]. During this superinflationary phase both the Hubble parameter and its time derivative are positive; thus the acceleration of the universe is faster than in a de-Sitter spacetime. Furthermore, the resolution of strong singularities has also been shown throughout this effective analysis [68]. A particularly important question is whether LQC may leave any testable footprints in the CMB. This issue has been analyzed also by effective means in several articles (see, for instance, [73,45,66,43,57,20]) with a wide variety of assumptions, but yet no final consensus has been achieved.
On the other hand, effective equations have been applied to the analysis of more complicated models, for which a complete quantum treatment is still out of reach. Regarding the anisotropic models, in [23] it has been shown that the chaotic behavior of the classical Bianchi IX model near the singularity is reduced by quantum effects for a µ 0 scheme. Effective equations for this model in theμ scheme have already been proposed [75], but a numerical analysis is still missing.
The effective dynamics of the vacuum Bianchi I model was first analyzed in [41] in the context of µ 0 scheme. This case is simple enough so that the analytical solution can be achieved, which clearly shows the resolution of the classical singularity since the volume is not allowed to vanish. For generic matter content, and in the context of theμ scheme, the effective study was carried out in [36]. In particular, the vacuum, scalar field and radiation matter contents were analyzed. For the particular case of a scalar matter field, in [35] a comparison between theμ andμ schemes was presented. This analysis showed that in both cases the classical singularities are resolved and replaced by big bounces (one for each spatial direction) but certain key details differ. In theμ scheme these bounces occur when the corresponding directional densities ρ I := p 2 φ /p 3 I , which depend on the fiducial cell, reach certain critical value. Whereas, for theμ case, the characterization of the bounce is given by the critical value of the matter density. Thus, in the last case, all three bounces take place almost simultaneously. The classical collapsing universe that is obtained "on the other side of the bounces" also depends on the quantization scheme.
There are also several effective analysis of the interior of the Schwarzschild black hole and even if, as far as we know, none of them is based on a rigorous quantum theory, appealing results have been obtained. For instance, in [13] the singularity resolution was already obtained for both µ 0 andμ quantization schemes. Constructing on these results, in a manner parallel to the above-commented Bianchi I case [35], Chiou performed a systematic comparison between the physical consequences ofμ andμ schemes for this model [37]. In both cases the classical singularity is resolved but the final picture is quite different. Taking into account an extended Schwarzschild spacetime, in theμ scheme quantum bounces replace the classical white and black hole singularities and work as a bridge between their corresponding regions. On the other hand, in theμ case, the classical singularity is resolved, which implies that the even horizon is diffused. This produces a baby black hole with a very reduced mass in the consecutive classical cycle. A process that continues until the spacetime enters a highly quantum regime.
Regarding cosmological inhomogeneous models, the vacuum Gowdy spacetime (with a threetorus spatial topology) has been numerically analyzed in [29] for theμ scheme. This effective treatment is based on the hybrid quantization of the mentioned model [50], where different degrees of freedom are treated on a distinct footing. The homogeneous degrees of freedom are quantized polymerically, whereas a regular Fock quantization is developed for the inhomogeneities. The main result of this analysis is that the bounce picture is robust in the presence of inhomogeneities. In addition, a Monte Carlo analysis showed that the amplitude of the inhomogeneities is statistically conserved through the bounce for highly inhomogeneous universes, whereas it is amplified for almost homogeneous cases.
In order to finalize, we would like to briefly mention that a different effective method (the so-called truncation method [8]) has been used in quantum cosmology [26,25]. Following the geometric quantization mentioned above, a state can be completely parameterized by the expectation values of the basic operators, likeq andp in quantum mechanics, and an infinite number of moments (expectation values ofq apb for all a and b). The equations of motion for all these variables, which constitute an infinite system of coupled differential equations, is then completely equivalent to the quantum dynamics of the corresponding wave function. Nonetheless, in certain (semiclassical) situations, this infinite system can be truncated by keeping terms up to a given order (defined as the sum a + b). In this way, it is possible to analyze the back-reaction of the high-order moments on the trajectories of the expectation values. For instance, this approach has been used in [24] to obtain and analyze the effective equations (for both WdW and loop quantizations) of cosmological homogeneous models with a massive or interacting scalar field at second order. In contrast to the free scalar field case, the back reaction is very relevant in this model since the classical equations for the expectation values receive corrections due to the back-reaction of fluctuations. The LRS Bianchi I model was also analyzed in [27]. In this case, the matter content was chosen as an isotropic scalar field with a negative energy density given in terms of the anisotropic scale factors. This exotic matter field makes the system treatable in terms of physical coherent states. Therefore, a comparison between the effective and the purely quantum system was performed and both approaches showed an excellent agreement. Finally, in [19], in order to systematize the computation of equations of motion for the moments at very high order, powerful algebraic computational tools have been constructed. Furthermore, these tools have then been applied to the study the quantum back reaction on a homogeneous universe with positive cosmological constant.

Summary
In this article we reviewed a variety of numerical techniques and their application to several relevant LQC models. We commented on the features of the quantum Hamiltonian constraint in the context these models that make these techniques readily applicable and extremely useful, especially for exploring different quantization schemes and examining the semi-classical limit. More specifically, inspired by its role in numerical analysis, we reviewed the concept and application of von Neumann stability to several LQC models. We also reviewed the significant role played by lattice refinement in the context of the semi-classical behavior of these models. This article also includes a survey of a wide variety of isotropic and anisotropic LQC models and how these numerical techniques have contributed towards a deeper understanding of their features, and a discussion of open issues that require further work to advance this emerging field.

A Stability analysis of the closed FRW model
As an example of stability analysis, we consider here the case of the k = 1 FRW model, with two incarnations of the quantization method -the work of Green and Unruh [44] with the old quantization scheme, and Ashtekar, Pawlowski, Singh and Vandersloot [7] with the improved quantization method. In the original Green and Unruh paper, the authors found that their Hamiltonian constraint generically led to solutions that increased in size without bound for large values of the triad parameter p, obviously not in tune with the classical recollapse expected for a closed isotropic model. Their discussion of these results notes the absence of an inner product in their framework to eliminate such unphysical states, but they posited that the lack of a recollapse indicated the LQC program has less contact with the full LQG theory than previous results had indicated. However, the results of Ashtekar et al., using a Hamiltonian constraint with an improved quantization scheme with lattice refinement, put to rest those doubts by finding both a bounce in the quantum regime near the classical singularity, and a recollapse at the appropriate classical scale. In their comparison of the results in the two papers (Section VIIB in [7]), they point out two differences between the constraints -the self-adjoint nature of the improved quantum constraint, and the use of a scalar field by Ashtekar et al. as a time parameter for the evolution.
However, we show here that the Green and Unruh constraint suffers from another problem beyond this -namely the presence of an unstable "numerical" mode that grows without bound, leading to the lack of classical recollapse. We now derive this using von Neumann stability analysis, and show this mode does not exist for the improved quantization scheme of Ashtekar et al. For both models, we make the ansatz for the wavefunction Ψ(m, φ) = ψ m e iωφ .
Here, m is a placeholder for the appropriate variable of the two models (scaled as appropriate), while φ is the value of the scalar field. For the Green-Unruh model, this is µ, which is proportional to the triad component p; for the improved model, this is the volume v. With both equations in this ansatz, the quantum Hamiltonian constraint is of the form C + (m)ψ m+4 + C 0 (m)ψ m + C − (m)ψ m−4 = ω 2 B(m)ψ m .
If in addition, we assume ψ m = λ m/4 in the asymptotic limit m → ∞, where λ is a constant, then this leads to a quadratic equation in λ, namely where C + 0 , C 0 0 , B 0 and C − 0 are the limits of the appropriate function for m → ∞. Using the standard formula, this gives a solution If |λ| > 1, then we have a solution corresponding to this constant which increases in magnitude without bound; this can be easily seen, since in our ansatz, lim m→∞ ψ m+4 ψ m = λ. Note that the scalar field eigenvalue ω drops out, because the inverse volume dependence B(m) m −3/2 is lower order than that of C 0 (m) m 1/2 . This gives values to leading order of λ ± = e iµ 0 1 + 1 2 γ 2 µ 2 0 ± γµ 0 2 4 + γ 2 µ 2 0 .