Spectral analysis of non-equilibrium molecular dynamics: spectral phonon temperature and phonon local non-equilibrium in thin films and across interfaces

TL;DR

This paper introduces a molecular dynamics-based spectral phonon temperature method to observe local thermal nonequilibrium in nanomaterials, revealing mode-specific temperature differences and interfacial thermal resistance mechanisms.

Contribution

It develops a novel spectral phonon temperature calculation method within molecular dynamics, enabling direct observation of phonon mode temperatures and nonequilibrium in nanostructures.

Findings

Spectral phonon temperature reveals local thermal nonequilibrium in thin films and interfaces.

Interfacial thermal resistance is influenced by phonon mode coupling and mixing.

Mode-specific phonon interactions affect thermal transport across interfaces.

Abstract

Although extensive experimental and theoretical works have been conducted to understand the ballistic and diffusive phonon transport in nanomaterials recently, direct observation of temperature and thermal nonequilibrium of different phonon modes has not been realized. Herein, we have developed a method within the framework of molecular dynamics to calculate the temperatures of phonon in both real and phase spaces. Taking silicon thin film and graphene as examples, we directly obtained the spectral phonon temperature (SPT) and observed the local thermal nonequilibrium between the ballistic and diffusive phonons. Such nonequilibrium also generally exists across interfaces and is surprisingly large, and it provides an additional thermal interfacial resistance mechanism. Our SPT results directly show that the vertical thermal transport across the dimensionally mismatched graphene/substrate interface is through the coupling between flexural acoustic phonons of graphene and the longitudinal phonons in the substrate with mode conversion. In the dimensionally matched interfaces, e.g. graphene/graphene junction and graphene/boron nitride planar interfaces, strong coupling occurs between the acoustic phonon modes on both sides, and the coupling decreases with interfacial mixing. The SPT method together with the spectral heat flux can eliminate the size effect of the thermal conductivity prediction induced from ballistic transport. Our work shows that in thin films and across interfaces, phonons are in local thermal nonequilibrium.

Figures

Equations

$$n_{\lambda }=\frac{k_B{T}_{\lambda }}{\hslash {\omega }_{\lambda }}.$$

$$E_{\lambda }=n_{\lambda }ϵ=\frac{k_B{T}_{\lambda }}{\hslash {\omega }_{\lambda }}\hslash {\omega }_{\lambda }=k_B{T}_{\lambda }.$$

$$⟨E_{K,\lambda }⟩=\frac{1}{2}{k}_B{T}_{\lambda }.$$

$$E_{K,\lambda }=\frac{1}{2}{\dot{Q}}_{\lambda }^{\ast }{\dot{Q}}_{\lambda },$$

$${\dot{Q}}_{\lambda }(t)=\frac{1}{\sqrt{N_c}}\sum _{l,b}^{N_c,n}\sqrt{m_b}\exp (-i\mathbf{k}\cdot {\mathbf{r}}_{l,b}){\mathbf{e}}_{b,\lambda }^{\ast }\cdot {\dot{\mathbf{u}}}_{l,b;t}.$$

$$T_{\lambda }=⟨{\dot{Q}}_{\lambda }^{\ast }(t){\dot{Q}}_{\lambda }(t)⟩/k_B,$$

$$V(r_{ij})=4ϵ\left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6\right]$$

$${\bar{T}}_{\nu }=\frac{1}{N_{\mathbf{k}}}\sum _{\mathbf{k}}^{N_{\mathbf{k}}}{T}_{\mathbf{k},\nu }=\frac{\int T_{\mathbf{k},\nu }d\mathbf{k}}{\int d\mathbf{k}},$$

$${\bar{T}}_{\lambda }=\frac{1}{N_{\lambda }}\sum _{\lambda }^{N_{\lambda }}{T}_{\lambda }=\frac{1}{3n}\sum _{\nu }^{3n}{\bar{T}}_{\nu },$$

$${\kappa }_{\lambda }=\frac{q_{\lambda }}{\nabla T_{\lambda }},$$

$$q_{\lambda }=\sum _{l,b}^{N_c,n}⟨\frac{1}{\sqrt{N_c{m}_b}}\left[E_{l,b}(t)-{\mathbf{S}}_{l,b}(t)\right]{\mathbf{e}}_{b,\lambda }\exp (i\mathbf{k}\cdot {\mathbf{r}}_{l,b}){\dot{Q}}_{\lambda }⟩.$$

$$\kappa =\frac{q_{\mathrm{MD}}}{\nabla T_{\mathrm{MD}}}=\frac{{\sum }_{\lambda }{q}_{\lambda }}{\nabla T_{\mathrm{MD}}},$$

$$\kappa =\sum _{\lambda }\frac{q_{\lambda }}{\nabla T_{\lambda }},$$

$$E_{K,atom}=3N{k}_{B}T/2.$$

$$E_{K,phonon}=\frac{1}{2}\sum _{\lambda }^{3N}{n}_{\lambda }\hslash {\omega }_{\lambda }=\frac{1}{2}\sum _{\lambda }^{3N}{E}_{\lambda }=\frac{1}{2}\sum _{\lambda }^{3N}{K}_B{T}_{\lambda }.$$

$${\bar{T}}_{\lambda }=\frac{1}{3N}\sum _{\lambda }^{3N}{T}_{\lambda }=\frac{1}{3n}\sum _{\nu }^{3n}{\bar{T}}_{\nu },$$

$${\bar{T}}_{\nu }=\frac{1}{N_c}\sum _{\mathbf{k}}^{N_c}{T}_{\mathbf{k},\nu }=\frac{\int T_{\mathbf{k},\nu }d\mathbf{k}}{\int d\mathbf{k}}.$$

$${\bar{T}}_{\nu }=\frac{\int k^2{T}_{\mathbf{k},\nu }dk}{\int k^{2}dk}\approx \frac{{\sum }_m{k}_m^2{T}_{\mathbf{k},\nu }}{k_m^2},3D$$

$${\bar{T}}_{\nu }=\frac{\int k{T}_{\mathbf{k},\nu }dk}{\int kdk}\approx \frac{{\sum }_m{k}_m{T}_{\mathbf{k},\nu }}{k_m}.2D$$

$$E_{\lambda }(\omega )={\left|\mathcal{F}[{\dot{Q}}_{\lambda }(t)]\right|}^2=\frac{C_{\lambda }}{(\omega -{\omega }_{\lambda }^A{)}^2+({\tau }_{\lambda }^{-1}{)}^2/4}.$$

$$\Phi (\mathbf{k},\omega )=\sum _{\nu }^{3n}{E}_{\lambda }(\omega )$$

$$=\frac{1}{4\pi t_0}\sum _{\alpha }^3\sum _b^n\frac{m_b}{N_c}{\left|\sum _l^{N_c}{\int }_0^{t_0}{\dot{u}}_{\alpha }^{l,b}(t)\exp (i\mathbf{k}\cdot {\mathbf{r}}_0^l-i\omega t)dt\right|}^2.$$

Full text

I Introduction

The transport of phonons drives the heat dissipation and shielding, thermoelectric energy conversion, energy storage and saving, etc. With the advances of nanotechnology, nowadays devices evolve towards smaller size which can be even smaller than phonon mean free path (MFP). In this case, phonons become ballistic and can travel without scattering or energy dissipation. Therefore, recently extensive experimental and theoretical work were focused on the study of ballistic and diffusive phonon transport in nanomaterials such as silicon Minnich (2012); Regner et al. (2013); Hu et al. (2015); Zeng et al. (2015); Cuffe et al. (2015), holey silicon Lee et al. (2015), silicon nanomesh Zen et al. (2014); Maldovan (2015), superlattices Ravichandran et al. (2013); Maldovan (2015); Luckyanova et al. (2012), SiGe alloy nanowire Hsiao et al. (2013), etc. By examining the size dependent thermal conductivity, researchers have found a considerable portion of ballistic phonons in the materials of several hundred nanometers to micrometer. Even though the transport properties have been extensively investigated Regner et al. (2013); Hu et al. (2015); Zeng et al. (2015); Lee et al. (2015); Zen et al. (2014); Ravichandran et al. (2013); Maldovan (2015); Luckyanova et al. (2012); Hsiao et al. (2013); Feng and Ruan (2014), the study particularly focused on an indirect demonstration, leaving the information of phonon temperature unknown and a direct observation of the local nonequilibrium among phonons unrealized. On the other hand, the phonon transport across interface is crucially important in real devices. Even though extensive works have been conducted Minnich (2015); Tian et al. (2012); Dechaumphai et al. (2014); Seol et al. (2010); Chen et al. (2014); Li and Yang (2012); Vermeersch et al. (2014); Siemens et al. (2010); Zhou et al. (2016), the understanding of spectral phonon interfacial conduction are still based on simple models and approximations. Thus, a direct and accurate understanding of the interfacial phonon transport is urgently needed. Recently Dunn et al. Dunn et al. (2016) have studied the non-equilibrium of lattice vibrations in nonequilibrium molecular dynamics (NEMD). Their results are on the frequency level while not on the phonon modal level. Another recent study Rurali et al. (2016) has concluded that interfaces are at local thermal equilibrium, but the spectral nature of phonons was not considered. Here, we report a spectral phonon temperature (SPT) method that directly calculate the spatial temperatures of all the phonon modes in real systems within the framework of NEMD, and we have directly observed the local temperature nonequilibrium among phonons in nanomaterials and across interfaces. It is demonstrated as an effective way to probe the spectral phonon thermal transport mechanisms across interface.

II Spectral Phonon temperature formalism

In MD, the phonon population is described by the Boltzmann distribution

[TABLE]

The total energy of the phonon mode $\lambda$ is the per phonon energy $\epsilon=\hbar\omega_{\lambda}$ multiplied by its population

[TABLE]

Here $\lambda$ is short for $(\mathbf{k},\nu)$ with $\mathbf{k}$ and $\nu$ representing the phonon wave vector and dispersion branch, respectively. Based on the energy equipartition theorem, the time averaged kinetic energy $\langle E_{K,\lambda}\rangle$ and potential energy $\langle E_{V,\lambda}\rangle$ are both half of the total energy, e.g.,

[TABLE]

Based on the lattice dynamics Dove (1993), the kinetic energy of the mode $\lambda$ is

[TABLE]

where $\dot{Q}_{\lambda}(t)$ is the time derivative of normal mode amplitude, which is given by the Fourier transform of atomic displacement in the real space

[TABLE]

$l$ and $b$ label the indexes of the primitive cells and basis atoms with the total numbers represented by $N_{c}$ and $n$, respectively. $m$, $\mathbf{r}$, and $\dot{\mathbf{u}}$ are the mass, equilibrium position, and velocity vector, respectively. $\mathbf{e}_{b,\lambda}^{*}$ is the complex conjugate of the eigenvector component at the basis $b$ for the mode $\lambda$. By comparing Eq.(3) and (4), we can get the temperature of the phonon mode $\lambda$

[TABLE]

where $\langle\rangle$ denotes the time average. To eliminate the fluctuation MD, Eq.(6) needs to be averaged over an enough long time.

Note that concept of phonon normal mode is still valid under non-equilibrium Zhou and Hu (2015). In Eq.(2), the temperature of a phonon mode is defined as a convenient representation of the carrier energy density, as commonly done in literature for both experimental and theoretical studies, such as Ref.Mann et al. (2006); Sullivan et al. ; Maassen and Lundstrom (2016); Vallabhaneni et al. (2016a). Another theoretical approach, spectral Boltzmann transport equation, also uses the concept of non-equilibrium temperatures of carriers such as electrons and phonons Ni et al. (2012). Note that the temperature in MD is defined by the atomic kinetic energy ($E_{K}=k_{B}T/2$) rather than the total energy $E_{K}+E_{V}=k_{B}T$, although theoretically $\langle E_{K}\rangle=\langle E_{V}\rangle$.

III Simulation Setups

Before applied to NEMD, our SPT method was firstly validated in equilibrium MD, in which all the phonons are in equilibrium and thus should have the same temperature. This is verified by our results in both 3D system Si and 2D system graphene (see the Appendix A).

In the following, we use our SPT method based on NEMD to study the ballistic and diffusive nature of the phonon transport in nanomaterials, e.g, Si thin film and graphene, as well as across interfaces, e.g., graphene/Si, graphene/graphene, and graphene/boron nitride. These systems cover the 2D, 3D, 2D-3D, 2D/2D junction and 2D-2D in-plane phonon transport.

The MD simulations were performed using the LAMMPS package Plimpton (1995). The interatomic potentials used for silicon, graphene, and BN are the original Tersoff (1989, 1990), optimized Lindsay and Broido (2010), and the modified Kınacı et al. (2012) Tersoff potentials, respectively. The interactions between graphene and silicon, and between graphene layers are van der Waals (vdW) forces modeled by the Lennard-Jones (LJ) potential

[TABLE]

with the parameters $\epsilon_{C-Si}=8.909$ meV, $\sigma_{C-Si}=3.629$ ÅRappé et al. (1992), $\epsilon_{C-C}=4.6$ meV and $\sigma_{C-C}=3.276$ ÅNicklow et al. (1972); Lindsay et al. (2011). The cutoffs for C-Si and C-C are $5\sigma$ and $3\sigma$, respectively. Here $r_{ij}$ is the distance between the atoms $i$ and $j$. The parameters can best fit the $c$-axis phonon dispersion and layer separation of 0.335 nm for graphite. The lattice constants at room temperature are obtained by relaxing the atomic structures in MD. The time step for simulations is set as 0.5 fs which is short enough to resolve all the phonon modes. For silicon, the time step can be set as 1 fs as well, and the results do not change. The two ends of the systems are fixed during the simulations as shown in Fig. 1 (a)-(e), except for the graphene/Si system in which only the bottom boundary of silicon is fixed as shown in Fig. 1 (d). We first relax the geometries under constant pressure and temperature (NPT) for 10 ns, and then change the ensemble to constant energy and temperature (NVE) except for the fix boundaries and heat reservoirs. After that, we apply the temperature differences on the two reservoirs and stabilize the heat current as well as temperature gradient for 10 ns (20 million time steps). Finally the simulations are run for 20 ns to extract the atomic velocities every 10 time steps. Totally, the atomic velocities of 2 million time steps are stored and used to calculate the phonon temperature. These setups are found to be able to give stable results. The entire simulation domain is divided into many cells, and a temperature is calculated for each cell using atomic velocities in it. One such cell is the sample volume drawn in Fig. 1, which indicates that the atoms in that volume are used to calculate the phonon temperatures at that position. The sample volume contains $8\times 8\times 8$ conventional cells with 4096 atoms for silicon, and $20\times 50$ with 2000 atoms for graphene and BN. To study a subtler spatial distribution of phonon temperature, a smaller sample volume can be used. Although some phonon modes are delocalized, we can still calculate their energy density inside the sample volume by transforming the real periodic vibration of the atoms to phase space. For the ballistic phonons, we indeed see that their temperature is almost a constant throughout the system, which indicates that they are indeed delocalized phonons. While for modes with short mean fee path, we can see appreciable temperature gradient. Periodic boundary conditions are applied to the lateral directions to model infinite large dimensions. The temperatures applied on the two reservoirs are 350 K and 250 K for Fig. 1 (a) Si; 310 K and 290 K for Fig. 1 (b) SLG; 350 K and 250 K for Fig. 1 (c) graphene layers; 500 K and 300 K for Fig. 1 (a) graphene on substrate; and 325 K and 275 K for Fig. 1 (e) C/BN. Different simulation sizes, durations, temperature differences, heat reservoir time constants, and boundary conditions have been tried in our benchmark process, and the conclusions are the same. The main results shown in this paper were obtained using the simulation domain sizes and boundary conditions labeled in Fig. 1. The time constant for the thermal reservoirs are set as 0.5 ps, which is suitable for a stable heat current. The time constant can vary to some extent but not too much. A time constant that is too short (e.g. $<$ 0.01 ps) may lead to unstable temperature profile, while a time constant that is too long (e.g. $>$ 10 ps) is not able to build up a suitable temperature gradient.

IV Spectral phonon temperature in silicon thin films and graphene ribbons

For Si thin film with the thickness $L\approx 152$ nm, the left and right sides are applied with the temperatures $T_{L}=350$ K and $T_{R}=250$ K, respectively, as shown in Fig. 2. Periodic boundary conditions are applied in the $y$ and $z$ directions to mimic an infinitely large cross section. Figure 2 (a) shows the overall MD temperature profile $T_{MD}$, linear in the middle and with two jumps near the thermal reservoirs, which is seen in general NEMD simulations. Then, we use the SPT method to convert $T_{MD}$ to the spectral phonon temperatures as shown in Fig. 2 (b), in which we plot a total of 240 modal temperatures covering 5 representative directions, i.e., [1,0,0], [1,1,0], [1,1,1], [0,1,0] and [0,0,1]. It is seen that most phonons exhibit similar T profile as $T_{MD}$. These phonons are diffusive. Some phonon modes, however, have flatter T profiles as indicated by the light green color, and, they are ballistic phonons. To validate our method and results, we have calculated the branch-averaged temperature $\bar{T}_{\nu}$ and overall-averaged temperature $\bar{T}_{\lambda}$, using

[TABLE]

respectively (See the Appendix B for the derivation). Here $n$ is the number of basis atoms, respectively. We can find that they agree well with $T_{MD}$, as shown in Figs. 2 (c) and (d). To gain more insight, the temperatures of the phonons traveling in the 5 representative directions are separately shown in Figs. 2 (e)-(i). For the phonons traveling along the x direction, nearly half of the modes are (quasi)ballistic. We calculate the MFP of these phonons using a separate independent method – spectral energy density (SED) analysis based on EMD (Appendix C), and find that their MFP are generally longer than the Si thin film thickness 152 nm. The other modes, which are diffusive, are found to have considerably shorter MFP, 60-130 nm, as indicated in Fig. 2 (e). Such an agreement between our NEMD-based SPT results and the EMD-based SED results strengthens the validity of the SPT method. These results agree qualitatively well with recent experimentally reconstructed MFP spectra of silicon Minnich (2012); Cuffe et al. (2015); Hu et al. (2015); Zeng et al. (2015); Regner et al. (2013), while providing more direct observation of ballistic and diffusive phonons whose temperatures have not been measured yet. For the phonons traveling in other directions, their effective phonon MFP is that projected in the x direction $\Lambda_{x}$. Therefore, oblique the traveling directions are, the more diffusive the phonon modes become. For example, in the [1,1,0] direction, $\Lambda_{x}=\Lambda/\sqrt{2}$, and more diffusive modes appear. In [0,1,0] and [0,0,1] directions, $\Lambda_{x}=0$ and all modes are diffusive. Note that in Figs. 2 (e)-(i), a few long-wavelength modes have large uncertainty due to the limitation of MD. To see clearly which phonon modes are ballistic, we designate the modal temperature gradient in the phonon dispersion relation as shown in Figs. 2 (j)-(n). It is seen that the fraction of ballistic modes decreases as phonon traveling direction becomes more oblique. We also tried the Berendsen thermostat and got similar results as shown in Appendix D.

The SPT method also provides a way to calculate the phonon MFP spectra. According to the phonon Boltzmann transport equation (BTE), the phonons’ temperature jumps $\Delta T_{\lambda}$ near the contacts are related to their MFP $\Lambda$ by Jeong et al. (2010); Kaiser et al. (2016) $\frac{T_{L}-T_{R}}{2\Delta T_{\lambda}}=1+L/(\frac{4}{3}\Lambda_{\lambda})$.

We have also obtained the modal temperatures of a 2D system, single-layer graphene, with length of 65 nm. As shown in Fig. 3, along $\Gamma$-M the acoustic phonons are ballistic while most optical phonons are diffusive. As the direction deviates from the heat flow direction, the portion of ballistic phonons reduces. Along $\Gamma$-K’ which is perpendicular to the heat flow, nearly no phonon is ballistic. Since most conclusions are similar to those of Si thin film, we will not reiterate them in the main text (See Appendix E for details).

We note that in both Si thin films and graphene ribbons, the assigned temperature can be well maintained deep inside the reservoir, while cannot be maintained in the region of the reservoir near the channel. As a result, the modal temperature is in equilibrium deep inside the reservoirs while in non-equilibrium near the channel. Dunn et al. Dunn et al. (2016) found phonon non-equilibrium in the reservoirs as they sampled the whole reservoirs to obtain the phonon information. This result is consistent with ours.

V Spectral Phonon temperature and phonon local non-equilibrium across Si/graphene interface

We now investigate the thermal transport across interfaces. Due to the 2D nature, graphene is typically supported on a substrate in real devices. The presence of the substrate can affect the electronic structure, phonon dispersion Zhou et al. (2007); Ong and Pop (2011); Qiu and Ruan (2012), and transport properties Fratini and Guinea (2008); Seol et al. (2010); Chen et al. (2008); Morozov et al. (2008); Qiu and Ruan (2012). For example, the mobility can be reduced by a factor of 10 Fratini and Guinea (2008) and the thermal conductivity reduced by a factor of 5 Seol et al. (2010). The heat dissipation from graphene to substrate is an important thermal management issue in devices Freitag et al. (2009); Yan et al. (2012). Although it is well known that the flexural modes couple most strongly with the substrate Qiu and Ruan (2012); Seol et al. (2010), it has been an open question, however, on which phonon modes in substrate they couple to Sadeghi et al. (2013). We tackle this problem using our spectral phonon temperature method. Figure 4 (a) shows the schematic of the heat conduction from graphene at 500 K to the Si substrate with bottom temperature maintained at 300 K. Periodic boundary conditions are applied to the in-plane directions to eliminate nano-size effect. As a control case, a suspended graphene with T = 500 K is also studied. The wavevector-resolved spectral phonon temperatures in suspended and supported graphene are shown in Fig.  4 (b) panels. Clearly the phonons are in thermal equilibrium in suspended graphene as expected, surprisingly they are in strong nonequilibrium in supported graphene. The branch temperatures of supported graphene and the substrate are shown in Fig.  4 (c) with comparison to $T_{MD}$. In graphene, the out-of-plane acoustic (ZA), transverse acoustic (TA) and longitudinal acoustic (LA) phonons have much lower temperatures than the transverse optical (TO), longitudinal optical (LO) and out-of-plane optical (ZO) phonons. ZA is the lowest among them. While on the surface of the Si substrate, the LA and LO branches have much higher temperature than the other modes. This indicates the most efficient interfacial transport channel is that the acoustic branches especially the ZA phonons in graphene transmit into the LA and LO branches in the substrate. The ZA mode in graphene couples strongly with the substrate because its out-of-plane atomic vibrations directly ”press” the substrate and excite the substrate LA and LO modes, which have the same out-of-plane atomic motion. Phonon mode conversion behavior is surprising and cannot be obtained using other methods such as acoustic or diffuse mismatch methods. Interestingly, the ZO mode, which was believed to couple strongly with the substrate for the in-plane thermal transportQiu and Ruan (2012); Seol et al. (2010), does not contribute much to the vertical thermal transport. This is probably due to that its frequency is higher than the phonon cutoff frequency in silicon, and the elastic transmission process does not allow its transmission. Not only different branches behave differently, but inside a branch the phonon with different wavevectors are in quite large nonequilibrium, as shown in Fig. 4 (b). The lower-energy acoustic phonons typically have much lower temperature, and thus couple more strongly with the substrate, which can be due to their higher transmission to the substrate and weaker coupling to the other hot phonon modes in graphene Vallabhaneni et al. (2016b).

We find that the substrate also breaks the in-plane symmetry of graphene, in addition to the breakdown of its out-of-plane symmetry that was discovered previously Qiu and Ruan (2012). For example, the spectral phonon temperatures of graphene along $\Gamma$-$M$ become quite different from those along $\Gamma$-$M^{\prime}$ when a substrate is applied as shown in Fig. 4 (b). This is because graphene has a sixfold axis of rotation while Si’s is fourfold. The atoms’ vibrations along $\Gamma-M$ in graphene align with those along [0,0,1] in silicon, while those along $\Gamma-M^{\prime}$ in graphene aligns with [0,2,$\sqrt{3}$] in silicon. This finding may be used to engineer the thermal transport between graphene and substrate using alignment.

Moreover, the phonon thermal nonequilibrium near interface creates a new resistance mechanism. Heat needs to be transferred from the optical modes to the acoustic modes in graphene first; then to the longitudinal modes in the substrate; and finally back to other modes in the substrate. This is in analogy to electron-phonon coupled thermal transport across metal-dielectric interfaces, where electron-phonon nonequilibrium in metals introduces an additional resistance Majumdar and Reddy (2004); Wang et al. (2016).

VI Phonon local non-equilibrium in graphene-graphene 2D junction

Here, we investigate a 2D-2D junction, the graphene-graphene junction, which is common in graphene foams that are promising for thermal management and energy storage Zhao et al. (2014); Zhou et al. (2013); Han et al. (2014); Huang et al. (2014); Thiyagarajan et al. (2014). As indicated in Fig. 5 (a), the heat flows from the left side of the top graphene to the bottom graphene with the cross-plane distance 0.335 nm (AB stack), same as in graphite. Inside the graphene layers away from the junction, the MD temperature profile is quite smooth and flat due to the strong sp2 bond and ballistic transport nature. Near the junction, at $x=0$ for the top layer and $x=100$ Å for the bottom layer, there are small jumps, indicating that the interlayer van der Waals field indeed drives the heat flow.

To gain insight of the spectral phonon heat conduction channel, we calculated the modal temperatures as shown in Fig. 5 (b)-(d). First, we inspect the phonons traveling along the $x$ (heat flow) direction Fig. 5 (b). Surprisingly, the phonon local thermal nonequilibrium is quite large throughout the whole system. At the junction the spectral phonon non-equilibrium is due to the different spectral heat transfer efficiency, while away from the junction the spectral non-equilibrium is due to the ballistic effect since the length is much shorter than phonon MFP in graphene. Clearly, the three acoustic phonon branches have much smaller jumps than the optical phonon branches, and they dominate the heat transfer. Among them, ZA mode is the most efficient one, whose the temperature jump is about only 10% of the overall $\Delta T_{MD}$. The out-of-plane vibration of the ZA mode can directly transfer energy across the vertical junction. The optical branches, including the ZO mode, however, cannot transfer heat efficiently, probably due to the negative group velocity with respect to the heat flow direction.

Now we inspect the phonons traveling in the $y$ direction, perpendicular to the heat flow, as shown in Fig. 5 (c). These phonons do not contribute to the thermal transport which is along the $x$ direction. We find a completely different phenomenon compared to that in the thermal transport direction. Here the phonon temperatures match well with the MD temperature except for the junction region. This is because that these phonons are diffusive as we discussed previously, and the phonon non-equilibrium away from the junction does not exist anymore. While at the junction, the spectral non-equilibrium still exist due to the different spectral heat transfer efficiency. Among these phonons, ZA and ZO branches have the smallest temperature jumps, probably due to their out-of-plane vibration’s governing the interlayer interaction. The branch temperatures averaged among 4 representative directions are plotted in Fig. 5 (d).

Furthermore, based on the $\mathbf{k}$-resolved phonon temperature shown in Fig. 6, we find that the longer-wavelength acoustic phonons have lower temperature in the upper layer and higher temperature in the lower layer, while the optical phonons are opposite. This indicates that the longer-wavelength acoustic phonons contribute more to the junction heat transfer.

VII Phonon local non-equilibrium in graphene-boron nitride 2D interface

After inspecting the dimensionally mismatched interfaces and matched junctions, now we take the graphene/boron nitride interface Fig. 7 (a) as an example to examine the 2D planner interfacial thermal transport. Surprisingly, even across such a nearly perfect-lattice-matched interface, the phonon local thermal nonequilibrium is quite large as shown in Figs. 7 (b) and (c). The coupling between the acoustic phonon modes ZA, TA and LA is quite strong, and the temperature jump of the ZA mode is only 10-20% of the overall $\Delta T_{MD}$. One interesting phenomenon is that the temperature jump of the ZO mode is as high as 220% of the $\Delta T_{MD}$ although it vibrates in the out-of-plane direction as the ZA mode. Such large jump and thus inefficient heat conduction are due to the negligible phonon band overlap between graphene and BN as shown in Fig. 8, which leads to negligible elastic transmission. Furthermore, if interface mixing is introduced, $\Delta T_{MD}$ as well as all the phonon modal temperature jumps $\Delta T_{\lambda}$ increase. Interestingly, relatively the jump of the ZA mode increases the most, from 19% to 41%, due to the break of the out-of-plane symmetry, while the relative jumps of the in-plane modes decreases, demonstrating that the heat transfer channel shifts a certain amount from ZA mode to the TA and LA modes.

VIII Spectral phonon thermal conductivity

With our SPT method, the spectral phonon thermal conductivity $\kappa_{\lambda}$ that is defined by Fourier’s law

[TABLE]

can be obtained, where the spectral heat flux $q_{\lambda}$ can be calculated by the method developed by Zhou and Hu Zhou and Hu (2015) recently.

[TABLE]

Here $E_{l,b}$ and $\mathbf{S}_{l,b}$ are the total energy and stress of the atom (l,b), respectively. We take 26-nm and 152-nm Si thin films as examples, and calculate the spectral phonon heat flux and temperature gradient along $\Gamma$-X as shown in Fig. 9. It is clearly seen that the TA and LA branches carry the most amount of heat (Fig. 9 a, c) although their temperature gradients are the smallest (Fig. 9 b, d). This certainly indicates that their spectral thermal conductivities are the largest. Generally, the heat fluxes of all the acoustic and optical branches decrease with increasing phonon frequency, which agree well with decreasing trend of phonon relaxation time with frequency. At X point, the modal heat flux $q_{\rm LA}$ merges with $q_{\rm LO}$, and $q_{\rm TA}$ merges with $q_{\rm TO}$. The former is due to the degeneration of the LA and LO branches at X point, and the latter is due to the zero group velocity of both TA and TO branches at X point. Overall, the error bars of the results are small, considering that molecular dynamics is statistical. Relatively, the uncertainty of optical phonons are larger than that of acoustic phonons due to the diffusive nature. But this is tolerable since optical phonons only contribute a small portion to the total thermal conductivity.

To further verify the accuracy of our results, we calculated the modal heat flux and temperature gradient in the whole octant of the Brillouin Zone with $16\times 16\times 16$ k-mesh as shown in Fig. 10. We have found that the total heat flux $\sum_{\lambda}q_{\lambda}$ and the average temperature gradient $\nabla\bar{T}_{\lambda}$ agree excellently with the MD heat flux $q_{\rm MD}$ and temperature gradient $\nabla T_{\rm MD}$, respectively.

The conventional thermal conductivity calculated from the grey temperature $T_{\rm MD}$

[TABLE]

misses the spectral phonon temperature information, and thus underestimates the thermal conductivity contribution of the ballistic modes. For the 152-nm Si thin film, Eq. (12) yields 148 W/mK based on the modal heat flux$\sum_{\lambda}q_{\lambda}$ which agree well with the 151 W/mK calculated from direct MD heat flux $q_{\rm MD}$. After considering the spectral phonon temperature, the thermal conductivity

[TABLE]

is obtained as 178 W/mK which eliminates the size effect induced from ballistic transport. We note that it is still lower than the bulk value 235 W/mK determined by the NEMD length extrapolation based on Tersoff potential. This is probably due to the reflection at the boundaries near the reservoirs.

IX Conclusions

In summary, we have developed a powerful SPT method to extract the spectral phonon temperature directly from atomistic simulations, and demonstrated its capability of resolving ballistic and diffusive phonon modes in nanomaterials and spectral phonon mode coupling across interfaces. In the benchmark materials, nano-size Si and graphene, the thermal local nonequilibrium between the ballistic and diffusive phonons has been clearly observed. Such nonequilibrium exists in interfaces as well and is surprisingly large, based on our study of graphene/substrate, graphene/graphene junction and graphene/BN planar interface. The phonon local thermal nonequilibrium introduces a new mechanism of thermal interfacial resistance. In particular, the most efficient thermal transport channel across the dimensionally-mismatch graphene-substrate interface is the transmission of in-plane acoustic modes in graphene into the cross-plane LA and LO modes in the substrate. The interface roughness is shown to substantially affect such nonequilibrium and coupling, and shift the heat transfer channel. The SPT method together with the spectral heat flux can predict the thermal conductivity with the size effect induced from ballistic transport eliminated. Since our method can directly extract the spectral and spatial phonon temperature from MD simulations, it is expected to have broad applicability to many thermal applications such as thermal diodes, thermal interface materials (TIMs), coherent phononic meta-materials, etc.

Acknowledgements.

Simulations were preformed at the Rosen Center for Advanced Computing (RCAC) of Purdue University. The work was partially supported by the National Science Foundation (Award No. 1150948) and the Defense Advanced Research Projects Agency (Award No. HR0011-15-2-0037). We would like to thank Yanguang Zhou for helpful discussions on the spectral heat flux decomposition.

Appendix A Benchmark the method in equilibrium molecular dynamics

This section is to (1) determine the total simulation time to obtain good modal temperature results, (2) determine a better approach between the long-time simulation and the ensemble average of short-time simulations since the latter is more time-efficient, and (3) verify that our spectral phonon temperature method does not require periodic boundary conditions that is required for the normal mode analysis to calculate phonon relaxation time.

Appendix B Formalism of the phonon temperature averaging

The kinetic energy of the atoms is the $3N$ degrees of freedom times the average energy each degree of freedom $k_{B}T/2$

[TABLE]

The kinetic energy of the 3N phonon modes is

[TABLE]

Since the two total energies should be equal to each other, the averaged phonon temperature by combining Eqs. (14) and (15) is

[TABLE]

Since practically only finite number of phonon modes are calculated, and thus the integral is evaluated by a discrete form:

[TABLE]

Appendix C Phonon mean free path spectra obtained by spectral energy density analysis

The spectral phonon mean free path is obtained by $\Lambda_{\lambda}=v_{\lambda}\tau_{\lambda}$, where $v_{\lambda}$ is the group velocity, and $\tau_{\lambda}$ is the phonon relaxation time. $\tau_{\lambda}$ is obtained by performing the following NMA Dove (1993); Qiu and Ruan (2011); Feng and Ruan (2014); Feng et al. (2015) based on equilibrium MD simulations,

[TABLE]

$\mathcal{F}$ denotes the Fourier transformation. The spectral energy density $E_{\lambda}(\omega)$ of the phonon mode $\lambda$ is obtained by substituting $\dot{u}_{\alpha}^{l,b}(t)$ extracted from MD trajectory into Eq. (20), where $C_{\lambda}$ is a constant for a given $\lambda$. By fitting the spectral energy density as a Lorentzian function, the peak position $\omega_{\lambda}^{A}$ and full linewidth $\tau_{NMA,\lambda}^{-1}$ at half maximum are obtained. Our former work Feng et al. (2015) has shown that Eq. (20) are equivalent to another version of frequency-domain NMA that does not include phonon eigenvectors de Koker (2009); Thomas et al. (2010)

[TABLE]

A full discussion about the methods of predicting phonon relaxation time was given in Ref. Feng and Ruan (2014).

Figure C.1 (a) shows a sample plot of the spectral energy density of Si. The position of the peaks indicate phonon frequencies $\omega_{\lambda}^{A}$, and the full width at the half maximum (FWHM) indicates the reciprocal of relaxation time $\tau_{\lambda}^{-1}$. The Figs. C.1 (b) and (c) are the phonon MFP as a function of frequency in silicon and graphene, respectively.

Appendix D Spectral Phonon Temperature of Silicon in Non-equilibrium molecular dynamics under Berendsen thermostat

Appendix E Spectral Phonon Temperature of Graphene under Langevin thermostat