Kanka Ghosh

Andrzej Kusiak

Jean-Luc Battaglia

Significant four-phonon scattering and its heat transfer implications in crystalline Ge 2 Sb 2 Te 5

We systematically demonstrate the temperature-dependent thermal transport properties in crystalline Ge2Sb2Te5 via first-principles density functional theory-informed linearized Boltzmann transport equation. The investigation, covering a wide temperature range (30 K-600 K), reports the emergence of an unusual optical phonon-dominated thermal transport in crystalline Ge2Sb2Te5. Further, a significant contribution of four-phonon scattering is recorded which markedly alters the lattice thermal conductivity. Therefore, the combined effect of cubic and quartic phonon anharmonicity is seen to navigate the underlying physical mechanism and open up intriguing phononic interactions in Ge2Sb2Te5 at high temperature. Irrespective of three and four-phonon processes, Umklapp is seen to prevail over normal scattering events. Consequently, four-phonon scattering is found to notably reduce the lattice thermal conductivity of Ge2Sb2Te5 to 28 % at room temperature and 42 % at higher temperature. This quartic anharmonicity further manifests in the breakdown of T -1 scaling of thermal conductivity and challenges the idea of a universal lower bound to phononic thermal diffusivity at high temperature. The faster decay of thermal diffusivity compared to T -1 is rationalized encompassing the quartic anharmonicity via a modified time scale. These results invoke better understanding and precision to the theoretical prediction of thermal transport properties of Ge2Sb2Te5. Concomitantly, this also triggers the possibility to explore the manifestations of the lower bound of thermal diffusivity in materials possessing pronounced four-phonon scattering.

I. INTRODUCTION

Phase change memory devices (PCM) represent an efficient nonvolatile memory storage mechanism with thermally induced rapid transitions between crystalline and amorphous phases, featuring a wide range of applicability [1][2][3]. Amongst these PCM materials, Ge 2 Sb 2 Te 5 (GST) is one of the most suitable and promising candidates due to its fast switching speed and notable contrast in electrical and optical properties during phase transition [1,[4][5][6][7][8]. Recently, Ge-rich GST materials were also being investigated to study their potency in PCM applications [9][10][11]. Being such a convenient PCM material, the thermoelectric properties of Ge 2 Sb 2 Te 5 [12][13][14][15][16][17][18] have been the topic of intense investigations over the years. Ge 2 Sb 2 Te 5 exists in one amorphous [19] and two crystalline phases (cubic and hexagonal). Owing to its low thermal conductivity (κ), crystalline Ge 2 Sb 2 Te 5 was found to be suitable for high temperature thermoelectric applications due to its enhanced thermoelectric figure of merit (ZT = S 2 σT /κ, σ: electrical conductivity, S: Seebeck coefficient) [14]. Therefore, scrutinizing high temperature thermal transport is essential to provide better understanding of its thermoelectric behavior. Over the years, experimental [20][21][22][23][24][25] as well as state-of-the-art ab initio theoretical methods [15,[26][27][28][29][30] were successfully employed to predict and measure thermal conductivity of crystalline Ge 2 Sb 2 Te 5 . While Wiedemann-Franz law was quite successful in reproducing the electronic part [20,31], phonon contri- * kanka.ghosh@lspm.cnrs.fr bution to thermal conductivity of hexagonal Ge 2 Sb 2 Te 5 needs special attention as an unusual domination of optical phonons was observed [28]. Contrary to the conventional heat transport mechanism via acoustic phonons, optical phonon-dominated thermal transport is less common. Also, this feature seems to challenge the conventional Slack model [32,[START_REF] Kaviany | Heat transfer physics[END_REF] which takes only acoustic phonons into account. Optical phonons gain importance due to various features like large group velocities, anharmonic effects (2H MoS 2 [START_REF] Dong | [END_REF]), and large acoustic-optical band gap (AlSb [35]), among a few. The large acousticoptical band gap in AlSb and GaN imposes restriction to the acoustic-optic phonon scattering, resulting in a large thermal conductivity [36]. At higher temperature, optical phonons populate without significant resistances due to forbidden acoustic-optical channels. This is pivotal in sustaining the optical mode-dominance in κ L in these high-κ L materials. On the other hand, for lowκ L materials, moderate contributions of optical phonons (22 % of total κ L for T = 300-700 K for PbTe [37] and ∼ 30 % at T = 300 K for Mg 2 Si [38]) were seen to stem from relatively dispersive optical branches and their presence at lower frequencies [35]. The softening of the optical modes in PbTe was identified to be crucial to trigger acoustic-optical scattering channel to produce low κ L [37]. Although featuring acoustic phonon-dominated thermal conduction, non-negligible contributions of optical phonons were observed in a temperature range of 4K-600K (∼ 80 % acoustic and ∼ 20 % optical modes) for GeTe [39,40]and at T = 300 K (∼ 65 % acoustic and ∼ 35 % optical modes) for Sb 2 Te 3 [29]. First principles study by Mukhopadhyay et al. [28] revealed that around 75 % of κ L was contributed by the opti-cal modes in stable hexagonal crystalline Ge 2 Sb 2 Te 5 at room temperature. Also, a recent investigation focused on the coupling between acoustic and optical branches in the phonon spectrum of crystalline Ge 2 Sb 2 Te 5 [17]. This coupling between acoustic and optical phonon branches and the dominance of optical phonons at lower frequencies can effectively trigger higher order phonon scattering channels which is our particular interest in this article. The strong presence of four-phonon scattering explains the strong temperature dependent thermal conductivity, although the anharmonic renormalization method was also employed [41] to account for strong temperature dependent κ L . Including four-phonon processes was found to affect the optical phonon contributions to κ L for AlSb and BAs [36,42], diminishing the optical phonon dominance in thermal transport. In low-κ L material, the presence of lower frequency optical phonons, their strong contribution to κ L as well as coupling between acoustic-optical branches can incite four-phonon processes as temperature is raised. Recently, low-κ L ionic conductors AgI [43] and AgCrSe 2 [44] were found to display strong quartic anharmonicity characterized by flat phonon dispersion and rattling dynamics of weakly bonding atoms. CsPbBr 3 , a prominent low-κ L chalcogenide perovskite material, was also shown to feature reasonable four-phonon scattering [45]. However, the implication of the presence of four-phonon scattering on thermal transport in a broad variety of low-κ L materials is intriguing and not well studied. In this spirit, we investigate the phononic behavior of low-κ L chalcogenide Ge 2 Sb 2 Te 5 to test its higher order scattering features. Controlling thickness [46][47][48] as well as vacancies [49] or stoichiometric defects [50] was proved to be efficient in modulating thermal transport in crystalline Ge 2 Sb 2 Te 5 . A small amount of Ge, Sb vacancies and disorder were found to lower κ L around 70-80 % [29]. However, as a standard practice, computing κ L and comparing it with experimental measurements involve calculating threephonon scattering along with approximate models to incorporate disorder, vacancies, and grain-boundary scattering perturbatively. Hence, to reproduce experimental results, often higher order phonon scattering is ignored which, if significant, can overestimate the κ L coming from three-phonon scattering alone. Moreover, at smaller length scales of the sample, crystals can be found to exist in a defect-starved pristine state. Consequently, it is important to precisely predict κ L that solely comes from phonon scattering. In this context, we focus our investigation on the phonon transport in crystalline hexagonal Ge 2 Sb 2 Te 5 to decipher the role of higher order phonon scattering mechanisms in a wide temperature range. Therefore we have several goals in this paper: (a) To check, validate and introspect the dominance of optical phonons in the thermal transport of crystalline Ge 2 Sb 2 Te 5 as a function of temperature. (b) To study whether four-phonon scattering is significant to describe the thermal transport in crystalline Ge 2 Sb 2 Te 5 at high temperature. (c) If significant, to investigate whether different mode contributions are being influenced due to the inclusion of four-phonon scattering. (d) To gain precise predictions of the phonon-phonon scattering term in κ L for a better estimation of other scattering terms (effect of disorder, grains) in crystalline Ge 2 Sb 2 Te 5 . To answer these issues, we explore the temperature dependent thermal transport in crystalline Ge 2 Sb 2 Te 5 , ranging from T = 30 K to 600 K using first-principles density functional method coupled with the solution of the linearized Boltzmann transport equation (LBTE). A systematic study is carried out to investigate harmonic and anharmonic phononic properties of GST. The presence of significant four-phonon scattering has been identified and the relative competitions between different scattering channels are discussed. Lattice thermal conductivity and its mode-wise contributions are evaluated with and without four-phonon scattering to mark the contrast. The effect of higher order phonon scattering on the anisotropy in thermal conduction has also been studied. Finally, the deviation from high temperature T -1 scaling of lattice thermal conductivity is analyzed. This resolved picture of thermal conduction attempts to bring more precision to the microscopic understanding for Ge 2 Sb 2 Te 5 , which can improve thermal management of related chalcogenide PCM materials.

II. METHOD: DENSITY FUNCTIONAL THEORY-INFORMED BOLTZMANN TRANSPORT EQUATION

The structural parameters of the crystalline, hexagonal phase of Ge 2 Sb 2 Te 5 (space group P 3m 1 ) are optimized via first-principles density functional theoretical (DFT) calculations, and the obtained lattice parameters in comparison to other related works are presented in Table I. Phonon dispersion and phonon density of states (PDOS) of two different crystalline, hexagonal Ge 2 Sb 2 Te 5 structures (Kooi and Petrov (see supplementary Fig S4 for PDOS and dispersion relation for Petrov structure)) have been calculated using the density functional perturbation theory (DFPT) [51] as employed in the quantum-espresso [52] suite of programs. Self-consistent field calculations, within the framework of density functional theory (DFT), are carried out to compute the total ground state energy of the crystalline Ge 2 Sb 2 Te 5 . For this purpose, Perdew-Burke-Ernzerhof (PBE) [53] generalized gradient approximation (GGA) is used as the exchange-correlation functional. Electronion interactions are represented by pseudopotentials using the framework of projector-augmented-wave (PAW) method [54]. The Kohn-Sham (KS) orbitals are expanded in a plane-wave (PW) basis with a kinetic cutoff of 35 Ry (476.2 eV) and a charge density cutoff of 140 Ry whereas for relax calculation the plane-wave (PW) kinetic cutoff is taken as 80 Ry (1088.45 eV). The semiempirical correction due to Grimme (D2) is considered to incorporate the van der Waals (vdW) interactions. The Brillouin zone (BZ) integration for self consistent electron density calculations are performed using a 12×12×12 Monkhorst-Pack (MP) [55] q-point grid.

We compute the lattice thermal conductivity and associated phonon scattering rates employing the Boltzmann transport equation (BTE). A linearized version of the BTE is solved for phonons using an iterative selfconsistent method as well as relaxation-time approximation (RTA) approach via fourphonon [56], an updated version of the shengbte code [57]. The numerical solution to the phonon-BTE demands accurate estimates of harmonic and anharmonic force constants. Therefore, DFT-based self-consistent field calculations are carried out to output harmonic and anharmonic force constants. Harmonic force constants are obtained using the DFPT method. We use a 4×4×2 supercell of Ge 2 Sb 2 Te 5 within the finite displacement approach to calculate the thirdorder interatomic force constants using finite differences method via thirdorder package [57], allowing interactions up to the fourth-nearest neighbors. We adopt same parameter values as that of the work by Ibarra et al. [15] for getting the third order interatomic force constants (IFC). 2×2×1 supercells with interactions up to second-nearest neighbors are used to obtain fourthorder IFCs from 1728 different configurations. Within the relaxation-time approximation (RTA), phonon lifetimes (τ RT A λ ) take the following form

1 τ RT A λ = 1 τ 3ph λ + 1 τ 4ph λ + 1 τ I λ (1)
with 1

τ 3ph λ = 1 N + λ ′ λ ′′ Γ + λλ ′ λ ′′ + - λ ′ λ ′′ 1 2 Γ - λλ ′ λ ′′ (2)
1

τ 4ph λ = 1 N ( ++ λ ′ λ ′′ λ ′′′ 1 2 Γ ++ λλ ′ λ ′′ λ ′′′ + +- λ ′ λ ′′ λ ′′′ 1 2 Γ +- λλ ′ λ ′′ λ ′′′ + -- λ ′ λ ′′ λ ′′′ 1 6 Γ -- λλ ′ λ ′′ λ ′′′ ) (3)
and

1 τ I λ = 1 N I λ ′ Γ I λλ ′ ( 4
)
where N is total number of q points, + and -denote absorption and emission respectively corresponding to three-phonon scattering. ++, +-and --represent four-phonon recombination, redistribution and splitting processes respectively. Three-phonon, four-phonon and isotope scattering rates are represented via Γ λλ ′ λ ′′ , Γ λλ ′ λ ′′ λ ′′′ and Γ I λλ ′ respectively. While fourth order IFCs are only implemented in RTA, third-order IFCs are computed using both RTA and iterative method. For a specific mode λ, linearized phonon-BTE can be expressed as

F λ = τ RT A λ (v λ + ∆ λ ) (5)
with

∆ β λ = 1 N [ + λ ′ λ ′′ Γ + λλ ′ λ ′′ (ξ β λλ ′′ F β λ ′′ -ξ β λλ ′ F β λ ′ )+ - λ ′ λ ′′ 1 2 Γ - λλ ′ λ ′′ (ξ β λλ ′′ F β λ ′′ + ξ β λλ ′ F β λ ′ )] (6)
Here, ∆ λ accounts for the iterative process and retrieves RTA with ∆ λ = 0.

ξ β λλ ′ = ω λ ′ v β λ ′ / ω λ v β λ
, where v λ denotes modal group velocity and F λ represents the linear deviation of the phonon population (f λ ) from the equilibrium Bose-Einstein distribution (f 0 λ ) as

f λ = f 0 λ -F λ • ∂f λ
∂T ∇T -... . Finally, using third-order and fourth-order IFCs and knowing the phonon distribution function, lattice thermal conductivity tensor (κ αβ L ) is computed from the following expression

κ αβ L = 1 k B T 2 N V 0 λ f 0 (f 0 + 1)(ℏω λ ) 2 v α λ F β λ (7)
where λ defines every phonon mode, N denotes the total number of modes in BZ, V 0 is the unit cell volume, f 0 = (e ℏω/k B T -1) -1 represents equilibrium Bose-Einstein phonon distribution function. In the RTA approach, F β λ = τ β λ v β λ , where τ β λ stands for the relaxation time for phonons of polarization λ propagating along the β direction. Phonon-BTE for both three and four-phonon processes are solved in a 16 × 16 × 4 q-grid over the Brillouin zone (see convergence of three and four-phonon processes in Supplementary Fig S15).

III. RESULTS AND DISCUSSIONS

The stable hexagonal phase of crystalline Ge 2 Sb 2 Te 5 possesses a crystal structure of P 3m 1 symmetry with nine atoms per unit cell stacked along the c axis. However, the stacking of these nine atoms along c is still being debated [15,29,[58][59][60][61]. The stacking sequence along the c axis by Petrov et al. [59] (defined as Petrov structure hereafter) was proposed as Te-Sb-Te-Ge-Te-Te-Ge-Te-Sb-while Kooi and De Hosson [58](defined as Kooi structure hereafter) observed the stacking to be Te-Ge-Te-Sb-Te-Te-Sb-Te-Ge-via electron diffraction experiments. Experiments also suggested a disordered phase in Ge/Sb sublattice with layers randomly occupied by Ge and Sb [60]. However, as we focus only on the phononic part of the heat transport, investigating the mentioned disordered structure is out of the scope of this paper. Lately, another stacking arrangement called inverted Petrov structure [61] was also proposed for hexagonal Ge 2 Sb 2 Te 5 , featuring interchanging adjacent Te and Ge atoms. It was noted in earlier studies [15,29] that Kooi structure exhibits lowest energy amongst these differently stacked structures. Therefore, although we characterize both Petrov and Kooi structures (Table I), hereafter all relevant studies are carried out using the Kooi I. The small differences can be attributed to the numerical uncertainties due to the variation of exchange-correlation functionals used in these studies. For example, PBE generalized gradient approximation (GGA) was employed in both Campi et al. [29] and this work, whereas Ibarra et al. [15] used GGA with a revised version of PBE for solids and surfaces [62] (PBEsol) as exchange-correlation functional. We note that the experimental values of lattice parameters of the Kooi [58] structure were measured as a = 4.25 Å and c = 17.27 Å, which are also consistent with our ab initio findings depicted above. After bench-marking our theoretically relaxed structure, we proceed in describing the phonon transport mechanisms via harmonic and anharmonic properties of Ge 2 Sb 2 Te 5 in the following subsections.

A. Harmonic properties and optical-mode dominance

To understand phonon transport, we start our discussion with the phonon dispersion and phonon density of states (PDOS) realizations for hexagonal Ge zone (BZ) portraying a stable structure with phonon frequencies going to zero at the Γ point in the BZ. The phonon frequencies are seen to span a large frequency regime ranging from 0 to 174 cm -1 . We mark some specific features in the phonon bands. Firstly, we observe the presence of low-frequency optical phonon bands starting from ω = 30 cm -1 . These low frequency optical bands are seen to couple with acoustic bands, indicating strong scattering events to reduce thermal conduction. We will get to this point later. As acoustic and optical phonons feature a strong coupling in an intermediate frequency range (30 cm -1 < ω < 60 cm -1 ), distinguishing individual contribution of the phononic bands to the thermal transport seems unrealistic. However, we can define frequency bands that consist of only acoustic, a mixture of acoustic and optical modes and purely optical modes to comprehend the contributions to thermal transport coming from these frequency bandwidths. Figure 1(C) distinguishes these frequency regimes that comprise of purely acoustic (0 < ω < 30 cm -1 ), purely optical (60 cm -1 < ω < 174 cm -1 ) and the mixed modes (30 cm -1 < ω < 60 cm -1 ) and hereafter will be denoted by [28]. At small frequency range, heavier atoms (Sb, Te) are seen to dominate the PDOS. Around ω = 60 cm -1 , the contribution of Ge atoms starts to be significant. From ω = 90 cm -1 till 120 cm -1 , the enhancement of PDOS is mostly contributed by the Te atoms. At frequencies higher than 130 cm -1 , this enhancement dies out and contributions from Ge, Sb, Te become comparable. Relating the elementwise contributions to PDOS with the frequency regimes A, M and O, we observe the following: (a) In regime A (0 cm -1 < ω < 30 cm -1 ), all three elements contribute comparably. (b) Partial DOS of Te atoms starts dominating with a competition between Sb and Ge PDOS in regime M (30 cm -1 < ω < 60 cm -1 ). (c) In regime O (60 cm -1 < ω < 174 cm -1 ), PDOS for Te mostly dominates and controls the total PDOS, except at higher frequencies (ω > 90 cm -1 ), where all contributions become comparable again. At extreme lower and extreme higher frequencies, all three elements show comparable number of excited vibrational modes. In standard Debye theory, frequency-dependent PDOS can be expressed as

PDOS (ω) = k 2 2π 2
1 dω/dk with the assumption of the linear relation between speed of sound (c s ) and the phonon frequencies (ω = c s k). Consequently, lower group velocities (v g = dω/dk) of phonons enable higher density of states as a function of phonon frequencies. On the contrary, dispersive phonon branches enhance the phonon group velocities causing a lower PDOS. Near the center of BZ (Γ point) in regime A, acoustic branches are seen to behave almost in a linear fashion (Fig 1(C)) which is caused by the vibrational modes of mostly heavy Te and Sb atoms along with little less vibrational contribution from comparably light Ge atoms. In the regime M (30 cm -1 < ω < 60 cm -1 ), dispersion relation shows the coupling between acoustic and optic branches as well as the onset of the comparatively flat bands. These flat bands give rise to an enhancement of PDOS in this regime. Sb and Te atoms seem to contribute equally until ω = 50 cm -1 , after which Te modes start dominating and the vibrational modes of Ge atoms become comparable or even greater than that of the Sb atoms (Fig 1(B)). In regime O (60 cm -1 < ω < 174 cm -1 ), two distinct features are observed. For ω < 130 cm -1 , larger contribution of flat phonon optical bands of low v g (Fig 1(C)) is seen to develop a larger PDOS with several peaks, mostly contributed by the vibrations of Te atoms (Fig 1(B)). For ω > 130 cm -1 , optical branches becomes quite dispersive warranting a higher v g and a decrease in PDOS. Small mass difference between Sb, Te (m Sb = 121.75 a.m.u, m T e = 127.6 a.m.u) makes it difficult to feature separated high and low frequency bands as the evolution of the PDOS goes through competitive vibrations between the elements. In contrast, another chalcogenide material GeTe, owing to the large mass difference between its constituents (m Ge = 72.64 a.m.u, m T e = 127.6 a.m.u), was seen [29,39] to possess PDOS with a separated low frequency acoustic bands (due to the vibrations of heavy Te atoms) and high frequency optical bands (due to light Ge atoms). Also, due to their similar masses, distinct contributions from Sb and Te to PDOS mostly rely on their vibrational modes due to the bonding arrangement in Ge 2 Sb 2 Te 5 . In this context, Sosso et al. [64] studied the Raman spectra of hexagonal Ge 2 Sb 2 Te 5 (Kooi structure) and A-type Raman-active modes (specifically A 1g modes) were identified at high frequencies. These A-type modes indicate vibration along the c-axis due to modulation of the outermost Te-Sb bond [64]. The high frequency vibrational modes involving both Sb and Te manifest themselves in comparable and similar phonon DOS at high frequencies (> 125 cm -1 ) for Te and Sb (see Fig 1(B)). To investigate thermal transport, we employ harmonic and third-order anharmonic force constants and therefore the solution of Linearized Boltzmann transport equation (LBTE) takes the form

κ αβ L = 1 k B T 2 N V 0 λ f 0 (f 0 + 1)(ℏω λ ) 2 v α λ F β λ ( 8
)
where κ αβ L is lattice thermal conductivity tensor, λ defines each phonon mode, N is the total number of modes in BZ, V 0 is unit cell volume, f 0 = (e ℏω/k B T -1) -1 denotes Bose-Einstein phonon distribution function, and F is the term in the expansion of the phonon distribution function which is related to the group velocity (v g ). We use both relaxation time approximation (RTA) as well as iterative approach for the solution of LBTE to obtain κ L . Scattering of phonons by randomly distributed isotopes is also incorporated. The isotope scattering rate (1/τ I λ (ω)), was given by Tamura [65] via second-order perturbation theory as

1/τ I λ (ω) = πω 2 λ 2N λ ′ δ (ω -ω ′ λ ) k g k | α W α (k, λ) W * α (k, λ) | 2 (9
) where g k is the mass variance parameter, defined as

g k = i n i 1 - m ik m k 2 ( 10
)
n i is the mole fraction, m ik is the relative atomic mass of ith isotope, m k is the average mass = i n i m ik , and W is a polarization vector. The natural abundance data for different elements [66] are used to obtain the mass variance parameters. Figure 2 , optical modes become the driving force to contribute κ L with higher contributions than both acoustic and mixed phonon modes. To invoke transparency to this feature, we recall the RTA approach, where the solution of LBTE is reduced to

κ L = 1 N V 0 λ C λ v 2 g λ τ λ (11)
which is obtained via using Eq.8 with zeroth-order approximation to F as F λ = τ λ v gλ and expressing modal heat capacity

C λ = 1 k B T 2 λ f 0 (f 0 +1)(ℏω λ ) 2 .
Increasing temperature does not seem to change the group velocity (v g ). Similar behavior is also seen for heat capacity as it is known to follow Dulong-Petit law at high temperature (see Supplementary Fig S3). However, the variation of phonon group velocities with frequency, shown in Fig 2(c), indicates significant contribution from the optical phonons with high group velocities comparable to that of the acoustic phonon modes. Figure 2(c) also presents v g coming out of three qualitatively different frequency regimes (Region I (A), Region II (M), Region III (O)). To elucidate this point further, probability distributions of v g , corresponding to optical as well as all the modes are presented in Fig 2(d). We find the probability distribution (P (v g )) of v g for all the phonon modes are almost merging with that of the optical modes, asserting the connection between the significant group velocities of optical phonons and the optical phonon-dominated κ L for Ge 2 Sb 2 Te 5 . P (v g ) shows a higher probability of lower group velocities with the maximum is found to be around 500 m/s. However, it also features a long tail distribution of higher group velocities extending up to ∼ 5000 m/s. The presence of both high as well as low v g -optical phonons corresponds to the presence of significantly dispersive and comparatively flat optical bands respectively.

As both group velocity and heat capacities are temperature insensitive, probing the temperature dependent implications in thermal transport in Ge 2 Sb 2 Te 5 warrants investigation on the variation of the phonon scattering rates with temperature. The presence of lower frequency optic phonons, their strong coupling with acoustic modes as well as the presence of flat bands, observed via sharp PDOS peaks, lead us to investigate the possibility of an existence of four-phonon scattering processes. We discuss this feature in the next subsection.

B. Three and four-phonon processes: competing effects and scattering channels

We start our investigation on phonon scattering rates by comparing the relative importance between threephonon (third-order anharmonicity, denoted 3ph hereafter) and four-phonon scattering processes (fourth-order anharmonicity, denoted 4ph hereafter) in a wide temperature range (30 K -600 K). At this point, we mention two important factors that usually imprint the signature of the four-phonon processes: (a) anharmonicity and (b) scattering phase space distribution. The extent of anharmonicity can be captured via the Gruneisen parameter (γ) which is defined as the rate of change of the vibrational frequency of a particular mode with volume. Thus, γ features the departure of crystal from harmonicity. For few modes within 'M' regime, corresponding to the strong overlap between acoustic and optical bands, γ becomes anomalously large (see supplementary Fig S1(a)). However, mode-averaged Gruneisen parameter smooths out this feature at a particular temperature and therefore, the average γ at fixed T behaves in less anomalous way (Figure S1(b)). This average γ is found to be 2.49 at T = 300 K which is similar for other low-κ L materials with comparable total thermal conductivities like PbTe (γ = 2.18) and PbS (γ = 2.5) at room temperature [67].

Having recorded the strong anharmonicity, we seek to explore the relative dominance between 3ph and 4ph processes via computing their respective weighted phase space as described by Li et al. [68]. This weighted phase space is realized via the sum over frequency dependent terms in the 3ph and 4ph transition probabilities without the scattering matrix elements. Higher phase space suggests larger number of available scattering channels, leading to stronger scattering effect and consequently lower lattice thermal conductivity. Figure 3 represents the relative weight-age of 3ph and 4ph scattering phase spaces as a function of temperature ranging from 30 K to 600 K. As we increase the temperature, a clear indication of overwhelming 4ph scattering phase space compared to its 3ph counterpart is observed. However, we must mention that the weighted phase space (WP) is a qualitative descriptor and their units are different for 3ph and 4ph. Moreover, it also includes the phonon distribution function (f 0 ) while summing over all the frequency-dependent terms. This involves the effect of temperature as in the high temperature limit f 0 ∼ k B T /ℏω. However, it is evident that as the temperature is raised, 4ph phase space becomes prominent with a significant increment due to the gradual population of high frequency optical modes at higher temperature ( processes in thermal transport of Ge 2 Sb 2 Te 5 .

Figure 4 describes the frequency variation of both 3ph and 4ph scattering rates for different temperature. Phonon scattering rates have been found to vary roughly 100 and 1000 times for 3ph and 4ph scattering rates respectively within the studied temperature range (30 K -600 K). The difference between 3ph and 4ph scattering rates are seen to decrease and strong competition between these two scattering rates is realized via the gradual collapse of two scattering rates as we increase the temperature beyond T = 100 K (Fig 4(d), (e), (f)). We attribute this trend to the significant 4ph scattering phase space beyond T = 100 K, following the weighted phase space analyzes (Fig 3). The 4ph scattering rate is seen to gradually eclipse the 3ph scattering rate at high temperature, mostly in the O frequency regime (Fig 4). In the M regime, both of the scattering rates are comparable as 4ph scattering is triggered due to the coupling between acoustic and optical bands (Fig 4). In the regime O, till ω ∼ 122 cm -1 , the dominance of 4ph scattering over 3ph scattering is seen to prevail, especially at higher temperatures. This is the regime in the phonon dispersion which consists of many flat optical bands (see Appendix A for a detailed description of flat and optical bands and the corresponding allowed 3ph processes). These flat bands are shown to enable significant four phonon scattering with optical bands, consistent with an earlier study [44].

At frequencies beyond 122 cm -1 , dispersive optical bands are seen to reappear, influencing the 4ph scattering to become weaker compared to the 3ph scattering rates. To elucidate this relative behavior of 3ph and 4ph scattering rates more transparently, we study the phase space evolution of both 3ph and 4ph scattering processes, separating out the temperature effect to make a one to one correspondence between the dispersion and the phase space variation with frequency (Fig 5 ). To further understand the 3ph-minima and 4ph-maxima in their respective phase space, we deconstruct the separate roles of different scattering channels in Fig 6 . 3ph scattering consists of absorption (q + q 1 = q 2 ) and emission (q = q 1 + q 2 ) channels ( ) comprises of three scattering channels namely (a) recombination (q + q 1 + q 2 = q 3 ), (b) redistribution (q + q 1 = q 2 + q 3 ) and (c) splitting (q = q 1 + q 2 + q 3 ) along with their momentum destroying counterparts via the reciprocal lattice vector G. Here q, q 1 , q 2 , q 3 define phonon wave vectors. Figure 6(a) depicts that the acoustic phonons with low frequency modes (regime 'A') enable absorption processes which seems to gradually decrease as optical phonons overlap with the acoustic ones in the regime 'M'. In regime 'O' (60 cm -1 < ω < 174 cm -1 ), optical modes with higher frequencies are more prone to undergo emission processes which dominates at higher frequencies while absorption phase space going to zero (Fig 6(a)). The aforementioned 3phminima are found to originate as a result of the crossing between absorption and emission phase space around 90-100 cm -1 . Also, Fig 6(b) shows that 4ph redistribution phase space peaks around the same frequency where 3ph phase space shows minima. This is a result of the occurrence of flat optical bands in this frequency regime (86 cm -1 < ω < 122 cm -1 ) that severely restrict OOO and weaken AOO processes but do not prevent optical phonons to participate in four-phonon processes. Rather flat bands enhance the redistribution processes as energy conservation between four phonons is preferable in this scenario. Similar feature was observed in AlSb and BAs [36]. This causes a simultaneous dip and peak in 3ph and 4ph phase space respectively. Recombination process is seen to dominate the whole 4ph phase space spectrum of Ge 2 Sb 2 Te 5 along with a small contribution from recombination and splitting processes at low and high frequencies respectively (Fig 6 (b)). These 3ph and 4ph

phase space distributions get mildly modified with temperature as can be seen via their weighted phase space distributions (see Supplementary Fig S7). We end this section by studying the importance of phonon-isotope scattering in crystalline Ge 2 Sb 2 Te 5 compared to both 3ph and 4ph scattering rates (see Fig S2

in the supplementary information). As mentioned in an earlier section, phonon-isotope scattering rates can be computed following the work by Tamura [65], where the mode dependent isotope scattering rate varies as ω 2 but remains independent of temperature. Therefore, at low temperature (30 K < T < 60 K), phonon-isotope scattering contributes significantly to increase the resistive scattering processes and helps reducing the lattice thermal conductivity. However, as temperature is further raised (see Fig S2(d),(e),(f) in the supplementary), both 3ph and 4ph scattering start dominating the total phonon scattering processes due to their strong temperature dependencies (τ

-1 3ph ∼ T , τ -1 4ph ∼ T 2 [42]).
C. The role of normal and Umklapp scattering: three and four-phonon processes

To further investigate the nature of 3ph and 4ph scattering and their competitive roles in dictating thermal transport, we study each of the 3ph and 4ph scattering channels distinguished as separate normal (N) and Umklapp (U) scattering processes. We compare the relative weights between N and U processes by calculating the thermodynamic average of scattering rates for N and U processes separately for both 3ph and 4ph scattering via

Γ i = ⟨τ -1 i ⟩ ave = λ C λ τ -1 λi λ C λ (12)
Here, λ defines phonon modes (q, j) comprising wave vector q and branch j. Index i denotes either N or U scattering rates used. C λ is the modal heat capacity, given by

C λ = k B ℏω λ k B T 2 exp(ℏω λ /k B T ) [exp(ℏω λ /k B T ) -1] 2 (13)
where, T denotes temperature, ℏ is the reduced Planck constant and k B is the Boltzmann constant. N scattering enables momentum conserving phonon scattering, while U and I scattering bring resistivity in thermal conduction. Therefore, the relative strength of N and U scattering is pivotal to understand the underlying phonon thermal conduction as well as to capture peculiar phenomena like phonon hydrodynamics [69,70]. Although U scattering is seen dominating over N in both 3ph and 4ph scattering events, it is instructive to check their contributions to different scattering channels to understand the phonon transport mechanism. Figure 8 represents the temperature variation of thermodynamically averaged 3ph absorption and emission scattering rates (Γ 3ph ) for both N and U scattering events. N and U absorption scattering rates are seen to overlap and merge together (Fig 8). Consequently, the emission process determines the relative strength between N and U 3ph scattering. The emission processes corresponding to U scattering events are seen to dominate over the N scattering events throughout the temperature range studied. This stems from the fact that emission processes are mostly activated at high frequency regime, where optical phonons with higher momentum and large wave vectors are more prone to satisfy q = q 1 + q 2 + G, where G is the reciprocal lattice vector. Figure 9(a), (b) present the variation of average 4ph scattering rates with T for different scattering channels for N and U scattering respectively. We clearly observe that both N and U scattering rates corresponding to redistribution 4ph processes are strong and they dominate over the two other scattering channels. Moreover, the collapse of the total average scattering rates (black dashed lines D. Thermal conductivity and its variation with temperature for three and four-phonon processes:

Mode analyses

After analyzing the phonon scattering rates and their relative importance, we calculate the variation of lattice thermal conductivity with temperature incorporating both three as well as four-phonon scattering rates for Ge 2 Sb 2 Te 5 . Isotope scattering rates, which are temperature independent, are also included within the scattering terms via Matthiessen's rule [START_REF] Kaviany | Heat transfer physics[END_REF].

1 τ λ = 1 τ 3ph λ + 1 τ 4ph λ + 1 τ I λ ( 14
)
where λ defines each mode and τ 3ph λ , τ 4ph λ and τ I λ denote 3ph, 4ph and phonon-isotope scattering lifetimes respectively. Figure 10( Figure 10(b) records a value of 1.03 W/mK at 300 K and 0.52 W/mK at 600 K, for the average κ L , calculated using only three-phonon and isotope scattering rates (red solid line for iterative and blue dashed line for RTA), consistent with earlier works by Ibarra et al. [15]. Campi et al.

[29] and Mukhopadhyay et al. [28] obtained a value of 1.2 W/mK and 1.4 W/mK respectively, considering an ideal crystal with the inclusion of only 3ph scattering for the Kooi structure at 300 K. Further, the significant implication of incorporating four-phonon scattering in κ L can also be seen in Fig 10(b). By including four-phonon processes with three-phonon and isotope scattering events, a significant drop in κ L is observed at high temperature. At 100 K, we record a reduction of ∼ 13 % (κ L = 2.69 W/mK), followed by a ∼ 28 % (κ L = 0.74 W/mK) and a striking ∼ 42 % (κ L = 0.3 W/mK) decrement in κ L at 300 K and 600 K respectively.

Further, we study the alteration of relative contributions coming from acoustic (A), mixed (M) and optical (O) frequency regimes with and without 4ph scattering processes. Figure 11 represents the cumulative lattice . At T = 30 K, acoustic phonons seem to contribute mostly (47.3 %) to κ L followed by the M regime (42.6 %). As the temperature is increased this trend seems to alter with optical phonons dominating κ L (40 %). We note that the M regime also contributes significantly (∼ 33 %) to κ L . At temperatures beyond 300 K, this trend continues with the contributions from different frequency regimes to the total κ L are found to be saturated (regime A: 26.4 %, regime M: 33.2 %, regime O: 40.4 %). To explore the modification of this picture after adding 4ph scattering, modified κ c L is also presented (red dashed lines in Fig 11(a)-(f)). Initially, at low temperature (T = 30 K), we see a similar trend of acoustic phonon dominated phonon transport for Ge 2 Sb 2 Te 5 (regime A: 46.7 %, regime M: 43.3 %, regime O: 10 %). As we increase the temperature, optical modes seem to dominate κ L , similar to that of the cubic anharmonicity. However, the contributions from M and O regimes are found to be slightly affected after the inclusion of 4ph scattering events. When temperature is raised beyond 300 K, the contributions from A, M and O frequency regimes are found to be saturated around 26.2 %, 35.8 % and 38 % respectively (see red dashed lines in Fig 11(d), (e), (f)). We observe that at high temperature, including 4ph scattering does not alter the contributions from the regime A. However, the relative contributions of M and O regimes are marginally interchanged. These results indicate that including 4ph scattering processes does not alter the optical mode-dominance in κ L for Ge 2 Sb 2 Te 5 . Rather, including 4ph scattering affects both acoustic and optic modes non-selectively, retaining the significant contribution of optical modes in κ L of Ge 2 Sb 2 Te 5 . This is in sharp contrast to high-κ L material like AlSb [36], where optical phonon-dominated κ L gets almost washed out while incorporating 4ph scattering processes.

Next, to investigate the size dependence of the thermal conductivity of Ge 2 Sb 2 Te 5 and its implications due to the presence of 4ph scattering, we present the variation of the maximum mean free path (λ max ) as a function of temperature (Fig 12). We note that grain boundaries, their scattering with phonons as well as the thickness were found to influence κ L for Ge 2 Sb 2 Te 5 [46,47]. Therefore, studying size-dependence of thermal conduction serves not only the purpose of understanding the MFP distribution and their contributions to κ L , but also of designing nanostructures by assessing the maximum mean free path inside the material (see supplementary path along a-b plane (κ x L ) and along c (κ z L ), with and without 4ph scattering). Figure 12 shows that λ max gradually decreases with temperature. Over the temperature range studied here, λ max is seen to vary from 3.1 µm to 75.3 nm for 3ph+I and from 2.1 µm to 24.7 nm for 3ph+4ph+I. Further, λ max is found to vary roughly as T -1 (λ max ∼ T -1.1 ) while including 3ph and I scattering (red circles and their fitted dashed line in the inset of Fig 12). This trend emerges due to the fact that higher temperature populates more phonons according to phonon distribution (BE distribution), inducing more scattering between phonons to reduce the mean free path between them. However, including 4ph scattering is seen to modify this negative slope with even a higher value of the exponent (λ max ∼ T -1.4 , blue diamonds and their fitted dashed line in the inset of Fig 12). We recall that mean free paths and phonon group velocities are connected as λ = vτ . Therefore, this is consistent with a higher order scattering scenario where the scattering rate varies as τ -1 ∼ T α where α > 1.

To explain this size effect in a more comprehensive way, we present the percentage of cumulative contribution to κ L as a function of maximum allowed MFPs for different temperature in ). This finding relates to a recent discussion on the room temperature experimental thermal conductivity of both cubic and hexagonal Ge 2 Sb 2 Te 5 [25]. Time-domain thermoreflectance measurements with grains of 10-15 nm were found to influence negligibly the room-temperature κ L of hexagonal Ge 2 Sb 2 Te 5 [25]. This behavior can be rationalized using the fact that almost 80 % of the κ L of hexagonal Ge 2 Sb 2 Te 5 are seen to come from phonons with mean free paths below ∼ 17 nm (20 nm if considering only 3ph+I) as shown in Fig 13. This suggests that the experimental κ L (0.47 W/mK) at 300 K has negligible effects due to phonon-grain boundary scattering. To realize serious effects of grain boundary scattering in hexagonal Ge 2 Sb 2 Te 5 , grains need to be smaller to enhance grain-boundary scattering (scattering rate τ -1 B = v g /L). For example, grains of size 3 nm can effectively reduce 60 % of the κ L at room temperature as can be seen via Fig 13(b). Here, we note that the obtained κ L at high T (≥ 600 K) is sensitive to the validity of phonon BTE as we include the 4ph scattering. Figure 13 shows growing contributions of phonons with MFPs comparable to the lattice parameter (a) from 3.4 % to 12.4 %, while including 4ph scattering (blue shaded regimes in Fig 13). Though λ max (24.7 nm) is well above a, non-negligible contributions from MFPs (≤ a) hint at the closeness of κ L to the amorphous limit (κ min ∼ 0.27 W/mK [9]) where the usage of phonon BTE seems questionable.

The negligible effect of grain-boundary scattering as well as the significant presence of four-phonon scattering trigger the discussion concerning the gap between experiments and theory to estimate the thermal conductivity of Ge 2 Sb 2 Te 5 . Experimental lattice thermal conductivity had been found to be around 0.42-0.48 W/mK [20,25,26,71], which is very close to the minimal thermal conductivity using the Cahill model [72], considering disorder-dominated structure of Ge 2 Sb 2 Te 5 . This raises concern about the Boltzmann transport equation approach to compute the thermal conductivity as was mentioned earlier [29]. However, our findings suggest that phonon-phonon scattering has a strong influence with the presence of quartic anharmonicity that can bring down κ L to 28 % at 300 K and 42 % at 600 K. More importantly, 4ph and disorder scattering (both phononvacancy [73] and phonon-isotope scattering [65]) are two completely different mechanisms following T 2 and Tindependent behavior respectively. Therefore, excluding quartic anharmonicity at room temperature and beyond, can lead to a crucial overestimation of the predicted disorder content at the cost of an underestimation of phonon-phonon contribution to κ L of Ge 2 Sb 2 Te 5 . Further, experiments [74] had also shown to produce fairly ordered structure for Ge 2 Sb 2 Te 5 by weakening the randomness by ordering the vacancies. Therefore, higher order anharmonicity in Ge 2 Sb 2 Te 5 is crucial to estimate thermal behavior of these highly ordered materials.

E. Anisotropy under three and four-phonon scattering

The anisotropy in the thermal conductivity for GST materials can propose efficient designing of PCM devices that hinders the cross talks between neighboring PCM cells [47]. Indeed, significant differences between x and z component of κ L are recorded while assessing its variation with phonon frequency and mean free path (see

supplementary Fig S11 -Fig S14)
. Figure 14 presents the temperature variation of κ L -anisotropy by using the ratio of κ L along a-axis to that of the c-axis. We see a marginal increment in the anisotropy factor (κ xx L /κ zz L ) from 2.04 to 2.1 with temperature, employing RTA and iterative method respectively within the three-phonon picture. We notice that the factor remains almost constant (∼ 2) at temperature beyond 300 K, consistent with the work of Mukhopadhyay et al. [28]. However, various experimental and theoretical studies, having consistent κ L , feature a wide variation of the anisotropy factor from 1.9 to 4.6 [17,28,29,47,48] for Ge 2 Sb 2 Te 5 at 300 K. Further, including 4ph scattering is found to increase the anisotropy factor from 2.04 to 2.21 (Fig 14) which is a marginal (∼ 8 %) increment. This is an indication of the negligible influence of scattering rates in the κ L -anisotropy. We evaluate the changes in anisotropy as we include 4ph scattering in terms of frequency evolution of the cumulative κ L . Figure 15 describes the variation of the anisotropic factor (κ c L (xx)/κ c L (zz)) as a function of frequencies both including (black solid lines) and excluding 4ph scattering events (blue dashed lines). We quantify the anisotropic behavior by the increment of the factor beyond 1 (shown via dashed green lines in Fig 15). A marginally higher value of anisotropy factor is seen to evolve with frequency while including 4ph scattering to the existing 3ph and isotope scattering processes. The increment of this factor towards a final value comprises of contributions from all three frequency regimes. At low temperature of 30 K (Fig 15(a)), A, M and O regimes contribute to ∼ 37 %, ∼ 55 % and ∼ 8 % respectively to the total anisotropy factor, while including 4ph scattering makes these contributions to ∼ 49 %, ∼ 46 % and ∼ 5 % respectively. At T = 300 K (Fig 15(d)), these contributions are found to modify as ∼ 24 %, ∼ 48 % and ∼ 28 % respectively for 3ph+I scattering and ∼ 39 %, ∼ 39 % and ∼ 22 % respectively for 3ph+4ph+I scattering events. At even higher temperature (T = 600 K, Fig 15(f)), with and without four-phonon scattering, the anisotropy factor is found to be distributed as follows: (i) excluding four-phonon scattering: ∼ 25 %, ∼ 48 % and FIG. 16. The variation of inverse lattice thermal conductivity (1/κL) with temperature for crystalline Ge2Sb2Te5 considering only three-phonon and isotope scattering using iterative (red circles) and RTA method (green hexagons). T -linear high temperature regime is obtained with a slope a1 = 0.0032 and a2 = 0.0034 for iterative and RTA methods respectively. The coefficient from the Slack model is also shown (0.026).

∼ 27 % contributions are coming from A, M and O regimes respectively. (ii) Including four-phonon scattering: ∼ 37 %, ∼ 40 % and ∼ 23 % contributions are coming from A, M and O regimes respectively. Irrespective of the included scattering processes, the contribution from optical phonons (regime O) is seen to significantly increase at the cost of a reduced contribution from the A regime. The variation of the relative contributions from different frequency regimes is found to be weakly-dependent on temperature as κ L -anisotropy in Ge 2 Sb 2 Te 5 is predominantly dictated by the anisotropy in phonon group velocity. Mukhopadhyay et al. [28] used a density regression method using PDOS to capture the anisotropy in phonon group velocity along the a and c axes. We observe a consistent behavior with a little steeper acoustic and optical branches along Γ-M or Γ-K compared to that of the Γ-A in the dispersion relation for Ge Throughout this study we unfold two significant features of Kooi-Ge 2 Sb 2 Te 5 : (a) significant domination of optical phonon in thermal transport and (b) considerable effect of quartic anharmonicity especially at high temperature. Crucial effect of these factors on the thermal transport leads us to introspect their further implications in the scaling properties of thermal conductiv-ity with temperature for crystalline Ge 2 Sb 2 Te 5 . At high temperature, in the kinetic thermal transport regime, κ L follows a T -1 scaling [START_REF] Kaviany | Heat transfer physics[END_REF] as a result of dominating Umklapp scattering between phonons. This high temperature behavior in nonmetallic materials had been analytically derived by Slack and Morelli [32,[START_REF] Morelli | High thermal conductivity materials[END_REF][START_REF] Morelli | [END_REF]. This model was quite successful in describing the high temperature T -1 behavior of κ L for a wide range of systems for which acoustic phonons in general and Umklapp scattering in particular, play a dominant role to κ L [32]. Therefore, often T -1 behavior of nonmetallic systems at high temperature range is attributed to this model. Considering only acoustic phonons as heat transport carriers, Slack and Morelli [START_REF] Morelli | [END_REF] developed the intrinsic κ L of a solid in a temperature range where Umklapp scattering is dominant as

κ L = A M θ 3 D δ γ 2 n 2/3 T ( 15
)
with

A = 2.43 × 10 -6 1 -0.514/γ + 0.228/γ 2 ( 16
)
where n is the number of atoms per unit cell, δ 3 is volume per atom (δ is in the unit of Å), M is the average atomic mass of the unit cell of the solid, γ is Gruneisen parameter and θ D is the Debye temperature of the material. Firstly, we analyze the high temperature T -1 behavior of our first-principles result of κ L considering only 3ph and isotope scattering for Ge 2 Sb 2 Te 5 . We recall the optical phonon dominance in κ L at high temperature, revealed by our ab initio studies. However, simultaneously we have also found the dominance of Umklapp scattering at high temperature regime. This abundance of Umklapp resistance is seen to manifest itself via a T -1 scaling of κ L at high temperature as shown in Figure 16. The linear fit between 1/κ L and T at high temperature regime yields the coefficients a 1 (= 0.0032) or a 2 (= 0.0034) considering iterative and RTA method respectively (Figure 16). To compare these coefficients to that of the Slack-Morelli model, we compute the coefficient from Eq.15 using different parameters presented in Table II. However, the Slack-Morelli coefficient yields a value 0.026, which is higher than the fitted ab initio results by an order of magnitude. Although the abundance of Umklapp scattering dictates T -1 scaling of κ L for both cases, the mismatch between the coefficients can be attributed to the fact that Slack model is analytically derived for acoustic branches while our results show optical phonon dominated thermal transport in Kooi-Ge 2 Sb 2 Te 5 .

The dominance of optical modes in thermal transport was observed to feature a deviation of κ L from its T -1 scaling at higher temperature range in monolayer GaN [77]. Although, κ L of Ge 2 Sb 2 Te 5 is seen to possess T -1 scaling, including four-phonon can radically change the scaling behavior. As both specific heat and group velocity are temperature independent, the temperature scaling of κ L at high temperature relies on the temperature dependencies of τ -1 3ph and τ -1 4ph . Figure 17(a) presents the temperature variation of the total three-phonon and four-phonon scattering rates of all modes. At high temperature regime, an almost T -linear scaling is observed for τ -1 3ph (∼ T 0.97 ) while τ -1 4ph retrieves a ∼ T 2.16 scaling. Further, we present 3ph and 4ph scattering rates for each mode as a function of temperature to more clearly visualize the T and T 2 variation of τ -1 3ph (Figure 17(b)) and τ -1 4ph (Figure 17(c)) respectively, featuring a generic trend observed in other materials [42]. Figure 17(d) shows that regardless of acoustic or optical phonons, specific scaling relations are followed by 3ph and 4ph scattering processes. 3ph and 4ph scattering rates for the first six phonon bands are seen to feature linear (dashed lines in Fig 17(d)) and quadratic (solid lines in Fig 17(d)) temperature variation respectively, regardless of acoustic (blue lines in Fig 17(d)) or optical modes (red lines in Fig 17(d)). The nonlinear temperature scaling of the four-phonon scattering rate indicates an altering temperature dependence of the lattice thermal conductivity in crystalline Ge 2 Sb 2 Te 5 . Figure 18(a) shows that the inclusion of 4ph scattering processes leads to a deviation of the typical T -1 scaling of κ L at high temperature (T > θ D ). The variation of 1/κ L with T including 4ph scattering yields an exponent 1.3 (black dashed line in Fig 18(a)). As the temperature variation of κ L is directly dependent on the T dependence of the phonon scattering rates, the departure from κ L ∼ 1/T can be understood as κ L ∝ 1/τ -1 ∝ 1/(τ -1 3ph + τ -1 4ph ) ∼ 1/(aT +bT 2 ). Indeed, Figure 18(a) shows that 1/κ L varies as aT +bT 2 , where a and b are constants. The departure of κ L from typical 1/T dependence for Ge 2 Sb 2 Te 5 at high T evokes complexities to the physical understanding of the universality of the lower bound to the thermal diffusivity of materials, as discussed recently [78]. Behnia and Kapitulnik [78] extended the idea of a universal lower bound to thermal diffusivity in electronic transport in metals [79,80] to the non-metallic phonon transport in materials. At high temperature, phononic thermal diffusivity D was found to exhibit a lower bound, governed by the speed of sound and a Planckian scattering time (τ p ) via [78]

D = sv 2 s τ p ( 17
)
where, thermal diffusion coefficient or diffusivity D = κ L /C, C denotes heat capacity, τ p = ℏ k B T , v s is average sound speed obtained from the longitudinal acoustic (LA) mode and s (> 1) denotes a dimensionless parameter that is constant for a specific material. This bound had also been realized to have a quantum-mechanical origin which needs further investigations [78]. We find a consistent D -1 ∝ T behavior at higher temperature . At temperature higher than θ D (=100 K for Ge 2 Sb 2 Te 5 ), heat capacity is almost constant. Therefore, the lower bound of D, featuring a relaxation time which is inverse of temperature, stems from the T -1 decay of κ L in the kinetic regime of phonon transport. This suggests a physical picture of phonons performing random walk to create a completely diffusive dynamics under dominating Umklapp scattering in this kinetic regime. However, this theory faces a serious challenge while in-corporating quartic anharmonicity for Ge 2 Sb 2 Te 5 . The deviation from T -1 behavior of κ L for Ge 2 Sb 2 Te 5 at high temperature severely limits to obtain a universal lower bound of thermal diffusivity using Eq.17. Within the domain of electronic transport, electron-doped cuprates were found to violate this Planckian bound owing to their high temperature T 2 dependence of resistivity in lieu of a typical T -linear resistivity [81]. This T 2 variation of resistivity helps surpassing the Planckian bound featuring a faster decay with τ -1 ≫ τ -1 p , a representative of super-Planckian behavior [80] in electronic transport. Borrowing the idea of 'super-Planckian' behavior to our phononbased description, we connect the departure of thermal diffusivity from T -1 scaling and the violation of Planckian bound to the combined effect of T -linear (τ -1 3ph ∼ T ) and T 2 (τ -1 4ph ∼ T 2 ) variation of phonon resistance at high temperature thermal transport in crystalline Ge 2 Sb 2 Te 5 . Figure 18(b) (red circles) shows the deviation of D -1 from the T -linear behavior while incorporating 4ph scattering processes. This deviation was rationalized for insulators via a velocity bound v M ≳ v s [82], that leads to super-Planckian scattering. Here v M is melting velocity defined as v M = k B T M ℏ a, T M is crystal melting temperature, and a denotes lattice spacing. Using T M = 900 K, v s = 3187.65 m/s, and lattice parameter from Table I, we get v M /v s ≈ 15.5 for Ge 2 Sb 2 Te 5 , confirming a super-Planckian type phononic resistive behavior. Further, to capture this deviation as well as to extend Eq.17 compatible for phonon scattering with quartic anharmonicity, we use the form

τ ′ p = ℏ k B (T + bT 2 ) ( 18
)
where b is a constant fitted parameter. Employing Eq.18 in D = sv 2 s τ ′ p and using this expression to fit the thermal diffusivity at high temperature of Ge 2 Sb 2 Te 5 yields a value of s = 2.76. We recall that s was originally defined [78] as a material specific constant. Consequently, thermal diffusivity with and without the inclusion of fourphonon scattering gives consistent s (2.76 and 2.79 respectively) for Ge 2 Sb 2 Te 5 (Fig 18(b)). Hence, thermal diffusivity of crystalline Ge 2 Sb 2 Te 5 with its significant higher order phonon anharmonicity seems to exhibit dissipation of super-Planckian type and the corresponding time scale needs to be modified to capture its high temperature behavior. This feature definitely demands more introspection in the future.

IV. SUMMARY AND CONCLUSIONS

First-principles theoretical investigation has been systematically carried out to unveil the phonon thermal transport in crystalline Ge 2 Sb 2 Te 5 , a supremely useful PCM and thermoelectric material possessing low lattice thermal conductivity. The variation of harmonic and anharmonic thermal transport properties of Ge 2 Sb 2 Te 5 with phonon frequencies is realized in a wide range of temperature from 30 K to 600 K. Using the Kooi structure of Ge 2 Sb 2 Te 5 , the study reveals two notable thermal transport features: (a) significant optical modedominance of phonon transport and (b) crucial contributions of quartic anharmonicity to the lattice thermal conductivity, especially at high temperature. Phonon dispersion relations are found to exhibit three broadly distinguished frequency regimes comprised of acoustic bands (A), mixed acoustic and optical bands (M) and optical bands (O). Low-frequency optical bands are found to lie close to the acoustic branches in the 'M' regime, enhancing the scattering of acoustic phonons by the optical ones to decrease κ L of Ge 2 Sb 2 Te 5 . Phase space analysis reveals that four-phonon (4ph) scattering becomes gradually important at higher temperature compared to that of the three-phonon (3ph) scattering phase space. Moreover, 4ph phase space is found to display a maxima simultaneously with a minima in its 3ph counterpart in a frequency domain 86 cm -1 < ω < 122 cm -1 within the 'O' regime. This regime hosts flat bands that prevent some 3ph scattering (especially OOO type) but enable 4ph redistribution processes. Moreover, this 3ph minima results from a crossover between absorption and emission phase space. Absorption dominates at low frequencies whereas emission becomes significant at higher frequencies. Amongst the 4ph scattering channels, redistribution dominates throughout the whole frequency regime. Both 3ph and 4ph scattering processes are found to be predominantly of Umklapp type in the temperature regime studied here. Eventually, this assists in imparting a strong resistance to phonon transport in Ge 2 Sb 2 Te 5 . The closeness between iterative and RTA solutions of BTE for κ L reinstates this feature. Inevitably, including 4ph processes with 3ph and isotope scattering events are found to significantly reduce κ L of Ge 2 Sb 2 Te 5 , starting from a drop of ∼ 13 % at 100 K (κ L = 2.69 W/mK), followed by ∼ 28 % (κ L = 0.74 W/mK) and a significant ∼ 42 % decrement (κ L = 0.3 W/mK) at 300K and 600 K respectively. The domination of optical phonons in κ L is retained in the presence of 4ph scattering with 'M' and 'O' regimes contributing ∼ 36 % and ∼ 38 % respectively to κ L , at high temperature beyond 300 K. Including 4ph scattering has also been found to systematically decrease the maximum allowed mean free paths with temperature with a different scaling (∼ T -1.4 ) compared to a typical ∼ T -1 relation. Almost 80 % of the κ L of hexagonal Ge 2 Sb 2 Te 5 are borne out of phonons with mean free paths below ∼ 17 nm. Recent experiments [25] used similar grain size to experimentally obtain κ L for Ge 2 Sb 2 Te 5 which suggests that at room temperature, grain boundary scattering has negligible contribution to κ L . Our investigation suggests that the damping of κ L in Ge 2 Sb 2 Te 5 is strongly influenced by four-phonon processes, where a significant drop in κ L is found to occur at 300 K and beyond. Therefore, the earlier theoretical estimates on the κ L of Ge 2 Sb 2 Te 5 , which did not include 4ph scattering but used different disorder concentration to match the experimental κ L , need to be scrutinized. This modification is crucial since explaining the experimental results with theoretical predictions without 4ph scattering in κ L , can lead to misleading physical understanding. 4ph scattering is also seen to marginally enhance the anisotropy in thermal transport of Ge 2 Sb 2 Te 5 compared to the 3ph only framework. The anisotropy has a very weak temperature dependence as it is mostly dictated by the anisotropy in phonon group velocities along a -b and c axes. However, we must mention some limitations of this study. Firstly, we note here the importance and differences of the theoretical approaches to analyze the Kooi and the Matsunaga [60] structure of Ge 2 Sb 2 Te 5 . While theoretical studies on Kooi or Petrov structure assume a pristine initial configuration with disorder included perturbatively, Matsunaga structure [29] is stabilized with the presence of Sb/Ge disorder within the initial configuration. Therefore, this disorder-rich Ge 2 Sb 2 Te 5 was seen to possess lower v g leading to lower κ L compared to that of the Kooi structure [29]. This was also seen to affect the optical phonon contributions to κ L which goes negligibly small. Four-phonon processes may be of negligible importance in that case. This meager effect can also possibly mask the faster decay of κ L compared to T -1 scaling at high temperature. However, ordering vacancies [74] or producing disorder-starved crystals via recent advanced experimental methods still demands a better understanding of phonon dynamics realized via Kooi and Petrov structures of Ge 2 Sb 2 Te 5 . Secondly, with current stateof-the-art, four-phonon processes could only be realized in RTA to reduce a serious overburden of computational cost. This feature, at least for Ge 2 Sb 2 Te 5 , seems to be less critical as at high temperature regime four-phonon Umklapp scattering is found to dominate over normal scattering, which suffices the usage of RTA at the fourphonon level. Inclusion of four-phonon processes has further serious implications on the temperature-scaling properties of κ L in the high temperature regime. Optical mode-dominated κ L , observed in our study, is found to show its signature through the disagreement with the Slack-Morelli description, captured via the coefficient of the T -1 scaling of κ L . Including 4ph scattering, a departure of κ L from the conventional T -1 scaling is observed and κ L is found to vary as 1/(aT +bT 2 ). This is rationalized via the T dependencies of the combined effect of 3ph and 4ph scattering. This deviation from T -1 scaling further questions a recently conceived universal lower bound of thermal diffusivity in phononic materials. While 3ph-dominated κ L maintains this bound, we recover a super-Planckian type phononic resistance. Consequently, a modified time scale that sets the lower bound in thermal diffusivity in crystalline Ge 2 Sb 2 Te 5 is obtained. These results will help predicting lattice thermal conductivity of Ge 2 Sb 2 Te 5 in a more precise fashion as well as can offer better predictions to disorder-driven phonon scattering which can be controlled and tuned for thermoelectric and PCM appli-cations. Also, further studies in the future should dedicate in elucidating the high temperature lower bound of thermal diffusivity in phononic materials having significant four-phonon scattering. Even more introspection is needed to fully decipher its quantum mechanical underpinnings.

1. AAO: These 3ph processes can occur if there is an involvement of at least one acoustic phonon whose frequency is higher than the acoustic-optical band gap. As dispersion relation (Fig 1(C)) suggests that there is no apparent A-O band gap, rather an overlap between A and O phonons, this criterion is satisfied. However, there is also a cut off frequency [83] for O phonons (2∆ω A ), beyond which O phonons can not engage in AAO processes. We note that three acoustic branches have bandwidths (∆ω A ) of 44.4 cm -1 (TA1), 49.3 cm -1 (TA2), and 59.1 cm -1 (LA). Therefore bands 21-24 (122 cm -1 < ω < 174 cm -1 ) only have weak AAO scattering due to the fact that only LA branch can participate.

2. AOO: It is possible if ω A ≤ ∆ω O . Therefore, few optical bands will always allow this process as shown in Column 3 of Table III. Obviously, due to extremely small ∆ω O , flat bands will allow lesser AOO processes than the highly dispersive bands.

3. OOO: These 3ph processes will activate if ∆ω O ≥ ω min O . Therefore, these processes will be absent in flat bands (bands from 8-20) and will be present marginally for dispersive bands at high frequencies (bands from [21][22][23][24]. III. Distinguishing optical phonon bands in Ge2Sb2Te5 as dispersive and flat bands and their corresponding selection rules due to energy conservation [83] to activate different three-phonon (3ph) processes. ωA and ωO denote acoustic and optical mode frequencies respectively. Here ω min O ∼ 30 cm -1 , denotes lowest available optical phonon frequency. The three acoustic branches have bandwidths (∆ωA) of 44.4 cm -1 (TA1), 49.3 cm -1 (TA2), and 59.1 cm -1 (LA). Optical band ids and band widths corresponding to flat bands are designated via blue in Table III We observe very few overdamped [84,85] modes in the 'M' frequency regime of Ge 2 Sb 2 Te 5 via the variation of 3ph scattering rates with frequency (Fig 4). This feature was also present in an earlier first-principles study of hexagonal Ge 2 Sb 2 Te 5 [28] although not explained. At high temperature beyond 100 K, very few sporadic overdamped phonon modes (ω 0 τ 3ph < 1) are seen in the 'M' regime, while at low temperature the modes become underdamped (ω 0 τ 3ph > 1). We link this feature to the strong anharmonicity in hexagonal Ge Although, the sporadic overdamped modes seem to pose a question on the quasiparticle picture, its validity also depends on the order of the phonon mean free paths.

As was mentioned by Allen [84], quasiparticle theory requires larger mean free path compared to the lattice constant. The variation of the maximum mean free path (λ max ) as a function of temperature is shown in Fig 12 . λ max is seen to vary from 3 µm to 75 nm for 3ph processes which is several order larger compared to the lattice parameters of Ge 2 Sb 2 Te 5 . Moreover, within the three-phonon framework, the variation of λ max shows 1/T trend (λ max ∼ T -1.1 , inset of Fig 12). This comes from the fact that higher temperature populates more phonons according to phonon distribution (BE distribution), inducing more scattering between phonons to reduce the mean free path. Thus the variation of λ max with T also indicates that phonon quasiparticle picture is valid. We note that Campi et al. [29] also mentioned about this validity of quasiparticle picture and the usage phonon BTE in the context of Ge 2 Sb 2 Te 5 . However, the validity of perturbation techniques beyond harmonic approximation needs more attention in systems with overdamped modes [84] and it definitely demands more understanding. Ge 2 Sb 2 Te 5 could be a good candidate to test this issue in future.

Appendix C: Supplementary material 1. Gruneisen parameter of Kooi-Ge2Sb2Te5 λI , beyond which 3ph and 4ph scattering dominate the scattering processes due to their strong temperature dependencies in contrast to the T -independent phononisotope scattering.

The temperature variation of the specific heat in

Kooi-Ge2Sb2Te5

The Debye temperature (θ D ) of the material is 100 K which is noted in Table II in However, unlike Kooi structure, Petrov structure is seen to feature a frequency gap between the optical bands around 90 cm -1 . This feature is reflected in the PDOS (Fig S4(a)) via a drop of PDOS to almost zero around 90 cm -1 , indicating almost absence of phonon modes at this frequency. This quasibipartite phonon density of states had also been discussed in earlier studies [15,29] with the same onset frequency of 90 cm -1 [15] in Petrov-Ge 2 Sb 2 Te 5 . Figure S6 represents the relative strength between 3ph and 4ph weighted scattering phase spaces in Petrov-Ge 2 Sb 2 Te 5 as a function of temperature ranging from 30 K to 600 K. As the temperature is raised, we find a gradually increased relevance of 4ph scattering phase space compared to that of the 3ph phase space. This stems from the fact that the evaluation of weighted phase space takes the phonon distribution function (f 0 ) into account which varies linearly with temperature at high temperature limit (f 0 ∼ k B T /ℏω). This indicates that the thermal transport of Petrov-Ge 2 Sb 2 Te 5 can also be influenced by incorporating the 4ph scattering processes.

6. Three and four-phonon weighted phase space for Kooi-Ge2Sb2Te5 decreasing 4ph phase space indicates restriction to 4ph scattering channels due to fewer sparsely located highly dispersive optical bands. Also, the crossover between the redistribution-dominated and splitting-dominated scattering processes is marginally postponed to higher frequencies at higher temperatures (Fig S7 (d)). For both three and four-phonon processes in Ge 2 Sb 2 Te 5 , Umklapp scattering is seen to dominate over normal scattering processes as realized through the temperaturevariation of their relative strength via the ratio Γ U /Γ N . Γ is average scattering rate defined as Γ = λ C λ τ -1 λ / λ C λ and is calculated for both Umklapp (Γ U ) and normal scattering (Γ N ) events. In the studied temperature range, this ratio is almost constant with a value 1.3 for 3ph scattering which increases to 2.5 for 4ph scattering. . At high T (600 K), these few sporadic N absorption modes feature overdamped phonon dynamics (ωτ < 1) as discussed in Appendix B. We also observe that τ -1 3ph is always larger for emission processes than that of the absorption processes regardless the temperature. Further a more quantitative picture of the relative importance between N and U scattering is obtained via the thermodynamic average of scattering rates as described in the article, which suggests U scattering dominance in 3ph scattering. Fig S10 presents the frequency distribution of fourphonon scattering rates in different scattering channels, distinguished as N and U processes. Redistribution processes in general, resistive U redistribution to be particular (q + q 1 = q 2 + q 3 + G), are seen to overpower at both high and low temperature. At low temperature (Fig S10(a), (c), (e)), the M frequency regime is seen to be populated by mostly redistribution as well as recombination processes that exhibit U scattering. In O frequency regime, recombination channels are mostly blocked while splitting processes are found to gain importance along with the dominating redistribution scattering rates. At very high frequencies splitting processes aggravate due to the increasing weighted scattering phase space. This frequency distribution remains almost similar at higher temperature (Fig S10(b), (d), (f)). We realize that both 3ph and 4ph scattering processes are predominantly of Umklapp type in Ge 2 Sb 2 Te 5 which is more evident from the temperature variation of the average scattering rate, discussed in the main text. We use the q-grid 16×16×4 (mentioned in the manuscript in the Methods section) for calculating the lattice thermal conductivity both including and excluding 4ph with the existing 3ph processes, to compare the effect of three and four-phonon scattering consistently. Figure S15(a) shows the variation of average lattice thermal conductivity with q-grid at T = 300 K using 3ph and isotope scattering. Beyond q-grid 16×16×4, thermal conductivity is seen to increase marginally (∼ 3.8% and ∼ 4% differences between 16×16×4 and 24×24×4 for iterative and RTA method respectively), suggesting a convergence of κ L in three-phonon framework. For convergence of 4ph processes, we present the variation of total four-phonon phase space volume with q-grid in Fig S15 (b). A tiny difference of 1.2 % between the calculated total 4ph phase space with q-grid 12×12×4 and 16×16×4 ensures convergence in 4ph processes. We must mention that computation costs involving 4ph scattering are very high as Ge 2 Sb 2 Te 5 consists of 9 atoms in the unit cell. Therefore, 16×16×4 q-grid optimally fit for this investigation.

  structure. Both Kooi and Petrov stacking structures for hexagonal Ge 2 Sb 2 Te 5 are shown in Fig 1(A)(a) and Fig 1(A)(b) respectively, using two formula units stacked along c and the periodic copies of atoms at the edges in the ab planes. A fairly good agreement on the lattice parameters between earlier first-principles calculations and this work is observed from Table

2 2

 22 Sb Te 5 Kooi structure, depicted in Fig 1(C) and Fig 1(B) respectively. Phonon dispersion relation is shown along the Γ-M -K-Γ-A high symmetry direction in the Brillouin

FIG. 1 .

 1 FIG. 1. (A) The schematic hexagonal structures of Ge2Sb2Te5 in two different stacking arrangements (a) Kooi and (b) Petrov, are presented. C axes are designed to accommodate two formula units of the nine atom hexagonal unit cell. The corresponding nine atom unit cell stacking orders for Ge, Sb, and Te atoms in Kooi and Petrov structures have also been shown in (a) and (b) respectively. The representations are realized through VESTA [63]. (B) Phonon density of states (PDOS) is presented as a function of frequency for crystalline hexagonal Ge2Sb2Te5 stacked in Kooi structure. Corresponding partial DOS for Ge, Sb and Te are also superimposed. (C) Phonon dispersion relations of Ge2Sb2Te5 Kooi structure along high symmetry path Γ-M -K-Γ-A. The frequency regimes corresponding to acoustic (A), mixed (M) and optical (O) modes are also distinguished.

FIG. 2 .

 2 FIG. 2. (a) The variation of average lattice thermal conductivity (κL) with temperature realized via both iterative and relaxation time approximation (RTA). (b) The variation of average lattice thermal conductivity, obtained via iterative method, as a function of temperature with distinguished contributions from acoustic (A), mixed (M) and optical frequency regimes (O) as described in the text. The regimes with different colors mark two separate crossings between acousticoptic contributions and optic-mixed contributions. (c) The variation of phonon group velocity (vg) for both acoustic and optical modes with frequency. Acoustic, mixed and optical regimes are denoted via region I, II, III respectively. (d) Probability distribution (P (vg)) of phonon group velocities (vg) for optical as well as all modes.

FIG. 3 .

 3 FIG. 3. The variation of three-phonon (3ph) and four-phonon (4ph) weighted phase space with phonon frequencies at different temperatures: (a) T = 30 K, (b) T = 60 K, (c) T = 100 K, (d) T = 300 K, (e) T = 500 K, and (f) T = 600 K for Ge2Sb2Te5.

  FIG. 4. The variation of three-phonon (3ph) and four-phonon (4ph) scattering rates with phonon frequencies at different temperatures: (a) T = 30 K, (b) 60 K, (c) 100 K, (d) 300 K, (e) 500 K, and (f) 600 K. Different phonon frequency regimes are marked using white (A), grey (M) and light red (O).

FIG. 5 .

 5 FIG. 5. The variation of three-phonon (3ph) and four-phonon (4ph) scattering phase spaces with phonon frequencies. Different phonon frequency regimes are marked using white (A) and grey (M) and light red (O). Flat bands within regime 'O' are distinguished via the dashed regime with 86 cm -1 < ω < 122 cm -1 .

Figure 5

 5 FIG. 6. The variation of (a) 3ph and (b) 4ph phase spaces with phonon frequencies for different scattering channels for Ge2Sb2Te5.beyond 122 cm -1 , dispersive optical bands reappear which rejuvenate the 3ph processes. As a result 3ph phase space is seen to dominate again over 4ph. The Petrov structure of Ge 2 Sb 2 Te 5 is also found to possess a dominating 4ph weighted phase space over its 3ph counterpart although in a less pronounced fashion compared to that of the Kooi structure (see Supplementary Fig S5 and Fig S6 for details).To further understand the 3ph-minima and 4ph-maxima in their respective phase space, we deconstruct the separate roles of different scattering channels in Fig 6.3ph scattering consists of absorption (q + q 1 = q 2 ) and emission (q = q 1 + q 2 ) channels (Fig 6(a)) while 4ph scattering (Fig 6(b)) comprises of three scattering channels namely (a) recombination (q + q 1 + q 2 = q 3 ), (b) redistribution (q + q 1 = q 2 + q 3 ) and (c) splitting (q = q 1 + q 2 + q 3 ) along with their momentum destroying counterparts via the reciprocal lattice vector G. Here q, q 1 , q 2 , q 3 define phonon wave vectors. Figure6(a) depicts that the acoustic phonons with low frequency modes (regime 'A') enable absorption processes which seems to gradually decrease as optical phonons overlap with the acoustic ones in the regime 'M'. In regime 'O' (60 cm -1 < ω < 174 cm -1 ), optical modes with higher frequencies are more prone to undergo emission processes which dominates at higher frequencies while absorption phase space going to zero (Fig6(a)). The aforementioned 3phminima are found to originate as a result of the crossing between absorption and emission phase space around 90-100 cm -1 . Also, Fig6(b)shows that 4ph redistribution phase space peaks around the same frequency where 3ph phase space shows minima. This is a result of the occurrence of flat optical bands in this frequency regime (86 cm -1 < ω < 122 cm -1 ) that severely restrict OOO and weaken AOO processes but do not prevent optical phonons to participate in four-phonon processes. Rather flat bands enhance the redistribution processes as energy conservation between four phonons is preferable in this scenario. Similar feature was observed in AlSb and BAs[36]. This causes a simultaneous dip and peak in 3ph and 4ph phase space respectively. Recombination process is seen to dominate the whole 4ph phase space spectrum of Ge 2 Sb 2 Te 5 along with a small contribution from recombination and splitting processes at low and high frequencies respectively (Fig6(b)). These 3ph and 4ph

  Fig 6(a)) while 4ph scattering (Fig 6(b)

Figure 7

 7 FIG. 7. Thermodynamic average of normal (N) and Umklapp (U) phonon scattering rates (Γ i = ⟨τ -1 i ⟩ave) are presented as a function of temperature for both (a) three-phonon (Γ 3ph ) and (b) four-phonon (Γ 4ph ) scattering events. For each of the cases, R (black dashed line) denotes the total resistive scattering rate involving Umklapp (U) and isotope (I) scattering.

FIG. 9 .

 9 FIG.8. The variation of 3ph average scattering rates (Γ 3ph ) with temperature for different scattering channels (absorption and emission) for both normal and Umklapp scattering events.

  observed in Fig 6(b). Looking at the phonon-mode contributions at different frequencies, we find that the 'M' frequency regime is seen to be mostly populated by redistribution as well as recombination processes that exhibit U scattering (supplementary Fig S10). In 'O' frequency regime, recombination channels are mostly blocked while splitting processes are found to gain importance along with the dominating redistribution scattering rates. At very high frequencies splitting processes aggravate due to the increasing weighted scattering phase space (see supplementary Fig S7(c), (d)). At higher temperature, marginal change is observed on this frequency distribution of scattering processes (Fig S7(d) and Fig S10(b), (d), (f) in the Supplementary information).

  a) represents the variation of aver-age lattice thermal conductivity (κ ave L ) with temperature, distinguishing different contributions coming from acoustic (A), optical (O) and mixed (M) frequency regimes defined via dispersion bands. We see a similar optical phonon dominance in κ L as that of the earlier case considering only 3ph and isotope scattering (shown in Fig 2(b)). At low temperature, as expected, acoustic contribution (red circles in Fig 10(a)) is large and the contributions from acoustic and mixed modes (regime M, green triangles in Fig 10(a)) dominate. Beyond 100 K, the dominance of acoustic modes decays while optical contribution (blue diamonds in Fig 10(a)) takes the lead with almost equal contributions from the mixed regimes.

FIG. 10

 10 FIG. 10. (a) The variation of average lattice thermal conductivity with temperature including 3ph, 4ph as well as isotope scattering processes, distinguished via contributions from acoustic, mixed and optical frequency regimes for Ge2Sb2Te5. (b) Average lattice thermal conductivity as a function of temperature for Ge2Sb2Te5. Three different curves correspond to the inclusion of different scattering processes; Red solid line: Inclusion of three-phonon (3ph) iteratively and the phonon-isotope scattering, blue dashed line: three-phonon using RTA and isotope scattering, black solid line with diamonds: three-phonon in RTA, four-phonon in RTA and isotope scattering. Experimental values of κL, obtained by Li et al. [25](κL = 0.47 W/mK) and Lee et al. [71](κL = 0.42 W/mK) are also presented via orange circle and green pentagon respectively.

  FIG. 11. Frequency variation of cumulative lattice thermal conductivity (κ c L ) using 3ph+I (blue solid lines) as well as 3ph+4ph+I (red dashed lines) for (a) T = 30 K, (b) T = 60 K, (c) T = 100 K, (d) T = 300 K, (e) T = 500 K and (f) T = 600 K. Blue, green and grey shaded regions denote A, M and O frequency regimes. The percentage contributions to κL coming from these regimes, are presented via blue (3ph+I) and red (3ph+4ph+I) for each T .

  Fig 13. Also, we represent the effect of

FIG. 12 .

 12 FIG.12. The variation of maximum mean free path (λ max ) as a function of temperature for Ge2Sb2Te5 both including and excluding 4ph scattering. Inset shows the scaling relation of λ max with T with (3ph+4ph+I) or without (3ph+I) 4ph scattering.

FIG. 13 .

 13 FIG. 13. The percentage contribution of cumulative lattice thermal conductivity (κ c L ) to that of the total κL is presented as a function of maximum allowed phonon mean free paths (MFP) including (a) 3ph, I scattering and (b) 3ph, 4ph and I scattering for different temperature. The dashed lines in both representations correspond to the 80 % contribution to κL. The corresponding maximum allowed MFPs for T = 30 K, 60 K, 100 K, 300 K, 500 and 600 K are denoted via a1, a2, a3, a4, a5, a6 in (a) and b1, b2, b3, b4, b5, b6 in (b). For comparison, the values correspond to T = 300 K and 600 K are shown via a4, a6 respectively in (a) and via b4, b6 respectively in (b). The blue shaded regimes in (a) and (b) represent the contributions from MFPs less than or equal to the lattice parameter a (4.193 Å).

FIG. 14 .

 14 FIG. 14. The variation of the anisotropy factor (κ xxL /κ zz L ) with temperature for crystalline Ge2Sb2Te5 using various scattering processes to calculate κL. Squares: 3ph (iterative) + isotope scattering, diamonds: 3ph (RTA) + isotope scattering, circles: 3ph (RTA) + 4ph (RTA) + isotope scattering.

FIG. 15 .

 15 FIG. 15. The variation of the anisotropy factor (κ c L (xx)/κ c L (zz)) with frequency for crystalline Ge2Sb2Te5 at different temperature. Blue dashed lines and black solid lines denote the scattering processes by 3ph+I and 3ph+4ph+I respectively. Different frequency regimes A, M, O are also shown via shaded regions. The dashed line in each case refers to the isotropic behavior in κ c L (xx)/κ c L (zz) beyond which anisotropy is measured.

2 2

 22 Sb Te 5 (Fig 1(C)), indicating an anisotropy in v g between its in-plane (a-b plane) and out-of-plane (along c) components in the hexagonal Brillouin zone.

F.

  Thermal conductivity scaling at high temperature: Relation between scattering rates and fundamental constants

  FIG. 17. (a) Total three-phonon and four-phonon scattering rates by summing up the contributions from all the modes (τ -1 all modes ), as a function of temperature. T -scaling for 3ph and 4ph scattering rates at high temperature are also shown via blue (T 0.97 ) and black fitted curves (T 2.16 ) respectively. (b) The variation of 3ph scattering rates (τ -1 3ph ) with temperature for each phonon mode shows ∼ T scaling. (c) The variation of 4ph scattering rates (τ -1 4ph ) with temperature for each phonon mode shows ∼ T 2 scaling. (d) 3ph (dashed lines) and 4ph (solid lines) scattering rates are presented for first six phonon bands. Scattering rates for both acoustic and optical bands are seen to follow ∼ T for 3ph and ∼ T 2 for 4ph scattering processes. Blue: acoustic phonons, Red: Optical phonons. For the lowest lying acoustic bands for 3ph and 4ph scattering, dashed and solid lines are merged and are indistinguishable.

FIG. 18 .

 18 FIG. 18. (a) The variation of inverse lattice thermal conductivity (1/κL) with temperature both including (red circles) and excluding four-phonon scattering processes (grey pentagons). The linear and non-linear trends of 1/κL as a function of temperature are represented via brown solid (∼ T ), blue solid (∼ T + T 2 ) and black dashed lines (∼ T 1.3 ). (b) The variation of inverse lattice thermal diffusivity (D -1 ) with temperature both including (red circles) and excluding four-phonon scattering processes (grey pentagons). The linear and non-linear trends of D -1 as a function of temperature are represented via blue (∼ T ) and black solid (∼ T + T 2 ) lines. in the kinetic transport regime of Ge 2 Sb 2 Te 5 considering three-phonon and isotope scattering as shown in Fig 18(b). Fitting the lower bound approximation using v s (= 3187.65 m/s) and τ p as known parameters for Ge 2 Sb 2 Te 5 yields s = 2.79 (see blue solid line in Fig 18(b)). At temperature higher than θ D (=100 K for Ge 2 Sb 2 Te 5 ), heat capacity is almost constant. Therefore, the lower bound of D, featuring a relaxation time which is inverse of temperature, stems from the T -1 decay of κ L in the kinetic regime of phonon transport. This suggests a physical picture of phonons performing random walk to create a completely diffusive dynamics under dominating Umklapp scattering in this kinetic regime. However, this theory faces a serious challenge while in-

FIG. 19 .

 19 FIG.19. Dispersion relation for optical bands in hexagonal-Kooi Ge2Sb2Te5. The optical bands are designated from 1 to 24 gradually from lowest to highest frequency.

2 2

 22 Sb Te 5 , especially in the 'M' regime, as is evident from our frequency variation of mode Gruneisen parameter (γ) (shown in supplementary FigS1(a)). For these few modes within 'M' regime, owing to the strong overlap between acoustic and optical bands, γ becomes anomalously large. However, mode-averaged γ smooths out this feature at a particular temperature and therefore, the average γ at fixed T behaves in less anomalous way (supplementary FigS1(b)). This anharmonicity seems to originate from the weak Te-Te bonds in Ge 2 Sb 2 Te 5 , which allow substantial frequency modulation of Ge 2 Sb 2 Te 5 with volume. Indeed, rattling dynamics of loosely bounded atoms were discussed by Xie et al.[44] for low-κ material AgCrSe 2 , where four-phonon scattering rates due to TA modes were found to satisfy the overdamped limit. Further, scattering channel distinction reveals that these modes result from the 3ph absorption process via normal scattering (see supplementary FigS9(b)).

FIG

  FIG. S1. (a) Gruneisen parameter (γ) for the Kooi-Ge2Sb2Te5 as a function of phonon frequencies. Acoustic, mixed and optical regimes, discussed in the main text, are also shown via A, M and O respectively. (b) The temperature variation of γ, calculated as the weighted sum of the mode contributions.

Figure

  Figure S1(a) presents the Gruneisen parameter (γ) for the Kooi-Ge 2 Sb 2 Te 5 as a function of phonon frequencies.FigureS1(b)shows the temperature variation of γ, calculated via the weighted sum of the mode contributions, at each temperature. γ varies in a wide range from -9 to ∼ 51, however the weighted sum over the contribution from different modes reveals that γ increases from 2.37 at T = 30 K to 2.49 at T = 100 K and gets saturated to a value of 2.49 throughout the high temperature regime beyond 100 K. Some sporadic points with higher values of γ are observed in the 'M' regime, corresponding to strong anharmonicity. This anharmonicity seems to come from the weak Te-Te bonds in Ge 2 Sb 2 Te 5 , which allows substantial frequency modulation of Ge 2 Sb 2 Te 5 with volume. These few strongly anharmonic modes feature overdamped phonon dynamics as described in Ap-

  Figure S1(a) presents the Gruneisen parameter (γ) for the Kooi-Ge 2 Sb 2 Te 5 as a function of phonon frequencies.FigureS1(b)shows the temperature variation of γ, calculated via the weighted sum of the mode contributions, at each temperature. γ varies in a wide range from -9 to ∼ 51, however the weighted sum over the contribution from different modes reveals that γ increases from 2.37 at T = 30 K to 2.49 at T = 100 K and gets saturated to a value of 2.49 throughout the high temperature regime beyond 100 K. Some sporadic points with higher values of γ are observed in the 'M' regime, corresponding to strong anharmonicity. This anharmonicity seems to come from the weak Te-Te bonds in Ge 2 Sb 2 Te 5 , which allows substantial frequency modulation of Ge 2 Sb 2 Te 5 with volume. These few strongly anharmonic modes feature overdamped phonon dynamics as described in Ap-

  FIG. S2. The variation of three-phonon (3ph), four-phonon (4ph) and phonon-isotope scattering (I) rates with phonon frequencies at different temperatures, denoted via blue, red and grey respectively: (a) T = 30 K, (b) 60 K, (c) 100 K, (d) 300 K, (e) 500 K, and (f) 600 K. Different phonon frequency regimes are marked using white (regime 'A'), light blue (regime 'M') and light red (regime 'O').

  FIG.S3. The variation of specific heat as a function of temperature for Ge2Sb2Te5.

  the manuscript. We refer to high temperature when T is considerably higher than θ D . The variation of specific heat as a function of temperature is shown below. Almost constant C v at temperature beyond 300 K is found (FigS3), satisfying the Dulong-Petit law. This result also matches considerably well with C v obtained via recent TDTR (time-domain thermoreflectance) experiments on hexagonal Ge 2 Sb 2 Te 5 by Li et al.[25].

4 .

 4 FIG. S4. (a) Phonon density of states (PDOS) is presented as a function of phonon frequencies for crystalline hexagonal Ge2Sb2Te5 stacked in Petrov structure. (b) Phonon dispersion relations of Petrov-Ge2Sb2Te5 along high symmetry path Γ-M -K-Γ-A.

5 .

 5 FigureS5presents the evolution of three-phonon (3ph) and four-phonon (4ph) phase space with frequency in Petrov-Ge 2 Sb 2 Te 5 . Although 3ph scattering phase space is dominant over 4ph scattering phase space, in a small frequency range within optical bands (shown via the grey shaded regime in Fig S5), 4ph phase space exceeds its 3ph counterpart. Also, this frequency window records the minima of the 3ph scattering phase space. This feature of Petrov-Ge 2 Sb 2 Te 5 is similar to that of the Kooi structure (discussed in the main text) although in a less pronounced fashion. We also observe similar but fewer flat bands in this regime via the dispersion relation compared to that of the Kooi-Ge 2 Sb 2 Te 5 . FigureS6represents the relative strength between 3ph and 4ph weighted scattering phase spaces in Petrov-Ge 2 Sb 2 Te 5 as a function of temperature ranging from 30 K to 600 K. As the temperature is raised, we find a gradually increased relevance of 4ph scattering phase space compared to that of the 3ph phase space. This stems from the fact that the evaluation of weighted phase space takes the phonon distribution function (f 0 ) into account which varies linearly with temperature at high temperature limit (f 0 ∼ k B T /ℏω). This indicates that the thermal transport of Petrov-Ge 2 Sb 2 Te 5 can also be influenced by incorporating the 4ph scattering processes.

Figure

  Figure S7 represents the variation of weighted phase space coming from distinct scattering channels for both 3ph (WP 3 , Fig S7(a), (b)) and 4ph scattering processes (WP 4 , Fig S7(c),(d)) at T = 30 K and 600 K. We observe that the distributions of WP 3 and WP 4 almost remain unchanged in the whole temperature range studied. In the frequency regime A, WP 3 is seen to be dominated by the absorption processes which soon starts decaying in the regime M, after which emission seems to dominate at higher frequencies in the O regime. At higher temperature, WP 3 increases and the crossover between absorption and emission processes is seen to be marginally postponed at higher frequencies. Four-phonon scattering is found to be dominated mostly by the redistribution process throughout the frequency range (Fig S7(c), (d)). Only at frequencies beyond 130 cm -1 , splitting process takes over. We note that this is the same regime where optical bands become markedly dispersive and sparsely populated. In frequency regime O beyond ω > 130 cm -1 ,

7 .

 7 Relative strength between Umklapp and Normal scattering in

8 .

 8 Figure S9 explains three-phonon (3ph) absorption (Fig S9(a), (b)) and emission processes (Fig S9(c), (d)) for two extreme temperature points (T = 30 K and 600 K).As temperature is raised from lowest (30 K) to highest value (600 K), for all the cases, scattering rates increase ∼ 60 and ∼ 11 folds for absorption and emission respectively. In the absorption process, few N scattering modes are found to contribute higher 3ph-scattering rates in the M regime (Fig S9(a), (b)). At high T (600 K), these few sporadic N absorption modes feature overdamped phonon dynamics (ωτ < 1) as discussed in Appendix B. We also observe that τ -1 3ph is always larger for emission processes than that of the absorption processes regardless the temperature. Further a more quantitative picture of the relative importance between N and U scattering is obtained via the thermodynamic average of scattering rates as described in the article, which suggests U scattering dominance in 3ph scattering. FigS10presents the frequency distribution of fourphonon scattering rates in different scattering channels, distinguished as N and U processes. Redistribution processes in general, resistive U redistribution to be particular (q + q 1 = q 2 + q 3 + G), are seen to overpower at both high and low temperature. At low temperature (FigS10(a), (c), (e)), the M frequency regime is seen to be populated by mostly redistribution as well as recombination processes that exhibit U scattering. In O frequency regime, recombination channels are mostly blocked while splitting processes are found to gain importance along with the dominating redistribution scattering rates. At very high frequencies splitting processes aggravate due to the increasing weighted scattering phase space. This frequency distribution remains almost similar at higher temperature (Fig S10(b), (d), (f)). We realize that both 3ph and 4ph scattering processes are predominantly of Umklapp type in Ge 2 Sb 2 Te 5 which is more evident from the temperature variation of the average scattering rate, discussed in the main text.

  FIG. S9. The frequency variation of 3ph normal (N) and Umklapp (U) absorption scattering rates are presented for (a) T = 30 K and (b) T = 600 K. Similarly, frequency variation of 3ph normal (N) and Umklapp (U) emission scattering rates are shown for (c) T = 30 K and (d) T = 600 K. A, M and O frequency regimes are also distinguished.

9 .

 9 Figure S11 presents the cumulative κ L as a function of frequency for crystalline Ge 2 Sb 2 Te 5 , considering 3ph and isotope scattering events at different temperature. The evolution of κ L along a-b plane (κ x L , shown via red dashed lines) and c axis (κ z L , shown via yellow dashed lines) suggests a strong anisotropic thermal conduction behavior in Ge 2 Sb 2 Te 5 . The average κ L (black lines in Fig S11(a)-(f)) shows a gradual domination of M and O regimes to dictate κ L as T is increased. At high temperature beyond 300 K, the contributions to κ L are saturated and A, M, and O regimes contribute to 26.4 %, 33.2 % and 40.4 % of κ L . Similar feature of anisotropic κ L is retained in the four-phonon framework where we employ 4ph, 3ph and I scattering to get the cumulative κ L for Ge 2 Sb 2 Te 5 (Fig S12). The optical phonon dominance in κ L is also retained at high temperature with a marginal modification of contributions coming from different frequency regimes as follows: regime A: 26.2 %, regime M: 35.8 % and regime O: 38 % of κ L (Fig S12).

10 .

 10 Figure S13 presents cumulative κ L as a function of maximum allowed phonon mean free path (MFP) in Ge 2 Sb 2 Te 5 , including 3ph and I scattering processes. Similar to that of the Fig S11, anisotropic thermal conduction is reflected via the evolution of cumulative κ L with mean free paths along a-b (blue dashed lines in Fig S13(a)-(f)) and along c (green dashed lines in Fig S13(a)-(f)). The maximum allowed mean free paths are identified at each T , which is found to decrease with T (Fig S13). Similar features are retained while including 4ph scattering Fig S14(a)-(f)). We observe that increasing T from 30 K to 600 K decreases λ max from 3.13 µm to 0.08 µm (80 nm) with 3ph and I scattering (FigS13). This effect is even more severe while including 4ph scattering where λ max decreases from 2.15 µm to 0.025 µm (25 nm) while increasing T from 30 K to 600 K (FigS14).

11 .

 11 Convergence of three-phonon and four-phonon processes in

  FIG. S11. Frequency variation of cumulative lattice thermal conductivity including 3ph and phonon-isotope scattering for (a) T = 30 K, (b) T = 60 K, (c) T = 100 K, (d) T = 300 K, (e) T = 500 K and (f) T = 600 K. For all the cases, black lines denote average κL, while red and yellow dashed lines represent κL along x and z respectively. The contributions of different frequency regimes (A, M, O) to the cumulative lattice thermal conductivities (in percentage) are presented for each temperature via the shaded zones.

  FIG. S12. Cumulative lattice thermal conductivity including 3ph, 4ph and phonon-isotope scattering is presented as a function of frequency for (a) T = 30 K, (b) T = 60 K, (c) T = 100 K, (d) T = 300 K, (e) T = 500 K and (f) T = 600 K. For all the cases, black lines denote average κL, while blue and green dashed lines represent κL along x and z respectively. The shaded zones for each temperature represent the contributions of different frequency regimes (A, M, O) to the cumulative lattice thermal conductivities (in percentage).

  FIG. S13. Cumulative lattice thermal conductivity including 3ph and phonon-isotope scattering processes as a function of maximum allowed mean free path for (a) T = 30 K, (b) T = 60 K, (c) T = 100 K, (d) T = 300 K, (e) T = 500 K and (f) T = 600 K. For each temperature, black lines denote average κL, while blue and green dashed lines represent κL along x and z respectively. The red vertical lines at each temperature indicate the maximum mean free path (λmax) that contributes to κL.

  FIG. S14. Cumulative lattice thermal conductivity including 3ph, 4ph and phonon-isotope scattering processes as a function of maximum allowed mean free path, for (a) T = 30 K, (b) T = 60 K, (c) T = 100 K, (d) T = 300 K, (e) T = 500 K and (f) T = 600 K. For each temperature, black lines denote average κL, while red and yellow dashed lines show κL along x and z respectively. The blue vertical lines at each temperature indicate the maximum mean free path (λmax) that contributes to κL.

  FIG. S15. (a)The variation of average lattice thermal conductivity with mesh size, considering 3ph and isotope-phonon scattering. The comparison between iterative and RTA method is also shown. (b) The total 4ph phase space volume is presented as a function of q-mesh grid.

TABLE I .

 I Comparison between equilibrium lattice parameters of differently stacked hexagonal Ge2Sb2Te5 obtained in this work and in earlier Ab initio studies.

		Kooi
	Lattice parameters ( Å) Ref. [29] Ref. [15] This work
	a	4.191	4.23	4.193
	c	17.062	16.88	17.031
		Petrov
	Lattice parameters ( Å) Ref. [29] Ref. [15] This work
	a	4.178	4.2	4.17
	c	17.41	17.14	17.37

TABLE II .

 II Theoretical and experimental values of various parameters used to calculate lattice thermal conductivity using the model proposed by Slack and Morelli[32,[START_REF] Morelli | High thermal conductivity materials[END_REF].

	Ge2Sb2Te5 (Kooi)	Parameter details	Value(s)
	V0 ( Å3 )	Volume of the unit cell	259.28 [this study]
	δ( Å)	(V0/atom) 1/3	3.066 [this study]
	θD(K)	Debye temperature	100 [47]
	γ	Gruneisen parameter	2.49 [this study]
	n	Number of atoms per unit cell	9
	MGe (g/mol)	Atomic mass of Ge	72.64
	M Sb (g/mol)	Atomic mass of Sb	121.76
	MT e (g/mol)	Atomic mass of Te	127.6
	M (g/mol)	Average atomic mass	114.09
		of Ge2Sb2Te5
	FIG. 17

  .

	Optical bands id	Band width(∆ωO (cm -1 ))	Allowed 3ph processes
	1	34.09	AAO
	2	34.35
	3	33.84	AOO (only for ωA ≤ ∆ωO)
	4	35.64
	5	33.84
	6	28.15	OOO
	7	35.64
	8	20.14
	9	13.43	AAO (only for ωO ≤ 2∆ωA)
	10	22.74
	11	14.73
	12	11.37
	13	9.04
	14	7.23	AOO (weak, only for ωA ≤ ∆ωO)
	15	7.49
	16	10.59
	17	9.82
	18	10.07
	19	7.86
	20	19.31
			AAO (weak)
	21	28.73
	22	25.27	AOO (only for ωA ≤ ∆ωO)
	23	29.71
	24	27.56	OOO(only for ∆ωO ≥ ω min O	)
	Appendix B: Overdamped phonon modes in
	acoustic-optic mixed ('M') regime

ACKNOWLEDGMENTS

We acknowledge the help of Cluster Curta at the University of Bordeaux for high-performance computing. This project had received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 824957 ("Before-Hand:" Boosting Performance of Phase Change Devices by Hetero-and Nanostructure Material Design).

 III). Consequently, following [83], we present the selection rules for 3ph processes concerning optical phonons (AAO, AOO, OOO) below: