Molecular simulations and lattice dynamics determination of Stillinger-Weber GaN thermal conductivity

The bulk thermal conductivity of Stillinger-Weber (SW) wurtzite GaN in the [0001] direction at a temperature of 300 K is calculated using equilibrium molecular dynamics (EMD), non-equilibrium MD (NEMD), and lattice dynamics (LD) methods. While the NEMD method predicts a thermal conductivity of 166 6 11 W/m ! K, both the EMD and LD methods predict thermal conductivities that are an order of magnitude greater. We attribute the discrepancy to signiﬁcant contributions to thermal conductivity from long-mean free path phonons. We propose that the Gr € uneisen parameter for low-frequency phonons is a good predictor of the severity of the size effects in NEMD thermal conductivity prediction. For weakly anharmonic crystals characterized by small Gr € uneisen parameters, accurate determination of thermal conductivity by NEMD is computationally impractical. The simulation results also indicate the GaN SW potential, which was originally developed for studying the atomic-level structure of dislocations, is not suitable for prediction of its thermal conductivity.


I. INTRODUCTION
Due to its electronic and thermal transport properties, GaN is an important component in microelectronic devices such as LEDs, laser diodes, and high electron mobility transistors. [1][2][3][4][5][6] Its high thermal conductivity is crucial for efficient heat dissipation in high-power GaN-based microelectronics. The room-temperature thermal conductivity, k, along the [0001] direction of single-crystal wurtzite GaN reported in experiments [7][8][9][10][11][12] varies from 170 to 250 W/mÁK. The GaN crystals used in experiments contain a variety of defects such as dislocations, impurities, and isotopes. Although the defect density in GaN crystals has been significantly reduced in the past decade, [7][8][9][10][11][12] it is still several orders of magnitude greater than that in typical silicon wafers. 13 Measured thermal conduciveness are thus smaller than that characterizing a perfect, single-crystal sample.
Atomic-level simulations and calculations have been applied to calculate the k of the perfect (i.e., single crystal and isotopically pure) wurtzite GaN crystal in the [0001] direction at a temperature of 300 K. [14][15][16] Using first principles calculations and lattice dynamics (LD) methods, Lindsay et al. predict a thermal conductivity of 385 W/mÁK. 16 Their prediction when natural isotopes are included is 239 W/mÁK, in agreement with the experimental measurements. Zhou et al. calculated the k of Stillinger-Weber (SW) GaN using the nonequilibrium molecular dynamics (NEMD) method. 14 By linear extrapolation of 1/k vs. 1/L to the 1/L!0 limit, where L is the sample length in the [0001] direction, they obtained a value of 185 W/mÁK, which falls within the range of the experimental values, although their simulations did not contain the defects found in experimental samples. Subsequent work by Sellan et al., however, showed that the validity of the linear extrapolation procedure in the work of Zhou et al. may be questionable since phonons with mean free paths (MFPs) much longer than the sample lengths in NEMD simulations can contribute significantly to bulk k of high-modulus semiconductors. 17 Using equilibrium MD (EMD) simulations and the Green-Kubo (GK) formula, Kawamura et al. found k (averaged over three directions) of SW GaN to be between 310 and 380 W/mÁK, depending on the data range selected in the evaluation of the integral of the heat flux autocorrelation function. 15 While this value is significantly higher than the experimental values, the heat flux autocorrelation function in the EMD simulations of Kawamura  In this work, we determine k of the perfect GaN wurtzite crystal along the [0001] direction using EMD, NEMD, and LD methods. By comparing the results, we evaluate the validity of the linear extrapolation procedure for NEMD prediction of k and assess the suitability of the existing SW potential 18 for quantitative modeling of thermal transport in GaN.

A. GaN model
We study the GaN wurtzite hexagonal crystal structure. The atomic interactions are modeled using the SW potential with parameters developed by Bere and Serra. 18 Although this potential was originally developed for studying the atomic-level structure of dislocations in GaN, Zhou et al. [This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP: used NEMD method to show that it predicts GaN thermal conductivities in agreement with experimental measurements. This agreement, however, could be a coincidence. Sellan et al. 17 found that while the thermal conductivity of SW Si predicted by the NEMD method agrees with experimental measurements, it underestimates the contribution from phonons with long MFPs and thus underpredicts the EMD thermal conductivity by a factor of 2.5 at a temperature of 500 K. As GaN is more thermally conductive than silicon and our simulations are conducted at a temperature of 300 K, the difference between EMD and NEMD simulation results could be even greater.
In our simulations, only nearest neighbor interactions are considered, which is appropriate at room temperature where the bonding structure does not change. Using the SW potential, our MD simulations yield the GaN lattice constants at a temperature of 300 K and a pressure of 1 atm to be a ¼ 3.194 Å and c ¼ 5.217 Å , which agree well with experimental data. 19 The ½1 100, ½11 20, and ½0001 directions are aligned with the x, y, and z directions of the simulation box. An orthogonal unit cell is defined with dimensions in the x, y, and z directions of 2acos(p/6), a, and c.

B. NEMD determination of k
For the NEMD simulations, the simulation box has a cross sectional area of 3 unit cells in the x direction by 5 unit cells in the y direction. With periodic boundary conditions (PBCs) applied in the x and y directions, this cross sectional area was found to be large enough to remove size effects related to the finite cross sectional area. 14 In the z direction, the simulation box is bordered by 10-unit-cell long heat source and heat sink regions and free boundaries.
The system is first equilibrated at a temperature of 300 K for 1 ns using the Berendsen thermostat. 21 The velocity Verlet algorithm with a time step size of 1 fs is used to integrate the equations of motions in all MD simulations. 20 After the system reaches thermal equilibrium, the global thermostat is turned off and a heat flux of q ¼ 15 GW/m 2 is applied by adding a constant amount of energy to the heat source and removing the same amount of energy from the heat sink at each time step using the velocity rescaling method. 22 We allow 6 ns for the system to reach a steady state and then run for 50 ns for data collection and averaging. We calculate the temperature profile in bins with width of 5 unit cells in the heat flow direction. A typical temperature profile at steady state is shown in Fig. 1. The temperature gradient is obtained by fitting to the linear region of the temperature profile. The thermal conductivity is obtained from the Fourier law as The uncertainty in k is determined by analyzing four independent simulations.
To determine the dependence of k on L, the number of unit cells along the z direction is varied from 75 to 300, which corresponds to L varying from 38 to 150 nm. This range of sample lengths is similar to that in the work of Zhou et al. 14 In Fig. 2, we show that, for this range of L, the inverse of thermal conductivity, 1/k, is a linear function of the inverse of the simulation cell length, 1/L. Using the linear extrapolation method proposed by Schelling et al., 23 we obtain k NEMD ¼ 166 6 11 W/mÁK, which is consistent with the value of 185 W/mÁK reported by Zhou et al., 14 who used the same SW potential for GaN but with PBCs applied in all three directions. We note that Zhou et al. 14 did not provide the uncertainty associated with their result.
Our NEMD prediction of the room-temperature thermal conductivity and that of Zhou et al. 14 are a factor of two lower than that obtained by Kawamura et al. using EMD. 15 Sellan et al. 17 pointed out that when low-frequency, long-MFP phonons contribute significantly to the thermal conductivity, the linear extrapolation method used with NEMD for small and moderate simulation cell lengths can severely underestimate the thermal conductivity. Increasing the simulation cell length further while keeping the cross section constant, however, may not solve the problem as the thermal conductivity might diverge with length, as typical of onedimensional systems. 24 To determine if this problem also exists for GaN, we increased the system length to 800 unit cells, which corresponds to L ¼ 418 nm, and calculated the thermal conductivity for this size to be 177 6 9 W/mÁK. This data point, also shown in Fig. 2, does not fall on the linear fit of 1/k vs. 1/L obtained for data points with smaller L. A similar nonlinear behavior was also found by Zhou et al. for SW GaN at a temperature of 800 K. 14 Such a nonlinear trend could reflect the true contributions of low-frequency, long-MFP phonons to k, 14 but could also be caused by one-dimensional phonon effects. 24 This artifact can be possibly mitigated by increasing the cross sectional area. Such an approach, however, significantly increases the total number of atoms in the system and thus the computational cost, which makes an accurate NEMD prediction of k challenging.
Based on the results and analyses presented in this section, we find it difficult to eliminate system-size effects on the NEMD prediction of k of bulk GaN. In Sec. II C, we present calculations of k of GaN using the EMD method, which, for silicon, has been shown to be able to eliminate size effects in systems with less than 10 000 atoms. 17,25 C. EMD determination of k The EMD method evaluates the thermal conductivity via the Green-Kubo formula 20 where V is the volume of the system, T is the temperature, k B is the Boltzmann constant, t represents time, hÁ Á Ái denotes ensemble average, and q z is the heat flux in the z direction (i.e., [0001]), which can be computed from 26,27 In the above equation, v i and E i are the velocity and total energy of atom i, r ij is the position vector from atom j to atom i, f ij is the two-body force acting on atom i by atom j, and f jik is the three-body force acting on atom i due to atoms j and k.
PBCs are applied in all three directions in our EMD simulations. The system is first equilibrated at a temperature of 300 K and a pressure of 1 atm for 1 ns using the Berendsen thermostat. 23 After the system reaches thermal equilibrium, the thermostat is turned off and data are collected from simulations run in the microcanonical ensemble for 300 ns. Such a long simulation time is needed to limit statistical error in the heat flux autocorrelation function (HFACF).
The HFACF and its running integral for a system built from 10 Â 10 Â 10 unit cells are shown in Fig. 3(a). The HFACF obtained using Eq. (3) definition of q z exhibits fast oscillations that are associated with high-frequency optical modes. Landry et al. 28 pointed out that these fast oscillations are related to the first sum in Eq. (3). Since the atoms in our perfect GaN crystal do not diffuse during the simulations, the heat flux can be written as 28 where r ij,0 is the separation vector from atom j to atom i based on their equilibrium positions. As shown in Fig. 3(b), the oscillations in the HFACF are significantly reduced when Eq. (4) is used to calculate q z and the HFACF decays to zero when the correlation length reaches 600 ps. The integrals of the two HFACFs converge to the same value of k EMD ¼ 1190 6 85 W/mÁK. The uncertainties of k are determined from the analysis of four independent simulation runs.
To investigate size effects on the EMD prediction of k of bulk GaN, we carried out similar simulations with systems containing 4 Â 7 Â 4 and 6 Â 10 Â 6 unit cells. As shown in Table I, there are no discernible size effects on k for the three system sizes considered. Most strikingly, the EMD determined thermal conductivity is eight times greater than that determined in NEMD. A similar behavior was observed in the prediction of k of SW silicon at a temperature of 500 K, where the NEMD result obtained by the linear extrapolation procedure underpredicted the EMD result by a factor of 2.

D. Prediction of k using LD calculations
To assess the validity of our EMD-determined thermal conductivity, we performed lattice dynamics calculations to obtain an independent prediction. The phonon contribution to the thermal conductivity in the n direction of a crystalline material can be obtained by solving the Boltzmann transport equation and using the Fourier law as The summation in Eq. (5) is over all the phonon modes in the first Brillouin zone. On the right-hand side, c ph;i is the volumetric specific heat, v g;n;i is the n-component of the phonon group velocity vector v g;i , and s i is the phonon lifetime, which we calculate under the relaxation time approximation. For GaN modeled from first principles calculations, Lindsay et al. found that the relaxation time approximation underpredicts the thermal conductivity from a full solution 29-31 of the Boltzmann transport equation by 7%. 16 The specific heat in classical statistics, as realized in the MD simulations, is k B =V. The group velocity vector is v g;i ¼ @x i =@j, where x i is the phonon frequency and j is the phonon wave-vector. We obtain phonon frequencies and phonon lifetimes from harmonic and anharmonic LD calculations. The details of the LD calculations can be found elsewhere. 32 We discretized the first Brillouin zone using an N o Â N o Â N z uniformly spaced phonon wave-vector grid. To obtain a converged value of the [0001] bulk thermal conductivity, we varied the grid density while keeping the ratio N z =N o close to a=c. The variation of 1=k z BTE with 1=N z is plotted in Fig. 4. The thermal conductivity decreases with increasing grid resolution as more phonon modes become available to satisfy the energy and momentum conservation rules. A similar trend was found by Jain and McGaughey for a family of model compound semiconductors. 33 The bulk thermal conductivity of 1; 8256173 W/mÁK is obtained by fitting a linear curve to 1=k z BTE vs 1=N z and extrapolating it to N z ! 1. This value is one order of magnitude larger than the NEMD thermal conductivity and 38% larger than the size-averaged EMD value, which may be due to only including up to three-phonon processes in the LD calculation. Nevertheless, the LD results combined with the EMD simulation results indicate that the thermal conductivity of SW GaN is in the range of 1500 W/mÁK, which is an order of magnitude above the value obtained from the NEMD simulations.
We note that while the EMD thermal conductivity prediction is converged for relatively small system sizes (Table  I), the LD prediction for equivalent system sizes is not yet converged (Fig. 4). A similar behavior was observed by Sellan et al. 17 for Lennard-Jones (LJ) argon and SW silicon. The origin of these different size effects is unknown and warrants further study.

E. Gr€ uneisen parameter
If a crystal contains no defects, free electrons, or internal interfaces, its thermal resistance is only due to intrinsic phonon-phonon scattering that results from anharmonicity. The magnitude of the anharmonicity can be quantified by the phonon-mode dependent Gr€ uneisen parameters c, 34 which are essentially a measure of the magnitudes of the thirdorder derivatives of the energies of the atomic interactions. When c is close to zero, the mode anharmonicity is small and the associated phonon lifetime is large. 35,36 Therefore, small values of c indicate that contributions from long-MFP phonons can be important. If the system length L in an NEMD simulation is not much longer than the MFPs of the majority of the phonons that significantly contribute to k, the NEMD method can severely underestimate thermal conductivity. Such an underestimation is observed in the current work for SW GaN and was observed for SW silicon by Sellan et al., 17 but was not observed for LJ argon. 32 As such, we will evaluate c's for low-frequency modes for LJ argon, SW GaN, and SW silicon to determine if they are good predictors of a discrepancy between EMD and NEMD thermal conductivity values.
Based on Debye theory, the low-frequency Gr€ uneisen parameter c can be calculated by 34 where P is the system pressure. To evaluate the derivatives in Eq. (6), we first minimize the system energy at a temperature of 0 K and a pressure of 1 atm. Then, we hydrostatically strain the simulation cell from À0.5% to 0.5% in volume, equilibrate the system to a temperature of 0 K using the Berendsen thermostat, 23 and then determine P at each volume from MD simulations. Using Eq. (6), we obtain c values of 1.32 for SW GaN, 1.30 for SW silicon, and 3.83 for LJ argon. This result suggests that the low-frequency Gr€ uneisen parameter is a good indicator for the accuracy of the bulk k predicted using NEMD simulations via the linear extrapolation procedure. If the low-frequency Gr€ uneisen parameter is small ($less than 2), there is a significant contribution to the bulk k from long-MFP phonons and accurate prediction of bulk k by NEMD method is difficult.

III. DISCUSSION AND CONCLUSIONS
The determination of the thermal conductivity of SW GaN via EMD simulations and LD calculation yielded $1500 W/mÁK, which is an order of magnitude larger than that obtained via NEMD and the linear extrapolation procedure. The underlying origin of this discrepancy is significant contribution of low frequency/long mean free path phonons to k and is associated with the highly harmonic nature of the GaN SW potential. Similar behavior was found by Sellan et al. for SW silicon. 17 While in principle NEMD simulations should provide convergent results, increasing just the length of the simulation cell is not sufficient, as one-dimensional effects associated with the sampling of the Brillouin zone may lead to artificially divergent thermal conductivity. 24 The proper application of NEMD to determine the bulk thermal conductivity of a highly harmonic solid would require increasing both the length and the cross-sectional area of the simulation cell. Such a procedure, however, would be very costly and to our knowledge has not been reported for high thermal conductivity materials.
We suggest that a relatively simple assessment of the applicability of the standard NEMD method for thermal conductivity determination is via evaluation of the lowfrequency Gr€ uneisen parameter. For values $<2, the EMD method should be used, as this method is characterized by much less severe finite size effects. Even for EMD, however, size effects should be examined to ensure convergence.
Finally, we comment on comparison with the experimental thermal conductivity values. The proper evaluation of the thermal conductivity of SW GaN at a temperature of 300 K yields a value an order of magnitude larger than that measured in experiment. This discrepancy can be partially attributed to the fact that a real GaN sample has defects, which reduce thermal conductivity, while in our MD simulations and LD calculations perfect crystal structures are used. This difference, however, is probably a minor effect. It is most likely that the SW GaN potential underestimates anharmonicity and thus predicts a larger thermal conductivity. This conclusion is corroborated by the agreement between the first-principles based predictions of Lindsay et al. 16 with the experimental measurements. This result indicates that the GaN SW potential, which was originally developed for studying the atomic-level structure of dislocations, 18 is not suitable for prediction of its thermal conductivity.