Charge order from structured coupling in VSe$_2$

Charge order is ubiquitous in the phase diagrams of correlated materials. It is customarily described in the weak-coupling limit, where the electron-phonon coupling is implicitly assumed to be featureless. Detailed comparison between the resulting theoretical predictions and modern, high-resolution experimental data, however, regularly yields paradoxical results. As a best-case scenario for the weak-coupling approach, we consider $1T$-VSe$_2$, whose charge density wave (CDW) phase has a temperature-dependent, three-dimensional wave vector. We show that the electronic response in the weak-coupling limit cannot explain the observed thermal evolution, whereas consideration of a structured, momentum-dependent electron-phonon coupling does reproduce the measured result. The temperature dependence of the structured susceptibility moreover clarifies the apparent transition between two consecutive CDW phases observed in x-ray diffraction experiments, and suggests that a particle-hole asymmetry in the CDW gap results in its apparent absence in spectroscopic experiments. These results agree with recent findings in more strongly coupled compounds, and showcase both the use and the need of a structured coupling to obtain quantitative agreement with experiments in materials that would traditionally be considered unconventional.


I. INTRODUCTION
Density waves of charge, spin, or orbital occupation play a central role in determining the physical properties of many materials, ranging from elements [1,2], to cuprate high-T c superconductors [3][4][5], pnictides [6,7], complex oxides [8][9][10], and (multi)ferroics [11][12][13]. Understanding the mechanisms driving the formation of density waves is important for understanding their interplay with other types of order, and is key in artificially tuning phases of matter to obtain ideal properties for applications [14][15][16]. Here, we focus on the most common mechanism underlying charge order. Thus far, no broadly applicable theoretical framework quantitatively describes the properties of experimentally observed charge density wave (CDW) order [17,18].
In ideal single-band, one-dimensional (1D) metals, it is well-known that a Peierls transition will lead to CDW order accompanied by periodic lattice distortions [17,19]. This instability is signalled for example in the Lindhard function, which describes the (bare) electronic susceptibility and which diverges in 1D metals at the wave vector Q = 2k F . Despite the fact that this instability occurs for any strength of electron-phonon coupling (EPC), it is often referred to as the "weak coupling" scenario for CDW formation, implying that a perturbative field theory description is appropriate and the mean-field approximation is valid. In direct analogy with the Bardeen-Cooper- * vanwezel@uva.nl Schrieffer (BCS) theory for superconductivity, the meanfield relation between the zero-temperature energy gap and the transition temperature in weak coupling CDWs is 2∆(T = 0) = 3.52 k B T c [17]. Charge ordered materials in which the zero temperature gap is much larger than 3.5 k B T c are considered "strongly coupled". For those, higher-order scattering processes must be included in a field-theoretical description, and fluctuations typically suppress the transition temperature and give rise to a short-range ordered phase without long-range coherence above T c [17].
Both in weak and in strong coupling CDW materials, the instability towards charge order is signalled by a divergence of the Lindhard function. In the idealised Peierls scenario, this can be understood as the result of Fermi surface (FS) nesting: shifting electronic bands by the ordering wave vector maps (significant portions of) the Fermi surface onto itself. In real materials, however, perfect FS nesting never occurs, and inspection of the Lindhard function commonly indicates either no clear peak, or a dominant peak at a wave vector that does not correspond with the observed CDW wave vector [20]. Examples of supposedly nesting-driven density wave materials for which this has been explicitly demonstrated are 2H-NbSe 2 and 2H-TaSe 2 [20][21][22][23], TbTe 3 and other rare-earth tellurides [20,[24][25][26], blue bronze (K 0.3 MoO 3 ) [27], and chromium (which hosts spin density waves) [28]. Even for materials in which the electronic band structure is commonly considered well-nested, the assumption that their density waves arise purely from FS nesting is thus demonstrably incomplete.
One reason that the Lindhard function, and its strongcoupling corrections, commonly fail to predict the CDW ordering vector, is the neglect of any dependence of the electron-phonon coupling on crystal momentum, orbital character, or spin [17,20]. Here, we argue that a structured coupling is a primary component of any quantitative model for charge order, in both the strong and weak coupling regimes. The need for including the coupling structure has been previously argued for, but is yet to become standard practice [17,20,24,26,27]. Furthermore, this need has commonly been associated with the strong-coupling nature of specific materials, whereas here, we show it to be equally relevant for obtaining a quantitative understanding of weak-coupling materials. This is consistent with earlier approaches. For example, the concept of "hidden nesting" was introduced to account for the role of the real-space shape and orientation of valence orbitals in CDW formation, which effectively corresponds to including an orbital dependence in the EPC [25,[29][30][31]. Similarly, in the case of the prototypical (strong-coupling) CDW compound 2H-NbSe 2 , the momentum-and orbital-dependence of the EPC has been explicitly computed to determine the "structured" electronic susceptibility, and was shown to correctly predict the wave vector of the electronic instability [22,23,32,33].
To showcase the need for considering a structured coupling, even in weak-coupling CDW materials, we consider in this paper the CDW compound 1T -VSe 2 , which at first sight appears to be an ideal weak-coupling system. Its CDW gap was found in recent scanning tunnelling spectroscopy (STS) measurements to be 2∆ = 24 ± 6 meV at a temerpature of 5 K [16], while the critical temperature is approximately T c ≈ 110 K [14,[34][35][36]. This is close to the BCS ratio, placing VSe 2 firmly within the weak-coupling limit. Electronically, a single band of predominantly a single orbital character makes up the FS, consistent with the model Peierls scenario. Despite this apparent simplicity, the highly three-dimensional nature of the electronic bandstructure in VSe 2 make FS nesting weak at best, and the charge order in VSe 2 shows some unusual features.
Unlike other transition metal dichalcogenides (TMDC) such as 2H-NbSe 2 , 2H-TaSe 2 and 1T -TaS 2 , angleresolved photo-emission spectroscopy (ARPES) studies of 1T -VSe 2 show minimal or no gap openings in the spectral function when going from high to low temperature, such that the CDW gap structure remains unclear [37][38][39][40]. Furthermore, while the in-plane components of the three simultaneous CDW wave vectors Q i (i = 1, 2, 3) in VSe 2 are commensurate (periodicity 4a with lattice parameter a), it was determined via X-ray diffraction that their common out-of-plane component is incommensurate and varies from q z = 0.314 c * at 105 K to q z = 0.307 c * below 83 K with c * the length of the reciprocal lattice vector along k z [34]. The measured thermal evolution of the ordering wave vector was deemed anomalous, and led to the suggestion that this material might host two distinct CDW phases [34,41]. This is unusual, since both phases remain incommensurate, and the transition between them therefore is not of the common lock-in type. Phase contrast in satellite dark field images have led to the suggestion that the transition may be between a high-T , three-component CDW to a low-T phase with only two of the three symmetry-related Q i present, a so-called 2Q phase [41]. Although theoretically allowed, such a 2Q phase would itself be unusual, as it can only be stable in a small region of phase space and requires fine-tuned contributions from sixth order terms in a Landau expansion of the free energy [23,41]. Various scanning tunnelling microscopy (STM) experiments (down to 4.2 K) report a 3Q CDW phase [14,16,42,43], while others observe an enhanced intensity of just one or two CDW wave vectors [44,45]. It has been suggested that the latter could be the effect of an anisotropic STM tip [16] and/or spatial variations in the relative phases of the three CDW components [43]. Concrete evidence for a phase transition around 83 K from thermodynamic probes such as resistivity or specific heat is thus far lacking [14,35,36].
To understand the thermal evolution of the out-ofplane component of the Q i , and to determine whether VSe 2 is naturally susceptible to two distinct CDW ordering wave vectors, we compute the electronic response function, including the contribution of the momentumdependent EPC. We also use a computationally inexpensive method to find where gaps will open in the spectral function, and use it to predict the structure of the elusive CDW gap near the Fermi level.

II. STRUCTURED ELECTRONIC SUSCEPTIBILITY
To obtain the electronic structure of 1T -VSe 2 , we performed an ab-initio calculation within the local density approximation (LDA), using the all-electron fullpotential linear augmented plane wave (FLAPW) Elk code [46]. The resulting band structure is shown in Figure 1. We used the experimental lattice parameters of a = 3.356Å, c = 6.104Å, and the relative distance of the Se planes from the V planes z Se = 0.25 in the calculations [47]. Relaxation of the Se position did not significantly affect the FS or band structure near the Fermi level (E F ). We used a mesh of 32 × 32 × 24 k-points in the full Brillouin zone (> 2000 in the irreducible wedge) to achieve convergence. Unlike most TMDCs, the valence electrons in VSe 2 significantly disperse out-of-plane (along k z ). Our ab initio calculations show that only a single band crosses the Fermi level, of primarily Vanadium, 3dorbital character (see Fig. 1). We evaluated the energies for this band on a 100 × 100 × 400 k-point mesh for subsequent calculations of the Lindhard response function, as well as the structured susceptibility, which includes a momentum-dependent EPC. Phonons are expected to soften at those wave vectors for which the structured susceptibility, rather than the Lindhard function, shows the highest peaks [20,22].
For the EPC matrix elements, we use the expression derived by Varma et al., which has been well-tested for other transition metal compounds with d-orbital character at E F [23,48]. In the case of just a single band crossing E F , the expression simplifies to: Here, ξ k is the electronic dispersion taken from the density-functional theory calculation. The direction of the vector g k,k+q indicates the polarisation of the phonons coupled to. The displacements of Vanadium atoms associated with the CDW transition are known from earlier first-principles calculations to be purely inplane and longitudinal [49]. From here on, we therefore consider only the component of the EPC vector parallel to the (normalised) in-plane phonon momentum: g k,k+q = g k,k+q · q /|q |. The structured electronic susceptibility can be derived from a perturbative expansion of the phonon propagator. Since VSe 2 falls in the weak-coupling regime, it is sufficient to consider uncorrelated virtual electron-hole excitations. That is, we use the random phase approximation (RPA), and neglect vertex corrections, which should be small [50]. The renormalised phonon propagator is then described by D RP A = (D −1 0 − D 2 ) −1 , with bare phonon propagator D 0 and structured electronic susceptibility D 2 , given by [22,32]: Here, f (ξ) is the Fermi-Dirac distribution function and in our numerical computations we use a small regulator δ = 0.1 meV. The Lindhard function χ(q) is the same as the above, but taking g k,k+q = 1. In terms of the electronic susceptibility, the renormalised phonon dispersion is given by: with Ω 0 (q) the bare (high-temperature) phonon dispersion [23]. At the mean-field T c , the phonons will exhibit a Kohn anomaly such that Ω RP A (Q i ) = 0. This implies that, as long as Ω 0 (q ≈ Q i ) has no sharp features, the maximum of D 2 (q) close to the transition temperature should lie at q = Q i . In Fig. 2 we show χ(q) and D 2 (q) for three values of q z , while q x and q y span one reciprocal lattice vector each. Notice that D 2 (q) is not periodic across Brillouin zones, because of the projection of g k,k+q onto the inplane radial direction of q. We find that the highest peak in both χ and D 2 at q z = 0.31 lies at (q x , q y ) = (0, 0.25) (in rlu), in exact agreement with the in-plane value for the CDW wave vector found experimentally by Tsutsumi [34]. Comparing χ and D 2 directly, it is clear that the former disperses significantly less than the latter, and χ is far-removed from showing a divergence. This is indicative of the small degree of nesting in VSe 2 (see also the Appendix).
To better assess the exact out-of-plane position of the peak in D 2 , and to compare this to the experimentally determined values for the CDW wave vector, we compute D 2 for set values of q x and q y , while varying q z . We then track the peak position as a function of temperature, as shown in Fig. 3. This is the equivalent of finding the CDW ordering wave vector of VSe 2 if we were to quench the system from its high-temperature state directly to the chosen temperature.
A remarkable agreement between the thermal evolution of the peak position of D 2 from 50-105 K and the q z values observed by Tsutsumi [34] is found if we add a small downward shift of the chemical potential to our computed electronic structure. The shift of µ = −20 meV is within the error margin obtained by comparing the computed electronic structure to available ARPES spectra [37][38][39][40]. Moreover, it is well-documented that self-intercalation during the growth process of TMDCs often leads to small variations in the chemical potential [51]. Thermal evolutions for different values of the chemical potential shift are shown in the Appendix for completeness.
Although we reproduce the observed values and thermal evolution of Q z with a downward chemical shift of µ = −20 meV, the position of the peak in the structured susceptibility at T = 100 K along q y is moved to 0.26 rlu, away from the observed value of Q y = 0.25 rlu. (Notice that that we are limited by the resolution of the band-structure calculation, δq x = δq y = 0.01 rlu, and δq z = 0.0025 rlu.) This discrepancy is likely to be the result of lock-in effects, which are not taken into account in the present computations. Because the in-plane position of the peak in the electronic susceptibility is always close to the commensurate value of q y = 0.25 rlu, the coupling between the CDW order parameter and the lattice is likely to lock the in-plane component of the CDW wave vector into the lattice-preferred commensurate value [52,53]. In the out-of-plane direction, however, we expect no lock-in effect to take place because of the peak in susceptibility being further from commensurate values, the inter-layer coupling being weak [16], and the atomic displacements being purely in-plane [49]. This agrees with the observed values of q z remaining incommensurate at all temperatures [34].

III. THE CDW GAP
The CDW gap structure in VSe 2 has so far remained elusive. Two low-temperature ARPES measurements are suggestive of a gap around k z = 0.5c * [37,39,40]. The reported gap size of 80-100 meV in Ref. [37], however, refers to the shift in peak position of energy dispersion curves, while the CDW gap is known to be more closely related to the leading edge shift [54]. In Ref. [39,40] on the other hand, the sides of the Fermi pockets around the L-point lack some spectral weight already at temperatures above T c , while Strocov et al. report no gaps in their data [38]. The exact location and shape of the CDW gaps in VSe 2 thus remain to be determined conclusively.
Based on the energy dispersion of the band that makes up the Fermi surface, we can compute the spectral function A(k, ) probed by photoemission experiments. The electron propagator in the normal state, G 0 , is given by Here, iω n are Matsubara frequencies, while Σ k = Σ k + iΣ k is the complex-valued electron self-energy, and µ is the chemical potential. Wick rotating iω n to an energy and an infinitesimal imaginary part iδ, we obtain the spectral function, given by: We assume Σ k is independent of k, and set Σ + µ = −20 meV. Σ k describes the experimental resolution, for which we use a constant value of 10 meV. We obtain the normalstate spectral function at = 0 as shown on the left-hand side of the plots in Fig. 4. Knowing the normal-state electron propagator, we can go one step further and predict where gaps will open up in the spectral function as the CDW order sets in. To do so, we use a similar method to that developed in the context of superconductivity by Nambu [55] and Gor'kov [56]. Rather than constructing a new field theory starting from propagators with the symmetries of the ordered state, we complement the disordered propagators with additional, anomalous electron propagators F that do not conserve momentum up to one CDW wave vector: Here ψ † k creates an electron with momentum k, and k 1 − k 2 = ±Q i . Since the CDWs in VSe 2 are incommensurate in the out-of-plane direction, an infinite number of different F k1 k2 could be constructed. For the sake of computation, we approximate Q z i = 1 3 c * . Further noting that 12Q i = 0 and Q 1 + Q 2 + Q 3 = c * , we can represent all possible k 1 and k 2 by k + mQ 1 + nQ 2 (with m, n ∈ [0, 11]), as long as the gap function ∆(k) is k z -independent. Doing so, we generate a matrixĜ whose elements are the renormalised electron propagators G(k + mQ 1 + nQ 2 ) on the diagonal, and anomalous propagators F k1 k2 and (F k1 k2 ) † for any off-diagonal element with k 1 and k 2 differing by exactly one Q i . The integers m and n thus represent the matrix indices ofĜ. We then construct a matrix Dyson equation: Here,Ĝ 0 is a diagonal matrix of bare propagators evaluated at momenta k + mQ 1 + nQ 2 ;Σ is a matrix with self-energies on the diagonal, and gaps ∆ in off-diagonal elements connected by one Q i ; andÎ is the identity matrix. Solving this equation for the top left element of G, we find G(k), and hence the spectral function in the presence of a CDW gap. This method provides an inexpensive way to predict the gap structure of any CDW system, given only the bandstructure and the ordering wave vectors.
Starting from the computed electronic structure in the presence of a chemical potential shift and assuming a constant imaginary self-energy and a constant gap size ∆ = 15 meV we find a clear CDW gap opening at approximately 20 meV above E F . That the CDW gap is offset from the Fermi level is not surprising, since it is allowed by the absence of particle-hole symmetry (in contrast to for example a superconducting gap) and asymmetric CDW gaps have been observed in various TMDCs before [57]. The value of the offset found here is consistent with the scanning-tunnelling spectroscopy results of Jolie et al. [16], which show a suppression of the density of states with width 2∆ = (24 ± 6) meV centred at 10 meV above E F .
The clearest gaps in our computed spectral function appear on the long sides of the oval-shaped lobes around the L point (k z = 0.5c * , see Fig. 4). These gaps are connected to others at k z = 0.17c * by precisely one Q i , as indicated by the blue circles in Fig. 4. Planes at other k z values are either unaffected by the CDW order setting in, such as the k z = 0 plane, or they experience some moderate loss in spectral weight, like in the sides of the lobes at k z = 0.33c * .
It is probable that the true gap function ∆(k) in VSe 2 has a non-trivial k-dependence. Obtaining a selfconsistent solution for the full, momentum-dependent gap function is possible in principle [23], but the threedimensional nature of the electron dispersion implies a significant computational cost in the present case. Moreover, including any momentum dependence in the CDW gap will not allow gaps to open other than those already observed in Fig. 4. Rather, the k-dependence can only change the relative size of the gaps at different momenta, and the gap locations at 20 meV above E F highlighted by blue circles in Fig. 4 are the primary candidates for observing the elusive CDW gap in VSe 2 .

IV. CONCLUSION
In conclusion, we have shown that consideration of the structured electronic susceptibility is necessary for obtaining a quantitative understanding of (charge) density wave order, even for materials well within the weakcoupling limit. We showcase this for the case of 1T -VSe 2 , which is a three-dimensional CDW material with a ratio of the zero-temperature gap to the CDW transition temperature close to the BCS ideal.
We find that Lindhard response function or bare electronic susceptibility does not indicate any noticeable nesting in this material. On the other hand, the full structured susceptibility, including the momentumdependence of the electron-phonon coupling, shows a sharp peak at the observed CDW ordering vector. Moreover, its temperature dependence reproduces the thermal evolution of the CDW wave vectors observed in 1T -VSe 2 by X-ray diffraction experiments. Our results indicate that the thermal variation of the CDW wave vector in VSe 2 is a purely intrinsic effect, whose observation does not necessitate a description in terms of multiple consecutive CDW phases. Additionally considering the error margins of the reported x-ray diffraction data (indicated in Fig. 3), the lack of indicators for a second transition in bulk probes like resistivity [14,35,36], and the fact that satellite dark-field phase contrast may be due to a natural phase variation of the CDWs [43], our results suggest that VSe 2 hosts just a single CDW phase. The resolution of this paradox by the effect of having a strongly structured electronic susceptibility brings VSe 2 in line with more strongly coupled CDW materials in which a temperature dependence of the incommensurate CDW ordering wave vectors is commonly observed [27,52].
Based on a computation of the spectral function in the ordered phase, we predict the elusive CDW gap in VSe 2 to be centrered above E F , and to be most pronounced on the sides of the Fermi surface lobes around k z = 0.17c * and k z = 0.5c * . The offset in energy from the Fermi level is again owing to the momentum-dependence of the structured electron-phonon coupling. High-resolution, k z -resolved ARPES measured at varying temperatures should be able to resolve the predicted gaps, and verify its size and position in both momentum and energy.
The results presented here for the weakly coupled material VSe 2 highlight the need for considering the structure of the electron-phonon coupling in quantitative models for any density wave material, be it strongly coupled or weakly coupled. Besides the momentum-dependence considered here, the coupling in different materials may also have orbital and even spin dependence. All of these contribute to the detailed physical properties of density wave materials, and understanding their quantitative impact is indispensable in understanding the emergence of and interaction between charge, orbital, and magnetic order throughout (unconventional) superconductors, magnets, and (multi-)ferroics.