van der Waals driven anharmonic melting of the 3D charge density wave in VSe2

Understanding of charge-density wave (CDW) phases is a main challenge in condensed matter due to their presence in high-Tc superconductors or transition metal dichalcogenides (TMDs). Among TMDs, the origin of the CDW in VSe2 remains highly debated. Here, by means of inelastic x-ray scattering and first-principles calculations, we show that the CDW transition is driven by the collapse at 110 K of an acoustic mode at qCDW = (2.25 0 0.7) r.l.u. The softening starts below 225 K and expands over a wide region of the Brillouin zone, identifying the electron-phonon interaction as the driving force of the CDW. This is supported by our calculations that determine a large momentum-dependence of the electron-phonon matrix-elements that peak at the CDW wave vector. Our first-principles anharmonic calculations reproduce the temperature dependence of the soft mode and the TCDW onset only when considering the out-of-plane van der Waals interactions, which reveal crucial for the melting of the CDW phase.

T he study of electronic ordering and charge-density-wave (CDW) formation is attracting massive efforts in condensed matter physics 1 . In particular, its dynamical nature is the focus of a strong debate in correlated oxides and high-Tc superconducting cuprates 2 , where fluctuations of the charge order parameter 3 , dispersive CDW excitations 4 , and phonon anomalies 5 are observed. Microscopically, the subtle balance between electron-phonon interaction (EPI) and nested portions of the Fermi surface (singularities in the electronic dielectric function, χ q , at q CDW = 2k F ) determines the origin and stabilization of the charge periodicities 6 . While the Fermi surface nesting scenario survives for 1D and quasi-1D systems (Peierls transition), its role in higher dimensions remains largely questioned 7,8 .
Among the solids showing electronic charge ordering, layered transition metal dichalcogenides (TMDs) represent the first crystalline structures where 3D CDWs were discovered 9 . 1T-VSe 2 (space group P3m1) belongs to the series of layered TMDs that develops a 3D-CDW as a function of temperature, T CDW = 110 K. However, unlike the isostructural 1T-TiSe 2 , which adopts a commensurate 2 × 2 × 2 CDW ordering with q CDW = (0.5 0 0.5) r.l.u 10 , 1T-VSe 2 develops a more complex temperature dependence 3D incommensurate pattern in its CDW phase with a q CDW = (0.25 0 −0.3) r.l.u CDW wave vector 11 , modulating the interlayer distances. 1T-VSe 2 is rather unique among the 1T-polytypes because it develops anomalies in its transport properties and magnetic susceptibility 12 that more closely resemble those of 2H-polytypes (T CDW [2H-NbSe 2 ] = 33 K, T CDW [2H-TaSe 2 ] = 122 K) and presents the lowest onset temperature among them, i.e., T CDW [1T-TiSe 2 ] = 200 K, T CDW [1T-TaS 2 ] = 550 K 11 . The sizable difference between T CDW [1T-VSe 2 ] and its 1T counterparts can be attributed to the occurrence of large fluctuation effects that lower the mean-field transition temperature 13 or to the out-of-plane coupling 14 between neighboring VSe 2 layers assisted by the weak short-range van der Waals interactions 15 . Moreover, the theoretical input based on ab initio calculations is also limited for all these TMDs undergoing CDW transitions due to the breakdown of the standard harmonic approximation for phonons, which cannot explain the stability of the high-temperature undistorted phases 16 . This hinders the study of both the origin and the melting of the electronically modulated state, complicating the comprehensive understanding of the CDW formation.
From the electronic point of view, angle-resolved photoemission (ARPES) experiments in VSe 2 reported asymmetric dogbone electron pockets centered at M(L) 17 that follow the threefold symmetry of the Brillouin zone (BZ) interior, with nesting vectors closely matching those observed by x-ray scattering 18 . The formation of the CDW results from the 3D warping of the Fermi surface in the ML plane ( Fig. 1a shows the high-symmetry points of the Brillouin zone of the hexagonal lattice of VSe 2 ). Moreover, photoemission data also find a partial suppression of the density of states near E F on the nested portion below 180 K, indicating that a pseudogap opens at the Fermi surface 19 . However, a detailed investigation of the electronic structure is complicated by the 3D nature of the CDW order, and the momentum dependence of the EPI and the response of the lattice to the opening of the gap at E F remains unsolved. In fact, inelastic x-ray scattering (IXS) and theoretical calculations discarded the Fermi surface nesting scenario proposed for 2H-NbSe 2 20,21 and 1T-TiSe 2 22,23 and emphasized the critical role of the momentum dependence of the EPI. In addition, it has been recently demonstrated that large anharmonic effects are required to suppress the CDW phases in TMDs and understand their phase diagrams, both in the bulk and in the monolayer limit 16,[24][25][26][27] . Indeed, an evolution from the (4 × 4) CDW in bulk VSe 2 to a ( ffiffi ffi 7 p ffiffi ffi 3 p ) electronic reconstruction has been reported by means of scanning tunneling microscopy 28 , imperatively calling for a comprehensive description of the nature of the 3D CDW in VSe 2 .

Results
Quasi-elastic central peak. Figure 1b displays the temperature dependence of the elastic signal at the critical wave vector q CDW = (2.25 0 0.7) r.l.u upon cooling from 300 K. The elastic line due to incoherent scattering is barely visible at high temperature and is temperature independent down to 150 K, implying low structural disorder. Below 150 K ( Supplementary  Fig. 4), a smooth increase of the quasi-elastic intensity is observed at q = q CDW due to low-energy critical fluctuations and displays a sharp onset at the CDW transition T ≈ 110 K. No indications of charge instabilities were observed along the Γ → M and Γ → L directions. The mean-field critical exponent obtained in the disordered phase at T > T CDW , γ = 1.303 ± 0.004, is consistent with the existence of a 3D regime of critical fluctuations of an order parameter of dimensions n = 2, as expected for a classical XY universality class 29 . A similar critical exponent has been observed in the quasi-1D conductor blue bronze K 0.3 MoO 3 30 and ZrTe 3 31 , which develops a giant Kohn anomaly at the CDW transition.
Experimental and theoretical phonons. Figure 1c displays the momentum dependence of the inelastic spectra at (2 + h 0 0.7) r.l.u. for 0.15 < h < 0.45 at 300 K. Optical phonons appear above 17 meV and do not overlap with the acoustic branches. At all momentum transfers, 0 < h < 0.5, the spectrum consists of 2 phonons, labeled as ω 1 and ω 2 in Fig. 1d, in good agreement with the results of the theoretical calculations (see Supplementary Fig. 5 for a precise description and assignment of the 2 branches). The third acoustic mode is silent in IXS as its polarization vector is perpendicular to the wave vector. Both ω 1 and ω 2 belong to the same irreducible representation and, thus, do not cross. For h < 0.2, ω 1 develops more spectral weight than ω 2 and, for h > 0.2, the intensity of ω 2 increases and ω 1 leads an apparent asymmetric broadening of ω 2 , as depicted in Fig. 1d. To obtain quantitative information of the frequency and the phonon lifetime, the experimental scans were fitted using standard damped harmonic oscillator functions convoluted with the experimental resolution of ≈1.5 meV (see Fig. 1d and Supplementary Fig. 6 for a detailed analysis of the fitting). The frequencies of the low-energy acoustic branches ω 1 and ω 2 start around 4 and 8 meV, respectively, and end at ≈13 meV. Remarkably, the results of our ab initio anharmonic phonon calculations with the stochastic self-consistent harmonic approximation (SSCHA) [32][33][34] , which are performed with forces calculated within density-funcitonal theory (DFT) and including van der Waals interactions, show that both ω 1 and ω 2 do not follow a sinusoidal dispersion, but develop a dip at h ≈ 0.25 r.l.u. The theoretical dispersion nicely matches the experimental data from the zone center to the border of the Brillouin zone (BZ), as shown in Fig. 1e. In fact, the results of the harmonic phonon calculations indicate that the high-temperature structure of 1T-VSe 2 is unstable towards a CDW transition. It is clear, thus, that anharmonicity stabilizes 1T-VSe 2 at high temperatures. On the other hand, the linewidth extracted from the analysis (Fig. 1f, symbols) of the ω 2 mode is resolution limited across the whole BZ. Nevertheless, the linewidth of the ω 1 branch is no longer resolution limited between 0.2 < h < 0.3 r.l.u. and develops an anomalously large broadening of ≈4 meV at h = 0.25 r.l.u. Again, the experimental broadening is well captured by our calculations (dashed lines in Fig. 1f), indicating that the large enhancement of the broadening is mainly due to the EPI even if the anharmonic contribution to the linewidth also peaks at h = 0.25 r.l.u. (Supplementary Fig. 10).
Given the observation of the phonon broadening at room temperature and the good agreement between theory and experiment, we proceed with the analysis of the lattice dynamics at lower temperatures. At 250 K, the phonon with energy ≈7 meV (ω 2 ) shows a clear asymmetric broadening at q CDW , i.e, the corresponding branch ω 1 appears to develop a redshift as a function of temperature (Fig. 2b). The dispersion of ω 2 at 150 K is similar to the one at 300 K. Contrarily, ω 1 lowers its energy, softening from room temperature down to 110 K. The softening extends over a wide region of momentum space 0.225 < h < 0.3 r.l.u. (0.15 Å −1 ) at 150 K, see green dotted line in Fig. 2a (and Supplementary Fig. 8). The pronounced instability of this acoustic mode and its broad extension in momentum space are consistent with the results of our anharmonic phonon calculations (solid lines in Fig. 2c). The momentum space spread of the softening indicates a substantial localization of the phonon fluctuations in real space due to the EPI, questioning the pure nesting mechanism suggested by ARPES 17 . More importantly, the softening of this branch represents the first indication of the lattice response to the formation of the 3D-CDW in VSe 2 . The analysis of the linewidth reveals that the lifetime of ω 2 remains nearly constant across the BZ and is resolution limited (Fig. 2d). On the other hand, the softening of the ω 1 mode at 150 K is accompanied by an enhancement of the linewidth, as shown in Fig. 2d (6 meV linewidth at 120 K, Fig. 3f) and, again, well modeled by the ab initio calculations (dashed lines in Fig. 2d).
At the critical temperature, T CDW = 110 K and q ≈ q CDW , the spectrum is dominated by an elastic central peak at zero energy loss (FWHM = 0.05 r.l.u. and ΔE = 1.6 meV), thus, the soft mode is no longer resolvable (see Fig. 3a-d and Supplementary  Fig. 9). Figure 3e displays the temperature dependence of the soft mode, ω 1 , as well as the frequency of the phonon obtained ab initio with and without including van der Waals corrections. As plotted in Fig. 3e, the phonon frequency obtained ab initio follows the temperature dependence of the experimental acoustic branch. Moreover, the high temperature 1T structure of VSe 2 remains unstable at all temperatures after withdrawing the van der Waals functional from DFT calculation (see blue triangle in Fig. 3e). The softening of the acoustic phonon is accompanied by a linewidth broadening at T CDW (Fig. 3f).

Role of EPI
Having achieved a comprehensive description of the CDW and its temperature dependence, we address the crucial role of the EPI and nesting mechanism in the formation of the charge modulated state. In Fig. 4, we plot the calculated harmonic phonon frequency together with the electron-phonon linewidth of the three acoustic modes along q = (h 0 −1/3) r.l.u calculated within densityfunctional perturbation theory (DFPT). As it can be seen, the harmonic phonon instability of ω 1 coincides with a huge increase of its linewidth associated with the EPI. The softening and the increase of the electron-phonon linewidth specially affect the ω 1 mode, which suggests that the electron-phonon matrix elements are strongly mode and momentum dependent and have a strong impact on the real part of the phonon self-energy, which determines the harmonic phonon frequencies 8,21 . This behavior is similar to the one reported for 1T-TiSe 2 and 2H-NbSe 2 20,22 . The real part of the non-interacting susceptibility χ 0 (q), which captures the full Fermi surface topology and also affects the real part of the phonon self-energy (see Supplementary Information), has a softening of around 4% at q CDW , which seems insufficient to explain the large softening of the ω 1 mode. This suggests that the electron-phonon matrix elements are crucial to induce the harmonic softening and that the topology of the Fermi surface is not the driving mechanism. In order to further clarify the point, we calculate the so-called nesting function ζ(q) which measures the topology of the Fermi surface and peaks at the nesting q points (Supplementary Information). As shown in Fig. 4c, it peaks at q CDW , which indicates that the CDW vector coincides with a nested region of the Fermi surface. Considering that for constant electron-phonon matrix elements the nesting function coincides with the phonon linewidth given by the EPI, it is illustrative to compare them. Clearly, the phonon linewidth of the ω 1 mode coming from the EPI depends much more drastically on momentum than the nesting function: it changes by orders of magnitude as a function of q while the nesting function only by less than a factor of two. This is highlighted in the ratio between the linewidth and the nesting function plotted in Fig. 4d, which measures the momentum dependence of the electron-phonon matrix elements and should be flat if the electron-phonon matrix elements were constant. This ratio depends much more strongly on momentum than the nesting function itself and resembles the linewidth dependence, underlining again that the momentum dependence of the electron-phonon matrix elements plays a crucial role here. In conclusion, the EPI is the main driving force of the CDW transition in 1T-VSe 2 despite the presence of nesting at q CDW . Nevertheless, the q-range over which the phonon softens, Δq ≈ 0.075 r.l.u., even if it coincides with an increase of the electron-phonon linewidth, is a factor of 3 less than in 1T-TiSe 2 22 , where EPI and excitonic correlations are responsible for the structural instability and the CDW order, pointing to an intricate relationship between EPI and Fermi surface nesting scenarios in VSe 2 .

Discussion
Our anharmonic calculations, which predict that the ω 1 frequency vanishes between 75 and 110 K, are in good agreement with the experimentally measured phonon frequencies and the CDW temperature onset, T CDW = 110 K. When the SSCHA anharmonic calculation is repeated without including the van der Waals corrections (blue triangles in Fig. 3e), the softest acoustic mode at q CDW remains unstable even at room temperature. Remarkably, the weak van der Waals forces (of the order of 1mRy/a 0 for a typical SSCHA supercell calculation) are responsible for the stabilization of the 1T structure of VSe 2 and play a crucial role in melting the CDW. On the other hand, the damping ratio, Γ=ω q , increases upon cooling and the phonon becomes critically overdamped at q CDW and 110 K. The damping ratio Γ=ω q is given by ω 0 ¼ ðω 2 q À Γ 2 Þ 1=2 , where Γ is the linewidth ω q is the phonon energy renormalized by the real part of the susceptibility and ω 0 is the energy of the phonon fitted to damped harmonic oscillator. The critical exponent derived from the fitting of the phonon frequency vs reduced temperature ((T-T CDW )/ T CDW ), β = 0.52 ± 0.04, agrees with the square-root power law   expected from the mean-field theory (inset of Fig. 3f) and, therefore, fluctuation corrections are unnecessary to invoke the low T CDW of VSe 2 as compared to its 1T counterparts. The critical role of the EPI has been recently suggested by Raman scattering 35 and DFT calculations 36 . Indeed, revisited ARPES experiments 37 in NbSe 2 revealed a pronounced dispersion along k z discarding the nesting-driven CDW formation and leaving EPI as the major contributor 20 . Although our results indicate an EPI-driven CDW instability, nesting is present and, thus, the charge modulated ground state of VSe 2 has to be understood as an interplay between EPI and Fermi surface nesting scenarios.
In conclusion, we have observed with high-resolution IXS that the CDW transition in 1T-VSe 2 is driven by the collapse of an acoustic mode at q CDW = (0.25 0 −0.3) at T CDW = 110 K. The high-temperature 1T-VSe 2 phase is stable thanks to anharmonic effects. The observed wide softening in momentum space, the calculated strongly momentum dependent electron-phonon linewidth that peaks at q CDW , and the weaker dependence on the wave vector of the susceptibility suggest that the EPI is the main driving force of the CDW transition despite the presence of nesting. Moreover, the results show that van der Waals forces are responsible for the melting of the CDW. The dominant role of van der Waals forces here may be attributed to the out-of-plane nature of the CDW, which modulates the interlayer distance. This is not the case in 2H-NbSe 2 , where the bulk and monolayer transition temperatures seem to be similar 26,38 . This line of thinking is consistent with the enhancement of the CDW in monolayer VSe 2 , T CDW = 220 K 28 , since the out-of-plane van der Waals interactions are absent in this case. The critical role of out-of-plane coupling of layers has also been highlighted in the development of the 3D CDW in high-Tc cuprate superconductors [39][40][41] .

Methods
Sample growth and characterization. High-quality single crystals of VSe 2 with dimensions 2 × 2 × 0.05 mm 3 were grown by chemical vapor transport (CVT) using iodine as transport agent (see Supplementary Figs. 1, 2 for their structural, magnetic 42 and electronic characterization).
Inelastic x-ray scattering (IXS) measurements. The high-resolution IXS experiments were carried out using the HERIX spectrometer at the 30-ID beamline of the Advanced Photon Source (APS), Argonne National Laboratory. The incident beam energy was 23.72 keV and the energy and momentum resolution was 1.5 meV and 0.7 nm −1 , respectively, 43 . The components (hkl) of the scattering vector are expressed in reciprocal lattice units (r.l.u.), (hkl) = ha * + kb * + lc * , where a * , b * , and c * are the reciprocal lattice vectors. The experimental lattice constants of the hexagonal unit cell at room temperature are a = 3.346 Å, c = 6.096 Å, and γ = 120 ∘ . Here, we focus on the low-energy acoustic phonon branches dispersing along the (0 < h < 0.5 0 −0.3) r.l.u direction in the Brillouin zone near the reciprocal lattice vector G 201 , thus, in the range (2 + h 0 −0.3) r.l.u with 0 < h < 0.5.
First-principles calculations. The variational SSCHA 32-34 method was used to calculate temperature-dependent phonons fully accounting for non-perturbative anharmonic effects. The variational free energy minimization of the SSCHA was performed by calculating forces on 4 × 4 × 3 supercells (commensurate with q CDW ) making use of DFT within the Perdew-Burke-Ernzerhof (PBE) 44 parametrization of the exchange-correlation functional. Van der Waals corrections were included within Grimme's semiempirical approach 45 . Harmonic phonon frequencies and electron-phonon matrix elements were calculated within density-functional perturbation theory (DFPT) 46 . The force calculations in supercells needed for the SSCHA as well as the DFPT calculations were performed within the QUANTUM ESPRESSO package 47,48 (See Supplementary Information for further details on the calculations, which includes citations to refs. [49][50][51] ).

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. See author contributions for specific data sets.