{'fraction_non_alphanumeric': 0.032067925057890675, 'fraction_numerical': 0.021121324819310925, 'mean_word_length': 4.837190251894327, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
How does the accuracy of interatomic force constants affect the prediction of lattice thermal conductivity

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\"uneisen parameters cannot be used as a necessary and sufficient criterion for the accuracy of third-order IFCs in the aspect of predicting thermal conductivity.


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.  The good agreement between these calculated values and measured experimental data proves its accuracy and reliability. [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][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][29][30][31] silicene, [32][33][34][35] phosphorene, [36][37][38][39][40] MoS 2 , [41][42][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][37][38][39][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 2 , 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 3 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, 4 discretization errors, etc. [49] For example, the attainable fractional precision [49] in forces using VASP is 4 10  . 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.

Single mode relaxation time approximation method
For a periodic crystal structure, the potential energy can be expanded as the Taylor series [50,51] 5 0 ,, 11 ...
where 0 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.
where  denotes different phonons that can be distinguished by wave vector q and phonon branch  .
, ph c  is the volumetric phonon specific heat. For more details about the method we refer the reader to Refs. [52][53][54]. 6 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 3 10  , 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.

Simulation details
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 7 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.

Silicon
Tersoff potential [59] was used to describe the Si-Si interactions in silicon. The face-centered cubic primitive unit cell and 333  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 , T , (see the Appendix for the details). 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 0 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. 8 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.   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 8 10  and 6 10  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 5 10  eV/Å force accuracy case, modified IFC3 9 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 5 10  eV/Å. 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 10 serve as a circumstantial evidence that the third-order IFCs of the benchmark case are quite accurate.
The modified IFC3 for 5 10  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 4 10  eV/Å accuracy case and IFC2. Since 5 10  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
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 5 5 1  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  In order to calculate third-order IFCs, second-nearest neighbors were considered. 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. 5 10  eV/Å case. 12 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.  Figure 4 shows our thermal conductivity result of graphene for different force accuracy. Please note that there is a break in the 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 5 10  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 4 10  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 13 out of the reasonable range we pre-defined. For the worst accuracy 3 10  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 4 10  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 4 10  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. 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 4 10  eV/Å, original IFC3 will 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 4 10  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. 15 Fig. 6 Grüneisen parameters of graphene from optimized Tersoff potential

Silicene
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 5 5 1  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 that 16 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. 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)  gray line without symbols also shows our benchmark result. It can be seen that when the force accuracy is on the order of or smaller than 8 10  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 6 10  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 6 10  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 6 10  eV/Å for silicene, original third-order IFCs will yield translational invariance condition . 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 5 10  eV/Å does not yield a thermal 18 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

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 19 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 5 10  / 6 10  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 4 10  eV/Å. Notably, when the force accuracy is 4 10  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.

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 SS 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) T0 i ij j 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 must be symmetric with respect to and 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 , T ,

B. Grüneisen parameters
Grüneisen parameters can be calculated from second-order IFCs with the following equation [60] (2) 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]   where e is the phonon eigenvector and m is the atomic mass.
with 0 2 1 2  1   1  8  8  1  , , 0, , 12 12 12 12 Using fourth-order accuracy finite difference 23 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.