How does the accuracy of interatomic force constants affect the prediction of lattice thermal conductivity
11 How does the accuracy of interatomic force constants affect the prediction of lattice thermal conductivity
Han Xie, Xiaokun Gu, and Hua Bao * University of Michigan-Shanghai Jiao Tong University Joint Institute,
Shanghai Jiao Tong University, Shanghai 200240, China Institute of Engineering Thermophysics, School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
ABSTRACT
Solving Peierls-Boltzmann transport equation with interatomic force constants (IFCs) from first-principles calculations has been a widely used method for predicting lattice thermal conductivity of three-dimensional materials. With the increasing research interests in two-dimensional materials, this method is directly applied to them but different works show quite different results. In this work, classical potential was used to investigate the effect of the accuracy of IFCs on the predicted thermal conductivity. Inaccuracies were introduced to the third-order IFCs by generating errors in the input forces. When the force error lies in the typical value from first-principles calculations, the calculated thermal conductivity would be quite different from the benchmark value. It is found that imposing translational invariance conditions cannot always guarantee a better thermal conductivity result. It is also shown that Grüneisen parameters cannot be used as a necessary and sufficient criterion for the accuracy of third-order IFCs in the aspect of predicting thermal conductivity.
KEYWORDS:
Thermal conductivity, Peierls-Boltzmann transport equation, Interatomic force constants, Grüneisen parameters * Corresponding author: [email protected]
1. INTRODUCTION
In crystalline semiconductors and insulators, phonons ( i.e. lattice vibrations) are the major heat carriers, so the thermal conductivity can be calculated with the knowledge of phonon properties. Recently there are growing interests to predict the lattice thermal conductivity by solving Peierls-Boltzmann transport equation (PBTE), either under single mode relaxation time approximation (SMRTA) or with full iterative solution.[1,2] With interatomic force constants (IFCs) extracted from first-principles calculations as the input, this method has been widely used to calculate the thermal conductivity of three-dimensional (3D) materials.[3–26] The good agreement between these calculated values and measured experimental data proves its accuracy and reliability.[3–21] The discovery of ultrahigh thermal conductivity of graphene[27] has stimulated a growing research interest in the thermal conductivity of two-dimensional (2D) materials. The first-principles PBTE method for predicting thermal conductivity has been directly applied to 2D materials like graphene,[28–31] silicene,[32–35] phosphorene,[36–40] MoS ,[41–43] borophene,[44] etc. However, previous calculations show quite different results in different works. For example, the predicted thermal conductivity for black phosphorene from SMRTA[36–40] at room temperature along armchair direction varies a lot from 5.46 to 33 W/mK, and the result along zigzag direction spans a rather large range between 15.33 and 83.5 W/mK. Similar discrepancy is also observed in the case of MoS , whose thermal conductivity at room temperature was predicted from SMRTA to be 83 W/mK by Li et al. [41] as well as Gu and Yang.[42] In comparison, Yan et al. [43] got a result of 35.5 W/mK with the same method. The lack of consistency in the predicted thermal conductivity may arise from using iterative method instead of SMRTA,[30,31] imposing translational invariance conditions to third-order IFCs,[15,35] different exchange-correlation functionals used in first-principles calculations,[45] enforcing a quadratic branch in the dispersion of 2D materials,[44] etc . In materials where resistive Umklapp scattering processes dominate in the phonon-phonon scattering, both SMRTA and iterative method should give similar results.[1,2] For graphene, iterative method yields a much larger thermal conductivity value than SMRTA because the momentum-conserving normal scattering processes also play an important role in phonon transport,[30,31] which has been explained by the hydrodynamic phonon transport at room temperature.[46,47] Regarding the translational invariance conditions ( i.e. acoustic sum rules), Lindsay et al. [15] demonstrated that imposing it to third-order IFCs would play an important role in determining thermal conductivity. Our previous results for silicene also showed its importance.[32,35] Translational invariance conditions will affect the calculated phonon scattering rate, especially for long-wavelength phonons near the Brillouin zone center. Concerning the exchange-correlation functionals, Jain and McGaughey[45] studied their effects on predicting the thermal conductivity of crystalline silicon from first-principles calculations. Different exchange-correlation functionals could lead to an under-prediction or over-prediction within 20% of the experimental value. Recently, Carrete et al. [44] showed that the second-order IFCs directly extracted from first-principles calculations might yield problematic phonon dispersion curve for 2D materials. The dispersion curve of unstrained 2D materials can be demonstrated to contain a quadratic branch near Brillouin zone center[48] but the second-order IFCs directly extracted from first-principles calculations would yield linear ones. By generating physically sound IFCs, they showed that the thermal conductivity result could be quite different from the value calculated with the raw IFCs from first-principles calculations.[44] Previous discussions explained part of the inconsistency in the predicted thermal conductivity of 2D materials from first-principles calculations. However, with all the above mentioned reasons considered, based on our own testing (unpublished), discrepancy in predicted thermal conductivity could still exist for some 2D materials. As will be shown later, the accuracy of third-order IFCs also plays an important role on the predicted thermal conductivity. For first-principles calculations, the accuracy might be affected by energy cutoff, k -point grid, reciprocal space projection technique, aliasing errors, discretization errors, etc .[49] For example, the attainable fractional precision[49] in forces using VASP is . Currently, the most widely used method to extract interatomic force constants generally takes force as the input parameter, and uses finite-difference method to calculate the IFCs. The raw IFCs generally do not satisfy the translational invariance conditions, so these conditions are artificially imposed by adding small compensation to each term. All these processes will induce additional uncertainty to the IFCs, especially the third-order IFCs. In this work, we will discuss how the accuracy of IFCs could affect the predicted lattice thermal conductivity values. Due to the large uncertainty in first-principles calculations, we used classical potential to do such an investigation. Classical potential has the advantage that it has an explicit analytical form, so that the error only comes from numerical computation and can be reduced to a negligible amount to get an “accurate” thermal conductivity result as benchmark. Based on the benchmark case from classical potential, inaccuracies are artificially introduced to third-order IFCs and effects of these inaccuracies on the predicted thermal conductivity are investigated. SMRTA is used to calculate thermal conductivity due to the computational cost consideration and also because phonon relaxation time is well defined with this approach. Previously, Grüneisen parameters have been used as a simple test for the accuracy of the third-order IFCs.[28] The applicability of this criterion is also examined. Our result will shed some light on predicting thermal conductivity from first-principles calculations. In what follows, we describe the simulation methods and details in Sec. 2. Simulation results for silicon, graphene, and silicene are shown in Sec. 3. Our conclusions are summarized and discussed in Sec. 4.
2. SIMULATION METHODS AND DETAILS 2.1 Single mode relaxation time approximation method
For a periodic crystal structure, the potential energy can be expanded as the Taylor series[50,51] ij i j ijk i j kij ijk
U U u u u u u , (1) where U is the equilibrium potential energy. i u , j u , and k u are the displacements of i -th atom in direction, j -th atom in direction, and k -th atom in direction, respectively. ij is the second-order IFC and ijk is the third-order IFC. Even higher-order IFCs are neglected in this equation. Physically correct second-order and third-order IFCs have to satisfy the point/space group symmetry relations, translational invariance conditions, and rotational invariance conditions, which are shown in the Appendix. Comparison of Grüneisen parameters from second-order and third-order IFCs has been used as a simple test for the accuracy of third-order IFCs. The calculation of Grüneisen parameters and the definition of relative difference are also shown in the Appendix. The force acting on each atom is ii U F and for a structure under equilibrium state, i F 0 . With the IFCs as the input, the thermal conductivity of semiconducting or insulating materials can be calculated from SMRTA with the following equation , l ph c v v , (2) where denotes different phonons that can be distinguished by wave vector q and phonon branch . , ph c is the volumetric phonon specific heat. v and v are the phonon group velocity in and directions, respectively. , ph c and v can be calculated with the second-order IFCs as the input.[52] is the phonon relaxation time, i.e. the inverse of phonon scattering rate. In our calculation of , only phonon-phonon scattering is considered. More specifically, only three-phonon scattering is considered. can be calculated with both second-order and third-order IFCs as the input. For more details about the method we refer the reader to Refs. [52–54]. The second-order and third-order IFCs as the input for SMRTA method were calculated from classical potential. Careful tests were carried out to reduce the numerical error. GULP package[55] was first used to optimize the primitive unit cell of silicon, graphene, or silicene. After that, a supercell was constructed and the lattice constant was re-optimized with LAMMPS package.[56] In this re-optimization process, the bisection method was used to find the lattice constant corresponding to the lowest energy state. Forces acting on each atom were computed from LAMMPS through the analytical derivatives of the potential function,[57] which would be free of truncation error coming from numerical differentiation. As the input for calculating IFCs, forces were output with sixteen significant digits to retain accuracy. In order to reduce the truncation error, fourth order accuracy method was used instead of central difference method (see the Appendix for the details) to compute the second-order and third-order IFCs with our own in-house code and revised THIRDORDER.PY,[54] respectively. These third-order IFCs were not modified by adding small compensation because they already satisfied translational invariance conditions to a reasonable extent. Point/Space group symmetry conditions were enforced and utilized to reduce computational cost. This result was then used to obtain the benchmark thermal conductivity. Inaccuracies in the third-order IFCs were simulated by either truncating digits or adding random numbers to the forces since forces were the input for calculating third-order IFCs. For truncation case, the force digits after a prescribed force accuracy value of each term were thrown. For example, when 0.1234567890123456 was truncated with a force accuracy of , the value would be 0.123. For random addition case, random numbers between the negative and positive force accuracy values were added to the force terms, which were generated with the rand function in standard C library. Five such random addition cases were conducted and the thermal conductivity was calculated for each of them. With the inaccuracies introduced, both original and modified third-order IFCs were used to calculate thermal conductivity. Modified third-order IFCs were produced by adding small compensation to each term using Lagrange multiplier method[12,54,58] in order to satisfy translational invariance conditions. In all the above-mentioned simulations, ShengBTE package[54] was used to calculate thermal conductivity with SMRTA. More simulation details that are related to different materials can be found in Sec. 3.
3. SIMULATION RESULTS 3.1 Silicon
Tersoff potential[59] was used to describe the Si-Si interactions in silicon. The face-centered cubic primitive unit cell and supercell were constructed and optimized successively. The lattice constant for equilibrium configuration with the lowest energy was found to be 5.4321 Å. Because a cutoff distance exists in the potential function, long-range IFCs are exactly zero. Accordingly, only second-nearest neighbors were considered in computing third-order IFCs.[60] Table 1 shows the maximum absolute values of translational invariance and rotational invariance tensors (2) T , (3) T , (2) 2 T( ) I I , (3) ) T(3 I I (see the Appendix for the details). (2) T and (2) 2 T( ) I I are the translational and rotational invariance tensors for second-order IFCs, respectively. Meanwhile, (3) T and (3) ) T(3 I I are for third-order IFCs. “Original IFC3” in the first column refers to the case where nothing was added to the third-order IFCs and “Modified IFC3” refers to the case where small terms were added to make the third-order IFCs satisfy the translational invariance conditions. Modified IFC3 case is shown in the table merely for comparison purpose. Theoretically, these tensors should be equal to if invariance conditions are strictly satisfied. However, the theoretical result is unattainable because error always exists in numerical computation. According to our experience, the results in Table 1 are much smaller than the commonly attainable values of first-principles calculations and we believe from these values that translational and rotational invariance conditions are well satisfied. By comparing the results between the two cases we find that both invariance conditions are better satisfied for modified third-order IFCs.
36 36 36 q mesh was used to calculate the thermal conductivity at 300 K. Our results were about 319 W/mK for both cases and later we would take this value as the benchmark result. We noticed that previous calculation[60] based on the same method yielded a value about 15% higher than our result. Such a difference may arise from the better accuracy of our calculated IFCs and the denser q mesh we have used. Furthermore, when comparing the results between the original IFC3 case and modified IFC3 case, we find that the difference is so small that it should be negligible, which might serve as the circumstantial evidence that our calculated thermal conductivity values are accurate. Table 1 Translational, rotational invariance conditions and thermal conductivity for silicon from Tersoff potential (benchmark case) Case (2) max T i (eV/Å ) (3) max T ij (eV/Å ) (2), , max(2) I I i i (eV/Å) (3) (3), , max I I ij ij (eV/Å ) Thermal conductivity (W/mK) Original IFC3 x axis means better force accuracy. The horizontal solid gray line without symbols shows our benchmark result. Truncation cases are plotted as diamond and circular symbols connected with dashed line. The results from five random addition cases are plotted as floating bars in the figure. It can be seen from Fig. 1 that when the force accuracy value is on the order of or smaller than eV/Å, all the calculated thermal conductivity values agree well with the benchmark result. The maximum error for these cases is less than 2% and therefore modified IFC3 is not quite necessary for getting an accurate result. When the force accuracy value is between and eV/Å, original IFC3 would yield 4.5%-20% error in thermal conductivity result but modified IFC3 can correct this error to less than 2%. For silicon, we take ±15% of the benchmark thermal conductivity as an acceptable range, and later we also use this criterion for graphene and silicene. For eV/Å force accuracy case, modified IFC3 can reduce the error to this acceptable range for both truncation case and random addition cases. When the force accuracy is even worse, the absolute value of the error becomes larger than 15% and looks quite obvious in the figure. Therefore, we believe a reasonable thermal conductivity can be obtained for silicon from modified IFC3 with force accuracy at least on the order of eV/Å. Fig. 1 Thermal conductivity of silicon from Tersoff potential vs. force accuracy We use the results from truncation cases for further discussions. Otherwise, five random addition cases need to be discussed for each force accuracy. In addition, random addition cases show similar trend as truncation cases. When the force accuracy is on the order of eV/Å, original IFC3 will yield translational invariance condition (3) max 3 T 0.468 e /ÅV ij and modified IFC3 can reduce this value to
10 3 . Modification of third-order IFCs can make translational invariance conditions much better satisfied. Since it has been shown in Fig. 1 that modified IFC3, not original IFC3, gives reasonable thermal conductivity for this force accuracy, Grüneisen parameters from modified IFC3 cases are plotted in Fig. 2, together with those from second-order IFCs (labelled as IFC2). The x axis indicates the reduced coordinates[61] of reciprocal lattice vectors. It can be seen from the figure that Grüneisen parameters calculated from IFC2 and benchmark IFC3 agree quite well, with the relative difference e =0.02% (see the Appendix for the definition of this relative difference). This can also serve as a circumstantial evidence that the third-order IFCs of the benchmark case are quite accurate. The modified IFC3 for eV/Å accuracy still gives Grüneisen parameters that agree reasonably with IFC2, with a relative difference of 2.87%. It can be obviously seen from the figure that Grüneisen parameters are different between IFC3 from eV/Å accuracy case and IFC2. Since eV/Å force accuracy gives reasonable thermal conductivity, Grüneisen parameters may serve as one of the good criteria for the accuracy of modified third-order IFCs for silicon. Fig. 2 Grüneisen parameters of silicon from Tersoff potential Graphene has attracted lots of scientific research interests since its discovery and has stimulated a growing research interest in 2D materials. The optimized Tersoff potential for graphene[62] was used to describe the C-C interactions. A hexagonal primitive unit cell structure was first constructed and optimized. Then a supercell was constructed and re-optimized. The lattice constant for the lowest energy state was finally found to be 2.492049 Å. For 2D materials, we used the original second-order IFCs to compute the phonon dispersion curve and check whether it would be consistent with the theoretical demonstration that a quadratic branch exists near Brillouin zone center.[48] The phonon frequencies for 1000 q points from Γ to M were calculated and part of it was plotted in Fig. 3. Please note that the unit for x axis is 2π/ a , where a is the lattice constant. We can see from the figure that the dispersion curve for out-of-plane acoustic (ZA) mode looks like quadratic while the longitudinal and transverse acoustic (LA/TA) modes have linear dispersions. In order to further confirm the trend for ZA mode, the group velocities along Γ to M direction (derivatives of frequency) of ZA branch were calculated with finite difference method and were plotted as red circular dots in the inset of Fig. 3. The solid black line in the inset is a straight line crossing the origin point and the last point in the figure. It can be seen from the inset that all the dots for ZA group velocities sit on the straight line. Therefore, we believe the ZA phonon mode has a quadratic dispersion curve near Brillouin zone center, which agrees with the theoretical demonstration.[48] Fig. 3 Phonon dispersion curve of graphene calculated from optimized Tersoff potential for acoustic modes from Γ to M direction around Γ point. Inset: Group velocity of ZA mode along Γ to M direction around Γ point. In order to calculate third-order IFCs, second-nearest neighbors were considered. Table 2 shows the maximum absolute values of translational and rotational invariance tensors. The values of (3) max T ij and (3) (3), , max I I ij ij are larger than those of silicon, but we should notice that the C-C bonding in graphene is stronger than Si-Si bonding in silicon. These values are still very small either according to our experience or by comparing with the critical force accuracy case for silicon, i.e. eV/Å case. Modification of third-order IFCs can increase the accuracy.
201 201 1 q mesh was used to calculate the thermal conductivity at 300 K. The results for original IFC3 and modified IFC3 cases are almost the same and we take the benchmark result as 633 W/mK. For graphene, using iterative method instead of SMRTA should give a much larger thermal conductivity[30] and it can be seen that the result from iterative method in Ref. [62] is much larger than our result here. Due to the reasons explained before, we just consistently use SMRTA method to do the accuracy analysis for silicon, graphene, and silicene. Table 2 Translational, rotational invariance conditions and thermal conductivity for graphene from optimized Tersoff potential (benchmark case) Case (2) max T i (eV/Å ) (3) max T ij (eV/Å ) (2), , max(2) I I i i (eV/Å) (3) (3), , max I I ij ij (eV/Å ) Thermal conductivity (W/mK) Original IFC3 y axis because some results span a very large range. The horizontal gray line without symbols is the benchmark result. It can be seen that when the force accuracy value is on the order of or smaller than eV/Å, all the calculated thermal conductivity values agree well with the benchmark result. Even the largest error is less than 5%. When the force accuracy is eV/Å, the results from original IFC3 are within ±15% of the benchmark thermal conductivity, no matter for truncation cases or for random addition cases. In accordance with silicon, these results are taken as reasonable values. Anomalously, the results from modified IFC3 cases deviate more from the benchmark result than those from original IFC3. For truncation case, the error grows from 13% to 28%. For random addition case, modified IFC3 will yield a larger span of thermal conductivity from 555 to 812 W/mK than original IFC3 in the range of 551 to 671 W/mK. The results from modified IFC3 are out of the reasonable range we pre-defined. For the worst accuracy eV/Å case studied here, it can be obviously seen from the figure that the discrepancy is quite large. Finally, we believe a reasonable thermal conductivity can be obtained for graphene from original IFC3 with force accuracy at least on the order of eV/Å. This critical value of force accuracy for graphene is larger than that for silicon, which might arise from the stronger C-C bonding in graphene than the Si-Si bonding in silicon. The accumulated thermal conductivity as a function of phonon frequency is helpful to understand the contributions of different phonons to the total thermal conductivity. We also show this plot for the benchmark case and eV/Å truncation cases in Fig. 5. It can be seen that the three curves in the figure are smooth and show similar trend. It can be clearly seen that the inaccuracy of thermal conductivity is not due to a few phonon modes, but rather from the entire phonon spectrum. Fig. 4 Thermal conductivity of graphene from optimized Tersoff potential vs. force accuracy Fig. 5 Accumulated thermal conductivity of graphene as a function of phonon frequency Results from truncation cases are used for further discussions. From the Grüneisen parameters plotted in Fig. 6, we find that the result from benchmark IFC3 agrees quite well with that from IFC2, with a relative difference of 2.01%. It is found that when the force accuracy is eV/Å, original IFC3 will yield (3) max 3 T 0.806 e /ÅV ij and modified IFC3 can reduce this value to . The translational invariance conditions are much better satisfied by modified third-order IFCs. Compared with original IFC3, Grüneisen parameters from modified IFC3 also show better agreement with those from IFC2. The relative difference of Grüneisen parameters is 22% between modified IFC3 and IFC2 but much larger than that between original IFC3 and IFC2. For graphene with eV/Å accuracy, modification of third-order IFCs can reduce the error in translation invariance conditions and correct the trend of Grüneisen parameters, but does not yield a thermal conductivity that agrees better with benchmark result. Fig. 6 Grüneisen parameters of graphene from optimized Tersoff potential
Silicene is the counterpart material of graphene with a buckled structure, and therefore it was investigated with the same method as graphene. We used the Stillinger-Weber (SW) potential[63] to describe the Si-Si interactions and the optimized SW2 potential parameters for silicene from Ref. [64] was used. A buckled hexagonal primitive unit cell structure was first constructed for silicene and optimized. And then a supercell was constructed and re-optimized. The lattice constant for the optimized structure was 3.812425 Å and the buckling distance was 0.426948 Å. For silicene, we also calculated the phonon frequencies for 1000 q points from Γ to M with original second-order IFCs. Our structure was well optimized and it could guarantee that no negative frequencies existed near Brillouin zone center even when we used so dense q points, as shown in Fig. 7. The dispersion curves of acoustic modes for silicene show similar trends as those of graphene. It should be noted that the flexural acoustic (FA) mode in silicene is similar to the ZA mode in graphene but not purely out-of-plane.[32,35] It can be seen that LA/TA dispersion curves seem linear and FA dispersion curve looks like quadratic. We further plotted the group velocities of FA branch along Γ to M direction near Brillouin zone center as red circular dots in the inset of Fig. 7. The solid black line in the inset is also a straight line crossing the origin point and the last point in the figure. It can be seen from the inset that6
Silicene is the counterpart material of graphene with a buckled structure, and therefore it was investigated with the same method as graphene. We used the Stillinger-Weber (SW) potential[63] to describe the Si-Si interactions and the optimized SW2 potential parameters for silicene from Ref. [64] was used. A buckled hexagonal primitive unit cell structure was first constructed for silicene and optimized. And then a supercell was constructed and re-optimized. The lattice constant for the optimized structure was 3.812425 Å and the buckling distance was 0.426948 Å. For silicene, we also calculated the phonon frequencies for 1000 q points from Γ to M with original second-order IFCs. Our structure was well optimized and it could guarantee that no negative frequencies existed near Brillouin zone center even when we used so dense q points, as shown in Fig. 7. The dispersion curves of acoustic modes for silicene show similar trends as those of graphene. It should be noted that the flexural acoustic (FA) mode in silicene is similar to the ZA mode in graphene but not purely out-of-plane.[32,35] It can be seen that LA/TA dispersion curves seem linear and FA dispersion curve looks like quadratic. We further plotted the group velocities of FA branch along Γ to M direction near Brillouin zone center as red circular dots in the inset of Fig. 7. The solid black line in the inset is also a straight line crossing the origin point and the last point in the figure. It can be seen from the inset that6 nearly all the red dots sit on the straight line, except that the first two points have a negligible deviation. These deviations are believed to arise from numerical uncertainty that brings very small residual strain. This result shows that the FA phonon for silicene also has a quadratic dispersion curve. Our previous results based on first-principles calculations[32,33,35] showed three linear dispersion curves for acoustic modes. Later Kuang et al. [34] showed that it would be quadratic from first-principles calculations. Actually, from our previous calculation based on SW potential, a quadratic branch was observed.[64] Fig. 7 Phonon dispersion curve of silicene calculated from optimized SW potential for acoustic modes from Γ to M direction around Γ point. Inset: Group velocity of FA mode along Γ to M direction around Γ point. The second-nearest neighbors were considered in computing third-order IFCs. Table 3 shows the maximum absolute values of translational and rotational invariance tensors. These values are comparable to those of silicon and are also very small. Modified third-order IFCs better satisfy rotational and translational invariance conditions than original third-order IFCs.
201 201 1 q mesh was used to calculate the thermal conductivity at 300 K. The two results are exactly the same and we take the benchmark result as 11.94 W/mK. This result is higher than our previous calculation[64] because we used a much denser q mesh compared with the previous
16 16 1 one. Table 3 Translational, rotational invariance conditions and thermal conductivity for silicene from optimized SW potential (benchmark case) Case (2) max T i (eV/Å ) (3) max T ij (eV/Å ) (2), , max(2) I I i i (eV/Å) (3) (3), , max I I ij ij (eV/Å ) Thermal conductivity (W/mK) Original IFC3 eV/Å, all the calculated thermal conductivity values fully agree with the benchmark result. Even the maximum error is less than 0.2%. When the accuracy is eV/Å, the error in thermal conductivity will be larger but modification of third-order IFCs can reduce it. For either truncation case or random addition case, thermal conductivity calculated from modified third-order IFCs is within ±15% of the benchmark result. When the force accuracy is even worse, the error for all the cases becomes very large and can be obviously seen in Fig. 8. Therefore, we believe a reasonable thermal conductivity can be obtained for silicene from modified IFC3 with force accuracy at least on the order of eV/Å. Truncation cases are used for further discussion about the accuracy of third-order IFCs. It is found that when the force accuracy is around eV/Å for silicene, original third-order IFCs will yield translational invariance condition T 0.0224 V Åe / ij and modified third-order IFCs can reduce this value to . Grüneisen parameters are also plotted for IFC2 and modified IFC3 in Fig. 9. The Grüneisen parameters for all the cases shown in the figure seem to agree reasonably. The relative differences of Grüneisen parameters between IFC2 and IFC3 are less than 2% for the three different IFC3 sets. Considering that force accuracy on the order of eV/Å does not yield a thermal conductivity result that agrees well with the benchmark result, we do not recommend using Grüneisen parameters as the criterion for the accuracy of modified third-order IFCs of silicene. Fig. 8 Thermal conductivity of silicene from optimized SW potential vs. force accuracy Fig. 9 Grüneisen parameters of silicene from optimized SW potential
4. SUMMARY AND DISCUSSION
In summary, we used classical potential to investigate the effect of the accuracy of third-order IFCs on the predicted thermal conductivity based on SMRTA. The benchmark thermal conductivity was calculated very carefully until the translational and rotational invariance conditions were well satisfied and Grüneisen parameters from second-order and third-order IFCs agreed quite well. Inaccuracies were introduced to the third-order IFCs by truncating digits or adding random numbers to the input forces. For classical potential, a reasonable thermal conductivity can be obtained from modified IFC3 for silicon/silicene with force accuracy at least on the order of / eV/Å. We can see that calculating the thermal conductivity of silicene requires higher accuracy of forces (and therefore third-order IFCs) than silicon, despite that the Si-Si interaction in silicon is stronger than that in silicene. For graphene, the required force accuracy is eV/Å. Notably, when the force accuracy is eV/Å for graphene, imposing translational invariance conditions to third-order IFCs does not give thermal conductivity that agrees better with benchmark value than using the original third-order IFCs. Regarding the first-principles calculations, it was noticed that large supercell and exceedingly high precision on forces were required to obtain accurate phonon properties.[49] For instance, discretization errors could be introduced from exchange correlation energy to forces and could be as large as 0.001 eV/Å per atom.[49] This would also cause the translational invariance conditions for forces to be destroyed[49] and therefore make the third-order IFCs inaccurate. The error in first-principles calculations may be affected by energy cutoff, k -point grid, reciprocal space projection technique, aliasing errors, discretization errors, and choice of control parameters.[49] When the force error introduced to classical-potential calculations lies in the typical value from first-principles calculations, our calculated thermal conductivity would be quite different from the benchmark value. This can explain the large differences of thermal conductivity values calculated from first-principles for 2D materials in literature. Especially, we have shown that very high accuracy of forces must be obtained in order to get accurate thermal conductivity values. Grüneisen parameters cannot be used as a necessary and sufficient criterion for the accuracy of third-order IFCs in the aspect of predicting thermal conductivity. Finally, unfortunately we are not able to provide more instructive guidance on how to more accurately predict lattice thermal conductivity using the first-principles PBTE method. It is left for future investigation. ACKNOWLEDGEMENT
This work was supported by the National Natural Science Foundation of China (Grant No. 51676121), the Materials Genome Initiative Center of Shanghai Jiao Tong University, and High Performance Computing Center of Shanghai Jiao Tong University. We thank Guangzhao Qin (RWTH Aachen University) for his helpful discussions about the thermal conductivity of black phosphorene and Fuman Xie (Shanghai Jiao Tong University) for her help in writing the code.
APPENDIX A.
Symmetry and invariance conditions for interatomic force constants
For a crystal with point/space group symmetries, the second-order and third-order IFCs must satisfy the following equations[50,51] i j ij
S S , (3) i j k ijk S S S , (4) where S is a point/space group symmetry operation matrix. i , j , and k are the atom indices corresponding to atom i , j , and k after symmetry operation S . These point and space group symmetries can be enforced by calculating an irreducible set of IFCs and then generating the other values according to Eq. (3) and (4), which can also reduce the computational cost. Under the fact that the potential energy would stay unchanged if the structure is translated as a whole under an arbitrary displacement, we can derive the translational invariance conditions ( i.e. acoustic sum rules)[12,15,50,58,65,66] for second-order and third-order IFCs (2) T 0 i ijj , (5) (3) T 0 ij ijkk . (6) Similarly, rotational invariance conditions[50,51] can be derived under the fact that the potential energy would not change if the structure is rotated as a whole under an arbitrary angle, as shown below (2),
I must be symmetric with respect to and i ij jj r , (7) ,(3) I must be symmetric with respect to and ij ijk k ij ijk r , (8) where j r is the Cartesian coordinate in direction of atom j under equilibrium state. denotes the Kronecker delta function. (2) (2) T i i I I and (3) (3) T ij ij I I are implemented in my calculation and they should be equal to if rotational invariance conditions are strictly satisfied. On condition that Eq.s (3)-(8) are not strictly satisfied, the IFCs cannot be physically accurate in the strictest sense. However, numerical error is unavoidable in computing IFCs and these analytical equations can never be strictly satisfied, i.e. satisfied with 0% error. Therefore the maximum absolute value of tensors (2) T , (3) T , (2) 2 T( ) I I , (3) ) T(3 I I are used to examine to what extent translational and rotational invariance conditions are satisfied. B. Grüneisen parameters
Grüneisen parameters can be calculated from second-order IFCs with the following equation[60] (2)
V V , (9) where is the phonon frequency and can be obtained from second-order IFCs. V is the crystal volume. Alternatively, Grüneisen parameters can also be calculated from third-order IFCs[60] (3) , *2 ji ji jijk kijk e e i rm m R q , (10) where e is the phonon eigenvector and m is the atomic mass. j R is the unit cell vector of the unit cell where j -th atom locates. The relative difference between the Grüneisen parameters calculated from second-order IFCs (2) and those calculated from third-order IFCs (3) is defined as (2) (3)(2) 22 e . C. Numerical error in finite difference method
For classical potential, IFCs are computed with finite difference method and the error mainly comes from rounding error and truncation error. Rounding error results from the fact that floating-point numbers are represented in a computer by a finite number of digits of precision.[67] Using double-precision floating-point numbers in C/C++ programming would yield negligible rounding error. Truncation error comes from the finite difference method. The second-order and third-order IFCs can be calculated with[68] (a) central difference method
21 21 ij i jm m C F R mh O hh , (11) mijk i j kn nm C C F R mh R nh O hh , (12) with C C C . (b) finite difference method with fourth-order accuracy O h
22 42 ij i jm m C F R mh O hh , (13) mijk i j kn nm C C F R mh R nh O hh , (14) with
02 1 21
C C C C C . Using fourth-order accuracy finite difference method instead of central difference method can reduce the truncation error but requires calculating more cases for force and thus would lead to more computational cost. References
Ab initio theory of the lattice thermal conductivity in diamond, Phys. Rev. B. 80 (2009) 125203. doi:10.1103/PhysRevB.80.125203. [6] A. Ward, D.A. Broido, Intrinsic phonon relaxation times from first-principles studies of the thermal conductivities of Si and Ge, Phys. Rev. B. 81 (2010) 085205. doi:10.1103/PhysRevB.81.085205. [7] X. Tang, J. Dong, Lattice thermal conductivity of MgO at conditions of Earth’s interior, Proc. Natl. Acad. Sci. 107 (2010) 4539–4543. doi:10.1073/pnas.0907194107. [8] K. Esfarjani, G. Chen, H.T. Stokes, Heat transport in silicon from first-principles calculations, Phys. Rev. B. 84 (2011) 085204. doi:10.1103/PhysRevB.84.085204. [9] J. Shiomi, K. Esfarjani, G. Chen, Thermal conductivity of half-Heusler compounds from first-principles calculations, Phys. Rev. B. 84 (2011) 104302. doi:10.1103/PhysRevB.84.104302. [10] J. Garg, N. Bonini, B. Kozinsky, N. Marzari, Role of Disorder and Anharmonicity in the Thermal Conductivity of Silicon-Germanium Alloys: A First-Principles Study, Phys. Rev. Lett. 106 (2011) 045901. doi:10.1103/PhysRevLett.106.045901. [11] Z. Tian, J. Garg, K. Esfarjani, T. Shiga, J. Shiomi, G. Chen, Phonon conduction in PbSe, PbTe, and PbTe${}_{1\ensuremath{-}x}$Se${}_{x}$ from first-principles calculations, Phys. Rev. B. 85 (2012) 184303. doi:10.1103/PhysRevB.85.184303. [12] W. Li, L. Lindsay, D.A. Broido, D.A. Stewart, N. Mingo, Thermal conductivity of bulk and nanowire Mg 2 Si x Sn 1 − x alloys from first principles, Phys. Rev. B. 86 (2012) 174307. doi:10.1103/PhysRevB.86.174307. [13] L. Lindsay, D.A. Broido, T.L. Reinecke, Thermal Conductivity and Large Isotope Effect in GaN from First Principles, Phys. Rev. Lett. 109 (2012) 095901. doi:10.1103/PhysRevLett.109.095901. [14] T. Shiga, J. Shiomi, J. Ma, O. Delaire, T. Radzynski, A. Lusakowski, K. Esfarjani, G. Chen, Microscopic mechanism of low thermal conductivity in lead telluride, Phys. Rev. B. 85 (2012) 155203. doi:10.1103/PhysRevB.85.155203. [15] L. Lindsay, D.A. Broido, T.L. Reinecke,
Ab initio thermal transport in compound semiconductors, Phys. Rev. B. 87 (2013) 165201. doi:10.1103/PhysRevB.87.165201. [16] L. Lindsay, D.A. Broido, T.L. Reinecke, Phonon-isotope scattering and thermal conductivity in materials with a large isotope effect: A first-principles study, Phys. Rev. B. 88 (2013) 144306. doi:10.1103/PhysRevB.88.144306. [17] H. Dekura, T. Tsuchiya, J. Tsuchiya, \textit{Ab initio} Lattice Thermal Conductivity of ${\mathrm{MgSiO}}_{3}$ Perovskite as Found in Earth’s Lower Mantle, Phys. Rev. Lett. 110 (2013) 025904. doi:10.1103/PhysRevLett.110.025904. [18] J.W.L. Pang, W.J.L. Buyers, A. Chernatynskiy, M.D. Lumsden, B.C. Larson, S.R. Phillpot, Phonon Lifetime Investigation of Anharmonicity and Thermal Conductivity of ${\mathrm{UO}}_{2}$ by Neutron Scattering and Theory, Phys. Rev. Lett. 110 (2013) 157401. doi:10.1103/PhysRevLett.110.157401. [19] T. Luo, J. Garg, J. Shiomi, K. Esfarjani, G. Chen, Gallium arsenide thermal conductivity and optical phonon relaxation times from first-principles calculations, EPL Europhys. Lett. 101 (2013) 16001. doi:10.1209/0295-5075/101/16001. [20] W. Li, N. Mingo, Thermal conductivity of fully filled skutterudites: Role of the filler, Phys. Rev. B. 89 (2014) 184304. doi:10.1103/PhysRevB.89.184304. [21] A.H. Romero, E.K.U. Gross, M.J. Verstraete, O. Hellman, Thermal conductivity in PbTe from first principles, Phys. Rev. B. 91 (2015) 214310. doi:10.1103/PhysRevB.91.214310. [22] J. Garg, N. Bonini, N. Marzari, High Thermal Conductivity in Short-Period Superlattices, Nano Lett. 11 (2011) 5135–5141. doi:10.1021/nl202186y. [23] A. Kundu, N. Mingo, D.A. Broido, D.A. Stewart, Role of light and heavy embedded nanoparticles on the thermal conductivity of SiGe alloys, Phys. Rev. B. 84 (2011) 125426. doi:10.1103/PhysRevB.84.125426. [24] L. Lindsay, D.A. Broido, T.L. Reinecke, First-Principles Determination of Ultrahigh Thermal Conductivity of Boron Arsenide: A Competitor for Diamond?, Phys. Rev. Lett. 111 (2013) 025901. doi:10.1103/PhysRevLett.111.025901. [25] W. Li, N. Mingo, Ultralow lattice thermal conductivity of the fully filled skutterudite YbFe 4 Sb 12 due to the flat avoided-crossing filler modes, Phys. Rev. B. 91 (2015) 144304. doi:10.1103/PhysRevB.91.144304. [26] W. Li, J. Carrete, G.K.H. Madsen, N. Mingo, Influence of the optical-acoustic phonon hybridization on phonon scattering and thermal conductivity, Phys. Rev. B. 93 (2016) 205203. doi:10.1103/PhysRevB.93.205203. [27] A.A. Balandin, Thermal properties of graphene and nanostructured carbon materials, Nat. Mater. 10 (2011) 569–581. doi:10.1038/nmat3064. [28] N. Bonini, J. Garg, N. Marzari, Acoustic Phonon Lifetimes and Thermal Transport in Free-Standing and Strained Graphene, Nano Lett. 12 (2012) 2673–2678. doi:10.1021/nl202694m. [29] L. Paulatto, F. Mauri, M. Lazzeri, Anharmonic properties from a generalized third-order ab initio approach: Theory and applications to graphite and graphene, Phys. Rev. B. 87 (2013) 214303. doi:10.1103/PhysRevB.87.214303. [30] G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, F. Mauri, Thermal Conductivity of Graphene and Graphite: Collective Excitations and Mean Free Paths, Nano Lett. 14 (2014) 6109–6114. doi:10.1021/nl502059f. [31] L. Lindsay, W. Li, J. Carrete, N. Mingo, D.A. Broido, T.L. Reinecke, Phonon thermal transport in strained and unstrained graphene from first principles, Phys. Rev. B. 89 (2014). doi:10.1103/PhysRevB.89.155426. [32] H. Xie, M. Hu, H. Bao, Thermal conductivity of silicene from first-principles, Appl. Phys. Lett. 104 (2014) 131906. doi:10.1063/1.4870586. [33] X. Gu, R. Yang, First-principles prediction of phononic thermal conductivity of silicene: A comparison with graphene, J. Appl. Phys. 117 (2015) 025102. doi:10.1063/1.4905540. [34] Y.D. Kuang, L. Lindsay, S.Q. Shi, G.P. Zheng, Tensile strains give rise to strong size effects for thermal conductivities of silicene, germanene and stanene, Nanoscale. 8 (2016) 3760–3767. doi:10.1039/C5NR08231E. [35] H. Xie, T. Ouyang, É. Germaneau, G. Qin, M. Hu, H. Bao, Large tunability of lattice thermal conductivity of monolayer silicene via mechanical strain, Phys. Rev. B. 93 (2016) 075404. doi:10.1103/PhysRevB.93.075404. [36] L. Zhu, G. Zhang, B. Li, Coexistence of size-dependent and size-independent thermal conductivities in phosphorene, Phys. Rev. B. 90 (2014) 214302. doi:10.1103/PhysRevB.90.214302. [37] A. Jain, A.J.H. McGaughey, Strongly anisotropic in-plane thermal transport in single-layer black phosphorene, Sci. Rep. 5 (2015) 8501. doi:10.1038/srep08501. [38] G. Qin, Q.-B. Yan, Z. Qin, S.-Y. Yue, M. Hu, G. Su, Anisotropic intrinsic lattice thermal conductivity of phosphorene from first principles, Phys Chem Chem Phys. 17 (2015) 4854–4858. doi:10.1039/C4CP04858J. [39] G. Qin, X. Zhang, S.-Y. Yue, Z. Qin, H. Wang, Y. Han, M. Hu, Resonant bonding driven giant phonon anharmonicity and low thermal conductivity of phosphorene, Phys. Rev. B. 94 (2016) 165445. doi:10.1103/PhysRevB.94.165445. [40] J. Zhu, H. Park, J.-Y. Chen, X. Gu, H. Zhang, S. Karthikeyan, N. Wendel, S.A. Campbell, M. Dawber, X. Du, M. Li, J.-P. Wang, R. Yang, X. Wang, Revealing the Origins of 3D Anisotropic Thermal Conductivities of Black Phosphorus, Adv. Electron. Mater. 2 (2016) 1600040. doi:10.1002/aelm.201600040. [41] W. Li, J. Carrete, N. Mingo, Thermal conductivity and phonon linewidths of monolayer MoS2 from first principles, Appl. Phys. Lett. 103 (2013) 253103. doi:10.1063/1.4850995. [42] X. Gu, R. Yang, Phonon transport in single-layer transition metal dichalcogenides: A first-principles study, Appl. Phys. Lett. 105 (2014) 131903. doi:10.1063/1.4896685. [43] R. Yan, J.R. Simpson, S. Bertolazzi, J. Brivio, M. Watson, X. Wu, A. Kis, T. Luo, A.R. Hight Walker, H.G. Xing, Thermal Conductivity of Monolayer Molybdenum Disulfide Obtained from Temperature-Dependent Raman Spectroscopy, ACS Nano. 8 (2014) 986–993. doi:10.1021/nn405826k. [44] J. Carrete, W. Li, L. Lindsay, D.A. Broido, L.J. Gallego, N. Mingo, Physically founded phonon dispersions of few-layer materials and the case of borophene, Mater. Res. Lett. 4 (2016) 204–211. doi:10.1080/21663831.2016.1174163. [45] A. Jain, A.J.H. McGaughey, Effect of exchange–correlation on first-principles-driven lattice thermal conductivity predictions of crystalline silicon, Comput. Mater. Sci. 110 (2015) 115–120. doi:10.1016/j.commatsci.2015.08.014. [46] S. Lee, D. Broido, K. Esfarjani, G. Chen, Hydrodynamic phonon transport in suspended graphene, Nat. Commun. 6 (2015) 6290. doi:10.1038/ncomms7290. [47] A. Cepellotti, G. Fugallo, L. Paulatto, M. Lazzeri, F. Mauri, N. Marzari, Phonon hydrodynamics in two-dimensional materials, Nat. Commun. 6 (2015) 6400. doi:10.1038/ncomms7400. [48] L.N. Virgin, Vibration of axially loaded structures, Cambridge University Press, 2007. [49] G. Kresse, VASP: Accurate force calculations and VASP datasets, (n.d.). http://wolf.ifj.edu.pl/workshop/work2004/files/talks/kresse_vasp_phonon.pdf (accessed January 16, 2017). [50] L. G., L. W., Theory of Anharmonic Effects in Crystals, n.d. [51] K. Esfarjani, H.T. Stokes, Method to extract anharmonic force constants from first principles calculations, Phys. Rev. B. 77 (2008). doi:10.1103/PhysRevB.77.144112. [52] J. Turney, E. Landry, A. McGaughey, C. Amon, Predicting phonon properties and thermal conductivity from anharmonic lattice dynamics calculations and molecular dynamics simulations, Phys. Rev. B. 79 (2009) 064301. doi:10.1103/PhysRevB.79.064301. [53] G.P. Srivastava, The Physics of Phonons, Adam Hilger, Bristol, 1990. [54] W. Li, J. Carrete, N. A. Katcho, N. Mingo, ShengBTE: A solver of the Boltzmann transport equation for phonons, Comput. Phys. Commun. 185 (2014) 1747–1758. doi:10.1016/j.cpc.2014.02.015.6