Origins of thermal conductivity changes in strained crystals

Kevin D. Parrish,1 Ankit Jain,1 Jason M. Larkin,1 Wissam A. Saidi,2 and Alan J. H. McGaughey1,* 1Department of Mechanical Engineering, Carnegie Mellon University Pittsburgh, Pennsylvania 15213, USA 2Department of Mechanical Engineering and Materials Science, University of Pittsburgh, Pittsburgh, Pennsylvania 15213, USA (Received 19 June 2014; revised manuscript received 23 October 2014; published 5 December 2014)


I. INTRODUCTION
To simplify modeling efforts, the assumption that a material or device is at zero strain is often made.Under realistic conditions in many systems, however, zero strain is not the case.For example, the extreme pressures in the Earth's core (∼300 GPa) can increase the thermal conductivity of the constituent materials (e.g., iron) by a factor of three [1].Strain can also be induced by epitaxial interfaces [2][3][4], mechanical loading, and the thermal expansion that results from an imposed temperature gradient or Joule heating.Strain can be used to engineer the transport properties of electrons and phonons [5], to tune the thermal conductance of interfaces [6,7], and can modify the thermal conductivities of nanostructures such as nanowires and one-and twodimensional materials [8][9][10][11][12].
There have been previous efforts to quantify the effect of strain on thermal conductivity in bulk crystalline systems.Ross et al. presented a comprehensive review of experimental measurements of the thermal conductivity of strained materials, including covalent and semiconducting structures [13].They concluded that compression increases the thermal conductivity of this class of materials.Picu et al. studied strain effects on a model Lennard-Jones (LJ) crystal using molecular dynamics (MD) simulations for isotropic strains ranging from −0.03 to 0.03, where a negative strain indicates compression [14].Using the Green-Kubo (GK) MD-based method, they found that the thermal conductivity increased under compression and decreased under tension.They attributed this behavior to changes in "stiffness."Tretiakov and Scandolo explored thermal transport in crystalline argon at high temperatures and pressures using MD simulations, the GK method, and the exponential-6 potential [15].As the compressive stress increased from 2 to 50 GPa, they found a power-law increase of thermal conductivity.Bhowmick and Shenoy studied strained, crystalline LJ argon using MD simulations and the GK method, also finding a power-law scaling of thermal conductivity with * mcgaughey@cmu.eduincreasing compressive strain [16].They predicted the same power-law scaling from Fermi's golden rule (exponent within 5%) based on analytical scalings of the phonon lifetime and sound speed.Goncharov et al. experimentally investigated thermal transport in high-temperature, compressed argon using a transient heating technique [17].Using a diamond anvil cell, argon was compressed in the stress range of 10 to 50 GPa.They found a power-law dependence of thermal conductivity with increasing compression at a temperature of 300 K. Chernatynskiy and Phillpot explored high-temperature (400 to 1200 K) argon using harmonic and anharmonic lattice dynamics calculations with force constants obtained from density functional theory (DFT) calculations over a compressive stress range of 10 to 50 GPa [18].Their thermal conductivity predictions fit the power-law functional form proposed by Tretiakov and Scandolo [15].Li et al. studied thermal transport in strained crystalline Tersoff silicon and diamond using MD simulations and the GK method [8].They found thermal conductivity to monotonically decrease from compression to tension over a strain range of −0.09 to 0.12 and justified this trend using changes in the phonon group velocities.While these previous studies focused on quantifying how thermal conductivity is altered by strain, modifications to the properties of individual phonon modes due to strain have yet to be explored.
In this work, we apply atomistic calculations to investigate how isotropic strain affects phonon mode properties in bulk systems where electrons make a negligible contribution to thermal conductivity.From a solution to the Boltzmann transport equation, the lattice (i.e., phonon) thermal conductivity in direction n, k n , can be expressed as [19] Here, the summation is over all phonon modes with wave vector κ κ κ and polarization branch ν.Inside the summation, C ph ( κ κ κ ν ) is the heat capacity, v g,n ( κ κ κ ν ) is the component of the group velocity vector v g in the direction n, and τ ( κ κ κ ν ) is the lifetime.We will predict thermal conductivity in the [100] direction, although all directions are equivalent through the cubic symmetry of the materials considered.
We first consider an LJ argon crystal as a representative soft system.Thermal transport in LJ argon is well-studied [14,16,[20][21][22][23][24], primarily because it uses a computationally-inexpensive two-body potential.As typical of a soft material, LJ argon has a low zero-strain thermal conductivity [O(1) W/m K] that results from contributions from the full spectrum of phonon modes [25].Molecular dynamics simulations and the GK method are first used to predict thermal conductivity for benchmarking purposes.Phonon properties are then predicted from harmonic lattice dynamics (HLD) [26,27] and anharmonic lattice dynamics (ALD) [22,[28][29][30] calculations and are used to evaluate Eq. (1).
We then consider silicon as a representative stiff material.Stiff materials typically exhibit a high thermal conductivity that is dominated by low-frequency acoustic phonon modes [25,31].The phonon properties required to evaluate Eq. ( 1) are predicted from HLD and ALD calculations with force constants obtained from DFT and density functional perturbation theory (DFPT) calculations.
The rest of the paper is organized as follows.In Sec.II, we describe the computational methods and details.The thermal conductivities and phonon mode properties are then presented and discussed for LJ argon (Sec.III A) and silicon (Sec.III B).To explain differences in the phonon lifetime trends for these two materials, the potential energy wells experienced by the atoms in each are compared in Sec.III C. We conclude by summarizing our results in Sec.IV.

A. Lennard-Jones Argon
Argon is an insulating face-centered cubic crystal that can be described by the two-body LJ interatomic potential [32] where i and j denote a pair of atoms, r ij is the distance between the two atoms, LJ is the energy scale (1.67×10 −21 J), and σ LJ is the length scale (3.40×10 −10 m).Argon has an atomic mass m, of 6.63×10 −26 kg.All calculations are made at a temperature of 20 K, where the zero pressure lattice constant, a • , is 5.315 Å [23].We use a cutoff of 2.5σ LJ with a shifted potential well for the MD simulations and the lattice dynamics (LD) calculations.Isotropic strains, ε, of −0.06 to 0.06 are considered, such that the lattice constant, a, is a • (1 + ε).The strains, lattice constants, and resulting stresses are provided in Table I.
The primitive one-atom basis is used to perform the HLD and ALD calculations with in-house codes.The phonon wave-vector grid is sampled using N • uniformly spaced wave vectors in each direction.Ten values of N • are considered ranging from 32 to 50 in increments of two.The phonon frequencies ω are obtained from HLD calculations, which require the second-order force constants [27].Using the frequencies, we calculate the group velocities using finite differencing from v g,n = dω/dκ n .The second-order force constants are calculated up to the potential cutoff of 2.5σ LJ .Since harmonic plane waves do not interact, HLD calculations do not provide any information about the phonon scattering and lifetimes.The phonon lifetimes are predicted using ALD calculations, which use the third-order force constants to account for three-phonon interactions [22,33].The third-order force constants are calculated up to first-nearest neighbor.The force constants are calculated by displacing selected atoms by 10 −4 a from their equilibrium position in a supercell and finite differencing of the resulting forces.
Using the group velocities and lifetimes, we evaluate thermal conductivity from Eq. (1).Because the MD simulations to be used for benchmarking are classical, we use classical occupation statistics in the ALD calculations and the classical harmonic value of heat capacity where k B is the Boltzmann constant and V is the system volume a 3 N 3 • .The inverse of thermal conductivity is plotted verses 1/N • and a line is fit to the data, whose extrapolated value at 1/N • = 0 gives the infinite system-size (i.e., bulk) value [22].
Turney et al. observed a difference between GK and ALD predictions of thermal conductivity for LJ argon at zero strain for temperatures above 40 K, where increasing anharmonicity cannot be captured with ALD.Because strain affects the anharmonicity, we also use the MD-based GK method as a self-consistent check of the thermal conductivities calculated from ALD.
The MD simulations are performed using the open-source LAMMPS package [34].The time step is 4.28 fs and the cubic simulation cell contains 256 atoms.We did not observe changes in thermal conductivity due to size effects over all strains.At each strain, a temperature rescale is performed for 250,000 time steps to bring the system temperature to 20 K.The system is then evolved in an NV E ensemble (a constant number of atoms, volume, and energy) for 250,000 time steps.A further 1,000,000 time steps are run where the heat current vector is recorded every five time steps.Each strain case is run with ten seeds (randomized initial velocities) and the heat current autocorrelation functions are averaged.The heat current autocorrelation function is the input to the GK method, which is an equilibrium technique for extracting the  thermal conductivity k GK based on the fluctuation-dissipation theorem [35].There is inevitably noise in the heat current autocorrelation function due to a limited sampling of the phase space from the finite sampling time.To find the converged value of the integral of the heat current autocorrelation function, we applied the first-avalanche method [36].We found k GK to be converged to within 4% with 10 seeds.This error was estimated by removing one seed and calculating the difference in the predicted thermal conductivity when all ten seeds are used.

B. Silicon
We consider isotopically pure silicon at a temperature of 300 K.At this temperature, silicon undergoes a phase change from the diamond lattice structure to the β-Sn lattice structure at compressive stresses greater than 12 GPa [37].To avoid this phase transition, we explored strains between −0.03 and 0.03, as provided in Table II.HLD and ALD calculations using in-house codes are performed to obtain phonon frequencies and lifetimes with force constants obtained from DFPT and DFT calculations [29].A norm-conserving pseudo-potential in the local density approximation is employed in the plane-wave package QUANTUM ESPRESSO [38].The total energy is converged to within 1.6 meV for a Monkhorst-Pack electronic wave-vector grid of 6×6×6 and a plane-wave energy cutoff of 80 Ryd.The harmonic force constants are obtained using DFPT calculations with a phonon wave-vector grid of 8×8×8 using the primitive (i.e., two-atom) unit cell.The third-order force constants are obtained using finite differences of DFT forces on a 64-atom supercell with a cutoff of third-nearest neighbor.The forces are obtained by applying displacements of 0.006a to selected atoms.The translational invariance constraint (i.e., the acoustic sum rule) for the second-and third-order force constants is enforced using the Lagrangian approach discussed by Lindsay et al. [39].For the HLD and ALD calculations, the primitive basis and a 20×20×20 sampling of the phonon wave-vector grid are used.
We are studying silicon at a temperature of 300 K, which is less than half of its Debye temperature of 625 K [32].Quantum effects on the heat capacity are important and are included using Bose-Einstein statistics.The thermal conductivity of unstrained silicon calculated using lifetimes from an iterative solution to the Boltzmann transport equation [30] is 151 W/m K and is converged to within 2% for the above choice of parameters.The experimental value for isotopically enriched silicon is 158 W/m K [40].This 5% difference is a large improvement over empirical potential-based predictions of thermal conductivity, which can differ from the experimental value by as much as a factor of six at a temperature of 300 K [41].All reported thermal conductivities are from the iterative solution to the Boltzmann transport equation.For the analysis to follow in Sec.III B, the lifetimes obtained from the relaxation time approximation will be used because the iterative lifetimes have larger noise that makes direct comparison across strains unfeasible.The relaxation time approximation lifetimes lead to a thermal conductivity prediction of 147 W/m K, such that we do not believe that this choice will affect our conclusions.

Thermal conductivity
The strain-dependence of thermal conductivity for LJ argon predicted from MD simulations and LD calculations is presented in Table I and is plotted in Fig. 1.The thermal conductivities span two orders of magnitude.Because the MDbased GK method includes all anharmonic effects (compared to ALD, where we only include up to third-order phonon scattering effects), k GK is taken as the benchmark.The two sets of thermal conductivity data follow the same trend: an exponential decrease when moving from compression to tension (as evidenced by a straight line on the semilog plot).This trend for LJ argon at 20 K was also observed by Bhowmick and Shenoy using the GK method and by Fermi's golden rule-based analytical scalings [16].The difference between k LD and k GK is less than 25% for all strains considered.

Group velocities and lifetimes
We now seek to understand the thermal conductivity trend in Fig. 1 by considering the phonon mode properties.The changes in group velocity can be visualized by considering the [100] phonon dispersion plotted in Fig. 2 for strains of −0.06, 0, and 0.06.The maximum phonon frequency decreases as the system moves from compression to tension, causing the group velocities to decrease.Since the dispersion is calculated using harmonic force constants, it only captures the harmonic effect of strain.The phonon lifetimes are plotted in Fig. 3.As with the group velocities, the lifetimes decrease as the system moves from compression to tension.
We quantify the mode-level changes by dividing the strainand mode-dependent squared group velocities (|v g | 2 ) and lifetimes by their zero-strain counterparts.The results for strains of −0.06 and 0.06 are plotted in Figs.4(a [Fig.4(a)] and lifetimes [Fig.4(b)] are relatively flat, with effectively no spectral dependence.The relative change in the squared group velocities with strain is larger than that of the lifetimes.
To get a system-level measure of the mode scalings, we average the normalized properties at each strain.The results are plotted in Fig. 5.The normalized squared group velocities and lifetimes both scale exponentially and the least-squares fits are plotted in Fig. 5 [42].The exponential constant is −17 for the squared group velocity and −9 for the lifetime.Therefore the group velocity scaling is exp (−8 ε) stronger than that of the lifetime, making it the dominant factor in thermal conductivity changes with strain.

Thermal conductivity
The strain-dependence of the thermal conductivity of silicon predicted from first-principles is plotted in Fig. 6.Thermal conductivity is relatively constant under compression and then decreases with increasing tension.This trend is in contrast to that for LJ argon, where the thermal conductivity monotonically decreased in moving from compression into tension (Fig. 1).The overall change in thermal conductivity is smaller in silicon (∼15%) as compared to LJ argon (two orders of magnitude).

Heat capacity, group velocities, and lifetimes
In interpreting the silicon data at the phonon mode level, we consider changes in heat capacity, squared group velocity, and lifetime.In Fig. 7, we plot the silicon dispersion.The maximum frequency decreases as the system moves from compression into tension, but the difference is less than that observed in LJ argon (Fig. 2).The stretching of the dispersion results in a decrease in group velocity under tension with the exception of the transverse acoustic branch, which reduces in frequency near the Brillouin zone edge.The silicon lifetimes plotted in Fig. 8 display an increase of the lifetimes in tension, opposite to the trend for LJ argon (Fig. 3).
The mode-dependent property normalizations to the zero strain values are plotted in Figs.9(a)-9(c).Because the thermal conductivity of silicon is dominated by low-frequency acoustic phonons [25,29], we do not include modes with frequencies past the point where 90% of the thermal conductivity is realized.This choice eliminates modes that do not meaningfully contribute to the thermal conductivity, but can show large scatter in the normalizations.The changes in mode properties with strain are smaller than the changes observed in LJ argon, which is consistent with the smaller changes in thermal conductivity.The heat capacity follows a predictable scaling from Bose-Einstein statistics.As the phonon mode's frequency (i.e., energy) decreases in moving from compression to tension, its heat capacity decreases.No trend in the normalized squared group velocities is discernible from Fig. 9

(b).
There is an overall reduction of lifetimes with compressive strain and an overall increase in lifetimes with tensile strain.The squared group velocities and lifetimes show more spectral dependence than in LJ argon.
In Fig. 10, we plot the mean values of the normalized properties versus strain, as plotted in Fig. 5 for LJ argon.Before performing the average, we exclude any group velocity normalizations that are greater than a value of five.In doing so, we eliminate nonphysical values due to the effects of van Hove singularities (i.e., where the group velocity approaches zero).The phonon frequencies decrease as the system moves from compression to tension, and thereby the heat capacity  and squared group velocity also decrease.These changes are related to harmonic effects.There is an increase in the lifetimes when moving from compression to tension over the entire strain range explored, opposite to what is observed in LJ argon.These competing effects balance in compression, leading to a constant thermal conductivity.In tension, the heat capacity and squared group velocity effects dominate, leading to a decrease in thermal conductivity.Using MD simulations, Li et al. found a monotonic decrease in the thermal conductivity of Tersoff silicon in moving from compression to tension.They attribute this behavior to a reduction in the group velocities and heat capacities [8].We note that they explored a larger strain range (−0.09 to 0.12) and used an empirical potential with the classical statistics inherent to MD simulation.

Root mean squared displacement
The strain-dependence of the phonon lifetimes is dissimilar in LJ argon and silicon.In LJ argon, the lifetimes decrease as the system moves from compression to tension (Figs. 3 and 5), while in silicon they increase (Figs. 8 and 10).To explore these opposite trends, we analyze the local potential well felt by atoms in the two systems.First, the root mean-square (RMS) displacement is calculated to quantify how much of the unit cell is explored by an atom.Second, the local potential energy well is mapped to decompose the potential energy into harmonic and anharmonic components.
The RMS displacement |u| 2 1/2 can be calculated through MD simulations by time-averaging the positions or from the harmonic phonon properties using the analytical expression [43] Here, u is the displacement of an atom from its equilibrium position, is the reduced Planck constant, n is the number of atoms in the unit cell, e( κ κ κ b ν α ) is the component of the eigenvector for mode (κ κ κ,ν) for atom b in direction α, e * ( κ κ κ b ν α ) is the complex conjugate eigenvector component, and f ( κ κ κ ν ) is the distribution function.
For the LJ argon crystal, the RMS displacement is calculated directly using the MD simulations and using Eq.(3).We use classical-limit value of f ( κ κ κ ν ) = k B T / ω( κ κ κ ν ) in Eq. ( 3) so as to make the prediction comparable to the MD result.The value is normalized by the nearest-neighbor distance, a/ √ 2, for each strain and the results are plotted in Fig. 11(a).The RMS displacement for LJ argon from MD is greater than that calculated from Eq. ( 3), which we attribute to the inclusion of the full anharmonicity in the MD simulations.We note that Turney et al. [22] erroneously included the zero-point energy in the phonon occupation in Eq. ( 3) yielding accidental agreement with MD simulation results.As it is computationally prohibitive to run DFT-driven MD simulations, the RMS displacement for silicon is calculated only using Eq.(3).The results, normalized by the nearest-neighbor distance, √ 3a/4, are plotted in Fig. 11(b).We use the Bose-Einstein distribution for the distribution function to be consistent with the thermal conductivity prediction.
The normalized RMS displacement for LJ argon more than doubles from compression to tension.The percentage of the unit cell that is accessible by the atom is reduced in compression and increased in tension.In silicon, the opposite trend is observed: the normalized RMS displacement decreases from compression to tension.Since the atom is only allowed to displace from equilibrium through the natural energy of the thermal fluctuations, k B T , the further the atom deviates from the equilibrium position, the more phonon-phonon scattering occurs (i.e., the system is more anharmonic).The RMS data for LJ argon and silicon are therefore consistent with the lifetime data in Figs. 3, 5, 8, and 10.

Potential snergy surface mapping
To further understand the difference in the lifetime trends, we next performed molecular statics (i.e., single-point energy) calculations for LJ argon and silicon to map the potential energy surface experienced by an atom.A single atom is displaced in small increments along the [100] direction and the total potential energy of the system is calculated.This analysis is performed up to the RMS displacement.
In Figs.12(a) and 12(b), the displacement energy (i.e., the energy minus the zero-displacement energy) is plotted for LJ argon and silicon.To capture the harmonic part of the curve, a parabola is fit to the first four displacements and is plotted over the full range.The curvature of the harmonic potential decreases in both argon and silicon as the strain is increased.For both systems, the harmonic curve deviates from the full potential well as the atom is displaced further from its equilibrium position, an indication of increasing anharmonicity.The deviation from the harmonic curve is smaller in silicon as compared to argon, which is expected because silicon is stiffer.
We next subtract the harmonic fits from the potential wells and the results are plotted in Figs.12(c) and 12(d).For LJ argon [Fig.12(c)], the anharmonic contribution is positive for all strains.The anharmonicity at the RMS displacement increases as the strain is increased, thereby increasing phonon-phonon scattering and reducing the lifetimes (Fig. 5).In silicon [Fig.12(d)], the anharmonic contribution is negative for all strains and the greatest amount of anharmonicity (i.e., deviation from the harmonic potential) is felt under compression.This decreasing anharmonicity from compression to tension increases the lifetimes (Fig. 10).
In silicon, we observe competing effects between harmonic (i.e., heat capacity and group velocity) and anharmonic (i.e., lifetime) effects.These two effects are thus decoupled and need not scale together under strain, as seen in LJ argon.Compression does not just increase the "stiffness" [14], but can induce positive or negative changes in the anharmonic contribution to the energy potential.The opposite behaviors of the LJ argon and silicon systems demonstrate the need to separately inspect harmonic and anharmonic effects when interpreting strain-dependent properties.

IV. SUMMARY
We applied atomistic calculations to predict the straindependent phonon properties and thermal conductivities of LJ argon and silicon.LJ argon undergoes an exponential decrease in thermal conductivity with increasing strain (Fig. 1).For silicon, as shown in Fig. 6, the thermal conductivity remains constant under compressive strain and decreases with increasing tensile strain.
For LJ argon, the mode-averaged lifetimes and squared group velocities both decrease exponentially as the system moves from compression to tension, with the squared group velocity effect dominating the thermal conductivity (Fig. 5).For silicon, the squared group velocities decrease from compression to tension, while the lifetimes increase anomalously, as shown in Fig. 10.
To explain this behavior, we examined the local potential well.The normalized RMS displacement, plotted in Figs.11(a increased.These trends are consistent with the lifetime trends, which is physically justified as larger displacements will lead to more anharmonicity and reduced lifetimes. We performed molecular statics to map the local potential well in the [100] direction.As shown in Figs.12(a) and 12(b), for LJ argon, the anharmonic contribution increases with increasing strain, similar to how the harmonic contribution increases.For silicon, however, the anharmonic contribution decreases with increasing strain, opposite to the argon trend [Figs.12(c) and 12(d)].Strain can thus affect the anharmonic contribution to the potential differently than it affects the harmonic contribution, qualitatively changing the scaling of the phonon lifetimes.
FIG. 5. (Color online) Normalized mean lifetimes and mean squared group velocities for LJ argon for N • = 50.Each data set is fit to an exponential function (solid lines).

FIG. 12 .
FIG. 12. (Color online) Potential energy changes for displacement in the [100] direction for multiple strains.(a) LJ argon displacement energy and harmonic fit.(b) Silicon displacement energy and harmonic fit.(c) LJ argon anharmonic energy.(d) Silicon anharmonic energy.