Thermal conductivity prediction for GaN nanowires from atomistic potential

A model is developed to evaluate the thermal conductivity of semiconducting compounds as a function of their atomistic structures during phase transformations induced by mechanical loading. The approach uses atomistic configurational information and interatomic interactions as input. The harmonic and anharmonic behaviors of phonons are captured through force constants which are sensitive to structural changes. The calculations focus on changes in thermal conductivity of GaN nanowires in response to deformation and phase transformation. Results show that the model yields results consistent with data obtained using the Green-Kubo method and is 50 times more efficient than calculations based on molecular dynamics.


I. INTRODUCTION
For many nanostructures, external loading can lead to structural transitions, resulting in changes in mechanical, thermal, electrical, and optical properties.3][4] In particular, such structural changes can have significant impact on thermal transport which is an important consideration in the design of nanodevices that use nanomaterials.Molecular dynamics (MD) simulations are widely used to study the mechanical behavior of materials, including structural changes under loading.8][19][20][21] Among these methods, the MD method is often used for the calculation of thermal conductivity in nanowires because relatively large numbers of atoms (>5000) can be considered to capture both structural changes and associated thermal response changes simultaneously.Models based on the BTE are not usually used because they do not distinctly consider the effects of mechanical strain and phase transformations.
In our previous work, 22,23 the thermomechanical behaviors of [0001]-oriented GaN nanowires are investigated using MD simulations.Phase transformations are observed during tensile loading and unloading.Specifically as shown in Fig. 1 using the equilibrium Green-Kubo method.The thermal conductivity of the nanowires is found to be dependent on the strain and the crystalline structure state during the tensile loading and unloading.
In order to evaluate the thermal conductivity using the Green-Kubo approach, the heat flux should be calculated over periods of at least 500 ps.This usually corresponds to one million time steps in MD simulations. 24,25 n contrast, to resolve a mechanical loading step, 20,000 time steps are sufficient to reach equilibrium in stress, temperature, and energy and capture the corresponding mechanical response.For Cu nanowires 1 and ZnO nanowires, 3 which show a lattice reorientation and the phase transformation respectively under tensile loading, the equilibrium times for each loading increment are less than 20,000 time steps.To put it in perspective, the computational cost for the calculation of thermal conductivity is fifty times that for the mechanical response.This imbalance makes coupled thermal-mechanical simulations quite inefficient.
Bhowmick and Shenoy 26 presented a theoretical model for the dependence of thermal conductivity on mechanical strain.This model allows the thermal conductivity of solids to be calculated using empirical atomistic potentials.This model does not account for the effect of phase transformations because the model is based on the effect of strain only.Thermal conductivity predicted by this model simply decreases with mechanical strain.This trend is similar to the result for GaN nanowires reported by us previously. 23However, the dependence of thermal conductivity on the crystalline structures cannot be described properly.This model fails to explain the thermomechanical behaviors of [0001]-oriented GaN nanowires in Ref. 23 because thermal conductivities are different at the same strain in the loading and unloading processes due to the difference in crystal structures.This model is valid only for materials with constant structures.Here, we propose a new approach to allow the dependence of thermal conductivity on both strain and crystal structure of deformed nanowires to be accounted for.The focus is on thermal response changes resulting from both mechanical deformation and structural transformation.The account of the effect of structural changes on thermal conductivity is through the use of deformed atomistic configurations obtained from separate MD calculations at different stages of deformation.Force constants are calculated first using the results of MD simulations of mechanical response.We then compute changes in thermal conductivity in response to deformation and phase transformation.

II. ANALYTICAL APPROACH
In the Green-Kubo approach, the thermal conductivity is evaluated using the ensemble average of the autocorrelation of heat flux.The calculation of heat flux concerns the positions, velocities, and forces of atoms.In MD simulations, such atomistic quantities are determined by the atomistic potential.This means that the atomistic configuration and potential are essential factors for the evaluation of thermal conductivity.As illustrated in Fig. 2, the new model uses atomistic configurational information and interatomic potential to calculate the thermal conductivity as a function of material structure.The atomistic configuration is obtained from MD simulations used to calculate the mechanical properties.
The potential energy of the material can be expanded as a Taylor series in powers of atomic displacement u relative to the current macroscopic equilibrium state such that where 0 is the potential energy associated with the macroscopic equilibrium state at each stage of deformation; the i, j, and k sums are over all the atoms in the system; the α, β, and γ sums are over the x-, y-, and z-directions, respectively; and u i,α is the α-th component of the displacement vector of atom i.In the harmonic approximation, 26 the first derivative of the potential energy with respect to atomic displacement is zero because the energy is the lowest at equilibrium.The third-and higher order terms are truncated.Following the approximation, the above expression is rewritten as Under this approximation, the equation of motion can be written as where m i is the mass of i-th atom.The displacement of atom i can be written as a superposition of the normal modes of the system, i.e., 27 where k is the wave vector, v is the mode for the wave vector, e (k, v) is the mode eigen-vector, r i is the position of atom i, ω (k, v) is the angular frequency, and t is time.Normal mode coordinate Q is obtained by using the inverse Fourier transformation of the displacement as 27 where N is number of atoms in the system.The normal mode system is equivalent to the system of harmonic oscillators.The frequency of the system ω is obtained as In this work, the potential energy can be regarded as a summation of the interatomic energy φ r i j between atom i and j, i.e., because the atomistic potentials considered in this work are pairwise.From ( 6) and ( 7), the frequency is dependent on the second-order force constant V 2 such that The harmonic scattering rate of phonons due to mass disorder is dependent on the frequency. 15,28  reflect the effect of structural deformation, V 2 is assumed to be dominant for changes in the harmonic scattering rate and the frequency of phonons.The scattering rate is dependent on where g is a measure of the mass disorder and D (ω) is the phonon density of states. 15,28 n the above relation, τ harmonic can be reduced to be a function of V 2 because τ harmonic is a function of the phonon frequency.The anharmonicity of phonon behavior mostly originates from the third-order term in (1).The anharmonic scattering rate is computed by considering the three-phonon scattering processes in the single-mode relaxation time approximation 11,15,29 via 1 τ anharmonic = π q , j , j V j, −q; j , q ; j , q 2 × [(1 + n j q + n j q )δ(ω j q + ω j q − ω j,−q ) + 2(n j q − n j q )δ(ω j q − ω j q − ω j,−q )], (10)   where V j, −q; j , q ; j , q are the three-phonon coupling matrix elements, 29 n is the Bose-Einstein phonon population, q is the phonon momentum, and j is the branch index.The formula for the anharmonic scattering rate contains terms of phonon frequencies and three-phonon coupling matrix elements, which are related to the third-order interatomic force constant. 14,16,29 Fr systems described by pairwise potentials, the third-order force constant V 3 is obtained from the third-order derivative of the potential energy in a way similar to the manner in which V 2 is obtained, i.e., The anharmonic scattering rate is dependent on V 2 and V 3 through From the kinetic theory of heat carrier, the thermal conductivity can be evaluated using the phonon frequency and scattering rates.The phonon frequency ω is dependent on V 2 as shown in (8), and the scattering rates are functions of V 2 and V 3 in ( 9) and (12).Since V 2 and V 3 are considered to be dominant for changes in phonon behavior, the thermal conductivity can be expressed as a function of V 2 and V 3 as where κ 0 is the thermal conductivity of the undeformed structure.Function F represents normalized thermal conductivity as a function of strain.At zero strain, the value of F is 1, so that κ = κ 0 .The function F reflects the effects of second-order and third-order force constants.A simple form of F is chosen as where a and b are parameters and V 2 (ε) and V 3 (ε) are second-order and third-order force constants at strain ε, respectively.This equation is not obtained through a derivation.However, advantages of using such a simple expression are convenience and efficiency in the prediction of thermal response.The relation between the thermal response and function F is similar to the relation between the physical properties and the empirical potential.For example, the empirical potential is obtained by fitting against ab initio calculations or experimental properties.The MD calculation using empirical potentials is an approximation of the quantum mechanical calculation.In spite of this limitation, empirical potentials are widely used for prediction of physical properties because of its efficiency.Similarly, this expression of thermal conductivity in Eqn. ( 14) is useful in the prediction of the structure dependence of thermal conductivity, as described in the next section.In this paper, the parameters are chosen such that a = −1 and b = 2 from the relation between phonon group velocity and phonon relaxation time. 26The strain dependence of the thermal conductivity is obtained by following the Peierls-Boltzmann formulation. 26Each parameter contains the effects of phonon relaxation time and phonon group velocity.The anharmonic scattering is dependent on both V 2 and V 3 in Eqns.( 10) and (12).Using this model, we can predict the rate of change of thermal conductivity of deformed nanowires.The values for parameters a and b are the same regardless of wire diameter, strain, and temperature.
For the model of Bhowmick and Shenoy, 26 the thermal conductivity is expected to scale with strain as κ ∼ ε −μ , where μ is a potential dependent parameter.Because the phonon frequency and the third-order force constant are considered to scale with temperature and strain only, the effect of phase transformations is not considered.This model does not use the MD simulation as input, rather, data obtained from MD simulations are fitted to the function of κ = AT −λ ε −μ , where A, λ, and μ are constants.As shown in Fig. 4 of Ref. 26, the thermal conductivity decreases under tensile stress and increases under compression monotonically.This model quantifies the effect of strain on thermal conductivity.However, it does not have predictive power because the constants are obtained by fitting MD results to the function.To calculate the thermal conductivity at each strain, the model requires MD simulations using the equilibrium Green-Kubo method or the non-equilibrium direct method.For this reason, this model is not appropriate for prediction of thermal responses of the structures that undergo phase transformations in response to mechanical loading.
In the model presented here, the effects of strain and structural change are included in the prediction of thermal conductivity.The calculation of the second-order and third-order force constants uses the configurational information obtained from MD simulations.This new model accounts for the effects of deformation (due to both elastic straining and structure changes) on thermal conductivity at a fundamental, atomistic level.Specifically, the differences in thermal conductivity between different crystal structures 23 are captured because the atomistic configuration from MD simulations contains information from both elastic straining and crystal structure.The dependence of the thermal conductivity on structural changes results from the dependence of phonon behavior on material structure and is not described in the previous model which is based on strain. 26Since the thermal conductivity of unstressed nanowires has been investigated extensively, 20,21,[30][31][32][33][34] the focus here is on the thermal conductivity of deformed nanowires.The analysis combines the use of MD simulations of mechanical response to applied loading and the use of the model for the thermal evaluation, with the former providing input for the later in the form of structural quantification of the deformed wires.The model predicts the thermal conductivity of deformed nanowires by using the thermal conductivity of the unstressed nanowire and MD simulations of mechanical response.This approach requires only 2% of the computational cost of the classical Green-Kubo calculations.

III. RESULTS
The thermal conductivity of [0001]-oriented GaN nanowires is evaluated using the model.The atomistic configurations of deformed nanowires come from MD simulations of the mechanical responses in our previous work in Ref. 23.The MD simulations in the previous study are carried out at T = 300 K using the LAMMPS code.The nanowires are oriented in the [0001] direction with hexagonal cross-sections as shown in Fig. 1 of Ref. 23.The Buckingham potential 35 and the Wolf summation 36,37 are used to describe the interactions between atoms.The mechanical responses are analyzed using a quasi-static loading increment of 20,000 time steps with an NVT ensemble (canonical ensemble).To evaluate the thermal conductivity using the Green-Kubo method, the heat flux of each deformed nanowire is calculated over one million time steps with an NVE ensemble (microcanonical ensemble).Details of the MD simulations are described in Ref. 23.The computational cost for the mechanical response is only 2% of that for the Green-Kubo calculation of thermal conductivity.For each deformed configuration, V 2 and V 3 are calculated using the same interatomic potential as used in the previous study.The force constants associated with the thermal conductivity in the axial direction are calculated by summing the second and third-derivatives of the potential (Buckingham and columbic) over all possible pairs.Equations ( 8) and ( 11) can be used to calculate V 2 and V 3 because the Buckingham potential 35 and the Wolf summation 36,37 are both pairwise.Finally, the thermal conductivity of the deformed nanowire is calculated based on the configurational data and the force constants.
Figure 3 compares the thermal conductivity predicted by the model with that calculated using the Green-Kubo method in Ref. 23.The diameters of the wires considered are 2.91 and 3.55 nm and the temperature is 300 K.The solid lines indicate model prediction and the symbols with error bars indicate Green-Kubo results.The dependence of thermal conductivity on strain and structural change is successfully captured by the model prediction which is in agreement with the Green-Kubo results.As shown in Fig. 3(a) (diameter d = 2.91 nm), the thermal conductivity decreases by 15% as the strain increases from zero to 0.06 under tensile loading, and increases by 30% as the strain decreases from zero to -0.06 under compression.The model captures the dependence of thermal conductivity on structural deformation quite well in both loading and unloading.The thermal conductivity of the nanowires at a strain of 0.01 in the unloading process (when both TS-and WZ-structured regions coexist) is 25% lower than that of the WZ-structured nanowires at the same strain under loading.The same trend is seen for the nanowire with a diameter of d = 3.55 nm in Fig. 3(b), with good agreement between the model prediction and the Green-Kubo results.The thermal conductivity of the nanowires that contain the inversion domain boundaries 23 at zero strain in the unloading process is 20% lower than that of the unstressed WZ wires.
In Fig. 3, the thermal conductivity predicted by the model does not follow a smooth line.This phenomenon comes from the uncertainty in MD simulations.The thermal conductivity values from  The thermal responses of the same nanowire (d = 3.55 nm) during loading and unloading at T = 1500 K are also analyzed.The phase transformation from WZ to TS is observed, as shown in Fig. 4(b).At 1500 K, the thermal conductivity predicted by the model is similar to that given by the Green-Kubo method, as shown in Fig. 5.As the strain increases from zero to 0.12 under tensile loading, the thermal conductivity decreases by 30% according to Green-Kubo and 40% according to the model.In the unloading process, the interior region of the nanowire transforms from TS to WZ, as shown in Figs.4(c) and 4(d).The thermal conductivity of the nanowire at a strain of 0.02 in the unloading process is 18% lower than that for the WZ-structured wires at the same strain under tensile loading.The model predicts the thermal conductivity of the same nanowire as 25% lower than that of the WZ wires under tensile loading.The model prediction is in good agreement with the Green-Kubo results at high temperatures as well.The percent error between the model prediction and the Green-Kubo results for the nanowires in Fig. 5 is 11.7%.

IV. SUMMARY
An analytical model is developed to evaluate the thermal conductivity of semiconducting compounds during deformations that involve structural changes.This model reflects the fact that the ensemble averages of atomistic behaviors are determined by the atomistic potentials of materials.The evaluation of thermal conductivity assumes that the second-order and third-order force constants dominate phonon behaviors during structural deformation.The force constants are calculated by the interatomic potential and the atomistic configuration obtained in the MD calculations for mechanical response.The dependence of thermal conductivity on strain and crystal structure is described through the use of deformed atomistic configurations.Application of this approach to GaN nanowires 2.91-3.55nm in size over the temperature range of 300-1500 K has yielded results that are in good agreement with data obtained using MD-based Green-Kubo method.This model is 50 times more efficient than MD calculations.
FIG. 1.A schematic view of phase transformation of [0001]-oriented GaN nanowires during tensile loading and unloading.

FIG. 2 .
FIG.2.Schematic illustration of the approach for predicting thermal conductivity of deformed structures using atomistic configurations from MD simulations and an atomistic potential.

FIG. 3 .
FIG. 3. Thermal conductivities of GaN nanowires at 300 K calculated using the Green-Kubo method and the model: (a) nanowire with diameter d = 2.91 nm; and (b) nanowire with diameter d = 3.55 nm.Error bars denote standard deviation of conductivity.

FIG. 4 .
FIG. 4. Atomistic configurations of a nanowire with diameter d = 3.55 nm at T = 1500 K during loading and unloading: (a) WZ-structured nanowire at zero strain; (b) TS-structured nanowire at a strain of 0.1; (c) WZ-TS structured nanowire at a strain of 0.02; and (d) WZ-TS structured nanowire at a strain of -0.02.

FIG. 5 .
FIG. 5. Thermal conductivities of GaN nanowires with diameter d = 3.55 nm at 1500 K calculated using the Green-Kubo method and the model.Error bars denote standard deviation of conductivity.