Dynamical Stability of Face Centered Cubic Lithium at 25 GPa

1Centro de Física de Materiales CFM, CSIC-UPV/EHU, Paseo Manuel de Lardizabal 5, 20018 Donostia/San Sebastián,Basque Country, Spain 2Donostia International Physics Center (DIPC), Manuel Lardizabal pasealekua 4, 20018 Donostia/San Sebastián, Basque Country, Spain 3Fisika Aplikatua 1 Saila, Bilboko Ingeniaritza Eskola, University of the Basque Country (UPV/EHU), Rafael Moreno “Pitxitxi” Pasealekua 3, 48013 Bilbao, Basque Country, Spain 4Departamento de Física de la Materia Condensada, University of the Basque Country (UPV/EHU), 48080 Bilbao, Basque Country, Spain


Introduction
Lithium, the lightest metal on the periodic table, is a nearly-free-electron material which adopts a bcc structure at ambient conditions [1].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 [2,3,4], but also gives rise to a plethora of fascinating physical properties [5] normally induced by nonlocality in the core pseudopotential with pressure [6,7].For instance, lithium becomes a semiconductor near 80 GPa [8], it melts below ambient temperature (190 K) at around 50 GPa [2] and displays a periodic undamped plasmon according to theoretical calculations [9].Moreover, it presents one of the highest superconducting critical temperatures (Tc) for an element, reaching values as high as 15 K at around 30 GPa [10][11][12][13].Recently, the measurement of an inverse isotope effect in the 20-26 GPa range [10] has put this element back under the spotlight.
Experimental evidence [2,3,4,8] shows that Tc rapidly soars in the face centered cubic (fcc) phase, which is adopted by lithium roughly between 7 and 40 GPa, pressure after which it transforms to the rhombohedral hR1 phase.This phase is really similar to fcc as it consists of just a distortion along the c axis if one switches to a hexagonal representation.The transformation to the cubic cI16 phase occurs shortly after, at around 43 GPa.
Theoretical calculations within the harmonic approximation in fcc lithium show a highly softened transverse acoustic branch in the Γ-K high-symmetry line [12,14,15,16,17].At qinst=2π/a(2/3,2/3,0), this anomalous mode presents a huge electron-phonon coupling, becoming a key factor to explain the high Tc observed in lithium [12,14,16].This softening is associated to a well-defined Fermi surface nesting [14,17] 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, and at even lower pressures if one uses the generalized gradient approximation.As seen in other systems, such as simple cubic calcium [18], anharmonicity is expected to have a significant role stabilizing this structure and, due to phonon frequency renormalization, also determining its superconducting properties.
In this work we analyze the mentioned anomalous vibrational mode in terms of ab-initio total energy calculations.We study the convergence of the harmonic phonon frequency of the unstable mode and we estimate the impact of anharmonicity assuming it does not interact with the rest of the modes in the crystal, by solving the time independent Schrödinger equation for the potential associated to the atomic displacements along the mode of interest.

Computational details
Our density functional theory (DFT) calculations were done within the Perdew-Burke-Ernzerhof parametrization of the generalized-gradient approximation [19].The electron-ion interaction was considered making use of an ultrasoft pseudopotential [20] as implemented in Quantum ESPRESSO [21].An energy cutoff of 885 eV was necessary for expanding the wave-functions in the plane-wave basis.Converging the total electronic energy within 0.1 meV.required a 30x30x30 k-point mesh and a 0.136 eV (0.01 Ry) Hermite-Gaussian electronic smearing for electronic integrations in the first Brillouin zone (BZ).Harmonic phonon frequencies were also calculated within density functional perturbation theory (DFPT) [22,23].

Numerical convergence of harmonic phonons
In Fig. 1 we show the DFPT harmonic phonon spectrum of fcc lithium at 25 GPa.The dynamical matrices have been explicitly calculated in a 9x9x9 q-point grid and Fourier interpolated afterwards to the chosen path in the BZ.The transverse acoustic branch T1 is clearly softened in the Γ-K path, displaying even imaginary (negative) frequencies at qinst.This point is part of the 9x9x9 grid in which an explicit DFPT dynamical matrix calculation was performed, which discards it being an artifact of the Fourier interpolation.
In order to better understand the nature of this anomaly, we obtained its associated potential energy profile by calculating the total electronic energy for different displacements following the polarization vector of the T1 mode at qinst.In Fig. 2 we plot the total energy along the path described by the anomalous mode for different choices of k-point grids and smearing widths for electronic integrations.While for displacements larger than 0.05 Å the curves do not vary much among the different calculations, for small displacements the chosen convergence criteria of 0.1 meV is clearly not enough for describing the profile properly.As one increases the BZ sampling and decreases the smearing width, the double-well shape of the energy surface sharpens.As the harmonic phonon frequency is obtained by calculating the second derivative at the equilibrium position, we expect the convergence of this value to be harsh.This is confirmed by DFPT calculations.In Fig. 3 we can see that, while the frequency of the L longitudinal and T2 transverse modes converge within 1 cm⁻¹, the frequency of the T1 transverse anomalous mode decreases monotonically as smearing values are reduced for the densest samplings of the BZ.A proper convergence would require an even denser grid and a smaller smearing width.However, such small energies and displacements are physically meaningless considering the zero point fluctuations.Thus, we hold to the choice of a 30x30x30 k-point mesh and a 0.01 Ry smearing width.
We estimated the harmonic frequency of the T1 mode at qinst by calculating the second derivative at the equilibrium position numerically.By fitting a parabola to the three data-points closest to the equilibrium position (see Fig. 4), we have estimated a phonon frequency of -28.6 cm -1 .The value obtained within DFPT for the same calculation parameters is -22.1 cm -1 , in quite good agreement with the previous frozen phonon result.

Solving the Schrödinger equation
Due to the highly anharmonic shape of the energy profile of the mode under analysis, we solved the one-dimensional time-independent Schrödinger equation for the potential associated to the vibrational mode.This way, we are neglecting the interaction of the analyzed mode with the rest of the normal modes of the lattice.We fitted a 10 th degree polynomial to the energy vs. displacement data points and solved the eigenvalue problem using ParametricNDSolve in Mathematica [24].We assumed that for displacements larger than 0.6 Å the eigenfunctions vanish, discarding tunneling to adjacent positions.
In Fig. 4 we show the data points, the fitted polynomial and the obtained three lowest energy eigenfunctions and eigenvalues, with values 6.3, 20.6 and 36.3 meV respectively.The possible transitions between states yield phonon frequencies of 14.3, 15.7 and 30.0 meV (115.3,121.8 and 241 cm -1 ) respectively.If we compare the value of the most probable transition, 115.3 cm -1 , which is the transition between the ground state and the first excited state, with the harmonic value of -28.6 cm -1 , we see a huge difference.First, the magnitude of the anharmonic correction is even larger than the harmonic value itself.Second, the shape and width of the wavefunctions and the magnitude of the eigenvalues confirms not only the lack of physical relevance of the convergence problems observed when computing the second derivative, but also the absolute failure of the harmonic approach for describing this vibrational mode.Finally, anharmonicity changes the sign of the frequency, and, thus, guarantees the dynamical stability of the system within this frozen-phonon approach.

Conclusions
We estimated the anharmonic correction of the anomalous transverse phonon at 2/3 ΓK in fcc lithium at 25 GPa within the frozen-phonon approximation by solving the Schrödinger equation.

011103-4
While within the harmonic approach one obtains a frequency of -28.6 cm -1 , anharmonicity renormalizes dramatically this value up to 115.3 cm⁻¹.Thus, it seems that anharmonicity is dynamically stabilizing this system.Plus, this strong phonon renormalization could potentially cause the inverse isotope effect [10] by making the electron-phonon coupling constant differ for the two stable lithium isotopes considering that in PdH anharmonicity explains the inverse isotope effect [25].Anyway, our frozen-phonon calculations neglect the interaction of the analyzed vibrational mode with the rest of the modes in the crystal.Therefore, a proper quantitative analysis would require of a method taking into account the phonon-phonon interaction also among different modes.

Fig. 2 .
Fig. 2. Electronic energy vs. displacement along the T1 mode at qinst for different choices of k-point grids and smearing width values.In red (30x30x30) and blue (39x39x39) curves the same smearing width values have been chosen as in the black ones (21x21x21).The inset shows a zoomed view of the small displacements region.The origin of energy is set at the zero displacement energy for a 30x30x30 grid and 0.01 Ry smearing width.b) Convergence analysis of the DFPT harmonic frequency of the three normal modes at qinst.

Fig. 3 .
Fig. 3. Convergence analysis of the DFPT harmonic frequency of the three modes at qinst.

Fig. 4 .
Fig. 4. Total energy (black circles), the fitted 10 th order polynomial (black curve), the parabolic fit at equilibrium position (red), and the eigenfunctions (dashed blue) obtained by solving the Schrödinger equation for the polynomic potential.The amplitude of the eigenfunctions is arbitrary for better visualization, and the respective values are shifted by the corresponding eigenvalues.The inset shows the data points and the parabolic fit at the equilibrium position in the small displacements region.