Anharmonicity and the isotope effect in superconducting lithium at high pressures: a first-principles approach

Recent experiments [Schaeffer 2015] have shown that lithium presents an extremely anomalous isotope effect in the 15-25 GPa pressure range. In this article we have calculated the anharmonic phonon dispersion of $\mathrm{^7Li}$ and $\mathrm{^6Li}$ under pressure, their superconducting transition temperatures, and the associated isotope effect. We have found a huge anharmonic renormalization of a transverse acoustic soft mode along $\Gamma$K in the fcc phase, the expected structure at the pressure range of interest. In fact, the anharmonic correction dynamically stabilizes the fcc phase above 25 GPa. However, we have not found any anomalous scaling of the superconducting temperature with the isotopic mass. Additionally, we have also analyzed whether the two lithium isotopes adopting different structures could explain the observed anomalous behavior. According to our enthalpy calculations including zero-point motion and anharmonicity it would not be possible in a stable regime.


I. INTRODUCTION
The strongly anomalous isotope effect recently measured in lithium in the 15-25 GPa pressure range 1 brought this element back under the spotlight. The reported superconducting critical temperatures (T c ) contrast starkly with the BCS theory, where T c is expected to scale as ∝ 1/M α , with M being the atomic mass and α the isotope coefficient (0.5 within the BCS theory). Actually, for most phonon mediated superconductors, α does not deviate much from 0.5. However, the above mentioned experiment shows a highly erratic behavior of α as a function of pressure, with values ranging from 1 to 4 from 15 to 21 GPa, decreasing sharply between 21 and 25 GPa, where it even becomes negative, with values as low as -2.
It is just another fascinating example of the rich and exotic phenomena emerging in lithium under pressure. The lightest metal on the periodic table shows a nearly free-electron bcc structure at ambient conditions 2 . Although it could be expected to evolve to an even more free-electron like system with increasing pressure, it has been shown that pressure not only induces several structural transformations [3][4][5][6][7] , but also gives rise to a plethora of fascinating physical properties 8 . For instance, lithium becomes a semiconductor near 80 GPa 9 , it shows a maximum in the melting line 10 and melts below ambient temperature (190 K) at around 50 GPa 3 . It also presents one of the highest T c for an element 1,[11][12][13][14][15] and it is expected to display a periodic undamped plasmon 16,17 . Additionally, according to a recent experiment lithium shows quantum and isotope effects in its low temperature and pressure phase transformations 7 .
Theoretical calculations within the harmonic approximation in fcc lithium show a highly softened transverse acoustic mode in the ΓK high-symmetry line 13,[18][19][20][21] . Around q inst = 2π/a(2/3, 2/3, 0), where a is the lattice parameter, this anomalous mode presents a huge electron-phonon coupling, becoming a key factor to explain the high T c observed in lithium [18][19][20] . This softening is associated to a well defined Fermi surface nesting 13,[18][19][20][21][22] and even yields imaginary phonon frequencies at pressures where fcc is known to be stable; the instability emerges at pressures higher than 30 GPa in the local density approximation (LDA), and at even lower pressures if one uses the generalized gradient approximation (GGA). As seen in other systems, such as simple cubic Ca 23 , PdH 24 , the record superconductor H 3 S 25 and NbSe 2 26 , anharmonicity is expected to have a significant role stabilizing this structure and, due to phonon frequency renormalization, also determining its superconducting properties 27 . As it has been measured at lower pressures of the phase diagram of lithium 7 , zeropoint vibrational energy could strongly impact the phase transitions of lithium in the 15-25 GPa pressure range, specially considering the small enthalpy differences between the most competitive candidates according to previous calculations 4,28,29 . In fact, the anharmonic correction to the vibrational energy could be significant as well.
The origin of the observed unconventional isotope effect in high pressure lithium remains unclear. Here we consider the following two hypothesis to explain this behavior. (i) Phonon frequencies scale with the atomic mass differently as expected within the harmonic approximation. Therefore, while in the harmonic approach the electron phonon coupling constant λ is independent of the isotopic mass, anharmonicity could make it differ from one isotope to the other, as it happens in palladium hydride 24 . (ii) 6 Li and 7 Li isotopes adopt different crystal structures due to the significant role of the vibrational energy in the phase diagram. Experimental evidence and previous theoretical calculations claim Li adopts the fcc phase from as low as 7 GPa to as high as 40 GPa in the temperature regime where superconductivity has been measured 3,5,6,9 . However, there is a considerable lack of experimental data in the mentioned region of the phase diagram and all previous calculations have been done in the static approach.
In this work we present an exhaustive analysis of the superconducting properties of fcc and cI16 structures of lithium in the 15-45 GPa pressure range, with vibrational degrees of freedom treated at the anharmonic level. We also analyze the possible existence of the hR1 phase in the pressure range of interest.

II. COMPUTATIONAL DETAILS
Our density functional theory (DFT) calculations were done within the Perdew-Burke-Ernzerhof (PBE) parametrization of the GGA 30 . Harmonic phonon fre- quencies and the electron-phonon deformation potential were calculated within density functional perturbation theory (DFPT) 31 as implemented in Quantum ESPRESSO 32 . The electron-proton interaction was considered making use of an ultrasoft pseudopotential 33 which includes 1s and 2s electrons. Anharmonic calculations, including the vibrational contribution to the enthalpy, were performed using the stochastic selfconsistent harmonic approximation (SSCHA) 34 . Anharmonic force constant matrices of fcc lithium were obtained by calculating forces in 3×3×3 supercells. Therefore, anharmonic dynamical matrices were obtained in the respective commensurate q-point grids and interpolated to a finer 9 × 9 × 9 mesh afterwards. These were combined with DFPT electron-phonon calculations obtained in the fine 9 × 9 × 9 mesh to calculate the anharmonic Eliashberg function α 2 F (ω). The same procedure was used for the cI16 structure, being 2 × 2 × 2 and 6 × 6 × 6 the coarse and fine grids respectively. The vibrational contribution to the enthalpy of hR1, which is a distortion of the fcc phase, was calculated using a 2 × 2 × 2 grid for obtaining anharmonic force constant matrices and interpolating the differences with respect to the undistorted fcc structure. More details and the convergence parameters are given in the Supplementary Material. Fig. 1 shows the DFPT harmonic phonon dispersion of fcc 7 Li at 26 GPa and the anharmonic corrections calculated within the SSCHA. Anharmonic force constant matrices were obtained by calculating forces in 3 × 3 × 3 and 4 × 4 × 4 supercells. Consequently, anharmonic dynamical matrices were obtained in the respective commensu-rate q-point grids. We see that anharmonicity is primarily localized around the phonon softening at the transverse acoustic T 1 branch at q inst , where the frequency is strongly shifted up by anharmonic effects. This well known phonon softening has been widely analyzed and explained in terms of Fermi surface nesting 13,[18][19][20][21][22] and, as shown in Fig. 2, it even yields imaginary frequencies at pressures higher than 25 GPa; a considerably lower pressure than the 30 GPa obtained within the LDA. In the same graph we also show the anharmonic frequency of the same mode, confirming fcc lithium is dynamically stabilized by anharmonicity above 25 GPa. However, as it is shown in the inset and even though this soft mode shows huge anharmonic effects, its frequency scales practically as in the harmonic case ( ω ∝ 1/M ). Despite the large anharmonicity, a similar harmonic scaling was previously calculated for high pressure simple cubic calcium 23 .

III. RESULTS AND DISCUSSION
Our DFPT electron-phonon coupling calculations displayed in Fig. 3 show the total coupling constant λ rises abruptly with increasing pressure in the fcc phase. Starting from an already high value of 0.85 at 15 GPa and reaching a value as high as 2.6 at 36 GPa, this dramatic growth is directly related to the also rapid increase of the electron-phonon linewidth γ of the T 1 mode at q inst , which doubles its value in the mentioned pressure range. The remarkable peak in the Eliashberg function α 2 F (ω) and the associated abrupt growth of the integrated electron-phonon coupling constant λ(ω) around the frequency of the anomaly is another indicator of how relevant this softening is in the superconducting properties of fcc lithium. However, while the phonon renormalization of the mentioned mode due to anharmonicity is huge, λ is nearly identical for both isotopes at every pressure except at 35 GPa, where the difference is just 7%, even if anharmonicity is already really strong. As mentioned above, this is due to the fact that the frequency of the anomalous mode scales harmonically. Our λ values are slightly larger than the ones by Maheswari et al. 20 and Profeta et al. 18 and quite larger than the ones by Akashi et al. 19 and Bazhirov et al. 13 . We attribute these disagreements to the large dependence of λ with the qpoint grid. While we used a 9 × 9 × 9 sampling of the BZ for the electron-phonon and lattice dynamics calculations, where q inst is explicitly taken into account, the mentioned works use 8 × 8 × 8 grids (7 × 7 × 7 in the case of Maheswari et al.), where it is not. According to our convergence tests, those grids clearly underestimate λ due to the absence of q inst in the grid (see Supplementary Material). Including this extremely anharmonic anomalous point is crucial for estimating the impact of anharmonicity in the electron-phonon coupling and, as a consequence, the superconducting T c .
Considering that for large electron-phonon coupling constants the McMillan equation underestimates the superconducting T c 35 , we solved the isotropic Migdal-Eliashberg equations 36,37 . We estimated a µ * value of Li fcc-Ref. [18] fcc-Ref. [13] fcc-Ref. [19] fcc-Ref. [20] fcc-Ref. [29] cI16-Ref. [29] Figure 3. Total electron-phonon coupling constant λ of fcc and cI16 lithium calculated for its two isotopes at different pressures. The inset shows the phonon linewidth of the T1 mode of fcc Li at qinst multiplied by the atomic mass, the product being independent of the phonon frequency and the isotopic mass. the calculated λ is compared to previous calculations 13, 18-20,29 . 0.17 using the Morel-Anderson formula 38 : . (1) The average electron-electron Coulomb repulsion term µ was obtained from Thomas-Fermi screening theory, a free-electron Fermi energy ε f was chosen, and the Debye cutoff phonon frequency ω D was taken as the highest frequency of the longitudinal acoustic modes 39 . Changes in phonon frequencies and electronic density for different pressures and isotopes only alter the fourth significant digit of µ * , so that differences in µ * cannot explain the isotope effect anomalies and we assume the same value for both isotopes. Fig. 4 shows the superconducting critical temperature of fcc lithium for both isotopes at 15, 20, 26 and 36 GPa. We find T c increases monotonically with pressure the same way λ does, ranging from 11.2 K (10.7 K) at 15 GPa to 34.8 K (32.5 K) at 36 GPa for 6 Li ( 7 Li). As in the case of λ, we do not see any anomalous scaling of the superconducting temperature with the isotopic mass; as it can be seen in fcc-SCDFT Ref. [18] fcc-µ*=0.13 MM Ref. [13] fcc-µ*=0.2 MM Ref. [13] fcc-SCDFT+Plasmons Ref. [19] fcc-µ*=0.13 MM Ref. [20] fcc-µ*=0.23 ME Ref. [29] cI16-µ*=0.12 ME Ref. [29] fcc 6 Li-µ*=0.17 ME (our work) cI16 6 Li-µ*=0.17 ME (our work) fcc 7 Li-µ*=0.17 ME (our work)  Li Ref. [1] natural Li Ref. [11] natural Li Ref. [12] natural Li Ref. [14] natural Li Ref. [5] fcc 6 Li-µ*=0.17 ME (our work) cI16 6 Li-µ*=0.17 ME (our work) fcc 7 Li-µ*=0.17 ME (our work) cI16 7 Li-µ*=0.17 ME (our work) (b) Estimated Tc of fcc and cI16 lithium for its two stable isotopes at different pressures (lines with symbols) and comparison with experimental values (characters) 1,5,11,12,14 . fects should not be isotope dependent and, due to the harmonic scaling of phonon frequencies, we do not expect vertex corrections to yield any anomalous isotope effect either. Therefore, we discard hypothesis (i).
After discarding that the anomalous isotope effect comes from strong anharmonicity in the fcc phase, we analyzed the possibility of the two isotopes showing different structures at the same pressure in a thermodynamically stable way. Fig. 6 shows the enthalpies of the competing phases cI16 and hR1 relative to their respective fcc ones for the two isotopes. Our static calculations, i.e. not including zero-point energy (ZPE), compare well with literature (there are no previous works including ZPE) 28 and just show the fcc to cI16 transi- Li: fcc,  Figure 5. The isotope coefficient α against pressure. Lines with symbols show the coefficients obtained for the cases in which the two isotopes adopt the same crystal structure (either cI16 or fcc). Curves without symbols show the coefficients for the cases in which the isotopes adopt different structures.
tion. No important changes are shown for both isotopes when anharmonic ZPE is included and, although in the pressure range where this phase transition happens the enthalpy difference with the hR1 is less than 1 meV per atom, that is, roughly the same as the error one assumes when converging total energy calculations within DFT, it remains metastable. Therefore, small changes in the calculation parameters or the choice of exchange and correlation potential might cause modifications in the transition pressures and phase sequence. Accordingly, when ZPE is included the fcc to cI16 transition pressure shifts from 37 GPa to 33 GPa for both isotopes, as the enthalpy difference is reduced by around 3 meV due to lattice vibrations. Additionally, in the 21-25 GPa pressure range, where the inverse isotope effect was observed, the enthalpy difference between cI16 and fcc structures is really small (around 4-6 meV/atom). In conclusion, our results do not support hypothesis (ii) as 6 Li and 7 Li isotopes are not expected to adopt different stable crystal structures.
Due to the extremely small enthalpy differences metastable coexistence of phases can not be discarded as it happens at ambient pressure for its martensitic transition 7 . In order to see if 6 Li and 7 Li adopting different structures could lead to the observed anomalous isotope effect, we have also made lattice dynamics and electron-phonon coupling calculations in the cI16 structure. We do not further consider hR1 as a candidate because, according to our calculations, the local minimum in the total energy surface associated to hR1 disappears for pressures lower than 28 GPa (see Supplementary Material). In Fig. 3 Figure 6. Relative enthalpies of cI16 and hR1 6 Li and 7 Li isotopes with respect to their fcc counterparts. In solid and dashed lines ZPE has been included, while in dotted curves only electronic energy has been considered. The low pressure limit for the hR1 curves has been set at the pressure which corresponds, in each case, to the maximum volume at which the phase shows a local minimum in the total energy surface (see Supplementary Material).
varies only between 0.9 and 1.2 in the 15-44 GPa pressure range. λ is fairly similar for both isotopes, so that anharmonicity does not have almost any impact. Actually, at the lowest pressures, cI16 values differ more than the fcc ones from one isotope to the other. This is due to the fact that, while the overall phonon spectrum is very slightly modified by anharmonicity in the cI16 phase, anharmonic corrections occur mostly at the lowest frequencies, which are the ones that contribute most to the total electronphonon coupling. In fact, our λ and T c estimations, with µ * =0.17 obtained with the Morel-Anderson formula as in the fcc case, shown in Figs. 3 and 4 yield values higher than in fcc below 20 GPa, being the opposite at higher pressures. The isotope effect coefficient is close to the harmonic value at 27 and 44 GPa, with α =0.42 and 0.57, respectively, while it deviates considerably from 0.5 at 15 and 19 GPa as it yields α =0.77 and 0.34, respectively. All this agrees with the higher anharmonicity we found at lower pressures. Although our enthalpy calculations do not predict both isotopes can stabilize in different structures, we have also analyzed this metastability driven hypothetical scenario: 6 Li stabilizing in the fcc phase and 7 Li in the cI16, and viceversa. As shown in Fig. 5, in the pressure range where the inverse isotope effect was experimentally observed (21-25 GPa), experimental values would only be qualitatively reproduced if 6 Li adopted the cI16 structure while 7 Li were in the fcc phase. This qualitative picture does not vary much if one uses the McMillan formula with µ * =0.22, but it could notably change if we used different µ * values for the different phases.

IV. CONCLUSIONS
According to our calculations, even though anharmonicity is crucial to stabilize the fcc phase in lithium under pressure, its λ remains almost the same for both isotopes and yields a conventional scaling of T c with isotopic mass and, therefore, it does not explain the experimentally observed anomalous isotope effect. On the other hand, including anharmonic ZPE in the enthalpy curve does not modify lithium phase diagram in the pressure range of interest, so that it is unexpected to have both isotopes in different structures. The anomalous isotope effect could only be qualitatively explained if 7 Li adopted the fcc structure while 6 Li adopted the cI16 one in a metastable way. All these, added to the large error bars and quite chaotic behavior of T c with pressure in Ref. 1-with considerably different temperature values for the same pressure-puts in question the experimental observation of an anomalous isotope effect in lithium at high pressure. This way, our work encourages further research to determine the phase sequence and superconducting properties of the two stable isotopes of lithium.

Supplementary Material
Anharmonicity and the isotope effect in superconducting lithium at high pressures: a first-principles approach  1 Phonon spectra and electron-phonon coupling 1

.1 fcc structure
Harmonic dynamical matrices of fcc Li have been obtained in a 9 × 9 × 9 q-grid for every analyzed pressure and isotope. A proper convergence of phonon frequencies required a 30×30×30 k-point grid and a Methfessel-Paxton smearing width of 0.01 Ry for electronic integrations in the first BZ. An energy cutoff of 65 Ry was necessary for expanding the wave-functions in the plane-wave basis. The electron-proton interaction was considered making use of an ultrasoft pseudopotential [1], in which 1s 2 core electrons where also included (the same pseudopotential has been used for the whole work). In Fig. 1 we show the harmonic phonon spectra obtained at 15,21,26 and 36 GPa for 7 Li after Fourier interpolating the dynamical matrices from the 9 × 9 × 9 grid to the desired path.
Anharmonic dynamical matrices where obtained in a 3 × 3 × 3 q-grid commensurate to the supercell size for our SSCHA calculations, in which we calculate forces acting on atoms. The difference between the harmonic and anharmonic dynamical matrices was interpolated to the finer 9 × 9 × 9 grid, so that anharmonic dynamical matrices were obtained in 9 × 9 × 9 grid after adding the harmonic dynamical matrix to the interpolation.
Electron-phonon matrix elements where calculated within DFPT, where converging the double Dirac delta in the equation for the phonon linewidth required a denser 80 × 80 × 80 k-point mesh. Superconducting T c was calculated solving isotropic Migdal-Eliashberg equations considering that for large electron-phonon coupling constants McMil- lan's equation underestimates T c . In Fig. 3a we can see how converging T c and λ with the q-point grid becomes tedious due to the large contribution of q inst to the total electron-phonon coupling. Our chosen 9 × 9 × 9 grid overestimate T c by 2.5 K comparing to the 12 × 12 × 12 case while λ is 0.2 larger. However, increasing the grid size would make the calculation really demanding and, our goal being to check whether anharmonicity could explain the anomalous isotope effect, this overestimation would only make anharmonic effects to be more visible. However, even in this case anharmonic effects are not big enough to explain the anomalous isotope effect. Moreover, we clearly see grids not containing q inst (dimensions not multiple of 3) yield smaller T c and λ values that the ones they do, and using such grids would obviously neglect how anharmonicity affects the electron-phonon coupling and superconductivity.

cI16 structure
Harmonic dynamical matrices have been obtained in a 6× 6× 6 q-grid for every analyzed pressure and isotope. A proper convergence of phonon frequencies required a 16×16×16 k-point grid and Methfessel-Paxton smearing width of 0.01 Ry for electronic integrations in the first BZ. An energy cutoff of 65 Ry was necessary for expanding the wave-functions in the plane-wave basis. In Fig. 2 we show the phonon spectra obtained at 18, 27 and 44 GPa for both 6 Li and 7 Li after Fourier interpolating the dynamical matrices from the 6 × 6 × 6 grid to the desired path.
Anharmonic dynamical matrices where obtained in a 2 × 2 × 2 q-grid, commensurate to the supercell in which the SSCHA was perfomed. We interpolated the results to the finer 6 × 6 × 6 grid with the same method as in the fcc case. In this case anharmonicity has practically no influence on phonon frequencies at 27 and 44 GPa, while at 15 and 19 GPa low frequency modes are more noticeably affected. We can see this in Fig. 2 for 6 Li (we do not show the result for 7 Li as they are practically identical).
Converging the double Dirac delta in the equation for the phonon linewidth required a 32 × 32 × 32 k-point mesh. Superconducting T c was calculated solving isotropic Migdal-Eliashberg equations. Converging T c within 1 K required to calculate the electronphonon matrix elements in a 6 × 6 × 6 q-point grid (see Fig. 3b).
The cI16 structure (Space Group I-43d) has all the Li atoms placed in the Wyckoff 16c positions (conventional coordinates (x, x, x) and all symmetry equivalent) , which has a free parameter x. As the SSCHA miniminization of the free energy is also performed with respect to x, final average atomic positions are different from the harmonic or static ones. In principle, one should perform the electron-phonon coupling calculations in the new anharmonic atomic positions for each isotope and pressure. However, we checked that the impact on λ and T c for 6 Li at 19 GPa, where the change in x is the greatest (∆x =0.004) is within the convergence criteria. Therefore, we use the electron-phonon coupling calculation calculated at the static equilibrium positions at each pressure for both isotopes.

Enthalpy curves
For obtaining the enthalpy H = E T + P V of the different structures, we calculated each contribution to the total energy E T = E el + E v at several unit-cell volumes and fitted them separately, due to the fact that the computational cost of a data point differs significantly from one contribution to another, as electronic energy E el is faster to compute than the vibrational E v one.
We calculated E el for fcc and ci16 Li for volumes per atom ranging from 50 to 100 a 3 0 with a step size of approximately 1.5 a 3 0 . We fitted the data using a Birch-Murnaghan equation of state. Due to the different properties of the phonon spectra, the vibrational contribution required a different treatment for each crystal structure. We find convenient to write the total vibrational contribution as E v = E f req + < V − V > [2], where E f req comes from the sum of the SSCHA frequencies over all the modes of the crystal and < V − V > comes from the difference of the actual anharmonic energy surface and the SSCHA harmonic one. E f req can be further splitted into the harmonic contribution and the anharmonic correction, E f req = E har + E anh , where E har is the energy coming from the harmonic frequencies.
For cI16 calculating harmonic dynamical matrices in a 2 × 2 × 2 q-grid was enough to converge E har within 0.5 meV/atom. We calculated E har at seven different volumes, from 50 to 100 a 3 0 , and fitted a fourth order polynomial to the data points. We calculated E anh and < V − V > at four different volumes per atom (60,70,80 and 84 a 3 0 ) by perfoming SSCHA calculations in 2 × 2 × 2 supercells, and fitted the data with a second order polynomial.
Fcc Li presents a more complex situation due to the anomaly in the Γ-K path. We computed E har using a 8 × 8 × 8 grid, which does not show any imaginary frequency down to at least 65 a 3 0 /atom (which corresponds to around 35 GPa), and converges E har within 0.2 meV/atom (see Fig. 4). We calculated E har at seven different volumes, from 50 to 100 a 3 0 , and fitted a fourth order polynomial to the data points.  Figure 5: Phonon spectra of hR1 7 Li at around 40 GPa. Harmonic dynamical matrices have been explicitly calculated in a 6 × 6 × 6 q-grid, while the anharmonic ones have been calculated in a 2 × 2 × 2 grid and interpolated to the finer 6 × 6 × 6 q-grid afterwards. The harmonic spectrum shows phonon instabilities in large regions of the BZ. Anharmonicity renormalizes strongly those instabilities, yielding real frequencies for every q-point in the 2 × 2 × 2 grid. However, after interpolating to the 6 × 6 × 6 grid some modes remain unstable.
the anharmonic contribution, we performed SSCHA calculations to obtain anharmonic dynamical matrices and < V − V > in a 3 × 3 × 3 grid for four different volumes (66,72,77 and 84 a 3 0 ). To overcome the situation of using different grids for each contribution of the vibrational energy, we needed to treat E anh carefully. We interpolated our anharmonic dynamical matrices from the 3 × 3 × 3 grid to a finer 9 × 9 × 9 one to obtain E f req , and substracted the harmonic contribution in a 8 × 8 × 8 grid as E anh = E f req − E har , as imaginary frequencies prevent us obtaining E har in a 9 × 9 × 9 grid. Finally, we fitted these four data points with a second order polynomial.
For hR1 we have proceeded in a different way due to the fact that it shows plenty of imaginary frequencies in the harmonic phonon spectra (see Fig. 5). These imaginary frequencies are strongly renormalized by anharmonicity and become real after applying the SSCHA in a 2 × 2 × 2 q-grid. However, the interpolation method was not useful in this case as some of the interpolated anharmonic matrices in a 6 × 6 × 6 remained yielding imaginary frequencies. We overcame this situation making use of the similarity of hR1 with the fcc phase. If one chooses a rhombohedral unit cell, hR1 differs from fcc only by the rhombohedral angle α. Thus, taking α and the unit cell volume V as variables, we can focus our attention to their associated potential energy surface. We define the total energy as E T (α, V ) = E T,f cc (V )+∆E(α, V ), where E T,f cc (V ) is the total energy of the fcc phase (α =60 o ) and ∆E(α, V ) is the difference in energy due to the change in rhombohedral angle. We only need to calculate ∆E(α, V ) in this case as we had previously calculated E T,f cc (V ). ∆E(α, V ) is the sum of electronic and vibrational contributions. The electronic contribution ∆E el (α, V ) is easily obtained by DFT total  energy calculations. For obtaining the vibrational contribution ∆E v (α) we assumed that it is independent of the unit cell volume. This way, we performed SSCHA calculations in 2 × 2 × 2 supercells at four different cos(α) values (0.25,0.35,0.412 and 0.5) for a single volume (60 a 3 0 ) and fitted it with a 3 rd order polynomial. In Fig. 6 we show ∆E(V, α) against α for different choices of the unit cell volume V, which is kept constant in each curve. Two relative minima can be distinguished below 70 a 3 0 : one at cos α = 0.5, which corresponds to the fcc structure, and another one corresponding to the hR1 phase, which even has a lower energy than the previous one for volumes smaller than 63 a 3 0 . Plus, the angle at which this minimum occurs increases with decreasing volume. Above 70 a 3 0 hR1 could not exist as it lacks of a local energy minimum.
In Fig. 7 we show the pressure vs. volume curves for each isotope and structure, obtained by taking the first derivative of E T with respect to the volume.