Introduction

Thermoelectric materials provide a platform to study phonon–phonon interactions, which is a central topic in condensed matter physics1,2,3,4. The dimensionless thermoelectric figure of merit (\(ZT = \frac{{S^2\sigma }}{\kappa }T\)) implies that good thermoelectrics require low thermal conductivity, κ, and high electrical conductivity, σ, along with a large Seebeck coefficient S5,6. To optimize the thermal transport in materials, it is crucial to understand the correlations among the phonon modes and how they impact the thermodynamics. A good example is thermoelectric tin selenide (SnSe), which has recently attracted tremendous attention due to its ultrahigh ZT of ~2.6 at 923 K7,8. SnSe is a semiconductor with an indirect electron band gap ~0.86 eV9,10,11,12, therefore, phonons are the dominant heat carrier13,14,15. Its high ZT value has been associated with soft, anharmonic distortions off the lattice, which cause strong phonon–phonon scattering16,17,18. The SnSe crystal structure features Sn-Se zigzag-chains along the b axis and armchair-chains along the c axis19. These b–c planes are stacked along the a axis to form an orthorhombic unit cell with strong anisotropy in electronic conductivity, phonon dispersion relation, and optical properties9,16,20,21. Around 800 K, SnSe undergoes a second-order phase transition between two different symmetry orthorhombic phases from Pmna (#62) phase to Cmcm (#63) phase, which is associated with the condensation of a soft optical phonon mode with an Ag symmetry22,23.

Here we report the use of an ultrafast pump-probe approach to directly unveil the soft, anharmonic lattice distortions involved in the strong phonon-phonon scattering in SnSe. These are identified by exploiting the fact that these soft, anharmonic distortions are associated with long recovery timescales after photoexcitation. The ensuing motions are tracked by MeV electron diffraction through its access to large regions of momentum space24,25. We demonstrate a photoinduced domain formation and evolution in the b–c plane in the transient state at 90 K. Combining quantitative ultrafast diffraction analysis and first-principles calculations, we find that a soft phonon mode, transverse optical (TO) Ag phonon mode at the Brillouin zone center in the Pnma phase, is highly excited23. The strong anharmonicity of the phonon mode induces atomic displacements from their equilibrium positions mainly along the c axis. The lattice distortion corresponding to the Ag phonon modes introduces different strain effects to the system. Their competition gives rise to a shear strain between the layers, which significantly reduces the correlation length along the c axis via lattice distortion, yielding in-plane domain nucleation. The domain boundaries likely introduce additional diverse phonon scatterings to reduce thermal conductivity. Furthermore, the contributions from these phonon modes are sensitive to the sample temperature as the strain generation and the domain formation in SnSe are suppressed at 300 K.

Results and discussion

Bragg peak intensity evolution

Figure 1a shows a schematic of the MeV-UED experimental setup, where we record diffraction patterns along [100] zone from a single-crystal SnSe as a function of time delay during the pump-probe process. The normalized diffraction intensities I/I0 (I0 being the peak intensity before time-zero) of (0k0) and (00l) Bragg peaks as a function of time delay are shown in Fig. 1b. Upon photoexcitation, the intensities of all the Bragg peaks were found to decrease. The intensity of {040} reflections, averaged over its equivalent reflections, starts to bounce back around 50 ps, while the intensities of other Bragg peaks, such as, {020}, {002}, {004}, remain low. In Fig. 1c, d we plot the intensities of the {002} and {004} reflections as a function of random thermal vibration (Debye–Waller factors) and atomic displacement, respectively, based on structure factor calculations (for details see Supplementary Note 1). The calculation illustrates that with the increase of the thermal vibration \(\left\langle {u^2} \right\rangle\) the diffraction intensity decreases more rapidly for high-index q {004} than for low-index q {002}. For atomic displacement, however, the intensity decreases with the increase of the displacement Δz of Sn for {002} but increases for {004}. If there only exist the changes of thermal vibrations, as described by thermal mean square displacement \(\left\langle {u^2} \right\rangle\), during the relaxation process after photoexcitation, the normalized intensity change (ΔI/I0) of {040} and {004} reflections would be larger than that of {020} and {002}, which is inconsistent with our experiments, as shown in Fig. 1b. According to these calculations, specific atomic displacements should be triggered by the incident photons. During the first 50 ps after laser pumping, the relative peak intensity can be fitted with a single exponential decay function with a time constant of ~18.3 ± 3.1 ps. The observed timescales support the idea that the atomic displacements correspond to the low energy-cost lattice distortion through electron–phonon interaction and phonon–phonon scattering (See Supplementary Fig. 1). The slow thermalization process is mainly due to the low thermal conductivity in SnSe. Therefore, a combination of thermal vibrations and atomic displacement may be responsible for the observed intensity variations with a rather slow time scale compared with other systems26,27. Additional intensity measurements for the Bragg peaks can be found in Supplementary Fig. 2.

Fig. 1: UED experimental setup and observations.
figure 1

a Schematic of the ultrafast laser-pump electrons-probe setup to study photoinduced lattice dynamics in SnSe. The crystallographic axes are labeled. b Temporal evolution of the normalized diffraction intensity I/I0 of Bragg peaks {020}, {040}, {002} and {004} at 90 K. The solid lines are fits of {020}, {002} and {004} using a single-exponential function; the solid line for the {040} peak intensity is fitted to a double-exponential function. The time constant for each plot is inserted. c, d Normalized intensity calculated based on structure factors as a function of thermal vibration <u2> (c) and atomic displacement Δz for {002} and {004} reflections (d).

Full size image

Asymmetrical intensity distribution of Bragg peaks

The full effect of the phonon contributions to the electron diffraction is also encoded in the diffuse scattering surrounding the Bragg peaks. To monitor the subtle changes in diffraction intensities upon photoexcitation, we constructed experimental intensity difference maps, ΔI (q, t) = I (q, t) – I (q, t ≤ 0), where t ≤ 0 represents the times before the pump arrives and q represents the in-plane momentum transfer. Figure 2a shows a typical difference map of the diffraction patterns for the time-delay of t = 1.7 ps. The figure shows the Brillouin zones of the (020) Bragg peak of SnSe in the Pnma phase (the full diffraction pattern is shown in Supplementary Fig. 3). We observe an angular dependence of the intensity distribution around the Γ points with elongation along the Γ–Z direction. The line profile of the Bragg peaks in the intensity difference map along Γ–Z and Γ–Y directions are shown in Fig. 2b. Away from the Γ point center, the intensity sharply rises along the Γ–Z direction, indicative of a drastic intensity distribution change after photoexcitation. In contrast, along the Γ–Y direction there is little change. The asymmetrical intensity distributions around the Bragg peaks were only observed during the first ~20 ps at 90 K. After that the intensity distribution in the difference maps becomes symmetrical again and decreases equally in all directions. Full details are provided in Supplementary Fig. 4, showing that the asymmetrical intensity distribution of the Bragg peaks changes with the time delay.

Fig. 2: Photoinduced domain evolution observed in UED patterns.
figure 2

a An intensity difference map between the UED patterns taken before time-zero and at 1.7 ps. The Brillouin zones around the (020) Bragg peak (Γ point) with the 1st, 2nd, and 3rd Brillouin zone boundaries marked by white, yellow and black lines, respectively. The (000) beam position is indicated by the circle. b The dots present the line profiles of the diffraction intensity along the Γ–Z (purple dots) and Γ–Y directions (green dots) in a, showing significant asymmetry of the intensity profiles. The solid lines are the intensity profiles from the simulated pattern in c. c Simulated intensity difference map calculated using Monte Carlo based multislice method with large supercell containing domain structure between t ≤ 0 (initial) and t > 0. d Line profiles along Γ-Z direction of the Bragg peak in the UED intensity difference pattern as a function of time. The squares are experimental data, the solid lines are simulation result. e Schematic of the domain formation and evolution after photoexcitation, showing the change of domain size (ξc) along the c direction with time, while the domain size (ξb) along the b direction remains unchanged. The domain size in b and c directions becomes equal in the metastable state. f Magnified drawing from the area marked by the blue rectangle in e, showing lattice distortion after forming domains. g Projected potential calculated based on the distorted lattice in f. The lattice distortion is exaggerated for clarity in e–g. h [100] projected electric potential calculated from the undistorted SnSe structure. The [100] projection of SnSe structure model is embedded with blue and yellow spheres representing Sn and Se, respectively.

Full size image

The origin of the asymmetrical intensity distribution

The asymmetrical intensity distribution in the difference maps of the Bragg peaks in time domain is unexpected and cannot be explained by short-wavelength phonons since there is no observable diffuse scattering in the rest of the Brillouin zone of the UED patterns, further away from the Bragg peaks. Inelastic scattering arising from specific phonon modes as the cause can also be excluded28.

Instead, the variation of reflection peak width can be related to the change of the correlation length (ξ) of a mosaic domain structure in real space. The Equation 2 in Supplementary Note 5 describe the scattering intensity for a finite-sized domain, i.e., correlation length, with dimension of ξa, ξb, and ξc along a, b, and c direction, which may also be applied to the situation where there are a large number of domains scattered incoherently by incident electrons29. Supplementary Fig. 4 shows the shapes and peak width of the reflection for different domain sizes along the b and c directions when viewed along the a direction. For domains with ξb = ξc, the reflection intensity profile is symmetrical. Domains with ξb ≠ ξc lead to an asymmetrical peak broadening, e.g., elongated along the Γ-Z direction if ξb > ξc. The resulting difference map (Supplementary Fig. 5c) captures exactly what we observed in the experiment. We therefore attribute the observed asymmetrical intensity distribution of the Bragg peaks to the temporal formation of domains with specific aspect ratios in response to the photoexcitation.

Our results are consistent with the formation of domains that have various shapes and sizes, arising from the photoinduced lattice distortion in SnSe. Considering the scattering among these photoinduced domains with small lattice distortions, we perform diffraction simulations based on the fact that in our UED measurements the Bragg peak positions remain unchanged before and after photoexcitation and there is no detectable lattice parameters or volume change in the system. We use Monte Carlo and dynamic multislice method by constructing a large supercell consisting of nb × nc × 1 (nb, nc ≥ 64) SnSe unit-cells along b, c, a direction with sampling points of 4096 × 4096 for each slice with thickness of 0.29 nm. The supercell is composed of domains with averaged domain size of m × n × 1 unit-cells, where m and n determine the domain dimension along the b and c directions. We consider the lattices inside each domain are rigid or uniformly distorted based on the corner positions of the domains that deviate from their equilibrium state (Fig. 2f). The small deviation (≤1%) of the domain corner follows a statistically random distribution both in space and in time, accounting for the domains with different distortions in the probing area and the frozen supercells with various distorted domains at a particular time, i.e., at each single-shot acquisition, respectively (see “Methods” section). This approach is consistent with our experimental setup where a UED pattern is recorded at certain time delay through an accumulation of many single-shots. The electrostatic potentials of the supercell were then used to calculate the diffraction patterns based on the multislice method, as shown in Fig. 2g and Supplementary Fig. 6. We refer the details on the Monte Carlo and multislice simulations and domain configurations to Supplementary Note 6.

Figure 2c shows simulated intensity distribution in a difference map, where there is a reduction in domain dimension by a factor of four along the c direction and no change along the b direction after photoexcitation (Fig. 2e). The intensity features in the simulated diffraction pattern agree well with the experimental observations, as shown in the corresponding line profiles along the Γ-Z and Γ-Y directions in Fig. 2b. Figure 2e depicts the domain formation, evolution, and annihilation dynamics based on the Bragg peaks analysis from UED patterns at different time delays (for details, see Supplementary Fig. 4). Figure 2f, g show the lattice distortion (exaggerated for clarity) of the domain and the corresponding projected electric potential. For comparison, the potential for undistorted SnSe crystal is shown in Fig. 2h. Upon photoexcitation, the initial crystal structure of SnSe is fragmented into small mosaic domains by mainly reducing the domain size along the c direction. For instance, the domain size along the c axis is reduced by factor of two at around 1 ps and by a factor of four at around 8 ps, at which the shoulder of the peak intensity reaches the maximum. Subsequently, the small domains annihilate, and eventually the system enters a metastable state with comparable domain sizes along the b and c directions after 20 ps.

We note that the elongation direction of the Bragg peaks does not always point to the Γ-Z direction (Supplementary Fig. 4), but varies in direction at different time delays over a small angular range. This implies that the photoinduced domains may rotate slightly around the [100] direction (Γ-Z direction in reciprocal space). As described in Supplementary Note 7, the effect of domain rotation can be simulated using the same approach by considering both distortion and different rotations of the domains in the supercell (Supplementary Fig. 8d). Similar to the Γ-Z oriented unrotated domains, here the averaged SnSe lattice parameters and orientation remain unchanged, being consistent with the unchanged Bragg peak positions observed in UED experiments. When the domain rotation away from the [100] direction increases, the lattice misfit strain and interfacial energy at the domain boundaries increases. This is likely the reason that the elongation direction of the reflections in UED experiments exhibits a small angular distribution around Γ-Z direction.

Temperature-dependent behavior

The intensity variations of Bragg peaks at 300 K show a similar trend compared with those observed at 90 K (Fig. 3a), implying that atomic displacements also take place at 300 K during the photoexcitation process. However, the Bragg peak intensities drop faster at 300 K (time constant τ300K ~1.7 ps), comparing to that at 90 K (τ90K ~20 ps) after the laser pump (Fig. 3b and Supplementary Fig. 9). The thermal diffuse scattering (TDS) intensity in Supplementary Fig. 1 shows that the saturated TDS intensity at 300 K is smaller than that at 90 K, which demonstrates the photoinduced increase of lattice temperature at 300 K (ΔTlattice_300K) is less than that at 90 K (ΔTlattice_90K). The smaller lattice temperature variation is related to the larger heat capacity in SnSe at 300 K. Additionally, the TDS intensity at 300 K saturates at a relatively shorter time delay compared with that at 90 K, which could be associated with the small temperature jump (ΔTlattice_300K) at 300 K. In terms of the faster response (τ300K ~1.7 ps) at 300 K, apart from the thermalization impact, atomic displacements following the phonon modes also play an important role in the intensity variation. Knoop et al. proposed “anharmonic force” associated with the degree of phonon anharmonicity, which governs the atomic displacements30. In our case, we infer that the anharmonic force of the phonon modes is larger at 300 K in SnSe, which would drive the atoms move faster. Importantly, the asymmetrical intensity distribution, not observed at 300 K (Fig. 3c and Supplementary Fig. 10), only exists for the photoexcitation process at low temperatures and likely not thermally driven. The absence of asymmetrical intensity distribution infers that there is no domain structure with a specific shape and size generated at and above room temperature.

Fig. 3: UED experimental observation at 300 K.
figure 3

a Temporal evolution of the normalized diffraction intensity I/I0 of Bragg peaks {020}, {040}, {002} and {004} at 90 K. b Comparison of the intensity changes of Bragg peaks {002} at 90 and 300 K. The intensity drops faster at 300 K, i.e., time constant τ300K ~1.7 ps < τ90K ~20 ps. c The isotropic distribution around each Bragg peak in the intensity difference map at 300 K at 1.7 ps after laser pump.

Full size image

Possible driving force for the domain formation

To understand the origin of the domain formation, we performed the DFT calculations (Supplementary Note 11). After tens to hundreds of femtoseconds, the “hot” electrons relax to thermal distribution described by electron temperature Te, which is much higher than sample’s base temperature. Due to their larger predicted electron–phonon coupling the higher-energy optical phonons are more likely to be excited than the acoustic modes, yielding a nonthermal distribution of phonons31,32,33,34. Figure 4a shows the calculated partial density of states for Te = 0 K (light color) and Te = 6000 K (dark color), where the 1.55 eV photons excite the electrons just below the Fermi surface to the conduction band. The changes in the electron density of state manifest the simultaneous hole doping into the valence bands and electron doping into the conduction bands.

Fig. 4: DFT calculations.
figure 4

a DFT calculations of partial density of states for the structures relaxed with Te = 0 K (light lines) and Te = 6000 K (thick lines), showing that the 1.55 eV incident photons excite the electrons near the Fermi surface to high energy levels. b Photoinduced atomic displacements and relaxation that can be attributed to the transverse optical Ag phonon modes. The Sn-Se movement increases the intralayer lattice symmetry, while induces interlayer shear strain along c axis, which yields lattice distortions and formation of domains. c Projected weights of the atomic displacements into the four Ag modes of the low-temperature Pnma phase for different electron temperature Te. d Lattice distortion models corresponding to Ag phonon modes with different energy.

Full size image

The DFT calculations predicted the anharmonic lattice distortions from the ground state to a few finite-Te structures, are shown in Fig. 4b. They are characterized by a decrease in the intralayer lattice distortion but an increase in the interlayer lattice distortion, instigating a shear motion along the c axis. Since the structural relaxation with a finite Te does not change the lattice symmetry, the Te induced atomic displacements can be expressed as a superposition of all the four Ag modes referred to Ag-5, Ag-12, Ag-19, Ag-23 based on their energy ranking, as shown in Fig. 4c, d. The calculated results indicate that atomic displacements are mainly governed by phonon modes Ag-5 and Ag-23. While the soft phonon mode Ag-5 represents the thermal phase transition, the shear strain is tied to the highest-energy Ag-23 mode. The projected atomic distortions induced by these two strongly anharmonic phonon modes are in the opposite direction in the b–c plane. Hence, the atomic displacements following these two competitive phonon modes introduce a shear strain along the c direction as shown in Fig. 4b, and the similar photoinduced strain formation has been reported in ref. 35. The extra shear strain can be relieved by inducing lattice deformation in b–c plane and reducing the corresponding correlation length, which results in domain formation in the b–c plane36. Our domain structure model for the diffraction simulation is based on this domain formation mechanism. Due to the introduced shear strain, the atomic position, i.e., corners in each block, can be deviated from the original positions to new positions. And the nearby atoms distort following the displacements of the corner atoms. The atom deviation is random, which reduces the correlation length, i.e., domain formation. Since the atomic displacements are mainly along c direction in SnSe during the photoexcitation, the correlation length in c axis is highly reduced, that is, the domain size along c axis is smaller.

In the domain structure model, the domain boundary is defined by the corners of each block, which are outlined by red and green lines in Supplementary Fig. 6a, d. Additionally, the dense domain boundaries contribute to phonon scatterings, including Rayleigh scattering, Umklapp scattering and displacement scattering, and suppress the lattice thermal conductivity. Upon subsequent release of the shear strain, the opposite pathway of atomic displacements leads to a recovery of the c-axis domain size, or correlation length. The electron diffraction calculations incorporating domain structures and atomic displacements are shown in Supplementary Figs. 6–8. The results illustrate that the domain size dominates the asymmetrical intensity distribution around the Bragg peaks. Based on our experimental observation, the photoinduced domain size starts to increase after t2 ~ 8 ps in Fig. 2d. At longer time delay, t > 20 ps, the asymmetrical intensity distribution disappears, suggesting the correlation length difference in the b and c direction vanishes. Then the atomic displacement following Ag-5 phonon mode and the Debye–Waller factor are dominant in the intensity variation. The calculated diffraction intensities induced by the atomic displacement corresponding to the Ag-5 phonon mode are shown in Supplementary Note 12. Regarding the UED observation at 300 K and the DFT calculation result, the absence of asymmetrical intensity distribution infers that the light-induced atomic displacement is likely dominated by the Ag-5 mode at higher temperatures, which is related to the structure phase transition. The absence of the competition among these phonon modes suppresses the shear strain formation and the domain generation at 300 K.

In SnSe, the Sn-Se bonds along a, b and c axes are anisotropic and anharmonic, which leads to a strong phonon scattering and low thermal conductivity7,8,16. Among these bonds, it has been pointed out that c-axis Sn-Se bond strength is notably weak, which enables the atoms to move preferentially along c axis16,23. Additionally, there is a negative thermal expansion (NTE) in this direction, which could be related to the electronic instability and its coupling to the anharmonic vibration of atoms along the c axis22,37. Additionally, the lattice distortion involved in the NTE is strongly coupled to the anharmonic Ag optical phonon mode17,37. During the photoexcitation process, the atomic displacements corresponding to the selected Ag-5 and Ag-23 optical phonon modes have been realized via electron–phonon coupling at 90 K. Due to the far-from-equilibrium state, the interlayer incompatibility from the lattice displacements introduces a shear strain, reduces correlation length and generates domain structures. The photoinduced strain evolution is incorporated into tens of ps time scale, which is comparable with the timescale observed in35. It is well known that phonon, as the heat carrier, plays a key role in the thermal transport. In laser pumped SnSe, with the induced domain structures discovered here, the enhanced phonon scatterings at the domain boundaries could shorten the mean free path of phonon propagation in the crystal. In this case, the involved phonons range could be from low to medium frequency38, and further reduce thermal conductivity of SnSe, leading to an improved thermoelectric performance.

Bring all pieces together as shown in Fig. 5, we can understand the emergence of domain formation after photoexcitation as following. The electrons are excited on the short (tens of femtoseconds) time scale and the corresponding electron temperature could reach thousands of Kelvin. Then, the energy of the electrons is transferred to the lattice through electron–phonon coupling. Due to the strong anharmonicity of phonon modes in SnSe, atomic displacements coupled to anharmonic phonon modes are realized. Based on our experimental observations and DFT calculations, the atomic displacement along c axis driven by competing Ag phonon modes, dominates in the first 20 ps after photoexcitation at 90 K. The relative Sn-Se atomic displacement introduces a shear strain in the interlayers, which causes more lattice deformations and reduces the correlation length mainly along the c direction. Namely, domains with a specific dimension are generated in the transient state. Additionally, the domain may slightly rotate in the b–c plane to accommodate strain, which impacts on the domain corner positions and the overall lattice distortion in the domain structure. The release of strain enables the atoms to move back to the original positions following other Ag phonon modes, which gradually increases the domain size along the c direction. The absence of domain formation at 300 K infers that the phase-transition related Ag-5 phonon mode dominants after photoexcitation, which significantly reduces the shear strain in the b–c plane. Summarizing our findings, the relatively long-time scale, i.e., tens of ps, observed in our UED experiment is a combination result from the low thermal conductivity in SnSe, lattice displacements induced by the anharmonic phonon modes and the shear strain. Given the close connection to anharmonicity, a non-linear laser-fluence dependence would be expected, which will be tested in future experiments.

Fig. 5: Photoinduced domain formation in SnSe.
figure 5

After photoexcitation, the electrons are ‘heated’ into excited states. During electron excitation and relaxation processes, via energy flowing from electrons to phonons, in the harmonic energy potential, the atom vibration near the equilibrium position is enhanced, which increases the instability of the lattice periodicity. Additionally, the transverse optical (TO) Ag phonon modes are excited through electron–phonon coupling, atoms deviate from the original positions along c axis (Δdc). There are four Ag phonon modes in the b–c plane in SnSe. At 90 K, the direction of atomic displacements dominates at different delay times. The incompatibility reduces the correlation length \(\left( {\xi _{t1}^c} \right)\) generates domains driven by the shear strain between layers. The strain release drives the atoms move to the opposite direction, which gradually increases the correlation length \(\xi _{t2}^c\).

Full size image

In conclusion, using MeV-UED we have observed unexpected anharmonic lattice distortion, which launches a transient domain formation in single-crystal SnSe following photoexcitation. We reveal that this domain formation in the non-equilibrium state is due to an interlayer shear strain, induced by atomic displacements along the c axis. DFT calculation indicates that the observed atomic displacements correspond to TO Ag phonon modes with strong anharmonicity and phonon–phonon scattering in SnSe. Additionally, the transient state of domain structure state only survives for ~20 ps at 90 K. The domains are created by incident photons, and then start to merge and grow during a strain release process. Our observations demonstrate that ultrafast optical pulses are capable of not only manipulating transient and metastable material properties but also creating lattice structures that are not known to exist in the thermal equilibrium state. We expect the domains and their boundaries we observed can be important sources for phonon scattering in reducing thermal conductivity and further improving thermoelectric performance in SnSe.

Methods

MeV UED

The MeV-UED experimental setup at BNL with 3 MeV electron pulses of <180 fs duration transmitted through the sample at the normal incidence, similar to that reported previously39,40. A laser pulse with a duration of 180 fs at 800 nm wavelength and pump fluence of 1.0 mJ cm−1 was used in the experiment which is sufficient to produce a photoinduced lattice dynamic process but prevent photoinduced sample damage. Overall temporal resolution considering RF synchronization jittering (100–150 fs), pump laser duration (180 fs), probe duration (<180 fs), horizontal probe beam size (100–300 µm) and angle between pump and probe (15°), we have a temporal resolution between 250 and 350 fs. The SnSe single crystal with layered orthorhombic crystal structure was exfoliated along the (100) plane and suspended on a 3 mm Mo grid. Measurements were performed at liquid nitrogen temperature (~90 K) and room temperature. Since the beam direction in the UED experiment is parallel to the a axis, our probe is more sensitive to the lattice distortion in the b–c plane. A total of 54 Bragg peaks are captured and analyzed. A UED pattern is recorded by accumulating 72 single-shots. We repeat each measurement ten times and the results were averaged to ensure the data is statistically meaningful. The intensities of equivalent (hkl) and \((\bar h\bar k\bar l)\) reflections are further averaged as I{hkl} to improve the signal noise ratio.

UED simulations

To understand the origin of the intensity elongation in the difference maps, we developed computer codes to calculate UED patterns based on Monte Carlo and dynamical multislice methods. To account for the photoinduced strain and domain formation we built a large supercell which consists of many domains. We consider the lattices inside the domain are rigid and distorted along with the corner positions of the domains when they are out-of-equilibrium (Fig. 2f and Supplementary Figs. 6 and 8). The deviation of the domain corner from their equilibrium positions exhibits statistical randomness following a Gaussian distribution with its mean (or equilibrium position) at the center of the bell-shaped curve. The supercell is used to calculate the UED pattern with the multislice program developed in-house41. More than 20 of such a UED pattern, each corresponding to a different supercell configuration with a random corner position, were then calculated and averaged to account for the multi-shots UED experiment when a diffraction pattern is recorded.

First-principles calculations

First principles simulations were performed by density functional theory in generalized gradient approximation implemented in Vienna Ab-initio Simulation Package. The projected augmented wave method was adopted for the potentials of Sn (5s25p2) and Se (4s24p4) atoms. The finite-displacement method was used to compute the phonon dispersion as implemented in Phonopy code. (See details in Supplementary Note 11).