Elastic moduli and thermal conductivity of quantum materials at finite temperature

The viability of materials for various technologies, including electronics, and energy storage, harvesting, and conversion, depends on their mechanical and thermal properties. At the microscopic level, these properties are dictated by lattice vibrations, also called phonons. Specifically, phonons are the main thermal energy carriers in insulators and semiconductors, whereas their interaction with electrons is the primary limiting factor of thermal transport in metals.

Lattice dynamics (LD) is the standard numerical approach to compute the phonon properties both in crystalline materials and in amorphous solids because it naturally includes statistical quantum effects.^1–8 LD relies on the Taylor expansion of the interatomic potential, which can be either semi-empirical or computed by first principles, usually by the density functional theory (DFT). The second-order term in the expansion gives the normal modes and their frequencies, thus giving access to the harmonic properties of the system, such as phonon dispersion relations, group velocities, heat capacity, and elastic constants. Higher-order terms account for anharmonic phonon–phonon scattering events, and they are necessary to calculate the spectral linewidths and the thermal conductivity (⁠ $κ$⁠), which is obtained by solving the Boltzmann transport equation (BTE).⁹ These terms are normally truncated at the third order, but four-phonon processes have been recently found essential to accurately estimate the thermal conductivity of certain materials, such as $BAs$⁠.^10–14 

The second and higher-order terms in the expansion of the potential energy, also identified as interatomic force constants (IFCs) are computed as zero-temperature derivatives of the potential either by finite differences or by the perturbation theory.^15,16 The harmonic IFC can be exploited in the standard quasiharmonic approximation to compute thermodynamic properties.^17,18

However, this implementation of LD is viable only for systems that are stable or metastable at zero temperature, i.e., their harmonic IFC matrix does not exhibit negative eigenvalues which would correspond to imaginary frequencies. Additionally, even for stable systems temperature effects are only accounted for in the solution of the phonon population terms appearing in the solution of the BTE.

The need to overcome these fundamental issues has led to the development of various methods in which IFCs at finite temperatures are inferred from molecular dynamics simulations^19–21 or self-consistent stochastic sampling of atomic displacements.^22–26 

These methods have been implemented in several open-source codes, including SSCHA,²⁷ QSCAILD,²⁸ HiPHive,²⁹ ALAMODE,³⁰ TDEP,³¹ and the “user-phonon” package in LAMMPS.²⁰ 

In this Tutorial Article, we illustrate a state-of-the-art workflow to compute the temperature-renormalized vibrational, elastic, and thermal transport properties of crystals by anharmonic LD at finite temperatures integrating nuclear quantum effects through Feynman’s path integral. The application of this workflow to crystalline silicon shows that even for such a simple system, nuclear quantum effects and finite-temperature displacements must be considered to obtain accurate estimates of the thermal conductivity, especially at low temperatures. To achieve convergence in the sampling of the forces used to compute the IFCs, we exploit a machine-learning neuroevolution potential (NEP), consisting of a single-layer neural feed-forward network model fitted on a DFT database with a natural evolutionary algorithm.³² The forces extracted from path-integral molecular dynamics (PIMD) simulations implemented in the GPUMD code,³³ are used to compute IFCs by a multi-order force constant fitting through the temperature-dependent effective potential (TDEP) method.³¹ Finally, the IFCs are used to compute the phonon properties, elastic constants, and thermal conductivity at each temperature as implemented in $κ$ALDo.⁸ Each step of this workflow takes full advantage of GPU acceleration for enhanced efficiency.

In Sec. II, we illustrate the theory of anharmonic lattice dynamics, the background of the methods used to compute IFCs at finite temperatures, and the calculation of mechanical and thermal properties. In Sec. III, we describe the computational workflow, and report the results for bulk crystalline silicon, with a critical discussion of how they differ from those obtained by standard approaches.

The harmonic approximation, i.e., the truncation of Eq. (2) at the second order, not only provides the phonon dispersion relations $ω( qq)$ but also macroscopic and thermodynamic properties of materials, such as the elastic moduli, thermal expansion, and lattice heat capacity.³⁵ Conversely, to compute the thermal conductivity, it is necessary to consider phonon–phonon scattering mechanisms that are described as a first approximation by the cubic anharmonic term in Eq. (2).

An accurate estimate of IFC2 and IFC3 tensors is essential to characterize the elastic and thermal properties of materials. In Eq. (2), these tensors are defined as the second and third derivatives in the Taylor expansion of the potential energy. This perturbative approach implies that the system is considered in its potential energy minimum at zero temperature. The temperature dependence enters only in the statistical distribution of phonons occupations in Eq. (4) when one computes thermodynamic or transport properties. This approach limits the application of LD to systems that are stable or metastable at zero temperature, thus excluding temperature-stabilized polymorphs, and cannot account for the temperature-induced renormalization of the vibrational frequencies, which can be substantial in strongly anharmonic materials.^36,37 To overcome these limitations, approaches have been developed stemming from the self-consistent harmonic approximation (SCHA).^{4,22,23,27,38,39}

In this class of methods, the normalized IFC2 tensor is represented as the curvature of the variational Gibbs free energy mapped on an effective quadratic Hamiltonian, sampled on a multivariate Gaussian distribution. When the variational minimum is obtained, the same sampling provides higher-order renormalized force constants. An alternative class of methods relies on sampling the canonical distribution through classical molecular dynamics (MD) or Monte Carlo.^19,21 In the temperature-dependent effective potential (TDEP) approach,²¹ the harmonic and higher-order IFCs are fitted over a set of frames from MD trajectories at the desired temperature. Higher-order IFCs are fitted on the residual forces from the previous order. While TDEP was originally developed as an empirical method, it recently found a solid foundation in the Mori–Zwanzig projector formalism of linear response and the mode-coupling theory.⁴⁰ In the same work, Castellano et al. generalize the use of TDEP to path-integral molecular dynamics (PIMD) to obtain a quantum mechanical sampling of the system without assuming multivariate Gaussian distributions. As the theory is derived in Ref. 40, here we outline the practical aspects of the calculation of the IFCs.

It is important to observe that this framework does not require PIMD when nuclear quantum effects are negligible. At a high enough temperature, the sample can be taken from a classical MD calculation, equivalent to a 1-bead PIMD simulation, using classical correlation functions as the classical limit of Kubo correlations.⁴⁴ Alternatively, sampling can be achieved through iterative multivariate Gaussian sampling and fitting in a SCHA framework, which is mathematically equivalent to the SSCHA method.^23,27,45 However, Gaussian sampling can cause artificial biasing toward a more harmonic distribution. In general, while computationally expensive, PIMD sampling offers the most general way to obtain the correct quantum distributions and is the best method to fit the IFCs gauging nuclear quantum effects.^31,40 These effects often have an important role in determining the vibrational properties of materials at low temperatures, with an impact not only on thermal transport but also on superconductivity and, in general, the phase diagram of quantum materials.^45–47 

For fitting IFCs at the quantum level, the advantage of using PIMD over other sampling techniques that assume a Gaussian probability distribution of displacements is that it can naturally account for large deviations from harmonic potentials, as in the case of solid helium.⁴⁰ Additionally, PIMD can be performed at constant pressure by coupling the dynamics with a barostat.⁵⁵ This allows us to obtain the equilibrium density at any temperature and compute the thermal expansion without any further approximation.

The isomorphism represented in Eq. (14) is exact in the limit of $n→∞$⁠. In practical applications, the number of beads needed to achieve convergence depends on the material and the temperature. Materials with light atoms and low temperatures require larger $n$⁠. Convergence in the number of beads must be tested carefully. All PIMD simulations were conducted using GPUMD-3.8.³³ 

The high costs of PIMD sampling and the need to simulate sufficiently large supercells make the direct use of first-principles DFT impractical. To efficiently sample the free energy surface at different temperatures, we employ a machine-learning potential (MLP) trained on a suitable dataset of energies and forces computed by DFT. MLPs have been used extensively to predict vibrational and thermal properties of both crystalline and disordered materials with very promising results since the early applications on phase change and thermoelectric compounds.^56–59 An exhaustive list of thermal transport studies with MLPs is provided in Ref. 60.

Here, we chose the neuroevolution potential (NEP) approach,³² which consists of a feed-forward neural network with 2-body $( d n i)$⁠, 3-body $( d n l i)$⁠, 4-body $( d n l 1 l 2 l 3 i)$⁠, and 5-body descriptors $( d n l 1 l 2 l 3 l 4 i)$⁠. The native implementation of NEP on graphical processing units (GPUs) in the GPUMD code³³ provides a substantial computational advantage, making this approach very efficient for running PIMD with large $n$⁠.⁶¹ 

The dataset used to fit the NEP model contained 2475 structures for crystalline silicon with varying deformations, vacancies, bond numbers, and various polymorphs, including amorphous silicon.⁶² DFT calculations were performed using the Perdew–Wang (PW91) generalized gradient approximation functional.⁶³ 

On the practical side, a QHA calculation consists of computing the zero-temperature energy and phonon frequencies for a grid of volumes around the equilibrium volume and using them to evaluate Eqs. (20) and (23). Since all the calculations are performed at zero temperature, QHA cannot be applied to temperature-stabilized crystals and does not account for the temperature renormalization of phonon frequencies, thus leading to substantial errors for strongly anharmonic systems. In this paper, our QHA simulation utilizes seven volumes (0 $%$⁠, $±0.25%$⁠, $±0.5%$⁠, $±1.0%$⁠.) with each QHA calculation performed on a $4×4×4$ supercell using the 8-atom conventional cell and q-point mesh of $30×30×30$⁠.

As such, accurate modeling of the elastic constants of a material is necessary for understanding how the material responds to strain and material strength.

The standard implementation of anharmonic lattice dynamics to compute the vibrational properties and the thermal conductivity of crystals relies on the calculation of IFCs at zero temperature, either by finite differences⁶⁷ or by perturbation theory.^15,16 The details of this workflow are given in Ref. 5. Solutions of the linearized phonon BTE are implemented in several open-source codes, including Phono3py,⁸⁴ ShengBTE,⁸⁵ TDEP,⁸⁶ D3Q-Thermal2,^16,79 almaBTE,⁸⁷ ALAMODE,³⁰ and $κ$ALDo.⁸ 

As for the solution of the BTE, Phono3py, almaBTE, and $κ$ALDo implement the direct inversion or Eq. (40), ShengBTE and TDEP use the iterative self-consistent approach, and Thermal2 the variational method.

All these codes have the option to read or calculate finite differences IFCs, except D3Q-Thermal2, which is fully integrated with the Quantum-Espresso suite,⁸⁸ and exploits density functional perturbation theory (DFPT) to compute both second and third-order IFCs. ShengBTE and AlmaBTE can read second-order IFCs computed by both finite differences and DFPT, while third-order IFCs are normally computed by finite differences. $κ$ALDo features full compatibility with Quantum-Espresso and D3Q force constants, finite-difference calculators, and interfaces with several force calculators, including LAMMPS and TDEP.

The procedure proposed in this article for calculating the elastic and thermal properties of materials at finite temperatures is outlined in the flow chart in Fig. 1. The proposed procedure consists of the following steps:

1. Identify or fit from scratch a suitable MLP for the system of interest, in our case, bulk silicon. Here, we use a pre-trained NEP that was formerly shown to reproduce the vibrational and thermal properties of bulk crystalline silicon.³² If necessary, NEP fitting can be performed with GPUMD.³³ 

2. Generate a sufficiently large supercell replicating the unit cell of the crystal of interest and run PIMD. We suggest running a shorter constant pressure simulation for equilibration and eventually a longer production run at constant volume for sampling displacements and forces in the canonical ensemble.

3. Constant pressure PIMD simulations give the equilibrium lattice parameter at each temperature and allow one to infer directly the thermal expansion coefficients.

4. Fit second- and third-order IFCs using TDEP.³¹ 

5. Calculate the elastic constants from the second-order IFCs.

6. Calculate the thermal conductivity in the ALD-BTE framework. For these last two steps, we decided to use $κ$ALDo, which is our in-house lattice dynamics code.

Each step of the workflow requires careful evaluations of the numerical results and convergence tests. The following parameters need to be considered: accuracy of the MLP, number of the PIMD beads, size of the supercell and spacial range of the TDEP fit, PIMD sampling time, Brillouin zone sampling (number of q-points) for the calculations of phonon properties, elastic constants, and thermal conductivity. Some of these checks may be performed for a single step, but for others, it is necessary to complete the full workflow to verify convergence.

Here, we assume that the NEP model is sufficiently accurate as it was developed and tested in a previous work.³² As the lowest temperature exhibits the largest nuclear quantum effects, we test the convergence PIMD total energy with the number of beads, $n$⁠, at 50 K, and we can safely assume that $n=24$ is sufficiently well-converged at 50 K and higher temperatures (Fig. 2).

To compute the IFCs by TDEP, we need to run PIMD on a supercell, i.e., a replica of the unit cell. The second-order IFCs are computed up to a cutoff radius which can be at most the largest sphere that can be included in the supercell. Third-order cutoffs were chosen such that the coefficient of determination is maximized, indicating the best fit. To obtain the largest cutoff radius with the smallest number of atoms it is convenient to use supercells as close as possible to cubic. In the case of bulk silicon, we can refer to the 8-atom conventional cubic cell as the unit cell and use $n r e p× n r e p× n r e p$ replicas of it. The convergence of $κ$ with $n r e p$ needs to be tested, which unfortunately implies running the full workflow in Fig. 1. IFCs are computed up to the maximum radius available within the supercell, corresponding to half the length of the supercell edge. $κ$ as a function of the size of the supercell at $T=300$ K is shown in Fig. 3, with error bars computed by block averages over five blocks of the PIMD trajectories. $n r e p=3$ already provides a reasonable estimate of $κ$⁠, which can be well converged for $n r e p=4$ as the value is well within one standard deviation compared to $n r e p=5$⁠. We then use $n r e p=4$ for the calculations at the other temperatures. The second-order IFC cutoff corresponds to twice the lattice parameter (⁠ $∼11$ Å), while the optimal third-order IFC cutoff is approximately 6 Å. We note that $κ$ for silicon is not very sensitive to the range of the IFCs and the size of the supercell, as also observed in previous works.⁵ This would be different for polar materials where long-range interactions have a large impact on the thermal conductivity.^89–91 

Finally, we need to test the convergence of the thermal conductivity against the density of the q-point mesh used to sample the first Brillouin zone. We conducted this test at the lowest temperature considered, 50 K, as converging $κ$ at low temperature is usually more challenging. Figure 4 shows that $κ$ plateaus for a $11×11×11$ q-point grid. However, we choose to use a $12×12×12$ grid since even grids are often computationally preferable for cubic systems. Note that this test was also conducted using the conventional 8-atom cell. The $12×12×12$ grid corresponds to a $19×19×19$ grid for the 2-atom primitive cell, which in our implementation of the BTE is well converged.⁸ 

In the standard approach to computing $κ$ by first principles BTE, thermal expansion is not considered.³ In principle, one should calculate the lattice parameter and then the IFCs at each temperature to obtain a more accurate estimate of $κ$⁠. As described in Sec. II E, the QHA approach allows one to calculate the temperature-dependent lattice parameters at the cost of a few more phonon calculations. Here, we compare the accuracy of QHA to the direct calculation of the lattice parameter of bulk silicon by constant pressure PIMD. At low temperatures, in PIMD simulations, it is important to converge very carefully the lattice parameter with the number of ring-polymer beads [ $n$ in Eq. (14)]. To this scope, at temperatures below 200 K, we have computed the equilibrium lattice parameter as a function of $1/n$ up to $n=64$⁠, extrapolating for $n→∞$⁠.

Figure 5 provides a comparison of these two approaches with the experimental measurements.⁹² Theoretical results suffer from the systematic overestimate of the lattice parameter of the underlying DFT functional used to fit the NEP model, a well-known shortcoming of GGA.⁹³ However, it is meaningful to compare the relative thermal expansion correcting for this systematic shift. The temperature dependence of the lattice parameter agrees with former quantum-mechanical calculations and experiments.^92,94–96 In particular, at temperatures below 100 K, both QHA and PIMD calculations support previous evidence of negative thermal expansion, as a direct consequence of nuclear quantum effects.⁹⁴ The $T→0$ limit of the lattice parameter for QHA and PIMD, obtained by minimizing the free energy including quantum zero-point motion, is about 0.01 Å larger than the lattice parameter that minimizes the total energy of the system. This correction corresponds to 0.2 $%$ of the lattice parameter, in line with previous calculations.⁹⁷ QHA provides estimates of the lattice parameter that are qualitatively consistent with PIMD up to room temperature. QHA and PIMD results start to diverge at temperatures larger than 600 K due to the increasing impact of anharmonicity. Figure 5 shows that overall PIMD simulations provide a better agreement with experiments than QHA. It is worth pointing out that above the Debye temperature (645 K for silicon) one can use classical MD, which is much more computationally efficient than PIMD.

To calculate the linear thermal expansion coefficient (⁠ $α L$⁠), we have interpolated the calculated temperature-dependent lattice parameters with a cubic spline and computed the derivative at room temperature. Compared to the experimental value of $2.72× 10 − 6$ K $− 1$⁠, we obtain a larger $α L$ by PIMD (⁠ $4.03× 10 − 6$ K $− 1$⁠) and even larger with QHA (⁠ $4.59× 10 − 6$  K $− 1$⁠).^95,96 This indicates that constant pressure PIMD (or MD) provides a better estimate of $α L$ than QHA, but the main source of discrepancies with experiments comes from the machine-learning potential or, most likely from the parent DFT functional used to construct the training set.

As elastic constants can be directly inferred from second-order IFCs using the Born long-wave method, the proposed PIMD+TDEP scheme provides a viable tool to compute the elastic properties of solids at finite temperatures alternative to molecular dynamics or Monte Carlo^72,99,100 and potentially more efficient. The elastic constants reported in Table I are in reasonable agreement with the experiments overall. The main reason for the observed discrepancies with experiments is the approximations in the DFT functional, but there are also differences between the DFT reference and the NEP model, size convergence issues inherent to the Born long-wave method, and statistical uncertainties in the TDEP fit. We notice that the reference DFT functional⁶² systematically underestimates the elastic constants compared to experiments (see first two rows of Table I), which is not uncommon for GGA functionals. To compare DFT and NEP, we computed the elastic constants at zero temperature by directly fitting the stress/strain relation through calculations of the stress tensor of strained supercells.³³ The NEP model further underestimates $c 11$ and $c 44$ but gives a slightly higher $c 12$⁠. These calculations also serve as a reference to check the accuracy and the convergence of the Born long-wave method. Computing the second-order IFC matrix by finite differences at zero temperature, we find that elastic constants converge for a cutoff radius corresponding to half the edge of a $5×5×5$ supercell. Finite temperature calculations are performed with this supercell size.

The finite temperature elastic constants exhibit the expected dependence on the temperature, as the material gets mechanically softer as the temperature increases. Softening is more pronounced for the diagonal components of the elastic tensor (⁠ $c 11$ and $c 44$⁠), whereas the off diagonal $c 12$ has a weaker temperature dependence. The decrease of the elastic constants becomes more prominent above room temperature and the trends are reflected in the bulk, Young, and shear moduli also reported in the table. It is worth noting that there is a crossover between $c 12$ and $c 44$ at about 600 K, as the latter becomes smaller than the former. This determines a change of sign of the Cauchy pressure, defined as $c 12− c 44$ which turns from negative to positive. A negative Cauchy pressure is characteristic of covalent directional bonding, whereas a positive value is typical of metallic bonding. This also marks a transition from brittle to ductile that takes place when the Pugh ratio $B/G$ exceeds 1.75.¹⁰¹ In our calculations, this happens at 900 K. However, whereas this trend is physically meaningful, the NEP model’s error on the shear modulus is larger than that over the bulk modulus, thus leading to an overestimate of the ductility of the material. More quantitative estimates of the elastic moduli would require better MLPs, possibly trained on datasets with a higher level of theory.

Anharmonicity in the potential causes a temperature renormalization of phonon frequency, which affects the thermal properties of materials and stands at the base of experimental techniques such as Raman thermometry. Our workflow allows the calculation of phonon dispersions at finite temperatures (Fig. 6). Temperature renormalization leads to a monotonic shift toward lower frequencies and mostly affects optical modes. This shift can be substantial, exceeding 1 THz at the highest temperature (1200 K).

Temperature-dependent Raman measurements allow us to assess the accuracy of our calculations. In Fig. 7, we compare the temperature dependence of the Raman peak to the frequency of the longitudinal optical mode at the center of the Brillouin zone. The absolute difference between the two data sets stems from approximations in the DFT functional and further inaccuracy in the NEP model. However, the temperature trends are in good agreement, with a slightly stronger temperature dependence in the calculations at temperatures higher than 300 K. It is worth noting that, in the framework of TDEP with MD, the eigenvalues of the mass-normalized second-order IFC matrix correspond to the finite-temperature frequency spectrum.⁴⁰ This is not exactly the case when sampling a Gaussian distribution as in the sTDEP or SCHA approaches.^40,103

Temperature renormalization of acoustic modes is more subtle but may have a significant impact on the mechanical and thermal properties of materials. For example, the decrease of the elastic constants discussed in the previous section is strictly related to the softening of acoustic modes with temperature. This effect can be quantified by computing the group velocities of the three acoustic branches at the $Γ$-point corresponding to the speed of sound in the material. Table II shows that the group velocities of both longitudinal and transverse acoustic modes are considerably reduced at high temperatures (900 and 1200 K).

The main goal of this study is to develop and assess a workflow to calculate thermal conductivity including temperature effects more accurately than in standard anharmonic lattice dynamics approaches. $κ(T)$ computed by inverting the BTE with second and third-order IFCs obtained from the PIMD+TDEP workflow is plotted in Fig. 8 alongside experimental measurements.^104–106 Our results are compared to calculations performed with IFCs obtained by finite differences (FD) at zero temperature either with the lattice parameter corresponding to the minimum of the equation of state or considering thermal expansion in the QHA. Figure 8(a) shows that PIMD+TDEP results are closer to the experimental reference, especially at low temperatures. The percent difference between calculations and experiments is reported in Fig. 8(b). We notice that sampling the IFCs at 50 K with classical MD in the isobaric-isothermal ensemble provides an estimate of $κ$ within the statistical uncertainty of that obtained using PIMD. This means that including thermal motion effects in the calculation of both second and third-order IFCs substantially improves the accuracy of low-temperature calculations, but nuclear quantum effects have on the IFCs have a small impact on the determination of $κ$ for silicon, whereas they are essential to describe thermal expansion. We expect nuclear quantum effects to have a stronger impact on materials with light atoms. Conversely, for silicon, the renormalization of the IFCs at temperatures above 600 K has a relatively smaller effect. The correction produced by computing IFCs with PIMD+TDEP compared to FD [Fig. 8(c)] is positive, except for $T=1200$ K and the largest at the lowest temperatures. The QHA correction, in turn, is always negative and remains similar at all temperatures.

To analyze the source of the differences between PIMD+TDEP and FD results across the broad temperature range considered, we computed anharmonic phonon properties, including phonon linewidths (which are the inverse of phonon lifetimes), and the cumulative thermal conductivity as a function of phonon mean free paths at 50, 300, and 1200 K (Fig. 9). Phonon linewidths provide a measure of the anharmonicity of vibrational states, and mean free paths, together with phonon group velocities, provide phonon-resolved contributions to the total lattice thermal conductivity, according to Eq. (39). In general, we observe that the linewidths of low-frequency acoustic modes computed from PIMD+TDEP are significantly smaller than those obtained from FD. Up to 300 K, we do not see significant temperature effects on phonon group velocities, hence we can conclude that the large thermal conductivity corrections at low temperatures stem from nuclear quantum effects on the anharmonic phonon properties, i.e., third-order IFCs. Indeed, we observe substantially longer phonon mean free paths providing a large contribution to $κ$ at 50 K. At room temperature, PIMD+TDEP and FD give similar $κ$⁠. The slightly larger $κ$ obtained by PIMD+TDEP is due to a slightly lower anharmonicity of the low-frequency acoustic modes that exhibit longer lifetimes leading to longer mean free paths than those predicted by FD. At 1200 K, $κ$ from TDEP+PIMD is lower than from FD. This is mostly due to the renormalization of phonon group velocities as anharmonic effects make acoustic modes softer, whereas lifetimes computed from temperature-renormalized IFCs still result slightly longer than from FD. It should be stressed that four-phonon and higher-order scattering processes become more impactful at high temperatures and would further reduce the thermal conductivity, potentially bringing it more in line with the experiment.¹⁰⁷ 

Silicon is usually considered a simple benchmark system for thermal conductivity calculations as it is made of relatively heavy atoms and does not exhibit temperature-stabilized phases. In fact, the standard FD approach usually provides reasonable results, especially around room temperature.^3,5 Instead, we found that accounting for nuclear quantum effects on third-order IFCs is essential to compute $κ$ at low temperatures and provides a substantial correction to $κ$ even at room temperature. This is probably the case for several other crystalline solids with high Debye temperatures that are weakly anharmonic and commonly considered suitable for treating with a perturbative FD approach. Generally, we see that the temperature renormalization of third-order IFCs produces a positive correction of $κ$⁠. As expected, at high temperatures, above 600 K, frequency renormalization plays a role and results in a decrease of phonon group velocities, thus providing a negative correction to $κ$⁠. In our silicon model, the former effect is predominant up to 600 K, whereas the latter becomes predominant at 1200 K.

In summary, we introduced a workflow that uses PIMD and TDEP to account for nuclear quantum effects and temperature renormalization of the interatomic force constants, providing a more accurate way to compute the elastic and thermal properties of materials at finite temperatures. We have illustrated the theory and the step-by-step procedure, including convergence tests and finite-size scaling, to apply this approach to a benchmark crystalline system, bulk silicon. The results highlight the importance of considering nuclear quantum effects in the calculation of anharmonic force constants to obtain an accurate estimate of the thermal conductivity at low temperatures. These modern numerical methods, when combined, achieve closer alignment with experimental results providing a basis for studying more novel or poorly understood materials that exhibit temperature-induced phase transitions, strong anharmonicity, or substantial nuclear quantum effects.

We are grateful to Olle Hellman and Matthieu J. Verstraete for fruitful discussions on theories and implementations of temperature-dependent effective potentials. This work was partially supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering (Grant No. DE-SC0022288). Florian Knoop acknowledges support from the Swedish Research Council (VR) Program No. 2020-04630 and the Swedish e-Science Research Centre (SeRC).

The authors have no conflicts to disclose.

Dylan A. Folkner: Investigation (lead); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (equal). Zekun Chen: Data curation (equal); Software (equal); Writing – review & editing (supporting). Giuseppe Barbalinardo: Software (lead); Writing – review & editing (supporting). Florian Knoop: Conceptualization (equal); Methodology (supporting); Software (equal); Writing – review & editing (supporting). Davide Donadio: Conceptualization (lead); Methodology (equal); Project administration (lead); Resources (lead); Writing – original draft (supporting); Writing – review & editing (lead).

The complete input files of PIMD simulations, TDEP force constant fitting, BTE calculations with $κ$ALDo, alongside the supplement data are publicly available in GitHub Repository at https://github.com/dafolkner//Silicon_project, Ref. 108.