Non-Equilibrium Thermal Transport: A Review of Applications and Simulation Approaches

Zexi Lu; Xiulin Ruan

ES Energy & Environment

ISSN: 2578-0646 (Print Version) ISSN: 2578-0654 (Online Version)

Review Paper

Citation: ES Energy & Environment, 2019, 4, 5-14

DOI:10.30919/esee8c155

Received: 22 Jan 2019
Revised: 22 Feb 2019
Accepted: 01 Apr 2019
Published online: 02 Apr 2019

View PDF Version

View volume articles

Non-Equilibrium Thermal Transport: A Review of Applications and Simulation Approaches

Zexi Lu and Xiulin Ruan

School of Mechanical Engineering, Purdue University, West Lafayette, Indiana, 47907 USA

*E-mail: lu267@purdue.edu (Z. Lu); ruan@purdue.edu (X. Ruan)

Abstract

As modern electronics shrink into nano-scale with increasing power consumption, thermal management inevitably becomes an important issue. Contrary to the conventional view of a single “apparent” temperature profile, electrons and different modes of phonons usually have different temperatures under many circumstances. This local thermal non-equilibrium appears widely in applications such as laser-matter interaction, hot electron cooling in solar cells, thermal transport across metal-insulator and semiconductor-semiconductor interfaces, etc. The resolved temperature profiles usually deviate significantly from the apparent temperature and have non-negligible effects on the measurement and interpretation of thermal properties and performance. Regardless, this phenomenon has been largely overlooked in early studies and has only started to attract attentions in the recent decades. In this work, we will review the background and research efforts on non-equilibrium thermal transport. The applications where non-equilibrium thermal transport is significant, and the corresponding simulation approaches, focused on two-temperature model, molecular dynamics, Boltzmann transport equations, and multi-temperature model, are reviewed.

Table of Content

We for the first time present a comprehensive review of applications and simulation approaches on non-equilibrium thermal transport.

Keywords: Non-equilibrium thermal transport; electron-phonon coupling; spectral phonon temperature

1. Introduction

Nowadays as electronics shrink into nano-scale and operate with increasing volumetric power consumption, thermal management in devices becomes an important issue in advancing the technology.¹ The thermal transport processes in metallic systems and across interfaces draw researchers’ attention as they appear widely in modern electronics, such as central processing unit in computers, hard drives, heat-assisted magnetic recording (HAMR) devices, laser diodes and thermoelectrics.^2–10 A very important feature in such systems is the coupled thermal transport among electrons and different phonons. In many applications, such as laser-matter interaction and hot electron cooling, electrons and phonons can be driven into strong non-equilibrium due to selective electron-phonon (e-p) coupling and phonon-phonon (p-p) coupling, and their respective temperature profiles significantly deviate from each other.^11–14 This phenomenon can significantly affect the heat transfer processes such as hot electron relaxation and thermal measurements in experiments, and has been drawing attention among the research community in both experimental and theoretical areas in the recent decade.

In thermodynamics, temperature is defined at locations which are in local thermal equilibrium. Phonons, which are bosons and the particle representation of lattice vibration, follow the Bose-Einstein distribution:¹⁵

where ω is the phonon frequency, k_Bis the Boltzmann constant and T is the temperature. Electrons, which are fermions, follow the Fermi-Dirac distribution:

where ∈ is the electron energy state. A local thermal equilibrium means there is a common temperature T that can be applied to Eqs. ([eq:BE_dist]) & ([eq:FD_dist]), which is also the apparent temperature we can observe through measurements. Many classical thermal theories and techniques, such as the acoustic mismatch model (AMM), diffuse mismatch model (DMM) and heat diffusion equations, are based on the assumption of local thermal equilibrium, and have been successful in analyzing and predicting the thermal transport process in many materials.¹⁶

Non-equilibrium thermal transport, on the other hand, involves thermal states not in local equilibrium, which cannot be described by Eq. ([eq:BE_dist]) or ([eq:FD_dist]). All materials have multiple phonon branches. Most metals with a face-centered cubic (FCC) unit cell have 3 phonon branches, while many semiconductors such as Si, Ge and semi-metals such as graphene and graphite have 6 phonon branches, and compounds can have 9 or more phonon branches. Phonons belonging to different branches have different vibration patterns and physical properties. As a result different phonon modes and electrons in the material can be driven into non-equilibrium under certain conditions. It has been shown that in many applications, such as laser-matter interactions, heat transfer across metal-dielectric interfaces, high-energy irradiation in nuclear reactors or space and hot electron relaxation in semiconductors and solar cells, local thermal non-equilibrium is observed to dominate the physical process.^17–22 For example, due to the different e-p coupling strength of different branches, during an ultrafast laser heating process, the amounts of energy received by different phonons from the electrons vary and they will not have a common temperature, as is shown in Fig. 1. Under such conditions, electrons and different phonon branches have different temperatures and Eqs. ([eq:FD_dist]) & ([eq:BE_dist]) have to be applied to them separately. If an analysis based on local thermal equilibrium is still applied, one may get inaccurate or even wrong results. For example, in a device under non-equilibrium thermal transport, some of the phonons can be much hotter than the others, leading to non-thermal failure although the average lattice temperature is still within the damage threshold.^23,24 Therefore a fundamental understanding of the non-equilibrium thermal transport is necessary for accurately predicting the physics and performance of devices, thus providing guidance in designing such nano-structures.

Illustration of the ultrafast laser heating and hot electron cooling processes in aluminum in Ref. . The temperature profiles are predicted from calculations with fitted parameters matching the experimental data. The 3 different phonon branches have different temperature due to their different coupling strength to the electrons.

Fig. 1 Illustration of the ultrafast laser heating and hot electron cooling processes in aluminum in Ref. 11. The temperature profiles are predicted from calculations with fitted parameters matching the experimental data. The 3 different phonon branches have different temperature due to their different coupling strength to the electrons.

While experiments are a solid approach to investigate these phenomena, they also have certain limitations. The direct measurement of the temperatures of electrons and individual phonon branches can be extremely difficult, hindering us from acquiring a complete picture of the non-equilibrium thermal transport from experiments alone. Therefore theoretical works are also needed for complement. In this article, we will present a review of simulation approaches on non-equilibrium thermal transport in different applications. This article is organized as follows: applications where non-equilibrium thermal transport dominates are divided into four categories along with their commonly-used simulation approaches: 1) e-p and 2) p-p non-equilibrium in a single material, and 3) e-p and 4) p-p non-equilibrium at interfaces between two different materials. It is noteworthy that during ballistic energy transport, local thermal equilibrium is not established either.^25–29 When the system size is so small that it becomes comparable with the energy carrier’s mean free path (MFP), electrons and phonons travel almost without scattering with each other. As a result, temperature cannot be defined using Eq. ([eq:BE_dist]) & ([eq:FD_dist]). However, this is a case where local thermal equilibrium cannot be established, which is different from our discussion of non-equilibrium thermal transport. Therefore it will not be discussed in our work.

Non-equilibrium thermal transport: applications and simulation approaches

Electron-phonon non-equilibrium in metals

The observation of thermal local non-equilibrium dates back to decades ago during experiments on laser excitation of crystalline metals, where electrons and phonons were found to be at different temperatures and couple to each other^30–33 Electrons are first heated up by the incident photons and subsequently transfer heat to the lattice. As the laser technique advanced, more time-resolved observation of the relaxation process at femtosecond-scales have been presented,^34–37 and efforts have been delivered to the measurement of e-p coupling.^38–40 Similar phenomena have also been observed in metals and their alloys under other forms of irradiation, e.g. ion irradiation in nuclear reactors and space.^41–47 When the nuclear fuel and its coating go through an ion cascade, the incident ions strike into the materials. Unlike laser irradiation, incident ions not only heat up electrons due to the electronic stopping effect, part of their kinetic energy is also transferred to the lattice which at the same time produce point defects in the structure. Due to the large difference in heat capacities between electrons and the lattice, a temperature difference is induced between electrons and phonons, and the subsequent defect evolution in the material is significantly affected by the e-p coupling.^48–50

While the direct measurement of electron’s temperature is difficult, researchers rely on calculations and simulations to predict and interpret experimental results. Theoretical models have been actively developed. The relaxation between excited electrons and the crystalline lattice used to be described by a classical model, where they follow Fermi and and Bose distributions with their own temperatures respectively.³⁰ With the application of the e-p coupling factor G_ep, the calculation approach was able to bypass going into the detailed scattering of all the particles in the system,⁵¹ and used the following simplified equations to describe the relaxation process:

where e and p refer to electron and phonon respectively. T is the temperature while is time. C refers to the volumetric heat capacity of the energy carriers. A simple mathematical derivation can straightforwardly show that the temperature evolution follow exponential curves. The reliability of Eq. [eq:P_TTM_GE] strongly depends on the accuracy of G_ep , which varies significantly with the electronic temperature.^52,53 Theoretical approaches, especially ab initio calculations, have been developed to predict its value.⁵⁴

Although Eq. (3) has been successful in predicting the results from pump-probe and time-domain thermal reflectance experiments,^1,55–57 practical calculation methods that can resolve the spatial temperature distribution had been absent for decades. In principle, the theoretical prediction could be achieved through rigorously solving the Boltzmann transport equations (BTE), but the cumbersome calculations hinder its applications. It was not until the introduction of the Fourier two-temperature model (TTM) by Majumdar and Reddy that a pratical method which can do calculations in both time and space domain was presented.⁵⁸

Two-temperature model

In metals, electrons contribute to most of the thermal transport over phonons. The two types of energy carriers can easily be driven into non-equilibrium through their different interaction with external excitations or different boundary conditions. In TTM, electrons and phonons are described as two interacting subsystems with their own temperatures.⁵⁸The governing equations are:

here k refers to the thermal conductivity of the energy carriers. Compared with Eq. (3), Eq. (4) adds the particle conduction term in the space domain and allows us to calculate the analytic solution for the spatial temperature profiles. It hence became a popular tool in modeling the two-channel thermal transport in metals, especially across metal-dielectric interfaces, which will be introduced in the later sections.

Two-temperature molecular dynamics

While TTM presents a powerful tool in analyzing the e-p coupled thermal transport, its relability heavily depends on the accuracy of the input parameters. In addition, the change of lattice structure such as the defect dynamics cannot be modeled by the Fourier equations alone. To complement its limitations, two-temperature molecular dynamics (TTM-MD) was developed.

MD simulations provides atomic scale calculations that can be used to investigate the properties of materials.^59–62 Using interatomic potentials and the lattice structure as inputs, the simulation provides a classical prediction of the atom behaviors. MD is advantageous in that by simulating the movements of all atoms in the system, it automatically includes many sophisticated physics such as inelastic phonon scattering, etc. Many properties including physical, mechanical and chemical properties can be directly acquired or derived from the simulation results. However, despite the multiple advantages MD also has its limitations. First and foremost is that MD simulation does not include electronic effects, therefore the electrical properties, including electrical thermal conductivity and e-p coupling etc. cannot be modeled. Moreover, a valid prediction of the material properties requires accurate interatomic potentials, but developing potentials, especially for materials with complex lattice structures, is very difficult and the accuracy of the result is still limited by the approximated form of the potential expressions. The first-principles MD can bypass this issue, but comes with the cost of high computation demand.⁶³ Last but not least, the classical simulation with no quantum effects also makes the result deviate from reality, e.g. phonons in thermal equilibrium do not follow Bose-Einstein distribution. Nevertheless, MD still proves to be a useful tool in many cases and has drawn a lot of attention among researchers.

By combining MD with TTM, one is able to include the effect of electrons while keeping MD’s advantages in simulating phonons. Originally developed by Duffy and Rutherford,⁶⁴ the basic idea of TTM-MD is to represent electrons’ scattering with phonons as an “friction force” added to the atoms in the simulation domain.⁶⁵ The system is first discretized into finite volume grids. Electron “gas” is added to each grid with proper electronic properties. The electronic gas has its own thermal diffusion, and atoms will interact with the gas in the same grid. The motion equation for atom i becomes:

where m and υ are the atomic mass and velocity, F_i is the total force exerted on atom i. γ is a friction factor representing the electron-phonon interaction, and it is related to G_ep as follows:

where n is the atom density. is a random force term commonly seen in Langevin dynamics and has the following expression:

where Δt is the time step in the MD simulation and R is a random vector [R₁, R₂, R₃] with each element ranging from -0.5 to 0.5. Then the coupling term in Eq. ([eq:TTM_GE]) can be expressed as:

Then the electronic thermal diffusion and e-p coupling mechanism has been successfully added to MD simulations. TTM-MD has proved to be successful in predicting the thermal transport in metallic systems.^66–69

Phonon-phonon non-equilibrium in single materials

In addition to local thermal non-equilibrium between electrons and phonons, different phonon branches in a single material can also be driven into thermal non-equilibrium. The observation of p-p non-equilibrium dates back to 1980’s when the “hot phonon bottleneck” was observed for photoexcited semiconductors.^70–72Similar with crystalline metals, laser irradiation also causes electron excitation and the subsequent hot electron cooling process in semiconductors. However, due to the fact that semiconductors have optical phonons, while most metals only have acoustic phonons, the cooling processes are significantly different. Experiments have shown that the hot electron cooling rates in photoexcited GaAs are significantly slower than the TTM prediction of Eq. ([eq:P_TTM_GE]), which has been ascribed to the presence of non-equilibrium phonons generated during the thermalization process.^73–77Recent works on solar cell materials, especailly perovskites, have observed similar results.^78,79 Due to selective e-p coupling between electrons and different phonon branches [80], acoustic and optical phonons are heated up by electrons at different rates. Optical phonons, which have stronger coupling to electrons, have a faster and larger temperature rise than acoustic phonons. Their temperature rapidly approaches that of electrons which prevents electrons from cooling, forming the hot phonon bottleneck.

Besides semiconductors, p-p non-equilibrium has also been observed in 2D materials and metals. For example, single-layer graphene (SLG) has exceptional electrical and thermal properties, and is often treated as a metallic material. The reported values of the thermal conductivity of single-layer graphene (SLG) ranges from 600 W/mK to as high as 5800 W/mK.^81–85 While no consensus has been reached on the explanation of this issue,^86–89 non-equilibrium thermal transport is supposed to be one of the reasons due to its existance in Raman spectroscopy measurements which involves laser heating. Experiments on short-pulse laser excitation of noble metals such as Al and Au has also been performed, along with proposed theoretical models to describe the p-p non-equilibrium.^11,90

Similar with the e-p non-equilibrium case, the direct acquisition of the temperatures of electrons and different phonon branches are extremely difficult, therefore researchers turn to theoretical modeling approaches to predict the detailed process. As is mentioned, the generalization of TTM failed when experimental results showed that the actual cooling rate can be more than 4 times slower than the TTM prediction,⁷⁹ as is shown in Fig. 2. Therefore a more sophisticated method is needed to describe the p-p non-equilibrium. Refined models based on Eq. ([eq:P_TTM_GE]) has been proposed to take into account the p-p non-equilibrium, where acoustic and optical phonons are considered as separated particles and together with electrons they form a “three temperature model”.⁷⁵ The refined model can capture the local thermal non-equilibrium among phonons and give results that better match the experiments.

$Illustration of hot electron cooling process in lead-iodide perovskites under different incident photon intensity measured in experiment compared with TTM predictions in Ref. . Symbols represent the experimental electron temperature cooling curve while the dash line represents the TTM prediction. n_{0} denotes the incoming photon density determined from laser power. It is clearly seen that the actual cooling rate is much slower than prediction at high power laser input.$

Fig. 2 Illustration of hot electron cooling process in lead-iodide perovskites under different incident photon intensity measured in experiment compared with TTM predictions in Ref. 79. Symbols represent the experimental electron temperature cooling curve while the dash line represents the TTM prediction. n₀ denotes the incoming photon density determined from laser power. It is clearly seen that the actual cooling rate is much slower than prediction at high power laser input.

Boltzmann transport equations

Spectral BTE can provide rigorous solutions as a semi-classical approach to describe the statistical behavior of thermodynamic systems.^25,26 A general form of the governing equation describes how the energy carrier’s population changes with the scattering, its momentum and external forces. For example the governing equation for phonons is:

where n denotes the Bose-Einstein distribution function and is the index for modes, F is external force, and are the trajectory and momentum vector respectively. The rigorous solution of Eq. ([eq:bte]) of all phonons is very cumbersome, and generally Monte Carlo (MC) simulation is used.^91–94 Spectral phonon transport based on MC calculation has been presented,⁹⁵ still due to the high computation demand it is still limited to small systems with low dimensions, which hinders its application. Therefore for more practical applications several simplification are often made and corresponding models are developed. One of the practical simplifications is the application of relaxation time approximation (RTA), which assumes the scattering between particles does not affect the their distributions. It is commonly used to simplify the calculation of the scattering term on the right side^96,97

where n₀ is the equilibrium distribution, and τ is the relaxation time of the particle. Then Eq. ([eq:bte]) can be integrated over the frequency range to get its energy density form:

Eq. ([eq:bte_e]) is the more commonly-used form of BTE, which focuses on a physical quantity that is more straightforward to acquire the heat fluxes and temperatures in the calculation domain. Based on RTA and an analogy between phonon and photon radiative transport, EPRT has been developed which expressed the phonon energy flow in terms of radiation intensity.^25,27 The temperature profile calculated from EPRT is usually non-linear, which deviates from Fourier’s Law prediction.²⁹ Another approach is the lattice Boltzmann method (LBM) which is based on a discretized form of Eq. ([eq:bte_e]).^98–100 LBM makes it more convenient for BTE to be applied in numerical simulations, but the accuracy also becomes highly dependent on how the system is meshed, especially in high-dimension cases. These simplified BTE approaches have since been widely applied in investigating the phononic heat transport. It has proved to be successful in predicting the size effect of thermal conductivity in thin films which the classical methods such as Fourier’s Law fails to capture.

When applied to non-equilibrium thermal transport, Eq. ([eq:bte]), or Eq. ([eq:bte_e]) when RTA is still assumed, is applied to electrons and each phonon polarization. However, caution should be taken as RTA requires that the scatterings have negligible effect on the carrier distributions, which will be weakened in materials where the p-p scattering is strong.¹⁰¹

Based on RTA, Chen et. al. developed the spectral BTE in dielectric materials which greatly simplifies the phonon scattering calculation and can resolve the steady-state temperature profiles of different phonon modes in thin films with fixed boundary temperatures.^96,102 A more recent work by Sadasivam et. al. utilized first-principles BTE to model the hot electron cooling process in several semiconductors,¹⁰³ as is shown in Fig. 3. Their results showed that the hot electron cooling time is 1-15 ps longer than the simple TTM prediction, which is consistent with experimental observations. A generalized TTM was proposed to take into account the phonon thermalization.

Comparison of calculated electron cooling and lattice heating curves in BN and BAs after electronic excitation acquired from TTM and first-principles spectral BTE from Ref. . The hot electron relaxation process predicted by BTE is more than 10 ps slower than TTM prediction in both materials.

Fig. 3 Comparison of calculated electron cooling and lattice heating curves in BN and BAs after electronic excitation acquired from TTM and first-principles spectral BTE from Ref. 103. The hot electron relaxation process predicted by BTE is more than 10 ps slower than TTM prediction in both materials.

Multi-temperature model

Besides BTE, models beyond TTM have also been developed. By extending TTM to also include p-p coupling, Waldecker et. al. proposed a model similar with RTA BTE but has the same calculation simplicity as TTM.¹¹ Symmetric p-p coupling strength is assumed for all phonons and the parameter is acquired by fitting the temperature cooling curve to their experimental data. Though based on hypotheses, their model is successful in manifesting the local non-equilibrium among acoustic phonons in metal, which matches their corresponding pulse laser heating experiment. A contemporary work by Vallabhaneni et. al. implemented RTA in a similar model, and extended TTM to a general multi-temperature model (MTM), which can model non-equilibrium thermal transport with thermal properties acquired from first-principles calculations.¹² The governing equations of MTM are:

Compared with Eq. ([eq:TTM_GE]), the e-p coupling term is modified to a branch-resolved form and a phonon-phonon coupling term is added for each phonon branch. Here i is the index for phonon branches, and G_ep,i is the coupling factor between electrons and phonon branch i. The summation of G_ep,i over all the phonon branches will lead to the e-p coupling factor G_ep in TTM:

The phonon-phonon scattering is represented by the coupling between each phonon branch and the “scattering lattice reservoir" which is represented by an averaged T_lat. It is an analogy of the e-p coupling using RTA based on the assumption that the phonon-phonon scattering has negligible effect on phonon distribution and phase space. G_pp,i is the phonon-phonon coupling factor between phonon branch i and the scattering lattice reservoir, which is calculated using RTA from the following equation:^104–107

where τ_i is the relaxation time of phonon branch i, and the scattering lattice reservoir’s temperature T_lat is defined to ensure the energy transfer among phonon branches is conserved:

Then by solving Eq. ([eq:MTM_GE]) numerically, we can obtain the transient and steady state temperature profile of the system with resolved electron and phonon branch temperatures. It can be seen that our MTM in Eq. ([eq:MTM_GE]) is essentially a spectral treatment of e-p coupling rather than the gray treatment in TTM in Eq. ([eq:TTM_GE]). A constant laser heating process of SLG which represents the Raman Spectroscopy measurement was simulated, as is shown in Fig. 4. Results showed great p-p non-equilibrium due to selective e-p coupling. The predicted k vary as one chooses different modal temperatures to perform the calculation, with the average value almost 2 times as large as that of a simple TTM.^12,13 A later work by Lu et. al. extended Vallabhaneni’s work by including the investigation into transient laser heating in SLG, and it was shown that due to the hot phonon bottleneck induced by phonon non-equilibrium, MTM predicted hot electron cooling rates 60 times slower than the original TTM.¹⁴

Fig. 4 a) Illustration of SLG sample irradiated by laser and b) corresponding MTM prediction of the resolved temperature profiles showing e-p and p-p non-equilibrium induced by selective e-p coupling during a Raman spectroscopy experiment Ref. 12. c) Experimental measurement of the temperatures of different phonon modes in Ref. 13.

Electron-phonon non-equilibrium at metal-dielectric interfaces

Thermal transport across interfaces between different materials is a core topic that concerns the research community nowadays, especially when the system dimension is very small such as the nano-size electronic devices and multi-layer structures. At such scale, interfacial thermal transport dominates over the conduction inside each material, and the prediction of interfacial thermal conductance is crucial for the design and evaluation of these nano-structures.

When heat goes across the interface between two different materials, it will experience temperature discontinuity due to the lattice mismatch between the two materials as well as surface conditions, etc. The interfacial thermal conductance (ITC) h_Bd is defined as the ratio of the heat flux over the temperature jump ΔT, and the thermal boundary resistance (TBR) R_Bd is defined as the inverse of h_Bd:

R_Bd is determined by many factors as mentioned above. If we consider an ideal interface while isolating the surface conditions and external forces, then R_Bd is only determined by the lattice and electronic properties of the two materials.

Traditional AMM and DMM approaches have been widely used due to their simplicity, and have proved to be successful in the prediction of ITC at low temperatures (<100 K).¹⁶However, the measured results often do not agree with theoretical predictions, and the discrepancy is not well understood.^55,108–113 Thermal non-equilibrium may play an important role, so in this section, we will review the works done on this topic.

Metal-dielectric interfaces appear widely in modern devices. While the thermal transport in dielectric is mostly through lattice vibration, or the phonons, the thermal transport in metals is complicated by the coupled transport between electrons and phonons. At these interfaces, usually electrons are considered as not transmittable to the dielectric while phonons can transfer energy between the two materials. This will cause e-p thermal non-equilibrium near the interface even if the system is in steady state. Experiments have shown that the thermal conduction in metal-semiconductor multilayers with period below 10 nm has been found to be dominated by phonons regardless of the presence of electrons.^114–116 Combined with theoretical calculations, it is demonstrated that e-p interactions have a strong impact at the interfacial thermal transports.

At a 1D metal-dielectric interface, usually electrons are considered as insulated while phonons can transfer energy between the two materials. The different boundary conditions are as follows:

x = 0 denotes the position of the interface. The different boundary conditions will cause the two temperature profiles to deviate from each other near the interface, resulting in the so-called e-p non-equilibrium, as is shown in Fig. 5. This presents an extra temperature jump in addition to the one introduced by the p-p coupling R_Bd , and can be conceived as an additional TBR R_ep in series with the phonon coupling resistance. Wang et al has derived an expression for R_ep by solving Eq. ([eq:TTM_GE]) while taking the extreme limit when the film is infinite large so that electrons and phonons are in thermal equilibrium away from the interface:¹⁹

Illustration of e-p non-equilibrium near a metal-dielectric interface derived from TTM in Ref. . Electrons are insulated from entering the non-metal therefore have a flat curve near the interface. Their energy is dumped to phonons and transferred to the other side. The e-p non-equilibrium results in an additional TBR in series with the p-p TBR.

Fig. 5 Illustration of e-p non-equilibrium near a metal-dielectric interface derived from TTM in Ref. 58. Electrons are insulated from entering the non-metal therefore have a flat curve near the interface. Their energy is dumped to phonons and transferred to the other side. The e-p non-equilibrium results in an additional TBR in series with the p-p TBR.

TTM has since been widely applied in both experimental and modeling studies of metal-dielectric interfaces.^117,118 Generally the studies have consistent results with the theoretical TTM prediction, that the e-p non-equilibrium introduces thermal resistance at metal-dielectric interfaces. This decreases the measured interfacial conductance as well as the effective thermal conductivity, as is shown in Fig. 6. Li et. al. applied TTM to interpret their experimental results of thermal conduction in Mo-Si multilayers with period below 10 nm.¹¹⁴ The e-p non-equilibrium is reported to be responsible for the low thermal conductivities k measured, which is almost identical to the materials’ lattice portion k_p. A later study of Lu et. al. applied TTM in metal nanoparticle-polymer matrix, and derived a more general form of R_ep.¹¹⁹ Their analysis shows that due to e-p non-equilibrium the thermal conduction in metal thin films sandwiched between dielectrics is dominated by phonons, which is consistent with Li’s work. TTM has also been combined with effective medium theory (EMA) in modeling random size-metal nanoparticle-composite materials in higher dimensions.^120–123 Consistently the addition of e-p non-equilibrium to the original p-p TBR will further reduce the prediction of the effective k of the composite material. Simplified TTM-BTE in 1D has also been solved practically using the equation of phonon radiative transport (EPRT) and lattice Boltzmann methods.^124,125 In addition, TTM is also extended to “three temperature model" with fitted parameters from experiments to describe the cross-interface e-p coupling at Au-Si and Au-glass interfaces in the studies of Hopkins et. al.^52,126. While the classical theoretical basis of TTM determines that it cannot include anharmonicity lattice effects which has been demonstrated to be significant at interfacial transports,^127–131 impletation of TTM in MD enables the investigation of electronic properties in metals while still including this complex physics in a practical way.^66,132,133 In a series of their works, Wang and Lu implemented TTM in molecular dynamics, adding the electronic participation in the classical atomic simulation tool while keeping its advantage in modeling complicated lattice vibrations.^19,134 By modifying the model of Duffy and Rutherford,⁶⁴ they enabled TTM-MD simulations for thermal interfacial transport for the first time. The simulation results of ITC at Cu-Si and Cu-graphene interfaces show that R_ep usually have comparable value with the p-p TBR, accounting for 30% to 70% of the total TBR.

$Results of interfacial thermal resistance R_{I} from classical MD and TTM-MD simulations with different e-p coupling strength from Ref. . The trends of R_{ep}, which is the difference between the R_{I}’s from MD and TTM-MD respectively, are consistent with Eq. ([eq:Rep_Wang]).$

Fig. 6 Results of interfacial thermal resistance R₁ from classical MD and TTM-MD simulations with different e-p coupling strength from Ref. 19. The trends of R_ep, which is the difference between the R₁’s from MD and TTM-MD respectively, are consistent with Eq. ([eq:Rep_Wang]).

Phonon-phonon non-equilibrium at interfaces

Besides e-p non-equilibrium at metal-dielectric interfaces, p-p non-equilibrium also exists. When heat transfers across the interface between two different materials, different phonon modes have different transmission coefficients, which results in a thermal non-equilibrium similar with the e-p non-equilibrium at metal-dielectric interfaces. This phenomenon is especially significant in several 2D materials due to the large difference in properties between the in-plane (IP) and out-of-plane (OP) phonons, as has been observed in Ref. 135. As is shown in Fig. 8, the temperatures of different phonon branches deviate significantly from each other at an SLG-BN interface.

Temperature profile of acoustic and optical phonons at a solid-solid interface acquired from MD simulations from Ref. . The local non-equilibrium between the two branches in the solid on the right is clearly shown. The solid on the left only has acoustic phonons.

Fig. 7 Temperature profile of acoustic and optical phonons at a solid-solid interface acquired from MD simulations from Ref. 127. The local non-equilibrium between the two branches in the solid on the right is clearly shown. The solid on the left only has acoustic phonons.

MD has been applied to investigate p-p non-equilibrium at dielectric interfaces. However, despite MD’s advantage in practically modeling complicated lattice vibration physics, acquiring the resolved temperature profiles needs additional post-processing of simulation data. In the study of Wu et. al., non-equilibrium MD simulation was performed on an imaginary solid-solid interface.¹²⁷ By calculating the frequency-resolved energy of the phonons, they separated the temperatures of acoustic and optical phonons. As is shown in Fig. 7, a local thermal non-equilibrium is clearly observed in the material on the right, which has a diatomic strucuture and 6 phonon branches while the material on the left only has acoustic phonons. The more recent study of Feng et. al. utilized MD simulation on several systems made of different materials including an SLG-BN interface.¹³⁵ Spectral analysis was then performed to calculate the temperature of each phonon branch,^136,137 as is shown in Fig. 8. The optical phonon modes have significantly larger interfacial temperature jump than the acoustic phonons, especially the ZO phonons whose band overlap in the phonon dispersions of SLG and BN is negligible. The results are consistent with the expectation that p-p non-equilibrium can be induced by their different transmission at the interface. Their later following works further presented the branch-wise condutance of different phonons, along with a modified Landauer approach which includes inelastic phonon scattering.¹³¹ In a later MD study conducted by An et. al., simulation is performed across an SLG-SLG interface.¹⁰¹ Phonons are divided into two groups: IP phonons (LA, TA, LO and TO phonons) and OP phonons (ZA and ZO phonons), and in one SLG only in-plane vibration is allowed thus removing ZA and ZO phonons in that side. A spectral analysis revealed similar IP-OP non-equilibrium with the e-p non-equilibrium at metal-dielectric interfaces, which is as expected since an analogy of IP-OP coupling to e-p coupling can be made with similar boundary conditions. Their results also showed that IP and OP phonons have weak coupling, so that their distributions will not be significantly affected by the scattering and the application of RTA is validated in SLG.

Illustration of phonon thermal non-equilibrium induced by different transmission at an SLG-BN interface from Ref. .

Fig. 8 Illustration of phonon thermal non-equilibrium induced by different transmission at an SLG-BN interface from Ref. 135.

Spectral BTE based on RTA, which greatly reduces the calculation needed, has been applied in 1D Si-Ge superlattice.¹⁰⁵ The p-p transmission at the interfaces are calculated by phonon branch-resolved DMM. Their results clearly demonstrate the p-p non-equilibrium at Si-Ge interfaces induced by different phonon transmission, which can be seen in Fig. 9. Similar with e-p non-equilibrium, the p-p non-equilibrium also introduces thermal resistance. Incorporating p-p non-equilibrium will predict a lower ITC than a simple DMM or AMM result based on local equilibrium assumption. In analogy to TTM which can be applied at interfaces when combined with boundary conditions, MTM can also be extended. A recent work by Lu et. al. combined AMM with MTM and simulated the thermal transport across Si-Ge interfaces, as is shown in Fig. 10. Their results of the reduced interfacial conductance as a consequence of the phonon local non-equilibrium in the bulk materials near the interface, showed that non-equilibrium thermal transports in the bulk and across the interface are not decoupled.¹³⁸

Fig. 9 Temperature profile predicted by spectral BTE based on RTA across a Si-Ge interface from Ref. 105. Different line styles represent different system dimensions, while different colors represent temperature of the lattice, acoustic phonons or optical phonons. Due to different transmission coefficients at the interface, the temperatures of the two phonon groups deviate from each other, resulting in local thermal non-equilibrium.

$Temperature profile predicted by MTM across a Si-Ge interface from Ref. . Temperatures of different phonon branches deviate from each other near the interface due to different interfacial transmission coefficients. T_{fit} refers to the fitted temperature acquired by extrapolating the linear (thermal equilibrium) part of the temperature profile to the interface. As a result the interfacial conductance is decreased due to the thermal non-equilibrium.$

Fig. 10 Temperature profile predicted by MTM across a Si-Ge interface from Ref. 138. Temperatures of different phonon branches deviate from each other near the interface due to different interfacial transmission coefficients. T_fit refers to the fitted temperature acquired by extrapolating the linear (thermal equilibrium) part of the temperature profile to the interface. As a result the interfacial conductance is decreased due to the thermal non-equilibrium.

Summary and outlook

In this review, we have discussed the physical significance of non-equilibrium thermal transport and recent numerical simulation studies to investigate this feature. Non-equilibrium thermal transport consists of e-p non-equilibrium and p-p non-equilibrium. They have significant effect and are manifested during thermal transport across the interface between two different materials, as well as in materials under irradiation. Generally thermal non-equilibrium introduces resistance to the system, e.g. adding extra TBR at an interface, or slows down the hot electron cooling rate in pulse laser-irradiated metals. If treated improperly it can derive misleading results and predictions, such as the under-estimation of thermal conductivity in Raman spectroscopy measurements. While experimental observations are limited to detecting their impacts rather than directly measuring the resolved deviated temperatures, multiple simulation approaches compensate this by predicting the non-equilibrium thermal transport process under different conditions. TTM, BTE, MD and MTM are most widely used and have been successful in predicting the local thermal non-equilibrium, bridging the discrepancy between theory and measurements caused by misleading interpretations from simple methods based on local equilibrium assumption.

Numerical simulations provide excellent complement to the study of non-equilibrium thermal transport, as well as extending beyond the limitations of current experiment techniques. Their predictions have helped people better interpret the measurements as well as provided guidance in the design of systems and devices, which are expected to be more efficient and effective in the future. Not just limited to non-equilibrium thermal transport, but also in other fields concerning thermal transport and other physics. Still, the current limitations in simulation approaches should not be overlooked. As is mentioned, MD simulations are limited by their incapability of including quantum effects and its reliability strongly depends on the accuracy of the interatomic potentials deployed. Though spectral BTE can always provide rigorous solutions, the high computation cost of a full-scattering model has rendered it too complicated. Approaches such as the MC methods have been developed, but the calculation is still far from practical. TTM and MTM are both based on diffusive thermal transport equations, and may lose ballistic transport physics at nano-scales. Recently first-principles calculations have received increasing popularity, but similar to spectral BTE the high computation cost is a major issue. Unresolved questions also remain in the area of non-equilibrium thermal transport, such as how the apparent temperature is related to the multiple resolved temperatures of different energy carriers, and how significant the p-p non-equilibrium is at dielectric interfaces as mentioned above. There are also possible topics that emerge from the e-p non-equilibrium near the interface, such as the cross-interface electron-phonon coupling.^{126,134,139–142} These limitations, open questions and potential research aspirations certainly deserve investigating, and we expect more rigorous and efficient modeling approaches in the future.

References

D. G. Cahill, P. V. Braun, G. Chen, D. R. Clarke, S. Fan, K. E. Goodson, P. Keblinski, W. P. King, G. D. Mahan, A. Majumdar, H. J. Maris, S. R. Phillpot, E. Pop and L. Shi, Appl. Phys. Rev. 2014, 1, 011305.