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

Motivated by recent experiments [Lee et al. Nature 498, 209 (2013)], we present here a detailed theoretical analysis of the Joule heating in current-carrying single-molecule junctions. By combining the Landauer approach for quantum transport with ab initio calculations, we show how the heating in the electrodes of a molecular junction is determined by its electronic structure. In particular, we show that in general the heat is not equally dissipated in both electrodes of the junction and it depends on the bias polarity (or equivalently on the current direction). These heating asymmetries are intimately related to the thermopower of the junction as both these quantities are governed by very similar principles. We illustrate these ideas by analyzing single-molecule junctions based on benzene derivatives with different anchoring groups. The close relation between heat dissipation and thermopower provides general strategies for exploring fundamental phenomena such as the Peltier effect or the impact of quantum interference effects on the Joule heating of molecular transport junctions.


Introduction
The advent of experimental techniques like the scanning tunneling microscope and break junctions have enabled the creation of single-molecule junctions and the study of their transport properties. This has triggered the hope to use single molecules as building blocks in novel nanoscale electronic devices, devices that could take advantage of the unlimited physical properties of molecules. In turn, this has also given rise to the research field of molecular electronics [1]. Although there are still many basic experimental challenges to be resolved before molecular electronics can become a viable technology, it is clear by now that molecular junctions constitute an excellent playground to test basic concepts of quantum transport, which could then be applied to novel charge and energy nanoscale devices. Thus for instance, the exhaustive study of electronic transport in single-molecule junctions has led to remarkable progress in the understanding of the basic mechanisms that govern the electrical conduction at the molecular scale (for a recent review see [2]).
The next step toward quantitative exploration of heat dissipation and conduction at the molecular scale has been taken very recently [28]. In that work we demonstrated that it is possible to measure the heat dissipated in the electrodes of an atomic-scale junction as a result of the passage of an electrical current (joule heating). Further, we also highlighted the novel heat dissipation effects that arise in nanoscale systems. To elaborate, while in macroscopic wires the heat dissipation is volumetric, it is not obvious a priori how the heat is dissipated in devices of nanometric size, including molecular junctions. In most of these junctions the inelastic mean free path is larger than the characteristic device dimensions and the electrical resistance is dominated by elastic processes. Thus, almost all the heat dissipation must take place inside the electrodes at a distance from the junction on the order of the inelastic scattering length. Using novel scanning tunneling probes with integrated nanothermocouples, we have been able to show that the power dissipated in the electrodes of a molecular junction is controlled by its transmission characteristics. In general, heat is not equally dissipated in both electrodes and it also depends on the bias polarity, i.e. on the direction of the current.
The goal of this work is to explore in more detail the basic principles that govern the heat dissipation in molecular junctions. For this purpose, we present here a detailed theoretical analysis of the heat dissipation in these systems based on the combination of the Landauer approach with ab initio calculations. We pay special attention to the relation between the heating asymmetries in single-molecule junctions and their thermopower. As test systems, we study junctions based on benzene derivatives that bind to gold electrodes via different anchoring groups. For these junctions, we compare the results obtained with a transport method based on density functional theory (DFT) with those obtained with the so-called DFT+ approach [29], which has been recently introduced to correct some of the known deficiencies of DFT-based approaches applied to transport problems. The study presented here provides a clear guide for controlling the heat dissipation in nanoscale systems and it outlines strategies to explore fundamental issues like the Peltier effect in single-molecule junctions.
The rest of the paper is organized as follows. In section 2 we present a detailed discussion of the Landauer theory of heat dissipation in nanoscale conductors. In particular, we discuss the basic formulae and the physical picture underlying heat dissipation in these systems. We also present some general considerations and establish the connection between heat dissipation and thermopower. Moreover, we illustrate our ideas with the help of a toy model for molecular junctions. Section 3 is devoted to a description of the two ab initio methods that we have employed, in combination with the Landauer formalism, to describe the heat dissipation and the thermopower of single-molecule junctions. In section 4 we present a detailed analysis of the heat dissipation and thermopower of single-molecule junctions based on benzene derivatives with different anchoring groups. Finally, we summarize in section 5 the main conclusions of our work.

Heat dissipation: basic formulae and physical picture
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 welldefined 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][33][34][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 8 . 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 8 Let us stress that throughout this work, we define the source electrode as the lead from which electrons flow for a positive bias voltage (and zero temperature), while the drain electrode is meant to be the lead into which electrons are injected. Notice that with our convention, for a negative bias (and zero temperature) electrons flow from the drain to the source. 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] 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 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: Obviously, the sum of the powers dissipated in both electrodes must be equal to the total power injected in the junction Q Total = I V , as can be seen from equations (4) and (5): 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 The source and drain are described by two Fermi seas with electrons occupying states up to their respective electrochemical potentials. The molecular region is described by a set of discrete state occupied up the HOMO. When a bias voltage V is applied across the junction, the electrochemical potential of the source and drain shift such that µ S − µ D = eV . (b) When an electron of energy E tunnels from the source to the drain, it leaves a hole behind. The electron releases its excess energy E − µ D in the drain, while the hole is filled up dissipating an energy equal to µ S − E in the source. (c) If the transmission is energy-dependent and, in particular, if it is higher in the lower part of the transport window, then more power is dissipated in the source for positive bias. The transmission function is represented here by a solid line, while the arrows indicate the direction of the electron flow. The red shaded areas help to visualize the different amount of heating in the electrodes. (d) Alternatively, if the transmission is higher in the upper part of the transport window, more heat is dissipated in the drain for positive bias. energy dissipated in the electrodes is equal to µ S − µ D = eV , which is exactly the work done by the external battery on each charge.
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.

General considerations
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 where f (E) = 1/[1 + exp(E/k B T )] is the Fermi function. Using the change of variable E → −E, the previous expression becomes 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 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 energyindependent 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 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, . This relation implies the following one: 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: 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 Equation (10) can also be written in the familiar form: Q S,D (V ) = ± I [38], where = T S 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 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 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: 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 2 .) 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 2 , 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 −1 [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).

Some lessons from a single-level model
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 Here, 0 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): 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 0 = +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 −3 . 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| 0 |. 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 4 ) and the bias asymmetry in the power dissipated in the source, Q S (V ) − Q S (−V ), contains a significant cubic term (∝V 3 ) and even higher order ones.  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 0 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.
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)  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 0 = +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. 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 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 (| 0 | ). Thus, the total power can be written in this limit as while the heating asymmetry adopts the form Here, G is the linear conductance, which in the off-resonant limit reads G ≈ G 0 ( / 0 ) 2 , where G 0 = 2e 2 / 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 2 ), then the power dissipated in the source can be expressed as [ Total (for negative bias), 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.

Ab initio calculations of the transmission characteristics of molecular junctions
As we have seen in the previous section, within the Landauer approach all the transport properties are determined by the transmission characteristics of the junctions. Thus, in order to describe realistic systems, we need microscopic methods to compute these characteristics.
In this section we discuss the two ab initio methods that we have employed in this work to compute the transmission through molecular junctions and, in turn, to study the corresponding heat dissipation.

Density functional theory (DFT)
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 −6 atomic units, and structure optimizations were carried out until the maximum norm of the Cartesian gradient fell below 10 −4 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 Within the Landauer approach, the low-temperature conductance is given by G = G 0 τ (E F ). It is important to remark that we shall restrict ourselves to the case of zero bias throughout this  (27))], are determined in a separate calculation [40].
discussion. The energy-dependent transmission τ (E) can be expressed in terms of the Green's functions as [1] τ where the retarded Green's function is given by and G a CC = [G r CC ] † . The self-energies in the previous equation adopt the form r On the other hand, the scattering rate matrices that enter the expression of the transmission are given by X (E) = −2 Im r X (E) , and g r X X (E) = (ES X X − H X X ) −1 are the electrode Green's functions with X = L, R.
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 X X (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.

DFT+
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 DFTbased 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: 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, H mol c m = m S mol c m , to obtain the eigenenergies, m , and the eigenvectors, 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: The potential energy, m , of a point charge distribution ρ m ( r ) placed between two perfectly conducting planes located at z = ±L/2 is given by (in atomic units) where r is the component of r = r + 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,H 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.

Heat dissipation in benzene-based single-molecule junctions
In this section we discuss the results obtained with the ab initio methods described in the previous section for the heat dissipation and its relation to thermopower in single-molecule junctions based on benzene derivatives. To be precise, we shall analyze benzene molecules attached to gold electrodes via four common anchoring groups in molecular electronics: amine (-NH 2 ); isonitrile (-NC); nitrile (-CN); and thiol (-SH). These four groups have different character (electron-donating versus electron-withdrawing) [51,52] and they will help us to illustrate that the concepts used to tune the thermopower are very similar to those that govern the heating asymmetries. Moreover, two of these molecules were investigated in our heat dissipation experiments of [28], which allows us to gauge the quality of the theoretical results, and for some of them the conductance and the thermopower have been reported experimentally, as we discuss below.
Let us start our discussion with the case of benzenediamine (BDA). We have first explored different binding geometries for this molecule between gold electrodes and found that the amine group only binds to undercoordinated Au sites, in agreement with previous studies  [28]. [29,53,54]. We show two examples of atop geometries for the Au-BDA-Au junctions in figures 6(a) and (b), which will be referred to as TT1 and TT2. These geometries are potential candidates to describe the configurations that give rise to the peaks in the experimental conductance histograms [28,53]. The corresponding results for these two geometries for the zero-bias transmission as a function of energy are shown in figures 6(c) and (d). We show the results obtained with the two methods described in the previous section, namely DFT and DFT+ . The values of the DFT conductances are equal to 2.1 × 10 −2 G 0 (TT1) and 1.1 × 10 −2 G 0 (TT2), which are slightly higher than the preferential experimental value of 6.4 × 10 −3 G 0 reported in [53] or 5.0 × 10 −3 G 0 found in [28]. As mentioned in the previous section, this typical overestimation of the conductance is attributed to exchange and correlations effects that are not properly described within the DFT framework with conventional exchangecorrelation functionals. Such functionals tend to underestimate the HOMO-LUMO gap and therefore they typically overestimate the conductance. On the other hand, the transmission curves indicate that the HOMO of the molecule dominates the charge transport in these Table 1. Relevant energies in the DFT+ method for the different molecular junctions discussed in the text. Here, −IP − H and −EA − L correspond to the energy shifts of the occupied and unoccupied states, respectively, of the molecules in gas phase as determined from the HOMO and LUMO. occ and virt correspond to the image charge energies for the occupied and unoccupied states, respectively. occ ( virt ) is the total energy shift applied to the occupied (unoccupied) states of the molecule in the junction. All the energies are in units of eV.  junctions, which results in a negative slope at the Fermi energy and, in turn, in a positive Seebeck coefficient. Using equation (15), we obtained a value for the thermopower at room temperature (300 K) of 7.2 µV K −1 for TT1 and 5.2 µV K −1 for TT2, which are similar to other DFT calculations [27]. These values have to be compared with the experimental value of 2.3 µV K −1 reported in [7]. In the case of the DFT+ method, the transmission curves exhibit a much larger HOMO-LUMO gap, as compared to DFT, and the values at the Fermi energy are around an order of magnitude smaller. As explained in section 3.2, this is due to the significant shift introduced in this method in all the occupied and empty states, see On the other hand, making use of equations (4) and (6), which account for the nonlinear contributions in the bias voltage, we have computed the power dissipated in the source electrode for these two junctions as a function of the total power dissipated in the junction and the results for positive and negative bias are shown in figures 6(e) and (f). In these panels we have also included the experimental results of [28]. It is worth stressing that in these calculations we have approximated the transmission curves by the zero-bias functions shown in figures 6(c) and (d). This approximation is justified by the fact that the highest bias explored in the experiments is still low as compared to the HOMO-LUMO gap of the molecules. Moreover, in our relatively weakly coupled and symmetric molecular junctions the voltage is expected to drop mainly at the metal-molecule interfaces. This means that the relevant orbitals in the molecule are not significantly shifted by the bias and thus, the transmission is not expected to vary appreciably with voltage. This has been shown explicitly for similar molecules, for instance, in [55,56]. Let us remind that within our approximation, the power dissipated in the drain electrode is given by Q D (V ) = Q S (−V ). As one can see in figures 6(e) and (f), both methods, DFT and DFT+ , give very similar results for the relation between the power dissipations in spite of the different results that they produce for the conductance. Notice also that the results for both junction geometries are very similar and, more importantly, all the theoretical results are in good agreement with the experimental results of [28]. Furthermore, both the theoretical and the experimental results show that for this molecule there is more heat dissipation in the source electrode for positive bias, while the situation is reversed for negative bias, where a higher power is dissipated in the drain electrode. As explained in sections 2.2 and 2.3, this behavior is due to the fact that the slope of the transmission at the Fermi energy is negative, which results in a positive Seebeck coefficient. On the other hand, the fact that the DFT and DFT+ results agree and that the relation between the power dissipations does not depend significantly on the junction geometry is a nice illustration of one of the central conclusions in section 2.3, namely that this relation is quite insensitive to the level alignment in off-resonant situations like the one realized in these Au-BDA-Au junctions.
To further illustrate the robustness of the conclusions drawn above, we have also explored how the heat dissipation evolves with the stretching of the contacts. For this purpose, we have started with the optimized geometries of figures 6(a) and (b) and we have separated step-wise the electrodes by a distance z and in every step we have re-optimized the geometries. This process simulates the junction stretching in typical break-junction experiments. In figure 7 we show both the transmission curves and the corresponding heat dissipations, calculated with the DFT-based method, for a series of geometries obtained by elongating the original TT1 and TT2 geometries by up to 1.2 Å. The important thing to remark here is that the relation between the power dissipations is hardly affected by the elongation process in the range of total powers explored in the experiments of [28]. Again, such an insensitivity can be attributed to the fact that we are in an off-resonant situation, as explained in section 2.3.
The heating asymmetries found in the Au-BDA-Au junctions can, in simple terms, be attributed to the electron-donating [51,52] character of the amine anchoring group, which results in HOMO dominated transport (or hole transport). In order to prove that the character of the terminal group determines both the thermopower and the heat dissipation, we now discuss the case of benzenediisonitrile (BDNC), since the isonitrile group has a well-known electron-withdrawing character [51,52]. Following the same analysis as for BDA, we have first investigated the most probable geometries of the Au-BDNC-Au junctions and we have found that the isonitrile group binds also preferentially to single low-coordinated gold atoms in atop positions. This is in agreement with previous studies of adsorption of isocyanides on gold surfaces [57]. Two representative examples of the atop geometries found in our analysis are shown in figures 8(a) and (b), from now on referred to as TT1 and TT2. Notice that in both cases the C atom is directly bound to a single Au atom and the main difference between the two geometries lies in the shape of the electrodes. The results for the zero-bias transmission as a function of the energy for these two junctions are shown in figures 8(c) and (d). For both geometries, and irrespective of the employed method, the low-bias transport is dominated by the LUMO. This is in qualitative agreement with previous results [58,59]. The DFT conductance values in these two examples are 4.8 × 10 −2 G 0 (TT1) and 0.14G 0 (TT2). These values are clearly higher than the preferential value of 3 × 10 −3 G 0 [60] and 2 × 10 −3 G 0 [28] found experimentally. Again, we attribute this discrepancy to the intrinsic deficiencies of the existent DFT functionals that tend to underestimate the HOMO-LUMO gap. In this case, DFT+ produces much more satisfactory results for the low-bias conductance, 6.2 × 10 −3 G 0 for TT1 and 3.0 × 10 −3 G 0 for TT2. With respect to the thermopower, as one can see in table 2, in all cases we find a negative value as a consequence of the fact that the transport is dominated by the LUMO of the molecule (electron transport). Notice also that the DFT+ method produces smaller values for the thermopower due to the smaller slope of the transmission at the Fermi energy, which is due to the fact that within this method the LUMO lies around 2 eV higher than in the DFT method, where the LUMO is located relatively close to the Fermi energy. We are not aware of measurements of the thermopower for BDNC. With respect to the heat dissipation, we show the corresponding results for BDNC in figures 8(e) and (f). Let us emphasize that these results were obtained using the zero-bias results of [28]. The main conclusions from these results are the following. Firstly, a higher power is dissipated in the source electrode for negative bias, in strong contrast to the BDA case. This is due to the fact that the slope of the transmission at the Fermi energy is positive, which leads to a negative value for the thermopower. Secondly, both methods produce similar results for the relation between the power dissipations, in spite of the big differences in their low-bias conductances. Again, this is a consequence of the insensitivity of this relation to the exact level alignment in an off-resonant situation. Thirdly, the theoretical results reproduce qualitatively the experimental findings. Fourthly, the heating asymmetries are more pronounced than in the BDA case. This is due to the larger Seebeck coefficient in the case of BDNC.
The calculations discussed so far, and their agreement with the experiments, provide strong evidence of the relationship between the sign of the thermopower and the heating asymmetries.
To provide further evidence, we have also investigated single-molecule junctions based on benzenedithiol (BDT) and benzenedinitrile (BDCN). Following the classification in terms of Hammett constants [51,52], the amine group is strongly electron-donating, thiol is weakly electron-withdrawing, and nitrile and isonitrile are both strongly electron-withdrawing. In this sense, the naive expectation is that the thermoelectricity and the heat dissipation in Au-BDT-Au and Au-BDCN-Au junctions may be similar to those in Au-BDA-Au and Au-BDNC-Au, respectively. Let us show that this is indeed the case. In figure 9 we depict the results for the transmission functions and the power dissipations for the two geometries of the upper panels. The main features of these results are the following. In the case of BDT, transport is dominated by the HOMO, the thermopower is positive, and there is higher power dissipation in the source electrode for positive bias than for negative bias. These results are indeed very similar to those obtained for the Au-BDA-Au junctions. The DFT transmission curve is similar to previous results in the literature [14,47,61]. As usual, the DFT+ method 9 [62] predicts a lower conductance value, which for this geometry is equal to 1.2 × 10 −2 G 0 , see table 2. This result is indeed very close to the experimental value of 1.1 × 10 −2 G 0 reported in [63]. With respect to the thermopower, both methods predict positive values that are close to the experimental one of 7.2µVK −1 reported in [6].
The results for the Au-BDCN-Au junction of figure 9(b) show that the transport in this case is dominated by the LUMO, which leads to a negative Seebeck coefficient and a higher power dissipation in the source electrode for negative bias, i.e. very much like in the case of Au-BDNC-Au junctions. Again, in this case the results obtained with DFT+ predict a much lower conductance than the DFT method, see table 2. In this case, the difference between these methods for the power dissipation is more significant, as one can see in figure 9(f). The difference is particularly clear at low bias, where the DFT method predicts a small negative power (cooling effect) for positive bias. This result is most likely an artifact of the position of the LUMO, which lies very close to the Fermi energy resulting in a very large and negative thermopower, see table 2. This cooling effect at relatively large total powers is absent in the DFT+ result.

Conclusions
We have presented a detailed theoretical analysis of the basic principles that govern heat dissipation in single-molecule junctions. With the help of the Landauer approach, we have shown how heat dissipation in the electrodes of a two-terminal molecular contact is determined by its transmission characteristics. In particular, we have shown how heat dissipation is, in general, different in the source and the drain electrodes and it depends on the bias polarity or current direction. We have especially emphasized the connection between these heating asymmetries and the thermopower of the junctions. As model systems, we have studied singlemolecule junctions based on benzene derivatives with different anchoring groups. Making use of DFT and DFT+ electronic structure methods, in combination with the Landauer approach, we have shown that both the heat dissipation and the thermopower can be chemically tuned by means of an appropriate choice of the terminal group. Moreover, we have shown that our results are in good agreement with our experiments of [28].
Our analysis also provides simple guidelines for future studies of related phenomena such as the Peltier effect in molecular junctions [3,64,65]. In all the junctions studied in this work, joule heating dominates over Peltier cooling over a wide range of voltages, but this can change with an appropriate choice of molecules and electrodes. For instance, resonant situations where the thermopower can take very large values can lead to cooling. It would also be interesting to explore the heat dissipation in situations where quantum interferences are expected to give rise to anomalously large values of the thermopower [16,18,65]. Another interesting possibility is to use semiconducting electrodes where the gap may help to remove hot electrons and thus to cool down one of the electrodes, as it has been demonstrated in hybrid superconducting systems in mesoscopic physics [36].