# Generative deep-learning reveals collective variables of Fermionic systems

Raphaël-David Lasseri  Centre Borelli, ENS Paris-Saclay, 91190 Gif-sur-Yvette,
France  Magic LEMP, 94110 Arcueil, France  David Regnier<EMAIL_ADDRESS>CEA, DAM, DIF, 91297 Arpajon, France  Université Paris-Saclay, CEA,
Laboratoire Matière en Conditions Extrêmes, 91680 Bruyères-le-Châtel, France
Mikaël Frosini  CEA, DEN, IRESNE, DER, SPRC, 13108 Saint-Paul-lès-Durance,
France  Marc Verriere  Nuclear and Data Theory Group, Nuclear and Chemical
Science Division, Lawrence Livermore National Laboratory, Livermore, CA,
United States  Nicolas Schunck  Nuclear and Data Theory Group, Nuclear and
Chemical Science Division, Lawrence Livermore National Laboratory, Livermore,
CA, United States

Deep-Learning, Generative learning, Autoencoder, Quantum phase transition,
Nuclear structure, Nuclear deformation

Complex processes ranging from protein folding to nuclear fission often follow
a low-dimension reaction path parameterized in terms of a few collective
variables. In nuclear theory, variables related to the shape of the nuclear
density in a mean-field picture are key to describing the large amplitude
collective motion of the neutrons and protons. Exploring the adiabatic energy
landscape spanned by these degrees of freedom reveals the possible reaction
channels while simulating the dynamics in this reduced space yields their
respective probabilities. Unfortunately, this theoretical framework breaks
down whenever the systems encounters a quantum phase transition with respect
to the collective variables. Here we propose a generative-deep-learning
algorithm capable of building new collective variables highly representative
of a nuclear process while ensuring a differentiable mapping to its Fermionic
wave function. Within this collective space, the nucleus can evolve
continuously from one of its adiabatic quantum phase to the other at the price
of crossing a potential energy barrier. This approach applies to any Fermionic
system described by a single Slater determinant, which encompasses electronic
systems described within the density functional theory.

Describing the dynamics of Fermionic systems is key to understanding how
molecules react to a laser excitation [1], study quantum phase transitions in
ultracold gases [2] or simulate low energy heavy-ion collisions [3]. For most
mesoscopic systems, the exponential growth of the Hilbert space with the
number of particles involved makes an exact resolution of this quantum many-
body problem unreachable. Yet, various complex processes such as large
amplitude vibrations in atomic clusters or atomic nuclei emerge as collective
behaviors of their constituents [4]. In such cases, the system’s wave function
mostly remains within a low-dimensional sub-manifold of the complete phase-
space. A few collective variables are sufficient to parameterize its dynamics.
In this situation, a widespread theoretical approach consists in computing a
series of static low-energy Fermionic wave functions along the collective path
[5, 6, 7]. This yields the so-called adiabatic energy landscape that reveals
favored configurations or reaction channels. Going one step further, it is
possible to recover a classical or quantum mechanical equation of motion in
terms of the collective variables [8, 9]. Integrating this equation in time
gives access to the timescales of the phenomenon along with the probability of
populating different final states. The success of this theoretical framework
requires that (i) the collective subspace is representative of the phenomenon
under study (ii) the mapping from the collective variable to the quantum state
of the system is differentiable. However, computing a landscape of adiabatic
quantum states that minimizes the energy at a given collective coordinate
gives no guarantee that the latter condition is fulfilled.

Chemistry and biochemistry possess a large corpus of methods to build
collective variables describing the positions of the nuclei along a reaction
path [10, 11, 12, 13, 14]. A method such as meta-dynamics built on top of ab-
initio electronic states ensures the continuity of the positions of the nuclei
as a function of the collective variable [15, 6]. Yet brutal changes in the
electronic configurations may still appear and manifest, for instance, as
energy cusps at conical intersections [16]. A similar situation holds in
nuclear physics, where collective variables related to the intrinsic shape of
the nuclear density are commonly used to study low-energy structure, giant
resonances, fusion, and fission reactions [17]. Within such collective space,
the clusterization of some nucleons [18, 19] or the discontinuous changes in
the shape of a fissioning nucleus appear as first-order quantum phase
transitions (QPT) [20, 21]. In several studies, the adiabatic energy landscape
was pragmatically smoothed out in order to be able to simulate the dynamics
nonetheless [22, 23, 16, 24]. Another approach consists in going beyond the
adiabatic approximation at the price of greater complexity and numerical cost
[25, 26, 24, 27, 28, 29]. Recently, Lau et al. proposed a novel method that
avoids these QPT by locally modifying the collective path [30]. Unfortunately,
this approach requires additional mean-field calculations, and its success is
not guaranteed for all QPT.

In parallel with these efforts, the fast development of dimension reduction
algorithms in the field of artificial intelligence opens new alternatives to
the design of collective variables. In the context of molecular systems, small
dimensional representations of the positions of nuclei were obtained with the
help of genetic algorithms [31], principal component analysis, Isomaps,
Sketch-maps or diffusion maps [10] but also autoencoders and variational
autoencoders [32, 33]. An additional difficulty arises when the collective
variables directly map Fermionic wave functions that belong to a non-Euclidean
space and possess peculiar symmetries [34]. Yet both fields of neural networks
representation of Fermionic states [35, 36, 37, 38, 39] and neural networks
for Riemannian geometries [40, 41, 42] are quickly progressing.

Based on these achievements, we show that it is possible to build
differentiable collective variables directly mapping Slater determinants and
demonstrate its ability to better describe the dynamics close to a QPT in
nuclear physics.

##  I Quantum phase transition in the 16O atomic nucleus

We consider the case of the atomic nucleus modeled as an ensemble of point-
like neutrons and protons. Among the thousands of nuclear systems of interest,
we choose to demonstrate our method on a nucleus of 16O. Due to its closed
shell nature, its superfluidity is negligible, which makes it similar to a
large number of electronic systems encountered in chemistry and condensed
matter. The Hartree-Fock ground state of 16O, obtained with a Machleidt ab-
initio nucleon-nucleon interaction [43], possesses a spherical one-body
density. As standard in nuclear physics, we study the static response of this
nucleus to an external force field constraining its elongation [44]. Here, the
collective variable is the expectation value of the intrinsic quadrupole
moment $Q_{20}=\langle\hat{Q}_{20}\rangle$ of the nucleus. Solving the
constrained Hartree-Fock equation for 200 regularly spaced values of the
quadrupole moment yields the adiabatic collective path.

The Figure 1 illustrates this collective path in terms of the quadrupole and
octupole moments of the solutions along with their Hartree-Fock energy.

Figure 1: Left Panel: The quadrupole-constrained Hartree-Fock path of 16O
plotted in the ($Q_{20}$, $Q_{30}$, energy) space. The dark continuous lines
represent its projections on the ($Q_{20}$, Energy) and ($Q_{20}$,$Q_{30}$)
planes. The color of the path corresponds to its $Q_{30}$ values. The ground
state of the system is the spherical configuration A. The points B and C are
separated by a first-order QPT highlighted by dashed lines. The point D shows
a typical high quadrupole deformation. Right Panel: Local nucleonic densities
for the indicated points of the quadrupole-constrained Hartree-Fock path.

In spite of the continuous behavior exhibited by the energy as a function of
the quadrupole moment, a sudden transition occurs in the system, shifting
discontinuously from the pear-shaped configuration B to the axially-symmetric
configuration C. This abrupt change in the quantum state of the nucleus is a
first-order quantum phase transition. The Fubini-Study distance [45] between
the quantum states B and C is two orders of magnitude higher than typical
between other neighboring states. This phase transition can be understood as
having a geometric nature, where the order parameter is represented by the
nucleus’s octupole moment, denoted as $Q_{30}=\langle\hat{Q}_{30}\rangle$. In
addition to this first-order QPT, there are four other points in the
collective space where the derivative of the octupole moment also exhibits
discontinuities. These points correspond to second-order QPTs of the system,
where the first derivative of the quantum state with respect to the collective
variable $Q_{20}$ is discontinuous.

The presence of these QPTs in the collective space poses a significant
challenge when attempting to describe the system’s dynamics. The collective
inertia, which plays a crucial role in simulating large-amplitude collective
dynamics and is a key driver of spontaneous fission half-lives [46], becomes
ill-defined in the vicinity of the QPTs [9]. The practical approach to
regularize these inertia commonly employed in state-of-the-art fission
dynamics calculations would likely result in an uncontrolled overestimation of
the decay rate from the prolate metastable state at $Q_{20}\simeq 30$ fm2
towards the ground state.

##  II Building an alternative collective variable

In this study, we use a set of adiabatic $Q_{20}$-constrained quantum states
to learn a continuous representation of the nuclear collective space. The data
set represented by Fig. 1 consists of 200 quantum states for 16O. As usual in
machine learning, we split the data randomly into a 70% training set, 20%
validation set and 10% test set. The calibration of our dimension reduction
algorithm involves the training and validation sets only while its final
predictive power is estimated from the test set. As shown later, this sampling
has a low impact on the final results. Each data point stands for a tensor
product of one time-reversal-invariant Slater determinant for the 8 neutrons
and another one for the 8 protons. The single-particle orbitals are expanded
on a truncated basis of 330 eigenstates of an axially-symmetric harmonic
oscillator Hamiltonian. The coefficients of this expansion fully define a
nuclear state. Yet, this mapping is not bijective since a unitary transform of
these coefficients may parameterize the same Slater determinant. More
specifically, the classes of normalized Slater determinants which are equal up
to a phase, often named Slater rays, form a set that is isomorphic to the
Grassmannian manifold, a non-Euclidean geometry [47].

In order to leverage standard dimension reduction techniques on such objects,
we first need to establish a differential mapping from a vector of real
coordinates to the Slater determinants of our data set. In this paper, we
explore two alternatives choices for this mapping. A first approach
parameterizes the class of Slater determinants up to a phase by the elements
of the upper triangular part of their one-body density matrices expressed in
the harmonic oscillator basis. Within the one-body space considered, this
mapping requires 108240 real numbers for the parameterization of one 16O
state. It is simple to implement but does not yield a one-to-one mapping, as
some values for the density matrix may not represent any Slater determinant. A
more complex but more robust approach leverages the Thouless theorem to
parameterize our finite set of Slater determinants. In this case, we first
compute a pivotal Slater determinant as the Karcher mean of our training
dataset [48, 49]. This new state minimizes its average geodesic distance to
the training dataset in the Grassmannian manifold. We then parameterize any
Slater determinant relative to this pivotal point $|\phi_{0}\rangle$ as

$|\phi(N,\theta,Z)\rangle=Ne^{i\theta}\text{exp}\left(\sum_{mi}Z_{mi}a^{\dagger}_{m}a_{i}\right)|\phi_{0}\rangle.$
(1)

The positive real number $N$ and the angle $\theta$ parameterize the norm and
phase of the state while the real matrix elements $Z_{mi}$ contain the
information on the Slater ray itself. In most applications, only the $Z$
coordinates present interest as a phase and norm convention could be chosen
arbitrarily. The sum runs over particle (resp. hole) orbitals of the pivotal
Slater determinant for the index $m$ (resp. $i$). This method ensures a one-
to-one differentiable mapping as long as the pivotal point has a non-vanishing
overlap with all the Slater determinant in the set. Describing an 16O
elongated state requires two vectors of 1304 real parameters, one associated
with the neutron $Z$ matrix and the other one with the protons. In this study,
both the density matrix and Thouless representations lead to similar
conclusions. In the following we present only results obtained with the
Thouless representation.

The next step to building a collective variable consists in finding a low-
dimension representation of the training data set. To do so, we combine a
principal component analysis (PCA) with a variational autoencoder (VAE) in a
similar manner as in [32]. The former performs a nearly lossless
compression/decompression from the pairs of the neutrons and protons $Z$
coordinates to an intermediate space $\mathcal{S}_{\text{PCA}}$. Selecting the
first 20 principal components enables a root mean square error lower than
$10^{-7}$ for the $Z$ matrices in the training set. This linear transformation
already gives a significant dimension reduction at a relatively low numerical
cost. Finally, we train a variational auto-encoder to further reduce the
dimension from the training set in $\mathcal{S}_{\text{PCA}}$ to a final
latent space of dimension one. We specifically choose a VAE as opposed to an
autoencoder so to regularize the latent space and improve the generative
quality of the algorithm. The details of the VAE architecture selected and the
process of hyper-parameter tuning are reported in Ref.V.2. After training, the
PCA decomposition followed by the encoder of the VAE maps any Thouless
representation (1) into a single latent variable $\lambda$. In the same way,
any $\lambda$ value passed through the VAE decoder, the inverse PCA
transformation and the inverse of the Thouless transformation gives a well
defined Slater determinant of 16O. Finally, using differentiable activation
functions for the variational auto-encoder ensures the differentiable nature
of the encoder and decoder. This makes $\lambda$ a new differentiable
collective variable. Its ability to span at least the same states as the
previous $Q_{20}$ representation is directly related to the convergence of our
training step. Here, the root mean square reconstruction error on the energy
is below 0.3% of the ground state energy for the training, validation and test
sets. Similarly the root mean square errors on the quadrupole and octupole
moments are below 0.01% of their respective variation range. This demonstrates
the quality of the training.

##  III Deformation properties along the new collective variable

Figure 2 shows the Hartree-Fock energy and the quadrupole moment of 16O as a
function of the $\lambda$ collective variable for the quantum states of the
training set. The topology of the energy curve is very close to the $Q_{20}$
representation of Fig. 1 up to a change of sign in abscissa. The $\lambda$
coordinate decreases monotonously with $Q_{20}$ and can be interpreted in a
similar way as encoding an information about the elongation of the system. The
$\lambda$ representation clearly separates the peared-shape configuration $B$
from the axially-symmetric one $C$ even though they share the same quadrupole
moment. In addition, small gaps appear at points like $\lambda=-2,Q_{20}=45$
fm2 corresponding to the second order QPTs in Fig. 1. As a first result, the
$\lambda$ collective variable better discriminates the quantum phases of our
system.

Figure 2: Energy landscape in the latent space discovered by the machine
learning algorithm. Points are colored according to the $Q_{20}$ coordinate
revealing the monotonous evolution of $Q_{20}$ with $\lambda$. Arrows report
the configurations A, B, C, D of Fig. 1.

Going one step further, we can now connect in a differentiable manner the
adiabatic phases. For this, we choose a fine enough mesh in the latent space
spanned by the variable $\lambda$ and generate for each point of this mesh a
new quantum state. The numerical cost associated with this step is negligible
compared to the one required by a series of new Hartree-Fock calculations.
Figure 3 shows the one-body density and multipole moments along the new
collective coordinates.

Figure 3: Upper panel: Collective path generated by the machine learning
algorithm represented in the ($Q_{20},Q_{30}$) space. Lower panel: Local
nucleonic densities for the indicated points. The new configuration $O$
corresponds to a $\lambda$ value in between the ones of configurations $B$ and
$C$.

Their values continuously connect configurations B and C.

Finally, figure 4 shows the Hartree-Fock energy as a function of $\lambda$
obtained by (i) choosing a fine $\lambda$-grid (ii) generating quantum states
for each $\lambda$ of this grid. A potential barrier between configurations B
and C appears as a striking new feature. The height of this new barrier is
13.2 MeV which should significantly hinder a dynamical transition from the
prolate configurations to the ground state of the system. This value depends
on the detail of the many-body path from configuration B to C. Typically, we
found that training multiple times the same VAE with different initial weights
always yields a new potential barrier. Its height varies with a standard
deviation of 2.1 MeV (cf. V.3).

Figure 4: Energy landscape in the latent space discovered by the machine
learning algorithm. Points generated from the latent space connect the
configurations B and C lying in two different quantum phases and reveal a new
potential barrier.

##  IV Conclusion

In this study, we present a groundbreaking approach that harnesses the power
of generative machine learning for dimensional reduction techniques, enabling
the construction of a novel collective variable that accurately represents
nuclear processes. By ensuring a differentiable parametrization of quantum
states, our method generates continuous trajectories between distinct quantum
phases, a crucial aspect for estimating the energy required for quantum phase
transitions.

Our findings open up new avenues for exploring nuclear dynamics within a
comprehensive quantum framework, such as the time-dependent generator
coordinate method. This advancement holds the potential to predict half-lives
of metastable states and the relative probabilities of reaction outcomes.

Moreover, our approach is applicable to a broad range of Fermionic systems
represented by Slater determinant states, including the vast domain of
electronic systems described through density functional theory. The
adaptability and versatility of neural networks ensure that our method can be
extended to accommodate even more complex Fermionic states in the future.

In summary, this study paves the way for a deeper understanding of nuclear
processes, electronic systems, and Fermionic states, demonstrating the
remarkable potential of integrating machine learning techniques in the realm
of quantum physics.

## References

  * González [2020] L. González,  _Quantum Chemistry and Dynamics of Excited States: Methods and Applications_ (Wiley-Blackwell, Hoboken, NJ, 2020).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Simenel and Umar [2018] C. Simenel and A. S. Umar, Progress in Particle and Nuclear Physics 103, 19 (2018).
  * Dinh  _et al._ [2017] P. M. Dinh, P.-G. Reinhard, and E. Suraud, “Similarities and differences between nuclei and metal clusters,” in  _Nuclear Particle Correlations and Cluster Physics_ (World Scientific, 2017) pp. 3–27.
  * Manzhos and Carrington [2021] S. Manzhos and T. J. Carrington, Chem. Rev. 121, 10187 (2021).
  * Schlegel [2003] H. B. Schlegel, Journal of Computational Chemistry 24, 1514 (2003).
  * Schunck and Regnier [2022] N. Schunck and D. Regnier, Progress in Particle and Nuclear Physics , 103963 (2022).
  * Hollingsworth and Dror [2018] S. A. Hollingsworth and R. O. Dror, Neuron 99, 1129 (2018).
  * Verriere and Regnier [2020] M. Verriere and D. Regnier, Front. Phys. 8 (2020), 10.3389/fphy.2020.00233.
  * Rogal [2021] J. Rogal, Eur. Phys. J. B 94, 223 (2021).
  * Maeda  _et al._ [2015] S. Maeda, Y. Harabuchi, Y. Ono, T. Taketsugu, and K. Morokuma, International Journal of Quantum Chemistry 115, 258 (2015).
  * Stewart  _et al._ [1987] J. J. Stewart, L. P. Davis, and L. W. Burggraf, Journal of Computational Chemistry 8, 1117 (1987).
  * Denzel  _et al._ [2019] A. Denzel, B. Haasdonk, and J. Kästner, J. Phys. Chem. A 123, 9600 (2019).
  * Tiwary and Parrinello [2013] P. Tiwary and M. Parrinello, Phys. Rev. Lett. 111, 230602 (2013).
  * Bussi and Laio [2020] G. Bussi and A. Laio, Nat. Rev. Phys. 2, 200 (2020).
  * Evenhuis and Martínez [2011] C. Evenhuis and T. J. Martínez, J. Chem. Phys. 135, 224110 (2011).
  * Bender [2008] M. Bender, The European Physical Journal Special Topics 156, 217 (2008).
  * Girod and Schuck [2013] M. Girod and P. Schuck, Phys. Rev. Lett. 111, 132503 (2013).
  * Ebran  _et al._ [2020] J.-P. Ebran, M. Girod, E. Khan, R. D. Lasseri, and P. Schuck, Phys. Rev. C 102, 014305 (2020).
  * Dubray and Regnier [2012] N. Dubray and D. Regnier, Computer Physics Communications 183, 2035 (2012).
  * Regnier  _et al._ [2019] D. Regnier, N. Dubray, and N. Schunck, Phys. Rev. C 99, 024611 (2019).
  * Ho and Rabitz [1996] T.-S. Ho and H. Rabitz, J. Chem. Phys. 104, 2584 (1996).
  * Hurd  _et al._ [2010] P. Hurd, T. Cusati, and M. Persico, Journal of Computational Physics 229, 2109 (2010).
  * Galvão  _et al._ [2016] B. R. L. Galvão, V. C. Mota, and A. J. C. Varandas, Chemical Physics Letters 660, 55 (2016).
  * Curchod and Martínez [2018] B. F. E. Curchod and T. J. Martínez, Chem. Rev. 118, 3305 (2018).
  * Varandas [2009] A. J. C. Varandas, Chemical Physics Letters 471, 315 (2009).
  * Bernard  _et al._ [2011] R. Bernard, H. Goutte, D. Gogny, and W. Younes, Phys. Rev. C 84, 044308 (2011).
  * Dietrich  _et al._ [2010] K. Dietrich, J.-J. Niez, and J.-F. Berger, Int. J. Mod. Phys. E 19, 521 (2010).
  * Zhao  _et al._ [2019] J. Zhao, T. Nikšić, D. Vretenar, and S.-G. Zhou, Phys. Rev. C 99, 014618 (2019).
  * Lau  _et al._ [2022] N.-W. T. Lau, R. N. Bernard, and C. Simenel, Phys. Rev. C 105, 034617 (2022).
  * Ma and Dinner [2005] A. Ma and A. R. Dinner, J. Phys. Chem. B 109, 6769 (2005).
  * Baima  _et al._ [2022] J. Baima, A. M. Goryaeva, T. D. Swinburne, J.-B. Maillet, M. Nastar, and M.-C. Marinica, Phys. Chem. Chem. Phys. 24, 23152 (2022).
  * Belkacemi  _et al._ [2022] Z. Belkacemi, P. Gkeka, T. Lelièvre, and G. Stoltz, J. Chem. Theory Comput. 18, 59 (2022).
  * Verriere  _et al._ [2022] M. Verriere, N. Schunck, I. Kim, P. Marević, K. Quinlan, M. N. Ngo, D. Regnier, and R. D. Lasseri, Frontiers in Physics 10 (2022).
  * Carleo  _et al._ [2018] G. Carleo, Y. Nomura, and M. Imada, Nat Commun 9, 1 (2018).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019).
  * Schütt  _et al._ [2019] K. T. Schütt, M. Gastegger, A. Tkatchenko, K.-R. Müller, and R. J. Maurer, Nat Commun 10, 5024 (2019).
  * Companys Franzke  _et al._ [2022] M. Companys Franzke, A. Tichai, K. Hebeler, and A. Schwenk, Physics Letters B 830, 137101 (2022).
  * Cassella  _et al._ [2023] G. Cassella, H. Sutterud, S. Azadi, N. D. Drummond, D. Pfau, J. S. Spencer, and W. M. C. Foulkes, Phys. Rev. Lett. 130, 036401 (2023).
  * Asif  _et al._ [2021] N. A. Asif, Y. Sarker, R. K. Chakrabortty, M. J. Ryan, M. H. Ahamed, D. K. Saha, F. R. Badal, S. K. Das, M. F. Ali, S. I. Moyeen, M. R. Islam, and Z. Tasneem, IEEE Access 9, 60588 (2021).
  * Zhang  _et al._ [2018] J. Zhang, G. Zhu, R. W. Heath Jr., and K. Huang, arXiv:1808.02229 [cs, eess, math, stat] (2018).
  * Chakraborty  _et al._ [2022] R. Chakraborty, J. Bouza, J. H. Manton, and B. C. Vemuri, IEEE Transactions on Pattern Analysis and Machine Intelligence 44, 799 (2022).
  * Machleidt and Sammarruca [2016] R. Machleidt and F. Sammarruca, Phys. Scr. 91, 083007 (2016).
  * Bender  _et al._ [2003] M. Bender, P. H. Heenen, and P. G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * Carollo  _et al._ [2020] A. Carollo, D. Valenti, and B. Spagnolo, Physics Reports Geometry of Quantum Phase Transitions, 838, 1 (2020).
  * Sadhukhan [2020] J. Sadhukhan, Front. Phys. 8, 567171 (2020).
  * Aoto and da Silva [2020] Y. A. Aoto and M. F. da Silva, Phys. Rev. A 102, 052803 (2020).
  * Hüper  _et al._ [2013] K. Hüper, M. Kleinsteuber, and H. Shen, Signal Processing 93, 459 (2013).
  * Townsend  _et al._ [2016] J. Townsend, N. Koep, and S. Weichwald, Journal of Machine Learning Research 17, 1–5 (2016).
  * Kingma and Welling [2019] D. P. Kingma and M. Welling, MAL 12, 307 (2019).
  * Snoek  _et al._ [2012] J. Snoek, H. Larochelle, and R. P. Adams, in  _Advances in Neural Information Processing Systems_, Vol. 25, edited by F. Pereira, C. Burges, L. Bottou, and K. Weinberger (Curran Associates, Inc., 2012).
  * Smith and Topin [2019] L. N. Smith and N. Topin,  _Artificial Intelligence and Machine Learning for Multi-Domain Operations Applications_ ,  11006, 1100612 (2019).
  * Lasseri  _et al._ [2020] R.-D. Lasseri, D. Regnier, J.-P. Ebran, and A. Penon, Phys. Rev. Lett. 124, 162502 (2020).
  * Schran  _et al._ [2020] C. Schran, K. Brezina, and O. Marsalek, The Journal of Chemical Physics 153 (2020), 10.1063/5.0016004, 104105.

##  V Methods

###  V.1 Dataset Generation and Splitting

We generated an initial dataset consisting of 200 pairs of Slater
determinants, which we randomly partitioned into training, validation, and
test sets. Each Slater determinant is represented by its $(330\times 330)$
one-body density matrix , generated using a Hartree-Fock code with an ab-
initio interaction. Each data point corresponds to an Oxygen-16 isotope,
axially constrained by the $Q_{20}$ collective variable. To check the
stability of the results to the sampling of the training, validation, and
testing sets, we performed 5-fold cross-validation. Each fold correspond to a
different sampling of these data sets while keeping the partition 70% for the
training set, 20% for the validation set and 10% for the test set. Table 1
presents the root mean square (RMS) of the Hartree-Fock energy for each fold.

RMS$(E_{HF})$ (keV) | train (70%) | valid (20%) | test (10%)
---|---|---|---
Fold 1 | 88 | 291 | 113
Fold 2 | 83 | 287 | 108
Fold 3 | 95 | 298 | 131
Fold 4 | 101 | 299 | 102
Fold 5 | 98 | 308 | 122
Table 1: Root mean square error of the Hartree-Fock energy for 5 different
folds of the dataset, for a fixed partition of the train/validation/test sets.
These values are all below 0.3% of the ground state energy.

The stability of the method across different dataset splits indicates that our
training procedure and evaluation of generalization are suitable for this
problem.

###  V.2 Variational Autoencoder implementation, architecture and training:

The non-linear transformation involved in building the new collective variable
$\lambda$ relies on a variational autoencoder (VAE) [50]. This type of neural
generative model learns a probabilistic mapping between an input data space
and a lower-dimensional latent space. The VAE is composed of two main
components: an encoder, which maps the input data to a latent space, and a
decoder, which reconstructs the original data from the latent space
representation. The primary objective of a VAE is to learn a continuous and
smooth latent space, enabling the generation of new data points by sampling
the latent space. The objective function of the VAE is the sum of a
reconstruction loss and a regularization term, typically the Kullback-Leibler
(KL) divergence, which encourages the latent space to have a specific
structure, such as following a Gaussian distribution. The total loss of the
VAE can be written as a sum of two terms:

$\mathcal{L}(\mu,\sigma^{2},x,\hat{x})=\mathcal{L}_{\text{rec.}}(x,\hat{x})+\mathcal{L}_{\text{KL}}(\mu,\sigma^{2}).$
(2)

The reconstruction loss $\mathcal{L}_{\text{rec.}}(x,\hat{x})$ measures the
difference between the original input data $x$ and the reconstructed data
$\hat{x}$ generated by the decoder. Assuming Gaussian likelihood with a
variance of one in each $x$ dimension the reconstruction loss for a single
data point $x$ and its reconstructed output $\hat{x}$ can be expressed as:

$\mathcal{L}_{\text{rec.}}(x,\hat{x})=\frac{1}{2}\sum_{i=1}^{D}\left(x_{i}-\hat{x}_{i}\right)^{2}$
(3)

where $D$ is the dimension of the data.

The KL divergence term of the loss measures the difference between the
approximate posterior distribution $q_{\phi}(z|x)$, represented by the
encoder’s output, and the true prior distribution $p(z)$, typically assumed to
be a standard normal distribution, i.e., $\mathcal{N}(0,I)$. The KL divergence
for a single data point can be expressed as:

$\mathcal{L}_{\text{KL}}(\mu,\sigma^{2})=-\frac{1}{2}\sum_{i=1}^{K}\left(1+\log(\sigma_{i}^{2})-\mu_{i}^{2}-\sigma_{i}^{2}\right),$
(4)

where $K$ is the dimension of the latent space, $\mu_{i}$ is the mean of the
approximate posterior distribution for the $i$-th latent variable, and
$\sigma_{i}^{2}$ is its variance.

Our VAE was implemented using PyTorch and the training process was completed
in around 10k epochs, taking approximately 20 minutes on a single RTX 3090
GPU. Early stopping was implemented based on the validation set.

The VAE architecture was optimized through Bayesian optimization of its
hyperparameters [51], including the number of hidden layers, density of each
layer, and the optimizer. To assess the quality of the reconstruction, the
mean squared error of the energy on the validation set was used. To accelerate
the training procedure and consequently reduce the time footprint of the
Bayesian optimization algorithm, we employed an implementation of the
superconvergence method [52] to accelerate the training. Superconvergence is a
phenomenon observed in the training of deep learning models, where the
convergence rate during optimization is significantly accelerated compared to
traditional training methods. This rapid convergence is achieved by tuning
hyperparameters such as the learning rate, batch size, and weight decay, along
with employing techniques like cyclical learning rate schedules, learning rate
warm-up, and regularization methods.

In our case, we used a cyclical learning rate, which adjusts the learning rate
periodically throughout training, alternating between higher and lower values.
This approach allows the model to escape local minima and explore the loss
landscape more effectively, leading to faster convergence and improved
generalization. We also used the learning rate warm-up technique, where the
learning rate is gradually increased during the initial phase of training,
allowing the model to adapt more smoothly to the optimization process. No
additional regularization was employed.

The final architecture consists of symmetric Encoder and Decoder, with their
weights trained independently from each other, and four hidden layers each.

  * •

Encoder:

    1. 1.

Input layer: 20 neurons

    2. 2.

Hidden layer 1: 100 neurons

    3. 3.

Hidden layer 2: 80 neurons

    4. 4.

Hidden layer 3: 60 neurons

    5. 5.

Hidden layer 4: 10 neurons

    6. 6.

Latent space layer ($\lambda$): 1 neuron

  * •

Decoder:

    1. 1.

Hidden layer 1: 10 neurons

    2. 2.

Hidden layer 2: 60 neurons

    3. 3.

Hidden layer 3: 80 neurons

    4. 4.

Hidden layer 4: 100 neurons

    5. 5.

Output layer: 20 neurons

###  V.3 Uncertainty of the potential barrier height

One way to check the robustness of our approach is to use the concept of
committees of neural networks. The idea described for instance in [53] and
[54] is to clone several time the neural model. The clones are identical in
terms of architecture and the only difference in between the members of this
neural committee consists on the initial weights initialization. Theses
differences in the initial weights allow the probing of the variational space
explored during the training phase for a fixed dataset. Each model of the
committee will then have a prediction on its own. If the modelisation is
robust enough the predictions will be relatively close to each other, if it is
not the case one can see this result as an hint that the specific datapoint
predicted is an outlayer. Moreover in the first case the spread of the
predictions can be directly linked to the uncertainty of the prediction and as
such is a good tool to estimate the reliability of the model. Our specific
case is quite interesting: from a physical point of view the topology of the
latent space does not matter at all since it describe a non-observable
representation. The only characteristic that can be used to quantify the
performance of our approach is the one which is directly linked to an
experimental observable e.g. the height of the energy barrier predicted. As
such by using a committee of VAE, one can probe the sensitivity and robustness
of our model. In our case for a committee of 10 VAEs we can estimate the
height of the barrier by computing the mean and the variance over the
different VAE, obtaining $E_{barrier}=13.21\pm 2.10$ MeV as such the results
are coherent with each other and reinforce the robustness claim of our
approach.

Figure 5: For a committee of 10 VAEs, (7 displayed here) the sampling of the
latent space at the level of the discontinuity. The presence of a new
potential barrier is a robust feature while the height this barrier varies
from a member of the committee to another. The standard deviation of this
variation is 2.10 MeV which represents 16% of its mean value.

##  VI Acknowledgements

The authors would like to thank S. Hilaire for his contribution in testing
this approach with a relativistic energy density functional structure code,
J.-P. Ebran for his insight on relativist energy density functional formalism,
N. Dubray for his support with the HFB3 solver and A. Penon for his advices on
training the VAE. This work was partly performed under the auspices of the
U.S. Department of Energy by Lawrence Livermore National Laboratory under
Contract No. DE-AC52-07NA27344 and by the Office of LDRD.

##  VII Author contributions

R.-D. Lasseri and D. Regnier developed and implemented the core ideas to build
new collective variables from a set of existing Slater determinants. M.
Frosini contributed to the production of the ab-inito Hartree-Fock series of
16O states. M. Verriere and N. Schunck provided technical consulting on the
parameteization of Slater determinant states for use with neural networks. D.
Regnier and R. Lasseri wrote the first version of this article that was
further read, amended, corrected and approved by all authors.


# Cosmological implications of generalized holonomy corrections

Cyril Renevey  Laboratoire de Physique Subatomique et de Cosmologie,
Université Grenoble-Alpes, CNRS/IN2P3
53, avenue des Martyrs, 38026 Grenoble cedex, France  Killian Martineau
Laboratoire de Physique Subatomique et de Cosmologie, Université Grenoble-
Alpes, CNRS/IN2P3
53, avenue des Martyrs, 38026 Grenoble cedex, France  Aurélien Barrau
Laboratoire de Physique Subatomique et de Cosmologie, Université Grenoble-
Alpes, CNRS/IN2P3
53, avenue des Martyrs, 38026 Grenoble cedex, France

###### Abstract

Most of the phenomenology of loop quantum gravity in the cosmological sector
is based on the so-called holonomy correction to the Hamiltonian constraint.
It straightforwardly modifies the Friedmann equations. In this work, we
investigate the influence of corrections generalizing the one usually used in
loop quantum cosmology. We find that a long enough inflation phase can be
generated by purely quantum geometrical effects but we also underline the
limitations of this scenario. In addition, we study the effects of generalized
holonomy corrections on an inflationary phase generated by a massive scalar
field. At the level of perturbations, we investigate in detail the
consequences on the primordial scalar power spectrum. The results are actually
quite general and can be used beyond the “loop” framework.

##  I Introduction

Loop quantum gravity (LQG) is a nonperturbative framework Ashtekar and Bianchi
(2021) providing a tentative quantization of general relativity (GR). It has
been expressed both in the canonical form Rovelli (2004) and in a covariant
way Rovelli and Vidotto (2014). As for all speculative theories the challenge
is twofold. On the one hand, one has to check internal consistency. This is
far from being a trivial requirement, especially in quantum gravity. From
gauge issues to infrared corrections, quite a lot of questions remain – at
least partially – open (see, e.g., Kiefer (2007)). On the other hand, it is
mandatory to face the outstanding question of phenomenological consequences
Barrau et al. (2014). Making links with observations is the key missing
ingredient for all attempts in quantum gravity, including string theory
Quevedo (2016).

In this article, we address the question of the robustness of some predictions
of LQG in the cosmological sector. Many different aspects have already been
investigated, taking into account in particular (see, e.g. Ashtekar et al.
(2008); Ashtekar and Singh (2011); Agullo et al. (2012); Diener et al. (2014);
Ashtekar and Barrau (2015); Bolliet et al. (2015); Alesci and Cianfrani
(2015); Martineau et al. (2017); Gielen and Oriti (2018))

  * •

the way initial conditions are set,

  * •

the validity of the minisuperspace approximation,

  * •

the backreaction effects,

  * •

the deformation of the algebra of constraints,

  * •

the inclusion of shear and curvature,

  * •

the quantization of operators associated with negative powers of the volume
operator,

  * •

the inclusion of effects inferred from quantum reduced loop gravity or group
field theory,

  * •

numerical results beyond the semiclassical approximation, etc.

Here, we tackle a different and somehow underestimated question: the
consequences of a generalized holonomy correction. The point is not to invent
what would be a superexotic theory, with new free parameters, to boost the
phenomenological richness. Just the other way round, the aim is to investigate
how reliable are the predictions made so far, taking into account implicit
assumptions that went mostly unnoticed and may play an important role.

The issue of quantization ambiguities in this framework was pointed out in
Perez (2006). Those associated with the quantization of the connection-based
holonomy variable might deeply influence the dynamics and constitute the
subject of this article. New theoretical arguments are being given in Amadei
et al. (2022), while the present work focuses on potential observable effects.
The question is especially important and meaningful when considered from a
renormalization point of view.

In the following, the basics of loop quantum cosmology (LQC) are first briefly
reminded. We then go into the details of generalized holonomy corrections. In
the next section, we show that a long period of inflation can be generated
using only a modified holonomy correction, without any massive scalar field.
We also highlight the limits of such a model. The consequences of generalized
holonomy corrections on both the background inflationary dynamics generated by
a massive inflaton field and the scalar primordial power spectrum are finally
exposed.

##  II FLRW loop quantum cosmology

In order to set the notations and remind the basics to the unfamiliar reader,
we summarize the main ideas behind LQC. This also allows the article to be
self-contained. In the fully constrained Ashtekar-Barbero formulation of GR,
the canonical variables are

$\displaystyle A^{i}_{a}\equiv\Gamma^{i}_{a}+\gamma
K_{a}^{i}\quad\textrm{and}\quad
E_{i}^{a}\equiv\frac{1}{2}\varepsilon^{abc}\varepsilon_{ijk}e^{j}_{b}e^{k}_{c}~{},$
(1)

where $\Gamma_{a}^{i}$ is the $su(2)$ spin connection, $\gamma$ the Barbero-
Immirzi parameter, $K_{a}^{i}$ the extrinsic curvature, $\varepsilon^{abc}$
the totally anti-symmetric tensor, and $e^{i}_{a}$ the triads. In this work,
we use $a,b,c$ as spacetime indices and $i,j,k$ as internal $su(2)$ algebra
indices. Both sets run from 1 to 3. These canonical variables satisfy the
relations

$\displaystyle\left\\{A^{i}_{a}(x),E^{b}_{j}(y)\right\\}=\kappa\gamma\delta^{i}_{j}\delta^{b}_{a}\delta^{3}(x-y)~{},$
(2)

with $\kappa=8\pi G$.

As GR is a fully constrained theory, its Hamiltonian is written as a sum of
constraints, respectively called scalar, vector and Gauss constraints. When
homogeneity is assumed, the scalar contribution $C_{g}$ is the only one to
remain. The Hamiltonian can then be written using the lapse function $N$ as

$\displaystyle C_{g}^{N}$ $\displaystyle=\int_{\Sigma}\differential
x^{3}NC_{g}$ (3) $\displaystyle=\frac{1}{2\kappa}\int_{\Sigma}\differential
x^{3}\frac{N}{\sqrt{q}}E^{a}_{i}E^{b}_{j}\left(\varepsilon^{ij}_{k}F^{k}_{ab}-2(1+\gamma^{2})K^{i}_{[a}K^{j}_{b]}\right)~{},$
(4)

with $\Sigma$ a compact hypersurface. In a homogeneous, isotropic, and flat
space, the metric reduces to the form

$\displaystyle\differential s^{2}=-\differential
t^{2}+a^{2}(t)\delta_{ab}\differential x^{a}\differential x^{b},$ (5)

where the cosmic time is related to the $0$-coordinate by $\differential
t=N\differential x^{0}$. To avoid divergent integrals and an ill-defined
symplectic geometry, we perform the integration on an arbitrary cubic fiducial
cell of comoving volume $V_{0}$. As the homogeneity assumption also implies
that spatial derivatives vanish, the spin connection disappears and the
canonical variables become simply

$\displaystyle
A^{i}_{a}(t)=\gamma\dot{a}(t)\delta^{i}_{a}\equiv\frac{c(t)}{V_{0}^{1/3}}\delta^{i}_{a}~{},$
(6)

where the dot represents a derivative with respect to the cosmic time $t$, and

$\displaystyle
E^{a}_{i}(t)=a^{2}(t)\delta^{a}_{i}\equiv\frac{p(t)}{V_{0}^{2/3}}\delta^{a}_{i}~{},$
(7)

with the relation

$\displaystyle\\{c,p\\}=\frac{\kappa\gamma}{3}~{}.$ (8)

Finally, the scalar constraint in this setting reduces to the simple form

$\displaystyle C_{g}^{N}=-\frac{3}{\kappa\gamma^{2}}N\sqrt{p}c^{2}~{}.$ (9)

The lapse function $N$ represents a gauge freedom.

In addition to the gravitational sector, we introduce a scalar field $\phi$
with an arbitrary potential $V(\phi)$ to investigate an early inflationary
period. Using the canonical variables for the scalar field, namely $\phi$ and
$\pi_{\phi}=p^{3/2}\dot{\phi}$, such that $\\{\phi,\pi_{\phi}\\}=1$, the total
Hamiltonian describing the coupled system can be written as

$\displaystyle
C^{N}=C^{N}_{g}+C^{N}_{m}=N\left(-\frac{3}{\kappa\gamma^{2}}\sqrt{p}c^{2}+p^{3/2}\rho\right)~{},$
(10)

where $\rho=\pi_{\phi}^{2}/(2p^{3})+V(\phi)$. The first Friedmann equation can
easily be recovered using the evolution equation – that is
$\dot{p}=\\{p,C^{N}\\}$ – with the choice $N=1$ and the Hamiltonian
constraint. In its usual form, it is written as

$\displaystyle
H^{2}\equiv\left(\frac{\dot{a}}{a}\right)^{2}=\frac{\kappa}{3}\rho~{}.$ (11)

Up to now, we have simply recovered the usual GR result within a specific
framework. An effective Hamiltonian including corrections from LQG is yet to
be constructed. The first step toward canonical quantization with well-defined
operators in the quantum theory is to rewrite the Hamiltonian constraint using
the holonomy of the connection. In other words, instead of deriving the
curvature operator $F_{ab}^{i}$ directly from the connection $c$, we use the
holonomy $h_{\Box_{ij}}$ of the connection on a fiducial square curve
$\Box_{ij}$ of length $\mu V_{0}^{1/3}$, with edges in the directions $i,j$.
The holonomy measuring the extent to which the parallel transport of a vector
around closed loops fails to preserve the transported vector, its form depends
on the chosen $SU(2)$ representation for the parallel transport along the
curves. This is known as the spin ambiguity. In standard LQC, the holonomy is
calculated using the fundamental 2D representation of $SU(2)$, but the
holonomy correction has been calculated for arbitrary representations in Ben
Achour et al. (2017). In this section, we describe the procedure for the
fundamental spin $1/2$ representation and we recall the procedure to follow
for general corrections.

In the harmonic gauge, where $N=p^{3/2}$, the Hamiltonian constraint in terms
of the curvature operator reads

$\displaystyle
C^{h}_{g}=-\frac{1}{2\kappa\gamma^{2}}p^{2}V_{0}^{2/3}\bar{e}^{a}_{i}\bar{e}^{b}_{j}F^{k}_{ab}~{},$
(12)

where the superscript $h$ stands for the harmonic gauge and $\bar{e}^{a}_{i}$
are the cotriads such that
$q_{ab}=a^{2}(t)\bar{e}^{i}_{a}\bar{e}^{j}_{b}\delta_{ij}$. The holonomy
corrected curvature operator in the fundamental representation is

$\displaystyle F^{k}_{ab}=\lim\limits_{\mu\rightarrow
0}\frac{-2}{\mu^{2}V_{0}^{2/3}}\Tr{h_{\Box_{ij}}\tau^{k}}\bar{e}^{i}_{a}\bar{e}^{j}_{b}~{},$
(13)

where $\tau^{k}$ are the generators of the $su(2)$ algebra and, in this case,
are represented by the Pauli matrices. By taking the limit $\mu\rightarrow 0$,
the usual definition of the curvature operator and of the Hamiltonian of
standard GR are recovered. However, LQG teaches us that the lowest nonzero
eigenvalue of the quantum area operator is $\lambda^{2}=4\sqrt{3}\pi\gamma
l_{pl}^{2}$, hence taking the limit down to zero is in principle not allowed.
One therefore chooses $\mu\rightarrow\bar{\mu}=\lambda/\sqrt{p}$ as an
estimator of the smallest possible length for the edge of the square curve.
The holonomy can then be calculated along a closed curve defined by

$\displaystyle h_{\Box_{ij}}$ $\displaystyle=h_{i}\circ h_{j}\circ h_{-i}\circ
h_{-j}~{},$ (14) $\displaystyle\textrm{where}\quad h_{\epsilon i}:$
$\displaystyle=\exp{\epsilon\mu c\tau_{i}}~{},$ (15)

with $\epsilon=\pm 1$. Putting everything together, one obtains the
Hamiltonian constraint of LQC coupled to a scalar field, that is

$\displaystyle
C^{h}=-\frac{3}{\kappa\gamma^{2}\bar{\mu}^{2}}p^{2}\sin^{2}(\bar{\mu}c)+p^{3}\rho~{}.$
(16)

This modified Hamiltonian can be recovered from the Hamiltonian of GR using
the substitution

$\displaystyle c^{2}\rightarrow\frac{\sin^{2}(\bar{\mu}c)}{\bar{\mu}^{2}}~{},$
(17)

usually called “the holonomy correction”. Finally, it is possible to derive
the modified Friedmann equation of LQC using $\dot{p}=\\{p,C^{h}\\}$, together
with the Hamiltonian constraint $C^{h}=0$. This leads to:

$\displaystyle
H^{2}=\frac{\kappa}{3}\rho\left(1-\frac{\rho}{\rho_{c}}\right)~{},$ (18)

with $\rho_{c}=3/(\kappa\gamma^{2}\lambda^{2})$. The remarkable feature of
this new equation is the resolution of the Big Bang singularity. When
$\rho\rightarrow\rho_{c}$ the Hubble parameter vanishes, as obvious from Eq.
(18), and a bounce occurs instead of the GR singularity111We however want to
emphasize that contrary to what is often believed, a past singularity is not
unavoidable in GR, even without exotic matter contents Barrau (2020); Renevey
et al. (2021).. When choosing the usual value for the Barbero-Immirzi
parameter $\gamma=0.2375$, the critical density is of the order of
$\rho_{c}\approx 0.41\rho_{\text{Pl}}$.

The general case for the curvature operator, calculated using an arbitrary
$d$-dimensional irreducible representation of $SU(2)$, was studied in
Vandersloot (2005); Chiou and Li (2009). A new closed formula for the
Hamiltonian of flat FLRW models regularized with arbitrary spins was found in
Ben Achour et al. (2017) and happens to be polynomial in the basic variables,
which corresponds to well-defined operators in the quantum theory (taking also
into account the inverse-volume corrections). The key-point lies in the fact
that the computation in a representation of spin $j$ of the trace of an
$SU(2)$ group element does not require the explicit knowledge of all its
matrix elements and can be reduced to an expression involving only the trace
in the fundamental representation and the class angle. The curvature operator
can then be written as

$\displaystyle
F_{ab}^{k}=\frac{-3}{d(d^{2}-1)}\frac{1}{\bar{\mu}^{2}V_{0}^{2/3}}\frac{\sin^{2}(\bar{\mu}c)}{\sin\theta}\partialderivative{\theta}\left(\frac{\sin(d\cdot\theta)}{\sin\theta}\right)\varepsilon^{k}_{ij}\bar{e}^{i}_{a}\bar{e}^{j}_{b}~{},$
(19)

with

$\displaystyle\theta=\arccos\left(\cos(\bar{\mu}c)+\frac{1}{2}\sin^{2}(\bar{\mu}c)\right)~{}.$
(20)

This was derived with a technique quite similar to the one described
previously for the holonomy regularization. Another technique to find the
curvature operator, called connection regularization, can also be effectively
considered Ben Achour et al. (2017). In this approach, a new definition of the
curvature operator, only valid in homogeneous space, is used and the result
for $F_{ab}^{k}$ is slightly different. In the literature, higher order
holonomy corrections were also investigated in details Mielczarek and
Szydlowski (2008); Hrycyna et al. (2009); Chiou and Li (2009). They arise when
higher order terms in powers of $\bar{\mu}$, usually neglected, are taken into
account in the expression for the regularized curvature. There could exist a
link between these higher order holonomy corrections and the contribution of
higher spin representations. However, it was shown in Ben Achour et al. (2017)
that these effects have actually very different physical consequences. At any
order in holonomy corrections, a physical Hilbert space can be rigorously
constructed and a complete family of Dirac observable can be identified Chiou
and Li (2009).

It is therefore mandatory to understand the cosmological implications of more
general holonomy corrections.

##  III Cosmology with arbitrary holonomy corrections

Instead of focusing on specific cases within the LQC framework with either
different spin representations or higher order terms, we remain as general as
possible for the expression for the holonomy correction. This can be studied
by the substitution

$\displaystyle c^{2}\rightarrow g^{2}(c,p)~{},$ (21)

where $g(c,p)$ is an arbitrary function such that, in the low energy limit,
standard GR is recovered, that is $g(c,p)\rightarrow c$. It should be
emphasized that this is not only a way of taking into account the lessons from
specific situations in LQC, but that this also make sense from a fully generic
quantization ambiguity/renormalization viewpoint. Furthermore, in the totally
constrained Hamiltonian formalism of GR, all the constraints are first class.
We do not relax this requirement so that the evolution operator keeps the
subspace of physical states invariant. As shown in Han and Liu (2018), keeping
only first class constraints and following the usual Dirac prescriptions adds
an extra condition on the function $g(c,p)$:

$\displaystyle g(c,p)=\frac{1}{\bar{\mu}}f(b)~{},$ (22)

where $f(b)$ is an arbitrary function of $b=\bar{\mu}c$, which must behave as
$f(b)\approx b$ at low energies (i.e. when $b\ll 1$). As opposed to Han and
Liu (2018), we used the parameter $\lambda$ to respect the units of length of
$g(c,p)$ and be consistent with LQC. Fortunately, this is fully compatible
with the correction given in Eq. (19), ensuring that any holonomy modification
coming from an arbitrary spin-representation will keep the algebra of
constraints consistent. It is quite remarkable that the “anomaly freedom”
requirement (see Wu et al. (2018); Cailleteau et al. (2012); Barrau et al.
(2015) for general considerations) allows one to sharpen the general
expression, in a way precisely compatible with known corrections expected in
the loop framework. The simple – and mandatory – fact that the evolution
vector is asked to be parallel to the submanifold of constraints severely
reduces the a priori freedom.

The modified equations of motion for the canonical variables $c$ and $p$ are
calculated using Hamilton’s equations. It is more natural to write them in
terms of $p$ and $b$. Together with the Hamiltonian constraint, they take the
form

$\displaystyle\dot{b}$
$\displaystyle=-\frac{\lambda\kappa\gamma}{2}\rho(1+w)~{},$ (23)
$\displaystyle\dot{p}$
$\displaystyle=\frac{2}{\gamma\lambda}pf(b)f^{\prime}(b)~{},$ (24)
$\displaystyle\frac{\rho}{\rho_{c}}$ $\displaystyle=f^{2}(b)~{},$ (25)

where $f^{\prime}(b)$ should be understood as $\differential
f(b)/\differential b$, $w=P/\rho$ and $P=\pi_{\phi}^{2}/(2p^{3})-V(\phi)$. To
derive Eq. (23) we used the continuity equation

$\displaystyle\dot{\rho}=-3\frac{\dot{p}}{2p}\rho(1+w)~{}.$ (26)

A general modified Friedmann equation can be found using $H=\dot{p}/(2p)$
together with the constraint (25) and can be written as

$\displaystyle H^{2}=\frac{\kappa}{3}\rho\left(f^{\prime}(b)\right)^{2}~{}.$
(27)

Since, by construction, $f(b)\rightarrow b$ when $b\rightarrow 0$, we indeed
recover the usual Friedmann equation $H^{2}=\kappa\rho/3$ in this limit. Let
us now investigate the behaviour of the model starting in a regime where GR is
valid and going backward in cosmic time $t$, toward a classical singularity
associated with $\rho\rightarrow\infty$ in GR. Using the null energy
condition, $w\geq-1$ together with Eqs. (23) and (24) one can easily show that
$\differential b/\differential t<0$ $\forall t$. Hence, if we start with $b>0$
in the GR regime, $b$ is always positive and increasing when going backward in
time. This is expected as, in the GR limit, the proportionality relation
$\rho\propto b$ is satisfied and the density increases in the past direction.
Furthermore, since one has $f(b)\approx b$ in the GR regime, the function $f$
is monotonic and strictly increasing with $b$ around $b=0$. In the case where
there exists a local maximum $b_{bounce}>0$, implying
$f^{\prime}(b_{bounce})=0$, one can see with Eqs. (25) and (27) that the
density has to reach a critical value $\rho_{b}$, where the Hubble parameter
vanishes. It is therefore meaningful to conclude that the Big Bang singularity
is resolved by a bouncing scenario of geometrical origin if and only if the
holonomy correction reaches a local maximum. If the function $f$ is strictly
monotonic, two different scenarios have to be considered. Either
$\lim\limits_{b\rightarrow\infty}f(b)=\infty$, in which case the singularity
is not resolved, or $\lim\limits_{b\rightarrow\infty}f(b)<\infty$ and the
situation is similar to eternal inflation where $\rho$ tends to a final
constant value behaving as a positive cosmological constant.

##  IV First remarks on inflation from the holonomy correction

A natural question arising in this framework is to wonder whether it is
possible to describe a long-lasting phase of inflation using only a
modification of the holonomy correction without the need for a fluid
satisfying the equation of state $w<-1/3$. In the usual LQC framework, the
quantum geometrical super-inflation occurring after the bounce cannot account
for more than a few e-folds and most of the known inflationary features are
due to a hypothetical massive scalar field filling the Universe. Furthermore,
one should also investigate if the inflation associated with generalized
holonomies could explain the quasi scale invariance of the power spectrum
observed in the cosmological microwave background (CMB). The answer turns out
to be positive. However, important drawbacks inherent to the construction will
be mentioned in this section. We assume here that the content of the Universe
is a massless scalar field, that is $w=1$ at all scales. We also restrict
ourselves to holonomy corrections such that there exists a bounce so as to
keep the huge benefit of the singularity resolution.

First of all, the correction $f(b)$ has to be chosen so as to ensure an
exponential growth of the scale factor $a$. This is achieved by $\dot{H}\simeq
0$: the Hubble parameter is nearly constant during inflation. From the
modified Friedmann equation (27) and the Hamiltonian constraint
$\eqref{eq:hamiltonian_constraint}$, one can easily conclude that the holonomy
correction should take the following form:

$\displaystyle f(b)\sim\sqrt{b}~{}.$ (28)

We assume that $f(\pi)=0$ and call $b_{bi}$ the value of $b$ at the beginning
of inflation. We set $b=b_{rh}$ at the reheating, when classical cosmology is
recovered. It is therefore necessary to find a function $f(b)$ such that : i)
$f(b)\sim b$ around $b=0$, to recover GR, ii) $f(b)\sim\sqrt{b}$ for
$b_{rh}<b<b_{bi}$ during the inflationary period, and iii) $f(b)$ reaches a
maximum between $b_{bi}$ and $b=\pi$ to induce a bounce. In order to make
explicit that such a function can be constructed, we give an arbitrary
example:

$\displaystyle f^{2}(b)$
$\displaystyle=(1+C_{1}b)^{1-\alpha}\sin^{2}(b)\frac{1}{C_{1}+1}\sum\limits_{n=0}^{C_{1}}\cos^{2n}(b)$
(29)
$\displaystyle=\frac{(1+C_{1}b)^{1-\alpha}}{1+C_{1}}\left(1-\cos^{2(C_{1}+1)}(b)\right)~{},$
(30)

where $C_{1}$ can be chosen in accordance with the parameters $b_{bi}$ and
$b_{rh}$ so as to obtain the desired inflationary duration. The parameter
$\alpha$ can be chosen to make the power spectrum slightly red. To illustrate
the behaviour of this holonomy correction, Eq. (30) is plotted on Fig. 1 for
different values of $C_{1}$, assuming $\alpha=0$. The higher the value of
$C_{1}$, the lower the value of $b$ at the end of inflation (denoted $b_{rh}$)
and the closer $b_{bi}$ to $\pi$. Let us take a closer look at Eq. (30) in the
case $C_{1}\gg 1$. If $b\gg 1/\sqrt{C_{1}}$, one indeed recovers $f^{2}(b)\sim
b$, because the cosine term is small compared to $1$ and $C_{1}b\gg 1$. On the
other hand, when $b\ll 1/C_{1}$ one has $(1+C_{1}b)\rightarrow 1$ while
$(1-\cos^{2C_{1}}(b))\sim C_{1}b^{2}$, leading to $f^{2}(b)\sim b^{2}$.
However, in the intermediate case, $1/C_{1}\ll b\ll 1/\sqrt{C_{1}}$,
$(1+C_{1}b)\sim C_{1}b$ while $(1-\cos^{2C_{1}}(b))$ still behaves as
$C_{1}b^{2}$, leading to $f^{2}(b)\sim b^{3}$. Hence, the transition at
reheating does not straightforwardly go from $f^{2}(b)\sim b$ to $f^{2}(b)\sim
b^{2}$, but goes through transition phase, behaving as $f^{2}(b)\sim b^{3}$.
This is not problematic for our purpose, but it is worth being pointed out.
Equation (25) shows that in order to have $N$ e-folds of inflation in a matter
dominated universe, one needs $f^{2}(b_{rh})e^{3N}=f^{2}(b_{bi})\sim 1$.
Therefore, if one chooses $C_{1}\sim e^{6N}$, the inflationary period ends
when $b_{rh}\sim e^{-3N}\implies f^{2}(b_{rh})\sim e^{-3N}$, leading to the
desired $N$ e-folds. Then, the transition phase takes place while
$e^{-6N}<b<e^{-3N}$ – that can call reheating – and finally we recover
classical cosmology for $b<e^{-6N}$. Obviously, a quite strong fine-tuning is
needed but solutions matching all the requirements can be found. We shall
discuss this issue later.

FIG. 1: Holonomy corrections as described by Eq. (30) for different values of
the parameter $C_{1}$.

Another potential problem to consider is related with the evolution of the
energy density. In usual models of inflation, where the exponentially
accelerating expansion is produced by a negative pressure fluid, the density
of the latter stays roughly constant around (slightly above) the density of
reheating $\rho\sim\rho_{rh}$. However, the density of radiation, dust or a
massless scalar field evolve as $a^{-4}$, $a^{-3}$ and $a^{-6}$, respectively,
thus increase exponentially during inflation as well (when thinking backward
in time). If inflation is of quantum geometrical origin and not caused by the
standard mechanism of a field slowly rolling on its potential, there is no
reason to assume that the usual “contents” are not present and possibly
dominant. This means that if one requires a reheating around the GUT scale
$T_{GUT}\sim 10^{15}$GeV and imposes roughly $N=65$ e-folds of exponential
expansion, the energy at the beginning of inflation $T_{bi}$ would be vastly
trans-Planckian. On the other hand, one can impose the maximum energy at the
Planck scale $T_{bi}\sim T_{pl}$ and simply ask for more e-folds before the
reheating than after, to solve the usual cosmological problems. However this
translates into $T_{rh}\leq 1$ TeV. Although unusual, this value is not
strictly ruled out by observations.

Finally, it is worth emphasizing that the scalar power spectrum, as observed
in the CMB, cannot be easily recovered for other contents than a massive
scalar field. To see this, let us consider a massless scalar field and a
holonomy correction of the form (30) with a sufficiently long inflationary
phase. In this example, we choose $T_{bi}\sim T_{pl}$ and roughly $N=35$
e-folds. The gauge invariant scalar perturbation can, in this context, be
described by the Mukhanov-Sasaki (MS) variables $v_{k}(\eta)$ and
$z(\eta)=\dot{\phi}a/H$, where $\eta$ is the conformal time, satisfying the
equation

$\displaystyle
v_{k}^{\prime\prime}(\eta)+\left(k^{2}-\frac{z^{\prime\prime}(\eta)}{z(\eta)}\right)v_{k}(\eta)=0~{}.$
(31)

It should be noticed that this classical equation is also the one fulfilled by
perturbations in the dressed metric/hybrid quantization approach to LQC Agullo
et al. (2013); Fernandez-Mendez et al. (2012). The scalar power spectrum is
obtained by

$\displaystyle\mathcal{P}_{S}(k)=\left.\frac{k^{3}}{2\pi^{2}}\absolutevalue{\frac{v_{k}}{z}}^{2}\right|_{k=aH}\equiv\left.\frac{k^{3}}{2\pi^{2}}\absolutevalue{\mathcal{R}_{k}}^{2}\right|_{k=aH},$
(32)

where, in the standard inflation theory, the curvature perturbation
$\mathcal{R}_{k}$ is constant for super-Hubble modes $k\leq aH$. The scalar
power spectrum therefore keeps a low amplitude matching the observed value. In
the case of the massless scalar field with an inflationary stage due to the
holonomy correction, one can show using Eq. (31) that $\mathcal{R}_{k}$ is not
constant for super-Hubble modes anymore. During the inflationary period, the
Hubble parameter $H$ is constant, $a\propto e^{Ht}$ and
$\dot{\phi}=2\sqrt{\rho}$, meaning that the second MS variable behaves as

$\displaystyle z(a)\sim a^{-1/2(1+3w)}.$ (33)

In the considered example $w=1$ but we keep the equation of state arbitrary so
that the conclusion remains general. Rewriting the Mukhanov-Sasaki equation
(31) with respect to $a$ and taking the limit of large values of $a$, as one
might expect for super-Hubble modes, we get

$\displaystyle\derivative[2]{v_{k}}{a}+\frac{1}{a}\derivative{v_{k}}{a}-\frac{2}{a^{2}}v_{k}=0,$
(34)

with solutions $v(a)\sim a$ and $v(a)\sim 1/a^{2}$. In standard inflation with
a massive scalar field, where $w\approx-1$, choosing the Bunch-Davies vacuum
comes down to selecting the behaviour $v(a)\sim a$ in quasi-de Sitter space,
hence one recovers $\mathcal{R}_{k}=v_{k}/z\sim a^{0}$. In the case of a
massless field ($w=1$), such a choice of vacuum leads to the behaviour
$\mathcal{R}_{k}\sim a^{3}$ and therefore the scalar perturbations with super-
Hubble modes exponentially increase during inflation. To illustrate this, the
scalar power spectrum computed from Eq. (31) for an inflation induced by a
holonomy correction with massless scalar field is shown in Fig. 2. As one can
see, the spectrum is scale invariant in the UV and keeps its usual shape but
the amplitude is meaningless. In principle, there might exist another vacuum
selecting $v(a)\sim 1/a^{2}$ and one would recover $\mathcal{R}_{k}\sim
a^{0}$. However, in the case of other matter contents, such as radiation
($w=1/3$) or dust ($w=0$), the scalar perturbations $\mathcal{R}_{k}$ cannot
be frozen for $k\geq aH$ in geometric inflation. This is an important point,
often forgotten, which should be taken into account in phenomenological
studies.

FIG. 2: Scalar power spectrum obtained with a massless scalar field from a
stage of “geometrical” inflation entirely produced by the modified holonomy
correction.

It should however be emphasized that perturbations were, more than a decade
ago, reunderstood in the framework of the effective theory of inflation Cheung
et al. (2008). This “new” approach somehow disentangles the question of
perturbations themselves from the question of the process which generates
them. From this point of view, this work is anyway interesting as the
perturbations can be treated even without any use of a scalar field. All that
is needed is a background and an effective “clock”.

##  V Focus on the duration of inflation

Detailed studies of the duration of inflation as predicted by LQC have already
been made in Linsefors and Barrau (2015); Martineau et al. (2017); Bolliet et
al. (2017) for initial conditions set in the remote past. In the case of a
massive scalar field playing the role of the inflaton, it was shown that the
number of e-folds, as predicted by LQC, is around $N=140$. Interestingly, the
number of e-folds varies only slightly with respect to most contingent
parameters. This can easily be seen by studying trajectories in the potential
energy $x(t)=\sqrt{V(\phi(t))/\rho_{c}}$ versus kinetic energy
$y(t)=\sqrt{\dot{\phi}/(2\rho_{c})}$ plane. The number $N=140$ corresponds to
a quadratic potential $V(\phi)=m^{2}/2\phi^{2}$ with a mass parameter $m\sim
1.2\cdot 10^{-6}$. For other potentials, the number of e-folds can be
different, but remains of this order of magnitude as long as one deals with
confining potentials. Although anisotropies can slightly decrease this value,
the duration of inflation is a robust property of the background in LQC (a
different proposal was however suggested in Ashtekar and Sloan (2011) but
relies on conditions set at the bounce, a path that we do not follow here).
The fact that this number is way smaller than naïve expectations is a key
feature of bouncing scenarios. One can therefore naturally wonder if this
prediction changes for different shapes of the holonomy correction. The set of
initial prebounce conditions $(x_{0},y_{0})$ is described using a phase
parameter $\delta\in[0,2\pi[$ such that

$\displaystyle
x_{0}=\sqrt{\frac{\rho_{0}}{\rho_{c}}}\cos{\delta}\quad\textrm{and}\quad
y_{0}=\sqrt{\frac{\rho_{0}}{\rho_{c}}}\sin{\delta},$ (35)

where $\rho_{0}$ is the initial energy density and defines how far away in the
past initial conditions are set. It has been checked that the chosen
distribution for $\delta$ is conserved over time. For completeness, we
reproduce the probability distribution function (PDF) using the settings of
Linsefors and Barrau (2015) in Fig.3 with a uniform distribution of the phase
parameter $\delta$. In the following, the mean value of the number of e-folds
for different holonomy corrections is calculated using such a PDF.

FIG. 3: Probability distribution function of the number of e-folds of
inflation for a quadratic potential with $m=1.2\cdot 10^{-6}$ in LQC.

In this section, we investigate the effect of the form of the holonomy
correction on the duration of inflation. For this purpose, we follow the
strategy described in Martineau et al. (2017); Bolliet et al. (2017): we
assume a uniformly distributed set of initial conditions for $\delta$ and the
initial time $t=0$ sufficiently far in the contracting branch so as to remain
away from the quantum gravitational regime. In order to stay conservative
about the holonomy correction shape, we study different general properties
such as an asymmetry, an increase of the maximum of $f^{2}$, and a flattening
of the correction on the length of inflation. Overall, we want to keep the
asymptotic behaviour of LQC, i.e. $f(0+\delta b)=f(\pi+\delta b)=\delta
b+\mathcal{O}(\delta b^{2})$, hence we choose a correction of the form

$\displaystyle
f^{2}(b)=\sin^{2}(b)\left(1+A_{1}b^{n_{1}}(b-\pi)^{n_{2}}\right),$ (36)

where $n_{i}\geq 1$, $i=1,2$ and $A_{1}\geq 0$, or

$\displaystyle
f^{2}(b)=\sin^{2}(b)\frac{1}{C_{1}+1}\sum\limits_{n=0}^{C_{1}}\cos^{2n}(b),$
(37)

where $C_{1}\geq 1$. Such parametrizations ensure the desired behaviour at the
fix points $b=0,\pi$.

FIG. 4: Left: asymmetric holonomy corrections. Right: mean value of the PDF of
the number of inflationary e-folds.

The results of the numerical simulations for different left and right
asymmetries are shown in Figs. 4 and 5. The left panels represent the chosen
holonomy corrections whereas the right panels display the mean value of the
PDF of the number of e-folds when scanning the full range of initial phases.
It should be pointed out that the PDF is narrow enough ($\sigma\sim 10$
e-folds) so that its first moment gives the relevant information. The strength
of the asymmetry is measured by the ratio $I_{2}/I_{1}$ between the integral
of $f^{2}(b)$ on $[0,b_{max}]$ and the integral on $[b_{max},\pi]$, where
$b_{max}$ is the value that maximizes $f^{2}(b)$. The general trend is a
decrease of the number of e-folds (although this is not true for very small
deformations).

FIG. 5: Left: asymmetric holonomy corrections. Right: mean value of the PDF of
the number of inflationary e-folds.

The results of the simulations for different amplitudes and plateaulike
functions are displayed in Figs. 6 and 7, respectively. In general, increasing
the energy density at the bounce decreases the number of e-folds and
flattening the correction (or reducing the energy of the bounce) increases the
length of inflation. Since the bounce energy density grows up with
$\max{f^{2}}$, as shown by Eq.(25), this result seems counter-intuitive at
first sight. It is true that the inflation duration generally increases with
the bounce energy scale, but this is not the only important ingredient
relevant for the number of inflationary e-folds. Another key parameter is the
ratio between potential and kinetic energies at the bounce
$x(t_{b})/y(t_{b})$. For kinetic energy dominated bounces, more time will
indeed be required before the potential energy finally dominates and the slow-
roll conditions are met. This leads to a lower inflation energy scale than for
potential energy dominated bounces. Thus, bounces dominated by kinetic energy
will actually lead to shorter inflation phases. This is the reason why the PDF
of the number of inflationary e-folds in LQC is peaked around low values of
$N$, close to the experimental lower bound. Setting initial conditions in the
remote past of the contracting branch with a flat PDF on the initial phase of
the scalar field leads to trajectories with kinetic energy dominated bounces,
with ratios $x(t_{b})/y(t_{b})$ typically of order $10^{-6}$. In the present
study, the higher $\max{f^{2}}$ , the more kinetic energy dominated the bounce
is. This phenomenon is shown in Fig.8. First, this indeed confirms in the top
right and middle left plots, representing respectively $\max{f^{2}}$ wrt.
$x(t_{b})/y(t_{b})$ and the number of e-folds wrt. $x(t_{b})/y(t_{b})$, that
$\max{f^{2}}$ has a significant effect on this ratio, which in turn affects
the number of e-folds as expected. Using the middle and bottom right plots,
one can see that it takes longer for a kinetic energy dominated bounce to
start the inflationary period, leading to a lower energy density at inflation,
thus lowering the number of e-folds.

FIG. 6: Left: amplitude-varying holonomy corrections. Right: mean value of the
PDF of the number of inflationary e-folds.

There are several lessons to be learnt from this analysis. There is an obvious
loss of predictivity associated with a possible generalization of the holonomy
correction. This can hardly come as a surprise. Whether this should be
considered as good or bad news depends on the point of view. This obviously
makes the theory less easily falsifiable. But this also opens up interesting
possibilities. It shows that desired cosmological behaviours can be obtained
by purely geometrical effects, relaxing the need for exotic matter contents.
This also means that, in principle, a good knowledge of the cosmological
dynamics can severely constrain the shape of the holonomy correction. In
particular, the minimum number of e-folds $N>60$ required to solve the horizon
problem sets a constraint on the maximum of $f^{2}(b)$ as $N$ decreases when
the maximum density at the bounce increases. Most arbitrary shapes are
problematic. It should however be underlined that as long as one deals with
moderate corrections to the usual $\sin^{2}(\bar{\mu}c)/\bar{\mu}^{2}$, the
effects on the background dynamics remain quite small and the core of the
known LQC predictions remains valid.

FIG. 7: Left: plateau-like holonomy corrections. Right: mean value of the PDF
of the number of inflationary e-folds.

The ambiguities considered in this work do not, by themselves, reveal a
theoretical failure. So as to regularize quantum operators associated with
nonlinear functionals of the fundamental fields, one relies on the
diffeomorphism invariant prescription of “point-splitting” Perez (2006). It
happens to be that the regulator can be removed without encountering UV
divergences. One is then left with a well defined quantum Hamiltonian
constraint, at the price of having many different quantum theories. This is
obviously reminiscent of the usual problem of renormalization of quantum
fields: the correct theory must be fixed by the renormalization conditions.

FIG. 8: Top left: amplitude-varying holonomy corrections. Top right: maximum
of $f^{2}(b)$ wrt. to the ratio of potential to kinetic energy of the scalar
field at the bounce, $x(t_{b})/y(t_{b})$. Middle left: mean value of the PDF
of the number of inflationary e-folds wrt. $x(t_{b})/y(t_{b})$. Middle right:
mean value of the elapsed time between the bounce and inflation wrt.
$x(t_{b})/y(t_{b})$. Bottom left: mean value of the scalar field at the
beginning of inflation wrt. $x(t_{b})/y(t_{b})$. Bottom right: mean value of
the ratio of the density at the beginning of inflation over the one at the
bounce wrt. $x(t_{b})/y(t_{b})$.

##  VI Effects on the primordial power spectra

We now analyze the effects of a modified holonomy correction on the scalar
power spectrum when considering an inflationary period generated by a massive
scalar field of mass $m=1.2\times 10^{-6}~{}m_{\text{Pl}}$. We keep the form
of the Mukhanov-Sasaki equation for the scalar perturbations unchanged. This
allows to specifically investigate the effects of the modification of the
background dynamics, due to the new holonomy correction, on the shape of the
spectra, independently of possible modifications of the perturbations
propagation equation (which is anyway beyond the scope of this study).

Figures 9 and 10 display comparisons of the usual LQC spectrum (in grey) with
spectra obtained using new background dynamics due to modified holonomy
corrections (in blue). Figures 9 represents a double bounce while a single
bounce at higher energy is shown in Fig. 10. The new features in the spectra
corresponding to the modified dynamics are hardly distinguishable from the
usual spectrum. The spectra remain (almost) scale invariant and consistent
with CMB data in the so-called ultraviolet region, that is the region
corresponding to the highest presented wave numbers (and even higher values of
$k_{c}$), in which the observable scales are located for almost all the
parameter space. The modifications to $f^{2}(b)$ only impact the details of
the oscillations in the intermediate regime. This was expected as those
oscillations are mainly associated with the bounce and their details depend on
the detailed dynamics at this time. The ultraviolet regime, i.e the one of
interest for comparison with data, is mostly independent of the bounce
dynamics as the wavenumbers of the corresponding modes are much greater than
the potential $z^{\prime\prime}/z$ during the bounce phase. In other words,
those modes do not even feel the presence of the bounce(s) It is however worth
emphasizing that, as highlighted in Martineau and Barrau (2018), the position
of the oscillations is dictated by the value of the potential
$z^{\prime\prime}/z$ evaluated at the bounce. Thus, sharper bounces would lead
to oscillations in the spectra located further away in the ultraviolet regime.
It is also worth to underline that modifications of $f^{2}(b)$ around
$b\approx 0$ or $\pi$ would modify the low energy behavior and have a deeper
impact on the spectra.

FIG. 9: Left: holonomy corrections for standard LQC (in gray) and for a double
bounce (in blue). Right: comparison between the associated scalar power
spectra (standard LQC in gray and with modified holonomy correction in blue).
FIG. 10: Left: holonomy corrections for standard LQC (in gray) and for a
higher energy bounce (in blue). Right: comparison between the associated
scalar power spectra (standard LQC in gray and with modified holonomy
correction in blue).

We have previously seen that an asymmetry (in $b$) of the holonomy correction
can have a significant influence on the evolution of the inflaton in the high
energy regime. It is therefore important to estimate the effects of
asymmetries in $f^{2}(b)$ on the spectra. Associated results are shown in
Figs. 11 and 12. The main noticeable effect is an extension of the oscillatory
regime toward higher wave numbers. This extension is however (for all the
functions tested) not sufficient to be game changing for the noticeable
features of primordial power spectra in loop inspired cosmologies. All usual
conclusions still hold. This effect only makes a possible
detection/invalidation of the bounce in CMB data slightly more probable as it
decreases a little the level of fine tuning required to bring the window of
comoving wave numbers associated with CMB observations toward the oscillatory
regime.

FIG. 11: Left: holonomy corrections for standard LQC (in grey) and for a left
asymetric bounce (in blue). Right: comparison between the associated scalar
power spectra (standard LQC in grey and with modified holonomy correction in
blue). FIG. 12: Left: holonomy corrections for standard LQC (in grey) and for
a right asymetric bounce (in blue). Right: comparison between the associated
scalar power spectra (standard LQC in grey and with modified holonomy
correction in blue).

It should be underlined that the situation might be totally different in the
case of the Deformed Algebra (DA) approach to perturbations, where the
equation of propagation of perturbations is usually given byBarrau et al.
(2018); Schander et al. (2016):

$v_{k}^{\prime\prime}(\eta)+\left(\Omega(\eta)k^{2}-\frac{z^{\prime\prime}(\eta)}{z(\eta)}\right)v_{k}(\eta)=0~{},$
(38)

with $\Omega=1-2\rho/\rho_{c}$. The change of sign of $\Omega$ at high
energies can be interpreted as a switch from a Lorentzian to an Euclidean
geometry Bojowald and Mielczarek (2015). The apparition of this $\Omega$
function in the MS equation comes from its presence in the anomaly-free
algebra of constraints, more precisely in the Poisson bracket between scalar
constraintsBarrau et al. (2015):

$\left\\{S^{\text{a.f}}[M],S^{\text{a.f}}[N]\right\\}\approx\Omega
D^{\text{a.f}}\left[q^{ab}\left(M\partial_{b}N-N\partial_{b}M\right)\right],$
(39)

where $S$ and $D$ correspond respectively to the scalar and diffeomorphism
constraints, the label a.f stands for ”anomaly-free”, and the $\approx$ symbol
means that this relation is satisfied on the hypersurface of constraints. But
the $\Omega$ function appears in this algebra as a consequence of the
modification of the kinetic term in the scalar constraint (by the holonomy
correction, that is Eq.(17) in the usual case). If we assume a 1+1 dimensional
toy model then $\Omega=(1/2)d^{2}f^{2}(b)/db^{2}$, which, in the case
$f^{2}(b)=\sin^{2}(b)$, and using the Hamiltonian constraint, gives back the
previous expression $\Omega=\cos(2b)=1-2\rho/\rho_{c}$. It is thus clear that
modifying the form of $f(b)$ would unavoidably modify the expression of
$\Omega$ and therefore the Mukhanov-Sasaki equation. As a consequence, the
shape of the primordial power spectra could be modified in a much more
important way in the Deformed Algebra scheme than in the approach presented in
this manuscript. This point definitely requires deeper investigations that are
not the purpose of this paper.

We have focused on the scalar spectrum which is the more interesting one from
the viewpoint of observations (and the more intricate to calculate). The
general trends however obviously remain true for the tensor spectrum as the
new features are due to the modified background dynamics.

##  VII Conclusion

Quantization ambiguities are unavoidable. The only serious requirement in this
framework is to recover general relativity (or the Wheeler-DeWitt equation) in
the low energy limit of the theory. This leaves an infinite dimensional set of
ambiguities in the choice of the function $f$. Investigating few of its
phenomenological properties was the goal of this study.

We have shown that a long enough inflationary stage can be generated by an
appropriate generalized holonomy correction without the need for matter
violating the energy conditions. This however requires a serious amount of
fine-tuning in the choice of the parameters entering the definition of the
function. Nevertheless this concept is ill-defined when no natural measure is
available. More profoundly, it is worthwhile emphasizing that the fine-tuning
is a real issue only when it is required to produce a situation which is a
priori singular or specific (in a Bayesian sense). Otherwise, it is just a
mere version of the nonproblematic tautology “if the laws had been different,
the World would be different”. In addition, we have shown that even though a
nearly scale invariant primordial power spectrum can be generated, there are
some (often forgotten) associated issues to face, making the spectrum hardly
compatible with data.

The influence of the shape and amplitude of the holonomy correction on the
background dynamics has also been investigated in details. The effects are
quite weak and do not change drastically the usual conclusions of LQC as long
as the used function remains close enough to the standard case. Interestingly
some new effects tend to decrease the number of e-folds. This is a good news
for phenomenology. With 140 inflationary e-folds, all the subtle quantum
gravitational effects are deeply super-Hubble today. Smaller values bring back
the hope that some footprints of the bounce might be seen in the CMB,
especially when taking into account nonlinear effects where different modes
fail to fully decouple one from the other.

##  VIII Acknowledgements

We would like to thank Alejandro Perez who was at the origin of this work and
who provided many important inputs and explanations. We also thank Salvatore
Ribisi and Lautaro Amedei for valuable discussions.

## References

  * Ashtekar and Bianchi (2021) A. Ashtekar and E. Bianchi, Rept. Prog. Phys. 84, 042001 (2021), eprint 2104.04394.
  * Rovelli (2004) C. Rovelli, _Quantum gravity_ , Cambridge Monographs on Mathematical Physics (Univ. Pr., Cambridge, UK, 2004).
  * Rovelli and Vidotto (2014) C. Rovelli and F. Vidotto, _Covariant Loop Quantum Gravity_ , Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2014), ISBN 1107069629, 9781107069626, 9781316147290\.
  * Kiefer (2007) C. Kiefer, in _Time and Matter 2007_ (University of Nova Gorica Press, 2007).
  * Barrau et al. (2014) A. Barrau, T. Cailleteau, J. Grain, and J. Mielczarek, Class.Quant.Grav. 31, 053001 (2014), eprint 1309.6896.
  * Quevedo (2016) F. Quevedo (2016), eprint 1612.01569.
  * Ashtekar et al. (2008) A. Ashtekar, A. Corichi, and P. Singh, Phys. Rev. D 77, 024046 (2008), eprint 0710.3565.
  * Ashtekar and Singh (2011) A. Ashtekar and P. Singh, Class. Quant. Grav. 28, 213001 (2011), eprint 1108.0893.
  * Agullo et al. (2012) I. Agullo, A. Ashtekar, and W. Nelson, Phys. Rev. Lett. 109, 251301 (2012), eprint 1209.1609.
  * Diener et al. (2014) P. Diener, B. Gupt, and P. Singh, Class. Quant. Grav. 31, 105015 (2014), eprint 1402.6613.
  * Ashtekar and Barrau (2015) A. Ashtekar and A. Barrau, Class. Quant. Grav. 32, 234001 (2015), eprint 1504.07559.
  * Bolliet et al. (2015) B. Bolliet, J. Grain, C. Stahl, L. Linsefors, and A. Barrau, Phys.Rev. D91, 084035 (2015), eprint 1502.02431.
  * Alesci and Cianfrani (2015) E. Alesci and F. Cianfrani, Phys. Rev. D92, 084065 (2015), eprint 1506.07835.
  * Martineau et al. (2017) K. Martineau, A. Barrau, and S. Schander, Phys. Rev. D95, 083507 (2017), eprint 1701.02703.
  * Gielen and Oriti (2018) S. Gielen and D. Oriti, Phys. Rev. D 98, 106019 (2018), eprint 1709.01095.
  * Perez (2006) A. Perez, Phys. Rev. D 73, 044007 (2006), eprint gr-qc/0509118.
  * Amadei et al. (2022) L. Amadei, A. Perez, and S. Ribisi (2022), eprint 2203.07044.
  * Ben Achour et al. (2017) J. Ben Achour, S. Brahma, and M. Geiller, Phys. Rev. D 95, 086015 (2017), eprint 1612.07615.
  * Barrau (2020) A. Barrau, Eur. Phys. J. C 80, 579 (2020), eprint 2005.04693.
  * Renevey et al. (2021) C. Renevey, A. Barrau, K. Martineau, and S. Touati, JCAP 01, 018 (2021), eprint 2010.13542.
  * Vandersloot (2005) K. Vandersloot, Phys. Rev. D 71, 103506 (2005), eprint gr-qc/0502082.
  * Chiou and Li (2009) D.-W. Chiou and L.-F. Li, Phys. Rev. D 80, 043512 (2009), eprint 0907.0640.
  * Mielczarek and Szydlowski (2008) J. Mielczarek and M. Szydlowski, Phys. Rev. D 77, 124008 (2008), eprint 0801.1073.
  * Hrycyna et al. (2009) O. Hrycyna, J. Mielczarek, and M. Szydlowski, Gen. Rel. Grav. 41, 1025 (2009), eprint 0804.2778.
  * Han and Liu (2018) Y. Han and M. Liu, Class. Quant. Grav. 35, 105017 (2018), eprint 1711.04991.
  * Wu et al. (2018) J.-P. Wu, M. Bojowald, and Y. Ma, Phys. Rev. D 98, 106009 (2018), eprint 1809.04465.
  * Cailleteau et al. (2012) T. Cailleteau, A. Barrau, J. Grain, and F. Vidotto, Phys. Rev. D86, 087301 (2012), eprint 1206.6736.
  * Barrau et al. (2015) A. Barrau, M. Bojowald, G. Calcagni, J. Grain, and M. Kagan, JCAP 1505, 051 (2015), eprint 1404.1018.
  * Agullo et al. (2013) I. Agullo, A. Ashtekar, and W. Nelson, Class.Quant.Grav. 30, 085014 (2013), eprint 1302.0254.
  * Fernandez-Mendez et al. (2012) M. Fernandez-Mendez, G. A. Mena Marugan, and J. Olmedo, Phys. Rev. D 86, 024003 (2012), eprint 1205.1917.
  * Cheung et al. (2008) C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, JHEP 03, 014 (2008), eprint 0709.0293.
  * Linsefors and Barrau (2015) L. Linsefors and A. Barrau, Class. Quant. Grav. 32, 035010 (2015), eprint 1405.1753.
  * Bolliet et al. (2017) B. Bolliet, A. Barrau, K. Martineau, and F. Moulin, Class. Quant. Grav. 34, 145003 (2017), eprint 1701.02282.
  * Ashtekar and Sloan (2011) A. Ashtekar and D. Sloan, Phys. Lett. B694, 108 (2011), eprint 0912.4093.
  * Martineau and Barrau (2018) K. Martineau and A. Barrau, Universe 4, 149 (2018), eprint 1812.05522.
  * Barrau et al. (2018) A. Barrau, P. Jamet, K. Martineau, and F. Moulin, Phys. Rev. D 98, 086003 (2018), eprint 1807.06047.
  * Schander et al. (2016) S. Schander, A. Barrau, B. Bolliet, L. Linsefors, J. Mielczarek, and J. Grain, Phys. Rev. D 93, 023531 (2016), eprint 1508.06786.
  * Bojowald and Mielczarek (2015) M. Bojowald and J. Mielczarek, JCAP 08, 052 (2015), eprint 1503.09154.