We theoretically explore the influence of end-group chemistry (bond stiffness and mass) on the interfacial thermal conductance at a gold–alkane interface. We accomplish this using the nonequilibrium Green's function (NEGF) coupled with first principle parameters in density functional theory (DFT) within the harmonic approximation. Our results indicate that the interfacial thermal conductance is not a monotonic function of either chemical parameters but instead maximizes at an optimal set of mass and bonding strength. This maximum is a result of the interplay between the overlap in local density of states (LDOS) of the device and that in the contacts, as well as the phonon group velocity. We also demonstrate the intrinsic relationship between the diffusive mismatch model (DMM) and the properties from NEGF, and provide an approach to get DMM from first principles NEGF. By comparing the NEGF-based DMM conductance and range of conductance while altering the mass and bonding strength, we show that DMM provides an upper bound for elastic transport in this dimension-mismatched system. We thus have a prescription to enhance the thermal conductance of systems at low temperatures or at low dimensions where inelastic scattering is considerably suppressed.
Hybrid junctions between inorganic components and organic molecules have been of great interest due to their potential applications in electronic devices [1–3], thermoelectrics [4,5], solar cells [6,7], and optoelectronics . In recent years, the hybrid junctions, particularly the self-assembled monolayers (SAMs) have been intensively explored for thermal properties because of the discovery that the interfacial thermal conductance can be enhanced by more than a factor of four at room temperature, simply by embedding the commonly assumed thermally insulating polymers (alkane chains) at the metal/dielectric interface [9,10]. This enhancement can even be tuned over one order of magnitude by carefully varying the interfacial properties (bond stiffness and mass) with end-group chemistry .
The observations of the influence of the interfacial properties on the interfacial thermal conductance diverge in the existing literature. One set of articles [11,12] reports that the stronger the bonding strength, the larger the thermal conductance. This result is commonly predicted by molecular dynamics simulations of interfaces between solids. A contrary viewpoint [13–17] is arrived at by lattice dynamics and Green's function on one-dimensional (1D) harmonic chains, showing that a maximum conductance exists when the interface properties help match the impedance from two sides. This maximum conductance can be obtained with an optimal interfacial bonding (the harmonic mean of the bonding from two sides) and the favored connecting mass (the geometric mean of the masses from two sides). A Klemens  like scattering treatment would be consistent with this viewpoint, since the transmission depends quadratically on the mass deviation, a sweet spot for the mass value can be chosen to produce the maximum conductance.
The discrepancy can arise from two sources. The first is the presence of additional inelastic channels created by phonon–phonon interactions across interfaces. Note that the theoretical observations showing the stronger the bonding, the larger the conductance involves molecular dynamic simulations, where anharmonicity plays a significant role [11,12,19]. On the contrary, the maximum conductance arises in the absence of phonon–phonon interaction, relevant at low temperatures and low dimensions [13–16]. A separate origin for the discrepancy is the dimensionality difference in the two sets of simulations. An increase of conductance with increasing bonding prediction is found in three-dimensional (3D)–3D interfaces (e.g., thin films) where the system can be decoupled into separate 1D chains for each transverse momentum, which is fundamentally different from 1D–1D interfaces with optimal conductance prediction. One illustrative example of the significance of dimensionality is our recent work  on 3D–3D diamond system (Figs. 2 and 4 in Ref. ) which shows that a bridging junction with masses between the masses of the two materials at the interfaces can only decrease the ballistic interfacial thermal conductance, because of the need to conserve the transverse momentum. This mode matching issue described above, however, is not relevant in 1D with a single mode.
Recent research shows that anharmonicity has limited influence on the metal–SAM systems, which leads the phonons across the metal–alkane junctions elastically [21–25]. Phonon interference in the metal–SAM systems was predicted both theoretically [21,23] and experimentally , demonstrating that phonons across the chains encounter only weak phonon–phonon interaction. Recently, experiments from Majumdar  show a decrease of the interfacial thermal conductance of this system as the mismatch of the Debye temperature of the metal leads increases, which is a sign for elastic scattering as the conductance inversely follows the overlap of the phonon spectra. In this paper, we use the nonequilibrium Green's function (NEGF) method under the harmonic approximation to explore the influence of bond stiffness and end-group mass on the interfacial thermal conductance of the Au–alkane junctions. The interatomic force constants required for the calculations are obtained from first principles and the harmonic approximation is supported by the negligible role of phonon–phonon interaction in metal–alkane junctions [21–25]. We show that a maximum conductance can be obtained with a set of optimal interfacial properties, and the maximum conductance is related to the maximum overlap between the density of states (DOS) at the interface and the DOS of the materials at two sides, which serves like a frequency window. As the interfacial properties deviate from the optimal ones, the DOS overlap decreases causing a nonmonotonic change in conductance. Our result of maximum conductance should work for any 3D–1D systems as long as phonon transport is elastic. And it should be well applicable since, in the 3D–1D systems, elastic transport was proven by recent experiments [21–25].
The rest of the paper is organized as follows: first, we briefly discuss the structures of Au–alkane surfaces with both thiol and amine as end groups. Then, we describe the way the bond stiffness is varied in this study. After that, we introduce the first principle NEGF method with the “scooping” technique for the dimensionality transition followed by a short description of NEGF-based diffusive mismatch model (DMM) method. Finally, we discuss and analyze our first principle interfacial thermal conductance results and generalize them using a simple 3D–1D model.
Model and Simulation Method
Au/Alkane Surface Structure.
We model the Au–alkane surface using the unit cell [26,27]. This structure corresponds to the highest density possible of alkane chains assembled on gold surface, which also leads to the highly ordered chains . Figure 1 shows the top view of the absorption position of alkanes on gold surface and the stacked view of junctions. In Fig. 1(b), the black parallelograms denote the configurations, which is the simulation cell for the periodic boundary condition in the direction parallel to the gold surface. We use the hexagonal close-packed (hcp) position as the adsorption position in this simulation (Fig. 1(b)). The alkane chains before geometry optimization are set to tilt along the hcp (hcp hollow site)–fcc (face-centered cubic hollow site) dash lines in Fig. 1(b), with tilt angles of 10 deg, 20 deg, and 30 deg. This geometry is very close to the optimized ones reported in the existing literature (20∼30 deg between the tilted alkanes and the normal to the gold surface). During the first principle simulations, the alkane length is set to be 6 C atoms in the chain with the end group as thiol or amine. In addition to the fcc–hcp bridge site for adsorption of alkanes on the clean gold surface, we also considered a reconstructed gold surface [28–31]. In this case, a gold adatom is initially centered on top of the hcp position with the alkane end group connected to it. We consider these two structures because the bonding of a thiol-terminated gold interface is commonly acknowledged as covalent bonding with thiol on the clean gold surface , whereas the amine-terminated surface can bind when the surface is rough, and with Au adatoms . After these two kinds of structures are set up in density functional theory (DFT), the whole system is then relaxed during the optimization process. Both the clean surface and the adatom surface are simulated as shown in Fig. 1(c).
The geometry optimization and the force constants are calculated using the first principles package QUANTUM ESPRESSO . Ultrasoft pseudopotentials with local density approximation exchange-correlation function are used in the simulation. The cutoff energy for wave functions is calibrated for convergence with a minimum of 35 Ry, while the cutoff energy for charge density is set to be 300 Ry. The threshold for the convergence of ground state energy is set to be 1.0 × 10−10 Ry.
Before the geometry relaxation of the alkane chains assembled on the gold surface, we performed the structure optimization of bulk gold and infinite polyethylene. An 8 × 8 × 8 Monkhorst–Pack grid is used for k sampling in our bulk gold ground state simulation, with the Methfessel–Paxton smearing parameter set at 0.03 Ry. A 4 × 4 × 4 grid is used for q space sampling for gold force constants calculation using density functional perturbation theory (DFPT) . In the polyethylene simulation, a 1 × 1 × 12 grid is used in both the ground state calculation and DFPT in which the chain orientation is along the z-direction. The distance between chains is set to be fairly large (10 Å) to avoid interactions between them. The optimized lattice constant for gold is 4.03 Å and the optimized length for one C2H4 unit is 2.5215 Å. The dispersions using calculated force constants are shown in Figs. 2(a) and 2(b), along with the existing experimental data in the literature [35–38].
The distance between the radicals and the gold surface before relaxation is set to be 2 Å. A supercell contains a unit cell, four (111) Au layers and 20 Å vacuum is used as the repeated unit in DFT simulations to follow the periodic boundary conditions in three dimensions. The 20 Å vacuum is enough to avoid interactions between the adsorbed alkane chains and the neighboring gold slab. A 3 × 3 × 1 grid is used for both the structure optimization process and the force constant calculation process. The threshold for force convergence is set to be 1.13 × 10−3 Ry/Å.
After the calculation of the lattice dynamic matrix by DFPT, a Fourier transformation of the dynamic matrix is performed to obtain the real space force constant parameters. The q grid in DFPT calculations corresponds to the atomic interaction range. Up to second neighboring unit, atomic interactions are calculated for gold; up to seventh neighboring unit, cell interactions are calculated for polyethylene. Although the 3 × 3 × 1 grid provides the interaction between chains, in the following transport calculations, we ignore those interactions for simplicity.
Nonequilibrium Green's Function and Landauer Formula.
In the above equation, Md and Kd are the mass matrix and the force constant matrix for the device, respectively, while Σl and Σr are the self-energies for the contacts whose anti-Hermitian parts give us the broadening matrices Γl,r.
Diffusive Mismatch Model From Ab Initio Nonequilibrium Green's Function.
The DMM is commonly used to predict the interfacial thermal conductance. In order to set a reference conductance for Au–SAM junctions, we calculate the DMM conductance from the first principle parameters. Instead of using the frequency-dependent DOS and velocity from the full nonequilibrium lattice dynamics, we develop a quick estimate using DMM–NEGF, using variables calculated from Green's function with first principle parameters.
M1 and M2 are the number of modes of the bulk materials at each side of the interfaces. We can get M1 and M2 similarly in Green's function by defining the contacts and junction as the same materials, using . In the bulk materials, for each mode, transmission probability is 1, so is just the number of modes. Hence, the equals the number of modes of materials at two sides in parallel.
In our simulations, we use a compact Au–SAM density with the surface geometry as . Each alkane chain connects with 3 Au atoms in an area of 2.1 × 10−19 m2. The simulated number of modes for Au and alkane chain in that area are shown in Fig. 2(c), with the inset graph showing the . The number of modes for the Au with that specific density is three times as the number of modes of Au in primitive unit cell in (111) direction.
As shown in Fig. 4, the DMM conductance is larger than the maximum available conductance we can get by tuning the interfacial bonding and masses. This result seems to contradict the prediction by DMM as DMM underestimates the conductance when the modes at two sides of the interface are similar (Eq. (10)). However, DMM does not take interface properties into consideration. Since interface properties can severely degrade the thermal conductance, the resultant conductance ends up much lower than that predicted by DMM. Our system in study is a dimension-mismatched system; not all the modes in the contact reservoir can connect with the modes in the chain. Consequently, based on our simulations with various bonding and masses at the interfaces, the DMM conductance, in which all the modes at two sides are diffusively connected, is larger than the maximum available conductance.
Influence of Bonding and End-Group Mass
Comparison of Different End Groups.
Figure 5 shows the force constants for different end groups obtained from first principles calculations and the corresponding interfacial thermal conductance extracted from NEGF. The calculated trend of interfacial thermal conductances of different end groups is consistent with the experimental data from Losego et al. . The Au–thiol–alkane junctions have a larger interfacial bonding strength compared to Au–amine–alkane junctions, and a correspondingly larger interfacial thermal conductance. In this case, the interfacial thermal conductance follows the bonding strength. However, it is not conclusive because the bonding strengths are not large enough to attain maximum conductance.
The Au–thiol–alkane junction and Au–amine–alkane junction show different trends relating to the roughness of the gold surface. In the Au–amine–alkane junction, the reconstructed surface (adatom case) shows a larger conductance than that of the clean surface junction. In the Au–thiol–alkane junction, the interfacial thermal conductance is larger for the clean gold surface. This discrepancy can be explained by analyzing the bond strength. The strength of all the bonds in all directions for Au–amine–alkane junction with adatom is larger than that of Au–amine–junction without adatom, hence the conductance is larger in the former case. In contrast, the reconstructed Au–alkane–thiol junction only has bond in z-direction larger than the bonds of clean surface Au junction, which limits the coupling of transverse vibrational modes in gold to those within the alkane chain.
Note that the experimentally reported interfacial thermal conductance of a single Au–thiol–alkane junction is around 220 MW m−2 K−1  and of Au–alkane–thiol–Au junction with surface roughness correction is around 115 MW m−2 K−1 . Our result of a clean covalent bonding interface is around 210 MW m−2 K−1, which is consistent with the experimental works.
Role of Interfacial Parameters.
We use two different mechanisms to vary the interfacial bonding strength: varying the adsorption distance, or directly varying the bonding strength by a ratio λ relative to the bonding strength of the Au–alkane thiol junction. In the case of varying adsorption distance, we assume that the alkane chains are normal to the gold surface, with the S atom centered on the top of the equilateral triangle formed by 3 gold atoms. The alkane chain and the end group are moved collectively, and the structure is not optimized during the ground state or during the lattice dynamics simulations. The force constants are calculated by the small displacement method [43–45] by only moving the atoms around the interfaces. Alternately, we directly vary the interfacial bonding strength between the gold and the end group relative to its calculated bond strength at an adsorption distance of 2 Å by a ratio of λ. For adjusting the end-group mass, we assume once again that the bond stiffness is fixed at its 2 Å value, and then change the corresponding mass matrix of the end group.
In order to focus on the influence of bonding strength and mass of the end group attached to gold, we assume that there is no impedance mismatch between the alkane end C2H4 unit and the thiol end group, and set their bonding strength equal to the one between two C2H4 units. This would be satisfied by enforcing the ASR for the Sulfur atom.
Figure 6 shows the variation of bond strength K and thermal conductance G with varying adsorption distance between 1.6 Å and 2.2 Å. The interfacial bond strength shows an inverse proportionality to the adsorption distance. Interestingly, the corresponding interfacial thermal conductance reaches a maximum at a distance of 1.9 Å, which is close to the relaxed distance for alkanethiols self-assembled on gold. The existence of a maximum is independently validated when we directly vary the bond strength, as shown in Fig. 7(a). A maximum conductance is seen to arise when the bond strength is around 1.3 times that of the initial 2 Å distance value for thiol. When the bond strength continues to increase, the interfacial thermal conductance reduces and ultimately seems to saturate.
A maximum is also observed at a mass around 10 amu when we plot the thermal conductance against the mass of the end group. This mass is lower than the mass of the C atom, indicating that the vibrational dynamics of the hydrogen atoms cannot be ignored. For an impedance matched to maximize thermal transport across the interfaces, we would normally expect the preferred bonding mass to lie in between the masses of the two contacts. The maximizing mass also suggests that the prevailing Hautman–Klein model  in molecular dynamics method is used to simulate the alkane chains, where we neglect hydrogen atoms and add their masses instead to the carbon atoms, might not be suitable in the harmonic limit.
We can use the NEGF extracted mode count in the DMM formula to directly estimate the thermal conductance (Fig. 4), as discussed in Sec. 3. The shaded area in Fig. 4 represents the tunable range of the interfacial thermal conductance obtained by varying the interfacial strength and the end-group atomic mass. The range is predicted by the two approaches we described earlier in this section. Both approaches, as we cautioned, ignore the impedance between the alkane end C2H4 unit and the thiol anchoring group. Ignoring this contribution implies that the upper limit of the tunable range could well be smaller than the shaded area under the harmonic assumption. The solid black line in Fig. 4 represents the contribution from a clean gold surface optimally bonded to an alkanethiol. Comparing with the tunable range, we can see that, in the harmonic limit the interfacial thermal conductance can at most be enhanced by less than 30%. This maximum is still less than the conductance that DMM predicts for this system. Since DMM arises in the limit of complete dephasing, it seems that such dephasing events must break restrictive symmetry selection rules that would otherwise limit the bandwidth of the transmitting phonon modes across the two dissimilar materials. In reality, we expect the answer to be slightly lower than DMM, because true dephasing events retain partial memory, and also tends to reduce the average transmission per mode.
Three-Dimensional to One-Dimensional Generalized Results and Discussions.
To explore whether our conclusions generalize to dimensional mismatch, we also studied a simplified one degree-of-freedom 3D–1D interface with only one degree of mechanical freedom (Fig. 8(a)). In this model, the force constants exist in three dimensions to account for the transverse momenta in the 3D system but the atoms are only allowed to move in the transport direction. The mass of the atoms on the 3D side is set to that of gold, and the bond strength is set to match its cutoff frequency. Similarly, the atoms in the 1D chain have the mass of carbon while their bond strength matches the cutoff frequency of the alkane chain acoustic branches. The alkane channel is represented by a 1D channel between a 1D contact representing its extension on one side, and the 3D contact representing the Au substrate.
Figure 8(b) shows a map of the thermal conductance to visualize its dependence on interfacial parameters. The thermal conductance is scaled by the conductance in the classical limit. The figure confirms the existence of a maximum in the interfacial thermal conductance for a set of optimal interfacial parameters. We also see that the conductance is most sensitive to mass variations when the bond energy hits its maximizing value. In contrast, bond variations are critical when the end-group masses are relatively high. Note that the optimal bonding and mass for the maximum conductance is 35 N/m and 50 amu, respectively. These optimal values are different from those predicted by one-dimensional interfaces [13,14,17]. For instance, the optimal bonding strength predicted by the harmonic mean of a one-dimensional interface [14,17] is 25 N/m, which differs by about 50% from our result. The difference arises mainly from the 3D dimensionality of our metal contact, which serves as a reservoir of channels whose impedances span over a wide range of values. Then, the optimal properties for the 3D–1D interface result from a weighted average of the optimal properties of many 1D–1D interfaces between the 1D molecule and each of the channels of the reservoir contact.
The existence of a maximum conductance with interfacial parameters can be explained by considering the dynamics of the DOS as shown in Fig. 9. Figures 9(a), 9(c), and 9(e) correspond to varying the bond strength with the end-group mass set to 50 amu (up triangle symbols in Fig. 8(b)), while Figs. 9(b), 9(d), and 9(f) correspond to varying the end-group mass, with the bonding strength set to 35 N/m (square symbols in Fig. 8(b)). Let us first look at varying bonding strength. Figure 9(a) shows the DOS of the 3D and 1D contacts and the average local density of states (LDOS) of the device with various interfacial parameters. As the bonding strength between gold and end-group atom increases, the LDOS shifts toward the higher frequency states spanned by the 1D contact. This can be explained by the ASR that once the bond strength increases, the on-site energy on the end-group atom will also increase, resulting in an up-shifted spectrum. These high-frequency states cannot be used as transport states within a harmonic assumption, so that elastic transport cannot arise above the 3D cutoff frequency, in this case 30 Trad/s.
where D3D(ω) and D1D(ω) are the DOS of 3D and 1D contacts, respectively. Within the harmonic assumption, the ODOS of contacts should serve as a weighted window for the channel states. As the carbon–carbon bonding strength increases, one of the LDOS peaks gets pushed out of the 3D band, resulting in an eventual decrease in transmission and its corresponding nonmonotonicity with a preferred sweet spot. Comparing the DOS in Fig. 9(c) and transmission in Fig. 9(e), we see that the transmission approximately follows the LDOS of the devices within the weighted window. As described earlier, the difference between LDOS and transmission is the phonon group velocity. High-frequency phonons have lower velocity, which explains the suppressed transmission of the LDOS peaks around the cutoff frequency. Weaker bonding strength also creates lower phonon velocity, so that the lowest transmission arises for the weakest bonding strengths. We also observe that the conductance start saturating after the maximum conductance as the bonding increases. This is also due to the overlap of LDOS with ODOS, which starts saturating after the high-frequency peak moves out of the weighted window.
The transmission dependence on end-group mass has an origin similar to above, seen by comparing the transmission in Fig. 9(f) and the DOS in Fig. 9(d). Within the weighted window of the ODOS, the DOS for the lighter end-group atom is much smaller, making its transmission correspondingly weak. On the other hand, heavier end-group mass reduces the phonon group velocity especially at higher frequency, reducing their corresponding transmissions as well. The end result is once again a maximum in the conductance when varying the mass of the end-group atom.
In summary, using first principles parameters coupled with NEGF, we studied the influence of interfacial chemistry (bond strength and mass of the connector atoms) on the interfacial thermal conductance of gold–alkane junctions. Within the harmonic approximation, we predict the existence of a maximum conductance that is mirrored by a simpler 3D–1D model. We attribute the maximum conductance to the interplay between LDOS and phonon group velocity in determining the phonon transmission. This result complements the optimal set of interfacial properties predicted by one-dimensional harmonic studies and shows that in a 3D–1D dimension-mismatched system, the maximum conductance also exists. However, the optimal mass being geometric mean of masses or the optimal bonding being the harmonic mean of bonding predicted in purely 1D systems are no longer valid in 3D–1D systems. Our results of maximum conductance should be valid in SAMs systems with two contacts and finite length chains as long as transport is elastic. However, the conductance can change slightly due to the added resonance in the transmission . Meanwhile, interchain interactions can be another important factor that can influence thermal transport in such systems [48,49]. In applying our results to the SAMs system, the end groups should be chosen attentively to obtain the optimal mass and bonding. To set up a reference for the conductance, we developed a DMM model based on the number of modes of the materials at each side of the interface extracted from the DFT–NEGF calculations. This approach goes beyond the isotropic approximation and can be easily used in any studies performed by NEGF simulations.
J. Z., C. A. P., and A. W. G. acknowledge the School of Engineering and Applied Sciences at the University of Virginia that covered the registration fees for a Quantum ESPRESSO workshop. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) (DMR 130123) , which is supported by National Science Foundation grant number ACI-1053575.
Division of Chemical, Bioengineering, Environmental, and Transport Systems (IDR-CBET 1134311).
Division of Electrical, Communications and Cyber Systems (CAREER-QMHP 1028883).
- A =
- D =
density of states
- g =
surface Green's function
- G =
interfacial thermal conductance, MW m−2 K−1
retarded Green's function
reduced Planck's constant
- K =
force constant, N/m
- LDOS =
local density of states
- M =
number of modes
- M =
mass matrix, kg
- N =
- ODOS =
overlap of density of states
- q =
- T =
mode-averaged transmission coefficient
- Γ =
- Σ =
- ω =