Semiclassical and quantum behavior of the Mixmaster model in the polymer approach for the isotropic Misner variable

We analyze the semi-classical and quantum behavior of the Bianchi IX Universe in the Polymer Quantum Mechanics framework, applied to the isotropic Misner variable, linked to the space volume of the model. The study is performed both in the Hamiltonian and field equations approaches, leading to the remarkable result of a still singular and chaotic cosmology, whose Poincaré return map asymptotically overlaps the standard Belinskii–Khalatnikov–Lifshitz one. In the quantum sector, we reproduce the original analysis due to Misner, within the revised Polymer approach and we arrive to demonstrate that the quantum numbers of the point-Universe still remain constants of motion. This issue confirms the possibility to have quasi-classical states up to the initial singularity. The present study clearly demonstrates that the asymptotic behavior of the Bianchi IX Universe towards the singularity is not significantly affected by the Polymer reformulation of the spatial volume dynamics both on a pure quantum and a semiclassical level.


Introduction
The Bianchi IX Universe [1,2] is the most interesting among the Bianchi models. In fact, like the Bianchi type VIII, it is the most general allowed by the homogeneity constraint, but unlike the former, it admits an isotropic limit, naturally reached during its evolution [3,4] (see also [5,6]), coinciding with the positively curved Robertson-Walker geometry. Furthermore, the evolution of the Bianchi IX model towards the initial singularity is characterized by a chaotic structure, first outlined in [7] in terms of the Einstein equations morphology and then re-analyzed in Hamiltonian formalism by [8,9]. a e-mail: crino.1547608@studenti.uniroma1.it b e-mail: giovanni.montani@frascati.enea.it c e-mail: pintaudi.1274530@studenti.uniroma1.it Actually, the Bianchi IX chaotic evolution towards the singularity constitutes, along with the same behavior recovered in Bianchi type VIII, the prototype for a more general feature, characterizing a generic inhomogeneous cosmological model [10] (see also [6,[11][12][13][14][15]).
The original interest toward the Bianchi IX chaotic cosmology was due to the perspective, de facto failed, to solve the horizon paradox via such a statistical evolution of the Universe scale factors, from which derives the name, due to C. W. Misner, of the so-called Mixmaster model. However, it was clear since from the very beginning, that the chaotic properties of the Mixmaster model had to be replaced, near enough to the initial singularity (essentially during that evolutionary stage dubbed Planck era of the Universe), by a quantum evolution of the primordial Universe. This issue was first addressed in [8], where the main features of the Wheeler-De Witt equation [16] are implemented in the scenario of the Bianchi IX cosmology. This analysis, besides its pioneering character (see also [6,17]), outlined the interesting feature that a quasi-classical state of the Universe, in the sense of large occupation numbers of the wave function, is allowed up to the initial singularity.
More recent approaches to Canonical Quantum Gravity, like the so-called Loop Quantum Gravity [18], suggested that the geometrical operators, areas and volumes, are actually characterized by a discrete spectrum [19]. The cosmological implementation of this approach led to the definition of the concept of "Big-Bounce" [20,21] (i.e. to a geometrical cutoff of the initial singularity, mainly due to the existence of a minimal value for the Universe Volume) and could transform the Mixmaster Model in a cyclic Universe [22,23]. The complete implementation of the Loop Quantum Gravity to the Mixmaster model is not yet available [24][25][26], but general considerations on the volume cut-off led to characterize the semi-classical dynamics as chaos-free [27]. A quantization procedure, able to mimic the cut-off physics contained in the Loop Quantum Cosmology, is the so-called Polymer quantum Mechanics, de facto a quantization on a discrete lattice of the generalized coordinates [28] (for cosmological implementations see [29][30][31][32][33][34][35][36]).
Here, we apply the Polymer Quantum Mechanics to the isotropic variable α of the Mixmaster Universe, described both in the Hamiltonian and field equations representation. We first analyze the Hamiltonian dynamics in terms of the so called semiclassical polymer equations, obtained in the limit of a finite lattice scale, but when the classical limit of the dynamics forh → 0 is taken. Such a semiclassical approach clearly offers a characterization for the behavior of the mean values of the quantum Universe, in the spirit of the Ehrenfest theorem [37]. This study demonstrates that the singularity is not removed and the chaoticity of the Mixmaster model is essentially preserved, and even enforced in some sense, in this cut-off approach.
Then, in order to better characterize the chaotic features of the obtained dynamics, we translate the Hamiltonian polymer-like dynamics, in terms of the modified Einstein equations, in order to calculate the morphology of the new Poincaré return map, i.e. the modified BKL map.
We stress that, both the reflection rule of the point-Universe against the potential walls and the modified BKL map acquire a readable form up to first order in the small lattice step parameter. Both these two analyses clarify that the chaotic properties of the Bianchi IX model survive in the Polymer formulation and are expected to preserve the same structure of the standard General Relativity case. In particular, when investigating the field equations, we numerically iterate the new map, showing how it asymptotically overlaps to the standard BKL one.
The main merit of the present study is to offer a detailed and accurate characterization of the Mixmaster dynamics when the Universe volume (i.e. the isotropic Misner variable) is treated in a semiclassical Polymer approach, demonstrating how this scenario does not alter, in the chosen representation, the existence of the initial singularity and the chaoticity of the model in the asymptotic neighborhoods of such a singular point.
Finally, we repeat the quantum Misner analysis in the Polymer revised quantization framework and we show, coherently with the semiclassical treatment, that the Misner conclusion about the states with high occupation numbers, still survives: such occupation numbers still behave as constants of motion in the asymptotic dynamics.
It is worth to be noted that imposing a discrete structure to the variable α does not ensure that the Bianchi IX Universe volume has a natural cut-off. In fact, such a volume behaves like e 3α and it does not take a minimal (non-zero) value when α is discretized. This fact could suggest that the surviving singularity is a consequence of the absence of a real geometrical cut-off, like in Loop Quantum Cosmology (of which our approach mimics the semi-classical features).
However, the situation is subtler, as clarified by the following two considerations.
(i) In [30], where the polymer approach is applied to the anisotropy variables β ± , the Mixmaster chaos is removed like in [27] even if no discretization is imposed on the isotropic (volume) variable. Furthermore, approaching the singularity, the anisotropies classically diverge and their discretization does not imply any intuitive regularization, as de facto it takes place there. The influence of the discretization induced by the Polymer procedure, can affect the dynamics in a subtle manner, non-necessarily predicted by the simple intuition of a cut-off physics, but depending on the details of the induced symplectic structure. (ii) In Loop Quantum Gravity, the spectrum of the geometrical spatial volume is discrete, but it must be emphasized that such a spectrum still contains the zero eigenvalue. This observation, when the classical limit is constructed including the suitably weighted zero value, allows to reproduce the continuum properties of the classical quantities. Such a point of view is well-illustrated in [38], where the preservation of a classical Lorentz transformation is discussed in the framework of Loop Quantum Gravity. This same consideration must hold in Loop Quantum Cosmology, for instance in [20], where the position of the bounce is determined by the scalar field parameter, i.e., on a classical level, it depends also on the initial conditions and not only on the quantum modification to geometrical properties of space-time.
Following the theoretical paradigm suggested by the point (ii), the analysis of the Polymer dynamics of the Misner isotropic variable α is extremely interesting because it corresponds to a discretization of the Universe volume, but allowing for the value zero of such a geometrical quantity. We note that picking, as the phase-space variable to quantize, the scale factor a = e α , or any other given power of this, the corresponding Polymer discretization is somehow forced to become singularity-free, i.e. we observe the onset of a Big Bounce. However, in the present approach, based on a logarithmic scale factor, the polymer dynamical scheme offers an intriguing arena to test the Polymer Quantum Mechanics in the Minisuperspace, even in comparison with the predictions of Loop Quantum Cosmology [18].
Finally, we are interested in determining the fate of the Mixmaster model chaoticity, which is suitably characterized by the Misner variables representation and whose precise description (i.e. ergodicity of the dynamics, form of the invariant measure) is reached in terms of the Misner-Chitré-like variables [13,39,40], which are double logarithmic scale factors. Although these variables mix the isotropic and anisotropic Misner variables to some extent, the Misner-Chitré time-like variable in the configurational scale, once discretized would not guarantee a minimal value of the universe volume. Thus the motivations for a Polymer treatment of the Mixmaster model in which only the α dynamics is deformed relies on different and converging requests, coming from the cosmological, the statistical and the quantum feature of the Bianchi IX universe.
The paper is organized as follows. In Sect. 2 we introduce the main kinematic and dynamic features of the Polymer Quantum Mechanics. Then we outline how this singular representation can be connected with the Schrödinger one through an appropriate Continuum Limit.
In Sect. 3 we describe the dynamics of the homogeneous Mixmaster model as it can be derived through the Einstein field equations. A particular attention is devoted to the Bianchi I and II models, whose analysis is useful for the understanding of the general Bianchi IX (Mixmaster) dynamics.
The Hamiltonian formalism is introduced in Sect. 4, where we review the semiclassical and quantum properties of the model as it was studied by Misner in [8].
Section 5 is devoted to the analysis of the polymer modified Mixmaster model in the Hamiltonian formalism, both from a semiclassical and a quantum point of view. Our results are then compared with the ones derived by some previous models.
In Sect. 6 the semiclassical behavior of the polymer modified Bianchi I and Bianchi II models is developed through the Einstein equations formalism, while the modified BKL map is derived and its properties discussed from an analytical and numerical point of view in Sect. 7.
In Sect. 8 we discuss two important physical issues, concerning the link of the polymer representation with Loop Quantum Cosmology and the implications of a polymer quantization of the whole Minisuperspace variables on the Mixmaster chaotic features.
Finally, in Sect. 9 brief concluding remarks follow.

Polymer quantum mechanics
The Polymer Quantum Mechanics (PQM) is a representation of the usual canonical commutation relations (CCR), unitarily nonequivalent to the Schrödinger one. It is a very useful tool to investigate the consequences of the assumption that one or more variables of the phase space are discretized. Physically, it accounts for the introduction of a cutoff. In certain cases where the Schrödinger representation is welldefined, the cutoff can be removed through a certain limiting procedure and the usual Schrödinger representation is recovered, as shown in [29] and summed up in Sect. 2.4. Many people in the Quantum Gravity community think that there is a maximum theoretical precision achievable when measuring space-time distances, this belief being backed up by valuable although heuristic arguments [41][42][43], so that the cutoff introduced in PQM is assumed to be a fundamental quantity. Some results of Loop Quantum Gravity [44] (the discrete spectrum of the area operator) and String Theory [45] (the minimum length of a string) point clearly in the direction of a minimal length scale scenario, too.
PQM was first developed by Ashtekar et al. [28,46,47] who also credit a previous work of Varadarajan [48] for some ideas. It was then further refined also by Corichi et al. [29,49]. They developed the PQM in the expectation to shed some light on the connection between the Planckianenergy Physics of Loop Quantum Gravity and the lowerenergy Physics of the Standard Model.

The Schrödinger representation
Let us consider a quantum one-dimensional free particle with phase space (q, p) ∈ Γ = R 2 . The standard CCR are summarized by the relation q,p = iÎ (1) whereq andp are the position and momentum operators respectively andÎ is the identity operator. These operators form a basis for the so called Heisenberg algebra [50]. At this stage we have not made any assumption regarding the Hilbert space of the theory, yet.
To introduce the Polymer representation in Sect. 2.2, it is convenient to consider also the Weyl algebra W , generated by exponentiation ofq andp where α and β are real parameters with the dimension of momentum and length, respectively. The CCR (1) becomê where · indicates the product operation of the algebra. All other product combinations commute. A generic element of the Weyl algebra W ∈ W is a finite linear combination where A i and B i are complex coefficients and the product (3) is extended by linearity. It can be shown that W has a structure of a C * -algebra.
The Schrödinger representation is obtained as soon as we introduce the Schrödinger Hilbert space of square Lebesgueintegrable functions H S = L 2 (R, dq) and the action of the bases operators on it: qψ(q) = qψ(q);pψ(q) = −i∂ q ψ(q) (5) where ψ ∈ L 2 (R, dq) and q ∈ R.

The polymer representation: kinematics
The Polymer representation of Quantum Mechanics can be introduced as follows. First, we consider the abstract states |μ of the Polymer Hilbert space H poly labeled by the real parameter μ. A generic state of H poly consists of a finite linear combination where N ∈ N and we assume the fundamentals kets to be orthonormal The Hilbert space H poly is the Cauchy completion of the vectors (6) with respect to the product (7). It can be shown that H poly is nonseparable [29]. Two fundamental operators can be defined on this Hilbert space. The "label" operatorε and the "displacement" operatorŝ(λ) with λ ∈ R, which act on the states |μ aŝ ε|μ := μ|μ ;ŝ(λ)|μ := |μ + λ (8) Theŝ operator is discontinuous in λ since for each value of λ the resulting states |μ + λ are orthogonal. Thus, there is no Hermitian operator that by exponentiation could generate the displacement operator. We denote the wave functions in the p-polarization as ψ( p) = p|ψ , where ψ μ ( p) = p|μ = e iμp . TheV (λ) operator defined in (2) shifts the plane waves by λ As expected,V (λ) can be identified with the displacement operatorŝ(λ) and thep operator cannot be defined rigorously.
On the other hand the operatorq is defined bŷ and thus can be identified with the abstract displacement operatorε. The reasonq is said to be discrete is that the eigenvalues of this operator are the labels μ, and even when μ can take value in a continuum of possible values, they can be regarded as a discrete set, because the states are orthonormal for all values of μ.
TheV (λ) operators form a C * -Algebra and the mathematical theory of C * -Algebras [51] provide us with the tools to characterize the Hilbert space which the wave functions ψ( p) belong to. It is given by the Bohr compactification R b of the real line [52]: H poly = L 2 (R b , dμ H ), where the Haar measure dμ H is the natural choice.
In terms of the fundamental functions ψ μ ( p) the inner product takes the form

Dynamics
Once a definition of the polymer Hilbert space is given, we need to know how to use it in order to analyze the dynamics of a system. The first problem to be faced is that, as it was said after equation (9), the polymer representation isn't able to describep andq at the same time. It is then necessary to choose the "discrete" variable and to approximate its conjugate momentum with a well defined and quantizable function. For the sake of simplicity let us investigate the simple case of a one-dimensional particle described by the Hamiltonian: in p-polarization. If we assume thatq is the discrete operator (in the sense explained in Sect. 2.2), we then need to approximate the kinetic term p 2 2m in an appropriate manner. The procedure outlined in [29] consists in the regularization of the theory, through the introduction of a regular graph γ μ 0 (that is a lattice with spacing μ 0 ), such as: It follows that the only states we shall consider, in order to remain in the graph, are that of the form |ψ = n b n |μ n ∈ H γ μ 0 , where the b n are coefficients such as n |b n | 2 < ∞ and the |μ n (with μ n = nμ 0 ) are the basic kets of the Hilbert space H γ μ 0 , which is a separable subspace of H poly . Moreover the action of all operators on these states will be restricted to the graph. Since the exponential operator V (λ) = e iλp can be identified with the displacement onê s (Eq. (9)), it is possible to use a regularized version of it (V (μ 0 ) such thatV (μ 0 )|μ n = |μ n + μ 0 = |μ n+1 ) in order to define an approximated version ofp: This definition is based on the fact that, for p << 1 μ 0 , one gets p 1 μ 0 sin(μ 0 p) = 1 2iμ 0 (e iμ 0 p − e −iμ 0 p ). From this result it is easy to derive the approximated kinetic termp 2 μ 0 by applying the same definition ofp μ 0 two times: It is now possible to introduce a regularized Hamiltonian operator, which turns out to be well-defined on H γ μ 0 and symmetric: In the p-polarization, then,p 2 μ 0 acts as a multiplication operator (p 2 μ 0 ψ( p) = 1 μ 2 0 sin 2 (μ 0 p)ψ( p)), whileq acts as a derivative operator (qψ( p) = i∂ p ψ( p)).

Continuum limit
As we have seen in the previous section, the condition p << 1 μ 0 is a fundamental hypothesis of the polymer approach. If this limit is valid for all the orbits of the system to be quantized, such a representation is expected to give results comparable to the ones of the standard quantum representation (based on the Schrödinger equation). This is the reason why a limit procedure from the polymer system to the Schrödinger one is necessary to understand how the two approaches are related. In other words we need to know if, given a polymer wave function defined on the graph γ 0 (with spacing a 0 ), it is possible to find a continuous function which is best approximated by the previous one in the limit when the graph becomes finer, that is: The first step to take in order to answer this question consists in the introduction of the scale C n which is a decomposition of the real line in the union of closed-open intervals with the lattice γ n points as extrema, that cover the whole line and do not intersect. For each scale we can then define an effective theory based on the fact that every continuous function can be approximated by another function that is constant on the intervals defined by the lattice. It follows that we can create a link between H poly and H S for each scale C n : where χ α m is the characteristic function on the interval α m = [ma n , (m + 1)a n ). Thus it is possible to define an effective HamiltonianĤ C n for each scale C n . The set of effective theories is then analyzed by the use of a renormalization procedure [29], and the result is that the existence of a continuum limit is equivalent to both the convergence of the energy spectrum of the effective theory to the continuous one and the existence of a complete set of renormalized eigenfunctions.
and the constants of structure λ l , λ m , λ n of the inequivalent isometry group that characterize each Bianchi model, where "a, b, c" are said cosmic scale factors and describe the behavior in the synchronous time t of the three independent spatial directions [1].
where the logarithmic variable q l (τ ), q m (τ ), q n (τ ) and logarithmic time τ are defined as such: q l = 2 ln a, q m = 2 ln b, q n = 2 ln c, dt = (abc)dτ (21) and the subscript τ is a shorthand for the derivative in τ . The isometry group characteristic of Bianchi IX is the SO (3) group. The Class A Bianchi models [13] (of which Bianchi I,II and IX are members) are set apart just by those three constants of structure.

Bianchi I
The Bianchi I model in void is also called Kasner solution [53]. By substituting in Eq. (19) the Bianchi I constants of structure λ l = 0, λ m = 0, λ n = 0, we get the simple equations of motion: whose solution can be given as a function of the logarithmic time τ and four parameters: where p l , p m , p n are said Kasner indices and Λ is a positive constant that will have much relevance in Sect. 6.1. It is needed to ensure that the sum of Kasner indices is always one: From (23) and (21), if we exploit the constraint (24), we obtain From relation (25) stems the appellative "logarithmic time" for the time variable τ . By substituting the Bianchi I equations of motion (22) and their solution (23) into (20), after applying the constraint (24), we find the additional constraint With the exception of the null measure set ( p l , p m , p n ) = (0, 0, 1) and ( p l , p m , p n ) = (−1/3, 2/3, 2/3), Kasner indices are always different from each other and one of them is always negative. The (0, 0, 1) case can be shown to be equivalent to the Minkowski space-time. It is customary to order the Kasner indices from the smallest to the greatest p 1 < p 2 < p 3 . Since the three variables p 1 , p 2 , p 3 are constrained by Eqs. (24) and (26), they can be expressed as functions of a unique parameter u as The range and the values of the Kasner indices are portrayed in Fig. 1. Excluding the Minkowskian case (0, 0, 1), for any choice of the Kasner indices, the spatial metric η = diag(t 2 p l , t 2 p m , t 2 p n ) has a physical naked singularity, called also Big Bang, when t = 0 (or τ → −∞).

Bianchi II
A period in the evolution of the Universe when the r.h.s of (19) can be neglected, and the dynamics is Bianchi I-like, is called Kasner regime or Kasner epoch. In this section we show how Bianchi II links together two different Kasner epochs at τ → −∞ and τ → ∞. A series of successive Kasner epochs, where one cosmic scale factor increases monotonically, is called a Kasner era. Again, after substituting the Bianchi II structure constants It should be noted that the conditions hold for every τ .
Let us consider the explicit solutions of Eq. (29) where c 1 , . . . , c 6 are integration constants. We now analyze how (31) behave as the time variable τ approaches +∞ and −∞ (remembering also Eq. (25)). By taking the asymptotic limit to ±∞, we note that a Kasner regime at +∞ is "mapped" to another Kasner regime at −∞, but with different Kasner indices.
where we have exploited the notation of Eq. (23). Both these sets of indices satisfy the two Kasner relations (24) and (26).
The old (the unprimed ones) and the new indices (the primed ones) are related by the BKL map, that can be given in the form of the new primed Kasner indices p l , p m , p n and Λ as functions of the old ones. To calculate it we need at least four relations between the new and old Kasner indices Then we can invert these relations to find the BKL map.
One f i is simply the sum of the primed Kasner indices [that is the primed version of (24)]. Other two relations can be obtained from (30): The fourth relation is obtained by direct comparison of the asymptotic limits (32), and for this reason we will call it asymptotic relation.
All other relations that can be obtained from (32) are equivalent to (35). Finally, by inverting the four relations (24), (33), (34), (35) and assuming again that p l is the negative Kasner index p 1 , we finally get the classical BKL map: The main feature of the BKL map is the exchange of the negative index between two different directions. Therefore, in the new epoch, the negative power is no longer related to the l-direction and the perturbation to the new Kasner regime [which is linked to the l-terms on the r.h.s. of Eq. (19)] is damped and vanishes towards the singularity: the new (primed) Kasner regime is stable towards the singularity. If we then insert the parametrization (27) and (28) into the BKL map (36) just found, we get that the map in the u parameter is simple-looking.
We emphasize that the role of the negative index is swapped between the l and the m-direction.

Bianchi IX
Here we briefly explain how the BKL map can be used to characterize the Bianchi IX dynamics. This is possible because it can be shown [13] that the Bianchi IX dynamical evolution is made up piece-wise by Bianchi I and II-like "blocks".
First, we take a look again at the Einstein equations (19): but now we substitute the Bianchi IX structure constants λ l = 1, λ m = 1, λ n = 1 into them to get the Einstein equations for Bianchi IX: Equation (20) stays valid for every Bianchi model. Again, we assume an initial Kasner regime, where the negative Kasner index p 1 is associated with the l-direction. Thus, the "exploding" cosmic scale factor in the r.h.s. of (38) is identified again with a. By retaining in the r.h.s. of (38) only the terms that grow towards the singularity, we obtain a system of ordinary differential equations in time completely similar to the one just encountered for Bianchi II (29), with the only caveat that now there is no initial condition (i.e. no choice of the initial Kasner indices) such that the initial or final Kasner regimes are stable towards the singularity.
This because all the cosmic scale factors {a, b, c} are treated on an equal footing in the r.h.s. of (38). After every Kasner epoch or era, no matter on which axis the negative index is swapped, there will always be a perturbation in the Einstein equations (38) that will make the Kasner regime unstable. We thus have an infinite series of Kasner eras alternating and accumulating towards the singularity.
Given an initial u parameter, the BKL map that takes it to the next one u is Let us assume an initial x 0 is irrational, as it is statistically and more general to infer, after the k 0 Kasner epochs that make up that Kasner era, its inverse 1/x 0 would be irrational, too. This inversion happens infinitely many times until the singularity is reached and is the source of the chaotic features of the BKL map [39,[55][56][57][58][59].

Hamiltonian formulation of the Mixmaster model
In this section we summarize some important results obtained applying the Arnowitt, Deser and Misner (ADM) Hamiltonian methods to the Bianchi models. The main benefit of this approach is that the resulting Hamiltonian system resembles the simple one of a pinpoint particle moving in a potential well. Moreover it provides a method to quantize the system through the canonical formalism. The line element for the Bianchi IX model is: where N is the Lapse Function, characteristic of the canonical formalism and σ i are 1-forms depending on the Euler angles of the SO(3) group of symmetry. β i j is a diagonal, traceless matrix and so it can be parameterized in terms of the two independent variables β ± as The factor 1 4 is chosen (following [8]) so that when β i j = 0, we obtain just the standard metric for a three-sphere of radius r = e α . Thus for β = 0 this metric is the Robertson-Walker positive curvature metric. The variables (α, β ± ) were introduced by Misner [9]; α describes the expansion of the Universe, i.e. its volume changes (V ∝ e 3α ), while β ± are related to its anisotropies (shape deformations). The kinetic term of the Hamiltonian of the Bianchi models, written in the Misner variables, results to be diagonal; we can then write, following [60], the super-Hamiltonian for the Bianchi models simply as: where ( p α , p ± ) are the conjugate momenta to (α, β ± ) respectively and N is the lapse function. It should be noted that the corresponding super-Hamiltonian constraint H B = 0 describes the evolution of a point β = (β + , β − ), that we call β-point, in function of a new time coordinate α (that is the shape of the Universe in function of its volume). Such point is subject to the anisotropy potential V B (β ± ) which is a non-negative function of the anisotropy coordinates, whose explicit form depends on the Bianchi model considered. In particular Bianchi type I corresponds to the V B = 0 case (free particle), Bianchi type II potential corresponds to a single, infinite exponential wall V B = e −8β + , while for Bianchi type IX we have a closed domain expressed by: The potential (42) (41), moreover, causes the potential walls to move outward as we approach the cosmological singularity (α → −∞).
Following the ADM reduction procedure, we can solve the super-Hamiltonian constraint (41) with respect to the momentum conjugate to the new time coordinate of the phase-space (i.e. α). The result is the so-called reduced Hamiltonian H ADM : Here we can see that far from the origin the function has the symmetry of an equilateral triangle, while near the origin (V B < 1) the equipotentials are closed curves from which we can derive the classical dynamics of the model by solving the corresponding Hamilton's equations: It should be noted that the choiceα = 1 fixes the temporal gauge N ADM = 6(4π) 2 H ADM κ e 3α . Since we are interested in the dynamics of the model near the singularity (α → −∞) and because of the steepness of the potential walls, we can consider the β-point to spend most of its time moving as a free particle (far from the potential is non-negligible only during the short time of a bounce against one of the three walls. It follows that we can separate the analysis into two different phases. As far as the free motion phase (Bianchi type I approximation) is concerned, we can derive the velocity of the β-point from the first of Eq. (44) (with V B = 0), obtaining: It remains to study the bounce against the potential. Following [8] we can summarize the results of this analysis in three points: (i) The position of the potential wall is defined as the position of the equipotential in the β-plane bounding the region in which the potential terms are significant. The wall turns out to have velocity |β wall | ≡ dβ wall dα = 1 2 , that is one half of the β-point velocity. So a bounce is indeed possible. (ii) Every bounce occurs according to the reflection law: where θ i and θ f are the incidence and the reflection angles respectively. (iii) The maximum incidence angle for the βpoint to have a bounce against a specific potential wall turns out to be: Because of the triangular symmetry of the potential, this result confirms the fact that a bounce against one of the three walls always happens. It follows that the β-point undergoes an infinite series of bounces while approaching the singularity and, sooner or later, it will assume all the possible directions, regardless the initial conditions. In short we are dealing with a pinpoint particle which undergoes an uniform rectilinear motion marked by the constants of motion ( p + , p − ) until it reach a potential wall. The bounce against this wall has the effect of changing the values of p ± . After that, a new phase of uniform rectilinear motion (with different constants of motion) starts. If we remember that the free particle case corresponds to the Bianchi type I model, whose solution is the Kasner solution (see Sect. 3.1), we recover the sequence of Kasner epochs that characterizes the BKL map (Sect. 3.3).
In conclusion, it is worth mentioning that, through the analysis of the geometrical properties of the scheme outlined so far near the cosmological singularity, we can obtain a sort of conservation law, that is: Here the symbol ... denotes the average value of H ADM α over a large number of runs and bounces. In fact, although it is not a constant, this quantity acquires the same constant value just before each bounce: H n ADM α n = H n+1 ADM α n+1 . In this sense quantities with this property will be named adiabatic invariants.

Quantization of the Mixmaster model
Near the Cosmological Singularity some quantum effects are expected to influence the classical dynamics of the model. The use of the Hamiltonian formalism enables us to quantize the cosmological theory in the canonical way with the aim of identify such effects.
Following the canonical formalism [8], we can introduce the basic commutation relations [β a ,p b ] = iδ ab , which can be satisfied by choosingp a = −i ∂ ∂β a . Hence all the variables become operators and the super-Hamiltonian constraint (41) is written in its quantum version (Ĥ Ψ (α, β ± ) = 0) i.e. the Wheeler-De Witt (WDW) equation: The function Ψ (α, β ± ) is the wave function of the Universe, which we can choose of the form: If we assume the validity of the adiabatic approximation (see [30]) we can solve the WDW equation by separation of variables.
In particular the eigenvalues E n of the reduced Hamiltonian (43) are obtained from those of the equationĤ ADM φ n = E 2 n φ n that is: The result, achieved by approximating the triangular potential with a quadratic box with infinite vertical walls and the same area of the triangular one (A = 3 4 √ 3α 2 ), turns out to be: where n, m are the integer and positive quantum numbers related to the anisotropies β ± . An important conclusion that Misner derived from this analysis is that quasi-classical states (i.e. quantum states with very high occupation numbers) are preserved during the evolution of the Universe towards the singularity. In fact if we substitute the quasi-classical approximation H E n in (48) and we study it in the limit α → −∞, we find: We can therefore conclude that if we assume that the present Universe is in a quasi-classical state of anisotropy (n 2 + m 2 >> 1), then, as we extrapolate back towards the Cosmological Singularity, the quantum state of the Universe remains always quasi-classical.

Mixmaster model in the polymer approach
Here the Hamiltonian formalism introduced in Sect. 4 is used for the analysis of the dynamics of the modified Mixmaster model. We choose to apply the polymer representation to the isotropic variable α of the system due to its deep connection with the volume of the Universe. As a consequence we will need to find an approximated operator for the conjugate momentum p α , while the anisotropy variables β ± will remain unchanged from the standard case.
The introduction of the polymer representation consists in the formal substitution, derived in Sect. 2.3: which we apply to the super-Hamiltonian constraint In what follows we will use the ADM reduction formalism in order to study the new model both from a semiclassical 1 and a quantum point of view.

Semiclassical analysis
First of all we derive the reduced polymer Hamiltonian (H α := − p α ) from (56): where the condition 0 ≤ μ 2 ( p 2 + + p 2 − + 3(4π) 4 κ 2 e 4α V B (β ± )) ≤ 1 has to be imposed due to the presence of the function arcsine. The dynamics of the system is then described by the deformed Hamilton's equations (that come from Eq. (44) with H ADM ≡ H α ): where the prime stands for the derivative with respect to α and p 2 = p 2 + + p 2 − .
As for the standard case, we start by studying the simplest case of V B = 0 (Bianchi I approximation), which corresponds to the phase when the β-point is far enough from the potential walls. Here the persistence of a cosmological singularity can be easily demonstrated since α results to be linked to the time variable t through a logarithmic relation [as we will see later in Eq. (90a)]: α ∼ ln(t) −→ t→0 −∞. Moreover the anisotropy velocity of the particle can be derived from Eq. (58a) (while p ± and H α are constants of motion) and the result is: Thus the modified anisotropy velocity depends on the constants of motion of the system, and it turns out to be always greater than the standard one (r α (μ, p ± ) > 1, ∀ p ± ∈ R). On the other side the velocity of each potential wall is the same of the standard case |β wall | = 1 2 . In fact the Bianchi II potential (i.e. our approximation of a single wall of the Bianchi IX potential) can be written in terms of the determinant of the metric ( √ η = e 3α ) as which, in the limit α → −∞, tends to the Θ-function: It follows that the position of the wall is defined by the equation 4 3 − 8β + 3α = 0, i.e. |β wall | = α 2 . The β-particle is hence faster than the potential, therefore a bounce is always possible also after the polymer deformation. For this reason, and from the fact that the singularity is not removed, we can state that there will certainly be infinite bounces onto the potential walls. Moreover, the more the point is moving fast with respect to the walls, the more the bounces frequency is high. All these facts, are hinting at the possibility that chaos is still present in the Bianchi IX model in our Polymer approach. It is worth noting that this is a quite interesting result if we consider the strong link between the Polymer and the Loop Quantum Cosmology, which tends to remove the chaos and the cosmological singularity [61,62]. In particular, when dealing with PQM, a spatial lattice of spacing μ must be introduced to regularize the theory (Sect. 2.4). Hence, one would naively think that the Universe cannot shrink past a certain threshold volume (singularity removal). Actually we know that Loop Quantum Cosmology applied both to FRW [61] and Bianchi I [62] predicts this Big Bounce-like dynamics for the primordial Universe. Anyway the persistence of the cosmological singularity in our model seems to be a consequence of our choice of the configurational variables, while its chaotic dynamics towards such a singularity is expected to be a more roboust physical feature. In order to analyze a single bounce we need to parameterize the anisotropy velocity in terms of the incident and reflection angles (θ i and θ f , shown in Fig. 4): As usual the subscripts i and f distinguish the initial quantities (just before the bounce) from the final ones (just after the bounce). Now we can derive the maximum incident angle for having a bounce against a given potential wall. If we consider, for example, the left wall of Fig. 2 we can observe that a bounce occurs only when (β + ) i > β wall = 1 2 , that is: Two important features should be underlined about this result: the first one is that θ α max is bigger than the standard maximum angle θ max = π 3 , as we can immediately see from its second order expansion. The second one is that when the anisotropy velocity tends to infinity one gets θ α max → π 2 (Fig. 3). Both these characteristics along with the triangular symmetry of the potential confirm the presence of an infinite series of bounces. Each bounce involves new values for the constants of motion of the free particle ( p ± , H α ) and a change in its direction. The relation between the directions before and after a bounce can be inferred by considering the Hamilton's equations (58) again, this time with V B = e −8β + . First of all we observe that V B depends only from β + , from which we deduce that p − is yet a constant of motion. A second constant of motion is K := H α − 1 2 p + , as we can verify from the Hamilton equations (58). An expression for p ± can be derived from Eq. (58a): Fig. 3 The maximum angle to have a bounce as a function of the anisotropy velocity r α ∈ [1, ∞). The two limit values π 3 (r α = 1) and π 2 (r α → ∞) are highlited where the compact notation μ 2 ( p 2 + 3(4π) 4 κ 2 e 4α−8β + ) = sin(μH α ) is used. By the substitution of this expression and of the parameterization (62) in the two constants of motion we obtain two modified conservation laws: The reflection law can be then derived after a series of algebraic passages which involve the substitution of (65a) in (65b), and the use of the explicit expression (59) for the anisotropy velocity r α : In order to have this law in a simpler form, comparable to the standard case, an expansion up to the second order in μp << 1 is required. The final result is: where we have defined Π 2 = 1 6 μ 2 p 2 .

Comparison with previous models
As we stressed in Sect. 5.1 our model turns out to be very different from Loop Quantum Cosmology [61,62], which predicts a Big Bounce-like dynamics. Moreover, if the polymer approach is applied to the anisotropies β ± , leaving the volume variable α classical [30], a singular but non-chaotic dynamics is recovered. It should be noted that the polymer approach in the perturbative limit can be interpreted as a modified commutation relation [63]: where μ > 0 is the deformation parameter. Another quantization method involving an analogue modified commutation relation is the one deriving by the Generalized Uncertainty Principle (GUP): through which a minimal length is introduced and which is linked to the commutation relation: In [64] the Mixmaster model in the GUP approach (applied to the anisotropy variables β ± ) is developed and the resulting dynamics is very different from the one deriving by the application of the polymer representation to the same variables [30]. The GUP Bianchi I anisotropy velocity turns out to be always greater than 1 (β 2 GU P = 1 + 6sp 2 + 9s 2 p 4 ) and, in particular, greater than the wall velocity β GU P wall (which in this case is different from 1 2 ). Moreover the maximum angle of incidence for having a bounce is always greater than π 3 . Therefore the occurrence of a bounce between the β-particle and the potential walls of Bianchi IX can be deduced. The GUP Bianchi II model is not analytic (no reflection law can be inferred), but some qualitative arguments lead to the conclusion that the deformed Mixmaster model can be considered a chaotic system. In conclusion if we apply the two modified commutation relations to the anisotropy (physical) variables of the system we recover very different behaviors. Instead the polymer applied to the geometrical variable α leads to a dynamics very similar to the GUP one described here: the anisotropy velocity (59) is always greater than 1 and the maximum incident angle (63) varies in the range π 3 < θ max < π 2 , so that the chaotic features of the two models can be compared.

Quantum analysis
The modified super-Hamiltonian constraint (56) can be easily quantized by promoting all the variables to operators and by approximating p α through the substitution (55): If we consider the potential walls as perfectly vertical, the motion of the quantum β-particle boils down to the motion of a free quantum particle with appropriate boundary conditions. The free motion solution is obtained by imposing V B = 0 as usual, and by separating the variables in the wave function: ψ( p α , p ± ) = χ( p α )φ( p ± ). It should be noted that, although we can no longer take advantage of the adiabatic approximation (51), due to the necessary use of the momentum polarization in Eq. (71), nonetheless we can suppose that the separation of variables is reasonable. In fact in the limit where the polymer correction vanishes (μ → 0), the Misner solution (for which the adiabatic hypothesis is true) is recovered. The anisotropy component φ( p ± ) is obtained by solving the eigenvalues problem: where k 2 = k 2 + + k 2 − (k ± are the eigenvalues ofp ± respectively). Therefore the eigenfunctions are φ( where A, B, C, D are integration constants. By substituting this result in (71) we find also the isotropic component χ( p α ) and the related eigenvaluep α : The potential term is approximated by a square box with vertical walls, whose length is expressed by L(α) = L 0 + α, so that they move outward with velocity |β wall | = 1 2 : The system is then solved by applying the boundary conditions: to the eigenfunctions (73), expressed in coordinate representation (i.e. φ(β ± ), which can be calculated through a standard Fourier transformation). The final solution is: with: Instead the isotropic component results to be: We emphasize the fact that, as we announced at the beginning of this section, the Misner solution (53) can be easily recovered from Eq. (80), by taking the limit μ → 0.

Adiabatic invariant
Because of the various hints about the presence of a chaotic behavior, outlined in Sect. 5.1 and that will be confirmed by the analysis of the Poincaré map in Sect. 7, we can reproduce a calculation similar to the Misner's one [8] in order to obtain information about the quasi-classical properties of the early universe. In Sect. 5.1 we found that the β-point undergoes an infinite series of bounces against the three potential walls, moving with constant velocity between a bounce and another. Every bounce implies a decreasing of H α due to the conservation of K = H α − 1 2 p + and the definition of a new direction of motion through Eq. (67). Because of the ergodic properties of the motion we can assume that it won't be influenced by the choice of the initial conditions. Thus we can choose them in order to simplify the subsequent analysis. A particularly convenient choice is: θ i + θ f = 60 • for the first bounce. In fact as a consequence of the geometrical properties of the system all the subsequent bounces will be characterized by the same angle of incidence: θ i = θ i (Fig. 4). Let us now analyze a single bounce, keeping in mind Fig. 5. At the time α in which we suppose the bounce to occur, the wall is in the position β wall = |α| 2 . Moving between two walls, the particle has the constant velocity |β | = r α (Eq. (59)), so that we can deduce the spatial distance between two bounces in terms of the "time" |α| necessary to cover such a distance: Now if we apply the Sine theorem to the triangles A BC and A B D, we find the relation: Moreover, from (65a) one gets: From these two relations we can then deduce a quantity which remains unaltered after the bounce:  Now we can generalize the argument to the n-th bounce (taking place at the "time" α n instead of α). By using again the Sine theorem (applied to the triangle A B D) we get the relation r α i |α i | = sin 120 • 2 sin θ f |α n | = C|α n |, where C takes the same value for every bounce. Correspondingly r α f |α f | = C|α n+1 |. If we substitute these new relations in Eq. (85), we can then conclude that: where the average value is taken over a large number of bounces. As we are interested in the behavior of quantum states with high occupation number, remembering that for such states the quasi-classical approximation H α p α is valid, we can substitute Eq. (80) and Eq. (59) in (86) and evaluate the correspondent quantity in the limit α → −∞. The result is that also in this case the occupation number is an adiabatic invariant: that is the same conclusion obtained by Misner without the polymer deformation. It follows that also in this case a quasiclassical state of the Universe is allowed up to the initial singularity.

Semiclassical solution for the Bianchi I and II models
As we have already outlined in Sect. 3, to calculate the BKL map it is necessary first to solve the Einstein's Eqs for the Bianchi I and Bianchi II models. Therefore, in Sect. 6.1 we solve the Hamilton's Eqs in Misner variables for the Polymer Bianchi I, in the semiclassical approximation introduced in Sect. 5.1.
Then in Sect. 6.2 we derive a parametrization for the Polymer Kasner indices. It is not possible anymore to parametrize the three Kasner indices using only one parameter u as in (27), but two parameters are needed. This is due to a modification to the Kasner constraint (26).
In Sect. 6.3 we calculate how the Bianchi II Einstein's Eqs are modified when the PQM is taken into account. Then we derive a solution for these Eqs, using the lowest order perturbation theory in μ. We are therefore assuming that μ is little with respect to the Universe volume cubic root.

Polymer Bianchi I
As already shown in Sect. 5.1, the Polymer Bianchi I Hamiltonian is obtained by simply setting V B = 0 in (56). Upon solving the Hamilton's Eqs derived from the Hamiltonian (56) in the case of Bianchi I, variables p α and p ± are recognized to be constants throughout the evolution. We then choose the time gauge so as to have a linear dependence of the volume V from the time t with unity slope, as in the standard solution (23): Moreover, to have a Universe expanding towards positive times, we must restrict the range of sin(2μp α ) < 0. This condition reflects on p α as The branch connected to the origin n = 0, will be called connected-branch. For μ → 0, it gives the correct continuum limit: p α < 0. Even though, for different choices of n = 0, we cannot provide a simple and clear explanation of their physical meaning, when we will define the parametrization of the Kasner indices in Sect. 6.2, we will be forced to consider another branch too, in addition to the connected-one. We choose arbitrarily the n = 1 branch and, by analogy, we will call it the disconnected-branch.
The solution to Hamilton's Eqs for Polymer Bianchi I in the (88) time gauge, reads as Since p α , p ± and μ in the (90) are all numerical constants, it is qualitatively the very same as the classical solution (23).
In particular, even in the Polymer case the volume V (t) ∝ e 3α(t) = t goes to zero when t → 0: i.e. the singularity is not removed and we have a Big Bang, as already pointed out in Sect. 5.1. By making a reverse canonical transformation from Misner variables (α, β ± ) to ordinary ones (q l , q m , q n ), we find out that (24) holds unchanged while (26) is modified into: If the continuum limit μ → 0 is taken, we get back the old constraint (26) as expected. Since the relation holds, the inequality p l 2 + p m 2 + p n 2 ≥ 1, ∀ p α ∈ R is directly linked to β-point velocity (59) being always equal or greater than one, as demonstrated in Sect. 5.1. This is a noteworthy difference with respect to the standard case of Eq. (45), where the speed was always one.
We can derive an alternative expression for the Kasner constraint (91) by exploiting the notation of Eq. (23): where the plus sign is for the connected-branch while the minus for the disconnected-branch. We have defined the dimensionless quantity Q as It is worth noticing that the lattice spacing μ has been incorporated in Q, so that the continuum limit μ → 0 is now equivalently reached when Q → 0. Condition (89) reflects on the allowed range for Q: both for the connected and disconnected branches.

Parametrization of the Kasner indices
In this section we define a parametrization of the Polymer Kasner indices. In the standard case of (27), because there are three Kasner indices and two Kasner constraints (24) and (26), only one parameter u is needed to parametrize the Kasner indices. On the other hand, in the Polymer case, even if there are two Kasner constraints (24) and (93) as well, constraint (93) already depends on another parameter Q. This means that any parametrization of the Polymer Kasner indices will inevitably depend on two parameters, that we arbitrarily choose as u and Q, where u is defined on the same range as the standard case (27) and Q was defined in (94). They are both dimensionless. We will refer to the following expressions as the (u, Q)-parametrization for the Polymer Kasner indices: where the plus sign is for the connected branch and the minus sign is for the disconnected branch. By construction, the standard u-parametrization (27) is recovered in the limit Q → 0 of the connected branch. We will see clearly in Sect. 7.2 how the Q parameter can be thought of as a measure of the quantization degree of the Universe. The more the Q of the connected-branch is big in absolute value, the more the deviations from the standard dynamics due to PQM are pronounced. The opposite is true for the disconnected-branch.
To gain insight on multiple-parameters parametrizations of the Kasner indices, the reader can refer to [65,66]. The (u, Q)-parametrization (96) is even in Q: p a (u, Q) = p a (u, −Q) with a = l, m, n. We can thus assume without loss of generality that 0 ≤ Q ≤ 1. Another interesting feature of the Polymer Bianchi I model, that is evident from Fig. 6, is that, for every u and Q in a non-null measure set, two Kasner indices can be simultaneously negative [instead of only one as is the case of the standard u-parametrization (27)].
In a similar fashion as the u parameter is "remapped" in the standard case [second line of (39)], also the (u, Q) parametrization is to be "remapped", if the u parameter happens to become less than one. Down below we list these remapping prescriptions explicitly: we show how the (u, Q)parametrization (96) can be recovered if the u parameter becomes smaller than 1, through a reordering of the Kasner indices and a remapping of the u parameter.
where the indices on the left are defined for the u shown after the bullet •, while the indices on the right are in the "correct" range u > 1. The arrow '→' is there to indicate a conceptual mapping. However, in value it reads as an equality '='. As far as the Q parameter is concerned, a "remapping" prescription for Q is not needed because values outside the boundaries (95) are not physical and should never be considered. In Fig. 6 the values of the ordered Kasner indices are displayed for the (u, Q)-parametrization (96), where Q ∈ [0, 1]. Because the range u > 1 is not easily plottable, the equivalent parametrization in u ∈ [0, 1] was used. We notice that the roles of p m and p n are exchanged for this range choice.

Polymer Bianchi II
Here we apply the method described Sect. in 6.1 to find an approximate solution to the Einstein's equations of the Polymer Bianchi II model. We start by selecting the V B potential appropriate for Bianchi II (60) and we substitute it in the Hamiltonian (56) (in the following we will always assume the time gauge N = 1).
Then, starting from Hamilton's equations, inverting the Misner canonical transformation and converting the synchronous time t-derivatives into logarithmic time τ -derivatives, the Polymer Bianchi II Einstein's equations are found to be: where v(τ ) := q l τ (τ ) + q m τ (τ ) + q n τ (τ ). In the Q 1 approximation it is enough to consider only the connectedbranch of the solutions, i.e. the branch that in the limit μ → 0 make Eq. (97) reduce to the correct classical Eq. (29). Now we find an approximate solution to the system (97), with the only assumption that μ is little compared to the cubic root of the Universe volume. We are entitled then to exploit the standard perturbation theory and expand the solution until the first non-zero order in μ. Because μ appears only squared μ 2 in the Einstein's Eq. (97), the perturbative expansion will only contain even powers of μ.
First we expand the solution at the first order in where q 0 l,m,n are the zeroth order terms and q 1 l,m,n are the first order terms in μ 2 . We recall the zeroth order solution q 0 l,m,n is just the classical solution (31), where the q l,m,n (τ ) there must be now appended a 0-apex, accordingly.
Considering the fact that q m τ τ = q n τ τ , we can simplify q n (τ ) from (97), by setting where the first six constants of integration c 1 , c 2 , . . . , c 6 play the same role at the zeroth order as in Eq. (31), while the successive six c 7 , . . . , c 12 are needed to parametrize the first order solution.
Then we substitute (99) in (97) and, as it is required by perturbations theory, we gather only the zeroth and first order terms in μ 2 and neglect all the higher order terms. The Einstein's Eqs for the first order terms q 1 l and q 1 m are then found to be: where we have defined C ≡ 2(4π) 2 κ 2 for brevity and we have substituted the zeroth order solution (31) where needed. Equation (100) are two almost uncoupled nonhomogeneous ODEs. Equation (100a) is a linear second order non-homogeneous ODE with non-constant coefficients and, being completely uncoupled, it can be solved straightly. Equation (100b) is solved by mere substitution of the solution of (100a) in it and subsequent double integration.
To solve (100a) we exploited a standard method that can be found, for example, in [67,Lesson 23]. This method is called reduction of order method and, in the case of a second order ODE, can be used to find a particular solution once any non trivial solution of the related homogeneous equation is known.
The homogeneous equation associated to (100a) is whose general solution reads as: (102) By applying the above-mentioned method, we obtain the following solution for (100a) The solution for (100b) is found by substituting (103) in it and then integrating two times: Now that we know q 0 l , q 1 l , q 0 m and q 1 m , the complete solution for q l , q m and q n is found through (98) and (99) by mere substitution.

Polymer BKL map
In this last section, we calculate the Polymer modified BKL map and study some of its properties.
In Sect. 7.1 the Polymer BKL map on the Kasner indices is derived while in Sect. 7.2 some noteworthy properties of the map are discussed.
Finally in Sect. 7.3 the results of a simple numerical simulation of the Polymer BKL map over many iterations are presented and discussed.

Polymer BKL map on the Kasner indices
Here we use the method outlined in Sects. 6.1 to 6.3 to directly calculate the Polymer BKL map on the Kasner indices.
First, we look at the asymptotic limit at ±∞ for the solution of Polymer Bianchi II (98). As in the standard case, the Polymer Bianchi II model "links together" two Kasner epochs at plus and minus infinity. In this sense, the dynamics of Polymer Bianchi II is not qualitatively different from the standard one. We will tell the quantities at plus and minus infinity apart by adding a prime to the quantities at minus infinity and leaving the quantities at plus infinity un-primed. The two Kasner solutions at plus and minus infinity can be still parametrized according to (23).
By summing together equations (23) and using the first Kasner constraint (24), we find that: By taking the limits at τ → ±∞ of the derivatives of the Polymer Bianchi II solution (98), and comparing (106) and (107) with (23) and (105), we can find the following expressions for the primed and unprimed Kasner indices and Λ: where the functions f l,m,n,Λ and f l,m,n,Λ correspond to the r.h.s. of (106) and (107) respectively and c is a shorthand for the set {c 1 , . . . , c 12 }. Now, in complete analogy with (35), we look for an asymptotic condition. We don't need to solve the whole system (108), but only to find a relation between the old and new Kasner indices and Λ at the first order in μ 2 . In practice, not all of the (108) relations are actually needed to find an asymptotic condition.
We recall from (36) that the standard BKL map on Λ is where we have assumed that p l < 0, as we will continue to do in the following. As everything until now is hinting to, we prescribe that the Polymer modified BKL map reduces to the standard BKL map when μ → 0, so that relation (109) is modified in the Polymer case only perturbatively. Since we are considering only the first order in μ 2 , we can write for the Polymer case: Solving for h(c), we find that We can invert the subsystem made up by the zeroth order of equations (108e), (108f) and (108h) to express h = h( p l , p m , Λ). Finally the asymptotic condition reads as: We stress that this is only one out of many equivalent ways to extract from the system (108) an asymptotic condition. Now, we need other three conditions to derive the Polymer BKL map. One is provided by the sum of the primed Kasner indices at minus infinity. This is the very same both in the standard case (24) and in the Polymer case.
Sadly, the two conditions (33) and (34) are not valid anymore. Instead, one condition can be derived by noticing that Lastly, we choose as the fourth condition the sum of the squares of the Kasner indices at minus infinity (93). Because of the assumption Q 1, we will consider here only the connected branch of (93). We gather now all the four conditions and put them in a system: where we have also used definition (94). Finding the polymer BKL map is now only a matter of solving the system (112) for the un-primed indices. The Polymer BKL map at the first order in μ 2 is then: p n ≈ 2 p l + p n 1 + 2 p l + 2 9 Q 2 p l × −3 + p l + 9 p l 2 + 8 p l 3 + p n 6 + 11 p l + 7 p l where all the terms of order O(Q 4 ) were neglected. It is worth noticing that, by taking the limit μ → 0, the standard BKL map (36) is immediately recovered. We stress that the form of the map is not unique. One can use the two Kasner constraints (24) and (93) to "rearrange" the Kasner indices as needed.

BKL map properties
First we derive how the Polymer BKL map can be expressed in terms of the u and Q parameters: Then we discuss some of its noteworthy properties.
We start by inserting in the polymer BKL map (113) the connected branch of the (u, Q) parametrization (96) and requiring the relations to be always satisfied. The resulting Polymer BKL map on (u, Q) is quite complex. We therefore split it in many terms: where we recall that the dominions of definition for u and Q are u ∈ [1, ∞) and Q ∈ [0, 1]. All the square-roots and denominators appearing in (116) are well behaved (always with positive argument or non zero respectively) for any (u, Q) in the intervals of definition. It is not clearly evident at first sight, but the polymer BKL map reduces to the standard one (39) if Q → 0. Now, we study the asymptotic behavior of the polymer BKL map (116a) for u going to infinity. This limit is relevant to us for two reasons. Firstly, In Sect. 5.1 we already pointed out that infinitely many bounces of the β-point against the Fig. 7 The Polymer BKL map for u (116a) is plotted for Q = 1/10. There a plateau u ≈ 1200 is reached at about u ≈ 10 4 . This means that for any u 10 4 the next u will inevitably be u ≈ 1200, no matter how big u is potential walls happen until the singularity is reached. This means that the Polymer BKL map is to be iterated infinitely many times.
Secondly, we suppose the Polymer BKL map [or at least its u portion (116a)] to be ergotic, so that every open set of the parameter space is visited with non null probability. This assumption is backed up by the observation that the Polymer BKL map (116a) tends asymptotically to the standard BKL map (39) (that is ergodic), as shown in the following Sect. 7.3 through a numerical simulation.
Hence, every open interval of the u > 1 line is to be visited eventually: we want therefore to assess how the map behaves for u → ∞. Physically, a big value for u means a long (in term of epochs) Kasner era, and in turn this means that the Universe is going deep inside one of the corners of Fig. 2. As we can appreciate from Fig. 7, for u → ∞, u reaches a plateau: The plateau (117) is higher and steeper the more Q is close to zero. This physically means that there is a sort of "centripetal potential" that is driving the β-point off the corners and towards the center of the triangle. We infer that this "centripetal potential" is somehow linked to the velocitylike quantity v(τ ) that appears in the polymer Bianchi II Einstein's Eq. (97). Because of this plateau, the mechanism that drives the Universe away from the corner, implicit in the standard BKL map (39), seems to be much more efficient in the quantum case. For Q → 0, the plateau tends to disappear: limu→∞ Q→0 u (u, Q) = ∞. As a side note, we remember that it has not been proved analytically that the standard BKL map is still valid deep Fig. 8 For any Q ∈ [0, 1] it is shown the corresponding u * (Q) such that for u < u * the Polymer BKL map on Q (116b) is monotonically decreasing. For a single iteraction it can be shown that for every Q < 0.96 the following Q < 1. The function u * (Q) cannot be given in closed form inside the corners. As a matter of fact, the BKL map can be derived analytically only for the very center of the edges of the triangle of Fig. 2, because those are the only points where the Bianchi IX potential is exactly equal to the Bianchi II potential. The farther we depart from the center of the edges, the more the map looses precision. At any rate, there are some numerical studies for the Bianchi IX model [68,69] that show how the standard BKL map is valid with good approximation even inside the corners.
The analysis of the behavior of the Polymer BKL map for the Q parameter (116b) is more convoluted. Little can be said analytically about the overall behavior across multiple iterations because of its evident complexity. For this reason, in Sect. 7.3 we discuss the results of a simple numerical simulation that probes the behavior of the map (116) over many iterations.
The most important point to check is if the dominion of definition (95) for Q is preserved by the map. We remember that for Q ≈ 1 the perturbation theory, by the means of which we have derived the map, is not valid anymore. So every result in that range is to be taken with a grain of salt.
In Fig. 8 it is shown the maximum value of u * = u * (Q) for which the polymer BKL map on Q is monotonically decreasing, i.e. for any Q ∈ [0, 1] it is shown the corresponding u * , such that for u < u * the polymer BKL map on Q is monotonically decreasing. It is clear how the more Q is little the more u can grow before coming to the point that Q stops decreasing and starts increasing.
That said, as soon as we consider the behavior of u, and particularly the plateau of Fig. 7, we notice that the Universe cannot "indulge" much time in "big-u regions": even if it happens to assume a huge value for u, this is immediately dampened to a much smaller value at the next iteration. The net result is that Q is almost always decreasing as it is also strongly suggested by the numerical simulation of Sect. 7.3.
Summarizing, the Polymer BKL map on Q (116b), apart from a small set of initial conditions in the region where the perturbation theory is failing Q ≈ 1, preserves the dominion of definition of Q (95) and is decreasing at almost any iteration.

Numerical simulation
In this section we present the results of a simple numerical simulation that unveils some interesting features of the Polymer BKL map (116). The main points of this simulation are very simple: • An initial couple of values for (u, Q) is chosen inside the dominion u ∈ [1, ∞) and Q ∈ [0, 0.96]. • We remember from Sects. 3.3 and 6.2 that, at the end of each Kasner era, the u parameter becomes smaller than 1 and needs to be remapped to values greater than 1. This marks the beginning of a new Kasner era. This remapping is performed through the relations listed on page 15. In the standard case, because the standard BKL map is just u → u − 1, the u parameter cannot, for any reason, become smaller than 0. For the Polymer Bianchi IX, however, it is possible for u to become less than zero. This is why we have derived many remapping relations, to cover the whole real line. • Many values (≈ 2 18 ) for the initial conditions, randomly chosen in the interval of definition, were tested. We didn't observe any "anomalous behavior" in any element of the sample. This meaning that all points converged asymptotically to the standard BKL map as is discussed in the following.
One sample of the simulation is displayed in Fig. 9 using a logarithmic scale for Q to show how effective is the map in damping high values of Q.
From the results of the numerical simulations the following conclusions can be drawn: • The polymer BKL map (116) is "well behaved" for any tested initial condition u ∈ [1, ∞) and Q ∈ [0, 0.96]: the dominion of definition for u and Q are preserved. • The line Q = 0 is an attractor for the Polymer BKL map: for any tested initial condition, the Polymer BKL map eventually evolved until becoming arbitrarily close to the standard BKL map. • The Polymer BKL map on Q is almost everywhere decreasing. It can happen that, especially for initial values Q ≈ 1, for very few iterations the map on Q is increasing, but the overall behavior is almost always decreasing. The probability to have an increasing behavior of Q gets smaller at every iteration. In the limit of infinitely many iterations this probability goes to zero.   (39), we expect that the notion of chaos for the standard BKL map, given for example in [56], can be applied, with little or no modifications, to the Polymer case, too (although this has not been proven rigorously).

Physical considerations
In this section we address two basic questions concerning the physical link between Polymer and Loop quantization methods and what happens when all the Minisuperspace variables are discrete, respectively. First of all, we observe that Polymer Quantum Mechanics is an independent approach from Loop Quantum Gravity. Using the Polymer procedure is equivalent to implement a sort of discretization of the considered configurational variables. Each variable is treated separately, by introducing a suitable graph (de facto a one-dimensional lattice structure): the group of spatial translations on this graph is a U (1) type, therefore the natural group of symmetry underlying such quantization method is U (1) too, differently from Loop Quantum Gravity, where the basic group of symmetry is SU (2).
It is important to stress that Polymer Quantum Mechanics is not unitary connected to the standard Quantum Mechanics, since the Stone -Von Neumann theorem is broken in the discretized representation. Even more subtle is that Polymer quantization procedures applied to different configurational representations of the same system are, in general, not unitary related. This is clear in the zero order WKB limit of the Poly-mer quantum dynamics, where the formulations in two sets of variables, which are canonically related in the standard Hamiltonian representation, are no longer canonically connected in the Polymer scenario, mainly due to the non-trivial implications of the prescription for the momentum operator.
When the Loop Quantum Gravity [18,70] is applied to the Primordial Universe, due to the homogeneity constraint underlying the Minisuperspace structure, it loses the morphology of a SU (2) gauge theory (this point is widely discussed in [71,72]) and the construction of a kinematical Hilbert space, as well as of the geometrical operators, is performed by an effective, although rigorous, procedure. A discrete scale is introduced in the Holonomy definition, taken on a square of given size and then the curvature term, associated to the Ashtekar-Barbero-Immirzi connection, (the so-called Euclidean term of the scalar constraint) is evaluated on such a square path. It seems that just in this step the Loop Quantum Cosmology acquires the features of a polymer graph, associated to an underlying U (1) symmetry. The real correspondence between the two approaches emerges in the semiclassical dynamics of the Loop procedure [46], which is isomorphic to the zero order WKB limit of the Polymer quantum approach. In this sense, the Loop Quantum Cosmology studies legitimate the implementation of the Polymer formalism to the cosmological Minisuperspace.
However, if on one hand, the Polymer quantum cosmology predictions are, to some extent, contained in the Loop Cosmology, on the other hand, the former is more general because it is applicable to a generic configurational representation, while the latter refers specifically to the Ashtekar-Berbero-Immirzi connection variable.
Thus, the subtle question arises about which is the proper set of variables in order that the implementation of the Polymer procedure mimics the Loop treatment, as well as which is the physical meaning of different Polymer dynamical behaviours in different sets of variables. In [73] it is argued, for the isotropic Universe quantization, that the searched correspondence holds only if the cubed scale factor is adopted as Polymer variable: in fact this choice leads to a critical density of the Universe which is independent of the scale factor and a direct link between the Polymer discretization step and the Immirzi parameter is found. This result assumes that the Polymer parameter is maintained independent of the scale factor, otherwise the correspondence above seems always possible. In this respect, different choices of the configuration variables when Polymer quantizing a cosmological system could be mapped into each other by suitable re-definition of the discretization step as a function of the variables themselves. Here we apply the Polymer procedure to the Misner isotropic variable and not to the cubed scale factor, so that different issues with respect to Loop Quantum Gravity can naturally emerge. The merit of the present choice is that we discretize the volume of the Universe, without preventing its vanishing behavior. This can be regarded as an effective procedure to include the zero-volume eigenvalue in the system dynamics, like it can happen in Loop Quantum Gravity, but it is no longer evident in its cosmological implementation.
Thus, no real contradiction exists between the present study and the Big-Bounce prediction of the Loop formulation, since they are expectedly two independent physical representations of the same quantum system. As discussed on the semi-classical level in [74], when using the cubed scale factor as isotropic dynamics, the Mixmaster model becomes non-singular and chaos free, just as predicted in the Loop Quantum Cosmology analysis presented in [27]. However, in such a representation, the vanishing behavior of the Universe volume is somewhat prevented a priori, by means of the Polymer discretization procedure.
Finally, we observe that, while in the present study, the Polymer dynamics of the isotropic variable α preserves (if not even enforces) the Mixmaster chaos, in [30] the Polymer analysis for the anisotropic variables β ± is associated to the chaos disappearance. This feature is not surprising since the two approaches have a different physical meaning: the discretization of α has to do with geometrical properties of the space-time (it can be thought as an embedding variable), while the implementation of the Polymer method to β ± really affects the gravitational degrees of freedom of the considered cosmological model.
Nonetheless, it becomes now interesting to understand what happens to the Mixmaster chaotic features when the two sets of variables (α and β ± ) are simultaneously Polymer quantized. In the following subsection, we provide an answer to such an intriguing question, at least on the base of the semi-classical dynamics.

The polymer approach applied to the whole Minisuperspace
According to the polymer prescription (55), the super-Hamiltonian constraint is now: where μ α is the polymer parameter associated to α and μ the one associated to β ± . The dynamics of the system can be derived through the Hamilton's equations (44). As usual, we start considering the simple Bianchi I case V B = 0. Following the procedure of Sect. 6.1, we find that the Universe is still singular. In fact, the solution to Hamilton's equations in the (88) time gauge (with μ α instead of μ) is: sin(2μp ± ) sin(2μ α p α ) ln(t) (119b) from which it follows that α(t) → −∞ as t → 0. Moreover, the sum of the squared Kasner indices, calculated by making a reverse canonical transformation from Misner to ordinary variables, reads as: from which is possible to derive the anisotropy velocity, according to Eq. (92). The anisotropy velocity can be also calculated through the ADM reduction method, as in Sect.
that turns out to be greater than one when μ α ≥ μ (see Fig. 10). Finally, it should be noted that since in the general V B = 0 case the wall velocity β wall is not influenced by the introduction of the polymer quantization, from Eq. (123) we find that the maximum incident angle for having a bounce against a given potential wall θ max = arccos β wall β is, also in this case, always greater than π/3 when μ α ≥ μ. We can then conclude that also when the polymer approach is applied to all the three configuration variables simultaneously, the Universe can be singular and chaotic just such as the one analyzed above.

Conclusions
In the present study, we analyzed the Mixmaster model in the framework of the semiclassical Polymer Quantum Mechanics, implemented on the isotropic Misner variable, according to the idea that the cut-off physics mainly concerns the Universe volume.
We developed a semiclassical and quantum dynamics, in order to properly characterize the structure of the singularity, still present in the model. The presence of the singularity is essentially due to the character of the isotropic variable conjugate momentum as a constant of motion.
On the semiclassical level we studied the system evolution both in the Hamiltonian and field equations representation, generalizing the two original analyses in [7,8], respectively. The two approaches are converging and complementary, describing the initial singularity of the Mixmaster model as reached by a chaotic dynamics, that is, in principle, more complex than the General Relativity one, but actually coincides with it in the asymptotic limit.
This issue is a notable feature, since in Loop Quantum Cosmology is expected that the Bianchi IX model is chaosfree [26,27] and it is well known [28] that the Polymer semiclassical dynamics closely resembles the average feature of a Loop treatment in the Minisuperspace. However, we stress that, while the existence of the singularity in Polymer Quantum Mechanics appears to be a feature depending on the nature of the adopted configuration variables, nonetheless the properties of the Poincaré map of the model is expected to be a solid physical issue, independent on the particular representation adopted for the system.
The canonical quantization of the Mixmaster Universe that we performed in the full Polymer quantum approach, i.e. writing down a Wheeler-DeWitt equation in the momentum representation, accordingly to the so-called Continuum Limit discussed in [29], is completely consistent with the semiclassical results. In fact, the Misner demonstration, for the standard canonical approach, that states with high occupation numbers can survive to the initial singularity, remains still valid in the Polymer formulation, here presented. This issue confirms that the cut-off we introduced in the configuration space on the isotropic Misner variable does not affect the cosmological character of the Mixmaster model. This result appears rather different from the analysis in [30], where the polymer approach has been addressed for the anisotropic Misner variables (the real degrees of freedom of the cosmological gravitational field) with the emergence of a non-chaotic cosmology. Such a discrepancy suggests that the Polymer regularization of the asymptotic evolution to the singularity produces more profound modifications when it touches physical properties than geometrical features. Actually, the isotropic Misner variable can be suitably interpreted as a good time variable for the system (an embedding variable in the language of [75,76]), while the Universe anisotropies provide a precise physical information on the considered cosmological model.
Despite this consideration about the gauge-like nature of the Misner isotropic variable, which shed light on the physics of our dynamical results, nonetheless we regard as important to perform further investigations on the nature of the singularity when other variables are considered to characterize the Universe volume, since we expect that, for some specific choice the regularization of the Big-Bang to the Big-Bounce must take place (see for instance [74]). However, even on the basis of the present analysis, we suggest that the features of the Poincaré map of the Bianchi IX model and then of the generic cosmological singularity (locally mimicked by the same Bianchi IX-like time evolution) is a very general and robust property of the primordial Universe, not necessarily connected with the existence of a cut-off physics on the singularity.