# Electron hopping heat transport in molecules

Galen T. Craven  Theoretical Division, Los Alamos National Laboratory, Los
Alamos, New Mexico 87544  Abraham Nitzan  Department of Chemistry, University
of Pennsylvania, Philadelphia, PA 19104, USA  School of Chemistry, Tel Aviv
University, Tel Aviv 69978, Israel

###### Abstract

The realization of single-molecule thermal conductance measurements has driven
the need for theoretical tools to describe conduction processes that occur
over atomistic length scales. In macroscale systems, the principle that is
typically used to understand thermal conductivity is Fourier’s law. At
molecular length scales, however, deviations from Fourier’s law are common in
part because microscale thermal transport properties typically depend on the
complex interplay between multiple heat conduction mechanisms. Here, the
thermal transport properties that arise from electron transfer across a
thermal gradient in a molecular conduction junction are examined
theoretically. We illustrate how transport in a model junction is affected by
varying the electronic structure and length of the molecular bridge in the
junction as well as the strength of the coupling between the bridge and its
surrounding environment. Three findings are of note: First, the transport
properties can vary significantly depending on the characteristics of the
molecular bridge and its environment; second, the system’s thermal conductance
commonly deviates from Fourier’s law; and third, in properly engineered
systems, the magnitude of electron hopping thermal conductance is similar to
what has been measured in single-molecule devices.

##  I Introduction

The development of experimental set-ups capable of single-molecule thermal
conductance measurements has opened the possibility for advanced understanding
of the interplay between electronic and heat transport at the molecular level.
Cui et al. (2019); Mosso et al. (2019); Reddy et al. (2007); Wang et al.
(2020); Gehring et al. (2021) The thermal conductance of molecular devices is
typically investigated in order to optimize operational stability of nano-
electronics or to control the functionality in thermal and thermoelectric
devices. Reddy et al. (2007); Wang et al. (2020); Gehring et al. (2021); Malen
et al. (2009); Ke et al. (2009); Tan et al. (2011); Dubi and Di Ventra (2011);
Li et al. (2012); Maldovan (2013); Luo and Chen (2013); Lee et al. (2013); Kim
et al. (2014); Lee et al. (2014); Leitner (2015); Segal and Agarwalla (2016);
Pandey and Leitner (2016); Cui et al. (2019); Mosso et al. (2019) One of the
major problems in examining thermal conductance at the nanoscale is that
thermal fluctuations can significantly affect the stability and functionality
of small systems, making single-molecule experimental measurements difficult
to perform. A system that is broadly used in the analysis of nanoscale thermal
and electrical conductivity is a molecular junction—a device consisting of a
molecular bridge seated between two electrodes. Aviram and Ratner (1974);
Nitzan and Ratner (2003); Ratner (2013) Molecular junctions are central in the
design of many nanoscale technologies, particularly molecular electronics, and
have therefore been studied extensively both theoretically and experimentally.
Reddy et al. (2007); Tan et al. (2011); Lee et al. (2013); Kim et al. (2014);
Capozzi et al. (2015); Garner et al. (2018); Cui et al. (2019); Mosso et al.
(2019) Applying either a temperature bias or a voltage bias, or both, across a
molecular junction induces an electric and/or heat current through the
molecular bridge which enables the control of conductivity at the molecular
level, a significant step toward the realization of technological devices at
the smallest level of miniaturization.

At the macroscale, a fundamental result that is broadly applied to analyze
heat transport is Fourier’s law:

$\mathbf{J}_{\mathcal{Q}}=-\kappa\nabla T,$ (1)

which states that the heat current density $\mathbf{J}_{\mathcal{Q}}$
generated in a system in response to a temperature gradient $\nabla T$ is
proportional to the system’s thermal conductivity $\kappa$. Dubi and Di Ventra
(2011); Segal and Agarwalla (2016) Fourier’s law is applied in the heat
transport analysis of myriad and diverse macroscopic systems. Despite its
importance, however, the question of how Fourier’s law arises from atomistic
interactions at the microscale is currently unresolved.Bonetto et al. (2000,
2004); Segal (2009); Chang et al. (2008) Fourier’s law is a limiting behavior
of atomistic dynamics and the physics underlying the way this limit is reached
under different environmental conditions is not universal. While progress has
been made toward a full understanding of the molecular origins of Fourier’s
law, the methods and models used in these studies typically rely on imposed
environmental conditions that are designed specifically to generate Fourier’s
law behavior, and thus the fundamental question of how Eq. (1) manifests at
the microscale remains open.

In purely molecular systems, phononic (vibrational) heat transport is
traditionally thought to be the sole thermal conductance mechanism.Velizhanin
et al. (2015); Dhar and Dandekar (2015); Li et al. (2012, 2006); Chang et al.
(2006); He et al. (2016); Leitner (2008, 2013); Pandey and Leitner (2016)
Recently, however, another heat transport mechanism, termed electron-transfer-
induced heat transport (ETIHT) has been identified in molecular systems.Craven
and Nitzan (2016); Cui et al. (2021) This heat transport is generated by the
transfer of electrons across a thermal gradient between molecules. Theoretical
studies have suggested that the magnitude of ETIHT can be similar to phononic
heat current and thus contribute significantly to the overall thermal
conductance in molecular systems. Chen et al. (2017)

In this article, we examine the electron-transfer-induced (ETI), i.e.,
electron hopping, thermal conductance properties of a model molecular junction
consisting of a molecular bridge of multiple charge transfer sites seated
between two metal electrodes, where each site in the bridge is characterized
by a different local temperature. This article advances our previous work on
electron hopping thermal and thermoelectric transport Craven and Nitzan
(2016); Matyushov (2016); Craven and Nitzan (2017a, b); Chen et al. (2017);
Craven et al. (2018a); Craven and Nitzan (2020) by examining how the magnitude
of ETIHT depends on the length of the molecular bridge in a molecular junction
and by performing a detailed analysis of how variations in electronic
properties and electron-phonon interactions affect ETIHT. The dominant
transport channels in molecular junctions often have quasi one-dimensional
topology Segal et al. (2003) and therefore, in this model, thermal conductance
$\mathcal{K}$ replaces thermal conductivity as the pertinent transport
property. Under linear response conditions, the total thermal conductance can
be defined as

$\mathcal{K}=\\!\lim_{\Delta T\to 0}\frac{\mathcal{J}_{\mathcal{Q}}}{\Delta
T}\bigg{|}_{\mathcal{J}_{\text{el}}=0},$ (2)

where $\mathcal{J}_{\mathcal{Q}}$ is the total heat current through the
junction calculated under zero electronic current conditions
$\mathcal{J}_{\text{el}}=0$.

In molecular junctions, a simple picture of the total thermal conductance
$\mathcal{K}$ can be constructed by taking a sum over four major
contributions:
$\mathcal{K}=\mathcal{K}_{\text{el}}+\mathcal{K}_{\text{rad}}+\mathcal{K}_{\text{ph}}+\mathcal{K}_{\text{ETI}}$,
where $\mathcal{K}_{\text{el}}$ is the electronic thermal conductance,
$\mathcal{K}_{\text{rad}}$ is the thermal conductance due to radiative heat
transfer between electrodes, and $\mathcal{K}_{\text{ph}}$ and
$\mathcal{K}_{\text{ETI}}$ are, respectively, the thermal conductances due to
phononic and ETI heat transport. What contribution each mechanisms makes to
the total thermal conductance will varying depending on the properties of the
molecular bridge, its surrounding environment, and the composition and
geometry of the electrodes. In this article, we disregard the phononic and
radiative contributions in our calculations in order to concentrate on the
other mechanisms. Although it is formally disregarded in our calculations,
comparisons will be made throughout between literature values for the phononic
conductance and the thermal conductances calculated here. The ETI contribution
to the total thermal conductance is examined by applying previously developed
theoretical formalisms.Craven and Nitzan (2016); Matyushov (2016); Craven and
Nitzan (2017a, b); Chen et al. (2017); Craven et al. (2018a); Craven and
Nitzan (2020)

The remainder of this article is organized as follows: Section II contains
details of both the model and the formalism that are applied to examine ETI
thermal transport in molecular junctions. Specifically, expressions are given
for the electric current and ETI heat current in a model molecular conduction
junction. These expressions are used to evaluate the system’s thermal
properties under different conditions. The numerical procedures used to
evaluate these currents are also discussed. Section III contains the
theoretical formalism used to examine the same transport properties generated
by ballistic electron transport. In Sec. IV, the results of this model are
presented, focusing on analyzing the thermal response trends as function of
temperature, electronic structure of the bridge, and molecular bridge length.
The thermal conductance values obtained in the hopping inelastic ETI heat
transport limit are compared with those obtained in the ballistic elastic
transport limit. Conclusions and future directions are discussed in Sec. V.

##  II Model details

We consider a physical model used previously by us in Refs. 39 and 42 that
consists of a molecular bridge seated between two metal electrodes. For
completeness of exposition in this work, we reexamine the details of this
model below.

The molecular bridge is comprised of $N$ charge transfer sites, labeled
$1,2,\ldots,N$. The electrodes are are denoted as $\text{M}_{\text{L}}$ (left
electrode) and $\text{M}_{\text{R}}$ (right electrode). Each molecular charge
transfer site $s$ in the bridge is associated with an electronic occupation
energy $E_{s}^{\prime}$, a local thermal environment at temperature $T_{s}$,
and a set of vibrational modes that are coupled to the ET process. These modes
are assumed to be in equilibrium at site’s local temperature. not  For
simplicity, herein, we consider the case of a single vibrational mode $x_{s}$
associated with each site $s$, although it is straightforward to generalize
this model to the case of multiple modes. Craven and Nitzan (2017b) This
assumption does not affect the general observations about thermal conductance
that will be made later. Each vibrational mode can be thought of as a
collective variable that parameterizes all of the vibrational motions
associated with the environment of its respective site—a picture that is
similar to the typical Marcus picture of ET. The left and right electrodes are
assumed to be in both thermal and electrochemical equilibrium. These
equilibria are characterized by respective temperatures $T_{\text{L}}=T+\Delta
T/2$ and $T_{\text{R}}=T-\Delta T/2$ and chemical potentials
$\mu_{\text{L}}=\mu-eV/2$ and $\mu_{\text{R}}=\mu+eV/2$, where $T$ is
temperature and $\mu$ is the Fermi level. We disregard any temperature
dependence of the Fermi levels. The potential bias across the device is
$\Phi=eV$. The spatial distribution of the energy levels $E^{\prime}_{s}$ in
the junction, i.e., the electronic structure of the molecular bridge, is
termed an energy landscape $\mathcal{L}$. Varying $\mathcal{L}$ by modifying
the electronic properties of the bridge alters the conduction properties of
the device. A schematic diagram of the molecular junction model is shown in
Fig. 1.

Figure 1:  Diagram of a representative molecular junction comprising three
charge transfer sites (gray circles) seated between two metal electrodes
(rectangles). Each molecular site $s$ has an electronic occupation energy
$E^{\prime}_{s}$ and is in contact with a local thermal reservoir (semi-
ellipse) with temperature $T_{s}$. The wavy gray lines represent molecule-
molecule and molecule-metal bonding. The wavy red lines denote heat currents.
The electrodes have respective temperatures $T_{\text{L}}$ and $T_{\text{R}}$.
In the diagram of each electrode, the dashed line denotes the Fermi level and
the colored region represents the corresponding Fermi-Dirac distribution. A
heat current $\mathcal{J}^{(s)}_{\mathcal{Q}}$ is associated with each thermal
reservoir. The heat currents associated with the molecule-metal interfaces are
$\mathcal{J}^{(\text{L})}_{\mathcal{Q}}$ and
$\mathcal{J}^{(\text{R})}_{\mathcal{Q}}$.

In the hopping limit of electron transport, the vibrational environment
associated with the molecular bridge is strongly coupled with the electronic
motion through the junction. In this limit, electronic conduction proceeds
through an inelastic mechanism in which the charge density of an electron
moving through the junction localizes on a particular site in the bridge. The
system can therefore be described by $N+1$ electronic states
$a\in\left\\{\text{M},1,\ldots,N\right\\}$ where each state $a\neq\text{M}$
corresponds to electron density being localized on molecular site $a$. The
state in which the electron occupies an energy level in one of the metal
electrodes is denoted by $a=\text{M}$ . We consider only singly-occupied
sites, the rationale being that a Coulomb blockade makes sites with multiple
electron occupancy too high in energy. For this reason, multiple-occupancy
states are usually disregarded when modeling molecular conduction junctions,
and we also disregard them here.

Electron hopping dominated by electron-phonon interaction that arises from
electron localization on each site can be described using the Marcus formalism
Craven and Nitzan (2016); Marcus (1956); Hush (1961); Kuznetsov and Ulstrup
(1999); Nitzan (2006) in which the energy of each state $a$ is described using
a paraboloid free energy surface:

$\displaystyle
E_{a}\left(x_{1},\ldots,x_{N}\right)=\sum\limits^{N}_{{s=1}}\frac{1}{2}k_{s}\left(x_{s}-\bar{x}^{(a)}_{s}\right)^{2}+E^{\prime}_{a},$
(3)

where $k_{s}$ is the force constant of mode $x_{s}$, $\bar{x}^{(a)}_{s}$ is a
shift in mode $x_{s}$ in response to electron localization on site $a$, and
$E_{a}^{{}^{\prime}}$ is the electronic energy origin of state $a$. The main
approximation involved in using the Marcus formalism to describe sequential
electron hopping is the assumption that thermal equilibration is fast relative
to the time between subsequent hopping events. A generalization to
multithermal situations when this approximation is not invoked was recently
been described. Kirchberg et al. (2020); Kirchberg and Nitzan  We apply the
displaced harmonic oscillator approximation and therefore take the force
constant to be independent of the electronic state of the system. In the
particular case of state $a=\text{M}$, the energy origin is
$E^{\prime}_{\text{M}}=\mu$. The reorganization energy of mode $s$ over a
specific state transition $a\to b$ is

$\lambda^{(a,b)}_{s}=\frac{1}{2}k_{s}\left(\bar{x}_{s}^{(a)}-\bar{x}_{s}^{(b)}\right)^{2},$
(4)

and the total reorganization energy of the entire system over the same
transition is

$\lambda^{(a,b)}=\sum_{s=1}^{N}\lambda^{(a,b)}_{s}.$ (5)

In all the numerical results we present (See Sec. IV), the total
reorganization energy and the reorganization of every site are taken to be the
same for each state transition. The total reorganization energy can therefore
be defined simply as $\lambda$.

###  II.1 Molecule-to-molecule electron transfer rates

The ET rate across the thermal gradient between sites in the molecular bridge
is Craven and Nitzan (2016)

$\displaystyle k_{a\to b}$
$\displaystyle=\frac{|V_{a,b}|^{2}}{\hbar}\sqrt{\frac{\pi}{k_{\text{B}}\left(T_{a}\lambda^{(a,b)}_{a}+T_{b}\lambda^{(a,b)}_{b}\right)}}$
$\displaystyle\qquad\times\exp{\left[-\frac{\left(\Delta
E_{ba}+\lambda^{(a,b)}\right)^{2}}{4k_{\text{B}}\left(T_{a}\lambda^{(a,b)}_{a}+T_{b}\lambda^{(a,b)}_{b}\right)}\right]},$
(6)

where $V_{a,b}$ is the electronic coupling between states $a$ and $b$ and
$\Delta E_{ba}=-\Delta E_{ab}=E^{\prime}_{b}-E^{\prime}_{a}$ is the reaction
free energy. Equation (II.1) is a bithermal reaction rate, previously derived
by us,Craven and Nitzan (2016); Cui et al. (2021) that describes ET in the
presence of a temperature difference between molecules. In the unithermal
limit where $T_{a}=T_{b}=T$, it reduces to the traditional Marcus ET rate
expression that describes phonon-assisted electron hopping between molecular
sites in equilibrium in the limit of strong electron-phonon coupling and high
temperature.

###  II.2 Molecule-to-metal electron transfer rates

In the high-temperature limit, the ET rates across the temperature gradients
between molecule and metal at the left and right electrodes are, Craven and
Nitzan (2017a)

$\displaystyle k_{1\to\text{M}_{\text{L}}}=\sqrt{\frac{1}{4\pi
k_{\text{B}}T_{1}\lambda^{(1,\text{M})}_{1}}}\int_{-\infty}^{\infty}\big{[}1-f\left(T_{\text{L}},\mu_{\text{L}},\epsilon\right)\big{]}$
$\displaystyle\qquad\times\Gamma_{\text{L}}(\epsilon)\exp\Bigg{[}-\frac{\big{(}\Delta
E_{\text{L}}+\epsilon+\lambda^{(1,\text{M})}_{1}\big{)}^{2}}{4k_{\text{B}}T_{1}\lambda^{(1,\text{M})}_{1}}\Bigg{]}\,d\epsilon,$
(7) $\displaystyle k_{N\to\text{M}_{\text{R}}}=\sqrt{\frac{1}{4\pi
k_{\text{B}}T_{N}\lambda^{(N,\text{M})}_{N}}}\int_{-\infty}^{\infty}\big{[}1-f\left(T_{\text{R}},\mu_{\text{R}},\epsilon\right)\big{]}$
$\displaystyle\qquad\times\Gamma_{\text{R}}(\epsilon)\exp\left[-\frac{(\Delta
E_{\text{R}}+\epsilon+\lambda^{(N,\text{M})}_{N})^{2}}{4k_{\text{B}}T_{N}\lambda^{(N,\text{M})}_{N}}\right]\,d\epsilon,$
(8)

and the corresponding rates of ET from metal to molecule are

$\displaystyle k_{\text{M}_{\text{L}}\to 1}=\sqrt{\frac{1}{4\pi
k_{\text{B}}T_{1}\lambda^{(\text{M},1)}_{1}}}\int_{-\infty}^{\infty}f\left(T_{\text{L}},\mu_{\text{L}},\epsilon\right)$
$\displaystyle\qquad\times\Gamma_{\text{L}}(\epsilon)\exp\left[-\frac{(-\Delta
E_{\text{L}}-\epsilon+\lambda^{(\text{M},1)}_{1})^{2}}{4k_{\text{B}}T_{1}\lambda^{(\text{M},1)}_{1}}\right]\,d\epsilon,$
(9) $\displaystyle k_{\text{M}_{\text{R}}\to N}=\sqrt{\frac{1}{4\pi
k_{\text{B}}T_{N}\lambda^{(\text{M},N)}_{N}}}\int_{-\infty}^{\infty}f\left(T_{\text{R}},\mu_{\text{R}},\epsilon\right)$
$\displaystyle\qquad\times\Gamma_{\text{R}}(\epsilon)\exp\left[-\frac{(-\Delta
E_{\text{R}}-\epsilon+\lambda^{(\text{M},N)}_{N})^{2}}{4k_{\text{B}}T_{N}\lambda^{(\text{M},N)}_{N}}\right]\,d\epsilon,$
(10)

where
$f(T,\mu,\epsilon)=(\exp\left[(\epsilon-\mu)/k_{\text{B}}T\right]+1)^{-1}$ is
the Fermi-Dirac distribution. The molecule-metal coupling strength at
electrode E is defined by

$\Gamma_{\text{E}}(\epsilon)=\bigg{(}\frac{2\pi}{\hbar}|V_{\text{E}}|^{2}\rho_{\text{E}}\bigg{)}_{\\!\epsilon}:\text{E}\in\left\\{\text{L},\text{R}\right\\}.$
(11)

These rates can be derived from a transition state theory approachNitzan
(2006); Craven and Hernandez (2015) using assumptions for the transmission
factor Craven and Nitzan (2017a) or from a Fermi’s golden rule formulation.
Cui et al. (2021) Note that the single electron density of states
$\rho_{\text{E}}$ and electronic tunneling coupling $V_{\text{E}}$ between
molecule and metal both formally depend on the energy level $\epsilon$,
however, in all numerical calculations we assume them to be independent of
energy using the so-called wide-band approximation. The reaction free energies
for state transitions $1\to\text{M}$ at the left electrode and $N\to\text{M}$
at the right electrode are

$\Delta E_{\text{L}}=E^{\prime}_{\text{M}}-E^{\prime}_{1},$ (12)

and

$\Delta E_{\text{R}}=E^{\prime}_{\text{M}}-E^{\prime}_{N},$ (13)

respectively, where $E^{\prime}_{\text{M}}=\mu$. In the limit in which the
temperature bias vanishes, the rate expressions in Eqs. (II.2) - (II.2) are
the Marcus-Hush-Chidsey rates for ET between a molecular and metal under
equilibrium conditions. Marcus (1965); Hush (1968); Chidsey (1991); Nitzan
(2006); Migliore and Nitzan (2011); Laborda et al. (2012)

###  II.3 Electric current

An expression for the electric current $\mathcal{J}_{\text{el}}$ through the
junction can be derived from the kinetic master equations describing the
probability $\mathcal{P}_{a}$ that the system is in each state
$a\in\left\\{\text{M},1,\ldots,N\right\\}$:

$\displaystyle\dot{\mathcal{P}}_{\text{M}}$
$\displaystyle=-k_{\text{M}_{\text{L}}\to
1}\mathcal{P}_{\text{M}}+k_{1\to\text{M}_{\text{L}}}\mathcal{P}_{1}$ (14)
$\displaystyle\qquad-k_{\text{M}_{\text{R}}\to
N}\mathcal{P}_{\text{M}}+k_{N\to\text{M}_{\text{R}}}\mathcal{P}_{N},$
$\displaystyle\dot{\mathcal{P}}_{1}$
$\displaystyle=-\left(k_{1\to\text{M}_{\text{L}}}+k_{1\to
2}\right)\mathcal{P}_{1}$ $\displaystyle\qquad+k_{{\text{M}_{\text{L}}}\to
1}\mathcal{P}_{\text{M}}+k_{2\to 1}\mathcal{P}_{2},$ $\displaystyle\ldots$
$\displaystyle\dot{\mathcal{P}}_{a}$ $\displaystyle=-\left(k_{a\to
a+1}+k_{{a\to a-1}}\right)\mathcal{P}_{a}$ $\displaystyle\qquad+k_{{a+1\to
a}}\mathcal{P}_{a+1}+k_{{a-1\to a}}\mathcal{P}_{a-1},$ $\displaystyle\ldots$
$\displaystyle\dot{\mathcal{P}}_{N}$
$\displaystyle=-\left(k_{N\to\text{M}_{\text{R}}}+k_{N\to
N-1}\right)\mathcal{P}_{N}$ $\displaystyle\qquad+k_{{\text{M}_{\text{R}}}\to
N}\mathcal{P}_{\text{M}}+k_{N-1\to N}\mathcal{P}_{N-1}.$

The physical interpretation of this system is that it describes the
probability that an electron is localized on each site in the bridge. Under
steady-state (ss) conditions, $\dot{\mathcal{P}}_{a}=0$ for each state.
Therefore, when the system is at steady-state the system of differential
equations (14) transforms into a system of algebraic equations. After
augmenting the system with the probability conservation equation

$\mathcal{P}_{\text{M}}+\sum^{N}_{a=1}\mathcal{P}_{a}=1,$ (15)

the steady-state solution of (14) can be readily obtained as
$\boldsymbol{\mathcal{P}}^{(\text{ss})}=\boldsymbol{K}^{-1}\boldsymbol{\mathcal{S}}$
with
$\boldsymbol{\mathcal{P}}^{(\text{ss})}=\\{\mathcal{P}^{(\text{ss})}_{\text{M}},\mathcal{P}^{(\text{ss})}_{1},\ldots,\mathcal{P}^{(\text{ss})}_{N}\\}$
where $\mathcal{P}^{(\text{ss})}_{a}$ is the steady-state probability of state
$a$, $\boldsymbol{K}$ is a matrix of rate coefficients extracted from the LHS
of (14), and $\boldsymbol{\mathcal{S}}$ is a column vector with element values
that arise from using Eq. (15). As stated before, we have ignored multiple
site occupancies. The steady-state electronic current can be evaluated by
taking the difference between the unidirectional flux of electrons hopping
from left to right and from right to left in the junction, for example,

$\mathcal{J}_{\text{el}}=e\left(k_{1\to\text{M}_{\text{L}}}\mathcal{P}^{(\text{ss})}_{1}-k_{\text{M}_{\text{L}}\to
1}\mathcal{P}^{(\text{ss})}_{\text{M}}\right).$ (16)

In the specific case of $N=1$, which describes a single charge transfer site
seated between two electrodes, the electric current is Migliore and Nitzan
(2011)

$\mathcal{J}_{\text{el}}=e\left(\frac{k_{1\to\text{M}_{\text{L}}}k_{\text{M}_{\text{R}}\to
1}-k_{1\to\text{M}_{\text{R}}}k_{\text{M}_{\text{L}}\to
1}}{k_{1\to\text{M}_{\text{R}}}+k_{\text{M}_{\text{R}}\to
1}+k_{1\to\text{M}_{\text{R}}}+k_{\text{M}_{\text{L}}\to 1}}\right).$ (17)

###  II.4 Electron-transfer-induced heat current

The system of equations that expresses the ETI heat current
$\mathcal{J}^{(s)}_{\mathcal{Q}}$ associated the thermal environment of each
site $s$ is (see Fig. 1)

$\displaystyle\mathcal{J}^{(1)}_{\mathcal{Q}}$
$\displaystyle=-\mathcal{J}^{(\text{L})}_{\mathcal{Q}}+k_{1\to
2}\mathcal{P}_{1}\mathcal{Q}^{(1)}_{1,2}+k_{2\to
1}\mathcal{P}_{2}\mathcal{Q}^{(1)}_{2,1},$ (18)
$\displaystyle\mathcal{J}^{(2)}_{\mathcal{Q}}$ $\displaystyle=k_{1\to
2}\mathcal{P}_{1}\mathcal{Q}^{(2)}_{1,2}+k_{2\to
1}\mathcal{P}_{2}\mathcal{Q}^{(2)}_{2,1}$ $\displaystyle\qquad+k_{2\to
3}\mathcal{P}_{2}\mathcal{Q}^{(2)}_{2,3}+k_{3\to
2}\mathcal{P}_{3}\mathcal{Q}^{(2)}_{3,2},$ $\displaystyle\ldots$
$\displaystyle\mathcal{J}^{(s)}_{\mathcal{Q}}$ $\displaystyle=k_{s-1\to
s}\mathcal{P}_{s-1}\mathcal{Q}^{(s)}_{s-1,s}+k_{s\to
s-1}\mathcal{P}_{s}\mathcal{Q}^{(s)}_{s,s-1}$ $\displaystyle\qquad+k_{s\to
s+1}\mathcal{P}_{s}\mathcal{Q}^{(s)}_{s,s+1}+k_{s+1\to
s}\mathcal{P}_{s+1}\mathcal{Q}^{(s)}_{s+1,s},$ $\displaystyle\ldots$
$\displaystyle\mathcal{J}^{(N)}_{\mathcal{Q}}$
$\displaystyle=-\mathcal{J}^{(\text{R})}_{\mathcal{Q}}+k_{N-1\to
N}\mathcal{P}_{N-1}\mathcal{Q}^{(N)}_{N-1,N}$ $\displaystyle\qquad+k_{N\to
N-1}\mathcal{P}_{N}\mathcal{Q}^{(N)}_{N,N-1},$

where $\mathcal{J}^{(\text{L})}_{\mathcal{Q}}$ and
$\mathcal{J}^{(\text{R})}_{\mathcal{Q}}$ are the molecule-to-metal ETI heat
currents at the left and right electrodes, $k_{a\to b}\mathcal{P}_{a}$ is a
unidirectional electron flux term, and Craven and Nitzan (2016)

$\mathcal{Q}^{(s)}_{a,b}=\frac{\displaystyle T_{s}\lambda^{(a,b)}_{s}\Delta
E_{ab}+\left(T_{\alpha}-T_{s}\right)\lambda^{(a,b)}_{a}\lambda^{(a,b)}_{b}}{\displaystyle
T_{a}\lambda^{(a,b)}_{a}+T_{b}\lambda^{(a,b)}_{b}},$ (19)

is the heat change in the thermal environment of site $s$ during the
electronic transition $a\to b$ with $\alpha$ being the element in
$\left\\{a,b\right\\}$ that is not equal to $s$. Note that Eq. (19) is the
heat change per electron. We are interested in the steady-state limit of the
system of equations (18) which is obtained by replacing all the
$\mathcal{P}_{a}$ terms with the corresponding $\mathcal{P}^{(\text{ss})}_{a}$
terms in each heat current expression.

The molecule-to-metal heat current expressions can be constructed by examining
how ET across the thermal gradient between the terminal molecules in the
bridge (1 or $N$) and the corresponding electrode affects heat transport.
Craven and Nitzan (2017b) Using the approximation that there is only a single
vibrational mode associated with each molecular site, ET at the left electrode
only affects mode $x_{1}$ and ET at the right electrode only affects mode
$x_{N}$. Therefore, the ETI heat current between molecule and metal at the
left electrode is

$\displaystyle\mathcal{J}^{(\text{L})}_{\mathcal{Q}}$
$\displaystyle=\int_{-\infty}^{\infty}\mathcal{P}_{1}\big{[}1-f\big{(}T_{\text{L}},\mu_{\text{L}},E_{\text{L}}(x_{1})\big{)}\big{]}\Gamma\big{(}E_{\text{L}}(x_{1})\big{)}$
(20)
$\displaystyle\quad\times\big{(}E_{\text{L}}(x_{1})-\mu_{\text{L}}\big{)}\frac{\exp\big{[}-E_{1,1}^{\ddagger}(x_{1})\big{/}k_{\text{B}}T_{1}\big{]}}{Z_{1,1}^{\ddagger}}\,dx_{1}$
$\displaystyle\quad+\int_{-\infty}^{\infty}\mathcal{P}_{\text{M}}\,f\big{(}T_{\text{L}},\mu_{\text{L}},E_{\text{L}}(x_{1})\big{)}\Gamma\big{(}E_{\text{L}}(x_{1})\big{)}$
$\displaystyle\quad\times\big{(}\mu_{\text{L}}-E_{\text{L}}(x_{1})\big{)}\frac{\exp\big{[}-E_{\text{M},1}^{\ddagger}(x_{1})\big{/}k_{\text{B}}T_{1}\big{]}}{Z_{\text{M},1}^{\ddagger}}\,dx_{1},$

and at the right electrode

$\displaystyle\mathcal{J}^{(\text{R})}_{\mathcal{Q}}$
$\displaystyle=\int_{-\infty}^{\infty}\mathcal{P}_{N}\big{[}1-f\big{(}T_{\text{R}},\mu_{\text{R}},E_{\text{R}}(x_{N})\big{)}\big{]}\Gamma_{\text{R}}\big{(}E_{\text{R}}(x_{N})\big{)}$
(21)
$\displaystyle\quad\times\big{(}E_{\text{R}}(x_{N})-\mu_{\text{R}}\big{)}\frac{\exp\big{[}-E_{N,N}^{\ddagger}(x_{N})\big{/}k_{\text{B}}T_{N}\big{]}}{Z_{N,N}^{\ddagger}}\,dx_{N}$
$\displaystyle\quad+\int_{-\infty}^{\infty}\mathcal{P}_{\text{M}}\,f\big{(}T_{\text{R}},\mu_{\text{R}},E_{\text{R}}(x_{N})\big{)}\Gamma_{\text{R}}\big{(}E_{\text{R}}(x_{N})\big{)}$
$\displaystyle\quad\times\big{(}\mu_{\text{R}}-E_{\text{R}}(x_{N})\big{)}\frac{\exp\big{[}-E_{\text{M},N}^{\ddagger}(x_{N})\big{/}k_{\text{B}}T_{N}\big{]}}{Z_{\text{M},N}^{\ddagger}}\,dx_{N},$

where

$E_{a,j}^{\ddagger}(x_{s})=\frac{1}{2}k_{j}\left(x_{s}-\bar{x}^{(a)}_{s}\right)^{2},\\\
$ (22)

is the energy of mode $x_{s}$ above the corresponding energy origin
$E^{\prime}_{a}$ when the system is in state $a$,

$\displaystyle E_{\text{L}}(x_{1})$
$\displaystyle=k_{1}\Big{(}\bar{x}^{(\text{M})}_{1}-\bar{x}^{(1)}_{1}\Big{)}x_{1}+\frac{1}{2}k_{1}\left(\bar{x}^{(1)}_{1}\right)^{2}$
$\displaystyle\quad-\frac{1}{2}k_{1}\left(\bar{x}^{(\text{M})}_{1}\right)^{2}-\Delta
E_{\text{L}},$ (23)

is the energy difference between surfaces $E_{1}$ and $E_{\text{M}}$ in mode
$x_{1}$,

$\displaystyle E_{\text{R}}(x_{N})$
$\displaystyle=k_{N}\Big{(}\bar{x}^{(\text{M})}_{N}-\bar{x}^{(N)}_{N}\Big{)}x_{N}+\frac{1}{2}k_{N}\left(\bar{x}^{(N)}_{N}\right)^{2}$
$\displaystyle\quad-\frac{1}{2}k_{N}\left(\bar{x}^{(\text{M})}_{N}\right)^{2}-\Delta
E_{\text{R}},$ (24)

is the energy difference between surfaces $E_{N}$ and $E_{\text{M}}$ in mode
$x_{N}$, and

$Z_{a,s}^{\ddagger}=\int_{-\infty}^{\infty}\exp\big{[}-E_{a,s}^{\ddagger}(x_{s})/k_{\text{B}}T_{s}\big{]}\,dx_{s}=\sqrt{\frac{2\pi
k_{\text{B}}T_{s}}{k_{s}}},$ (25)

is a normalization factor. The steady-state limit of
$\mathcal{J}^{(\text{L})}_{\mathcal{Q}}$ and
$\mathcal{J}^{(\text{R})}_{\mathcal{Q}}$ are obtained by replacing the
$\mathcal{P}_{1}$ and $\mathcal{P}_{N}$ terms with
$\mathcal{P}^{(\text{ss})}_{1}$ and $\mathcal{P}^{(\text{ss})}_{N}$,
respectively.

###  II.5  Numerical implementation and physical approximations

Several specific numerical procedures are employed to calculate and analyze
ETI thermal conductance:

(a) A variant of the self-consistent reservoir (SCR) method Bolsterli et al.
(1970); Bonetto et al. (2004); Segal (2009); Sääskilahti et al. (2013) is
applied to define the temperature profile in the molecular bridge. In this
modified SCR procedure, the bridge temperatures $T_{1},T_{2},\ldots,T_{N}$ and
the voltage bias $V$ are obtained subject to the constraints:

$\mathcal{J}^{(s)}_{\mathcal{Q}}=0\,\forall\,s\in\left\\{1,2,\ldots,N\right\\}\quad\text{and}\quad\mathcal{J}_{\text{el}}=0.$
(26)

These constraints state that the heat current into/out of the thermal bath
associated with each molecular site and also the electronic current through
the junction must vanish. This is a model that could describe an isolated
molecular bridge where the nuclear environment is intramolecular or the
situation in which there is a temperature gradient in the solvent environment
in a solvated molecular junction. The latter constraint in Eq. (26) does not
appear in the traditional SCR method that is normally applied to examine
phononic heat transport mechanisms. A multidimensional Newton-Raphson root-
finding procedure is used to satisfy (26). Using the SCR procedure ensures
that there is no net heat flow associated with the thermal environment of each
molecular site and therefore the magnitude of the heat current through the
junction $\mathcal{J}_{\mathcal{Q}}$ can be well defined as

$\mathcal{J}_{\mathcal{Q}}=\big{|}\mathcal{J}^{(\text{L})}_{\mathcal{Q}}\big{|}=\big{|}\mathcal{J}^{(\text{R})}_{\mathcal{Q}}\big{|}.$
(27)

(b) The expression for the thermal conductance on the RHS of Eq. (2) is
evaluated by calculating the heat current at a specific temperature $T$ with
$\Delta T=\Delta T_{+}=1\,\text{K}$ and again with $\Delta T=\Delta
T_{-}=-1\,\text{K}$, and then reporting the measured thermal conductance as

$\mathcal{K}=\frac{\mathcal{J}_{\mathcal{Q}}(\Delta
T_{+})-\mathcal{J}_{\mathcal{Q}}(\Delta T_{-})}{\Delta T_{+}-\Delta T_{-}}.$
(28)

All calculations are performed under zero electric current conditions.

(c) Four energy landscapes are used: $\mathcal{L}_{1}\equiv
E^{\prime}_{s}=0\,\forall\,s$, $\mathcal{L}_{2}\equiv
E^{\prime}_{s}=0.05\,\text{eV}\,\forall\,s$, $\mathcal{L}_{3}$ which is linear
ramp of energy levels from $0.05\,\text{eV}$ to $-0.05\,\text{eV}$, and
$\mathcal{L}_{4}\equiv E^{\prime}_{s}=0.1\,\text{eV}\,\forall\,s$. Note that
at least two sites are needed in a molecular bridge to construct the energy
gradient in $\mathcal{L}_{3}$.

We also apply several physical approximations in our model:

(a) The wide-band approximation is applied in all numerical calculations and
thus $\Gamma_{\text{L}}$ and $\Gamma_{\text{R}}$ are taken to be independent
of the energy level $\epsilon$. Furthermore, throughout, we take the coupling
at the left and right electrodes to be equal
$\Gamma_{\text{L}}=\Gamma_{\text{R}}=\Gamma$.

(b) We also assume that electron hopping between sites only affects the two
sites and the associated vibrational environments involved in that process.
This implies that the total reorganization energy for the transition $a\to b$
is given by $\lambda^{(a,b)}=\lambda^{(a,b)}_{a}+\lambda^{(a,b)}_{b}$. We take
the site reorganization energies to be equal
$\lambda^{(a,b)}_{a}=\lambda^{(a,b)}_{b}$. The total reorganization energy is
assumed to be the same for each state transition, and we denote this
reorganization energy as $\lambda$.

(c) We assume that the nuclear system is not transferring energy to the
electrodes directly by a phononic mechanism.

##  III Ballistic Electron Transport

The ETI heat current can be compared with the Landauer electronic heat current
to illustrate how hopping ET and ballistic electron transport respectively
manifest in the thermal transport properties of the junction. The ballistic
currents are evaluated using Landauer theory and a nonequilibrium Green’s
function approach.

In the Landauer limit of transport, electronic motion through the molecular
bridge is ballistic and electron-phonon interaction is disregarded. This is
the opposite limit of transport with respect to electron-phonon coupling
strength than the electron hopping mechanism examined in ETIHT. The ballistic
electric current through the junction can be expressed using the Landauer
formalism as Galperin et al. (2008)

$\mathcal{J}_{\text{el}}=\frac{e}{\pi\hbar}\int_{-\infty}^{\infty}\\!\\!\mathcal{T}(\epsilon)\Big{[}f(T_{\text{R}},\mu_{\text{R}},\epsilon)-f(T_{\text{L}},\mu_{\text{L}},\epsilon)\Big{]}d\epsilon,$
(29)

where $\mathcal{T}(\epsilon)$ is the electronic transmission function. The
transmission function is constructed using a nonequilibrium Green’s function
approach as

$\mathcal{T}(\epsilon)=\text{Tr}\left[\hat{\Gamma}_{\text{L}}(\epsilon)\hat{G}^{\dagger}(\epsilon)\hat{\Gamma}_{\text{R}}(\epsilon)\hat{G}(\epsilon)\right],$
(30)

where $\hat{G}$ is the Green’s function for the molecular bridge that
satisfies

$\hat{G}^{-1}(\epsilon)=\begin{pmatrix}\epsilon-E^{\prime}_{1}+\frac{1}{2}i\hbar\Gamma_{\text{L}}(\epsilon)&V_{1,2}&0&0&\ldots&0\\\
V_{2,1}&\epsilon-E^{\prime}_{2}&V_{2,3}&0&\ldots&0\\\
0&V_{3,2}&\epsilon-E^{\prime}_{3}&V_{3,4}&&\vdots\\\
0&0&V_{4,3}&\ddots&\ddots&0\\\
\vdots&\vdots&&\ddots&\epsilon-E^{\prime}_{N-1}&V_{N-1,N}\\\
0&0&\cdots&0&V_{N,N-1}&\epsilon-E^{\prime}_{N}+\frac{1}{2}i\hbar\Gamma_{\text{R}}(\epsilon)\\\
\end{pmatrix}.$ (31)

The molecule-metal coupling matrices

$\hat{\Gamma}_{\text{L}}(\epsilon)=\begin{pmatrix}\hbar\Gamma_{\text{L}}(\epsilon)&0&\cdots&0\\\
0&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&0\\\ \end{pmatrix},$
(32)

and

$\hat{\Gamma}_{\text{R}}(\epsilon)=\begin{pmatrix}0&0&\cdots&0\\\
0&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\hbar\Gamma_{\text{R}}(\epsilon)\\\ \end{pmatrix},$ (33)

define the interaction between the respective terminal molecular site ($1$ or
$N$) at each end of the bridge and the corresponding electrode
($\text{M}_{\text{L}}$ or $\text{M}_{\text{R}}$). We again apply the wide-band
approximation implying that $\Gamma_{\text{L}}$ and $\Gamma_{\text{R}}$ are
taken to be constant and equal ($\Gamma_{\text{L}}=\Gamma_{\text{R}}=\Gamma$)
in all numerical calculations. Ballistic electron transport through the
junction generates an electronic heat current that can be calculated using the
Landauer expression Cui et al. (2017)

$\mathcal{J}_{\mathcal{Q}}=\frac{1}{\pi\hbar}\int_{-\infty}^{\infty}\\!\\!\mathcal{T}(\epsilon)\left(\epsilon-\mu\right)\Big{[}f(T_{\text{R}},\mu_{\text{R}},\epsilon)-f(T_{\text{L}},\mu_{\text{L}},\epsilon)\Big{]}d\epsilon,$
(34)

which is evaluated under the condition of zero electronic current
$\mathcal{J}_{\mathcal{Q}}|_{\mathcal{J}_{\text{el}}=0}$. The zero electronic
current condition is enforced by varying the voltage bias $V$. The electronic
thermal conductance generated by ballistic electron transport is calculated
numerically by combining Eq. (34) with Eq. (28).

##  IV Results

Figure 2:  Thermal conductance $\mathcal{K}$ as a function of total
reorganization energy $\lambda$ for bridges with $N=1$ (red), $N=2$ (blue),
$N=5$ (magenta), and $N=10$ (black) sites. The electronic coupling parameters
are $\Gamma=100\,\text{ps}^{-1}$ and $V_{a,b}=0.1\,\text{eV}$ in (a) and (c)
and $\Gamma=10\,\text{ps}^{-1}$ and $V_{a,b}=0.01\,\text{eV}$ in (b). The
energy landscape in all calculations is $\mathcal{L}_{1}\equiv
E^{\prime}_{s}=0\,\forall\,s$. The temperature is $300\,\text{K}$ in (a) and
(b) and $1000\,\text{K}$ in (c). Other parameters are $\mu=0$.

The ETI thermal conductance values are shown in Fig. 2 as function of
reorganization energy $\lambda$. Results are shown for various values of $N$.
The electronic coupling parameters are $\Gamma=100\,\text{ps}^{-1}$ and
$V_{a,b}=0.1\,\text{eV}$ in Fig. 2(a) and (c) and $\Gamma=10\,\text{ps}^{-1}$
and $V_{a,b}=0.01\,\text{eV}$ in Fig. 2(b). These electronic couplings are
taken from previously estimated values. Kubas et al. (2014); Gajdos et al.
(2014); Gillet et al. (2016); Futera and Blumberger (2017); Giannini et al.
(2018) It is important to note that our choice of parameter values is
designed, in part, to probe ETHIT at the limits of validity of Marcus
kinetics. The coupling parameter values $\Gamma=100\,\text{ps}^{-1}$ and
$V_{a,b}=0.1\,\text{eV}$ can be considered a limiting case for Marcus theory,
in range of the transition from nonadiabatic Marcus hopping to strong coupling
and/or adiabatic ET mechanisms using the temperatures and reorganization
energies examined here.Jain and Subotnik (2015) At a temperature of
$300\,\text{K}$, shown in Fig. 2(a) and (b), the thermal conductance is
maximized at $\lambda\approx 0.1\,\text{eV}$. Comparing Fig. 2(a) and Fig.
2(b), it can be inferred that the shape of the thermal conductance curve with
respect to variation of $\lambda$ does not depend strongly on the electronic
coupling parameters, however its magnitude does depend on these parameters.
The observation that thermal conductance tends to zero in the $\lambda\to 0$
and $\lambda\to\infty$ limits and goes through a single maximum between those
limits agrees with the observations in Ref. 36. The cause of this behavior is
that the ETI thermal conductance curves follow the electron transfer rates,
which also exhibit this behavior. In Fig. 2(c) we use the same electronic
properties as in Fig. 2(a), but with the temperature increased to
$1000\,\text{K}$. The same general trends are observed, however, the maximum
conductance is reached at $\lambda\approx 0.3\,\text{eV}$ as opposed the
maximum at $\lambda\approx 0.1\,\text{eV}$ observed for 300K. This thermal
conductance behavior again follows the functional behavior of the electron
transfer rates, i.e., the maximum thermal conductance is obtained when the
electron transfer rates are maximized.

The magnitude of the calculated ETI thermal conductances are notable. In Fig.
2(a), the maximum conductance for $N=1$ is $\approx 175\,\text{pW/K}$, a value
that is approximately $10$ times greater than the $\approx 20\,\text{pW/K}$
value that has recently been measured in alkane molecular junctions using
single-molecule experiments. Cui et al. (2019); Mosso et al. (2019) For $N=2$,
the maximum conductance is $\approx 120\,\text{pW/K}$, also a significant
value compared to the measured single-molecule phononic thermal conductance
values. Figure 2(b) shows results with both the molecule-metal and molecule-
molecule electronic couplings reduced a factor of $10$. Making the electronic
couplings weaker significantly reduces the magnitude of the ETI thermal
conductance, however for molecular bridges of length $N=1$ and $N=2$, the
maximum thermal conductance is still of the same order to the measured single-
molecule values that are mostly attributed to phononic heat transport. It is
notable that increasing the system temperature to $1000\,\text{K}$, as shown
in Figure 2(c), changes the functional shape of the thermal conductance curves
but does not significantly change their maximum value. There are three
important implications of these predicted magnitudes: (a) the ETI thermal
conductances are in range of what is currently experimentally measurable, (b)
the ETI thermal conductance could make an appreciable contribution to the
overall thermal conductance and, (c) ETIHT could be the principle heat
transport mechanism in properly engineered systems.

Figure 3:  Thermal conductance $\mathcal{K}$ as a function of total
reorganization energy $\lambda$ for energy landscapes $\mathcal{L}_{1}\equiv
E^{\prime}_{s}=0\,\forall\,s$ (red), $\mathcal{L}_{2}\equiv
E^{\prime}_{s}=0.05\,\text{eV}\,\forall\,s$ (blue), and $\mathcal{L}_{3}$
(magenta) which is linear ramp of energy levels from $0.05\,\text{eV}$ to
$-0.05\,\text{eV}$. Parameters are (a) $N=2$, $\Gamma=100\,\text{ps}^{-1}$ and
$V_{a,b}=0.1\,\text{eV}$ and (b) $N=2$, $\Gamma=10\,\text{ps}^{-1}$ and
$V_{a,b}=0.01\,\text{eV}$. Other parameters in both panels are $\mu=0$ and
$T=300\,\text{K}$. The dashed lines in (b) correspond to the values obtained
in the ballistic electron transport limit using Landauer theory.

Figure 3 illustrates how variations in the energy landscape affect the thermal
conductance for different values of the reorganization energy. The red curves
in both Figure 3 (a) and (b) correspond to $\mathcal{L}_{1}\equiv
E^{\prime}_{s}=0\,\forall\,s$, the blue curves correspond to
$\mathcal{L}_{2}\equiv E^{\prime}_{s}=0.05\,\text{eV}\,\forall\,s$ and the
magenta curves correspond to $\mathcal{L}_{3}$ which is linear ramp of energy
levels from $0.05\,\text{eV}$ to $-0.05\,\text{eV}$. The $\mathcal{L}_{1}$
landscape represents the situation in which all of the electronic occupation
energies in the molecular bridge are aligned with the unbiased Fermi levels of
the electrodes. In Fig. 3(a) the electronic coupling parameters are
$\Gamma=100\,\text{ps}^{-1}$, $V_{a,b}=0.1\,\text{eV}$ and in Fig. 3(b) they
are $\Gamma=10\,\text{ps}^{-1}$, $V_{a,b}=0.01\,\text{eV}$. Similar trends are
observed for the different landscapes. For each $\mathcal{L}$, the thermal
conductance tends to zero for $\lambda\to 0$ and $\lambda\to\infty$ and goes
through a single maximum between these limits, following the same behavior in
the molecule-metal electron transfer rates. The $\lambda$ value at which
$\mathcal{K}$ is maximized is different for each landscape, again reflecting
the same behavior in the electron transfer rates. The dashed lines in Fig.
3(b) represent the values obtained in the ballistic electron transport limit
using Landauer theory for the respective energy landscape. The primary
observation that can be taken from these ballistic calculations is that the
ETI thermal conductance can exceed the ballistic thermal conductance in
certain regimes of electronic coupling strength. In general, however,
electron-phonon coupling that manifests through the electron transfer
reorganization energy reduces the thermal transport resulting in lower values
for the thermal conductance—the examples shown Fig. 3(b) are an exception.

Figure 4:  Thermal conductance $\mathcal{K}$ as a function of temperature $T$
for bridges with $N=1$ (red), $N=2$ (blue), $N=5$ (magenta), and $N=10$
(black) sites. Parameters are (a) $\lambda=0.2\,\text{eV}$ and (b)
$\lambda=0.5\,\text{eV}$. Other parameters are $\Gamma=100\,\text{ps}^{-1}$,
$V_{a,b}=0.1\,\text{eV}$, and $\mu=0$. The bottom panel (c) shows the results
calculated in the ballistic electron transport limit using Landauer theory
with the same parameter values. The dashed curves in (c) denote that the
calculations are performed in the ballistic limit. The energy landscape in all
calculations is $\mathcal{L}_{1}$.

The temperature dependence of the thermal conductance is shown in Fig. 4.
Results are shown for molecular bridges of various lengths in each panel. In
all cases, the ETI thermal conductance is greater the shorter the molecular
chain. It is notable that the thermal conductance does not increase linearly
and without bound as the temperature is increased. Instead, a maximum value is
reached when $T\approx\lambda/4k_{\text{B}}$ for the $\mathcal{L}_{1}$
landscape examined here. This maximum thermal conductance temperature
coincides with temperature that maximizes the electron transfer rates. The
ballistic limit results are shown in Fig. 4(c). The first point to note is
that the thermal conductance due to ballistic transport is typically greater
than the thermal conductance due to ETI transport. This implies that, in
general, the electron-phonon coupling, here expressed through the
reorganization energy, will reduce the amount of heat that is transported.
This is particularly pronounced for long chains (large values of $N$). Also
note that for ballistic transport the thermal conductance is not a
monotonically decreasing function with increasing $N$. This can be observed,
for example, by noting that the conductance values $N=2$ are greater than
$N=1$ values. This behavior is observed for ballistic transport and not for
ETI transport, an interesting distinguishing characteristic between the two
mechanisms.

Figure 5:  Thermal conductance $\mathcal{K}$ as a function of electronic
coupling $V_{a,b}$ for bridges with $\Gamma=10\,\text{ps}^{-1}$ (red),
$\Gamma=50\,\text{ps}^{-1}$ (blue), and $\Gamma=100\,\text{ps}^{-1}$
(magenta). Parameters are (a) $\lambda=0.2\,\text{eV}$ and (b)
$\lambda=0.5\,\text{eV}$. Other parameters are $N=2$, $T=300\,\text{K}$, and
$\mu=0$. The bottom panel (c) shows the results calculated in ballistic
electron transport limit using Landauer theory with the same parameter values.
The dashed curves in (c) denote that the calculations are performed in the
ballistic limit. The energy landscape in all calculations is
$\mathcal{L}_{1}$.

The dependence of the thermal conductance on electronic coupling between
molecular sites $V_{a,b}$ is shown in Fig. 5. A molecule bridge with $N=2$
sites is used in all calculations. In the limiting case of $V_{a,b}\to 0$ we
observe the expected behavior that the thermal conductance vanishes. This
simply implies the intuitive result that when the two sites in the bridge are
decoupled, thermal conductance due to electronic motion goes to zero. In each
plot, results are shown for $\Gamma=10\,\text{ps}^{-1}$ (red),
$\Gamma=50\,\text{ps}^{-1}$ (blue), and $\Gamma=100\,\text{ps}^{-1}$
(magenta). Figure 5(a) shows the results for $\lambda=0.2\,\text{eV}$. The
thermal conductance exhibits a nonlinear dependence on $V_{a,b}$, with the
results following functional forms similar to logistic curves. Figure 5(b)
shows results for $\lambda=0.5\,\text{eV}$. The thermal conductance again
exhibits nonlinear dependence on $V_{a,b}$. The corresponding ballistic
thermal conductance values are shown in Fig. 5(c). Comparing the ETI and
ballistic results, there are two important observations: (a) the ballistic
results goes through a maximum and then decreases as $V_{a,b}$ is increased,
behavior that is not observed for the ETI results and (b) the maximum
ballistic thermal conductance value is larger than the maximum ETI value.

Figure 6:  Thermal conductance $\mathcal{K}$ as a function of molecule-metal
electronic coupling strength $\Gamma$. Parameters are (a) $N=1$ (red) and
$N=2$ (blue), and $\lambda=0.2\,\text{eV}$ and (b) $\lambda=0.5\,\text{eV}$.
Other parameters are $T=300\,\text{K}$, $V_{a,b}=0.1\,\text{eV}$, and $\mu=0$.
The bottom panel (c) shows the results calculated for an in the ballistic
electron transport limit using Landauer theory with the same parameter values.
In all panels, the solid curves denote results calculated in the hopping limit
and the dashed curves denote the results calculated in the ballistic electron
transport limit using Landauer theory with the same parameter values. The
energy landscape in all calculations is $\mathcal{L}_{1}$.

The dependence of the thermal conductance on molecule-metal coupling strength
$\Gamma$ is shown in Fig. 6. Analogous to the the molecule-molecule electronic
coupling case, in the limit $\Gamma\to 0$ the thermal conductance vanishes due
to the molecular bridge being decoupled from the metal electrodes in this
limit. In Fig. 6(a) and (b), curves are shown for $N=1$ (red) and $N=2$ (blue)
sites. Note that for $N=1$, the thermal conductance does not depend on
$V_{a,b}$. In Fig. 6(a), with $\lambda=0.2\,\text{eV}$, the ETI thermal
conductance scales approximately linearly with $\Gamma$ and reaches a maximum
value of $\sim 100\,\text{pW/K}$ at the end of the examined range. The
corresponding ballistic results are shown by dashed curves. The magnitude of
the ballistic values for $\mathcal{K}$ are in general greater than the ETI
values. Figure 6(b) shows the results with the reorganization energy increased
to $\lambda=0.5\,\text{eV}$. The scaling is again approximately linear in
$\Gamma$ and reaches maximum values of $\approx 40\,\text{pW/K}$ for $N=1$ and
$\approx 20\,\text{pW/K}$ for $N=2$. The results for a bridge with $N=3$ sites
are shown in Fig. 6(c) with $V_{a,b}$ reduced by a factor 10 to
$0.01\,\text{eV}$. In this case, the thermal conductance rises sharply for
small $\Gamma$, but then shows an approximately linear increase. The
corresponding ballistic result (dashed curve) goes through a maximum where its
magnitude is approximately a factor $2$ greater than ETI thermal conductance,
but then decays with increasing $\Gamma$.

Fourier’s law predicts that heat current in molecule or molecular junction
with $N$ sites will scale as

$\mathcal{J}_{\mathcal{Q}}\propto\frac{\Delta T}{N},$ (35)

i.e, that the heat current $\mathcal{J}_{\mathcal{Q}}\propto N^{-1}$. In rest
of this Section, we examine under what conditions the ETI thermal conductance
obeys Fourier’s law and under what conditions deviations from Fourier’s law
are observed. To do so, we apply the the general scaling formula

$\mathcal{J}_{\mathcal{Q}}\propto N^{\alpha}$ (36)

and extract values of $\alpha$ for various physical conditions. When
$\alpha=-1$ Fourier’s law behavior is observed.

Figure 7:  Log-log plot showing heat current $\mathcal{J}_{\mathcal{Q}}$ as a
function of molecular length $N$. The parameters for each curve are:
$\Gamma=100\,\text{ps}^{-1}$, $V_{a,b}=0.1\,\text{eV}$, $T=300\,\text{K}$
(red), $\Gamma=100\,\text{ps}^{-1}$, $V_{a,b}=0.1\,\text{eV}$,
$T=1000\,\text{K}$ (blue), $\Gamma=10\,\text{ps}^{-1}$,
$V_{a,b}=0.1\,\text{eV}$, $T=300\,\text{K}$ (magenta) with
$\lambda=0.2\,\text{eV}$ in (a) and $\lambda=0.5\,\text{eV}$ in (b). The
energy landscape for all curves in (a) and (b) is $\mathcal{L}_{1}$. The
bottom panel (c) shows results for energy landscapes $\mathcal{L}_{2}$ and
$\mathcal{L}_{3}$, both of which are defined in the caption of Fig. 3. Other
parameters in all panels are $\mu=0$ and $\Delta T=1\,\text{K}$. The dashed
curves in each panel show the slope of $N^{-1}$ (Fourier’s law) behavior.

Figure 7 shows a log-log plot of the ETI heat current as a function of $N$ for
various parameter values. The dashed lines in each panel show representative
$N^{-1}$ scaling that would arise from Fourier’s law behavior. In Fig. 7 (a),
results are shown for $\lambda=0.2\,\text{eV}$ with various combinations of
$\Gamma$, $V_{a,b}$, and $T$. The principal observation is that, in general,
the ETI heat current does not follow Fourier’s law. Another important
observation is that, in all cases, and even for long molecules, the thermal
conductance magnitude is greater than a picowatt. In Fig. 7(b), results are
shown for $\lambda=0.5\,\text{eV}$ using the same combinations of $\Gamma$,
$V_{a,b}$, and $T$ from Fig. 7 (a). The scaling does scale as $\alpha=-1$ in
the regime $N\leq 2$ for the red and blue curves, but quickly deviates from
this trend for $N>2$. For this larger reorganization energy, the heat currents
can be sub-pW as $N$ becomes large in the linear response limit (here with
$\Delta T=1\,\text{K}$). Figure 7(c) illustrates the heat current behavior as
a function $N$ for energy landscapes $\mathcal{L}_{2}$ and $\mathcal{L}_{3}$
with various reorganization energies. Again, Fourier’s law is in general not
obeyed. This is especially pronounced as $N$ becomes large.

Figure 8:  Heat current $\mathcal{J}_{\mathcal{Q}}$ as a function of molecular
length $N$ for (a) $\Gamma=100\,\text{ps}^{-1}$, $V_{a,b}=0.1\,\text{eV}$,
$\lambda=0.2\,\text{eV}$, (b) $\Gamma=100\,\text{ps}^{-1}$,
$V_{a,b}=0.1\,\text{eV}$, $\lambda=0.5\,\text{eV}$, (c)
$\Gamma=100\,\text{ps}^{-1}$, $V_{a,b}=0.1\,\text{eV}$,
$\lambda=0.2\,\text{eV}$, $\mathcal{L}_{3}$ and (d)
$\Gamma=10\,\text{ps}^{-1}$, $V_{a,b}=0.01\,\text{eV}$,
$\lambda=0.2\,\text{eV}$. Other parameters in all panels are
$T=300\,\text{K}$, $\Delta T=1\,\text{K}$, and $\mu=0$. The red curve is the
calculated value and the black curve is a fit to the functional form
$\mathcal{J}_{\mathcal{Q}}=\mathcal{J}^{(1)}_{\mathcal{Q}}N^{\alpha}$. The
energy landscape is $\mathcal{L}_{1}$ in (a), (b), and (d) and
$\mathcal{L}_{3}$ in (c).

In order to asses validity of Fourier’s law to describe ETIHT in long
molecules (here defined as molecules with $N\geq 5$ sites), we fit the
calculated heat current data for $N\geq 5$ to a power law functional form
$\mathcal{J}_{\mathcal{Q}}=\mathcal{J}^{(1)}_{\mathcal{Q}}N^{\alpha}$ and
extracted values for $\alpha$ from these fits. The results for four
representative systems, varying $\Gamma$, $V_{a,b}$, and $\lambda$, are shown
in Fig. 8. The electronic coupling parameter values in Fig. 8 (a), (b), and
(c) are chosen to represent strong electronic coupling
($\Gamma=100\,\text{ps}^{-1}$, $V_{a,b}=0.1\,\text{eV}$) and the values in
Fig. 8 (d) are chosen to represent weak electronic coupling
($\Gamma=10\,\text{ps}^{-1}$, $V_{a,b}=0.01\,\text{eV}$). In all cases, it can
be observed that in the limit of large $N$, fitting the data to a power law
provides an acceptable approximation to the decay of the heat current with
$N$. The scaling factor $\alpha$ can be taken from these fits and compared to
the result predicted by Fourier’s law.

Figure 9:  Scatter plot showing the $\alpha$ values for various parameter sets
obtained by fitting the heat currents for molecular lengths of $N\geq 5$ to
the functional form
$\mathcal{J}_{\mathcal{Q}}=\mathcal{J}^{(1)}_{\mathcal{Q}}N^{\alpha}$. The
parameter values for each set are listed in Table 1 in the Appendix. The
dashed black line denotes behavior that follows Fourier’s law ($\alpha=-1$).

Figure 9 shows the calculated values of $\alpha$ for the fifteen parameter
sets $\left\\{\Gamma,V_{a,b},\lambda,T,\mathcal{L}\right\\}$ listed in Table 1
in the Appendix. In almost all cases, we observe $\alpha<-1$. This implies two
important points: (a) deviations from Fourier’s law ($\alpha=-1$) are common
and (b) because we observe deviations over a broad spectrum of parameter
values, this suggests that these deviations arise from the nature of the
examined physical mechanism, i.e., ETI heat transport. We have used parameter
values in our calculations that are in agreement with values that have been
either observed experimentally or calculated using electronic structure
methods and that represent typically experimental conditions. We therefore
expect that deviations from Fourier’s law will be observed in molecules and
molecular junctions in which the dominant heat transport mechanism is due to
electron hopping.

##  V Conclusions

A theoretical formalism has been developed to describe the thermal transport
properties generated by electron hopping in molecules and molecular junctions.
This heat transport mechanism has been termed electron-transfer-induced heat
transport. Using this formalism, we have illustrated how the magnitude of
electron-transfer-induced heat transport in a molecular junction is affected
by varying the electronic properties of the junction’s molecular bridge as
well as the strength of the coupling between the bridge and its surrounding
solvent environment. We have shown how varying the electronic properties of
the molecular bridge alters electron transport and, subsequently, the
magnitude of electron-transfer-induced thermal conductance. The dependence of
electron-transfer-induced thermal conductance on the length of the molecular
bridge has also been examined, and it has been shown that the heat currents
generated by electron hopping commonly deviate from Fourier’s law behavior.

Comparing the magnitude of the thermal conductance values calculated here with
recent experimental measurements in molecular junctions, we expect that
electron-transfer-induced thermal conductance can make an appreciable and
measurable contribution to the overall thermal conductance in properly
engineered systems. Experimental work in that direction is currently underway,
with possible applications being to infer values for electron transfer
reorganization energies using relationships between thermal and electrical
conductance. Kondratenko et al. (2022)

Recently, Sowa et al. have derived an expression for molecule-to-metal
electron transport rate that interpolates between the Landauer (ballistic) and
Marcus (hopping) limits of transport, and that also accounts for level
broadening. Sowa et al. (2018, 2019) The main advantage of this expression is
that it can be applied to understand electronic transport properties away from
the limiting cases. Although the validity of this expression outside of the
limiting behaviors is unknown, it offers an interesting way to access
intermediate transport regimes. However, it should be noted that in the
context of the thermal response properties examined here, this expression is
not applicable, principally because using it in heat transfer calculations
predicts an unphysically large amount of heat will be transported. This is
because using this expression gives rise to divergent heat transfer integrals.
Additionally, the formalism from Refs.68 and 69 can only describe the case of
a single charge transfer site in the molecular bridge, not the case of $N$
sites treated here, and therefore its applicability is limited in the context
of describing transport in long molecules and molecular junctions.

In the theoretical picture developed here, it is assumed that energy
relaxation into the solvent environment in response to an electron transfer
event occurs instantaneously, and therefore the solvent environment is always
in an equilibrium state. This physical picture relies on strong assumptions,
and recent work that goes beyond this picture to include nonequilibrium
effects in the solvent have been developed. Matyushov (2019); Kirchberg et al.
(2020); Kirchberg and Nitzan ; Sowa and Marcus (2021) Extensions of those
theories to include the role that soft and complex fluids Likos (2006);
Coslovich and Ikeda (2013); Craven et al. (2014); Wilding and Sollich (2014);
Singh and Hernandez (2018) play in electron transfer reactions may elucidate
how structural solvent effects, such as caging in soft matter fluids, affects
thermal transport in molecules. Work in this direction may lead to enhanced
understanding of how tuning the electronic and structural properties of a
solvent environment can enhance or suppress thermal transport in junctions for
specific technological applications. In future work, applying theoretical
formalisms that treat energy activation and energy relaxation events
separately may give more insight into the physical mechanisms giving rise to
electron-transfer-induced heat transport. Craven et al. (2018b); Craven and
Nitzan (2018)

##  VI Acknowledgments

The research of AN is supported by the U.S. National Science Foundation under
the Grant No. CHE1953701 and the University of Pennsylvania. The research of
GTC is supported by Los Alamos National Laboratory through a Laboratory
Directed Research and Development (LDRD) grant.

##  VII Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request and institutional approval.

*

##  Appendix A  Parameter Sets

Table 1: Parameter sets used in Fig. 9 set | $\lambda(\mathrm{eV})$ | $\Gamma(\mathrm{ps^{-1}})$ | $V_{a,b}(\mathrm{eV})$ | $T(\mathrm{K})$ | $\mathcal{L}$ | $\alpha$
---|---|---|---|---|---|---
1 | 0.2 | 10 | 0.01 | 300 | $\mathcal{L}_{1}$ | -1.86
2 | 0.2 | 10 | 0.01 | 1000 | $\mathcal{L}_{1}$ | -1.94
3 | 0.5 | 10 | 0.01 | 300 | $\mathcal{L}_{1}$ | -1.99
4 | 0.2 | 10 | 0.1 | 300 | $\mathcal{L}_{1}$ | -0.95
5 | 0.5 | 10 | 0.1 | 300 | $\mathcal{L}_{1}$ | -1.07
6 | 0.2 | 100 | 0.1 | 300 | $\mathcal{L}_{1}$ | -1.31
7 | 0.2 | 100 | 0.1 | 1000 | $\mathcal{L}_{1}$ | -1.47
8 | 0.2 | 100 | 0.1 | 300 | $\mathcal{L}_{2}$ | -1.20
9 | 0.2 | 100 | 0.1 | 300 | $\mathcal{L}_{3}$ | -1.22
10 | 1.0 | 100 | 0.1 | 1000 | $\mathcal{L}_{1}$ | -1.97
11 | 0.5 | 100 | 0.1 | 300 | $\mathcal{L}_{1}$ | -1.61
12 | 0.5 | 100 | 0.1 | 1000 | $\mathcal{L}_{1}$ | -1.58
13 | 0.5 | 100 | 0.1 | 300 | $\mathcal{L}_{4}$ | -1.21
14 | 0.5 | 100 | 0.1 | 1000 | $\mathcal{L}_{4}$ | -1.53
15 | 0.5 | 100 | 0.1 | 1000 | $\mathcal{L}_{2}$ | -1.46

## References

  * Cui et al. (2019) L. Cui, S. Hur, Z. A. Akbar, J. C. Klöckner, W. Jeong, F. Pauly, S.-Y. Jang, P. Reddy, and E. Meyhofer, Nature 572, 628 (2019), eprint doi:10.1038/s41586-019-1420-z.
  * Mosso et al. (2019) N. Mosso, H. Sadeghi, A. Gemma, S. Sangtarash, U. Drechsler, C. Lambert, and B. Gotsmann, Nano Letters 19, 7614 (2019), eprint doi:10.1021/acs.nanolett.9b02089.
  * Reddy et al. (2007) P. Reddy, S.-Y. Jang, R. A. Segalman, and A. Majumdar, Science 315, 1568 (2007), eprint doi:10.1126/science.1137149.
  * Wang et al. (2020) K. Wang, E. Meyhofer, and P. Reddy, Adv. Funct. Mater. 30, 1904534 (2020), eprint doi:10.1002/adfm.201904534.
  * Gehring et al. (2021) P. Gehring, J. K. Sowa, C. Hsu, J. de Bruijckere, M. van der Star, J. J. Le Roy, L. Bogani, E. M. Gauger, and H. S. van der Zant, Nature Nanotechnology 16, 426 (2021).
  * Malen et al. (2009) J. A. Malen, P. Doak, K. Baheti, T. D. Tilley, A. Majumdar, and R. A. Segalman, Nano Lett. 9, 3406 (2009).
  * Ke et al. (2009) S.-H. Ke, W. Yang, S. Curtarolo, and H. U. Baranger, Nano Lett. 9, 1011 (2009), eprint doi:10.1021/nl8031229.
  * Tan et al. (2011) A. Tan, J. Balachandran, S. Sadat, V. Gavini, B. D. Dunietz, S.-Y. Jang, and P. Reddy, J. Am. Chem. Soc. 133, 8838 (2011), eprint doi:10.1021/ja202178k.
  * Dubi and Di Ventra (2011) Y. Dubi and M. Di Ventra, Rev. Mod. Phys. 83, 131 (2011), eprint doi:10.1103/RevModPhys.83.131.
  * Li et al. (2012) N. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, and B. Li, Rev. Mod. Phys. 84, 1045 (2012), eprint doi:10.1103/RevModPhys.84.1045.
  * Maldovan (2013) M. Maldovan, Nature 503, 209 (2013), eprint doi:10.1038/nature12608.
  * Luo and Chen (2013) T. Luo and G. Chen, Phys. Chem. Chem. Phys. 15, 3389 (2013), eprint doi:10.1039/C2CP43771F.
  * Lee et al. (2013) W. Lee, K. Kim, W. Jeong, L. A. Zotti, F. Pauly, J. C. Cuevas, and P. Reddy, Nature 498, 209 (2013), eprint doi:10.1038/nature12183.
  * Kim et al. (2014) Y. Kim, W. Jeong, K. Kim, W. Lee, and P. Reddy, Nature Nanotech. 9, 881 (2014), eprint doi:10.1038/nnano.2014.209.
  * Lee et al. (2014) E.-S. Lee, S. Cho, H.-K. Lyeo, and Y.-H. Kim, Phys. Rev. Lett. 112, 136601 (2014), eprint doi:10.1103/PhysRevLett.112.136601.
  * Leitner (2015) D. M. Leitner, Adv. Phys. 64, 445 (2015), eprint doi:10.1080/00018732.2015.1109817.
  * Segal and Agarwalla (2016) D. Segal and B. K. Agarwalla, Annu. Rev. Phys. Chem. 67, 185 (2016), eprint doi:10.1146/annurev-physchem-040215-112103.
  * Pandey and Leitner (2016) H. D. Pandey and D. M. Leitner, J. Phys. Chem. Lett. 7, 5062 (2016), eprint doi:10.1021/acs.jpclett.6b02539.
  * Aviram and Ratner (1974) A. Aviram and M. A. Ratner, Chem. Phys. Lett. 29, 277 (1974), eprint doi:10.1016/0009-2614(74)85031-1.
  * Nitzan and Ratner (2003) A. Nitzan and M. A. Ratner, Science 300, 1384 (2003), eprint doi:10.1126/science.1081572.
  * Ratner (2013) M. A. Ratner, Nat. Nanotechnol. 8, 378 (2013), eprint doi:10.1038/nnano.2013.110.
  * Capozzi et al. (2015) B. Capozzi, J. Xia, O. Adak, E. J. Dell, Z.-F. Liu, J. C. Taylor, J. B. Neaton, L. M. Campos, and L. Venkataraman, Nature Nanotech. 10, 522 (2015), eprint doi:10.1038/nnano.2015.97.
  * Garner et al. (2018) M. H. Garner, H. Li, Y. Chen, T. A. Su, Z. Shangguan, D. W. Paley, T. Liu, F. Ng, H. Li, S. Xiao, et al., Nature 558, 415 (2018), eprint doi:10.1038/s41586-018-0197-9.
  * Bonetto et al. (2000) F. Bonetto, J. L. Lebowitz, and L. Rey-Bellet, _Mathematical Physics_ (Imperial College, London, 2000), pp. 128–150.
  * Bonetto et al. (2004) F. Bonetto, J. L. Lebowitz, and J. Lukkarinen, J. Stat. Phys. 116, 783 (2004), eprint doi:10.1023/B:JOSS.0000037232.14365.10.
  * Segal (2009) D. Segal, Phys. Rev. E 79, 012103 (2009), eprint doi:10.1103/PhysRevE.79.012103.
  * Chang et al. (2008) C. W. Chang, D. Okawa, H. Garcia, A. Majumdar, and A. Zettl, Phys. Rev. Lett. 101, 075903 (2008), eprint doi:10.1103/PhysRevLett.101.075903.
  * Velizhanin et al. (2015) K. A. Velizhanin, S. Sahu, C.-C. Chien, Y. Dubi, and M. Zwolak, Sci. Rep. 5 (2015), eprint doi:10.1038/srep17506.
  * Dhar and Dandekar (2015) A. Dhar and R. Dandekar, Physica A 418, 49 (2015), eprint doi:10.1016/j.physa.2014.06.002.
  * Li et al. (2006) B. Li, L. Wang, and G. Casati, Appl. Phys. Lett. 88, 143501 (2006), eprint doi:10.1063/1.2191730.
  * Chang et al. (2006) C. W. Chang, D. Okawa, A. Majumdar, and A. Zettl, Science 314, 1121 (2006), eprint doi:10.1126/science.1132898.
  * He et al. (2016) D. He, J. Thingna, J.-S. Wang, and B. Li, Phys. Rev. B 94, 155411 (2016), eprint doi:10.1103/PhysRevB.94.155411.
  * Leitner (2008) D. M. Leitner, Annu. Rev. Phys. Chem. 59, 233 (2008), eprint doi:10.1146/annurev.physchem.59.032607.093606.
  * Leitner (2013) D. M. Leitner, J. Phys. Chem. B 117, 12820 (2013), eprint doi:10.1021/jp402012z.
  * Craven and Nitzan (2016) G. T. Craven and A. Nitzan, Proc. Natl. Acad. Sci. 113, 9421 (2016), eprint doi:10.1073/pnas.1609141113.
  * Cui et al. (2021) B. Cui, G. T. Craven, and A. Nitzan, J. Chem. Phys. 155, 194104 (2021), eprint doi:10.1063/5.0068303.
  * Chen et al. (2017) R. Chen, G. T. Craven, and A. Nitzan, J. Chem. Phys. 147, 124101 (2017), eprint doi:10.1063/1.4990410.
  * Matyushov (2016) D. V. Matyushov, Proc. Natl. Acad. Sci. 113, 9401 (2016), eprint doi:10.1073/pnas.1610542113.
  * Craven and Nitzan (2017a) G. T. Craven and A. Nitzan, J. Chem. Phys. 146, 092305 (2017a), eprint doi:10.1063/1.4971293.
  * Craven and Nitzan (2017b) G. T. Craven and A. Nitzan, Phys. Rev. Lett. 118, 207201 (2017b), eprint doi:10.1103/PhysRevLett.118.207201.
  * Craven et al. (2018a) G. T. Craven, D. He, and A. Nitzan, Phys. Rev. Lett. 121, 247704 (2018a), eprint doi:10.1103/PhysRevLett.121.247704.
  * Craven and Nitzan (2020) G. T. Craven and A. Nitzan, Nano Lett. 20, 989 (2020), eprint doi:10.1021/acs.nanolett.9b04070.
  * Segal et al. (2003) D. Segal, A. Nitzan, and P. Hänggi, J. Chem. Phys. 119, 6840 (2003), eprint doi:10.1063/1.1603211.
  * (44) The molecular junction model considered here is a model we have used previously in Refs. 39 and 42.
  * Marcus (1956) R. A. Marcus, J. Chem. Phys. 24, 966 (1956), eprint doi:10.1063/1.1742723.
  * Hush (1961) N. Hush, Trans. Faraday Soc. 57, 557 (1961), eprint doi:10.1039/TF9615700557.
  * Kuznetsov and Ulstrup (1999) A. M. Kuznetsov and J. Ulstrup, _Electron Transfer in Chemistry and Biology: An Introduction to the Theory_ (John Wiley & Sons, Ltd.: Hoboken, NJ, 1999).
  * Nitzan (2006) A. Nitzan, _Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems_ (Oxford University Press: New York, 2006).
  * Kirchberg et al. (2020) H. Kirchberg, M. Thorwart, and A. Nitzan, J. Phys. Chem. Lett. 11, 1729 (2020), eprint doi:10.1021/acs.jpclett.0c00118.
  * Kirchberg and Nitzan (0) H. Kirchberg and A. Nitzan, J. Chem. Phys. 0, null (0), eprint doi:10.1063/5.0086319.
  * Craven and Hernandez (2015) G. T. Craven and R. Hernandez, Phys. Rev. Lett. 115, 148301 (2015), eprint doi:10.1103/PhysRevLett.115.148301.
  * Marcus (1965) R. A. Marcus, J. Chem. Phys. 43, 679 (1965), eprint doi:10.1063/1.1696792.
  * Hush (1968) N. Hush, Electrochim. Acta 13, 1005 (1968), eprint doi:10.1016/0013-4686(68)80032-5.
  * Chidsey (1991) C. E. D. Chidsey, Science 251, 919 (1991), eprint doi:10.1126/science.251.4996.919.
  * Migliore and Nitzan (2011) A. Migliore and A. Nitzan, ACS Nano 5, 6669 (2011), eprint doi:10.1021/nn202206e.
  * Laborda et al. (2012) E. Laborda, M. C. Henstridge, and R. G. Compton, J. Electroanal. Chem. 667, 48 (2012), ISSN 1572-6657, eprint doi:10.1016/j.jelechem.2011.12.011.
  * Bolsterli et al. (1970) M. Bolsterli, M. Rich, and W. M. Visscher, Phys. Rev. A 1, 1086 (1970), eprint doi:10.1103/PhysRevA.1.1086.
  * Sääskilahti et al. (2013) K. Sääskilahti, J. Oksanen, and J. Tulkki, Phys. Rev. E 88, 012128 (2013), eprint doi:10.1103/PhysRevE.88.012128.
  * Galperin et al. (2008) M. Galperin, A. Nitzan, and M. A. Ratner, Mol. Phys. 106, 397 (2008), eprint doi:10.1080/00268970701837784.
  * Cui et al. (2017) L. Cui, R. Miao, C. Jiang, E. Meyhofer, and P. Reddy, J. Chem. Phys. 146, 092201 (2017), eprint doi:10.1063/1.4976982.
  * Kubas et al. (2014) A. Kubas, F. Hoffmann, A. Heck, H. Oberhofer, M. Elstner, and J. Blumberger, J. Chem. Phys. 140, 104105 (2014), eprint doi:10.1063/1.4867077.
  * Gajdos et al. (2014) F. Gajdos, S. Valner, F. Hoffmann, J. Spencer, M. Breuer, A. Kubas, M. Dupuis, and J. Blumberger, J. Chem. Theory Comput. 10, 4653 (2014), eprint doi:10.1021/ct500527v.
  * Gillet et al. (2016) N. Gillet, L. Berstis, X. Wu, F. Gajdos, A. Heck, A. de la Lande, J. Blumberger, and M. Elstner, J. Chem. Theory Comput. 12, 4793 (2016), eprint doi:10.1021/acs.jctc.6b00564.
  * Futera and Blumberger (2017) Z. Futera and J. Blumberger, J. Phys. Chem. C 121, 19677 (2017), eprint doi:10.1021/acs.jpcc.7b06566.
  * Giannini et al. (2018) S. Giannini, A. Carof, and J. Blumberger, J. Phys. Chem. Lett. 9, 3116 (2018), eprint doi:10.1021/acs.jpclett.8b01112.
  * Jain and Subotnik (2015) A. Jain and J. E. Subotnik, J. Chem. Phys. 143, 134107 (2015), eprint doi:10.1063/1.4930549.
  * Kondratenko et al. (2022) K. Kondratenko, D. Guérin, X. Wallart, S. Lenfant, and D. Vuillaume, Nanoscale pp. – (2022), eprint doi:10.1039/D2NR00819J.
  * Sowa et al. (2018) J. K. Sowa, J. A. Mol, G. A. D. Briggs, and E. M. Gauger, J. Chem. Phys. 149, 154112 (2018), eprint doi:10.1063/1.5049537.
  * Sowa et al. (2019) J. K. Sowa, J. A. Mol, and E. M. Gauger, J. Phys. Chem. C 123, 4103 (2019), eprint doi:10.1021/acs.jpcc.8b12163.
  * Matyushov (2019) D. V. Matyushov, J. Chem. Phys. 150, 074504 (2019), eprint doi:10.1063/1.5082970.
  * Sowa and Marcus (2021) J. K. Sowa and R. A. Marcus, J. Chem. Phys. 154, 034110 (2021), eprint doi:10.1063/5.0034782.
  * Likos (2006) C. N. Likos, Soft Matter 2, 478 (2006), eprint doi:10.1039/B601916C.
  * Coslovich and Ikeda (2013) D. Coslovich and A. Ikeda, Soft Matter 9, 6786 (2013), eprint doi:10.1039/C3SM50368B.
  * Craven et al. (2014) G. T. Craven, A. V. Popov, and R. Hernandez, Soft Matter 10, 5350 (2014), eprint doi:10.1039/C4SM00751D.
  * Wilding and Sollich (2014) N. B. Wilding and P. Sollich, J. Chem. Phys. 141, 094903 (2014), eprint doi:10.1063/1.4894374.
  * Singh and Hernandez (2018) R. S. Singh and R. Hernandez, Chem. Phys. Lett. 708, 233 (2018), eprint doi:10.1016/j.cplett.2018.08.032.
  * Craven et al. (2018b) G. T. Craven, R. Chen, and A. Nitzan, J. Chem. Phys. 149, 104103 (2018b), eprint doi:10.1063/1.5007854.
  * Craven and Nitzan (2018) G. T. Craven and A. Nitzan, J. Chem. Phys. 148, 044101 (2018), eprint doi:10.1063/1.5007854.