The nature of photoinduced phase transition and metastable states in vanadium dioxide

Photoinduced threshold switching processes that lead to bistability and the formation of metastable phases in photoinduced phase transition of VO2 are elucidated through ultrafast electron diffraction and diffusive scattering techniques with varying excitation wavelengths. We uncover two distinct regimes of the dynamical phase change: a nearly instantaneous crossover into an intermediate state and its decay led by lattice instabilities over 10 ps timescales. The structure of this intermediate state is identified to be monoclinic, but more akin to M2 rather than M1 based on structure refinements. The extinction of all major monoclinic features within just a few picoseconds at the above-threshold-level (~20%) photoexcitations and the distinct dynamics in diffusive scattering that represents medium-range atomic fluctuations at two photon wavelengths strongly suggest a density-driven and nonthermal pathway for the initial process of the photoinduced phase transition. These results highlight the critical roles of electron correlations and lattice instabilities in driving and controlling phase transformations far from equilibrium.


I. Sample Synthesis and Characterization
: Sample Geometry. a) SEM image of the VO2 thin film. b) Histogram showing the distribution of the nano-grain size in the VO2 thin film determined from SEM images. The mean grain size is ~31 nm, with a 1σ inhomogeneous broadening of ~ 7nm.  2 ), whose intensity decay represents structural transformation of VO2. The red solid line represents the insulator-metal transition, which is measured on the substrate area via the four-probe method. The critical temperatures of structural and electronic transitions are near ~339 K, over a transition temperature window of ± 8 K. (b) Room-temperature ground-state electron diffraction pattern obtained using the ultrafast electron diffraction setup.

II. Experimental Setup
The in situ knife-edge measurements only allow us to measure the pump-laser spot size along x direction, while we determine the size in the other direction (y) by laser-electron cross-correlation approach, as illustrated in Fig. S3(a). In this measurement, we move the pump-laser spot on the sample and monitor the changes of electron diffraction patterns of the sample at a time delay ~+100 ps, which gives us the cross-correlation response between the IR pump laser and the electron beam.
Because the electron beam is focused to ~30 μm on the sample that is much smaller than the pump laser spot size, the cross-correlation measurement can give a reliable measurement of the pumplaser spot size. The intensity of selected Bragg reflection as a function of the displacement of pump laser spot along x direction ( pump x  ) and along y direction ( pump y  ) for both pump wavelengths (800 nm and 2000 nm) is plotted in Fig. S3(b). We note that the displacement values shown in Fig.  S3(b) correspond to the internal reading of the motor which we used to control the moving optics, so they do not directly give the real-space size of the laser pump spot. With proper calibrations, the measurements can provide us a measurement of the aspect ratio of the pump-laser spot for both wavelengths. The estimated error in pump-size is ~17% for 800nm and ~25% for 2000 nm pump, which is determined by both knife-edge measurements and the laser-electron cross-correlation measurements. Figure. S3: Laser-electron cross-correlation measurements. a) Experimental setup of the measurement. The electron beam is focused to ~30 μm on the sample. The laser spot is ~400 μm in size. The laser spot is moved by a motor controlled optics with well-defined steps and the variation of the diffraction pattern is measured after sample. b) The intensity of the selected Bragg reflection is plotted as a function of the displacement of the laser spot on sample along x and y directions. Note that the internal reading of the motor that is linearly related to the real-space displacement is used in the plots.

4
The repetition-rate dependent study has been conducted to verify the full recovery of the samples. Fig. S4 show the results conducted at the base temperature TB=290K. When the repetition rate is set at 1 kHz (1ms between pump-probe cycles), the intensity of (30 2 ) peak before zero-oftime (-100 ps) decreases rapidly when pump laser fluence increases. The transition fluence for the intensity at 150 ps is ~3.5 mJ/cm 2 , much lower than the results obtained at the repetition rates of 500 Hz and 200 Hz. So, the negative-time state of VO2 has an increased temperature than TB in the 1kHz repetition rate experiments due to the thermal build-up between the cycles, resulting in a lower than expected threshold energy. In contrast, at 500 Hz repetition rate, the intensity change at -100 ps is close to zero compared to pump-blocked intensity, even when laser fluence is increased to very high level. By further reducing repetition rate to 200 Hz, we observe the results of fluence scan are consistent with 500 Hz results, which means that at 500 Hz (2ms between cycles) VO2 samples can fully recover to its ground state as without applying the laser pulses. So, the full thermal recover time of our sample under ultrahigh vacuum is between 1ms to 2ms. By operating experiment at ≤500 Hz, we can eliminate the effects of thermal build-up at the negative times. Figure S4: Repetition rate dependent study. The intensity change of the Bragg reflection (30 2 ) at negative (-100 ps, hollow symbols) and positive (+150 ps, solid symbols) times over the fluence ranges covering M1-to-R transition were determined at repetition rate of 200, 500 Hz, and 1kHz. The negative time signals obtained at 1kHz show fluence dependence, indicating that the sample is not fully recovered. In contrast, the experiments conducted at 200 and 500 Hz show no fluence dependence in the negative time frames and have a similar transition curve, indicating that the VO2 samples pumped at these two repetition rate are fully reversible.

III. Calibration of Optical Constants
The optical constants of the employed VO2 samples were determined based on the broadband optical reflectance and transmittance measurements on an as-grown 50 nm VO2 films, deposited on a 45 nm a-Si layer sputtered on a quartz substrate. The optical transfer matrix 60 T 5 was constructed to formulate the correlation between optical transmission and reflection in each layer: (S1) Here, we denote (S2) Here,   For the specific multilayer VO2 film ( Fig. S5(b)) the overall reflectance R and transmittance T measured experimentally can be linked to 2VO T and Si a T  , representing the transfer matrix of the VO2 and a-Si films: where r and t represent reflected and transmitted wave field in the regions I and IV, determining R (r*r) and T (t*t). The measured wavelength-dependent optical transmittance and reflectance are depicted in Fig. S5 We calculate the absorption ratio , defined here as the ratio between the net absorbed fluence (Fabs) by the VO2 film and the applied fluence F, namely =Fabs/F. Assuming the incident light-wave amplitude to be unitary, the absorption ratio of the 50 nm VO2 film is directly given by where Si n is the refractive index of the silicon film, A ) are the amplitude of the reflected (incident) wave field in the a-Si layer and its complex conjugate, respectively. These parameters along with the regions they are defined can be found in Fig. S5(b). The calculation takes into account the difference in the supporting a-Si membranes. In the ultrafast electron diffraction experiment, the a-Si membrane is reduced to 9 nm to ensure transparency for transmission electron diffraction studies. The reported  is determined based on the transfer matrix using the measured optical constants applied to the 9nm a-Si substrate. Nonetheless, it was found that up to 98% of the absorbed optical energy is absorbed by the 50 nm VO2 film in the optical wavelength ranging from 800 nm to 2000 nm in the calibration experiment, hence the absorption by the 9 nm a-Si thin film is in fact negligible.
To compare with other results reported in the literature, it's interesting to see the difference in  due to the optical interference in the thin film geometry. Without considering interference, , and the absorption ratio ' in such a scenario is . Table S1 shows the comparison between  and ', for the 50nm VO2 films.
The main difference occurs for 2000nm, where the net absorption is reduced to half if the interference effect were to be ignored. Including such an effect,  differs by nearly 300% between 800 nm and 2000 nm. However, if the optical interference were to be ignored, the difference would further increase. Such suppression of optical absorption in mid-IR is well linked to the diminishing density of states near the optical gap edge.

IV. Determination of Energy and Optical Dose
With Table S1, we can now calculate the absorbed energy density H : where 2 VO a is the thickness of the VO2 film (50 nm). For consideration of optical doping effect, we calculate the optically generated electron and hole density, which is equal to the absorbed photon density n: To quantitatively compare with the optically induced phase transitions, as shown in Fig.  2(a), we calculate the absorbed thermal energy density in the temperature-induced phase transition at temperature T, based on integrating the heat capacity v C 42 and latent heat L0 of VO2 from the initial crystal temperature TB to T: The latent heat in single crystal VO2, 0 L , has been determined to be 235 6 10  J/m 3 42 , or 1.47 eV/nm 3 .

8
The local bonding changes associated with optically induced M1 to R transition can be inspected based on the time-dependent radial distribution function G(r; t), obtained via a Fourier transform of the normalized structure function sM(s) presented in Fig. 1(b) 62 . Fig. S6(b) depicts the difference radial distribution function G(r;t), obtained by taking the difference between G(r; t) and G(r; -2ps). G(r;t )can be used to identify the formation of new bonds and depletion of old bonds. In transition from the M1 state to the R state, the correlation density at the location of rR (2.85Å) will increase, whereas the correlation density near rS and rL will decrease in the nearest V-V bonding regionsee Fig. S6(a). Similar trend will also occur in the 2 nd nearest neighbor region around 3.52Å, representing the closest V-V bonding between the two sublattice chains, as indicated in Fig. S6(a). These phenomena are evidenced in Fig. S6(b). Whereas the radial distribution function analysis is informative of local bonding changes, for understanding the symmetry-related structural modes directly coupled to electronic excitations, it is more straightforward to inspect the integrated intensity of symmetry-related Bragg reflections.
The atomic movements within the unit cell lead to effects of constructive (intensity increase) or destructive (intensity decrease) interferences of the scattered electron waves, directly contributing the observed diffraction intensities. These changes can be determined by the change of the structure factor for a given direction and zone axis. Specifically, the structural factor F(hkl) is the sum of the scattered electron waves from every atom within the unit cell: .
(S9) Phase transition from the R to the M1 state can be well described by two symmetrybreaking structural modes, namely the pairing of V-V bonding along the c-axis, and the twisting of V-V bonding involving displacements of V atoms towards ab plane (see Fig. S6(b)). These two independent structural modes can manifest in diffraction intensity over two different sets of Bragg reflections. For tracking the pairing status, we can resort to (3 0 2 )M1 and (3 1 3 )M1, which arise only after the lattice dimerizes along the c-axis. In contrast, we can also identify (2 2 0)M1 and (2 3 1)M1, which carry strong contributions from twisting mode. The symmetry recovery from M1to-R transition will diminish the (3 0 2 ) and (3 1 3 ) intensities and increase the (2 2 0) and (2 3 1) intensities. This can occur without changing the octahedral motif, namely the unit cell constant does not instantaneously change over the ultrafast timescale 16 .
The atomic displacement vectors associated with the two structural modes are built on the fractional atomic coordinates based on a fixed unit cell. In the models, the O atoms are frozen during the processes of depairing and detwisting, based on our current understanding of electronically driven structure phase transition between M1 and R 8 . We also keep the volume constant as discussed, namely the unit cell constants are also frozen, only the fractional atomic positions are adjusted. From M1 to R, the unit cell size is reduced by a factor of two exactly. This is in contrast to thermally induced transition, where a lattice strain is transferred from the M1 structure to the R structure 40,63,64 , which leads the following ratios: Here, we set all three ratios to unity throughout the phase transition, justified by the fact that changes in unit cell distances will only affect the position of the Bragg reflections, not the integrated intensity.
In simulating the scattering intensity, we used the Mott-Bethe formula to convert the X-ray form factors to the electron ones 65 . The crystal unit cell is based on x-ray crystallography results of M1 40 . For fitting the dynamics as shown in the diffraction intensity evolution (Fig. S7a), we employed three different dynamical models pathways (see Fig. S7b). Model A considers continuous movements of bond distances and angles from M1 to R as a result of superposing the two structures, which emulates a thermal transition where the intermediate states are mixtures of M1 and R domains prior to establishing a uniform R state. Model B considers a simple two-step transformation from two separate modes with the step (i) predominately conducting depairing along cR, which is followed by the step (ii) changing course to de-twisting. These movements take place simultaneously on A and B chains. Model C considers a metastable M2 separating the two steps. In Model C, sublattices undergo separate de-pairing and detwisting towards untwisted A and unpaired B chains at the M2 stage. The detailed structural displacement steps associated with the three models are depicted in Fig. S7(b) along with the corresponding simulated intensity changes for the 5 Bragg reflections as depicted in Fig.  S7(a). The overall velocity of the movements was obtained by fitting (3 0 2 ) M1 and (3 1 3 ) M1 simulated intensities to the data and the intensity changes of the Bragg reflections (2 2 0)M1, (4 0 2 ) M1 and (2 3 1) M1 were generated accordingly. The results are presented in Fig.  S7(b). From the key features of evolutions, we clearly rule out Model A for lacking any visible step-wise changes in all the intensities. Model B shows steps in (220)M1, ( 1 23 )M1 (as expected), but not in [ 2 40 ]M1, which rises unceasingly with a disproportionally higher intensity than those of (220)M1, ( 1 23 )M1 within the step (i), which is inconsistent with the data. Mode C agrees best with the data, reproducing all the key features.
The Model C agrees well with the experimental results despite of the fact that the movements of individual atoms are highly constrained to only two different structural modes. We note, in reaching this conclusion we have examined the influence of the scale factors of the two structural modes in fitting the experimental results. Using the same scale factor for both structural modes, we found the changes associated with (2 2 0)M1, (4 0 2 )M1, and (2 3 1)M1 reflections, which are directly sensitive to the twisting mode, are greater than the experimental values. Simply reducing the initial detwisting angle by 50%, the agreement is drastically improved. This may suggest that the initial monoclinic state under the surface strain has a reduced lateral distortion, whereas the dimerized V-V bonding distance remains near that of pristine M1 state, putting such a strained initial state to be between the pristine M1 and M2 states. More sophisticated modeling using separate scaling factor to model the amplitude of the two structural modes may capture more details of changes. The model presented here represents the simplest self-consistent model, capturing the essential features of the experiments. We note that in various straincontrolled VO2 nanobeam experiments, the intermediate monoclinic insulating states, namely T or M2, may become metastable 38,59,66 . Such observations could provide the basis for understanding the ground state properties of VO2 nanocrystals. Nonetheless, our results stand on a transient meta-stable state well distinguished in the dynamics on the 1-2 ps timescales, which can only be understood as being driven by electronic excitations, rather than strain or temperature.