# Quantum symmetrization transition in superconducting sulfur hydride from quantum Monte Carlo and path integral molecular dynamics

Romain Taureau, Marco Cherubini, Tommaso Morresi, Michele Casula

► To cite this version:

Romain Taureau, Marco Cherubini, Tommaso Morresi, Michele Casula. Quantum symmetrization transition in superconducting sulfur hydride from quantum Monte Carlo and path integral molecular dynamics. npj Computational Materials, 2024, 10 (1), pp.56. 10.1038/s41524-024-01239-0. hal-04776561

HAL Id: hal-04776561

https://hal.science/hal-04776561v1

Submitted on 11 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Quantum symmetrization transition in superconducting sulfur hydride from quantum Monte Carlo and path integral molecular dynamics

<figure><img src="image_3.png" /></figure>

Check for updates

Romain Taureau<fnref n="1" />, Marco Cherubini,<fnref n="1" />, Tommaso Morresi ② & Michele Casula ①✉

We study the structural phase transition, originally associated with the highest superconducting critical temperature $T_c$ measured in high-pressure sulfur hydride. A quantitative description of its pressure dependence has been elusive for any ab initio theory attempted so far, raising questions on the actual mechanism leading to the maximum of $T_c$. Here, we estimate the critical pressure of the hydrogen bond symmetrization in the $Im\overline{3}m$ structure, by combining density functional theory and quantum Monte Carlo simulations for electrons with path integral molecular dynamics for quantum nuclei. We find that the $T_c$ maximum corresponds to pressures where local dipole moments dynamically form on the hydrogen sites, as precursors of the ferroelectric $Im\overline{3}m$-$R\overline{3}m$ transition, happening at lower pressures. For comparison, we also apply the self-consistent harmonic approximation, whose ferroelectric critical pressure lies in between the ferroelectric transition estimated by path integral molecular dynamics and the local dipole formation. Nuclear quantum effects play a major role in a significant reduction (≈50 GPa) of the classical ferroelectric transition pressure at 200 K and in a large isotope shift (≈25 GPa) upon hydrogen-to-deuterium substitution of the local dipole formation pressure, in agreement with the corresponding change in the $T_c$ maximum location.

Since its discovery in 1911¹, superconductivity has been one of the most investigated topics in both theoretical and experimental physics. While it was discovered that almost every conductor could reach zero resistance at low-enough temperatures ($T < 10$ K)<fnref n="2" />, the quest for higher critical temperature ($T_c$) superconductors became the new challenge. Until recently, cuprates were leading the race with a $T_c$ as large as 133 K for Hg-Ba-Ca-Cu-O systems³, although the pairing mechanism in these materials is considered unconventional and it is not explained by the standard Bardeen-Cooper-Schrieffer (BCS) theory<fnref n="2" />.

In 2015, the discovery of conventional superconductivity in $H_3S$ with a maximum $T_c$ of 203 K reached at a pressure $P_c$ as high as 150 GPa⁵ paved the way to a new era of high-$T_c$ materials. Indeed, hydrogen (H)-based systems are nowadays the most promising candidates to achieve room-temperature superconductivity. As a matter of fact, in 2019, the same team that

discovered $H_3S$ claimed to have measured an even higher $T_c$ in $LaH_{10}$, superconducting already at 250 K⁸, later followed by a similar discovery in the yttrium hydride⁷. In a rush toward room-temperature superconductivity, more recent claims of $T_c$ larger than the one found in $LaH_{10}$ did not meet the consensus of the whole community⁸,⁹. The main issue of these materials is the extreme pressure conditions, usually larger than 150 GPa, needed to obtain the high-$T_c$ superconducting phase. Indeed, while all the binary candidates involving hydrogen were theoretically investigated, none of them seems to sufficiently decrease the pressure of the superconducting state. Eyes are now turned toward ternary materials¹⁰.

In this work, we focus on the prototypical case of $H_3S$ and we study its structural phase transition generally associated with the maximum of the superconducting critical temperature, located at around 150 GPa¹¹⁻¹⁴. According to x-ray diffraction data¹⁵, at lower pressures the sulfur (S) sites

<figure><img src="image_4.png" /></figure>

<footnote n="1">¹Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie, Sorbonne Université, CNRS UMR 7590, MNHN, 4 Place Jussieu, 75252 Paris, France.</footnote>

<footnote n="2">²European Center for Theoretical Studies in Nuclear Physics and Related Areas—Fondazione Bruno Kessler*, 286 Strada delle Tabarelle, 38123 Trento, Italy.</footnote>

<footnote>e-mail: michele.casula@sorbonne-universite.fr</footnote>
are arranged in a geometry that is compatible with the trigonal R3m symmetry (Fig. 1b) and, upon compression, the system undergoes a phase transition toward a body-centered-cubic (bcc) Im3̄m structure (Fig. 1a).

After the first theoretical prediction of high-$T_c$ superconductivity in $H_3S^{16}$, several works tried to explain the origin of the maximum of $T_c$ found in experiments as a function of pressure. Even if the magnitude of the calculated $T_c$ is right, confirming the BCS origin of the superconducting state, a quantitative disagreement between various theoretical approaches was found, with estimated $T_c$ values fluctuating over a 50 K range for the high-pressure phase¹⁷. Moreover, theoretical studies more oriented to understand the underlying structural properties of $H_3S$, revealed a significant disagreement in the transition pressures between the predicted phases. In those works¹⁸,¹⁹, the structural phase transition is explained by a quantum proton symmetrization from the R3m phase, with displaced protons, to the Im3̄m one, where every hydrogen lies in the midpoint of the two neighboring sulfur atoms (S-S midpoint). This is also called ferroelectric transition, because the hydrogen atoms displaced from the S-S midpoint lead in the R3m phase to a long-range order of local dipole moments, created by the H-S bond asymmetry. In that context, the shuttling mode of hydrogen atoms, namely their vibrational mode along the direction linking two neighboring sulfur atoms, was thoroughly investigated. The phase transition was then identified by looking at the dynamical instability of the symmetric Im3̄m phase when the pressure is lowered and the shuttling mode softens. On general grounds, this reflects the sudden transformation of the free energy profile, leading to a sign change of its curvature across the transition between two different crystal structures, one with lower symmetry than the other.

These findings were obtained by solving the nuclear Hamiltonian within the Stochastic Self-Consistent Harmonic Approximation (SSCHA)²⁰⁻²², which has proven to be one of the best approximated theories to deal with nuclear quantum effects (NQE). Within this framework, the electronic part was solved by Density Functional Theory (DFT) using different parametrizations for the exchange-correlation functional, like the Perdew-Burke-Ernzerhof (PBE)²³ and the Becke-Lee-Yang-Parr (BLYP)²⁴,²⁵ ones. Independently of the DFT functional used, a sizable underestimation of the experimental critical pressure $P_c$ by ≈40 GPa was always observed, leaving open the question about the origin of this mismatch, and whether this should be attributed to the electronic or to the nuclear components.

Here, we go beyond the previous state-of-the-art calculations by treating the electronic problem not only at the DFT-BLYP, but also at Quantum Monte Carlo (QMC) level, which provides a benchmark for the DFT methods. QMC is known to yield very accurate total energies in both molecules and solids²⁶⁻²⁸, thanks to its stochastic Green's function algorithms²⁹,³⁰, such as the lattice regularized diffusion Monte Carlo³¹, projecting any initial trial wavefunction toward the ground state of the system within the fixed node approximation. Moreover, we solve the nuclear

Hamiltonian by using Path Integral Molecular Dynamics (PIMD), which is in principle exact, outperforming any other approximation for the nuclear degrees of freedom. Then, we analyze the resulting phase diagram by looking at the ferroelectric order parameter, at the hydrogen/deuterium density, focusing on its transformation from the unimodal to bimodal distribution, and finally at its quantum fluctuations, detecting when the associated local polarization freezes in a displaced geometry.

In this work, we have been able to track the evolution of the mode distribution with a high resolution in volume (and pressure), thanks to a three-dimensional (3D) model of the shuttling mode. The reliability of our model has been benchmarked using ab initio PIMD simulations with BLYP electrons, across the local moment formation. The advantage of the 3D model is that its potential energy surface (PES) can still be derived by much more expensive, although more accurate, QMC calculations, allowing us to check the impact of the electronic description on the occurrence of a local polarization.

In the model we developed, all hydrogen atoms in the system are allowed to move in the same way. However, only the spatial degrees of freedom of a single H site are retained. This feature induces some limitations, such as the lack of spatially disordered H configurations, and of correlations beyond a single-site description. In spite of this, we can accurately describe the local path from the symmetric proton arrangement to the asymmetric one, by detecting the local moment formation in the system, related to the shuttling mode softening. We have finally performed both SSCHA and PIMD simulations of the 3D model to investigate how NQE treated at different levels of approximation affect the final outcome.

## Results

### Harmonic and anharmonic phonons

At high pressure, above 150 GPa, the $H_3S$ crystal is expected to be in the cubic Im3̄m symmetric phase (Fig. 1a), where every hydrogen atom sits on the midpoint of two neighboring sulfur atoms. Upon pressure release, the lattice undergoes a trigonal distortion and the hydrogen atoms leave the aforementioned midpoint to move closer to one of the two flanking sulfur atoms, leading to the R3m asymmetric phase, depicted in Fig. 1b. In our description, we introduce a simplification by neglecting the trigonal distortion, which is however very weak (<0.06°)¹⁹. Thus, the R3m phase considered here differs from the Im3̄m one just by the hydrogen positions.

In Fig. 2, we report the analysis of the phonon dispersion for different volumes of the cubic Im3̄m unit cell, obtained at the ab initio level using the BLYP functional, either through the harmonic approximation via Density Functional Perturbation Theory (DFPT), or with the inclusion of quantum anharmonicity via PIMD simulations at the temperature $T = 200$ K. Hereafter, volumes and energies will be expressed per $H_3S$ unit, while the unit cell will be taken as cubic with S atoms arranged in a bcc lattice.

<figure><img src="image_5.png" /><figcaption>**Fig. 1 |** $H_3S$ crystal structures. **a** Crystal structure of the Im3̄m symmetric phase (smaller volume, higher-pressure phase). **b** The R3m asymmetric phase (larger volume, lower-pressure phase). $d_{SS}$ is the lattice parameter of the bcc crystal.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 2 | Phonon dispersion for different volumes of the cubic $Im\bar{3}m$ unit cell, whose electronic structure has been computed with the BLYP functional.** At $V = 83.6$ $a_0^3$ (volume per H$_3$S unit), only the harmonic dispersion is reported. For all the other volumes, we compare the harmonic dispersion (light-blue color) with the</figcaption></figure>

PIMD one (green color). Full dispersions are obtained by interpolating harmonic (anharmonic) dynamical matrices defined on a $6 \times 6 \times 6$ ($2 \times 2 \times 2$) q-grid. PIMD simulations are performed at $T = 200$ K.

<figure><img src="image_7.png" /><figcaption>**Fig. 3 | Shuttling mode frequencies as a function of volume by different approaches.** At the PIMD level, the phonon frequencies are computed using the S-S midpoint as reference position for the quantum displacement-displacement correlator$^{32}$. Within the SSCHA framework, phonons are computed using the centroid position obtained through the free energy minimization and the full self-energy dynamical correction is included$^{30}$. Imaginary phonon components are represented with negative values.</figcaption></figure>

At this point, it is important to underline that the DFPT and the PIMD phonons bear different information (see also section “Phonons”). The PIMD phonons are computed through the quantum displacement-displacement correlator recently developed in ref. 32. They describe the lowest vibrational excitations$^{33}$, that is the energy difference between the first excited state and the ground state of the nuclear Hamiltonian. This is the quantity normally measured by experimental probes, such as infrared or

Raman spectroscopies. Consequently, phonons computed in this way fully include anharmonic effects and are always positive definite, meaning that they cannot describe dynamical instabilities via the appearance of imaginary phonons. This is at variance with the harmonic case or with approximated theories devised to deal with NQE, such as the SSCHA$^{20}$, which instead provide information about the sign of the free energy curvature at the reference geometry.

While for $V = 83.6$ $a_0^3$ only the harmonic dispersion is reported in Fig. 2, for larger volumes we compare the PIMD phonons (green lines) obtained in a $2 \times 2 \times 2$ supercell with the harmonic ones (light-blue lines). For PIMD phonons, the spatial range of the force constant matrix is such that the $2 \times 2 \times 2$ supercell is large enough to allow for a q-interpolation of the phonon branches $\omega_m = \omega_m(q)$. The comparison between PIMD and harmonic phonons of Fig. 2 clearly shows how strong NQE are and how sizable is the softening of the most energetic phonons due to quantum anharmonicity, particularly at the largest volumes. In the harmonic framework, for $V > 85$ $a_0^3$ (see Figs. 2 and 3), the appearance of imaginary frequencies indicates the dynamical instability of the $Im\bar{3}m$ structure. More specifically, the softening of the shuttling hydrogen mode at $q = \Gamma$ signals the transition toward the asymmetric $R\bar{3}m$ phase$^{18}$. From Fig. 2, one can see that imaginary frequencies disappear in PIMD phonons and their evolution as a function of volume is much smoother than in the harmonic case. As expected from the definition of PIMD phonons, PIMD simulations never yield imaginary frequencies for the shuttling mode. In this regard, see also Fig. 3, where we report the shuttling mode frequency we obtained as a function of volume at different levels of theory. This analysis shows that the putative transition implied by the maximum in $T_c$ cannot be determined using solely the shuttling mode frequency as a proxy, because in the volume range corresponding to the experimental $T_c$ maximum$^5$, i.e., between $98$ $a_0^3$
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | 3D-PES landscapes (in meV per H₃S unit) for different volumes. Top panel: DFT-BLYP PES computed on the plane containing the S-S direction (x axis) and bisecting the y-z quadrant (y = z plane). The points on the plane are fully defined by their x and y coordinates, expressed in fractional units. Middle panel: 1D</figcaption></figure>

projection along the axis connecting the two minima. Bottom panel: same as Top panel with QMC PES. In the top and bottom panels the zero of energy is the PES value at (0.5,0,0). In the middle panel, the zero is the PES minimum.

and 100 $a_0^3$ (see equation of state (EOS) in Fig. 9a), there is no anomalous behavior of the proton shuttling mode frequency. We need to rely upon other observables in a framework describing nuclei as quantum particles.

So far, we have reported the structural behavior as a function of volume, fixed in our simulations. However, we can easily deduce the corresponding pressure by deriving the EOS P = P(V) using the Vinet relation³⁴ computed with the same functionals employed to calculate the phonon dispersions. Nevertheless, our goal is to go beyond DFT and reach a more accurate electronic description of the system using QMC methods (the details of our QMC calculations are reported in the “Electronic structure calculations for the PES model” section). A simple comparison of the EOS produced by the two approaches, shown in Fig. 9a, reveals visible differences, suggesting that a description of the electronic structure at the QMC level is crucial to estimate correctly the critical pressure. Unfortunately, QMC calculations are much more expensive than DFT, and coupling them with ab initio PIMD simulations to study the real crystalline system is out of reach. Therefore, we need a simplified PES describing the hydrogen shuttling mode that can be derived, after a fitting procedure, from QMC total energy calculations performed on a coarse grid of nuclear configurations. This model PES can then be used to compute the shuttling mode frequencies and to study the local polarization properties induced by the proton displacement at the PIMD level.

### Classical 3D model

The model PES is derived by considering the collective and coherent motion of all the hydrogen atoms along the direction connecting the two S atoms flanking each H (S-S direction), by allowing also hydrogen out-of-axis mobility, while the S atoms are pinned in their bcc positions. In this way, we aim at reproducing the shuttling mode dynamics that takes place at $q = \Gamma$, thus having the same modulation for all hydrogen atoms in the crystal. Therefore, we reduce the 3N dimensions of the ab initio potential (with N the number of atoms in the supercell) to only 3 dimensions. The PES is fitted over total energies generated either by DFT-BLYP or by QMC for nuclear configurations defined on a cylindrical grid. Further details about the

model description can be found in the “Potential energy surface parametrization” section.

In Fig. 4, we report the PES profiles obtained by solving the electronic problem within the DFT-BLYP (first row) and QMC (third row) methods. At the volumes taken into account here, both DFT-BLYP PES and QMC PES have two minima connected through the inversion symmetry with respect to the S-S midpoint ((0.5,0,0) in fractional units). The second row shows a comparison of both energy profiles cut along the line connecting these two PES minima, going through the S-S midpoint. For the smallest volume analyzed, $V = 90.9 a_0^3$, we found a good agreement between the DFT-BLYP PES and QMC PES, suggesting that electron correlation effects are reasonably well described at the DFT level at high-enough pressures. However, the discrepancy between the two approaches appears when we increase the volume and it grows continuously upon pressure release. For the largest volume considered, $V = 110.0 a_0^3$, the height of the double well barrier for QMC is ~270 meV per H₃S unit, 80% larger than the DFT-BLYP one.

The ferroelectric transition volume for classical nuclei can be estimated based on the PES by using the Landau theory for continuous phase transitions³⁵. This method relies on the sign change of the free energy curvature (the total energy curvature at $T = 0$ K) at the volume when the two displaced minima merge into a single one, in the symmetric configuration corresponding to the point (0.5,0,0) in Fig. 4. For DFT-BLYP, we found a critical volume $V_{ferro}$ around $85 a_0^3$ corresponding to a pressure of 263 GPa, while for QMC we found the same volume ($\approx 85 a_0^3$) which corresponds to 238 GPa in this case (see Fig. 9a). We note that the $V_{ferro}$ yielded by the BLYP 3D-PES is in nice agreement with the value at which the shuttling mode frequency vanishes, computed ab initio in the harmonic approximation (see Fig. 3). This is a signature of our model PES quality.

### Quantum 3D model

In order to have a reliable description of the structural phase transition based on our 3D-PES, we need to include nuclear quantum effects. We add them by performing PIMD calculations at $T = 200$ K as implemented in ref. 36. Numerical details of these simulations can be found in the “PIMD simulations” section. Here, we mention only that in the PIMD simulations of our
<figure><img src="image_9.png" /><figcaption>**Fig. 5** | **Quantum hydrogen distributions.** a Two-dimensional projection on the $y=z$ plane of hydrogen positions with our 3D model solved by PIMD at 200 K. The results of BLYP PES calculations are represented in the top panel and the QMC PES ones in the bottom. b Hydrogen distribution projected along the shuttling mode.</figcaption></figure>

Blue points are for the BLYP 3D-PES; green points for the ab initio BLYP simulations; red points for QMC 3D-PES. Lengths are expressed in Bohr. The origin of the reference frame is centered at the S-S midpoint. $|\delta x|$ is the distance from the origin.

3D model the hydrogen atom has an effective mass equal to three times the physical hydrogen mass, owing to the fact that the PES is expressed per H$_3$S unit and the 3D motion of all hydrogen atoms in the H$_3$S molecule is concerted by construction.

In Fig. 5, we report the projections of the resulting 3D proton density, which takes two distinct shapes depending on the volume. The density exhibits only one peak centered in the middle of the S-S axis for small volumes and the central peak splits into two lobes for the largest volumes. In the contour plot of Fig. 5a, one can clearly see that the doubling of the peak happens in QMC at smaller volumes (higher pressures) than in DFT-BLYP, as expected from the analysis of the classical PES, which shows deeper minima in the QMC PES at fixed volume. In Fig. 5b, we plot the distribution of the hydrogen position projected along the shuttling mode direction, by including also the data coming from the ab initio PIMD simulations driven by DFT-BLYP forces. Our 3D model has a similar behavior in comparison with the full 3N dimensional system. The mode distribution assumes a double peak shape at approximately the same volume for the model (blue lines) and the ab initio system (green lines), evaluated for the same BLYP functional. The main differences are the broadness of the distribution,

<table><caption>**Table 1** | Transition pressures and volumes for the local moment formation yielded by PIMD</caption><thead><tr><td>Theory</td><td>Isotope</td><td>Probe</td><td>V<sub>c</sub> [a<sub>0</sub><sup>3</sup>]</td><td>P<sub>c</sub> [GPa]</td></tr></thead><tbody><tr><td>DFT-BLYP</td><td>H<sub>3</sub>S</td><td>Fluctuations</td><td>104.2</td><td>130</td></tr><tr><td>DFT-BLYP</td><td>H<sub>3</sub>S</td><td>Density</td><td>99.6</td><td>153</td></tr><tr><td>DFT-BLYP</td><td>D<sub>3</sub>S</td><td>Fluctuations</td><td>97.6</td><td>165</td></tr><tr><td>DFT-BLYP</td><td>D<sub>3</sub>S</td><td>Density</td><td>96.6</td><td>171</td></tr><tr><td>QMC</td><td>H<sub>3</sub>S</td><td>Fluctuations</td><td>100.9</td><td>126</td></tr><tr><td>QMC</td><td>H<sub>3</sub>S</td><td>Density</td><td>96.3</td><td>152</td></tr><tr><td>QMC</td><td>D<sub>3</sub>S</td><td>Fluctuations</td><td>95.6</td><td>156</td></tr><tr><td>QMC</td><td>D<sub>3</sub>S</td><td>Density</td><td>93.6</td><td>169</td></tr></tbody></table>

underestimated by the model, and the position of the peak, which lies closer to the S atoms in the ab initio simulations. These differences can be understood based on the enhanced quantum-thermal fluctuations of the ab initio system compared to the one with a reduced number of degrees of freedom. Nevertheless, as far as the peak splitting is concerned, the ab initio and the model PIMD calculations are in agreement. This validates the accuracy of our 3D-PES model, which then allows one to compare directly BLYP and QMC results. The projected 1D distribution in Fig. 5b reveals that the QMC PES leads to a smaller volume for the peak splitting, as shown already in the contour plot of Fig. 5a.

By fitting the distribution in Fig. 5b and interpolating the parameters obtained for several volumes, it is possible to determine precisely the position of its maximum as a function of volume, and thus the occurrence of the bimodal distribution. We estimate the peak splitting to take place in H$_3$S at a volume of 99.6 a$_0^3$ for DFT-BLYP, and of 96.3 a$_0^3$ for QMC. According to the EOS of Fig. 9a, these volumes correspond to pressures of 153 GPa for DFT-BLYP and of 152 GPa for QMC. These values, reported in Table 1, are in good agreement with the position of the maximum $T_c$ measured in experiments$^5$, and they are strongly affected by NQE. Nevertheless, it is important to underline that the two similar pressures obtained by BLYP PES and QMC PES after inclusion of NQE originate from a compensation of errors in the BLYP values, if we take QMC as reference. Indeed, a volume overestimation found in the BLYP PES compensates with a pressure overestimation in the BLYP EOS to yield approximately the same pressure for the peak splitting as the one found in QMC (see Fig. 9a).

The occurrence of the bimodal distribution signals the proximity of a critical region where the quantum proton is more localized, either dynamically or statically, in one of the two wells. However, from a more rigorous point of view the instantaneous localization of the proton, namely the local dipole formation, is better characterized by the “local moment” susceptibility, as defined below. Indeed, quantum fluctuations are at work across the transition to make hydrogen atoms shuttle between the two PES minima. As
<figure><img src="image_10.png" /><figcaption>**Fig. 6 |** Imaginary time correlator $g(β/2)$ and normalized local moment susceptibility. Local moment susceptibility for H$_3$S (a) and D$_3$S (b) for quantum nuclei with the DFT-BLYP (blue) and QMC (red) 3D-PES. The vertical dashed lines indicate the susceptibility maxima. c Imaginary time correlator $g(β/2)$ for H$_3$S. Blue</figcaption></figure>

circles and green stars are DFT-BLYP results for the 3D-PES and ab initio simulations, respectively. For reference, vertical dashed lines indicate the position of the local moment susceptibility peak and the ferroelectric transition, as determined in Fig. 7.

the volume increases and the minima deepen, the fluctuations will start freezing, leading to the creation of a local electric dipole moment, generated by the statically displaced proton in the R3m phase, or generated dynamically, by instantaneous configurations where the whole path representing the quantum proton is fully localized in one of the two wells. PIMD fully accounts for quantum fluctuations, thanks to its imaginary time resolution. We can measure them by computing the imaginary time correlator $g(β/2) = ⟨δx(0)δx(β/2)⟩$, with $β = 1/(k_B T)$ the inverse temperature used in the PIMD simulations, and $δx(τ) = x(τ) - ⟨x⟩$, where ⟨x⟩ is the thermal quantum average of the $x$ coordinate, corresponding to the symmetric position in the 3D model. Quantum fluctuations reduce the value of $g(β/2)$. A non-zero value of $g(β/2)$ can be interpreted by the presence of a finite moment in the distribution. In our 3D-PES, this moment is by definition local, because by model construction the hydrogen dynamics is condensed in a single 3D site. Therefore, the local moment susceptibility $χ_g$ is the normalized variance of $g(β/2)$, namely $χ_g = Var[g(β/2)]/Var[g(0)]$. Within the local moment fluctuation picture, the occurrence of local polarization can then be estimated by evaluating the volume at which $χ_g$ is maximum, as shown in Fig. 6a, b. This quantity has already been used in a previous work³⁷ to identify the transition from a paraelectric phase to a disordered regime in an anharmonic oscillator chain, characterized by a tunable double well potential, where the symmetry is locally and instantaneously broken in favor of displaced configurations. If we describe the change of regime based on local moment fluctuations, we obtain $V_c = 104.2 a_0^3$ for DFT-BLYP, corresponding to $P_c = 131$ GPa, and $V_c = 100.9 a_0^3$ for QMC, corresponding to $P_c = 126$ GPa (Fig. 6a and Table 1). Also in this case, like for the density probe, cancellation of errors is at play and, by consequence, the two electronic descriptions provide almost the same critical pressure.

We evaluated $g(β/2)$ also in our ab initio PIMD simulations, and in Fig. 6c we compare it against the values of $g(β/2)$ coming from the PIMD solution of the 3D model. The ab initio and model results are in statistical agreement for this local quantity, confirming that the 3D model correctly captures the volume evolution of the local polarization.

The two probes we used in this work, the local moment susceptibility and the peak splitting, allow us to determine a lower and an upper bound for the pressure where the fluctuating local dipoles disappear in favor of a paraelectric phase, by squeezing the compound. Notice that this does not correspond to the ferroelectric transition pressure, associated instead with the global Im3m-R3m symmetry breaking, and long-range dipole order, which happens at lower values. The same analysis is carried out for both the hydrogen H$_3$S and deuterium D$_3$S crystals, to estimate the magnitude of isotope effects (see Fig. 6b). We summarize the results in Table 1, where we show that the hydrogen-to-deuterium substitution brings about an increase of the local polarization formation pressure that falls into the [17–35] GPa range.

## Full BLYP-PIMD solution of the H$_3$S phase diagram at 200 K and comparison with SSCHA

After having analyzed the local moment formation with the help of the 3D model, we turn now the attention to the ferroelectric transition, associated with the global R3m → Im3m transformation. The suitable order parameter to identify this transition is $Δ = ⟨1/N ∑_{i=1}^N δx_i⟩$, where the sum runs over all the N hydrogen atoms in the supercell, and $δx_i$ is the distance of the ith proton from the S-S midpoint at a given snapshot. An equivalent order parameter, showing usually less statistical fluctuations, is $Δ_{abs} = ⟨|1/N ∑_{i=1}^N δx_i|⟩$. The brackets indicate the average over the classical or quantum nuclear distribution. In PIMD simulations an additional average is then done over the beads positions. From these definitions, it is clear that this order parameter can only be computed in our ab initio simulations, being the 3D model local.

In Fig. 7, we plot the volume dependence of the order parameter $Δ$ and $Δ_{abs}$ and their susceptibilities. Their peak is located at the ferroelectric transition, occurring at $V_{ferro} = 117a_0^3$, $P_{ferro} \approx 82$ GPa, as found in our BLYP-PIMD simulations in the $2 \times 2 \times 2$ supercell at $T = 200$ K. We notice that the peak location is correlated with the jump in the shuttling mode frequency, reported in both Figs. 3 and 7. The agreement between the $Δ$ and $Δ_{abs}$ susceptibilities and the shuttling mode frequency jump strengthens the reliability of our estimate. Interestingly, at 200 K the ferroelectric transition takes place at a pressure much lower than the one where the local moments are suppressed. For quantum nuclei, between these two pressures the system is in a regime characterized by disordered local moments and Im3m symmetry (see Fig. 8 for the resulting phase diagram). Accounting for thermal and quantum effects leads to a strong reduction of the critical ferroelectric pressure observed in the classical framework, which is as large as 263 GPa at zero temperature. Furthermore, in order to distinguish between anharmonicity coming from thermal and NQE, we also performed classical ab initio MD simulations at the same temperature (see Supplementary Note V of the Supplementary Information (SI)), which yield a ferroelectric transition at ≈133 GPa (see Fig. 8). Thus, classical anharmonicity accounts for about 70% of the total pressure reduction of the ferroelectric transition at 200 K. The remaining 30% is due to NQE at the same temperature.

The ferroelectric transition cannot be estimated at the QMC level, due to its computational cost when applied to the dynamics of a real extended system. However, as we have seen, the 3D model, derived at both the DFT-BLYP and QMC levels, is enough to determine the formation of local electric dipole moments, whose pressure $P_c$ matches well the position of the experimental $T_c$ maximum.

At variance with previous state-of-the-art calculations based on a combination of DFT-BLYP and SSCHA frameworks, our PIMD results yield $P_c$ in a substantial agreement with the experimental finding, and this irrespective of the electronic theory used to generate the PES. One should notice here that the original SSCHA approximation is not able to capture the
<figure><img src="image_11.png" /><figcaption>**Fig. 7 | Ferroelectric transition from ab initio BLYP PIMD simulations at 200 K in a 2 × 2 × 2 supercell.** Black (blue) points are values of the ferroelectric order parameter $\Delta (\Delta_{\text{abs}})$ as a function of volume. The variance (i.e., susceptibility) of these order parameters is represented by green and red crosses, respectively. We also report the shuttling mode frequencies (gold pentagons). All quantities are in atomic units except for the phonon frequencies, expressed in cm$^{-1}$. The vertical red-dashed line locates $V_{\text{ferro}}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8 | Phase diagram at 200 K for classical and quantum nuclei as determined by BLYP-driven classical and path integral molecular dynamics.** Both ferroelectric transition (black dots) and the region where the local polarization vanishes into a paraelectric Im3m phase (bracketed by the blue star and the purple pentagon) are reported. Between ferro and para, an Im3m phase with disordered local moments is expected to take place$^{37}$. The disordered region is reported only for the quantum case. We also report the ferroelectric critical pressures predicted by SSCHA, from both ab initio$^{19}$ (white square) and our 3D-PES model (orange diamonds). For comparison, the pressure of the experimental $T_c$ maximum is shown by a teal diamond.</figcaption></figure>

disordered Im3m phase, being a mean-field theory with no disordered configurational entropy and no direct information of imaginary time correlations, key to detect the dynamical local moment formation. Thus, the critical pressure that SSCHA can normally compute is the ferroelectric one, $P_{\text{ferro}}$, and not $P_c$. However, time resolved extensions of SSCHA have been recently proposed, able in principle to access also retardation effects$^{38}$.

To investigate more deeply this mismatch, we carry out SSCHA calculations with our model PES (see section “SSCHA simulations” for details). In SSCHA, the occurrence of the asymmetric R3m phase is signaled by a centroid displaced with respect to the S-S midpoint. The SSCHA critical values are $V_{\text{ferro}} = 107.8a_0^3$, $P_{\text{ferro}} = 114$ GPa for the DFT-BLYP PES, and $V_{\text{ferro}} = 102.4a_0^3$, $P_{\text{ferro}} = 118$ GPa for the QMC PES. As in PIMD, there is no significant difference in the transition pressures between the electronic structure methods used to generate the PES. Our SSCHA results for the DFT-BLYP PES are in a very good agreement with the outcome of previous

SSCHA simulations for the full ab initio system$^{19}$, calculated with the same DFT-BLYP functional. We find that the SSCHA overestimates $P_{\text{ferro}}$ with respect to the one obtained in PIMD for the same BLYP functional, as expected from a mean-field theory, and underestimates $P_c$, which is however out of reach by the SSCHA, suggesting that the approximated description of the nuclear Hamiltonian is the source of disagreement with both PIMD and experimental results.

Let us look now at the predictions for the shuttling mode frequencies, plotted in Fig. 3 for various methods. It has to be noted that, within the SSCHA, the phonon frequency of the shuttling mode shows a jump at $V_{\text{ferro}}$. This is due to the hop of the SSCHA centroid from its symmetric position to a different minimum of the free energy, already “preformed”, which breaks the symmetry and becomes energetically more favorable at $V_{\text{ferro}}$. Moreover, we also observe an increase of the SSCHA phonon line-width across the transition of the order of $10 \text{ cm}^{-1}$. A similar jump in the shuttling phonon frequencies is detected by our ab initio BLYP simulations in correspondence with the ferroelectric transition (see Figs. 3 and 7).

Nevertheless, our PIMD phonon determination shows a progressive phonon softening without jumps across the volume region of local moment formation. This is not only true within our 3D model PES, but also for our PIMD calculations driven by ab initio forces computed at the DFT-BLYP level, as shown in Fig. 3. The agreement between shuttling mode frequencies yielded by the 3D model and the ones given by ab initio calculations in this volume region highlights once again the quality of our model PES. This supports the hypothesis of two different transitions. The first one is a smooth transition, or crossover, from the paraelectric Im3m to a phase sharing the same Im3m symmetry and characterized by the formation of local and spatially disordered local moments. This phase cannot be detected by looking at the phonon frequencies, and it is not accessible within the SSCHA formulation. The second one is the ferroelectric transition from the disordered Im3m to the asymmetric R3m phase, which happens at significantly lower pressure than the first one, where the shuttling phonon frequency shows a jump. The phase diagram deducible from our combined ab initio and 3D model results is drawn in Fig. 8.

## Discussion

In this work, starting from ab initio electronic structure calculations, we generated a model PES to describe the shuttling mode of hydrogen in H3S, responsible for the R3m → Im3m transition, which was originally associated with the $T_c$ maximum as a function of pressure. Despite the fact that such a hydrogen symmetrization is expected to happen in H3S upon compression, so far no theoretical method has been able to spot it at pressures near the one that maximizes $T_c$ in experiments. This raised doubts on the original association between superconductivity and structural transition$^{39,40}$, worsened by the fact that other competing symmetries could be stable in the same pressure range$^{41–44}$. The mismatch found between previous theoretical estimates of the critical pressure $P_c$ and the experimental values for the $T_c$ maximum is solved by applying state-of-the-art computational methods in both the electronic and nuclear Hamiltonians, namely using QMC calculations for electrons, and the PIMD approach for nuclei. Within our QMC+PIMD approach, the experimental pressure where $T_c$ is maximum is bracketed by the $P_c$ value estimated from the local fluctuations probe and the one determined by the transformation of the bimodal hydrogen distribution into a unimodal one. Consequently, these two probes provide a lower and an upper bound for the critical pressure, with a range between the two of ≈20 GPa. The range of transition pressures identified is consistent with the available experimental data for the $T_c$ maximum$^{5,11,13}$ for both H3S and D3S, as we can see in Fig. 9b.

We have thus shown that the occurrence of the $T_c$ maximum should be linked with the formation of the phase characterized by disordered local moments$^{37}$, and it cannot be associated with the ferroelectric R3m → Im3m transformation, which takes place at a lower pressure $P_{\text{ferro}}$ compared with $P_c$. According to our outcome, $T_c$ reaches its maximum when the local dipole moments melt upon compression, and protons become fully delocalized across the PES barrier.
**Fig. 9 | Transition values from theory and experiment.** a Equation of state $P = P(V)$ (see section “Equation of state” and Eq. (7)) for both DFT-BLYP (blue color) and QMC (red color) calculations of H$_3$S in the Im$\overline{3}$m phase. The transition volumes and the corresponding pressures identified at different levels of theory using the density probe (peak splitting) are displayed by dashed lines. **b** $T_c$ as a function of pressure in H$_3$S and D$_3$S from Drozdov et al.⁵, Einaga et al.¹¹ and from Minkov et al.¹³. The shaded areas on each data set represent the pressure range when the transition occurs according to our PIMD results for the QMC PES, where the lower limit is based on the local quantum fluctuations analysis and the upper one on the density evolution (see Table 1 for the actual values).

<figure><img src="image_13.png" /></figure>

Furthermore, we notice that the ab initio electronic structure computed at the DFT-BLYP level predicts very good results for the critical pressure, similar to those obtained by QMC. However, it is important to stress that the DFT-BLYP pressures are affected by error compensation, the overestimation of the critical volume being balanced by a different EOS if compared against QMC calculations. This aspect underlines the importance of using an accurate electronic description, beyond the DFT level. The generation of our model PES, built to describe the hydrogen shuttling mode, allowed us to exploit the QMC energies in a PIMD framework, otherwise unfeasible in the full 3N dimensional system.

We conclude by noting that the R$\overline{3}$m → Im$\overline{3}$m structural phase transition in sulfur hydride has strong analogies with the hydrogen bond symmetrization in other compounds such as high-pressure ice, where, upon compression, phase VII and VIII hosting displaced protons, stable at lower pressure, are expected to transform into the symmetric phase X$^{45,46}$. However, it is still a matter of debate whether the transformation is direct or whether other intermediate disordered structures appear, with protons only partially symmetrized. In this respect, further work is needed to extend our model beyond the collective path dynamics to treat non-local spatial correlations and disordered patterns. Machine learning schemes could then be useful to generate more extended PES from QMC data$^{47-49}$ with the aim at including a larger variety of hydrogen configurations in PIMD calculations by keeping the same QMC accuracy.

## Methods

### Electronic structure calculations for the PES model

For the DFT electronic structure calculations, we used the Quantum Espresso (QE) suite of codes$^{50,51}$, while for the QMC calculations, we employed the TurboRVB package$^{52}$. For sake of consistency, in both DFT and QMC calculations, we used the same set of pseudopotentials. Namely, we treated the sulfur atom with the ccECP neon-core pseudopotential$^{53}$ particularly suited for correlated calculations, available in both the QE-compatible Unified Pseudopotential Format (UPF) and in the TurboRVB-compatible Gaussian expansion format. For hydrogen, we used the bare Coulomb potential, with a very short-range cutoff for a QE usage within the plane-wave framework. In the QMC calculations instead, no short-range cutoff is needed for the bare Coulomb potential, because the nuclear cusp conditions are automatically fulfilled by our QMC wave function (see below). These pseudopotentials have been chosen after performing preliminary calculations at the DFT level to test their accuracy. We also tested other pseudopotentials (ultrasoft (US), projector augmented wave (PAW), and a combination of the above), by comparing the total energy profile obtained by moving the hydrogen atom away from the S-S midpoint, and constrained to stay on the S-S axis. This leads to a very crude one-dimensional (1D) PES, which is however useful for testing purposes, with the advantage that it is easily computable for its simplicity. We took as reference the total DFT energy computed with the all-electron LAPW

approach, as implemented in Elk$^{54}$. The ccECP pseudopotential for the sulfur atom and the bare Coulomb potential with short-range cutoff for the hydrogen atom turned out to be the most accurate choice (see Supplementary Note I).

For single-point calculations at selected nuclear configurations, we carried out DFT calculations with the Becke-Lee-Yang-Parr (BLYP) functional$^{24,25}$. The cutoff energy for plane waves is set to 200 Ry (due to the hardness of the H Coulomb pseudopotential), with the smearing parameter equal to 0.002 Ry and a k-points grid of 32 × 32 × 32.

For the QMC calculations, we used a Slater-Jastrow wavefunction $\Psi$, which reads as:

$$ \Psi = \Phi_S \cdot \exp(J), \quad (1) $$

where the term $\exp(J)$ is the Jastrow factor, symmetric under electron exchange, while $\Phi_S$ is the antisymmetric Slater determinant. The Slater orbitals in $\Phi_S$ are generated by DFT calculations within the Local Density Approximation (LDA)$^{55}$, performed in a Gaussian basis set by means of the DFT code built in TurboRVB. For the sulfur atom, we employed a modified cc-pVTZ primitive basis set with 6s6p2d1f components, contracted into 11 hybrid orbitals through the Geminal Embedded Orbitals (GEO) procedure$^{56}$. For hydrogen, we used a modified cc-pVTZ primitive basis set with 4s2p1d components contracted into 6 GEO hybrid orbitals.

The Jastrow exponent $J$ introduces explicitly electronic correlation in the wavefunction, and it can be decomposed into three terms, such that $J = J_1 + J_2 + J_3$.

$J_1$ is the so-called one-body term, which takes into account the interaction effects between the electrons $i$ and a nucleus $I$, and it depends on the relative electron-nucleus distances $r_{il}$. $J_2$ is the so-called two-body term, treating the correlations between electrons $i$ and $j$, and depending on their relative distance $r_{ij}$. Both $J_1$ and $J_2$ are designed to fulfill the electron-nucleus and electron-electron cusp conditions, respectively. They read as $J_1 = \sum_{i=1}^{N_e} \sum_{l=1}^{N} u_i(r_{il})$, and $J_2 = \sum_{i,j=1}^{N_e} v(r_{ij})$, where $N_e$ is the number of nuclei (electrons) in the supercell, and the functions $u$ and $v$ are defined as follows:

$$ u_i(r) = \frac{Z_1}{a} (1 - e^{-ar}) \quad (2) $$

$$ v(r) = \frac{r}{2(1 + br)}, \quad (3) $$

with $a$ and $b$ variational parameters, and $Z_1$ the charge of the $I$th pseudotatom. The coefficients in Eqs. (2) and (3) are set to fulfill the Kato cusp conditions for electron-nucleus and electron-electron coalescence, respectively$^{57}$.
$J_3$ is the three-body term that accounts for the electron-electron-nucleus interactions. As defined in TurboRVB, it is also intrinsically non-homogeneous, because it depends on the individual electron positions and not only on the relative distances, which is less accurate. Being non-homogeneous, it is expanded on a modified atomic Gaussian basis set of $2s2p1d$ atomic orbitals, for both sulfur and hydrogen atoms.

The $J_3$ parameters, together with a and b, are optimized by minimizing the variational energy of the many-body wavefunction in Eq. (1). The Slater part is instead kept frozen as determined by DFT-LDA. As stochastic minimization algorithm, we employed the linear method$^{58}$. We then carried out lattice regularized diffusion Monte Carlo (LRDMC) calculations$^{31}$, to stochastically project the initial wavefunction toward the ground state of the system, within the fixed node approximation. Within this approximation, the LDA nodes provide accurate results for this system, as verified in Supplementary Note II. In LRDMC, we used a lattice space of 0.25 $a_0$, which is known to produce converged energy differences. We started the projection from the best variational state optimized in the previous step, taken as trial wavefunction. Finite-size scaling has been performed on the $2 \times 2 \times 1$, $2 \times 2 \times 2$, $3 \times 2 \times 2$ and $3 \times 3 \times 2$ real-space supercells in order to extrapolate the LRDMC total energy to the thermodynamic limit, by also using Kwee-Zhang-Krakauer (KZK)$^{59}$ corrections to make its size dependence milder.

This workflow has been repeated for every point in the real-space grid used to interpolate the PES model from ab initio data (see section “Potential energy surface parametrization”).

## Potential energy surface parametrization

To derive an effective low-dimensional PES, we considered the collective and concerted motion of all the hydrogen atoms of the cubic unit cell, with sulfur atoms forming a bcc sublattice. The position of a hydrogen atom is described by the cylindrical coordinates $x, r$ and $\phi$, defined along the axis connecting the two flanking sulfur atoms (S-S axis): $r$ is the position of the hydrogen atom along the S-S axis, $r$ is the radial distance from the S-S axis, and $\phi$ is the azimuthal angle, wrapping around the same axis. We use fractional coordinates, where the lengths are expressed in $d_{SS}$ units, $d_{SS}$ being the lattice parameter of the bcc unit cell. Within this reference system, the S-S midpoint has coordinates $(x, r, \phi) \equiv (0.5, 0, 0)$. We assume that all hydrogen atoms in the unit cell move in the same way. This fixes the choice of a collective path connecting the $Im\bar{3}m$ symmetry (with all hydrogen atoms sitting at the S-S midpoints) to the $R\bar{3}m$ one (with all hydrogen atoms coherently displaced from the midpoint). In this way, we apply a dimensionality reduction of the full potential, depending on $3N$ dimensional coordinates, where $N$ is the number of atoms in the cell, to a much simpler 3D PES: $E = E(x, r, \phi)$.

The functional form of our 3D PES is constructed as follows:

$$E(x, r, \phi) = A(x, r) + B(x, r) \sin(\phi + 5\pi/4), \quad (4)$$

with:

$$A(x, r) = \frac{f_{\max}(x, r) + f_{\min}(x, r)}{2}, \quad B(x, r) = \frac{f_{\max}(x, r) - f_{\min}(x, r)}{2}, \quad (5)$$

and where $f_{\min}$ and $f_{\max}$ are defined as:

$$\begin{aligned} f_{\min, \max}(x, r) = & a + \frac{1}{2}b(x - 0.5)^2 + \frac{1}{4}c(x - 0.5)^4 \\ & + dr + \frac{1}{2}er^2 \\ & \pm f(x - 0.5)r \pm g(x - 0.5)r^2 \\ & \pm h_1(x - 0.5)^3 r \pm h_2(x - 0.5)^3 r^2 \\ & + h_3(x - 0.5)^2 r + h_4(x - 0.5)^2 r^2 \\ & + h_5(x - 0.5)^4 r + h_6(x - 0.5)^4 r^2 \end{aligned} \quad (6)$$

The choice of this functional form is motivated by the symmetries of the system. For a fixed {$x, r$}, the potential $E$ has an angular dependence that

varies following a sine curve with $2\pi$-periodicity. In particular, for $x < 0.5$, $E$ has a minimum given by $f_{\min}$ at $\phi = \pi/4$ and the maximum $f_{\max}$ at $\phi = 5\pi/4$. This dependence is built in Eqs. (4) and (5). The {$f_i(x, r)$}$_{i=\min,\max}$ functions in Eq. (6) are a composition of the following terms: a Landau-type potential that well describes the energy profile for $r=0$, a second-order polynomial function in $r$ for $x=0.5$, and mixed terms made of cross products of factors up to the fourth order in $(x-0.5)$ and up to the second order in $r$, which give enough flexibility in order to well reproduce the total PES. The signs in {$f_i(x, r)$}$_{i=\min,\max}$ ensure the symmetry: $E(1-x, r, \phi+\pi) = E(x, r, \phi)$, fulfilled by the system.

We sampled the PES by discretizing the 3D space according to the following grid defined in cylindrical coordinates: $x = [0.42, 0.44, 0.46, 0.48, 0.5]$ (in $d_{SS}$ units), $r = [0.00, 0.02, 0.05, 0.08]$ (in $d_{SS}$ units) and $\phi = [\pi/4, 5\pi/4]$. For these points we computed the ab initio total energies, given either by DFT-BLYP or by QMC calculations. We finally used the generated datasets to best fit the PES, parametrized according to Eqs. (4)-(6). The root mean square error of these fits amounts to $\approx 1$ meV per $H_3S$ unit.

## Equation of state

In order to get the pressure associated to each volume, $P = P(V)$, we use the Vinet EOS:

$$P(V) = 3B_0 \frac{1-\eta}{\eta^2} \exp\left(-\frac{3}{2}(B'_0-1)(1-\eta)\right). \quad (7)$$

with $\eta = (V/V_0)^{1/3}$. In Eq. (7), the parameters $V_0$, $B_0$ and $B'_0$ are the equilibrium volume, the isothermal bulk modulus, and the derivative of bulk modulus with respect to pressure, respectively. The Vinet EOS$^{34}$ is empirical and, despite having only a few parameters, it is very accurate to describe solids under extreme conditions. We obtained $V_0$, $B_0$ and $B'_0$ by fitting the $E=E(V)$ relation for the $Im\bar{3}m$ phase, where the total energy is computed from first principles, either by DFT-BLYP or by QMC, on a grid of volumes (see Fig. 9a). In the fit, we disregarded the Zero-Point Energy (ZPE) contribution, because we verified that the ZPE variation is very small (<1 mHa per $H_3S$ unit) in the range of pressures analyzed here within the same $Im\bar{3}m$ phase.

## PIMD simulations

The PIMD simulations are carried out at 200 K using 20 beads with ab initio DFT forces, while using 40 beads with 3D-PES forces, to take into account quantum effects. A convergence study of the PIMD results with respect to the number of beads is reported in Supplementary Note III. Nuclei are evolved in time using the PIOUD integrator$^{36}$ with a time step equal to 0.75 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. 36, where it is found to be optimal for both stochastic and deterministic forces. Simulations lasted around 6 ps, until the convergence of the vibrational modes at $\Gamma$ is reached. Forces are computed from the Born-Oppenheimer PES evaluated at DFT level within the QE package, or from the model PES defined in the “Potential energy surface parametrization” section. In case of ab initio PIMD, we used a BLYP functional for computing the PES. The wavefunction cut-off for the PES is set to 90 Ry (420 Ry for the charge density), while the Fermi smearing is Gaussian and set equal to 0.03 Ry. PIMD simulations are performed using $2 \times 2 \times 2$ real-space supercells, containing in each case 32 atoms, and the corresponding reciprocal-space mesh is always equal to $9 \times 9 \times 9$. We used a smaller plane-wave cutoff than the one used in single-point DFT calculations, because in PIMD we replaced the hard H Coulomb pseudopotential with a smoother PAW one. This has been necessary to speed up the PIMD calculations, which would otherwise have been too time consuming.

In our PIMD framework, non-adiabatic effects and permutation exchanges are not included. However, in this work we focused on the normal electronic state, where the system is purely metallic, with no relevant excitonic effects, and without superconducting gap opened. As far as quantum statistics is concerned$^{60}$, exchange effects should not play a major role in the range of volumes analyzed for $H_3S$. Indeed, by comparing the
density of hydrogen atoms in the largest volume we studied ($ρ_H$) against the typical density of pristine hydrogen where exchange effects are significant ($ρ < ρ_0 = 0.0435$ mol per cm$^3$)$^{61}$, we find $ρ_H ≥ 3ρ_0$. Furthermore, the average distance between hydrogen atoms ($d_{HH}$) in the range of volumes considered is always such that $d_{HH} ≥ 2$ Å and no pure roto-librations between two hydrogen atoms are present in the vibrational spectrum.

## SSCHA simulations

Besides the exact description of quantum nuclear motion provided by PIMD, one can also rely on approximated theories like the SSCHA$^{30}$, based on a variational principle on the free energy, which allows one to include quantum nuclear anharmonicity in a non-perturbative way. Here, we performed SSCHA simulations on the 3D H$_3$S (D$_3$S) model using up to 30,000 configurations. The average proton position (centroid) reported in Supplementary Fig. 5 (in Supplementary Note IV) are directly accessible through the SSCHA free energy minimization.

## Phonons

Harmonic phonons are obtained through DFPT simulations$^{62}$ as implemented within QE$^{51}$. The same set of DFT parameters and pseudopotentials employed for PIMD simulations were used to compute harmonic phonon dispersions, except for the k-space grid that was chosen equal to 18 × 18 × 18, as in this case it is referred to the unit cell. The results of these calculations are shown in Fig. 2. We specify that the DFPT shuttling mode frequency at q = I, reported in Fig. 3 has been computed with higher precision by employing the more accurate H Coulomb pseudopotential, requiring a plane-wave cutoff of 200 Ry.

Anharmonic phonon frequencies at PIMD level are evaluated by computing the zero frequency component of the phonon Matsubara Green's function from PIMD simulations. This method has been recently implemented in ref. 32 and it has been shown to describe accurately the vibron frequencies of solid phases of hydrogen. Conversely, within the SSCHA, auxiliary phonons are a byproduct of the free energy minimization. However, to get the physical phonons of Fig. 3, probed by spectroscopies, we apply the full self-energy dynamical corrections to the auxiliary dynamical matrix, described in detail in ref. 20, including both the third and fourth-order terms.

## Data availability

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

## Code availability

The TurboRVB code$^{52}$ used to carry out our QMC calculations is available under the GPLv3 license. It can be downloaded from the following open source repository: https://github.com/sissaschool/turborvb. In this work we have used Sandro Sorella's legacy version (v1.0.0), which can be found at this link: https://github.com/sissaschool/turborvb/releases/tag/v1.0.0.

Received: 28 July 2023; Accepted: 8 March 2024;
Published online: 15 March 2024

## References

1. Onnes, H. K. The superconductivity of mercury. *Comm. Phys. Lab. Univ. Leiden* **122**, 122–124 (1911).

2. Tresca, C. et al. Why mercury is a superconductor. *Phys. Rev. B* **106**, 180501 (2022).

3. Schilling, A., Cantoni, M., Guo, J. D. & Ott, H. R. Superconductivity above 130 K in the Hg-Ba-Ca-Cu-O system. *Nature* **363**, 56–58 (1993).

4. Bardeen, J., Cooper, L. N. & Schrieffer, J. R. Theory of superconductivity. *Phys. Rev.* **108**, 1175–1204 (1957).

5. Drozdov, A. P., Eremets, M. I., Troyan, I. A., Ksenofontov, V. & Shylin, S. I. Conventional superconductivity at 203 kelvin at high pressures in the sulfur hydride system. *Nature* **525**, 73–76 (2015).

6. Drozdov, A. et al. Superconductivity at 250 K in lanthanum hydride under high pressures. *Nature* **569**, 528–531 (2019).

7. Kong, P. et al. Superconductivity up to 243 k in the yttrium-hydrogen system under high pressure. *Nat. Commun.* **12**, 5075 (2021).

8. Service, R. F. At last, room temperature superconductivity achieved. *Science* **370**, 273–274 (2020).

9. Ferreira, P. P. et al. Search for ambient superconductivity in the Lu-N-H system. *Nat. Commun.* **14**, 5367 (2023).

10. Cataldo, S. D., von der Linden, W. & Boeri, L. First-principles search of hot superconductivity in La-X-H ternary hydrides. *npj Comput. Mater.* **8** https://doi.org/10.1038/s41524-021-00691-6 (2022).

11. Einaga, M. et al. Crystal structure of the superconducting phase of sulfur hydride. *Nat. Phys.* **12**, 835–838 (2016).

12. Mozaffari, S. et al. Superconducting phase diagram of H$_3$S under high magnetic fields. *Nat. Commun.* **10**, 2522 (2019).

13. Minkov, V. S., Prakapenka, V. B., Greenberg, E. & Eremets, M. I. A boosted critical temperature of 166 K in superconducting D$_3$S synthesized from elemental sulfur and hydrogen. *Angew. Chem. Int. Ed.* **59**, 18970–18974 (2020).

14. Osmond, I. et al. Clean-limit superconductivity in *Im$\overline{3}$mH$_3$S synthesized from sulfur and hydrogen donor ammonia borane. *Phys. Rev. B* **105**, 220502 (2022).

15. Goncharov, A. F., Lobanov, S. S., Prakapenka, V. B. & Greenberg, E. Stable high-pressure phases in the H-S system determined by chemically reacting hydrogen and sulfur. *Phys. Rev. B* **95**, 140101 (2017).

16. Duan, D. et al. Pressure-induced metallization of dense (H$_2$S)$_2$H$_2$ with high-T$_c$ superconductivity. *Sci. Rep.* **4**, 6968 (2014).

17. Sano, W., Koretsune, T., Tadano, T., Akashi, R. & Arita, R. Effect of van Hove singularities on high-T$_c$ superconductivity in h$_3$S. *Phys. Rev. B* **93**, 094525 (2016).

18. Errea, I. et al. Quantum hydrogen-bond symmetrization in the superconducting hydrogen sulfide system. *Nature* **532**, 81–84 (2016).

19. Bianco, R., Errea, I., Calandra, M. & Mauri, F. High-pressure phase diagram of hydrogen and deuterium sulfides from first principles: structural and vibrational properties including quantum and anharmonic effects. *Phys. Rev. B* **97**, 214101 (2018).

20. Monacelli, L. et al. The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects. *J. Phys. Condens. Matter* **33**, 363001 (2021).

21. Errea, I., Calandra, M. & Mauri, F. 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).

22. Bianco, R., Errea, I., Paulatto, L., Calandra, M. & Mauri, F. 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).

23. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. *Phys. Rev. Lett.* **77**, 3865–3868 (1996).

24. Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. *Phys. Rev. A* **38**, 3098–3100 (1988).

25. Lee, C., Yang, W. & Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. *Phys. Rev. B* **37**, 785–789 (1988).

26. Lester, W. A., Mitas, L. & Hammond, B. Quantum Monte Carlo for atoms, molecules and solids. *Chem. Phys. Lett.* **478**, 1–10 (2009).

27. Saritas, K., Mueller, T., Wagner, L. & Grossman, J. C. Investigation of a quantum Monte Carlo protocol to achieve high accuracy and high-throughput materials formation energies. *J. Chem. Theory Comput.* **13**, 1943–1951 (2017).

28. Raghav, A., Maezono, R., Hongo, K., Sorella, S. & Nakano, K. Toward chemical accuracy using the Jastrow correlated antisymmetrized
geminal power ansatz. *J. Chem. Theory Comput.* **19**, 2222–2229 (2023).

29. Foulkes, W. M. C., Mitas, L., Needs, R. J. & Rajagopal, G. Quantum Monte Carlo simulations of solids. *Rev. Mod. Phys.* **73**, 33–83 (2001).

30. Wagner, L. K. & Ceperley, D. M. Discovering correlated fermions using quantum Monte Carlo. *Rep. Prog. Phys.* **79**, 094501 (2016).

31. Casula, M., Filippi, C. & Sorella, S. Diffusion Monte Carlo method with lattice regularization. *Phys. Rev. Lett.* **95**, 100201 (2005).

32. Morresi, T., Paulatto, L., Vuilleumier, R. & Casula, M. Probing anharmonic phonons by quantum correlators: a path integral approach. *J. Chem. Phys.* **154**, 224108 (2021).

33. Morresi, T., Vuilleumier, R. & Casula, M. Hydrogen phase-IV characterization by full account of quantum anharmonicity. *Phys. Rev. B* **106**, 054109 (2022).

34. Vinet, P., Smith, J. R., Ferrante, J. & Rose, J. H. Temperature effects on the universal equation of state of solids. *Phys. Rev. B* **35**, 1945–1953 (1987).

35. Landau, L. D. On the theory of phase transitions. *I. Phys. Z. Sowjet.* **11**, 26 (1937).

36. Mouhat, F., Sorella, S., Vuilleumier, R., Saitta, A. M. & Casula, M. Fully quantum description of the Zundel ion: combining variational quantum Monte Carlo with path integral Langevin dynamics. *J. Chem. Theory Comput.* **13**, 2400–2417 (2017).

37. Srdinšek, M., Casula, M. & Vuilleumier, R. Rényi entropy of quantum anharmonic chain at non-zero temperature. *Phys. Rev. B* **108**, 245121 (2023).

38. Monacelli, L. & Mauri, F. Time-dependent self-consistent harmonic approximation: anharmonic nuclear quantum dynamics and time correlation functions. *Phys. Rev. B* **103**, 104305 (2021).

39. Akashi, R., Sano, W., Arita, R. & Tsuneyuki, S. Possible “Magnéli” phases and self-alloying in the superconducting sulfur hydride. *Phys. Rev. Lett.* **117**, 075503 (2016).

40. Azadi, S. & Kühne, T. D. High-pressure hydrogen sulfide by diffusion quantum Monte Carlo. *J. Chem. Phys.* **146**, 084503 (2017).

41. Goncharov, A. et al. Hydrogen sulfide at high pressure: change in stoichiometry. *Phys. Rev. B* **93**, 174105 (2016).

42. Li, Y. et al. Dissociation products and structures of solid H₂S at strong compression. *Phys. Rev. B* **93**, 020103 (2016).

43. Guigue, B., Marizy, A. & Loubeyre, P. Direct synthesis of pure H₃S from S and H elements: no evidence of the cubic superconducting phase up to 160 GPa. *Phys. Rev. B* **95**, 020104 (2017).

44. Cui, T. T., Chen, D., Li, J. C., Gao, W. & Jiang, Q. Favored decomposition paths of hydrogen sulfide at high pressure. *N. J. Phys.* **21**, 033023 (2019).

45. Pruzan, P. et al. Phase diagram of ice in the VII-VIII-X domain. Vibrational and structural data for strongly compressed ice VIII. *J. Raman Spectrosc.* **34**, 591–610 (2003).

46. Benoit, M., Marx, D. & Parrinello, M. Quantum effects on phase transitions in high-pressure ice. *Comput. Mater. Sci.* **10**, 88–93 (1998).

47. Tirelli, A., Tenti, G., Nakano, K. & Sorella, S. High-pressure hydrogen by machine learning and quantum Monte Carlo. *Phys. Rev. B* **106**, 041105 (2022).

48. Huang, C. & Rubenstein, B. M. Machine learning diffusion Monte Carlo forces. *J. Phys. Chem. A* **127**, 339–355 (2022).

49. Niu, H. et al. Stable solid molecular hydrogen above 900 K from a machine-learned potential trained with diffusion quantum Monte Carlo. *Phys. Rev. Lett.* **130**, 076102 (2023).

50. Giannozzi, P. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. *J. Phys. Condens. Matter* **21**, 395502 (2009).

51. Giannozzi, P. et al. Advanced capabilities for materials modelling with quantum ESPRESSO. *J. Phys. Condens. Matter* **29**, 465901 (2017).

52. Nakano, K. et al. TurboRVB: a many-body toolkit for ab initio electronic simulations by quantum Monte Carlo. *J. Chem. Phys.* **152**, 204121 (2020).

53. Bennett, M. C. et al. A new generation of effective core potentials for correlated calculations. *J. Chem. Phys.* **147**, 224106 (2017).

54. The Elk Code. https://elk.sourceforge.io/#documentation.

55. Kohn, W. & Sham, L. J. Self-consistent equations including exchange and correlation effects. *Phys. Rev.* **140**, 1133–1138 (1965).

56. Sorella, S., Devaux, N., Dagrada, M., Mazzola, G. & Casula, M. Geminal embedding scheme for optimal atomic basis set construction in correlated calculations. *J. Chem. Phys.* **143**, 244112 (2015).

57. Kato, T. On the eigenfunctions of many-particle systems in quantum mechanics. *Commun. Pure Appl. Math.* **10**, 151–177 (1957).

58. Umrigar, C. J., Toulouse, J., Filippi, C., Sorella, S. & Hennig, R. G. Alleviation of the fermion-sign problem by optimization of many-body wave functions. *Phys. Rev. Lett.* **98**, 110201 (2007).

59. Kwee, H., Zhang, S. & Krakauer, H. Finite-size correction in many-body electronic structure calculations. *Phys. Rev. Lett.* **100**, 126404 (2008).

60. Hirshberg, B., Rizzi, V. & Parrinello, M. Path integral molecular dynamics for bosons. *Proc. Natl Acad. Sci.* **116**, 21445–21449 (2019).

61. Mao, H.-k & Hemley, R. J. Ultrahigh-pressure transitions in solid hydrogen. *Rev. Mod. Phys.* **66**, 671–692 (1994).

62. Baroni, S., de Gironcoli, S., Dal Corso, A. & Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. *Rev. Mod. Phys.* **73**, 515–562 (2001).

## Acknowledgements

The authors thank M. Calandra, I. Errea, F. Mauri and L. Monacelli for useful discussions. They acknowledge computational resources provided by GENCI under the allocation number 0906493, which granted access to the HPC resources of IDRIS and TGCC. They also thank RIKEN for providing computational resources of the supercomputer Fugaku through the HPCI System Research Project ID hp220060. The authors are grateful to the European Centre of Excellence in Exascale Computing TREX-Targeting Real Chemical Accuracy at the Exascale, which supported this work. This project has received funding from the European Union's Horizon 2020 Research and Innovation program under Grant Agreement No. 952165.

## Author contributions

R.T. carried out the DFT and QMC calculations. M.Cherubini performed the SSCHA calculations. The PIMD phonons were computed by T.M. and M.Cherubini. All authors analyzed and interpreted the results. M.Cherubini, T.M. and M.Casula wrote the paper, with contributions from R.T. M.Casula conceived the project.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41524-024-01239-0.

**Correspondence** and requests for materials should be addressed to Michele Casula.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024