A size-consistent Grüneisen-quasiharmonic approach for lattice thermal conductivity

Thermal conductivity of a material is of fundamental interest and crucial for the performance of thermal management, heat dissipation, thermoelectrics, and thermal barrier coatings[1, 2, 3, 4, 5, 6, 7, 8, 9]. Determining lattice thermal conductivity $\kappa_{l}$ from a truly first-principles approach requires handling of phonon-phonon interactions[10, 11] that is extremely compute intensive and therefore cannot be done in a routine manner. To treat many systems in a reliable and computationally efficient manner a few feasible methods have been proposed. These include the quasiharmonic Debye model (QDM)[12], the relaxation-time approximation (RTA) method[13], and a modified Debye-Callaway model[14]. Among these the QDM is interesting due to its high computing efficiency. The workhorse of the QDM for calculating $\kappa_{l}$ is the Slack formulae[15] that require two key input parameters, i.e., the acoustic Debye temperature $\theta_{a}$ and the Grüneisen parameters. These two quantities could be deduced approximately by an equation-of-state analysis[16] which requires the total energy as a function of volume that may be easily obtained from only first-principles total energy calculations, thereby avoiding the moderately expensive phonon calculations. Even though the QDM somewhat consistently underestimates $\kappa_{l}$ due to overestimations of Grüneisen parameters, it can still reliably predict the ordinal ranking of $\kappa_{l}$ for several classes of semiconductor materials. It was shown[12] that the QDM is able to deliver high Pearson correlation between the predicted $\kappa_{l}$ and experimental $\kappa_{l}$. In view of the advantages of the QDM, it is interesting to ask if the prediction of $\kappa_{l}$ could be further improved when $\theta_{a}$ and the Grüneisen parameters are obtained from the standard first-principles phonon calculations. One may argue that phonon calculations may be still be expensive for a routine calculation, but they are considered to be very inexpensive compared to a full density-functional theory (DFT) anharmonic treatment[17] that accounts for multiple-phonon scatterings. Indeed, phonon databases[18, 19, 20] have been established that show the feasibility of standard phonon calculation route with today’s computing resources. We note that systematic Grüneisen parameter determinations have also been feasibly carried out for the evaluation of the thermal expansion coefficients within the Grüneisen formalism.[21, 22, 23, 24, 25, 26, 27, 28, 29]

The outline of this paper is as follows. In Section II we review the existing Slack formulae in estimating $\kappa_{l}$. The size consistency issue of the original Slack formulae is then diagnosed and discussed. Subsequently we propose the necessary modifications, using the specific heat argument, to the existing Slack formulae to preserve size consistency. In Section III we present the results obtained using the proposed GQA. Section IV contains conclusions and a summary.

The lattice thermal conductivity $\kappa_{l}(T)$ at temperature $T$ according to the Slack formula [15] is given by

Here $\theta_{a}$ is the acoustic Debye temperature[15] that can be deduced from the Debye temperature $\theta_{D}$ and the number of atoms in a primitive cell $n$ according to

We note there are different phononic properties, one related to the heating and another to heat conduction, that are appropriately described by, respectively, the Debye temperature $\theta_{D}$, and the acoustic Debye temperature $\theta_{a}$ where the optical modes contribute little to the heat conduction. At $T=\theta_{a}$, $k_{l}(\theta_{a})$ is given by[30, 31, 32, 15]

where $\gamma_{a}$ is the macroscopic Grüneisen parameter $\gamma_{m}(T)$ evaluated at $T=\theta_{a}$ (we shall use Eq. 8 to evaluate $\gamma_{m}(T)$). The primitive cell has a mass $M$ and a volume $V$. $h=2\pi\hbar$ and $k_{B}$ are the Planck and Boltzmann constants, respectively.

In this work, we obtain $\theta_{D}=h\nu_{c}^{\rm min}/k_{B}$ by finding the Debye cutoff frequency $\nu_{c}^{\rm min}$ that gives the best match[25] between the constant-volume heat capacity $C_{v}^{D}(\nu_{c},T)$ from the Debye model and the constant-volume heat capacity $C_{v}^{DFT}(T)$ from the DFT. This is achieved through minimizing the absolute error

According to the Debye model[33], the Debye phonon density of states $\rho_{D}(\nu)=A\nu^{2}$ for $0\leq\nu\leq\nu_{c}$ where $A=9n/\nu_{c}^{3}$ to ensure normalization, and $\rho_{D}(\nu)$ vanishes outside the $0\leq\nu\leq\nu_{c}$ range. This allows the calculation of

Here the heat capacity contributed by a phonon mode of frequency $\nu$ is $c(\nu,T)=k_{B}r^{2}/\sinh^{2}r$, $r=h\nu/2k_{B}T$. From the DFT calculations, we obtain

The summation is over all phonon mode indices $\lambda$ and wavevector ${\boldsymbol{q}}$ in the Brillouin zone (BZ). We note that $C_{v}^{DFT}(T)$ may also be determined using the DFT phonon density of states information $\rho_{DFT}(\nu)$ where

where $\nu_{\rm max}$ is the maximum frequency in the phonon spectrum.

The mode Grüneisen parameter is defined as $\gamma_{\lambda{\boldsymbol{q}}}=-(V/\nu_{\lambda{\boldsymbol{q}}})\partial\nu_{\lambda{\boldsymbol{q}}}/\partial V$. With $\gamma_{\lambda{\boldsymbol{q}}}$, the macroscopic Grüneisen parameter[25, 34] as a function of $T$ is

Now we discuss the size-consistency issue. Eq. 3 was originally derived for the monoatomic case, where $n=1$. While the efficacy of the Slack formulae (Eqs. 1, 2, and 3) have been well established[15] when $n=2$, we note that Eq. 2 may not be appropriate for cells with large $n$ since $\theta_{D}$ depends on the atomic mass and bond strength but is independent[15] of $n$. In fact we expect the expression for $\theta_{a}$ is to be independent of $n$ if we assume there are two crystals that are similar physically but differing only in $n$. Take the rocksalt NaCl crystal as an example. We could choose a primitive cell with two atoms ($n=2$) to do a first calculation. However, we could also use the conventional cubic cell with eight atoms ($n=8$) to do a second calculation and we expect $\theta_{a}$ to be the same as that from the first calculation. In Section  III we will show the heat capacities at $T=\theta_{a}$ deduced from Eq. 2 for some wurtzite compounds with $n=4$ is indeed smaller compared with that of the zincblende and rocksalt compounds with $n=2$. Next, because of Eqs. 2 and  3, $\kappa_{l}(\theta_{a})$ is seen to scale as $n^{-\frac{1}{3}}$, which is undesirable since we would expect $\kappa_{l}(\theta_{a})$ to be also independent of $n$. To overcome this size-consistency issue while retaining the good efficacy[15] of the method when $n=2$, we propose to define

with $\Lambda=2^{\frac{1}{3}}$. This choice of $\Lambda$ is made so that we can retain the same numerical results as that obtained with the original Slack formulae in order to make a direct comparison between the GQA (that uses the DFT input parameters) and QDM for the case of $n=2$. The optimal choice of $\Lambda$ may be further investigated in a future study. To be consistent with Eq. 9 we restore the size-consistency of Eq. 3 by using the following modified equation

where $\gamma_{a}^{\prime}=\gamma_{m}(\theta_{a}^{\prime})$. Notice that the exponent of $n$ in the above equation has been changed from $1$ in Eq. 3 to $\frac{4}{3}$ to achieve size consistency. The size consistency of Eq. 10 is seen by observing that both $M$ and $V$ scale linearly with $n$, therefore the appearance of $MV^{1/3}$ in the numerator cancels out the effect of $n^{\frac{4}{3}}$ in the denominator. With this modification, we obtain a modified equation for the lattice thermal conductivity where

In this work, since the calculations for each compound involve many steps, it is important to design a workflow that will automatically generate the structures required for atomic relaxations and phonon calculations. It also needs to handle automatic job submissions and post-processing of DFT results such as scanning of Debye cutoff frequencies to find the Debye temperature.

We use the optimized structures from a phonon database [19] as a starting point to perform the necessary atomic and cell relaxations. Density-functional theory (DFT) calculations are carried out as implemented in the Vienna Ab initio Simulation Package (VASP)[35] with projected augmented-wave pseudopotentials and a PBEsol[36] functional. A relatively high cutoff energy of $600$ eV is used throughout. Atomic relaxation is stopped when the maximum force on all atoms is less than $10^{-3}$ eV/Å. Phonon calculations are performed with a small-displacement method[37, 38]. We have used a supercell of $2\times 2\times 2$ for the cubic based compounds (i.e., diamond, zincblende, and rocksalt) and $3\times 3\times 2$ for the wurtzite compounds.

From the phonon density of states obtained from the DFT calculations we carry out the error minimization as outlined in Section II to find the Debye frequency $\nu_{c}^{\rm min}$ and then determine the Debye temperature $\theta_{D}=h\nu_{c}^{\rm min}/k_{B}$. Note that most numerical values related to the GQA have been tabulated in the Supplementary Material. Fig. 1, shows the Debye temperatures $\theta_{D}$ obtained from the VASP and the ABINIT codes that are very similar even though one is based on the small-displacement method (VASP) and the other on the density-functional perturbation theory (ABINIT)[19]. As expected, the compounds with large $\theta_{D}$ tend to consist of only very light elements.

Fig. 2 shows the results of $\theta_{a}$ from the GQA and QDM[12] (that are deduced from Eq 2), as well as that from the experiment. We note that the estimation of $\theta_{a}$ from the QDM is rather impressive even though it uses only the energy vs. volume information (i.e., without any direct phonon information). We find that the GQA delivers $\theta_{a}$ that is somewhat larger than that from experiment for the zincblende compounds. However, the agreement is better for the rocksalt compounds. The agreement between $\theta_{a}$ from the GQA and experiment is poor for compounds involving very light elements (such as LiH), presumably due to the presence of very high-frequency phonons. We note that better estimates of $\theta_{a}$ may be achieved by a proper modification of Eq 2 for each crystal class. However, since we aim to report the results obtained from the original Slack formulae, we leave the improvement of $\theta_{a}$ estimation in a future study.

The Grüneisen parameters are calculated using a finite central-difference scheme with strains of $\pm 1\%$ and a phonon-band connectivity technique.[39]. Fig. 3 shows the results of the macroscopic Grüneisen parameters $\gamma_{a}=\gamma_{m}(\theta_{a})$ from the GQA, QDM, and experiment. Here we see that the GQA gives a consistent trend as the experiment where diamond, zincblende, and wurtzite compounds have $\gamma_{a}$ mostly ranges between $0.5$ and $1.0$. For the rocksalt compounds, $\gamma_{a}$ is larger than $1.5$. It is noticed that the QDM consistently overestimates the Grüneisen parameters, especially for the diamond and zincblende compounds, and this leads to low $\kappa_{l}$ that is seen in Fig. 4. For the diamond and zincblende compounds, the GQA somewhat overestimates $\kappa_{l}$ which is attributed the fact that $\theta_{a}$ are generally overestimated as shown in Fig. 2. In contrast to these compounds, the rocksalt compounds show better agreement between the GQA and experiment, where $\kappa_{l}$ are generally smaller than the diamond and zincblende compounds due to larger phonon anharmonicities as reflected in larger Grüneisen parameters for the rocksalt compounds (see Fig. 3).

In Section II we have discussed the possible improper estimations of $\theta_{a}$ that is based on Eq. 2, especially when the number of atoms in the primitive cell $n$ is large. In our study all diamond, zincblende, and rocksalt compounds have $n=2$, while all six wurtzite compounds SiC, AlN, GaN, ZnO, BeO, and CdS have $n=4$. Fig. 5 shows the normalized heat capacity $C^{\prime}=C_{v}/3nk_{B}$ evaluated at $\theta_{D}$ fluctuates around a mean of $\sim 0.94$. However, it is seen that $C^{\prime}$ evaluated at $\theta_{a}$ fluctuates around $\sim 0.91$ and exhibits the same variation as that of $C^{\prime}$ at $\theta_{D}$ except for the wurtzite compounds where there is a sudden drop to lower values of $\sim 0.86$. However, $C^{\prime}$ restores its expected variation when it is evaluated using Eq. 9. This shows that $\theta_{a}$ may be more appropriately estimated by the same consistent fractional change of the $\theta_{D}$ as suggested by Eq. 9. We adopt a simple modification that results in Eq. 9 because the original Slack formulae made very good predictions[15] of $\kappa_{l}$ for compounds with $n=2$.

For the wurtzite compounds, Fig. 6 shows the results of the acoustic Debye temperatures $\theta_{a}$, $\theta^{\prime}_{a}$ (top panel), the macroscopic Grüneisen parameters $\gamma_{a}=\gamma_{m}(\theta_{a})$ and $\gamma_{a}^{\prime}=\gamma_{m}(\theta^{\prime}_{a})$ (middle panel), the lattice thermal conductivities $\kappa_{l}$ and $\kappa^{\prime}_{a}$ obtained from the original and modified Slack formulae (bottom panel) at $300$ K. For each of these compounds, it is not surprising that $\theta^{\prime}_{a}$ calculated based on the modified Slack formulae (Eq. 9) is predicted to be larger than that based on the original Slack formulae (Eq. 2). However, these values of $\theta_{a}$ and $\theta_{a}^{\prime}$ make almost the same prediction of the macroscopic Grüneisen parameters. Finally, we find that $\kappa_{l}$ obtained with the modified equation has a closer agreement with the experimental values for SiC, AlN, and BeO. These results show that the prediction of $\kappa_{l}$ is rather sensitive on the way in which $\theta_{a}$ is estimated, which is already noticeable even with a relatively small value of $n=4$. Future studies may include conducting comparative studies with and without the modified Slack formulae in studying compounds with large $n$ such as the special quasirandom structures[40] or half-Heusler crystals[41].

In this paper we have calculated the lattice thermal conductivity $\kappa_{l}$ that is based on the original Slack formulae[15] where various equation inputs such as the Debye temperature $\theta_{D}$, the acoustic Debye temperature $\theta_{a}$, the Grüneisen parameters are determined from the density-functional theory (DFT) phonon calculations. This approach, called the Grüneisen-quasiharmonic approach (GQA), is able to give $\theta_{a}$, Grüneisen parameters, and $\kappa_{l}$ that follow the experimental trends well, despite the empirical nature of the Slack formulae. We have identified a possible size consistency problem of the original Slack formulae that involves the number of atoms in the primitive cell $n$. We then proposed a simple way to rectify it. The modified Slack formulae preserved the high efficacy of $\kappa_{l}$ prediction when $n=2$, while it was shown that the effect of size consistency is noticeable even for $n$ as small as $4$ when we considered the wurtzite structures. It is expected the GQA with the modified Slack formulae could be used as an effective and practical predictor of $\kappa_{l}$ for a large range of materials.

We acknowledge fruitful discussions with D. Ding. We thank the National Supercomputing Center, Singapore (NSCC) and A*STAR Computational Resource Center, Singapore (ACRC) for computing resources. This work is supported by RIE2020 Advanced Manufacturing and Engineering (AME) Programmatic Grant No A1898b0043.