Dynamic disorder phonon scattering mediated by Cu atomic hopping and diffusion in Cu3SbSe3

Cu3SbSe3 that exhibits distinct liquid-like sublattice due to the heterogeneous bonding environment has emerged as a promising low cost superionic semiconductor with intrinsic ultralow thermal conductivity. However, the relationship between atomic dynamics resulting in liquid-like diffusion and anomalous phonon transport properties remains poorly understood. Herein, combing ab initio molecular dynamics with temperature-dependent Raman measurements, we have performed a thorough investigation on the lattice dynamics of Cu3SbSe3. Superionic transition is unveiled for both structurally inequivalent Cu atoms at elevated temperatures, while the Se-formed tetrahedral framework can simultaneously maintain. An intermediate state of Cu3SbSe3 through the mixture of quasi-1D/2D Cu nearest-neighbor vacancy hopping is discovered below the superionic transition temperature. Our results also manifest that phonons predominately involved with Cu contributions along diffusion channels have been strongly scattered during the superionic transition, whereas the liquid-like diffusion of Cu is too slow to completely breakdown the propagation of all transverse phonon modes. The insight provided by this work into the atomic dynamics and phonon scattering relationship may pave the way for further phonon engineering of Cu3SbSe3 and related superionic materials.


INTRODUCTION
Thermal conduction inside condensed system can be considered as internal energy transfer through microscopic collisions of particles and the capability to support propagation of transverse modes often serves as a fundamental criterion to dynamically distinguish solid from liquid 1,2 . Despite the thermal conductivity can be tremendously suppressed by disrupting lattice phonon transport with extrinsic defects or nanostructures 3 , single crystals with intrinsic ultralow thermal conductivity are rare due to their well-defined lattice to conduct both transverse and longitudinal waves. Meanwhile, minimizing the thermal energy flow via earthabundant low thermal conductivity materials is of vital technological importance in many disciplines such as rechargeable batteries 4 , thermal barrier coating 5 , and high-performance thermoelectrics [6][7][8] . As the thermoelectric conversion efficiency is inversely proportional to thermal conductivity, a desired strategy to achieve advancements in thermoelectrics follows the recent paradigm of phonon-liquid electron-crystal (PLEC) 9 , which evolves from phonon-glass electron-crystal 10 , and formulates that the decoupled lattice thermal conductivity can be decreased even below glasses because the propagation of transverse phonons are suppressed by the liquid-like fluctuating sublattice. Superionic conductors that exhibit long-range liquid-like ionic behavior while simultaneously maintaining solid crystalline sublattice for electronic conduction, therefore, have emerged as a class of promising candidate to realize intrinsic ultralow thermal conductivity and represent potential PLEC materials [11][12][13] .
Although structural analysis has demonstrated that the liquidlike thermal transport property could be closely related to the presence of energetically accessible interstitial voids which enable intensive vibration of loosely bonded rattlers 9,14-16 , the thermal structural evolution is yet to be fully understood in respect of lattice dynamics in non-caged superionic systems 12,17 . More fundamentally, the relationship between phonon scattering and liquid-like stochastic diffusivity also remains elusive due to the lack of nonperturbative treatment of anharmonic phonon-phonon interactions 18 . It has been argued that whether the transverse phonons are completely suppressed by the liquid-like diffusion depends critically on the jumping time scales involved 19 . Therefore, a comprehensive study of the relation between the atomic dynamics and the corresponding anharmonic phonon properties at elevated temperatures is fundamental to understand the phonon scattering evolution across the superionic transition.
Despite possessing identical elements and similar stoichiometries, the ternary Cu 3 SbSe 3 and Cu 3 SbSe 4 compounds, manifest distinctive thermal transport properties. Cu 3 SbSe 4 demonstrates a classical thermal conductivity 20 , whereas Cu 3 SbSe 3 exhibits a liquid-like thermal conductivity of~0.7-1.0 Wm −1 K −1 at 300 K and a strong temperature dependence deviation from characteristic T −115,21 . Thus, these compounds provide us an excellent opportunity to dynamically probe the thermal structure and phonon scattering evolution from well-defined solid vibration to continuous superionic diffusion. Although diverse frameworks have been proposed to theoretically account for the anomalous phonon transport in Cu 3 SbSe 3 , including strong lattice anharmonicity 20 , lone s 2 pair electrons 22 and chemical bond hierarchy 15 , a detailed atomistic understanding has remained elusive due to the lack of sufficiently long timescale simulations that can fully capture the dynamics of the collective liquid-like ions and the corresponding anharmonic phonon calculations.
Herein, an effective hybrid scheme 23,24 based on ab initio molecular dynamics (AIMD) simulations, which is capable of treating the anharmonicity of the atomic interactions at any order, has been applied to depict the phonon quasiparticles. Combined with theoretical approaches, temperature-dependent Raman spectra have also been experimentally measured on Cu 3  and Cu 3 SbSe 4 pellets to further elucidate their distinct phonon behaviors. Our results reveal the activation of superionic transition for both structurally inequivalent Cu atoms in Cu 3 SbSe 3 , resulting in a collective broadband phonon scattering. A distinct temperature-dependent diffusion mechanism has been found for asymmetrically bonded Cu atoms across the superionic transition, which is initiated by normal solid-like vibration, followed by the thermally activated Cu hopping and then continuous liquid-like Brownian diffusion. In addition, our results indicate that phonon quasiparticles dominated by Cu atoms are strongly scattered across the superionic transition along diffusion channels, while the liquid-like diffusion of Cu is too slow to completely collapse the propagation of all transverse phonon modes of Cu 3 SbSe 3 .

RESULTS
Thermal structure evolution of the inhomogeneous Cu sublattice The crystal structure of Cu 3 SbSe 4 is close-packed diamond-like with Cu atoms tetrahedrally coordinated by four equivalent Se atoms, as shown in Fig. 1a. On the other hand, as illustrated in Fig. 1b, Cu 3 SbSe 3 crystallizes in a more complicated orthorhombic structure with inequivalent Cu atoms located inside the Sb distorted tetrahedrons, which provide an inhomogeneous sublattice and enable Cu1 and Cu2 to preferentially oscillate along z direction (1D) and x-z plane (2D), respectively. Cu 3 SbSe 3 manifests liquid-like dynamic characteristics above the critical transition temperature of around 380 K 25 . The thermal structure evolution for mobile sublattice has been first tracked by evaluating the root mean square displacement (RMSD) Δ(r) during AIMD simulations, which often serves as an indicator to quantitatively distinguish solid from liquid 26 .
As shown in Fig. 1c, d, Δ(r) Cu1 and Δ(r) Cu2 plateau at T = 100 K, corresponding to a solid vibrational feature originated from the crystalline bonding. Upon heating to 200 K, however, a sudden increase of Δ(r) Cu1 is observed at around 95 ps, which indicates a thermally activated vacancy hopping of Cu1 atom from a shallow potential well to the next equilibrium position; whereas, only slight deviation from the classic solid characteristic is found for Δ(r) Cu2 . At 300 K, a similar hopping pattern with shorter hopping residence time has been observed for both Cu1 and Cu2 atoms, which can be attributed to the low bonding strength to the inhomogeneous sublattice. After the thermally induced occupational transition, Cu atoms are trapped locally in the neighboring equilibrium positions, thus leaving a plateau in RMSD curve as illustrated in Fig. 1c. Furthermore, warming across 400 K, Δ(r) Cu1 and Δ(r) Cu2 increase approximately linearly with ffiffi t p and exhibit distinct disorder effect and continuous liquid-like diffusion, which is consistent with previous temperature-dependent atomic dynamics studies 15,21 and high temperature X-ray diffraction experiment 27,28 for Cu 3 SbSe 3 .
To further examine the thermal structure evolution, anisotropic AIMD trajectories of Cu1 atom have been shown in Fig. 1e, which clearly demonstrates the aforementioned temperature-dependent diffusion process and is compatible with the pair distribution function, as shown in Supplementary Fig. 2. A significant broadening of the pair distribution function is observed for Curelated peaks after the activation of the thermal hopping, indicating Cu atoms lose their initial coordinations and gradually become more randomly distributed at elevated temperatures. In addition, projection plots in Fig. 1e and Supplementary Fig. 1 indicate the liquid-like mobility of Cu1 and Cu2 atoms arises overwhelmingly from the motion along the interstitial void in z and x-z directions, respectively, consistent with previous structure analysis based on atomic displacement parameters 21 .
Atom-projected phonon dynamics across the superionic transition To illustrate correlated phonon evolution from normal solid-like vibration to continuous liquid-like diffusion in the superionic regime, atom-projected power spectra of Cu 3 SbSe 4 and Cu 3 SbSe 3 are calculated through the Fourier transform of the projected velocity autocorrelation function, as shown in Fig. 2a, b. The slightly progressive broadening and softening of atomic resolved power spectra for Cu 3 SbSe 4 at elevated temperatures manifests relatively weak lattice anharmonicity and strong bonding within the symmetric Se-formed tetrahedrons, which is consistent with previous theoretical calculation based on quasi-harmonic approximation and three phonon scattering process 20,21 . The atomprojected power spectra of Cu 3 SbSe 3 , however, reveal a prominent broadband phonon scattering with a dramatic loss of Cu spectral weight accompanying the activation of atomic vacancy hopping. As one can see, Cu dominated power spectra show pronounced peaks at low T, which then undergo considerably broadening upon heating to 200 K, indicating the dominance of Cu hopping phonon scattering. In addition, the phonon peaks (4.0-6.5 THz) mainly involved with Sb and Se contributions also undergo sudden changes from 100 to 200 K, associated with the distortion of the tetrahedral framework resulted from the hopping of Cu atoms; sudden increases of RMSD of Se and Sb atoms are also observed as shown in Supplementary Fig. 3. After the superionic transition, the Cu power spectra at low frequency range become extremely broad and almost disappear with no sharp peaks remaining, in accordance with alike superionic materials involving inhomogeneous bond strength 9,16 . In contrast, phonon peaks mainly involve Sb and Se atomic contributions are relatively less affected across the superionic transition, signifying the coexistence of the solid supporting tetrahedral framework. The collective large amplitude anharmonic vibration of Cu atoms, therefore, plays a similar role as the rattlers within cage compounds or the diffusive atoms in intercalated layered compounds 16,29,30 , but involving a much broader characteristic size for heat carrying lattice phonon scattering due to the widely distributed Cu dominated phonon modes. To further clarify the physical interpretation of the phonon peaks and the structure evolution, the corresponding AIMD trajectories of Cu 3 SbSe 3 supercell at different temperatures in the y-z plane have been depicted in the right panel of Fig. 2b. Upon heating from 100 to 300 K, Cu1 atoms hop predominantly along the path from the tetrahedral center toward the tetrahedral face and finally trap locally after passing through the boundary face, resulting in an intermediate state before the superionic transition of Cu 3 SbSe 3 .

Raman spectra of Cu 3 SbSe 3 and Cu 3 SbSe 4
The temperature-dependent Raman spectra were recorded for pellets of Cu 3 SbSe 4 in a temperature range from 198 to 398 K to investigate the phonon behavior, as shown in Fig. 3c. A factor group analysis of Cu 3 SbSe 4 structure (I42m) indicates the main peak located at 5.5 THz can be assigned to the stretching vibrations of the symmetric CuSe 4 /SbSe 4 tetrahedron; this peak demonstrates moderate broadening and softening as temperature increases, in agreement with the typical homogeneous bonded crystalline characteristic and previous Raman spectrum measurements based on thin film samples 31,32 . The gradually disappeared Raman peak that locates near 6.9 THz as temperature increases appears to be related to the intermediate Sb 2 Se 3 deposition and the peak presents near 6.6 THz may originate from Sb precipitation 32 , which is consistent with the XRD and scanning electronic microscopy (SEM) images of Cu 3 SbSe 4 samples as illustrated in Fig. 3a and Supplementary Fig. 4. Shown in Fig. 3b is the XRD pattern recorded for the Cu 3 SbSe 3 pellet. In contrast, the Raman spectra for Cu 3 SbSe 3 broaden noticeably even below the superionic transition temperature, which could be attributed to the strong anharmonic vibration of the weakly bonded Cu atoms. The most significant phonon peak locates at 6.0 THz, which is mainly contributed by the collective motions of Cu and Se atoms, gradually disappears at elevated temperatures; the Raman spectrum also exhibits discrete changes in frequency and linewidth across the superionic transition, consistent with previous experimental findings on signature order-disorder transition of Cu 3 SbSe 3 25,33 .
Differences from previous Raman spectra on peak quantity 25 , however, are likely related to the different sample synthesis methods 34 . The distinct phonon behaviors between Cu 3 SbSe 4 and Cu 3 SbSe 3 , therefore, indicate the pivotal role of inhomogeneous Cu sublattice in dynamic disorder phonon scattering.
Mode-projected phonon dynamics across the superionic transition To further analyze the response of different phonons across the superionic transition, the intrinsic mode-projected phonon frequencies and linewidths at different q-vectors have been calculated as shown in Fig. 4a, b. For both q = (0.5, 0.0, 0.0) and q = (0.0, 0.5, 0.0), a dramatic broadening of phonon linewidth, which is inversely proportional to the phonon lifetime, is observed upon heating from 100 to 400 K, indicating the giant dynamic phonon scattering induced by the liquid-like Brownian diffusion. After transforming into the superionic state, the diffusive dynamic scattering has demonstrated mode selectivity; although the phonons predominately involved with Cu contributions in the middle frequency range collapse, phonons arise from the collective vibrations of all atoms at lower end (<1 THz) are relatively less affected. Furthermore, by analyzing the polarization vectors of the three lowest energy phonon modes, as illustrated in Fig. 4c, d, we find that the most strongly affected transverse phonons are predominated by Cu1 vibrations along the z diffusion channel, while the longitudinal phonons that involve Cu1 vibration in planes perpendicular to the flow direction are less affected after the superionic transition. The participation of Cu1 vibration to those less influenced optical phonons is also found to be insignificant, which is consistent with the persistence of non-Ag dominated phonons in the superionic phase of AgCrSe 2 35 . The survived phonon modes at low frequency range indicate the liquid-like Cu diffusion may not be fast enough to completely breakdown the propagation of all transverse phonon modes in Cu 3 SbSe 3 , which is different from typical phonon-liquid superionic materials, for example, Cu 2 Se 9 . The Cu1 self diffusion constant estimated from Einstein's Brownian motion equation is D = 4.05 × 10 −7 cm 2 s −1 at 400 K and increases to D = 2.87 × 10 −6 cm 2 s −1 at 500 K, which is relatively small compared with that of multi-layered superionic CuCrSe 2 12 and AgCrSe 2 36 . The strong dynamic phonon scattering with broad characteristic frequencies due to the weakly bonded Cu atoms rather than the collapse of all transverse phonons, therefore, is more likely to be responsible for the suppression of phonon transport in heterogeneous bonded Cu 3 SbSe 3 37 .

DISCUSSION
In summary, we have explored the structure evolution and the corresponding phonon scattering mechanisms for Cu 3 SbSe 3 and Cu 3 SbSe 4 through AIMD simulations and temperature-dependent Raman spectrum measurements. Our results reveal the existence of atomic level structural heterogeneity facilitated by the activation of superionic transition of Cu 3 SbSe 3 . An intermediate state of Cu 3 SbSe 3 has been discovered through the Cu nearestneighbor vacancy hopping along multiple diffusion channels below the superionic transition temperature. Phonons predominately involved with Cu contributions along the diffusion channels are strongly scattered, whereas the liquid-like diffusion of Cu may not be fast enough to completely breakdown the propagation of all transverse phonon modes. The insight provided by this work into the atomic dynamics and phonon scattering relationship may pave the way for further phonon engineering of Cu 3 SbSe 3 and related superionic materials.

Computational details
Phonon calculations for both Cu 3 SbSe 3 and Cu 3 SbSe 4 were performed using the plane-wave basis Vienna ab initio simulation package 38 , implementing the collinear spin-polarized density functional theory and the Perdew-Burke-Ernzerhof generalized gradient approximation functional 39 . We adopted the projector-augmented wave 40   lattice parameters and atomic positions were fully optimized until atomic forces were smaller than 1 meV Å −1 .
Harmonic force constants were obtained using the small displacement method as implemented in Phonopy 42 , with an atomic displacement amplitude of 0.01 Å. We used supercells of 2 × 2 × 1 (112 atoms) and 2 × 2 × 2 (128 atoms) for Cu 3 SbSe 3 and Cu 3 SbSe 4 , respectively. The duration of MD simulations was 240 ps with a 4 fs time step. A Γ point mesh combined with an energy convergence of 10 −5 eV was applied to all AIMD simulations. Before data collection, all systems were initially relaxed by minimizing the potential energy and then equilibrated for at least 20 ps using NVT ensemble with Nose-Hoover thermostat. The first half of the AIMD trajectory of Cu 3 SbSe 3 had been discarded during the power spectra calculation due to the extremely slow relaxation process mediated by Cu hopping. Unless otherwise noted, remaining AIMD parameters were kept identical to those of the harmonic phonon calculations. The modeprojected power spectra were calculated by projecting the atomic velocities onto the phonon eigenvectors using the normal-modedecomposition technique 24 , as implemented in DynaPhoPy 43 . The power spectra were then convolved with Lorentzian function to obtain phonon peak positions and linewidths.

Experimental details
Polycrystalline Cu 3 SbSe 3 and Cu 3 SbSe 4 samples were synthesized by melting the stoichiometric amount of high purity elements (>99.99%) at 923 K for 8 h, quenching in cold water followed by a further annealing at 723 K for 3 days. The annealed ingots were hand ground into powders for identifying the phase compositions. The powders of Cu 3 SbSe 3 and Cu 3 SbSe 4 were densified for 30 min by hot pressing under a uniaxial pressure of 50 MPa at 673 K and 80 MPa at 643 K, respectively. The obtained pellets (>98% of the theoretical densities) were used for Raman measurements. The phase compositions of the samples were identified by X-ray diffraction (DX2700) and SEM equipped with energy dispersive spectrometer. Raman spectra were excited by the 532 nm lines of an argon laser in the back-scattering geometry, using a Jobin Yvon model U-1000 monochromator equipped with a conventional photo-counting system. The sintered samples (pellet in geometry) were held onto a glass slide for Raman measurements. The sample temperature was varied from 173 to 398 K by using a continuous flow LN 2 optical cryostat. To ensure a better signal-to-noise ratio, four sets of measurements were counted for averaging.

DATA AVALIABILITY
All data generated or analyzed during this study are included in this published article (and its supplementary information files).