Correlating the Energetics and Atomic Motions of the Metal-Insulator Transition of M1 Vanadium Dioxide

Materials that undergo reversible metal-insulator transitions are obvious candidates for new generations of devices. For such potential to be realised, the underlying microscopic mechanisms of such transitions must be fully determined. In this work we probe the correlation between the energy landscape and electronic structure of the metal-insulator transition of vanadium dioxide and the atomic motions occurring using first principles calculations and high resolution X-ray diffraction. Calculations find an energy barrier between the high and low temperature phases corresponding to contraction followed by expansion of the distances between vanadium atoms on neighbouring sub-lattices. X-ray diffraction reveals anisotropic strain broadening in the low temperature structure’s crystal planes, however only for those with spacings affected by this compression/expansion. GW calculations reveal that traversing this barrier destabilises the bonding/anti-bonding splitting of the low temperature phase. This precise atomic description of the origin of the energy barrier separating the two structures will facilitate more precise control over the transition characteristics for new applications and devices.

Beginning with the work of Cavalleri et al. 5 pump-probe measurements have shed considerable light on the lattice dynamics occurring across the transition. A structural "bottleneck" associated with the phonon connecting the monoclinic and tetragonal structures was observed upon hole photo-doping 26 , suggesting that the insulating phase depends significantly on the lattice potential, indicating band-like character. Kim et al. 27 used pump-probe measurements in conjunction with X-ray diffraction and found that the sharp resonance corresponding to the monoclinic A g peak disappears at the transition with lower energy and less intense tetragonal B g resonances replacing them. Wall et al. 28 also used pumping to modify the lattice local potential and examine its effects on the coherent phonon spectrum of VO 2 , as an example of the general applicability of the use of pumping to induce a change in lattice potential, which can be used to study relaxation processes. However, a theoretical description of the interplay between the lattice potential and the atomic and electronic structure has proven elusive.
The computational study of Zheng and Wagner 29 utilised the Quantum Monte Carlo approach to show that the MIT is a direct consequence of the change in structure, that is the monoclinic structure is insulating, and the tetragonal form exhibits metallic behaviour. While sounding trivial, there has been some conjecture over the coincidence of the structural and electronic phase transitions 27,30 , which Zheng and Wagner, and also this work resolve. Chen et al. 31 explored the properties of the parameter space spanned by the β angle and the tetragonal c-axis using the DFT+ U approach 32 , and suggested that changing orbital occupancy is initially responsible for opening the band gap as a result of dimerisation, which is widened by a subsequent increase in the antiferroelectric distortion.
Thus what is missing from the literature as it currently stands is an exploration of the energy landscape of the structural phase transition with respect to the metal-insulator transition in terms of exactly what constitutes the separation between the two structures. The intent of this work is to utilize a comprehensive computational approach to determining the processes occurring during the metal-insulator and structural phase transitions, and to combine it with experimental data to confirm the predictions of our calculations. Specifically, the outstanding questions we seek to answer are: i) literature data suggests a latent heat of ~40 J/g for the transition 33 , to what does this energy barrier correspond? Which particular atomic motions give rise to this barrier, ii) if a minimum energy path can be mapped between the structures, and the aforementioned atomic displacements determined, what are the effects of these displacements on the electronic structure? Are the structural phase transitions and metal-insulator transitions necessarily coincident as suggested by Zheng and Wagner 29 ?
We start by computing the lowest energy path between the structures using the nudged elastic band technique 34 and density functional theory to determine this energy landscape. The DFT data reveal that in order for the structural transition to occur, the inter-vanadium spacing along the [110] or [110] directions must be compressed, generating electronic repulsion and thus an energy barrier. High resolution X-ray diffraction measurements reveal anisotropic strain related to the atomic spacing in these directions in the monoclinic structure, which is not present in the tetragonal form. Frequency-dependent GW calculations reveal that the top of the barrier corresponds to the opening of the gap due to bonding/anti-bonding splitting as the vanadium atoms dimerise. The data indicate that the most efficient modulations of the transition temperature involve stress input along [110] or [110] of the tetragonal structure or the [011] or [011] directions of the monoclinic structure, consistent with the action of doping with tungsten 35 .

Results
Structural Rearrangements. The relevant structural characteristics of tetragonal and M 1 VO 2 are presented in Fig. 1. Figure 1a,b compare the view of the tetragonal structure down 〈 001〉 T (the subscript T or M refers to tetragonal or monoclinic respectively) with the view of the M 1 structure down 〈 100〉 M . The comparison illustrates that the structural rearrangements occurring in the transition from the tetragonal to the monoclinic structure orthogonal to the monoclinic a-axis can be summarized as an alternating off-set of the vanadium atoms from the centers of the oxygen octahedra. This off-set occurs in the long axis of the octahedra. The (110) T and (011) M planes of the tetragonal and monoclinic structure are also indicated with black lines, which reveals that they correspond to equivalent atomic spacings in each structure (although due to a slight expansion of the tetragonal structure its diffraction peak manifests at slightly lower angle). Figure 1c,d present a comparison of the tetragonal and M 1 structures down 〈 010〉 (this axis is coincident for the tetragonal and monoclinic structures), which indicates that the changes occurring across the structural phase transition parallel to the monoclinic a-axis consist of the evenly spaced vanadium atoms (the "vanadium chains") of the tetragonal structure pairing up (the so-called Peierls pairing), forming an alternating long-short pattern of inter-vanadium spacing. The (101) T and (202) M planes are indicated in the tetragonal and monoclinic structures respectively, the (202) M plane has been shifted by one half of its spacing to illustrate that it is the equivalent distance in the monoclinic structure of the (101) T plane. Figure 1e illustrates the four characteristic distances of interest in this study which are orthogonal to the monoclinic a-axis. The V-O Apical Long V-O Apical Short distances describe the amount to which the vanadium atoms are off-set from the center of the octahedron; if the numbers are equal then the atom sits at the center of the oxygen octahedron. The V-V Corner Long and Short distances describe the two shortest distances between the vanadium atoms on neighboring chains, these distances lie parallel to the (101) T and (200) M , (201) M planes respectively. Figure 1f defines the V-V Chain Long (indicated by the letter "L") and short ("S") distances. These are the distances between the vanadium atoms in the chain which undergoes Peierls pairing. If both distances are equal, as in the tetragonal structure, the vanadium atoms are evenly spaced. As the tetragonal structure transforms into the monoclinic form, the atoms pair up and one of these distances decreases, while the other increases.
Nudged Elastic Band Calculations. The total free energies of the structures obtained along the minimum energy path between the monoclinic and tetragonal structures determined by the nudged elastic band method are plotted in Fig. 2. While the energies of the monoclinic and tetragonal structures are almost degenerate, there is a clear energy barrier between the two structures. This corresponds to an energy of 18.6 J/g, which compares favorably with the experimentally observed specific heat of the phase transition of ~40 J/g 33 . However, this result view down 〈 100〉 M 1 structure, the (011) planes are indicated by black lines, illustrating that they correspond to the same distance as the tetragonal (110) planes, the off-centre positioning of the vanadium atoms in the M 1 structure from the anti-ferroelectric twist is also apparent (c) view down 〈 010〉 of the tetragonal structure, the evenly spaced vanadium chains are visible, running parallel to the c-axis, and the (101) planes are indicated with black lines, (d) view down 〈 010〉 of the M 1 structure, with the Peierls paired vanadium chains visible running parallel to the a-axis and the (202) planes are marked with black lines. These are shifted by one half of a lattice spacing for better comparison to the tetragonal (101) planes, illustrating that both correspond to the same distance, (e) same perspective as (b) but with the "V-V Corner Long", "V-V Corner Short", "V-O Apical Long" and "V-O Apical Short" distances marked and (f) same perspective as (d) but with the "V-V Chain Short" and "V-V Chain Long" distances indicated by the letters "S" and "L" respectively.  Figure 3a,b indicate that the Peierls pairing distance increases continuously across the transition from monoclinic to tetragonal, while unsurprisingly, the corresponding long inter-vanadium distance decreases monotonically. This data simply expresses the fact that the evenly spaced vanadium atom chains running along the tetragonal c-axis (monoclinic a-axis) experience a Peierls distortion and adopt a long-short internuclear spacing configuration. The monotonicities of the plots do not suggest an origin for the energy barrier of Fig. 2. Figure 3c,d however tell a different story. Figure 3c indicates that the short distance between the central vanadium atoms and the corner vanadium atoms (see Fig. 1c) increases monotonically across the transition, however the longer distance (plotted in Fig. 3d) initially contracts, and then expands. Comparison of Figs 2 and 3 suggest that the peak of the total energy corresponds approximately to the minimum in the long V-V corner distance. Figure 3e,f illustrate the trends of the apical vanadium-oxygen distances, and the shorter distance again displays a monotonic increase across the metal-insulator transition, however the longer distance initially plateaus, before decreasing significantly.
Thus the only behavior occurring across the metal-insulator transition consistent with an energy barrier, i.e. an initial increase and subsequent decrease, is the compression and expansion of the long V-V corner distance. This indicates that the force needed to effect the transition between the structures is directed approximately along the diagonals of the unit cell, perpendicular to the vanadium chains. This corresponds to the [110] T and [110] T directions of the tetragonal structure. These directions describe the spacing of the {110} planes of the tetragonal structure, and the {011} planes of the monoclinic structure. Such an effect closely mirrors the observations of Pouget et al. 36 who found that inputting a uniaxial stress along the [110] T direction resulting in the appearance of the M 2 monoclinic form of vanadium dioxide. X-ray absorption also revealed that changes in this distance in tungsten-doped VO 2 correlated with the amount of tungsten doped into the lattice and therefore the degree to which the transition temperature was depressed 35 . Thus this direction seems to be of significance in the structural phase transition. Investigation of any changes in these spacings occurring may therefore confirm the prediction of the computational approach.
High Resolution X-ray Diffraction. To determine if there was any manifestation of the effects of this distortion in experimental data, high resolution X-ray diffraction was performed using the powder diffraction beamline at the Australian Synchrotron. Diffraction data of a sample of pure VO 2 were recorded above and below the structural phase transition temperature of ~340 K, and Fig. 4 illustrates the most significant properties. Figure 4a contrasts the (110) T and (011) M peaks respectively (which correspond to the same planes, (see Fig. 1a,b). The data shows clearly that the transition from tetragonal to monoclinic results in slight splitting and considerable broadening of the peaks corresponding to the inter-vanadium distances along the [110] T and [110] T directions.
However, this broadening is not uniform. Figure 4b illustrates the (101) T peak, and a triplet corresponding to the (202) M , (211) M , and (200) M peaks (from low to high angle). In this case, the (101) T and (202) M peaks correspond to the same distance, and despite a difference in amplitude, the peak shapes are very similar. Therefore, the data indicates that these distances do not experience any broadening as the structure transforms from tetragonal to monoclinic. Thus the data indicates that peaks of the form (0xx) or (0xy) experience significantly more broadening than other orientations. Such spacings describe distances with the same orientation as that of Fig. 3d; directed toward the neighboring vanadium chain. This disorder may be reconciled with the NEB data by taking into account the effects of defects and grain boundaries in the structure of the experimental sample. Figure 3d indicates that a distance is initially compressed, and subsequently extends. Figure 2 suggests that this compression costs energy. Thus, if structural defects are present which allow dissipation of this energy, the transformation to the monoclinic structure may be incomplete. Obviously individual grain boundaries will place limits on the extent of the propagation of this, and therefore it is possible that due to this, the strain broadening of the individual grains comprise a distribution which is anisotropic in nature.
To investigate the possible manifestation of this, anisotropic strain broadening parameters were extracted using the Stephens method 37 . Figure 5a plots the magnitudes of the S 022 and S 040 contributions to the broadening for the M 1 structure for five temperatures below the critical point, and contrasts them with the equivalent Tetragonal S 400 and S 220 contributions respectively (the mismatch in indices between the monoclinic and tetragonal parameters is due to the aforementioned different naming conventions of crystallographic axes in the cells, thus S 040M = S 400T and S 022M = S 220T ). The co-existence of the monoclinic and tetragonal structures near the critical  temperature creates issues for fitting, and thus the data of Fig. 5 is limited to those points near the transition temperature which exhibited the best fitting parameters. While the S 040M data is independent of temperature, and approximately equivalent in magnitude to the S 400T data, the S 022M data diverges as the temperature approaches the critical point. This contrasts sharply with the S 220T data which is almost un-correlated, as S 220T = − 1.4, which is very close to zero. This data therefore indicates that as the critical temperature is approached from below, the contributions to peak broadening from variations in the h and k spacing become increasingly correlated, and the magnitudes of the variances increase rapidly. In other words, the broadening observed contains components along both the crystallographic a and b axes. Figure 5b illustrates that the Stephens parameters corresponding to the only the b-axis (S 040 ), and those with no contribution at all from the b-axis are low in magnitude, and show no temperature dependencies. Figure 5c illustrates that while there is a correlation between the monoclinic a-and b-axes in the behaviour of the S 220M parameter, it is dwarfed by the behaviour of the parameters in which contributions from the b-and c-axes are present: S 121M and S 022M .
Figures 5a-c thus reveal that the diffraction data contains contributions to the broadening of the peaks which is anisotropic in nature, and that the most significant contributions are those of the S 022 and S 121 parameters, while the tetragonal structure shows no significant correlation in the equivalent parameters. This data is therefore in line with the data of Fig. 4, and Tables 1 and 2, and supports the hypothesis that the energy barrier between the structures corresponds to the compression and expansion of the characteristic distance of Fig. 3d.
The observed behaviour of the S 220 data is in some respects not surprising, as when plotted as a function of temperature, it is basically an un-normalized temperature correlation function of the variances of the a* and b* axes. The divergence of this correlation at the T c point is in line with critical behaviour expected at a phase transition, however, we do not attempt to explore this aspect in this work.   Electronic Structure. What remains to be determined however, is the effect of these structural rearrangements on the electronic structure. There is little question that the electronic structures of the tetragonal and monoclinic forms are metallic and insulating respectively 29 , however of interest in this study is the behaviour of the band structure across this structural phase transition, in order to determine whether the structural and electronic phase transitions are intrinsically related, or in fact merely coincident. The next section explores this in detail. The band structures of the monoclinic ground state, and structures "Step 1", " Step 2" and "Step 3" calculated using the GW approximation (black filled circles, fitted with blue splines) are presented in Fig. 6. The densities of states are plotted next to each band structure, on the same energy scale (blue filled curve). As expected, the monoclinic structure is insulating, with a band gap of ~0.70 eV, in excellent agreement with experiments (0.70 eV) 38 .
However as the structure transitions to the slightly more symmetric forms of Step 1 and Step 2 the densities of states indicate that the gap closes, and the structure becomes metallic. The corresponding band structures indicate that this occurs via two simultaneous mechanisms. The dispersions of the valence bands in the Γ → A direction suggests that in comparison to the band minima at Γ , the higher energy states near E F are shifting upwards, most significantly near A and D. At the same time, the conduction band minima at Γ shift downwards. In structure "Step 2" the conduction and valence bands overlap (this was determined by inspection of the charge density), and the indirect gap closes. The band structure of Fig. 6c, when compared with the total energy data of Fig. 2, indicates that the electronic band gap destabilises before the top of the energy barrier is crested.
We can gain a better idea of how the structural transitions are affecting these states by transforming them to real space charge densities and comparing them. Figure 7a,b present charge density isosurfaces of the valence and conduction band states at Γ in the (011) plane of the ground state monoclinic structure respectively, while Fig. 7c-d present charge density isosurfaces in the (011) plane of the valence and lowest energy conduction band states at the D point of the M 1 structure.
From comparison of Fig. 7a,b, it is obvious that the valence band state at Γ consists of shared charge density between the vanadium (grey spheres) atoms, while the conduction band state corresponds to isolated density on each vanadium atom. This suggests that the gap between the valence and conduction bands arises from bonding/ antibonding splitting. The same story is repeated at D. The valence band state contains density linking the Peierls paired vanadium atoms, while the conduction band state consists again of isolated density on each vanadium atom. The band structure data of Fig. 6 indicate that as the structure transitions away from the M 1 form and towards the tetragonal, the splitting between these states decreases considerably, with the eigenvalues at D crossing over by Step 3. This charge density, combined with the eigenvalues, and the inter-vanadium spacing data of Fig. 3 indicate that as the structural phase transition progresses, starting from the M 1 structure, the inter-vanadium distance of the Peierls pairs increases, destabilising bonding states with respect to anti-bonding states, narrowing the gap between the conduction and valence bands, until it disappears completely and the structure becomes metallic. Figure 8 plots the value of the local potential along a line segment connecting the Peierls paired vanadium atoms of the M 1 , "Step 3", "Step 6" and tetragonal structures, and from the data, it is obvious that as the inter-vanadium distance increases, the height of the potential barrier between the nuclei also increases. This increase will significantly affect the wavefunctions of the highest energy electrons which are obviously less tightly bound, resulting in less electron sharing between the paired atoms, consequently raising the energy of bonding configurations with respect to anti-bonding.   input of strain along the [110] T direction 36 . The chains consist of a linear, Pe These displacements result in the destabilisation of the bond/anti-bonding splitting of the monoclinic structure, due to the increase in the potential barrier between the Peierls paired vanadium atoms. This results in the conduction and valence bands overlapping; an insulator-metal transition.
If the energy barrier is indeed a consequence of the atomic motion of Fig. 3d, then attempts to manipulate the transition temperature which involve modulation of the stress or strain along the [110] T and [011] T directions would produce the most significant effects. This is consistent with X-ray absorption studies 35 which indicate that the depression of the transition temperature by tungsten-doping correlates with an increase in the V-V corner spacing. If the energy increase of Fig. 2 is a consequence of electronic repulsion, then increasing the V-V corner distances will lower the repulsion due to contraction as the increased internuclear separation will result in a lowering of the electrostatic potential between the atoms, reducing the barrier height.

Methods
Variable temperature Synchrotron X-ray Powder Diffraction was conducted at the Australian Synchrotron Powder Diffraction Beamline. Samples were sealed in 0.3 mm borosilicate glass capillaries. Prior to data collection, the wavelength was set at 0.82732 Å using a Si (111) double crystal monochromator. The exact wavelength was refined using the NIST 660b LaB 6 standard reference material. A Cyberstar hot-air blower was used to control the temperature to within 0.1-0.2 °C at each data collection temperature. Traces were recorded for 5 minutes at each of the two detector settings after the sample had reached the set point temperature and equilibrated for 10 minutes.
Quantitative Rietveld analysis was performed on the data using the Bruker TOPASTM V4.2 program to determine the weight percentage of phases present. Background signal was described using a Chebyshev polynomial linear interpolation function. A broad pseudo-Voight function was also used to model the background contribution form the capillary. Cell parameters, atom positions, (tightly constrained) isotropic thermal parameters, Gaussian and Lorentzian contributions to peak full widths at half maximum and scale factor were all refined.
The anisotropic broadening of the diffraction peaks observed in Fig. 4 was hypothesized to originate from one of two sources. Either there was some anisotropy in the crystallite shapes, leading to broadening of lines corresponding to directions in which fewer planes are stacked, or the crystal grains exhibit a distribution of residual strains originating from the transition from tetragonal to monoclinic. Thermal effects were deemed an unlikely origin, as refinements produced similar thermal parameters at all temperatures, and the anisotropic broadening observed is smaller at higher temperature. In addition, the Debye-Waller factor 39,40 tends to reduce the scattered intensity, however in comparisons between the equivalent monoclinic and tetragonal peaks, such as (011) M and (110) T the peaks integrate to the same total intensity. Employing Jarvinen's method 41 to account for anisotropic broadening led to rather poor fits in comparison to those generated by the Stephens method 37 for strain broadening, indicating that while some crystallite size ansiotropy may exist, the broadening is dominated by the strain distribution.
Strain analysis of the X-ray data was performed using the method developed by Stephens 37 , which is a phenomenological approach to determining the contributions to broadening induced by anisotropic variations in plane spacings. We repeat the central thesis of this approach here, but for a more complete treatment the reader is referred to the original work 37 .
The spacing of planes with Miller indices hkl is given by: Re-labeling the metric parameters {A, … , F} as {α i } and assuming that they have Gaussian distributions characterised by a covariance matrix , the variance of M hkl can be written: etc. can be re-written: This Γ A is combined with the usual parameters for Gaussian and Lorentizian line-widths to give expressions for anisotropically broadened line-shapes which are fitted to the experimental data, and the S hkl are extracted from the fit. The Gaussian and Lorentzian broadening parameters of Tables 1 and 2 were extracted from fits to the individual peaks presented in the data of Fig. 4.
Force calculations. The monoclinic 42 and tetragonal 43 structural parameters were input to DFT geometry relaxations using the VASP code 44  Perdew et al. 45 , on 6 × 6 × 6 and 8 × 8 × 6 Monkhorst-Pack 46 k-space grids. The structures were then relaxed to their respective ground states using Methfessel and Paxton smearing 47 and the conjugate gradient algorithm. Upon reaching the desired ground states, a 1 × 1 × 2 supercell of the tetragonal structure was constructed in order to have the same dimensions as the monoclinic form.
The Cartesian atomic positions of the tetragonal structure were then subtracted from those of the monoclinic structure, which generated vectors describing the movement of the atoms across the transition. Vectors describing the changes in unit cell dimensions were obtained in the same manner. These vectors were then divided such that 10 structures were generated, with the monoclinic structure being the first, and the tetragonal being the last and the intermediate structures are labelled "Step 1" to "Step 8". The elastic band technique 34 was then applied to these structures, in order to find the minimum energy path between them. The use of DFT to determine the total energies of Fig. 2 from the structures optimised by the elastic band method, rather than DFT+ U or hybrid functionals stems from the requirement to maintain a consistent Hamiltonian for the calculation of the energies along the minimum energy path as the energy landscape, by definition, will be Hamiltonian dependent.
Electronic Structure. The relaxed structural parameters were used as input to Density Functional Theory 44,48 calculations on 6 × 6 × 6 Monkhorst-Pack k-space grids, again using the Generalized Gradient Approximation (GGA) approach to exchange and correlation of Perdew et al. 45 with the Brillouin zone integration approach of Bloechl et al. 49 . Frequency-dependent GW calculations 50 were performed using the implementation of Shishkin and Kresse 51 in VASP 44 . The GW calculations were performed using a grid of 50 frequency points and an energy cutoff of 200 eV.