Local structure elucidation of tungsten-substituted vanadium dioxide (V1-x\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{1-x}$$\end{document}Wx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_x$$\end{document}O2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document})

Initially, vanadium dioxide seems to be an ideal first-order phase transition case study due to its deceptively simple structure and composition, but upon closer inspection there are nuances to the driving mechanism of the metal-insulator transition (MIT) that are still unexplained. In this study, a local structure analysis across a bulk powder tungsten-substitution series is utilized to tease out the nuances of this first-order phase transition. A comparison of the average structure to the local structure using synchrotron x-ray diffraction and total scattering pair-distribution function methods, respectively, is discussed as well as comparison to bright field transmission electron microscopy imaging through a similar temperature-series as the local structure characterization. Extended x-ray absorption fine structure fitting of thin film data across the substitution-series is also presented and compared to bulk. Machine learning technique, non-negative matrix factorization, is applied to analyze the total scattering data. The bulk MIT is probed through magnetic susceptibility as well as differential scanning calorimetry. The findings indicate the local transition temperature (Tc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c$$\end{document}) is less than the average Tc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c$$\end{document} supporting the Peierls-Mott MIT mechanism, and demonstrate that in bulk powder and thin-films, increasing tungsten-substitution instigates local V-oxidation through the phase pathway VO2→\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2\, \rightarrow$$\end{document} V6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_6$$\end{document}O13→\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{13} \, \rightarrow$$\end{document} V2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}O5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_5$$\end{document}.

Once this occurs and the structure is transitioned to the tetragonal phase, the a-direction expansion decreases by a factor of three. The slower lattice expansion after phase transformation is due to the 12.9% ionic radii difference between V (0.58 Å) and W (0.66 Å), assuming they are both in the 4+-oxidation state in an 6-coordinate environment 79 .
While in the monoclinic phase, the lattice expansion leads to a decrease in the c/a parameter and an increase in c/a, when in the tetragonal phase. The structural transition occurs when the c/a ratio reaches 0.62663 (4), in good agreement with the 0.625 value reported previously 80 for the structural transformation in the rutile phase. As W is introduced into the structure, V−V bonding strength is increased as the atoms are driven closer together due to the larger W-ions. This continues until the structure fully transitions into tetragonal with continuous metal-metal bonding along the rutile c-axis or monoclinic a-axis. Average structure data will now be compared to similar local structure information.
"Boxcar" refinements 81 were performed on synchrotron total x-ray scattering PDF where varying r-ranges were fit separately, to create distinct regions of "local" structure, "intermediate" structure, and "long-range" structure ( Fig. S4). The refinements produced goodness-of-fit parameters below 10% supporting the phase purity conclusion drawn from the previous Rietveld refinements. Due to the anisotropic nature of the monoclinic structure ( a = b = c ), three boxcar fits best represented the data if there was any phase fraction of monoclinic structure in either the PDF or the XRD (Tables S11-S18). However, once the structure was fully tetragonal the best fits were achieved with two boxcars corresponding to a decrease in degrees of freedom with a = b � = c . These refinements produced local structure unit cell information that was directly comparable to the average structure unit cell information from XRD.
The data was grouped into two boxcars labeled as the short-range and long-range data corresponding approximately to r-values lying between 0-19.2 Å and 19.2-30 Å, respectively. The r-ranges of each boxcar does depend on the dataset however, and more precise ranges can be found in Fig. S4). Lattice parameters were followed similarly to the Rietveld refinements but for each r-range (Fig. S5). Similarly to the average structure unit cell, there is an increase of the a-axis in both the short-range and long-range, when the structure is mostly monoclinic (Table S19). The hypothesis that W introduction drives V atoms closer together and increases the metal-metal bonding is supported by the almost 5.5 times increase in the a-axis expansion when comparing the local structure to the average structure. This corresponds to a local decrease in c/a of − 0.003(2) at% −1 (Table S20) compared to an increase of 0.00027(1) at% −1 (Table S21) in the long range with the average structure lying between the two at − 0.00071 (6) at%. Most of the recent W-substitution literature 23,61,[82][83][84][85][86][87][88][89][90][91][92][93] uses average structure information to draw conclusions from, but it has been demonstrated that locally the structure changes much more drastically than the average structure indicates.
Approximate strain energies were calculated from the rates of lattice parameter change across the W-substitution series. Energies were derived from typical strain energy 94 , �G S ≈ 4µδ 2 V , where µ is the shear modulus, δ is the unconstrained misfit, and V is the volume of the inclusion atom. If the following assumptions are made: µ correlates with B, δ correlates with V V , and the inclusion atom is W therefore Theoretical strain energies were calculated assuming that the volumetric lattice expansion depends on differences in the atomic radii of W and V, as well as the W-substitution amount such that V = V 0 V W V V x where V 0 is the unsubstituted lattice volume, V W is the effective radii of W in a 6-coordinate environment, V V is the effective radii of V also in a 6-coordinate environment, and x is the fractional amount of W within the lattice. The resulting strain energies with the corresponding 95% confidence interval based on error propagation for the average structure, the local structure (1.50 Å ≤ r ≤ 19.2 Å) from PDF, the intermediate structure (19.2 Å ≤ r ≤ 30.0 Å) from PDF, and the theoretical strain energy (Fig. 2) demonstrate a significantly higher localized strain as well as increased strain overall compared to the theoretical calculations. These localized strains have been previously suggested as nucleation centers in thin-film 96 and single-crystal 36 VO 2 systems.
The PDF data shows an increasing probability of finding V−O as the structure becomes more symmetric with increasing W-substitution demonstrated by the increasing peak intensities of the V−O peaks. The V−V interatomic distance also increases from 2.65 Å to 2.92 Å across the W-substitution series demonstrated by the peak maxima shifting with increasing W concentration (Fig. S6). The V-V distance quickly increases at a rate of 0.07(1) Å/at% when W-substitution is below 3.6 at% but this rate decreases by an order of magnitude at higher substitution amounts (3.6-15 at%) (Table S22). This is reminiscent of the fast increase in the c/a ratio locally followed by a sharp decrease. This increase in the V−V interatomic distance happens similarly in Cr x V 1−x O 2 but conversely Sc-substitution decreases the V−V interatomic distance (Fig. S6).
For the Cr-substitution system, the V−V interatomic distance increases at a rate of 0.010(1) Å/at%. Whereas, for the Sc-substitution system, the V−V interatomic distance decreases at a rate of − 0.0013(9) Å/at%, an order of magnitude lower. The Cr and Sc systems have a significantly slower, almost 1-2 orders of magnitude, respectfully, rate of V−V interatomic distance change (Fig. S6) compared to W-substitution which has a rate of 0.07(1) Å/at% for similar substitution ranges. Phase purity of the Cr-and Sc-substitution systems was verified through Rietveld refinement of synchrotron XRD (Fig. S7). This corresponds well with the MIT differences between the three systems, where Cr increases the MIT by 3 °C/at% 45,47 , Sc decreases the MIT by − 0.7 °C/at% 40 , and W decreases the MIT by − 28 °C/at% 41,61 . The rate of V−V interatomic distance change correlates well ( R 2 = 0.9286 ) with the change in the MIT upon increasing substitution amount of the relevant element. This may indicate that V-V dimerization occurs faster in systems that create a larger MIT change, supporting the Peierls hypothesis of structural distortions driving the MIT 1 . This needs further support through experimentation into oxidation states and electronic correlations, however.
A direct comparison of the residuals from the Rietveld refinement and average real-space Rietveld refinement across the multiple boxcars ( Fig. 3) illustrates that the local structure begins transitioning from monoclinic to tetragonal sooner than the average structure based on tungsten-substitution amount. Even after 0.8 at% inclusion of tungsten, locally the structure fits best to a co-refinement to both the monoclinic, and tetragonal structures.
Local approximate strain is greater than long-range strain from PDF and Rietveld refinements. Both local and long-range strain are greater than theoretical strain calculated from effective ionic radii differences between W and V. www.nature.com/scientificreports/ But this is not seen until 2.5 at% is reached in the average structure. The local structure is then fully tetragonal at 3.6 at% while the average structure needs to be at 6.3 at% to have a complete phase transformation. This is direct tracking of the nucleation and growth event of the phase transformation. This nucleation and growth is also exhibited through in-situ TEM BF imaging during heating and cooling (Fig. 4). Contrast differences distinguish the high-scattering, lower-symmetry monoclinic phase from the lowscattering, higher-symmetry tetragonal phase. Even at temperatures well below the expected T c regions of brighter contrast can be seen amidst regions of darker contrast for both the unsubstituted VO 2 as well as the W 0.008 V 0.992 O 2 (Figs. S8 and S9) illustrating that nuclei of P4 2 /mnm are already within the structure prior to the MIT.
Using the in-situ imaging data, we also extracted approximate rates of transformation by tracking the contrast, or phase, boundary migration while cooling. This was then compared to the rate of V-V interatomic distance Figure 3. Comparison of the goodness-of-fit parameter, R wp , from both XRD (> 30 Å) and PDF data (1.5-30 Å) illustrates how the local structure begins to transition from the monoclinic to the tetragonal structure before the average structure based on tungsten substitution amount. This exemplifies the need for greater local structure studies for predicting and explaining the properties if substituted VO 2 is to be engineered and implemented in a functional capacity.  (21) Å/ • C for TEM and PDF, respectfully, upon cooling. The large errors are due to the large temperature steps taken during TEM cooling (1 °C) and PDF (2 °C). However, upon comparison of the two samples, the W-substituted phase transformation rate is significantly greater than pure VO 2 . This suggests W-substitution promotes the local SPT depression.
In-situ PDF's during heating and cooling were refined similarly to room-temperature PDF. The full r-range was used to construct the phase diagram and the boxcar fits were used to track the local lattice parameters (Fig. S5). These results (Figs. S10-S13) were used to validate the accuracy of the unsupervised machine learning results and compare to the TEM BF rate of transformation data.
Non-negative matrix factorization (NMF) is an unsupervised machine learning technique previously applied to XRD 97 and PDF 98 analyses. It is similar to the principle component analysis (PCA) method 99 for encompassing a whole as a sum of its parts. However, NMF is different from PCA in that the parts of the whole are more intuitive for positive-valued data than PCA 100 . Mathematically, NMF decomposes a compressed data set into non-negative components whereas, PCA utilizes an orthogonality constraint.
NMF approximates an n × m data set matrix, V, by two non-negative matrices W ( n × r ) and H ( r × m ), V ≈ WH . In this in-situ PDF case, NMF was used to cluster the full r-range, 30 Å PDF into pre-transition and post-transition regions. In the case of in-situ PDF, m is the number of NMF components used to form the model, n the number of data sets which in this case corresponds to the number of different temperatures G(r) data was collected at, and r, in the context of NMF, is the number of data points in each G(r). Usually the number of components chosen, m, is less than both m and r, compressing the data set. After compressing the data set, the unsupervised part of NMF quantifies the quality of the approximation by calculating a cost function 101 . This tracks the divergence of the data from a linear mixing of the two end members of the data set 97 . The minimization of this deviation was taken as the SPT transition temperature.
The coefficients of linear mixing of these two components tracked across the temperature gradient, resulted in a sigmoidal curve which was a measure of how closely the two NMF-determined components matched the data against the temperature series. These approximately relate to the amount of the component; therefore, the coefficient of linear mixing for component one was normalized between 0 and 100 for simplicity. The inflection point was tracked through the first derivative, and the peak maximum was taken as the transition temperature (Fig. S14). The NMF used the full r-range of the data (1.5-30 Å) instead of broken down into boxcars as was done when fitting. This did not significantly change the results from the fitting where the local SPT occurred prior to the long-range SPT.
From the in-situ PDF refinements, a phase diagram (Fig. 5) was produced that shows deviations between the local SPT temperature compared to the average SPT temperature supporting the findings from the roomtemperature data that a nucleation and growth phase transformation occurs as observed by others 75 . The rates of P2 1 /c → P2 1 /c + P4 2 /mnm structural phase transformation rate was found to be −23(1) • C/at% given by the NMF analysis. However, a phase transition temperature could not be extracted from in-situ PDF. The rate of www.nature.com/scientificreports/ structural phase transformation between the P2 1 /c + P4 2 /mnm → P4 2 /mnm phases was found to be −21(1) • C/at% and −20.7(2) • C/at% for the PDF fitting and NMF analysis, respectively. Further comparison to the bulk properties ( Fig. 6) confirms the local SPT not only occurs prior to the bulk SPT but also prior to the bulk MIT. The local SPT happening prior to the bulk SPT and the bulk MIT is critical support of the Peierls transformation occurring prior to the Mott-Hubbard mechanism as hypothesized in c-DMFT studies 69 and observed by previous characterization techniques 75 .
Local SPT occurs prior to the bulk MIT. The bulk MIT was characterized using DSC and magnetization experiments (Fig. 6). There is expected hysteresis 37,93,102 in the heating and cooling T c 's. T c was determined from the peak maxima positions of the DSC data, and the inflection point of the magnetization data. To determine the inflection point when the W-substitution was greater than 2.5 at%, fitting the Curie-Weiss law ( χ = C T+θ W + χ 0 ) to the low-temperature ( −200 • C to −263 • C) magnetic susceptibility data was performed ( Fig. S15 and Table S33). This fitting was extrapolated to the full temperature range, and subtracted from the overall magnetic susceptibility to emphasize the inflection point at T c .
Magnetization experiments showed that below the MIT, the system follows Curie-Weiss paramagnetism, and after the MIT it follows Pauli paramagnetism and is temperature independent. A similar study on the V 1−x Mo x O 2 system has been performed 58 . Holman et. al. indicated that the substitution of Mo into the system introduces the Curie-Weiss moment seen at temperatures below −175 • C, which can be extended to the W-substitution system as well. The Weiss temperatures were low, within approximately 10 °C for all samples, demonstrating weak interactions among the magnetic species similar to the Mo-system 58 . The correction factor χ 0 is also very small, on the order of 10 −5 , for all samples. Previous literature found similar results for the Weiss temperature in the W-system 3 .
Using the SPT's gleaned from refinements of the short-range PDF and long-range PDF data as well as the NMF analysis, linear regressions proffered the similarities between the two techniques and the difference when compared against the rate of MIT. All linear regressions fit well producing R 2 values and negative correlation coefficients of ± 0.99 or greater (Table S34).
The NMF asymmetric Gaussian, correlates to the observed MIT phenomenon. For example, in the DSC (Fig. 6a) the full-width half maximum also increases from 7.95(2) to 8.46(3) with increasing substitution amount from 0 at% to 3.6 at%, respectively. The peak intensity decreases as well, also mimicking the NMF-derived SPT behavior. Also, in comparison to the magnetization data (Fig. 6b) the magnetization differences between the pre-MIT and post-MIT decreases as W substitution amount increases, similar to how the changes in heat flow and SPT magnitude decreases upon increasing W-substitution.
Disagreement between the SPT and MIT can been seen in comparison of the transition temperature depression upon cooling in DSC, − 18(2) °C/at%, and from magnetization experiments, − 22.046(3) °C/at%. The magnetization experiment MIT T c agrees with the P2 1 /c → P2 1 /c + P4 2 /mnm SPT, −23(1) • C/at%. However, the DSC MIT T c agrees better with the P2 1 /c + P4 2 /mnm → P4 2 /mnm SPT's, −21(1) • C/at% and −20.7(2) • C/at% for the PDF fitting and NMF analysis, respectively. Since magnetization experiments are more sensitive than DSC measurements, it is not necessarily surprising to observe this. www.nature.com/scientificreports/ Localized vanadium oxidation. Locally, 1.5-3.6 Å as W-substitution increased, three unidentified peaks emerged, while R w p stayed below 10%. By comparing these unidentified local structure peaks between 1.5 and 3.5 Å to simulated PDF patterns, it was deduced that V is undergoing oxidation upon W-substitution. In the V 6 O 13 structure, V is in the 4+ and 5+ oxidation states. The best simulated representation of the data corresponded to increasing amounts of V 6 O 13 in the C2/m phase along with the main VO 2 phase (Fig. 7a) captured by the XRD pattern. As the W-substitution amount increases, the ratio of V 6 O 13 (C2/m) to the main VO 2 phase is best described as a linear relationship from the PDF simulation data (Fig. S16). This was corroborated (Table S35) by EXAFS fitting of thin-film W x V 1−x O 2 (Fig. 7b). According to the Ellingham diagram, V will oxidize W to form V 2 O 5 and WO 2 or WO 3 103 . Due to the lowsubstitution amounts, neither crystalline WO 2 or WO 3 can be seen in the PDF or XRD. But, it seems that there is an intermediate in this redox process is V 6 O 13 . This oxidation process of VO 2 has previously been demonstrated in hydrothermal VO 2 syntheses 104-106 but has not yet been captured in bulk. A potential mechanism for V oxidation is O diffusion through WO x species. The enhanced O diffusivity in WO 3 107 could provide a pathway for V oxidation in this system.
Refined room-temperature EXAFS data of the V K-edge of thin-film W x V 1−x O 2 (x = 0, 0.006, and 0.01) gave similar results to the PDF simulations. The best fits were determined through minimizing the goodness-of-fit parameter, χ 2 ν . All three samples refined best to a combination of paths from VO 2 ( P2 1 /c ), V 6 O 13 (C2/m), V 2 O 5 (Pnma), and in the case of all but the unsubstituted sample, VO 2 ( P4 2 /mnm ) was also needed along with the previous phases to produce the best fit (Table S36). This further supports the oxidation pathway VO 2 → V 6 O 13 → V 2 O 5 . These findings support that the local SPT occurs prior to the average MIT.
It was found previously from magnetic susceptibility experiments that W-substitution forms bonding pairs of V 3+ -W 6+ , V 3+ -V 4+ . This was deduced by comparing the effective magnetic moment to a theoretical effective magnetic moment derived from S 1 = 1 and S 2 = 1 2 3 . Using the Curie constant found from fitting, the effective magnetic moment was calculated from µ eff = ( 3k B C N ) 1/2 µ B where k B is Boltzmann's constant, C is the Curie constant from fitting, and N is the amount of W-ions in moles, and µ B is the Bohr magneton. The theoretical effective magnetic moment was calculated using µ eff = gµ B [S(S + 1)] 1/2 with g as the Landé g-factor taken as two arising from solely spin contributions to the angular orbital momentum 108 . The number of unpaired electrons, n can be approximated by the spin-only magnetic moment, µ SO = µ eff = √ n(n + 2) ( Table 1). Given that the number of unpaired spins is decreasing as the W-substitution amount is increasing, the conclusion is drawn that oxidation is occurring of V and/or W. Adding the results from the PDF and EXAFS fitting, the most likely scenario is the W reducing from 6+ to 4+ while the V is oxidizing from 4+ to 5+ , supported by the Ellingham diagram 103 .

Conclusions
A W-substitution series of VO 2 was analyzed through a slew of diffraction analysis techniques and compared to bulk property measurements from DSC and magnetization. Commonly used average structure identification is not enough to fully capture the complex first-order phase transformation occurring. XRD indeed indicates This local structure is crucial to determine the driving mechanism of this phase transition. More robust analysis of the local structure of phase transitions is needed to comprehensively discuss the origin of these fascinating transformations and ultimately predict them. XRD indicated phase purity, as well as average lattice expansion upon substitution. The local structure transformation occurs more gradually and prior to the average structure phase transformation. This has not been studied before in regards to the W x V 1−x O 2 system. The PDF structural phase transformation was analyzed through conventional real-space Rietveld fitting techniques as well as NMF modeling. The NMF analysis was able to expediently extract the same and more information as the fitting analysis. PDF fitting was unable to uncover the short-range SPT of highly W-substituted ( > 3.6 at%) samples, but NMF was able to extract this information. NMF analyses of other complex phase transitions would be beneficial to demonstrate differences in the local structure with less direct user involvement. Caution should always be taken when using these techniques as to avoid incorrect conclusions from the data but that does not draw away from the power of more machine based data analysis techniques.
It was found that the average structural phase transformation temperature correlates well with the bulk property metal-insulator transition, with all T c depression rates being approximately − 20 °C/at% agreeing with previous literature 1 . The local structural phase transformation temperature however, occurs prior to the MIT and average SPT, as illustrated both through fitting of room-temperature PDF, fitting of in-situ PDF, and NMF analysis of the in-situ PDF. This supports the Peierls-Mott hypothesis of the origin of the MIT as proposed by other in-situ experiments 75 , as well as cluster dynamical mean field theory calculations 69 . Thin-film EXAFS fitting and inspection of highly local (1.5 Å-3.6 Å) PDF data as well as Curie-Weiss fitting of the magnetic susceptibility indirectly uncovered V oxidation due to W-substitution.

Methods
Bulk powder synthesis. Reduction of V 2 O 5 to V 2 O 3 was performed under of flow of 5% H 2 /95% N 2 at 800 °C for 24 h. Molar equivalents of V 2 O 3 , V 2 O 5 , and WO 2 or Cr powder or Sc 2 O 3 were then dry-ground for 15 min in an agate mortar and pestle, and sealed under vacuum using conventional Schlenk line techniques into a 9 mm inner diameter fused silica ampoule (mass loading ≈ 330 mg. Ampoules were then loaded into a Thermolyne FD1540M Box Furnace equipped with a Eurotherm 2614 temperature control unit and annealed at 1050 °C for 216 h (9 days) with a heating rate of 2 °C/min and a cooling rate of 15 °C/min, except for the Sc-substituted system which was annealed at 1050 °C for 264 h (11 days) with a heating rate of 4 °C/min and a cooling rate of 15 °C/min. Thin film synthesis. The procedure for synthesis of VO x used for generating VO 2 thin films was reported by Paik et al. 109 . In brief, oleic acid, 1-octadecene, and dopants were evacuated at 100 °C for 30 min using a Schlenk line technique for purging and refilling with N 2 . After degassing, the solution was exposed to air and the temperature was increased to 200 °C. At this point, VOCl 3 was injected and the solution was aged for an additional 20 min. VO x nanocrystals were collected, washed with ethanol, and centrifuged at 7000 RPM to remove the supernatant. After centrifugation, supernatant was discarded, and nanoparticles were redispersed in hexanes. The dispersion was then dropcast onto 2 cm diameter quartz discs spinning at around 1000 rpm until a homogenous coating was formed. W-substituted VO x thin films were rapidly thermally annealed to crystalline thin films in a MILA-5000 series RTA (Advance Riko Inc.) at 500 °C for 5 min after the chamber was evacuated to < 1 mTorr.

Inductively coupled plasma optical emission spectroscopy. Inductively Coupled Plasma Optical
Emission Spectroscopy (ICP-OES) was executed using a Perkin Elmer Optima 4300DV spectrometer equipped with a Meinhard concentric glass nebulizer. Samples were digested in a 1:1 mixture of 1% HF:HNO 3   was conducted on a TA Instruments Discovery DSC 2500 under a heating and cooling rate of 10 °C/min. Magnetization measurements were performed using a Quantum Design MPMS3 at the Ohio State University Nanosystems Lab with an applied magnetic field of 70 kOe and a cooling rate of 10 °C/min. The sample was packed into a pill capsule, and loaded into a straw to avoid aberrant signal. Centering of the capsule was performed prior to the experiment. Thermal equilibration of the sample at each temperature was ensured prior to the collection of data. Fitting parameters to the Curie-Weiss law included the Curie constant, C, the Weiss temperature, θ W , and a correction factor, χ 0 for the temperature independent component of the magnetic susceptibility. The fitting parameters are presented in the supporting information (Table S33). All fitting was completed in OriginPro.
Synchrotron X-ray diffraction and total scattering. W In-situ PDF experiments were performed at a sample-to-detector distance of 185 mm to accommodate the Oxford Cryosystems Cryostream 700 Plus situated above the sample capillary. The in-situ experimental configuration is shown in our previous literature 35 . A cooling rate of 6 °C/min preceding and proceeding the phase transformation, and 2 °C/min through the phase transformation ( T c ± 10 • C) was performed. T c was determined using DSC and magnetization data prior to the experiments. Rietveld refinement analysis of the synchrotron XRD data was performed using GSAS-II 78 . The following parameters were refined: (i) lattice parameters, (ii) site occupancies for W (Cr, or Sc) and V, (iii) atomic displacement parameters with the cations held equivalent and the anions held equivalent, (iv) fractional atomic coordinates with W, (Cr, or Sc) and V held equivalent, (v) peak shape, (vi) background Chebyshev coefficients of degree eight, and (vii) scale factor. Refinements utilized the following VO 2 CIFs: ICSD-34033 ( P2 1 /c ) and ICSD-1504 ( P4 2 /mnm ). The site occupancy of the V was altered using VESTA to include the ICP-OES substitution amounts.
The V K-edge EXAFS data was normalized in ATHENA 113 . For all samples the ionization energy was set to 5482.03 eV, and R bkg to 1.2 Å. Ṫhe pre-edge range was −150 to −30 eV, the normalization order was three, and the normalization range was, 150-866.446 eV. The spline clamps were strong for both the low-and high-energy data.
The normalized EXAFS data was then refined in ARTEMIS 113 . All refinements occurred over the r-range 1.2-3.6 Å. Ṫhe following parameters were refined: (i) E not , a correctional energy shift, (ii) S 2 0 , the electronic core-hole relaxation was kept equal for all paths given that it depends on the core element which in this case was V, (iii) phase fraction which was represented as a coefficient to S 2 0 and (iv) α × R eff , represented isotropic lattice expansion of the effective interatomic distance from the FEFF calculation Artemis performed. The mean square displacement about the path length, σ 2 , was not refined as it is contingent upon S 2 0 which was already being modified by the phase fraction coefficient.
The phase fraction coefficient was constrained to be between 0.001 and 1.000 for each phase, and the sum of all phase fraction coefficients was restrained to be between 0.999 and 1.000. The paths chosen as the best fit were ones that minimized χ 2 ν while maintaining − 10 eV ≤ E not ≤ 10 eV, and S 2 0 ≈ 0.7 . The Fourier transform range was chosen based off of Ifeffit's suggestion, and k weights two and three were fit for all data sets. The following .cif 's were used in FEFF for path generation: ICSD-34033 (VO 2 , P2 1 /c ), ICSD-1504 (VO 2 , P4 2 /mnm ), ICSD-15028 (V 6 O 13 , C/2m), and ICSD-267175 (V 2 O 5 , Pnma). For structures with more than one V-site, V 6  www.nature.com/scientificreports/ and V 2 O 5 , the structure was aggregated with FEFF prior to path generation. The final fitting parameter results can be found in Supplemental Table S36.
Transmission electron microscopy. The BF TEM imaging was collected using a FEI Tecnai F20 TEM, and the heating experiment was performed using a DENSsolutions MEMS-based Wildfire heating holder with heating/cooling rate of 2 °C/min. Powder samples were prepared by dispersal in ethanol and drop-cast onto the electron-transparent windows silicon nitride, SiN x , in the MEMS-based device. Particles directly attached the edge of the electron transparent windows were chosen to ensure temperature homogeneity across the whole particle.
Non-negative matrix factorization analysis. The dimensionality reduction with non-negative matrix factorization was done using an in-house Python code compiled in Jupyter notebook. The dimensionality of the data is reduced into a 2D space using scikitlearn's NMF module in the decomposition learning class 114 . Dimensionality was reduced to two components, which were compared to the G(r) at each temperature. The first derivative was taken of the sigmoid produced when analyzing the linear coefficient of one of the components as a function of temperature. The resulting Gaussian peaks were fit to an asymmetric BiGaussian using OriginPro. The width on either side of the peak half-maximum was fit and the resulting temperature was used as the local SPT or the SPT-onset temperature and the long-range SPT or the SPT-termination temperature.