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

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 \emph{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$\bar{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$\bar{3}$m-R3m 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 ($\approx$ 50 GPa) of the classical ferroelectric transition pressure at 200K and in a large isotope shift ($\approx$ 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.


I. INTRODUCTION
Since its discovery in 1911 [1], 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) [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 [3], although the pairing mechanism in these materials is considered unconventional and it is not explained by the standard Bardeen-Cooper-Schrieffer (BCS) theory [4].
In 2015, the discovery of conventional superconductivity in H 3 S with a maximum T c of 203 K reached at a pressure P c as high as 150 GPa [5] 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 3 S claimed to have measured an even higher T c in LaH 10 , superconducting already at 250 K [6], later followed by a similar discovery in the yttrium hydride [7].In a rush towards 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 [8,9].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 * michele.casula@sorbonne-universite.frall 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 towards ternary materials [10].
In this work, we focus on the prototypical case of H 3 S and we study its structural phase transition generally associated with the maximum of the superconducting critical temperature, located at around 150 GPa [11][12][13][14].According to x-ray diffraction data [15], at lower pressures the sulfur (S) sites are arranged in a geometry that is compatible with the trigonal R3m symmetry (Fig. 1(b)) and, upon compression, the system undergoes a phase transition towards a body-centered-cubic (bcc) Im 3m structure (Fig. 1(a)).conductivity in H 3 S [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 [17].Moreover, theoretical studies more oriented to understand the underlying structural properties of H 3 S, revealed a significant disagreement in the transition pressures between the predicted phases.
In those works [18,19], the structural phase transition is explained by a quantum proton symmetrization from the R3m phase, with displaced protons, to the Im 3m 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 Im 3m 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) [20][21][22], 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) [23] and the Becke-Lee-Yang-Parr (BLYP) [24,25] 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 [26][27][28], thanks to its stochastic Green's function algorithms [29,30], such as the lattice regularized diffusion Monte Carlo [31], projecting any initial trial wavefunction towards 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 ex-act, 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.

A. Harmonic and anharmonic phonons
At high pressure, above 150 GPa, the H 3 S crystal is expected to be in the cubic Im 3m symmetric phase (Fig. 1(a)), 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. 1(b).In our description, we introduce a simplification by neglecting the trigonal distortion, which is however very weak (<0.6°) [19].Thus, the R3m phase considered here differs from the Im 3m one just by the hydrogen positions.
In Fig. 2, we report the analysis of the phonon dispersion for different volumes of the cubic Im 3m 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.Hereafter, volumes and energies will be expressed per H 3 S unit, while the unit cell will be taken as cubic with V=90.9 a 0 3 Harmonic (6 x 6 x 6 q-grid) PIMD (2 x 2 x 2 q-grid) V=94.9 a 0 3 Harmonic (6 x 6 x 6 q-grid) PIMD (2 x 2 x 2 q-grid) Harmonic (6 x 6 x 6 q-grid) PIMD (2 x 2 x 2 q-grid)

Γ
V=110.0 a 0 3 Harmonic (6 x 6 x 6 q-grid) PIMD (2 x 2 x 2 q-grid) Phonon dispersion for different volumes of the cubic Im 3m unit cell, whose electronic structure has been computed with the BLYP functional.At V=83.6 a 3 0 (volume per H3S unit), only the harmonic dispersion is reported.For all the other volumes, we compare the harmonic dispersion (light-blue color) with the PIMD one (green color).Full dispersions are obtained by interpolating harmonic (anharmonic) dynamical matrices defined on a 6x6x6 (2x2x2) q-grid.PIMD simulations are performed at T =200 K.

S atoms arranged in a bcc lattice.
At this point, it is important to underline that the DFPT and the PIMD phonons bear different information (see also Sec.IV F).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 3 0 only the harmonic dispersion is reported in Fig. 2, for larger volumes we compare the PIMD phonons (green lines) obtained in a 2x2x2 supercell with the harmonic ones (light-blue lines).We notice that for PIMD phonons, the spatial range of the force constant matrix is such that the 2x2x2 supercell is large enough to allow for a q-interpolation of the phonon branches ω m = ω 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 3 0 (see Figs. 2 and 3), the appearance of imaginary frequencies indicates the dynamical instability of the Im 3m structure.More specifically, the softening of the shuttling hydrogen mode at q = Γ signals the transition towards the asymmetric R3m 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 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 [20].Imaginary phonon components are represented with negative values.
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 3 0 and 100 a 3 0 , [34] 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 equation of state (EOS) P = P (V ) using the Vinet relation [35] 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 Sec.IV A).A simple comparison of the EOS produced by the two approaches, shown in Fig. 9(a), 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.

B. 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 = Γ, 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 Sec.IV B.
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 3 0 , 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 3 0 , the height of the double well barrier for QMC is ∼270 meV/H 3 S, 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 [36].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 3 0 corresponding to a pressure of 263 GPa, while for QMC we found the same volume (≈ 85 a 3 0 ) which corresponds to 238 GPa in this case (see Fig. 9(a)).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.

C. 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 as implemented in Ref. [37].Numerical details of these simulations can be found in Sec.IV D. Here, we mention only that in the PIMD simulations of our 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. 5(a), 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. 5(b), 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, 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 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. 5(b) reveals that the QMC PES leads to a smaller volume for the peak splitting, as shown already in the contour plot of Fig. 5(a).
By fitting the distribution in Fig. 5(b) 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.Moreover, in PIMD we can easily quantify isotope effects by replacing the hydrogen with the deuterium mass.
We estimate the peak splitting to take place in H 3 S at a volume of 99.6 a 3 0 for DFT-BLYP, and of 96. of 152 GPa for QMC.These values, reported in Tab.I, 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. 9(a)).
The occurrance 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 rigourous 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 the hydrogen shuttle between the two PES minima.As 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. 6(a-b).This quantity has already been used in a previous work [38] 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 3 0 for DFT-BLYP, corresponding to P c = 131 GPa, and V c = 100.9 a 3 0 for QMC, corresponding to P c = 126 GPa (Tab.I).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. 6(c) 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 favour of a paraelectric phase, by squeezing the compound.Notice that this does not correspond to the ferroelectric transition pressure, associated instead with the global Im 3m-R3m symmetry breaking, and long-range dipole order, which happens at lower values.The same analysis is carried out for both the H 3 S and D 3 S crystals, to estimate the magnitude of isotope effects.We summarize the results in Tab.I, where we show that the hydrogen-to-deuterium substitution brings about an increase of the local polar- ization formation pressure that falls into the [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] GPa range.

D. Full BLYP-PIMD solution of the H3S 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 → Im 3m transformation.The suitable order parameter to identify this transition is ∆ = ⟨ 1 N N i=1 δx i ⟩, where the sum runs over all the N hydrogen atoms in the supercell, and δx i is the distance of the i-th proton from the S-S midpoint at a given snapshot.An equivalent order parameter, showing usually less statistical fluctuations, is ∆ abs = ⟨| 1 N N i=1 δ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 3 0 , P ferro ≃ 82 GPa, as found in our BLYP-PIMD simulations in the 2×2×2 supercell.We notice that the peak location is correlated with the jump in the shuttling mode frequency, reported in both Fig. 3 and 7.The agreement between the ∆ and ∆ abs susceptibilities and the shuttling mode frequency jump strengthens the reliability of our estimate.Interestingly, at 200K 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 mo- ments and Im 3m 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 200K (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 200K.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-theart 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 SS-CHA approximation is not able to capture the disordered Im 3m phase, being a mean-field theory with no configurational entropy and no direct information of imaginary time correlations, key to detect the dynamical local moment formation [39].Thus, the critical pressure SSCHA can normally compute is the ferroelectric one, P ferro , and not P c .To investigate more deeply this mismatch, we Between ferro and para, an Im 3m phase with disordered local moments is expected to take place [38].We also report the ferroelectric critical pressures predicted by SSCHA, from both ab initio [19] (white square) and our 3D-PES model (green squares).For comparison, the pressure of the experimental Tc maximum is shown by a teal diamond.
carry out SSCHA calculations with our model PES (see Sec. IV E for details).In SSCHA, the occurrence of the asymmetric R3m phase is signalled by a centroid displaced with respect to the S-S midpoint.The SSCHA critical values are V ferro = 107.8 a 3 0 , P ferro = 114 GPa for the DFT-BLYP PES, and V ferro = 102.4 a 3 0 , P ferro = 118 GPa for the QMC PES.As in PIMD, there is no significant difference between the electronic structure method 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 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 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 ferro .Moreover, we also observe an increase of the SSCHA phonon line-width across the transition of the order of 10 cm −1 .A similar jump in the shuttling phonon fre-quencies 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 Im 3m to a phase sharing the same Im 3m 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 Im 3m 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.

III. CONCLUSIONS
In this work, starting from ab initio electronic structure calculations, we generated a model PES to describe the shuttling mode of hydrogen in H 3 S, responsible for the R3m → Im 3m 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 H 3 S 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 [40,41], worsened by the fact that other competing symmetries could be stable in the same pressure range [42][43][44][45].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 H 3 S and D 3 S, as we can see in Fig. 9(b).
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 [38], and it cannot be associated with the ferroelectric R3m −→ Im 3m transformation, which takes place at a lower pressure P 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.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 the R3m → Im 3m 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 [46,47].However, it is still a matter of debate whether the transformation is direct or whether another intermediate disordered structure appears, 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 [48][49][50] with the aim at including a larger variety of hydrogen configurations in PIMD calculations by keeping the same QMC accuracy.

A. Electronic structure calculations for the PES model
For the DFT electronic structure calculations, we used the Quantum Espresso (QE) suite of codes [51,52], while for the QMC calculations, we employed the TurboRVB package [53].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 [54] 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 calcula-  [11] and from Minkov et al. in [13].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.
tions 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 [55].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 Fig.

S.1 of the SI).
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 32x32x32.
For the QMC calculations, we used a Slater-Jastrow wavefunction Ψ, which reads as: where the term exp (J) is the Jastrow factor, symmetric under electron exchange, while Φ S is the antisymmetric Slater determinant.The Slater orbitals in Φ S are generated by DFT calculations within the Local Density Approximation (LDA) [56], 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 [57].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 electronnucleus distances r iI .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 electronelectron cusp conditions, respectively.They read as J 1 = Ne i=1 N I=1 u I (r iI ), and J 2 = Ne i<j=1 v(r ij ), where N (N e ) is the number of nuclei (electrons) in the supercell, and the functions u and v are defined as follows: with a and b variational parameters, and Z I the charge of the I-th pseudoatom.The coefficients in Eqs. 2 and 3 are set to fulfill the Kato cusp conditions for electron-nucleus and electron-electron coalescence, respectively [58].
J 3 is the three-body term that accounts for the electron-electron-nucleus interactions.As defined in Tur-boRVB, 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 manybody 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 [59].We then carried out lattice regularized diffusion Monte Carlo (LRDMC) calculations [31], to stochastically project the initial wavefunction towards 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 of the SI.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 2x2x1, 2x2x2, 3x2x2 and 3x3x2 real-space supercells in order to extrapolate the LRDMC total energy to the thermodynamic limit, by also using Kwee-Zhang-Krakauer (KZK) [60] 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 Sec. IV B).

B. 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 ϕ, defined along the axis connecting the two flanking sulfur atoms (S-S axis): x is the position of the hydrogen atom along the S-S axis, r is the radial distance from the S-S axis, and ϕ 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, ϕ) ≡ (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 3m symmetry (with all hydrogen atoms sitting at the S-S midpoints) to the R3m 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: The functional form of our 3D PES is constructed as follows: E(x, r, ϕ) = A(x, r) + B(x, r) sin(ϕ + 5π/4), (4) with: and where f min and f max are defined as: 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π-periodicity.In particular, for x < 0.5, E has a minimum given by f min at ϕ = π/4 and the maximum f max at ϕ = 5π/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, ϕ + π) = E(x, r, ϕ), 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 ϕ = [π/4, 5π/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, 5 and 6.The root mean square error of these fits amounts to ≈ 1 meV/H 3 S.

C. Equation of state
In order to get the pressure associated to each volume, P = P (V ), we use the Vinet EOS: with η = (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 [35] 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 3m phase, where the total energy is computed from first principles, either by DFT-BLYP or by QMC, on a grid of volumes (see Fig. 9(a)).In the fit, we disregarded the Zero-Point Energy (ZPE) contribution, because we verified that the ZPE variation is very small (< 1 mHa/H 3 S) in the range of pressures analyzed here within the same Im 3m phase.

D. 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 of the SI.Nuclei are evolved in time using the PIOUD integrator [37] with a time step equal to 0.75 fs and a friction parameter of the Langevin thermostat equal to 1.46•10 −3 atomic units.The latter value is the same as in Ref. [37], 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 Γ 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 Sec.IV B. In case of ab initio PIMD, we used a BLYP functional for computing the PES.The wavefunction cutoff 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 2x2x2 real-space supercells, containing in each case 32 atoms, and the corresponding reciprocal-space mesh is always equal to 9x9x9.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.

E. SSCHA simulations
Besides the exact description of quantum nuclear motion provided by PIMD, one can also rely on approximated theories like the SSCHA [20], based on a variational principle on the free energy, which allows one to include quantum nuclear anharmonicity in a nonperturbative way.Here, we performed SSCHA simulations on the 3D H 3 S (D 3 S) model using up to 30000 configurations.The average proton position (centroid) reported in Supplementary Fig. S.5 of the SI are directly accessible through the SSCHA free energy minimization.

F. Phonons
Harmonic phonons are obtained through DFPT simulations [61] as implemented within QE [52].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 18x18x18, 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 = Γ, 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 Ryd.
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 fourthorder terms.

ACKNOWLEDGMENTS
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 partially supported this work.This project has received funding from the European Unions Horizon 2020 Research and Innovation program under Grant Agreement No. 952165.

FIG. 1 .
FIG. 1.(a) Crystal structure of the Im 3m symmetric phase of H3S (smaller volume, higher-pressure phase).(b) The R3m asymmetric phase (larger volume, lower-pressure phase).dSS is the lattice parameter of the bcc crystal.
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[20].Imaginary phonon components are represented with negative values.

58 QMCFIG. 4 .
FIG.4.3D-PES landscapes (in meV/H3S units) for different volumes, 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.Top panel: DFT-BLYP PES.Middle panel: 1D projection along the axis connecting the two minima.Bottom panel: 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.

FIG. 5 .
FIG. 5. (a) Two-dimensional projection on the y = z plane of hydrogen positions with our 3D model solved by PIMD.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.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.|δx| is the distance from the origin.

FIG. 6 .
FIG. 6. Local moment susceptibility for H3S (panel (a)) and D3S (panel (b)) for quantum nuclei with the DFT-BLYP (blue) and QMC (red) 3D-PES.The vertical dashed lines indicates the susceptibility maxima.Panel (c): Imaginary time correlator g(β/2) for H3S.Blue 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.

FIG. 7 .
FIG. 7. Ferroelectric transition from ab initio BLYP PIMD simulations at 200K in a 2×2×2 supercell.Black (blue) points are values of the ferroelectric order parameter ∆ (∆ 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).
FIG.8.Phase diagram at 200K 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 Im 3m phase (bracketed by the blue star and the purple pentagon) are reported.Between ferro and para, an Im 3m phase with disordered local moments is expected to take place[38].We also report the ferroelectric critical pressures predicted by SSCHA, from both ab initio[19] (white square) and our 3D-PES model (green squares).For comparison, the pressure of the experimental Tc maximum is shown by a teal diamond.

FIG. 9 .
FIG. 9. (a) Equation of state P = P (V ) (see Sec. IV C and Eq.7) for both DFT-BLYP (blue color) and QMC (red color) calculations of H3S in the Im 3m phase.The transition volumes and the corresponding pressures identified at different levels of theory using the density probe are displayed by dashed lines.(b) Tc as a function of pressure in H3S and D3S from Drozdov et al. in [5], Einaga et al. in[11] and from Minkov et al. in[13].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.