Methodology Perspective of Computing Thermal Transport in Low-Dimensional Materials and Nanostructures: The Old and the New.

Demands for engineering thermal transport properties are ever increasing for a wide range of modern micro- and nanodevices and materials-based energy technologies. In particular, there is a severe situation due to the rapid progress in the synthesis and processing of materials and devices with structural characteristic length on the nanometer scales, which are comparable or even smaller than the intrinsic length scales (such as mean free path and wavelength) of basic energy carriers (such as phonons, electrons, and photons). Although advanced approaches for controlling the electronic and photonic transport have been proposed in the past decades, progress on controlling lattice vibrations (i.e., the phonons) is still far behind. Gaps between the fundamental understandings of the behavior of the basic energy carriers at small scales and the technological demands still remain, particularly from a computer modeling point of view. Herewith, we give a perspective of the computational approaches for predicting the thermal transport properties of low-dimensional materials and nanostructures, which are mainly sorted into three categories: empirical molecular dynamics, anharmonic lattice dynamics based Boltzmann transport equation, and Landauer theory. The advantage and disadvantage of each method are discussed and some possible solutions are suggested. The discussion is focused on fully and accurately characterizing the mode-level phonon behavior, possible all-order phonon scattering process, and incorporation of realistic nanostructures. Moreover, emerging challenges of phonon coupling effects, such as electron-phonon, phonon-photon, and phonon-magnon coupling, are also discussed. We expect that this perspective will stimulate future research in computer modeling of micro-/nanoscale heat transfer beyond traditional phonons.

Introduction

In recent years, increasing attention has been focused on the thermal transport properties of nanomaterials and nanostructures due to the rapid development of nanotechnology. On the one hand, one expects to have a high thermal conductivity for cooling the nanoelectronics and optoelectronics or for efficient interfacial thermal transport, where the heat can dissipate rapidly and the temperature of the systems can be effectively controlled, and therefore their functional performance can be maintained. On the other hand, a low thermal conductivity is preferred for thermoelectric energy conversion, where the low thermal conductivity is the key point to achieve a high figure of merit (ZT) because ZT is inversely proportional to thermal conductivity. Based on the fact that in most thermoelectrics the phonon mean free path is much longer than that of electrons, one strategy to improve ZT is to reduce the lattice thermal conductivity by nanostructuring. Therefore, developing computational methods to model heat transfer in nanostructures will provide an insight into the understanding of physical aspects, and therefore guide the nanofabrication techniques.

In the past few decades, many computational methods have been developed and applied to predict the thermal transport properties of nanostructures. These methods can be largely grouped into three categories. The first category includes methods based on molecular dynamics (MD) simulations. These include the equilibrium MD (EMD) method based on the Green–Kubo formula,[1] the nonequilibrium MD (NEMD) method,[2−4] directly based on the Fourier’s law, and the approach-to-equilibrium MD (AEMD) method.[5,6] The second category includes methods based on the Boltzmann transport equation (BTE)[7−10] and the lattice dynamics (LD).[11] The third category includes methods based on the Landauer formula or more generally the nonequilibrium Green’s functions (NEGF) formalism.[12,13] For all of the methods, the dynamics of the studied systems can be described by either the empirical potentials or the quantum mechanics. For the BTE and the NEGF methods, mode-by-mode spectral properties can be naturally obtained, whereas only the total thermal conductivity (or conductance) is calculated in the MD-based methods. However, many spectral decomposition methods in the framework of MD simulations have also been recently developed, including the time-domain normal mode analysis,[14] the spectrum energy density method,[15] and the spectrum heat current method or the frequency domain direct decomposed method.[16,17] In the following, we discuss the advantages and challenges for each method.

Molecular Dynamics: Quantum Corrections, Improving the Accuracy and Speed

A crucial advantage of the MD-based methods is that the anharmonicity of the lattice vibration is fully accounted for. Furthermore, MD simulations can also be used to predict the thermal stability of the systems. Despite these advantages, there are a few challenges in these methods, as discussed below. The most fundamental limitation of the MD simulations is the classical distribution of phonons, in terms of correctly capturing the phonon behavior at different temperatures, i.e., the trajectories of atoms are integrated according to the classical Newtonian mechanics and therefore quantum effects are ignored. Due to the classical assumptions, phonons in the MD simulations follow the Boltzmann distribution, i.e., f(ω) = 1/eℏω, instead of the physically correct Bose–Einstein distribution, i.e., f(ω) = (eℏω – 1)−1. As a result, the energy equipartition theorem is applied in the MD simulations, which means that at any temperature T all of the phonons with whatever frequency ω are activated and each phonon mode has the same energy of kBT. Physically speaking, this assumption is only valid when the temperature of the system is much higher than the Debye temperature TD (T ≫ TD), at which the Bose–Einstein distribution can be simplified to the Boltzmann distribution. At temperatures well below the Debye temperature, the high-frequency phonons are at most partially activated rather than fully occupied as in MD simulations. Therefore, MD simulations are usually conducted at relatively high temperatures, in which case, the de Broglie wavelength of the atoms is smaller than the interatomic distance and the atoms can be treated classically. Figure shows an example of the thermal conductivity of bulk silicon, germanium, and diamond. We can easily see that the thermal conductivity calculated using the MD simulations is inaccurate at low temperatures due to the wrong distribution function of phonons. To overcome this fundamental limitation, different quantum corrections in MD simulations have been proposed, such as temperature correction[18,19] and quantum thermal bath.[20] The temperature correction is made by equating the total energies of the classical and quantum systems, which means thatwhere the summations are taken over all of the phonon modes denoted by N wave vectors k⃗ and 3n dispersion branches ν, and EC and EQ are the classical and quantum phonon energies, respectively. By assuming that the quantum energy is harmonic and the classical energy is equipartitioned, the relation between TC and TQ can be determined and the quantum-corrected thermal conductivity can be written aswhere κC and κQ are the classical and quantum thermal conductivities, respectively. However, the above equation is questionable, as it is not derived from the fundamental relations and has not been rigorously proved. Seeking a reliable quantum correction method in the MD method is significantly important and still unsolved.[21,22]

Figure 1

Comparison of thermal conductivity of bulk silicon (a), germanium (b), and diamond (c) between EMD simulations (black dots) and BTE (lines). The red and black line is the BTE result from the single relaxation time approximation and iterative solution, respectively.

The second challenge for the MD-based methods is that the results sensitively depend on the accuracy of the interatomic potential used. For instance, the thermal conductivity of bulk silicon can differ by two times[23,24] when using different potentials to describe the interactions among the Si atoms. Generating accurate potentials, particularly for complex systems, is quite challenging and of great importance for the MD-based methods. Furthermore, transferability of the empirical potentials is also crucial. Most of the potentials are developed for describing specific bulklike materials. Applying them directly to low-dimensional materials and nanostructures is usually questionable. To improve the transferability, one may need to construct different potential models for different situations, and the atomicrex package[25] provides many features to facilitate this task. Alternatively, with the development of computational technologies, these problems may be solved by using ab initio MD (AIMD) simulations, in which there is no requirement of any empirical inputs. Recently, Marcolongo et al.[26] and Stackhouse et al.[27] have demonstrated that this idea is possible to be realized in AIMD. However, more effort is needed to improve this method due to the quite expensive computational resources required by AIMD.

The last challenge for the MD-based methods concerns the space and time scales achievable using the computational resources available currently and in the future. The MD-based methods are linear scaling with respect to the number of atoms in terms of both memory usage and computation time, allowing for relatively large simulation cell sizes. The number of atoms has reached 5 million using the NEMD method[28] and 9 million using the AEMD method,[29] which are indeed impressive. However, it is still computationally challenging to use the MD-based methods to simulate systems with experimentally comparable large-scale structures. On the other hand, the maximum correlation time for observing a converged thermal conductivity in the EMD method can sometimes be very long, and it also requires many (of the order of 100) independent simulations to obtain accurate results. For example, it has been demonstrated that the thermal conductivity of pristine graphene only converges up to a correlation time of about 10 ns,[28] and it even takes about 100 ns for the convergence of the thermal conductivity of ultrathin silicon nanowires (NWs).[30] Therefore, high-performance computing via parallelization is desirable for the MD-based methods. Apart from the conventional parallelization with CPUs, parallel graphics processing units (GPUs) computing provides an emerging way to achieve a high price–performance ratio. An efficient GPU-based MD code called GPUMD[31] has been built from the ground up to improve the state-of-the-art in terms of the performance for simulations with many-body potentials on GPUs by an order of magnitude (Figure ). The current single-GPU version of GPUMD can simulate systems with up to tens of millions of atoms, and we expect that the continuous improvement in the computational power of GPUs and further parallelization using multiple GPUs can bring this limit up to about 1 billion, allowing for modeling many interesting nanostructures.

Figure 2

Computational speed of the GPUMD package and a GPU version of the LAMMPS package using a Tesla K40 GPU, with both double and single precisions. The test system is bulk silicon with 256 000 atoms (using a cubic simulation cell and periodic boundary conditions in all of the three directions) with the many-body Tersoff potential. This figure is adapted from ref (31).

Boltzmann Transport Equation: Considering the Reality of Structures and Higher-Order Scattering Process

When the wavelength (mean free path) is much smaller (larger) than the size of the wave packet, the lattice vibration, which is the main heat carrier in dielectric solids, can be treated as a quantum particle: phonon. The thermal conductivity tensor can then be computed aswhere c(k⃗, ν) is the volumetric phonon specific heat, v⃗g(k⃗, ν) is the phonon group velocity, and τ(k⃗, ν) is the phonon relaxation time. The group velocity can be calculated using LD, whereas the relaxation time can be obtained by solving the BTE. In recent years, the anharmonic LD-based BTE method[7−9] has become one of the most popular methods to calculate thermal conductivity due to its easy input (only the force constants) and accuracy (can be combined with first principles). This method has been proved to be able to accurately calculate the thermal conductivity for most bulk materials[32,33] when the three-phonon scattering process is dominant in the system and the force constants are evaluated directly by the first principles. To calculate the thermal conductivity of nanostructures using the BTE method, one needs to assume that phonon confinement can be ignored in the system and only the relaxation time is affected by the so called “nano effect” by using the Matthiessen’s rule.[10] However, it is well known that the phonon dispersion will be modified considerably when the size of the nanostructure is below 20 nm,[34] which means that the popular BTE method used to calculate the thermal conductivity of nanostructures with a grain size smaller than 20 nm is questionable. One possible way to solve this problem is to calculate the thermal conductivity of the nanostructures directly, e.g., treating the nanostructure as a new system with a new primitive cell or unit cell. In this case, the phonon dispersions could be quite different from the bulk counterpart. Due to the computational complexity, the maximum number of atoms in the system is limited to up to several hundreds of atoms,[30] and only some special systems, such as nanowires[30] or nanotubes,[7] can be treated. Another concern for the BTE method is the higher-order phonon scattering process. Recently, researchers have found that only considering the three-phonon scattering process is far from being enough for some systems at high temperatures.[35] For instance, the thermal conductivity of diamond at 1000 K computed by considering three-phonon process is overestimated by about 20% compared to the EMD result using the same empirical potential (Figure c). Introducing a higher-order phonon process to the BTE method is important for strongly anharmonic systems[35] and systems with large phonon band gaps.[36] However, explicitly calculating the higher-order phonon process is computationally expensive. It will be helpful to develop an efficient approach to evaluate the importance of the higher-order phonon processes without doing expensive calculations.

Finally, the static method (i.e., finite displacement method and density functional perturbation theory) used to calculate the force constants as inputs for solving the BTE is inaccurate or even totally wrong for some systems, such as PbTe[37] and SnSe.[38] For example, one would obtain imaginary frequencies using the static method for SnSe,[38] which means that the system is “unstable”. Such a conclusion is wrong because experiments evidently show that bulk SnSe is stable.[39] Therefore, it is necessary to include the anharmonicity into the force constants. A promising way to solve this problem is to calculate the temperature-dependent force constants using a fitting procedure,[40,41] which has been demonstrated to successfully account for the intrinsic anharmonicity.

Landauer Theory: Beyond Ballistic and Small-Size System

The phonon conduction can be represented on the concept of wave scattering, which is similar to the electron just as Landauer did in 1957. In the Landauer theory, the spectral heat current from the left lead to the right lead through a junction connecting the two leads at two different equilibrium heat-bath temperatures (TL and TR, respectively) is

Then, the spectral thermal conductance G(ω) can be written aswhere T(ω) is the frequency-dependent transmission and A is the cross-sectional area of the system. The transmission can be calculated by using the NEGF method.[12,13,42] The main drawback of this method is that it is computationally challenging to account for phonon anharmonicity, which is particularly important for systems with short phonon mean free paths.[43,44] To include inelastic phonon–phonon scattering, one has to calculate the many-body self-energy and the device Green’s function self-consistently,[45,46] which severely limits the system size that can be treated. An alternative way to solve this problem is to calculate the transmission coefficient from NEMD.[17,47] In NEMD, the heat current spectrum Q(ω) can be obtained viawhere U is the potential and ΔT is the temperature difference between the left lead and the right lead, which are separated by a true or imaginary interface.[47] The angular brackets here represent the time average in MD simulations. Combining eqs and 6, it is easy to know

Figure shows the phonon transmission at K or Cl twin boundary. The results inherently include the anharmonicity due to the intrinsic anharmonicity in the atomic velocity and interatomic potential used. Using this method, the system one can deal with can be as large as ten millions of atoms. In addition, by using the force constants to replace the potential in eq , it is easy to extend this method to ab initio NEMD, which is being done by Hu’s group. Another bottleneck of the Landauer theory is the lack of phonon mode information. Generally speaking, there always exist several different phonon modes at one particular frequency, which cannot be distinguished via the traditional Green’s function method. Although these phonon modes hold the same frequency, they might play different contributions to the thermal transport. As a result, to thoroughly reveal the phonon transport process, we must extend the atomistic Green’s function and obtain the mode-resolved phonon transmission. Such an attempt has recently been done by Latour et al.[48] and Sadasivam et al.[49]

Figure 3

(a) KCl nanotwinned structure. (b) The corresponding transmission coefficient calculated via eq . The black and red line denotes the transmission for K and Cl at the twin boundary, respectively.

Thermal Transport Beyond Traditional Phonons: Coupling Effects

Although phonons are the main heat carrier in most semiconductor structures, other energy carriers, such as electrons, magnons, and photons, can also carry heat especially when there is high carrier density, strong external electric or magnetic field, and light exposure, and their quantum behaviors are not independent but interdependent on each other[50] (Figure ). In addition, for metallic systems, electrons dominate the thermal transport and the key quantity of electron relaxation time is determined by the electron–phonon interaction.[50] In polar semiconductors, such as GaN,[51] the Fröhlich electron–phonon coupling is strong and its influence on phononic thermal transport is of great significance. Generally speaking, all of these quasi-particles can be adopted in BTE because they are either bosons or fermions. For instance, the electrical thermal conductivity can be calculated by the semiclassical BTE under the constant carrier relaxation time assumption.[52] By applying the Fermi’s golden rule, the scattering rate matrix can be expressed, and then, the relaxation time of the charge carriers, magnons and photons, can be computed by solving the BTE. The robust theory and computing algorithms for carriers has been realized recently for bulk materials[53] and interfaces.[54,55] However, such a method is only valid for small cell sizes with no more than 50 atoms, and thus can only deal with some simple systems. Also, to accurately determine the phononic thermal transport in magnetic materials, the phonon–magnon interaction plays a nontrivial role. However, currently, there is no unified theory to quantify it, despite the experiments that have verified its existence.[56] For nanostructures, the low-dimension confinement can induce relatively strong magnetic properties, and then the phonon–magnon interaction would be critical to determine the phonon transport. To quantify the contribution of phonon–magnon interaction to thermal transport, one feasible way is to compute phonons, such as dispersion curve and density of states, under the influence of external magnetic field and then calculate the thermal conductivity within empirical formula. Actually, all of the three approaches mentioned above have the potential to compute the thermal transport in such a system. For the first-principles MD method, the heat current contributed from phonon, photon, and electron can be possibly calculated,[26] and then the thermal transport properties considering such coupling effects can be computed by Green–Kubo EMD or Fourier’s law (NEMD). In the framework of BTE, the scattering rates of phonon, electron–phonon, and photon–phonon for the bulk materials can be calculated first, and the total scattering rate can then be obtained via the Matthiessen’s rule by combining with the boundary scattering. At the same time, under the promise of computational power, the system shown in Figure can be treated as a single unit cell directly, and thus the assumption of boundary scattering can be avoided. In the framework of Landauer theory, similar to the BTE method, one also first needs to calculate the relaxation time due to phonon–phonon, phonon–photon, and phonon–electron scattering, etc., and then the total thermal transport properties of such a system can be obtained under the assumption of boundary scattering. Meanwhile, it is also possible to treat the system as a single unit cell and then calculate the relevant thermal transport properties.

Figure 4

Schematic of phonons, electrons, magnons, photons, and their possible interactions in nanostructures.

Discussion and Conclusions

Although there is still large room to improve the methodologies mentioned above, some of them have been already successfully applied to predict the thermal transport properties of low-dimensional materials and nanostructures. For instance, the BTE method has been widely used to calculate the lattice thermal conductivity of nanowires (NWs) and nanostructures, and the predictions are in good agreement with the experimental measurements.[57,58] MD simulations have promised to provide the correct trend and comparable results with respect to experimental investigations in the polycrystalline nanostructures.[59,60] The Landauer theory is a powerful method to give a valid prediction of the thermal conductance at the interface compared to experimental results. For more details, we recommend reading the review article by Wang et al.[13]

To summarize, in this perspective, we discuss the major computational methods to compute the thermal transport properties of low-dimensional materials and nanostructures. Three types of methods, namely molecular dynamics, anharmonic lattice dynamics based Boltzmann transport equation, and Landauer theory, are introduced. The disadvantages of these methods in dealing with the low-dimensional materials and nanostructures are presented and some possible solutions are given. (1) For classical molecular dynamics simulations, quantum effects and inaccurate potentials are the two main issues. The probable solution for the first difficulty is to develop an effective quantum correction method. The second challenge can be solved via using ab initio molecular dynamics and developing more adequate potentials, such as using machine learning. Meanwhile, the computational efficiency can be improved significantly by using GPU computing. (2) In anharmonic lattice dynamics based BTE, structural reality, higher-order phonon scattering process, and temperature-dependent interatomic force constants are three tremendous challenges. Future work should focus on these three aspects to make this method more robust. (3) Running NEMD instead of the atomistic Green’s function to introduce anharmonicity into the Landauer theory could be an efficient way to enhance the accuracy of thermal transport calculations. Finally, we also give the perspective for studying the coupling effect between traditional phonons and other energy carriers and/or quasi-particles, such as electrons, magnons, and photons, which is expected to be one of the future research direction for thermal or more broad energy transport process in some emerging nanotechnologies.