Giant Phonon Anharmonicity and Anomalous Pressure Dependence of Lattice Thermal Conductivity in Y2Si2O7 silicate

Modification of lattice thermal conductivity (κL) of a solid by means of hydrostatic pressure (P) has been a crucially interesting approach that targets a broad range of advanced materials from thermoelectrics and thermal insulators to minerals in mantle. Although it is well documented knowledge that thermal conductivity of bulk materials normally increase upon hydrostatic pressure, such positive relationship is seriously challenged when it comes to ceramics with complex crystal structure and heterogeneous chemical bonds. In this paper, we predict an abnormally negative trend dκL/dP < 0 in Y2Si2O7 silicate using density functional theoretical calculations. The mechanism is disclosed as combined effects of slightly decreased group velocity and significantly augmented scattering of heat-carrying acoustic phonons in pressured lattice, which is originated from pressure-induced downward shift of low-lying optic and acoustic phonons. The structural origin of low-lying optic phonons as well as the induced phonon anharmonicity is also qualitatively elucidated with respect to intrinsic bonding heterogeneity of Y2Si2O7. The present results are expected to bring deeper insights for phonon engineering and modulation of thermal conductivity in complex solids with diverging structural flexibility, enormous bonding heterogeneity, and giant phonon anharmonicity.

Phonon engineering, in terms of controlling heat flow within the lattice by manipulating phonon behavior, has became a powerful paradigm in the design and optimization of advanced materials with optimal thermal conductivity. The key concern is to modulate phonon transport by modifying different phonon scattering mechanism, including intrinsic phonon-phonon scatterings, and extrinsic scatterings by boundary and point defects 1 . In fact, significant reduction of thermal conductivity has been realized in a wide range of energy conservation technologies via material design strategies such as doping, nanostructuring and synthesis of superlattice [2][3][4][5] . However, the intricacy and uncertainty of these methods leaves open the necessity of searching for effective phonon engineering methodology based on manipulation of intrinsic structural characteristics.
Over the years, tailoring thermal transport in solids by applying hydrostatic pressure has been an attractive subject of experimental and theoretical investigations in a broad variety of materials from thermoelectrics and thermal insulators to minerals in mantle, because the phonon behaviors, which provide the detailed profile of lattice thermal conduction, is sensitive to volume change of the unit cell 6 . In a majority of materials, including diamond 7 , alkali halides 8,9 , ice 10,11 and mantle materials [12][13][14][15][16][17] , thermal conductivity has been found to increase with compressive pressure. Well documented knowledge in the literature states that, elevated frequency scale of phonons caused by strengthened atomic bonding tends to facilitate thermal conduction by increasing the acoustic phonon velocity and reducing the phonon scattering rate, which has long been approved as an universal feature of most solids. However, negative pressure dependence of lattice thermal conductivity (κ L ) has also been reported in recent years, which enlightens the novel strategy of bringing down κ L by applying compressive hydrostatic pressure with respect to manipulating phonon behaviors. By using first principles calculation combined with iterative solution of Boltzmann equation, Lindsay investigated the pressure dependence of κ L in some III-V, II-VI cubic compounds (such as BeTe) 18 , with a characteristic of significant acoustic-optic frequency gap created by large atomic mass ratio. In this case, the separation of longitudinal and transverse acoustic phonon branches Scientific RepoRts | 6:29801 | DOI: 10.1038/srep29801 under hydrostatic pressure raises the scattering rate among them and thus leads to decreased κ L . Nevertheless, the reported upward shift of each phonon branch to higher frequency range under pressure still complies with the traditional cognition, mainly due to the simple crystal structure and weak bonding heterogeneity of these studied binary compounds. Using similar methodology, Ouyang 19 observed a dκ L /dP < 0 relationship in HgTe with partially ionic atomic bonds, in which the reduction of relaxation time of transverse acoustic phonons considerably overwhelms the less enhancement of group velocity of longitudinal acoustic and optic modes. Such mechanism, however, cannot directly guide the investigations of complicated-structured ceramics with multiple chemical compositions and covalent and heterogeneous bonding features.
Recently, our density functional theoretical (DFT) calculations on high-pressure Raman spectroscopy of γ -Y 2 Si 2 O 7 (referred as Y 2 Si 2 O 7 for brevity) have identified obvious softening of low-frequency optic phonon modes. And this is much likely to account for significant acoustic-optic coupling, which has been recognized in other complex-structured materials with similar structural characteristics [20][21][22][23][24] . Intuitively, it is necessary to explore how this abnormal pressure-induced softening of optic phonons would affect κ L . The unexpected results might cast light on tuning κ L with respect to modulating the scattering mechanism among phonons, as well as the group velocity of acoustic phonons. Besides, considering the promising application of Y 2 Si 2 O 7 as high-temperature environment/thermal barrier coatings (ETBC) 25 , a comprehensive understanding of κ L on the basis of phonon engineering bares its values in both scientific understanding and technological application. In a broader perspective, the interlaced stacking of rigid [SiO 3 -O-SiO 3 ] units (or referred to as Si 2 O 7 where two adjacent SiO 4 tetrahedra are bridged by an central O atom in a Si-O-Si bond) and soft YO 6 octahedra in Y 2 Si 2 O 7 makes it a desirable model system for materials consisting of structural units in diverging flexibility, which covers large varieties of ternary/quaternary oxides.
In this study, we investigate the phonon dispersion and anharmonicity, and lattice thermal conductivity of Y 2 Si 2 O 7 at various compressive hydrostatic pressure using DFT calculations. The signature of a giant acoustic-optic coupling is identified with the character of some low-lying optic phonons interacting with acoustic phonons. Unexpected anomalous decrement of κ L is observed in pressured lattice; and the physical mechanism is disclosed as combined effects of 1) slightly decreased acoustic phonon group velocity and 2) significantly augmented scattering of acoustic phonons owing to pressure-induced enhancement of phonon anharmonicity. Additionally, we find that longitudinal acoustic phonons are more effective heat-carriers in Y 2 Si 2 O 7 and dominate the magnitude of decrement of total κ L under pressure; whereas transverse acoustic phonons are more intensely scattered by low-lying optic phonons and yield negligible contribution to thermal conduction.

Results
The crystal structure of Y 2 Si 2 O 7 remains stable up to the hydrostatic pressure of P = 3.06 GPa, equivalent to a volume contraction of 2%. In practice, the structural stability is guaranteed by steady changes of lattice parameters, bonding characteristics, configuration of structural units and elastic parameters. Pressure dependences of lattice constants and atomic bond lengths (See Supplementary Fig. S1 for details) indicate that deformation pattern of Y 2 Si 2 O 7 under hydrostatic pressure involves coordinated contraction of both Si-O and Y-O structural units, leading to packing densification. Besides, based on a statistical analysis on the configuration of specific Si 2 O 7 and YO 6 polyhedra, we find that the average length of Y-O bonds changes from 2.240 Å at P = 0 GPa to 2.225 Å at P = 3.06 GPa, whilst that of Si-O bonds changes from 1.626 Å to 1.621 Å. Thus, the variation percentage is larger for Y-O (0.67%) than Si-O (0.31%) bonds. Similarly, the angle of intra-octahedral O-Y-O bonds changes by as much as 1.01%, larger than that of intra-tetrahedral O-Si-O (no more than 0.64%). Evidently, YO 6 octahedron deforms in higher magnitude under pressure as if it is more "flexible" than Si 2 O 7 unit, and actively undergoes localized distortion to accommodate the loaded strain. Such unique response of structural feature to hydrostatic pressure will be further manifested in elastic properties, which will be discussed later.
Phonon dispersions -Giant acoustic-optic coupling in Y 2 Si 2 O 7 . Phonon dispersion curves along high-symmetry directions in Brillouin zone (BZ) are calculated at different pressure, and the results of two typical cases (P = 0 and 3.06 GPa) are plotted in Fig. 1(a). A general feature is the highly mixing of low-frequency optic phonons and longitudinal acoustic (LA) phonons at multiple BZ points, indicative of a giant acoustic-optic coupling in Y 2 Si 2 O 7 due to the close-up of acoustic-optic frequency gap. Particularly, a striking observation is the avoided crossing 26 between the lowest optic phonons (ν ~ 1.9 THz at Г point) and transverse acoustic (TA) phonons, seen as a repulsive trend between the two groups of curves. To explicitly elucidate this phenomenon, phonon dispersions (below 6 THz) of the equilibrium structure along the principal directions of reciprocal lattice are presented in Fig. 1(b), together with corresponding mode Grüneisen constant (γ i (q)) curves zoomed in for acoustic and three low-frequency optic phonons. As is depicted, the lowest optic phonon branch exhibits an apparent drop into the acoustic phonon regime at around the BZ center, and results in either flattening out or curving back of TA modes at approximately midway to the BZ boundary, leading to abrupt changes of calculated γ TA (q) due to avoided crossing between TA and low-lying optic phonons. On comparison, LA mode could better retain its feature towards the BZ boundary without drastic softening; and accordingly, calculated γ LA (q) curve exhibits as some small bumps, originated when LA cross the low-lowing optic phonons with slight fluctuation of its slope. Typically, optic phonons do not directly participate in thermal transport process due to the low group velocity, yet they provide essential scattering channels for heat-carrying acoustic phonons. Therefore, one could expect that the coupling between low-lying optic phonons and TA/LA phonons should play an important role in determining the thermal conduction in Y 2 Si 2 O 7 . One thing to notice is the much higher absolute values of γ for TA (above ~10) than those for LA (within ~3) phonons especially along Г→ Z path; which, based on the 1/γ 2 dependency of phonon relaxation time 27,28 , indicates that TA phonons have considerably higher anharmonicity and undergo stronger scattering than LA phonons. It is shown in Fig. 1(a) that under pressure of P = 3.06 GPa, the majority of high-lying optic phonons are raised to higher frequencies as a result of strengthened atomic bonding, analogous to the general trend in most solids. In contrast, the collective of low-lying optic and acoustic phonons shift downwards across entire BZ, which is an unconventionally interesting phenomenon; accordingly, γ i (q) (< 0 for most points) for all three acoustic modes moves to higher absolute values, indicative of enhanced acoustic phonon anharmonicity. Further comparison among acoustic phonons indicates that TA phonons are more susceptible to hydrostatic pressure than LA phonons, reflected by the higher magnitude of downward frequency shift as well as higher enhancement of γ TA (q) values. Here, we note that the softening of low-frequency phonons under hydrostatic pressure is a direct manifestation of the localized accommodation of applied strain, rather than the softening of entire lattice that implies structural instability.
Pressure dependence of intrinsic lattice thermal conductivity. Intrinsic lattice thermal conductivity of Y 2 Si 2 O 7 is calculated as a function of hydrostatic pressure based on modified Debye-Callaway model 27,28 . As is presented in Fig. 2(a), κ L of Y 2 Si 2 O 7 monotonously decreases upon hydrostatic pressure over the calculated temperature range; and the remarkable decrement (55.2%) of κ L at 600 K (around Debye temperature of Y 2 Si 2 O 7 , 575 K) sharply contradicts that of other referenced materials in Fig. 2(b) [7][8][9][10][11][12][13][14][15][16][17] . We also investigate the partial contributions of TA and LA branches to the total κ L and results are plotted in Fig. 2(c) for the end-pressure cases (P = 0 and 3.06 GPa). Of great interest, LA phonons contribute as much as 90% to the total κ L within studied temperature (pressure) range, declaring its overwhelmingly effective role over TA phonons in thermal conduction of Y 2 Si 2 O 7 . In fact, the trivial contribution of TA phonons could be fundamentally explained by stronger scattering with low-lying optic phonons, as we have observed from phonon dispersion curves. And such phenomenon has also been found in rare-earth pyrochlores with stacked polyhedra units 29 . Consequently, even though the partial κ L i of TA phonons decreases in higher proportion than LA (~71% decrement for κ L TA1 , 58% for κ L TA2 and 55% for κ L LA at P = 3.06 GPa), it is the pressure dependence of LA phonons that dominates the abnormal tendency as well as the magnitude of dκ L /dP < 0 relationship in Y 2 Si 2 O 7 .
In order to validate the reliability of modified Debye-Callaway model used in the present study, we herein compare our predicted κ L of Y 2 Si 2 O 7 with experimental data 30 measured for polycrystalline samples at ambient pressure. As is shown in Fig. 2(d), κ L calculated from present model manages to capture the κ L ~ 1/T dependency at around and above Debye temperature where Umklapp phonon scattering process dominates, and generally agrees well with experimental data. The slight discrepancy is understandable considering the fact that experimental data normally renders an over-flattened κ L -T relationship, i.e. defects contained in polycrystalline experimental samples might serve as effective phonon scatters at relatively low temperature regime in addition to the intrinsic Umklapp phonon-phonon scattering; while high-temperature measurements are susceptible to thermal radiation, which brings extra heat and thus result in higher or even ascending experimental values at high temperature region. Besides, we also calculate κ L of Y 2 Si 2 O 7 at a typical case (P = 0 GPa) from the combination of single-mode relaxation-time approximation and a full solution of linearized phonon Boltzmann equation (SMRT-LBTE), using DFT-calculated harmonic and anharmonic interatomic force constant (IFC) as inputs 31 . However, calculated κ L is seriously underestimated, which is a widely observed trend of this method originated from its basic methodology by taking phonon-phonon Normal scattering process as an independent channel for thermal resistance 32 . From above comparison, we could find that modified Debye-Callaway model is sufficiently reliable to reproduce the temperature dependence of κ L for Y 2 Si 2 O 7 .
Elucidating the driving force of the abnormal dκ L /dP < 0 in Y 2 Si 2 O 7 requires explicit examination of decisive phonon transport parameters, as listed in Table 1. Here, parameters are directly extracted from phonon  30 . Minimum thermal conductivity (κ min ) 52 at equilibrium geometry is presented in dash-dot line as a guide for the eye.
, and yields a 9.5% decrement up to P = 3.06 GPa. These results are consistent with the observed acoustic phonon softening upon hydrostatic pressure. Another key factor entering the calculation of κ L i is the frequency-dependent phonon relaxation time τ i U , and results for three acoustic phonons calculated at a typical temperature of 600 K are plotted in Fig. 3. As we could see, TA phonons in overall have overwhelmingly smaller relaxation time than LA phonons, with orders of difference in magnitude. Plus, values of τ i U for each acoustic phonons decrease under hydrostatic pressure, up to approximately 72% for TA1, 67% for TA2 and 54% for LA phonons. According to the inverse relation between phonon scattering rate and phonon relaxation time under relaxation time approximation (RTA), we come to a vivid screenshot of the phonon transport details in Y 2 Si 2 O 7 : 1) TA phonons are strongly scattered by low-lying optic phonons, which fundamentally explains its minor contribution to thermal conduction of Y 2 Si 2 O 7 ; 2) anharmonic phonon scattering is augmented under applied hydrostatic pressure, which could also be derived from tremendous increment of γ i for all the three acoustic phonons. To sum up, the predicted abnormal decrement of κ L in Y 2 Si 2 O 7 is originated from significantly augmented scattering of acoustic phonons induced by enhanced phonon anharmonicity, in corporation with slightly decreased acoustic phonon group velocity.
Pressure-induced augment of phonon scattering. In order to better understand the underlying physical mechanism of the abnormal dκ L /dP < 0 dependency in Y 2 Si 2 O 7 , we move on to specifically examine the key factor, i.e. scatterings of acoustic phonons. In dielectric crystals, the intrinsic κ L is governed by phonon-phonon scatterings rising from anharmonicity of interatomic potential 33 , and the most important resistive three-phonon scatterings involves two acoustic phonons combining with one optic phonon (aao) or with another acoustic phonon (aaa) 34,35 . Naturally, one could expect that the frequency change of both low-lying optic and acoustic phonons might bring significant modulation to different three-phonon scattering options, and thus determines heat transport process. It is shown in Fig. 1(a) that under P = 3.06 GPa, the overall decreased frequency scale within the lower part of phonon spectrum (ν ≦ 5 THz) increases the populations of low-lying optic and acoustic phonons 36 ; and meanwhile, the severe softening of TA phonons leads to a lower degree of "acoustic bunching" seen as the separation of TA further away from LA phonons, illustrated by grey arrows in Fig. 1(a). These features, to a qualitative extent, indicates an enlarged phase space for both (aao) and (aaa) scattering options 18,[34][35][36][37] , and thus points to a higher phonon scattering rate. Here we note that, such speculation based on analysis of phonon dispersion could explain the enhanced phonon scattering under pressure; nevertheless, future investigations (e.g. on energy-momentum surfaces, phonon scattering channels or phase space) 38 are required to fully clarify the interplay of different scattering options involving low-lying optic and acoustic phonons in complex ceramics, as well as its evolution upon modulation of structural features.
Structural origin of low-lying optic phonons and giant phonon anharmonicity. Clearly evidenced from above results, the low-lying optic phonons play a significant role in determining the thermal conductivity of Y 2 Si 2 O 7 as well as its pressure dependency, through severe coupling with heat-carrying acoustic phonons. Analysis of mode eigenvectors reveals that the low-lying optic phonons (with ν ≦ 5 THz) correspond with low-energy translation of Y atoms, coordinated by tilting or bending of rigid [O 3 Si-O-SiO 3 ] units. Distinguishment of the eigenvectors of these optic phonons from those lying in higher frequency range is presented in Supplementary Discussion and Supplementary Fig. S2. In fact, such unique vibration pattern is originated from the specific heterogeneity of interatomic environment in Y 2 Si 2 O 7 . To better elucidate this point, potential energy curves of Y and Si atoms are calculated by shifting the atoms away from their static equilibrium positions and calculate the energy difference with respect to their equilibrium state. As is shown in Fig. 4(a), Y atom (in heavier mass) sits in a very flat potential well and is loosely bonded with surrounding O atoms. In contrast, the potential well of Si atom (in lighter mass), which forms much stronger chemical bonds with the surrounding atoms, is over twice steeper than that of Y. As a result, the translation of Y atom could be easily exited with lower energy, coordinated by low-energy tilting or bending of [O 3 Si-O-SiO 3 ] bridged units, during which the configuration of end SiO 4 tetrahedra shows negligible changes. What worthy of mention is that, the inherent correlation between interatomic bonding heterogeneity and low-frequency optic phonons, contributed by coupled vibration of a set of atoms and thus run over a wide range in the lower part of phonon spectrum, has also been observed in other ternary materials showing a part-crystalline part-liquid state 39 .
Plus, it would be illuminating to dig out the structural origin of the pronounced enhancement of phonon anharmonicity in pressured Y 2 Si 2 O 7 . According to an empirical model treating crystal lattice as spring network 40 , the extent of lattice anharmonicity of materials has a positive correlation with the mismatch of force constants among neighboring bonds. There are four (six) nonequivalent Si-O (Y-O) bonds in the unit cell of Y 2 Si 2 O 7 , constituting a sample of ten different chemical bonds in our study. As is shown in Fig. 4(b), IFCs of Si-O bonds are approximately three times higher than those of Y-O bonds, indicative of higher strength of Si-O than Y-O bonds; and accordingly, significant bonding heterogeneity in Y 2 Si 2 O 7 . Besides, as hydrostatic pressure increases to P = 3.06 GPa, IFCs for both types of bonds are raised to higher values, which is a direct reflection of the enhanced overall atomic bonding strength due to volume contraction. In order to quantify the mismatch of force constant of atomic bonding in each crystal geometry, a parameter M IFCs is defined as (n = 10 in our statistical analysis) based on the average squared deviation of IFC for i th individual bonds (IFC i ) from their arithmetic mean value (IFC). Results are summarized in the inset of Fig. 4(b). Clearly evident, the value of M IFCs is gradually raised by 2.2% at P = 3.06 GPa. Therefore, the increased lattice anharmonicity in pressured crystal is a reasonable consequence of the enlarged mismatch of force constants between neighboring bonds.

Correlation between elastic properties and acoustic phonon softening.
It is well known that elastic moduli of a material are closely related with the long-wavelength part of the acoustic phonon branches near BZ center, and thus are critical parameters in tracing the physical properties of acoustic phonons. Calculation of elastic parameters reveals that under hydrostatic pressure, the overall strength of the structure (measured by bulk modulus B) is enhanced by Δ B/Δ P = 3.00, as expected from higher average atomic bonding strength upon lattice contraction; however, the resistance to shear deformation (measured by shear modulus G) shows a weak decrease  Supplementary Fig. S3 for details). This abnormity could be attributed to active accommodation of shear strain by the localized deformation of soft and flexible YO 6 polyhedral units, which has been reported as the typical shear deformation mechanism in yttrium silicates 41 , and again, is originated from the intrinsic bonding heterogeneity of Y 2 Si 2 O 7 . In fact, it is fundamental knowledge that average sound velocity (v m ) could be approximated as acoustic phonon group velocity by assuming linearity of phonon dispersion, written in the form of = 3/2 3/2 1/3 4 3 (ρ being density of crystal); and hence calculated Δ v m /Δ P < 0 dependency qualitatively echoes the negative pressure dependence of v ave , as we have covered before. Therefore, the weakly decreased (obviously increased) G (B) value is a mirror of the pressure-induced softening of acoustic phonons in Y 2 Si 2 O 7 , and has also been observed in other materials with the complicated interpenetrating of both soft and rigid (usually tetrahedrally coordinated) structural units 42 . Accordingly, one might be able to make a primary deduction on the acoustic phonon behavior, and furthermore, the positive or negative κ L (P) relationship of a novel material by monitoring the pressure dependence of elastic moduli, which could be easily obtained from either theoretical or experimental approaches.

Discussions
In this study, we predict an anomalous decrement (55.2%) of lattice thermal conductivity (κ L ) in Y 2 Si 2 O 7 under hydrostatic pressure up to P = 3.06 GPa. The mechanism of the dκ L /dP < 0 relationship is disclosed as the combination of slightly decreased group velocity and significantly augmented scattering of heat-carrying acoustic phonons. The coordinated vibration of flexible YO 6 and rigid [SiO 3 -O-SiO 3 ] structural units with highly diverse interatomic force constant renders the system a collective of low-lying optic phonons that strongly interact with acoustic phonons through acoustic-optic coupling; and the particularly severer interaction with TA than LA phonons makes the former (latter) negligible (dominant) contributors to thermal conduction. Under pressure, the anomalous downward shift of low-lying optic phonons and heat-carrying acoustic phonons tend to significantly raise acoustic phonon anharmonicity, and eventually block the phonon transport by slightly suppressing the acoustic phonon group velocity and facilitating three-phonon scattering process. Besides, upon tracing to experimentally-accessible elastic properties of Y 2 Si 2 O 7 , we find that the crystal undergoes an enhancement of the overall lattice strength under hydrostatic pressure; and at the same time maintains the capability to positively accommodate shear strain via localized structural deformation of flexible YO 6 octahedra. This unique elastic deformation mechanism, governed by the intrinsically significant bonding heterogeneity, leads to a positive (weak) pressure-dependence of bulk (shear) modulus. Based on above discussions, a guidance could be reached for the identification of solids in which loading hydrostatic pressure is an initiator to abnormally decreased lattice thermal conductivity. Such material should have diverging structural flexibility, enormous bonding heterogeneity, and giant phonon anharmonicity. And these criteria are predicted to cover numerous candidates of complex-structured oxide ceramics with interlaced stacking of rigid/soft polyhedral units, including various members of silicates, phosphates, zirconates and tungstates, to name a few.
From the perspective of phonon engineering, we would like to point out that the motivation of this paper is to open up the possibility of negatively tailoring lattice thermal conductivity of complex ceramics via hydrostatic pressure, based on modulation of intrinsic structural characteristics. And hence, it would be interesting to compare the phonon transport mechanisms in Y 2 Si 2 O 7 to those reported in HgTe having simple crystal structure 19 . One point in common is that both work emphasize the significance of augmented phonon scattering in determining the dκ L /dP < 0 relationship due to enhanced phonon anharmonicity. However, major discrepancies expose if one look into specific phonon behaviors of the two systems by taking into consideration the diverging bonding characteristics. Firstly, the highly mixing of low-lying optic phonons in acoustic regime of Y 2 Si 2 O 7 results in substantial scattering of TA phonons, leading to negligible contribution of TA phonons to thermal conduction as compared with LA phonons; In HgTe, however, the sizable acoustic-optic frequency gap seems to freeze out most (aao) process involving low-frequency TA phonons due to the restriction of phonon energy/momentum conservation; whereas LA phonons could decently participate in both (aaa) and (aao) scattering options, resulting in shortened relaxation time and attenuated contribution to thermal conduction. Consequently, it is the behaviors of LA (TA) phonons that dominate thermal transport process in Y 2 Si 2 O 7 (HgTe). Secondly, a wide range of low-frequency optic and acoustic phonons in Y 2 Si 2 O 7 undergo downward frequency shift upon hydrostatic pressure; while conventional phonon stiffening is observed in HgTe with minor exception for several TA phonon modes near BZ boundary. As a result, the decrement (increment) of phonon group velocity induced by phonon softening (stiffening) serves as a contributive (compensated) factor to the dκ L /dP < 0 relationship in Y 2 Si 2 O 7 (HgTe). In essence, the role of acoustic-optic coupling is accentuated in Y 2 Si 2 O 7 , contributed by a collective of low-lying optic phonons, i.e. coupled vibration of a set of atoms in hierarchical neighboring atomic environment, which runs over respectively wide range of phonon spectrum. And such phenomenon is inclined to be observed in ceramics with multiple chemical composition, highly complicated crystal structure and significant bonding heterogeneity, rather than in the simple crystals, such as HgTe (in cubic lattice with two non-equivalent atoms in primitive unit cell). Therefore, future work on various material systems is encouraged to elucidate the interplay of different scattering options involving low-lying optic and acoustic phonons as well as its behavior upon modulation of structural characteristics, so as to better guide the practicability of phonon engineering. phonons (κ L TA1 and κ L TA2 ); whereas the direct participation of optic phonons is neglected due to their low group velocity. And the reliability of such approximation in Y 2 Si 2 O 7 has been verified by previous work 30 , in which relaxation time of optic phonons are found to be lower than that of acoustic phonons by 3 ~ 4 orders of magnitude.

Intrinsic
The partial conductivities κ L i are given by: Here, i denotes phonon branches (LA, TA1 and TA2), k B is the Boltzmann constant,  is the reduced Planck constant, T is absolute temperature, v i is corresponding acoustic phonon group velocity, τ C is the total relaxation time of phonon scattering process, is a dimensionless quantity derived from phonon frequency ω; and θ i is the longitudinal and transverse acoustic Debye temperature obtained from phonon dispersion as , where ω i max is the maximum frequency of the acoustic phonon branch i in each BZ direction 39 . Here, phonon group velocity v i is calculated from the slope of acoustic phonon dispersion around Г point, averaged along each BZ direction.
The total phonon scattering rate (τ C In perfect lattice, such as pure Y 2 Si 2 O 7 in our study, the contribution from latter terms could be ignored since Umklapp phonon-phonon scattering effect dominates at above Debye temperature, and thus τ C ≈ τ U . Although many expressions exist for τ U , they all take the form as quadratic inverse of γ with similar pre-factors 1 . Based on the work of Slack and Galginaitis 43 , we employ the following expression for the Umklapp scattering rate for phonon branch i 27 : where M is the average mass of an atom in the crystal.
In the present calculation, mode Grüneisen constant is obtained by slightly perturbing the crystal away from the original volume V and calculate the resulted frequency change of each mode: where γ i (q) and ω i (q) are mode Grüneisen constant and phonon frequency, respectively, at wave vector q for phonon branch i. Accordingly, the average Grüneisen constant γ i for each acoustic branch could be written as the average of γ i (q) weighed by mode contribution of specific heat capacity 29 : Here, the absolute value of γ q ( ) i is emphasized to avoid the cancelation between modes with positive or negative values. All the parameters presented above are directly extracted from lattice dynamic calculations, allowing estimation of thermal conductivity from phonon behaviors.
As an alternative, κ L could also be obtained from anharmonic lattice dynamics calculation, by solving linearized phonon Boltzmann transport equation under single-mode relaxation-time approximation (SMRT-LBTE) and employing DFT-calculated anharmonic interatomic force constants as inputs. In this method, κ L takes the form of 44 : where V 0 is the volume of a unit cell; N is the number of unit cells in the crystal; λ is abbreviated for phonon mode (i, q); C λ , v λ and τ λ are the specific heat capacity, group velocity and relaxation time of mode λ, in regular definition as stressed before. Phonon relaxation time (τ λ ) is computed from imaginary part of phonon self-energy, by using many-body perturbation theory with the crystal potential expanded to include third-order terms. We should note that ceramics with low thermal conductivity (such as Y 2 Si 2 O 7 in our study) typically has complicated crystal structure, i.e. in low crystal symmetry and contains as much as tens of atoms in the unit cell, which makes anharmonic dynamics calculation computationally formidable and time-wise unaffordable. In this situation, RTA-based Debye-Callaway model has long been approved as a convenient methodology, considering its reliability in reproducing experimental data of various ceramic systems 29,45 .

Computation details. Density functional theory (DFT) calculation is performed using Vienna Ab Initio
Simulation Package (VASP) program 46,47 , where localized density approximation (LDA) 48 and projected augmented wave (PAW) method 49 are employed, with a cutoff energy of 600 eV. The valence shells and electronic configurations for pseudo-atoms are 4s 2 4p 6 4d 1 5s 2 for Y, 3s 2 3p 2 for Si and 2s 2 2p 4 for O. The special points sampling integration over Brillouin zone (BZ) is realized by using Γ -centered 5 × 3 × 5 Monkhorst-Pack 50 k-points method. Lattice parameters and internal atomic positions are optimized until the total energy and force converge to 10 −9 eV and 10 −4 eV/Å, respectively. Calculated lattice parameters at equilibrium structure of Y 2 Si 2 O 7 are a = 4.667 Å, b = 10.730 Å, c = 5.510 Å, β = 95.84°, which agrees well with experimental data 51 and previous calculations 52 . The effect of hydrostatic pressure is simulated by imposing a volume contraction on the equilibrium geometry of Y 2 Si 2 O 7 up to 2%, with an interval of 0.5%. The unit-cell vectors and atomic positions are fully relaxed via geometry optimization under each volume-fixed lattice until convergence criteria is satisfied; and the structural parameters including lattice constants, crystal symmetry, and atomic configuration are carefully monitored to ensure structural stability upon volume modulation. Lattice-dynamics calculations are carried out using the Phonopy package 53 . Second-order interatomic force constants (IFC) are calculated within the framework of density function perturbation theory (DFPT) via linear-response approach 54 . Phonon frequencies and eigenvectors are obtained from diagonalization of the dynamical matrix built by IFCs. During post-process, phonon frequencies are sampled on a Γ -centered 3 × 3 × 3 q mesh, which could converge the thermal dynamic properties to a sufficient extent.
SMRT-LBTE calculation is conducted using Phono3py 31 , with harmonic and anharmonic interatomic force constant extracted from DFT calculations. Calculation parameters are chosen identical as above, although 2 ndand 3 rd -order IFCs are calculated from supercell approach (2 × 1 × 2 supercell containing 88 atoms) with finite atomic displacement of 0.03 Å; and translational invariance constant is enforced when dealing with 3 rd IFC. Atomic interaction cut-off of 3 rd IFC is chosen as sixth-and fourth-nearest neighbors for Y and Si atoms, respectively; whereas increasing the cut-off to include twelfth-/eighth-nearest neighbors for Y/Si renders a negligible correction (0.5%) to κ L . Thermal conductivity is calculated with q-mesh of 3 × 3 × 3; and a tetrahedron method 55 is employed in calculating imaginary part of self-energy and thus phonon lifetime.
Calculation of elastic parameters. Second-order elastic coefficients (c ij ) are determined from the linear fit of stress as a function of applied homogeneous elastic strains 56 . We apply four homogeneous strain patterns with finite values and calculate the generated stress after optimizing the internal degrees of freedom. For each strain pattern, six amplitudes ε (three positive and three negative) are applied with |ε | ≦ 0.5%. The polycrystalline bulk modulus B and shear modulus G are determined by Voigt-Reuss-Hill (VRH) approximations [57][58][59] , where values calculated from second-order stiffness constants (c ij in C) and second-order compliance constants (s ij in S = C −1 ) are averaged to represent the real elastic moduli: