Matryoshka phonon twinning in α-GaN

Understanding lattice dynamics is crucial for effective thermal management in electronic devices because phonons dominate thermal transport in most semiconductors. α-GaN has become a focus of interest as one of the most important third-generation power semiconductors, however, the knowledge on its phonon dynamics remains limited. Here we show a Matryoshka phonon dispersion of α-GaN with the complementary inelastic X-ray and neutron scattering techniques and the first-principles calculations. Such Matryoshka twinning throughout the basal plane of the reciprocal space is demonstrated to amplify the anharmonicity of the related phonons through creating abundant three-phonon scattering channels and cutting the lifetime of affected modes by more than 50%. Such phonon topology contributes to reducing the in-plane thermal transport, thus the anisotropic thermal conductivity of α-GaN. The results not only have implications for engineering the thermal performance of α-GaN, but also offer valuable insights on the role of anomalous phonon topology in thermal transport of other technically semiconductors. Lattice dynamics can govern the thermal conductivity solids and understanding the underlying mechanisms may enable enhanced performance of devices via judicious engineering. Here, insight into the anisotropic thermal conductivity of semiconducting α-GaN is gained by measuring the Matryoshka phonon dispersions.


Results and discussion
Matryoshka phonon twinning behavior in α-GaN. The phonon dispersions of α-GaN were measured by IXS below the phonon gap around 45 meV, as shown in Fig. 1a for 300 K (the dispersion measurements at 50 and 175 K can be found in Supplementary Fig. 3 and Supplementary Note 3). The measurements are in excellent agreement with first principles calculations and previous IXS results 9 . Most importantly, two parallel sections were found along Γ-M and Γ-K directions in the Brillouin zone (BZ), showing dispersion nesting behavior 35 : the near-constant energy difference is around 16 meV along Γ-M direction, while it is around 15 meV along Γ-K direction. The upper nested branch consists of low-energy optical (LEO) phonons in the lower q range and longitudinal acoustic (LA) phonons in the higher q range (referred to as the "Arc" branch in the following text), while the lower nested branch is one of the transverse acoustic branches (TA 2 ). and Supplementary Note 3). It can be clearly observed that, along both Γ-M and Γ-K directions, the energies of these phonon modes exhibit similar dispersive behaviors with increasing q, leading to a nearly parallel section between the branches.
An INS experiment was performed to measure the full phonon lattice dynamics in the reciprocal space to complement the IXS data, which only cover selected high symmetry directions. The acquired dynamical susceptibility, χ″(Q, E), is shown in Fig. 2 as twodimension slices along Γ-M, Γ-K, and Γ-A directions at 14 (Fig. 2a, d, g) and 300 K (Fig. 2b, e, h), overlaid with the first principles phonon dispersion calculation. The χ″(Q, E) calculations (Fig. 2c, f, i) based on first-principle phonon eigenvectors agree well with the INS measurements (Fig. 2b, e, h). The nesting behavior of the TA 2 and the Arc branches can be clearly observed along Γ-M direction in both the experimental and the calculated χ″(Q, E). The nesting behavior of the TA 2 and the Arc branches can also be observed along Γ-K, although the intensity of χ″(Q, E) is weaker due to the scattering structure factor.
Jointly analyzing the IXS and INS results, it is remarkable that the phonon dispersion nesting behavior is exhibited throughout the basal plane in the BZ of α-GaN (Supplementary Fig. 1b and Supplementary Note 1), as visualized by the volume and cross-section views of the measured χ″(Q, E) at 300 K in Fig. 3a-d (the results at 14 K are shown in Supplementary Fig. 6 and Supplementary Note 4). Similar to the phonons along the high symmetry directions in Fig. 2, the phonons along the non-high symmetry directions exhibit similar nesting behavior (additional cuts and visual angles can be found in Supplementary Fig. 7 and Supplementary Note 4). Such phenomenon can be vividly illustrated by the phonon dispersion surfaces in Fig. 3e from the first principles calculation (also shown by the calculated phonon dispersion in Supplementary Fig. 8 and Supplementary Note 5). In the basal plane of the reciprocal space, the nested TA 2 and Arc phonon dispersion surfaces are similarly cone-shaped and nest together in a Matryoshka-like twinning behavior (Fig. 3e). Such Matryoshka twinning is reflected by both the first principles calculations and INS projection contours (Fig. 3f, g and Supplementary Fig. 8 and Supplementary Note 5). Similar to the case of CuCl where the phonon nesting arises from the ionic radius mismatch 36 , the large ionic radius mismatch of N (0.75 Å) and Ga (1.62 Å) atoms are also responsible for the Matryoshka phonon twinning in α-GaN. In addition, such Matryoshka behavior may also be attributed to the two [GaN4] tetrahedrons along c-axis in the crystal structure. Unlike the monolayer GaN 3 , the double tetrahedral units along c-axis induce the TA branch folding from the zone boundary (A-point), and thus lead to a low-energy of the Arc at Γpoint. Such Matryoshka twinning may provide vast three-phonon scattering channels in α-GaN by reducing the momentum transfer constraint, thus expanding the phonon scattering phase space. For example, in a three-phonon emission process where one phonon mode on the Arc branch (q Arc , ω Arc ) decays into two phonon modes (q 1 , ω 1 ) and (q 2 , ω 2 ), for the mode (q 1 , ω 1 ) on the Arc branch near Γ point, we can always find a TA 2 mode (q 2 = q Arc − q 1 , ω 1 = ω Arc − ω 1 ) in the basal plane that enables such scattering (Supplementary Fig. 9 and Supplementary Note 6). This scattering behavior is consistent with recent report of weighted phonon scattering space in α-GaN, which shows the emission process is significantly larger than that of the absorption process in the energy range from 20 to 40 meV (~5 THz to~10 THz), (Supplementary Fig. 10a and Supplementary Note 7) 6 in line with the energy range of the Arch branch in the Matryoshka twinning.
The impact of Matryoshka phonon twinning on the anharmonic scattering rate. Phonon linewidths of α-GaN at 300 K were extracted from IXS by de-convoluting the energy and momentum resolution functions, INS by de-convoluting the energy resolution function, and Raman spectroscopy (at Γ-point) without considering the instrument resolution. The INS results are moderately in consistent with the IXS measurement (Supplementary Fig. 10b and Supplementary Note 7). The Raman spectroscopy (https://chenwli.com/facilites/) gives a linewidth of 0.8 meV for the Arc branch at the Γ-point (Fig. 4d), where the linewidth is also found to be narrow by INS. Although the measured phonon linewidths vary between different techniques because of their different resolution functions, it is reasonable to compare them within one technique to quantify the linewidth difference between the phonons related to Matryoshka twinning and other branches. Phonon linewidth is proportional to the total scattering rate from all processes 35,37 and dominated by anharmonic phonon-phonon interactions in α-GaN. As shown in Fig. 4a, while the linewidths of most phonon modes are between 0.5 and 1.5 meV, the linewidths of the phonon modes on the TA 2 and the Arc branches are much broader and double of some other phonon modes. These large linewidths indicate that the nested phonons are scattered more strongly than the others, leading to the phonon lifetime of nested modes is cut by more than 50%. In Fig. 4b, c, it is found that the linewidths of phonon modes on the TA 2 and Arc branches (low-q in TO and high-q in LA) are much larger than those on the high energy optical (HEO) branch along Γ-M direction and those on the HEO and LEO branches along Γ-A direction. The anharmonic phonon scattering rate is expressed as 38 where V is the third-order interatomic potential, ω qs is the phonon frequency of the phonon mode of wavevector q and branch index s, and n is Bose occupation. This anharmonic scattering rate is determined by both the third-order interatomic potentials and the phonon scattering phase space: the former represents the scattering strength, and the latter represents the phase space of available scattering channels 39 . The large linewidths of the phonons involved in the Matryoshka twinning suggest that such behavior significantly promotes three-phonon scatterings in α-GaN.
To elucidate the anharmonicity of the phonon modes on the Arc branch in α-GaN, the evolution of the phonon energy at q = 0.1 along Γ-M direction was fitted over a wide temperature range with the following expression 40 where A and ω(0) are fitting constants, and ω A (T) is the temperature-dependent phonon energy. The result in Fig. 5a shows a moderate anharmonicity below the Debye temperature (~636 K) 11 , supported by the moderate difference between experimental isobaric 41 and calculated volumetric mode Grüneisen parameters 42 (Fig. 5b), see the details in Supplementary Note 8. Additionally, the frozen phonon potential for the phonon mode on the Arc branch at q = 0.1 along Γ-M direction (Fig. 5c) only slightly departs from the harmonic behavior, indicating moderate anharmonicity (the phonon eigenvectors are shown in the insert in Fig. 5c). Considering such moderate anharmonicity, it is confirmed that the large linewidths of the phonon modes in the TA 2 and Arc branches are dominated by the enlarged scattering phase space from the vast scattering channels induced by the Matryoshka twinning.
The impact of Matryoshka phonon twinning on the anisotropic thermal conductivity. The anisotropic thermal conductivity of α-GaN recently attracted intensive attention because it is critical for heat management in power electronics, where α-GaN is a key component. However, the answer remains controversial: due to the difficulty of the anisotropic thermal conductivity measurement 7 , previous theoretical calculations reported that the in-plane thermal conductivity is greater than the out-of-plane one while others showed opposite results 3,5,6,19,20 . Among these reports, three recent theoretical calculations showed that the in-plane thermal conductivity is lower. Qin et al. reported that the in-plane and the out-of-plane thermal conductivities are 280 and 325 W m −1 K −1 , with an anisotropy ratio of 14% 3 . Wu et al. 19 and Yang et al. 20 reported that the inplane and the out-of-plane thermal conductivities are 228 and 260 W m −1 K −1 , with an anisotropy ratio of 12%. These two reports both show good agreement with the measured value around 230 W m -1 K −18 and thus provide relatively reliable reference for the anisotropy of the thermal conductivity of α-GaN.
Understanding the temperature dependence of the anisotropic thermal conductivity requires considering both changes of the anharmonic phonon scattering rates and phonon dispersion/ group velocities 26 , as shown in the following equation: where C qs is the heat capacity, v qs the group velocity, and τ qs the phonon lifetime of wavevector q and branch index s. Usually, acoustic and low-energy optical phonons contribute more to the thermal conductivity than other less dispersive modes so only their contribution is considered here. The temperature-dependent acoustic phonon group velocities were extracted from the measured phonon spectra (Fig. 2), as shown in Table 1. The acoustic phonon group velocities along the out-of-plane (Γ-A) direction are slightly larger than those along the in-plane (Γ-M and Γ-K) directions. For the LA (TA 2 ) branch, the phonon group velocities along Γ-A are 1%, 1.7%, and 0.8% (3.7%, 3%, and 4%) larger than those along Γ-M at 14, 50, and 300 K, respectively. This behavior indicates that the anisotropy of acoustic phonon group velocity is minor and may not induce a sizable thermal conductivity anisotropy in α-GaN.
As discussed previously, the effects of in-plane Matryoshka twinning on the scattering rates of in-plane phonons are significant, shown by the anomalously large linewidths of the modes involved at 300 K. In Fig. 4, the Arc, TA 2 , and other in-plane branches have phonon linewidths around 2.3, 2, and 1.3 meV, respectively. The out- of-plane HEO and other branches (LEO, TA, and LA) have phonon linewidths around 1.5 and 1.2 meV, respectively. Considering the instrument resolution, the phonon lifetimes of TA 2 and the high-q part of LA are estimated to be <50% of their counterparts in the basal plane. Assuming similar group velocities for in-plane and out-ofplane phonon branches and using Eq. (3), the thermal conductivity contribution by TA 1 /LA, TA 2 , and Arc (TO) branches is roughly estimated to be 8%, 40%, and 48% lower, respectively, along the inplane direction than the out-of-plane direction. The Matryoshka twinning amplifies the acoustic-optical three-phonon scattering in α-GaN and thus reduces the in-plane thermal conductivity, especially at temperatures above the Debye temperature, when the lattice thermal  transport is dominated by the phonon-phonon interactions. Therefore, the anisotropy of thermal conductivity is mostly a result of the twinning behavior. If such Matryoshka phonon twinning could be suppressed by phonon engineering through strain or doping, the inplane thermal transport properties of α-GaN may potentially be enhanced for better thermal management at high temperatures.

Conclusion
Through a combination of inelastic X-ray and neutron scattering experiments and first principles calculations, a Matryoshka phonon twinning is found throughout the basal plane of the BZ in α-GaN. This phenomenon creates a vast number of three-phonon scattering channels and leads to enhanced anharmonic phonon scattering, shown by anomalously large phonon linewidths. The strong in-plane phonon scatterings due to the Matryoshka phonon twinning suppress the in-plane thermal conductivity of α-GaN and enhance the thermal transport anisotropy, while its acoustic phonon group velocities remain fairly isotropic. Amplification and suppression of such Matryoshka phonon twinning by phonon engineering could provide a valuable means to control lattice thermal transport in many related materials, such as electronic devices, thermoelectrics, and thermal barriers.

Methods
Sample preparation. High-quality GaN single crystals (un-doped, n-type, MTI Corporation https://www.mtixtl.com/GaN-C-50D025C2.aspx) used in this work were grown by the hydride vapor phase epitaxy (HVPE)-based method with low dislocation density (<1 × 10 7 cm −2 ) (Supplementary Fig. 1c and Supplementary Note 1). The quality of the crystals was checked with both X-ray and neutron diffraction. The full width at half maximum (FWHM) of the X-ray diffraction peak at (002) plane is around 0.10 ± 0.02°( Supplementary Fig. 1d and Supplementary Note 1).

Inelastic scattering experiments
Inelastic X-ray scattering measurements. High-resolution IXS experiment was performed to measure the phonon dispersions of an α-GaN single crystal (250 μm thickness, the lower panel of Supplementary Fig. 1c and Supplementary Note 1) at 50, 175, and 300 K. The measurements were conducted at 30-ID-C (the High-resolution Inelastic X-ray Scattering beamline, HERIX) at the Advanced Photon Source (APS) 43,44 . The incident photon energy was~23.7 keV with an energy resolution ΔE of 1.2 meV and a momentum resolution of 0.65 nm −1 . The single crystal was attached to a copper post by varnish and the copper post was mounted in a closed-cycle cryostat. The IXS measurements were accomplished at constant wavevector mode in reflection geometry. The orientation matrix was defined by using Bragg peaks at (4 0 0), (0 0 4), (0 0 5), and (2 2 0) respectively.
Inelastic neutron scattering measurements. The INS measurements were carried out on a bigger single crystal (φ = 5 cm, upper panel of Supplementary Fig. 1c and Supplementary Note 1) on the Wide Angular-Range Chopper Spectrometer, timeof-flight neutron spectrometer at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory 45 . An incident energy of E i = 50 meV and a Fermi chopper frequency of 420 Hz were used to optimize instrument energy resolution. At the elastic scattering, the energy resolution is 2 meV. For INS, the [110] axis was set vertical, and scattering plane was selected to obtain both the in-plane and outplane phonons at 14, 50, 300, and 630 K. When collecting the data, 1°step was used with a rotation range from −90°to 90°for 300 K, and 2°step was used with a rotation range from −70°to 50°for other temperatures. Because of the geometry of the substrate sample, the neutron absorption is very small when incident beam is near perpendicular to the substrate surface. The neutron absorption significantly reduces the counting statistics when incident and scattered neutron beams are near parallel to the substrate surface, but this do not affect the dispersion measurement. Data reduction was performed using the Mantid program 46 . The data were normalized by the accumulated incident neutron flux, and the detector efficiency correction was applied based on the incoherent scattering of the vanadium standard.
Scattering data processing. INS data shown in the present work were corrected for the sample temperature to obtain the dynamical susceptibility, χ″(Q, E) = S(Q, E)/[n T (E) + 1], where n T (E) stands for the Bose distribution, S(Q, E) is the fourdimensional scattering dynamical structure factor 47 where b d is the neutron scattering length, Q = k − k′ is the wavevector transfer, and k′ and k are the final and incident wavevector of the scattered particle. q is the phonon wavevector, ω s is the eigenvalue of the phonon corresponding to the branch index s, τ is the reciprocal lattice vector, d is the atom index in the unit cell, r d is the atom position, W d is the corresponding Debye-Waller factor, and e ds is the phonon eigenvectors. For given Q, the measured scattering spectra were fitted with a damped-harmonic-oscillator model 48 where M is the effective mass, ω 0 is the bare phonon energy in the absence of damping forces, and 2γ is a damping factor that describes the phonon scattering rates. n is the Bose distribution. The bare phonon energy and phonon linewidths were extracted by de-convoluting with both instrument energy and momentum resolution functions, the latter of which eliminates the slope effect of highly dispersive phonons on linewidths.
First principles calculations. The first principles calculations were performed based on the density functional theory as implemented in the Vienna Ab Initio Simulation Package 49 . The exchange-correlation energy was computed using the local-density approximation functional, and the projector-augmented-wave potentials were used (4s 2 4p 1 for Ga and 2s 2 2p 3 for N). For plane-wave expansion in reciprocal space, we used 600 eV as the kinetic energy cutoff value. The q-mesh was chosen as 7 × 7 × 5. The convergence criteria for total energy was set to 1 × 10 −8 eV, and that for atomic force was 10 −3 eV/Å. The structure was fully relaxed and the lattice constants, a = b = 3.188 Å, c = 5.190 Å, are slightly smaller than the experimental values (a = b = 3.191 Å, c = 5.191 Å obtained from HERIX at 300 K). The phonon dispersion of α-GaN at the level of harmonic approximation was calculated using the Phonopy package 50 with 4 × 4 × 3 supercell of α-GaN containing 192 atoms.

Data availability
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Figs. 1-10. Additional data related to this paper may be requested from the authors.