Unraveling the role of V-V dimer on the vibrational properties of VO₂ by first-principles simulations and Raman spectroscopic analysis

Vanadium oxide (VO₂) has been a focus of intense research since its discovery in 1959Morin (1959). Besides academic interest, it underpins a plethora of applications including gas sensorsStrelcov et al. (2009), window coatingsManning et al. (2004), resistive random access memory (RRAM) devices Kim et al. (2013a), electronic switchesShin et al. (1977); Lim et al. (2012); Yang et al. (2011), and so on. The versatile phase diagram of VO₂ showcase well known first order metal-insulator transition (MIT) at almost room temperature (340 K) and ambient pressure accompanied by a structural transition from a high-temperature rutile (R phase) to a low-temperature monoclinic (M₁ phase). The intricate concomitant nature of these two transitions has fueled a long-standing debate on the plausible mechanism for the MIT in VO₂. Several scenarios have been adopted including Peierls typeCavalleri et al. (2004); Wentzcovitch et al. (1994); Goodenough (1971), Mott-Hubbard-type and even combination of these twoBiermann et al. (2005); Haverkort et al. (2005); Weber et al. (2012)utilizing the role of lattice instabilities, electron-electron interactions, electron-phonon interactions, and so on. Though extensively studied, a consensus for a exact mechanism of MIT in VO₂ is still elusive.

The role of lattice vibrations in this MIT is still not completely understood. Studying lattice vibrational properties of the VO₂ turns out to be a non-trivial task from both the experimental and theoretical point of view. For example, the bottleneck of single-crystal inelastic neutron scattering (INS) is the incoherent vanadium neutron scattering cross-section. Even, some of the existing Raman spectroscopic studies of the structural phase transition (SPT) in epitaxial VO₂ films have been found to be inconsistent due to the large thicknesses of studied films (80- 100 nm) or even misinterpreted due to the strong signals from the substratePetrov et al. (2002); Shibuya and Sawa (2014); Okimura et al. (2014); Yang et al. (2016).

Besides these experimental challenges, achieving an exact theoretical description of the lattice dynamics of VO₂ has been a tougher nut to crack. Since lattice vibrational properties are strongly related with the structure of the material, first-principle based methods that take into account the much needed material specific information can be an ideal tool for this task. However, a reliable exchange functional which can correctly describe the vibrational properties of all the phases of VO₂ is not readily available and still an going topic of research. It has been found that the standard functional such as local density approximation (LDA) and generalized gradient approximations (GGA) can not describe the structural, electronic and magnetic features of the R and M₁ phases simultaneously Wentzcovitch et al. (1994); Zhu and Schwingenschlögl (2012). To ameliorate the inaccuracy of the band structure calculations based on the standard functional, several attempts have been already made to go beyond the LDA or GGA functional. For example, the modified Becke-Johnson (mBJ) exchange potential has been demonstrated to be very efficient for describing the MIT in VO₂Zhu and Schwingenschlögl (2012). Unfortunately, the Hellmann–Feynman force calculations are not possible using the mBJ ruling out the possibility of phonon calculations within this exchange approximation. The GW method and its variants successfully predict the electronic feature of the MITContinenza et al. (1999); Gatti et al. (2007); Sakuma et al. (2008). However, equivalent GW calculations for phonons are still not clear. The dynamical mean-field theory and its cluster extensions correctly describe the electronic properties of both the R and M₁ phases of VO₂Laad et al. (2006); Biermann et al. (2005). Cluster extension of DMFT such as dynamical cluster approximations (DCA) and its typical medium extension named as typical medium DCA (TMDCA) for disordered phononic systems has been recently developedMondal et al. (2017, 2019). Unfortunately, such DCA and TMDCA method for phononic systems are currently available for parameter-based model calculations. Moreover, hybrid functionals efficiently describe the electronic properties of VO₂Eyert (2011) and is found to be also successful in describing the structural properties of all the phases of VO₂Xu et al. (2017). However, phonon calculations using hybrid functionals turn out to be a very expensive computational task. First-principle based Quantum Monte Carlo based calculations have been also performed for VO₂Zheng and Wagner (2015). Again, computational expenses discard such possibility for applying the method in the lattice dynamics calculations.

Considering all these limitations, the GGA+U method Korotin et al. (2002); Kim et al. (2013b); Mellan et al. (2019); Budai et al. (2014) has been extensively used for phonon calculations of VO₂. However, the GGA+U method remains controversial in producing an accurate description of VO₂. Anisimov and co-authorsKorotin et al. (2002) find insulating nature of the R phase which strongly disagree with experimental findings. Mellan et al.Mellan et al. (2019) predict electronic structure of VO₂ that is consistent with experiment using non-spin polarised GGA+U calculations. But, U value within the LDA+U method has been treated so far as a free parameter and the results are also very sensitive with the choice of U value. For example, Budai et al.Budai et al. (2014) point out that the U value less than 3.4 eV is more appropriate in determining energy difference between the R and M₁ phase that could be more consistent with experiment, but such lower U value can produce more smaller bandgap. Hence, an accurate prediction of the U value rather using it as a free parameter or bare atomic value, may have a paramount importance for reliable estimation of the vibrational properties of VO₂. To best of our knowledge, calculations for determining U value of VO₂ via the constrained random phase approximation (cRPA) or linear response method are still missing. Another important thing is to include temperature effect and anharmonicity. Such calculationsBudai et al. (2014) have been already performed and can help to validate further theoretical development.

In this paper, we present a comprehensive study of the lattice dynamics and electronic structures using various DFT functionals and benchmark our calculations with experiments. By comparing with Raman spectroscopy, we investigate the interplay between strain effect and local environment (e.g. spectral shifts of the V-V and V-O phonon modes) in epitaxial VO₂ thin films deposited on MgF₂ (001) and (110) substrates. We show that our calculations are in a good agreement with the observed Raman shifts of the characteristic phonon modes induced by the epitaxial strains imposed on VO₂ thin films. Our work demonstrates that Raman spectroscopy could serve as a reliable way to detect the subtle change of V-V dimer in the strained VO₂.

We have organized this paper as follows: we give details of our first-principles simulations in section II. We provide information of our experimental set up in section III. We present and discuss our results in section IV. Finally , we conclude our discussions in section V.

We consider a plane-wave basis and the projector-augmented- wave (PAW) method as implemented in the Vienna Ab initio Simulation Package (VASP)Kresse and Furthmüller (1996, 1996). We have used several different exchange correlation functional including the LDA, the PBE GGAPerdew et al. (1996), and the SCANSun et al. (2015). We note that the SCAN has been implemented for self-consistent calculations in the VASP and such self-consistent implementation is still not available in other first-principles packages like WIEN2K. We choose PBE GGA pseudopotential in our calculations using SCAN. For electronic structure calculations, we use the plane-wave cutoff energy as 500 eV. We also verify that there is no appreciable changes in the results using a cutoff energy of 600 eV. For all Brillouin-zone sampling, $\Gamma$-centered k-point grids is considered as $16\times 16\times 16$. The self-consistent calculations are considered to be converged by considering $10^{-4}$ eV between successive iterations, and the convergence of structural relaxation is decided by the total energy difference between two successive ionic steps as $10^{-3}$ eV. We have also cross-checked our electronic structure calculations by using mBJ calculations as available in WIEN2K. The results obtained from the WIEN2K are not shown here. For the phonon calculations, we have employed the pseudopotential band method and the supercell approach that are implemented in VASPKresse and Furthmüller (1996, 1996) and PHONOPYTogo and Tanaka (2015). As a standard supercell based phonon calculations, force-constants are computed by means of the Hellmann-Feynman theoremParlinski et al. (1997).

High quality 10-nm-thick VO₂ films were grown on rutile MgF₂ (001) and (110) single crystal substrates by reactive molecular-beam epitaxy (MBE) via codeposition method under a distilled ozone background pressure at the PARADIM Thin Film Growth Facility at Cornell university.Paik et al. (2015) Raman microscopy analysis was conducted at the Center for Nanoscale Materials at the Argonne National Laboratory. The spectra were recorded using a Raman microscope (inVia Reflex, Renishaw, Inc.) with spectral resolution of 0.5 cm^-1 using 532 nm excitation from a diode pumped solid state laser (RL532C50, Renishaw, Inc.). Samples were held in a nitrogen-purged temperature-controlled stage (THMS600/TMS94/LNP94, Linkam Scientific Instruments Ltd). Excitation and collection of scattered light occurred through a 50x objective (Leica, NA=0.50). The laser power was set to 0.64 mW to exclude local heating effect. Collected spectra are typically consisted of the average of 20 integrations, where each integration was collected for 1 min.

We begin our discussions by understanding the well-known experimental crystal structure of VO₂Eyert (2002) as described in Fig. 1. The unit cell contains 6 atoms for the R phase and 12 atoms for the M₁ phase in our considered experimental structure. As shown in Fig. 1(a), the V-V atoms form a periodic V chain with fixed V-V bond distance as 2.86 $\AA$ in the rutile structure, whereas there are significance differences in the arrangement of the V atoms along the rutile $C_{R}$ axis in the M₁ phase, as demonstrated in Fig. 1(b). In this M1 phase, the V atoms are arranged such a way that they form dimmer alternatively and tilt along the $C_{R}$ axis, which leads to doubling of the unit cell volume of that for the R phase with V-V distance as 2.62 $\AA$(bonding) and 3.16 $\AA$ (anti-bonding). With this understanding of available experimental crystal structure, we perform structure optimization of both the R and M1 phases of VO₂.

For structure optimization, we consider the ionic optimization process, in which the lattice parameters are fixed and atoms are allowed to move to search for minimum energy configurations. First, we optimize the crystal structure of the R phase using different functional and compare with experimentEyert (2002) as well as previously reported valuesXu et al. (2017) as given in Table 1. For SCAN meta-GGA functional, the V-V pair distance changes from the experimental value 2.86 $\AA$ to 2.76 $\AA$, whereas it is 2.77 $\AA$ for the GGA functional. We note that there is no significant changes in the results obtained from the SCAN meta-GGA and GGA functional. However, this is certainly not true for the case of M₁ phase. Our calculations suggest that the optimization process is extremely sensitive to the functional used in the calculations as summarized in Table 2. We find that the LDA functional completely fails to optimize the M1 structure. We observe that the long V-V dimer and short V-V dimer distance becomes almost equal as 2.72 $\AA$ which imply that the M1 structure reduces to the rutile like structure through the LDA optimization process. Hence, like electronic properties, our calculations reveal that the LDA is not applicable to determine the structural properties of the M1 phase. The situation becomes better when the GGA is applied in our calculations. In the optimized M1 structure based on the GGA functional, the long V-V dimer distance decreases to 3.0 $\AA$ from the experimental value $3.16\AA$ and short V-V dimer separation decreases to $2.59\AA$ from the experimental value $2.62\AA$. Such drastic change in the dimmer distance clearly indicates that the M1 structure again tends to go to the rutile phase through the optimization process. As has been already foundXu et al. (2017), this failure of the LDA and GGA can be cured by the Hybrid functional. To overcome this barrier, we have utilized SCAN meta-GGA functional for optimizing the M₁ structure. As summarized in Table 2, in the SCAN optimized M1 structure, the long V-V dimer nature survive as $3.16\AA$ which is very close to the value $3.14$ as found in hybrid calculationsXu et al. (2017) and exact as $3.16\AA$ as found in experimentEyert (2002). The short V-V dimer using SCAN functional is found to be $2.47\AA$ which is also to $2.44\AA$ as found from the hybrid functional calculations. Thus, like hybrid functional, the optimized crystal structure using SCAN has all the expected features of the experimentally derived structure: All V-V chains of M1 are found to be dimerized and canted along the rutile axis, whereas all V-V chains form straight undimerized one in the rutile phase.

Here we want to investigate the effect of the structural optimization on the electronic structure, in order to draw an accurate band theory picture of the M1 phase. Fig. 2 presents the electronic density of states (DOS) based on our optimized structures. To treat the comparison on equal footing, we use SCAN functional in the calculation of electronic DOS from the structures that are obtained from different functional. As shown in the top panel of Fig. 2, there is no gap in the DOS corresponding to the LDA optimized structure. However, a pseudogap feature tends to emerge in the DOS of the GGA optimized structure as shown in the middle panel of the Fig. 2. Such small gap in the DOS with the GGA optimized structure was also previously reported and our calculations are consistent with these previous calculationsKim et al. (2013b); Budai et al. (2014). Interestingly, we find a fully gapped spectrum with the SCAN optimized structure as displayed in the bottom panel of Fig. 2.

Next, we study the effect of V-V dimerization on the electronic properties of M₁ phase. Fig. 3 displays the variation of the density of states as a function of various values of short V-V dimer separation. As it is clearly seen from the top to bottom of this figure, there is no gap in the density of states for the V-V dimer distance as $2.59\AA$. Interestingly, a finite gap emerges as the V-V dimer separation is varied from $2.59\AA$ to $2.48\AA$. Thus, our calculations suggest that a gap can be opened up for the M₁ phase depending only on the structure property without considering any Hubbard on-site coulomb interactions. From our calculations, it is certainly clear that the role of Peierl’s distortion cannot be neglected in determining the insulating nature of the M₁ phase. Such observations strongly suggest that the M₁ is not pure a conventional Mott insulator.

Besides successfully describing $M_{1}$ phase, determining the correct electronic structure of VO₂ including both rutile and M₁ phases using a single-functional is a considerable difficult task from band-theory point of view. To test the performance of SCAN against other functional in describing both the rutile and M1 phases of VO₂, we calculate the density of states as presented in Fig. 4. From the optimized structure, SCAN based calculations predict the high temperature rutile phase as a metal and the low temperature M1 phase as an insulator. These findings are consistent with experiments. We note that previous calculations based different functional fail to correctly describe both the phases of VO₂ simultaneously. For example, the LDA+U method find insulating nature of the rutile phase of VO₂Korotin et al. (2002); Liu et al. (2010). However, mBJ and hybrid functional as we discussed before can successfully describe both VO₂ phases. As we have already mentioned, mBJ cannot be utilized for calculating the forces and hybrid functional is computationally very expensive. Hence, they are not highly suitable particularly for phonon calculations of VO₂.

After establishing the correctness of our optimized structure and producing consistent electronic properties with the SCAN optimized structure, now we are ready to discuss our phonon results. Fig. 5 shows the phonon dispersion curves of R-VO₂ using SCAN meta-GGA functional. We have also performed same calculations using GGA functional. We find that the phonon dispersion for the rutile phase using the SCAN does not deviate significantly from the results using the GGA functional and observe the phonon softening instabilities in both the SCAN and GGA calculations. Such phonon softening of the rutile phase has already been reported before Kim et al. (2013b). We note that zero temperature DFT calculations within harmonic approximations fail to predict stable rutile phase which is also expected since the rutile phase is stable only at high temperature. We obtain phonon softening at $\mathbf{q}=\mathbf{\Gamma},\mathbf{M}$ and $\mathbf{X}$ which are excellent agreement with previous GGA based calculationsKim et al. (2013b). Existing studies Brews (1970); Hearn (1972); McWhan et al. (1974); Terauchi and Cohen (1978); Gervais and Kress (1985); Woodley (2008) suggest that the phonon softening at $\mathbf{q=R}$ leads structural transition from the high temperature R structure to low temperature $M_{1}$ structure. However, we have not observed any phonon softening at $\mathbf{q}=\mathbf{R},\mathbf{A}$ and $\mathbf{Z}$. The phonon softening at $\mathbf{q}=R$ is observed in the previous LDA+U calculationsKim et al. (2013b). It has been argued electron correlation is the driving force for such instability. However, recent experimental studyBudai et al. (2014) shows that experimental dispersion curve with no negative frequencies of the rutile phase can be reproduced by including temperature and anharmonicity in the theoretical calculations. However, the role of structural instability in driving the MIT is still not well understood.

To test the limitations of harmonic approximations and zero temperature calculations, we calculate phonon density of states of rutile phase as shown in Fig. 7. Followed by the dispersion, there are negative frequencies in the low energy region of the spectrum which is not shown here. Interestingly, we obtain a good agreement with the experimental phonon density of states. We note that the experimental phonon spectrum has been collected at temperature T = 381 K, whereas our calculations are limited to T = 0 K within harmonic approximations. Despite these limitations, we are able to capture qualitatively peak features of the experimental spectrum around 18 meV, 33 meV, 45 meV, and 72 meV. We also observe that the height of the peak around 18 meV and 33 meV is more compared to the experimental one. This could be consequence of temperature effect. In addition, we also compare our results against previous calculationsBudai et al. (2014). As we expected, we notice that the height of the peak decreases significantly particularly around 18 meV and 33 meV as temperature increases from 0 to 425 K. We also observe that our calculations are able to capture some features of the spectra that agree better with the experimental spectrum compared to previous theoretical calculationsBudai et al. (2014). For example, the peak around 55 meV corresponding to the experimental DOS is shifted to around 65 meV in the DOS obtained from the previous calculationsBudai et al. (2014). Also, previous calculations find DOS at the high energy region remain around 83 meV, whereas the bandwidth of the experimental spectrum is larger. We note that the bandwidth of our calculated spectrum is much closer to the band width of the experimental DOS compared to previous finite temperature calculations. Although the high temperature rutile phase is believed to be anharmonic, the low temperature M1 phase is found to be highly harmonic. So, we expect our theoretical framework is more suitable for the M1 phase, which is our next topic of discussion.

We have computed phonon dispersion of the M₁ phase as displayed in Fig. 6. Our calculated dispersion curve for the M1 phase is in good agreement with previous calculationsLee et al. (2017). As shown in Fig. 6, unlike rutile phase, there is no phonon softening which is in agreement with experimental observationsBudai et al. (2014). It confirms that this M1 phase is the stable structure of VO₂ at low temperature. Such structural information can be directly probed in Raman spectra. To understand that, a careful analysis of the phonon density of states can be good a starting point.

The phonon density of states of M₁ phase is presented in Fig. 8. Followed by dispersion, there is no negative frequency in the low frequency region of our calculated phonon density of states. We compare our calculated phonon density of states at T =0 K with the phonon density states corresponding to T=298 K obtained from experiment using neutron spectrometer as reported in RefBudai et al. (2014). We find our zero temperature calculations using SCAN functional successfully capture the experimentally observed phonon density of states at room temperature. To check the performance of SCAN against other functional, we also compare our results with previous LDA+U calculationsBudai et al. (2014) as shown in the lower panel of Fig. 8. Comparison between the results obtained from SCAN and LDA+U clearly shows that the results obtained from SCAN serves better than the LDA+U. In particular, the peak features in the low frequency region around 20 meV is more intense compared to the LDA+U results as observed in experiments.

We take a more closer look on the phonon density of states and study the effect of the V-V dimer on it. The role of V-V dimerization on the insulating behavior of the M1 phase is well understood now, whereas its effect on the lattice vibrations is still not extensively studied. Here, we aim to explore such physics. As we have already discussed that the V-V dimer distance deviates significantly from the experimental one, we compare phonon density of states of M₁ phase using both optimized and experimental structure as shown in Fig. 9. We observe that there are large number of negative frequencies in the low frequency region of the phonon density of states derived from experimental structure as displayed in upper panel of Fig. 9, whereas there are no such negative frequencies in the low frequency region of the phonon density of states considering optimized structure as shown in the lower panel of Fig. 9. The absence of such negative frequencies in the low frequency region of the phonon density of states corresponding to optimized structure is well expected. Moreover, there are also other noticeable differences. There is a peak in the low frequency region around 190 cm^-1 and a peak in the high frequency region around 609 cm ^-1 in the phonon density of states obtained from with out optimizing the structure. The low frequency peak is mainly contributed from V-V pair and high frequency peak is originated from V-O pair which can be understood from Fig. 10. The positions of the peak and nature of the vibrations are in excellent agreement with the observed Raman spectra. However, the position of the peaks are shifted in the phonon density of states corresponding to optimized structure compared to without optimization one. We attribute this to the drastic change in the short V-V dimer distance that occurred in the optimization process Table 2.

Modulation of the V-V dimer length can be accomplished with strain engineering. It has previously been shown in VO₂/TiO₂ epitaxial thin films that tensile strain on the c_R-axis shifts the MIT mechanism towards bulk-like characteristics largely driven by the Peierls physics while elongation of the c_R-axis induces a larger influence from electron correlation or Mott physics.Quackenbush et al. (2016); Paez et al. (2020); Lee et al. (2019). However, due to the strong signal from the TiO₂ in Raman spectroscopy, the strain effect on the phonon modes has been misinterpreted.Yang et al. (2016) Moreover, although strain has been accomplished on numerous substrates such as ZnOWong et al. (2013), Al₂O₃Fan et al. (2013), LSATLiu et al. (2018), and RuO₂Fischer et al. (2020), the formation of rotational domains may complicate fundamental structural analysis in the investigation of the MIT. Like TiO₂, MgF₂ is isostrucutral to VO₂ in its high temperature phase, and exhibits weaker Raman response than TiO₂ substrateShibuya et al. (2013) and thus we selected it as the substrate for the epitaxial growth of VO₂ thin films. Fabrication of VO₂ on MgF₂ (001) induces the tensile strain of the c_R-axis while MgF₂ (110) films elongate the c_R-axis. Indeed, using Raman spectroscopy we observed spectral shift of the characteristic Raman peaks which represent V-V and V-O phonon modes. As shown in the right panel of Fig. 11 and Fig. 12, low frequency peak is shifted from 196 cm ^-1 to 197 cm ^-1, whereas high frequency peak is shifted from 609 cm ${-1}$ to 618 cm ^-1 with the application of strain. Such spectral displacement is shown in our calculated phonon density of states by manipulating V-V dimmer. As presented in the left panel of Fig. 11 and Fig. 12, the low frequency peak has been shifted from 190 $cm^{-1}$ to 230 cm^-1 as short dimer V-V distance is varied from 2.62 $\AA$ to 2.47 $\AA$. Similarly, the high frequency peak is deviated from 609 cm ^-1 to 628 cm ^-1. Thus, our calculations qualitatively explain the effect of strain on the lattice vibrations as observed in our Raman experiments.

In this paper, we have computed vibrational properties of the R and M₁ phases of VO₂ using first-principles calculations. A crucial step for calculating such vibrational properties is the structure optimization which is found to be non trivial for the M₁-VO₂. We have employed SCAN meta-GGA functional and successfully optimized the M₁ structure. The performance of the SCAN has been compared with various other functional in optimizing the M₁ structure. Particularly, the SCAN has been found to be as accurate as the Hybrid functional. The SCAN optimized structure has been shown to produce consistent electronic properties including both the R and M₁ phases of VO₂. We have shown that the R phase is a metal, whereas the M₁ phase is found to be an insulator with finite band gap value. Moreover, we numerically demonstrate that a gap can be opened up as the V-V dimer length is varied without incorporating electronic correlations. Hence, our results strongly suggest that the M₁ phase is not a conventional Mott insulator given the sensitivity to the dimer distance to its electronic phase. Thus, we have validated the SCAN for phonon calculations of VO₂. The phonon softening related with the R phase and the phonon stiffening of the M₁ phase explain the structural phase transition from the R to M₁ as temperature is lowered. We argue that such structural phase transition can play a vital role in the MIT of VO₂ as found in recent experiment. The V-V dimer has significant impact on both the electronic and vibrational properties of the M₁ phase. We show that the vibrational density of states significantly changes with the variation of the V-V dimer distance. Our calculated vibrational spectra are corroborated by previous neutron scattering experiment as well as our Raman experiment. Moreover, our first principles based calculations of bulk VO₂ enable us to gain insight into the structural properties of strained VO₂ as found in our Raman experiment. Our work also suggests that the SCAN can be utilized for investigating phonon properties of related materials and such work is in progress.