Atomic‐Scale Visualization and Quantification of Configurational Entropy in Relation to Thermal Conductivity: A Proof‐of‐Principle Study in t‐GeSb2Te4

Abstract It remains a daunting task to quantify the configurational entropy of a material from atom‐revolved electron microscopy images and correlate the results with the material's lattice thermal conductivity, which strides across statics, dynamics, and thermal transport of crystal lattice over orders of magnitudes in length and time. Here, a proof‐of‐principle study of atomic‐scale visualization and quantification of configurational entropy in relation to thermal conductivity in single crystalline trigonal GeSb2Te4 (aka t‐GeSb2Te4) with native atomic site disorder is reported. A concerted effort of large t‐GeSb2Te4 single crystal growth, in‐lab developed analysis procedure of atomic column intensity, the visualization and quantification of configurational entropy including corresponding modulation, and thermal transport measurements enable an entropic “bottom‐up” perspective to the lattice thermal conductivity of t‐GeSb2Te4. It is uncovered that the configurational entropy increases phonon scattering and reduces phonon mean free path as well as promotes anharmonicity, thereby giving rise to low lattice thermal conductivity and promising thermoelectric performance. The current study sheds lights on an atomic scale bottom‐up configurational entropy design in diverse regimes of structural and functional materials research and applications.

Where I 0 is the background intensity, A i is the amplitude of each peak. All of the Gaussian peaks shared the same width . The experimental images are pixels × pixels, corresponding to nm×nm.
Theoretically, the incident probe intensity and the HAADF detector should be taken into account in the quantitative HAADF simulation process. [1][2][3] We measured the intensity of the incoming beam and obtained the electron probe on an absolute scale by applying the method of LeBeau et al. [1][2][3] , but without external amplifier. Figure S1 shows the color map of the intensity of the incident electron probe while scanning over the detector and the radially averaged sensitivity of the detector as a function of the scattering angle. And then the GST image intensities relative to the intensity of the incoming electron beam allowed us to direct compare with simulated image intensities. The absolute normalized intensity can be expressed as , where is the raw intensity, is the intensity in the vacuum region, is the averaged intensity on the detector. Figure S1 (a) Color map of the intensity of the incident electron probe while scanning over the detector, (b) the radially averaged sensitivity of the detector as a function of the scattering angle.
In general, the crystal structure of Ge 1 Sb 2 Te 4 (lattice system: rhombohedral, crystal system: trigonal, space group 3 ̅ ) is derived from the cubic ABC stacking sequence along c axis consisting of three septuple layers interpenetrated by van der Waals gaps, which can be represented in terms of the stacking of the close-packed planes described by AβCαBγA BγAβCαB CαBγAβC. The capital Roman letters represent Te layers and the Greek letters represent the Ge/Sb layers. Thus, the stacking sequence of Ge 1 Sb 2 Te 4 is -Te1-GS2-Te2-GS1-Te2-GS2-Te1-, where Te atoms occupy two non-equivalent sites (denoted as Te1 and Te2), and similarly, GS1 and GS2 are non-equivalent octahedral cation sites for Ge and Sb atoms, thereby forming Te sublattice and Ge/Sb sublattice. The atomic number difference between Te and Ge/Sb naturally bring about the contrast in HAADF-STEM image. Figure S2 shows the HAADF-STEM image of GST along the [ ̅ ] direction, and the inset is the corresponding fast Fourier transform (FFT) image. The intensity fluctuations of Ge/Sb sublattice (the weak peak between two adjacent strong peaks from Te sublattices) in the intensity line-profile that correspond to the septuple layers marked by the red box indicate cation site occupancy disorder.

Figure S2
The filtered HAADF-STEM image of GST along the [ ̅ ] direction, and the inset is the corresponding FFT pattern. The intensity line-profile (right panel) of septuple layers marked by the red boxes.
First, we obtain the absolute intensity and the position of each atomic column in a HAADF image based on the Gaussian distribution as shown in Figure S3. Then, we apply the IDA method to better visualize the change of spot brightness in the HAADF image. Specifically, the absolute intensity of each Te column in Figure S3(b) (cf. the solid circles in the middle in Figure S4(b)) is first divided by the average intensity of its "nearest" Te neighbors (cf the dotted circles in the same figure) as the normalized intensity of Te sublattice, the result is shown in Figure S4(c). The intensity of a Ge/Sb column is further divided by the average intensity of adjacent Te columns (taken as the reference, cf. Figure  S4(b)) to obtain the normalized intensity of Ge/Sb column, which indirectly reflects the cation atoms distribution in the columns as shown in Figure S4(d). The IDA process can be described using the following two equations.
where the denotes the absolute intensity from the HAADF image of the neighboring Te columns (cf. the red circles in Figure S4(b)), and are the intensities of center anion Te and those of center cation Ge/Sb (cf. the yellow circles in Figure S4(b)), respectively. The and are the normalized intensity of anion column and cation column, respectively; , are the numbers of the neighboring Te columns for the center anion Te and cation Ge/Sb columns, respectively.
Specifically, the number of the neighboring Te columns for the center anion Te is six in GST, and the number of the neighboring Te columns of the center cation Ge/Sb columns is four. (n=6, m=4). So the normalized intensity equations of GST show as following: The normalized intensities are represented by different colors (see the colorbar) in the normalized maps (Figure S4(c) and (d)). In ⅣⅤ-V 2 VI 3 pseudo-binary compounds, intrinsic (native) point defects refer to vacancies and antisite defects, the formation of which intimately depends on the stoichiometry [5][6] . Excess Te in the starting materials of synthesis facilitates the formation Ge and Sb vacancies on random cationic sites, whereas deficient Te would facilitate anion antisites. Indeed, the color gradient of Te 1 layers indicates a large amount of anionic antisites resulting from the off-stoichiometry (Table S2). Meanwhile, Sb Te " and Ge Te "" anionic antisites responsible for the observed high carrier concentration. We calculated the defect concentration from the measured carrier concentration of 4.3 × 10 20 cm -3 to be about 2.25% Sb Te antisites defects or 1.13% Ge Te antisites defects as tolerable and energy favorable ( Figure S18). From the cation normalized map we conclude that the Ge is rich in the middle layer than the outer layers. Beyond the normalized intensity, we further carry the quantitative analysis of the cation distribution down to the atomic level. As known, the brightness of the HAADF-STEM image is roughly proportional to the exponential function of the atomic number. Here we suppose Ge occupancy + Sb occupancy = 1, and the cation vacancy is left out due to the small amount. Then we calculate the quantitative results of the cation distribution using Equation 6: where a is the Ge percentage in each cation column, , and are the intensity, , and are the atom mumber of Ge, Sb and Te, respectively. The / ratio is the normalized intensity value. The normalized intensities obey a Gaussian distribution, so the mean value is identified as 50% Ge-occupancy and 50% Sb-occupancy based on the literature. [7][8][9][10][11][12] The exponential α is found to be 1.74, which is different from the square of atomic numbers as commonly used. The scattering cross-sections are highly sensitive to collection angle [13] , surface strain fields [14] , the number of atoms in a column and the chemical composition [15] , the real exponent will be lower than expected. The darkest atomic column corresponds to about 82.0% Ge-occupancy. In our IDA method, the ratio of the normalized intensity was used to calculate the quantitative results of the cation distribution as following: . So, the intensity of incident electron probe did not influence the results in the quantitative process.
The curve of the Ge percentage as a function of normalized intensity is shown in Figure  S5. The point of intersection with vertical axis represents the normalized intensity as 100% Ge occupancy. Similarly, the point of intersection with horizontal axis represents 100% Sb. Based on the normalized intensity value of each atom column, we can map the cation distribution (Figure 1i and Figure 1j) in the main text. The Ge occupancy of each atom column is thus derived. Figure S6 shows the statistics data of atomic column Ge occupancy of GS1 and GS2 cationic sublattice sites. Finally, we obtain the quantitative analysis of the cation distribution to reach an unprecedented atomic level. The average model of our results is -Te-Ge 26.9 /Sb 73.1 -Te-Ge 49.9 /Sb 50.1 -Te-Ge 26.9 /Sb 73.1 -Te-, a statistical result over 7000 atoms columns. Likewise, we can quantify the anion distribution as well, assuming 2.25% Sb Te or 1.13% Ge Te antisites defects.

Table S1
The atomic number of Ge, Sb and Te, respectively.

Table S2
The chemical content of the ideal Ge 1 Sb 2 Te 4 , as well as the real single crystal GST measured by the super X-ray energy dispersive spectrometer (EDS).

Table S3
The concentration of the defects when the carrier concentration reaches 4.3×10 20 cm -3 .

Figure S5
The Ge percentage in the cation (Ge/Sb) column as a function of normalized intensity.

Table S4
The fitting values of Gauss distribution of both GS1 and GS2 cationic columns.

Figure S6
The statistics data of atomic column Ge occupancy of GS1 and GS2 cationic sublattice sites.

The source-of-error analysis the IDA method
The statistical histograms of Te normalized intensity over 4000 atom columns show that Te 2 obeys the Gaussian distribution while Te 1 deviate that in Figure S7. The full width at half maximum (FWHM) of the intensity distribution is composed of the intrinsic atom disorder distribution fluctuation and instrument resolution errors, F(i, q) = D(i)  E(q), (7) where the symbol "" represents convolution, F(i), D(i) and E(q) present the fluctuation sum, disorder part and instrument part, respectively. The instrument resolution error bar of TEM and IDA in quantification analysis is estimated to be about 4%, evaluated on the FWHM of the Te2 peak shown in Figure S7.

Figure S7
The statistical histogram of Te normalized intensity over 4000 atom columns. The error bar of TEM and IDA in quantification analysis, evaluated on the full width at half maximum (FWHM) of the Te2 peak, is estimated less than 4%.

Simulated HAADF STEM images
HAADF imaging gives highly intuitive images and can often be directly interpreted on a semi-quantitative level (Z-contrast images), however, the combined elastic and thermal diffuse scattering (TDS) make the background somewhat complicated, especially when we try to extract the quantitative information of the sample at the atomic scale. [16] Image simulations are applied for a comparison of experiment and simulation on the same absolute intensity scale with atomic column resolution. The effects of dynamical scattering of the probe and also the contribution of thermal diffuse scattering must be included in the simulation. Meanwhile, the effects of channeling can be significant, because atomic resolution HAADF STEM is often performed on crystals aligned along high-symmetry directions where the atoms form resolvable columns. [16] So far, the approaches to including both channeling and TDS in the calculations fall into two main categories: those that use absorptive potentials and those that make use of the frozen phonon approach.
In this section, some of the detailed simulation results are presented, which used the Dr. Probe® simulation software [17] , providing quantitative HAADF-STEM image simulations at the atomic scale. Five atomic structure models are input for the simulations, example images are displayed in Figure S8 for different thicknesses of a GeSb 2 Figure S8. The HAADF images have obvious difference in the different models, and gradually change with thickness. The sample thickness and Ge/Sb site disorder have an effect on the intensity. Through careful analysis, Figure S9 shows a comparison of experiment and simulation on the same absolute intensity scale with atomic column resolution. We randomly pick and extract the intensities of three septuple-structures using the IDA method, the results agree with the Matsunaga model (-Te-Ge 25 /Sb 75 -Te-Ge 50 /Sb 50 -Te-Ge 25 /Sb 75 -Te-) (cf. the first section of Supporting Materials, the "nearest" neighbor atomic column intensity differential analysis).

Figure S9
Comparison of experiment and simulation on the same absolute intensity scale with atomic column resolution.
Besides, we analyze the exponent α of simulation images by IDA. HAADF imaging gives highly intuitive images and can often be directly interpreted on a qualitative level (Z-contrast images). The intensity of scattering would follow the Rutherford scattering model of being proportional to Z 2 where Z is the atomic number of the illuminated atom. The Rutherford scattering model assumes an unscreened Coulomb potential, but in practice screening will reduce the Z-dependence to some extent. Simulations of scattering intensities using realistic atomic scattering factors suggest that a slightly-lower-than-2 exponent (1.6-1.9) is more reasonable for thin samples. Due to the effects of dynamical scattering of the probe, the channeling effect (Atomic resolution HAADF STEM is often performed on crystals aligned along high symmetry directions where the atoms form resolvable columns), thermal diffuse scattering (TDS) and sample thickness, the real exponent will be much lower than expected. [18,19] We randomly pick and extract the intensities of three septuple-structures using the IDA method, the results agree with the Matsunaga model (-Te-Ge 25 /Sb 75 -Te-Ge 50 /Sb 50 -Te-Ge 25 /Sb 75 -Te-) (cf. the first section of Supplementary Materials, the "nearest" neighbor atomic column intensity differential analysis). Besides, we also analyze the exponent α of simulation images by IDA. In our simulation of GeSb 2 Te 4 , the exponent α of simulated images can be as high as about 1.98 for thin samples, and as low as about 1.89 for thicker sample (about 50 nm), which is very close to the value used in our IDA result (1.74). Above all, the IDA method proposed in Section 1 allows us to statistically extract and quantify the site disorder of HAADF-STEM images, which are comparable with the reference of multislice simulations given by Dr. Probe software.  Figure S11 XRD measurement of GST at different temperature (from room temperature to 773 K) using PANalytical Empryean with CuKα radiation.

Calculations of configurational entropy
The Boltzmann's formula of configurational entropy stemming from the occupation of any microstate is assumed to be equally probable for an isolated system at a global thermodynamic equilibrium, [20] which is in contrast to our findings that Ge/Sb presents Gauss distribution in microscopic scale as shown in Figure 1 (cf. main text). The data of atomic site disorder presented in the preceding section allow us to statistically extract and quantify the configurational entropy. Here, we define the configurational entropy of i th atomic column as based on the concept of entropy and expressed in terms of a discrete set of probabilities [21] − ∑ log , ∑ , We can now substitute Equation 11 into Equation 9 to calculate the configurational entropy. variance, respectively, 0<f<1, its configurational entropy will become The statistical histograms of configurational entropy of t-GST A wave-like vibration mode of the CE along the a-axis of GS1 and GS2 lattices are revealed. The pair correlation function p(r) can be used to describe the fluctuation period. [22] The discrete data points of GS1 layers and GS2 layers are used to calculate the pair correlation function by the self-developed MATLAB code, the defined correlation function can be expressed as ∑ , where N is the total number of pairs of and .
present the configurational entropy values at the distance range of , where r is the spacing between any two CE peaks in a particular atomic layer. The smooth curves of GS1 layers reveal a high frequency CE modulation (1/3 Å -1 ) and the short-range order of Ge/Sb, and the characteristic length of modulation of atomic column CE of GS1 layers is estimated to be close to the interatomic spacing ~3 Å. By contrast, a sharp distinctive feature occurs in the GS2 layers, and the characteristic length for GS2 layers are found to be at about 9 Å with a low frequency CE modulation (1/9 Å -1 ). Extended to the data containing 60  50 (layers) atom columns, the CE of GS1 layers and GS2 layers are shown in Figure S12

Figure S12
The entropy statistical histograms of GS1 and GS2 cationic sublattice sites over 3000 cation columns.

Configurational entropy induced phonon scattering
The disorder scattering parameter Γ, where the scattering parameters Γ and Γ are due to mass and strain field fluctuations, respectively. Γ Γ Γ , (14) The chemical composition of a material can be expressed as A 1c1 A 2 c 2 A 3c3 A 4c4 …A ncn , where the A i are crystallographic sublattices in the structure and the c i are the relative degeneracies of the respective sites. In general, there will be several different types of atoms that occupy each sublattice, and the k th atom of the i th sublattice has mass , radius , and fractional occupation . The average mass and radius of atoms on the i th sublattice are: The average atomic mass of the compound is: The mass fluctuation scattering parameter is then given by: where the mass fluctuation scattering parameter for the i th sublattice is: For two different atoms on each of the i th sublattices, i.e., k=1, 2. In the case of t-GST, n=2 because of 2 sublattices. We have masses and , fractional concentrations and , and and ̅ , Equation 18 becomes The strain field fluctuation scattering parameter is then given by: Where ̅ , and is a phenomenological adjustable parameter for the i th sublattice. The parameter is a function of the Grüneisen parameter γ, which characterizes the anharmonicity of the lattice. [23][24][25][26][27][28] Sound velocity is another important parameter to understand the thermal transport in a solid. The measured sound velocities and series calculation values present in Table S5. The longitudinal (ν l ), transverse (ν t ), and mean (ν m ) sound velocities for the GST along c axis are 2700 m s ￚ 1 , 1773 m s ￚ 1 and 1941 m s ￚ 1 , respectively. The phonon mean-free path of GST is estimated to be 0.36 nm according the measured sound velocities. The experimental values are close to the amorphous limit (κ Lmin ) estimated by the Cahill approximation model [29] .

The sound velocities and related mechanical properties
Measurements of sound velocities are carried out on the pellet samples at room temperature. Longitudinal and transverse sound velocities are measured using pulse-receiver (®Olympus-NDT) equipped with an oscilloscope (®Keysight). Shear gel (®Olympus) and water are used as couplants between the sample and the ultrasonic transducers for transverse and longitudinal sound velocity measurements, respectively. Average sound velocity ( ) is calculated from the longitudinal sound velocity ( ) and the transverse sound velocity ( ). ) are calculated by the Cahill approximation model [24] .

Density function theory calculations Structure models
Here we compare two systems in terms their thermodynamic stability. (1) The pristine ordered GST: Ge ions occupy the GS1 layer while Sb ions occupy the GS2 layer. The other system is disordered GST, there is cross-layer site occupational disorder, sometimes called "self-solid-solution". The disordered GST is simulated using a Special Quasi-random Structure (SQS) [29] and the Alloy Theoretic Automated Toolkit (ATAT) [31] , which satisfies the mixing ratios of Ge/Sb, (2) Ge : Sb = 1 : 1 at the GS1 layer, and Ge : Sb = 1 : 3 at the GS2 layer, which is a high entropy state.

Figure S13
The structure models of GST in the pristine order and high entropy disorder system, Te (cyan), Sb (purple), Ge (pink).

Figure S14
The electron localization function (ELF) and bonding lengths of GST in the pristine order and high entropy disorder system.
Since the Ge-Sb-Te system is a layered structure and contains heavy elements (Sb and Te), we need to test the effects of van der Waals (VdW) interactions and spin-orbital (SO) coupling on the energetic trend. The VdW corrections is used the Tkatchenko-Scheffler (TS) method [32] . We calculate the energy differences of the two compounds (pristine order and high entropy) with the VdW, SO or VdW+SO corrections, and compare them with the simple PBE results (Figure S15). It turns out both the stability trend and the energy differences using the corrections are similar to those using PBE. This is due to that the three Ge-Sb-Te compounds share the same layered structure and the same elements. This confirms that the PBE calculations are suitable to investigate the phase stability trend of the two compounds.

Figure S15
Theoretically calculated energy differences of the pristine order and high-entropy Ge-Sb-Te compounds using the PBE functional, SO, VdW and VdW+SO corrections.
The geometry structure parameters (aka the lattice constants and interlayer distances) of the pristine GeSb 2 Te 4 using the PBE functional and the VdW corrections are presented in Table S6. We find that the in-plane lattice constant is the same with or without the VdW corrections. While along the c direction, the VdW corrections shrink the interlayer distance and the lattice constant along c axis. Smaller interlayer distance implies stronger interlayer interactions and thus higher lattice thermal conductivity along the c direction. Nevertheless, we notice that the pristine and high entropy configurations undergo similar trend of the interlayer distance change after including the VdW corrections (Table S6). Therefore, whether including VdW corrections or not won"t alter our conclusions on the lattice thermal conductivity.  Table S6 The geometry structures of the pristine and high-S GeSb 2 Te 4 compounds using the PBE functional and VdW corrections.
We use the Deybe model to evaluate the vibrational properties [33,34] . The enthalpy (H vib ) and vibrational entropy (S vib ) can be written as, where, , v, ħ and k B are the frequency, phonon velocity, Boltzmann constant, and Plank constant, respectively.  D is the Debye frequency ( D = k B /ħ,  is the Debye temperature).
Further including the configurational entropy (S conf ), we can evaluate the Gibbs free energy as formula − − , (29) The configurational entropy plays a role of stabilizing the t-GST phase by theoretical calculations of Gibbs free energy in Figure 4c of main text.

The electronic band structure
The calculated band gap of pristine order model and high entropy model are ~ 0.33 eV and ~ 0.14 eV, respectively.

Figure S16
The calculated density of states of the pristine order and high entropy models of GST.

Formation energy of various point defects
In Ge 1 Sb 2 Te 4 alloys, inherent point defects occur during the crystal growth from the stoichiometric melt. A large amount of anion antisite defects result from the off-stoichiometric feature (Table S2). Negatively charged Sb Te " and Ge Te "" anion antisite defects is reasonable upon the deficiency of Te atoms, which accounts for the p-type conductivity. The smaller the differences in electronegativity and atomic size between cation and anion atoms, the less resistance to the formation of antisite defects. [35] Figure S17 The formation energy of various point defects calculated by the density functional theory (DFT), which indicates the formation of anion antisite defects and cation disorder is more thermodynamically favorable than the vacancy in GST.

Lattice distortion effects
The lattice distortion of the local octahedron motif of (Ge)Te 6 in the disorder GST system is more severe than that in order GST. The details of bonding lengths and bonding angles present in Figure S17. in the pristine order and high entropy GST system.

Anharmonicity
Since the Grüneisen parameter characters the lattice anharmonicity (the relationship between phonon frequency and crystal volume change), low Young"s modulus, shear modulus, and large Gruneisen parameter usually imply a low lattice thermal conductivity. Table S7 shows the calculated bulk and shear moduli (B and G in GPa), phonon velocities (ν in m/s), Grüneisen parameters (γ), and acoustic Debye temperatures (Θ D in K) of both pristine order and high entropy GST systems. The calculated results show that the high entropy (GST-SQS) has lower phonon velocity and larger Grüneisen parameter than that of the pristine order structures.

The thermal conductivity calculation
The relaxation time of the cationic location disorder scattering. [36,37] Where is the phonon frequency, is the average volume per atom, is the acoustic phonon velocity, is the mass fluctuation and and is the strain fluctuation due to the Ge/Sb disorder.
The relaxation time of the cationic location disorder scattering only considering the mass fluctuation is given by [37] : (31) Where k B is the Boltzmann constant, ħ is the Plank constant, T is the absolute temperature, ħ / , is the mass of the doping atom, is the atomic occupation ratio, and ̅ is the average mass of the unit cell.
The r-GeSb 2 Te 4 has high site occupational disorder and thus high configurational entropy. It is challenging to derive the 2nd and 3rd order force constants required for the ShengBTE® code [38] , not to mention the computation cost of computing anisotropic transport matrices. The Debye-Callaway model has drawbacks, e.g., the negligence of optical modes, may result in deviations of the absolute magnitude of lattice thermal conductivity. Nonetheless, acoustic modes are of utmost importance in the thermal transport. The efficacy of the Debye-Callaway model has been proved in numerous studies, e.g., Cu-Sb-Se [39] , SnSe [40] and a layered structure Ge 2 Sb 2 Te 5 [41] . The Debye-Callaway model herein adopted include the effect of mass fluctuation (cf. Eq. 4 in main text and Eq. 31 in SI). For the purpose of unraveling the relation between entropy and lattice thermal conductivity in r-GeSb 2 Te 4 with high native site disorder, the Debye-Callaway model correctly explains the trend of variation of lattice thermal conductivity and gives semi-quantitative yet microscopic insights, compared to other existing models.
We additionally provide the geometry structures (lattice constants and interlayer distances) of the pristine GeSb 2 Te 4 using the PBE functional and the VdW corrections (Table S6). We find that the lattice constant in plane is the same with or without the VdW corrections. While along the c direction, the VdW corrections shrink the c axis and the interlayer distance. The smaller interlayer distance might indicate the stronger interatomic interactions and higher lattice thermal conductivity. Nevertheless, we notice that the three compounds (pristine and high-S) show the similar trend of the interlayer distance change with including the VdW corrections (Table S6). Therefore, this obviously will not change our conclusions on the lattice thermal conductivity.

Figure S19
The calculated lattice thermal conductivity of order and disordered GeSb 2 Te 4 using the elastic properties. (a) The red line presents the lattice thermal conductivity of disordered GeSb 2 Te 4 only considering the ideal mass fluctuation. (b) The lattice thermal conductivity of pristine order and high entropy state of GeSb 2 Te 4 .

Differential thermal analysis (DTA)
Figure S20 Differential thermal analysis (DTA) results indicate that GST appears to melt at 889 K with one endothermic peak observed on the heating curve.

Hall mobility
The carrier concentration is high at a level of about 4.3×10 20 cm ￚ 3 at room temperature, which is relevant to the large density of native point defects. The Hall mobility exhibits an exponent behavior with temperature. A relatively weak temperature dependence of mobility presents at low temperature, which may be relative to configurational entropy. And acoustic phonon scattering is the dominant at higher temperatures.

Figure S21
Hall mobility as a function of temperature for GST single crystals along a-axis.

Lorenz number and thermal diffusivity
For the Ge 1 Sb 2 Te 4 sample, the Lorenz number is calculated based on a single-band Kane model (SKB) [42,43] for the high carrier concentration, and we assume the charge carriers are mainly scattered by acoustic phonons. A room temperature κ L as low as 0.4 W m -1 K -1 is attained along c-axis of value of t-GST, comparable to other state-of-the-art single crystalline thermoelectric materials that are claimed to have ultralow thermal conductivity, e.g. SnSe [40] , SnS [44,45] and SnSe 0.9 S 0.1 [43] (cf.

Figure S25
The lattice thermal conductivity (κ L ) of GST compared with SnSe, SnS and SnSe 0.9 S 0.1 thermoelectric single crystals.
We compare the lattice thermal conductivity of t-GST, Sb 2 Te 3 [46,47] and GeTe in Figure   S26 and Figure S27, respectively. It is found that the lattice thermal conductivity of t-GST is lower than that of Sb 2 Te 3 , which had a sharp peak in temperature dependence. However, the peaks of t-GST are flatter, illustrating the significant influence of configurational entropy on thermal conductivity.

Figure S26
The total thermal conductivity and lattice thermal conductivity for Sb 2 Te 3 [46,47] and GST single crystals.

Figure S27
The reported available data of thermal conductivities of GeTe in the literature of the past decade. The solid black squares and red squares represent the total thermal conductivity and lattice thermal conductivity, respectively.

Electrical transport and thermoelectric performance
In general, entropy is intimately related to thermoelectrics: a thermoelectric process is a coupled process between charge and entropy flow (including heat flow and entropy creation); the nature of Seebeck coefficient is the average entropy carried per each charge carrier. The performance of thermoelectric materials can be expressed by the dimensionless figure of merit zT=α 2 σT/κ, where α, σ, T and κ are the Seebeck coefficient, the electrical conductivity, the absolute temperature and the total thermal conductivity, respectively. In this section, we present the electrical transport and thermoelectric performance of single crystal t-GST in details, as shown in Figure S28. Electrical conductivities (~10 5 S m -1 ) show a decrease tendency with rising temperature in Figure S28a, suggesting a degenerate semiconductor behavior, and σ along a-axis is almost two times higher than that along c-axis. The maximum S value reaches up to 151.5 μV K -1 at 773 K along c-axis in Figure S28b, because only high energy holes could pass over the energy barrier induced by van der Waals gaps. The strong anisotropy in electrical conductivity and Seebeck coefficient is probably due to the difference in carrier mobility along two directions, which is analogous to V 2 VI 3 compounds [48] . The bipolar effect does not explicitly appear in the curves of temperature-dependent σ and α due to the high hole concentration. The power factor (α 2 σ) is displayed in Figure S28c. The highest power factor up to 1.34 mW m -1 K -2 is observed along c-axis at 723 K, which is higher than that along a-axis (1.1 mW m -1 K -2 at 773 K). The zT value as high as 0.8 at 723 K could be achieved along c-axis in Figure S28d, which is 8 times larger than the previously reported value of polycrystalline GST samples. [49,50] We also compare the thermal conductivity and thermoelectric zT of our single crystals t-GST with literature [50][51][52] in Figure S29. These results should be of interest to explore new TE materials in ⅣⅤ-V 2 VI 3 pseudo-binary compounds, such as Ge 2 Sb 2 Te 5 , Sn 1 Sb 2 Te 4 , Sn 1 Bi 2 Te 4 and Pb 1 Bi 2 Te 4 et al.  Figure S29. Comparison of the thermal conductivity (a) and thermoelectric performance zT (b) as a function of temperature for t-GST crystals. [50][51][52]