Multiphonon diffuse scattering in solids from first principles: Application to layered crystals and two-dimensional materials

Time-resolved diffuse scattering experiments have gained increasing attention due to their potential to reveal nonequilibrium dynamics of crystal lattice vibrations with full momentum resolution. Although progress has been made in interpreting experimental data on the basis of one-phonon scattering, understanding the role of individual phonons can be sometimes hindered by multiphonon excitations. In Ref. [M. Zacharias, H. Seiler, F. Caruso, D. Zahn, F. Giustino, P. C. Kelires, and R. Ernstorfer, Phys. Rev. Lett. 127, 207401 (2021)], we have introduced a rigorous approach for the calculation of the all-phonon inelastic scattering intensity of solids from first-principles. In the present work, we describe our implementation in detail and show that multiphonon interactions are captured efficiently by exploiting translational and time-reversal symmetries of the crystal. We demonstrate its predictive power by calculating the scattering patterns of monolayer molybdenum disulfide $\left({\mathrm{MoS}}_2\right)$, bulk ${\mathrm{MoS}}_2$, and black phosphorus (bP), and we obtain excellent agreement with our measurements of thermal electron diffuse scattering. Remarkably, our results show that multiphonon excitations dominate in bP across multiple Brillouin zones, while in ${\mathrm{MoS}}_2$ they play a less pronounced role. We expand our analysis for each system and examine the effect of individual atomic and interatomic vibrational motion on the diffuse scattering signals. We further demonstrate the high-throughput capability of our approach by reporting all-phonon scattering maps of two-dimensional ${\mathrm{MoSe}}_2$, ${\mathrm{WSe}}_2$, ${\mathrm{WS}}_2$, graphene, and ${\mathrm{CdI}}_2$, rationalizing in each case the effect of multiphonon processes. As a side point, we show that the special displacement method reproduces the thermally distorted configuration that generates precisely the all-phonon diffuse pattern. The present methodology opens the way for systematic calculations of the scattering intensity in crystals and the accurate interpretation of static and time-resolved inelastic scattering experiments.

Figure 1

(a) Schematic illustration of FEDS experiment on bulk ${\mathrm{MoS}}_2$. More details about the setup can be found in Sec. 3a and Ref. [16]. (b) Typical scattering pattern of bulk ${\mathrm{MoS}}_2$. The left subplot shows the raw pattern as collected on the detector. $\mathbf{G}$ indicates the Bragg peak vector (110), $\mathbf{q}$ the reduced phonon wave vector, and $\mathbf{Q}=\mathbf{G}+\mathbf{q}$ the scattering vector. A hexagonal Brillouin zone with the high-symmetry points $\mathrm{\Gamma }$, $K$, and $M$ is also shown. The right subplot shows the difference scattering pattern $\mathrm{\Delta }I(\mathbf{Q},100\mathrm{ps})=I(\mathbf{Q},100\mathrm{ps})-I(\mathbf{Q},t<t_0)$, where $I(\mathbf{Q},t<t_0)$ is the average intensity prior to photoexcitation.

Figure 2

(a) Zero-plus-one-phonon, (b) all-phonon, (c) multiphonon, and (d) Einstein model scattering intensity of monolayer ${\mathrm{MoS}}_2$ calculated for $T=300$ K. In plot (a) we show the fundamental Brillouin zone together with the high-symmetry points $\mathrm{\Gamma }$, $K$, and $M$. We also show the (1 0) and (0 1) Bragg peaks. The blue circle indicates a rapid decrease in the diffuse scattering intensity. The sampling of the Brillouin zone was performed using a $50\times 50$ q-grid. For all plots, the scattering intensity is divided by the maximum Bragg intensity, i.e., with $I_0(\mathbf{Q}=\mathbf{0},T)$.

Figure 3

Percentage contribution of multiphonon interactions to diffuse scattering of (a),(b) monolayer ${\mathrm{MoS}}_2$, (c) bulk ${\mathrm{MoS}}_2$, and (d) bulk black phosphorus (bP) calculated as $\mathcal{P}=I_{\mathrm{multi}}/(I_1+I_{\mathrm{multi}})\times 100$ at $T=300$ K. Parts (a), (c), and (d) represent calculations within the LBJ theory and (b) using ZG displacements. Maps are separated into Brillouin zones to highlight the extent of multiphonon interactions.

Figure 4

Individual and distinct atomic contributions to the all-phonon scattering intensity of monolayer ${\mathrm{MoS}}_2$ calculated for $T=300$ K. Parts (a) and (b) are for the Mo and S individual scattering terms. Parts (c) and (d) are for the MoS and ${\mathrm{S}}_1{\mathrm{S}}_2$ distinct scattering terms. The Brillouin zone sampling was performed using a $50\times 50$ q-grid. Data are divided by the Bragg intensity at the center of the Brillouin zone, i.e., with $I_0(\mathbf{Q}=\mathbf{0},T)$.

Figure 5

(a)–(d) Convergence of the ZG scattering intensity of monolayer ${\mathrm{MoS}}_2$ at $T=300$ K with respect to the Brillouin zone sampling. (e) Exact all-phonon scattering intensity calculated using Eq. (10). All data are divided by the Bragg peak at the center of the Brillouin zone, i.e., with $I_0(\mathbf{Q}=\mathbf{0},T)$.

Figure 6

(a) Zero-plus-one-phonon, (b) all-phonon, (c) ZG (all-phonon), and (d) multiphonon scattering intensity of bulk ${\mathrm{MoS}}_2$ calculated for $T=300$ K. The calculated intensities are divided by the Bragg intensity at the center of the Brillouin zone, i.e., with $I_0(\mathbf{Q}=\mathbf{0},T)$. In plot (a) we show the fundamental Brillouin zone together with the high-symmetry points $\mathrm{\Gamma }$, $K$, and $M$. We also show the (1 0) and (0 1) Bragg peaks. In plot (a) we indicate the (1 0) and (0 1) Bragg peaks, as well as regions of diffuse and Bragg scattering. (e) Zero-plus-one-phonon, (f) all-phonon, and (g) ZG (all-phonon) diffuse maps of bulk ${\mathrm{MoS}}_2$ calculated as $\mathrm{\Delta }I(\mathbf{Q},500\mathrm{K})=I(\mathbf{Q},500\mathrm{K})-I(\mathbf{Q},300\mathrm{K})$, corresponding to the temperature difference estimated from the experiments. (h) Experimental scattering signals of bulk ${\mathrm{MoS}}_2$ measured at 100 ps. Signals are divided by the maximum count due to elastic scattering. For comparison purposes, the simulated data are multiplied by $500000$ to match the experiment. The sampling of the Brillouin zone was performed using a $50\times 50\times 50$ q-grid.

Figure 7

(a) Zero-plus-one-phonon, (b) all-phonon, (c) ZG (all-phonon), and (d) multiphonon scattering intensity of bulk black phosphorus (bP) calculated for $T=300$ K. The calculated intensities are divided by the elastic scattering at the central Bragg peak, $I_0(\mathbf{Q}=\mathbf{0},T)$. In plot (a) we show the fundamental Brillouin zone together with the high-symmetry points $\mathrm{\Gamma }$, $A$, and $X$. We also show the (2 0) and (0 2) Bragg peaks. (e) Zero-plus-one-phonon, (f) all-phonon, and (g) ZG (all-phonon) difference scattering maps of bulk black phosphorus calculated as $\mathrm{\Delta }I(\mathbf{Q},300\mathrm{K})=I(\mathbf{Q},300\mathrm{K})-I(\mathbf{Q},100\mathrm{K})$, compatible with experimental conditions. (h) Experimental difference scattering signals measured at 50 ps. Signals are divided by the maximum count due to elastic scattering. Simulations are multiplied by $400000$ to match the experimental maximum intensity [22]. The sampling of the Brillouin zone was performed using a $50\times 50\times 50$ q-grid.

Figure 8

Individual and distinct atomic contributions to the all-phonon scattering intensity of bulk black phosphorus calculated for $T=300$ K. Plots (a) and (b) are for the individual ${\mathrm{P}}_1$, ${\mathrm{P}}_2$, ${\mathrm{P}}_3$, and ${\mathrm{P}}_4$ contributions. Plots (c), (d), (e), and (f) are for the distinct (and inequivalent) ${\mathrm{P}}_i{\mathrm{P}}_j$ contributions. We also report a ball-and-stick model of bP. The Brillouin zone sampling was performed using a $50\times 50\times 50$ q-grid, and data are divided by the Bragg intensity at the center of the Brillouin zone, i.e., with $I_0(\mathbf{Q}=\mathbf{0},T)$.

Figure 9

Zero-plus-one-phonon, all-phonon, and multiphonon scattering intensities, and percentage contribution of multiphonon processes to diffuse scattering intensity $\mathcal{P}$ calculated for 2D (a) ${\mathrm{MoS}}_2$, (b) ${\mathrm{MoSe}}_2$, (c) ${\mathrm{WSe}}_2$, (d) ${\mathrm{WS}}_2$, (e) graphene, and (f) ${\mathrm{CdI}}_2$ all at $T=300$ K. The energy transfer to the crystal from multiphonon scattering, $\mathrm{\Delta }\mathcal{E}$, the total atomic mass per unit cell in atomic mass units (amu), $M_T$, and the mean phonon frequency, ${\omega }_{\mathrm{E}}$, are indicated on each plot. The sampling of the Brillouin zone was performed using a $50\times 50$ q-grid, and data are divided by the maximum Bragg intensity, i.e., with $I_0(\mathbf{Q}=\mathbf{0},T)$. We also provide the ball and stick model, primitive-cell, and optimized lattice parameter $a$ of each structure.

Figure 10

Individual and distinct atomic contributions to the all-phonon scattering intensity of monolayer ${\mathrm{CdI}}_2$ calculated for $T=300$ K. Parts (a) and (b) are for the Cd and I individual contributions. Parts (c) and (d) are for the CdI and ${\mathrm{I}}_1{\mathrm{I}}_2$ distinct terms. The Brillouin zone was sampled using a $50\times 50$ q-grid. Data are divided by the maximum Bragg intensity, i.e., with $I_0(\mathbf{Q}=\mathbf{0},T)$.

Figure 11

Percentage contribution of multiphonon processes to the all-phonon scattering intensity, $I_{\mathrm{multi}}/I_{\mathrm{all}}$, of monolayer ${\mathrm{MoS}}_2$ as a function of the scattering vector's distance $\left|\mathbf{Q}\right|$ from the zone center. Black circles represent calculations using the exact formula in Eq. (10) and a Brillouin zone q-grid of size $50\times 50$. The green disks represent calculations using ZG displacements and a Brillouin zone q-grid of size $300\times 300$. Red arrows indicate distances for which $\left|\mathbf{Q}\right|=\left|\mathbf{G}\right|$, i.e., when Bragg scattering occurs. The horizontal dashed line intersects the vertical axis at 50%.