# Structural Motifs and Bonding in Two Families of Boron Structures Predicted
at Megabar Pressures

Katerina P. Hilleke  Department of Chemistry, State University of New York at
Buffalo, Buffalo, NY 14260-3000, USA  Tadashi Ogitsu  Lawrence Livermore
National Laboratory, Livermore, CA, 94550, USA  Shuai Zhang  Laboratory for
Laser Energetics, University of Rochester, Rochester, NY, 14623, USA  Eva
Zurek<EMAIL_ADDRESS>Department of Chemistry, State University of New
York at Buffalo, Buffalo, NY 14260-3000, USA

###### Abstract

The complex crystal chemistry of elemental boron has led to numerous proposed
structures with distinctive motifs as well as contradictory findings. Herein,
evolutionary structure searches performed at 100 GPa have uncovered a series
of potential new metastable phases of boron, and bonding analyses were carried
out to elucidate their electronic structure. These polymorphs, dynamically
stable at 100 GPa, were grouped into two families. The first was derived from
the thermodynamic minimum at these conditions, $\alpha$-Ga, whereas channels
comprised the second. Two additional intergrowth structures were uncovered,
and it was shown they could be constructed by stacking layers of $\alpha$-Ga-
like and channel-like allotropes on top of each other. A detailed bonding
analysis revealed networks of four-center $\sigma$-bonding functions linked by
two-center B-B bonds in the $\alpha$-Ga based structures, and networks that
were largely composed of three-center $\sigma$-bonding functions in the
channel-based structures. Seven of these high pressure phases were found to be
metastable at atmospheric conditions, and their Vickers hardnesses were
estimated to be $\sim$36 GPa.

##  I Introduction

The potential for structural diversity in boron clusters [1, 2], sheets [3, 4]
and crystalline lattices [5, 6, 7], arises from boron’s propensity to adopt
novel multicenter bonding schemes to accomodate its three valence electrons.
Various stable and metastable polymorphs of elemental boron have been proposed
as superhard [8, 9, 10], superconducting [11, 12, 13, 14], and even
topological materials [15]. Furthermore, because of its higher density and
tensile strength as compared to plastics, this low-$Z$ material offers an
option as an ablator in inertial confinement fusion and high energy density
experiments, and its phase behavior as a function of temperature and pressure
is of extreme interest [16, 17]. The structural chemistry of boron, including
metastable phases that could be created in experiment, must be understood in
order to accurately model its behavior under such conditions.

Some of the synthesized forms of elemental boron are illustrated in Figure 1.
At ambient pressures two elemental phases, both based on packings of B12
icosahedra, are known. The unit cell of the $\alpha$-rhombohedral polymorph,
which can be synthesized at $T<1300$ K, contains a single icosahedron linked
to neighbors by classic two-center two-electron (2c-2e) bonds between polar
atoms, and three-center two-electron (3c-2e) bonds between equatorial atoms.
The idealized 105-atom unit cell of the $\beta$-rhombohedral polymorph, which
can be accessed at higher temperatures, contains an outer layer of 8
icosahedra at the corners and 12 icosahedra along the edges, with two fused
icosahedra in the middle connected by an interstitial boron atom. Multiple
computational studies have explored the relative energetics of these two
polymorphs [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. Much emphasis
has been placed on unravelling the role of defects on the stability of the
$\beta$ phase [6, 27, 28]. With these defects taken into account, first-
principles calculations identify it as being the most stable boron allotrope
at ambient pressure [24, 26, 27], although the $\alpha$ structure is close
enough to make experimental differentiation of their energetics intractable
[30]. Tetragonal boron allotropes have also been synthesized. It has been
questioned whether one of these, T-50 [31, 32, 33, 34], can be stable in its
pure form, whereas it is generally accepted that the other, T-192 [35, 36,
37], is the only allotrope other than $\alpha$\- or $\beta$\- rhombohedral
boron that is synthesizable at atmospheric pressures, although others have
been reported [38, 39, 40].

Figure 1: Select experimentally observed allotropes of boron. At ambient
pressures the (a) $\alpha$\- and (b) $\beta$-rhombohedral forms predominate,
with transitions to the (c) $\gamma$-orthorhombic, and (d) $\alpha$-Ga phases
at 19 and 89 GPa, respectively. The characteristic icosahedra of boron-based
structures (green: free standing; blue: face- and vertex-sharing) are evident
in the phases that are preferred at lower pressures.

Computations have predicted the following sets of structural phase transitions
for boron:
$\alpha\text{-B}_{12}\rightarrow\gamma\text{-B}_{28}\rightarrow\alpha\text{-Ga}$
at 19 and 89 GPa [41]. The orthorhombic unit cell of $\gamma\text{-B}_{28}$
features B12 icosahedra and B2 dumbells arranged in a rock-salt pattern. The
structure of this phase, which was first synthesized in 1965 by Wentorf [42],
was elucidated via experiment and theory more than 40 years after it was made
[43, 41]. The disputed question of electron transfer between the dumbbells and
icosahedra, and the possibility of ionic behavior in an elemental phase have
led to numerous studies of the electronic structure and bonding in this
allotrope [41, 9, 44, 45].

von Schnering and Nesper were the first to propose that boron could adopt a
$Cmca$ symmetry $\alpha$-Ga structure that does not contain any icosahedral
motifs [46]. First-principles calculations have verified the stability of this
phase under pressure [47, 48], and predicted its superconducting behavior
[13]. It is suspected that the synthesized superconducting boron phase whose
critical temperature, $T_{c}$, was measured to be as high as 11 K at 250 GPa
[12] corresponds to the $\alpha$-Ga structure. Recently, the synthesis of this
polymorph has been reported at 115 GPa and 2100 K, and the phase has been
dubbed $\zeta$-B [49]. We note that static compression of $\beta$-B105 leads
to amorphization at $\sim$100 GPa, not to the lower enthalpy $\alpha$-Ga
phase. This may be a result of the slow kinetics of boron phase transitions
[50]. Such kinetic effects are believed to be important in $\alpha$-B12 where
static compression of both doped and undoped samples showed a superconducting
phase transition without a structural transformation [51, 52, 53]. Theoretical
calculations predict a hitherto unobserved transformation to a non-icosahedral
$P6_{3}/mcm$ symmetry structure above 375 GPa, whose $T_{c}$ was estimated to
be 44 K at 400 GPa [14].

The Thomas-Fermi-Dirac model predicts that compression increasingly leads to
bonding delocalization and the formation of simpler crystal structures with
large coordination numbers. However, recent theoretical and experimental
discoveries have shown that the behavior of elements under pressure may be
much more complex [54], as in the case of lithium [55, 56], sodium [57, 58]
and aluminum [59]. Moreover, experimental variables such as the temperature
and compression mechanism (static vs. dynamic) can affect the phases that are
formed. For boron, in particular, discrepancies between the measured and
calculated shock Hugoniot, as well as the computed melting temperatures
suggest that new post-$\alpha$-Ga phases can be formed at megabar pressures
[17].

With this in mind, herein we employ evolutionary algorithms coupled with
first-principles calculations to predict metastable phases of boron at 100
GPa. Two families of allotropes with distinct structural and bonding motifs,
one based on derivatives of the $\alpha$-Ga structure and the other on phases
containing channels, whose enthalpies are within 100 meV/atom of the
thermodynamic minimum are found. The results suggest that many more related
allotropes of boron can be constructed based on simple building blocks
arranged and rearranged with different periodicity and alignments, as well as
intergrowths of these two structure types. In addition to $\alpha$-Ga, six of
the discovered polymorphs are metastable at atmospheric pressures, and their
estimated Vickers hardnesses are comparable to rhombohedral $\alpha$-B and
orthorhombic $\gamma$-B28[8, 60, 61].

##  II Computational Details

Structure searches were carried using the XtalOpt evolutionary algorithm (EA)
version r11 [62, 63] on unit cells containing 8, 12 and 20 boron atoms at 100
GPa. The initial set of random symmetric structures was produced via the
RandSpg algorithm [64], and duplicates were identified and removed from the
gene pool using the XtalComp algorithm [65]. Structures within 150 meV/atom of
the most stable phase found, $\alpha$-Ga, were fully optimized, and those
within 100 meV/atom (a total of 10 phases) were chosen for further analysis.
The structural coordinates, relative enthalpies, equations of states (EOS),
and other computed parameters are provided in the Supporting Information (SI).

Geometry optimizations and electronic structure calculations (densities of
states, band structures, electron localization functions [66]) were performed
via density functional theory (DFT) as implemented in the Vienna _Ab Initio_
Simulation Package (VASP) [67, 68] using the Perdew-Burke-Ernzerhof (PBE)
gradient-corrected exchange and correlation functional [69]. The projector
augmented wave (PAW) method [70] described the core states, and the boron
2s22p1 electrons were treated explicitly. In the structure searches the
geometries were optimized with the accurate precision setting in a two-step
process. First, the atomic positions were allowed to relax within a set unit
cell, followed by a full relaxation. A plane-wave basis set with an energy
cutoff of 600 eV, and a $\Gamma$-centered Monkhorst-Pack scheme [71] where the
number of divisions along each reciprocal lattice vector was chosen such that
its product with the real lattice constant was 30 Å was employed. These values
were increased to 800 eV and 50 Å for precise reoptimizations and electronic
structure calculations. Phonon densities of states and band structures were
obtained using the supercell approach as implemented in the PHONOPY code [72]
in concert with VASP.

Helmholtz free energy data under finite temperature conditions for the
$\alpha$-Ga and Channel-I phases, including the electron thermal and
vibrational contributions, were calculated under the quasiharmonic
approximation as implemented in the phonopy-qha script. The data at various
volume were fitted to a volume integrated third-order Birch-Murnaghan EOS
[73], which then gave the pressure and the Gibbs free energies were obtained
thereafter.

The LOBSTER package [74, 75] was used to calculate the crystal orbital
Hamilton population (COHP) [76], and the negative of the COHP integrated to
the Fermi level (-iCOHP) between select atom pairs. A bonding analysis was
performed using the Solid State Adaptive Natural Density Partitioning
(SSAdNDP) method [77], an extension of the AdNDP method [78] for molecules to
solids. A periodic implementation of the Natural Bond Orbital (NBO) technique
[79, 80] interfaced with SSAdNDP was used to project the plane-wave density
onto the 6-31G basis set. Calculations performed with the cc-PVTZ basis [81]
on a subset of the structures yielded similar results (see the SI). The
visualization module of the SSAdNDP code was employed to generate Gaussian
cube files for each bond, which were visualized using the VESTA software
package [82]. Powder X-ray diffraction patterns for each structure were
generated using the Diamond software package [83].

##  III Two Families of Metastable Boron Phases

The evolutionary searches carried out at 100 GPa identified $\alpha$-Ga as the
most stable phase, in agreement with previous theoretical studies [47, 48, 84,
41], and the recent synthesis of this phase in a diamond anvil cell at 115 GPa
and 2100 K [49]. Provided the energy or enthalpy of a metastable phase is not
too high, it could potentially be made by a judicious choice of the
temperature, heating method, pressure, and starting material, or by using
synthesis techniques to access phases far away from equilibrium [85]. A data
mining study of the Materials Project high throughput database, which contains
DFT energies computed at $T=0$ K of structures found in the Inorganic Crystal
Structure Database (ICSD), was recently carried out to quantify the
thermodynamic scale of metastability [86]. The 90th percentile of the DFT-
calculated metastability of all of the known inorganic crystalline materials
was found to be $\sim$70 meV/atom. However, examples of known polymorphs whose
energy was $\sim$150 meV/atom above the ground state illustrated that a low
energy of metastability did not necessarily correspond with synthesizability.
With this in mind, we chose to analyze the phases identified by the EA search
whose enthalpies, computed within the static lattice approximation, were
within 100 meV/atom of the $\alpha$-Ga phase. Several recent noteable examples
of metastable materials that have been synthesized under pressure include
superconducting PHn [87, 88, 89, 90] and CSHx [91, 92, 93] phases.

Table I lists the space-groups and enthalpies (with and without zero-point-
energy, ZPE, corrections) with respect to the $\alpha$-Ga phase of the
structures chosen for further analysis. Phonon calculations confirmed they
were dynamically stable at 100 GPa. Generally speaking, phases that were
related to $\alpha$-Ga had a lower ZPE than the channel structures, with the
$\alpha$-Ga having the second lowest ZPE and the Channel-I the highest. The
quasiharmonic approximation was used to calculate the relative Gibbs free
energies of the $\alpha$-Ga and Channel-I phases (see Fig S17-S20). These
calculations revealed that finite temperature conditions stabilize the
Channel-I structure more than the $\alpha$-Ga phase, although the $\alpha$-Ga
phase remains thermodynamically preferred at 100 GPa up to 1000 K. However,
the Channel-I phase (to be described below) is within $k_{B}T$, where $k_{B}$
is the Boltzmann constant, of $\alpha$-Ga at room temperature. While distinct
from one another, a number of morphological connections could be drawn between
these phases, resulting in two main families – those based on $\alpha$-Ga, and
those featuring channels. Additionally, two phases arising from the
intergrowth of the $\alpha$-Ga and channel structures were identified.

Table 1: Enthalpies (meV/atom) of the ten most stable boron allotropes predicted at 100 GPa with and without zero-point-energy (ZPE) corrections relative to the enthalpy of the $\alpha$-Ga phase.  Structure | Space group | $\Delta H$ |  $\Delta H$+ZPE
---|---|---|---
$\alpha$-Ga | $Cmca$ | 0.0 | 0.0
Channel-I | $C2/m$ | 12.5 | 15.1
Channel-II | $C2/m$ | 45.8 | 47.5
Intergrowth-I | $C2/m$ | 51.1 | 52.7
Channel-III | $C2/m$ | 61.9 | 63.2
$\alpha$-Ga’-I | $P\bar{1}$ | 63.9 | 65.7
$\alpha$-Ga’-II | $P\bar{1}$ | 70.8 | 71.9
$\alpha$-Ga’-III |  $P4_{2}2_{1}$2 | 77.1 | 76.6
Intergrowth-II | $P\bar{1}$ | 85.3 | 85.9
Channel-IV | $Ima2$ | 89.2 | 90.2

The first group of structures is derived from the $\alpha$-Ga phase, which is
the thermodynamic minimum. As shown in Figure 2(a) this $Cmca$ symmetry
structure is composed of layers of sheared honeycomb nets stacked in an ABAB…
pattern along the $a$ axis. Figure 2(b) highlights another viewpoint: layers
of buckled triangular sheets of boron atoms that are stacked along the $c$
axis and connected by short B-B contacts measuring 1.63 Å comprise this phase.
Within the buckled sheets, the B-B distances are slightly longer, ranging from
1.67-1.75 Å.

Figure 2: Illustration of the $\alpha$-Ga structure of elemental boron
($\zeta$-B) [49]. It can be viewed as (a) a layered structure consisting of
sheared honeycomb nets, with thicker contacts along short B-B distances, and
(b) as layers of buckled trigonal nets connected by the short B-B contacts
highlighted in (a). The light green and dark green boron atoms do not lie in
the same plane.

The structural features of three of the metastable phases found, denoted as
$\alpha$-Ga’-I, II, and III, place them within the same family as $\alpha$-Ga.
As shown in Figure 3, the first of these polymorphs, with P$\bar{1}$ symmetry,
is constructed via an ABCABC… stacking of sheared honeycomb nets. The shearing
pattern differs from that of the parent $\alpha$-Ga structure, resulting in a
changed repeat period and a much lower symmetry. However, rotating the crystal
by 90∘ still reveals stacked layers of buckled triangular nets of boron, a
consistent feature in these $\alpha$-Ga-based structural derivatives. The
$\alpha$-Ga’-II phase, also of P$\bar{1}$ symmetry, consists of yet another
ABCABC… stacking pattern, this time not only of hexagons but also of pentagons
and octagons. Lastly, the $\alpha$-Ga’-III phase, with P42212 symmetry, is an
ABAB… stacking of sheared hexagonal layers very similar to that of the parent
$\alpha$-Ga structure. However, viewing it as layers of buckled triangular
boron nets stacked along the $c$ axis reveals that every other net is rotated
by 90∘ relative to the ones above and below, resulting in a 42 screw axis.

Figure 3: Boron allotropes related to $\alpha$-Ga. The (a) $\alpha$-Ga’-I, and
(b) $\alpha$-Ga’-II structures are based on stacks of sheared honeycomb nets
with different stacking orders and shearing patterns. The (c,d)
$\alpha$-Ga’-III phase can be directly derived from $\alpha$-Ga by rotating
every buckled layer by 90∘ from the layer above it.

One could imagine multiple other derivatives of the $\alpha$-Ga structure
based on different shearing or stacking patterns of hexagonal nets; indeed,
several more such polymorphs were produced by XtalOpt, although out of the
range of enthalpies we considered. Moreover, the $Imma$ symmetry $o$-B16 phase
proposed by Fan et al. [84] is built from sheared honeycomb nets. At 100 GPa
we calculate its enthalpy as being 33 meV/atom higher than that of
$\alpha$-Ga, in-line with the structures presented here, although it was not
found in our EA searches.

The second family of phases found in the EA searches features large open
channels. The three most stable members of this group, dubbed Channel-I, II,
and III are shown in Figure 4. As in the case of $\alpha$-Ga and its
derivatives, they can be viewed as stacked layers of buckled boron nets.
However, unlike the triangular nets in the $\alpha$-Ga derivatives,
parallelograms and pairs of triangles join to form diamond-like motifs that
comprise the nets. Instead of being connected by B-B bonds, as in the
$\alpha$-Ga family, the buckled layers in the channel structures are joined by
strips that also contain these parallelogram and diamond motifs, as shown in
the insets. In the Channel-I phase, two open parallelogram units are adjacent
in these connecting strips (Figure 4(a), inset). The Channel-II phase can be
viewed as a distorted version of the Al network comprising the $I4/mmm$ BaAl4
structure type [94]. In fact, a recent theoretical study suggested that a
boron allotrope with this symmetry could be made by removing the Na atoms from
an $I4/mmm$ NaB4 phase that was predicted to be stable under pressure, and
quenching to 1 atm [11]. Placing this open-channel boron framework under
pressure might yield the channel-based structures highlighted here. Another
predicted compound that features similar, albeit larger channels based on
networks of B4 diamond and square motifs is a $Pnnm$-B16 phase whose
electronic structure reveals exotic topological states [95].

Figure 4: Channel-based family of boron structures that contain diamond-like
B4 motifs made from two edge-sharing triangles, along with more open
parallelograms (strips connecting buckled layers are shown in the insets),
with different arrangements producing different phases: (a) Channel-I, (b)
Channel-II, and (c) Channel-III. The buckled layers comprising these phases
are made from: two adjacent rows of diamond-like motifs separated by a single
row of parallelogram units in (d) Channel-I, or alternating strips of diamond-
like B4 units and parallelograms as in the (e) Channel-II, and (f) Channel-III
structures.

Changing patterns and repeat periods of strips of diamonds and squares in the
buckled nets give rise to the Channel-I, II, and III structures, which all
possess $C2/m$ symmetry. In comparison to these, the Channel-IV phase shown in
Figure 5, which assumes the $Ima2$ spacegroup, is more complex. Whereas the
buckled layers in the three $C2/m$ symmetry phases are connected by planar
strips composed of diamonds and squares, the connections in the Channel-IV
structure are crimped, resulting in a visually less open channel. As
illustrated in Figure 5(b), the Channel-IV phase can also be viewed as a
derivative of $\alpha$-Ga that is constructed from layers of hexagonal nets in
such a way that leads to the formation of channels.

Figure 5: (a) The Channel-IV structure is built up of corrugated strips of B
triangles (inset) connecting buckled layers of B atoms, leading to open
channels in the phase. (b) It can be alternately viewed as another derivative
of the $\alpha$-Ga structure built up of sheared hexagonal nets. (c) The
buckled layers of this structure reveal a complex pattern of diamond-like B4
and open parallelogram motifs.

Finally, the two structures shown in Figure 6, dubbed Intergrowth-I and
Intergrowth-II, which contained distinct $\alpha$-Ga and channel-like regions
were identified. In the first of these, with $C2/m$ symmetry, the two regions
of near equal size are intergrown along the $c$ axis. In Figure 6(a) regions
of the $\alpha$-Ga structure are drawn with thicker connections between atoms,
highlighting the sheared honeycomb motifs, whereas thinner lines connect boron
atoms within the channel-based regions (highlighted in grey; see inset). In
Intergrowth-I, a view of the buckled layers reveals that each boron atom is
linked in a purely triangular network. The open parallelograms seen in purely
channel-based structures are thus not present, due to the small spatial extent
of the channel-based regions in the intergrown phase. In Intergrowth-II, with
P$\bar{1}$ symmetry, the two regions are again intergrown along the $c$ axis,
although with the channel structure-based regions (highlighted in grey) being
much larger than in Intergrowth-I. Consequently, the distortions of the
buckled layers present in Intergrowth-II are similar to those seen in the
polymorphs that only contain channels. Additionally, as with the multiple
possible adaptations that could be imagined for the $\alpha$-Ga derivatives,
one can anticipate a plethora of possible intergrowths between the $\alpha$-Ga
and channel structures, with differing repeat periods and relative amounts of
each phase being incorporated. Indeed, one might imagine blurring the line
between an intergrowth structure as shown here and two crystals with
$\alpha$-Ga and channel morphology interfaced with one another.

Figure 6: (a) The Intergrowth-I structure, and one of its (b) buckled layers.
(c) The Intergrowth-II structure, and (d) one of its buckled layers. Both
phases are based on the $\alpha$-Ga and channel structures, with thicker bonds
highlighting $\alpha$-Ga-based regions. Grey backgrounds span channel-based
regions (shown in the insets), which are larger in the Intergrowth-II
structure.

##  IV Factors Affecting Stability

In the static lattice approximation at 0 K the thermodynamic variable
governing relative phase stability is the enthalpy, which is the sum of the
internal energy and the pressure volume term, given as $H=U+PV$. At high
pressures the magnitude of the $PV$ term, which favors denser structures, can
be greater than that of the $U$ term, which is a reflection of the bonding
present in the system. We therefore wondered if the geometric pecularities
between the two families of phases could be traced to different factors
affecting their stability? To explore the driving forces for stabilization we
calculated the aforementioned quantities for the nine metastable phases
predicted here, and compared them with those obtained for the $\alpha$-B12,
$\beta$-B105, and $\gamma$-B28 phases up to 300 GPa.

Figure 7: Enthalpies of select phases of boron, including those listed in
Table 1, with respect to that of $\alpha$-B12 as a function of pressure. (b)
The internal energy of these boron phases relative to $\alpha$-B12. (c) The
$PV$ contribution to the enthalpy of these boron phases relative to that of
the $\alpha$-Ga polymorph. Channel-based phases tend to have a lower internal
energy, while those based on $\alpha$-Ga have smaller $PV$ contributions to
the enthalpy.

For ease of visualization, the enthalpies are plotted relative to that of the
$\alpha$-B12 phase in Figure 7(a). At 0 GPa $\alpha$-B12 and a model of the
$\beta$-B105 structure wherein all of the atomic positions are fully occupied
are within 26 meV/atom of each other. $\alpha$-B12 remains the thermodynamic
minimum until $\sim$19 GPa, at which point $\gamma$-B28 becomes the
thermodynamic minimum. The $\gamma$-B28 polymorph persists until 89 GPa, when
the $\alpha$-Ga phase becomes lower in enthalpy [41]. Remarkably, until
$\sim$70 GPa, the Channel-I phase is actually lower in enthalpy than
$\alpha$-Ga, although neither is more stable than $\gamma$-B28. At higher
pressures the $\gamma$-B28 phase becomes progressively more destabilized as
compared to the polymorphs found within the EA search, a reflection of the
poor space-filling abilities of icosahedra.

Examination of Figure 7(b) illustrates that the stability of the Channel-I
structure can be linked to its low internal energy. With increasing pressure
the internal energies of $\gamma$-B28, $\alpha$-Ga, as well as the phases
generated in the EA search approach and even begin to overtake $\alpha$-B12.
Generally speaking, the channel polymorphs tend to have lower internal
energies than the $\alpha$-Ga derivatives, with those of the intergrowths
being intermediate between the two. Meanwhile, as shown in Figure 7(c), the
$\alpha$-Ga structure has the lowest $PV$ contribution to the enthalpy, with
the $\alpha$-B12, $\gamma$-B28, and $\beta$-B105 having the highest values,
respectively, reflecting the low packing density achieved by icosahedra.
Unsurprisingly, the $\alpha$-Ga derivatives tend to be denser than the channel
phases, with the intergrowth phases being intermediate.

Thus, the low enthalpies of the channel structures can be traced back to their
large negative internal energies, which are offset by the large positive $PV$
terms. In comparison, $\alpha$-Ga and its derivatives have smaller volumes,
but the magnitude of their internal energies are also smaller. The Channel-IV
allotrope, whose $PV$ term is less than that of both the Intergrowth-II and
$\alpha$-Ga’-II phases, is an exception. As previously discussed, although
this polymorph resembles the channel-based structures, it can also be
considered as an $\alpha$-Ga derivative. Its relatively low $PV$ term and high
internal energy are in-line with the latter interpretation, however a closer
look (Figure 5(a)) shows that, rather than the planar chains of square-and-
diamond boron motifs connecting buckled layers of boron, the chains running
along the $c$ axis are themselves corrugated, leading to a smaller relative
volume.

##  V Bonding Analysis

The prevalence of peculiar bonding schemes in boron-based materials [96, 46,
21, 97, 98, 99, 100] including the $\gamma$-B28 phase [44, 45, 41], which
frequently involve multicenter bonds, led us to wonder what sorts of curious
features were present in the discovered polymorphs, and if their sorting into
structural families heralded similar bonding arrangements? To investigate this
we began by calculating the Electron Localization Function (ELF) for each
polymorph; a summary of the results obtained for the $\alpha$-Ga, Channel-I,
and Channel-II phases at 100 GPa is presented in Figure 8. For $\alpha$-Ga,
ELF maxima are observed along the short 1.63 Å B-B contacts connecting the
buckled nets that lie in the $ac$ plane, suggestive of classic two-center two-
electron (2c-2e) bonding regions. Such maxima are also present in the ELF plot
for the Channel-I phase along B-B contacts 1.62 Å apart, on the edges shared
by two parallelograms within the chain of parallelograms and diamonds that
connect the buckled nets in the $bc$ plane. ELF maxima are also present within
the triangles comprising the B4 diamond-like motifs, and along the buckled
nets, although these are spread over a larger region. Finally, the ELF maxima
within the Channel-II structure are located within the B4 diamonds, similar to
the Channel-I phase, although no signatures of 2c-2e bonding interactions are
detected. As shown in the plots provided in the SI (Fig. S8-S14) these general
features are present in the remaining structures, with the $\alpha$-Ga
derivatives displaying ELF maxima indicative of 2c-2e bonds formed along the
short contacts connecting buckled nets, and the maxima in the channel-based
compounds arising within the diamond motifs.

Figure 8: Isosurface (isovalue = 0.8) plots of the ELF for the (a)
$\alpha$-Ga, (b) Channel-I, and (c) Channel-II structures at 100 GPa. Plots of
the ELF in the (d) (100) plane of $\alpha$-Ga, (e) (010) plane of Channel-I,
and (f) (010) plane of Channel-II are also provided.

To gain a deeper understanding of the bonding features present in these
systems, we used the SSAdNDP bonding analysis [77], an extension of the NBO
technique. Within this method, the plane-wave DFT electronic structure is
projected onto a Natural Atomic Orbital basis and the resulting density matrix
is used to locate putative bonding functions with occupancies above a user-
defined cutoff. In this way, a chemically meaningful, localized bonding scheme
can be derived containing lone pairs and classical 2c-2e bonds, as well as the
more complex multicenter bonding motifs that may be required for an accurate
description of the bonding in extended systems. Indeed, a hypothetical planar
boron allotrope, the $\alpha$-sheet, was used as a proof of concept for the
applicability of the SSAdNDP analysis, with bonding functions ranging from
3c-2e to 7c-2e being required to account for all valence electrons of the
system [77]. In addition, a novel phase consisting of stacked wiggle
$\alpha$-sheets, 3D-$\alpha$’ boron, was proposed to be the first three-
dimensional topological boron structure, with a spindle nodal chain arising
from the intersection of nodal lines and rings in momentum space [15]. Its
curious electronic structure was linked to the same sorts of $\pi$-bonds
present in the 2D $\alpha$-sheet analogue.

Let us begin with the simplest boron allotrope considered herein, the $Cmca$
$\alpha$-Ga phase whose unit cell contains eight boron atoms, resulting in a
total of 24 valence electrons to be accounted for. The SSAdNDP analysis
assigns eight of these to four 2c-2e $\sigma$-bonding functions with
occupation numbers (ONs) of 1.77 $|$e$|$ that connect the buckled nets, as
shown in Figure 9(a). These functions span the same region of space as the ELF
maxima, reflecting the high degree of localization in these classical bonds.
The remainder of the electrons are apportioned to eight 4c-2e $\sigma$-bonds
with ON = 1.92 $|$e$|$, see Figure 9(b), that span the buckled nets present in
the $\alpha$-Ga phase. Thus, each boron atom takes part in one 2c-2e
$\sigma$-bond with a boron atom on a neighboring buckled net, and four 4c-2e
$\sigma$-bonds that are shared among the surrounding eight boron atoms in the
buckled net. Each boron atom can then be assigned one electron from the 2c-2e
function and half an electron from each of the 4c-2e functions, yielding three
electrons in total, in-line with what is expected for a phase of boron.

Figure 9: Isosurfaces (isovalue=0.08 a.u.) of the bonding functions obtained
with the SSAdNDP method for $\alpha$-Ga at 100 GPa. (a) Four 2c-2e
$\sigma$-bonds between buckled layers of boron atoms, and (b) eight 4c-2e
$\sigma$-bonds within buckled layers of boron atoms account for all 24 valence
electrons in the structure.

Figure 10 plots the difference between the converged electron density from a
self-consistent calculation and that of the noninteracting atoms for the
$\alpha$-Ga, Channel-I and Channel-II structures. The resulting map is
positive in regions where electron density accumulates in proceeding from the
non-interacting picture to the density in the compound, that is, where bonds
might be expected. The results for $\alpha$-Ga are in-line with the ELF, and
also with the SSAdNDP analysis: significant electron accumulation is observed
along the B-B contacts that correspond to 2c-2e bonds, whereas more
delocalized regions are found along the B4 diamond like motifs, reflective of
the 4c-4e bonds found to span the buckled net.

Figure 10: Electron density difference isosurface plots (positive isovalue of
0.015 a.u. is shown) for the (a) $\alpha$-Ga, (b) Channel-I, and (c) Channel-
II structures. Bottom: buckled layers; top: chains connecting them.

Next we turn our attention to the bonding framework of the Channel-I and
Channel-II polymorphs. In the former, the ELF and density difference analysis
both indicate the presence of 2c-2e bonds in the chains running perpendicular
to and connecting the buckled layers. The SSAdNDP analysis confirms this
characterization, with two such bonds of ON = 1.7 $|$e$|$ per unit cell, as
seen in Figure 11. The diamond-like fragments comprising the remainder of
these chains are spanned by four 3c-2e bonds with ON = 1.9 $|$e$|$, with two
of these found in each B4 diamond motif, so that twelve electrons in total are
allocated to the chains. In the buckled layer, strips of triangular nets
alternate with strips of corrugated open squares. The bonding can be
decomposed into eight 3c-2e bonds per unit cell with ON = 1.9 $|$e$|$ that lie
along the edges of the triangular net strip, and four 4c-2e bonds per unit
cell with ON = 1.8 $|$e$|$ that lie in the middle, totalling 24 electrons.

The Channel-I structure is peculiar among the channel phases in the placement
of two open parallelogram units adjacent to one another in the chains
connecting the buckled layers, which is then spanned by a two-center bond.
Consequently, the local coordination of the boron atoms in the two-center bond
is very similar to that of the boron atoms in the $\alpha$-Ga structure, which
each participate in one two-center bond connecting buckled layers and four
four-center bonds in the buckled layers. The atoms involved in the two-center
bond in the Channel-I phase participate in four- and three-center bonds in the
buckled layers. Additional repetition of adjacent open parallelogram units
could extend this network of boron atoms involved in two-center bonds between
layers and four-center bonds within layers, leading to the bonding scheme
exhibited by the $\alpha$-Ga based phases, rather than the three-center bonds
found in the channel-based structures.

Figure 11: Isosurfaces (isovalue=0.08 a.u.) of the bonding functions obtained
for the Channel-I structure of boron at 100 GPa. (a) Two 2c-2e and four 3c-2e
$\sigma$-bonds between buckled layers of boron atoms, and (b) eight 3c-2e and
four 4c-2e $\sigma$-bonds within buckled layers of boron atoms account for all
36 valence electrons in the phase.

In the less stable Channel-II phase the SSAdNDP analysis, whose results are
plotted in Figure 12, reveals that the chains linking buckled layers composed
of diamond-like fragments can be described by 3c-2e bonds with ON = 1.9
$|$e$|$, and do not contain any 2c-2e bonds. Within the buckled layer, the
strips of triangular nets are narrower than in the Channel-I structure and are
spanned entirely by 3c-2e functions with ON = 1.8 $|$e$|$. With eight 3c-2e
functions in the buckled layers and four in the connecting chains, all 24
electrons in the unit cell are accounted for. In essence, this bonding scheme
matches the one derived for Channel-I but with the 2c-2e and 4c-2e functions
deleted.

Figure 12: Isosurfaces (isovalue=0.08 a.u.) of the bonding functions obtained
for the Channel-II structure of boron at 100 GPa. (a) Four 3c-2e
$\sigma$-bonds between buckled layers of boron atoms, and (b) eight 3c-2e
$\sigma$-bonds within buckled layers of boron atoms account for all 24 valence
electrons in the phase.

The prevalence of 3c-2e bonds in both channel phases can be traced to features
in their electron difference density maps (bottom of Figure 10). The map for
the $\alpha$-Ga phase reveals regions of electron accumulation in the buckled
layers that are spread across four B atoms in diamond motifs. However, the
difference density in the triangular nets of the channel structures, where
three-center bonding functions were detected by the SSAdNDP analysis, displays
greater concentration along the edges of the triangular regions rather than
being contained within B4 units as in the $\alpha$-Ga phase. This may reflect
the greater amount of electron localization within a region constrained by
three atoms rather than by four.

To further support the assignment of the bonding to 2c-2e, 3c-3e and 4c-4e
functions in the $\alpha$-Ga and two channel phases, we calculated the
negative of the crystal orbital Hamilton populations integrated to the Fermi
level, -iCOHP, which reflects the strength of a bond between select atom
pairs. The results, provided in Table 2, highlight the relative strengths of
interaction between pairs of boron atoms involved in 2c, 3c, or 4c bonds as
identified by the SSAdNDP analysis. B-B contacts of more than 2.0 Å were
neglected, as their -iCOHP values were an order of magnitude smaller than the
remaining pairs of atoms.

By far the largest -iCOHP value, at 7.59 eV/bond, corresponds to the 2c-2e
bond in the $\alpha$-Ga structure, with the second-largest value of 6.81
eV/bond arising from the two-center contact identified in the Channel-I phase.
Between pairs of boron atoms involved in 3-center bonding functions, the
-iCOHP values range from 6.58 eV/bond (in the functions connecting buckled
layers in the Channel-I phase) to 4.68 eV/bond (within the buckled layers in
Channel-I), with most of them between 5 to 6 eV/bond. While robust, these
interactions are not as strong as those between the pairs of boron atoms
involved in two-center bonds. Finally, regions containing four-center bonds
tended to have the weakest bond strengths, with the highest value of 5.38
eV/bond being found within the buckled layer of the $\alpha$-Ga structure and
most lying between 4 and 5 eV/bond. In this way, the -iCOHP analysis provides
additional delineation between the bonding motifs identified by SSAdNDP and
the difference density analyses.

Table 2: The crystal orbital Hamilton populations integrated to the Fermi level (-iCOHP) for select B-B contacts in the $\alpha$-Ga, Channel-I, and Channel-II phases at 100 GPa. The B-B contacts are labeled according to the type of bonding function they comprise according to the SSAdNDP analysis (2c, 3c, 4c, or between 3c and 4c regions). Location | Distance (Å) | -iCOHP (eV/bond)
---|---|---
$\alpha$-Ga |  |
Interlayer (2c) | 1.63 | 7.59
Buckled (4c) | 1.67 | 5.38
Buckled (4c) | 1.72 | 4.77
Buckled (4c) | 1.75 | 4.09
Channel-I |  |
Interlayer (2c) | 1.62 | 6.81
Interlayer (3c) | 1.58 | 6.58
Interlayer (3c) | 1.70 | 5.18
Interlayer (3c) | 1.72 | 4.68
Buckled (3c) | 1.66 | 5.84
Buckled (3c) | 1.63 | 5.32
Buckled (3c/4c) | 1.68 | 4.88
Buckled (4c) | 1.67 | 4.88
Buckled (4c) | 1.71 | 4.48
Channel-II |  |
Interlayer (3c) | 1.59 | 6.14
Interlayer (3c) | 1.68 | 5.01
Interlayer (3c) | 1.70 | 5.11
Interlayer (3c) | 1.64 | 5.11
Buckled (3c) | 1.66 | 5.66
Buckled (3c) | 1.67 | 4.77

The SSAdNDP analysis was carried out for the remaining structures, and the
results are summarized in Table 3 with further details provided in the SI
(Fig. S8-S14). The bonding in the $\alpha$-Ga’-I, II, and III phases can be
described by 4c-2e functions spanning the buckled triangular nets of boron,
with 2c-2e bonds bridging the buckled nets. The diamond-like motifs that
comprise the Channel-III structure are filled with 3c-2e bonds.The Channel-IV
polymorph, also interpretable as an $\alpha$-Ga derivative, contains two-,
three-, and four-center bonds in arrangements reflecting its relationships to
both $\alpha$-Ga and channel polymorphs. Finally, the bonding within the two
intergrowths of alternating $\alpha$-Ga and channel-like layers is
characterized by two- and four-center networks in the $\alpha$-Ga based
regions, while the channel-based regions of the Intergrowth-I and -II
polymorphs contain four- and three-center bonds, respectively. The three-
center bonds in the channel-like regions of the Intergrowth-II phase are
analogous to the 3c bonding networks in the parent channel-based structures,
while the four-center bonds found in the Intergrowth-I structure reflect the
narrowness of the channel-like regions in that structure in comparison to
those based on $\alpha$-Ga, which contain four-center bonding functions in the
buckled nets.

Table 3: Summary of SSAdNDP analyses of boron phases at 100 GPa. Structure | Formula units | 2c-2e | 3c-2e | 4c-2e
---|---|---|---|---
$\alpha$-Ga | 8 | 4 | 0 | 8
Channel-I | 12 | 2 | 12 | 4
Channel-II | 8 | 0 | 12 | 0
Intergrowth I | 12 | 6 | 0 | 12
Channel-III | 24 | 0 | 36 | 0
$\alpha$-Ga’-I | 12 | 6 | 0 | 12
$\alpha$-Ga’-II | 8 | 4 | 0 | 8
$\alpha$-Ga’-III | 8 | 4 | 0 | 8
Intergrowth-II | 12 | 4 | 4 | 10
Channel-IV | 16 | 4 | 16 | 4

##  VI Properties at 1 Atmosphere

Phonon calculations were carried out to determine which boron allotropes could
conceivably be quenched to atmospheric pressures (effectively 0 GPa in our
calculations). In addition to the $\alpha$-Ga structure, the Channel-I and
-II, $\alpha$-Ga’-I and -II, and Intergrowth-I and -II phases were found to be
dynamically stable. Their enthalpies range from 199 meV/atom (Channel-I) to
287 meV/atom (Intergrowth-I) relative to the $\alpha$-B12 phase. Band
structure calculations within the PBE functional suggested that Channel-I and
Intergrowth-I were both small gap semiconductors, whereas the remainder were
metallic. However, the PBE functional is known to underestimate band gaps, and
inclusion of Hartree-Fock exchange would likely open up the band gap in some
of these.

Covalent compounds that contain boron atoms are known to behave as hard
materials. Therefore, the Vickers Hardness, $H_{\text{v}}$, of these phases at
ambient pressure were obtained using the Teter equation [101, 102] combined
with shear moduli obtained via a machine-learning model [103] trained on the
materials within the AFLOW repository [104, 105]. This method was shown to
provide hardness estimates that correlated well with experiment for a broad
range of compounds [106]. The estimated $H_{\text{v}}$s were 39.6 GPa for
$\alpha$-Ga, 32.8 GPa for the $\alpha$-Ga’-I structure, and between 35.5 and
36.6 GPa for the remainder, all of which are slightly higher than the
$H_{\text{v}}$s of 30.5, 31.5 and 28.4 GPa obtained for the known
$\alpha$-B12, $\gamma$-B28 and $\beta$-B105 phases. The $m$-B16 and $o$-B16
allotropes predicted in Ref. [84], were estimated to be somewhat harder,
whereas the $I4/mmm$-B4 and $Pm$-B17 phases from Ref. [11] were not quite as
hard. However, it is well known that both the computational and experimental
methods employed to guage hardness are imperfect, with the former often
yielding values that can differ by more than 10% [107]. Therefore, we suggest
that the $H_{\text{v}}$s of these various allotropes of boron at 1 atm do not
differ significantly from each other.

##  VII Conclusions

Evolutionary crystal structure prediction searches were employed to identify a
series of metastable polymorphs of elemental boron at 100 GPa. Potentially
synthesizable phases, which were within 100 meV/atom of the thermodynamic
minimum at this pressure, the $\alpha$-Ga structure, were further analyzed.
These could be grouped into families based upon prevalant structural motifs:
three were clearly derived from the $\alpha$-Ga phase, three were based on
structures that contained channels of boron atoms, and two were layered
intergrowths of the $\alpha$-Ga structure with the channel-based phases. One
additional structure bears a clear resemblance with the channel-based
compounds, but can also be interpreted as another derivative of $\alpha$-Ga
based on layers of sheared honeycomb nets.

An examination of the factors stabilizing these allotropes revealed that the
$PV$ contribution to the enthalpy favored the denser $\alpha$-Ga-derived
structures, while the internal energies of the channel-based phases were
lower. We further found that finite temperature contributions stabilize the
Channel-I phase more than the $\alpha$-Ga structure. This effect is magnified
at lower pressures owing to the changes in the relative $PV$ and Helmholtz
free energy terms, suggesting that a viable route to synthetically accessing
the Channel-I structure could be a high-pressure high-temperature synthesis
followed by quenching to lower-pressure conditions.

To further understand the bonding in these phases we employed the Solid State
Adaptive Natural Density Partitioning (SSAdNDP) method, coupled with
calculations of the electron localization functions (ELF), crystal orbital
Hamilton populations (COHP), and density difference analysis. In the
$\alpha$-Ga structure, four-center two-electron (4c-2e) bonds were found to
span the buckled layers with classic two-center two-electron (2c-2e) bonds
connecting these layers. All of the $\alpha$-Ga derived phases were
characterized by these same types of 4c-2e and 2c-2e bonds. The channel-based
polymorphs, meanwhile, were largely characterized by three-center two-electron
(3c-2e) bonds distributed across the diamond-like B4 motifs present in the
phases. The intergrowth structures contained features from both $\alpha$-Ga
and channel-based parents. The $\alpha$-Ga structure, along with six of the
predicted structures remained dynamically stable at atmospheric pressures,
with estimated Vickers hardnesses of $\sim$36 GPa.

The plethora of known and hypothesized allotropes of boron reflect the
complexity of its crystal chemistry, suggesting a wealth of materials for
further study under appropriate experimental conditions. Furthermore, the
structural relations revealed in this study indicate a possibility of numerous
additional metastable variants of elemental boron based on minor modifications
or intergrowths of a few basic parent structures. We look forward to further
experimental and theoretical explorations of the complex crystal chemistry of
elemental boron, in particular, the twin roles of pressure and temperature in
the stability of various allotropes, including the effects of static versus
dynamic compression.

##  VIII Acknowledgments

K.H. is thankful to the U.S. Department of Energy, National Nuclear Security
Administration, through the Capital-DOE Alliance Center under Cooperative
Agreement DE-NA0003975 for financial support. Calculations were performed at
the Center for Computational Research at SUNY Buffalo [108]. Part of this work
was performed under the auspices of the U.S. Department of Energy by Lawrence
Livermore National Laboratory under Contract DE-AC52-07NA27344. T.O.notes
LLNL-JRNL-817818.

## References

  * Zhai  _et al._ [2003] H. J. Zhai, A. N. Alexandrova, K. A. Birch, A. I. Boldyrev, and L. S. Wang, Hepta- and octacoordinate boron in molecular wheels of eight- and nine-atom boron clusters; observation and confirmation, Angew. Chem. Int. Ed. 42, 6004 (2003).
  * Alexandrova  _et al._ [2006] A. N. Alexandrova, A. I. Boldyrev, H. J. Zhai, and L. S. Wang, All-boron aromatic clusters as potential new inorganic ligands and building blocks in chemistry, Coord. Chem. Rev. 250, 2811 (2006).
  * Li  _et al._ [2019] D. Li, J. Gao, P. Cheng, J. He, Y. Yin, Y. Hu, L. Chen, Y. Cheng, and J. Zhao, 2d boron sheets: Structure, growth, and electronic and thermal transport properties, Adv. Funct. Mater. , 1904349 (2019).
  * Sun  _et al._ [2017] X. Sun, X. Liu, J. Yin, J. Yu, Y. Li, Y. Hang, X. Zhou, M. Yu, J. Li, G. Tai, and W. Guo, Two-dimensional boron crystals: Structural stabilitiy, tunable properties, fabrications and applications, Adv. Funct. Mater. , 1603300 (2017).
  * Albert and Hillebrecht [2009] B. Albert and H. Hillebrecht, Boron: Elementary challenge for experimenters and theoreticians, Angew. Chem. Int. Ed. 48, 8640 (2009).
  * Ogitsu  _et al._ [2013] T. Ogitsu, E. Schwegler, and G. Galli, $\beta$-rhombohedral boron: At the crossraods of the chemistry of boron and the physics of frustration, Chem. Rev. 113, 3425 (2013).
  * Shirai [2017] K. Shirai, Phase diagram of boron crystals, Jpn. J. Appl. Phys. 56, 05FA06 (2017).
  * Veprek  _et al._ [2011] S. Veprek, R. Zhang, and A. Argon, Mechanical properties and hardness of boron and boron-rich solids, J. Superhard Mater. 33, 409 (2011).
  * Zarechnaya  _et al._ [2009a] E. Y. Zarechnaya, L. Dubrovinsky, N. Dubrovinskaia, Y. Filinchuk, D. Chernyshov, V. Dmitriev, N. Miyajima, A. El Goresy, H. F. Braun, S. Van Smaalen, I. Kantor, A. Kantor, V. Prakapenka, M. Hanfland, A. S. Mikhaylushkin, I. A. Abrikosov, and S. I. Simak, Superhard semiconducting optically transparent high pressure phase of boron, Phys. Rev. Lett 102, 185501 (2009a).
  * Zhou  _et al._ [2010] W. Zhou, H. Sun, and C. Chen, Soft bond-deformation paths in superhard $\gamma$-boron, Phys. Rev. Lett. 105, 215503 (2010).
  * Zhang  _et al._ [2020a] S. Zhang, X. Du, J. Lin, A. Bergara, X. Chen, X. Liu, X. Zhang, and G. Yang, Superconducting boron allotropes, Phys. Rev. B 101, 174507 (2020a).
  * Eremets  _et al._ [2001] M. Eremets, V. Struzhkin, H.-k. Mao, and R. Hemley, Superconductivity in boron, Science 293, 272 (2001).
  * Ma  _et al._ [2004] Y. Ma, J. Tse, D. Klug, and R. Ahuja, Electron-phonon coupling of $\alpha$-ga boron, Phys. Rev. B 70, 214107 (2004).
  * Li  _et al._ [2014] D. Li, K. Bao, F. Tian, X. Jin, D. Duan, Z. He, B. Liu, and T. Cui, Hih-pressure close-packed structure of boron, RSC Adv. 4, 203 (2014).
  * Gao  _et al._ [2018] Y. Gao, Y. Xie, Y. Chen, J. Gu, and Z. Chen, Spindle nodal chain in three-dimensional $\alpha$’ boron, Phys. Chem. Chem. Phys. 20, 23500 (2018).
  * Zhang  _et al._ [2018a] S. Zhang, B. Militzer, M. C. Gregor, K. Caspersen, L. H. Yang, J. Gaffney, T. Ogitsu, D. Swift, A. Lazicki, D. Erskine, R. A. London, P. M. Celliers, J. Nilsen, P. A. Sterne, and H. D. Whitley, Theoretical and experimental investigation of the equation of state of boron plasmas, Phys. Rev. E 98, 023205 (2018a).
  * Zhang  _et al._ [2020b] S. Zhang, H. Whitley, and T. Ogitsu, Phase transformation in boron under shock compression, Solid State Sci. 108, 106376 (2020b).
  * Masago  _et al._ [2006] A. Masago, K. Shirai, and H. Katayama-Yoshida, Crystal stability of $\alpha$\- and $\beta$-boron, Phys. Rev. B. 73, 104102 (2006).
  * Masago  _et al._ [2004] A. Masago, K. Shirai, and H. Katayama-Yoshida, The pressure dependence of solid boron, Molecular Simulation 30, 935 (2004).
  * Masago  _et al._ [2005] A. Masago, K. Shirai, and H. Katayama-Yoshida, Crystal stability of $\alpha$\- and $\beta$-boron, Physics of Semiconductors, Pts A and B 772, 87 (2005).
  * Prasad  _et al._ [2005] D. L. V. K. Prasad, M. M. Balakrishnarajan, and E. Jemmis, Electronic structure and bonding of $\beta$-rhombohedral boron using cluster fragment approach, Phys. Rev. B 72, 195102 (2005).
  * Shang  _et al._ [2007] S. Shang, Y. Wang, R. Arroyave, and Z.-K. Liu, Phase stability in alpha- and beta-rhombohedral boron, Phys. Rev. B 75, 092101 (2007).
  * Siberchicot [2009] B. Siberchicot, Ab initio equation of state of alpha- and beta-boron: Possible amorphization of beta-boron under high pressure, Phys. Rev. B 79, 224101 (2009).
  * van Setten  _et al._ [2007] M. J. van Setten, M. A. Uijttewaal, G. A. de Wijs, and R. A. de Groot, Thermodynamic stability of boron: The role of defects and zero point motion, J. Am. Chem. Soc. 129, 2458 (2007).
  * M. Widom [2008a] M. M. M. Widom, Thermodynamic stability of boron: The role of defects and zero point motion, 16th International Symposium on Boron, Borides, and Related Materials 176, 012024 (2008a).
  * M. Widom [2008b] M. M. M. Widom, Symmetry-broken crystal structure of elemental boron at low temperature, Phys. Rev. B 77, 064113 (2008b).
  * Ogitsu  _et al._ [2009] T. Ogitsu, F. Gygi, J. Reed, Y. Motome, E. Schwegler, and G. Galli, Imperfect crystal and unusual semiconductor: Boron, a frustrated element, J. Am. Chem. Soc. 131, 1903 (2009).
  * Ogitsu  _et al._ [2010] T. Ogitsu, F. Gygi, J. Reed, M. Udagawa, Y. Motome, E. Schwegler, and G. Galli, Geometrical frustration in an elemental solid: An ising model to explain the defect structure of beta -rhombohedral boron, Phys. Rev. B 81, 020102 (2010).
  * Ogitsu and Schwegler [2012] T. Ogitsu and E. Schwegler, The $\alpha$–$\beta$ phase boundary of elemental boron, Solid State Sci. 14, 1598 (2012).
  * White  _et al._ [2015] M. A. White, A. B. Cerqueira, C. A. Whitman, M. B. Johnson, and T. Ogitsu, Determination of phase stability of elemental boron, Angew. Chem. Int. Ed. 127, 3697 (2015).
  * Laubengayer  _et al._ [1943] A. W. Laubengayer, D. T. Hurd, A. E. Newkirk, and J. L. Hoard, Boron. i. preparation and properties of pure crystalline boron, J. Am. Chem. Soc. 65, 1924 (1943).
  * Hoard  _et al._ [1951] J. L. Hoard, S. Geller, and R. E. Hughes, On the structure of elementary boron, J. Am. Chem. Soc. 73, 1892 (1951).
  * Hoard  _et al._ [1958] J. L. Hoard, R. E. Hughes, and D. E. Sands, The structure of tetragonal boron, J. Am. Chem. Soc. 80, 4507 (1958).
  * Hoard and Newkirk [1960] J. L. Hoard and A. E. Newkirk, An analysis of polymorphism in boron based upon x-ray diffraction results, J. Am. Chem. Soc. 82, 70 (1960).
  * Talley [1960] C. P. Talley, A new polymorph of boron, Acta Crystallogr. 13, 271 (1960).
  * Vlasse  _et al._ [1979a] M. Vlasse, R. Naslain, J. S. Kasper, and K. Ploog, Crystal-structure of tetragonal boron, J. Less-Common Met. 67, 1 (1979a).
  * Vlasse  _et al._ [1979b] M. Vlasse, R. Naslain, J. S. Kasper, and K. Ploog, Crystal-structure of tetragonal boron related to alpha-b12, J. Solid State Chem. 28, 289 (1979b).
  * Parakhonskiy  _et al._ [2013] G. Parakhonskiy, N. Dubrovinskaia, E. Bykova, R. Wirth, and L. Dubrovinsky, High pressure synthesis and investigation of single crystals of metastable boron phases, High Pressure Res. 33, 673 (2013).
  * Karmodak and Jemmis [2017] N. Karmodak and E. D. Jemmis, Fragment approach to the electronic structure of $\tau$-boron allotrope, Phys. Rev. B 95, 165128 (2017).
  * An  _et al._ [2016] Q. An, M. Reddy, K. Y. Xie, K. J. Hemker, and W. A. Goddard III, New ground state crystal structure of elemental boron, Phys. Rev. Lett. 117, 085501 (2016).
  * Oganov  _et al._ [2009] A. Oganov, J. Chen, C. Gatti, Y. Ma, Y. Ma, C. Glass, Z. Liu, T. Yu, O. Kurakevych, and V. Solozhenko, Ionic high-pressure form of elemental boron, Nature 457, 863 (2009).
  * Wentorf [1965] R. Wentorf, Boron: Another form, Science 147, 49 (1965).
  * Zarechnaya  _et al._ [2009b] E. Zarechnaya, L. Dubrovinsky, N. Dubrovinskaia, N. Miyajima, Y. Filinchuk, D. Chernyshov, and V. Dmitriev, Synthesis of an orthorhombic high pressure boron phase, Sci. Technol. Adv. Mater. 9, 044209 (2009b).
  * Häussermann and Mikhaylushkin [2010] U. Häussermann and A. Mikhaylushkin, Structure and bonding of $\gamma$-b28: Is the high pressure form of elemental boron ionic?, Inorg. Chem. 49, 11270 (2010).
  * Mondal  _et al._ [2011] S. Mondal, S. van Smaalen, A. Schönleber, Y. Filinchuk, D. Chernyshov, S. I. Simak, A. S. Mikhaylushkin, I. A. Abrikosov, E. Zarechnaya, L. Dubrovinsky, and N. Dubrovinskaia, Electron-deficient and polycenter bonds in the high-pressure $\gamma$-b28 phase of boron, Phys. Rev. Lett. 106, 215502 (2011).
  * von Schnering and Nesper [1991] H. von Schnering and R. Nesper, $\alpha$-gallium: An alternative to the boron structure, Acta Chem. Scand. 45, 870 (1991).
  * Häussermann  _et al._ [2003] U. Häussermann, S. Simak, R. Ahuja, and B. Johansson, Metal-nonmetal transition in the boron group elements, Phys. Rev. Lett. 90, 065701 (2003).
  * Segall and Arias [2003] D. Segall and T. Arias, Ab initio approach for high-pressure systems with application to high-pressure phases of boron: Perturbative momentum-space potentials, Phys. Rev. B 67, 064105 (2003).
  * Chuvashova  _et al._ [2017] I. Chuvashova, E. Bykova, M. Bykov, V. Prakapenka, K. Glazyrin, M. Mezour, L. Dubrovinsky, and N. Dubrovinskaia, Nonicosahedral boron allotrope synthesized at high pressure and high temperature, Phys. Rev. B 95, 180102 (2017).
  * Sanz  _et al._ [2002] D. N. Sanz, P. Loubeyre, and M. Mezouar, Equation of state and pressure induced amorphization of beta -boron from x-ray measurements up to 100 gpa, Phys. Rev. Lett 89, 245501 (2002).
  * Shirai  _et al._ [2009] K. Shirai, H. Dekura, and A. Yanase, Electronic structure and electrical resistivity of alpha-boron under high pressure, J.Phys. Soc. Jpn. 78 (2009).
  * Shirai  _et al._ [2011] K. Shirai, H. Dekura, Y. Mori, Y. Fujii, H. Hyodo, and K. Kimura, Structural study of alpha-rhombohedral boron at high pressures, J.Phys. Soc. Jpn. 80 (2011).
  * Nagatochi  _et al._ [2011] T. Nagatochi, H. Hyodo, A. Sumiyoshi, K. Soga, Y. Sato, M. Terauchi, F. Esaka, and K. Kimura, Superconductivity in li-doped $\alpha$-rhombohedral boron, Phys. Rev. B 83, 184507 (2011).
  * Miao  _et al._ [2020] M. Miao, Y. Sun, E. Zurek, and H. Lin, Chemistry under high pressure, Nat. Rev. Chem. 508-527, 4 (2020).
  * Marques  _et al._ [2011] M. Marques, M. I. McMahon, E. Gregoryanz, M. Hanfland, C. L. Guillaume, C. J. Pickard, G. J. Ackland, and R. J. Nelmes, Crystal structures of dense lithium: A metal–semiconductor-metal transition, Phys. Rev. Lett. 106, 095502 (1 (2011).
  * Tsuppayakorn-aek  _et al._ [2018] P. Tsuppayakorn-aek, W. Luo, T. Watcharatharapong, R. Ahuja, and T. Bovornratanaraks, Structural prediction of host-guest structure in lithium at high pressure, Sci. Rep. 8, 5278 (2018).
  * Gregoryanz  _et al._ [2008] E. Gregoryanz, L. F. Lundegaard, M. I. McMahon, C. Guillaume, R. J. Nelmes, and M. Mezouar, Structural diversity of sodium, Science 320, 1054 (2008).
  * Ma  _et al._ [2009] Y. M. Ma, M. Eremets, A. R. Oganov, Y. Xie, I. Trojan, S. Medvedev, A. O. Lyakhov, M. Valle, and V. Prakapenka, Transparent dense sodium, Nature 458, 182 (2009).
  * Pickard and Needs [2010] C. J. Pickard and R. J. Needs, Aluminum at terapascal pressures, Nat. Mater. 9, 624 (2010).
  * Solozhenko  _et al._ [2008] V. Solozhenko, O. Kurakevych, and A. Oganov, On the hardness of a new boron phase, orthorhombic $\gamma$-b28, J. Superhard Mater. 30, 428 (2008).
  * Zhang  _et al._ [2018b] S. H. Zhang, X. Zheng, Q. Q. Jin, S. J. Zheng, D. Legut, X. H. Yu, H. Y. Gou, Z. H. Fu, Y. Q. Guo, B. M. Yan, C. Peng, C. Q. Jin, T. C. Germann, and R. F. Zhang, Unprecedented plastic flow channel in $\gamma$-b28 through ultrasoft bonds: A challenge to superhardness, Phys. Rev. Mater. 2, 123602 (2018b).
  * Lonie and Zurek [2011] D. Lonie and E. Zurek, Xtalopt: An open-souce evolutionary algorithm for crystal structure prediction, Comp. Phys. Commun. 182, 372 (2011).
  * Avery  _et al._ [2018] P. Avery, Z. Falls, and E. Zurek, Xtalopt version r11: An open-source evolutionary algorithm for crystal structure prediction, Comput. Phys. Commun. 222, 418 (2018).
  * Avery and Zurek [2017] P. Avery and E. Zurek, Randspg: An open-source program for generating atomistic crystal structures with specific spacegroups, Comput. Phys. Commun. 213, 208 (2017).
  * Lonie and Zurek [2012] D. C. Lonie and E. Zurek, Identifying duplicate crystal structures: Xtalcomp, an open-source solution, Comput. Phys. Commun. 183, 690 (2012).
  * Savin [2005] A. Savin, The electron localization function (elf) and its relatives: Interpretations and difficulties, J. Mol. Struc.-Theochem. 727, 127 (2005).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for Ab Initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13, 5188 (1976).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Birch [1952] F. Birch, J. Geophys Res. 57, 227 (1952).
  * Deringer  _et al._ [2011] V. L. Deringer, A. L. Tchougréeff, and R. Dronskowski, Crystal orbital hamilton population (cohp) analysis as projected from plane-wave basis sets, J. Phys. Chem. A 115, 5461 (2011).
  * Maintz  _et al._ [2016] S. Maintz, V. L. Deringer, A. L. Tchougréeff, and R. Dronskowski, Lobster: A tool to extract chemical bonding from plane-wave based dft, J. Comp. Chem. 37, 1030 (2016).
  * Dronskowski and Blöchl [1993] R. Dronskowski and P. E. Blöchl, Crystal orbital hamilton populations (cohp): energy-resolved visualization of chemical bonding in solids based on density-functional calculations, J. Phys. Chem. 97, 8617 (1993).
  * Galeev  _et al._ [2013] T. R. Galeev, B. D. Dunnington, J. R. Schmidt, and A. I. Boldyrev, Solid state adaptive natural density partitioning: a tool for deciphering multi-center bonding in periodic systems, Phys. Chem. Chem. Phys. 15, 5022 (2013).
  * Zubarev and Boldyrev [2008] D. Y. Zubarev and A. I. Boldyrev, Developing paradigms of chemical bonding: adaptive natural density partitioning, Phys. Chem. Chem. Phys. 10, 5207 (2008).
  * Weinhold and Landis [2005] F. Weinhold and C. R. Landis,  _Valency and Bonding: A Natural Bond Orbital Donor-Acceptor Perspective_ (Cambridge University Press, 2005).
  * Dunnington and Schmidt [2012] B. D. Dunnington and J. R. Schmidt, Generalization of natural bond orbital analysis to periodic systems: Applications to solids and surfaces via plane-wave density functional theory, J. Chem. Theory. Comput. 8, 1902 (2012).
  * Dunning [1989] J. T. H. Dunning, Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen, J. Chem. Phys. 90, 1007 (1989).
  * Momma and Izumi [2011] K. Momma and F. Izumi, Vesta 3 for three-dimensional visualization of crystal, volumetric, and morphology data, J. Appl. Crystallogr. 44, 1272 (2011).
  * [83] Diamond - crystal and molecular structure visualization.
  * Fan  _et al._ [2015] C. Fan, J. Li, and L. Wang, Phase transitions, mechanical properties, and electronic structures of novel boron phases under high-pressure: A first-principles study, Sci. Rep. 4, 6786 (2015).
  * Parija  _et al._ [2018] A. Parija, G. R. Waetzig, J. L. Andrews, and S. Banerjee, Traversing energy landscapes away from equilibrium: Strategies for accessing and utilizing metastable phase space, J. Phys. Chem. C 122, 25709 (2018).
  * Sun  _et al._ [2016] W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain, W. D. Richards, A. C. Gamst, K. A. Persson, and G. Ceder, The thermodynamic scale of inorganic crystalline metastability, Sci. Adv. 2, e1600225 (2016).
  * Drozdov  _et al._ [2015] A. Drozdov, M. Eremets, and I. Troyan, Superconductivity above 100 k in ph3 at high pressures, arXiv preprint arXiv:1508.06224 (2015).
  * Shamp  _et al._ [2016] A. Shamp, T. Terpstra, T. Bi, Z. Falls, P. Avery, and E. Zurek, Decomposition products of phosphine under pressure: Ph2 stable and superconducting?, J. Am. Chem. Soc. 138, 1884 (2016).
  * Flores-Livas  _et al._ [2016] J. A. Flores-Livas, M. Amsler, C. Heil, A. Sanna, L. Boeri, G. Profeta, C. Wolverton, S. Goedecker, and E. K. U. Gross, Superconductivity in metastable phases of phosphorus-hydride compounds under high pressure, Phys. Rev. B 93, 020508 (2016).
  * Bi  _et al._ [2017] T. Bi, D. P. Miller, A. Shamp, and E. Zurek, Superconducting phases of phosphorus hydride under pressure: Stabilization via mobile molecular hydrogen, Angew. Chem. Int. Ed. 56, 10192 (2017).
  * Snider  _et al._ [2020] E. Snider, N. Dasenbrock-Gammon, R. McBride, M. Debessai, H. Vindana, K. Vencatasamy, K. V. Lawler, A. Salamat, and R. P. Dias, Room-temperature superconductivity in a carbonaceous sulfur hydride, Nature 586, 373 (2020).
  * Cui  _et al._ [2020] W. Cui, T. Bi, J. Shi, Y. Li, H. Liu, E. Zurek, and R. J. Hemley, Route to high-_Tc_ superconductivity via ch4 intercalated h3s hydride perovskites, Phys. Rev. B 101, 134504 (2020).
  * Sun  _et al._ [2020] Y. Sun, Y. Tian, B. Jiang, X. Li, H. Li, T. Iitaka, X. Zhong, and Y. Xie, Computational discovery of a dynamically stable cubic sh3-like high-temperature superconductor at 100 gpa via ch4 intercalation, Phys. Rev. B 101, 174102 (2020).
  * Shatruk [2019] M. Shatruk, Thcr2si2 structure type: The “perovskite” of intermetallics, J. Solid State Chem. 272, 198 (2019).
  * Dong  _et al._ [2018] X. Dong, Q. Wu, O. V. Yazyev, X.-L. He, Y. Tian, X.-F. Zhou, and H.-T. Wang, Condensed matter realization of fermion quasiparticles in minkowski space (2018), arXiv:1804.00570 [cond-mat.mtrl-sci] .
  * Longuet-Higgins and d. V. Roberts [1955] H. C. Longuet-Higgins and M. d. V. Roberts, The electronic structure of an icosahedron of boron atoms, Proc. Royal Soc. London A230, 110 (1955).
  * Jemmis and M. [2001] E. D. Jemmis and B. M. M., Polyhedral boranes and elemental boron: Direct structural relations and diverse electronic requirements, J. Am. Chem. Soc. 123, 4324 (2001).
  * Fujimori  _et al._ [1999] M. Fujimori, T. Nakata, T. Nakayama, E. Nishibori, K. Kimura, M. Takata, and M. Sakata, Peculiar covalent bonds in $\alpha$-rhombohedral boron, Phys. Rev. Lett. 82, 4452 (1999).
  * Wade [1976] K. Wade, Structural and bonding patterns in cluster chemistry, Adv. Inorg. Chem. Radiochem. 18, 1 (1976).
  * Mingos [1984] D. M. P. Mingos, Polyhedral skeletal electron pair approach, Acc. Chem. Res. 7, 311 (1984).
  * Teter [1998] D. M. Teter, Computational alchemy: The search for new superhard materials, MRS Bull. 23, 22 (1998).
  * Chen  _et al._ [2011] X.-Q. Chen, H. Niu, D. Li, and Y. Li, Modeling hardness of polycrystalline materials and bulk metallic glasses, Intermetallics 19, 1275 (2011).
  * Isayev  _et al._ [2017] O. Isayev, C. Oses, C. Toher, E. Gossett, S. Curtarolo, and A. Tropsha, Universal fragment descriptors for predicting properties of inorganic crystals, Nat. Commun. 8, 15679 (2017).
  * Curtarolo  _et al._ [2012] S. Curtarolo, W. Setyawan, S. Wang, J. Xue, K. Yang, R. H. Taylor, L. J. Nelson, G. L. Hart, S. Sanvito, M. Buongiorno-Nardelli, N. Mingo, and O. Levy, AFLOWLIB.ORG: A distributed materials properties repository from high-throughput ab initio calculations, Comput. Mater. Sci. 58, 227 (2012).
  * Toher  _et al._ [2018] C. Toher, C. Oses, D. Hicks, E. Gossett, F. Rose, P. Nath, D. Usanmaz, D. C. Ford, E. Perim, C. E. Calderon, J. J. Plata, Y. Lederer, M. Jahnatek, W. Setyawan, S. Wang, J. Xue, K. Rasch, R. V. Chepulskii, R. H. Taylor, G. Gomez, H. Shi, A. R. Supka, R. A. R. Al Orabi, P. Gopal, F. T. Cerasoli, L. Liyanage, H. Wang, I. Siloi, L. A. Agapito, C. Syshadham, G. L. Hart, J. Carrete, F. Legrain, N. Mingo, E. Zurek, O. Isayev, A. Tropsha, S. Sanvito, R. M. Hanson, I. Takeuchi, M. J. Mehl, A. N. Kolmogorov, K. Yang, P. D’Amico, A. Calzolari, M. Costa, R. De Gennaro, M. B. Nardelli, M. Fornari, O. Levy, and S. Curtarolo, The aflow fleet for materials discovery, in  _Handbook of Materials Modeling. Volume 1 Methods: Theory and Modeling_ , Vol. 1, edited by W. Andreoni and S. Yip (Springer, Cham, Switzerland, 2018) pp. 1–28.
  * Avery  _et al._ [2019] P. Avery, X. Wang, C. Oses, E. Gossett, D. Proserpio, C. Toher, S. Curtarolo, and E. Zurek, Predicting superhard materials via a machine learning informed evolutionary structure search, Npj Comput. Mater. 5, 89 (1 (2019).
  * Gao and Gao [2010] F. M. Gao and L. H. Gao, Microscopic models of hardness, J. Superhard Mater. 32, 148 (2010).
  * [108] Http://hdl.handle.net/10477/79221, July 05, 2020\.


11institutetext:  † Centre de Physique Théorique de Luminy, Université de la
Méditerranée, F-13288 Marseille FR, and “Sapienza” Università di Roma,
Dipartimento di Fisica and ICRA,
P.le A. Moro 5, 00185 Rome IT<EMAIL_ADDRESS>
‡ Institut des Hautes Etudes Scientifiques, 91440 Bures-sur-Yvette FR, and
“Sapienza” Università di Roma, Dipartimento di Fisica and ICRA, P.le A. Moro
5, 00185 Rome IT<EMAIL_ADDRESS>
§ ICRA, ICRANet, ENEA and “Sapienza” Università di Roma, Dipartimento di
Fisica, P.le A. Moro 5, 00185 Rome IT<EMAIL_ADDRESS>
# GUP vs polymer quantum cosmology: the Taub model

Marco Valerio Battisti† Orchidea Maria Lecian‡ Giovanni Montani§

###### Abstract

The fate of the cosmological singularity in the Taub model is discussed within
the two frameworks. An internal time variable is ruled out and the only
remaining degree of freedom (the anisotropy) of the Universe is quantized
according to such schemes. The resulting GUP Taub Universe is singularity-
free, differently from the second case, where the classical singularity is not
tamed by the polymer-loop quantum effects.

##  1 Introduction

Two different quantum cosmology approaches are applied to the Taub model. The
study is performed at the classical and at the quantum level in both schemes.
In particular, the generalized uncertainty principle (GUP) and the polymer
(loop) frameworks are implemented to this system. In the first case [1], the
cosmological singularity appears to be probabilistically suppressed, while, in
the second one [2], the Universe is still singular. Such a feature then allows
us to better understand the avoidance of the cosmological singularity in other
different quantum gravity toy models.

The Taub Universe arises as a particular case of the Bianchi IX model, i.e.
the most general scheme allowed by the homogeneity constraint (for reviews see
[3]). It is obtained by restricting the dynamics to that of a one-dimensional
particle bouncing against a wall, when only one degree of freedom (the
Universe anisotropy) is taken into account. The relevance of the Taub Universe
in quantum cosmology is then due to the fact that it is a necessary step
towards the Bianchi IX model, being a generalization of other isotropic
models. In particular, it has been used to test the validity of the
minisuperspace scheme [4] and to explore the application of the extrinsic
cosmological time [5].

The paper is organized as follows. In Section 2 the Taub model and the
formalisms are reviewed. Section 3 and 4 are devoted to the classical and
quantum analysis respectively. Comparisons with other approaches follow in
Section 5. We adopt $\hbar=c=16\pi G=1$ units.

##  2 Taub Universe and formalisms

The Taub cosmological model is a particular case of the Bianchi IX Universe,
of which only one anisotropy is taken into account. The Bianchi IX model,
together with Bianchi VIII, is the most general spatially-homogeneous model
and is described by the line element
$ds^{2}=N^{2}dt^{2}-e^{2\alpha}\left(e^{2\gamma}\right)_{ij}\omega^{i}\otimes\omega^{j}$
[3]. Here $N=N(t)$ is the lapse function, $\omega^{i}=\omega^{i}_{a}dx^{a}$
are the $SO(3)$ left-invariant 1-forms, $\alpha=\alpha(t)$ describes the
isotropic expansion of the Universe and $\gamma_{ij}=\gamma_{ij}(t)$ is a
traceless symmetric matrix, which determines the anisotropies via
$\gamma_{\pm}$. The dynamics of this Universe towards the singularity is
described by the motion of a two-dimensional particle (the two physical degree
of freedom of the gravitational field) in a dynamically-closed domain [3]. In
the Misner picture, such a domain depends on the time variable $\alpha$,
while, in the Misner-Chitré one, it becomes stationary in time. Performing the
ADM reduction of the dynamics (according to which the classical constraints
are solved with respect to the given momenta before implementing any
quantization algorithm), an effective Hamiltonian is obtained, which depends
only on the physical degrees of freedom of the system. In particular, the
scalar constraint is solved with respect to the momentum conjugated to the
time variable $p_{\tau}$ (we adopt the time gauge $\dot{\tau}=1$) and,
performing another change of variables, we obtain

$-p_{\tau}\equiv H^{IX}_{ADM}=v\sqrt{p_{u}^{2}+p_{v}^{2}}.$ (1)

The dynamics of such a system is equivalent to a billiard ball on a
Lobatchevsky plane and the three corners of the Misner scheme are replaced by
the points $(0,0)$, $(-1,0)$ and $v\rightarrow\infty$ in the $(u,v)$-plane, as
in Fig. 1.

Figure 1: The dynamical-allowed domain in the $(u,v)$-plane where the dynamics
is restricted.

The Taub cosmological model is described by $\gamma_{-}=0$. The dynamics of
this Universe is equivalent to the motion of a particle in a one-dimensional
closed domain. Such a domain corresponds the choice of only one of the three
equivalent potential walls of the Bianchi IX model. The ADM Hamiltonian (1)
rewrites $H_{ADM}^{T}=vp_{v}$, where $v\in[1/2,\infty)$. This Hamiltonian can
be further simplified defining the new variable $x=\ln v$, and becomes

$H_{ADM}^{T}=p_{x}\equiv p$ (2)

where the configuration variable $x$ is related to the Universe anisotropy
$\gamma_{+}$ by the equation $\gamma_{+}=e^{\tau-x}(e^{2x}-3/4)/\sqrt{3}$. The
Hamiltonian (2) will be the starting point of our analysis. It is worth noting
that the classical singularity now appears for $\tau\rightarrow\infty$.

###  2.1 GUP quantum mechanics

Some issues and results of a non-relativistic quantum mechanics with non-zero
minimal uncertainties in position are briefly reviewed [6]. In one dimension,
we consider the Heisenberg algebra generated by $q$ and $p$ obeying the
commutation relation

$[q,p]=i(1+\beta p^{2}),$ (3)

where $\beta>0$ is a deformation parameter. This commutation relation leads to
the generalized uncertainty relation $\Delta q\Delta
p\geq\frac{1}{2}\left(1+\beta(\Delta p)^{2}+\beta\langle p\rangle^{2}\right)$,
which appears in string theory [7]. The canonical Heisenberg algebra can be
recovered in the limit $\beta=0$, and the generalization to more dimension is
straightforward, leading naturally to a “noncommutative geometry” for the
space coordinates.

It is immediate to verify that a finite minimal uncertainty in position
$\Delta q_{min}=\sqrt{\beta}$ is predicted. The existence of a non-zero
uncertainty in position is relevant since it implies that there cannot be any
physical state that is a position eigenstate. In fact, an eigenstate of an
observable necessarily has a vanishing uncertainty on it. Although it is
possible to construct position eigenvectors, they are only formal eigenvectors
and not physical states. The deformed Heisenberg algebra (3) can be
represented in the momentum space as

$p\psi(p)=p\psi(p),\qquad q\psi(p)=i(1+\beta p^{2})\partial_{p}\psi(p),$ (4)

on a dense domain $S$ of smooth functions. To recover information on position
we have to study those states, which realize the maximally-allowed
localization. Such states $|\psi^{ml}_{\zeta}\rangle$, which are proper
physical states around a position $\zeta$, have the proprieties
$\langle\psi^{ml}_{\zeta}|q|\psi^{ml}_{\zeta}\rangle=\zeta$ and $(\Delta
q)_{|\psi^{ml}_{\zeta}\rangle}=\Delta q_{min}$. We can project an arbitrary
state $|\psi\rangle$ on the maximally-localized states
$|\psi^{ml}_{\zeta}\rangle$ to obtain the probability amplitude for a particle
being maximally localized around the position $\zeta$ (i.e. with standard
deviation $\Delta q_{min}$). We call these projections the quasiposition wave
function $\psi(\zeta)=\langle\psi^{ml}_{\zeta}|\psi\rangle$, and explicitly we
have the generalized Fouries transformation

$\psi(\zeta)\sim\int^{+\infty}_{-\infty}\frac{dp}{(1+\beta
p^{2})^{3/2}}\exp\left(i\frac{\zeta}{\sqrt{\beta}}\tan^{-1}(\sqrt{\beta}p)\right)\psi(p).$
(5)

As $\beta\rightarrow 0$, the ordinary position wave function
$\psi(\zeta)=\langle\zeta|\psi\rangle$ is recovered.

###  2.2 Polymer quantum mechanics

The polymer representation of quantum mechanics consists in defining abstract
kets, labeled by a real number and assumed to form an orthonormal basis, and
then considering a suitable finite subset of them, whose Hilbert space is
defined by the corresponding inner product [8]. This procedure helps one gain
insight onto some particular features of quantum mechanics, when an underlying
discrete structure is somehow hypothesized. The request that the Hamiltonian
associated to the system be of direct physical interpretation defines the
polymer phase space, and the continuum limit can be recovered by the
introduction of the concept of the scale [9].

In the particular case of a discrete position variable in the momentum
polarization, the Hamiltonian variable $p$ cannot be implemented as an
operator, so that some restrictions on the model have to be required. If the
set of kets is restricted by the introduction of a regular graph
$\gamma_{\mu_{0}}$, the kynetic term of the Hamiltonian is approximated by the
polymer substitution

$p\rightarrow\frac{1}{\mu_{0}}\sin(\mu_{0}p),$ (6)

where the incremental ratio is evluated for an exponentiated operator. The
Hamiltonian operator $H_{\mu_{0}}$, which lives in
$\mathcal{H}_{\gamma_{\mu_{0}}}$, reads

$H_{\mu_{0}}=\frac{\hat{p}_{\mu_{0}}^{2}}{2m}+V(\hat{q}).$ (7)

The definition of a scale, $C_{n}$, eables one to approximate continuous
functions with functions that are constant on the intervals. As a result, at
any given scale $C_{n}$, the kinetic term of the Hamiltonian operator can be
approximated, and effective theories at given scales are related by coarse-
graining maps.

##  3 Deformed classical dynamics

The ordinary Taub model can be interpreted as a massless scalar relativisitic
particle moving in the Lorentzian minisuperspace ($\tau,x$)-plane, whose the
classical trajectory is its light-cone. More precisely, the incoming particle
($\tau<0$) bounces on the wall ($x=x_{0}=\ln(1/2)$) and falls into the
classical cosmological singularity ($\tau\rightarrow\infty$). Investigating
the modification of the dynamics within the two frameworks will show that the
two behaviors can be interpreted as complementary.

###  3.1 GUP framework

The GUP-classical dynamics is contained in the modified symplectic geometry
arising from the classical limit of (3), as soon as the parameter $\beta$ is
regarded as an independent constant with respect $\hbar$. It is then possible
to replace the quantum-mechanical commutator (3) via its Poisson brackets,
i.e. $-i[q,p]\Longrightarrow\\{q,p\\}=(1+\beta p^{2})$. The Poisson brackets
for any two-dimensional phase space function are

$\\{F,G\\}=\left(\frac{\partial F}{\partial q}\frac{\partial G}{\partial
p}-\frac{\partial F}{\partial p}\frac{\partial G}{\partial q}\right)(1+\beta
p^{2}).$ (8)

Applying this scheme to the Hamiltonian (2), we immediately obtain the
equations of motion for the model [1],

$x(\tau)=(1+\beta A^{2})\tau+cost,\qquad p(\tau)=cost=A,$ (9)

where $x\in[x_{0},\infty)$. Therefore, at the classical level, the effects of
the deformed Heisenberg algebra (3) on the Taub Universe are as follows. The
angular coefficient is $(1+\beta A^{2})>1$ for $\beta\neq 0$, and thus the
angle between the two straight lines $x(\tau)$, for $\tau<0$ and $\tau>0$,
becomes smaller as the values of $\beta$ grows. The trajectories of the
particle (Universe), before and after the bounce on the potential wall at
$x=x_{0}\equiv\ln(1/2)$, are closer to each other then in the canonical case
($\beta=0$).

###  3.2 Polymer framework

The polymer-classical dynamics relies on the substitution (6) in the
Hamiltonian of the model (2). This way, the equations of motion rewrite [2]

$\dot{x}=\left\\{x,H\right\\}=\cos(ap),\qquad\dot{p}=\left\\{p,H\right\\}=0,$
(10)

where a dot denotes differentiation with respect to the time variable $\tau$.
The equations of motion are immediately solved as

$x(\tau)=\cos(ap)\tau,\qquad p(\tau)=A,$ (11)

where $A$ is a constant. In the discretized (polymer) case, i.e. for $a\neq
0$, the one-parameter family of trajectories flattens. In fact, the angle
between the incoming trajectory and the outgoing one is greater than $\pi/2$,
since $p\in\left(-\pi/a,\pi/a\right)$. As these trajectories diverge rather
than converge, we expect the polymer quantum effects to be reduced with
respect to the classical case, as we will verify below.

##  4 Deformed quantum dynamics

The quantum dynamics of the Taub Universe is here investigated according to
the two different approaches. Particular attention is paid to the wave-packet
evolution and the consequential fate of the classical cosmological
singularity. In both frameworks, the variable $\tau$ is regarded as a time
coordinate and therefore ($\tau,p_{\tau}$) are treated in the canonical way.
The deformed quantization (GUP or polymer) is then implemented only to the
submanifold describing the only degree of freedom of the Universe, i.e. the
phase space spanned by ($x,p$). We then deal with a Schrödinger-like equation

$i\partial_{\tau}\Psi(\tau,p)=\hat{H}_{ADM}^{T}\Psi(\tau,p),$ (12)

where the operator $\hat{H}_{ADM}^{T}$ accounts for the modifications due to
the two frameworks. We have to square the eigenvalue problem in order to
correctly impose the boundary condition: in agreement with the analysis
developed in [10], we make the well-grounded hypothesis that the
eigenfunctions form be independent of the presence of the square root, since
its removal implies the square of the eigenvalues only. The wave packets,
which are superposition of the eigenfunctions
$\Psi(\tau,x)=\int_{0}^{\infty}dkA(k)\psi_{k}(x)e^{-ik\tau}$, are then
constructed for both models, taking $A(k)$ as a Gaussian-like weighting
function. The differences between the two approaches are due to the features
of the eigenfunctions $\psi_{k}(x)$. Analyzing such an evolution, we show that
the GUP Taub Universe appears to be probabilistically singularity-free,
differently from the polymer case, where the singularity is not tamed by the
cut-off-scale effects.

###  4.1 GUP framework

We now analyze the model in the GUP approach [1]. As explained before, we lost
all informations on the position itself, so that the boundary conditions have
to be imposed on the quasiposition wave function (5), i.e. $\psi(\zeta_{0})=0$
(where $\zeta_{0}=\langle\psi^{ml}_{\zeta}|x_{0}|\psi^{ml}_{\zeta}\rangle$, in
agreement with the previous discussion). The solution of the eigenvalue
problem is the Dirac $\delta$-distribution $\psi_{k}(p)=\delta(p^{2}-k^{2})$,
and therefore the quasiposition wave function (5) reads

$\psi_{k}(\zeta)=\frac{A}{k(1+\beta
k^{2})^{3/2}}\left[\exp\left(i\frac{\zeta}{\sqrt{\beta}}\tan^{-1}(\sqrt{\beta}k)\right)-\exp\left(i\frac{(2\zeta_{0}-\zeta)}{\sqrt{\beta}}\tan^{-1}(\sqrt{\beta}k)\right)\right],$
(13)

where $A$ is a constant and the boundary condition $\psi(\zeta_{0})=0$ has
been imposed. The deformation parameter $\beta$, i.e. the presence of a non-
zero minimal uncertainty for the configuration variable, is responsible for
the GUP effects on the dynamics. The physical interpretation of $\beta$ is
then a non-zero minimal uncertainty in the anisotropy of the Universe. To
better understand the modifications induced by the deformed Heisenberg algebra
on the canonical Universe dynamics, we have to analyze different
$\beta$-regions. In fact, when $\beta$ becomes more and more important, i.e.
when we are at some scale that allows us to appreciate the GUP effects, the
evolution of the wave packets is different from the canonical case. More
precisely, these effects are present when the product $k_{0}\sqrt{\beta}$
becomes remarkable, i.e. when $k_{0}\sqrt{\beta}\sim\mathcal{O}(1)$, and
therefore when $\beta$ is comparable to $1/k_{0}^{2}$. In fact, the correct
semiclassical behaviors of the model far away from the singularity is
described by wave packets peaked at energies much smaller then
$1/\sqrt{\beta}$ [1]. In particular, for $k_{0}=1$, we can distinguish between
three different $\beta$-regimes:

  * •

$\beta\sim\mathcal{O}(10^{-2})$ regime. The wave packets begin to spread and a
constructive and destructive interference between the incoming and outgoing
wave appears. The probability amplitude to find the Universe is still peaked
around the classical trajectory.

  * •

$\beta\sim\mathcal{O}(10^{-1})$ regime. It is no more possible to distinguish
an incoming or outgoing wave packet and, at this level, the notion of a wave
packet following a classical trajectory becomes meaningless.

  * •

$\beta\sim\mathcal{O}(1)$ regime. A dominant probability peak “near” the
potential wall appears. There are also other small peaks for growing values of
$\zeta$, but they are widely suppressed for bigger $\beta$. The motion of wave
packets shows a stationary behavior, i.e. these are independent of $\tau$. See
Fig. 2.

Following this picture we are able to learn the GUP modifications to the WDW
wave packets evolution. In fact, from small to big values of $\beta$, we can
see how the wave packets escape from the classical trajectories and approach a
stationary state close to the potential wall.

Figure 2: Wave packets $|\Psi(\tau,\zeta)|$ in the GUP framework as $\beta
k_{0}^{2}=1$ ($k_{0}=1$ and $\sigma=4$).

Such a behavior is, in some sense, expected from a classical point of view. In
fact, at classical level the ingoing and the outgoing trajectories shrink each
other. So a quantum probability interference is a fortiori predicted. On the
other hand, the stationarity feature exhibited by the Universe in the
($\beta\sim\mathcal{O}(1)$)-region is a purely quantum GUP effect. Such a
behavior cannot be inferred from a deformed classical analysis. From this
point of view, the classical singularity ($\tau\rightarrow\infty$) is widely
probabilistically suppressed, because the probability to find the Universe is
peaked just around the potential wall. This way we claim that the GUP-Taub
Universe is singularity-free.

###  4.2 Polymer framework

We now analyze the model in the polymer approach [2]. For the quantum analysis
of the model, we choose a discretized $x$ space, and solve the corresponding
eigenvalue problem in the $p$ polarization. Considering the time evolution for
the wave function $\Psi$ as given by $\Psi_{k}(p,\tau)=e^{-ik\tau}\psi_{k}(p)$
and the results of [10], we obtain the following eigenvalue problem

$(p^{2}-k^{2})\psi_{k}(p)=\left[\frac{2}{a^{2}}\left(1-\cos(ap)\right)-k^{2}\right]\psi_{k}(p),$
(14)

solved by

$\displaystyle k^{2}=k^{2}(a)=\frac{2}{a^{2}}\left(1-\cos(ap)\right)\leq
k^{2}_{max}=\frac{4}{a^{2}}$ (15a)
$\displaystyle\psi_{k,a}(p)=A\delta(p-p_{k,a})+B\delta(p+p_{k,a})$ (15b)
$\displaystyle\psi_{k,a}(x)=A\left[\exp(ip_{k,a}x)-\exp(ip_{k,a}(2x_{0}-x))\right]:$
(15c)

(15b) is the momentum wave function, with $A$ and $B$ two arbitrary
integration constant, and (15c) is the coordinate wave function, where an
integration constant has been eliminated by imposing suitable boundary
conditions. Moreover, we have defined the modified dispersion relation

$p_{k,a}\equiv\frac{1}{a}\arccos\left(1-\frac{k^{2}a^{2}}{2}\right)$ (16)

from (15a). Furthermore, we stress that $k^{2}$ is bounded from above, as
illustrated in (15a), but it is its square root, considered for its positive
determination, which accounts for the time evolution of the wave function.

We now construct suitable wave packets $\Psi(x,\tau)$ taking into account the
previous discussion (note that a maximum energy $k_{max}$ is now predicted).
Three relevant cases can be distinguished:

Figure 3: The spread polymer wave packet $|\Psi(x,\tau)|$ as $k_{0}a=1/2$
($a=50$, $k_{0}=0.01$, $\sigma=0.125$).

  * •

$k_{0}a\sim\mathcal{O}(1)$ and peaked weighting function. The resulting wave
packet is well approximated by a purely monochromatic wave. A small
interference phenomenon between the wave and the wall is then predicted.

  * •

$k_{0}a\sim\mathcal{O}(1)$ and spread weighting function. A strong
interference phenomenon between the wave and the wall now appears.
Nevertheless, this interference phenomenon is not able to probabilistically
tame the singularity, as it takes place in the ’outer’ region, in a way
complementary to that of the GUP approach (see Fig. 3). The polymer-Taub
Universe is then still a singular cosmological model.

  * •

$k_{0}a\ll\mathcal{O}(1)$ regime. This can be considered as the semi-classical
limit of the model. In fact, differently from the other cases, the value of
$k_{0}$ around which the wave packet is peaked is not arbitrary, but
constrained by the characteristic scale $a$ under investigation. The ordinary
WDW behavior is therefore recasted.

##  5 Comparison with other approaches

The Taub cosmological model offers a suitable scenario, where different
quantization techniques can be applied. In fact, it is possible to single out
a time variable, so that the anisotropy describes the real degree of freedom
of the Universe. It is therefore reasonable to investigate the fate of the
cosmological singularity without modifying the time variable. The comparison
with analysis of the cosmological singularity in other cosmological models
outlines how the features of the Taub model allow one to pick the cut off
effects out of those due to the choice of the Hamiltonian variables.

In the cosmological isotropic sector of GR, i.e. the FRW models, the
singularity is removed by loop quantum effects. The wave function of the
universe exhibits a non-singular behavior at the classical singularity, and
the big-bang is replaced by a big-bounce, when a free scalar field is taken as
the relational time [11]. The scale factor of the universe is directly
quantized by the use of the polymer (loop) techniques, so that the evolution
itself of the wave packet of the universe is deeply modified by such an
approach. The Hamiltonian constraint does not allow for a constant solution of
the variable conjugated to the scale factor, so that it is not possible to
choose a scale, such that the polymer modifications are negligible throughout
the whole evolution, so that the comparison with the ordinary representation
is not always possible. Anyhow, we stress that, for the Taub model, the
cosmological singularity is probabilistically suppressed, regardless to the
fact whether the system can appreciated or not the cut off during the whole
evolution.

In [12], all the degrees of freedom of the Bianchi cosmological models in the
ADM reduction of the dynamics are quantized by loop techniques. In particular,
also the time variable, i.e. the Universe volume, is treated at the same level
as the others. In most cases, the time variable is defined by a phase space
variable, i.e. it is an internal one. As a result, also the Bianchi Universes
are singularity-free [13]. In this respect, our analysis is based on
considering the time variable as an ordinary Heisenberg variable, while the
cut off is imposed on the anisotropy only.

The GUP dynamics of other cosmological models has been investigated in
different approaches. In particular, the big-bang singularity appears to be
tamed by GUP effects showing a stationary behavior of the wave packets [14].
Such a prediction is in agreement with those achieved in a noncommutative
quantum cosmology [15]. However, in order to predict a big bounce à la LQC, a
Snyder-deformed quantum cosmology has to be addressed [16]. As the last point,
it is interesting to notice that the GUP-Mixmaster Universe is still a chaotic
model [17], as opposite to the LQC one [18], the difference being essentially
based on the application of the deformed scheme to the time variable too.

Acknowledgments. The “Angelo Della Riccia” Fellowship and the “Sapienza CUN 2”
Fellowship are gratefully acknowledged.

## References

  * [1] M.V.Battisti and G.Montani, Phys.Rev.D 77 (2008) 023518.
  * [2] M.V.Battisti, O.M.Lecian and G.Montani, Phys.Rev.D 78 (2008) 103514.
  * [3] G.Montani, M.V.Battisti, R.Benini and G.Imponente, Int.J.Mod.Phys.A 23 (2008) 2353; J.M.Heinzle and C.Uggla, arXiv:0901.0776.
  * [4] K.Kuchar and M.P.Ryan, Phys.Rev.D 40 (1989) 3982.
  * [5] G.Catren and R.Ferraro, Phys.Rev.D 63 (2001) 023502.
  * [6] A.Kempf, G.Mangano and R.B.Mann, Phys.Rev.D 52 (1995) 1108; A.Kempf, J.Math.Phys. 38 (1997) 1347.
  * [7] D.J.Gross and P.F.Mendle, Nucl.Phys.B 303 (1988) 407; D.Amati, M.Ciafaloni and G.Veneziano, Phys.Lett.B 216 (1989) 41.
  * [8] A.Ashtekar, S.Fairhurst and J.L.Willis, Class.Quant.Grav. 20 (2003) 1031; A.Corichi, T.Vukasinac and J.A.Zapata, Phys.Rev.D 76 (2007) 0440163.
  * [9] A.Corichi, T.Vukasinac and J.A.Zapata, Class.Quant.Grav. 24 (2007) 1495.
  * [10] R.Puzio, Class.Quant.Grav. 11 (1994) 609.
  * [11] A.Ashtekar, T.Pawlowski and P.Singh, Phys.Rev.Lett. 96 (2006) 141301; Phys.Rev.D 73 (2006) 124038.
  * [12] M.Bojowald, G.Date and K.Vandersloot, Class.Quant.Grav. 21 (2004) 1253; M.Bojowald, Class.Quant.Grav. 20 (2003) 2595.
  * [13] D.W.Chiou, Phys.Rev.D 76 (2007) 124037; G.Date, Phys.Rev.D 71 (2005) 127502.
  * [14] M.V.Battisti and G.Montani, Phys.Lett.B 656 (2007) 96.
  * [15] B.Vakili and H.R.Sepangi, Phys.Lett.B 651 (2007) 79; H.R.Sepangi, B.Shakerin and B.Vakili, Class.Quant.Grav. 26 (2009) 065003.
  * [16] M.V.Battisti, arXiv:0805.1178; J.Phys.Conf.Ser. (at press), arXiv:0810.5039.
  * [17] M.V.Battisti and G.Montani, arXiv:0808.0831.
  * [18] M.Bojowald and G.Date, Phys.Rev.Lett. 92 (2004) 071302.