Weak Dimensionality Dependence and Dominant Role of Ionic Fluctuations in the Charge-Density-Wave Transition of NbSe_{2}.

Contradictory experiments have been reported about the dimensionality effect on the charge-density-wave transition in 2H NbSe_{2}. While scanning tunneling experiments on single layers grown by molecular beam epitaxy measure a charge-density-wave transition temperature in the monolayer similar to the bulk, around 33 K, Raman experiments on exfoliated samples observe a large enhancement of the transition temperature up to 145 K. By employing a nonperturbative approach to deal with anharmonicity, we calculate from first principles the temperature dependence of the phonon spectra both for bulk and monolayer. In both cases, the charge-density-wave transition temperature is estimated as the temperature at which the phonon energy of the mode driving the structural instability vanishes. The obtained transition temperature in the bulk is around 59 K, in rather good agreement with experiments, and it is just slightly increased in the single-layer limit to 73 K, showing the weak dependence of the transition on dimensionality. Environmental factors could motivate the disagreement between the transition temperatures reported by experiments. Our analysis also demonstrates the predominance of ionic fluctuations over electronic ones in the melting of the charge-density-wave order.

A charge-density wave (CDW) is a structural distortion of the crystal lattice that induces a modulation of the electronic density [1]. CDWs in transition metal dichalcogenides (TMDs) [2] are particularly relevant because they seem to compete with superconductivity [3][4][5][6][7][8][9] and form a phase diagram similar to the high-temperature superconductors [10,11]. The origin of the CDW instability and its interaction with superconductivity in TMDs continues to be a topic of intense debate [12][13][14][15], as it might enlighten the hidden physics in strongly correlated materials where different phases compete.
TMDs tend to adopt layered crystal structures. Each layer is formed by transition metal atoms sandwiched by covalently bonded chalcogens (see Fig. 1). Due to the weak van der Waals interaction that holds together the layers in the bulk, TMDs can be exfoliated down to the monolayer [16]. Monolayer TMDs can currently also be synthesized by chemical means [17]. This has opened the possibility to study the effect of dimensionality on CDW transitions. The results obtained thus far, however, do not show a clear trend and it is generally not clear whether experimental results are affected by environmental factors. In monolayer TaS 2 the CDW present in the bulk disappears [18], while in the isolectronic and isostructural TaSe 2 it remains unchanged [19]. Even if NbS 2 has no CDW transition in the bulk [20], a CDW transition has been observed in the monolayer grown on top of bilayer graphene [21]. This result seems substrate dependent, as a monolayer grown on Au(111) does not show any CDW feature down to 30 K [22]. Similarly, in monolayer TiSe 2 the CDW temperature (T CDW ) is enhanced with respect to the bulk, but the value of T CDW is strongly substrate dependent [23,24].
The effect of thickness on the T CDW in NbSe 2 is even more controversial. The most stable polytype of NbSe 2 is the 2H shown in Fig. 1, which undergoes a CDW transition at 33 K to an incommensurate structure very close to a 3 × 3 × 1 ordering [3,31]. Inelastic x-ray experiments evidence that the CDW transition is second order. In fact, a longitudinal acoustic (LA) phonon mode collapses exactly at 33 K at q CDW ∼ 2=3ΓM [32], a momentum consistent with the periodicity of the CDW phase. In other TMDs as well, the phonon frequency of the mode that drives the CDW is strongly reduced by cooling and eventually vanishes at T CDW [33,34]. The CDW transition temperature in the NbSe 2 monolayer has been determined by Raman measurements in exfoliated samples on sapphire substrates [35,36], as well as by scanning tunneling microscopy (STM) in single layers grown by molecular beam epitaxy (MBE) on bilayer graphene, which confirmed that the CDW order remains 3 × 3 [37]. The problem is that while the Raman experiments on exfoliated samples estimate a huge enhancement of T CDW up to 145 K [35,36], STM experiments determine that dimensionality does not affect T CDW as the CDW occurs between 25 and 45 K [37].
In this Letter, we present first-principles calculations of T CDW both in bulk and monolayer NbSe 2 . We determine that the intrinsic CDW transition temperature is barely affected by dimensionality. The theoretically calculated phonon spectrum and T CDW in the bulk are in good agreement with inelastic x-ray experiments. Since the value obtained for the monolayer in a completely comparable calculation is very similar, it is confirmed that bulk and suspended monolayer NbSe 2 are expected to have a similar CDW transition temperature as suggested by the STM experiments. Our study also demonstrates that, when anharmonicity is fully taken into account, the contribution to the melting of the CDW order given by the electronic thermal fluctuations is totally irrelevant compared to the contribution of the ionic thermal fluctuations.
In Fig. 2 we show the (0 K) harmonic phonon spectra calculated for the bulk and the monolayer (obtained with the Methfessel and Paxton [38] cold smearing technique, see the Supplemental Material [25]). As it has been already pointed out [32,39,40], in both cases the harmonic phonons show many unstable modes. Following the displacement pattern of any of them, the ionic potential energy surface VðRÞ is lowered. Even if in both cases the LA mode is unstable close to q CDW ¼ 2=3ΓM, there are some differences in the harmonic phonon spectra. The most unstable mode along ΓM in the monolayer is shifted to smaller momentum with respect to the bulk, around 0.56ΓM, in agreement with previous calculations [39]. Remarkably, there is no instability at the M point in the monolayer, but a deep instability emerges along MK. Even if it has been argued that the CDW spatial modulation can be inferred from the q point where the deepest instability occurs in the (0 K) harmonic calculation [39,40], this can only be understood as a first hint and it may yield to a wrong interpretation. As a matter of fact, the modulation of the CDW can only be determined theoretically by calculating the phonon spectra as a function of temperature and seeing at which q a phonon becomes unstable on cooling.
Calculating temperature-dependent phonon frequencies in systems that undergo a second-order structural phase  [25]). The gray areas represent imaginary phonon frequencies, which are given with negative values.
FIG. 1. Crystal structure of bulk NbSe 2 in the 2H polytype. The unit cell is depicted, which is formed by two nonequivalent layers. The crystal structure of the monolayer is marked in the figure, which only contains one single layer. The lattice paremeters used are given in the Supplemental Material [25]. In the right panel the Brillouin zone is shown, which is restricted to the ΓMK plane in the monolayer. transition upon cooling is not a trivial task. A first, elementary, approach is to perform a temperature-dependent analysis within the harmonic approximation [41][42][43]. We perform these calculations within Mermin's grandcanonical extension to density-functional theory (DFT) [44]. We use the Fermi-Dirac distribution function to obtain the finite-temperature occupation of the Khon-Sham states [45], and we adopt the natural approximation of using the ground-state exchange-correlation functional also at finite temperature [46,47]. In this way, the finitetemperature electronic free energy F el ðRÞ ¼ E el ðRÞ− TS el ðRÞ, as a function of the ionic configuration R, is obtained [with E el ðRÞ we denote the average electronic energy at the considered temperature, which coincides with the electronic ground state energy VðRÞ in the zerotemperature limit]. The Hessian (divided by the square root of the masses) of the electronic free energy around its minimum gives the finite-temperature generalization of the standard 0 K harmonic dynamical matrix [for which the Hessian of VðRÞ around its minimum is considered]. An imaginary phonon thus corresponds to an atomic displacement pattern that lowers the electronic free energy and drives a displacive second-order phase transition.
It is worthwhile to stress that, at this level, thermal (and quantum) fluctuations of ions are not taken into account, thus only thermal (and quantum) fluctuations of electrons can play a role in the melting of the CDW order. As shown in Figs. 2 and 4, the finite-temperature electronic contribution to the free energy melts the CDWat around 872 K for the bulk and at around 917 K for the monolayer, which are both very far from the T CDW s reported experimentally. On the other hand, it is interesting to observe that the differences in the q order of the instabilities between bulk and monolayer found in the harmonic results at 0 K become irrelevant: the temperature-dependent harmonic analysis shows that, in agreement with experiments, in both bulk and monolayer the instability appears close to q CDW ¼ 2=3ΓM.
The harmonic results show a very weak dimensionality dependence in both q CDW and T CDW (the monolayer and bulk temperature differ by 5%). However, while the ordering vector q CDW is in good agreement with experiments, the poor agreement of the calculated T CDW s makes it very difficult to draw any reliable conclusion at the harmonic level, and suggests that anharmonicity could be the key ingredient in the CDW melting, as it happens in other chalcogenides [13,[48][49][50][51][52]. Since the high-temperature undistorted phase is not a minimum of the Born-Oppenheimer potential VðRÞ, anharmonicity cannot be included through standard perturbative approaches on top of the harmonic result. We overcome this problem by making use of the stochastic self-consistent harmonic approximation (SSCHA) [53][54][55], a quantum free energy variational method that allows us to perform full anharmonic nonperturbative calculations. Within the Born-Oppenheimer approximation, if the electronic temperature fluctuations are neglected, the nuclei Hamiltonian is given by K þ VðRÞ, where K is the kinetic energy operator. With the SSCHA we evaluate the free energy F ¼ E − TS of the system by minimizing the free energy functional where ρ R;Φ is the ionic density matrix of an auxiliary harmonic potential parametrized with average ionic positions R and effective force constants Φ (related to the amplitude of the ionic fluctuations around R), h□i ρ R;Φ ¼ Tr½ρ R;Φ □, and S ion ½ρ R;Φ ¼ −hln ρ R;Φ i ρ R;Φ is the ionic entropy. Therefore, at this level, the system's entropy S includes only the ionic contribution S ion . The application of the SSCHA requires the calculation of atomic forces on supercells. We used 6 × 6 × 1 supercells for both the bulk and the monolayer, which are commensurate with q CDW ¼ 2=3ΓM, and the forces were calculated on these supercells using DFT (see the Supplemental Material [25] for further details on the calculations). In Fig. 3 we show the phonon spectra along a path, computed by Fourier interpolating the SSCHA free energy Hessian obtained from the 6 × 6 × 1 supercell calculations. The anharmonic spectra display a huge anharmonic renormalization for the low-energy modes, with a remarkable temperature dependence. This highlights the relevant role of ionic fluctuations, which strongly feel the anharmonic part of the potential, in the CDW order melting. In the bulk case, the obtained phonon spectra at 250 K is in very good agreement with inelastic x-ray experiments [32]. At 50 K our calculations are in rather good agreement with experiments as well, but underestimate the frequency of the LA mode at q CDW ¼ 2=3ΓM, which is still imaginary. The LA mode stabilizes between 50 and 100 K in our calculations. Similarly, this mode also stabilizes in a similar temperature range in the monolayer. Interestingly, the deep instability along MK in the monolayer is washed out by anharmonicity at all temperatures. Now we focus our attention on the determination of the CDW modulation and transition temperature in this framework. The anharmonic phonon spectra in Fig. 3 are obtained by Fourier interpolating the SSCHA free energy Hessian obtained in the 6 × 6 × 1 supercell. In this Fourier interpolated spectra, the CDW instability occurs close to q CDW ¼ 2=3ΓM in both bulk and monolayer, thus driving a 3 × 3 CDW ordering in agreement with experiments. However, calculations in a larger 9 × 9 × 1 supercell for the monolayer show that the Fourier interpolated spectra are not converged with the supercell size over the whole Brillouin zone (see the Supplemental Material [25]). As a consequence, inferring the CDW ordering vector from the sole Fourier interpolated spectra would be risky. Nevertheless, the 3 × 3 ordering of the CDW is an experimental fact in both bulk and monolayer, consistent PHYSICAL REVIEW LETTERS 125, 106101 (2020) as well with the temperature-dependent harmonic calculation, where we are not limited to a coarse interpolation grid. Moreover, the phonon frequency of the LA mode at q CDW ¼ 2=3ΓM is well converged with the supercell size (as confirmed in the monolayer case by the calculations on the 9 × 9 × 1 supercell [25]). For this reason, we can confidently estimate T CDW within SSCHA, and readily compare it between the bulk and the monolayer, by studying the temperature dependence of the obtained anharmonic phonon frequency of the LA mode at q CDW ¼ 2=3ΓM.
As shown in Fig. 4, the square of the calculated anharmonic phonon frequency of the LA mode at q CDW ¼ 2=3ΓM shows a temperature dependence in agreement with the experimental trend [32]. The frequencies are slightly underestimated and, consequently, the value of the theoretical T CDW is around ∼59 K, close to the experimental value of 33 K, but slightly overestimated due to the presence of systematic temperature-independent DFT-related errors. In the monolayer the frequency of the LA mode at q CDW , as well as its temperature dependence, is practically on top of the bulk result. The CDW temperature in the monolayer is consequently very close to the bulk result, ∼73 K. Considering that the SSCHA calculations in the bulk and in the monolayer are performed with the same supercell, with consistent DFT parameters, the comparison between the results is perfectly justified. We can thus conclude that the CDW temperature in NbSe 2 is weakly dependent on the dimensionality, supporting the results obtained with STM experiments [37].
The weak dimensionality dependence of the CDW in NbSe 2 might suggest a weak interlayer interaction. However, the large energy difference between the two softened phonon modes in the bulk at q CDW , which are imaginary in the 0 K harmonic calculations and have a similar but not identical distortion pattern for each layer (see the Supplemental Material [25]), shows that there is a non-negligible interlayer interaction. Therefore, the electronic screening, the electron-phonon coupling, and the electronic and ionic fluctuations are expected to play a different role in monolayer and bulk, and it is the complicate interplay between these effects that yields a very similar CDW phonon branch regardless of the thickness.
A CDW phase, as any order in condensed matter, melts with increasing temperature due to entropy or, in other words, fluctuations. The presented SSCHA calculations  included the entropic effects coming from the ionic fluctuations only, with the electrons kept at 0 K. The good agreement between the experiments and the T CDW obtained with the SSCHA, compared to the extremely high T CDW found at harmonic level (which, on the contrary, included only the electronic fluctuations), suggests that the ionic fluctuations play the dominant role in suppressing the CDW order. However, since in bulk and monolayer NbSe 2 the electronic density of states at the Fermi level is quite sizeable [39,56], mainly due to localized d states, the electronic contribution to entropy is not guaranteed to be negligible when anharmonicty is taken into account. In order to carefully estimate the electronic contribution to the fluctuations also at the anharmonic level, so as to give a quantitative comparison between the ionic and electronic thermal effects, we perform SSCHA calculations including finite-temperature electronic fluctuations. We consider the Hamiltonian K þ F el ðRÞ for the ionic system (if electrons have finite temperature, in the adiabatic approximation forces and equilibrium position of the ions are ruled by the electronic free energy), and we minimize the functional where S½ρ R;Φ ¼ hS el ðRÞi ρ R;Φ þ S ion ½ρ R;Φ . Therefore, the SSCHA estimation of the system's entropy now includes contributions from both electrons (averaged through the ionic distribution) and ions. For the monolayer at 200 K, the difference between the SSCHA phonon frequencies obtained with and without the finite-temperature electronic contribution is around 2.6%, which has a negligible impact on the estimation of T CDW . At lower (electronic) temperatures, since the harmonic results remain identical (see the Supplemental Material [25]), the corrections to the SSCHA frequencies coming from the inclusion of finite-temperature electronic effects are not expected to be larger than that. Therefore, this test confirms that the ionic fluctuations totally overshadow the electronic fluctuations when anharmonicity is taken into account.
In conclusion, making use of first-principles calculations within the stochastic self-consistent harmonic approximation, we show that anharmonicity melts the CDW in both bulk and monolayer NbSe 2 , and that ionic thermal fluctuations predominate over electronic thermal effects. Moreover, we show that, in spite of interlayer interaction, the CDW transition temperature and spatial modulation weakly depend on thickness. Our result supports the STM measurements on single layers grown by MBE [37], but it questions the Raman results on exfoliated samples that estimate a large enhancement of T CDW in the monolayer [35,36]. It is to be seen whether the enhanced T CDW observed in Raman experiments is a consequence of the sample deterioration, e.g., oxidation, during the exfoliation process or if it is affected by the substrate. Indeed, similar theoretical calculations as those presented here have recently shown that strain or charge doping from the substrate can affect the CDW transition in other TMDs [48,49].