We proceed by discretizing the areas *pi* and imposing in (60) the polymer substitution for the connections *ci* :

$$c\_i \to \frac{1}{b\_i} \sin(b\_i c\_i) \,\tag{110}$$

where *bi* (*i* = 1, 2, 3) refer to the three independent polymer lattices and {*ci* , *pj*} = *γδij*. Then, the polymer Hamiltonian takes the following form:

$$\mathcal{L}\_{\text{BI}}^{\text{poly}} = -\frac{1}{\gamma^2 V} \sum\_{i \neq j} \frac{\sin(b\_i c\_i) p\_i \sin(b\_j c\_j) p\_j}{b\_i b\_j} + \frac{P\_\phi^2}{2V} = 0,\tag{111}$$

where *i*, *j* = 1, 2, 3 and *V* = √*p*1 *p*2 *p*3.

We derive the dynamics after choosing *φ* as the internal time, i.e., imposing the gauge *N* = √*p*1 *p*2 *p*3 *Pφ* , so the equations of motion take the following form:

$$\frac{dp\_i}{d\phi} = -\frac{p\_i \cos(b\_i c\_i)}{\gamma P\_\phi} \left[ \frac{p\_j}{b\_j} \sin(b\_j c\_j) + \frac{p\_k}{b\_k} \sin(b\_k c\_k) \right] \tag{112a}$$

$$\frac{dc\_i}{d\phi} = \frac{\sin(b\_i c\_i)}{\gamma b\_i P\_\phi} \left[ \frac{p\_j}{b\_j} \sin(b\_j c\_j) + \frac{p\_k}{b\_k} \sin(b\_k c\_k) \right] \tag{112b}$$

for *i*, *j*, *k* = 1, 2, 3, *i* 6= *j* 6= *k*. It is possible to solve this system by establishing the initial conditions on the variables (*ci* , *pi*), which also satisfy the Hamiltonian constraint (111). In this respect, we make the following choice:

$$\begin{aligned} c\_1(0) &= c\_2(0) = c\_3(0) = \frac{\pi}{2b\_0}, \\ p\_1(0) &= \vec{p}\_1, \quad p\_2(0) = \vec{p}\_2, \\ p\_3(0) &= \vec{p}\_3 = \frac{P\_\phi^2 \gamma^2 b\_0 - 2\vec{p}\_1 \vec{p}\_2}{2(\vec{p}\_1 + \vec{p}\_2)}, \end{aligned} \tag{113}$$

where *b*1 = *b*2 = *b*3 = *b*0 without loss of generality. Moreover, it can be easily seen that the momentum conjugate to the scalar field is a first integral since the variable *φ* is cyclic in (111). Other constants of motion can be obtained by combining the Hamilton Equation (112), so we obtain the following:

$$\frac{p\_i \sin(b\_0 c\_i)}{b\_0} = \mathcal{K}\_{i\prime} \qquad P\_\Phi = \mathcal{K}\_\Phi \, , \tag{114}$$

where the considered values of K*φ* and K*i* depend on the initial conditions. Identifying these first integrals allows to transform the six-equations system shown in (112) in three closed systems along the three spatial directions as follows:

$$\frac{dp\_i}{d\phi} = -\frac{p\_i \cos(b\_0 c\_i)}{\gamma P\_\phi} \left[\mathcal{K}\_j + \mathcal{K}\_k\right],\tag{115a}$$

$$\frac{dc\_i}{d\phi} = \frac{\sin(b\_0 c\_i)}{\gamma b\_0 P\_\phi} \left[\mathcal{K}\_j + \mathcal{K}\_k\right],\tag{115b}$$

where *i* 6= *j* 6= *k*.

Thanks to this procedure, the equations of motion can be solved analytically, leading to the following solutions:

$$\begin{split} p\_1(\phi) &= \bar{p}\_1 \cosh\left[\frac{(\gamma^2 P\_\phi^2 b\_0^2 + 2\bar{p}\_2^2)\phi}{2\gamma P\_\phi b\_0(\bar{p}\_1 + \bar{p}\_2)}\right], \\ p\_2(\phi) &= \bar{p}\_2 \cosh\left[\frac{(\gamma^2 P\_\phi^2 b\_0^2 + 2\bar{p}\_1^2)\phi}{2\gamma P\_\phi b\_0(\bar{p}\_1 + \bar{p}\_2)}\right], \\ p\_3(\phi) &= \bar{p}\_3 \cosh\left[\frac{(\bar{p}\_1 + \bar{p}\_2)\phi}{\gamma P\_\phi b\_0}\right]. \end{split} \tag{116}$$

We have reported only the explicit expression of the variables *pi* as functions of *φ* since we are interested in the Universe volume behavior *V*(*φ*) = p *p*1(*φ*)*p*2(*φ*)*p*3(*φ*) that is shown in Figure 7. The resulting trajectory highlights that a quantum Big Bounce replaces the classical Big Bang thanks to the polymer effects, which are expected to become dominant near the Planckian region.

**Figure 7.** Polymer trajectory of the Universe volume *V* = √*p*1 *p*2 *p*3 as function of *φ*: the Big Bounce replaces the classical singularity of the Bianchi I model.

Now, we focus our attention on the analysis of the critical energy density of matter at the Bounce:

$$\rho\_{\rm crit}^{\phi} = \frac{P\_{\phi}^{2}}{2V(\phi)^{2}}\bigg|\_{\phi\_{\rm B}} = \frac{P\_{\phi}^{2}}{2\vec{p}\_{1}\vec{p}\_{2}\vec{p}\_{3}} = \frac{P\_{\phi}^{2}(\vec{p}\_{1} + \vec{p}\_{2})}{\vec{p}\_{1}\vec{p}\_{2}(\gamma^{2}P\_{\phi}^{2}b\_{0}^{2} - 2\vec{p}\_{1}\vec{p}\_{2})},\tag{117}$$

i.e., the energy density of the matter scalar field in correspondence of the minimum Universe volume at the time of the Bounce *φB* = 0. As we can see from (117), *ρ φ* crit clearly depends on the initial conditions on the motion. Moreover, in the simplest case K1 = K2 = K3 it reduces to the following:

$$
\rho\_{\rm crit}^{\Phi} = \left(\frac{3}{\gamma^2 b\_0^2}\right)^{\frac{3}{2}} \frac{\sqrt{2}}{P\_{\Phi}}.\tag{118}
$$

reproducing a consistent behavior with that one obtained in the FLRW model in the same variables; see (100). Indeed, it is possible to see that for *Pφ* 1 the critical matter energy density increases until it diverges, while on the other hand, it approaches zero when *Pφ* 1, highlighting that in the Ashtekar formulation of the Bianchi I model, the Big Bounce has no universal features.

Let us now study this model on a quantum level by implementing the Dirac quantization. In the momentum representation of the polymer formulation, the fundamental operators act as follows:

$$\mathfrak{P}\_{i} := -i\gamma \frac{d}{d\mathfrak{c}\_{i}} , \quad \mathfrak{E}\_{i} := \frac{\sin(b\_{i}c\_{i})}{b\_{i}} , \quad \mathfrak{P}\_{\Phi} := -i\frac{d}{d\phi} . \tag{119}$$

Before quantization, we rewrite the WDW equation in the form of a Schrödinger one in the attempt of defining a positive and conserved probability density. Therefore, we recall the semiclassical scalar constraint and perform an ADM reduction of the variational principle:

$$
\epsilon\_{\phi}^{2} - \Theta = 0,\tag{120}
$$

where

$$\Theta = \frac{2}{\gamma^2} \left[ \frac{\sin(b\_1 c\_1) p\_1 \sin(b\_2 c\_2) p\_2}{b\_1 b\_2} + \frac{\sin(b\_1 c\_1) p\_1 \sin(b\_3 c\_3) p\_3}{b\_1 b\_3} + \frac{\sin(b\_2 c\_2) p\_2 \sin(b\_3 c\_3) p\_3}{b\_2 b\_3} \right]. \tag{121}$$

*P*

After choosing the scalar field *φ* as the temporal parameter, we derive the ADM Hamiltonian by solving the scalar constraint (120) with respect to the momentum associated to the scalar field:

$$P\_{\Phi} \equiv \mathcal{C}\_{\text{BI}}^{\text{ADM-poly}} = \sqrt{\Theta}. \tag{122}$$

where we choose the positive root in order to guarantee the positive character of the lapse function. Thanks to this procedure, the WDW equation can be rewritten in the form of a Schrödinger one by promoting the ADM Hamiltonian to a quantum operator:

$$-i\partial\_{\phi}\Psi = \sqrt{\bigoplus} \Psi,\tag{123}$$

where the operator √ Θˆ , that we assume well-defined, can be written as follows:

$$\sqrt{\Theta} = \left[\frac{2}{\gamma^2} \left(\partial\_{\mathbf{x}\_1}\partial\_{\mathbf{x}\_2} + \partial\_{\mathbf{x}\_1}\partial\_{\mathbf{x}\_3} + \partial\_{\mathbf{x}\_2}\partial\_{\mathbf{x}\_3}\right)\right]^{1/2},\tag{124}$$

where the new variables *xi* are defined from the connections *ci* as follows:

$$\chi\_i = \ln \left[ \tan \left( \frac{\mu c\_i}{2} \right) \right] + \bar{\chi}\_i. \tag{125}$$

Now, we are allowed to introduce the probability density as follows:

$$\mathcal{P}(\vec{\mathbf{x}},\phi) = \Psi^\*(\vec{\mathbf{x}},\phi)\Psi(\vec{\mathbf{x}},\phi),\tag{126}$$

where

$$\Psi(\vec{x},\phi) = \int\_{-\infty}^{\infty} dk\_1 \, dk\_2 \, dk\_3 \, A(k\_1, k\_2, k\_3) \, e^{i(k\_1\mathbf{x}\_1 + k\_2\mathbf{x}\_2 + k\_3\mathbf{x}\_3 + \sqrt{2|k\_1k\_2 + k\_1k\_3 + k\_2k\_3|} \cdot \phi)}\,,\tag{127}$$

$$A(k\_1, k\_2, k\_3) = \exp\left(\sum\_{i=1}^3 -\frac{(k\_i - \overline{k}\_i)^2}{2\sigma\_{k\_i}^2}\right). \tag{128}$$

This quantity is positive everywhere by definition, and its spatial integral remains constant through time.

Hence, we have discarded the covariant formulation of the WDW equation in favor of a Schrödinger-like one in order to define a probability density through which we analyze the quantum dynamics of the Bianchi I wave packet presented in (127). This way, we have avoided all the issues regarding the sign of the probability density that would be encountered in the interpretation of the WDW formulation as a Klein–Gordon-like theory.

In Figure 8, some different sections of the probability density P are shown at different times in order to verify how its shape and its maximum evolve. The present sections were obtained by fixing two of the three coordinates through the values that they assume in the semiclassical trajectories1 . As we can see from Figure 8, the normalized quantum distributions of *x*1, *x*2, *x*3 are shown in sequence, and their spreading behavior over time is evident. Additionally, in Figure 9, the position of the peaks of P is represented by the red dots that are fitted by means of a linear interpolation and compared with the semiclassical trajectories. We can affirm that there is a good correspondence between the quantum behavior of the wave packet and the solutions of the semiclassical dynamics since all the three slopes of the functions resulting from the fit of the red dots are consistent with the semiclassical ones with a confidence level of 3 standard deviations (we have supposed a relative error of 10%, accounting for quantum fluctuations and numerical integration errors). Moreover, we want to highlight that this quantum analysis based on the sections of P is justified by the semiclassical decoupling of the equations of motion. In conclusion, we remark that this feature of our analysis in the Ashtekar variables suggests the presence of a bouncing dynamics with non-universal properties also at a quantum level.

**Figure 8.** The normalized sections P(*xi* , *φ*) are shown in sequence for *i* = 1, 2, 3 respectively at different times. Their spreading behavior over time is evident together with the Gaussian-like shape.

**Figure 9.** The three pictures show the position of the peaks of P(*xi* , *φ*) for *i* = 1, 2, 3 respectively in function of time *φ* (red dots). The resulting fitting functions (red dashed straight lines) overlap the semiclassical trajectories (continuous lines) with a confidence level of the order of 3 standard deviations.

#### *4.3. Polymer Cosmology in the Volume-like Variables*

In this section, we describe the FLRW dynamics in the polymer framework when the polymer lattice is implemented on the Universe volume *v*. Then, we generalize the analysis by taking under consideration the Bianchi I model, which admits two different sets of volume-like variables thanks to its anisotropic structure. In particular, when the Universe volume itself is considered to be one of the configurational variables, it is possible to derive a polymer modified Friedmann equation for the Bianchi I model. We will see that the total critical energy density has universal features, as demonstrated in [37] for the isotropic case.

#### 4.3.1. The FLRW Universe in the Volume Variables

Following [44], we outline the semiclassical dynamics of the model and then we perform a full quantum analysis by analyzing the properties of the Universe wave packet. In order to apply the polymer semiclassical framework to the flat FLRW model in the volume variables, we introduce the Hamiltonian constraint as follows:

$$\mathcal{L}\_{\text{FLRW}} = -\frac{27}{4\gamma^2} v \ \vec{c}^2 + \frac{P\_\phi^2}{2v} \, , \tag{129}$$

where the phase space variables *c*˜, *v* are linked to the Ashtekar variables *c*, *p* through the following canonical transformation:

$$
v = |p|^{\frac{3}{2}} = a^3, \quad \pounds = \frac{2}{3} \frac{c}{\sqrt{|p|}} \approx \frac{d}{a} \tag{130}$$

that preserves the Poisson brackets {*c*˜, *v*} = *γ* 3 . Additionally, a massless scalar field has been added to the dynamics with the role of a relational time. Thus, the modified polymer Hamiltonian is as follows:

$$\mathcal{L}\_{\text{FLRW}}^{\text{poly}} = -\frac{27}{4\gamma^2 b\_0^2} v \sin^2(b\_0 \tilde{c}) + \frac{P\_\phi^2}{2v} = 0,\tag{131}$$

where the variable *v* is defined as discrete and so the semiclassical polymer substitution (87) *c*˜ → sin(*b*0*c*˜) *b*0 is used for the generalized coordinate *c*˜.

The equations of motion for these variables are as follows:

$$\psi = -\frac{9\text{N}}{2\gamma b\_0} v \sin(b\_0 \mathfrak{x}) \cos(b\_0 \mathfrak{x}) \,\text{.}\tag{132a}$$

$$\mathcal{E} = \frac{\mathcal{N}}{3} \left( \frac{27}{4 \gamma b\_0^2} \sin^2(b\_0 \mathcal{E}) + \gamma \frac{P\_\phi^2}{2v^2} \right). \tag{132b}$$

As expected in the polymer framework, a modified Friedmann equation appears; in particular, we can obtain its analytical expression by using (132a) and the vanishing Hamiltonian constraint (131):

$$H^2 = \left(\frac{\dot{v}}{3v}\right)^2 = \frac{\rho}{3}\left(1 - \frac{\rho}{\rho\_{\rm crit}}\right), \quad \rho\_{\rm crit} = \frac{27}{4\gamma^2 b\_0^2}.\tag{133}$$

The presence of a bouncing point for the Hubble parameter represents the mechanism through which the initial singularity is regularized. In addition, the critical energy density at which the Big Bounce occurs is fixed by the spacing *b*0 associated to the polymer lattice, giving the dynamics a universal character. This result is in strong contact with the analogous considerations made in the context of the *µ*¯ approach in the LQC formulation in Section 3.3.

Considering now the scalar field *φ* as the internal time of the dynamics, we fix the time gauge, i.e.,

$$1 = \phi = \mathbf{N} \frac{\partial \mathcal{C}\_{\text{FLRW}}^{\text{poly}}}{\partial P\_{\phi}} = \mathbf{N} \frac{P\_{\phi}}{v}, \quad \mathbf{N} = \frac{v}{P\_{\phi}} = \frac{1}{\sqrt{2\rho}};\tag{134}$$

thus the effective Friedmann equation in the (*v*, *φ*) plane reads as follows:

$$\left(\frac{1}{v}\frac{dv}{d\phi}\right)^2 = \frac{3}{2}\left(1 - \frac{4\gamma^2 b\_0^2}{54}\frac{P\_\phi^2}{v^2}\right),\tag{135}$$

whose analytical solution is as follows:

$$v(\phi) = \sqrt{\frac{4\gamma^2 b\_0^2}{54}} \, P\_\phi \, \cosh\left(\sqrt{\frac{3}{2}}\,\phi\right). \tag{136}$$

As shown in Figure 10 the Bounce is clearly visible since the Universe possesses a minimum volume.

Now, we promote the system to a full quantum level by applying the Dirac procedure [67]. In particular, the Hamiltonian operator annihilates and selects the physical wave functions when applied to the generic quantum states, yielding the WDW equation:

$$
\hat{\mathcal{L}}\_{\text{FLRW}}^{\text{poly}} |\Psi\rangle = 0. \tag{137}
$$

If we promote the classical variables to a quantum level in the momentum representation, we obtain the following:

$$\mathcal{O} = -\frac{i\gamma}{3} \frac{d}{d\vec{\mathcal{E}}} , \quad \hat{\mathcal{E}} = \frac{1}{b\_0} \text{sin}(b\_0 \vec{\mathcal{E}}) , \quad \hat{P}\_\Phi = -i \frac{d}{d\phi} , \tag{138}$$

recalling that *v* is discrete and therefore *c*˜ must be regularized. So, the quantum version of the Hamiltonian constraint (131) is as follows:

$$\left[-\frac{3}{2b\_0^2} \left(\sin\left(b\_0\mathfrak{z}\right)\frac{d}{d\mathfrak{z}}\right)^2 + \frac{d^2}{d\phi^2}\right] \Psi(\mathfrak{z}, \phi) = 0;\tag{139}$$

Using again a mixed factor ordering and the substitution

$$\mathfrak{X} = \sqrt{\frac{2}{3}} \ln \left[ \tan \left( \frac{b\_0 \mathfrak{X}}{2} \right) \right] + \tilde{\mathfrak{X}},\tag{140}$$

we can recast recast the expression (139) in the form of a massless Klein–Gordon-like equation:

$$\frac{d^2}{d\tilde{\mathfrak{X}}^2} \Psi(\tilde{\mathfrak{x}}, \phi) = \frac{d^2}{d\phi^2} \Psi(\tilde{\mathfrak{x}}, \phi). \tag{141}$$

This way, the wave packet representing the wave function of the Universe can be written as follows: 2

$$\Psi(\pounds,\phi) = \int\_0^\infty dk\_\phi \, \frac{e^{-\frac{\left|k\_\phi - k\_\phi\right|}{2\sigma^2}}}{\sqrt{4\pi\sigma^2}} \, k\_\phi \, e^{ik\_\phi \circ} e^{-ik\_\phi \phi} \tag{142}$$

in the *x*˜-representation. We have constructed a superposition of plane waves in *φ* by means of Gaussian-like coefficients and we have restricted the analysis only to the positive energy-like eigenvalues *kφ*.

**Figure 10.** The polymer trajectory of the volume in (*v*, *c*˜) for the flat FLRW Universe. The volume presents a minimum, showing that a Big Bounce takes place.

In this case, to study the bouncing behavior of the model, we can evaluate the expectation values of both the volume and the energy density operators:

$$
\hat{\mathcal{V}} = \mathfrak{d}, \quad \mathfrak{d} = \frac{\hat{P}\_{\Phi}^{2}}{2\hat{\mathfrak{d}}^{2}}.\tag{143}
$$

In what follows, we express all the quantities as functions of the new variable *x*˜ using (140), so that we can calculate the expectation values through the Klein–Gordon scalar product as follows:

$$
\langle \Psi | \hat{O} | \Psi \rangle = \int\_{-\infty}^{\infty} d\tilde{x} \,\, i \left( \Psi^\* \,\partial\_{\Phi} (\hat{O} \Psi) - (\hat{O} \Psi) \,\partial\_{\Phi} \Psi^\* \right), \tag{144}
$$

where we assume normalized wave functions. The action of the operators (143) is derived from the basic operators shown in (138). In Figures 11 and 12, we show, respectively, the expectation values h*V*ˆ(*φ*)i and h*ρ*ˆ(*φ*)i as functions of time. In Figure 13, the value *V*ˆ(*φB*) of the volume at the Bounce is shown as a function of the initial value *kφ*; it is clear how the minimum volume scales linearly with the energy-like eigenvalue, in accordance with (136). Then, in Figure 14, we show the Bounce density h*ρ*ˆ(*φB*)i for different values of *kφ*. In accordance with the semiclassical critical density (133), choosing the volume itself as the configurational variable for the quantization of the system makes the density at the Bounce independent from the initial conditions on the quantum solution (142). This quantum analysis allows a direct comparison with the *µ*¯ scheme of LQC, and it shows how polymer quantum cosmology is unable to reproduce the inverse triad corrections. We notice that in Figures 11–14, the numerically integrated points are fitted with the continuous lines to check for consistency with the semiclassical solutions.

**Figure 11.** The expectation value of the Universe volume as a function of time (blue dots), fitted with a function in accordance with the semiclassical evolution (continuous red line).

**Figure 12.** The expectation value of the energy density as a function of time (blue dots), fitted with a function in accordance with the semiclassical evolution (continuous red line).

**Figure 13.** The expectation value of the Universe volume at the time *φB* of the Bounce as a function of *kφ* (blue dots), fitted with a function in accordance with the semiclassical evolution (continuous red line).

**Figure 14.** The expectation value of the energy density at the time *φB* of the Bounce as a function of *kφ* (blue dots), fitted with a function in accordance with the semiclassical evolution (continuous red line).

4.3.2. The Bianchi I Model in the Anisotropic Volume Variables: (*v*1, *v*2, *v*3)

In this subsection, we study the polymer semiclassical dynamics of the Bianchi I model in complete analogy with the analysis performed for the FLRW model in the volume variables, as presented in [45]. More specifically, the anisotropic character of the Bianchi I model leads to the possibility of taking into account two different sets of volume-like variables that coincide in the case of the isotropic model. The first that we take under consideration consists of three equivalent generalized coordinates (see [129]):

$$w\_i = \text{sgn}(p\_i) |p\_i|^{\frac{3}{2}}, \quad \eta\_i = \frac{c\_i}{\sqrt{|p\_i|}} \tag{145}$$

for *i* = 1, 2, 3, where *ηi* are the conjugate momenta and the new symplectic structure for the system is characterized by the following Poisson brackets:

$$\{\eta\_{i\nu}v\_{j}\} = \frac{3}{2}\gamma\delta\_{ij}.\tag{146}$$

In this case, we are not promoting one of the configurational variables to represent the Universe volume. On the contrary, we are imposing that the three independent coordinates are isomorphic to the isotropic volume for each direction so that *V* = (*v*1*v*2*v*3) 1 3 . In particular, each coordinate *vi* reduces to the Universe volume in the isotropic limit.

The Hamiltonian constraint for this framework in the polymer representation reads as follows:

$$\mathcal{L}\_{\text{BI}}^{\text{poly}} = -\frac{1}{4\gamma^2 V} \sum\_{i \neq j} \frac{v\_i \sin(b\_i \eta\_i) v\_j \sin(b\_j \eta\_j)}{b\_i b\_j} + \frac{P\_\phi^2}{2V} = 0,\tag{147}$$

where *i*, *j* = 1, 2, 3 and the substitution *ηi* → sin(*biηi* ) *bi* is performed.

If we fix the gauge through the choice *N* = *V Pφ* , the Hamilton equations will describe the dynamics of the model with respect to *φ* playing the role of the time variable:

$$\frac{dv\_i}{d\phi} = -\frac{v\_i \cos(b\_i \eta\_i)}{4\gamma P\_\phi} \left[ \frac{v\_j}{b\_j} \sin(b\_j \eta\_j) + \frac{v\_k}{b\_k} \sin(b\_k \eta\_k) \right] \tag{148a}$$

$$\frac{d\eta\_i}{d\phi} = \frac{\sin(b\_i \eta\_i)}{4\gamma b\_i P\_\Phi} \left[ \frac{v\_j}{b\_j} \sin(b\_j \eta\_j) + \frac{v\_k}{b\_k} \sin(b\_k \eta\_k) \right] \tag{148b}$$

for *i*, *j*, *k* = 1, 2, 3 and *i* 6= *j* 6= *k*. In analogy with the previous treatment, we can identity the following constants of motion:

$$\frac{v\_l \sin(b\_0 \eta\_l)}{b\_0} = \mathcal{K}\_{i\prime} \qquad P\_\phi = \mathcal{K}\_\phi \tag{149}$$

where *b*1 = *b*2 = *b*3 = *b*0. Taking general initial conditions that satisfy the constraint (147), the analytical solutions for the anisotropic volume coordinates read as follows:

$$\begin{aligned} v\_1(\phi) &= \bar{v}\_1 \cosh\left[\frac{(2\gamma^2 P\_\phi^2 b\_0^2 + \bar{v}\_2^2)\phi}{4\gamma P\_\phi b\_0 (\bar{v}\_1 + \bar{v}\_2)}\right], \\ v\_2(\phi) &= \bar{v}\_2 \cosh\left[\frac{(2\gamma^2 P\_\phi^2 b\_0^2 + \bar{v}\_1^2)\phi}{4\gamma P\_\Phi b\_0 (\bar{v}\_1 + \bar{v}\_2)}\right], \\ v\_3(\phi) &= \bar{v}\_3 \cosh\left[\frac{(\bar{v}\_1 + \bar{v}\_2)\phi}{4\gamma P\_\phi b\_0}\right], \end{aligned} \tag{150}$$

where

$$\begin{aligned} v\_1(0) &= \overline{v}\_1 & v\_2(0) &= \overline{v}\_2. \\ v\_3(0) = \overline{v}\_3 &= \frac{(2\gamma^2 b\_0^2 P\_\phi^2 - \overline{v}\_1 \overline{v}\_2)}{(\overline{v}\_1 + \overline{v}\_2)}. \end{aligned} \tag{151}$$

By combining the solutions (150), we can find the Universe volume behavior in function of *φ* as *V*(*φ*) = (*v*1(*φ*)*v*2(*φ*)*v*3(*φ*))1/3. As shown in Figure 15, the Big Bounce appears as a polymer regularization effect in place of the classical Big Bang. Analogously, it is possible to investigate the properties of the critical energy density of matter by computing the energy density of the scalar field in correspondence of the minimum volume at the time *φB* = 0 of the Bounce, i.e., the following:

$$\rho\_{\rm crit}^{\phi} = \frac{P\_{\phi}^{2}}{2V(\phi)^{2}}\bigg|\_{\phi\_{\rm B}} = \frac{P\_{\phi}^{2}}{2(\bar{v}\_{1}\bar{v}\_{2}\bar{v}\_{3})^{2/3}} = \frac{P\_{\phi}^{2}(\upsilon\_{1} + \upsilon\_{2})^{2/3}}{2[\upsilon\_{1}\upsilon\_{2}(2\gamma^{2}P\_{\phi}^{2}b\_{0}^{2} - \upsilon\_{1}\upsilon\_{2})]^{2/3}};\tag{152}$$

this results to be dependent on the initial conditions also in this set of anisotropic volume variables. In the next section, we implement directly the Universe volume as a generalized coordinate, trying to overcome this issue.

**Figure 15.** Bouncing trajectory of the Universe volume *V* = (*v*1*v*2*v*3) 1/3 as function of *φ* in the semiclassical polymer framework.

4.3.3. The Bianchi I Model in the Volume Variables: (*v*, *q*1, *q*2)

Now, we consider the set of volume variables introduced in [45]:

$$q\_i = \text{sgn}(p\_i)\sqrt{|p\_i|}, \qquad v = q\_1 q\_2 q\_3 \tag{153}$$

where *ηi* (with *i* = 1, 2, 3) are the conjugate momenta and *v* represents the Universe volume. The Poisson brackets are as follows:

$$\{\eta\_{\dot{\nu}}\eta\_{\dot{\jmath}}\} = \gamma \delta\_{\dot{\imath}\dot{}}.\tag{154}$$

and the Hamiltonian takes the following form:

$$\mathcal{L}\_{\rm BI} = -\frac{1}{4\gamma^2 V} (q\_1 \eta\_1 q\_2 \eta\_2 + q\_1 \eta\_1 v \eta\_3 + q\_2 \eta\_2 v \eta\_3) + \frac{P\_\phi^2}{2V} = 0. \tag{155}$$

When we rewrite this expression using the polymer substitution *ηi* → sin(*biηi* ) *bi* , we obtain the following:

$$\mathcal{C}\_{\text{BI}}^{\text{poly}} = -\frac{1}{4\gamma^2 V} \left( \sum\_{i=1,2} \frac{q\_i \sin(b\_i \eta\_i) \nu \sin(b\_i \eta\_3)}{b\_i b\_3} + \frac{q\_1 \sin(b\_1 \eta\_1) q\_2 \sin(b\_2 \eta\_2)}{b\_1 b\_2} \right) + \frac{P\_\phi^2}{2V} = 0. \tag{156}$$

In order to derive the dynamics of the model in function of the relational time *φ*, we impose the following:

$$\dot{\phi} := \mathbf{N} \frac{\partial \mathcal{C}\_{\mathrm{BI}}^{\mathrm{poly}}}{\partial P\_{\Phi}} = 1, \quad \mathbf{N} = \frac{V}{P\_{\Phi}}.\tag{157}$$

so the Hamilton equations for the couple (*v*, *η*3) take the following expressions:

$$\frac{dv}{d\phi} = -\frac{v\cos(b\_3\eta\_3)}{4\gamma P\_\phi} \left[ \frac{q\_1}{b\_1} \sin(b\_1\eta\_1) + \frac{q\_2}{b\_2} \sin(b\_2\eta\_2) \right] \tag{158a}$$

$$\frac{d\eta\_3}{d\phi} = \frac{\sin(b\_3 \eta\_3)}{4\gamma b\_3 P\_\phi} \left[ \frac{q\_1}{b\_1} \sin(b\_1 \eta\_1) + \frac{q\_2}{b\_2} \sin(b\_2 \eta\_2) \right] \tag{158b}$$

while for the conjugate variables (*q*1, *η*1), (*q*2, *η*2) we have the following:

$$\frac{dq\_i}{d\phi} = -\frac{q\_i \cos(b\_i \eta\_i)}{4\gamma P\_\phi} \left[ \frac{v}{b\_3} \sin(b\_3 \eta\_3) + \frac{q\_j}{b\_j} \sin(b\_j \eta\_j) \right],\tag{159a}$$

$$\frac{d\eta\_i}{d\phi} = \frac{\sin(b\_i \eta\_i)}{4\gamma b\_i P\_\Phi} \left[ \frac{v}{b\_3} \sin(b\_3 \eta\_3) + \frac{q\_j}{b\_j} \sin(b\_j \eta\_j) \right] \tag{159b}$$

for *i*, *j* = 1, 2 and *i* 6= *j*.

Once fixed the initial conditions on the variables (*q*1, *η*1), (*q*2, *η*2), (*v*, *η*3) according to (156), we can solve this system analytically since the 3D motion is decoupled in three one-dimensional trajectories, thanks to the use of the constants of motion as follows:

$$\frac{q\_l \sin(b\_0 \eta\_i)}{b\_0} = \mathcal{K}\_{i\prime} \quad P\_\phi = \mathcal{K}\_\phi \tag{160}$$

where *i* = 1, 2, 3 and *b*1 = *b*2 = *b*3 = *b*0. Differently from the previous analyses, we fix the constants of motion as follows:

$$\begin{aligned} \mathcal{K}\_1 &= \sqrt{\frac{2\gamma^2 P\_\phi^2 + \mathcal{K}^2}{3}} + \mathcal{K}, \\ \mathcal{K}\_2 &= \sqrt{\frac{2\gamma^2 P\_\phi^2 + \mathcal{K}^2}{3}} - \mathcal{K}, \\ \mathcal{K}\_3 &= \sqrt{\frac{2\gamma^2 P\_\phi^2 + \mathcal{K}^2}{3}}. \end{aligned} \tag{161}$$

where K is an arbitrary constant. As we will see below, this choice allows to write a convenient form for the Friedmann equation in order to develop a more rigorous analysis of the bouncing behavior of the model. Indeed, the existence of a non-trivial solution to the equation *H*2 = 0 identifies the expression of the critical energy density for which the scale factor velocity becomes null. More precisely, this information allows to identify also the anisotropy contribution to the total critical energy density added to the standard one associated to the matter fields. This way, it is possible to rigorously analyze the physical properties of the critical point.

In this set of volume variables, the Hubble parameter can be written as follows:

$$\begin{split}H^2 &= \left(\frac{1}{3v}\frac{dv}{dt}\right)^2 = \frac{(\mathcal{K}\_1 + \mathcal{K}\_2)^2}{144\gamma^2 v^2} \cos^2(b\_0\eta\_3) = \\ &= \frac{(\mathcal{K}\_1 + \mathcal{K}\_2)^2}{144\gamma^2 v^2} [1 - \sin^2(b\_0\eta\_3)] = \frac{(\mathcal{K}\_1 + \mathcal{K}\_2)^2}{144\gamma^2 v^2} (1 - \frac{b\_0^2}{v^2}\mathcal{K}\_3^2) \end{split} \tag{162}$$

where we restored the synchronous time-gauge N = 1 in the equation for the volume written in (158a). Now, if we substitute the conditions expressed above in (161), we obtain the following:

$$\bar{H}^2 = \frac{P\_\phi^2 + \bar{\mathcal{K}}^2}{54v^2} \left[ 1 - \frac{4\gamma^2 b\_0^2}{3} \left( \frac{P\_\phi^2 + \bar{\mathcal{K}}^2}{2v^2} \right) \right], \quad \bar{\mathcal{K}} = \frac{\mathcal{K}}{\sqrt{2\gamma^2}}.\tag{163}$$

This expression represents the polymer modified Friedmann equation for the Bianchi I model in the volume variables that, written in this convenient form, allows to derive the total critical energy density of the model. In particular, the additional term K¯ 2/2*v* 2 reasonably mimics the anisotropic contribution *ρ* aniso, so that we can compute *ρ* tot crit as follows:

$$
\rho\_{\rm crit}^{\rm tot} = \rho\_{\rm crit}^{\Phi} + \rho\_{\rm crit}^{\rm aniso} \,\,\,\,\,\tag{164}
$$

where the term regarding the matter scalar field takes the usual expression *P* 2 *φ*/(2*v* 2 ). Moreover, from (163) the total critical energy density results to be as follows:

$$
\rho\_{\rm crit}^{\rm tot} = \frac{\mathfrak{Z}}{4\gamma^2 b\_0^2} \,\mathrm{.}\tag{165}
$$

The solution for the Universe volume *v*(*φ*) when the initial conditions on the motion satisfy (161) is as follows:

$$w(\phi) = \sqrt{\frac{2\gamma^2 P\_{\phi}^2 + \mathcal{K}^2}{3}} \; b\_0 \cosh\left(\frac{\sqrt{2\gamma^2 P\_{\phi}^2 + \mathcal{K}^2}}{2\sqrt{3}} \frac{\phi}{\gamma P\_{\phi}}\right) . \tag{166}$$

clearly resembling a bouncing behavior as shown in Figure 16.

**Figure 16.** Semiclassical polymer trajectory of the Universe volume *v*(*φ*).

Now, we can verify whether the total critical energy density computed as (164) acquires the expression (165). In particular, we have the following:

$$\rho\_{\rm crit}^{\rm tot} = \frac{P\_{\Phi}^2}{2v(\phi)^2} \Bigg|\_{\phi\_{\rm B}} + \frac{\mathcal{K}^2}{2v(\phi)^2} \Bigg|\_{\phi\_{\rm B}} = \frac{3P\_{\Phi}^2}{2b\_0^2(2\gamma^2 P\_{\phi}^2 + \mathcal{K}^2)} + \frac{3\mathcal{K}^2}{4\gamma^2 b\_0^2(2\gamma^2 P\_{\phi}^2 + \mathcal{K}^2)} \,, \tag{167}$$

which clearly reduces to (165); therefore, the total critical density computed from the laws of motion for the Universe volume is consistent with the expression derived from the Friedmann Equation (163). Imposing a cut-off on the volume variables makes the critical energy density independent from the initial conditions on the scalar field, and therefore, produces a Big Bounce with universal properties. This result shows that taking the Universe volume itself as a configurational variable makes the Big Bounce acquire universal physical properties, in agreement with the behavior obtained for the set (*v*, *c*˜) in the FLRW model.

We notice that, even if the polymer-modified Friedmann equation in this convenient form is derived by considering the particular choice (161) for the constants of motion, the physical properties of the bouncing behavior of our model as derived from (163) have a general meaning since they are independent from the value assigned to the constant K.

#### *4.4. Polymer Cosmology in the Misner-like Variables*

In this section, we treat the semiclassical polymer dynamics of the most general anisotropic Universe, i.e., the Bianchi IX model. First, we implement the polymer paradigm in the Misner variables [42], and then we generalize the analysis by considering a inhomogeneous extension of the model in Misner-like variables, i.e., the anisotropies together with the Universe volume *v* ∝ *e* 3*α* [38].

#### 4.4.1. The Bianchi IX Model in the Misner Variables

The aim of this section is to discuss the main features of the Mixmaster semiclassical dynamics in the polymer representation [42]. We choose to define the Misner variables (*α*, *β*±) as discrete, so we perform the following formal substitutions:

$$P\_{\pm}^2 \rightarrow \frac{1}{b^2} \sin^2(bP\_{\pm}) \, , \tag{168a}$$

$$P\_a^2 \to \frac{1}{b\_a^2} \sin^2(b\_a P\_\hbar) \, , \tag{168b}$$

where *b* is the polymer parameter for the anisotropies, while *bα* is that one related to the isotropic variable *α*.

The super Hamiltonian constraint (14) becomes the following:

$$\mathcal{L}\_{\text{EIX}}^{\text{poly}} = \frac{\text{N}}{3(8\pi)^2} e^{-3a} \left[ -\frac{1}{b\_a^2} \sin^2(b\_a P\_a) + \frac{1}{b^2} \sin^2(b P\_+) + \frac{1}{b^2} \sin^2(b P\_-) + \frac{\mathcal{U}\_{\text{EIX}}(\mathcal{S}\_\pm)}{b\_a^2} \right] = 0,\tag{169}$$

while the ADM Hamiltonian has the following expression:

$$\mathcal{C}^{\text{ADM-poly}}\_{\text{BIX}} = \frac{\arcsin\sqrt{\frac{b\_a^2}{b^2} [\sin^2(bP\_+) + \sin^2(bP\_-)] + \mathcal{U}\_{\text{BIX}}}}{b\_a},\tag{170}$$

with the condition 0 ≤ *b* 2 *α b* 2 [sin2 (*bP*+) + sin2 (*bP*−)] + UBIX(*β*±) ≤ 1, due to the presence of the inverse sine function. In both (169) and (170), we have performed the substitution UBIX(*β*±) = *b* 2 *α*3(4*π*) 4 *e* 4*αU*BIX(*β*±).

The dynamics of the model is described by the following Hamilton equations:

$$\beta'\_{\pm} = \frac{\partial \mathcal{C}^{\text{ADM-poly}}\_{\text{BIX}}}{\partial P\_{\pm}} = \frac{b\_{a}}{b} \frac{\sin(2bP\_{\pm})}{\sin\left(2b\_{a}\mathcal{C}^{\text{ADM-poly}}\_{\text{BIX}}\right)} \text{ }\tag{171a}$$

$$P'\_{\pm} = -\frac{\partial \mathcal{C}^{\text{ADM-poly}}\_{\text{BIX}}}{\partial \beta\_{\pm}} = -\frac{3b\_{a}(4\pi)^{4}e^{4a}}{\sin\left(2b\_{a}\mathcal{C}^{\text{ADM-poly}}\_{\text{BIX}}\right)}\frac{\partial \mathcal{U}(\beta\_{\pm})}{\partial \beta\_{\pm}},\tag{171b}$$

$$\left(\mathcal{C}\_{\text{BIX}}^{\text{ADM-poly}}\right)' = 4e^{4a - 8\mathcal{G}\_{+}} \frac{\mathfrak{Z}b\_{a}(4\pi)^{4}}{\sin\left(2b\_{a}\mathcal{C}\_{\text{BIX}}^{\text{ADM-poly}}\right)'} \tag{171c}$$

where the symbol 0 denotes the derivative with respect to *α*.

Due to the steepness of the potential walls, we initially study the dynamics in the regime *U* ∼ 0. Under this condition, it can be demonstrated that there is a logarithmic relation between the time variable *t* and the isotropic one *α*, so we have that *α* ∼ ln(*t*) → −∞ for *t* → 0, even if *α* is described in the polymer formulation. This result points out that the discrete nature of the isotropic variable *α* does not prevent the Universe volume from vanishing, so the cosmological singularity is still present. We remark that the hypothesis *U* ∼ 0 becomes more reliable near the singularity because for *α* → −∞, the walls move outwards, and the region where the approximation is valid asymptotically covers the whole (*β*+, *β*−) plane.

Regarding the study of the chaotic features in the polymer modified picture, we have to analyze the relative motion between the particle-Universe and the potential walls, whose velocity is still |*β* 0 wall| = 1 2 since the polymer representation leaves the potential unchanged. On the other hand, if we study the anisotropy velocity of the particle while varying the values of the polymer parameters, we can see that the following holds:

$$\beta'(b\_{\mathfrak{a}}, b\_{\prime}P\_{\pm}) \ge 1 \,\,\forall \, P\_{\pm} \in \mathbb{R} \Leftrightarrow \frac{b\_{\mathfrak{a}}}{b} \ge 1 \,\,\, \tag{172}$$

where

$$\beta' \equiv \sqrt{\beta\_+^2 + \beta\_-^2} = \sqrt{\frac{\sin^2(bP\_+)\cos^2(bP\_+) + \sin^2(bP\_-)\cos^2(bP\_-)}{r^2(bP\_+,bP\_-)\left[1 - \frac{b\_a^2}{b^2}r^2(bP\_+,bP\_-)\right]}},\tag{173}$$

and *r*(*x*, *y*) = q sin2 (*x*) + sin2 (*y*). We notice that (*P*+, *P*−) (and consequently C ADM-poly BIX ) are constants of motion in the free particle regime. In addition, it is easy to verify that the relation (19) is recovered for *b*, *bα* → 0.

The figures represented in Figure 17 show that the anisotropy velocity represented on the vertical axis is always greater than 1 only if we choose the ratio *bα*/*b* to be greater than or equal to one. Therefore, the dynamics of the Mixmaster model is expected to be still chaotic because of the existence of the singularity and the presence of a never-ending series of rebounds against the potential walls. Instead, if we choose the polymer parameters such that *bα*/*b* < 1 (Figure 18), the series of rebounds occurs until the particle velocity becomes smaller than the velocity of the potential walls; then, the point-Universe reaches the singularity with no other rebound.

**Figure 17.** The 3D-profiles of the anisotropy velocity (173) with *bα*/*b* = 1 (left panel) and *bα*/*b* = 1.5 (right panel).

**Figure 18.** The 3D-profile of the anisotropy velocity (173) with *bα*/*b* = 0.1.

Moreover, it can be demonstrated that for *bα*/*b* < 1, we have *θ* max poly < *π* 3 , and this implies the absence of rebounds, even if the particle moves towards a specific wall. On the other hand, if *bα*/*b* ≥ 1, we have *π* 3 ≤ *θ* max poly < *π* 2 ; this means that a rebound is always possible, given the triangular symmetry of the system.

Regarding the polymer modified reflection law, by using analogous constants of motion with respect to the standard Misner case, we are able to derive the following law:

$$\begin{split} &\frac{1}{4b} \left[ \arcsin \left( \cos \theta^f \frac{\sin \theta^i}{\sin \theta^f} r (2bP^i\_+, 2bP^i\_-) \right) + \arcsin \left( \cos \theta^i r (2bP^i\_+, 2bP^i\_-) \right) \right] = \\ & \quad = \frac{1}{b\_\mathbb{R}} \left[ \arcsin \left( \frac{b\_\mathbb{R}}{b} r (bP^f\_+, bP^f\_-) \right) - \arcsin \left( \frac{b\_\mathbb{R}}{b} r (bP^i\_+, bP^i\_-) \right) \right], \end{split} \tag{174}$$

where *r*(*x*, *y*) = q sin2 (*x*) + sin2 (*y*).

In order to make a comparison with the standard case (20), an expansion up to the second order for *b* and *bα* is required:

$$\frac{\sin\left(\theta^{i} + \theta^{f}\right)}{2} = \sin\theta^{i}(1 + \Pi\_{f}^{2})(1 + \mathcal{R}\_{f}) - \sin\theta^{f}(1 + \Pi\_{i}^{2})(1 + \mathcal{R}\_{i})\tag{175}$$

where

$$
\Pi^2 = \frac{1}{6} b\_\alpha^2 (P\_+^2 + P\_-^2) \,\text{.}\tag{176a}
$$

$$R = \frac{2}{3}b^2 \frac{(P\_+)^4 + (P\_-)^4}{(P\_+)^2 + (P\_-)^2} \,\prime \tag{176b}$$

and it is easy to show that in the limit *b*, *bα* → 0, we find the standard reflection law (20) obtained by Misner.

Recovering this limit leads us to infer that the map above still admits stochastic properties; therefore, when *bα* < *b*, sooner or later the parameter region where no rebound takes place is reached, as in [27].

Given the results of this analysis and that of previous sections, we can infer that not only the physical nature of the Big Bounce, but also its very existence depends on the geometrical dimension of the variable chosen as discrete: we have seen here that using a logarithmic variable does not indeed avoid the singularity. Therefore, in the last model that we present in the next subsection, we choose to use a mixed representation involving the anisotropies as they are defined by Misner, together with an isotropic volume variable instead of *α* ∝ ln *v*.

#### 4.4.2. The Inhomogeneous Mixmaster Model in the Volume Variable

In [38], the semiclassical polymer dynamics of the inhomogeneous Mixmaster model is analyzed by choosing the cubed scale factor (i.e., the Universe volume) alone as the discretized configurational variable. In addition, a massless scalar field and a cosmological constant are included in the dynamics, accounting respectively for a quasi-isotropization and inflationary-like mechanisms. The resulting dynamics is a singularity-free Kasnerlike phase that is linked with a homogeneous and isotropic de Sitter evolution. Moreover, the chaotic character of the Mixmaster dynamics is absent in this framework. Thus, the study presented in [38] demonstrates that the generic cosmological solution is singularity-free and non-chaotic once the polymer discretization of the Universe volume variable is performed at a semiclassical level. This result is alike to that achieved also in LQC, as presented in Section 3.5.

In order to summarize more in detail all the main results presented in [38], let us start by characterizing the generic cosmological problem and its relationship with the homogeneous Mixmaster model. The most general line element can be written as follows:

$$ds^2 = \mathbf{N}^2 dt^2 - h\_{ab} (d\mathbf{x}^a + \mathbf{N}^a dt) (d\mathbf{x}^b + \mathbf{N}^b dt) \, , \tag{177}$$

where N*a*,*b* (*t*,~*x*) is the shift vector and the following holds:

$$h\_{ab}(t, \vec{x}) = e^{m\_1(t, \vec{x})} l\_a^1 l\_b^1 + e^{m\_2(t, \vec{x})} l\_a^2 l\_b^2 + e^{m\_3(t, \vec{x})} l\_a^3 l\_b^3 \,. \tag{178}$$

with *a*, *b* = 1, 2, 3. Rigorously, this is not the most general picture since we are not considering the rotation of the Kasner vectors *l* 1,2,3 *a*,*b* (*t*,~*x*). Differently from the homogeneous Bianchi line element (12), the lapse function N(*t*,~*x*) depends explicitly on the spatial coordinates, as well as the Misner-like variables, that in their generalized version are defined as follows:

$$m\_d(t, \vec{\mathfrak{x}}) = \frac{2}{3} \ln[v(t, \vec{\mathfrak{x}})] + 2\beta\_d(t, \vec{\mathfrak{x}}) \,. \tag{179}$$

where *a* = 1, 2, 3 and *βa* = (*β*+ + √ 3 *β*−, *β*+ − √ 3 *β*−, −2*β*+). Here, we are considering *v* as the isotropic variable proportional to the spatial volume of the Universe (instead of *α* as in the proper Misner variables), while *β*+, *β*− are the anisotropies. Let us include a massless, self-interacting scalar field *φ* in the dynamics, taking into account the slow-roll condition *φ*˙ 2 *U*(*φ*) as well as the following hypothesis:

$$\left|\nabla\phi\right|^{2}\ll\mathcal{U}[\phi(t,\vec{x})]\sim\Lambda(\vec{x})\,.\tag{180}$$

both typical of the inflationary paradigm.

The action in the inhomogeneous case reads as follows:

$$S\_{\mathbf{B}}^{\text{inhom}} = \int d^3 \mathbf{x} \, dt \left( P\_{\upsilon} \psi + \sum\_{j} P\_{j} \not{p}\_{j} - \mathbf{N} \mathcal{C}\_{\mathbf{B}}^{\text{inhom}} - \mathbf{N}^{j} \mathcal{C}\_{i} \right) \,. \tag{181}$$

where *j* = +, −, *φ* and *βφ* ≡ *φ*. The super Hamiltonian constraint is as follows:

$$\mathcal{L}\_{\mathsf{B}}^{\text{inhom}} = \frac{3}{4} \left[ -\upsilon P\_{\upsilon}^{2} + \sum\_{j} \frac{P\_{j}^{2}}{9\upsilon^{2}} + \frac{\upsilon^{1/3}}{3} \mathcal{U}\_{\mathsf{B}}^{\text{inhom}} + \upsilon \Lambda(\vec{x}) \right],\tag{182}$$

while C*i* represents the supermomentum one. We remark how the self-interacting scalar field, under the slow-roll condition, is equivalent to a free scalar field plus a cosmological constant Λ(~*x*).

The potential term is due to the three-dimensional scalar curvature and can be split as *U*inhom = *W* + *U*inhom B , where the contribution due to the spatial gradients of the configurational variables is encoded in *W*, and the term *U*inhom B represents the inhomogeneous generalization of the Bianchi model's potential (the Bianchi VIII and Bianchi IX models are the most general choices). By using the Belinskii–Khalatnikov–Lifshitz (BKL) conjecture [130], the term *W* is considered to be negligible and so the super Hamiltonian (182) reduces to that one of the Bianchi IX model, where the spatial coordinates appear only as parameters. The reliability of this assumption can be verified by explicitly evaluating the term *W* by means of the supermomentum constraint. Moreover, by means of the gauge choice N*i* = 0, each point of space evolves independently and can be described, using a homogeneous Bianchi IX model. This way, point by point, the inhomogeneous cosmological evolution is approximated by a Mixmaster-like one. We notice that the corresponding solutions of the dynamics must satisfy also the supermomentum constraint, which identically vanishes for a homogeneous spacetime. In addition, with the hypothesis that the physical scale of the inhomogeneities is much bigger than the average Hubble horizon, each causal connected region verifies a Mixmaster-like evolution, instead of each point of space [131].

Let us now analyze the semiclassical polymer dynamics of the model. Using the substitution (89) for the momentum conjugate to the variable *v*, i.e., *Pv* → sin(*b*0*Pv*) *b*0 , the Hamiltonian constraint becomes the following:

$$\mathcal{C}\_{\text{BIX}}^{\text{poly}} = \frac{3}{4} \left[ -v \frac{\sin^2(b\_0 P\_v)}{b\_0^2} + \sum\_j \frac{P\_j^2}{9v^2} + \frac{v^{1/3}}{3} \mathcal{U}\_{\text{BIX}} + v\Lambda \right],\tag{183}$$

where the spatial gradients are neglected and the potential term is that one of the Bianchi IX model. Firstly, we can notice that the following condition is imposed by the form of the Hamiltonian when the potential term is negligible:

$$\sin^2(b\_0 P\_v) = b\_0^2 \left(\frac{\sum\_j P\_j^2}{9v^2} + \Lambda\right). \tag{184}$$

As a consequence, the right-hand side of (184) must be smaller than 1, and therefore, the Universe volume acquires a lower bound:

$$v > v\_{\mathbf{B}} \equiv \frac{b\_0}{3} \sqrt{\frac{\sum\_{j} P\_j^2}{b\_0^2 \Lambda}}\,,\tag{185}$$

with the condition *b* 2 0Λ < 1, which is physically reasonable in the Planck units. So, in the polymer semiclassical dynamics, the initial singularity is replaced by a Big Bounce. We notice that the minimum value of the Universe volume depends not only on the polymer scale *b*0, but also on the cosmological constant Λ and on the conjugate momenta *Pj* , which are constants of motion when the potential term is negligible. This feature is in accordance with the results obtained in LQC [49]. The absence of singularity thanks to the semiclassical polymer formalism can be shown more in general by solving the Hamilton equation for the Universe volume, from which we obtain the following:

$$w(t) = \sqrt{\frac{\sum\_{j} P\_{j}^{2} \left[ \cosh \left( 3\sqrt{\Lambda} \sqrt{1 - b\_{0}^{2} \Lambda} t \right) - 1 + 2b\_{0}^{2} \Lambda \right]}{18 \Lambda (1 - b\_{0}^{2} \Lambda)}}. \tag{186}$$

In particular, as shown in Figure 19 the Bounce occurs when the volume reaches its minimum, which is given by Equation (185).

**Figure 19.** Semiclassical polymer trajectory of the Universe volume *v* in the synchronous time *t*. Image from [38].

Moreover, by performing the ADM reduction and choosing *v* as the relational time, the Hamilton equations for the anisotropies and the scalar field can be solved, showing that the former do not diverge approaching the Big Bounce, as it happens in the standard case. Additionally, far from the Bounce, the anisotropies tend to a constant value that can be absorbed by an appropriate and local rescaling of the coordinate system. In this sense, the quasi-isotropization mechanism described in [62,131] results to be preserved under the polymer modification.

In conclusion, the semiclassical polymer dynamics is well described by a bridge solution connecting an initial non-singular, inhomogeneous and anisotropic Kasner-like regime with a later homogeneous and isotropic de Sitter phase. We notice that this result is local, i.e., valid in a nearly homogeneous region. Nonetheless, with the hypothesis of having an inhomogeneity scale much larger than the average Hubble horizon, each causal connected region undergoes such an isotropization phase and can be expanded by the inflationary mechanism on a scale much larger than the Hubble radius. Thus, if the inflationary phase is long enough, it is possible to obtain a homogeneous and isotropic region larger than the Hubble horizon today.

Finally, in [38], it is also demonstrated that the chaotic properties of the dynamics are suppressed when the polymer representation is implemented at a semiclassical level. This property is a consequence of the following considerations. As we have seen in Section 2.3, the chaoticity of the Bianchi IX model is derived by analyzing the relative velocity between the point-Universe and the potential walls. It is possible to verify that, similar to the

analysis in Section 4.4.1, in this model, the velocity of the walls is finite at the Bounce, while that of the point-Universe diverges in the polymer framework. So, rebounds always occur, even in the presence of a scalar field, even though in the classical framework, its presence is able to remove the chaos [113]. However, due to the polymer modifications, the following condition holds:

$$\sum\_{j} P\_j^2 + 9\Lambda v^2 > \max\left(3v^{4/3}\mathcal{U}\_{\text{BIX}}\right). \tag{187}$$

Therefore, there are two possible scenarios for the removal of chaos, similar to the case of the Loop quantized homogeneous Bianchi IX model described in Section 3.5: either the condition (187) is satisfied before the Bounce, resulting in a final stable Kasner-like epoch (similarly to the AVTD regime reached classically when a scalar field is present), or chaos is removed by the fact that the presence of a Big Bounce makes the number of rebounds on the potential walls finite.

#### *4.5. The Link between Polymer Quantum Mechanics and Loop Quantum Cosmology: Canonical Equivalence*

In Sections 4.2–4.4, we have presented the main results gained thanks to the application of the polymer representation at the level of the semiclassical dynamics, with some extensions to the full quantum regime in Sections 4.2 and 4.3. Along the presentation, we have proceeded by considering more and more general cosmological models, starting from the FLRW one and then extending the analysis to the Bianchi I and IX models afterwards. For each model, we have compared the obtained results in different sets of variables, always performing the canonical transformation between the phase space variables before implementing the polymer approximation for the variables conjugate to those represented as discrete on the polymer lattices. If, on one hand, the polymer formulation permits to successfully overcome the GR shortcoming of the initial singularity, on the other hand, it reproduces different physical pictures for the same cosmological model when represented in different sets of variables canonically linked in the classical regime. In this respect, the polymer behavior of the Bianchi IX cosmological model in the standard Misner variables is of particular relevance since it shows a singular dynamics in contrast with the bouncing cosmology and the chaos removal ensured by the use of a volume-like variable (note that in the Misner variables the presence of chaos depends on the ratio of the different discretization parameters associated to the isotropic variable and the anisotropy coordinates *β*± [42]).

This fact raises an issue about what it really means to do a canonical transformation in the polymer framework. For a proposal about recovering the equivalence between different sets of variables after the polymer is implemented on the Hamiltonian constraint, see [44]. In this work, it is demonstrated how, by fixing the preferred system of generalized coordinates and assigning a constant lattice step (actually in the considered case, there is only the scale factor, but a generalization of the analysis can be easily inferred), the dynamical features of the model remain the same in any other set of coordinates, at the price of considering the discretization step as a function of the adopted variables. For example, starting from the volume coordinates and then changing to the Ashtekar connections maintains the canonical equivalence between the dynamics thanks to the dependence of the polymer parameter on *p*. Therefore, the dynamics is ruled by the set of variables in which the polymer spacing results in being constant. We note that this study was performed only at a semiclassical level since the implementation of a translational operator depending on the coordinates is not trivial, so the equivalence in the full quantum theory is still an open question.

Nevertheless, this result achieved in [44] about the recovery of an equivalence class of configurational variables could have relevant implications for the *µ*¯ scheme of LQC. In particular, in the *µ*¯ scheme of LQC the translational operator acting on the physical states acquires a parameter depending on the momentum variable (due to the physical rescaling of the area element) that, in the polymer framework, plays exactly the role of *p*. This analogy suggests that a change of coordinates is at the basis of the improved LQC approach [49–51] in a way that a translational operator of constant step is restored and the theory is made technically viable.

However, such a change of variables opens a possible criticism about the quantization procedure, which is no longer referred to the natural *SU*(2) connection so that LQC would be departing from full LQG. While on a quantum level, this question is somewhat an open issue, on the level of an effective LQC theory, the result mentioned above can somewhat support the validity of the change in variables in view of the link between PQM and LQC. Indeed, since the physical picture of the model is dictated by the scheme with a constant polymerization parameter, we can conclude that the semiclassical dynamics of the *µ*¯ formulation is the same of that obtained in the Ashtekar variables if the lattice parameter is canonically transformed so that no ambiguity arises about the universal character of the Big Bounce, i.e., about its independence on the initial conditions.

#### **5. Concluding Remarks**

In this review, we analyzed basic aspects of LQC and presented polymer quantum cosmology in depth in order to trace some solid implications on the early Universe evolution precisely regarding the regularization of the singularity.

In this respect, we can firmly say that a bouncing cosmology always emerges in LQC and in polymer quantum cosmology, when the Universe volume is considered one of the configurational variables. Furthermore, the Big Bounce has also the features of a universal cut-off since the critical energy density of the Universe takes a maximum value depending only on fundamental constants and the Immirzi parameter. This result is obtained in the improved scheme of LQC as shown in Section 3.3 and also demonstrated in the context of PQM for the isotropic model in [37]. Moreover, similar results to those obtained in the *µ*¯ scheme are outlined in [44,45] when PQM is implemented on volume-like variables.

The situation is slightly different when we polymerize different sets of generalized variables in PQM. In particular, when we consider the polymerization of the Ashtekar– Barbero–Immirzi connections, the Big Bounce still emerges in the semiclassical dynamics, but the critical energy density now depends on the initial conditions on the motion, or equivalently on the wave packets when the full quantum approach is implemented (see [44,45]). In these analyses, the polymerization of the natural connection seems to reproduce a physical picture very similar to the *µ*0 scheme of LQC presented in Section 3.2.

The importance of the choice of the adopted variable for the polymer quantization of the Universe is outlined by observing that the polymerization of the Bianchi IX model (both in the homogeneous and inhomogeneous cases) expressed in the standard Misner variables yields a model that is still singular [34,42], while discretizing the Universe volume in the metric formulation does avoid the singularity [38]. This result highlights the link between the dimensionality of the discretized variable and the regularization of the singularity in PQM; also, it opens the non-trivial question about which equivalence class among all possible choices of canonically related coordinates leads to the same physical picture.

An interesting result toward the solution of this puzzling question is elucidated from the analysis in [44]. We have discussed how the equivalence is ensured for the isotropic model on a semiclassical level by the possibility to restate the problem in terms of a new polymerization step dependent on the configurational coordinates. This result can have positive implications in stating the equivalence of the effective dynamics in the *µ*0 and *µ*¯ schemes of LQC [20,50], as discussed in Section 4.5, at least on a semiclassical level.

On the base of the results illustrated here, we can conclude that, if on one hand, some basic features of cut-off quantum cosmology are well traced, many other detailed points must be addressed. For instance, it is necessary to better characterize the physical properties of the Big Bounce, both by providing a pure quantum description instead of just a semiclassical one, and by further investigating the thermodynamical nature of the quantum cosmological fluid. Moreover, it is crucial to understand how a collapsing Universe before our expanding branch is generated and how it can influence the morphology of the present

Universe (see some pioneering works [37,39,132–138]). A typical example of this issue affecting bouncing cosmologies is provided in [37], where it is argued that the flatness paradox is still present, even though the horizon paradox is naturally solved, thanks to the pre-existing collapsing Universe. These kinds of physical questions become even more meaningful if we postulate that the Universe dynamics is associated to a cyclical evolution between a Big Bounce and a later classical turning point [139,140]: in this context, the behavior of the Universe entropy becomes a delicate feature to be addressed.

**Author Contributions:** All authors contributed to all parts of the research and manuscript. G.M. has been slightly more active than the others in the conceptualization while G.B. and E.G. have contributed slightly more on the software and visualization. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Note**

1 In [45] it is shown that by combining the semiclassical solutions for *ci* (*φ*) with (125) we obtain that *xi* (*φ*) ∝ *φ* with slopes depending on the initial conditions on the motion.

#### **References**

## *Review* **A Brief Overview of Results about Uniqueness of the Quantization in Cosmology**

**Jerónimo Cortez 1 , Guillermo A. Mena Marugán 2 and José M. Velhinho 3,\***

**Abstract:** The purpose of this review is to provide a brief overview of recent conceptual developments regarding possible criteria to guarantee the uniqueness of the quantization in a variety of situations that are found in cosmological systems. These criteria impose certain conditions on the representation of a group of physically relevant linear transformations. Generally, this group contains any existing symmetry of the spatial sections. These symmetries may or may not be sufficient for the purpose of uniqueness and may have to be complemented with other remaining symmetries that affect the time direction or with dynamical transformations that are, in fact, not symmetries. We discuss the extent to which a unitary implementation of the resulting group suffices to fix the quantization—a demand that can be seen as a weaker version of the requirement of invariance. In particular, a strict invariance under certain transformations may eliminate some physically interesting possibilities in the passage to the quantum theory. This is the first review in which this unified perspective is adopted to discuss otherwise different uniqueness criteria proposed either in homogeneous loop quantum cosmology or in the Fock quantization of inhomogeneous cosmologies.

**Keywords:** quantum cosmology; uniqueness of the quantization

#### **1. Introduction**

Quantization is the process of constructing a description that incorporates the principles of Quantum Mechanics starting with a given classical system. In the present work, we consider exclusively the so-called canonical quantization process. This means that the classical system can be described in canonical form (e.g., in terms of variables that form canonical pairs) and that one aims at promoting classical variables to operators in a Hilbert space, preserving the canonical structure as much as possible. The prototypical system is, of course, the phase space R2*n* , with coordinates {(*qi* , *pi*), *i* = 1, · · · , *n*}, equipped with the Poisson bracket {*qi* , *pj*} = *δij*. The standard quantization is realized in the Hilbert space *L* 2 (*R n* ) of square integrable functions, with configuration variables promoted to multiplicative operators

$$
\mathfrak{q}\_i \mathfrak{\psi} = \mathfrak{q}\_i \mathfrak{\psi}\_{\prime} \tag{1}
$$

and momentum variables acting as derivative operators

$$
\not p\_i \psi = -i \frac{\partial}{\partial q\_i} \psi. \tag{2}
$$

Here, and in the following, we set the reduced Planck constant ¯*h* equal to one. These operators satisfy the canonical commutation relations (CCRs), thus, implementing Dirac's quantization rule that *Poisson brackets go to commutators* [1]. This quantization is irreducible, in the sense that any operator that commutes with all the *q*ˆ*i* 's and *p*ˆ*i* 's is necessarily

**Citation:** Cortez, J.; Mena Marugán, G.A.; Velhinho, J.M. A Brief Overview of Results about Uniqueness of the Quantization in Cosmology. *Universe* **2021**, *7*, 299. https://doi.org/ 10.3390/universe7080299

Academic Editor: Paulo Vargas Moniz

Received: 19 July 2021 Accepted: 8 August 2021 Published: 13 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

proportional to the identity operator. In addition, this quantization satisfies a technical continuity condition, since it provides a continuous representation of the Weyl relations associated with the CCRs (see Section 2 for details).

These two conditions uniquely determine the quantization. This is the celebrated Stone–von Neumann uniqueness theorem (see, e.g., [2]): every irreducible representation of the CCRs coming from a continuous representation of the Weyl relations is unitarily equivalent to the aforementioned quantization, meaning that, given operators *q*ˆ 0 *i* and *p*ˆ 0 *i* with the above properties in a Hilbert space H, there exists a unitary operator *U* : H → *L* 2 (*R n* ) relating the two sets of operators *q*ˆ*i* , *p*ˆ*i* and *q*ˆ 0 *i* , *p*ˆ 0 *i* (see Equation (6) below).

On the other hand, it has been known since Dirac's work (and it was rigorously proved by Groenewold and van Hove [3,4]) that imposing Dirac's quantization rule on a large set of observables is not viable in the sense that it is impossible to satisfy the relations

$$\hat{i}\left[\hat{f},\hat{\mathfrak{g}}\right] = \hat{i}\widehat{\{f,\mathfrak{g}\}}\tag{3}$$

for a large set of classical observables, and certainly not for the whole algebra of classical observables (see [5] for a thorough discussion).

Nevertheless, in any given physical system, there are certainly observables of interest, other than the coordinate variables *q*ˆ*i* and *p*ˆ*i* , which need to be quantized as well. Certain canonical transformations typically stand out in a given classical system, e.g., dynamics or symmetries, and these require, as well, a proper quantum treatment (these two aspects are, in fact, related since observables of physical interest often emerge as generators of some special groups of canonical transformations). This poses no problem as far as one considers the standard quantization of linear canonical transformations and the corresponding generators in R2*n* , precisely because of the above uniqueness theorem. One can easily illustrate this issue by considering a canonical transformation

$$
\begin{pmatrix} q \\ p \end{pmatrix} \rightarrow \begin{pmatrix} q' \\ p' \end{pmatrix} = \mathcal{S} \begin{pmatrix} q \\ p \end{pmatrix} \tag{4}
$$

in R2*n* , where *S* is a symplectic matrix. Then, the operators *q*ˆ 0 *i* and *p*ˆ 0 *i* , defined as

$$
\begin{pmatrix} \hat{\mathfrak{q}}' \\ \mathfrak{p}' \end{pmatrix} = \mathcal{S} \begin{pmatrix} \hat{\mathfrak{q}} \\ \hat{\mathfrak{p}} \end{pmatrix} . \tag{5}
$$

provide a new representation of the CCRs (in this case, in the same Hilbert space) with the same properties of irreducibility and continuity. Thus, it is guaranteed that there exists a unitary operator *U* such that

$$\mathfrak{q}'\_{\dot{i}} = \mathcal{U}^{-1} \mathfrak{q}\_{\dot{i}} \mathcal{U}, \quad \mathfrak{p}'\_{\dot{i}} = \mathcal{U}^{-1} \mathfrak{p}\_{\dot{i}} \mathcal{U}. \tag{6}$$

The unitary operator *U* is naturally interpreted as the quantization of the symplectic transformation *S*. We will also refer to it as the unitary implementation (at the quantum level) of the canonical transformation *S*. If, instead of a single linear transformation *S*, one has a 1-parameter group, generated, e.g., by a quadratic Hamiltonian function *H*, one obtains a corresponding 1-parameter group of unitary operators *U*(*t*), which is typically continuous, so that a self-adjoint generator *H*ˆ such that *U*(*t*) = *e* −*iHt* ˆ can be extracted. Non-quadratic classical Hamiltonians of the type *H* = (1/2) ∑ *p* 2 *i* + *V*(*qi*) do not fall in this last category; however, there is nevertheless a standard and well-defined procedure to obtain their quantum version (see [6] for details), which is simply to define *H*ˆ as *H*ˆ = (1/2) ∑ *p*ˆ 2 *i* + *V*(*q*ˆ*i*).

Some ambiguities may occur in the quantization of more general functions, involving e.g., products of *q*'s and *p*'s; however, these ambiguities are typically not overly severe. These are precisely the type of ambiguities that may happen in standard homogeneous quantum cosmology (QC). In fact, in a homogeneous cosmological model, the number

of both gravitational and matter degrees of freedom (DoF) is hugely reduced, ending up with just a finite number of global DoF, precisely due to homogeneity. In this socalled minisuperspace setup, the configuration variables on the gravitational side are typically given by the different scale factors, the number of which depend on the degree of anisotropy.

The object of interest here is the Hamiltonian constraint, the quantization of which leads to what is often called the Wheeler–de Witt equation *H*ˆ Ψ = 0. This quantization may not be entirely trivial, owing to a possibly complicated dependence of the constraint with respect to the basic configuration and momentum variables. Nevertheless, the ambiguity that may emerge from the quantization of the Hamiltonian constraint typically involves a choice of factor ordering in *H*ˆ with respect to the basic quantum operators. Although different choices may lead to different versions of the Hamiltonian constraint, it is often the case that this does not affect the physical predictions substantially, in the sense that the predictions remain qualitatively the same.

Thus, the formalism of standard homogeneous QC, based on the standard quantization for a finite number of DoF, is to a large extent free of major ambiguities, as follows from the Stone–von Neumann theorem.

This last paradigm can be broken in two different ways, for very distinct reasons. First, the quantization—even of the set of basic configuration and momentum variables of systems with an infinite number of DoF escapes the conclusions of the Stone–von Neumann theorem. On the contrary, in that case, there are representations of the CCRs leading to physically inequivalent descriptions of the same system. This occurs when local DoF are considered, of which one can mention two distinct situations of interest in cosmology: (i) the quantization of gravitational DoF (and possibly also of matter fields) in inhomogeneous cosmologies (such as in the example of Gowdy models [7]), and (ii) the quantum treatment of fields propagating in a non-stationary curved spacetime (e.g., of the FRW or de Sitter type), which is considered as a classical background.

While the first case embodies genuine applications to QC, i.e., a full quantum treatment of the gravitational DoF (in cases with considerable symmetry, like the Gowdy model, in the so-called midisuperspace setup), the second situation finds applications in the treatment of quantum perturbations in cosmology (both of gravitational and matter DoF), with the homogeneous background kept as a classical entity, such as, e.g., in inflationary scenarios. There is, thus, the need of selecting physically relevant quantizations corresponding to a given system containing an infinite number of DoF. Note that available selection criteria (leading to uniqueness) typically rely on stationarity and are, therefore, not applicable in the above described situations.

Section 4 is devoted to review selection criteria that were recently introduced and proved viable, leading to unique and well defined quantizations in the aforementioned cases. Such criteria are based on remaining symmetries, present in the cosmological system, and crucially on the unitary implementation of the dynamics, which can be seen as a weaker version of the requirement of invariance under time-translations, which can be applied only in stationary settings.

The other avenue for departure from the conclusions of the Stone–von Neumann theorem, which is relevant in the cosmological context, is exemplified by the quantization approach for homogeneous cosmologies known as loop quantum cosmology (LQC)1 . Although the same models with a finite number of DoF are considered as in standard homogeneous cosmology, the obtained quantizations are not physically equivalent.

The type of quantization used in LQC is not unitarily equivalent to that of standard QC, and the reason why this is possible is that, in LQC, one of the conditions of the Stone– von Neumann theorem is broken in a hard way. The LQC-type of quantization starts from the Weyl relations, which are the exponentiated version of the CCRs and considers representations of the Weyl relations that are not continuous, thus, violating one of the conditions of the Stone–von Neumann theorem. In particular, it is the configuration part of the representation that is not continuous.

As a result, and although the correspondent of the unitary group *e itq*ˆ is well defined in the LQC-type of representation as a unitary group *U*(*t*), the would-be generator *q*ˆ cannot be defined, owing to the lack of continuity. This, in turn, is at the heart of the emergence of the discretization (in the canonically conjugate variable) that is so characteristic of LQC. Together with quantization methods adapted from those of loop quantum gravity (LQG) [8,9], this discretization is responsible, at the end of the day, for the results about singularity avoidance for which LQC is known.

One question that naturally arises is the following: are there other representations of the Weyl relations with different physical properties? Or is this particular LQC-type of representation naturally selected in some way? In Section 3, we discuss and comment on a uniqueness result for isotropic LQC recently put forward by Engle, Hanusch, and Thiemann [10], following a previous discussion concerning the Bianchi I case [11].

We would like to stress that, to the best of our knowledge, this is the first time that the results reviewed in Section 3 and those mentioned in Section 4 have been considered and discussed together. In particular, this joint and integrated review brings about a discussion on the two possible ways to use relevant transformations in order to select a unique quantization. In fact, results like those described in Section 3 are rooted on a requirement of strict invariance, while the results of Section 4 relax that condition (regarding what dynamical transformations are concerned), requiring only the weaker condition of unitary implementation of the transformations in question. These two approaches are discussed, providing a better understanding on the results achieved thus far in cosmology.

For completeness, we will start with a very brief review of the formalism for the study of Weyl algebras and their representations.

Also, we include an Appendix A sketching the proof of the uniqueness of the representation results mentioned in Section 4, in the simplest case of a scalar field in *S* 1 with time-dependent mass, with the purpose of providing the main steps and typical arguments of the proofs to the interested reader, without overloading the main text.

#### **2. Weyl Algebra and Standard Representations**

Let U and V be a pair of unitary representations of the commutative group R, in the same Hilbert space H, i.e., U(*a*) and V(*b*) are unitary operators for all real values of *a* and *b* and such that U(*a* + *a* 0 ) = U(*a*)U(*a* 0 ) and V(*b* + *b* 0 ) = V(*b*)V(*b* 0 ). The pair U, V is said to satisfy the Weyl relations if

$$\mathcal{V}(b)\mathcal{U}(a) = e^{iab}\mathcal{U}(a)\mathcal{V}(b). \tag{7}$$

The standard representation of the Weyl relations, corresponding to the usual Schrödinger representation of the CCRs, is obtained as follows. Consider the Hilbert space *L* 2 (R) of square integrable functions *ψ*(*q*) with respect to the usual Lebesgue measure *dq*. The expressions

$$(\mathcal{U}(a)\psi)(q) = e^{iaq}\psi(q) \tag{8}$$

and

$$(\mathcal{V}(b)\psi)(q) = \psi(q+b) \tag{9}$$

define unitary representations of R, which clearly satisfy the Weyl relations (7). These representations are, moreover, jointly irreducible and continuous, i.e., *a* 7→ U(*a*) and *b* 7→ V(*b*) are continuous functions. The appropriate notion of continuity of these operator valued functions is that of strong continuity, and irreducibility means that no proper subspace of H supports the action of both U(*a*) and V(*b*) ∀*a*, *b*.

It is precisely due to continuity that Stone's theorem guarantees that it is possible to define infinitesimal generators *q*ˆ and *p*ˆ such that U(*a*) = *e iaq*ˆ and V(*b*) = *e ibp*ˆ . In this case, it turns out that *q*ˆ is the multiplication operator *q* and *p*ˆ = −*i d dq* .

The celebrated Stone–von Neumann Theorem ensures that any other representation of the Weyl relations on a separable Hilbert space, with the same properties of irreducibility and continuity, is unitarily equivalent to the one above.

In order to make contact with the language of ?-algebras, we now introduce the so-called Weyl algebra. This is the algebra of formal products of objects U(*a*) and V(*b*), subjected to the Weyl relations (7). Note that, thanks to the Weyl relations, a generic element of the Weyl algebra can always be written as a finite linear combination of elements of the form U(*a*)V(*b*), *a*, *b* ∈ R, or equivalently of elements

$$\mathcal{W}(a,b) = e^{iab/2} \mathcal{U}(a) \mathcal{V}(b), \tag{10}$$

known as Weyl operators.

A representation of the Weyl relations is, thus, tantamount to a representation of the Weyl algebra, and, since this is a ?-algebra with identity, its representations can be discussed in terms of states of the algebra2 . In particular, given a state, one can construct a cyclic representation of the algebra, by means of the so-called GNS construction [12]. Taking into account the above remarks, it follows that states of the Weyl algebra, and therefore the corresponding representations, are uniquely determined by the values assigned to the Weyl operators.

Concerning the unitary implementation of automorphisms of ?-algebras, it is a well known fact that, if a state *ω* is invariant under a given transformation, then a unitary implementation of that transformation is ensured to exist in the GNS representation defined by *ω*. The above constructions related to the Weyl algebra are straightforwardly generalized to any finite number of DoF, and also without difficulties to field theories. In this respect, let us consider, for instance, a scalar field in R3 .

The starting point in the canonical quantization process is the choice of properly defined variables in phase space, that are going to play, e.g., the role of the *q*'s and the *p*'s. The integration of the field *φ* and of the canonically conjugate momentum *π* against smooth and fast decaying *test* functions provides just those amenable variables. Thus, given test functions *f* and *g* with the above properties, hence, belonging to the so-called Schwartz space S, one defines linear functions in phase space by

$$\Phi(\phi, \pi) \mapsto \int \phi f \, d^3 \mathbf{x} + \int \pi g \, d^3 \mathbf{x} =: \phi(f) + \pi(g). \tag{11}$$

In particular, the variables *φ*(*f*) and *π*(*g*), with Poisson bracket

$$\{\Phi(f), \pi(\mathbf{g})\} = \int f \mathbf{g} \, d^3 \mathbf{x} \tag{12}$$

replace in this context the familiar *q*'s and *p*'s. The corresponding Weyl relations are

$$\mathcal{V}(\mathbf{g})\mathcal{U}(f) = e^{i\int f\,\mathbf{g} d^3 \mathbf{x}} \mathcal{U}(f)\mathcal{V}(\mathbf{g}), \qquad f, \mathbf{g} \in \mathcal{S}, \tag{13}$$

with seemingly defined Weyl operators

$$\mathcal{W}(f,\mathfrak{g}) = e^{\frac{i}{2}\int fg d^3x} \mathcal{U}(f)\mathcal{V}(\mathfrak{g}).\tag{14}$$

Let us focus on a particular type of representations of the Weyl relations (or equivalently of the associated Weyl algebra), namely representations of the Fock type. These are representations defined by complex structures on the phase space or equivalently on the space S ⊕ S of pairs of test functions (*f* , *g*). In this respect, the pairs (*f* , *g*) define linear functionals in the phase space, and thus S ⊕ S is naturally dual to the phase space inheriting, therefore, a symplectic structure that is induced from that originally considered on phase space. We also recall that a complex structure *J* on a linear space with symplectic form Ω is a linear symplectic transformation such that *J* 2 = −**1** and is compatible with Ω in the sense that the bilinear form defined by Ω(*J*·, ·) is positive definite.

Let *J* be a complex structure on the symplectic space S ⊕ S. As mentioned above, a state of the Weyl algebra is defined by the values on the Weyl operators, and, therefore, the assignment

$$\mathcal{W}(f,\emptyset) \mapsto e^{-\frac{1}{4}\Omega(f(f,\emptyset),(f,\emptyset))}\tag{15}$$

defines a state and an associated cyclic representation of the Weyl algebra. The cyclic vector is, here, physically interpreted as the vacuum of the Fock representation, and therefore the above expression coincides precisely with the expectation values of the Weyl operators on the vacuum of the Fock representation defined by *J*.

Let us now discuss the question of unitary implementation of symplectic transformations in the specific context of Fock representations. We consider unitary operators W(*f* , *g*) providing a representation of the Weyl algebra and a linear canonical transformation *A*. There is, then, a new representation W*A*, defined (in the same Hilbert space) by W*A*(*f* , *g*) := W(*A* −1 (*f* , *g*)). If the representation W is defined by a complex structure *J*, it follows that W*A* corresponds to a new complex structure *JA* := *AJA*−1 . In general, *A* is not unitarily implementable, i.e., there is no unitary operator *UA* such that

$$\mathcal{U}\_A^{-1} \mathcal{W}(f, \mathcal{g}) \mathcal{U}\_A = \mathcal{W}(A^{-1}(f, \mathcal{g})). \tag{16}$$

In fact, the Fock representations defined by *J* and *JA* are unitarily equivalent if and only if the difference *JA* − *J* is an operator of a special type, namely a Hilbert–Schmidt operator [13]. Two notorious cases where that condition is automatically satisfied are the following. First, every operator in a finite dimensional linear space is of the Hilbert–Schmidt type, and therefore the unitary implementation of symplectic transformations comes for free in finite dimensional phase spaces, as expected from the Stone–von Neumann theorem. On the other hand, the null operator is always of the Hilbert–Schmidt type, regardless of the dimensionality, and therefore a transformation *A* that leaves *J* invariant is always unitarily implementable *in the Fock representation defined by J* .

With respect to previous remarks in our exposition, we note that invariance of *J* immediately translates into invariance of the associated Fock state defined by (15). Of course, the two situations that we have described correspond only to sufficient conditions for unitary implementation, which are by no means necessary. In particular, unitary implementation of a canonical transformation *A* can be achieved via a non-invariant complex structure *J*, provided that *JA* − *J* is Hilbert–Schmidt.

In any case, and in clear contrast with the situation found for a finite number of DoF, in field theory, no Fock representation supports the unitary implementation of the full group of linear canonical transformations. Fock representations are, therefore, distinguished by the class of transformations that are unitarily implementable. Now, in a particular theory, specified, e.g., by a given Hamiltonian, a particular set of canonical transformations stands out, namely transformations generated by the Hamiltonian and by possible symmetries.

The requirement of unitary implementation of relevant canonical transformations, therefore, provides a criterion guiding the selection of one representation over another when we are trying to quantize a field theory. In this respect, we notice that the case of the set of transformations corresponding to classical time evolution is particularly relevant, given the role that unitarity plays in the probabilistic interpretation of the quantum theory.

The simplest situation is that of a free field of mass *m* in Minkowski spacetime. In this case, the representation is completely fixed by the requirement of invariance under spatial symmetries and time evolution (or the full Poincaré group), in the sense that a unique complex structure is selected by that requirement, namely *Jm* defined by

$$J\_m(f, \mathcal{g}) = \left( (m^2 - \Delta)^{\frac{1}{2}} \mathcal{g}, -(m^2 - \Delta)^{-\frac{1}{2}} f \right),\tag{17}$$

where ∆ is the Laplacian.

In addition to the above free field in Minkowski spacetime, there are other known situations of linear dynamics where uniqueness results apply. In fact, provided that the Hamiltonian is time independent, the criterion of positivity of the energy, together with invariance under the 1-parameter group of canonical transformations generated by the Hamiltonian, is sufficient to select a unique complex structure.

Here, positivity means that the unitary group implementing the dynamics possesses a positive generator, i.e., the quantum Hamiltonian is a positive operator [14,15]. This result finds remarkable applications in the quantization of free fields in stationary curved spacetimes (i.e., with a timelike Killing vector) [16,17]. On the other hand, no general uniqueness results are available for the non-stationary situations typical in cosmology.

#### **3. Loop Quantum Cosmology**

Let us now consider the representation of the Weyl relations used in LQC, sometimes referred to as the polymer representation. We will restrict our attention to its simplest version, namely the one associated with the homogeneous and isotropic flat FLRW model3 . For convenience, we set the speed of light and Newton constant multiplied by 4*π* equal to the unit, and we make the Immirzi parameter [18] equal to 3/2 in order to simplify our equations, without loss of generality.

At the classical level, the system is described by a pair of canonically conjugate variables, typically denoted by *c* and *p*, with Poisson bracket {*c*, *p*} = 1. The variables *c* and *p* parametrize, respectively, the (homogeneous) Ashtekar connection and the densitized triad (see [10] for details in the context of the current uniqueness discussion, and [19,20] for more general introductions to LQC). In particular, *p* is proportional to the square scale factor of the FLRW spacetime.

Let us consider the Hilbert space HP defined by the discrete measure in R, i.e., the space of complex functions *ψ*(*p*) such that

$$\sum\_{p \in \mathbb{R}} |\psi(p)|^2 < \infty,\tag{18}$$

with the inner product given by

$$
\langle \psi, \psi' \rangle = \sum\_{p \in \mathbb{R}} \bar{\psi}(p) \psi'(p),
\tag{19}
$$

where the overbar denotes complex conjugation. We note that this Hilbert space, also referred to as the polymer Hilbert space, is very different from the standard one, *L* 2 (R). In particular, HP is non-separable4 . An orthonormal basis of HP is formed, e.g., by the uncountable set of functions Ψ*p*0 , for all *p*0 ∈ R, where

$$\Psi\_{p\_0}(p) = \delta\_{pp\_0} \tag{20}$$

with *δpp*0 as the Kronecker delta.

We, then, define the operators UP(*a*) and VP(*b*), for *a*, *b* ∈ R, acting on HP by

$$(\mathcal{U}\_{\mathcal{P}}(a)\psi)(p) = \psi(p-a) \tag{21}$$

and

$$(\mathcal{V}\_{\mathcal{P}}(b)\psi)(p) = e^{ibp}\psi(p). \tag{22}$$

It is clear that these operators satisfy the Weyl relations (7) and that the representation is irreducible. The map *b* 7→ VP(*b*) is continuous, and thus one can define the infinitesimal generator. We will denote it by *π*(*p*) and not *p*ˆ, to distinguish it from the standard Schrödinger representation in *L* 2 (R). It follows that

$$(\pi(p)\psi)(p) = p\psi(p). \tag{23}$$

On the other hand, UP(*a*) is not continuous. To see this, it suffices to note that, for arbitrarily small *a*, the vector Ψ*p*0 is mapped by UP(*a*) to an orthogonal one Ψ*p*0+*a*. The would-be generator of the unitary group UP(*a*) cannot, therefore, be defined. Taking into account the commutator [*π*(*p*), UP(*a*)], one can see that the operators UP(*a*) can be regarded as providing a quantization of the classical variables *e iac*. We note that the reality conditions are properly satisfied, since U † (*a*) = U(−*a*), with the dagger denoting the adjoint. Thus, the set of operators UP(*a*), together with *π*(*p*), provide a quantization, in the usual Dirac's sense, of the Poisson algebra of phase space functions made of finite linear combinations of the functions *p* and *e iac*, with *a* ∈ R (see [23] for details).

This Poison algebra is, of course, different from the kinematical algebra usually associated with the linear phase space R2 with coordinates *c* and *p*, which is simply the Heisenberg algebra of linear (non-homogeneous) functions in *c* and *p*. At the foundation of LQC, there is, thus, a situation similar to that of LQG: a non-standard choice of basic variables and a representation of the associated algebra that is non-continuous (at least in part of the algebra), thus, obstructing the quantization of the connection itself [8]. However, we note that, contrary to the situation in LQG, the standard Schrödinger quantization gives us a different representation of the very same Poisson algebra, since the variables *e iac* are trivially quantized in *L* 2 (R) by multiplication operators, e.g., by operators *e*c*iac* such that

$$\left(\widehat{e^{iac}}\psi\right)(q) = e^{iaq}\psi(q). \tag{24}$$

In the case of this Schrödinger quantization, the continuity of the representation U(*a*) allows us to define the operator *c*ˆ itself, and therefore to extend the quantization of configuration variables to a much larger set of functions *f*(*c*), simply by defining *f*d(*c*) = *f*(*c*ˆ); whereas, in the polymer quantization, one is restricted to quantize configuration variables of the type *e iac* (and linear combinations thereof).

In LQG, there is a celebrated result about the uniqueness of the quantization that gives robustness to the loop representation [24,25]. In [10], the authors proved a similar uniqueness result for LQC, that we now discuss. In order to make contact with that work, which uses the language of ?-algebras and corresponding states, we first introduce the LQC analogue of the LQG holonomy-flux algebra, which is again denoted in [10] as the quantum holonomy-flux ?-algebra U.

The LQC ?-algebra U is constructed in [10] as the algebra of formal products of operators corresponding to the variables *p* and *e iac* subjected to the conditions coming from the commutator [*π*(*p*), UP(*a*)]. However, since we have already introduced the Weyl algebra, it is more natural here to follow an alternative procedure, and identify instead the ?-algebra U with the Weyl algebra, i.e., the algebra of formal products of objects U(*a*) and V(*b*), subjected to the Weyl relations (7).

Recall now that the group of spatial diffeomorphisms is a "gauge" symmetry in General Relativity (GR). Physical states in quantum gravity should, therefore, be invariant under the quantum operators representing these diffeomorphisms (or annihilated by the quantum diffemorphism constraint). Thus, a unitary implementation of the group of spatial diffeomorphisms in the quantum Hilbert space is required. Since the LQG holonomy-flux algebra is again a ?-algebra with identity, it follows from previous comments that, in order to achieve the required unitary implementation of the group of diffeomorphisms, it is sufficient that the quantization be defined by a diffeomorphism invariant state.

The LQG result of uniqueness of the quantization [24,25] guarantees, precisely, that there exists a unique diffeomorphism invariant state of the LQG holonomy-flux algebra. Moreover, the GNS representation defined by such a unique invariant state is unitarily equivalent to the LQG representation that was previously known. The analogous result for LQC starts from the observation, explained in detail in [10], that a residual gauge group still remains, when descending from full GR to homogeneous and isotropic flat models. In fact, although almost all of the diffeomorphism gauge symmetry is automatically fixed, in

homogeneous and isotropic flat models, one is left with a small gauge group, namely the group of isotropic dilations, acting on phase space as

$$(p, c) \mapsto (\lambda p, c/\lambda), \qquad \lambda \in \mathbb{R}.\tag{25}$$

Thus, there is the possibility of exploring this residual gauge symmetry in order to select a state in LQC, very much like in the above mentioned LQG uniqueness result. The authors of [10] succeeded in proving that there is a unique dilation invariant state of the LQC holonomy-flux algebra, for which the associated GNS representation is unitarily equivalent to the polymer representation described above.

Although the result of [10] is interesting and rigorous, we argue here that, in a certain sense, its status is not as strong as that of the LQG uniqueness result. From the physical viewpoint, what is really required is a unitary implementation of the group of interest (which may correspond to dynamics, to symmetries, or to gauge transformations), and not necessarily an invariant state. There are even situations (see, e.g., Section 4) where invariant states are not available, and nevertheless a unitary implementation of the relevant group exists.

It is true that constructing the quantization by means of an invariant state is sufficient to achieve unitary implementation—and it is perhaps the "Kings way" of doing it—however, it is by no means necessary. In the present case, the Schrödinger representation, although not (unitarily equivalent to a representation) defined by an invariant state, carries a unitary implementation of the group of dilations (25), which is physically just as good as the one provided by the polymer representation. The existence of this unitary implementation actually follows from the Stone–von Neumann theorem, but let us show it explicitly. We consider the transformations *Uλ* : *L* 2 (R) → *L* 2 (R) defined by

$$
\psi(q) \stackrel{\mathcal{U}\_{\lambda}}{\longmapsto} \lambda^{1/2} \psi(\lambda q), \qquad \lambda \in \mathbb{R}, \tag{26}
$$

which are clearly unitary ∀*λ*, with respect to the standard inner product defined by the measure *dq*. We consider also the standard operator *p*ˆ = −*i d dq* and the Schrödinger quantization of configuration variables

$$\left(\widehat{f(c)}\psi\right)(q) = f(q)\psi(q), \qquad \psi \in L^2(\mathbb{R}), \tag{27}$$

which clearly includes the LQC variables *e iac*. A straightforward computation shows that

$$\left(\mathcal{U}\_{\lambda}^{-1}\mathcal{Y}\mathcal{U}\_{\lambda}\psi\right)(q) = -i\lambda \frac{d}{dq}\psi(q),\tag{28}$$

$$\left(\mathcal{U}\_{\lambda}^{-1}\widehat{f(c)}\mathcal{U}\_{\lambda}\psi\right)(q) = f(q/\lambda)\psi(q),\tag{29}$$

or

$$
\mathcal{U}\_{\lambda}^{-1} \mathcal{Y} \mathcal{U}\_{\lambda} = \widehat{\lambda p}, \qquad \mathcal{U}\_{\lambda}^{-1} \widehat{f(c)} \mathcal{U}\_{\lambda} = \widehat{f(c/\lambda)}, \tag{30}
$$

which is the announced unitary implementation of the group of dilations (25) in the Schrödinger representation.

From this perspective, we conclude that the physical criterion of a unitary implementation of the residual group of dilations in homogeneous and isotropic flat cosmology does not fully succeed in selecting a unique quantization, since both the polymer and the Schrödinger representations are viable from this viewpoint. Only the more mathematically stringent requirement of strict invariance selects a unique state. This is perhaps a reminder about the fact that strict invariance is not an unavoidable requirement and that the quantization of groups of interest via unitary implementations that are not *necessarily* based on invariant states is worth exploring.

Let us end with a brief comment regarding the analogous uniqueness result in LQG. In that case, the uniqueness is also proved by requiring the strict invariance of a state

of the holonomy-flux algebra under the action of spatial diffeomorphisms. There is, however, a key difference with respect to the above described situation in LQC: no other (irreducible) representation is known, of any kind, admitting a unitary implementation of the diffeomorphism group5 , and thus there is no alternative route that may cast any shadow on the uniqueness.

The situation remains, however, somewhat open, until a stronger uniqueness result is demonstrated that is based exclusively on a unitary implementation of the spatial diffeomorphisms and not just in strict invariance or otherwise until a new representation of the LQG algebra admitting a unitary non-invariant implementation of the diffeomorphism group is constructed.

#### **4. Fock Quantization in Non-Stationary Cosmological Settings**

As we have already mentioned, the Stone–von Neumann uniqueness result fails for an infinite number of DoF, and no general result on the uniqueness of the quantization is available for non-stationary situations. This includes, of course, cases of interest in QC. In some of those cases, the underlying theory can be recast in the form of a linear scalar field with a time-dependent mass, propagating in an auxiliary background spacetime, which is both static and spatially compact.

One such situation is the linearly polarized Gowdy model with the spatial topology of a three-torus, where the gravitational degrees of freedom are encoded by a scalar field on *S* 1 , evolving in time precisely as a linear field with a time-dependent mass of the type *m*(*t*) = 1/(4*t* 2 ) [31,32]. Other situations of interest include free scalar fields in cosmological scenarios, e.g., propagating in (compact) FLRW or the Sitter spacetimes. In those cases, the non-stationarity is transferred from the background to the field efective mass, by means of a simple transformation.

#### *4.1. Gowdy Models*

Midisuperspace models are symmetry reductions of full GR that retain an infinite number of degrees of freedom. Typically, these are local degrees of freedom, and thus they often describe inhomogeneous scenarios. Therefore, these midisuperspace models must face the inherent ambiguity that affects the quantization of fields.

One of the simplest inhomogeneous cosmologies obtained with a symmetry reduction is the linearly polarized Gowdy model on the three-torus, *T* 3 [7]. This model describes vacuum spacetimes with spatial sections of *T* 3 -topology containing linearly polarized gravitational waves, with a symmetry group generated by two commuting, spacelike, and hypersurface orthogonal Killing vector fields. As a consequence, the local physical degrees of freedom can be parametrized by a scalar field corresponding to those waves and effectively living in *S* 1 .

After a partial gauge fixing, the line element of the linearly polarized Gowdy *T* 3 model can be written as [33]

$$ds^2 = e^{\gamma - \phi/\sqrt{p}} \left( -dt^2 + d\theta^2 \right) + e^{-\phi/\sqrt{p}} t^2 p^2 d\sigma^2 + e^{\phi/\sqrt{p}} d\delta^2,\tag{31}$$

where (*∂*/*∂σ*) *a* and (*∂*/*∂δ*) *a* are the two Killing vector fields. The true dynamical field DoF are encoded in *φ*(*θ*, *t*), where *t* > 0 and *θ* ∈ *S* 1 . On the other hand, *p* is a homogeneous non-dynamical variable, and the field *γ* is completely determined by *p* and *φ* as a result of the gauge-fixing process [33]. There remains a global spatial constraint on the system, giving rise to the symmetry group of (constant) translations in *S* 1 . The time evolution is dictated by the field equation

$$
\ddot{\phi} + \frac{1}{t}\dot{\phi} - \phi'' = 0.\tag{32}
$$

Here, the prime denotes the derivative with respect to *θ*, and the dot denotes the time derivative.

A complete quantization of the system was obtained in [34]. Nevertheless, it was soon realized [33,35] that the classical dynamics could not be implemented as a unitary transformation in such quantization. With the purpose of achieving unitarity, and restoring, in particular, the standard probabilistic interpretation of quantum physics within the quantum Gowdy model, an alternative quantization was introduced by Corichi, Cortez, and Mena Marugán [31,32]. A crucial step toward a quantization with unitary dynamics is the following time-dependent transformation performed at the level of the classical phase space. Instead of working with the original field *φ* and its corresponding conjugate momentum *Pφ*, the authors introduced a new canonical pair *χ* and *Pχ* related to the first one by means of the canonical transformation

$$\chi = \sqrt{t}\phi, \quad P\_{\chi} = \frac{1}{\sqrt{t}} \left( P\_{\phi} + \frac{\phi}{2} \right), \tag{33}$$

taking advantage in this way of the freedom in the scalar field parametrization of the metric of the Gowdy model.

The evolution of the new canonical pair was found to be governed by a time-dependent Hamiltonian that, in our system of units, adopts the expression

$$H\_{\chi} = \frac{1}{2} \oint d\theta \left[ P\_{\chi}^2 + \chi'^2 + \frac{\chi^2}{4t^2} \right]. \tag{34}$$

The corresponding Hamiltonian equations are

$$
\dot{\chi} = P\_{\chi'} \quad \dot{P}\_{\chi} = \chi'' - \frac{\chi}{4t^2},
\tag{35}
$$

that, combined, give the second-order field equation

$$
\ddot{\chi} - \chi'' + \frac{\chi}{4t^2} = 0.\tag{36}
$$

We can, therefore, view the system as a linear scalar field with a time-dependent mass of the form *m*(*t*) = 1/(4*t* 2 ), evolving in an effective static spacetime with one-dimensional spatial sections with the topology of the circle6 . Another relevant aspect of the model is the invariance of the dynamics under (constant) *S* 1 -translations:

$$T\_{\mathfrak{a}} : \theta \mapsto \theta + \mathfrak{a}, \qquad \mathfrak{a} \in \mathbb{S}^1. \tag{37}$$

These translations are, moreover, symmetries generated by a remaining constraint, as we already mentioned.

The quantization of the system put forward by Corichi, Cortez, and Mena Marugán in [31,32] starts from the CCRs satisfied by the canonical pair *χ* and *Pχ*, or rather by the corresponding Fourier modes. The advantage of using Fourier components, say *χn* = (1/√ 2*π*) H *dθe* −*inθχ*(*θ*), instead of the field *χ*(*θ*) itself is clear: since the spatial manifold is compact, the set of Fourier modes is discrete, and one, therefore, avoids the issues of dealing with operator valued distributions, such as, e.g., *χ*ˆ(*θ*). The aforementioned quantization is of the standard Fock type, with the following remarkable properties. To begin with, the complex structure on phase space that effectively defines the Fock representation is invariant under *S* 1 -translations.

Thus, the corresponding state of the Weyl algebra is *S* 1 -invariant, leading to a natural unitary implementation of these gauge transformations. Secondly and most importantly, the classical dynamics in phase space defined by Equation (35) (i.e., generated by the Hamiltonian (34)) are unitarily implemented at the quantum level. In other words, let *t*0 be an arbitrary but fixed initial time, and *S*(*t*, *t*0) be the linear symplectic transformation corresponding to the classical evolution in phase space from the time *t*0 to the arbitrary time *t*. Then, to each transformation *S*(*t*, *t*0), there corresponds a unitary operator *U*(*t*, *t*0), that

intertwines between the quantum operators *χ*ˆ and *P*ˆ *χ* defined at the initial time *t*0 and those obtained from *χ*ˆ and *P*ˆ *χ* by application of *S*(*t*, *t*0), as exemplified in Equations (5) and (6).

This last transformation corresponds to the usual evolution in the Heisenberg picture, which can always be formally defined, once canonical operators are given at some initial time. The key difference with respect to systems with a finite number of DoF is that, whereas, in those cases, the relation between the "initial" and the evolved operators in principle is always unitary, the existence of such unitary operators is far from being guaranteed in an arbitrary representation of the CCRs for field theory (or generally with an infinite number of DoF).

We note that, although a unitary implementation of all the transformations *S*(*t*, *t*0) is achieved in the Corichi, Cortez, and Mena Marugán representation, the corresponding state is not invariant under these transformations. In fact, no state exists such that it remains invariant under all the transformations *S*(*t*, *t*0), ∀*t*.

All in all, the quantization of the linearly polarized Gowdy model proposed in [31,32] is one of the few available examples of a rigorous and fully consistent quantization of an inhomogeneous cosmological model. Nonetheless, the eventual robustness of its physical predictions might be affected by the possible existence of major ambiguities in the quantization process. Fortunately, a quantization with the aforementioned properties is, indeed, unique as shown in [36,37], where the uniqueness result that we now discuss was derived.

A source of ambiguity in the process leading to the Corichi, Cortez, and Mena Marugán quantization is the choice of representation for the canonical pair *χ* and *Pχ*. However, it was shown in [36] that any other Fock representation of the CCRs that (*i*) is defined by a *S* 1 invariant complex structure (or equivalently by a *S* 1 -invariant state of the Weyl algebra) and (*ii*) allows a unitary implementation of the dynamics defined by Equation (35), is unitarily equivalent to the considered representation (and therefore physically indistinguishable). We note that there actually is an infinite number of *S* 1 -invariant states, leading to many inequivalent representations. It is only after the requirement of unitary dynamics that a unique unitary equivalence class of representations is selected.

Another possible source of ambiguity concerns the choice of the "preferred" canonical pair (*χ*, *Pχ*). In this respect, note that a time-dependent transformation different from Equation (33) would lead to classical dynamics that would not reproduce Equation (35), and it is, in principle, conceivable that the new dynamics could be unitarily implemented in a different representation, thus, leading to a distinct quantization. This is, however, not the case.

It was shown in [37] that any other canonical transformation of the type (33) modifies the equations of motion in such a way as to render impossible the unitary implementation of the dynamics, with respect to any Fock representation defined by a *S* 1 -invariant complex structure. The kind of transformations considered in [37] is restricted by the natural requirements of locality, linearity, and preservation of *S* 1 -invariance. In configuration space, these are contact transformations that produce a time-dependent scaling of the field *φ* (see [37] for a detailed discussion). Such scalings can always be completed into a canonical transformation in phase space, of the general form

$$(\phi, P\_{\Phi}) \mapsto \left( f(t)\phi, \frac{P\_{\Phi}}{f(t)} + g(t)\phi \right). \tag{38}$$

which includes a contribution to the new momentum, which is linear in the field *φ*.

Finally, let us mention that completely analogous results were obtained for the remaining linearly polarized Gowdy models, namely those with spatial sections with topologies *S* 1 × *S* 2 and *S* 3 . This analysis was performed in the following independent steps. First, the classical models were addressed in [38], showing that, in these cases, the local gravitational DoF can also be parametrized by a single scalar field, namely an axisymmetric field in *S* 2 . Then, following a procedure similar to the one introduced by Corichi, Cortez, and Mena Marugán, a Fock quantization with unitary dynamics was obtained [39]. In particular, a

time-dependent scaling of the original field is again involved, now of the form *φ* 7→ √ sin *tφ*. Finally, the uniqueness of the quantization obtained in this way was proved in [40].

#### *4.2. Quantum Field Theory in Cosmological Settings*

A common feature of the Gowdy models mentioned in the previous section is that the local DoF are parametrized by a scalar field effectively living in a compact spatial manifold. Moreover, after the crucial scaling of the field, the dynamics are those of a linear field with time-dependent mass, i.e., they obeys a second-order equation of the type

$$
\ddot{\chi} - \Delta \, \chi + s(t)\chi = 0,\tag{39}
$$

where ∆ is the Laplace–Beltrami (LB) operator for the spatial sections in question, e.g., *S* 1 for the Gowdy model on *T* 3 and *S* 2 for the remaining two models.

Remarkably, a whole different type of situations in cosmology can also be described by an equation of the form (39). Let us consider, e.g., a free scalar field in a homogeneous and isotropic FLRW spacetime, with the line element

$$ds^2 = a^2(t) \left[ -dt^2 + h\_{ab} dx^a dx^b \right],\tag{40}$$

where *t* is the conformal time, *a*(*t*) is the scale factor, and *hab* (*a*, *b* = 1, 2, 3) is the Riemannian metric, for either flat Euclidean space or the 3-sphere. A minimally coupled scalar field of mass *m* obeys, in this cosmological spacetime, the equation

$$
\ddot{\phi} + 2\frac{d}{a}\phi - \Delta\ \phi + m^2 a^2 \phi = 0,\tag{41}
$$

where ∆ is the LB operator defined by the metric *hab* of the spatial sections.

The most obvious situation described by this setup is the propagation of an actual (test) scalar matter field (disregarding the backreaction) in an FLRW background. Nonetheless, the treatment of quantum perturbations, both of matter and of gravitational DoF, also fits the above description. In fact, in the context of cosmological perturbations, the leadingorder approximation in the action, together with a neglected backreaction, amounts to keep the homogeneous classical cosmology as the background and treats both matter and gravitational perturbations as fields propagating on that background [41–45].

The quantum treatment of the situations described above, therefore, faces the ambiguity of the choice of quantum representation. In particular, the criterion based on stationarity, mentioned in Section 2, is not available, since all these situations are inherently non-stationary.

There is, however, a natural avenue to address this issue, arriving hopefully at a unique quantum theory, with unitary dynamics. The way forward is actually suggested by a standard procedure, commonly found precisely in QFT in curved spacetimes (see, e.g., [46,47]) and in the treatment of cosmological perturbations (see, for example, [41–44]). This consists in the scaling of the field variable *φ* by means of the scale factor, thus, introducing the rescaled field *χ* = *a*(*t*)*φ*, which now obeys an equation of the type (39), with *s*(*t*) = *m*2 *a* 2 − (*a*¨/*a*).

Thus, the combined effect of the use of conformal time and the scaling *φ* 7→ *a*(*t*)*φ* is to recast the field equation in the form (39), which is effectively the field equation of a linear field with a time-dependent quadratic potential *V*(*χ*) = *s*(*t*)*χ* 2/2, propagating in a static background with metric

$$ds^2 = -dt^2 + h\_{ab}d\mathbf{x}^a d\mathbf{x}^b.\tag{42}$$

In this manner, we see that a generalization of the uniqueness result obtained for the Gowdy models would provide a useful criterion to select a unique quantization for fields in non-stationary backgrounds, such as FLRW universes, typically with compact spatial sections with the topology of *S* 3 or *T* 3 , also allowing applications to the usual flat universe case.

Note that the restriction to spatial compactness is most convenient from the viewpoint of mathematical rigor, as, otherwise, infrared issues would plague the analysis. Nonetheless, the physical effects of the artificially imposed compactness, e.g., in the spatially flat case, should be irrelevant when the physical problem at hand does not involve arbitrarily larges scales, e.g., going beyond the Hubble radius. The case of continuous scales, which corresponds to non-compactness, can be reached in a suitable limit for flat universes, after completing all the demonstrations of uniqueness in the framework of compact spatial sections [48].

The desired generalization of the result about the uniqueness of the representation in the aforementioned context of test fields and perturbations on cosmological backgrounds was obtained in [49–53]. In the rest of this Section, we very briefly explain the implications of these results.

Let I × Σ be a globally hyperbolic spacetime, where I ∈ R is an interval and Σ is a compact Riemannian manifold of dimension *d* ≤ 3 (for cosmological applications, one can think of Σ as being either *S* 1 , *S* 2 , *S* 3 , or *T* 3 ). The spacetime is assumed to be static, with a metric given by Equation (42), where *hab* is the time-independent, Riemannian metric on Σ.

Consider a linear scalar field in I × Σ obeying a field equation of the type (39), where *s*(*t*) is essentially7 an arbitrary function, and ∆ is the LB operator associated with the metric *hab*. Note that any symmetry of this metric is transmitted to the LB operator and, therefore, to the equations of motion. Let us consider the Weyl algebra associated with the field *χ* (and of course its canonical conjugate momentum) and their Fock representations. Then,

Remarkably, it again follows that the rescaling of the field *φ* 7→ *a*(*t*)*φ* is quite rigid and uniquely determined: no other time-dependent scaling can lead to a(n invariant Fock) quantization with unitary dynamics, and thus no ambiguity remains in the choice of the preferred field configuration variable. Furthermore, the unitarity requirement essentially selects, as well, the canonical momentum field. Physically, this can interpreted as a unique splitting between the time-dependence assigned to the background and that corresponding to the evolution of the scalar field DoF.

Finally, extensions of these uniqueness results were obtained in several directions. First, analogous results were obtained for scalar fields in homogeneous backgrounds of the Bianchi I type [54]. These spacetimes are not isotropic, and therefore the conformal symmetry, which was a common characteristic of the previous cases (at least asymptotically for large frequencies) is no longer present in an obvious way. Even more remarkable are the extensions of the uniqueness results attained for fermions, since they mark the transition to a largely unexplored territory [55–57].

A recent account of these results8 can be found in [59].

#### **5. Conclusions**

We reviewed and discussed several results concerning the quantization of systems with relevant applications in cosmology. In particular, we have focused our attention on results ascertaining the uniqueness of the quantization process. This uniqueness is crucial in order to provide physical robustness to the eventual cosmological predictions of the quantum models in question, which would otherwise be fundamentally affected by ambiguities.

Starting with homogeneous models in this cosmological context, we discussed a uniqueness result from Engle, Hanusch, and Thiemann, concerning the representation of the Weyl relations commonly used in LQC. Turning to recent investigations carried

out by our group and collaborators, the quantization of the family of inhomogeneous cosmologies, known as the linearly polarized Gowdy models, was reviewed next. Crucial in this discussion is the requirement of unitary implementation (at the quantum level) of the dynamics. Together with invariance under spatial symmetries, the criterion of unitary dynamics proved very effective in the selection of a unique and physically meaningful quantization in a variety of situations.

A common general mathematical model embodying all these cases is that of a scalar field with a time-dependent mass (with arbitrary time dependence except for some very mild conditions), propagating in a static spacetime with compact spatial sections. The existence and uniqueness of a Fock quantization with unitary dynamics has been proven for such a general model, thus, providing unique quantizations in cosmological systems, ranging from quantum fields in FLRW backgrounds to the quantization of (perturbative) gravitational DoF.

In particular, it follows from our analysis and the aforementioned discussions that the quantization of linear transformations of physical interest via a unitary implementation does not necessarily require the existence of an invariant state, and that it is worthwhile pursuing alternatives that are not based on such invariance.

In fact, what is really required in physical terms is a unitary implementation and not necessarily an invariant state, i.e., unitary implementations via non-invariant states are still physically acceptable and cannot be discarded. Clearly, a proof on the uniqueness of the quantization based solely on the unitary implementation of those transformations, rather than on the existence of an invariant state, is a stronger result inasmuch as the requirement of unitarity is weaker than invariance. Whereas, in some circumstances, there are good reasons to restrict attention to invariant states (for instance, when there is a time-independent Hamiltonian giving rise to a group of transformations), in the general case of arbitrary transformations, the requirement of an invariant state is not so compelling and non-invariant states cannot be simply disregarded.

In this sense, invariance can actually be an overly rigid demand in certain situations, especially if the considered transformations provide a notion of evolution that is crucial to describe the dynamics. For instance, the unitary implementation of the evolution via (dynamical-)invariant states in typically non-stationary scenarios is of doubtful physical use and, in fact, seems hopeless, whereas physically viable states, still leading to unitary dynamics and ensuring uniqueness, are available and have direct cosmological applications. In this respect, we conclude with the following remark.

The uniqueness of the Fock quantization of the scalar field with time-dependent mass was obtained by restricting the attention to states that remain invariant under spatial symmetries. Although the removal of this restriction seems unlikely to lead to physically new representations with unitary dynamics, it is, nevertheless, an open possibility. More precisely, it remains to be disproved the existence of representations with unitary dynamics and a unitary implementation of the spatial symmetries, which are, nevertheless, not unitarily equivalent to the representation defined by an invariant vacuum.

Thus, a conceivable line of future research in this area is to consider also Fock representations that, while not possessing an invariant vacuum, still allow a unitary implementation of the spatial symmetries. The uniqueness result would be strengthened if those representations were shown to be equivalent to the previous one; otherwise, new and potentially interesting representations could emerge. Another possibility, with a clearly higher potential to produce physically inequivalent results, is to use polymer-inspired representations for the scalar field, instead of Fock representations.

**Author Contributions:** Conceptualization: J.C., G.A.M.M. and J.M.V. Original Draft Preparation: J.C., G.A.M.M. and J.M.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Spanish MINECO grant number FIS2017-86497-C2-2-P, the Spanish MICINN grant number PD2020-118159GB-C41, and the European COST (European Cooperation in Science and Technology) Action number CA16104 GWverse. J.M.V. is grateful for the

support given by research unit Fiber Materials and Environmental Technologies (FibEnTech-UBI), on the extent of the project reference UIDB/00195/2020, funded by the Fundação para a Ciência e a Tecnologia (FCT).

#### **Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors are grateful to B. Elizaga Navascués for discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Sketchof the Proof of Uniqueness of the Representation for the Scalar Field with Time-Dependent Mass in** *S* **1**

In this appendix, we consider the simplest example of a linear scalar field in a static spacetime with compact spatial sections, namely, the case where the spatial sections are 1-dimensional with the topology of the circle. The field equation reads

$$
\ddot{\chi} - \Delta \chi + s(t)\chi = 0,\tag{A1}
$$

where the mass term *s*(*t*) can be an (essentially) arbitrary function of time.

Since the spatial manifold is compact and the field equation is linear, a Fourier decomposition gives us a discrete set of independent modes:

$$\chi(\theta, t) = \frac{1}{\sqrt{2\pi}} \sum\_{n=-\infty}^{\infty} \chi\_n(t) e^{in\theta} = \frac{1}{\sqrt{\pi}} \sum\_{n=1}^{\infty} \left( q\_n \cos(n\theta) + \chi\_n \sin(n\theta) \right) + \frac{q\_0}{\sqrt{2\pi}}.\tag{A2}$$

The configuration space for the scalar field is then described by the set of real variables *qn*, *n* ≥ 0, and *xn*, *n* > 0, which are completely decoupled. For simplicity, we drop all the modes *xn* (which can be treated like their cosine counterparts *qn* for *n* > 0) and *q*0, and continue with the infinite set {*qn*, *n* > 0}. The variable *q*0 is dropped to avoid introducing a special treatment in the case *s*(*t*) = 0, since *n* = 0 corresponds to a zero frequency oscillator or a free particle instead of a regular harmonic oscillator. In any case, it describes a single degree of freedom, which cannot affect the considered matters of unitary implementation.

The equations of motion for the modes are

$$
\ddot{q}\_n + [n^2 + s(t)]q\_n = 0.\tag{A3}
$$

The corresponding Hamiltonian equations are

$$
\phi\_n = p\_{n\nu} \quad \phi\_n = -[n^2 + s(t)]q\_{n\nu} \tag{A4}
$$

where *pn* is the momentum canonically conjugate to *qn*, i.e., {*qn*, *pn* 0} = *δnn*0 .

There are many representations of the CCRs satisfied by the infinite set of pairs {(*qn*, *pn*), *n* > 0}, or of the associated Weyl algebra. For instance, every sequence {*µn*, *n* > 0} of (quasi-invariant9 ) *probability* measures in R gives a representation, since it defines a regular product measure in the set of all sequences (*q*1, *q*2, . . .), thus, providing a Schrödinger type of quantization in the Hilbert space of square integrable functions in the configuration space.

What is not available, however, is the straightforward generalization of the usual representation in finite dimensions obtained from the Lebesgue measure *dq*, since no mathematical sense can be made of the formal infinite product ∏ ∞ *n*=1 *dqn*. In such a context, Fock representations of the Weyl relations are given by (normalized) Gaussian measures, which still make perfect sense in infinite dimensions. Even after restricting our attention to Gaussian measures, there are endless possibilities, and many of them lead to inequivalent representations.

In the present case, a particularly important measure in our infinite dimensional configuration space is

$$d\mu = \prod\_{n=1}^{\infty} e^{-n\,q\_n^2} \frac{\sqrt{n}dq\_n}{\sqrt{\pi}}\,,\tag{A5}$$

which is associated with the quantum operators

*q*ˆ*n*Ψ = *qn*Ψ, *p*ˆ*n*Ψ = −*i ∂ ∂qn* Ψ + *inqn*Ψ. (A6)

This is, in fact, a particular realization of the Fock representation given by the complex structure *J*0, defined as in Equation (17) 10 with *m* equal to zero. Thus, it should be no surprise that this particular representation allows a unitary implementation of the free massless field dynamics (i.e., with *s*(*t*) = 0). The quantization is the most natural one for the massless field, since the wave functional Ψ = 1 is invariant under the unitary group *U*(*t*) implementing the dynamics. Putting it differently, there is a well-defined quantum Hamiltonian, and Ψ = 1 is the zero-energy state of the free massless field.

The same representation also allows a unitary implementation of the time-dependent mass case. To see this, let us introduce the annihilation and creation-like variables *an* and *a*¯*n*, defined by

$$a\_n = \frac{nq\_n + ip\_n}{\sqrt{2n}},\tag{A7}$$

which are precisely the ones associated with the complex structure *J*0. In particular, this means that *J*0 takes the diagonal matrix form diag(*i*, −*i*) when written in terms of the basis in phase space made of the pairs (*an*, *a*¯*n*). With respect to this basis of (complex) variables, the classical evolution (from time *t*0 to time *t*) determined by Equation (A4) is given by non-vanishing 2 × 2 matrix blocks of the form

$$\mathcal{U}\_{\boldsymbol{n}}(\boldsymbol{t},\boldsymbol{t}\_{0}) = \begin{pmatrix} \boldsymbol{\alpha}\_{\boldsymbol{n}}(\boldsymbol{t},\boldsymbol{t}\_{0}) & \boldsymbol{\beta}\_{\boldsymbol{n}}(\boldsymbol{t},\boldsymbol{t}\_{0})\\ \boldsymbol{\beta}\_{\boldsymbol{n}}(\boldsymbol{t},\boldsymbol{t}\_{0}) & \boldsymbol{\alpha}\_{\boldsymbol{n}}(\boldsymbol{t},\boldsymbol{t}\_{0}) \end{pmatrix} \tag{A8}$$

with

$$|a\_{\boldsymbol{\pi}}(t, t\_0)|^2 - |\beta\_{\boldsymbol{\pi}}(t, t\_0)|^2 = 1, \qquad \forall \boldsymbol{n} > 0, \,\forall t, t\_0. \tag{A9}$$

Note that, in the case *s*(*t*) = 0, precisely because of the way in which the variables *an* are defined, the above parameters *βn* are all vanishing (and the parameters *αn* are only phases).

For quite general functions *s*(*t*) in Equation (A4), the following condition is satisfied [49]:

$$\sum\_{n}^{\infty} |\beta\_{n}(t, t\_{0})|^{2} < \infty, \quad \forall t, t\_{0}. \tag{A10}$$

This is precisely the necessary and sufficient conditions for unitary implementability of the dynamics, in the *J*0 representation. Thus, there is a Fock representation, defined by a complex structure that remains invariant under the spatial symmetries, such that a unitary quantum dynamics can be achieved.

The proof that a quantization with the above characteristics is unique is as follows. To begin with, any other invariant (under spatial isometries) complex structure *J* is related to *J*0 by *J* = *KJ*0*K* −1 , where *K* is a symplectic transformation given by a block diagonal matrix with 2 × 2 blocks of the form

$$\left(\begin{array}{c} \kappa\_{\hbar} \\ \overline{\lambda}\_{\hbar} \end{array} \begin{array}{c} \lambda\_{\hbar} \\ \mathbb{R}\_{\hbar} \end{array} \right)\_{\prime} \qquad |\kappa\_{\hbar}|^{2} - |\lambda\_{\hbar}|^{2} = 1, \ \forall n > 0. \tag{A11}$$

Suppose now that the dynamics are unitary in the Fock representation defined by *J*. This is equivalent to the unitary implementability in the *J*0-representation of a modified dynamics, obtained precisely by applying the transformation *K* to the canonical transformations corresponding to time evolution. A simple computation shows that, for this modified dynamics, the coefficients *βn* in Equation (A8) are replaced with

$$\mathcal{B}\_{n}^{\vec{I}}(t, t\_{0}) = 2i\mathbb{k}\_{n}\lambda\_{n}Im[a\_{n}(t, t\_{0})] + (\mathbb{k}\_{n})^{2}\beta\_{n}(t, t\_{0}) - \lambda\_{n}^{2}\vec{\beta}\_{n}(t, t\_{0}) \tag{A12}$$

where the notation *Im* denotes the imaginary part. Then, it follows from the hypothesis of unitary dynamics that the following condition holds:

$$\sum\_{n}^{\infty} |\beta\_n^I(t, t\_0)|^2 < \infty, \qquad \forall t, t\_0. \tag{A13}$$

Now, a detailed asymptotic analysis [49,50] shows that condition (A13) implies that

$$\sum\_{n}^{\infty} |\lambda\_n|^2 < \infty. \tag{A14}$$

Given the relation between *J* and *J*0, this last condition guarantees precisely that the operator *J* − *J*0 is of the Hilbert–Schmidt type, i.e., that the Fock representations defined by *J* and *J*0 are unitarily equivalent.

#### **Notes**

#### **References**

**João Marto 1,2**

1 Departamento de Física, Universidade da Beira Interior, Rua Marquês D'Ávila e Bolama, 6200-001 Covilhã, Portugal<EMAIL_ADDRESS>
2 Centro de Matemática e Aplicações da Universidade da Beira Interior, Rua Marquês D'Ávila e Bolama, 6200-001 Covilhã, Portugal