Heat dissipation and its relation to thermopower in single-molecule junctions

The molecular junctions based on short molecules that we investigate in this work are excellent examples of systems where the transport properties are dominated by elastic tunneling events. In this case, conduction electrons preserve their quantum mechanical coherence along the junction region. Such systems can be described theoretically within the framework of the Landauer approach [30] that has become one of the central pillars of theoretical nanoelectronics [31]. The key idea of this approach is that, if the electron transport in a nanojunction is dominated by elastic processes, it can be modeled as a scattering problem. In this problem the junction electrodes are assumed to be ideal reservoirs where electrons thermalize and acquire a well-defined temperature, while the central region plays the role of a scattering center. As a result, all the transport properties of these systems are completely determined by the transmission function τ(E,V ), which describes the total probability for electrons to cross the junction at a given energy E when a voltage bias V is applied across the system. Landauer's original approach was extended by several authors to describe both thermal transport and thermoelectricity in nanojunctions [32–35], and in this section we explain how this approach can be also used to describe heat dissipation in atomic-scale circuits.

Let us consider a two-terminal device with source S and drain D electrodes linked by a nanoscale junction⁸. If the electron transport is elastic, i.e. if electrons flow through the device without exchanging energy in the junction region, then obviously their excess energy—acquired by the application of a bias voltage—must be released in the electrodes. From simple thermodynamical arguments, the amount of heat released per unit of time in an electrode with electrochemical potential μ is given by [36]

$$\begin{equation} Q = \frac{\mu}{e} I - I_E , \end{equation} \tag{ 1 }$$

where I and I_E are the charge and energy current, respectively, and e > 0 is the electron charge. Within the Landauer theory, the electronic contribution to charge and energy currents in a two terminal device are expressed in terms of the transmission function as

$$\begin{equation} I(V) = \frac{2e}{h} \int^{\infty}_{-\infty} \tau(E,V) [ f_{\mathrm{S}}(E,\mu_{\mathrm{S}}) - f_{\mathrm{D}}(E,\mu_{\mathrm{D}}) ]\,\mathrm{d}E, \end{equation} \tag{ 2 }$$

$$\begin{equation} I_E(V) = \frac{2}{h} \int^{\infty}_{-\infty} E \tau(E,V) [ f_{\mathrm{S}}(E,\mu_{\mathrm{S}}) - f_{\mathrm{D}}(E,\mu_{\mathrm{D}}) ]\,\mathrm{d}E , \end{equation} \tag{ 3 }$$

where f_S,D are the Fermi functions of the source or the drain and μ_S,D are the corresponding electrochemical potentials such that μ_S − μ_D = eV . The prefactor 2 is due to the spin degeneracy that will be assumed throughout our discussion. Using equations (2) and (3) we arrive at the following expressions for the power (heat per unit of time) dissipated in the source and the drain electrodes:

$$\begin{equation} Q_{\mathrm{S}}(V) = \frac{2}{h} \int^{\infty}_{-\infty} (\mu_{\mathrm{S}} -E) \tau(E,V) [ f_{\mathrm{S}}(E,\mu_{\mathrm{S}}) - f_{\mathrm{D}}(E,\mu_{\mathrm{D}}) ]\,\mathrm{d}E, \end{equation} \tag{ 4 }$$

$$\begin{equation} Q_{\mathrm{D}}(V) = \frac{2}{h} \int^{\infty}_{-\infty} (E -\mu_{\mathrm{D}}) \tau(E,V) [ f_{\mathrm{S}}(E,\mu_{\mathrm{S}}) - f_{\mathrm{D}}(E,\mu_{\mathrm{D}}) ]\,\mathrm{d}E . \end{equation} \tag{ 5 }$$

Obviously, the sum of the powers dissipated in both electrodes must be equal to the total power injected in the junction Q_Total = IV , as can be seen from equations (4) and (5):

$$\begin{equation}\fl Q_{\mathrm{S}}(V) + Q_{\mathrm{D}}(V) = \frac{2eV}{h} \int^{\infty}_{-\infty} \tau(E,V) [ f_{\mathrm{S}}(E,\mu_{\mathrm{S}}) - f_{\mathrm{D}}(E,\mu_{\mathrm{D}}) ]\,\mathrm{d}E = IV = Q_{\rm Total}(V) . \end{equation} \tag{ 6 }$$

Equations (4) and (5) are the central results of the Landauer theory of heat dissipation and will be used throughout this work. Given their relevance, it is worth explaining the underlying physical picture. For this purpose, we shall follow here [37]. As shown in figure 1(a) for the case of a molecular contact, when a bias voltage is applied to a junction, the electrochemical potentials of the source and the drain are shifted. This opens up an energy window for electrons to cross the junction and it results in a net electron current in the system. Let us consider an electron with an energy E which is transmitted elastically through the junction region from the source to the drain. After reaching the drain, see figure 1(b), the electron undergoes inelastic processes (via electron–phonon interactions) and it decays to the electrochemical potential of the drain. In these processes, the electron releases an energy E − μ_D in the drain. On the source side, the original electron leaves behind a hole that is filled up by the same type of inelastic processes. In this way, an energy μ_S − E is dissipated in the source. Notice that the overall energy dissipated in the electrodes is equal to μ_S − μ_D = eV , which is exactly the work done by the external battery on each charge.

Zoom In Zoom Out Reset image size

Equations (4) and (5) describe in quantitative terms the argument just explained. The energy integrals take into account the contribution of all the electrons in the transport window, the transmission function accounts for the finite probability of an electron to tunnel through the junction, and the Fermi functions account for the occupation probabilities of the initial and final states in the tunneling events. In particular, the appearance of the difference of the Fermi functions is a result of the net balance between tunneling processes transferring electrons from the source to the drain and from the drain to the source.

The physical argument discussed above also allows us to answer in simple terms a central question in this work: Is the heat equally dissipated in both electrodes? In general, the answer is no, as we proceed to explain. Let us assume that the transmission function is energy-dependent and, for instance, that the transmission is higher for energies in the lower part of the transport window, as shown in figure 1(c). As we shall discuss below, this corresponds to a situation where the transport is dominated by the highest occupied molecular orbital (HOMO). In this case, for a positive bias, it is obvious that a larger portion of the energy is dissipated in the source because it is more probable that electrons with a smaller excess energy tunnel through the molecule and release their energy in the drain. Obviously, if we reversed the transmission landscape, see figure 1(d), the heat would be preferentially dissipated downstream following the electron flow, i.e. in the drain. To conclude this subsection, we would like to remark that these arguments also make clear the fact that, in general, heat dissipation in a given electrode depends on the bias polarity, i.e. it depends on the current direction.

As explained above, within the Landauer approach the heat dissipated in the electrodes of a nanoscale junction is fully determined by its transmission characteristics. From equations (4) and (5) one can draw some general conclusions about the symmetry of the heat dissipation between electrodes and with respect to the bias polarity. Let us first discuss under which circumstances heat is equally dissipated in both electrodes. For this purpose, we assume from now on (without loss of generality) that all the energies in the problem are measured with respect to the equilibrium electrochemical potential of the system, which we set to zero (μ = 0). Moreover, we assume that the electrochemical potentials are shifted symmetrically with the bias voltage, i.e. μ_S = eV/2 and μ_D = −eV/2. This means in practice that the energy E is measured with respect to the center of the transport window. With this choice, the power dissipated in the source electrode is given by

$$\begin{equation} Q_{\mathrm{S}}(V) = \frac{2}{h} \int^{\infty}_{-\infty} (eV/2 -E) \tau(E,V) [ f(E-eV/2) - f(E+eV/2) ]\,\mathrm{d}E, \end{equation} \tag{ 7 }$$

where $f(E) = 1/[1 + \exp (E/k_{\mathrm { B}}T)]$ is the Fermi function. Using the change of variable E → − E, the previous expression becomes

$$\begin{equation} Q_{\mathrm{S}}(V) = \frac{2}{h} \int^{\infty}_{-\infty} (E+eV/2) \tau(-E,V) [ f(E-eV/2) - f(E+eV/2) ]\,\mathrm{d}E, \end{equation} \tag{ 8 }$$

where we have used the relation f(− E) = 1 − f(E). This latter equation must be compared with the corresponding expression for the power dissipated in the drain, which now reads

$$\begin{equation} Q_{\mathrm{D}}(V) = \frac{2}{h} \int^{\infty}_{-\infty} (E + eV/2) \tau(E,V) [ f(E-eV/2) - f(E+eV/2) ]\,\mathrm{d}E . \end{equation} \tag{ 9 }$$

Thus, we conclude that a sufficient condition to have equal heat dissipation in both electrodes, i.e. Q_S(V ) = Q_D(V ), is that the transmission function fulfills that τ(E,V ) = τ(− E,V ). This means that if the transport is electron–hole symmetric, then heat is equally dissipated in the source and in the drain. A simple corollary of this condition is that if the transmission is energy-independent in the transport window, then the power dissipation is the same in both electrodes. In particular, if the junction is ballistic, i.e. if τ(E,V ) = M, M being an integer equal to the number of open conduction channels in the system, then heat is equally dissipated in both electrodes. Notice that these conclusions can be also drawn from the handwaving arguments explained at the end of the previous subsection.

Another crucial question that can be answered in general terms is: Does the power dissipated in a given electrode depend on the bias polarity (or current direction)? With manipulations similar to those of the previous paragraph, it is straightforward to show that Q_S,D(V ) = Q_S,D(− V ) if τ(E,V ) = τ(− E,−V ) is satisfied. This means, in particular, that if the transmission does not depend significantly on the bias voltage, the power dissipated in an electrode is independent of the bias polarity as long as there is electron–hole symmetry. Obviously, if the transmission is independent of both the energy and the bias, then the heating is symmetric with respect to the inversion of the bias polarity. Again, this is what occurs in any ballistic structure.

On the other hand, one can also show that the relation Q_S,D(V ) = Q_D,S(− V ) is satisfied if τ(E,V ) = τ(E,−V ). This relation implies the following one: Q_S,D(V ) + Q_S,D(− V ) = Q_S(V ) + Q_D(V ) = Q_Total(V ) = IV . This means that if the transmission is symmetric with respect to the inversion of the bias voltage, then the sum of the dissipated powers for positive and negative bias in a given electrode is equal to the total power dissipated in the junction. This relation is then expected to hold for left–right symmetric junctions, where the voltage profile is inversion symmetric with respect to the center of the junction.

From the general considerations above, it is clear that in order to have a heating asymmetry of the type Q_S(V ) ≠ Q_D(V ) or Q_S,D(V ) ≠ Q_S,D(− V ), one needs a certain degree of electron–hole asymmetry in the transmission function. This can be seen explicitly by expanding equations (7) and (9) to first order in the bias voltage, which leads to the following expression:

$$\begin{equation}\fl Q_{\mathrm{S,D}}(V) = \pm \frac{2eV}{h} \int^{\infty}_{-\infty} E \tau(E,V=0) \frac{\partial f(E,T)} {\partial E} \mathrm{d}E \, + \, O(V^2) = \pm GTSV \, + \, O(V^2) . \end{equation} \tag{ 10 }$$

Here, G is the linear electrical conductance of the junction, T is the absolute temperature and S is the thermopower or Seebeck coefficient of the junction. Let us remind that in the Landauer approach G and S can be expressed in terms of the zero-bias transmission as

$$\begin{equation} G = \frac{2e^2}{h} \int^{\infty}_{-\infty} \tau(E,V=0) [ - \partial f(E,T)/\partial E] \; \mathrm{d}E, \end{equation} \tag{ 11 }$$

$$\begin{equation} S = -\frac{1}{eT} \frac{\int^{\infty}_{-\infty} E \tau(E,V=0) [-\partial f(E,T)/\partial E] \; \mathrm{d}E}{\int^{\infty}_{-\infty} \tau(E,V=0) [-\partial f(E,T)/\partial E] \; \mathrm{d}E}. \end{equation} \tag{ 12 }$$

Equation (10) can also be written in the familiar form: Q_S,D(V ) = ± ΠI [38], where Π = TS is the Peltier coefficient of the junction and I = GV is the electrical current in the linear regime. Moreover, from equation (10) we can conclude that

$$\begin{equation} Q_{\mathrm{S}}(V) - Q\mathrm{_{D}}(V) = \pm 2GTSV + O(V^2), \end{equation} \tag{ 13 }$$

$$\begin{equation} Q_{\mathrm{S,D}}(V) - Q_{\mathrm{S,D}}(-V) = \pm 2GTSV + O(V^3) , \end{equation} \tag{ 14 }$$

i.e. the heating asymmetries between electrodes and with respect to the bias polarity are determined to first order in the bias by the Seebeck coefficient of the junction.

On the other hand, using the low temperature expansion of G and S

$$\begin{equation} G = \frac{2e^2}{h} \tau(E_{\rm F},V=0) \;\;\; \mathrm{and} \;\;\; S = - \frac{\pi^2 k^2_{\rm B}T}{3e} \frac{\tau^{\prime}(E_{\rm F}, V=0)}{\tau(E_{\rm F}, V=0)} , \end{equation} \tag{ 15 }$$

where τ'(E_F,V =0) is the energy derivative of the zero-bias transmission at the Fermi energy, E_F = 0, we arrive at the following expression for the linear voltage term of the power dissipations at low temperature:

$$\begin{equation} Q_{\mathrm{S,D}}(V) = \mp \left( \frac{2e}{h} \right) \frac{\pi^2}{3} (k_{\rm B}T)^2 \tau^{\prime}(E_{\rm F}, V=0) V \, + \, O(V^2). \end{equation} \tag{ 16 }$$

Thus, in the linear response regime the slope of the transmission function at the Fermi energy determines not only the thermopower but also the magnitude and the sign of the heating asymmetry.

Equation (10) implies that at sufficiently low bias, and if the thermopower is finite, one of the electrodes may be cooled down by the passage of an electrical current (Peltier effect). How low must the voltage be to observe this cooling effect? To answer this question we need to consider the quadratic term of Q_S,D(V ) in the bias. (Notice that a quadratic term must exist in order to satisfy energy conservation, Q_S(V ) + Q_D(V ) = Q_Total(V ), since at low bias Q_Total(V ) = GV².) If for simplicity we ignore the bias dependence of the transmission function, it is easy to see that the second-order term of Q_S,D(V ) is given by (1/2)GV², which is equal to half of the total power in the linear regime. This second-order term dominates over the linear one, see equation (10), when |V | > 2T|S| (above this bias, each electrode is heated up by the current). Thus for instance, if we assume room temperature (T = 300 K) and a typical value of |S| = 10 μV K⁻¹ [4], then the second-order term dominates the contribution to the heating for voltages |V | > 6 mV. This means in practice that in most molecular junctions we expect joule heating to dominate over Peltier cooling over a wide range of bias. This also means that, in general, in order to correctly describe the heat dissipation, one needs to go beyond the linear response theory. This can be easily done by using the full expressions for the power dissipated, see equations (7) and (9), once the transmission function is known. It is worth stressing that throughout this work we have used these full expressions to account for all the nonlinear contributions to the power dissipation in both electrodes. In particular, in order to correctly describe the experimental results of [28], it is absolutely necessary to go well beyond linear response (see section 4).

In order to illustrate some of the ideas discussed in the previous subsections and to gain some further insight, we analyze in this subsection the heat dissipation in molecular junctions with the help of a single-level model, sometimes referred to as resonant tunneling model [1, 39]. In this model one assumes that transport through a molecular junction is dominated by a single molecular orbital (typically the HOMO or the lowest unoccupied molecular orbital (LUMO)) and the transmission function is given by

$$\begin{equation} \tau(E) = \frac{\Gamma^2}{(E-\epsilon_0)^2 + \Gamma^2}. \end{equation} \tag{ 17 }$$

Here, ₀ is the energetic position of the molecular level measured with respect to the Fermi energy (set to zero) and Γ is the strength of the metal–molecule coupling (or equivalently the level broadening). For simplicity, we assume here that the junction is symmetric, i.e. the molecule is equally coupled to both electrodes. Let us recall that in this model the parameter Γ is taken as a constant (energy-independent), which means that local density of states in the metal electrodes is assumed to be constant in the energy range involved in the transport window. Let us also recall that in the symmetric case considered here, the level position does not depend on the bias and the transmission function is thus independent of the applied voltage.

Within this model, the power dissipated in the source and in the drain electrodes is obtained by substituting equation (17) into equations (7) and (9). Let us emphasize that these formulae account for any nonlinear contribution as a function of the bias voltage. Moreover, since the transmission does not depend explicitly on the bias voltage, the following two relations are satisfied (see section 2.2):

$$\begin{equation} Q_{\mathrm{S}}(V) = Q_{\mathrm{D}}(-V) \;\;\; \mathrm{and} \;\;\; Q_{\mathrm{S,D}}(V) + Q_{\mathrm{S,D}}(-V) = Q_{\rm Total}(V) = IV . \end{equation} \tag{ 18 }$$

The first relation tells us that the power dissipated in one of the electrodes can be obtained from the power dissipated in the other one by simply inverting the bias. The second relation implies that the sum of the power dissipated for positive and negative bias in a given electrode is equal to the total power dissipated in the junction.

Let us now use this model to simulate a junction where the transport is dominated by the LUMO. For this purpose, we choose ₀ = + 1 eV (level above the Fermi energy) and Γ = 40 meV. With these values the transmission at the Fermi energy is τ(E_F) = 1.59 × 10⁻³. In figure 2 we show for this example the transmission as a function of energy (panel (a)), the current–voltage (I–V ) characteristics (panel (b)), the total power and power dissipated in the source as a function of the bias (panel (c)), and the power dissipated in the source as a function of the total power for both positive and negative bias (panel (d)). The voltage range in panels (b) and (c) has been chosen so as to describe the typical regime explored in most experiments in which the highest bias is not sufficient to reach the resonant condition |eV | = 2|₀|. The results of panels (b)–(d) have been calculated at room temperature (300 K), but in this off-resonant situation the corresponding results at zero temperature are similar (not shown here). The main conclusion of these results is that, as it can be clearly seen in panel (d), the power dissipated in the source for negative bias is higher than for positive bias (the situation is the opposite in the drain). According to equation (14), this fact is due to the negative value of the Seebeck coefficient, or equivalently to the positive slope of the transmission at the Fermi energy, see equation (16). In this case, the nonlinearities of the I–V curve are reflected in the fact that the total power contains a significant quartic term (∝V⁴) and the bias asymmetry in the power dissipated in the source, Q_S(V ) − Q_S(− V ), contains a significant cubic term (∝V³) and even higher order ones.

Zoom In Zoom Out Reset image size

To illustrate how the sign of the Seebeck coefficient is closely related to the asymmetries in the heat dissipation, we now investigate a case where we just change the sign of ₀ with respect to the previous example. This case corresponds to a situation where the transport is dominated by the HOMO. The corresponding results are shown in figure 3. As one can see, while the I–V curve is exactly the same as in the previous case, now a higher power is dissipated in the source for positive bias contrary to the example in figure 2. These results nicely illustrate how the heat dissipation can be controlled by tuning the thermopower of the junctions: the power dissipation is always higher in the electrode with the highest electrochemical potential if the level lies below the Fermi level (in equilibrium), while it is higher in the electrode with the lowest electrochemical potential if the level is above it.

Zoom In Zoom Out Reset image size

The relation between Q_S and Q_Total is particularly useful to visualize the heating asymmetries (both between electrodes and with respect to the bias polarity). In this sense, it is interesting to analyze to what extent this relation depends on variations of the parameters of the model. In this analysis we focus on the case of off-resonant transport, which is the most common situation in the experiments. In figure 4(a) we show this relation for Γ = 40 meV and different values of the level position ranging from 0.4 to 1.4 eV. Notice that in all cases we are in the off-resonant regime (|₀| ≫ Γ). As one can see, the Q_S–Q_Total relation is rather insensitive to variations in the level position, in spite of the fact that the linear conductance changes by up to an order of magnitude. In figure 4(b) we explore how this relation depends on the coupling strength. For this purpose, we have fixed the level position to ₀ = + 1 eV and varied Γ between 20 and 60 meV. Notice that the relation between the powers is more sensitive to variations in Γ, but again the variations are relatively small taking into account that the conductance has been varied by a factor of 10.

Zoom In Zoom Out Reset image size

Now, let us try to shed some more light on the insensitivity of the Q_S–Q_Total relation to variations of the junction parameters in an off-resonant transport situation. For this purpose, we have carried out analytical calculations in the limit of zero temperature, which are also applicable to finite temperatures in an off-resonant situation. From equations (7), (9) and (17), it is straightforward to show that

$$\begin{eqnarray} Q_{\mathrm{S}}(V) & = & \frac{2}{h} \Gamma (eV/2-\epsilon_0) \left\{ \arctan \left[ \frac{eV/2-\epsilon_0}{\Gamma} \right] + \arctan \left[ \frac{eV/2+\epsilon_0}{\Gamma} \right] \right\}\nonumber\\ & & - \frac{\Gamma^2}{h} \left\{ \ln \left[ 1 + \left( \frac{eV/2-\epsilon_0}{\Gamma} \right)^2 \right] - \ln \left[ 1 + \left( \frac{eV/2+\epsilon_0}{\Gamma} \right)^2 \right] \right\}, \end{eqnarray} \tag{ 19 }$$

$$\begin{equation} Q_{\rm Total}(V) = \frac{2e \Gamma V}{h} \left\{ \arctan \left[ \frac{eV/2-\epsilon_0}{\Gamma} \right] + \arctan \left[ \frac{eV/2+\epsilon_0}{\Gamma} \right] \right\}. \end{equation} \tag{ 20 }$$

To establish the relation between these two powers, we first do a Taylor expansion in the bias of both expressions and focus on the off-resonant situation (|₀| ≫ Γ). Thus, the total power can be written in this limit as

$$\begin{equation} Q_{\rm Total}(V) \approx GV^2 + \frac{e^2 G}{4\epsilon^2_0} V^4 + O(V^6) , \end{equation} \tag{ 21 }$$

while the heating asymmetry adopts the form

$$\begin{equation} Q_{\mathrm{S}}(V) - Q_{\mathrm{D}}(V) = Q_{\mathrm{S}}(V) - Q_{\mathrm{S}}(-V) \approx \frac{eG}{3\epsilon_0} V^3 + \frac{e^3G}{10 \epsilon^3_0} V^5 + O(V^7). \end{equation} \tag{ 22 }$$

Here, G is the linear conductance, which in the off-resonant limit reads G ≈ G₀(Γ/₀)², where G₀ = 2e²/h is the conductance quantum. If we now restrict ourselves to the low-bias limit where the total power is quadratic (Q_Total(V ) ≈ GV²), then the power dissipated in the source can be expressed as [Q_S = Q_Total/2 + (Q_S − Q_D)/2]

$$\begin{equation} Q_{\mathrm{S}}(Q_{\rm Total}) \approx \left\{ \begin{array}{@{}cc@{}} \displaystyle\frac{1}{2} Q_{\rm Total} + \mathrm{sgn}(\epsilon_0) \displaystyle\frac{e}{6G^{1/2}_0} \displaystyle\frac{1}{\Gamma} Q^{3/2}_{\rm Total} & {\rm (for \ negative \ bias),} \\ \displaystyle\frac{1}{2} Q_{\rm Total} - \mathrm{sgn}(\epsilon_0) \displaystyle\frac{e}{6G^{1/2}_0} \displaystyle\frac{1}{\Gamma} Q^{3/2}_{\rm Total} & {\rm (for \ positive \ bias)} \end{array} \right. \end{equation} \tag{ 23 }$$

This expression shows that the relation between Q_S and Q_Total is independent of the level position in an off-resonant situation. Notice that the level position only enters via its sign (positive for LUMO-dominated cases and negative for HOMO-dominated ones). This analytical expression nicely illustrates the insensitivity of this relation to the level position found in the numerical results above.

The goal of this subsection is to describe the first ab initio method, which is based on density functional theory (DFT). This DFT-based transport method has been described in great detail in [40] and therefore, our discussion here will be rather brief. This method is built upon the quantum-chemistry code TURBOMOLE [41]. A central issue in our approach is the description of the electronic structure of the junctions within DFT. In all the calculations presented here we have used the BP86 exchange-correlation functional [42, 43] and the Gaussian basis set def-SVP [44]. The total energies were converged to a precision of better than 10⁻⁶ atomic units, and structure optimizations were carried out until the maximum norm of the Cartesian gradient fell below 10⁻⁴ atomic units. In what follows, we shall discuss the two main steps in our method: (i) construction of the geometries of the single-molecule junctions; and (ii) determination of the electronic structure of the junctions and calculation of the transmission function using Green's function techniques.

In the construction of the junction geometries we proceed as follows. We first place the relaxed molecule in between two gold clusters with 20 (or 19) atoms. Subsequently, we perform a new geometry optimization by relaxing the positions of all the atoms in the molecule as well as the four (or three) gold atoms on each side that are closest to it, while the other gold atoms are kept fixed. Afterwards, the size of the gold cluster is increased to about 63 atoms on each side in order to describe correctly both the metal–molecule charge transfer and the energy level alignment. Finally, the central region, consisting of the molecule and one or two Au layers on each side, is coupled to ideal gold surfaces, which serve as infinite electrodes and are treated consistently with the same functional and basis set within DFT.

To determine the electronic structure of the infinite junctions and to compute the transmission characteristics we make use of a combination of Green's function techniques and the Landauer formula expressed in a local nonorthogonal basis. Briefly, the local basis allows us to partition the basis states into L, C and R ones, according to the division of the contact geometry, see figure 5. Thus, the Hamiltonian (or single-particle Fock) matrix H, as well as the overlap matrix S, can be written in the block form

$$\begin{equation} {\bf H}=\left(\begin{array}{@{}ccc@{}} {\bf H}_{\rm LL} & {\bf H}_{\rm LC} & {\bf 0}\\ {\bf H}_{\rm CL} & {\bf H}_{\rm CC} & {\bf H}_{\rm CR}\\ {\bf 0} & {\bf H}_{\rm RC} & {\bf H}_{\rm RR}\end{array}\right). \end{equation} \tag{ 24 }$$

Within the Landauer approach, the low-temperature conductance is given by G = G₀τ(E_F). It is important to remark that we shall restrict ourselves to the case of zero bias throughout this discussion. The energy-dependent transmission τ(E) can be expressed in terms of the Green's functions as [1]

$$\begin{equation} \tau(E) = \mathrm{Tr} \left[{\boldsymbol{\Gamma}}_{\rm L} {\bf G}_{\rm CC}^{\rm r} {\boldsymbol{\Gamma}}_{\rm R} {\bf G}_{\rm CC}^{\rm a} \right], \end{equation} \tag{ 25 }$$

where the retarded Green's function is given by

$$\begin{equation} {\bf G}_{\rm CC}^{\rm r}(E) = \left[ E {\bf S}_{\rm CC} - {\bf H}_{\rm CC} - {\boldsymbol{\Sigma}}_{\rm L}^{\rm r}(E) - {\boldsymbol{\Sigma}}_{\rm R}^{\rm r}(E)\right]^{-1} , \end{equation} \tag{ 26 }$$

and ${\bf{ G}}_{\mathrm {CC}}^{\mathrm { a}} = [{\bf{ G}}_{{\rm CC}}^{\rm r}]^{\dagger }$. The self-energies in the previous equation adopt the form

$$\begin{equation} {\boldsymbol{\Sigma}}_{X}^{\rm r}(E) = \left({\bf H}_{{\rm C}X} - E {\bf S}_{{\rm C}X} \right) {\bf g}_{XX}^{\rm r}(E) \left({\bf H}_{X{\rm C}} - E {\bf S}_{X{\rm C}} \right). \end{equation} \tag{ 27 }$$

On the other hand, the scattering rate matrices that enter the expression of the transmission are given by ${\boldsymbol{\Gamma }}_{X}(E) = -2\,\mathrm {Im} \left [{\boldsymbol {\Sigma }}_{X}^{\mathrm { r}}(E) \right ]$, and g^r_XX(E) = (E S_XX −  H_XX)⁻¹ are the electrode Green's functions with X = L,R.

Zoom In Zoom Out Reset image size

In order to describe the transport through the contact shown in figure 5(a), we firstly extract H_CC and S_CC and the matrices H_CX and S_CX from a DFT calculation of the extended center cluster (ECC) in figure 5(b). The blue-shaded atoms in regions L and R of figure 5(a) are assumed to be those coupled to the C region. Thus, H_CX and S_CX, obtained from the ECC, serve as the couplings to the electrodes in the construction of Σ^r_X(E).

On the other hand, the electrode Green's functions g^r_XX(E) in equation (27) are modeled as surface Green's functions of ideal semi-infinite crystals. To obtain these Green's functions, we first compute separately the electronic structure of a spherical fcc gold cluster with 429 atoms. Secondly, we extract the Hamiltonian and overlap matrix elements connecting the atom in the origin of the cluster with all its neighbors. Then, we use these bulk parameters to model a semi-infinite crystal which is infinitely extended perpendicular to the transport direction. The surface Green's functions are then calculated from this crystal with the help of the decimation technique [40, 45]. In this way we describe the whole system consistently within DFT, using the same nonorthogonal basis set and exchange-correlation functional everywhere.

In order to cross-check the validity of our conclusions on heat dissipation in single-molecule junctions, we have employed a second method to compute the transmission characteristics. This method is known as self-energy corrected DFT + Σ and it is an extension of the DFT method that has been introduced to cure some of its known deficiencies [29]. It is well-known that due to self-interaction errors in the standard exchange-correlation functionals and image charge effects, DFT-based methods have difficulties to accurately describe the energy gap and level alignment of molecules at surfaces [46]. In particular, DFT tends to underestimate the HOMO–LUMO gap, which in the context of molecular transport junctions implies that this method tends to systematically overestimate the conductance. The DFT + Σ approach, which we describe in this subsection, has been shown to improve the agreement with the experiments for both conductance and thermopower [11, 21, 47, 48].

In our implementation of the DFT + Σ method we have followed [47] (see also  [27]). This method aims at constructing accurate quasiparticle energies and starts by correcting the DFT-based energy levels of the gas phase molecules. To be precise, the HOMO energy is shifted such that it corresponds to the negative ionization potential (IP), while the LUMO is shifted to agree with the negative electron affinity (EA). All other energies of occupied (unoccupied) levels are shifted uniformly by the same amount as the HOMO (LUMO). The IP and EA are computed within DFT from total energy calculations as follows:

$$\begin{equation} \mathrm{IP} = E(Q=+e) - E(Q=0) \;\;\; \mathrm{and} \;\;\; \mathrm{EA} = E(Q=0) - E(Q=-e) , \end{equation} \tag{ 28 }$$

where E(Q = 0) is the total energy of the neutral molecule and E(Q = + e) (E(Q = −e)) is the total energy of the molecule with one electron removed (or added).

These corrected levels are, in turn, shifted when the molecule is brought into the junction. In particular, image charge interactions shift the energy of the occupied states up in energy and the virtual (unoccupied) states down in energy [49]. The key idea in the DFT + Σ method is that the screening of the metallic electrodes can be described classically as the interaction of point charges (in the molecule) with two perfectly conducting infinite surfaces. The idea goes as follows. First, from the Hamiltonian H_CC of the ECC, see equation (24), and the corresponding overlap matrix S_CC, we extract the matrices H_mol and S_mol corresponding to the basis functions on the molecular atoms. Then, this molecular Hamiltonian is diagonalized, ${\bf{ H}}_{\mathrm {mol}} \skew3\vec {c}_m = \epsilon _m {\bf{ S}}_{\mathrm {mol}} \skew3\vec {c}_m$, to obtain the eigenenergies, _m, and the eigenvectors, $\skew3\vec {c}_m$, for the molecule in the junction. Now, this information is used to construct a point charge distribution associated to a given molecular orbital m: $\rho _m(\vec {r}) = \sum ^N_{i=1} Q^m_i \delta (\vec {r} - \vec {r}_i)$, where the N point charges Q^m_i are located at the positions $\vec {r}_i$ of the atoms of the molecule in the junction. The contribution Q^m_i for a given molecular orbital m can be obtained from the Mulliken populations for the selected molecular orbital at a given atom i as $Q^m_i = \sum _{\mu \in i} ({\bf{ D}}^m {\bf{ S}}_{\mathrm {mol}})_{\mu \mu }$, where the sum runs over the local atomic basis states |μ〉 belonging to atom i and D^m_μν = c_μmc_νm is the density matrix of molecular orbital m.

The potential energy, Δ_m, of a point charge distribution $\rho _m(\vec {r})$ placed between two perfectly conducting planes located at z = ± L/2 is given by (in atomic units)

$$\begin{eqnarray}&&\fl \Delta_m = \frac{1}{2} \sum^N_{i,j=1} Q^m_i Q^m_j \sum^{\infty}_{n=1} \left( \frac{-1}{\sqrt{ [z_i + z_j \!+\! (2n-1)L ]^2 + |\vec{r}^{\parallel}_i - \vec{r}^{\parallel}_j|^2}}\right.\nonumber\\ &&\left.+ \frac{-1}{\sqrt{[z_i + z_j - (2n-1)L ]^2 + |\vec{r}^{\parallel}_i - \vec{r}^{\parallel}_j|^2}} + \frac{1}{\sqrt{[z_i + z_j + 2nL ]^2 + |\vec{r}^{\parallel}_i - \vec{r}^{\parallel}_j|^2}} \right.\nonumber\\ &&\left.+ \frac{1}{\sqrt{[z_i + z_j - 2nL ]^2 + |\vec{r}^{\parallel}_i - \vec{r}^{\parallel}_j|^2}} \right), \end{eqnarray} \tag{ 29 }$$

where $\vec {r}^{\parallel }$ is the component of $\vec {r} = \vec {r}^{\parallel } + z \hat z$ parallel to the planes.

Following [27], we use the charge distribution of the HOMO to calculate the image charge correction, Δ_occ = −Δ_HOMO, for all the occupied states, whereas we use the LUMO charge distribution to determine the correction Δ_virt = Δ_LUMO for the virtual or unoccupied states. The position of the image plane is chosen to be 1.47 Å outside of the first unrelaxed gold layer of the left and right electrode. This, however, constitutes an approximation as we neglect the screening introduced by the apex Au atoms and the position of 1.47 Å is valid only for a perfectly flat Au surface [11, 50].

So finally, the energy shift applied to all the occupied states is given by Σ_occ = −IP − _H + Δ_occ, and the corresponding shift for the unoccupied ones by Σ_virt = −EA − _L + Δ_virt. Here, _H and _L are the Kohn–Sham HOMO and LUMO energies from the gas-phase calculation. With these shifts we obtain a modified molecular Hamiltonian, $\tilde {\bf{ H}}_{\mathrm {mol}}$, which replaces H_mol in H_CC. Then, the method proceeds exactly as in the DFT case to compute the transmission of the junctions.