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

We analyze the semiclassical 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 surprising result of a still singular and chaotic cosmology, whose Poincar\'e 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.


I. 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]. 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 Universe 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 cut-off 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].
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 for → 0 is taken. Such a semiclassical approach clearly offers a characterization for the behavior of the mean values of the quantum Uni-verse, in the spirit of the Ehrenfest theorem [37]. This study surprisingly 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 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.
The paper is organized as follows. In Section II 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 Section III 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 Section IV, where we review the semiclassical and quantum properties of the model as it was studied by Misner in [8].
Section V 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.
Finally, in Section VI the semiclassical behavior of the polymer modified Bianchi I and Bianchi II models is de-veloped 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 Section VII.

II. 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 well-defined, 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 Sec. II D.
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 [38][39][40], so that the cutoff introduced in PQM is assumed to be a fundamental quantity. Some results of Loop Quantum Gravity [41] (the discrete spectrum of the area operator) and String Theory [42] (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,43,44] who also credit a previous work of Varadarajan [45] for some ideas. It was then further refined also by Corichi et al. [29,46]. They developed the PQM in the expectation to shed some light on the connection between the Planckian-energy Physics of Loop Quantum Gravity and the lower-energy Physics of the Standard Model.

A. 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 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 [47]. At this stage we have not made any assumption regarding the Hilbert space of the theory, yet.
To introduce the Polymer representation in Sec. II B, 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 Lebesgue-integrable functions H S = L 2 (R, dq) and the action of the bases operators on it: where ψ ∈ L 2 (R, dq) and q ∈ R.

B. 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 |µ asˆ |µ := µ |µ ;ŝ(λ) |µ := |µ + λ 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 λ V (λ) · ψ µ (p) = e iλp e iµp = e i(λ+µ)p = ψ (λ+µ) (p) (9) 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 byq 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 [48] 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 [49]: where the Haar measure dµ H is the natural choice.
In terms of the fundamental functions ψ µ (p) the inner product takes the form

C. 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 Sec. II B), 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 operatorV (λ) = e iλp can be identified with the displacement oneŝ (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 . From this result it is easy to derive the approximated kinetic term p 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)).

D. 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: γ 0 → γ n = q k ∈ R | q k = ka n , with a n = a 0 2 n , ∀k ∈ Z (17) 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 : m ψ(ma n )δ man,q ∈ H poly → m ψ(ma n )χ αm (q) ∈ H S (18) 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Ĥ Cn 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.

III. THE HOMOGENEOUS MIXMASTER MODEL: CLASSICAL DYNAMICS
With the aim of investigating in Sec. VII the modifications to the Bianchi IX dynamics produced by the introduction of the Polymer cutoff, in this Section we briefly review some relevant results obtained by Belinskii, Khalatnikov and Lifshitz (BKL) regarding the classical Bianchi IX model. In particular they found that an approximate solution for the Bianchi IX dynamics can be given in the form of a Poincaré recursive map called BKL map. This map is obtained as soon as the exact solutions of Bianchi I and Bianchi II are known. Hence, we will first derive the dynamical solution to the classical Bianchi I and II models.
The Einstein equations of a generic Bianchi model can be expressed in terms of the spatial metric η(t) = diag[a(t), b(t), c(t)] 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]. 20) 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.
A. Bianchi I The Bianchi I model in void is also called Kasner solution [50]. By substituting in Eqs (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 Section VI A. 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 equations (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 2p l , t 2pm , t 2pn ) has a physical naked singularity, called also Big Bang, when t = 0 (or τ → −∞).

B. 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 Ilike, 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.
where we have exploited the notation of equation (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 f i (p l , p m , p n , Λ, p l , p m , p n , Λ ) = 0 with i = 1, . . . , 4. 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.

C. 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: Eq. (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 ldirection. 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 is the integer part and x 0 = u 0 − [u 0 ] is the fractional part. If 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 [52][53][54][55][56][57].

IV. HAMILTONIAN FORMULATION OF THE MIXMASTER MODEL
In this section we summarize some important results obtained applying the R. Arnowitt, S. Deser and C.W. 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. β ij 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 β ij = 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 C.W. Misner in [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 [58], 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: (42) 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 : 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 wall V B 0), while V B (β ± ) is nonnegligible 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 Eqs. (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 Sec. III A), we recover the sequence of Kasner epochs that characterizes the BKL map (Sec. III C).
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.

A. 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 equa-tionĤ 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 C.W. 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.

V. MIXMASTER MODEL IN THE POLYMER APPROACH
Here the Hamiltonian formalism introduced in Sec. IV 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 Sec. II C: which we apply to the super-Hamiltonian constraint H B = 0 (with H B defined in (41)): 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.

A. 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 Eqs. (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 1 In this case the word "semiclassical" means that our super-Hamiltonian constraint was obtained as the lowest order term of a WKB expansion for → 0.
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)): 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 surprising 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 [59,60]. In particular, when dealing with PQM, a spatial lattice of spacing µ must be introduced to regularize the theory (Sec II D). 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 [59] and Bianchi I [60] 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) 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): 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 .

B. Comparison with previous models
As we stressed in Sec. V A our model turns out to be very different from Loop Quantum Cosmology [59,60], which predicts a Big Bounce-like dynamics. Similarly if the polymer approach is applied to the anisotropies β ± , leaving the volume variable α classical [30], a nonsingular and 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 [61]: 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 [62] 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.

C. 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 φ(p ± ) = φ + (p + )φ − (p − ) with: 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.

D. Adiabatic Invariant
Because of the various hints about the presence of a chaotic behavior, outlined in Sec. V A and that will be confirmed by the analysis of the Poincaré map in Sec VII, 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 Sec. V A 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: [distance covered before the bounce] = r αi |α i | (82a) [distance covered after the bounce] = r α f |α f | (82b) Now if we apply the Sine theorem to the triangles A BC and A BD, 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 BD) 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 quasi-classical state of the Universe is allowed up to the initial singularity.

VI. SEMICLASSICAL SOLUTION FOR THE BIANCHI I AND II MODELS
As we have already outlined in Sec. III, 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 Sec. VI A we solve the Hamilton's Eqs in Misner variables for the Polymer Bianchi I, in the semiclassical approximation introduced in Sec. V A.
Then in Sec. VI B 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 Sec. VI C 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.

A. Polymer Bianchi I
As already shown in Sec. V A, 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 Sec. VI B, 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 disconnectedbranch.
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 Sec. V A. 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. 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 Sec. V A. 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 Eqs (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.

B. 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 (92) as well, constraint (92) 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 (93).
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 Sec. VII B 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 [63,64].
The (u, Q)-parametrization (95) 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 (95) 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 (94) are not physical and should never be considered. In Figure 6 the values of the ordered Kasner indices are displayed for the (u, Q)-parametrization (95), 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.

C. Polymer Bianchi II
Here we apply the method described in VI A to find an approximate solution to the Einstein's Eqs of the Polymer Bianchi II model. We start by selecting the V B potential appropriate for Bianchi II (60) and we substitute The connected-branch is the one containing the red (dark gray) lines, i.e. the standard parametrization shown in Fig. 1. The Q parameter of the connected-branch is increasing from 0 to 1 in the positive direction of the Q-axis, while for the disconnected-branch it is decreasing from 1 to 0. The black line delimits the part of the (u, Q)-plane where there is only one (above) or two (below) negative Kasner indices. The plot is cut on the bottom and on top but it is actually extending up to ±∞. it in the Hamiltonian (56) (in the following we will always assume the time gauge N = 1).
Then, starting from Hamilton's Eqs, inverting the Misner canonical transformation and converting the synchronous time t-derivatives into logarithmic time τderivatives, the Polymer Bianchi II Einstein's Eqs 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 Eqs (96) reduce to the correct classical Eqs (29). Now we find an approximate solution to the system (96), 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 Eqs (96), the perturbative expansion will only contain even powers of µ.
First we expand the solution at the first order in µ 2 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 (96), 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 (98) in (96) 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: for brevity and we have substituted the zeroth order solution (31) where needed.
Eqs (99) are two almost uncoupled non-homogeneous ODEs. (99a) is a linear second order non-homogeneous ODE with non-constant coefficients and, being completely uncoupled, it can be solved straightly. (99b) is solved by mere substitution of the solution of (99a) in it and subsequent double integration.
To solve (99a) we exploited a standard method that can be found, for example, in [65,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 (99a) is whose general solution reads as: By applying the above-mentioned method, we obtain the following solution for (99a) The solution for (99b) is found by substituting (102) 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 (97) and (98) by mere substitution.

VII. POLYMER BKL MAP
In this last Section, we calculate the Polymer modified BKL map and study some of its properties.
In Sec. VII A the Polymer BKL map on the Kasner indices is derived while in Sec. VII B some noteworthy properties of the map are discussed.
Finally in Sec VII C the results of a simple numerical simulation of the Polymer BKL map over many iterations are presented and discussed.

A. Polymer BKL map on the Kasner indices
Here we use the method outlined in Secs. from VI A to VI C 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 (97). 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 (97), and comparing (105) and (106) with (23) and (104), 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 (105) and (106) 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 (107), 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 (107) 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 (108) 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 (107e), (107f) and (107h) 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 (107) 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 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 (92) to "rearrange" the Kasner indices as needed.

B. 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 (112) the connected branch of the (u, Q) parametrization (95) and requiring the relations (114) to be always satisfied. The resulting Polymer BKL map on (u, Q) is quite complex. We therefore split it in many terms: +2398u 8 − 1350u 9 + 663u 10 − 162u 11 + 63u 12 F (u) = 72(−1 + u + 3u 3 + 3u 5 + u 6 + 2u 7 ) 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 (115) 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 (115a) for u going to infinity. This limit is relevant to us for two reasons. Firstly, In Sec. V A we already pointed out that infinitely many bounces of the β-point against the 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 (115a)) 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 (115a) tends asymptotically to the standard BKL map (39) (that is ergodic), as shown in the following Sec. VII C 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: lim u→∞ u (u, Q) = 24 + Q 2 + 3 (7Q 4 + 48Q 2 + 192) 4Q 2 (116) The plateau (116) 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 velocity-like quantity v(τ ) that appears in the polymer Bianchi II Einstein's Eqs (96). Because of this plateau, the mechanism that drives the Universe away from the corner, implicit in the standard BKL map (39), seems to  1] it is shown the corresponding u * (Q) such that for u < u * the Polymer BKL map on Q (115b) 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.
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 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 [66,67] 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 (115b) is more convoluted. Little can be said analytically about the overall behavior across multiple iterations because of its evident complexity. For this reason, in Sec. VII C we discuss the results of a simple numerical simulation that probes the behavior of the map (115) over many iterations.
The most important point to check is if the dominion of definition (94) 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 de-creasing as it is also strongly suggested by the numerical simulation of Sec. VII C.
Summarizing, the Polymer BKL map on Q (115b), 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 (94) and is decreasing at almost any iteration.

C. Numerical simulation
In this Section we present the results of a simple numerical simulation that unveils some interesting features of the Polymer BKL map (115). 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 Sec.III C and VI B 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 14. 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 (115) 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  (115) is just to take a single point and map it to the next one. A total of 14 Kasner eras are shown, each with a different color. The fact that, starting from the upmost point, the map on Q is almost always decreasing is very evident from the downwards-going ladder-like shape of the eras.
iteration. In the limit of infinitely many iterations this probability goes to zero.
• Since the Polymer BKL map (115) tends asymptotically to the standard BKL map (39), we expect that the notion of chaos for the standard BKL map, given for example in [53], can be applied, with little or no modifications, to the Polymer case, too (although this has not been proven rigorously).

VIII. 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 [8] and [7], 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 surprising feature, since in Loop Quantum Cosmology is expected that the Bianchi IX model is chaos-free [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 chaoticity 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 socalled 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 and it suggests that the chaoticity of the model is an intrinsic feature, rather robust in the asymptotic limit to the singularity. This result appears rather different from the analysis in [30], where the polymer approach has been addressed for the anisotropic Misner variables (the real gravitational 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 interpret as a good time variable for the system (an embedding variable in the language of [68,69]), while the Universe anisotropies provides 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. However, even on the basis of the present analysis, we suggest that the chaotic feature 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.