Using non-equilibrium molecular dynamics simulations, we investigate thermal transport in organic/inorganic Au/molecular nanolayer (MNL) nanolaminate. We examine the tunability of thermal conductivity via interfacial bonding by (i) homogeneous change of bonding strength and heterogeneous change of (ii) bond density and (iii) molecular coverage at the interface. By comparing the thermal conductivity of the nanolaminates with the interfacial thermal conductance of corresponding individual interfaces, we conclude that phenomenologically the thermal conductivity can be predicted from independent interfacial resistors connected in a series model, particularly at higher temperatures. However, interfacial thermal conductance shows a moderate increase with temperature, whereas the thermal conductivity of Au/MNL nanolaminates shows the opposite effect. We elucidate this apparent contradiction via phonon wave packet simulations at individual and multiple interface structures.

Nanolaminates are the class of composite materials assembled from alternating layers of two or more materials at the nanoscale. Nanolaminate composites display thermal, electrical, mechanical, and optical properties not exhibited by their bulk counterparts.1 Organic/inorganic nanolaminates in particular combine and exploit material properties across both the organic and the inorganic domains. Organic/inorganic nanolaminates such as polyelectrolytes/silicate display an approximately twofold increase in fracture toughness,2 while Al2O3/molecular nanolayer (MNL) nanolaminates display high (∼95%) optical transmittance and low water permeability making it an excellent candidate for organic light-emitting diodes (OLEDs),3 and SiO2/SiOCH nanolaminates are found to have improved scratch resistance tunable with interface density.4 Our prior work on viscoelastic damping characterization on the Au/MNL multilayer also revealed frequency-dependent viscoelastic damping bandgaps not seen either on the bulk crystal of Au or polyethylene crystals, bearing potential applications in nanomechanical dampening.5 

Interfacial thermal conductance and thermal transport across the organic/inorganic interface with monolayers have been well studied both computationally and experimentally in the last two decades.6–8 The interfacial thermal conductance in organic/inorganic interfaces are known to be tunable with interfacial properties such as bond stiffness, bond strengths, bond density, molecular coverage, and MNL chain lengths. Following the miniaturization of molecular electronics in recent years, the scientific community has been focusing on thermal transport across single molecular junctions (SMJs) extensively.6,7 The bulk thermal conductivity response in inorganic/inorganic interfaces9 and their superlattices10 are also well studied and therein thermal properties are often modulated by varying geometric parameters like superlattice periods (λ). In regard to thermal transport in organic/inorganic nanolaminates, however, there are only a few experimental studies in the literature1,11–14 and these works focus on the effect of interfacial density on the bulk thermal conductivity in the material. Consequently, there is still a knowledge gap in the understanding of the thermal transport mechanism in organic/inorganic nanolaminates and its dependence on interfacial bonding character. Thermal transport at individual inorganic/inorganic9,15–17 and organic/inorganic interfaces18–21 was well studied using molecular dynamics (MD) simulations. Also, there is a large body of work on inorganic–inorganic superlattices22–24 motivated by their thermoelectric applications and the development of ultralow thermal conductivity materials.25 However, thermal transport mechanisms in organic/inorganic nanolaminates are less explored. Here, we strive to understand the thermal transport response of organic/inorganic nanolaminates from the interfacial point of view and to assess how bulk thermal response can be tuned using interfacial properties.

Our approach is to tune interfacial bonding in nanolaminates both homogeneously and heterogeneously and determine the bonding-thermal conductivity relationship. In addition, we separately determine interfacial thermal conductance and its bonding dependence at corresponding individual interfaces. This allows us to examine the question of to what extent thermal transport in nanolaminates can be understood in terms of independent thermal interfacial resistors connected in a series model, or, alternatively, if collective interfacial response/coherent phonon transport across multiple interfaces is a significant factor in understanding thermal transport of these nanolaminates. To shed more light on this question, we also performed the phonon wave packet (WP) simulations on Au/MNL single and multiple interface junctions to determine phonon transmission and details of the interfacial phonon scattering.

The rest of the article is organized as follows: in Sec. II, we describe model structures and the computational details of thermal transport characterization by Non-Equilibrium Molecular Dynamics (NEMD) and phonon wave packet (WP) simulations. Section III presents results and discussion on the effect of interfacial bond strength (XIBS), interfacial bond density (XIBD), and interfacial molecular coverage (XIMC) on the thermal conductivity of the Au/MNL nanolaminates and the interfacial conductance of individual interfaces. In Sec. III, we also present the results on the effect of temperature on Au/MNL nanolaminate thermal conductivity and Au/MNL/Au junction interfacial thermal conductance. In addition, in Sec. III, we provide insight into our NEMD results with the results of phonon wave packet simulations. In Sec. IV, we summarize the major conclusions drawn in this study.

In the present study, the thermal transport properties of Au/MNL nanolaminates and Au/MNL/Au junction were examined. First, a two-bilayer unit-cell with each layer consisting of six atomic layers of Au(111) layer of thickness 11.75 Å covalently bonded with alkane dithiol [–S–(CH2)n–S–] MNL of thickness 9.53 Å was created using a Moltemplete builder.26 In the first layer, MNL molecules with the thiol end-group are anchored to ( 3 × 3 ) R30° surface of Au(111); the specific anchoring method is detailed elsewhere.18,19 In the second layer, alkane dithiol MNL molecules are rotated 120 ° anticlockwise relative to the MNL in the first layer to create a stable bilayer structure. We then replicated the bilayer unit-cell to create Au/MNL nanolaminates with a cross section of 51 × 48.9 Å 2 for studying thermal conductivity (see Fig. 1). For the determination of the interfacial thermal conductance, we used structures with only two-bilayer units but with significantly thicker Au and MNL layers [see Fig. 2(b) and details below].

FIG. 1.

Model structure of the nanolaminate consisting of b = 12 bilayers (∼307 Å long) of Au/MNL with a zoomed-in view of interfacial bonding between Au and S. The interfacial properties such as bond strength (XIBS), bond density (XIBD), and molecular coverage (XIMC) are highlighted.

FIG. 1.

Model structure of the nanolaminate consisting of b = 12 bilayers (∼307 Å long) of Au/MNL with a zoomed-in view of interfacial bonding between Au and S. The interfacial properties such as bond strength (XIBS), bond density (XIBD), and molecular coverage (XIMC) are highlighted.

Close modal
FIG. 2.

Steady-state temperature profile of (a) Au/MNL nanolaminates (b = 12) and (b) the Au/MNL/Au junction (∼54 Å thick Au/ ∼ 34 Å thick MNL) in a NEMD simulation.

FIG. 2.

Steady-state temperature profile of (a) Au/MNL nanolaminates (b = 12) and (b) the Au/MNL/Au junction (∼54 Å thick Au/ ∼ 34 Å thick MNL) in a NEMD simulation.

Close modal

In our modeling approach, the hydrogen atom was not modeled explicitly, but its mass was added to the carbon backbone, treating the methylene group (CH2) as a united atom. For simplicity, the methylene group will only be referred to as “C” for the rest of the article. The non-bonded van der Waals interactions were modeled using 6–12 Lennard-Jones (LJ) forcefield for C–C, S–S, Au–C, C–S, and Morse forcefield for Au–Au and Au–S pairs.18 For both LJ and Morse forcefield, cutoffs are typically chosen at 2.5ro (ro is the equilibrium separation), where the interactions are already very weak and interaction energy essentially goes to zero. In our case, the 2.5ro results in ∼11 and ∼8 Å as a cutoff distance for LJ and Morse forcefield, respectively. We used the harmonic pairwise bonded interactions between C–C, C–S; harmonic bond angle bending interactions between C–C–C, S–C–C;18 and Fourier-type dihedral angle interactions between C–C–C–C and S–C–C–C.27 To study the bonding effect on thermal transport, we (i) homogeneously tuned the Au–S interfacial bond strength (XIBS) via the bond energy parameter (Do) in the Morse forcefield from its full strength down to 10% of the full strength, or heterogeneously, (ii) either by turning off Morse forcefield bonding for a fraction of randomly selected bonds and keeping the rest of the bonds at full strength (interfacial bond density, XIBD) or by interfacial molecular coverage (XIMC) tuned by randomly removing the whole alkane dithiol MNL chain molecule sandwiched between the inorganic Au layers in both Au/MNL nanolaminates and Au/MNL/Au junctions (see Fig. 1).

We used non-equilibrium MD (NEMD) simulation with a constant heat flux method for computing thermal transport properties with the open-source Large-Scale Atomic/Molecular Massively Parallel Simulator (LAAMPS) package.28 At first, we performed an isothermal–isobaric (NPT) ensemble equilibration for 1 ns at constant pressure (0 bar) and constant targeted temperature followed by a microcanonical ensemble (NVE) equilibration for 1 ns. Then, we used the NVE ensemble for the whole system, and we added and removed heat at a constant rate from six atomic layers of Au to create the localized heat source and sink regions in the material by scaling the kinetic energy of the Au atoms (see Fig. 2). After ∼2 ns, we observed the development of a steady-state temperature profile in our model structures. The reported temperature profiles were obtained by averaging the kinetic energy of atoms from each atomic plane for the next 3 ns for the production run on both Au/MNL nanolaminates [see Fig. 2(a)] and Au/MNL/Au junction [see Fig. 2(b)]. All our NEMD simulations were performed using periodic boundary conditions (PBC) in all directions and with a time step of 0.5 fs. In all our simulations, heat flux was tuned such that we maintained the source/sink temperature difference of 10%–15% of the average absolute temperature of the whole structure.

The thermal conductivity of Au/MNL nanolaminates was calculated using the relation κ = ( d Q / d t ) / ( 2 A T ), where d Q / d t is the heart rate, T is the temperature gradient, A is the  cross-sectional area of the interface, and the factor of 2 indicates the fact that while using periodic boundary conditions, the total heat flux from the localized heat source crosses the two interfaces with a total cross section of 2A (see Fig. 2). The temperature gradient was estimated using linear regression on the temperature profile of Au/MNL nanolaminate between the heat source and the sink [see Fig. 2(a)]. The thermal conductance on the Au/MNL/Au junction was calculated as G = ( d Q / d t ) / ( 2 A Δ T ), where ΔT is the temperature drops at the interface. During the steady-state temperature gradient develops within the Au block while no obvious temperature gradient develops within the MNL chains, indicating ballistic phonon transport in MNLs [see Fig. 2(b)]. We implemented a linear regression on the temperature profiles of MNL and the Au blocks and computed the ΔT at the interface. Temperature data immediately adjacent to the interface were excluded from the linear regression fitting [see Fig. 2(b)].

Thermal conductivity and thermal conductance obtained from MD simulations can depend on the size of the model structures when the phonon-mean free path is comparable to the system size. We performed a size-effect study to identify the optimal system size for computing the thermal conductance for the Au/MNL/Au junction and thermal conductivity for Au/MNL nanolaminates (see Table I). For studying the size effects of thermal conductivity on Au/MNL nanolaminates, we used various numbers of bilayers b = 6, 8, 12, each bilayer consists of 11.75 Å thick Au layer and 9.53 Å thick MNL and 4.35 Å thick Au–S interface. To calculate thermal conductance, we used Au/MNL/Au junctions with varying thicknesses of the organic and inorganic phases in two bilayers for the structure. Both interfacial thermal conductance and thermal conductivity were computed for limiting bonding cases, i.e., weak interface (XIBS = 0.1) and strong interface (XIBS = 1) at 300 K. Based on the results presented in Table I, we used a bilayer of ∼54 Å thick Au/ ∼ 34 Å thick MNL junction for individual interface investigation and 12 bilayers thick (∼308 Å long nanolaminate) nanolaminate, which is sufficient to avoid significant finite-size effects in our estimates of thermal conductivity and interfacial thermal conductance for the rest of our study.

TABLE I.

Finite-size effects on interfacial thermal conductance (G) of the Au/MNL/Au junction and thermal conductivity (κ) of nanolaminates.

Interfacial thermal conductance (G) (MW m−2 K−1)
Au/MNL/Au junctionWeak interface (XIBS = 0.1)Strong interface (XIBS = 1)
∼26 Å thick Au/∼18 Å thick MNL 36.8 431.8 
∼40 Å thick Au/∼26 Å thick MNL 37.3 414.7 
∼54 Å thick Au/∼34 Å thick MNL 38.6 398.1 
Thermal conductivity (κ) (Wm−2 K−1
Nanolaminates b = no. of bilayers Weak interface Strong interface 
(1 b = 11.75 Å Au, 9.53 Å MNL, 4.35 Å Au–S interface) (XIBS = 0.1) (XIBS = 1) 
b = 6, Total length = ∼ 154 Å 0.061 0.667 
b = 8, Total length = ∼ 205 Å 0.063 0.698 
b = 12, Total length = ∼ 308 Å 0.068 0.711 
Interfacial thermal conductance (G) (MW m−2 K−1)
Au/MNL/Au junctionWeak interface (XIBS = 0.1)Strong interface (XIBS = 1)
∼26 Å thick Au/∼18 Å thick MNL 36.8 431.8 
∼40 Å thick Au/∼26 Å thick MNL 37.3 414.7 
∼54 Å thick Au/∼34 Å thick MNL 38.6 398.1 
Thermal conductivity (κ) (Wm−2 K−1
Nanolaminates b = no. of bilayers Weak interface Strong interface 
(1 b = 11.75 Å Au, 9.53 Å MNL, 4.35 Å Au–S interface) (XIBS = 0.1) (XIBS = 1) 
b = 6, Total length = ∼ 154 Å 0.061 0.667 
b = 8, Total length = ∼ 205 Å 0.063 0.698 
b = 12, Total length = ∼ 308 Å 0.068 0.711 

1. WP simulation method

Phonon wave packet simulations were first implemented by Schelling et al.29 to study interfacial phonon scattering directly via MD simulations. This method has been used to study phonon transport across all kinds of interfaces such as superlattices,30 crystalline/amorphous interfaces,20 solid/liquid interfaces,31 and heterostructures.32 In this work, we implement phonon wave packet simulations to gain insight into the thermal transport mechanism in the Au/MNL single and multiple interfaces at the individual phonon level.

In order to perform WP simulations, the phonon wave packet of well-defined frequency and polarization is injected into the relaxed gold lead and launched toward the interface. This procedure requires first to calculate phonon spectra and associated eigenvectors. Phonon dispersion of gold was computed using the lattice-dynamics (LD) calculation with the General Utility Lattice Program (GULP)33 package [see Fig. 3(a)]. The phonon dispersion of gold crystal with atomic interactions described by the Morse potential compares reasonably well with that obtained in the neutron scattering experiment on gold;34 however, the frequencies are overestimated by the factor of ∼1.4 throughout the Brillouin zone high symmetry direction. Although there are force fields that can capture dispersion relation accurately for gold, we used the Morse forcefield to be consistent with our earlier NEMD simulation on gold–alkane interfaces.5 In this work, the Γ–L phonon band of the face-centered cubic (FCC) Au crystal was used to launch longitudinal acoustic (LA) phonon wave packets, with phonons normally incident on the Au (111) surface. The normally incident phonons have group velocity orthogonal to the surface and have a higher energy-transmission coefficient compared to other phonon branches.31,35

FIG. 3.

(a) The phonon dispersion relations of FCC Au crystal with the Morse potential along the Γ–X–W–X–K–Γ–L path in the Brillouin zone calculated using lattice-dynamics calculation in GULP. L refers to longitudinal acoustic, and T refers to transverse acoustic modes in the phonon dispersion curve. (b) Snapshots of the spatial distribution of velocity for an LA mode phonon wave packet with ν = 3 THz on the Au/MNL single interface. (c) Snapshots of the spatial distribution of velocity for an LA mode phonon wave packet with ν = 3 THz on Au/MNL with four interfaces. Spectral probes located at ∼300–400 nm away from the interface indicate the positions where we monitor atomic velocities to provide data for frequency analysis for reflected and transmitted phonons.

FIG. 3.

(a) The phonon dispersion relations of FCC Au crystal with the Morse potential along the Γ–X–W–X–K–Γ–L path in the Brillouin zone calculated using lattice-dynamics calculation in GULP. L refers to longitudinal acoustic, and T refers to transverse acoustic modes in the phonon dispersion curve. (b) Snapshots of the spatial distribution of velocity for an LA mode phonon wave packet with ν = 3 THz on the Au/MNL single interface. (c) Snapshots of the spatial distribution of velocity for an LA mode phonon wave packet with ν = 3 THz on Au/MNL with four interfaces. Spectral probes located at ∼300–400 nm away from the interface indicate the positions where we monitor atomic velocities to provide data for frequency analysis for reflected and transmitted phonons.

Close modal
The wave packet with the wave vector ko and phonon branch γ was formulated and launched toward the interface by displacing the atoms in the gold crystal lead according to Ref. 29,
u j s μ = A ϵ j μ γ ( k o ) e ( i k o ( z s z 0 ) ) e ( z s z o ) 2 / ζ 2 ,
(1)
where u j s μ represents the μth component of displacement for atom j in the unit-cell s, A is the amplitude of the wave packet, ϵ j μ γ ( k o ) is the polarization vector along the band γ at ko, zo is the center of the wave packet, and 1/ζ is the spatial extent of the wave packet. The time-dependent displacement of atoms was computed as
u j s μ ( t ) = A ϵ j μ γ ( k o ) e ( i k 0 ( z s z o ) ) e ( z s z o ) 2 / ζ 2 e i ω t ,
(2)
where ω (=2πν) is the angular frequency of the phonon. The velocity of the wave packet was calculated using
v j s μ ( t ) = d u j s μ ( t ) d t = i ω u j s μ ( t ) .
(3)

Using Eq. (2) atoms are displaced on the gold lead to formulate the wave packet and using Eq. (3) initial velocity was assigned to the displaced atoms. These initial sets of atomic positions and velocities were utilized to launch phonon wave packets toward the individual or multiple Au/MNL interfaces (see below).

2. Model structures for WP simulations

The wave packet simulation on a single Au/MNL interface was performed on the model structure consisting of the 3600 Au (111) layers thick Au lead of dimension 17.3 × 20 × 8457 Å3 covalently bonded with –S–C7000–S– MNL chains [see Fig. 3(b) top panel]. Similarly, WP simulation on multiple interfaces was performed on Au/MNL nanolaminate with n = 4 interfaces (simply referred to as multiple interfaces henceforth). The 3600 Au (111) layers thick Au lead was attached to Au/MNL multiple interfaces of dimension 16.93 × 19.53 × 51.15 Å3 on top and bottom along the z-axis [see Fig. 3(c) top panel]. These model structures were first equilibrated with an NPT ensemble at 0.1 K for 0.8 ns using the Berendsen thermostat and Parrinello–Rahman barostat followed by potential energy minimization using the conjugate gradient (CG) algorithm with the Parrinello–Rahman barostat and damped dynamics minimization method within the LAMMPS platform. The combination of these minimization algorithms is performed on model structures with 10−15 eV Å−1 force tolerance and a maximum iteration of 6 × 105 steps. Following this protocol, we obtain zero temperature relaxed baseline structures. Then, the wave packet with a well-defined wave vector and polarization was introduced into the left Au lead and launched toward the Au/MNL single and multiple interfaces [see Figs. 3(b) and 3(c)]. We performed all our wave packet simulations in multiple LA modes along the Γ–L phonon band with displacement amplitude A = 0.01 Å, lattice parameter along Au (111) a o / 3 with a o = 4.079 Å, and spatial extent ζ = 80 a o for Au/MNL single interface and ζ = 160 a o for Au/MNL multiple interfaces.

Figures 3(b) and 3(c) show an example of the z-component of the spatial velocity distribution of a wave packet for single and multiple interfaces. The wave packet is launched toward the interface and upon scattering at the interface, part of the energy of the wave packet is reflected and transmitted [see Figs. 3(b) (blue curve) and 3(c) (red curve)]. We then monitored the transmitted and reflected energy to calculate the energy-transmission coefficient [α(t)] and reflection coefficient [1 − α(t)] on single and multiple interfaces. As we discuss in more detail in Sec. II C 3, we calculate the transmission and reflection coefficients once the interfacial transient scattering phase is completed. When calculating α(t), the baseline energy (without wave packet) of the model structures at 0 K was subtracted from the model structures with the wave packet to remove any artifacts resulting from the atomic relaxation. Transmission coefficient α(t) = 0 indicates total reflection of the wave packet energy, α(t) = 1 indicates complete transmission of energy, and 0 < α(t) < 1 indicates both transmission as well as the reflection of the wave packet.

3. Power spectrum

To analyze the wave packet scattering at individual and multiple interfaces, we determined the reflected and transmitted frequencies of the phonon wave packets after the interfacial scattering event is completed. We monitored the reflected (EReflected), transmitted (Etransmitted), and interface (EInterface) energy throughout the scattering process and concluded the completion of the scattering process after these energies do not change with time [see Fig. 4(a) at 130 ps for the case of ν = 0.67 THz], and there is no energy stored at the interface. For frequency analysis, the x, y, and z components of the velocity of selected atoms residing far away from the interface (300–400 nm) were monitored for reflected and transmitted waves. Examples of such velocities temporal dependence are shown in Figs. 4(b) and 4(c). Quite interestingly we observed cases where after scattering the incident LA wave can in part transform to a TA wave [Fig. 4(e)].

FIG. 4.

Longitudinal acoustic (LA) wave packet with ν = 0.67 THz scattering at the Au/MNL multiple interface structure. (a) Reflected (EReflected), transmitted (ETransmitted), and interface (Einterface) energy during the wave-packet scattering process. (b) and (c) Velocity of reflected and transmitted wave packet after the scattering event at ∼300 and ∼400 nm away from the interface for the reflected and transmitted wave, respectively. (d) and (e) Time Fourier transforms of reflected and transmitted wave-packet velocity signal in (b) and (c).

FIG. 4.

Longitudinal acoustic (LA) wave packet with ν = 0.67 THz scattering at the Au/MNL multiple interface structure. (a) Reflected (EReflected), transmitted (ETransmitted), and interface (Einterface) energy during the wave-packet scattering process. (b) and (c) Velocity of reflected and transmitted wave packet after the scattering event at ∼300 and ∼400 nm away from the interface for the reflected and transmitted wave, respectively. (d) and (e) Time Fourier transforms of reflected and transmitted wave-packet velocity signal in (b) and (c).

Close modal

In order to identify the frequencies of reflected and transmitted wave packets, we performed the Fast Fourier Transform (FFT) of the time-domain velocity signal [see Figs. 4(b) and 4(c)] with a sampling interval of 0.0125 ps, FFT decomposes the signal into its constituent frequency components. The Discrete Fourier Transform (DFT) routine implemented in the NumPy library37 was utilized in this work for performing the FFT of time-domain velocity signal, the details of the FFT algorithm implemented in the DFT routine can be found in Ref. 36. Power spectrum F(ν) was then calculated by squaring the absolute value of the FFT result. The frequencies of the reflected ν r e f l e c t e d i and the transmitted ν t r a n s m i t t e d i wave packets are identified from the peak positions on the power spectrum of reflected/transmitted wave packets.

As shown in Figs. 4(d) and 4(e), scattering leads to rich outcomes including conversion of the initial LA phonon into TA phonon, as well as clear anharmonic effects manifesting in frequency doubling and tripling. We further verified (not shown) that the anharmonic effects diminish with the lowering of the incident phonon amplitudes. We want to emphasize that our frequency analysis methodology implemented here identifies the scattered wave packet frequencies unencumbered from any transient phonon scattering effects at or near the interface due to spatial and temporal separation for the transient interfacial scattering.

First, we studied the effect of a homogeneously changing interfacial bonding by varying the Au–S interfacial bond strength (XIBS) with the thermal conductivity κIBS of Au/MNL nanolaminate. The computed κIBS is shown in Fig. 5 (blue curve). In our simulations, we controlled XIBS homogenously across all the interfaces by lowering the bond strength parameter (Do) of Au–S in the Morse potential. Experimentally, such variations are possible when the end-group (CH3, NH2, Br, SH) chemistry is altered across all the interfaces from the strong bonding (SH) to the weakly bonding (CH3) group.38 In our simulations, as we increased the strength 0.1 ≤ XIBS≤ 1, the interface became more conducting, and thermal conductivity increased from 0.07 ≤ κIBS ≤ 0.7 W/m K. Thermal conductivity κIBS increased with increasing XIBS displaying positive curvature. In the literature, organoclay nanolaminates with a nitrogen atom head group (weak interfacial bonding) are experimentally reported to have an ultralow thermal conductivity of ∼0.1 W m−1 K−1.11 Our simulation result at low bond strengths (XIBS = 0.1) in Au/MNL nanolaminate is comparable to the reported value.

FIG. 5.

Thermal conductivity, κ, of Au/MNL nanolaminates for different interfacial Au–S bond strengths (XIBS), bond density (XIBD), and molecular coverage (XIMC).

FIG. 5.

Thermal conductivity, κ, of Au/MNL nanolaminates for different interfacial Au–S bond strengths (XIBS), bond density (XIBD), and molecular coverage (XIMC).

Close modal

Second, we studied the effect of heterogeneously changing interface on the thermal conductivity κ of Au/MNL nanolaminate. The heterogeneous change was implemented by varying interfacial (i) bond densities, XIBD, and (ii) molecular coverage, XIMC. Interfacial bond density (XIBD) is manipulated by removing the bond between Au–S altogether for randomly selected bonds, which mimics the fraction of poor MNL contact at the interface. Experimentally, such variation is possible in a mixed MNL configuration where XIBD is controlled by the mole fractions of alkane dithiol and monothiol MNL chains sandwiched between Au substrates.39 Molecular coverage is simply varied by removing the MNL chain molecules in between the Au layers in a random fashion. Thermal conductivity as a function of bonding density within the 0.1 ≤ XIBD≤ 1 range is shown in Fig. 5 (red curve). As in the case of the homogenous bonding variation, the thermal conductivity is tunable from 0.07 ≤ κIBD ≤ 0.7 W m−1 K−1. However, simulated κIBD values increase parabolically with the negative curvature, while κIBS increases with the positive curvature. This result indicates that provided that the average bonding strength is the same, homogenously bonded interfaces are more conductive than heterogeneously bonded interfaces.

Finally, we computed the thermal conductivity of Au/MNL nanolaminate, κIMC with varying interfacial molecular coverage only in the range 0.3 ≤ XIMC ≤ 1 as shown in Fig. 5. We found that Au/MNL nanolaminate was structurally unstable for XIMC < 0.3 during the NPT equilibration. Figure 5 shows that κIMC follows κIBD very closely, indicating that the exact nature of the bonding heterogeneity is of secondary importance and thermal conductivity is mainly controlled by the average bonding strength.

While the thermal conductivity of nanolaminates increases with increasing XIBS/XIBD/XIMC, homogeneously bonded interfaces (XIBS) show a faster increase in thermal conductivity compared to heterogeneously bonded interfaces (XIBD/XIMC, see Fig. 5). This could potentially be due to higher phonon scattering at the heterogeneous interface, leading to more bond defects compared to the homogeneous interface.

In order to understand the thermal transport behavior of Au/MNL nanolaminate presented in Sec. III A, we can construct a simplest possible model based on thermal resistors (R) connected in series as follows:
d κ = R = A u / M N L ( R A u + R M N L + R i n t ) ,
(4)
where κ is the thermal conductivity of Au/MNL nanolaminate, d is the average interfacial distance in the Au/MNL nanolaminate, RAu is the thermal resistance of Au, RMNL is the thermal resistance of MNL, and Rintis the thermal resistance of individual Au/MNL/Au junction. In the case of Au, κAu = 318 W m−1 K−1 with both electronic and phononic contribution.40 For an nm thick Au block, RAu is of the order of ∼10−12 m2 K W−1, very negligible compared to Rint, which is of the order of ∼10−4 m2 K W−1. Alkanethiol MNL acts as ultrafast thermal energy transport channels conducting heat at the rate of 50 pW K−1 per chain,41 offering negligible resistance to the heat flow. So, we can neglect RAu and RMNL terms in our series resistance model presented in Eq. (4). Consequently, we can estimate thermal conductivity via
κ i n t d G ,
(5)
where G is interfacial thermal conductance (the inverse of interfacial thermal resistance). Equation (5) predicts the thermal conductivity of Au/MNL nanolaminates using only the interfacial thermal conductance, an outcome of the series thermal resistance model.

To evaluate the applicability of the model, we need to determine the interfacial thermal conductance, G, for the individual Au/MNL/Au junction, compute the interface predicted thermal conductivity κint of Au/MNL nanolaminate, and compare such values with those obtained in direct simulations on nanolaminates (Fig. 5).

Figure 6 shows the results on interfacial thermal conductance, G, in the Au/MNL/Au junction as a function of interfacial bonding strength, XIBS/XIBD/XIMC. Our results on interfacial thermal conductance for a homogeneously varying interfacial bonding, GIBS, agree well with the trends on MD simulation results of GIBS for Cu/MNL/SiO2 junction42 and Si/Polymer junction.20 The Au/MNL/Au junction thermal conductance GIBS for homogeneously varying interface shows an increase with XIBS with positive curvature that is qualitatively the same behavior as seen on κIBS in Au/MNL nanolaminate (see Fig. 5).

FIG. 6.

Interfacial thermal conductance, G, of the Au/MNL/Au junction for different interfacial Au–S bond strengths (XIBS), bond density (XIBD), and molecular coverage (XIMC).

FIG. 6.

Interfacial thermal conductance, G, of the Au/MNL/Au junction for different interfacial Au–S bond strengths (XIBS), bond density (XIBD), and molecular coverage (XIMC).

Close modal

Experimental work on similar interfaces was done by Majumdar et al.43 and Losego et al.38 Reference 43 reports the experimentally measured interfacial conductance of the Au/MNL/Au junction of 65 ± 7 MW m−2 K−1. This interface was characterized by rather poor Au/MNL bonding estimated to be at about 50% of the perfect interface value. Nevertheless, these values are significantly lower than those reported in Fig. 6 even at 50% bonding strength. Similarly, Losego et al.38 also determined the interfacial conductance of Au/X–C11–Si/Quartz interfaces, where X = CH3, SH, Br, NH2 to be in the range of 20–70 MW m−2 K−1 depending on the interfacial bonding strength.

There are potentially multiple reasons why the values reported in our simulations are significantly larger than in the experiments. First of all, in both the experimental studies the measured conductance characterizes two interfaces connected in series, whereas in Fig. 6, we report the interfacial thermal conductance for a single Au/MNL interface. Second, our model interfaces are perfectly planar, while Majumdar et al.43 indicated that experimentally prepared surfaces are rough and exhibit imperfect bonding as described above. Finally, we need to recognize that interatomic potential in our simulation is also a potential source of the “physical error” of the model. We, however, selected this model due to the comprehensive MD simulation study performed by Luo et al.18,19 with the same sets of interatomic potential.

Using Eq. (5) and the thermal conductance results presented in Fig. 6, we can quantitatively estimate the thermal conductivity of the Au/MNL nanolaminate. We calculated the average interfacial distance, d(XIBS/XIBD/XIMC) in the Au/MNL nanolaminate at 300 K for varying XIBS/XIBD/XIMC. The observed ranges as a function of bonding strengths were 12.76 Å d I B S 12.94 Å, 12.76 Å d I B D 14.29 Å, and 12.11 Å d I M C 12.76 Å. In the case of XIBS and XIBD, the average interfacial separation, d, was decreasing with increasing interfacial bonding strengths, while in the case of XIMC, removal of MNL chain molecules leads to the decrease of d. Figure 7 compares nanolaminate thermal conductivity obtained by direct MD simulation, κ, and that predicted by interfaces in the series model, κint. The first observation from Fig. 7 is that the simple resistors in the series model predict thermal conductivity quite well. This is quite surprising, considering the very close proximity of interfaces (∼1 nm). Detailed examination of the results presented in Fig. 7 reveals that thermal conductivities estimated from the interfacial model underpredict the actual values by about 10%–15% (see Fig. 7). This result indicates that when interfaces are close to each other, as in the case of nanolaminates, they have slightly lower effective interfacial resistance and are more conducting than individual interfaces. However, considering the small differences, we can consider thermal conductivity, κ in the Au/MNL nanolaminate well described by the resistances in the series model. This is in contrast to the behavior observed in semiconductor superlattices, where there are strong effects of coherent phonon propagation and scattering across multiple interfaces leading to an effect such as minimum thermal conductivity observed at a specific superlattice period (i.e., interfacial distance).10 

FIG. 7.

Thermal conductivity, κ, and interface predicted thermal conductivity, κint = G⋅d, of Au/MNL nanolaminate for different interfacial Au–S (a) bond strengths (XIBS), (b) bond density (XIBD), and (c) molecular coverage (XIMC).

FIG. 7.

Thermal conductivity, κ, and interface predicted thermal conductivity, κint = G⋅d, of Au/MNL nanolaminate for different interfacial Au–S (a) bond strengths (XIBS), (b) bond density (XIBD), and (c) molecular coverage (XIMC).

Close modal

All the results presented so far correspond to simulations carried out at 300 K. We have also explored the role of temperature on thermal conductivity and interfacial thermal conductance. Figure 8(a) shows that thermal conductivity decreases monotonically, albeit not significantly until 350 K, with increasing temperature, except for the weakest interface, which exhibits essentially no temperature effects. Majumdar et al.43 have also shown temperature independent behavior for imperfect Au/MNL/Au junction with G ∼ 100 MW m−2 K−1, which is qualitatively similar to our weakest interface results. A decrease in thermal conductivity in Au/MNL nanolaminates with increasing temperature as seen in the strongly bonded interfaces is typical behavior seen in the bulk crystalline materials and superlattices.44 This is attributed to increased phonon scattering with increasing temperature. For the weakest interface, interfacial scattering rather than phonon–phonon scattering likely dominates, thus the effect of temperature is not discernible.

FIG. 8.

Temperature dependence of thermal conductivity and interfacial thermal conductance under different interfacial Au–S bond strengths, XIBS. (a) Thermal conductivity of Au/MNL nanolaminates. (b) Interfacial thermal conductance of Au/MNL/Au junction. (c) Dihedral angle distribution of chains in Au/MNL/Au junction at XIBS= 1. (d) Snapshots of MNL chains at different temperatures at XIBS = 1.

FIG. 8.

Temperature dependence of thermal conductivity and interfacial thermal conductance under different interfacial Au–S bond strengths, XIBS. (a) Thermal conductivity of Au/MNL nanolaminates. (b) Interfacial thermal conductance of Au/MNL/Au junction. (c) Dihedral angle distribution of chains in Au/MNL/Au junction at XIBS= 1. (d) Snapshots of MNL chains at different temperatures at XIBS = 1.

Close modal

In contrast, the interfacial thermal conductance exhibits a moderate increase with the increasing temperature until 350 K, which agrees well with the literature.18,19 It is understood that interfacial thermal conductance increase with increasing temperature is associated with anharmonic phonon–phonon interactions that open new channels for heat flow across the interface, effectively lowering the junction resistance.45 Furthermore, we have provided numerical evidence of anharmonic interaction at the Au/MNL interface by performing the phonon wave packet simulations discussed later in Sec. III E 2.

While in Sec. III B, we argued that the thermal conductivity of the nanolaminates can be largely understood by the interfacial resistors in the series model, up to 350 K, the observed trend in thermal conductivity is opposite to that observed for interfacial thermal conductance. This difference is up to ∼50% at low temperature and bonding (discussed in Sec. III D) that is not accounted for by the resistors in the series model (see Fig. 7). This suggests that in addition to the serial independent thermal interfacial resistors, an explanation of complete thermal conductivity might require the employment of phonons propagating coherently across multiple interfaces.

The observed temperature dependence can be more formally understood within the recently proposed “Unified theory of thermal transport in crystals and glasses.”46 This theory expresses the total thermal conductivity as a sum of the two contributions: (i) propagating phonons that scatter due to thermal or structural fluctuations—this mechanism is dominant in crystalline materials and (ii) energy exchange between non-propagating vibrational modes—this contribution is significant in disordered and amorphous solids. When non-propagating vibrational energy dominates thermal transport, with increasing temperature thermal conductivity increases. On the other side, thermal scattering of propagating phonons decreases thermal conductivity with increasing temperature. Our observation of moderate thermal conductivity decrease with increasing temperature suggests that the propagating phonons temperature effect dominates the temperature dependence. Furthermore, the results of our phonon scattering simulations presented in Sec. III E suggest that indeed the vast majority of phonons propagate and scatter at the interfaces while preserving a significant fraction of their coherence.

Above 350 K, both thermal conductivity and interfacial thermal conductance follow a decreasing trend with increasing temperature. Prior studies on the Au/MNL/Au junction have not explored temperature effects beyond this point.18,19 We tracked down the origin of such a significant decrease in thermal conductivity and interfacial thermal conductance to the structural changes that occur in the MNL chains at high temperatures. As we increased the temperature beyond 350 K, we observed the increased density of kinks in the MNL chain configuration. The trans-straight chain conformations of MNL are characterized by a dihedral angle φ = 180 ° and gauche-rotational state is achieved when φ = 60 °. At a low temperature (T = 100 K), ∼95% of the MNL chains connecting the inorganic Au blocks are straight and have trans configurations as seen from the dihedral angle distribution in Fig. 8(c) and snapshots of MNL chains in Fig. 8(d). But, as we increase the temperature up to 350 K, dihedral angle distribution becomes broader and develops a second maximum at φ 70 ° [see Figs. 8(c) and 8(d)].

The fact that the resistance of the MNL phase increases with the introduction of kinks is consistent with numerous studies. Molecular simulations of Sasikumar et al.47 reported that heat flows ballistically in straight chains, while gauche conformations lead to an increase in chain resistance. Also, using MD simulations Duan et al.48 showed that the effective thermal conductivity of polyethylene chains decreases with an increasing number of kinks. Das et al.49 performed phonon wave packet simulations on folded polyethylene and showed interfacial thermal conductance decaying significantly with an increasing number of folds. A more intriguing question is why the presence of kinks affects the interfacial thermal conductance of Au/MNL interface. We will address this issue in Sec. III E 4 (see Fig. 13).

In Sec. III B, we discussed the applicability of the interfacial resistors in a series model regarding interfacial bonding at temperature 300 K. We showed that for XIBS/XIBD/XIMC, the resistors in the series model show thermal conductivity prediction on Au/MNL nanolaminate within the discrepancy of 10%–15% throughout the range of interfacial bonding. The temperature dependence of thermal conductivity, κ, interfacial thermal conductance, G, and interfacial separation, d, will allow us to further examine the applicability of the resistors in a series model, i.e., κint = G⋅d.

To this end, in Fig. 9(a), we show the variation of average spacing between interfaces, d with increasing temperature. It can be seen that d increases sharply with increasing temperature with a corresponding coefficient of thermal expansion (CTE) of ∼2 × 10−2 K−1. At first glance, this is very surprising considering the thermal expansion of gold is an order of magnitude lower, and the end-to-end distance of the MNL chains is expected to decrease with increasing temperature due to chain vibration and the development of kinks. In order to understand such a large thermal expansion of Au/MNL nanolaminate, we noticed a geometric change in MNL chain orientation in terms of MNL chains tilt angle. The average tilt angle of MNL chains decreases from 37.4° at 100 K to 16.3° at 500 K with a wider distribution [see Fig. 9(b)]. Such a large change in the tilt angle of MNL chains explains the large thermal expansion observed for Au/MNL nanolaminate.

FIG. 9.

Temperature and bonding effect on thermal conductivity in Au/MNL nanolaminate (a) Average interfacial distance, d, in Au/MNL nanolaminate. (b) Tilt angle distribution of MNLs in nanolaminate. (c) Thermal conductivity, κbulk, and interfacial resistors in the series model predicted thermal conductivity, κint = G⋅d for varying temperature and interfacial bonding strength. (d) The percentage difference of resistors in the series model predicted thermal conductivity from the actual one.

FIG. 9.

Temperature and bonding effect on thermal conductivity in Au/MNL nanolaminate (a) Average interfacial distance, d, in Au/MNL nanolaminate. (b) Tilt angle distribution of MNLs in nanolaminate. (c) Thermal conductivity, κbulk, and interfacial resistors in the series model predicted thermal conductivity, κint = G⋅d for varying temperature and interfacial bonding strength. (d) The percentage difference of resistors in the series model predicted thermal conductivity from the actual one.

Close modal

Both d(T) [Fig. 9(a)] and G(T) [Fig. 8(b)] increase with temperature up to 350 K, and consequently, the predicted thermal conductivity κint(T) = G(T)⋅d(T) of Au/MNL also increases with temperature—this is in contrast to observed thermal conductivity that decreases with temperature [see Fig. 9(c)]. Thus, the resistors in the series model do not explain temperature dependence. We attribute this behavior to the effects associated with coherent phonon transport across multiple interfaces and will discuss this further in Sec. III E 3. Additional credence to our speculations is provided in Figs. 9(c) and 9(d) showing that the resistance in the series model is most accurate at higher temperatures and bonding (∼6% error at 350 K) and has up to ∼50% discrepancy at the lowest temperature and bonding studied. Beyond 350 K, the resistor in the series model has to be modified to account for the non-negligible resistances of the organic chains due to the presence of kinks.

In summary, the resistors in the series model is more accurate for interfaces with strong bonding and at higher temperatures, while it is less predictive for weakly bonded interfaces and at lower temperatures, which we attribute to the effect of the coherent phonon transport across multiple interfaces.

All the results presented so far are the observations from direct MD simulations on Au/MNL nanolaminates and Au/MNL/Au junctions. We have observed simple interfacial resistors in the series model describing the majority of thermal transport contributions in Au/MNL nanolaminates with increasing temperature and bonding. In Secs. III A and III D, we speculated the origin of such difference being phonons propagating coherently across Au/MNL multiple interfaces. Similarly, we also attributed the increasing interfacial thermal conductance in the Au/MNL/Au junction (see Sec. III C) to anharmonic phonon–phonon interaction at the interface. In this section, we implement phonon wave packet simulations in the Au/MNL single and multiple interfaces at the individual phonon level to provide quantitative evidence of coherent phonon transport as well as anharmonicity at the individual Au/MNL interface. For this purpose, we have computed the frequency-dependent transmission coefficients, performed spectral analysis, and provided the WP-based insights into our direct MD simulation results presented earlier in Secs. III AIII D.

1. Transmission coefficient

Figure 10(a) shows the energy-transmission coefficient for several LA wave packets across the Au/MNL single interface. The average transmission coefficient α is calculated after the transmitted energy post-scattering event. We sampled a total of 80 frequency points in the LA mode. Figure 10(b) shows the frequency-dependent transmission coefficient for LA phonon wave packets, which exhibit a strong phonon frequency dependence with nonmonotonic behavior. The observed frequency-dependent transmission coefficient for the Au/MNL single interface is qualitatively similar to the Au/MNL/liquid system.31 The maximum transmission of α = 0.66 occurs at 4.29 THz and minimum α = ∼ 0 at 4.76 THz of Au/MNL single interface. Wei et al.31 showed maximum transmission at 4.4 THz for Au/MNL/liquid system using LA phonon mode and a larger frequency sampling interval of ∼1 THz. We also used the acoustic mismatch model (AMM)8 to predict the transmission coefficient of α = 0.34 for the Au/MNL model system using the relation α = 4ZAuZPE/(ZAu + ZPE)2, where ZAu is the acoustic impedance of Au and ZPE is the acoustic impedance of tilted polyethylene. AMM prediction simply uses acoustic impedance mismatch for predicting transmission coefficients without considering the exact nature of interfacial bonding. AMM prediction of α = 0.34 is close to the α = 0.27 observed directly from the WP simulation near the Γ point.

FIG. 10.

(a) Energy transmission across the single Au/MNL interface for LA phonon modes with varying frequency ν. (b) Frequency-dependent energy-transmission coefficient α(ν) from WP simulation in an Au/MNL single interface. (c) The energy transmission across the Au/MNL multiple interfaces for different phonon modes with varying frequency ν. (d) Frequency-dependent energy-transmission coefficient α(ν) from WP simulation in an Au/MNL multiple interfaces.

FIG. 10.

(a) Energy transmission across the single Au/MNL interface for LA phonon modes with varying frequency ν. (b) Frequency-dependent energy-transmission coefficient α(ν) from WP simulation in an Au/MNL single interface. (c) The energy transmission across the Au/MNL multiple interfaces for different phonon modes with varying frequency ν. (d) Frequency-dependent energy-transmission coefficient α(ν) from WP simulation in an Au/MNL multiple interfaces.

Close modal

Figures 10(c) and 10(d) show the LA wave packet scattering results from the Au/MNL multiple interfaces. Figure 10(d) shows the frequency-dependent transmission coefficient of Au/MNL multiple interfaces, which is dramatically different from the single interface. The wave packets reflected from multiple Au/S and S/Au interfaces can interfere constructively or destructively leading to multiple minima and maxima across the whole frequency range. This is indicative of the strong presence of coherent phonon interference effects in Au/MNL multiple interfaces and rather contradictory to the underpinnings of the resistor in the series model (see Secs. III B and III D ).

2. Spectral analysis

We analyzed the WP scattering mechanism in detail for the single and multiple interfaces by monitoring the reflected and transmitted phonon spectra obtained via power spectrum analysis described in Sec. II C 3. Figure 11 shows the analysis of phonon wave packet scattering across the single Au/MNL interface. Upon scattering from the interface, phonons scatter into the same frequencies upon reflection and transmission, indicated by a continuous array of points on a solid line (see Fig. 11 top and bottom panel). This can be attributed to the continuous presence of phonon modes in the MNL region across the whole range of phonon frequencies due to the overlapping of vibrational density of states (VDOSs) of Au and MNL. The presence of reflected and transmitted transverse acoustic components [see Figs. 11(a), 11(b), 11(d), and 11(e)] indicates that after scattering the incident longitudinal acoustic wave also transforms into transverse acoustic waves. In addition, at selected frequency locations wave packets also scatter into the higher frequency phonons upon reflection/transmission, resulting in phonon frequencies 2νincident (see data points in Fig. 11 dotted line). We have classified the low-intensity and high-intensity peaks in the power spectrum as (1–4)% of F(ν)max and (5–100)% of F(ν)max, respectively. In both cases of reflection and transmission, all the low-intensity peak contributions (+) of the power spectrum occur only at higher frequencies (2νincident) and all the high-intensity peak contributions (●) preserve the incident wave packet frequency νincident. This indicates strong harmonic phonon scattering at interfaces but also the presence, albeit weaker, anharmonic scattering.

FIG. 11.

Phonon wave packet scattering from Au/MNL single interface with varying LA mode frequencies in the Г–L phonon branch propagating in the z-direction. Top panels (a)–(c) are the reflected wave packet frequencies after phonon scattering from a single interface. Bottom panels (d)–(f) are the transmitted wave packet frequencies after phonon scattering from a single interface. The solid line indicates a slope of one while the dotted line indicates a slope of two. Data points shown with “+” markers are low-intensity peak contributions in the power spectrum representing (1–4)% of F(ν)max, and the rest of the data points shown with “●” markers are defined as high-intensity peak contributions.

FIG. 11.

Phonon wave packet scattering from Au/MNL single interface with varying LA mode frequencies in the Г–L phonon branch propagating in the z-direction. Top panels (a)–(c) are the reflected wave packet frequencies after phonon scattering from a single interface. Bottom panels (d)–(f) are the transmitted wave packet frequencies after phonon scattering from a single interface. The solid line indicates a slope of one while the dotted line indicates a slope of two. Data points shown with “+” markers are low-intensity peak contributions in the power spectrum representing (1–4)% of F(ν)max, and the rest of the data points shown with “●” markers are defined as high-intensity peak contributions.

Close modal

Figure 12 shows the analysis of phonon wave packet scattering across the multiple Au/MNL interface. Multiple interfaces have more complex scattering processes leading to a larger spectrum of scattered phonons with enhanced anharmonic phonon scattering. The complex wave packet scattering in the multiple interface results in multiple frequency phonons on reflection and transmission in contrast to single interface scattering characterized with scattering to the same frequency (νincident) and double frequency (2νincident). The reflection/transmission frequencies of the wave packet have the same frequency, νincident (see solid line Fig. 12), double the frequency, 2νincident (see dotted line Fig. 12), but also more than double and well as lower than the frequency of the incident wave packet.

FIG. 12.

Phonon wave packet scattering from Au/MNL multiple interfaces with varying LA mode frequencies in the Г–L phonon branch. Top panels (a)–(c) are the reflected wave packet frequencies after phonon scattering from a multiple interface. Bottom panels (d)–(f) are the transmitted wave packet frequencies after phonon scattering from a multiple interface. The solid line indicates a slope of one while the dotted line indicates a slope of two. Data points shown with “+” and “●” markers are the same as defined in Fig. 11.

FIG. 12.

Phonon wave packet scattering from Au/MNL multiple interfaces with varying LA mode frequencies in the Г–L phonon branch. Top panels (a)–(c) are the reflected wave packet frequencies after phonon scattering from a multiple interface. Bottom panels (d)–(f) are the transmitted wave packet frequencies after phonon scattering from a multiple interface. The solid line indicates a slope of one while the dotted line indicates a slope of two. Data points shown with “+” and “●” markers are the same as defined in Fig. 11.

Close modal
FIG. 13.

Phonon wave packet simulation for 3 THz LA phonon performed on Au/MNL single interface with and without the kink interface present in the MNL structure. (a) Snapshots of the spatial distribution of velocity for MNL without kink. (b) Snapshots of the spatial distribution of velocity for MNL with the kink interface located 100 nm away from the Au/MNL interface. (c) Energy transmission through the Au/MNL single interface with and without the kink interface.

FIG. 13.

Phonon wave packet simulation for 3 THz LA phonon performed on Au/MNL single interface with and without the kink interface present in the MNL structure. (a) Snapshots of the spatial distribution of velocity for MNL without kink. (b) Snapshots of the spatial distribution of velocity for MNL with the kink interface located 100 nm away from the Au/MNL interface. (c) Energy transmission through the Au/MNL single interface with and without the kink interface.

Close modal

3. Resistors in series model vs observation from WP simulations

It is clear from Fig. 10 that the phonon scattering mechanism is complex and interference effects are apparent at Au/MNL multiple interfaces. However, earlier [see Figs. 7 and 9(c)] we observed a simple thermal independent interfacial resistor in a series model that explains the majority of the thermal transport in Au/MNL nanolaminates more accurately at higher temperatures. This would suggest the lack of interference/coherent scattering effects.

A similar apparent inconsistency was observed by Shen et al.50 for the case of the Si/Graphene/Si interface, while the overall interfacial thermal conductance of a Si/Graphene/Si interface was independent of the number of graphene layers, the transmission coefficient of individual phonons was strongly dependent on the number of graphene layers. However, Shen et al.50 demonstrated that the interfacial thermal resistance estimated from an integral over all the incident phonons is virtually independent of the number of graphene layers. In other words, the actual scattering process shows coherent multi-interfacial scattering, which is effectively integrated out at the level of interfacial thermal transport coefficient.

To resolve a similar inconsistency observed in our simulations of organic/inorganic interfaces, we recall that the interfacial thermal conductance can be determined from a sum over all phonon modes by utilizing the relation35,51
G = 1 ( 2 π ) 2 k γ { L A , T A x , T A y } W k ω γ ( k ) v γ ( k ) α γ ( k ) n ( ω , T ) T
(6)
where G is interfacial thermal conductance, ℏ is the reduced Plank constant, ω γ ( k ) = phonon frequency, n(ω, T) is the Bose–Einstein distribution function at temperature T, Wk is the weight factor associated with k-space volume, v γ ( k ) = phonon group velocity. To compare this formula with results with classical MD simulation we use Boltzmann statistics n (ω, T) = kBT/ℏω, where kB is the Boltzmann constant. Integrating along the kx and ky direction in Eq. (6) within the classical limit reduces to31 
G = 4 k B a 2 γ { L A , T A x , T A y } v γ ( k z ) α γ ( k z ) Δ k z
(7)

Using the grid mesh Δ k z, the Г–L phonon band and taking only the γ { L A } phonon mode in Eq. (7) we have estimated the ratio G m u l t i p l e G s i n g l e = 0.345, where Gsingle is the interfacial thermal conductance of single interface. From the simple independent resistors in the series model, one gets G m u l t i p l e G s i n g l e = 0.25, i.e., the resistance of 4 interfaces is 4 times larger than the resistance of a single interface. Considering the approximate nature of the estimate provided by Eq. (7) using only LA modes, the agreement between the estimate given by Eq. (7) to that of the independent resistors in the series model is quite good. As was the for the Si/Graphene/Si interface50 discussed above, while the transmission coefficient behavior is drastically affected by the number of interfaces, the integrated effect is quantitively very similar to that predicted by the interfaces in the series model. Although the pattern of transmission coefficient is very different in single and multiple interfaces, an integral over all the phonon modes predicts that the overall resistance of multi-interfacial structure is more or less proportional to the number of interfaces present which explains the applicability of the resistors in series model here.

4. Temperature effect—insights from WP simulations

The coherent phonon propagation across the multiple interfaces observed from the WP simulation of Au/MNL [see Fig. 10(d)] suggests a reason for the thermal conductivity decrease with increasing temperature observed for Au/MNL nanolaminates [see Fig. 8(a)]. As shown in Fig. 12, phonon propagation across multiple interfaces is characterized by a significant role of anharmonic scattering processes. Such processes scatter phonons and thus limit their mean free path and in consequence, reduce thermal conductivity. This scattering is expected to increase with increasing temperature due to higher vibrational amplitudes.

According to Fig. 11, anharmonic scattering is also present at individual interfaces. This suggests a reason for increasing interfacial thermal conductance with increasing temperature [see Fig. 8(b)]. The reasoning here is based on the fact that anharmonic scattering channels increase overall transmission coefficient and according to Eq. (6), increase overall interfacial thermal conductance. Anharmonic scattering strength increases with increasing phonon vibrational amplitudes caused by increasing temperature thus explaining the increase of the interfacial thermal conductance as observed in Fig. 8(b).

Above 350 K, both thermal conductivity and interfacial thermal conductance decrease with increasing temperature. We already discussed the role of kinks on the thermal conductivity of the MNL phase (see Sec. III C), which contributes to the observed thermal conductivity decrease. A separate question is why interfacial thermal conductance decreases with increasing temperature above 350 K. To address this question, we study the effect of the kink interface in the MNL phase on the Au/MNL interfacial thermal conductance. We performed phonon wave packet simulation on a model structure with a plane of kinks at a distance of 100 nm away from the Au/MNL interface (see Fig. 13). According to the data presented in Fig. 13, the presence of the kink interface leads to a significant reduction (∼22%) of the phonon transmission coefficient at the Au/MNL interface. This is not a surprising result, considering that the phonon-mean free path along straight polyethylene chains is very large. This leads to the interdependence of the Au/MNL interface and kink interface. Such phonon-based interdependence of proximal interfaces was observed in simulations of nanoscopic adlayers,52 and the related effects were observed in experiments on thermal transport through point contacts between graphitic materials.53 

Using a combination of NEMD and phonon wave packet dynamics simulations, we have elucidated the thermal transport mechanism in an Au/MNL. We demonstrated that the thermal conductivity of organic/inorganic Au/MNL nanolaminates with an interfacial separation of ∼1 nm, can be tuned from low (∼0.7 W m−1 K−1) to ultralow (∼0.07 W m−1 K−1) values via various interfacial bonding manipulation strategies. We showed that the heterogeneous reduction of the bonding strength created by reduced bond density or reduced molecular coverage has a greater effect on thermal conductivity compared to a homogeneous reduction of bonding strength. We observed that the bulk thermal conductivity of the Au/MNL nanolaminate can be reasonably explained by the independent interfacial resistors connected in the series model more accurately for interfaces with strong bonding and at higher temperatures, while it is less predictive for weakly bonded interfaces and at lower temperatures. We attribute these effects to the coherent phonon transport across multiple interfaces. Such coherent transport is expected to be more pronounced at low temperatures, which explains large differences between model predicted and observed thermal conductivity values. At higher temperatures, coherence effects are expected to be diminished leading to more accurate predictions of resistors in series models where phonons are assumed to scatter independently at each interface.

Our wave packet simulations of phonon transmission across individual and multiple Au/MNL interfaces indicate that coherent phonon interference effects are present across the whole phonon frequency range. This is rather surprising considering that interfacial resistors in the series model provide a reasonable description of thermal transport in nanolaminates. We explain this contradiction by noticing that thermal conductivity as well the interfacial thermal conductance represents quantities that are integrated/averaged over all the phonon modes. Quite interestingly to us, it appears that while coherency is pronounced at the individual phonon level, its effect on integrated quantity, namely, interfacial thermal conductance is diminished due to averaging.

This work was funded by the Focus Interconnect Center at Rensselaer Polytechnic Institute (RPI) and supported by the New York State Department of Economic Development, Grant No. C180117. We are grateful to RPI's Center for Computational Innovations (CCI) for providing the computational time required for this work.

The authors have no conflicts to disclose.

Rajan Khadka: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Pawel Keblinski: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.

1.
J.
Azadmanjiri
,
C. C.
Berndt
,
J.
Wang
,
A.
Kapoor
, and
V. K.
Srivastava
,
RSC Adv.
6
,
109361
109385
(
2016
).
2.
B.
Yeom
,
A.
Jeong
,
J.
Lee
, and
K.
Char
,
Polymer
91
,
187
193
(
2016
).
3.
K. H.
Yoon
,
H. S.
Kim
,
K. S.
Han
,
S. H.
Kim
,
Y. E. K.
Lee
,
N. K.
Shrestha
,
S. Y.
Song
, and
M. M.
Sung
,
ACS Appl. Mater. Interfaces
9
,
5399
5408
(
2017
).
4.
T.
Poirié
,
T.
Schmitt
,
E.
Bousser
,
R.
Vernhes
,
L.
Martinu
, and
J. E.
Klemberg-Sapieha
,
Surf. Coatings Technol.
315
,
399
407
(
2017
).
5.
R.
Khadka
,
G.
Ramanath
, and
P.
Keblinski
,
Sci. Rep.
12
,
10788
(
2022
).
6.
S.
Park
,
J.
Jang
,
H.
Kim
,
D.
Il Park
,
K.
Kim
, and
H. J.
Yoon
,
J. Mater. Chem. A
8
,
19746
19767
(
2020
).
7.
K.
Wang
,
E.
Meyhofer
, and
P.
Reddy
,
Adv. Funct. Mater.
30
,
1904534
(
2020
).
8.
J.
Chen
,
X.
Xu
,
J.
Zhou
, and
B.
Li
,
Rev. Mod. Phys.
94
,
025002
(
2022
).
9.
A.
Giri
and
P. E.
Hopkins
,
Adv. Funct. Mater.
30
,
1903857
(
2020
).
10.
M.
Nomura
,
R.
Anufriev
,
Z.
Zhang
,
J.
Maire
,
Y.
Guo
,
R.
Yanagisawa
, and
S.
Volz
,
Mater. Today Phys.
22
,
100613
(
2022
).
11.
M. D.
Losego
,
I. P.
Blitz
,
R. A.
Vaia
,
D. G.
Cahill
, and
P. V.
Braun
,
Nano Lett.
13
,
2215
2219
(
2013
).
12.
C.
Wan
,
X.
Gu
,
F.
Dang
,
T.
Itoh
,
Y.
Wang
,
H.
Sasaki
,
M.
Kondo
,
K.
Koga
,
K.
Yabuki
,
G. J.
Snyder
,
R.
Yang
, and
K.
Koumoto
,
Nat. Mater.
14
,
622
627
(
2015
).
13.
A.
Giri
,
J.-P.
Niemelä
,
T.
Tynell
,
J. T.
Gaskins
,
B. F.
Donovan
,
M.
Karppinen
, and
P. E.
Hopkins
,
Phys. Rev. B
93
,
115310
(
2016
).
14.
A.
Giri
,
J.-P.
Niemelä
,
C. J.
Szwejkowski
,
M.
Karppinen
, and
P. E.
Hopkins
,
Phys. Rev. B
93
,
024201
(
2016
).
15.
E.
Lee
and
T.
Luo
,
Appl. Phys. Lett
112
,
011603
(
2018
).
16.
T.
Murakami
,
T.
Hori
,
T.
Shiga
, and
J.
Shiomi
,
Appl. Phys. Express
7
,
121801
(
2014
).
17.
K.
Gordiz
and
A.
Henry
,
Appl. Phys. Lett.
108
,
181606
(
2016
).
18.
T.
Luo
and
J. R.
Lloyd
,
Int. J. Heat Mass Transfer
53
,
1
11
(
2010
).
19.
T.
Luo
and
J. R.
Lloyd
,
J. Heat Transfer
132
,
032401
(
2010
).
20.
M.
Hu
,
P.
Keblinski
, and
P. K.
Schelling
,
Phys. Rev. B
79
,
104305
(
2009
).
21.
Y.
Wang
,
Y.
Cao
,
K.
Zhou
, and
Z.
Xu
,
Adv. Mater. Interfaces
4
,
1700355
(
2017
).
22.
X. K.
Chen
,
Z. X.
Xie
,
W. X.
Zhou
,
L. M.
Tang
, and
K. Q.
Chen
,
Appl. Phys. Lett.
109
,
023101
(
2016
).
23.
R.
Cheaito
,
C. A.
Polanco
,
S.
Addamane
,
J.
Zhang
,
A. W.
Ghosh
,
G.
Balakrishnan
, and
P. E.
Hopkins
,
Phys. Rev. B
97
,
085306
(
2018
).
24.
M.
Hu
and
D.
Poulikakos
,
Nano Lett.
12
,
5487
5494
(
2012
).
25.
Y.
Yin
,
K.
Baskaran
, and
A.
Tiwari
,
Phys. Status Solidi Appl. Mater. Sci.
216
,
1800904
(
2019
).
26.
A. I.
Jewett
,
D.
Stelter
,
J.
Lambert
,
S. M.
Saladi
,
O. M.
Roscioni
,
M.
Ricci
,
L.
Autin
,
M.
Maritan
,
S. M.
Bashusqeh
,
T.
Keyes
,
R. T.
Dame
,
J.-E.
Shea
,
G. J.
Jensen
, and
D. S.
Goodsell
,
J. Mol. Biol.
433
,
166841
(
2021
).
27.
X. P.
Liu
,
Y.
Ni
, and
L. H.
He
,
Nanotechnology
27
,
1135707
(
2016
).
28.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in ‘t Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
, and
S. J.
Plimpton
,
Comput. Phys. Commun.
271
,
108171
(
2022
).
29.
P. K.
Schelling
,
S. R.
Phillpot
, and
P.
Keblinski
,
Appl. Phys. Lett.
80
,
2484
(
2002
).
30.
P. K.
Schelling
and
S. R.
Phillpot
,
J. Appl. Phys.
93
,
5377
(
2003
).
31.
X.
Wei
and
T.
Luo
,
J. Appl. Phys.
126
,
015301
(
2019
).
32.
X.
Liu
,
Q.
Wang
,
R.
Wang
,
S.
Wang
, and
X.
Liu
,
J. Appl. Phys.
133
,
095101
(
2023
).
33.
J. D.
Gale
and
A. L.
Rohl
,
Mol. Simul.
29
,
291
341
(
2011
).
34.
J. W.
Lynn
,
H. G.
Smith
, and
R. M.
Nicklow
,
Phys. Rev. B
8
,
3493
(
1973
).
35.
C.
Kimmer
,
S.
Aubry
,
A.
Skye
, and
P. K.
Schelling
,
Phys. Rev. B
75
,
144105
(
2007
).
36.
J. W.
Cooley
and
J. W.
Tukey
,
Math. Comput.
19
,
297
(
1965
).
37.
C.R.
Harris
,
K.J.
Millman
,
S.J.
van der Walt
,
R.
Gommers
,
P.
Virtanen
,
D.
Cournapeau
,
E.
Wieser
,
J.
Taylor
,
S.
Berg
,
N.J.
Smith
,
R.
Kern
,
M.
Picus
,
S.
Hoyer
,
M.H.
van Kerkwijk
,
M.
Brett
,
A.
Haldane
,
J.F.
del Río
,
M.
Wiebe
,
P.
Peterson
,
P.
Gérard-Marchant
,
K.
Sheppard
,
T.
Reddy
,
W.
Weckesser
,
H.
Abbasi
,
C.
Gohlke
, and
T.E.
Oliphant
,
Nature
585
,
357
362
(
2020
).
38.
M. D.
Losego
,
M. E.
Grady
,
N. R.
Sottos
,
D. G.
Cahill
, and
P. V.
Braun
,
Nat. Mater.
11
,
502
506
(
2012
).
39.
S.
Majumdar
,
J. A.
Malen
, and
A. J. H.
McGaughey
,
Nano Lett.
17
,
220
227
(
2017
).
40.
C.
Uher
,
Thermal Conductivity: Theory, Properties, and Applications
(Kluwer Academic / Plenum Publishers, New York,
2004
), pp.
21
91
.
41.
Z.
Wang
,
J. A.
Carter
,
A.
Lagutchev
,
K. K.
Yee
,
N. H.
Seong
,
D. G.
Cahilll
, and
D. D.
Dlott
,
Science
317
,
787
790
(
2007
).
42.
P. J.
O’Brien
,
S.
Shenogin
,
J.
Liu
,
P. K.
Chow
,
D.
Laurencin
,
P. H.
Mutin
,
M.
Yamaguchi
,
P.
Keblinski
, and
G.
Ramanath
,
Nat. Mater.
12
,
118
122
(
2012
).
43.
S.
Majumdar
,
J. A.
Sierra-Suarez
,
S. N.
Schiffres
,
W. L.
Ong
,
C. F.
Higgs
,
A. J. H.
McGaughey
, and
J. A.
Malen
,
Nano Lett.
15
,
2985
2991
(
2015
).
44.
J.
Garg
,
N.
Bonini
, and
N.
Marzari
,
Nano Lett.
11
,
5135
5141
(
2011
).
45.
P. E.
Hopkins
,
J. C.
Duda
, and
P. M.
Norris
,
J. Heat Transfer
133
,
062401
(
2011
).
46.
M.
Simoncelli
,
N.
Marzari
, and
F.
Mauri
,
Nat. Phys.
15
,
809
813
(
2019
).
47.
K.
Sasikumar
and
P.
Keblinski
,
J. Appl. Phys
109
,
114307
(
2011
).
48.
X.
Duan
,
Z.
Li
,
J.
Liu
,
G.
Chen
, and
X.
Li
,
J. Appl. Phys.
125
,
164303
(
2019
).
49.
S.
Das
and
M.
Muthukumar
,
Macromolecules
56
,
393
403
(
2022
).
50.
M.
Shen
,
P. K.
Schelling
, and
P.
Keblinski
,
Phys. Rev. B
88
,
045444
(
2013
).
51.
B.
Deng
,
A.
Chernatynskiy
,
M.
Khafizov
,
D. H.
Hurley
, and
S. R.
Phillpot
,
J. Appl. Phys.
115
,
084910
(
2014
).
52.
Z.
Liang
,
K.
Sasikumar
, and
P.
Keblinski
,
Phys. Rev. Lett.
113
,
065901
(
2014
).
53.
J.
Yang
,
M.
Shen
,
Y.
Yang
,
W. J.
Evans
,
Z.
Wei
,
W.
Chen
,
A. A.
Zinn
,
Y.
Chen
,
R.
Prasher
,
T. T.
Xu
,
P.
Keblinski
, and
D.
Li
,
Phys. Rev. Lett.
112
,
205901
(
2014
).