# Hydrogen phase-IV characterization by full account of quantum anharmonicity

Tommaso Morresi,1∗ Rodolphe Vuilleumier,2 Michele Casula1∗

1Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC),
Sorbonne Université, CNRS UMR 7590, MNHN, 4 Place Jussieu, 75252 Paris, France
2PASTEUR, Département de chimie, École normale supérieure, PSL University,
Sorbonne Université, CNRS, 75005 Paris, France

###### Abstract

We devise a framework to compute accurate phonons in molecular crystals even
in case of strong quantum anharmonicity. Our approach is based on the
calculation of the static limit of the phononic Matsubara Green’s function
from path integral molecular dynamics simulations. Our method enjoys a
remarkably low variance, which allows one to compute accurate phonon
frequencies after a few picoseconds of nuclear dynamics, and it is further
stabilized by the use of appropriate constrained displacement operators. We
applied it to solid hydrogen at high pressure. For phase III, our predicted
infrared (IR) and Raman active vibrons agree very well with experiments. We
then characterize the crystalline symmetry of phase IV by direct comparison
with vibrational data and we determine the character and isotopic shift of its
Raman and IR vibron peaks.

##  I Introduction

One of the main structural features of hydrogen-rich materials is the presence
of large quantum fluctuations affecting hydrogen nuclear motion, and leading
to strong anharmonicity. Pristine hydrogen is intriguing because of its
unexpectedly rich phase diagram in the high-pressure range [1, 2, 3, 4, 5, 6,
7]. However, a precise experimental assessment of all its phases is difficult,
since hydrogen is a weak X-ray and neutron scatterer. In fact, the most direct
structural information comes from Raman and infrared (IR) techniques, probing
its vibrational properties. Hence, the underlying geometries could be deduced
by comparing experimental data with theoretical predictions. However, that
comparison is hindered by nuclear quantum fluctuations, which make harmonic
and perturbative methods unreliable.

In this work, we generalize and extend to molecular solids a framework
previously introduced [8] to compute vibrational properties of strongly
anharmonic systems. This goal is achieved by accessing the static limit of the
phononic Matsubara Green’s function and by imposing a constraint on the
quantum displacement operators defining the Green’s function. As application,
we determine the vibron modes of phase III and IV in solid hydrogen, entirely
from first principles and with unprecedented accuracy. We then address the
open issue of identifying the crystalline symmetry of phase IV.

Previous calculations of the vibrational properties of hydrogen were based on
classical molecular dynamics (MD) [9, 10], which cannot include quantum
anharmonicity. To overcome the MD limitations and include nuclear quantum
effects (NQE), more advanced methods have been developed, such as the
vibrational self-consistent field (VSCF) method [11, 12, 13] and the
stochastic self-consistent harmonic approximation (SSCHA) [14, 15, 16].
Although they account for NQE, these approximations heavily rely on the
quality of their variational ansatz wavefunction. In our framework instead, we
fully include NQE by sampling the exact quantum thermal distribution, for a
given electronic theory, through path integral MD (PIMD) simulations [17, 18]
driven by ab initio forces. PIMD gives direct access to imaginary-time
correlation functions and in particular to the exact phononic Matsubara
Green’s function [19], defined as

$G_{ij}(\tau)=-\sqrt{m_{i}m_{j}}\langle\mathcal{T}\delta\hat{x}_{i}(\tau)~{}\delta\hat{x}_{j}(0)\rangle,$
(1)

with $\tau$ the imaginary time, $\mathcal{T}$ the time-ordered operator,
$\delta x_{i}$ the displacement of $i$-th Cartesian degree of freedom with
respect to its equilibrium position, $m_{i}$ its mass, and the brackets
indicate the average over the PIMD distribution. However, while Eq. (1)
encloses all information on the vibrational properties of the system, the
analytic continuation to extract the phonon spectral function suffers from two
main issues: it is numerically ill-defined [20] and is plagued by sampling
errors inherent in the finite size of PI trajectories. So far, these
drawbacks, together with a usually large computational cost, have severely
limited practical PIMD-based phonon calculations in realistic Hamiltonians.

In our framework the above issues are solved by taking the Kubo-transform
$\tilde{\mathbf{G}}$ of the Green’s function in Eq. (1), obtained by averaging
the displacement operators over the whole ring-polymer chain representing a
quantum particle in the PIMD framework [21]. This corresponds precisely to the
static limit of the phononic Matsubara Green’s function. By calculating the
inverse $\tilde{\mathbf{G}}^{-1}$, we show that one can get an accurate
estimate of the NQE-renormalized phonon frequencies, without requiring any
analytic continuation. We also provide a computational scheme where
$\tilde{\mathbf{G}}^{-1}$ can be evaluated with low variance and low bias.
This goal is achieved by solving a generalized eigenvalue problem (GEV), where
the sampling errors of the displacement autocorrelation function are
compensated by those of the velocity autocorrelation.

Nevertheless, the calculation of vibrational properties of molecular solids is
hampered by the simultaneous presence of intramolecular and lattice degrees of
freedom.

Figure 1: 2D dimer model. (a) Convergence of the PIMD eigenvalue for the
vibron frequency ($\omega_{\delta\bar{x}\delta\bar{x},\ II}$) of the dimer
model in Eq. (4) with $c=$ 0.05 as a function of the PIMD evolution. The GEV
frequencies (dark-blue points) converge after a few picoseconds to the
solution obtained by solving numerically the corresponding Schrödinger
equation. (b) Potential shape for three different values of the parameter $c$
in Eq. (4): left column c=0.005; central column c=0.05; right column c=0.5.
(c) GEV frequencies (Eq. (2)) with exact unconstrained
$\delta\hat{\mathbf{x}}$ (top panel) and constrained displacements (lower
panel) as a function of the coupling strength $c$ (Eq. (4)). (d) Square
modulus of the ground state wave function ($|\Psi_{0}|^{2}$) for the three
values of $c$ in panel (b). (e) Imaginary part of the retarded Green’s
functions for the dimer model of Eq. (4) obtained by the numerically exact
solution (light green) and with the constrained (blue) displacement operators
for c=0.005, c=0.05 and c=0.5. The effect of constraints in Eq. (II.2) is to
suppress some of the rotational poles of the unconstrained Green’s function.
The small imaginary part in the denominator of the retarded Green’s function,
corresponding to the width of the peaks, is set equal to 2 cm-1.

Indeed, while it has been shown [8] that this approach accurately reproduces
the exact spectrum for various model Hamiltonians even in the strong
anharmonic limit, its reliability deteriorates for systems made of rotating
molecules, such as the molecular phases of solid hydrogen [22]. These systems
are characterized by the presence of intramolecular libration and vibron
modes. In case of large librations, it turns out that the vibron mode is
systematically underestimated, as the $\tilde{\mathbf{G}}^{-1}$ eigenvalues
are biased by the hybridization between these modes. In this paper, we
demonstrate how to fix this problem by introducing constraints on the
displacement operators in Eq. (1), that effectively reduce the phonon Fock
space to the most relevant sectors. After benchmarking the performances of our
constrained PIMD (cPIMD) framework (Sec. II.1) against the exact solution for
a rotor model, which represents a dimer molecule interacting with an external
field (Sec. II.2), we prove the generality of cPIMD by applying it to phases
III and IV of solid hydrogen (Secs. II.3 and II.4, respectively).

##  II Results

###  II.1 Constrained PIMD framework

While the evaluation of the Matsubara Green’s function (Eq. 1) is challenging,
PIMD allows one to compute also Kubo-transformed correlation functions [21],
such as
$\tilde{G}_{ij}=-\beta\sqrt{m_{i}m_{j}}\langle\delta\tilde{x}_{i}\delta\tilde{x}_{j}\rangle$,
where
$\delta\tilde{x}_{i}=\int_{0}^{\beta}\\!d\tau~{}\delta\hat{x}_{i}(\tau)/\beta$
is the time-averaged operator, $\beta=1/(k_{B}T)$ is the inverse temperature,
and $\tilde{G}_{ij}$ the Fourier transform of Eq. (1) at zero Matsubara
frequency. In other words,
$\tilde{\mathbf{G}}=-\frac{1}{\mathbf{D}^{\textrm{harm}}+\bm{\Pi}(0)}$, where
$\mathbf{D}^{\textrm{harm}}$ represents the harmonic dynamical matrix of
“bare” phonons, and $\bm{\Pi}(0)$ is the static limit of the phonon self-
energy. Therefore, by computing the inverse $\tilde{\mathbf{G}}^{-1}$, one can
have access to squared phonon frequencies renormalized by static correlations
and thus avoid the nuisance of analytically continuing Eq. (1) to real
frequencies. Moreover, dealing with $\tilde{\mathbf{G}}^{-1}$ cures also the
second drawback of the Green’s function approach, namely the impact of finite-
sampling errors. Indeed, there is an efficient way to drastically reduce its
statistical fluctuations. This is achieved by solving the GEV:

$\left[\langle\delta\tilde{\mathbf{x}}\delta\tilde{\mathbf{x}}^{T}\rangle^{-1}\right]_{ij}W_{jk}=\omega^{2}_{k}\left[\langle\tilde{\dot{\mathbf{x}}}\tilde{\dot{\mathbf{x}}}^{T}\rangle^{-1}\right]_{ij}W_{jk},$
(2)

where $\langle\tilde{\dot{x}}_{i}\tilde{\dot{x}}_{j}\rangle$ is the Kubo-
transformed velocity-velocity correlation function, $\omega^{2}_{i}$ is the
squared frequency of the $i$-th phonon mode, and $W_{ji}$ is the eigenvectors
matrix. By extending the maximum localization criterion for determining normal
modes in classical MD [23], the GEV in Eq. (2) was first proposed in Ref. [8]
within the PIMD framework. The advantage of using Eq. (2) is illustrated in
Fig. 1(a), where the GEV eigenvalues converge remarkably faster than the ones
obtained by direct inversion of $\tilde{\mathbf{G}}$. The GEV convergence
below a given target error requires one-order-of-magnitude shorter PIMD
trajectories, thus opening the avenue to realistic systems applications.

However, issues arise when treating rotational systems using this approach.
This is rationalized by rewriting the Matsubara Green’s function at zero
frequency by means of its Lehmann representation, such that it reads:

$\small\tilde{G}_{ij}=\beta\frac{\sqrt{m_{i}m_{j}}}{Z}\sum_{l,m}\langle
l|\delta\hat{x}_{i}|m\rangle\langle
m|\delta\hat{x}_{j}|l\rangle\frac{e^{-\beta E_{m}}-e^{-\beta
E_{l}}}{E_{m}-E_{l}},$ (3)

where $Z$ is the canonical partition function, and $E_{l}$ is the energy of
the $l$-th quantum eigenstate of the system. From Eq. (3), it is apparent that
the frequencies $\\{\omega_{i}\\}_{i=1,\ldots,N}$, resulting from the
diagonalization of $\tilde{\mathbf{G}}^{-1}$ for a system with $N$ degrees of
freedom, depend not only on the $N$ principal phonon excitations, but on the
full spectrum, with higher $m\rightarrow l$ excitations (overtones) weighted
by matrix elements given by $\langle l|\delta\hat{x}_{i}|m\rangle\langle
m|\delta\hat{x}_{j}|l\rangle(e^{-\beta E_{m}}-e^{-\beta E_{l}})$. By
definition, the principal phonon modes are associated with the poles of the
Green’s function with the largest spectral weights. When the weights of the
main poles are dominant, this framework is nearly exact. However, for strongly
anharmonic systems this is no longer the case, being the spectral weight
broadened over several excited states. In particular, for rotating molecules,
the anharmonicity of librations makes them hybridize with the vibrons.
Therefore, the vibron modes computed by diagonalizing $\tilde{G}_{ij}$ are
underestimated by the mixing with much softer librations.

In our approach, we solve this issue by defining appropriate displacement
operators $\delta\bar{\mathbf{x}}$ to partially decouple vibrational and
rotational modes. These operators are effectively restricted to a subspace of
the full Fock space. Henceforth, $\delta\bar{\mathbf{x}}$ will be defined as
constrained displacements. We express them in spherical coordinates around the
molecular equilibrium positions and we expand their angular components for
small librational angles. We show that in this way both librational and vibron
modes are accurately determined by our method even in the rotor limit, when
the librational amplitude is maximized.

###  II.2 Two-dimensional rotor model

As a first benchmark application, we take the model of a rotating dimer,
described by the two-dimensional (2D) potential that reads as:

$V_{c}(r,\phi)=d_{M}\cdot\left[1-e^{-a_{M}\cdot\left(r-r_{0}\right)}\right]^{2}+c\cdot\frac{\text{sin}^{2}\left(\frac{\phi-\phi_{M}}{2}\right)}{r^{2}},$
(4)

where $r_{0}$ is the classical bond length of the dimer, oriented such as to
form an angle $\phi_{M}$ with the $x$ axis. While the vibron mode is
controlled by the $d_{M}$ and $a_{M}$ parameters, the amplitude of the
libration is tuned by the strength $c$ of the external confining angular field
(see Fig. 1(b) for a plot of the potential as a function of $c$). The other
parameters are kept fixed. We set $d_{M}$=0.15 Hartree, $a_{M}$=1.1 Bohr-1 and
$r_{0}$=1.4 Bohr. Then, we numerically solved the corresponding 2D Schrödinger
equation. We note that the rotational symmetry of the problem allows one to
identify two quantum numbers, the radial and the angular one. In Fig. 1(c),
$E_{l,0}$ ($E_{0,l}$) indicates the $l$-th excited state in the vibrational
(rotational) channel.

We first obtain the GEV frequencies by solving Eq. (2) with the unconstrained
displacements $\delta\mathbf{x}=\mathbf{x}-\langle\mathbf{x}\rangle$, where
$\mathbf{x}$ is expressed in polar coordinates, such that its components read:

$\displaystyle x_{1}$ $\displaystyle=$ $\displaystyle r\
\text{cos}(\phi_{M}+\Delta\phi),$ $\displaystyle x_{2}$ $\displaystyle=$
$\displaystyle r\ \text{sin}(\phi_{M}+\Delta\phi),$ (5)

with $\Delta\phi=\phi-\phi_{M}$. From the upper panel of Fig. 1(c), one can
immediately see that the GEV eigenvalues of $\tilde{\mathbf{G}}^{-1}$ do not
reproduce the exact vibron frequency for small enough $c$. Indeed, already for
a moderate libration regime ($c<0.2$), the vibron mode is underestimated due
to the hybridization with the libration mode, which is instead very well
reproduced. To illustrate the origin of this failure in a more quantitative
way, we study the spectral function of the Green’s function in Eq. (1). For
this simple 2D rotor model, it can be computed exactly, and it is plotted in
Fig. 1(e) in light green lines. Fig. 1(e) shows that, as $c$ gets smaller,
more librational excitations become lower in frequency and larger in strength
than the first vibron peak, indicated by an arrow in the plot. The vibron peak
coincides with the transition from the ground state to the first vibrational
level (E1,0-E0,0). The corresponding vibron eigenvalue yielded by Eq. (2) is
thus a mix between the physical vibron mode and the rotational states whose
spectral weights are larger than (or comparable with) the vibron and for which
E${}_{0,J}<$ E1,0. This mixing is clearly seen also in the shape of the ground
state distribution, shown in Fig. 1(d). On the other hand, in the opposite
limit of $c\rightarrow 1$, describing the dimer strongly confined along a
preferred direction, one recovers the correct results for the fundamental
modes because, in this limit, the first two strongest peaks in the spectral
function are those related to $E_{1,0}$ and $E_{0,1}$. Thus, the mixing
between librational and vibrational sectors is suppressed and the result is
unbiased.

Figure 2: Geometries of the structures analysed in this work. (a) Top view of
the supercells. (b) Side view. All the supercells contain 96 hydrogen atoms.
The plots are generated with a cutoff for the bonds equal to 2.2 Bohr.

The failure for small coupling strength $c$ in Eq. (4) is fixed by the
_constrained_ displacement operators
$\delta\bar{\mathbf{x}}=\bar{\mathbf{x}}-\langle\bar{\mathbf{x}}\rangle$,
where $\bar{\mathbf{x}}$ in a 2D setting reads:

$\displaystyle\bar{x}_{1}$ $\displaystyle=$ $\displaystyle r\
\text{cos}(\phi_{M})-r\ \text{sin}(\phi_{M})\text{sin}(\Delta\phi),$
$\displaystyle\bar{x}_{2}$ $\displaystyle=$ $\displaystyle r\
\text{sin}(\phi_{M})+r\ \text{cos}(\phi_{M})\text{sin}(\Delta\phi).$ (6)

The displacements $\delta\bar{\mathbf{x}}$ defined by Eq. (II.2) have been
obtained from Eq. (II.2) by performing an expansion around the average angle
of rotation $\phi_{M}$ up to the second order in $\Delta\phi$ for the “slow”
cosine dependence, while keeping the full $\Delta\phi$ dependence in the sine
component. This is enough to suppress the hybridization in the relevant matrix
elements of Eq. (3), and get a very accurate phonon determination for both
librational and vibron modes through Eq.(2), as shown in the lower panel of
Fig. 1(c). The effect of the _constrained_ displacements is apparent in the
spectral functions of Fig. 1(e) (shown in blue color), where the secondary
rotational poles having larger strength than the vibron are suppressed. This
is related to the fact that, by using the constrained displacement operators,
the matrix elements that couple the vibron mode with these excited rotational
levels become negligible.

###  II.3 Phase III

After the 2D rotor model, we benchmark our method in a realistic setting, by
investigating the vibron modes of phase III of solid hydrogen. In the three-
dimensional molecular case, the strategy to constrain the displacement
operators is described by the following steps:

  1. 1.

we read the atomic Cartesian coordinates from the output trajectories of the
PIMD simulations;

  2. 2.

once the molecules are identified in the simulation cell, we express their
positions in spherical coordinates;

  3. 3.

we go back from spherical to Cartesian coordinates, by defining the new
constrained positions for each molecule in the following way:

$\bar{x}_{i}^{(\pm)}=x_{c.m.,i}\pm\frac{\bar{x}_{R,i}}{2},$ (7)

where in the above notation “$\pm$” identifies the two atoms of the dimer,
$c.m.$ is its center of mass and $\bar{x}_{R,i}$ are obtained from the
constraining relations:

$\begin{split}x_{R,1}&=r\cdot\text{sin}(\theta)\text{cos}(\phi)\\\ &\eqsim
r\cdot\left(\text{sin}(\theta_{M})+\text{cos}(\theta_{M})\text{sin}(\Delta\theta)\right)\cdot\\\
&\ \ \ \ \ \ \
\cdot\left(\text{cos}(\phi_{M})-\text{sin}(\phi_{M})\text{sin}(\Delta\phi)\right)\equiv\bar{x}_{R,1},\\\
x_{R,2}&=r\cdot\text{sin}(\theta)\text{sin}(\phi)\\\ &\eqsim
r\cdot\left(\text{sin}(\theta_{M})+\text{cos}(\theta_{M})\text{sin}(\Delta\theta)\right)\cdot\\\
&\ \ \ \ \ \ \
\cdot\left(\text{sin}(\phi_{M})+\text{cos}(\phi_{M})\text{sin}(\Delta\phi)\right)\equiv\bar{x}_{R,2},\\\
x_{R,3}&=r\cdot\text{cos}(\theta)\\\ &\eqsim
r\cdot\left(\text{cos}(\theta_{M})-\text{sin}(\theta_{M})\text{sin}(\Delta\theta)\right)\equiv\bar{x}_{R,3},\\\
\end{split}$ (8)

where $\Delta\theta=\theta-\theta_{M}$, $\Delta\phi=\phi-\phi_{M}$, being
$\theta_{M}$ and $\phi_{M}$ the polar and azimuthal angles of the molecular
axis at equilibrium.

Due to the nonlinear relations between Cartesian and spherical coordinates,
this operation should be performed explicitly over all the beads, and not on
the PIMD centroids. Therefore, the PIMD displacement operator is then computed
as
$\delta\bar{x}_{i}=\frac{1}{P}\sum_{j=1}^{P}\bar{x}_{i}^{(j)}-\langle\bar{x}_{i}\rangle$,
where $P$ is the number of beads in the PIMD simulation [18]. Although the
constraints in Eqs. (7) and (8) are defined for molecular dimers, it is worth
to note that this procedure can be easily generalized to an arbitrary
structure. Indeed, the constraints can be implemented in any situation once
the preferential vibrational axis, around which the slow cosine-dependence is
expanded, has been chosen. In the dimer case, this choice is dictated by the
direction of the bond linking the two atoms forming the dimer. In the most
general case instead, one should choose that axis according to the preferred
symmetries of the system.

Figure 3: Experimental measurements vs theoretical predictions for phase III.
(a) Experimental IR vibron peaks, taken from [24] at 190 GPa and from [25] at
247 GPa, versus theoretical IR anharmonic peaks of C2/c-24 (200 and 240 GPa)
and Cmca-12 (240 GPa) geometries. (b) Experimental Raman vibron peaks, taken
from [26], compared with theoretical Raman anharmonic peaks of C2/c-24 (200,
240 and 280 GPa) and Cmca-12 geometries (240 and 280 GPa). (c) Experimental
infrared [24] and (d) Raman [26] vibron peaks as function of pressure compared
with theoretical predictions of the C2/c-24 geometry at the harmonic level
(red), SSCHA (grey) and PIMD (green). SSCHA values are taken from [16]. Here
and in the following figures, literature data are digitalized using
WebPlotDigitizer [27].

Phase III is stable in a range of pressures between 160 GPa and 420 GPa at
temperatures lower than 200 K [24]. The well established symmetry for phase
III is the C2/c-24 geometry [28, 29, 16, 7] (see Fig. 2). Thus, we perform
three different simulations for the C2/c-24 at 200, 240 and 280 GPa, by
running PIMD with BLYP-driven forces. The BLYP functional [30, 31] has already
been proved to give results in reasonable agreement with diffusion Monte Carlo
calculations in this range of pressures [29, 32]. For testing purposes, we
examine also the Cmca-12 competing configuration [33, 29] at 240 and 280 GPa.
Computed vibron spectra with full inclusion of NQE are shown in Figs. 3(a) and
3(b), where we compare the anharmonic IR and Raman peaks with experiments. In
this comparison, we weighted the vibron modes using Born effective charges and
Raman tensor to reproduce the IR and Raman intensities, respectively.
Theoretical peaks are delta-like because, by means of our approach, we have
access to the static self-energy, without lifetime effects. We observe a nice
agreement between the C2/c-24 IR and Raman peaks and the experimental ones,
even if the frequencies of the computed spectra are still underestimated by
$\sim$100 cm-1. Furthermore, a large frequency difference is found between the
C2/c-24 and Cmca-12 geometries, with the latter in clear disagreement with the
experiments. This benchmark case shows that our PIMD phonons approach has the
predictive power to attribute the right symmetry to a given hydrogen phase,
from a genuine comparison with the experimental vibrons.

In Figs. 3(c) and 3(d) we report the pressure behaviour of the IR and Raman
experimental peaks, respectively, and we compare them with the C2/c-24 peaks
evaluated at different levels of theory. We observe a strong renormalization
of the anharmonic PIMD vibron frequencies with respect to the harmonic case.
Also, the vibron energy slope is correctly reproduced only when NQE are taken
into account. Furthermore, we find that our PIMD vibron frequencies are closer
to the experiment than the ones computed by the time dependent SSCHA framework
using the same BLYP functional [16]. This is due to the SSCHA Gaussian ansatz
limitations in presence of strong librational modes, and it is rationalized by
a detailed analysis of the 2D dimer model using the SCHA framework (see
Supplementary Information [34]).

###  II.4 Phase IV

We then addressed the open problem of determining the right symmetry for phase
IV [25, 35, 36, 37]. Experimentally, it is stable above 300 K and 220 GPa [38,
39]. A number of theoretical works [40, 41, 42, 43, 44] led to a list of
possible candidates for its structural characterization.

Figure 4: Experimental measurements vs theoretical predictions for phase IV.
(a) Top panel: Experimental IR and Raman vibron peaks for phase IV at 280 GPa.
Raman experimental data are taken from [38] while infrared experimental from
[35]. Lower panels: PIMD vibron modes for the four different geometries:
Pc-48; Pbcn-48; Ama2-24; Ibam-8. Continuous lines represent modes projected
over the molecular layers, while points are projections over atomic-like
layers. (b) Experimental infrared [35] and Raman [38] vibron peaks as function
of pressure compared with theoretical predictions of the Ama2-24 geometry at
the harmonic level (red) and PIMD (dark grey). (c) Experimental IR and Raman
vibron peaks for deuterium phase IV at 280 GPa, extrapolated by fitting the
slope of the vibron energy as a function of pressure. Raman experimental data
of Raman vibron are taken from [38] while infrared experimental from [36].
Lower panels: anharmonic theoretical vibron modes for Pc-48 and Ama2-24.

Here, we consider four different geometries, already introduced in previous
works, that are all characterized by alternating layers of strongly bonded H2
molecules and atomic-like layers (G-layers): Pc-48 [41], having a strongly
distorted honeycomb lattice as G-layer; Pbcn-48 [28], where the honeycomb
lattice is less distorted and slightly dimerized; Ama2-24 [44], where the
G-layer has a typical C2/c arrangement; Ibam-8 [28] with a perfectly hexagonal
G-layer. See Fig. 2 for their graphical representation. From previous PIMD
results [45], confirmed by our simulations, it turns out that these structures
are clearly quantum crystals [46], with spatial fluctuations even larger than
in phase III. At variance with phase III, where the H2 molecules librate but
never perform full $\pi$-rotations, the molecules belonging to the molecular
layers of phase IV are free to rotate around their center of mass. They behave
very similarly for all four geometries, with an average bond length of 1.36
Bohr in every case (see Tab. 2 in [34]). Moreover, high temperatures lead to
an additional drift of their center of mass. Therefore, we study their vibron
modes at 280 GPa and 20 K, thermodynamic conditions that allowed us to compute
less noisy PIMD phonons. Nevertheless, quantum fluctuations are strong enough
to place these structures in a regime where IR and Raman matrix elements are
no longer reliable within a perturbative approach. Thus, for phase IV we do
not weight the PIMD vibron modes by these matrix elements. We show them all in
the lower panels of Fig. 4(a) for the four considered geometries and we
finally compare their frequencies with corresponding experimental IR and Raman
values (top panel of Fig. 4(a)).

##  III Discussion

Experimentally, phase IV is believed to be made of layers of strongly bonded
molecules alternated by graphene-like sheets [41]. A signature of this
structure would be a strong softening of the G-layer vibron with increasing
pressure. Therefore, in order to study the character of the vibron modes of
the simulated structures, we project these modes over both molecular and
atomic-like layers for each system. A striking feature is the clear separation
between the vibron frequencies belonging to the molecular layers from the
atomic ones. Interestingly enough, quantum anharmonicity not only renormalizes
the full spectrum, but also widens the gap between these two sets of vibrons
with respect to the harmonic approximation. This effect is present in all
geometries analyzed, being the mildest in Ibam-8. From this analysis, we
deduce that both IR peaks measured in experiments and reported in the top
panel of Fig. 4(a) should necessarily be localized over molecular layers. This
is at variance with a former interpretation, which assigned the softer IR peak
to atomic-layer modes [36]. The rationale behind it is that strong NQE smooth
out the structural differences between competing arrangements, particularly in
the atomic layer, and make it closer to an undistorted hexagonal shape (see
Fig. 7 in [34]). This is in accordance to the fact that, for a perfect
graphene-like arrangement, the IR response is absent. From Fig. 4(a), by
comparing the predicted vibron frequency distribution with the experimental
one, it turns out that Ama2-24 is the best-matching geometry for phase IV.

Table 1: PIMD expectation value of the Hamiltonian ($\langle H\rangle$) vs static DFT energies. The energies with quantum nuclei are computed using the virial estimator [18] for the kinetic contribution. All the PIMD simulations are performed at 20 K and using 120 beads. Energies are given in eV and are normalized by the number of atoms in the simulation cells. Structure | Pressure (GPa) |  $\langle H\rangle$ PIMD | static DFT
---|---|---|---
$C2/c$–$24$ | 200 | -14.2816(1) | -14.5088
$C2/c$–$24$ | 240 | -14.0890(2) | -14.3254
$C2/c$–$24$ | 280 | -13.9095(1) | -14.1535
$Cmca$–$12$ | 240 | -14.0672(1) | -14.3091
$Cmca$–$12$ | 280 | -13.8912(2) | -14.1398
$Pbcn$–$48$ | 280 | -13.9054(2) | -14.1467
$Pc$–$48$ | 280 | -13.9089(2) | -14.1504
$Ibam$–$8$ | 280 | -13.8740(6) | -14.1070
$Ama2$–$24$ | 280 | -13.9103(6) | -14.1512

A further verification of the symmetry attribution for phase IV comes from
studying the static and fully anharmonic lattice energies in Tab. 1. Indeed,
for a pressure equal to 280 GPa, Ama2-24 is the most competitive among phase
IV candidates both at the static and anharmonic lattice level, followed by
Pc-48. We notice however that, from a structural perspective, NQEs make the
difference between these various symmetries much milder than the ones present
in static configurations.

Based on the above analysis, we chose the Ama2-24 geometry to perform another
simulation at different pressure (300 GPa) to inspect if, like for phase III,
the pressure dependence follows the slope of the experimental vibron energies.
We show the results in Fig. 4(b), where we compare the harmonic and PIMD
vibron frequencies (in red and black colours, respectively) with the
experimental values. We observe that the lowest Raman peak is caught only by
PIMD frequencies and that its slope is correctly described by anharmonic
vibron modes of atomic character. All PIMD vibron frequencies above are of
molecular character, and they are spread over an energy interval compatible
with the experimental findings over the full pressure range analyzed here.

To confirm the scenario, we also performed a simulation of Ama2-24 and Pc-48
using deuterium in place of hydrogen, at 280 GPa. At this pressure and with
the deuterium mass, both symmetries have the tendency to transform into the
C2/c-24 geometry for low enough temperatures. Ama2-24 is stable at 150 K,
while at 20 K we have been able to run Pc-48 simulations for 4.5 ps, a time
sufficiently long to estimate the PIMD phonons. As shown in Fig. 4(c), also in
this case Ama2-24 seems to best fit the experiments. The underestimation of
the softer Raman peak predicted in Ama2-24 is consistent with what found in
C2/c-24 for phase III.

In summary, in this work we introduced an accurate method for computing
phonons of molecular solids strongly affected by NQE. The framework proposed
here paves the way towards a robust assessment of crystalline symmetries, even
in cases where structural information comes primarily from vibrational
spectroscopy. This is particularly relevant in quantum crystals, where the
usual harmonic approximation badly fails, preventing a direct comparison with
experiments. The present approach is suitable to better quantify their unique
properties, such as polaron effects and superconductivity in hydrogen-rich
materials.

##  IV Methods

The PIMD simulations are carried out at 20 K using 120 beads to take into
account quantum effects. Nuclei are evolved in time using the PIOUD integrator
[18] with a time step equal to 0.5 fs and a friction parameter of the Langevin
thermostat equal to 1.46$\cdot$10-3 atomic units. The latter value is the same
as in Ref. [18], where it is found to be optimal for both stochastic and
deterministic forces. Simulations lasted around 10 ps, until the convergence
on vibron modes at Gamma is reached. Forces are computed from the Born
Oppenheimer potential energy surface (PES) evaluated at Density Functional
Theory level within the Quantum Espresso [47] engine. In particular we used a
BLYP functional for computing the PES and an LDA functional for the evaluation
of the Born effective charges and Raman tensors. The choice of the BLYP
functional comes from the benchmark tests using Quantum Monte Carlo
calculations performed in Ref. [48] and also from the fact that, between the
different available functionals, geometries computed at BLYP level give the
lowest Diffusion Monte Carlo energies [32, 29]. The wavefunction cut-off for
the PES is set to 60 Ry (300 Ry for the charge density), while the Fermi
smearing is Gaussian and set equal to 0.03 Ry. PIMD simulations are performed
using supercells containing in each case 96 hydrogen atoms and the
corresponding supercell reciprocal space mesh is always equal to 4x4x4. For
phase III, Raman active modes are evaluated from the response function in the
following way [16]:

$I_{Raman}(\omega)\propto-\sum_{\alpha=1}^{3}\sum_{a,b}\frac{A_{\alpha\alpha,a}A_{\alpha\alpha,b}}{\sqrt{m_{a}m_{b}}}\
\text{Im}\left[G_{a,b}(\omega)\right],$ (9)

where $a,b$ run over all the 3N degrees of freedom, $A$ is the Raman tensor
that we derived from static DFT calculations using the Quantum Espresso suite
[47], $m_{i}$ is the mass of the i-th degree of freedom, Im is the imaginary
part, $G^{-1}_{a,b}(\omega)=\omega^{2}-\Phi_{a,b}$ and $\Phi_{a,b}$ is the
force constant matrix at $\Gamma$ that we build from PIMD simulations,
therefore already including the renormalization due to anharmonic effects.
Infrared active modes are instead computed using the Born effective charge
tensors $Z^{*}$, with the expression:

$I_{IR}(\omega)=\sum_{\alpha,\beta=1}^{3}\sum_{a,b}\frac{Z^{*}_{\alpha,a}Z^{*}_{\beta,b}}{\sqrt{m_{a}m_{b}}}\
\text{Im}\left[G_{a,b}(\omega)\right].$ (10)

In particular, the points reported in Fig. 3 are computed by averaging the
frequencies with the spectra obtained from Eqs. (9) and (10), i.e:

$\begin{split}&\omega_{Raman,IR}=\frac{\int d\omega\ \omega\cdot
I_{Raman,IR}(\omega)}{\int d\omega\ I_{Raman,IR}(\omega)}.\end{split}$ (11)

For phase IV, we find that the above methods to estimate IR and Raman tensors
fail because the system is strongly anharmonic and the small displacement
around equilibrium positions employed to estimate tensors are no longer
reliable. In Figs. 4, phonons at $\Gamma$-point of the Brillouin zone are
projected over the layer of strongly bonded molecules (solid lines) and over
the atomic-like layer (points). The projector operator is built by
orthonormalizing the NP, where NP is the number of projections, 3N-dimensional
vectors $|v^{proj}\rangle$, representing the atoms of each H2 molecule over
which we project. Then, the reduced eigenvectors are found by:

$\phi^{proj}_{j,a}=\sum_{b}|v^{proj}_{j,b}\rangle\langle
v^{proj}_{j,b}|\Phi_{a,b}\rangle.$ (12)

The projected density of states (DOS) is finally obtained using the following
expression:

$PDOS(\omega)=\sum_{a=1}^{3N}\langle\phi^{proj}_{a}|\phi^{proj}_{a}\rangle\delta(\omega-\bar{\omega}_{a}).$
(13)

## Acknowledgements

T.M. and M.C. thank the French grand équipement national de calcul intensif
(GENCI) for the computational time provided under Project No. 0906493. All
authors thank Lorenzo Paulatto, Lorenzo Monacelli and Michele Lazzeri for
useful discussions.

## Funding

This work is supported by the European Centre of Excellence in Exascale
Computing TREX (Targeting Real Chemical Accuracy at the Exascale). This
project has received funding from the European Union’s Horizon 2020 (Research
and Innovation program) under grant agreement no. 952165.

## Author Contributions

TM coded the numerical algorithms and performed the calculations. All authors
developed the theory and analysed the data. TM wrote the paper with
contributions from MC and RV. MC proposed and led the project.

## Competing Interests

The authors declare no competing interests.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Code Availability

The codes implementing the calculations of this study are available from the
corresponding author upon request.

## References

  * Mao and Hemley [1994] H. K. Mao and R. J. Hemley, Ultrahigh-pressure transitions in solid hydrogen, Rev. Mod. Phys. 66, 671 (1994).
  * Loubeyre  _et al._ [1996] P. Loubeyre, R. LeToullec, D. Hausermann, M. Hanfland, R. J. Hemley, H. K. Mao, and L. W. Finger, X-ray diffraction and equation of state of hydrogen at megabar pressures, Nature 383, 702 (1996).
  * Loubeyre  _et al._ [2002] P. Loubeyre, F. Occelli, and R. LeToullec, Optical studies of solid hydrogen to 320 gpa and evidence for black hydrogen, Nature 416, 613 (2002).
  * Eremets and Troyan [2011] M. I. Eremets and I. A. Troyan, Conductive dense hydrogen, Nature Materials 10, 927 (2011).
  * McMinis  _et al._ [2015] J. McMinis, R. C. Clay, D. Lee, and M. A. Morales, Molecular to atomic phase transition in hydrogen under high pressure, Phys. Rev. Lett. 114, 105305 (2015).
  * Dias and Silvera [2017] R. P. Dias and I. F. Silvera, Observation of the wigner-huntington transition to metallic hydrogen, Science 355, 715 (2017), https://www.science.org/doi/pdf/10.1126/science.aal1579 .
  * Monacelli  _et al._ [2022] L. Monacelli, M. Casula, K. Nakano, S. Sorella, and F. Mauri, Quantum phase diagram of high-pressure hydrogen (2022), arXiv:2202.05740 [cond-mat.mtrl-sci] .
  * Morresi  _et al._ [2021] T. Morresi, L. Paulatto, R. Vuilleumier, and M. Casula, Probing anharmonic phonons by quantum correlators: A path integral approach, The Journal of Chemical Physics 154, 224108 (2021), https://doi.org/10.1063/5.0050450 .
  * Magdău and Ackland [2013] I. B. Magdău and G. J. Ackland, Identification of high-pressure phases iii and iv in hydrogen: Simulating raman spectra using molecular dynamics, Phys. Rev. B 87, 174110 (2013).
  * Zhang  _et al._ [2018] C. Zhang, C. Zhang, M. Chen, W. Kang, Z. Gu, J. Zhao, C. Liu, C. Sun, and P. Zhang, Finite-temperature infrared and raman spectra of high-pressure hydrogen from first-principles molecular dynamics, Phys. Rev. B 98, 144301 (2018).
  * Bowman [1978] J. M. Bowman, Self-consistent field energies and wavefunctions for coupled oscillators, The Journal of Chemical Physics 68, 608 (1978).
  * Monserrat  _et al._ [2013] B. Monserrat, N. Drummond, and R. Needs, Anharmonic vibrational properties in periodic systems: energy, electron-phonon coupling, and stress, Physical Review B 87, 144302 (2013).
  * Azadi  _et al._ [2014] S. Azadi, B. Monserrat, W. Foulkes, and R. Needs, Dissociation of high-pressure solid molecular hydrogen: A quantum monte carlo and anharmonic vibrational study, Physical review letters 112, 165501 (2014).
  * Errea  _et al._ [2014] I. Errea, M. Calandra, and F. Mauri, Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides, Phys. Rev. B 89, 064302 (2014).
  * Bianco  _et al._ [2017] R. Bianco, I. Errea, L. Paulatto, M. Calandra, and F. Mauri, Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation, Phys. Rev. B 96, 014111 (2017).
  * Monacelli  _et al._ [2021] L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Black metal hydrogen above 360 gpa driven by proton quantum fluctuations, Nature Physics 17, 63 (2021).
  * Ceriotti  _et al._ [2010] M. Ceriotti, M. Parrinello, T. Markland, and D. Manolopoulos, Efficient stochastic thermostatting of path integral molecular dynamics, The Journal of chemical physics 133, 124104 (2010).
  * Mouhat  _et al._ [2017] F. Mouhat, S. Sorella, R. Vuilleumier, A. M. Saitta, and M. Casula, Fully quantum description of the zundel ion: Combining variational quantum monte carlo with path integral langevin dynamics, Journal of Chemical Theory and Computation 13, 2400 (2017), pMID: 28441484, https://doi.org/10.1021/acs.jctc.7b00017 .
  * Mahan [2000] G. D. Mahan, Nonzero temperatures, in  _Many-Particle Physics_ (Springer US, Boston, MA, 2000) pp. 109–185.
  * Jarrell and Gubernatis [1996] M. Jarrell and J. Gubernatis, Bayesian inference and the analytic continuation of imaginary-time quantum monte carlo data, Physics Reports 269, 133 (1996).
  * Craig and Manolopoulos [2004] I. R. Craig and D. E. Manolopoulos, Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics, The Journal of Chemical Physics 121, 3368 (2004), https://doi.org/10.1063/1.1777575 .
  * Van Kranendonk and Kranendonk [1983] J. Van Kranendonk and J. Kranendonk,  _Solid Hydrogen: Theory of the Properties of Solid H2, HD, and D2_ (Springer US, 1983).
  * Martinez  _et al._ [2006] M. Martinez, M.-P. Gaigeot, D. Borgis, and R. Vuilleumier, Extracting effective normal modes from equilibrium dynamics at finite temperature, The Journal of Chemical Physics 125, 144106 (2006), https://doi.org/10.1063/1.2346678 .
  * Loubeyre  _et al._ [2020] P. Loubeyre, F. Occelli, and P. Dumas, Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen, Nature 577, 631 (2020).
  * Zha  _et al._ [2012] C.-S. Zha, Z. Liu, and R. J. Hemley, Synchrotron infrared measurements of dense hydrogen to 360 gpa, Phys. Rev. Lett. 108, 146402 (2012).
  * Eremets  _et al._ [2019] M. I. Eremets, A. P. Drozdov, P. P. Kong, and H. Wang, Semimetallic molecular hydrogen at pressure above 350 gpa, Nature Physics 15, 1246 (2019).
  * Rohatgi [2020] A. Rohatgi, Webplotdigitizer: Version 4.3 (2020).
  * Pickard and Needs [2007] C. J. Pickard and R. J. Needs, Structure of phase iii of solid hydrogen, Nature Physics 3, 473 (2007).
  * Drummond  _et al._ [2015] N. D. Drummond, B. Monserrat, J. H. Lloyd-Williams, P. L. Ríos, C. J. Pickard, and R. J. Needs, Quantum monte carlo study of the phase diagram of solid molecular hydrogen at extreme pressures, Nature Communications 6, 7794 (2015).
  * Lee  _et al._ [1988] C. Lee, W. Yang, and R. G. Parr, Development of the colle-salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 37, 785 (1988).
  * Becke [1988] A. D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 38, 3098 (1988).
  * Azadi  _et al._ [2017] S. Azadi, N. D. Drummond, and W. M. C. Foulkes, Nature of the metallization transition in solid hydrogen, Phys. Rev. B 95, 035142 (2017).
  * Singh  _et al._ [2014] R. Singh, S. Azadi, and T. D. Kühne, Anharmonicity and finite-temperature effects on the structure, stability, and vibrational spectrum of phase iii of solid molecular hydrogen, Phys. Rev. B 90, 014110 (2014).
  * [34] See Supplementary Information, which includes Ref. [49], for FIGS. [S.1–S.7] and for further information on: (i) Two dimensional roto-vibrational model; (ii) Details of the atomic structures; (iii) PIMD simulation analysis .
  * Zha  _et al._ [2013] C.-s. Zha, Z. Liu, M. Ahart, R. Boehler, and R. J. Hemley, High-pressure measurements of hydrogen phase iv using synchrotron infrared spectroscopy, Phys. Rev. Lett. 110, 217402 (2013).
  * Loubeyre  _et al._ [2013] P. Loubeyre, F. Occelli, and P. Dumas, Hydrogen phase iv revisited via synchrotron infrared measurements in h2 and d2 up to 290 gpa at 296 k, Phys. Rev. B 87, 134101 (2013).
  * Goncharov  _et al._ [2019] A. F. Goncharov, I. Chuvashova, C. Ji, and H.-k. Mao, Intermolecular coupling and fluxional behavior of hydrogen in phase iv, Proceedings of the National Academy of Sciences 116, 25512 (2019), https://www.pnas.org/content/116/51/25512.full.pdf .
  * Howie  _et al._ [2012a] R. T. Howie, C. L. Guillaume, T. Scheler, A. F. Goncharov, and E. Gregoryanz, Mixed molecular and atomic phase of dense hydrogen, Phys. Rev. Lett. 108, 125501 (2012a).
  * Howie  _et al._ [2012b] R. T. Howie, T. Scheler, C. L. Guillaume, and E. Gregoryanz, Proton tunneling in phase iv of hydrogen and deuterium, Phys. Rev. B 86, 214104 (2012b).
  * Liu  _et al._ [2012] H. Liu, L. Zhu, W. Cui, and Y. Ma, Room-temperature structures of solid hydrogen at high pressures, The Journal of Chemical Physics 137, 074501 (2012), https://doi.org/10.1063/1.4745186 .
  * Pickard  _et al._ [2012] C. J. Pickard, M. Martinez-Canales, and R. J. Needs, Density functional theory study of phase iv of solid hydrogen, Phys. Rev. B 85, 214114 (2012).
  * Goncharov  _et al._ [2013] A. F. Goncharov, J. S. Tse, H. Wang, J. Yang, V. V. Struzhkin, R. T. Howie, and E. Gregoryanz, Bonding, structures, and band gap closure of hydrogen at high pressures, Phys. Rev. B 87, 024101 (2013).
  * Azadi  _et al._ [2018] S. Azadi, R. Singh, and T. D. Kühne, Nuclear quantum effects induce metallization of dense solid molecular hydrogen, Journal of Computational Chemistry 39, 262 (2018), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.25104 .
  * Li  _et al._ [2020] G.-J. Li, Y.-J. Gu, Z.-G. Li, Q.-F. Chen, and X.-R. Chen, New possible candidate structure for phase iv of solid hydrogen, RSC Adv. 10, 26443 (2020).
  * Rillo  _et al._ [2018] G. Rillo, M. A. Morales, D. M. Ceperley, and C. Pierleoni, Coupled electron-ion monte carlo simulation of hydrogen molecular crystals, The Journal of Chemical Physics 148, 102314 (2018), https://doi.org/10.1063/1.5001387 .
  * Cazorla and Boronat [2017] C. Cazorla and J. Boronat, Simulation and understanding of atomic and molecular quantum crystals, Rev. Mod. Phys. 89, 035003 (2017).
  * Giannozzi et al. [2009] P. Giannozzi et al., Quantum espresso: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Clay  _et al._ [2014] R. C. Clay, J. Mcminis, J. M. McMahon, C. Pierleoni, D. M. Ceperley, and M. A. Morales, Benchmarking exchange-correlation functionals for hydrogen at high pressures using quantum monte carlo, Phys. Rev. B 89, 184106 (2014).
  * McMahon  _et al._ [2012] J. M. McMahon, M. A. Morales, C. Pierleoni, and D. M. Ceperley, The properties of hydrogen and helium under extreme conditions, Rev. Mod. Phys. 84, 1607 (2012).

## SUPPLEMENTARY INFORMATION

### Tommaso Morresi,1∗ Rodolphe Vuilleumier,2 Michele Casula1∗

1Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC),
Sorbonne Université, CNRS UMR 7590, MNHN, 4 Place Jussieu, 75252 Paris, France
2PASTEUR, Département de chimie, École normale supérieure, PSL University,
Sorbonne Université, CNRS, 75005 Paris, France

### Two-dimensional roto-vibrational model

The shape of the potential in Eq. (4) for three different values of $c$ is
presented in Fig. 1(b). In particular, the limit $c\to 0$ (c$\sim$0.005)
corresponds to the case in which the dimer is free to rotate, while the large
$c$ limit (c$\sim$0.5) describes the molecule confined along one preferred
direction. In the intermediate range (c$\sim$0.05), the dimer cannot perform a
full $\pi$-rotation but the amplitude of the librational motion is more or
less wide depending on the value taken by the $c$ parameter. As we have done
in Ref. [8], we carry out the analysis by solving numerically the Schrödinger
equation at low temperature (20 K) to minimize thermal effects and to focus
purely on quantum effects. The details of the numerical scheme can be found in
Appendix E of Ref. [8].

Figure 5: Eigenvalues of the Schrödinger equation for the potential in Eq.
(4). In particular we highlight the ground state level (red), the first
excited state (orange) which is a rotational state in the range of $c$
considered here, the first vibrational excited state (brown) and all the
remaining states which have rotational nature (purple and olive).

As seen in the main text (top panel Fig. 1(c)), using the _unconstrained_
displacement operators one is able to reproduce very well the rotational
frequency (orange line in Fig. 1(c)). However, for $c<$ 0.2 the computed
vibron frequency (brown line in Fig. 1(c)) is biased. Based on the eigenstates
decomposition (Eq. (3) of the main text and Eq. (26) of Ref. [8]), another
explanation for this failure can be given by looking at Figs. 5, reporting the
eigenvalues of the system.

Figure 6: (a) SCHA eigenvalues obtained with the auxiliary hessian matrix and
the second derivative of the SCHA free energy for different values of $c$ in
Eq. (4). (b) Density plot of the square modulus of the ground state
wavefunction ($|\Psi_{0}|^{2}$) versus SCHA wavefunction
($|\Psi^{\text{(SCHA)}}_{0}|^{2}$) for the points reported in panel (a).

Indeed, while one would expect to find the transition between the first
rotational state and the ground state (E0,1 \- E0,0) and the transition
between the first vibrational excited state and the ground state (E1,0 \-
E0,0) as fundamental modes, the presence of rotational states in between for
$c<$0.2 clearly biases the evaluation of the vibron mode (i.e., E1,0 \- E0,0).
By inspecting the matrix elements coupling different excitations, one can see
that, even if the transition from state $(1,0)$ to the state $(0,0)$ has the
largest weight, the matrix elements coupling higher rotational states (namely,
states $(0,J)$ with $J>1$) with the ground state yield a spurious contribution
to the calculated vibron frequency for $c<$0.2.

Finally, we would like to compare our framework with the Self Consistent
Harmonic Approximation (SCHA) [14, 15] results on the very same dimer model.
Indeed, the stochastic SCHA approach provides the main theoretical reference
for phase III, beside our calculations (Figs. 3(c) and 3(d)). This approach is
based on a variational principle for the free energy [15] and includes NQE as
well as PIMD. In Fig. 6 we report the behaviour of the SCHA eigenvalues for
the auxiliary Hessian matrix ($\omega^{aux}_{SCHA,I}$ and
$\omega^{aux}_{SCHA,II}$) and for the second derivative of the SCHA free
energy ($d^{2}F_{SCHA}/dR^{2}|_{I}$ and $d^{2}F_{SCHA}/dR^{2}|_{II}$) for a
few $c$ values in Eq. (4). In Fig. 6(a) we can observe a small underestimation
of the vibron mode using the SCHA free energy’s second derivative in the
region $c<$0.5, while for $c=$0.5 the vibron eigenvalue is correctly
reproduced. This behaviour is understood by the SCHA Gaussian ansatz
wavefunction, that cannot recover the full rotational character of the true
ground state wave function for small coupling parameters. This can be clearly
seen in Fig. 6(b). On the other hand, in the limit of the molecule strongly
confined along one direction ($c\geq$0.5 in our model), the SCHA wavefunction
reproduces (in a variational sense) very well the true ground state
wavefunction and the eigenvalues are accurately estimated.

### Details of the atomic structures

In Fig. 2 of the main text, we report a picture of the initial supercell
configurations for all the geometries studied in this work.

Figure 7: Top view snapshots of two different layers of mixed structures. Top
left: Pc-48. Top right: Pbcn-48. Bottom left: Ama2-24. Bottom right: Ibam-8.
The colours employed to draw the highlighted G-layers are the same of Fig.
4(a) of the main text. Apart from the different lattice symmetries, by
plotting all the beads an overall similarity between the geometries can be
recognized.

The details of their unit cells are given in the following. Purely molecular
structures (phase III):

  * •

C2/c-24: base-centered monoclinic structure and four different layers,
arranged in an ABCD fashion, of H2 molecules ordered in distorted hexagonal
rings;

  * •

Cmca-12: base-centered orthorhombic lattice and two different layers of H2
molecules arranged in an AB fashion;

Mixed structures (phase IV):

  * •

Pc-48: simple orthorombic lattice, four layers stacked in an ABAB fashion and
the G-layers made by strongly distorted hexagons (corresponding to weak bonded
molecules) [41];

  * •

Pbcn-48: simple orthorombic lattice, four layers stacked in an ABAB fashion
and the G-layers made by slightly distorted hexagons (corresponding to very
weak bonded molecules) [28];

  * •

Ama2-24: base-centered orthorombic lattice, four layers stacked in an ABCD
fashion and the G-layers very similar to the C2/c layers [44];

  * •

Ibam-8: body-centered orthorhombic lattice, four different layers stacked in
an ABAB fashion and perfectly hexagonal G-layers [28].

Finally, in Fig. 7 we exhibit a top-view plot of the mixed structures, with
all beads gathered together at given time (that we choose in the middle of the
trajectory) of the PIMD simulation. We highlight the G-layer atoms with the
same colours of Fig. 4(a) of the main text. This picture is interesting
because it clearly shows that, if one focus just on the instantaneous
configurations, quantum effects hide the differences between competing
geometries of phase IV.

### PIMD simulations analysis

In this section we add further comments to Tab. 1 of the main text, we analyze
the average bond length of the molecules and the pair distribution functions
for all the simulations carried out in this work. We then report the DOS for
the C2/c-24 and Cmca-12 geometries and the vibron modes of mixed structures to
analyse the impact of Nuclear Quantum Effects (NQE) with respect to the
harmonic theory.

Table 2: Average bond lengths of the H2 molecule in the different structures. Lengths are given in Bohr units. Structure | Pressure (GPa) | Bond length
---|---|---
$C2/c$–$24$ | 200 | 1.3968
$C2/c$–$24$ | 240 | 1.4053
$C2/c$–$24$ | 280 | 1.4145
$Cmca$–$12$ | 240 | 1.4357
$Cmca$–$12$ | 280 | 1.4443
$Pbcn$–$48$ | 280 | 1.3684
$Pc$–$48$ | 280 | 1.3685
$Ibam$–$8$ | 280 | 1.3693
$Ama2$–$24$ | 280 | 1.3690
Figure 8: Layered pair distribution functions for the structures studied in
this work. (a) Mixed structures (Pbcn-48, Pc-48, Ibam-8 and Ama2-24); (b)
C2/c-24; (c) Cmca-12. In (a) we plot both the g(r) belonging to atomic-like
planes (using points) and the g(r) belonging to molecular planes (solid
lines). A similar analysis but at different temperatures is done in Ref. [45].

In Tab. 1 of the main text we can observe that, for the C2/c-24 and Cmca-12
geometries, the energy expectation value decreases as the pressure increases
and that, at 280 GPa and 20 K, the lowest energy structure is the C2/c-24, as
already found in different works [28, 49, 29]. It is also worth to note that
we have attempted to simulate the Cmca-4 geometry [29] at the same
thermodynamic conditions, but after a few ns we observe a phase transition
into a mixed layered structure, very close to the Pbcn-48. Therefore, we
evince that the Cmca-4 structure is unstable at 280 GPa and 20 K using the
BLYP functional. The same finding is reported in Ref. [45] at 200 K. In Tab. 2
we report the average bond lengths of H2 molecules within the different
geometries (Fig. 2). For the mixed structures, i.e. Pbcn-48, Pc-48, Ibam-8 and
Ama2-24. This analysis is performed by taking into account only layers of
strongly bonded molecules. We observe that the bond length increases with
pressure in the C2/c-24 and Cmca-12 samples. Molecular bond lengths of the
mixed structures are very close to each other and significantly smaller than
the ones in purely molecular structures (C2/c-24 and Cmca-12). This can be
rationalized also by looking at the pair distribution function g(r) in Fig. 8.
In Fig. 8(b) and 8(c) we observe a small shift of the first peak to the left
by increasing pressure for the C2/c-24 and Cmca-12 geometries respectively. In
Fig. 8(a) we instead observe an almost exact overlap for the first peak,
denoting that the molecular layers of the four different mixed structures
behave similarly. In the latter case instead we denote a difference in the
atomic layer (points in Fig. 8(a)). Indeed, in that case while the Ibam-8 g(r)
has a single peak, meaning that during the simulation the hexagonal symmetry
of the plane is conserved, for the Pbcn-48, Pc-48 and Ama2-24 we observe a
double peak reflecting the presence of different distances between atoms in
the atomic-layer. Furthermore, the curves corresponding to the atomic planes
Pbcn-48 and Pc-48 overlap, making the two structures very hard to distinguish
during the PIMD simulation.

Figure 9: Theoretical phonon DOS of the C2/c-24 geometry for three different
pressures: (a) 200 GPa; (b) 240 GPa; (c) 280 GPa. In red we report the
harmonic DOS while in green the anharmonic DOS from PIMD simulations. All the
DOS are obtained by interpolating the four q-points sampled by the supercell
choice and smoothed with a normalized Lorentzian function having a linewidth
equal to 10 cm-1. Figure 10: Theoretical phonon DOS of the Cmca-12 geometry
for (a) 240 GPa and (b) 280 GPa. As in Fig. 9 we report the harmonic DOS in
red. The anharmonic DOS are given in dark-green. Figure 11: Anharmonic
theoretical vibron modes vs harmonic vibron modes (red) for the four different
geometries: (a) Pc-48; (b) Pbcn-48; (c) Ama2-24; (d) Ibam-8. Continuous lines
represent modes projected over the strongly bonded molecular layers, while
points projections over atomic-like layers.

The harmonic (red) and anharmonic DOS for the pure molecular structures
(C2/c-24 and Cmca-12) are reported in Figs. 9 and 10. We observe in both cases
that, while the lattice modes (modes at frequencies below 2500 cm-1) are not
heavily affected by NQE, vibron modes are strongly renormalized with respect
to the harmonic case. A similar finding was reported in [16] for the C2c/-24.

Finally, a remarkable softening of $\Gamma$-point PIMD phonons with respect to
the harmonic case is observed for the mixed-layers structures describing phase
IV in Fig. 11.


BA, Hiram College, Hiram, Ohio, 2012

Physics

Doctor of Philosophy

December

2020

# DATA ANALYSIS FOR A BALLOON BORNE COMPTON POLARIMETER

Sambid Kumar Wasti

(2020)

This dissertation has been examined and approved in partial fulfillment of the
requirements for the degree of Doctor of Philosophy in Physics by:

Dissertation Director, Mark McConnell,
Professor of Physics & Astronomy

Joe Dwyer,
Professor of Physics & Astronomy

Karsten Pohl,
Professor of Physics & Astronomy

Fabian Kislat,
Assistant Professor of Physics & Astronomy

Lynn Kistler,
Professor of Physics & Astronomy

on 11/25/2020.

Original approval signatures are on file with the University of New Hampshire
Graduate School.

This dissertation is dedicated to my parents.

First and foremost, I would like to thank my advisor, Dr. Mark McConnell, for
his guidance through out my research years at UNH. This would not have been
possible without his guidance, advice and patience. I would also like to
acknowledge and thank Camden Ertley, Peter Bloser, James Ryan and Jason Legere
for their help with the GRAPE and the balloon campaign. I would also like to
thank Taylor Connor, Chris Bancroft, Amanda Madden, and rest of the GRAPE
collaboration. I would also like to thank my committee members, Dr. Joseph
Dwyer, Dr. Karsten Pohl, Dr. Lynn Kister and Dr. Fabian Kislat for their
insightful remarks and guidance to make this thesis a reality.

Words cannot express how grateful I am to my family. My mom Sushama, who has
been more worried about my thesis completion than a graduate student like
myself. And my dad Sharat, whose wisdom has guided me during difficult times,
throughout my life. My brother Shashanka and my cousins, Reshmi, Asphota,
Ayush, Ashish, Abhijan, Prapti and Prafulla who have continuously reminded me
to be grounded and be myself. I would also like to thank my friends Jason
Shuster, Cristian Ferradas, Kate Luksha, Tejaswita Sharma, Karla Oñate, Nick
Lubinsky, Pranav Shakya, and Aaditya Raj Bhandari. The board game nights and
the philosophical discussions were an escape to life outside research and
classes. All of this would not have been possible without the support of my
friends and family.

Finally, I would like to acknowledge the funding that made my work possible:
NASA NNX13AB96G, NASA NNX13AB96G, and NASA 80NSSC19K0624. I would also like
thank CUA for supporting me during the last few months of completion of this
thesis.

The Gamma RAy Polarimeter Experiment (GRAPE), a balloon borne Compton
polarimeter for 50-500 keV gamma rays, was successfully flown for the second
time in 2014. GRAPE uses events in which a photon scatters between 2 detector
elements to measure the polarization of the measured photons. GRAPE consists
of 24 collimated polarimeter modules. Each module is made up of 64 rectangular
scintillators elements in a grid of 8 $\times$ 8 (36 plastic scintillators
surrounded by 28 CsI scintillators). GRAPE was flown from Fort Sumner, New
Mexico on the morning of September 26${}^{\text{th}}$ 2014\. The experiment
was at float altitude for 14.4 hours. The Crab (the pulsar and the nebula) was
the primary target of the experiment. A polarization measurement of the Crab
would not only validate the instrument design (for future missions observing
GRBs), but also improve our understanding of the emission mechanisms of the
Crab emission. The flight plan included 8 hours of Crab observation but the
flight was terminated before all data could be collected. The Crab was
observed for only 1.8 hours. Background dominates the observations at flight
altitudes, therefore a good background estimation is crucial for the analysis.
The background depends on many flight and instrument parameters (altitude,
instrument pointing, temperatures, etc). Estimating the Crab background using
these parameters were a challenging task. A technique based on the Principle
Component Analysis (PCA) was implemented which used the varying parameters to
estimate the background for the Crab. Our instrument design allows for an
optical crosstalk (a known issue) that arises from light exiting the
scintillator element and spilling to neighboring anodes. A model was developed
to represent the crosstalk, which was subsequently incorporated in our
simulations and the instrument response which further improved our analysis.
The analysis was focused on the phase-integrated data (due to limited
statistics of the shortened observation period). A power-law spectrum with
photon index of 1.70$\pm$0.24 and a normalization of 1.01$\pm$1.35
ph/keV/s/cm2 was measured. A polarization fraction of 0.43$\pm$0.4 and a
polarization angle of 56${}^{\circ}\pm$30∘ was measured for the phase
integrated Crab observation in the 70-200 keV energy range. This result was
not sufficiently significant enough to further our understanding of the
emission mechanism but was statistically consistent with other experiments.