Manipulation of Stacking Order in Td-WTe2 by Ultrafast Optical Excitation

Subtle changes in stacking order of layered transition metal dichalcogenides may have profound influence on the electronic and optical properties. The intriguing electronic properties of Td-WTe2 can be traced to the break of inversion symmetry resulting from the ground-state stacking sequence. Strategies for perturbation of the stacking order are actively pursued for intentional tuning of material properties, where optical excitation is of specific interest since it holds the potential for integration of ultrafast switches in future device designs. Here we investigate the structural response in Td-WTe2 following ultrafast photoexcitation by time-resolved electron diffraction. A 0.23 THz shear phonon, involving layer displacement along the b axis, was excited by a 515 nm laser pulse. Pump fluences in excess of a threshold of ∼1 mJ/cm2 result in formation, with an ∼5 ps time constant, of a new stacking order by layer displacement along the b axis in the direction toward the centrosymmetric 1T* phase. The shear displacement of the layers increases with pump fluence until saturation at ∼8 pm. We demonstrate that the excitation of the shear phonon and the stabilization of the metastable phase are decoupled when using an optical pump as evidenced by observation of a transition also in samples with a pinned shear phonon. The results are compared to dynamic first-principles simulations and the transition is interpreted in terms of a mechanism where transient local disorder is prominent before settling at the atomic positions of the metastable phase. This interpretation is corroborated by results from diffuse scattering. The correlation between excitation of intralayer vibrations and interlayer interaction demonstrates the importance of including both short- and long-range interactions in an accurate description of how optical fields can be employed to manipulate the stacking order in 2-dimensional transition metal dichalcogenides.


Phonon excitation at low temperature
The experiment was performed at 110 K using a liquid nitrogen cooling holder. A 41 nm thick sample was pumped with 1.9 mJ/cm 2 laser fluence at 515 nm. The sample was tilted to the [025] zone axis for selected area electron diffraction. The diffraction pattern is shown in Figure S1b with the corresponding simulated diffraction pattern (by SingleCrystal) overlayed. The experimental and simulated patterns show good agreement. Several indexed diffraction spots are shown in Figure S1b. The intensity evolution of two diffraction spots, following photo-excitation, is shown in Figure S1a. The temporal traces show oscillations with opposite phase. The extracted oscillation frequency from the fast Fourier transform shown in Figure S1c is 0.27 THz. A slightly higher frequency than the 0.23 THz obtained in the room temperature experiments shown in Figure 2 of the main text. The phonon is hardened at low temperatures. The initial phase of the oscillation at 110 K is best described by a cosine function, similar to what is observed at room temperature, and is in agreement with the displacive excitation of coherent phonons (DECP) mechanism.

Structure model of shear mode and structure factor calculation
The excited optical phonon which is called 'shear' mode here can be described as an antiphase motion of adjacent layers of WTe 2 along the b axis. The initial shear direction is depicted by the arrows in Figure S2b. For simplicity, we can assume that the middle layer shear towards the negative b direction relative to the top and bottom layer. Since the phonon can be explained by the DECP model, the initial stacking position for the middle layer is at the rightmost position of phonon. The metastable equilibrium position of the shear phonon exhibits a negative b offset for the middle layer of unit cell. Based on this, the shear phonon can be described by the following equation: where A is the amplitude of phonon, T is the period, and d is the displacement.
By assuming T = 4 ps and A = 2.5 pm we can calculate the modulation of the structure factor by the shear phonon according to the following equation: Where f j is the scattering factor and r j is the position of the jth atom.
The intensity of the (hkl) Bragg diffraction spot, I Bragg , is proportional to the square modulus of the structure factor . The atomic position for ground state Td-WTe 2 is obtained from Ref. [1] . When the ℎ shear phonon is excited, the atomic position of middle layer is modulated by d from eq. 1. The relative intensity modulation by the shear phonon were calculated for six diffraction spots and the results are shown in Figure S2a. The simulated results show that the expected oscillations in diffraction intensity from the shear phonon is in agreement with the experimental results shown in Figure 2 and in Figure  S1. Note that the Debye-Waller effect is not considered in the simulation. Figure S2. (a) Normalized modulation of the structure factor for several diffraction spots by a shear phonon along the b axis. The initial shear direction of the antiphase oscillation of the adjacent layers in WTe 2 is illustrated by the arrows in (b). The period of the phonon in the model is assumed to be 4 ps and the amplitude of the shear phonon displacement 2.5 pm.

Lifetime of photo-induced metastable phase
The photo driven metastable phase persists for at least 200 ps as shown by the intensity of the 130 Bragg diffraction spot in Figure S3.

Phonon amplitude as a function of pump fluence
The phonon amplitude is characterized by the amplitude of the intensity oscillation of the 130 Bragg spot in a sample with < 10 nm thickness. To reduce the influence of photo driven phase transition, the amplitude is extracted from the time range 6 -10 ps after time zero. The amplitude of the intensity oscillations as a function of pump fluence is shown in Figure S4. Below 0.5 mJ/cm 2 we observe a near linear increase in amplitude with increasing pump fluence. No significant increase in amplitude is observed after 0.5 mJ/cm 2 , indicating a complex behavior which may be influenced by the formation of a metastable phase. Figure S4. Intensity oscillation amplitude obtained from the 130 Bragg diffraction spot as a function of pump fluence in a thin sample (below 10 nm in thickness). The amplitude is extracted from the time range 6 -10 ps after time zero.

Image of samples with and without strong shear mode excitation
The presence of defects is arguably the reason behind the suppression of the shear phonon presented in Figure 4 of main text. TEM images of a sample with and without strong shear phonon excitation are shown in Figure S5. The flat sample in Figure S5a shows strong photo-excited phonon signal. While the sample in Figure S5b exhibit defect and small ripple domains. Figure S5. TEM images of two samples. In sample (a), the photoexcited shear phonon is strong. While in sample (b), the shear phonon is strongly suppressed. The presence of defects in (b) is arguably the reason behind the suppression of the shear phonon.

Calculation of shear displacement in the metastable phase
The formation of the metastable phase can be modeled by shearing the middle layer of the unit cell towards the negative b direction. The lattice arrives at the metastable equilibrium stacking position through a process with a time constant of ~5 ps. The shear displacement increases linearly with pump fluence until saturation at a pump fluence in excess of ~2.3 mJ/cm 2 (as shown in Figure 5 of the main text). In order to extract the saturation shear displacement, the normalized change in intensity as a function of shear displacement of the middle layer in the Td-WTe 2 unit cell along the negative b direction is calculated for several diffraction spots (solid lines in Figure S6). The intensity of the 130, 150, and 170 spots decrease with the increase of shear displacement while the intensity of the 120 spot increase.
From the intensity gap between 3 ps and 30 ps as shown in Figure 5a of main text the corresponding shear displacement of the metastable phase can be obtained from the solid curves in Figure S6. For conditions with saturated shear displacement, we extract normalized intensities for the 120, 130, 150, and 170 Bragg diffraction spots (indicated by square symbols with error bars in Figure S6). We extract an average saturation shear displacement of ~8 pm by considering the four diffraction spots. Figure S6. Simulated change in normalized intensity for several diffraction spots as a function of shear displacement of the middle layer in the Td-WTe 2 unit cell along the negative b direction. The square symbols with error bars are extracted from the maximized intensity gap observed between the 3 ps and 30 ps traces shown in Figure 5a. Influence of the Debye-Waller effect was considered for the intensity decay at 3 ps.

Diffuse scattering
The intensity of the diffuse scattering between Bragg diffraction spots was traced following photoexcitation. The results are summarized in Figure S7. As can be seen in Figure S7c, the diffuse intensity in Regin I and II increase within 2 ps after photo-excitation which is similar to the instrument response time. The average intensity before time zero is subtracted. Diffuse scattering appears quickly after the photo-excitation. Figure 2a and c.

Fitting method used in
The temporal trace of the diffraction intensity is fitted using a least squares method by a combination of an exponential decay with a sub-ps time constant and the 0.23 THz shear phonon. The thermal effect is considered by an exponential decay with sub-ps time constant as extracted from the change in intensity of the 400 diffraction spot shown in Figure 3b of main text. As presented in the 'Structure model of shear mode and structure factor calculation' section of the text, the intensity oscillation due to the shear phonon can be described by a cosine function. The functions used in the fitting procedure are listed below.
Function for fitting the evolution of the 120, , diffraction intensity: 064 152 (eq. 3) Function for fitting the evolution of the 130, diffraction intensity: 032 252 (eq. 4) Where A is the amplitude of phonon, T is the period of the shear phonon, is the damping constant 1 of the shear phonon, B is the intensity decay due to the Debye-Waller effect, is the decay time 2 constant of Debye-Waller effect, and t is time. was set to be 0.5 ps as obtained from fitting of 400 2 diffraction spot intensity evolution.  Figure 3a The fitting of 130 and 140 diffraction intensity of the 8 nm thickness sample use two exponential functions and a shear phonon, a fast exponential decay with a sub-ps time constant and a slower decay/rise time constant. Since the formation of metastable phase impacts the phonon amplitude and frequency, the first and second period of phonon was assumed to have a different frequency and amplitude. After the initial two periods, the phonon was assumed to have a damped amplitude with a certain frequency. The exact fitting functions are listed below.
Functions for fitting the evolution of the 140 diffraction intensity: , t < T 1 (eq. 5) , T 1 ≤ t < T 1 +T 2 (eq. 6) Functions for fitting the evolution of the 130 diffraction intensity: , t < T 1 (eq. 8) , T 1 ≤ t < T 1 + T 2 (eq. 9) , ] t ≥ T 1 + T 2 (eq. 10) Where A 0 is the amplitude of the shear phonon in the first oscillation, A is the amplitude of the shear phonon after the first period, B is the intensity decay due to the Debye-Waller effect, C is the intensity decay due to the formation of the metastable phase. T 1 is the first period of the shear phonon, T 2 is the second period of the shear phonon, T 3 is the period after the first two periods. is the damping 1 constant of the shear phonon, is the decay time constant due to the Debye-Waller effect. was set 2 2 to be 0.5 ps as obtained from the fitting of the 400 diffraction spots intensity evolution. is the time 3 constant for the formation of the metastable phase. t is time.

Simulated pulse shape vs experimental conditions
The experimental pulse profile is approximately described by a Gaussian function. The onset of the Gaussian pulse is slow, rendering the simulations exceedingly heavy in order to cover the time when the Gaussian profile has insignificant power density. In order to make simulations of the initial pump dynamics feasible, we have employed a rectangular pulse shape starting at ending at , in = 0 = 1 our simulations this is equivalent to continuous light conditions between and . It is clear from the 0 1 intensity profile in figure S8, panel a), that this representation captures approximately the rate of energy increase at the surface of the sample. Using the same peak power density as the Gaussian beam profile at 2.0 mJ/cm 2 , for a rectangular pulse shape, we reach 1.3 mJ/cm 2 after 200 fs, and 1.6 mJ/cm 2 after 250 fs, which is the maximum time reported in the main text. Our simulations use periodic boundary conditions; hence the field is the same in the entire crystal. Considering the slight attenuation of the field in the experimental setup, these conditions are representative of experimental situation. Table S3. Envelopes are described as G for Gaussian, with standard deviation  = 127 fs = FWHM/2.355, for a FHWM of 300 fs. R indicates a rectangular pulse shape, described by a boxcar function, ( ( -, in our simulations equivalent to a CW-illumination with amplitude , between t 0 and 0) -( -1)) t 1 .
Envelope (type, fs) G(=127) G(=127) G(=127) G(=127) R(t=200) R(t=250) R(t=300) Total Intensity (mJ/cm 2 )  Figure S8: Panel a) shows the intensity on the sample for different characteristics of the pump. Panel b) shows the power density profile. To facilitate comparison to the Gaussian profiles, the rectangular pulse is shifted so that intensity equals to 1 mJ at t = 0.