Supplementary Information for Photoinduced anisotropic lattice dynamic response and domain formation in thermoelectric SnSe

1. Condensed Matter Physics and Materials Science Division, Brookhaven National Laboratory, Upton, NY 11973 2. State Key Laboratory of Superhard Materials, College of Physics, Jilin University, Changchun 130012, China 3. Accelerator Science & Technology Initiative, Accelerator Test Facility, Brookhaven National Laboratory, Upton, NY 11973 4. London Centre for Nanotechnology, University College, London WC1E 6BT, UK # These authors contributed equally: Wei Wang, Lijun Wu # Current address: Los Alamos National Laboratory, MS K764, Los Alamos, NM 87545, USA *Corresponding author: zhu@bnl.gov


Supplementary Note 1. Structure factor calculation in SnSe
The diffraction intensity is directly related to the structure factor (SF) where is the scattering factor of atom μ; = 8 2 〈 2 〉 is the Debye-Waller factor with 〈 2 〉 being the mean square displacement of the thermal vibration; q is the scattering vector; and is the atomic position in the unit cell. The structure factor expression indicates that atomic thermal vibration 〈 2 〉 and atomic displacements Δr affect the scattering strength. To understand the role of these two parameters, we calculated the normalized kinematical-scattering-based intensities (structure factor squared), , with the increase of thermal vibration 〈 2 〉 and atomic displacement Δz (displacement along the c direction) of Sn atom, as shown in Figs. 1(c, d). The calculations show that the intensity decreases with the increase of the thermal vibration 〈 2 〉. The higher q, the lower intensity. In contrast to the thermal vibration, as the displacement increases, the diffraction intensity decreases for {002} reflection, but increases for {004} reflection. Comparing experimental UED ( Fig. 1(b)) with the calculations (Figs. 1(c, d)), we conclude that both the thermal vibrations and atomic displacements contribute to the intensity variation in the UED experiment.
Supplementary Figure 1. Top panel: Thermal diffuse scattering signal (TDS) as a function of time measured from 300 K and 90 K. The intensity saturates faster at 300 K than that at 90 K, which are labeled as t300K and t90K in the plot. Bottom panel: Intensity change of {022} reflection taken at 90 K and 300 K for comparison.

Supplementary Note 2. Intensity measurements results at 90 K
Supplementary Figure 1 compares the normalized intensity change as a function of time delay for a series of reflections measured at 90 K. The intensity variations of the Bragg peaks do not simply decrease with the increase of q, indicating that the Debye-Waller factor is not the dominant factor in the dynamic process after photoexcitation. Namely, atomic displacements and lattice distortions are triggered too.
Supplementary Figure 2. Temporal evolution of the normalized diffraction intensity I/I0 for a series of Bragg peaks at 90 K. The solid lines are fitted for eye guide.

Supplementary Note 3. Typical intensity difference map in UED patterns
The intensity difference map is ΔI (q, t) = I (q, t) -I (q, t ≤ 0). The UED pattern taken before time zero is used as the reference pattern, and it is subtracted from the UED patterns after time zero. A typical intensity difference map is shown in Supplementary Fig. 2.
Supplementary Figure 3. Typical intensity difference map at 90 K. (a) Electron diffraction patterns (ED) taken before time zero (t ≤ 0 ps, as a reference pattern), (b) after time zero (t = 1.7 ps). The streaking contrast in the middle is due to the dark-current of the photoelectron-gun in the UED system. (c) The intensity difference map of two ED (t = 1.7 ps) -ED (t ≤ 0 ps). The pattern shown in Fig. 2(a) is cropped from the intensity difference map marked by white frame.

Supplementary Note 4. Intensity difference pattern as a function of time at 90 K
To improve statistics, we repeated the measurement 10 times for each UED dataset at every designated time delay. Supplementary Figs. 3(a-h) show a series of intensity difference patterns as a function of time. These intensity changes can be attributed to the domain formation and evolution (see the main text). The intensity in the shoulder position increases at the beginning of laser pump shown in the top panel, then gradually decreases and disappears around 20 ps in the bottom panel.
The changes of the intensity indicate the domain size reduction and recovery process along c axis during the photoexcitation process. We note that the spots elongation direction sometimes deviates from the -Z direction, e.g., the direction is different between t = 13.3 ps and t = 16.7 ps (Supplementary Figs. 3(e) and (f)), but also from dataset to dataset, e.g., the spots elongation direction at t =1.7 ps in one data set (Fig. 2a) is different from another shown in Supplementary  Fig. 3(c). Statistically, the spots elongation direction shows a distribution around -Z direction with the reduced possibility for large deviation, indicating photoinduced domains may rotate in the b-c plane and their rotation angle follows the same distribution.

Supplementary Note 5. Bragg peak shape analysis
The Γ-Z and Γ-Y line profile measurement from the intensity difference map showing in Fig.  2(b), indicate that peak width along Γ-Z and Γ-Y directions before and after laser pump are asymmetrical. Namely, the correlation lengths in b and c axes in the real space are different upon photoexcitation, i.e., domain structure is formed. To figure out the domain size effect, we first consider the relation between domain size and the shape of diffraction spot g in the kinematic approximation: K = k'k = r* + sg, where K is the reflection vector, r* and sg are the reciprocal lattice vector and excitation error of the corresponding diffraction spot g, respectively. Supposing there are na, nb, and nc unit cells along a, b, and c directions in one domain structure, the amplitude of electron scattering from this domain is: where Fg is structure factor, sg = s1a*+s2b*+ s3c*. (2) The expressions above indicate that the intensity distribution of the Bragg peak depends on the values of na, nb, and nc, that is, domain size. Smaller domain size along one direction, the broader peak width in the direction in the reciprocal space 1 .
For UED pattern viewed along [100] direction, the reflection shape mainly depends on domain sizes along b and c directions. Supplementary Figs. 4(a), (b) show the calculation results based on domain size of 16 u.c. by 16 u.c. and 16 u.c. by 8 u.c., respectively. Supplementary Figure  4(a) presents a symmetrical peak due to the domain with the same correlation length. The peak width along c direction is elongated in Supplementary Fig. 4(b), resulting in a shoulder in the intensity difference map shown in Supplementary Fig. 4(c), which is similar to the observed intensity difference map ( Fig. 2(a) and Supplementary Fig. 3). The corresponding line profiles along Γ-Z is shown in Supplementary Fig. 4(d). We therefore attribute the observed asymmetrical intensity distribution of the Bragg peaks to the formation of domains with aspect ratios in response to the photoexcitation. Here, the domain dimension is equivalent to the correlation length. The above analysis is based on single domain. It may apply to single crystals with many domains if the scattering among domains is incoherently, i.e., the interference between domains is negligibly small. Since the peak width is inversely proportional to the correlation length, the wider peak width in the reciprocal space, the shorter correlation length in the real space. Namely, the correlation length along the c axis is significantly reduced, compared to that along the b axis.

Supplementary Note 6. Monte Carlo multislice electron diffraction simulation for domain structure
To understand the origin of the asymmetric peak shapes and confirm the domain structure in SnSe after photoexcitation, we performed dynamic electron diffraction simulations with various domain structures or configurations. Bloch wave method and multislice method are two main approaches for dynamical electron diffraction and image simulation, and the latter is particularly suitable for defect or disordered structure built onto a large supercell. In the multislice method, the multiple scattering is fully accounted by dividing the specimen into two-dimensional thin slices along the electron beam direction. The electron beam alternately gets transmitted through one slice and propagated to the next slice. Each slice is thin enough to be considered as a simple phase object and the propagation between slices is determined using Fresnel diffraction 2,3 . To include domain structure in the simulation, we use a large supercell consisting of at least 64  64 SnSe unit cells along b and c directions ( Supplementary Fig. 5(a)). We divide the supercell to domains based on their size and shape we are going to generate. Here block is defined as the undistorted and unrotated domain structure. For example, to generate a domain structure with a size of 16 u.c. in b direction and 8 u.c. in c direction, we divide the supercell into 4  8 blocks (the corners of each block are outlined by red circles in Supplementary Fig. 5(a)). The corners of the blocks are then randomly moved to the new positions (black circles) around its original position (red circles). The deviation of the corners from its original position depends on the amount of lattice distortion (≤ 0.1%), which is randomly generated and follows a Gaussian distribution. The lattices inside the block are uniformly distorted based on the corner positions of the block using a bilinear function, as outlined by green lines in Supplementary Fig. 5(d). The electrostatic potential of the supercell is calculated based on the atomic positions in the supercell using the multislice method, as shown in Supplementary Fig. 5(b) for undistorted supercell and Supplementary Fig. 5(e) for supercell with distorted domains. Supplementary Figures 5(c, f) show the Fast Fourier transformation (FFT) of (b, e). The spots of the FFT pattern are single-pixel sharp for undistorted supercell ( Supplementary  Fig. 5(c)). In contrast, the same spots with the distorted domains show spread in intensity along both b and c direction, especially along c direction due to its small domain size in this direction ( Supplementary Fig. 5(f)). In multislice calculations, we divide the supercell to four slices per unit cell along beam direction [100] with thickness of 2.875 Å for each slice. The sampling points of each slice is 4096  4096. The electron beam alternately propagated through these four slices until a desired thickness is reached to form UED pattern on the camera. A supercell generated by one set of corner positions of the blocks, e.g., black circles in Supplementary Fig. 5(d) using a set of random numbers that can be related to the percentage of the lattice distortion and follow a Gaussian distribution. To improve the statistics, we repeat UED simulations with more than 36 configurations. Moreover, account for the sample bending effect due to large beam size (~hundred m) in UED, we calculate the diffraction patterns with the beam along [100] zone as well as slightly deviating from [100] zone axis (deviation angle ≤ 0.5). The final UED pattern at one time-delay is the average of UED patterns calculated based on those domain configurations and beam directions. This is consistent with our experiment setup where UED pattern at each time delay is averaged by many single shots to improve the signal/noise ratio. We found averaging 32 domain configurations is sufficient, beyond that no improvement is observed. The approach is similar to the multislice simulations with frozen phonon approximation for simulating atomic resolution high angle annular dark field (HAADF) images in scanning transmission electron microscopy (STEM). This indicates that our UED simulations are stable and robust for the photoinduced domain structure. The simulated UED pattern is sensitive to domain size and averaged lattice distortion percentage, rather than local lattice distortion of the individual domains. Therefore, the simulation provides useful information about the average domain size and lattice distortion induced by photoexcitation.
The domain structure generated by above method is suitable for multislice-based dynamical electron diffraction calculations because the domain size and distortion can be easily adjusted to compare with UED patterns. As shown in Supplementary Fig. 5(e), there is no disconnect or void, nor atomic overlaps as long as the lattice distortion is small which is the case less than 1% in our simulation. Moreover, the averaged lattice parameters over the supercell remains unchanged, which is consistent with UED observations as we did not observe a noticeable Bragg peak position change in UED experiments. We realize that in the real situation the lattice distortions in the vicinity of the domain boundaries can be more complex due to the lattice relaxation after photoexcitation. For small lattice distortion, the electron scattering from the domain boundaries would remain small, thus will have little effect on the diffraction patterns.
To further demonstrate the effects of domain dimension on the intensity distribution, we did simulation based on the multislice method by changing domain size, the results are shown in Supplementary Fig. 6. In the simulation, we generate initial domains (before laser pump) consisting 32 u.c.× 32 u.c. in the b-c plane (u.c., i.e., unit cell), the sample thickness is 50 nm along beam direction. Such a domain structure is very similar to the misoriented domain in a mosaic crystal. After laser pump, the domain size is reduced to 32 u.c. × 24 u.c., 32 u.c. × 20 u.c., 32 u.c. × 18 u.c., 32 u.c. × 16 u.c., 32 u.c. × 14 u.c., 32 u.c. × 12 u.c., 32 u.c. × 10 u.c., 32 u.c. × 8 u.c., in b and c directions. The partial intensity difference maps between before and after photoexcitation with different domain sizes are shown in Supplementary Figs. 6(a)-(e). The corresponding line profiles along Γ-Z direction around (020) Bragg peak are shown in the plot. The simulation results indicate that the smaller size in the c direction, the stronger intensity distribution shown in the pattern. Based on the simulation results discussed above and the experimental measurement in Fig.2(d) and Supplementary Figs. 3, the simulation results shown in solid lines in Fig. 2(d) present a domain evolution during the photoexcitation process. The domain structures are 32 u.c. × 32 u.c. (t ≤ 0 ps) and 32 u.c. × 24 u.c., 32 u.c. × 18 u.c., 32 u.c. × 14 u.c., 32 u.c. × 12 u.c., and 32 u.c. × 8 u.c. (t > 0 ps) in the simulation, and Se Wykoff position in c axis is changed from 0.481 (t ≤ 0 ps) to 0.477 (t > 0 ps). The sample thickness is 50 nm, Debye Waller factors (t ≤ 0 ps) is 0.006, 0.006, 0.006; Debye Waller factors for Se (t > 0 ps) is 0.007, 0.007, 0.012; Debye Waller factors for Sn (t > 0 ps) is 0.006, 0.006, 0.010; the lattice distortion percentages along b and c direction are 0.1%; 0.1% (t ≤ 0 ps); 0.18%; 1% (t > 0 ps).

Supplementary Note 7. Domain rotation and lattice distortion
We note that the elongation direction of the Bragg peaks deviates from -Z direction, as shown in Supplementary Fig. 3. This implies that the photoinduced domain may rotate in the b-c plane. We can generate the rotated domain by defining the corner position of the original blocks in the undistorted supercell so the shortest side of the original block deviates from the c direction. As shown in Supplementary Fig. 7, we define the original blocks based on the corners marked by red squares, then move these corners to new positions (marked by green circles in Supplementary  Fig. 7(a)) based on a set of random number following Gaussian distribution with its standard deviation equal to the mean lattice distortion. The rotated domains are formed by performing rigid or uniformly distorted lattice inside the blocks. Supplementary Figs. 7(c, d) show electron diffraction pattern calculated based on the rotated domain structure shown in Supplementary Figs.  7(a, b). It is seen that the spots position remains unchanged, but the elongation direction of the spots rotates. Similar to non-rotated domains, the average lattice parameters and orientation over the supercell remains unchanged, and the lattice distortion only depends on the corner positions. However, we notice that the lattice misfit at the domain boundaries increases as with the increase of domain rotation.

Supplementary Note 8. Intensity variations of Bragg peaks at 90 K and 300 K
Intensity variations from the same reflection at 90 K and 300 K are shown below. All the measurement results present that the intensity drops faster after the photoexcitation at 300 K.

Supplementary Note 10. Crystal growth
SnSe single crystals were grown using chemical vapor transport method with iodine agent 4 . Phase purity was checked by X-ray diffraction (XRD) acquired on a Rigaku Miniflex powder diffractometer with Cu Kα (λ = 0.15418 nm) at the room temperature.

Supplementary Note 11. First-principles calculations
First principles simulations of Pnma SnSe were performed by density functional theory (DFT) in generalized gradient approximation (GGA) 5 implemented in Vienna Ab-initio Simulation Package (VASP) 6,7 . The projected augmented wave (PAW) method 6,8 was adopted for the potentials of Sn (5s25p2) and Se (4s24p4) atoms. The plane-wave cutoff energy was 700 eV and the k-point mesh with a grid spacing of 2π × 0.025 Å -1 was used in the calculations. The van der Waals (vdW) interaction with optB86b-vdW correlation functional 9-11 was employed to account for dispersion interactions in this layered material, resulting in the volume smaller by only 0.44% than the experimentally reported structure at 300 K 12 . The finite-displacement method was used to compute the phonon dispersion as implemented in Phonopy code 13 within the framework of harmonic approximation. The calculated phonon density of state (DOS) is shown in Supplementary Fig. 12. The phonon calculations used a 2 × 2 × 2 supercell of the primitive cell, which contained 64 atoms. The DFT relaxed crystal structure and electron band structure were verified by using the Quantum Espresso (QE) Package 14 , for which we used fully relativistic norm conserving pseudopotentials generated using the optimized norm-conserving Vanderbilt pseudopotentials 15 and Grimme's semi-empirical DFT-D3 vdW interaction 16 with the plane-wave cutoff energy of 950 eV and a 6 × 6 × 6 k-point mesh, resulting an increase in the volume by 3.92%.
Supplementary Figure 11. Total and partial phonon densities of state (DOS) for Pnma phase, shown in red (total), green (Sn) and blue (Se), respectively.

Supplementary Note 12. Diffraction calculation
In the first 50 ps at 90 K, the observed reflection intensity is a combination result from Debye-Waller factor, atomic displacements and domain structure evolution. After 50 ps, Since the domain structure has recovered, the atomic displacement following Ag-5 phonon mode and the Debye-Waller factor dominant in the intensity variation. The calculated diffraction intensities induced by the atomic displacement following Ag-5 phonon mode is shown in Supplementary Fig. 12. The calculation results demonstrate that the intensity of {040} reflection increases and {020} decreases with the atomic displacement, which also induces an intensity increase of {004} reflection. In the main text and Supplementary Note 1, we have discussed the intensity variations of {002} and {004} reflection after 50 ps considering atomic displacements and thermal vibration. In terms of intensities from {020} and {040}, apart from the atomic displacement following Ag-5 phonon mode, the thermal conductivity along b axis is larger than that along c axis in SnSe, thermal vibration along b axis decays faster, that is, the atom vibration amplitude along b axis reduces faster, which will increase the intensities of {020} and {040} reflections. The combined result leads to the plot shown in Fig. 1b. Correspondingly, I{033} increases in the diffraction calculation, and I{011} reduces faster than I{022} induced by the atomic displacements, which makes I{033} is higher as shown in the Supplementary Fig. 2. reflections. The intensity is normalized with that calculated based on the original structure, e.g., Δz = 0, for each reflection.